空间点模式分析 (#122)

* 空间点模式分析

* rats 频率派方法比较

* 简单介绍几个 R 包的作用

* rats 贝叶斯方法结果比较

* Test with Fedora 39 and Quarto 1.4.489

* INLA fix

* disable inla

* Quarto 1.4.504

_quarto.yml

@@ -64,6 +64,7 @@ book:
         - analyze-text-data.qmd
         - analyze-survival-data.qmd
         - analyze-time-series-data.qmd
+        - analyze-point-pattern.qmd
         - analyze-spatial-data.qmd
         - analyze-areal-data.qmd
     - part: "优化建模"

analyze-point-pattern.qmd

@@ -0,0 +1,170 @@
+# 空间点模式分析 {#sec-analyze-point-pattern}
+
+```{r}
+#| echo: false
+
+source("_common.R")
+```
+
+本章以斐济地震数据集 quakes 为例介绍空间点模式数据的操作、探索和分析。
+
+```{r}
+#| message: false
+
+library(spatstat)
+library(sf)
+library(ggplot2)
+```
+
+[**spatstat**](https://github.com/spatstat/spatstat/) 是一个伞包，囊括 8 个子包，构成一套完整的空间点模式分析工具。
+
+1.  spatstat.utils 基础的辅助分析函数
+2.  spatstat.data 点模式分析用到的示例数据集
+3.  spatstat.sparse 稀疏数组
+4.  spatstat.geom 空间数据类和几何操作
+5.  spatstat.random 生成空间随机模式
+6.  spatstat.explore 空间数据的探索分析
+7.  spatstat.model 空间数据的参数建模和推理
+8.  spatstat.linnet 线性网络上的空间分析
+
+**sf** 包是一个专门用于空间矢量数据操作的 R 包。**ggplot2** 包提供的几何图层函数 `geom_sf()` 和坐标参考系图层函数 `coord_sf()` 支持可视化空间点模式数据。
+
+## 数据操作
+
+### 类型转化
+
+先对斐济地震数据 quakes 数据集做一些数据类型转化，从 data.frame 转 Simple feature 对象。
+
+```{r}
+library(sf)
+quakes_sf <- st_as_sf(quakes, coords = c("long", "lat"), crs = st_crs(4326))
+quakes_sf
+```
+
+### 坐标转化
+
+快速简单绘图，可采用图层 `geom_sf()`，它相当于统计图层 `stat_sf()` 和坐标映射图层 `coord_sf()` 的叠加，`geom_sf()` 支持点、线和多边形等数据数据对象，可以混合叠加。 `coord_sf()` 有几个重要的参数：
+
+1.  `crs`：在绘图前将各个 `geom_sf()` 图层中的**数据**映射到该坐标参考系。
+
+2.  `default_crs`：将非 sf 图层（没有携带 CRS 信息）的数据映射到该坐标参考系，默认使用 `crs` 参数的值，常用设置 `default_crs = sf::st_crs(4326)` 将非 sf 图层中的横纵坐标转化为经纬度，采用 World Geodetic System 1984 (WGS84)。
+
+3.  `datum`：经纬网线的坐标参考系，默认值 `sf::st_crs(4326)`。
+
+下图的右子图将 quakes_sf 数据集投影到坐标参考系统[EPSG:3460](https://epsg.io/3460)。
+
+```{r}
+#| label: fig-quakes-ggplot2-grid
+#| fig-cap: 斐济地震的空间分布
+#| fig-subcap: 
+#| - 坐标参考系 4326（默认）
+#| - 坐标参考系 3460
+#| fig-width: 4
+#| fig-height: 4
+#| fig-showtext: true
+#| layout-ncol: 2
+
+library(ggplot2)
+ggplot() +
+  geom_sf(
+    data = quakes_sf, aes(color = mag)
+  )
+ggplot() +
+  geom_sf(
+    data = quakes_sf, aes(color = mag)
+  ) +
+  coord_sf(crs = 3460)
+```
+
+数据集 quakes_sf 已经准备了坐标参考系统，此时，`coord_sf()` 就会采用数据集相应的坐标参考系统，即 `sf::st_crs(4326)`。上图的左子图相当于：
+
+```{r}
+#| eval: false
+
+ggplot() +
+  geom_sf(
+    data = quakes_sf, aes(color = mag)
+  ) +
+  coord_sf(
+    crs = 4326,
+    datum = sf::st_crs(4326),
+    default_crs = sf::st_crs(4326)
+  )
+```
+
+### 凸包操作
+
+```{r}
+quakes_sf <- st_transform(quakes_sf, crs = 3460)
+# 组合 POINT 构造 POLYGON
+quakes_sfp <- st_cast(st_combine(st_geometry(quakes_sf)), "POLYGON")
+# 构造 POLYGON 的凸包
+quakes_sfp_hull <- st_convex_hull(st_geometry(quakes_sfp))
+```
+
+```{r}
+#| label: fig-convex-hull
+#| fig-cap: 凸包
+#| fig-subcap: 
+#| - 凸包（base R）
+#| - 凸包（ggplot2）
+#| fig-showtext: true
+#| fig-width: 4
+#| fig-height: 4
+#| layout-ncol: 2
+#| par: true
+
+# 绘制点及其包络
+plot(st_geometry(quakes_sf))
+# 添加凸包曲线
+plot(quakes_sfp_hull, add = TRUE)
+
+ggplot() +
+  geom_sf(data = quakes_sf) +
+  geom_sf(data = quakes_sfp_hull, fill = NA) +
+  coord_sf(crs = 3460, xlim = c(569061, 3008322), ylim = c(1603260, 4665206))
+```
+
+## 数据探索
+
+### 核密度估计
+
+给定边界内的[核密度估计与绘制热力图](https://stackoverflow.com/questions/68643517)
+
+```{r}
+# spatial point pattern ppp 类型
+quakes_ppp <- spatstat.geom::as.ppp(quakes_sf)
+# 限制散点在给定的窗口边界内平滑
+spatstat.geom::Window(quakes_ppp) <- spatstat.geom::as.owin(quakes_sfp_hull)
+# 密度估计
+density_spatstat <- spatstat.explore::density.ppp(quakes_ppp, dimyx = 256)
+# 转化为 stars 对象 栅格数据
+density_stars <- stars::st_as_stars(density_spatstat)
+# 设置坐标参考系
+density_sf <- st_set_crs(st_as_sf(density_stars), 3460)
+```
+
+### 绘制热力图
+
+```{r}
+#| label: fig-kernel-heatmap
+#| fig-cap: 热力图
+#| fig-subcap: 
+#| - 核密度估计
+#| - 核密度估计（原始数据）
+#| fig-showtext: true
+#| fig-width: 4
+#| fig-height: 4
+#| layout-ncol: 2
+
+ggplot() +
+  geom_sf(data = density_sf, aes(fill = v), col = NA) +
+  scale_fill_viridis_c() +
+  geom_sf(data = st_boundary(quakes_sfp_hull))
+
+ggplot() +
+  geom_sf(data = density_sf, aes(fill = v), col = NA) +
+  scale_fill_viridis_c() +
+  geom_sf(data = st_boundary(quakes_sfp_hull)) +
+  geom_sf(data = quakes_sf, size = 1, col = "black")
+```

hierarchical-normal-models.qmd

@@ -287,7 +287,7 @@ summary(fit_lme4)
 
 ### blme
 
-下面使用 **blme** 包 [@Chung2013] ，**blme** 包基于 **lme4** 包，截距项结果比较一致，随机效应（组间方差）和残差项（组内方差）完全不同。
+下面使用 **blme** 包 [@Chung2013] ，**blme** 包基于 **lme4** 包，参数估计结果完全一致。
 
 ```{r}
 # the mode should be at the boundary of the space.
@@ -713,13 +713,13 @@ rats_gam1 <- gam(
 summary(rats_gam1)
 ```
 
-输出结果中，固定效应部分的结果和 nlme 包完全一样。
+输出结果中，固定效应部分的结果和 **nlme** 包完全一样。
 
 ```{r}
 gam.vcomp(rats_gam1, rescale = TRUE)
 ```
 
-输出结果中，依次是随机效应的截距、斜率和残差的标准差（标准偏差），和 nlme 包给出的结果非常接近。
+输出结果中，依次是随机效应的截距、斜率和残差的标准差（标准偏差），和 **nlme** 包给出的结果非常接近。
 
 **mgcv** 包还提供函数 `gamm()`，它将混合效应和固定效应分开，在拟合 LMM 模型时，它类似 **nlme** 包的函数 `lme()`。返回一个含有 lme 和 gam 两个元素的列表，前者包含随机效应的估计，后者是固定效应的估计，固定效应中可以添加样条（或样条表示的简单随机效益，比如本节前面提及的模型）。实际上，函数 `gamm()` 分别调用 **nlme** 包和 **MASS** 包来拟合 LMM 模型和 GLMM 模型。
 
@@ -1065,7 +1065,7 @@ rats_samples <- coda.samples(rats_model,
 summary(rats_samples)
 ```
 
-输出结果与 rstan 十分一致，且采样速度极快。
+输出结果与 rstan 十分一致，且采样速度极快。类似地，`alpha0 = alpha_c - xbar * beta_c` 可得 alpha0 = 242.4752 - 22 \* 6.1878 = 106.3436。
 
 ### MCMCglmm
 
@@ -1077,6 +1077,7 @@ prior1 <- list(
   R = list(V = 1, nu = 0.002),
   G = list(G1 = list(V = 1, nu = 0.002))
 )
+set.seed(20232023)
 rats_mcmc1 <- MCMCglmm::MCMCglmm(
   weight ~ days, random = ~ rats,
   data = rats_data, verbose = FALSE, prior = prior1
@@ -1094,6 +1095,7 @@ prior2 <- list(
   R = list(V = 1, nu = 0.002),
   G = list(G1 = list(V = diag(2), nu = 0.002))
 )
+set.seed(20232023)
 rats_mcmc2 <- MCMCglmm::MCMCglmm(weight ~ days,
   random = ~ us(1 + days):rats,
   data = rats_data, verbose = FALSE, prior = prior2
@@ -1103,7 +1105,7 @@ summary(rats_mcmc2)
 
 G-structure 代表随机效应部分，R-structure 代表残差效应部分，Location effects 代表固定效应部分。**MCMCglmm** 包的这套模型表示术语源自商业软件 [ASReml](https://vsni.co.uk/software/asreml) 。
 
-随机截距的方差为 123.2867，标准差为 11.1034，随机斜率的方差 0.2746，标准差为 0.524，随机截距和随机斜率的协方差 -0.7634，相关系数为 -0.131，这与 **nlme** 包结果很接近。
+随机截距的方差为 124.1327，标准差为 11.1415，随机斜率的方差 0.2783，标准差为 0.5275，随机截距和随机斜率的协方差 -0.7457，相关系数为 -0.1268，这与 **nlme** 包结果很接近。
 
 ### INLA
 
@@ -1147,7 +1149,33 @@ summary(rats_inla2)
 
 基于 rats 数据建立变截距、变斜率的分层正态模型，也是线性混合效应模型的一种特殊情况，下表给出不同方法对模型各个参数的估计及置信区间。
 
-表中给出截距 $\beta_0$ 、斜率 $\beta_1$ 、随机截距 $\sigma_0$、随机斜率 $\sigma_1$、残差 $\sigma_{\epsilon}$ 等参数的估计及 95% 的置信区间，四舍五入保留 3 位小数。
+|                | $\beta_0$ | $\beta_1$ | $\sigma_0$ | $\sigma_1$ | $\rho_{\sigma}$ | $\sigma_{\epsilon}$ |
+|-----------------|-----------|-----------|-----------|-----------|-----------|-----------|
+| nlme (REML)    | 106.568   | 6.186     | 10.743     | 0.511      | -0.159          | 6.015               |
+| lme4 (REML)    | 106.568   | 6.186     | 10.744     | 0.511      | -0.16           | 6.015               |
+| glmmTMB (REML) | 106.568   | 6.186     | 10.743     | 0.511      | -0.16           | 6.015               |
+| MASS (PQL)     | 106.568   | 6.186     | 10.495     | 0.500      | -0.15           | 6.015               |
+| spaMM (ML)     | 106.568   | 6.186     | 10.49      | 0.499      | -0.15           | 6.015               |
+| hglm           | 106.568   | 6.186     | 10.171     | 0.491      | \-              | 6.091               |
+| mgcv (REML)    | 106.568   | 6.186     | 10.311     | 0.492      | \-              | 6.069               |
+
+: 频率派方法比较 {#tbl-rats-freqentist-compare}
+
+表中给出截距 $\beta_0$ 、斜率 $\beta_1$ 、随机截距 $\sigma_0$、随机斜率 $\sigma_1$、随机截距和斜率的相关系数 $\rho_{\sigma}$、残差 $\sigma_{\epsilon}$ 等参数的估计及 95% 的置信区间，四舍五入保留 3 位小数。固定效应部分的结果完全相同，随机效应部分略有不同。
+
+|                 | $\beta_0$ | $\beta_1$ | $\sigma_0$ | $\sigma_1$ | $\rho_{\sigma}$ | $\sigma_{\epsilon}$ |
+|-----------------|-----------|-----------|-----------|-----------|-----------|-----------|
+| rstan (NUTS)    | 106.4     | 6.2       | 14.6       | 0.5        | \-              | 6.1                 |
+| cmdstanr (NUTS) | 106       | 6.19      | 14.5       | 0.513      | \-              | 6.09                |
+| brms (NUTS)     | 106.47    | 6.18      | 11.27      | 0.54       | -0.11           | 6.15                |
+| rstanarm (NUTS) | 106.575   | 6.187     | 10.194     | 0.551      | -0.0969         | 6.219               |
+| blme (REML)     | 106.568   | 6.186     | 10.787     | 0.512      | -0.160          | 5.949               |
+| rjags (Gibbs)   | 106.344   | 6.188     | 14.623     | 0.518      | \-              | 6.073               |
+| MCMCglmm (MCMC) | 106.40    | 6.19      | 11.14      | 0.53       | -0.13           | 6.18                |
+
+: 贝叶斯方法比较 {#tbl-rats-bayesian-compare}
+
+其中，**INLA** 结果的转化未完成，表格中暂缺。**rstan** 、 **cmdstanr** 和 **rjags** 未考虑随机截距和随机斜率的相关性，因此，相关系数暂缺。MCMC 是一种随机优化算法，在不同的实现中，可重复性的要求不同，设置随机数种子仅是其中的一个必要条件，故而，每次运行程序结果可能略微不同，但不影响结论。Stan 相关的 R 包输出结果中，**rstan** 保留 1 位小数，**cmdstanr** 保留 3 位有效数字，**brms** 保留 2 位小数，**rstanarm** 小数点后保留 3 位有效数字，各不相同，暂未统一处理。
 
 ## 习题 {#sec-hierarchical-models-exercises}