From 01d361111795e16e0fa3751c85bdd4ed2d46a701 Mon Sep 17 00:00:00 2001 From: XiangyunHuang Date: Tue, 14 Nov 2023 21:50:40 +0800 Subject: [PATCH 1/7] =?UTF-8?q?=E7=A9=BA=E9=97=B4=E7=82=B9=E6=A8=A1?= =?UTF-8?q?=E5=BC=8F=E5=88=86=E6=9E=90?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- _quarto.yml | 1 + analyze-point-pattern.qmd | 158 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 159 insertions(+) create mode 100644 analyze-point-pattern.qmd diff --git a/_quarto.yml b/_quarto.yml index 61d7343f..600ffe3e 100755 --- 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 new file mode 100644 index 00000000..2d9d9f30 --- /dev/null +++ b/analyze-point-pattern.qmd @@ -0,0 +1,158 @@ +# 空间点模式分析 {#sec-analyze-point-pattern} + +```{r} +#| echo: false + +source("_common.R") +``` + +本章以斐济地震数据集 quakes 为例介绍空间点模式数据的操作、探索和分析。 + +```{r} +#| message: false + +library(spatstat) +library(sf) +library(ggplot2) +``` + +## 数据操作 + +### 类型转化 + +先对斐济地震数据 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") +``` + From 57f6cfe7a0e30f84feec5fe5862920e2ab57e2b9 Mon Sep 17 00:00:00 2001 From: XiangyunHuang Date: Wed, 15 Nov 2023 11:23:27 +0800 Subject: [PATCH 2/7] =?UTF-8?q?rats=20=E9=A2=91=E7=8E=87=E6=B4=BE=E6=96=B9?= =?UTF-8?q?=E6=B3=95=E6=AF=94=E8=BE=83?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- hierarchical-normal-models.qmd | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/hierarchical-normal-models.qmd b/hierarchical-normal-models.qmd index f0d68d96..d5fdb6f9 100644 --- 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 模型。 @@ -1147,7 +1147,19 @@ 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-compare} + +表中给出截距 $\beta_0$ 、斜率 $\beta_1$ 、随机截距 $\sigma_0$、随机斜率 $\sigma_1$、随机截距和斜率的相关系数 $\rho_{\sigma}$、残差 $\sigma_{\epsilon}$ 等参数的估计及 95% 的置信区间，四舍五入保留 3 位小数。 ## 习题 {#sec-hierarchical-models-exercises} From c2154dfb71a66725677b2d7784e7379d6ed1edae Mon Sep 17 00:00:00 2001 From: XiangyunHuang Date: Wed, 15 Nov 2023 11:27:53 +0800 Subject: [PATCH 3/7] =?UTF-8?q?=E7=AE=80=E5=8D=95=E4=BB=8B=E7=BB=8D?= =?UTF-8?q?=E5=87=A0=E4=B8=AA=20R=20=E5=8C=85=E7=9A=84=E4=BD=9C=E7=94=A8?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- analyze-point-pattern.qmd | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/analyze-point-pattern.qmd b/analyze-point-pattern.qmd index 2d9d9f30..04038313 100644 --- a/analyze-point-pattern.qmd +++ b/analyze-point-pattern.qmd @@ -16,6 +16,19 @@ 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()` 支持可视化空间点模式数据。 + ## 数据操作 ### 类型转化 @@ -155,4 +168,3 @@ ggplot() + geom_sf(data = st_boundary(quakes_sfp_hull)) + geom_sf(data = quakes_sf, size = 1, col = "black") ``` - From 1617bc38671aa1b54f36b2e5062d8d6eb39798ed Mon Sep 17 00:00:00 2001 From: XiangyunHuang Date: Wed, 15 Nov 2023 14:59:56 +0800 Subject: [PATCH 4/7] =?UTF-8?q?rats=20=E8=B4=9D=E5=8F=B6=E6=96=AF=E6=96=B9?= =?UTF-8?q?=E6=B3=95=E7=BB=93=E6=9E=9C=E6=AF=94=E8=BE=83?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- hierarchical-normal-models.qmd | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/hierarchical-normal-models.qmd b/hierarchical-normal-models.qmd index d5fdb6f9..6d73dd2e 100644 --- a/hierarchical-normal-models.qmd +++ b/hierarchical-normal-models.qmd @@ -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 @@ -1157,9 +1159,23 @@ summary(rats_inla2) | 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-compare} +: 频率派方法比较 {#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} -表中给出截距 $\beta_0$ 、斜率 $\beta_1$ 、随机截距 $\sigma_0$、随机斜率 $\sigma_1$、随机截距和斜率的相关系数 $\rho_{\sigma}$、残差 $\sigma_{\epsilon}$ 等参数的估计及 95% 的置信区间，四舍五入保留 3 位小数。 +其中，**INLA** 结果的转化未完成，表格中暂缺。**rstan** 、 **cmdstanr** 和 **rjags** 未考虑随机截距和随机斜率的相关性，因此，相关系数暂缺。MCMC 是一种随机优化算法，在不同的实现中，可重复性的要求不同，设置随机数种子仅是其中的一个必要条件，故而，每次运行程序结果可能略微不同，但不影响结论。Stan 相关的 R 包输出结果中，**rstan** 保留 1 位小数，**cmdstanr** 保留 3 位有效数字，**brms** 保留 2 位小数，**rstanarm** 小数点后保留 3 位有效数字，各不相同，暂未统一处理。 ## 习题 {#sec-hierarchical-models-exercises} From 524f7a630a206a076537048d65210a0374013a2b Mon Sep 17 00:00:00 2001 From: XiangyunHuang Date: Wed, 15 Nov 2023 18:04:11 +0800 Subject: [PATCH 5/7] [docker] Quarto 1.4.504 and Fix Matrix related issues https://github.com/XiangyunHuang/data-analysis-in-action/pull/121 --- .github/workflows/docker-build-push.yaml | 2 +- fedora_rstudio_pro.Dockerfile | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/docker-build-push.yaml b/.github/workflows/docker-build-push.yaml index 3b69f884..9fef1a08 100755 --- 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/fedora_rstudio_pro.Dockerfile b/fedora_rstudio_pro.Dockerfile index 98b2aab0..ee223174 100755 --- 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'))" \ From a4c0db69c20d51464f21a426a7ac509ddd15ae9d Mon Sep 17 00:00:00 2001 From: XiangyunHuang Date: Wed, 15 Nov 2023 18:14:34 +0800 Subject: [PATCH 6/7] =?UTF-8?q?[docker]=20ggstats=20=E5=B7=B2=E7=BB=8F?= =?UTF-8?q?=E4=BB=8E=20CRAN=20=E7=A7=BB=E9=99=A4=E4=BA=86=EF=BC=8C?= =?UTF-8?q?=E5=AE=89=E8=A3=85=20GitHub=20=E7=9A=84=E5=BC=80=E5=8F=91?= =?UTF-8?q?=E7=89=88?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- DESCRIPTION | 5 +++-- desc_pkgs.txt | 1 - 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index f37cd57e..00e0724c 100755 --- 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/desc_pkgs.txt b/desc_pkgs.txt index abb57825..27cfc5e0 100755 --- 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 From 1c823a4a5be9c29a1d9618d7d56c5ce58d5f3599 Mon Sep 17 00:00:00 2001 From: Xiangyun Huang Date: Wed, 15 Nov 2023 19:58:18 +0800 Subject: [PATCH 7/7] Test with Fedora 39 and Quarto 1.4.489 (#121) * Test with Fedora 39 and Quarto 1.4.489 * INLA fix * disable inla * Quarto 1.4.504 --- .github/workflows/quarto-book-fedora.yaml | 2 +- analyze-survival-data.qmd | 31 ++++++++++++++++++++--- 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/.github/workflows/quarto-book-fedora.yaml b/.github/workflows/quarto-book-fedora.yaml index 008f2693..4a776cc3 100755 --- 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/analyze-survival-data.qmd b/analyze-survival-data.qmd index 8012b2a4..700ead70 100755 --- 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) +```