XiangyunHuang · XiangyunHuang · Nov 15, 2023 · Nov 14, 2023 · Nov 15, 2023 · Nov 15, 2023
diff --git a/.github/workflows/docker-build-push.yaml b/.github/workflows/docker-build-push.yaml
@@ -11,7 +11,7 @@ jobs:
     runs-on: ubuntu-22.04
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
-      QUARTO_VERSION: "1.4.489"
+      QUARTO_VERSION: "1.4.504"
       CMDSTAN_VERSION: "2.33.1"
       REGISTRY: ghcr.io/xiangyunhuang
     steps:

diff --git a/.github/workflows/quarto-book-fedora.yaml b/.github/workflows/quarto-book-fedora.yaml
@@ -22,7 +22,7 @@ jobs:
     env:
       CMDSTAN_VERSION: "2.33.1"
     container:
-      image: ghcr.io/xiangyunhuang/fedora-rstudio-pro:1.4.455
+      image: ghcr.io/xiangyunhuang/fedora-rstudio-pro:1.4.504
       credentials:
          username: ${{ github.repository_owner }}
          password: ${{ secrets.GITHUB_TOKEN }}

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -54,7 +54,7 @@ Imports:
     ggrepel,
     ggridges,
     ggsignif,
-    ggstats (>= 0.3.0),
+    ggstats (>= 0.5.0),
     ggstream,
     ggsurvfit,
     ggTimeSeries,
@@ -133,7 +133,8 @@ Imports:
 LinkingTo: StanHeaders (>= 2.32.2), RcppParallel (>= 5.1.0)
 Remotes:
     davidsjoberg/ggbump,
-    davidsjoberg/ggstream
+    davidsjoberg/ggstream,
+    larmarange/ggstats
 Suggests:
     quarto (>= 1.2),
     rsconnect (>= 1.1.0)

diff --git a/_quarto.yml b/_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: "优化建模"

diff --git a/analyze-point-pattern.qmd b/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")
+```
diff --git a/analyze-survival-data.qmd b/analyze-survival-data.qmd
@@ -14,7 +14,7 @@ library(ggplot2)
 library(ggfortify) # autoplot
 library(ggsurvfit) # survfit2 and ggsurvfit
 library(glmnet)    # Cox Models
-# library(VGAM)      # R >= 4.4.0
+library(VGAM)      # R >= 4.4.0
 library(INLA)
 ```
 
@@ -91,6 +91,13 @@ autoplot(aml_survival, data = aml) +
   theme_minimal()
 ```
 
+参数化的生存分析模型（参数模型，相对于非参数模型而言）
+
+```{r}
+aml_surv_reg <- survreg(Surv(time, status) ~ x, data = aml, dist = "weibull")
+summary(aml_surv_reg)
+```
+
 下面 ggsurvfit 包再次拟合模型，并展示模型结果。
 
 ### ggsurvfit
@@ -146,7 +153,7 @@ INLA 包拟合 Cox 比例风险回归模型 [@gómez-rubio2020] 采用近似贝
 ```{r}
 library(INLA)
 inla.setOption(short.summary = TRUE)
-aml_inla <- inla(inla.surv(time, status) ~ x, data = aml, family = "exponential.surv")
+aml_inla <- inla(inla.surv(time, status) ~ x, data = aml, family = "exponential.surv", num.threads = "1:1")
 summary(aml_inla)
 ```
 
@@ -156,9 +163,25 @@ Tobit (Tobin's Probit) regression 起源于计量经济学中的 Tobit 模型，
 
 -   逻辑回归，响应变量是无序的分类变量，假定服从二项、多项分布，拟合函数 `glm()` 和 `nnet::multinom()`
 -   Probit 回归，响应变量是有序的分类变量，拟合函数 `MASS::polr()`
--   Tobit 回归，响应变量是有删失/截尾的，VGAM 包依赖少，稳定，推荐使用。
+-   Tobit 回归，响应变量是有删失/截尾的，VGAM 包依赖少，稳定，推荐使用。VGAM 包括了广义线性模型
 
 ```{r}
-# library(VGAM) # Vector Generalized Linear and Additive Models
+#| eval: false
+#| echo: false
+
+library(VGAM) # Vector Generalized Linear and Additive Models
 # VGAM::vglm(family = tobit(Upper = 800)) # Tobit regression
 ```
+
+```{r}
+library(VGAM)
+with(aml, SurvS4(time, status))
+```
+
+```{r}
+#| eval: false
+#| echo: false
+
+aml_vglm <- vglm(SurvS4(time, status) ~ x, data = aml, family = cens.poisson)
+summary(aml_vglm)
+```
diff --git a/desc_pkgs.txt b/desc_pkgs.txt
@@ -40,7 +40,6 @@ R-CRAN-ggraph
 R-CRAN-ggrepel
 R-CRAN-ggridges
 R-CRAN-ggsignif
-R-CRAN-ggstats
 R-CRAN-ggsurvfit
 R-CRAN-ggTimeSeries
 R-CRAN-ggVennDiagram

diff --git a/fedora_rstudio_pro.Dockerfile b/fedora_rstudio_pro.Dockerfile
@@ -137,8 +137,8 @@ RUN dnf -y copr enable iucar/cran \
 COPY DESCRIPTION DESCRIPTION
 # For Python
 COPY requirements.txt requirements.txt
-# Setup showtext V8 and INLA
-RUN install2.r showtextdb showtext \
+# Setup Matrix showtext V8 and INLA
+RUN install2.r MatrixModels TMB glmmTMB showtextdb showtext \
   && export GITHUB_PAT=${GITHUB_PAT} \
   && export DOWNLOAD_STATIC_LIBV8=1 \
   && Rscript -e "install.packages('INLA', repos = c(getOption('repos'), INLA = 'https://inla.r-inla-download.org/R/stable'))" \

diff --git a/hierarchical-normal-models.qmd b/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}