第 17 章 使用 R 进行生存分析
阅读提示
本章给出了许多代码块 (code chunk)，默认为全局打开。点击代码块的 “Hide” 可隐藏该代码块 ；若希望全局隐藏以方便阅读，需将页面划到最上方，右侧会显示一个额外的全局控制，将其设为 “Hide All Code” 即可（默认为 “Show All Code”）。
如无法显示/隐藏，可试着将全局控制的两个选项各点一次，原因不详。
同时，为尊重原著，未删除代码块中的提示符 > ，因此直接将代码复制粘贴到 R 中并运行是会报错的，需删除该提示符。
有许多统计软件包可用于实现本书中描述的技术，其中一些已在第 1 章 1.5 节中提及。其中，R (R Core Team, 2024) 可能是使用得最广泛的，至少有三个原因。首先，已经为各种各样的统计任务开发了 R 包，这些包由作者维护和更新。其次，新统计技术的开发和发布及其在 R 中的实施之间的周期相对较短，并且 R 代码通常伴随着新统计方法的发布。第三，R 是免费软件，可供所有人使用。然而，R 并不是一个容易学习的包，有一些编程经验是一个优势。其输出也往往相当简洁，但有时很难解释。就生存数据的分析而言，本章中的材料应该有助于解决这些问题。
本章首先介绍了 R 代码和命名法，概述了数据输入和编辑的基础知识。接着，本章介绍了使用第 2 章至第 16 章描述的方法建模生存数据的 R 包。每种情况下都包含足够的细节，使读者能够快速入门。尽管对于任何特定任务通常都有许多不同的 R 包可以使用，但本章主要关注用于展示本书所涉技术的那些包。在某些情况下，这里展示的 R 代码可能不是最优雅或最简洁的，但应该相对容易理解。本章并不试图对软件进行全面评述，随着新包的开发和其他包的弃用，本章中的部分材料不可避免地会过时。
17.1 R 的介绍
R 软件可以从 Comprehensive R Archive Network (CRAN) 下载，网址为 https://www.r-project.org。 该网站链接了许多国家的站点，从这些站点可以下载适用于 Linux, Windows 或 macOS 的 R code. R 的界面相当基础，但通过使用 RStudio Desktop 可以得到增强。RStudio Desktop 是 R 的一个集成开发环境，也是开源的，可从 https://posit.co/download/rstudio-desktop/ 下载。 除了具有源代码编辑器，可以在运行 R 命令之前输入命令外，RStudio 还简化了保存和检索文件以及访问命令语法帮助的操作。当运行代码时，输出会出现在单独的窗口中，还有用于图形的窗口、工作区中的变量列表以及已使用命令的历史记录。安装 R 和 RStudio 后，可以通过打开 RStudio 来启动 R.
与 R 的交互是通过带有 “>” 提示符的命令行进行的。在后续给出的代码中，R 命令将以 “>” 为前缀，但是该符号不是命令的一部分。R 命令直接或通过 RStudio 中的源编辑器输入到 R 工作区中。通常，使用 setwd(工作目录名称) 指定可以存储 R 代码的工作目录或文件夹比较方便，并且可以使用 getwd() 检查当前目录的名称。
17.2 数据输入和编辑
用于生存分析的数据集通常具有与生存时间、删失或状态指示符以及解释变量或因子相对应的列，而行则包含每个个体的变量值。要使用 R 进行分析的数据集称为数据框 (dataframe)，这是一个矩形数组，其中列对应于变量，行对应于观测。
变量的值可直接读入，例如
其中 <- 是赋值运算符。在许多情况下，可以使用等号 =。函数 c 将函数的自变量组合成一个向量。然后 R 可以用作计算器来操纵变量的值，如
用于计算 \(y = x \log(x) + \exp(x^3/20)\)，并且可以通过在命令提示符下键入 print(y) 或简单地输入 y 来输出结果。也可在计算中使用统计函数，例如
要创建包含变量 x 和 y 的 R 数据框，请将值分配给数据框名称，如下所示
可以使用前缀 # 将注释添加到 R 代码中。这一系列计算步骤产生的输出如下所示。
> y
[1] 1.491812e+01 4.665302e+01 1.448436e+09 1.726076e+04 1.491812e+01
> # 形成 x 值的均值和标准差
> mean(x)
[1] 4.94
> sd(x)
[1] 1.690562
> # 根据 x 和 y 构造数据框
> testdata<-data.frame(x,y)
> testdata
x y
1 3.6 1.491812e+01
2 4.2 4.665302e+01
3 7.5 1.448436e+09
4 5.8 1.726076e+04
5 3.6 1.491812e+01请注意，某些命令（例如 data.frame）由两个单词组成，并用点分隔它们。变量名称还可以包含点，如 newdata.x，并且所有变量名称和命令名称都区分大小写。有大量数学和统计函数可以这种方式操作数据，包括用于条件执行语句、数据取子集、矩阵代数等命令。任何命令的帮助文档，例如 data.frame，都可以使用 help(data.frame) 获得。
17.2.1 从文件中读取和操作数据
本书的数据集均可从出版商网站 https://www.routledge.com/
获得（译者注：具体链接为这里），并且可使用 read.table 函数读入 R 并形成数据框。假设本章中使用的数据集位于当前工作目录中。
示例 1.2 中有关女性确诊为乳腺癌后的生存时间数据可输入并分配给名为 bcancer 的数据框，如
这里，引号中的文件名是完整路径名，或者如果数据文件存储在当前目录中，则只是数据文件本身的名称。选项 header=TRUE 表示变量名称已包含在数据集中，并且这些变量名称将用作数据框中的变量名称。
请注意，为了使表 1.2 提供的数据适用于基于计算机的分析，time 和 status 这两个变量表示右删失数据，其中如果观测时间是事件，则 status=1，否则 status=0。可以使用 print 函数检查数据框，如 print(bcancer) ，或者只需在提示符处键入数据框名称，如
在前五个观测后，输出被截断。另一个有用的命令是 summary，它给出数据框中变量的描述性统计信息，如
> summary(bcancer)
stain time status
Min. :1.000 Min. : 5.00 Min. :0.0000
1st Qu.:1.000 1st Qu.: 40.00 1st Qu.:0.0000
Median :2.000 Median : 71.00 Median :1.0000
Mean :1.711 Mean : 96.24 Mean :0.5778
3rd Qu.:2.000 3rd Qu.:148.00 3rd Qu.:1.0000
Max. :2.000 Max. :225.00 Max. :1.0000 经常会发生具有相同名称的变量在不同数据框中使用的情况。为了引用特定数据框中的特定变量，使用 $ 符号来链接数据框名称和变量名称。例如，要通过将变量 stain 的值从 1, 2 转换为 0, 1 来形成新的变量 stain1，请使用
操作数据集时通常需要的另一个过程是创建因子的工具。在乳腺癌数据集中，stain 是值为 1, 2 的变量。要将其设置为具有两个水平的因子，并分别将它们标记为“阴性”和“阳性”，请使用
在该代码中，因子标签的分配顺序与数据中的顺序相同，因此当 bcancer$stain=1 时，因子水平被标记为 “Negative”.
要退出 R 会话，请使用 quit，这样您可以选择保存当前工作空间的内容以供将来使用。但是，要将工作区专门保存在 RData 类型的文件中，并在以后进行恢复，可以使用以下代码
在许多关于 R 的教科书中描述了可以进行的进一步操作，其中一些在本章末尾的延伸阅读部分提到。
17.2.2 R 包
R 程序具有许多操作、分析和展示数据的功能，其中一些功能已在本节中概述。然而，R 的强大之处是通过使用包 (packages) 来实现的。包由代码、示例数据和文档组成，可从 Comprehensive R Archive Network (CRAN) 下载。通常，一个包包含了许多用于包内特定任务的函数。现已有数千个包，涵盖统计技术、数据可视化、数据操作、机器学习、报告生成等。当需要用于某些特定目的时，可以通过首先安装包将包合并到R中。然后将已安装的包加载到当前工作区中，以使它们能够在 R 会话中使用。
对于生存分析来说，一个特别重要的包是 survival 包。 CRAN 网站上包含完整的文档，作者提供了描述该软件包 Therneau (2024) 附带的 vignettes. 该软件包实现了多种生存分析方法，包括描述性技术、Cox 回归建模、参数建模、时依变量建模、脆弱建模和模型检查技术。
base R 软件包含许多广泛使用的软件包，包括 survival，但是可以使用以下命令安装和加载该软件包以及任何其他软件包
可以使用 library() 检查库中的包列表。许多软件包需要其他软件包才能正常工作，这些软件包将与主软件包一起下载。
在以下各节中，按照本书的章节顺序，给出了用于实现生存数据建模特定技术的 R 代码。
17.3 非参数程序
分析生存数据的第一步通常是绘制两组或多组数据的生存函数的 Kaplan-Meier 估计。要在 R 中执行此操作，请使用 survival 包中的 survfit 函数。以下代码可用于估计乳腺癌预后数据中两组女性的生存函数
survfit 函数包含一个 survival 对象 Surv，用于描述模型公式中的响应变量。它的参数是生存时间和删失指示符，在本例中是 time 和 status。然后使用符号 ~ 将其链接到模型公式。如果模型中没有项，则模型公式可能仅包含 1，否则它是一串用符号 + 连接的变量和因子名称。还可以使用冒号 : 来连接项以构成交互作用项，以便 x 中变量的系数取决于 a 中因子水平，这由 a:x 表示。
在前面的代码中，指示了对两个染色组 Kaplan-Meier 估计的请求，但忽略组别的总体估计将使用 Surv(time,status) ~ 1 获得。summary 函数给出了每个 stain 值的生存函数估计的数值总结，该值可以是变量也可以是因子。stain=1 的 summary 输出如下
> km<-survfit(Surv(time,status) ~ stain,data=bcancer)
> summary(km)
Call: survfit(formula = Surv(time, status) ~ stain, data = bcancer)
stain=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
23 13 1 0.923 0.0739 0.789 1.000
47 12 1 0.846 0.1001 0.671 1.000
69 11 1 0.769 0.1169 0.571 1.000
148 6 1 0.641 0.1522 0.402 1.000
181 5 1 0.513 0.1673 0.271 0.972在这段代码中，survfit 的结果被保存到 survfit 对象 km 中。该对象包含可以提取或用于补充命令的信息。例如，使用带有 km 的 plot 函数可获得两个生存函数的图形，如
并得出与第 2 章图 2.9 所示非常相似的图。plot 函数有许多选项可用于自定义 CRAN 文档中描述的图形。还有其他的绘图函数，包括 lattice 包中的 xyplot 和 ggplot2 包。
图 2.9
survival 包中的 R 函数 surfdiff 可用于执行第 2 章 2.6 节描述的 log-rank 和 Peto-Peto 检验。该函数使用的检验统计量是 \(d_{1j}\) 和 \(e_{1j}\)（在 \(r\) 个事件时间的第 \(j\) 个事件时间，一组中观测的和预期的死亡人数）之差的加权和，由下式给出
\[\begin{aligned}U_T&=\sum_{j=1}^rw_j(d_{1j}-e_{1j})\end{aligned}\]
权重为 \(w_j=[\hat{S}_{KM}(t_{(j)})]\rho\)，其中 \(\hat{S}_{KM}(t_{(j)})\) 为在假设两组生存函数之间没有差异的情况下，在第 \(j\) 个有序事件时间 \(t_{(j)},j=1,2,\ldots,r\) 处的总体生存函数。当 \(\rho=0\) 时（这是 \(\rho\) 的默认值），得到 log-rank 检验，\(\rho=1\) 给出了 Peto-Peto 检验。
对于乳腺癌预后数据，其中数据包含在 R 数据框 bcancer 中，可以利用 survediff 使用以下代码进行 log-rank 检验
使用 print(lr) 得到以下输出，这与第 2 章示例 2.12 中的计算结果一致
Call:
survdiff(formula = Surv(time, status) ~ stain, data = bcancer)
N Observed Expected (O-E)^2/E (O-E)^2/V
stain=1 13 5 9.57 2.18 3.51
stain=2 32 21 16.43 1.27 3.51
Chisq= 3.5 on 1 degrees of freedom, p= 0.06 输出给出了每个染色组观测的和预期的事件数，最后一列中的检验统计量标题为 (O-E)^2/V。无论在计算中使用哪一组，检验统计量都是相同的。在这个例子中，默认的小数位数有时是不够的。输出中的位数可以通过使用 print 函数的 digits= 选项来更改，如 print(lr,digits=5)。结果输出的卡方值为 3.515（\(P=0.0608\)）。
Wilcoxon 检验未在 survival 包中实现，但由于结果通常与 Peto-Peto 检验的结果非常相似，因此可以使用它来代替。为此，我们在 survdiff 中包含选项 rho=1。进行 Peto-Peto 检验所需的代码和输出如下
> pp<-survdiff(Surv(time,status) ~ stain,rho=1,data=bcancer)
> print(pp)
Call:
survdiff(formula = Surv(time, status) ~ stain, data = bcancer,
rho = 1)
N Observed Expected (O-E)^2/E (O-E)^2/V
stain=1 13 3.0 6.6 1.97 4.11
stain=2 32 15.7 12.1 1.08 4.11
Chisq= 4.1 on 1 degrees of freedom, p= 0.04 结果与示例 2.13 中的计算一致。
如果特别需要 Wilcoxon 检验，则在首次安装 PHInfiniteEstimates 包后，可以按以下方式使用 PHInfiniteEstimates 包中的 R 函数 gehan.wilcoxon.test
> library("PHInfiniteEstimates")
> gehan.wilcoxon.test(Surv(time,status) ~ stain,data=bcancer)
Gehan-Wilcoxon
data:
= 4.18, p-value = 0.0409
alternative hypothesis: two-sided再次与示例 2.13 中的计算一致。
这些检验的 R 代码可以扩展为包含多个因子，并且可以通过在模型公式中包含 strata 项来实现第 2 章 2.8 节描述的分层检验。
17.4 Cox 回归模型
给出第 \(i\) 个个体（共 \(n\) 个）在时间 \(t\) 发生事件风险的 Cox 回归模型为
\[\begin{aligned}h_i(t)=\exp(\boldsymbol{\beta'x}_i)h_0(t)\end{aligned}\]
其中 \(\boldsymbol{\beta}^{\prime}\boldsymbol{x}_i\) 是第 \(i(i=1,2,\ldots,n)\) 个个体解释变量值的线性组合。 \(h_0(t)\) 是未指定的基线风险函数。该模型可以使用 survival 包中的 R 函数 coxph 进行拟合。指定模型的语法与函数 surfit 中使用的语法相同，因为 Surv 对象链接到模型公式。
如何处理结生存时间是一个重要的问题，这在第 3 章 3.3.2 节讨论过。许多 Cox 回归分析软件包使用 Breslow 法作为默认方法，主要是因为它在计算上很简单。然而，R生存函数 coxph 默认使用 Efron 法，因为该法更准确。这意味着 R 的默认输出可能与其他软件包的输出不匹配。为了更改 coxph 中处理结的方法，包含选项 ties="breslow"，如下所示。输出将与本书中使用 Breslow 法的示例相对应。
为了说明 coxph 函数的使用，将使用第 1 章示例 1.3 中给出的 48 名多发性骨髓瘤患者的生存时间数据。下面给出了用于输入数据并拟合包含变量 \(Hb\) 和 \(Bun\) 的 Cox 回归模型的 R 代码。
> myeloma<-read.table("Survival of multiple myeloma patients.dat", header=TRUE)
> coxregfit<-coxph(Surv(time,status) ~ hb+bun,ties="breslow", data=myeloma)使用 summary(coxregfit) 将得到如下输出
Call:
coxph(formula = Surv(time, status) ~ hb + bun, data = myeloma,
ties = "breslow")
n= 48, number of events= 36
coef exp(coef) se(coef) z Pr(>|z|)
hb -0.133639 0.874906 0.062051 -2.154 0.03126 *
bun 0.018571 1.018745 0.005674 3.273 0.00106 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
hb 0.8749 1.1430 0.7747 0.988
bun 1.0187 0.9816 1.0075 1.030
Concordance= 0.672 (se = 0.053 )
Likelihood ratio test= 13 on 2 df, p=0.002
Wald test = 14.86 on 2 df, p=6e-04
Score (logrank) test = 17.76 on 2 df, p=1e-04拟合的 Cox 回归模型中变量 \(Hb\) 和 \(Bun\) 的系数估计分别为 -0.1336 和 0.0186，标准误分别为 0.0620 和 0.0057。输出中包括了基于统计量 \(z=\hat{\beta}/\text{se}(\hat{\beta})\) 和相应 \(P\) 值的 Wald 检验。还给出了风险比 \(\exp(\hat\beta)\)，在标题为 lower .95 和 upper .95 的列中还包括了相应的 95% 置信限。然后是一致性度量，即第 3 章3.12.1 节中描述的 c 统计量，以及拟合模型中所有系数均为零的假设的似然比、Wald 和 得分检验结果。这些检验在附录 A 的 A.2 节中进行了定义，但很少使用 Wald 和得分检验。
coxph 的输出不包括 \(−2\log\hat L\) 统计量的值，但这可以使用从 coxph 对象的 coxregfit 中提取
该输出的两个值分别是没有解释变量的模型（即零模型）和拟合模型的统计量的值，这些与第 3 章表 3.9 中的值一致。
用于比较拟合模型与零模型的似然比检验的结果（如 coxph 的输出所示）是零模型和拟合模型的 \(−2\log\hat L\) 值之差异，即 215.940 − 202.938 = 13.00。这可以与两个自由度的卡方分布的百分位点进行比较，以检验模型中是否需要 \(Hb\) 和 \(Bun\)。如第 3 章所述，模型构建方法通常按顺序包含变量，并且通常不关注该检验的结果通常。相反，我们应为每个模型获得 \(−2\log\hat L\) 统计量，并将由于包含或省略项而引起的变化与卡方分布的百分位点进行比较。然而，当拟合了两个不同的模型时，R 输出中给出的似然比检验统计量之差将与 \(−2\log\hat L\) 值之差相同。
使用 coxph 拟合的模型的 AIC 或 BIC 统计量可以通过分别将 \(2q\) 或 \(q\log d\) 添加到 \(−2\log\hat L\) 的值来直接获得，其中 \(q\) 是模型中 \(\beta\) 参数的数量，\(d\) 是未删失观测的数量。或者，函数 extractAIC 可以应用于 coxph 对象。默认给出 AIC，但使用选项 k= 可用于设置 \(\log d\) 来给出 BIC。使用 k=2（默认值）将给出 AIC 的标准公式。
例如，确诊为多发性骨髓瘤患者的生存时间数据有 \(d=36\) 个事件。该值可利用 coxph 对象 coxregfit 的 coxregfit$nevent 获取，下面的代码给出包含 \(Hb\) 和 \(Bun\) 的模型的 AIC 和 BIC 统计量。
> # AIC
> extractAIC(coxregfit)
[1] 2.0000 206.9382
> # BIC
> extractAIC(coxregfit,k=log(coxregfit$nevent))
[1] 2.0000 210.1052在该输出中，第一个数字是模型中参数数量 \(q\) 的值，以及 AIC = 206.94，BIC = 210.10.
17.4.1 变量选择和 lasso
R 不包括用于变量选择的自动程序，因为它们有许多缺点，如第 3 章 3.6 节所述。然而，survival 包中的 anova 函数给出当按照模型公式中出现的顺序依次向模型中添加项时得到的 \(−2\log\hat L\) 的值。例如，anova(coxregfit) 的输出如下所示。
Analysis of Deviance Table
Cox model: response is Surv(time, status)
Terms added sequentially (first to last)
loglik Chisq Df Pr(>|Chi|)
NULL -107.97
hb -105.53 4.8717 1 0.027300 *
bun -101.47 8.1300 1 0.004354 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1这给出了最大化偏对数似然的值（从零模型开始，添加 \(Hb\)，然后添加 \(Bun\)），可以将其乘以 −2 以给出 \(−2\log\hat L\) 统计量。标题为 Pr(>|Chi|) 的列中的 \(P\) 值指的是序贯地添加每个项时 \(−2\log\hat L\) 的变化的显著性。将 \(Hb\) 添加到零模型中会显著降低 \(−2\log\hat L\)，将 \(Bun\) 添加到包含 \(Hb\) 的模型中也是如此。\(Bun\) 的 \(P\) 值 0.0043 对应于在考虑 \(Hb\) 的效应后，该变量对死亡风险的效应的显著性。
为了了解调整 \(Bun\) 后 \(Hb\) 的显著性，在 coxph 命令中以相反的顺序给出这两个变量，如下所示
第 3 章 3.7 节描述了使用 lasso 来进行 Cox 回归模型的变量选择，这可以使用 penalized 包中的 R 函数 penalized 来实现。现使用 lasso 程序来确定在拟合多发性骨髓瘤患者生存数据的 Cox 回归模型中需要七个潜在解释变量中的哪一个，如第 3 章示例 3.7 所示。安装并加载 penalized 包后（如 17.2.2 节所述），lasso 程序通过如下代码来实施
> penfit<-penalized(Surv(time,status),penalized= ~ hb+bun+age+sex+
ca+protein+pcells,lambda1=0,steps=20,standardize=TRUE,
data=myeloma)survival 对象 Surv(time,status) 再次用作响应，但模型公式在 penalized= 选项中给出。请注意，公式开头需要一个 ~ 符号。如果只惩罚某些参数，则可以在 unpenaltized= 选项中声明其余参数。选项 lambda1=0 指定使用 \(L_1\) 范数惩罚，它对应于 lasso，最小值为零。steps= 选项指定要使用 20 个 \(\lambda\) 值。选项 steps="Park" 使用 Park and Hastie (2007) 中描述的方法自动选择步长。若选项 standardize= 设置为 TRUE，则将解释变量设置在类似的测量尺度上，如 3.7.2 节所述。通过将 penalized 对象分配给 penfit，如下 R 命令
给出对应于每个 lasso 参数 \(\lambda\) 值的偏对数似然值，以及 \(\lambda\) 值的系数估计的轨迹，如图 3.2 所示。
图 3.2
17.4.2 预测能力和解释的变异的度量
第 3 章 3.12 节总结了预测能力和解释的变异的各种度量。Harrell’s \(c\) 统计量可使用 survival 包中的 concordance 函数计算，Gönen and Heller 的预测能力度量及其标准误可以使用 CPE 包中的函数 phcpe 计算，如下所示
> concordance(coxregfit) # Harrell's c 统计量
Call:
concordance.coxph(object = coxregfit)
n= 48
Concordance= 0.6716 se= 0.05308
concordant discordant tied.x tied.y tied.xy
585 286 0 20 0
> phcpe(coxregfit,CPE.SE=TRUE) # Gonen-Heller statistic
$CPE
[1] 0.6694815
$CPE.SE
[1] 0.04128793\(c\) 统计量也包含在 summary(coxregfit) 的输出中。
总结预测能力的 \(D\) 统计量以及表示为 \(R^2_P\) 和 \(R^2_D\) 的解释的变异的度量可以使用 survival 包中的 royston 函数获得。为多发性骨髓瘤数据拟合的包含 \(Hb\) 和 \(Bun\) 的模型的这些统计量如下所示。
> royston(coxregfit)
D se(D) R.D R.KO R.N C.GH
1.0076162 0.3505868 0.1950951 0.2728951 0.2399530 0.6694815 在该输出中，前两项是 \(D\) 统计量的值及其标准误，以第 3 章 3.12 节的表示法，R.D 为 \(R^2_D\)，R.KO 为 \(R^2_P\)。此外，R.N 为 Nagelkerke’s \(R^2\) 值 (Nagelkerk, 1991)，C.GH 是 Gönen and Heller 的一致性度量。这些度量值与第 3 章表 3.20 所示值一致。
17.4.3 时依 ROC 曲线
函数 survivalROC（在同名包中）可用于为生存数据拟合模型后估计时依 ROC 曲线。ROC 曲线由拟合模型的风险评分或线性预测器形成，默认情况下，风险评分的唯一值用作计算真阳性率和假阳性率的阈值。由于 survivalROC 仅需要风险评分来估计时依 ROC 曲线，因此该函数可与任何参数生存模型以及 Cox 回归模型一起使用。
为多发性骨髓瘤数据拟合 Cox 回归模型后，可以使用以下命令从 coxph 对象 coxregfit 中提取风险评分或线性预测器
lp1<-coxregfit$linear.predictors然而，lp1 中风险评分的结果值基于以均值为中心的变量 Hb 和 Bun 的值。因为 \(Hb\) 和 \(Bun\) 值的平均值分别为 10.2521 和 33.9167，这意味着线性预测器或风险评分是
\[\begin{aligned}-0.13364~(Hb-10.2521)+0.01857~(Bun-33.9167)\end{aligned}\]
基于未转换变量的得分为 \(-0.13364Hb+0.01857Bun\)，若要根据以零为中心的变量值计算，可以使用以下代码
这样做的效果是将每个变量的系数和样本均值的乘积之和加到由变量的中心值形成的线性预测器中。
安装并加载 survivalROC 后，可以通过以下方式使用 survivalROC 函数获得 12 个月时的时依 ROC 曲线的坐标。
mroc<-survivalROC(Stime=myeloma$time,status=myeloma$status,marker=lp0,
predict.time=12,method="KM")在该代码中，Stime= 和 status= 指定为骨髓瘤数据框中的生存时间和删失指示符，marker= 设置为 lp0 中的风险评分，时间点使用选项 predict.time=12 设置为 12 个月。使用 method="KM" 指定 Kaplan-Meier 估计方法。利用 survivalROC 对象 mroc，可使用 TP12<-mroc$TP 和 FP12<-mroc$FP 提取每个阈值的真阳性率和假阳性率 以及使用 mroc$AUC 获得 ROC 曲线下面积。之后，可以根据 FP12 绘制 TP12 的图形，得出 12 个月时的 ROC 曲线，如图 3.10 所示。
图 3.10
还可使用选项 method="NNE" 来实现最近邻估计方法。此外，使用选项 span=0.05，将式 (3.44) 中的平滑参数 \(\omega\) 设为 0.05.
17.5 Cox 回归模型中的模型检查
第 4 章描述的许多模型检查诊断都是基于各种类型的残差。为了说明如何使用 R 获得残差，将再次使用示例 1.3 中首次给出的多发性骨髓瘤患者生存时间数据。
17.5.1 残差分析
通过将 survival 包中的 residuals 函数应用于 survival 对象，可以轻松获得残差。在 17.4 节中，对多发性骨髓瘤数据拟合了 Cox 回归模型，并将结果存储在 survival 对象 coxregfit 中。然后可以得到鞅和偏差残差，并将其分配给变量 martres 和 devres，如下
> martres<-residuals(coxregfit,type=c("martingale"))
> devres<-residuals(coxregfit,type=c("deviance"))Cox-Snell 和修正的 Cox-Snell 残差可以使用第 4 章的式 (4.6) 和 (4.4) 根据鞅残差计算
> csres<-myeloma$status-residuals(coxregfit,type=c("martingale"))
> modcsres<-1-residuals(coxregfit,type=c("martingale"))得分、Schoenfeld 和缩放的 Schoenfeld 残差可分别使用 residuals 函数中的选项 type=c("score"), type=c("schoenfeld") 和 type=c("scaledsch") 获得。鞅和偏差残差对于每个个体都有一个单独的值，并且作为向量提供。然而，得分、Schoenfeld 和缩放的 Schoenfeld 残差是模型中每个个体和每个解释变量的值的矩阵。
对于多发性骨髓瘤数据，当 residuals 函数用于获得变量 \(Hb\) 和 \(Bun\) 的得分残差时，得分残差矩阵具有与数据框相同顺序的 48 行，每行对应每个患者，列对应于 \(Hb\) 和 \(Bun\) 的得分残差。这些列用数据框中的变量名（即 hb 和 bun）进行标记，但通常可以方便地更改这些名称，以避免与原始数据中的变量混淆。要获得得分残差并将列名更改为 score.hb 和 score.bun，请使用
将残差与原始数据一起置与数据框中也很有用，这样就可以很容易地对其进行操作。为此，使用 cbind 函数绑定包含残差的列。例如，若要将利用 residuals 函数计算的包含 Cox-Snell、鞅、偏差和得分残差的列绑定到名为 myeloma 的数据框，请使用
可以使用 print(myeloma.res) 检查名为 myloma.res 的对象中变量的值，其中一些变量的前 10 个值如下所示。
patient time status csres martres devres score.hb score.bun
1 1 13 1 0.23014540 0.76985460 1.18253068 3.2999546008 2.4140876
2 2 52 0 0.90494858 -0.90494858 -1.34532419 -1.9303609057 13.8875587
3 3 6 1 0.13748042 0.86251958 1.49783457 1.8290988299 -7.3828723
4 4 40 1 0.87666446 0.12333554 0.12880543 0.1792547576 1.9590824
5 5 10 1 0.20639730 0.79360270 1.25247724 2.9277335613 -2.6913044
6 6 7 0 0.15889219 -0.15889219 -0.56372367 -0.0772011450 6.0619565
7 7 66 1 1.34525433 -0.34525433 -0.31199757 -3.0274073993 7.0127113
8 8 10 0 0.27393472 -0.27393472 -0.74018204 -1.2483218431 -0.1676842
9 9 10 1 0.60622895 0.39377105 0.46200976 -0.8189308283 12.1545183
10 10 14 1 0.57177193 0.42822807 0.51144310 -0.0216079855 12.5389899对于一些诊断图，将生存时间的秩和风险评分与其他残差置于同一数据框中非常有用。通过将 rank 函数应用于生存时间，可以很容易地得出其秩，例如使用 ranktime<-rank(time)。线性预测器 \(\hat{\boldsymbol{\beta}}^{\prime}\boldsymbol{x}_i\) 可按照 17.4.3 节描述的方式从 coxph 对象获得。所得评分以变量均值为中心，但除了横轴上的这种偏移外，残差与风险分数的关系图将与第 4 章中所示的相同。然而，本书中使用的未中心化的风险评分可以使用如下命令获得
然后可将其与数据框 myeloma.res 绑定
然后可绘制残差图
它给出了如图 4.6 所示的偏差残差与风险评分的关系图。4.2.1 节中描述的 Cox-Snell 残差的累积风险图可以通过使用 survfit 函数并将残差指定为时间变量来获得。即
图 4.6
累计风险估计可从 survfitc 对象 cskm 中找到，即 cskm$cum，这给出了 2.3.3 节中定义的累积风险函数的 Nelson-Aalen 估计。累积风险函数也可根据 \(-\log\hat{S}(r_{C{i}})\) 估计，其中，\(\hat{S}(r_{C{i}})\) 是 Cox-Snell 残差的生存函数的 Kaplan-Meier 估计，这是使用 ch<- -log(cskm$surv) 计算的。然后，可绘制累积风险关于 Cox-Snell 残差（通过 cskm$time 获得）的图形，即 plot(ch ~ cskm$time)。所得图形与图 4.5 相似。
图 4.5
Schoenfeld 残差矩阵对于数据集中的每个事件时间有一行，每个变量有一列。行按与事件时间相关的变量排序。这些很难与其他残差组合在一起，但 Schoenfeld 残差和缩放的 Schoenfeld 残差的矩阵的行以事件时间作为标签。例如，多发性骨髓瘤数据的 Schoenfeld 残差可使用如下代码获得
矩阵 schoenres 的前 10 行为
> print(schoenres)
hb bun
1 -0.04266125 78.30817
1 -4.34266125 -66.69183
1 1.85733875 -30.69183
4 0.61182088 100.34294
4 -2.98817912 -23.65706
5 0.25208340 -28.41657
5 0.95208340 -15.41657
5 0.75208340 88.58343
5 -0.44791660 -18.41657
6 2.10211724 -12.21550同样，列可以重新命名，将行标签中的列（即事件时间）转换为变量也很有用。然后将矩阵转换为数据框以进行进一步操作，这通过如下代码完成
> colnames(schoenres)<-c("schoen.hb","schoen.bun")
> schoenres<-cbind(time=rownames(schoenres),
data.frame(schoenres,row.names=NULL))
> print(schoenres)并得到包含变量 time 的数据框 schoenres。
获取并操作缩放的 Schoenfeld 残差矩阵也是类似地。这得到了一个数据框，可以轻松地用于生成缩放的 Schoenfeld 残差与事件时间的关系图，以检查比例风险的假设，如第 4 章 4.4.2 节所述。
17.5.2 识别有影响的观测
delta-beta，即 \(\Delta_i\hat{\beta}_j\) 是模型拟合过程中省略第 \(i\) 个观测时第 \(j\) 个参数估计的近似变化。该量也可以缩放以给出标准化的 delta-beta。这两个量均在第 4 章 4.3.1 节进行了描述，用于探索特定观测如何影响参数估计及其显著性。这两个统计量都可从 survival 包中的 residuals 函数中获得。例如，对于为多发性骨髓瘤数据拟合的 Cox 回归模型中两个变量 \(Hb\) 和 \(Bun\)，其 delta-beta 和缩放的 delta-beta 可通过如下命令得到
这两个函数都会生成一个矩阵，其行对应于 48 名患者，列对应于两个解释变量 \(Hb\) 和 \(Bun\)。可以更改列名称，并将该矩阵与包含其他诊断的数据框 myeloma.res 合并，方法与 17.5.1 节中的得分残差相同。
在第 4 章 4.3.2 节介绍的两个统计量 \(LD_i\) 和 \(|\boldsymbol{l}_{\max}|_i\)，度量了个体观测对参数估计集的影响，可以很容易使用 R 的矩阵代数功能获得。具体而言，似然位移 \(LD_i\) 是式 (4.16) 中矩阵 \(\boldsymbol B\) 的第 \(i\) 个对角元，\(|\boldsymbol l_\max|_i\) 是与 \(\boldsymbol B\) 的最大特征值相对应的特征向量。对于多发性骨髓瘤数据，这些诊断可从拟合模型中的得分残差、得分和参数估计的方差-协方差阵中获得，该矩阵是从 coxph 对象 coxregfit 的 coxregfit$var 获得的。下面的 R 代码实现了这一点。
> B<-scoreres%*%coxregfit$var%*%t(scoreres)
> LD<-diag(B)
> e<-eigen(B)
# 按递减顺序给出特征值
> lmax<-abs(e$vectors[,1])在此代码中，使用 %*% 运算符执行矩阵乘法，函数 t 用于获得矩阵的转置。矩阵运算 diag 给出 \(\boldsymbol B\) 的对角元，函数 eigen 返回一个分量列表，包括按降序排列的特征值和相应的标准化特征向量。将该列表指定为 e，与最大特征值对应的特征向量是第一个特征向量，对应于特征向量形成的矩阵中的第一列，这可以使用 e$vectors[,1] 来访问。然后将这些值的绝对值指定为 lmax。可以使用 cbind 函数将变量 LD 和 lmax 添加到另一个数据框中，以进行进一步分析。
17.5.3 检验比例风险假设
第 4 章4.4.1 节描述的数据的对数累积风险图可以根据生存函数的 Kaplan-Meier 估计得出。在第 4 章的示例 4.9 中，与血清血红蛋白水平相关的变量 \(Hb\) 被重新定义为具有四个水平的因子，并检查了在这些水平上的比例风险假设。可通过如下命令在 myeloma 数据框中为因子 hbfactor 定义四个水平 1, 2, 3, 4
> myeloma$hbfactor<-as.factor(ifelse(myeloma$hb<=7,"1",
ifelse(myeloma$hb<=10,"2",
ifelse(myeloma$hb<=13,"3",
ifelse(myeloma$hb>13,"4",0)))))并使用如下命令生成 survfit 对象 km1
survfit 对象包括变量 km1$time 和生存函数 km1$surv 的 Kaplan-Meier 估计。要绘制对数累积风险，必须将 \(Hb\) 的因子水平附加给 km1 对象。由于该对象与原始数据框的长度不同，因此无法直接将 km1 对象与 hblevels 中的因子水平合并。然而，km1 包含了与血红蛋白水平相关的变量 km1$strata，该变量可用于生成 \(Hb\) 的因子水平，如下
该命令将值 1, 2, 3, 4 分配给 hblevels，并根据四个层中的观测数重复每个水平。对数累积风险图可使用对数累积风险的 Nelson-Aalen 估计或 Kaplan-Meier 估计来获得，通过如下的 plot 函数获得
> plot(log(km1$cum) ~ log(km1$time),pch=hblevels)
> plot(log(-log(km1$surv)) ~ log(km1$time),pch=hblevels)plot 中的选项 pch= 为每个 \(Hb\) 水平提供了不同的绘图符号。这两个图非常相似，但第二张图得到了图 4.13.
作者勘误
事实上第二张图与图 4.13 稍有不同。为复现图 4.13，生存函数估计和 km1 对象中的生存时间需要限制在事件时间内。一种方法是将生存函数估计值、生存时间、变量 hblevels 和每次事件的数量保存到数据框 lchplotdata 中。然后使用如下命令选择 km1$n.event>0 的行的子集
图 4.13
第 4 章 4.4.2 节描述的缩放的 Schoenfeld 残差与生存时间的关系图也可用于评估比例风险假设的有效性。该图可按照 17.5.1 节末尾描述的方式从缩放的 Schoenfeld 残差矩阵中获得。
第 4 章 4.4.3 节描述的检验统计量可使用 coxph 对象应用 survival 包中的函数 cox.zph 得到。例如，对于多发性骨髓瘤数据
给出了拟合模型中变量 \(Hb\) 和 \(Bun\) 的各个 Grambsch and Therneau 检验的卡方统计量及其 \(P\) 值以及总体 zph 检验。
17.6 参数生存模型
生存数据的参数模型，如第 5 章所述，可以使用 survival 包中的 survreg 函数进行拟合。该函数使用 5.11 节中描述的模型的对数线性表示来拟合加速失效时间模型，并且可用于对假定具有指数, Weibull, lognormal 或 loglogistic 等分布的生存时间进行建模。
在与第 \(i\) 个（共 \(n\) 个）生存时间相关的随机变量 \(T_i\) 的加速失效时间模型中，
\[\begin{align} \log T_i=\mu+\alpha_1x_i+\alpha_2x_{2i}+\cdots+\alpha_px_{pi}+\sigma\epsilon_i \tag{17.1} \end{align}\]
其中 \(\mu\) 和 \(\sigma\) 分别是截距和尺度参数，\(\alpha\) 是模型中 \(p\) 个解释变量的系数，\(\epsilon_i\) 是一个随机变量，其概率分布取决于 \(T_i\) 的分布。解释变量对生存的效应可总结为时间比 \(\exp(\alpha)\)，这是对应解释变量值每增加一个单位时，生存时间所乘的倍数。
survreg 函数的语法与 coxph 的语法非常相似，默认情况下拟合 Weibull 模型。例如，以下代码可用于为第 1 章示例 1.3 中首次给出的确诊多发性骨髓瘤后的生存时间数据拟合 Weibull 模型。
使用 summary(weibfit) 获得如下输出
Call:
survreg(formula = Surv(time, status) ~ hb + bun, data = myeloma)
Value Std. Error z p
(Intercept) 2.70425 0.57769 4.68 2.9e-06
hb 0.11456 0.05375 2.13 0.03306
bun -0.01580 0.00423 -3.74 0.00019
Log(scale) -0.10644 0.12746 -0.84 0.40365
Scale= 0.899
Weibull distribution
Loglik(model)= -153.5 Loglik(intercept only)= -159.8
Chisq= 12.57 on 2 degrees of freedom, p= 0.0019
Number of Newton-Raphson Iterations: 5
n= 48根据该输出，式 (17.1) 的加速失效时间模型中的参数估计为 \(\hat\mu = 2.704，\hat\sigma = 0.899\)，变量 \(Hb\) 和 \(Bun\) 的系数估计分别为 0.1146 和 -0.0158. 输出中的 Log(scale) 参数就是 \(\log \hat\sigma\)。这部分输出给出了每个参数估计的标准误、Wald 统计量（计算为 \(z=\hat{\alpha}/\mathrm{se}\left(\hat{\alpha}\right)\) ）以及相应的 \(P\) 值。\(Hb\) 和 \(Bun\) 均显著影响生存。\(Hb\) 值每增加一个单位，生存时间就会增加为原来的 \(\exp(0.1146) = 1.121\) 倍，即增加 12%，而 \(Bun\) 值每增加一个单位，生存时间就会减少 1.4%.
模型的对数似然 Loglik(model) 使用第 5 章的式 (5.46) 计算，并且 \(-2\log\hat{L}\) 统计量的值根据 -2*weibfit$loglik 计算，结果为 306.996. 该统计量用于比较替代模型。仅具有截距的模型的对数似然，标记为 Loglik(intercept only),，是零模型的对数似然，输出中的卡方统计量将拟合模型与该零模型进行比较，这种比较通常为意义不大。输出还给出了拟合模型所需的迭代次数以及数据框中的观测值数 n。
Weibull 加速失效时间模型也可表示为比例风险模型，其中第 \(i\) 个个体在时间 \(t\) 的死亡风险为
\[\begin{aligned}h_{\boldsymbol{i}}(t)=\exp(\boldsymbol{\beta'x}_{\boldsymbol{i}})h_{0}(t)\end{aligned}\]
假定所有解释变量为零的个体的生存时间具有尺度参数为 \(\lambda\) 和形状参数为 \(\gamma\) 的 Weibull 分布，基线风险函数由 \(h_0(t)=\lambda\gamma t^{\gamma-1}\) 给出。加速失效时间的参数与 Weibull 模型的比例风险版本之间的对应关系如第 5 章的式 (5.53) 所示。具体来说，模型比例风险版本中的对数风险比为\(\beta=-\alpha/\sigma\)。\(Hb\) 和 \(Bun\) 值单位变化的对数风险比估计分别为 \(−0.1146/0.899=−0.1275\) 和 \(0.0158/0.899=0.0176\)，从而风险比分别为 0.88 和 1.02.
可以使用 survreg 和选项 dist= 来拟合具有其他生存时间分布的加速失效时间模型，该选项具有 "exponential", "loglogistic" 或 "lognormal" 等值。在每种情况下，均会给出式 (17.1) 中对数线性模型中的截距、尺度和 \(\alpha\) 参数。这些模型和其他模型参数之间的对应关系在第 5 章中进行了描述。在每种情况下，输出都与 Weibull 模型获得的输出非常相似，并且使用 \(-2\log\hat{L}\) 统计量来比较替代模型。例如，使用以下代码为多发性骨髓瘤数据拟合 loglogistic 加速失效时间模型
summary(llfit) 得到
Call:
survreg(formula = Surv(time, status) ~ hb + bun, data = myeloma,
dist = "loglogistic")
Value Std. Error z p
(Intercept) 2.05628 0.68561 3.00 0.00271
hb 0.12619 0.06328 1.99 0.04612
bun -0.01392 0.00429 -3.25 0.00116
Log(scale) -0.45293 0.13562 -3.34 0.00084
Scale= 0.636
Log logistic distribution
Loglik(model)= -153.6 Loglik(intercept only)= -159.9
Chisq= 12.62 on 2 degrees of freedom, p= 0.0018
Number of Newton-Raphson Iterations: 4
n= 48 Weibull 模型和 loglogistic 模型的时间比和对数似然值非常相似。
17.7 灵活的参数模型
第 6 章中介绍的分段指数模型和基于样条函数的灵活参数模型可以很容易地使用 R 进行拟合。
17.7.1 分段指数模型
第 6 章 6.1 节描述的分段指数模型可使用 eha 包中的 phreg34 函数进行拟合。第 3 章示例 3.4 中关于 36 名肾癌患者生存时间的数据将用于说明该函数，第 6 章示例 6.1 中拟合了分段指数模型。首先输入关于生存时间、生存状态以及变量\(Age\) 和 \(Nephrectomy\) 的数据，并使用以下代码将变量 \(Age\) 定义为水平 <60, 60-70 和 >70 的因子。
> kidney<-read.table("Treatment of hypernephroma.dat",header=TRUE)
> kidney$agef<-factor(kidney$age)
> levels(kidney$agef)<-list("<60"=1,"60-70"=2,">70"=3)为拟合基线风险函数具有 22 个区间的模型，使其等价于 Cox 回归模型，使用
> pefit <- phreg(Surv(time,status) ~ agef+nephrectomy, dist="pch",
cuts=c(5,6,8,9,10,12,14,15,17,18,21,26,35,36,38,48,52,56, 68,72,84,108,115),data=kidney)函数 phreg 可用于拟合许多分布，默认为 Weibull. 使用选项 dist="pch" 拟合分段指数模型，选项 cuts= 指定基线风险函数的分割点。使用 print(pefit) 会生成以下输出
Call:
phreg(formula = Surv(time, status) ~ agef + nephrectomy, data = kidney,
dist = "pch", cuts = c(5, 6, 8, 9, 10, 12, 14, 15, 17, 18,
21, 26, 35, 36, 38, 48, 52, 56, 68, 72, 84, 108, 115))
Covariate W.mean Coef Exp(Coef) se(Coef) Wald p
agef
<60 0.643 0 1 (reference)
60-70 0.316 0.117 1.124 0.427 0.784
>70 0.041 1.508 4.518 0.608 0.013
nephrectomy 0.934 -1.536 0.215 0.530 0.004
Events 32
Total time at risk 1340
Max. log. likelihood -124.35
LR test statistic 13.95
Degrees of freedom 3
Overall p-value 0.00296981输出中标题为 W.mean 的列中的条目是变量的加权平均值或因子水平的相对频率。权重是每个因子水平或变量值的患者所占的总观测生存时间的比例。在这个数据集中，36 名患者的总生存时间是 1340 个月，这在输出中标记为 Total time at risk。对于因子 \(Age\) 的第 1 水平，有 19 名患者，该组的生存时间总和为 862. 第一年龄组患者所经历的总体生存时间的比例则为 \(862/1340 = 0.643\)，如输出所示。类似地，对于变量 \(Nephrectomy\)，\(Nephrectomy=0\) 的患者的总生存时间是 88 个月，而 \(Nephrectomy=1\) 的患者的总生存时间是 1252 个月。该变量的加权平均值则为 \((88 × 0 + 1252 × 1)/1340 = 0.934\)，如输出所示。
下一个输出块给出模型中每一项的系数和风险比，以及它们的标准误和从统计量 \(z=\hat{\beta}/\text{se}(\hat{\beta})\) 导出的 Wald \(P\) 值。因子的风险比是相对于第一水平的，并且与示例 6.1 中使用的一致。统计量 Max. log. likelihood 是拟合模型的最大对数似然，使用 -2*pefit$loglik 将其乘以 -2 即可获得 \(-2\log\hat{L}\) 统计量。LR test statistic 和相应的 Overall p-value 是拟合模型和零模型之间的比较。
phreg 对象 pefit 包含许多其他变量，包括使用 pefit$hazards 得出的 22 个区间中的基线风险函数值。
17.7.2 风险函数的模型
基线风险函数对数的 B 样条模型在贝叶斯生存分析中得到了更广泛的应用，稍后将在 17.15 节中进行说明。
17.7.3 对数累积风险函数的模型
如第 6 章 6.4 节所述，使用限制性立方样条对对数风险风险进行建模的 Royston and Parmar 模型可以使用软件包 flexsurv 中的函数 flexsurfsurf 进行拟合。一旦安装并加载了该包，就可以轻松拟合具有不同结数的模型。再次考虑第 1 章示例 1.3 中介绍的多发性骨髓瘤确诊后的生存数据，以及 17.4 - 17.6 节使用的数据。如第 6 章 6.4 节所述，没有内部结的 Royston and Parmar 模型对应于 Weibull 模型。选项 k= 用于指定结数，因此包含变量 \(Hb\) 和 \(Bun\) 的 Weibull 模型也可以使用如下命令拟合
这会生成如下输出
Call:
flexsurvspline(formula = Surv(time, status) ~ hb + bun, data = myeloma,
k = 0)
Estimates:
data mean est L95% U95% se exp(est) L95% U95%
gamma0 NA -3.00773 -4.57689 -1.43858 0.80060 NA NA NA
gamma1 NA 1.11229 0.83442 1.39015 0.14177 NA NA NA
hb 10.25208 -0.12744 -0.24453 -0.01036 0.05974 0.88034 0.78307 0.98970
bun 33.91667 0.01758 0.00730 0.02786 0.00525 1.01773 1.00733 1.02825
N = 48, Events: 36, Censored: 12
Total time at risk: 1122
Log-likelihood = -153.4981, df = 4
AIC = 314.9963输出的第一部分给出变量 \(Hb\) 和 \(Bun\) 的样本均值，后面是模型的参数估计。在第 6 章 6.4 节的表示法中，\(\hat{\gamma}_0=-3.0077,\hat{\gamma}_1=1.1123\)，模型线性部分中两个变量的 \(\beta\) 系数估计分别为 -0.1274 和 0.0176. 由于这里使用 Weibull 模型的比例风险版本，这些系数估计是对数风险比，它们的值与 17.4 节中拟合的 Cox 回归模型得到的值非常相似。\(\gamma\) 与 Weibull 基线风险中的参数 \(h_0(t)=\lambda\gamma t^{\boldsymbol{\gamma}-1}\) 相关，并且根据第 6 章 6.4 节，\(\gamma_0=\log\lambda\) 且 \(\gamma_1 = \gamma\)。输出展示了系数的标准误及其 95% 置信限。标题为 exp(est) 的列给出了系数的指数，即风险比及其 95% 置信限。在输出观测数、事件数和删失数的总结后，将生存时间的总和展示为 Total time at risk。然后是拟合模型的最大对数似然和自由度数，即未知模型参数的数量，以及模型比较中使用的 AIC 统计量。
使用一系列具有不同选项 k= 值的 flexsurvspline 函数可以轻松拟合结数不断增加的模型，并使用 AIC 进行比较。第 6 章 6.5 节中描述的具有限制性立方样条函数的比例几率模型也可以使用带有选项 scale="odds" 的 flexsurvspline 函数进行拟合。输出非常相似，但拟合模型的系数估计现在是对数几率比。
17.8 参数模型中的模型检查
第 7 章 7.1 节描述的参数模型的残差很容易获得。survreg 对象包含了加速失效时间模型的线性部分估计，由下式给出
\[\begin{aligned}\hat{\eta}_i=\hat{\mu}+\hat{\alpha}_1x_i+\hat{\alpha}_2x_{2i}+\cdots+\hat{\alpha}_px_{pi}\end{aligned}\]
从中可以得到标准化残差。为多发性骨髓瘤数据拟合后 Weibull 加速失效时间模型后，survreg 对象 weibfit 在 weibfit$linear.prpredictors 中包含了线性预测器 \(\hat{\eta}_i\)，以及在 weibfit$scale 中包含了尺度参数 \(\hat\sigma\)。然后使用如下命令得到标准化残差 \(r_{Si}\)
Cox-Snell 残差是使用第 7 章的式 (7.3) 和 (7.2) 获得的，其中对于 Weibull 分布的特定情况，\(S_{\epsilon_i}(\epsilon)=\exp(-e^{\epsilon})\)，如式 (7.7) 所示，因此 \(r_{Ci} = \exp(r_{Si})\)。这些残差可简单地使用 csres<-exp(stres) 计算。然后使用如下命令计算鞅和偏差残差
> martres<-myeloma$status-csres
> devres=sign(martres)*sqrt(-2*(martres+
myeloma$status*(log(myeloma$status-martres))))偏差残差也可以直接使用
来获得。然而，如此得到的偏差残差将与根据鞅残差计算出的偏差残差的符号相反。第 7 章 7.1.4 节解释了这种符号差异的原因。
与从拟合的 Cox 回归模型中获得的残差一样，残差列可以与数据框绑定，以便按照 @ref*sec17-5-1) 节所述的方式在绘图和进一步分析中使用。
17.8.1 有影响的值
可以使用第 7 章 7.4 节描述的 delta-beta 统计量来研究拟合参数生存模型中观测对参数估计 \(\hat\alpha\) 的影响。可以对 survreg 对象应用 residuals ，带上选项 type=c("dfbeta") 或 type=c("dfbetas") 以提取未缩放和缩放的 delta-beta. 第 7 章 7.4.2 节描述的似然位移统计量可使用选项 type=c("ldcase") 得到。
17.8.2 比较观测的和拟合的生存函数
将基于模型的生存函数估计与观测的 Kaplan-Meier 估计进行比较，提供了一种评估模型拟合的非正式方法，如第 7 章 7.3 节所述。需要进行一些操作，以获得具有不同风险评分水平的个体的生存函数，使用 26 名女性接受两种卵巢癌化疗后的生存时间数据来说明该过程，如第 5 章示例 5.11 所示。
数据首先被输入并且被指定为数据框 ovarian。然后为数据拟合包含变量 \(Age\) 和 \(Treat\) 的 Weibull 模型并保存为 survreg 对象 ofit。这是通过如下命令实现的。
> ovarian<-read.table("Chemotherapy in ovarian cancer patients.dat",
header=TRUE)
> ofit<-survreg(Surv(time,status) ~ age+treat,data=ovarian)然后从 Weibull 模型拟合的加速失效时间版本中提取每个患者的风险评分，并将其保存在变量 eta 中。此外，获得每个评分的秩，此外，还获得了每个评分的秩，以便根据评分的秩是否在 1-8, 9-17 或 18-26 的范围内，将患者分为三个风险组。这个变量被添加到数据框中，并得到了第 7 章表 7.2 中所示的分组。完成此操作的代码如下所示。
> ovarian$lp=ofit$linear.predictors > ovarian$ranklp<-rank(ovarian$lp)
> ovarian$group<-as.factor(ifelse(ovarian$ranklp<=8,"1",
ifelse(ovarian$ranklp<=17,"2",
ifelse(ovarian$ranklp>17,"3",0))))表 7.2
接下来，对于 26 名女性中的每一位，计算 Weibull 模型拟合的生存函数，其生存时间 t 的值从 0 到 1200，步长为 50. 这是通过如下代码完成的。
> ov1<-ovarian[rep(seq_len(26),each=25),]
> tt<-data.frame(tt=seq(0,1200,by=50))
> tt<-tt[rep(seq_len(nrow(tt)),26),]
> sdata<-cbind.data.frame(ov1,tt)结果是一个新的数据框 sdata，具有 \(25×26=650\) 行。最后，使用下式计算生存函数估计
\[\hat{S}(t_i)=\exp\left\{-\exp\left(\frac{\log t_i-\hat{\eta}_i}{\hat{\sigma}}\right)\right\}\]
其中 \(\hat{\eta}_i=\hat{\mu}+\hat{\alpha}_1Age_i+\hat{\alpha}_2\) 为线性预测器。\(Age_i\) 和 \(Treat_i\) 是第 \(i\) 个女性的 \(Age\) 和 \(Treat\) 的值。此 R 代码基于函数 aggregate 来形成 tt 和 group 中值的组合的组均值，如下所示。
> sdata$survest<-exp(-exp((log(sdata$tt)-sdata$lp)/ofit$scale))
> avsurv<-aggregate(survest ~ group+tt,FUN=mean,data=sdata)然后，数据框 avsurv 包含按风险组 group 和时间 tt 划分的 survest 中生存函数估计的均值。最后，可以使用以下命令轻松获得每个风险组的生存函数的 Kaplan-Meir 估计
然后，可以将平均生存函数叠加在每组的 Kaplan-Meier 估计图上，该图是从 survival 对象 KMest 使用如下命令获得的
> plot(KMest,xlab="Survival time",ylab="Survivor function")
> lines(survest ~ tt,pch=c("o","+","*")[as.numeric(group)],type="p", data=avsurv)给出的图形类似于图 7.5，但在 plot 函数中使用额外的选项来控制布局会得到一个更令人满意的图。
17.9 时依变量
通过以计数过程格式表示数据，可以将时依变量纳入 Cox 回归模型。第 8 章 8.3 节概述了该数据格式，第 13 章13.1 节对其进行了更详细的描述。为了实现这一点，通常需要进行一定量的数据操作，并且需要小心避免代码错误。
survival 包中的 R 函数 tmerge 可用于表示计数过程格式的数据集。该函数在关于时依变量的 R vignette中进行了描述 (Therneau, Crowson and Atkinson, 2022). 总之，该函数用于指定添加时依变量的基础数据集。然后，在函数的后续调用中定义时依变量，并与基础数据集合并。此过程最好用一个例子来说明，将使用骨髓移植治疗白血病的数据。该数据集在第 8 章的示例 8.2 中进行了描述，数据列于表 8.2 中。当患者在血小板水平恢复正常之前死亡时，这在变量 \(Ptime\) 中为缺失值，用 . 表示。因此，选项 na.strings=c(".") 附加于 read 语句中，从而为缺失值创建 NA。这里不需要变量 \(Precovery\)，因此从数据框中删除。使用以下代码还定义了与三个疾病组相关的因子，并删除了疾病组的原始变量。
> leukaemia<-read.table("Bone marrow transplantation in the treatment
of leukaemia.dat",header=TRUE,na.strings=c("."))
> leukaemia$fgroup<-as.factor(leukaemia$group)
> leukaemia<-subset(leukaemia,select=-c(group,precovery))表 8.2
然后使用 tmerge 函数以计数过程格式表示数据
这导致三个新变量被添加到基础数据框 leukaemia 中。变量 tstart 和 tstop 分别给出每个区间的开始和停止时间（这里都为零），而变量 \(Time\) 中的值代表实际的生存时间。选项 id= 指定用于合并使用 tmerge 创建的数据框的标识符，选项 death=event(time,status) 指定新变量 death 用于表示每个时间区间（tstart 到 tstop）结束时个体的状态。数据框 tdvdata 有 23 行，前 5 行如下所示。
patient time status page dage ptime fgroup tstart tstop death
1 1 1199 0 24 40 29 1 0 1199 0
2 2 1111 0 19 28 22 1 0 1111 0
3 3 530 0 17 28 34 1 0 530 0
4 4 1279 1 17 20 22 1 0 1279 1
5 5 110 1 28 25 49 1 0 110 1接下来，第二次调用 tmerge 函数用于指定 \(Ptime\) 是表示血小板恢复时间的变量。然后，从这个变量中导出的时变血小板水平指示符与先前创建的数据框 tdvdata 进行合并。选项 tdc= 指定 Plate 为新的时依变量，并将其添加到第一次调用 tmerge 得到的数据中。以下是执行此操作的代码。
现在得到的数据框 tdvdata 中，除了两名血小板没有恢复到正常水平的患者外，每位患者都有两行数据。数据现在处于计数过程格式，前六行如下所示。
patient time status page dage ptime fgroup tstart tstop death Plate
1 1 1199 0 24 40 29 1 0 29 0 0
2 1 1199 0 24 40 29 1 29 1199 0 1
3 2 1111 0 19 28 22 1 0 22 0 0
4 2 1111 0 19 28 22 1 22 1111 0 1
5 3 530 0 17 28 34 1 0 34 0 0
6 3 530 0 17 28 34 1 34 530 0 1在这种格式中，第一个患者的血小板水平指示符为零，当他的血小板水平恢复正常（第 29 天），血小板水平指示符变为一，其余患者以此类推。
现在，可以拟合包含新时依变量的 Cox 回归模型。为此，coxph 函数中的 Surv 对象用于指定开始和停止时间以及每个时间区间结束时患者的状态。使用以下函数调用拟合包含与疾病组相关的因子和时依变量的模型。
下面展示了输出结果的一部分，与前面 17.4 节描述的非常相似。
Call:
coxph(formula = Surv(tstart, tstop, death) ~ Plate + fgroup,
data = tdvdata, ties = "breslow")
coef exp(coef) se(coef) z p
Plate -2.31956 0.09832 1.22938 -1.887 0.0592
fgroup2 -2.07610 0.12542 1.15105 -1.804 0.0713
fgroup3 0.38943 1.47614 0.61825 0.630 0.5288根据该输出，时依变量的系数为 −2.3196. 因子 fgroup 的三个水平分别对应于 ALL、低风险 AML 和高风险 AML，并且由第一个水平 fgroup1 的效应设为零。ALL 患者在给定时间相对于低风险 AML 患者的风险比为 \(\exp(0+2.0761)=7.97\)，并且高风险 AML 患者相对于低风险患者在给定时间的风险比是 \(\exp(0.3894+2.0761)=11.77\)，正如示例 8.2 中发现的。
为了进一步说明函数 tmerge 的使用，给出了 R 代码用于分析示例 8.4 来自肝硬化治疗假设研究的数据。第一步是将表 8.4 中的基本数据与表 8.5 中的数据结合起来，给出时变的 log(bilirubin) 值。首先使用以下命令输入两个数据集
> liverbase<-read.table("Data from a cirrhosis study (baseline).dat",
header=TRUE)
> lbrdata0<-read.table("Data from a cirrhosis study (lbr data).dat",
header=TRUE)表 8.4
表 8.5
要创建的时变变量需要包括时间起点的 log(bilirubin) 值，即表 8.4 中的基线值。这可以通过构建一个包含患者编号、时间起点和基线 log(bilirubin) 值的数据框，然后将其与表 8.5 中的数据结合，并确保合并后的数据集按患者和随访时间顺序排列来实现。以下是实现这一目标的代码。
> baselbr<-liverbase[c("patient","time","lbr")]
> baselbr$time<-c(0)
> lbrdata1<-rbind(baselbr,lbrdata0)
> lbrdata<-lbrdata1[order(lbrdata1$patient,lbrdata1$time),]数据框 lbrdata 现在与基线数据合并，以生成计数过程格式的数据集。这可以使用下面两个 tmerge 函数调用来完成。
> tdvlbr<-tmerge(data1=liverbase,data2=liverbase,id=patient,
death=event(time,status))
> tdvlbr<-tmerge(data1=tdvlbr,data2=lbrdata,id=patient,
lbrt=tdc(time,lbr))患者 5 和 8 在其死亡的同一天记录了 log(bilirubin) 值。当使用时依变量时，在给定时间记录的值在该时间之后立即应用，因此，在第 341 天为患者 5 记录的 log(bilirubin) 的新值和在第 182 天为患者 8 记录的 log(bilirubin) 的新值将不会出现在分析中。合并后的数据框有 52 个观测，log(bilirubin) 水平的基线值在变量 lbr 中，时依值在 lbrt 中。使用如下命令
coxlbr<-coxph(formula = Surv(tstart, tstop, death) ~ lbrt + treat, data = tdvlbr,
ties = "breslow")拟合包括治疗和时依 log(bilirubin) 值的 Cox 回归模型后，使用 summary(coxlbr) 获得以下输出。
Call:
coxph(formula = Surv(tstart, tstop, death) ~ lbrt + treat, data = tdvlbr,
ties = "breslow")
coef exp(coef) se(coef) z p
lbrt 3.6048 36.7731 2.2449 1.606 0.108
treat -1.4791 0.2278 1.3403 -1.104 0.270
Likelihood ratio test=14.45 on 2 df, p=0.0007299
n= 52, number of events= 8使用 > -2*logLik(coxlbr) 获得 \(−2\log\hat L\) 统计量，为
这与第 8 章示例 8.4 中引用的值一致。
不包含时依变量的模型可以拟合到计数过程格式下的扩展数据集 lbrdata，或者拟合到数据框 liverbase 中的 12 个观测。两者都会得到相同的结果。
17.9.1 时变系数
具有时变系数的模型已在第 4 章的 4.4.4 节描述，这与检验比例风险假设有关，并在 11.2 节中描述了通过阶跃函数对治疗效应进行建模。这样的模型也可使用函数 coxph 进行拟合。
在第 4 章 4.4.4 节中，Cox 回归模型的系数是时间的函数，用于评估比例风险。coxph 函数的时间变换功能 (time
transform facility) 可用于在模型中包含项 \(\beta(t)=\beta t\)。为此，在 coxph 中定义一个函数 tt以形成变量与某个时间函数的乘积，命令为 tt=function(x,t,...) x*g(t)，其中 g(t) 是 t 的指定函数，如 t, log(t) 甚至是样条函数。
为了说明这一点，我们回到第 4 章示例 4.1 中给出的透析患者感染发生率数据。在示例 4.12 中，我们通过为变量 \(Age\) 和 \(Sex\) 添加时变系数，来扩展包含这两个变量的 Cox 回归模型，以检查年龄和性别效应是否随时间保持不变。将数据输入到名为 dialysis 的 R 数据框后（dialysis <- read.table("Infection in patients on dialysis.dat", header=TRUE)），可以使用以下代码来拟合包括 \(Age\) 和 \(Sex\) 以及 \(Sex\) 的时变系数的 Cox 回归模型。
> tdcfit<-coxph(Surv(time,status) ~ age+sex+tt(sex),ties="breslow",
tt=function(x,t,...)x*t,data=dialysis)下面是使用 summary(tdcfit) 的部分输出。
Call:
coxph(formula = Surv(time, status) ~ age + sex + tt(sex), data = dialysis,
ties = "breslow", tt = function(x, t, ...) x * t)
n= 13, number of events= 12
coef exp(coef) se(coef) z Pr(>|z|)
age 0.03318 1.03374 0.02686 1.235 0.217
sex -1.28076 0.27782 2.53323 -0.506 0.613
tt(sex) -0.07422 0.92847 0.12190 -0.609 0.543对于该模型，使用 -2*logLik(tdcfit) 得出 \(−2\log\hat L\) 的值为 34.104，并且当省略项 tt(sex) 时，\(−2\log\hat L\) 增加 0.364 至 34.468，如示例 4.12 所示。拟合的模型具有时变的性别系数，其估计为 \(\hat{\beta}(t)=-1.281-0.074t\)。这表明性别效应随着时间的推移而增加，但并不显著。
在连续时间区间内具有不同系数的项会产生分段 Cox 模型，用于对治疗之间的非比例性进行建模，如第 11 章 11.2.1 节所述。可以使用 survival 包中的 survSplit 函数将时变系数的阶跃函数纳入 Cox 回归模型。此函数用于将时间尺度划分为所需数量的区间，并使用所得数据框来拟合分层模型。
使用第 11 章示例 11.2 中给出的前列腺癌患者生存时间数据来说明该过程。在该示例中，假定治疗效应在 1–360, 361–720, 721–1080 天和 1080 天以上的时间区间内是恒定的，但允许它们之间有所不同。数据首先通过如下命令输入以形成 R 数据框 gcancer
在所得数据框中，治疗效应由变量 treat 表示，单独化疗为零，联合治疗为一。然后使用 surveSplit 函数对时间段进行分层
这段代码得到 survSplit 对象 gcancer1，时间段缩短为 360、720 和 1080 天。选项 episode= 指定用于标记结果数据框中的时间段的变量名称。原始数据框第 50–58 行中患者 32, 33 和 34 的数据如下。
patient treat tstart time status tperiod
32 0 0 360 0 1
32 0 360 720 0 2
32 0 720 778 1 3
33 0 0 360 0 1
33 0 360 720 0 2
33 0 720 786 1 3
34 0 0 360 0 1
34 0 360 720 0 2
34 0 720 797 1 3数据采用计数过程格式，并使用新变量 tstart 以及原始变量 time 和 status 定义时间段和状态。时变的治疗效应可以通过按时间段分层并包括治疗与时间段的交互作用来建模。治疗与分层变量之间的交互作用在模型公式中表示为 treat:strata(tperiod)，并且使用以下命令拟合包含 treat 和交互项的 Cox 回归模型
> gcfit<-coxph(Surv(tstart,time,status) ~ treat+treat:strata(tperiod),
ties="breslow",data=gcancer1)下面给出了部分输出。
coef exp(coef) se(coef) z Pr(>|z|)
treat 0.877 2.405 0.335 2.62 0.0088
treat:strata(tperiod)tperiod=2 -1.129 0.323 0.535 -2.11 0.0348
treat:strata(tperiod)tperiod=3 -1.978 0.138 0.870 -2.27 0.0229
treat:strata(tperiod)tperiod=4 -2.060 0.128 0.786 -2.62 0.0088根据该输出，在四个时间区间内，联合治疗相对于单独化疗的对数风险比分别为 \(0.877,0.877 − 1.129 = −0.25,0.877 − 1.978 = −1.10\) 和 \(0.877 − 2.060 = −1.18\)。此外，anova(gcfit) 还可以检验治疗效应在四个时期内是否恒定。该模型过度参数化，这就是为什么没有给出 treat:strata(tperiod)tperiod=1 的项。使用如下命令可获得更易解释的结果。
交互项 treat:strata(tperiod) 现在会在每个时间段产生单独的治疗效应，部分输出如下所示。
coef exp(coef) se(coef) z Pr(>|z|)
treat:strata(tperiod)tperiod=1 0.877 2.405 0.335 2.62 0.0088
treat:strata(tperiod)tperiod=2 -0.251 0.778 0.417 -0.60 0.5464
treat:strata(tperiod)tperiod=3 -1.101 0.333 0.802 -1.37 0.1701
treat:strata(tperiod)tperiod=4 -1.182 0.307 0.711 -1.66 0.0964风险比和区间估计可直接从此输出中读取，并且与表 11.3 中列出的一致。
表 11.3
17.9.2 纵向数据和生存数据的联合建模
可以使用 JM 包为纵向数据和生存时间拟合联合模型，但也可以使用其他包，比如 joineR 包。首先使用 nlme 包中的 lme 函数通过线性混合模型对纵向变量进行建模。该函数要求指定模型中的固定效应和随机效应。然后使用 coxph 对生存时间进行建模，以生成生存对象。最后，使用 jointModel 函数将这些单独的纵向模型和生存模型组合成一个联合模型。
使用第 8 章示例 8.5 给出的主动脉瓣置换术后患者生存数据来说明该过程。该数据集包括患者生存时间和与左心室质量指数 LVMI 相关的变量的纵向值，以及在研究开始时测量的解释变量。首先使用如下命令将数据输入到数据框 valve。
纵向变量将使用与时间的直线关系进行建模，每个患者具有随机截距和斜率。代表患者性别和术前糖尿病状态的变量也将包括在内，因为这些变量在第 8 章的示例 8.8 中被证明是显著的，拟合该模型的 lme 函数如下。
在该函数中，random= 选项指定随机效应的模型为 ~ time，分组变量 id 在 | 符号后指定。
为了对生存时间进行建模，需要一个数据框，其中每行代表一个患者。为了从包含每个患者多行数据的数据框 valve 中构建该数据框，可以使用以下代码。
! 符号用于选择数据框 valve 中变量 id 值不重复的行。使用如下命令删除变量 lvmi 和 time，因为它们不再需要。
这会产生一个数据框 lsurv，其前 5 行如下。
id futime status age sex redo emerg dm type
1 4.956 0 75 0 0 0 0 2
2 9.663 0 46 0 1 0 0 1
3 7.915 0 63 0 0 0 0 1
4 4.038 0 61 0 1 0 0 1
5 8.816 0 54 0 1 0 0 1在示例 8.8 中，患者年龄是显著影响生存的唯一变量。然后以通常的方式拟合仅包含该变量的 Cox 回归模型，但包括选项 x=TRUE，以在模型中保留解释变量的值矩阵以供以后使用。
然后，使用函数 jointModel 将 lme 对象 lmefit 与 coxph 对象 coxfit 相结合，并使用 timeVar= 选项指定与纵向测量时间相关的变量。在这个阶段，有许多模型可以用于联合模型中的生存时间，包括 Cox 回归模型、分段指数模型、全参数模型和基于样条的灵活模型。此处将使用基线风险函数对数的样条模型，如第 8 章示例 8.9 所示。生成以下代码35。
method= 选项指定使用默认为五个内部结的比例风险样条曲线模型，并使用 Gauss-Hermite 求积进行所需的数值积分。可以使用 lng.in.kn= 选项设置不同的内部结数。可以类似的方式指定替代的比例风险模型，例如 method="Cox-PH-GH", method="piecewise-PH-GH" 或 method="weibull-PH-GH"。
summary(jmod)的部分输出如下。
Call:
jointModel(lmeObject = lmefit, survObject = coxfit, timeVar = "time",
method = "spline-PH-GH")
Data Descriptives:
Longitudinal Process Event Process
Number of Observations: 988 Number of Events: 54 (21.1%)
Number of Groups: 256
Joint Model Summary:
Longitudinal Process: Linear mixed-effects model
Event Process: Relative risk model with spline-approximated
baseline risk function
Parameterization: Time-dependent
log.Lik AIC BIC
-5437.314 10912.63 10979.99
Variance Components:
StdDev Corr
(Intercept) 59.6043 (Intr)
time 7.4949 -0.5819
Residual 36.1661
Coefficients:
Longitudinal Process
Value Std.Err z-value p-value
(Intercept) 175.0903 5.1812 33.7932 <0.0001
sex -29.6409 10.0306 -2.9551 0.0031
dm -20.9350 10.2557 -2.0413 0.0412
time -1.4126 0.9104 -1.5516 0.1208
Event Process
Value Std.Err z-value p-value
age 0.0859 0.0160 5.3598 <0.0001
Assoct 0.0052 0.0030 1.7325 0.0832在给出最大化对数似然值 \(\log\hat{L}\)，AIC 和 BIC 统计量后，随机截距标准差估计 \(\hat u_{0i}\) 和随机斜率标准差估计 \(\hat u_{1i}\) 与残差项的标准差估计 \(\hat \sigma_\epsilon\) 在题为 StdDev 的列中给出。下一列给出了截距和斜率之间的相关系数为 −0.582.
纵向过程的系数是对 \(m_i(t)\) 的拟合模型中的 \(\beta\) 系数的估计，\(m_i\) 是第 \(i\) 名患者的基本响应曲线。这里，拟合模型在时间 \(t\) 上是线性的，包括第 \(i\) 名患者的性别 \(s_i\) 和糖尿病状况 \(d_i\)，并在式 (8.12) 中给出，其中根据第 8 章示例 8.9：
\[\begin{aligned}\hat{m}_i(t)=175.09-1.41t-29.64s_i-20.93d_i\end{aligned}\]
系数估计的标准误、\(z\) 值和 \(P\) 值得到了区间估计和假设检验。
联合模型的风险函数估计是根据标记为 Event Process 的部分获得的，其中联合模型中的系数估计 \(m_i(t)\) 标记为 Assect。与第 8 章的示例 8.9 一样，这是
\[\begin{aligned}\hat h_i(t)=\exp\{0.0857a_i+0.0051\hat m_i(t)\}\end{aligned}\]
该输出的 Event Process 部分还包含对数基线风险的样条模型指定所需的 9 个立方 B 样条基函数的系数，但这些已被省略。
接下来，可以对 jointModel 对象应用许多操作，包括 predict 以获得纵向变量的拟合值，plot 以获得模型的各种图形总结，以及 residuals 以导出拟合的联合模型的各种残差。
17.10 区间删失数据
如第 9 章 9.1 节所述，精确事件时间、左删失事件时间和右删失事件时间都可以以区间删失数据的形式表示。右删失事件时间具有无限上限的区间，左删失时间具有零下限，而精确事件时间的区间上下限相等。上限为 \(L\)、下限为 \(R\) 的区间删失数据可以表示为 survival 对象，以便在模型公式 Surv(L,R,type="interval2") 中使用。无限上限可以指定为 Inf 或 NA。
本节将使用第 9 章示例 9.1 给出的溃疡复发时间数据进行说明。第一步是输入表 9.1 中的数据并执行必要的操作将其表示为区间删失格式，表 9.2 总结了这一点。变量 \(Treat\) 和 \(Dur\) 的值也被转换，使得 \(Treat\) 对于治疗 \(B\) 为 1，否则为 0，并且如果疾病持续时间至少为五年，则 \(Dur\) 为 1，否则为 0. “基线”个体是指接受治疗 A 且疾病持续时间少于五年的个体。这可以通过以下方式来完成。
> ulcer <- read.table("Recurrence of an ulcer.dat",header=TRUE)
> ulcer$left<-ifelse(ulcer$result==1,ulcer$time, ifelse((ulcer$result==2 & ulcer$time<=6),0,
ifelse((ulcer$result==2 & ulcer$time>6),6,0)))
> ulcer$right<-ifelse(ulcer$result==1,NA,ulcer$time)
> ulcer$duration<-ulcer$duration-1
> ulcer$treatment<-ulcer$treatment-1表 9.1
表 9.2
数据框 ulcer 的前十行如下。
patient age duration treatment time result left right
1 1 48 1 1 7 2 6 7
2 2 73 0 1 12 1 12 NA
3 3 54 0 1 12 1 12 NA
4 4 58 1 1 12 1 12 NA
5 5 56 0 0 12 1 12 NA
6 6 49 1 0 12 1 12 NA
7 7 71 0 1 12 1 12 NA
8 8 41 0 0 12 1 12 NA
9 9 23 0 1 12 1 12 NA
10 10 37 0 1 5 2 0 517.10.1 生存函数的非参数最大似然估计
第 9 章 9.2 节描述并说明了区间删失数据生存函数的非参数最大似然估计。可以使用 icfit 函数获得该估计，该函数是 interval 包的一部分，代码如下。
译者注：可能需要如下通过命令安装 interval 所依赖的包 Icens：
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Icens")使用 summary(icsurv) 会得到以下输出。
treatment=0:
Interval Probability
1 (0,2] 0.1579
2 (6,12] 0.0648
3 (12,Inf) 0.7773
treatment=1:
Interval Probability
1 (0,5] 0.0833
2 (6,7] 0.2083
3 (12,Inf) 0.7083如果需要，可以对 Surv 对象 Surv(left,right,type="interval2") ~ 1 应用 icfit 来拟合总体生存函数。
summary 函数的输出给出了每个治疗组的 Turnbull 区间中的生存概率。据此，治疗 A（\(Treat=0\)）患者的生存函数估计在时间起点为 1，对于 2-6 个月之间的复发时间，为 \(1 − 0.1579 = 0.8421\)，在 12 个月时为 \(1 − 0.1579 − 0.0648 = 0.7773\)，而对于超过 12 个月是未定义的，如第 9 章示例 9.4 所示。
interval 包中的函数 getsurv 还可用于获取每种治疗在指定时间的生存函数估计。由于生存函数在 Turnbull 区间内不是唯一定义的，因此该函数使用线性插值从连接 Turnbull 区间开始和结束处的生存估计的线上的点获取估计。以下代码捕获了每个治疗组在事件时间为 \(0,1, 2,\ldots,12\) 时的估计，并将这些值形成一个矩阵 surv_ests，用于打印和进一步分析。
> times<-c(0:12)
> icsurvests<-getsurv(times,icsurv)
> surv_ests<-cbind(times,icsurvests[[1]]$S,icsurvests[[2]]$S)
> colnames(surv_ests)<-c("Time","S(0)","S(1)")
> print(surv_ests)这会产生以下输出，该输出可以与从 summary(icsurv) 输出中的生存概率获得的生存函数估计进行核对和统一。
Time S(0) S(1)
[1,] 0 1.0000000 1.0000000
[2,] 1 0.9210526 0.9833333
[3,] 2 0.8421053 0.9666667
[4,] 3 0.8421053 0.9500000
[5,] 4 0.8421053 0.9333333
[6,] 5 0.8421053 0.9166667
[7,] 6 0.8421053 0.9166667
[8,] 7 0.8313090 0.7083333
[9,] 8 0.8205128 0.7083333
[10,] 9 0.8097166 0.7083333
[11,] 10 0.7989204 0.7083333
[12,] 11 0.7881242 0.7083333
[13,] 12 0.7773279 0.7083333icfit 对象 icsurv 也可直接使用 plot(icsurv) 绘制，这会产生类似于图 9.2 所示的图形。
图 9.2
使用 icenReg 包中的 R 函数 ic_np 也可以得到区间删失数据生存函数的非参数最大似然估计。要使用此函数，必须按以下方式指定区间上下限的矩阵，而不是 Surv 对象。
summary 函数的输出内容较少，但 ic_np 对象 icsurv1 中包含了 icsurv1$scurves，可以通过 print(icsurv1$scurves) 进行展示，其产生的输出结果就包括了生存函数估计。
17.10.2 区间删失数据的半参数模型
第 9 章 9.3 节描述的半参数 Turnbull 模型可以使用 icenReg 包中的 ic_sp 函数拟合区间删失数据。例如，要为溃疡复发时间数据拟合仅包含治疗效应的模型，请使用36
选项 bs_samples= 指定在计算拟合模型中参数估计的标准误差的 bootstrap 估计时使用的 bootstrap 样本数。
第 9 章的示例 9.5 概述了此过程。可以使用 summary(icfit) 获得输出并给出
Model: Cox PH
Dependency structure assumed: Independence
Baseline: semi-parametric
Call: ic_sp(formula = Surv(left, right, type = "interval2") ~ treatment,
data = ulcer, bs_samples = 10000)
Estimate Exp(Est) Std.Error z-value p
treatment 0.1953 1.216 1.502 0.13 0.8966
final llk = -31.44204
Iterations = 6
Bootstrap Samples = 10000根据该输出，接受治疗 B 的患者相对于接受治疗 A 的患者的对数风险比估计为 0.195，标准误为 1.488. 该标准误差是使用 bootstrap 模拟获得的，因此不能完复现。\(-2\log\hat{L}\) 统计量的值根据 final llk 的值计算为 \(−2(−31.442) = 62.884\)，如示例 9.5 所示。
getSCurves 函数可用于获取每个治疗组的生存函数估计
这将产生以下的简略输出。
$Tbull_ints
lower upper
[1,] 1.0e-10 1e-10
[2,] 1.0e-10 2e+00
[3,] 6.0e+00 7e+00
[4,] 1.2e+01 Inf
$S_curves
$S_curves$`1`
[1] 1.0000000 0.8957260 0.7585258 0.0000000
$S_curves$`2`
[1] 1.0000000 0.8747077 0.7146427 0.0000000
attr(,"class")
[1] "sp_curves"Turnbull 区间和生存函数估计与第 9 章表 9.4 给出的值相对应。请注意，getSCurves 函数的默认设置是根据拟合模型中变量值的均值来估计生存函数。因此，getSCurves(icfit) 给出了 Treat=0.5881 的生存函数（因为 0.5881 是该变量在 43 名患者中的样本均值），但这没有意义。
表 9.4
17.10.3 区间删失数据的参数模型
利用 survival 包中的 survreg 函数，可以对区间删失数据拟合参数模型。该函数中使用的 Surv 对象与其他用于处理区间删失数据的函数所使用的相同，并且其输出结果与 17.6 节描述的使用survreg拟合参数模型时的情况相似。唯一需要注意的是，对于左删失时间为零的情况，必须将其替换为 NA。例如，要对溃疡复发时间的区间删失数据拟合 Weibull 模型，请使用以下命令：
> ulcer$left1<-ifelse(ulcer$left==0,NA,ulcer$left)
> icweib<-survreg(Surv(left1,right,type="interval2") ~ treatment,
data=ulcer)在计算逻辑表达式以给出 True 或 False 的结果时，ifelse 语句中需使用双等号。有关解释此输出的详细信息，请参见 17.6 节。结果与第 9 章示例 9.7 中给出的结果一致。
17.11 脆弱模型
包含个体随机效应或组内个体共享的随机效应的生存时间模型称为脆弱模型。这类模型可能是参数化的或半参数化的，并在第 10 章有详细描述。
17.11.1 拟合具有个体脆弱性的参数脆弱模型
时间 \(t\) 时事件风险的比例风险模型，其中第 \(i\) 个个体（共 \(n\) 个）的随机效应为 \(u_i\) ，为
\[h_i(t)=\exp(\boldsymbol{\beta'x}_i+u_i)h_0(t)\]
其中 \(\boldsymbol{\beta}^{\prime}\boldsymbol{x}_i\) 是解释变量 \(\boldsymbol x_i\) 值与系数 \(\boldsymbol\beta\) 的线性组合，\(h_0(t)\) 是基线风险函数。该模型也可写为
\[\begin{align} h_i(t)=z_i\exp(\boldsymbol{\beta'x}_i)h_0(t) \tag{17.2} \end{align}\]
其中 \(z_i = \exp(u_i)\) 称为脆弱性。尽管某些 R 包中使用脆弱性一词指代 \(u_i\) 或 \(z_i\)，但在后续讨论中，我们将保持随机效应 \(u_i\) 与其对应的脆弱项 \(z_i\) 之间的这种区分。
参数脆弱模型可以使用 parfm 包中的 R 函数 parfm 进行拟合。该函数可为假定有有指数, Weibull, loglogstic 或 lognormal 分布的数据进行建模，其中脆弱性分布包括 gamma 和 lognormal. parfm 函数中的基线分布以与 survreg 函数不同的方式进行参数化。具体来说，Weibull 基线风险表示为 \(\text{λρtρ-1}\)，因此 \(\rho\) 代替了第 5 章式 (5.4) 中的 \(\gamma\)，指数基线风险为 \(\rho=1\)。loglogstic 基线风险的参数化如第 5 章式 (5.8) 所示，但使用 \(\alpha\) 替代了 \(\theta\)，lognormal 分布如第 5 章 5.1.4 节所述。脆弱性分布的参数为 \(\theta\)，它是第 \(i\) 个个体的模型中脆弱项 \(z_i\) 的方差。因此，在 10.2 节描述的 gamma 脆弱性分布中的 \(\theta\) 替换为 \(\theta-1\)，而第 10 章 10.2.1 节中的 lognormal 脆弱性方差 \(\text{var }(z_i) = \theta\)，因此脆弱效应 \(u_i\) 的方差为 \(e^{\theta}(e^\theta − 1)\)。有关 parfm 中使用的模型范围和参数化的详细信息，请参阅与软件包关联的 R vignette (Munda, Rotolo and Legrand, 2017).
为了说明 parfm 函数的使用，将使用示例 10.3 中关于等待肺移植的患者的生存数据。首先将数据输入到 R 数据框，将性别声明为一个因子，其水平对应于男性和女性患者，而疾病为四个水平的因子，对应于 COPD、纤维化、化脓性疾病和其他疾病。基线水平设为患有化脓性疾病的男性患者。
> lung<-read.table("Survival of patients registered for a lung transplant.dat",
header=TRUE)
> lung$fdisease<-factor(lung$disease)
> levels(lung$fdisease)<-list("COPD"=1,"Fib"=2,"Supp"=3,"Other"=4)
> lung$fdisease<-relevel(lung$fdisease,ref=3)
> lung$fgender<-factor(lung$gender)
> levels(lung$fgender)<-list("male"=1,"female"=2)
> lung$fgender<-relevel(lung$fgender,ref=2)函数 parfm 的语法与生存数据建模中使用的其他 R 函数的语法相似，模型公式是使用 Surv 对象表示的。此外，dist= 选项用于指定生存时间的参数分布，cluster= 选项指定用于脆弱性的项，该项视为一个因子，而 frailty= 选项用于指定脆弱分布。为了拟合具有 gamma 脆弱性的 Weibull 模型，其中 196 名患者中每个患者都有单独的脆弱性值，使用以下代码。
> f_lung0<-parfm(Surv(time,status) ~ fgender+age+bmi+fdisease,
cluster="patient",dist="weibull",frailty="gamma",data=lung)不幸的是，模型拟合程序无法收敛，因为 Weibull 尺度参数的估计变得非常小。这个问题可以通过使用不同的生存时间时间尺度来解决。在这里，在除以 30.44 将生存时间从天转换为月后，模型拟合程序便收敛了。通过这种重新调整，原始 Weibull 基线风险变为 \(30.44h_0(t/30.44)\)，即 \(30\lambda\rho(t/30.44)^{\rho-1}\)，因此为以月为单位的生存时间拟合的模型的尺度参数估计为 \(\lambda^*=30.44^\rho\lambda\)。因此，原始尺度上的 Weibull 尺度参数为 \(\lambda^{*}/30.44^{\rho}\)。对数似然统计量 \(−2 \log \hat L\) 的值也会有所不同，但替代模型之间的该值之差不受生存时间的任何缩放的影响。
为了将时间尺度从天改为月，并拟合具有 gamma 脆弱性的 Weibull 模型，使用以下代码。
> lung$time1<-lung$time/30.44
> f_lung<-parfm(Surv(time1,status) ~ fgender+age+bmi+fdisease,
cluster="patient",dist="weibull",frailty="gamma",data=lung)print(f_lung) 得到如下输出。
Frailty distribution: gamma
Baseline hazard distribution: Weibull
Loglikelihood: -591.709
ESTIMATE SE p-val
theta 3.015 1.063
rho 1.304 0.229
lambda 0.023 0.025
fgendermale 0.433 0.412 0.294
age 0.003 0.020 0.867
bmi -0.014 0.047 0.763
fdiseaseCOPD -1.008 0.692 0.145
fdiseaseFib 1.383 0.722 0.055 .
fdiseaseOther -0.076 0.766 0.921
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1根据该输出，在比较嵌套模型时使用的 \(−2 \log \hat L\) 统计量为 \(−2(−591.709)=1183.418\)，这也可以通过 -2*logLik(f lung) 来计算。此外，theta 是指脆弱性方差，即式 (17.2) 中 \(z_i\) 的方差估计，rho 和 lambda 指定了 Weibull 基线风险中的形状和尺度参数。请注意，第 10 章表 10.3 中 \(\hat\lambda\) 的值为 \(0.023/(30.441.304)=0.00027\)。这些估计值后面是模型中项的系数估计，它与表 10.3 第四列中的条目一致。然而，由于生存时间进行了缩放，\(−2 \log \hat L\) 的值有所不同。可以使用 predict 函数获得脆弱效应估计，即 predict(f_lung) 给出 gamma 脆弱的 196 个单独估计 \(\hat z_i\)。通过将 frailty= 选项更改为 frailty="lognormal"，可得到表 10.3 第三列中的条目。
表 10.3
17.11.2 拟合具有共享脆弱性的参数脆弱模型
函数 parfm 也可用于拟合具有共享脆弱性的参数模型。其代码与个体脆弱性的代码非常相似，只是 cluster= 选项现在指定与具有共享脆弱性的个体组相关的变量。
作为说明，请考虑示例 10.4 中有关肾移植后的生存数据。在这里，一名肾脏供体可能会提供一或两个移植，因此两个受体可能有一个共同的供体，从而导致共享脆弱。以下代码用于输入数据、变换时间尺度并拟合具有 gamma 脆弱性的 Weibull 模型。
> tplant<-read.table("Survival following kidney transplantation.dat",
header=TRUE)
> tplant$time1<-tplant$time/30.44
> f_kidney<-parfm(Surv(time1,status) ~ age+cit+diabetes,
cluster="donor",dist="weibull",frailty="gamma",data=tplant)使用 print(f_kidney) 得到以下输出。
Frailty distribution: gamma
Baseline hazard distribution: Weibull
Loglikelihood: -433.31
ESTIMATE SE p-val
theta 0.421 0.620
rho 0.602 0.068
lambda 0.005 0.002
age 0.019 0.008 0.013 *
cit 0.024 0.020 0.229
diabetes -0.258 0.458 0.573
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Kendall's Tau: 0.174 该输出的解释与个体脆弱模型的解释相同。然而，由于生存时间的缩放，对数似然值再次与示例 10.4 中引用的值不同。gamma 脆弱的拟合模型的输出还包括 Kendall’s \(\tau\) 统计量，这是一种生存时间对之间关联的度量，参见第 10 章 10.8 节。类似的代码可以用于拟合参数分布和脆弱性分布的其他组合。
17.11.3 拟合具有个体 lognormal 脆弱性的半参数模型
使用同名包中的 coxme 函数可以为生存数据拟合具有 lognormal 脆弱项的模型。该函数使用惩罚偏似然法来拟合包含脆弱效应的模型，如第 10 章 10.5 节所述。
再次考虑等待肺移植期间生存时间的数据。第 17.11.1 节使用的 lung 数据框架包含以性别和疾病集为因子的数据。此分析将使用以天为单位的生存时间。为了拟合 Cox 回归模型，该模型包含 \(Age,Sex,BMI\) 和 \(Disease\)，以及假设具有 lognormal 分布的患者的个体脆弱效应 \(z_i\)，使用如下代码来完成。
在该代码中，通过在模型公式中添加 (1|patient) 来指定脆弱项，并且脆弱变量自动假定为一个因子。使用 summary(lnfcox) 可得到以下输出。
Mixed effects coxme model
Formula: Surv(time, status) ~ fgender + age + bmi + fdisease + (1 | patient)
Data: lung
events, n = 123, 196
Random effects:
group variable sd variance
1 patient Intercept 0.7227754 0.5224043
Chisq df p AIC BIC
Integrated loglik 17.75 7.00 1.314e-02 3.75 -15.93
Penalized loglik 116.62 49.59 2.532e-07 17.44 -122.02
Fixed effects:
coef exp(coef) se(coef) z p
fgendermale 0.11401 1.12076 0.23180 0.49 0.623
age 0.01209 1.01217 0.01312 0.92 0.357
bmi -0.02555 0.97477 0.02839 -0.90 0.368
fdiseaseCOPD -0.68092 0.50615 0.42471 -1.60 0.109
fdiseaseFib 0.55321 1.73882 0.39815 1.39 0.165
fdiseaseOther -0.19793 0.82043 0.44819 -0.44 0.659该输出的第一部分给出了三个对数似然值，展示为 NULL, Integrated 和 Fitted。零对数似然是指不包含解释变量或脆弱性的模型。这只是 Cox 回归模型的最大化偏对数似然，该模型对每个个体具有相同的基线风险函数 \(h_0(t)\)。积分对数似然是第 10 章式 (10.28) 给出的最大化边际对数似然 \(\log L_m(\hat{\boldsymbol{\beta}},\hat{\sigma}_u^2)\)。拟合对数似然是从式 (10.25) 获得的拟合模型的最大偏对数似然 \(\log L_p(\boldsymbol{\hat{\beta}},\hat{\boldsymbol{u}})\)。在这三个对数似然值中，积分对数似然是最有用的，因为它用于比较替代模型。
输出的下一部分中 Integrated loglik 所在行 给出了卡方值，用于将积分对数似然与零模型的对数似然进行比较。即，\(2(587.616−578.74)=17.75\)。标记为 Penalized loglik 的行将拟合或最大化偏对数似然与零模型的似然进行比较，即 \(2(587.616−529.304)=116.62\)。自由度是指拟合模型中独立参数的有效数。最大化偏对数似然包括随机效应的估计 \(\hat {\boldsymbol u}\)，而它们已经从边际对数似然性函数中积分出来。这解释了与两个检验相关的有效自由度数的差异。标记为 AIC 和 BIC 值的列中的 AIC 和 BIC 统计量是根据 Chisq 值计算通过减去自由度的两倍或 \(q \log d\)，其中 \(q\) 是自由度，\(d\) 是事件数。当关注于替代模型的比较时，这些统计量的用途有限。
输出的下一部分给出了参数估计、风险比、标准误、\(z\) 以及检验系数为零的假设的 \(P\) 值。接下来是对模型中随机效应项的标准差和方差的估计，即 \(\widehat{\text{var }}(u_i)\)。
一旦拟合了脆弱模型，就可以使用 coxme 对象来提取一些量的值。例如，如 10.4 节所述，lnfcox$wear 给出了个体脆弱性估计，可用于识别更脆弱的个体的群体。当使用 lnfcox$loglik 来产生最大化对数似然时，可获得以下输出。
这里，Penalized log-likelihood 是拟合模型的惩罚偏对数似然的最大值，即式 (10.27) 中的 \(\log\hat{L}_{pen}(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{u}},\hat{\sigma}_u^2)\)。这是通过从最大化偏对数似然 \(\log L_p(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{u}})\) 中减去惩罚项 \((2\hat{\sigma}_u^2)^{-1}\sum_i\hat{u}_i^2\) 来计算的。该惩罚函数的值可使用 lnfcox$penalty 得出，为 22.8272，因此，如所给出的，最大化惩罚对数似然为 \(−529.3037−22.8272=−552.1309\)。
17.11.4 拟合具有个体 gamma 脆弱性的半参数模型
具有 gamma 脆弱性分布的模型不能使用 coxme 拟合，但 survival 包中的 frailty 函数可与 coxph 一起使用，以纳入 gamma 脆弱效应。下面给出了为登记进行肺移植的患者的生存时间数据拟合具有 gamma 脆弱性的 Cox 回归模型的 R 代码。
> gfcox<-coxph(Surv(time,status) ~ fgender+age+bmi+fdisease+
frailty.gamma(patient),ties="breslow",data=lung)在这里，通过将 frailty 函数添加到模型公式中来指定脆弱项，summary(gfcox) 将产生以下输出。
Call:
coxph(formula = Surv(time, status) ~ fgender + age + bmi + fdisease +
frailty.gamma(patient), data = lung, ties = "breslow")
n= 196, number of events= 123
coef se(coef)
fgendermale 0.30599 0.35412
age 0.00557 0.01812
bmi -0.02430 0.04257
fdiseaseCOPD -0.88837 0.60317
fdiseaseFib 1.22698 0.60721
fdiseaseOther -0.08851 0.66758
frailty.gamma(patient)
se2 Chisq
fgendermale 0.23434 0.75
age 0.01198 0.09
bmi 0.02786 0.33
fdiseaseCOPD 0.40313 2.17
fdiseaseFib 0.39557 4.08
fdiseaseOther 0.43458 0.02
frailty.gamma(patient) 243.82
DF p
fgendermale 1.0 3.9e-01
age 1.0 7.6e-01
bmi 1.0 5.7e-01
fdiseaseCOPD 1.0 1.4e-01
fdiseaseFib 1.0 4.3e-02
fdiseaseOther 1.0 8.9e-01
frailty.gamma(patient) 120.5 2.2e-10
exp(coef) exp(-coef)
fgendermale 1.3580 0.7364
age 1.0056 0.9944
bmi 0.9760 1.0246
fdiseaseCOPD 0.4113 2.4312
fdiseaseFib 3.4109 0.2932
fdiseaseOther 0.9153 1.0925
lower .95 upper .95
fgendermale 0.6784 2.718
age 0.9705 1.042
bmi 0.8979 1.061
fdiseaseCOPD 0.1261 1.342
fdiseaseFib 1.0376 11.213
fdiseaseOther 0.2474 3.387
Iterations: 6 outer, 68 Newton-Raphson
Variance of random effect= 2.154654 I-likelihood = -576.7
Degrees of freedom for terms= 0.4 0.4 0.4 1.4 120.5
Concordance= 0.944 (se = 0.008 )
Likelihood ratio test= 304.5 on 123.2 df, p=<2e-16该输出给出了模型中项的系数估计及其标准误、计算为 \(\{(\hat{\beta}/\sec{(\hat{\beta})}\}^2\) 的卡方统计量以及相应的 \(P\) 值。标记为 se2 的量是对标准误的稳健估计，该标准误解释了模型指定中的不确定性。系数的指数，即风险比，也与 95% 的区间估计一起给出。
标记为 frailty.gamma(patient) 的行给出了无脆弱性假设的卡方检验统计量。概括地说，这基于计算为 \(\hat{\boldsymbol{u}}^{\prime}\boldsymbol{V}^{-1}(\hat{\boldsymbol{u}})\hat{\boldsymbol{u}}\) 的 Wald 检验统计量，其中 \(\boldsymbol{V}(\hat{\boldsymbol{u}})\) 为向量 \(\widehat{\boldsymbol{u}}\) 中脆弱性估计的方差-协方差阵，这是从惩罚偏对数似然的二阶导数矩阵中获得的，如附录 A 所述。
给出拟合模型所需迭代次数的信息后，给出了随机效应的方差估计，即 \(\widehat{\text{var}} \left ( u _ i \right )\)。同一行中是拟合模型的积分或边际对数似然，标记为 I-likelihood。将其乘以 −2 后，即为用于将替代模型与伽玛脆弱性进行比较的 \(−2 \log \hat L\) 统计量。似然函数的值包括第 10 章 10.5 节末尾提到的额外项 \(\sum_{i=1}^n\delta_i\)，因此可以更容易地将具有脆弱性的模型与没有脆弱性的模型进行比较。之后给出了模型中每个项（包括随机效应）的有效自由度数，然后是一致性估计，该一致性通过第 3 章 3.12.1 节中描述的 \(c\) 统计量来度量。标记为 Likelihood ratio 的统计量将拟合模型与零模型进行比较，其中拟合模型的似然是 在 \(\beta\) 和脆弱效应 \(u_i\) 估计处的最大化偏对数似然，但不包含惩罚项，即 \(\log L_p(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{u}})\)，以第 10.5 节表示法。
以 frailty.gaussian(patient) 形式的 frailty 函数也可以用于向模型中添加 lognormal 脆弱性。尽管通常认为 coxme 函数更可靠，但这两种方法往往产生相似的结果。
17.11.5 拟合具有共享脆弱性的半参数模型
R 函数 coxme 以及带 frailty.gamma 的 coxph 可分别用于对共享 lognormal 和 gamma 脆弱性进行建模。该代码与上一节中描述的个体对若行的代码非常相似。现用肾移植后的生存数据，以说明参数共享脆弱模型。以下代码使用函数 coxme 来拟合供体效应的 lognormal 脆弱性
summary(slnfcox) 得到以下输出。
Mixed effects coxme model
Formula: Surv(time, status) ~ age + cit + diabetes + (1 | donor)
Data: tplant
events, n = 71, 434
Random effects:
group variable sd variance
1 donor Intercept 0.5268985 0.277622
Chisq df p AIC BIC
Integrated loglik 5.44 4.00 0.245290 -2.56 -11.61
Penalized loglik 42.36 20.92 0.003677 0.51 -46.84
Fixed effects:
coef exp(coef) se(coef) z p
age 0.01876 1.01894 0.01037 1.81 0.0706
cit 0.02357 1.02385 0.02177 1.08 0.2791
diabetes -0.25418 0.77555 0.44938 -0.57 0.5716该输出与使用 coxme 拟合的个体脆弱性模型的格式相同，并且与第 10 章示例 10.4 中给出的结果一致。
带 frailty.gamma 函数的函数 coxph 用于拟合具有共享 gamma 脆弱性的模型，使用
> sgfcox<-coxph(Surv(time,status) ~ age+cit+diabetes+
frailty.gamma(donor),ties="breslow",data=tplant)模型中解释变量系数的估计结果与使用具有 lognormal 脆弱性的半参数模型时发现的结果非常相似。
17.12 竞争风险
竞争风险模型用于研究中的个体可能经历多个不同终点之一的情况，并在第 12 章进行了描述。在本节中，使用 R 来估计累积发生率函数。并拟合原因别风险和原因别发生率的模型。第 12 章示例 12.1 中描述的肝移植后移植失效时间的数据将用于说明 R 代码和输出。在此数据集中，许多患者的移植失效时间删失，有四个潜在的致命终点，即排异、血栓形成、疾病复发和其他原因。在数据集中，与失效原因相关的变量分别编码为 0, 1, 2, 3 和 4. 数据被输入到 R 数据框架中，性别和肝病都被设为因子，女性和酒精性肝病作为参考水平，这些使用如下代码完成。
> liver<-read.table("Survival of liver transplant recipients.dat", header=TRUE)
> liver$fdisease<-factor(liver$disease)
> levels(liver$fdisease)<-list("PBC"=1,"PSC"=2,"ALD"=3)
> liver$fdisease<-relevel(liver$fdisease,ref=3)
> liver$fgender<-factor(liver$gender)
> levels(liver$fgender)<-list("male"=1,"female"=2)
> liver$fgender<-relevel(liver$fgender,ref=2)17.12.1 原因别的风险函数的估计和建模
可以按照第 12 章 12.2 节所述的方式，使用每种事件类型的生存函数来总结竞争风险数据。这包括依次使用每个失效原因的事件时间，将所有其他事件的时间视为删失，以提供原因别数 据。这可以通过创建一个新的变量来实现，该变量指示给定肝衰竭原因的事件时间是一个事件，而所有其他原因的事件时间都是删失的。例如，如果死因是血栓，则定义一个变量，对于 liver 数据框中变量 cof 等于 2 的观测，该变量为 1，否则为 0：
然后可以在 survfit 函数中使用变量 status2 来获得生存函数估计。例如，当血栓形成是移植失效的原因时，每种肝病的生存函数的 Kaplan-Meier 估计使用如下代码获得。
随后，可以使用 coxph 函数拟合每种事件类型的风险的 Cox 回归模型，如 17.4 节所述。例如，要获得第 12 章表 12.5 中的原因别风险比，在以血栓形成引起的移植失败为事件拟合 Cox 回归模型后，使用
所得血栓形成的风险比如表 12.5 所示。移植失效的其他原因可以用类似的方式处理。
表 12.5
17.12.2 估计累计发生率函数
cmprsk 包中的函数 cuminc 可用于估计原因别累积发生率函数，其中选项 ftime= 和 fstatus= 用于指定包含生存时间和状态指示符的变量。例如，为了获得血栓形成的总体累积发生率，使用以下命令。
这可以扩展到为每个移植指示符提供累积发生率函数，使用 group= 选项指定分组变量，如
在此之后，例如，酒精性肝病的累积发病率和相应的时间点可以从变量 cuminc2$"ALD 1"$est 和 cuminc2$"ALD1"$time 中获得。这些可以分配变量名称，以便在进一步的分析或图形总结中使用。例如，可以使用plot(cuminc2) 获得累积发生率函数的图形，这类似于图 12.2.
17.12.3 累积发生率的 Fine and Gray 模型
第 12 章 12.5 节中描述的 Fine and Gray 模型用于在存在竞争风险的情况下对事件的累积发生率进行建模。该模型可以使用 cmprsk 包中的函数 crr 进行拟合。该函数不允许指定模型公式，相反，必须提供一个包含解释变量值的矩阵，该矩阵将被包含在模型中。对于因子，需要定义指示变量，但可使用 model.matrix 函数来辅助。例如，函数 model.matrix(~A+B)[,-1] 将为两个因子 A 和 B 生成指示变量。最后的 [,-1] 会从函数的输出中移除一个常数项，以避免对指示变量进行过度参数化。
例如，对于肝移植后移植失败时间的数据，性别和肝病因素的指标变量的水平在 liver 数据框的 fgender 和 fdisease 列中，使用如下命令。
患者的年龄也应包括在模型中，因此必须将其与 indvars 中的指标变量相结合，以给出包含年龄以及性别和疾病的指标变量的矩阵 xvars。生成的矩阵的列也被重命名。这些使用如下代码实现。
矩阵 xvars 的前 10 行如下所示。
age PBC PSC male
1 55 0 0 1
2 63 1 0 1
3 67 1 0 1
4 58 0 0 1
5 59 1 0 1
6 35 1 0 0
7 51 1 0 0
8 61 0 1 0
9 51 0 0 0
10 59 0 0 0由 model.matrix 函数生成的指标变量是这样的：如果疾病是 PBC（或 PSC），则 PBC（或 PSC）的值为 1，否则为 0；对于男性患者，male 的值为 1，否则为 0.
现在可用 crr 函数来拟合特定事件的累积发生率的模型。使用选项 ftime= 和 fstatus= 指定包含每个事件时间以及每种事件类型编码以及删失观测的变量。接下来，使用 cov1= 指定要包含在模型中的解释变量矩阵。可以使用 cencode= 选项指定删失时间的编码，failcode= 选项用于指定要建模的事件类型，所有其他事件都是竞争风险。在 liver 数据框中，变量 cof 包含事件编码，对于删失时间为 0，对于血栓形成为 2. 使用 failcode=2 指定要对血栓形成的发生率进行建模。选项 cencode=0 指定变量 cof 在删失时间为 0，但这是不必要的，因为这是默认值。
使用了以下 crr 函数调用，以对血栓形成发生率关于患者年龄、性别和原发性肝病的依赖性进行建模。
cifit2<-crr(ftime=liver$time,fstatus=liver$cof,cov1=xvars,cencode=0,
failcode=2)使用 summary(cifit2) 会生成以下输出。
Competing Risks Regression
Call:
crr(ftime = liver$time, fstatus = liver$cof, cov1 = xvars, failcode = 2,
cencode = 0)
coef exp(coef) se(coef) z p-value
age -0.0233 0.977 0.0104 -2.246 0.025
PBC 0.6057 1.833 0.2798 2.165 0.030
PSC 0.4648 1.592 0.3086 1.506 0.130
male 0.0760 1.079 0.2458 0.309 0.760
exp(coef) exp(-coef) 2.5% 97.5%
age 0.977 1.024 0.957 0.997
PBC 1.833 0.546 1.059 3.171
PSC 1.592 0.628 0.869 2.914
male 1.079 0.927 0.666 1.747
Num. cases = 1761
Pseudo Log-likelihood = -528
Pseudo likelihood ratio test = 9.28 on 4 df,输出在很大程度上是不言自明的，并给出了模型中解释变量系数的估计、相应的子风险比 exp(coeff)、系数的标准误、\(z\) 值和 \(P\) 值。然后重复输出了子风险比，标记为 2.5% 和 97.5% 的两列是它们的 95% 置信限。这些对应于表 12.6 中 “Thrombosis”（血栓形成）栏中的值。随后是拟合模型的对数似然值。Pseudo Log-likelihood 是根据第 12 章式 (12.12) 获得的最大化偏对数似然，Pseudo likelihood ratio test 将其与零模型（即没有解释变量的模型）进行比较。
可以从 crr 对象 cifit2 中提取各种量以用于进一步分析。使用值为 1, 3 和 4 的 failcode 将能够对其他三种事件类型的发生率进行建模。
17.13 多次事件和事件历史建模
复发事件模型，具体来说，Anderson and Gill (AG) 以及 Prentice, Williams and Peterson (PWP) 模型，可以在以计数过程格式表达相关数据后进行拟合，如第 13 章示例 13.2 所示。这可以使用 17.9 节中描述的 tmerge 和 surveSplit 函数直接完成。此外，tmerge 中的 cumtdc= 选项可用于枚举每个个体的复发事件数。当使用 coxph 函数拟合 PWP 模型时，需要通过复发次数进行分层。
事件历史模型是通过定义与不同转移相关的变量来拟合的，并使用该变量以 tmerge 的计数过程格式形成数据。将数据集转换为所需格式所需的数据操作可能很复杂，这取决于所提供数据的格式。因此，此处不再提供更多详细信息。
17.14 相依删失
14.3 节描述的具有相依删失的 Cox 回归模型可以通过以下方式进行拟合。首先，通过为删失时间拟合生存模型来估计删失概率。这是通过反转通常的事件时间指示符来完成的，使得它对删失时间取值 1，对事件时间取值 0. 如 14.3 节所述，拟合参数模型更方便，因此 survreg 函数用于对删失的风险进行建模。该模型中的参数估计产生线性预测器和尺度参数，它们被保存到数据框中以供以后使用。
然后使用 surveSplit 函数以计数过程格式表示事件时间数据。这是通过创建一个具有唯一事件时间的变量来实现的，该变量用于函数的 cuts= 选项。根据扩展的数据框，在每个区间终点计算删失时间的生存函数的估计。它们的倒数是相依删失建模中所需的权重。最后，使用具有选项 weights=（用于指定权重变量）的 coxph 函数为计数过程数据拟合 Cox 回归模型。
为了说明所需代码，将使用示例 14.1 中关于等待肝移植时死亡时间的数据来拟合考虑相依删失的 Cox 回归模型。首先输入数据，以给出与表 14.1 匹配的数据帧 livertx。
表 14.1
然后，使用删失状态变量 cstatus，通过为删失时间拟合 Weibull 分布来对删失概率进行建模。使用以下代码将线性预测器和尺度参数的估计分配给 lp 和 sigma 并保存在 livertx 数据框中。
> livertx$cstatus<-1-livertx$status
> wcens<-survreg(Surv(time,cstatus) ~ age+gender+bmi+ukeld,
data=livertx)
> livertx$lp<-wcens$linear.predictors
> livertx$sigma<-wcens$scale下一步是构造一个变量 survtime，其中包含数据集中事件时间的唯一值，并在 survSplit 函数中使用该变量以计数过程格式给出新的数据框 livertx1。
> survtime<-subset(livertx,status==1)
> survtime<-unique(survtime[,"time"])
> livertx1<-survSplit(Surv(time,status) ~ .,cut=survtime,data=livertx)变量 time 现在是每个时间区间的终点。现在可以根据式 (17.1) 中参数模型的对数线性形式找到在这些端点处的生存函数估计。根据第 14 章的式 (14.5)，超过时间 \(t\) 的删失概率估计为
\[\hat{S}_{ci}(t)=\exp\left\{-\exp\left[\frac{\log t-\hat{\eta}_{ci}}{\hat{\sigma}_{c}}\right]\right\}\]
其中，模型中的线性预测器是 \(\hat{\eta}_{ci}=\hat{\mu}_{c}+\hat{\boldsymbol{\alpha}}_{c}'\boldsymbol{x}_{i}\)，并且 \(\hat{\boldsymbol{\mu}}_{c}\) 是截距，\(\hat{\mathbf{O}}_{c}\) 是解释变量的系数向量，而 \(\hat{\sigma}_{c}\) 是 Weibull 模型对数线性形式的尺度参数。用于在计数过程数据中的每个区间的终点处获得的生存函数估计并计算权重的代码如下。
> livertx1$csurv<-exp(-exp((log(livertx1$time)-
livertx1$lp)/livertx1$sigma))
> livertx1$weight<-1/livertx1$csurv不再需要变量 cstatus, lp, sigma 和 csurv，可以使用如下代码从 livertx1 中删除。
livertx1 的前 8 行如下所示。
patient age gender bmi ukeld tstart time status weight
1 1 60 0 24.24 60 0 1 0 1.004505
2 2 66 0 30.53 67 0 2 1 1.011361
3 3 71 0 26.56 61 0 2 0 1.008534
4 3 71 0 26.56 61 2 3 0 1.012929
5 4 65 1 23.15 63 0 2 0 1.012837
6 4 65 1 23.15 63 2 3 0 1.019468
7 5 62 0 22.55 64 0 2 0 1.011611
8 5 62 0 22.55 64 2 3 0 1.017603现在可使用如下代码拟合加权 Cox 回归模型。
> wcoxfit<-coxph(Surv(tstart,time,status) ~ age+gender+bmi+ukeld,
weight=weight,ties="breslow",data=livertx1)下面是 summary(wcoxfit) 的部分输出。
Call:
coxph(formula = Surv(tstart, time, status) ~ age + gender + bmi +
ukeld, data = livertx1, weights = weight, ties = "breslow")
n= 8968, number of events= 64
coef exp(coef) se(coef) robust se z Pr(>|z|)
age 0.011805 1.011875 0.007758 0.031582 0.374 0.70858
gender -0.990764 0.371293 0.367368 0.655330 -1.512 0.13057
bmi -0.021822 0.978414 0.014422 0.049145 -0.444 0.65702
ukeld 0.155906 1.168716 0.018695 0.042688 3.652 0.00026 ***该输出中四个解释变量的系数估计及其稳健的标准误与表 14.4 中的加权估计一致。该模型的未加权版本可以使用如下代码实现。
表 14.4
如果需要，可以通过为每个单独的观测定义聚类标识符并将 cluster(id) 添加到模型公式中来获得该模型中系数标准误的稳健估计，如下所示。
> livertx1$id<-1:nrow(livertx1)
> coxfit<-coxph(Surv(tstart,time,status) ~ age+gender+bmi+ukeld+
cluster(id),ties="breslow",data=livertx1)该代码的输出给出了表 14.4 中的未加权估计。
17.15 贝叶斯生存分析
rstanarm 包支持使用 Markov
Chain Monte Carlo (MCMC) 方法和 R 建模语法来拟合各种回归模型。该包又是 rstan 包的补充，它提供了 Stan 的 R 接口，Stan 是一个用于贝叶斯建模和推理的 C++ 库。第 16 章描述的贝叶斯生存分析可以使用 rstanarm 中的函数 stan_surv 来实现，该函数在 Brilleman et al. (2020) 的文章中进行了描述。第 16 章示例 16.7 中描述的健康评估和与初级保健联系的 HELP 研究数据将用于展示如何在贝叶斯生存分析中使用函数 stan_surv，并说明结果输出。使用如下命令将数据输入到名为 help 的数据框。
安装 rstanarm 包后，可以使用 library("rstanarm") 加载它。
安装说明
rstanarm 最新发布的 CRAN 版本（version 2.32.1，发布日期：2024-01-18）尚无 feature/survival，即无法进行生存分析。要完成下列分析，需使用如下命令安装 rstanarm survival branch 的二进制版本（静态版本）：
install.packages("rstanarm", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))具体详见 rstanarm 的 Github页面 以及这里和这里的讨论。
17.15.1 贝叶斯参数建模
stan_surv 函数可用于拟合基线风险函数的贝叶斯版本的指数, Weibull 和 Gompertz 模型（如第 5 章所述）以及对数基线风险的 B 样条模型（如第 6 章的 6.3 节所述）。需要指定模型中未知参数的分布以及有关 MCMC 仿真过程的信息。贝叶斯 Weibull 比例风险模型包含数据集中的四个变量（\(Age, Gender, Housing\) 和 \(Help\)），可以使用以下代码为 HELP 研究的数据进行拟合。
> hweib<-stan_surv(Surv(days,status) ~ linkage+age+gender+housing,
basehaz="weibull",chains=1,iter=11000,warmup=1000,seed=123,
prior_intercept=normal(0,10),prior_aux=normal(0,10),
prior=normal(0,1),data=help)该模型是使用 Surv 对象以通常的方式指定的，而 Weibull 基线危险是使用 basehaz="weibull" 指定的。使用 stan_surv，Weibull 比例风险模型参数化为 \(h_i(t)~=~\exp(\boldsymbol{\beta'x}_i)\lambda\gamma t^{\gamma-1}\)，其中 \(\log \lambda\) 称为截距并在输出中表示为 (Intercept)，\(\gamma\) 为形状参数，在输出中标记为 weibull-shape。使用 algorithm=c("sampling") 启动 MCMC 采样过程，但由于这是 stan_surv 函数的默认设置，因此已省略此选项。关于要使用的链的数量、所需的模拟总数以及用于热身的模拟数的详细信息分别使用选项 chain=, iter= 和 warmup= 给出。seed= 选项可用于初始化模拟过程的伪随机数生成器，并使随机数序列能够再现。
接下来，指定先验分布。可使用广泛的先验分布，包括正态分布、柯西分布、指数分布和 Laplace 分布。Brilleman et al. (2020) 的 vignette 对这些进行了描述。在上述代码中，为 Weibull 模型的截距 \(\log \lambda\) 指定了具有均值为零和方差为 102 的正态分布。prior_aux= 选项用于指定形状参数的先验分布，该分布为非负分布。在这种情况下，prior_aux=normal(0,10) 指定 Weibull 形状参数 \(\gamma\) 为 half-normal 分布。prior= 选项用于指定模型线性部分中 \(\beta\) 系数的先验分布，这里它们都分配了 \(N(0,1)\) 分布。在第 16 章的示例 16.7 中使用了先验分布的这种选择。如果 \(\beta\) 系数具有不相同的先验分布，则可以使用不同的位置和尺度参数或不同的分布。
stan_surv 的函数调用得到以下输出。summary 和 print 函数提供不同的输出信息。summary(hweib) 的输出如下。
Model Info:
function: stan_surv
baseline hazard: weibull
formula: Surv(days, status) ~ linkage + age + gender + housing
algorithm: sampling
sample: 10000 (posterior sample size)
priors: see help('prior_summary')
observations: 447
events: 167 (37.4%)
right censored: 280 (62.6%)
delayed entry: no
Estimates:
mean sd 10% 50% 90%
(Intercept) -6.3613 0.5235 -7.0404 -6.3580 -5.6853
linkage 1.6422 0.1846 1.4043 1.6400 1.8778
age 0.0239 0.0100 0.0111 0.0239 0.0367
gender 0.3439 0.1945 0.0985 0.3390 0.5946
housing -0.1458 0.1530 -0.3403 -0.1456 0.0474
weibull-shape 0.6084 0.0414 0.5549 0.6074 0.6621
MCMC diagnostics
mcse Rhat n_eff
(Intercept) 0.0063 0.9999 6971
linkage 0.0021 0.9999 7710
age 0.0001 1.0000 9042
gender 0.0022 1.0000 8033
housing 0.0017 0.9999 8243
weibull-shape 0.0005 1.0004 6472
log-posterior 0.0244 1.0000 4865
For each parameter, mcse is Monte Carlo standard error, n_eff is a
crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence Rhat=1).该输出给出了 10000 次 MCMC 模拟中解释变量系数的均值和标准差，Weibull 截距和形状参数，以及它们分布的 10%, 50% 和 90% 百分位数。输出中对标记为 MCMC diagnostics 的输出内容进行了有益的描述。
使用 print(hweib) 的简略输出如下。
Median MAD_SD exp(Median)
(Intercept) -6.3580 0.5191 NA
linkage 1.6400 0.1835 5.1552
age 0.0239 0.0100 1.0242
gender 0.3390 0.1911 1.4036
housing -0.1456 0.1529 0.8645
weibull-shape 0.6074 0.0408 NA该输出给出了 MCMC 估计的中位数，以及中位绝对偏差 (median absolute deviation)，标记为 MAD_SD。这是对变异性估计的稳健度量，计算为每个模拟值与总体中位数之差绝对值的中位数。 exp(Median) 列表示中位数的指数，可用作模型中每个变量的风险比估计。
函数 prior_summary 提供了一种有用的方法来获得 stan_surv 函数中使用的先验总结，得到以下输出。
Priors for model 'hweib'
------
Intercept
~ normal(location = 0, scale = 10)
Coefficients
~ normal(location = [0,0,0,...], scale = [1,1,1,...])
Auxiliary (weibull-shape)
~ half-normal(location = 0, scale = 10)这证实了 Weibull 模型的截距和形状参数以及变量系数的所选先验。
17.15.2 贝叶斯半参数建模
stan_surv 函数也可用于拟合基于 B 样条的对数风险函数的模型（在第 6 章的 6.3 节中描述）的贝叶斯版本。在使用零次 B 样条的特殊情况下，该模型等价于分段指数模型，其中在第 \(j(j=1,2,\ldots,m+1)\) 个时间区间中，基线风险函数为 \(h_0(t)=\exp(\alpha_j)\)，其中 \(m\) 是内部结数，即最短和最长事件时间之间的结的数量。然而，R 函数 stan_surv 包含截距项，因此实际拟合的基线风险为 \(h_0(t)=\exp(\mu+\alpha_j)\)。由于该模型中现在有 \(m+2\) 个未知参数，与 \(m+1\) 个时间区间相关，所以该模型被过度参数化，\(\alpha_1\) 被设置为零。在第一个时间区间中，基线风险为 \(e^\mu\)，在其余 \(m\) 个时间区间中，则为 \(\exp(\mu+\alpha_{j-1}),j=2,3,\ldots,m+1\)。那么，对于具有 \(m\) 个结的模型，将只有 \(m\) 个B 样条系数。
所需的代码与为 HELP 研究数据拟合参数模型的代码非常相似，只是 basehaz= 选项设置为 bs，代表 B 样条曲线。为了拟合分段指数模型的贝叶斯版本，使用选项 basehaz_ops=list(degree=0,df=5) 来表示 B 样条函数的阶数和基线危险变化处的结的位置。这里，选项 df=5 意味着在最短 (0) 和最长 (456) 事件时间之间有五个内部结，并且时间区间的终点选择为在事件时间的等距百分位数处。这导致在第 18, 32, 42, 74 和 136 天出现内部结。选项 df=5 可以替换为用于内部结的事件时间列表，因此 basehaz_ops=list(degree=0,knots=c(18,32,42,74,136)) 将提供类似的输出。为每个事件时间区间设置不同基线风险的模型将非常接近 Cox 回归模型，这可以通过将结数设置为不同事件时间数，或者在 knots= 选项中列出唯一的事件时间来进行拟合。
为了说明贝叶斯半参数建模，以下代码用于为 HELP 研究的数据拟合具有 5 个内部结的贝叶斯分段指数模型。所有四个解释变量均已包含在模型中。\(N(0, 1)\) 先验用于解释变量的每个系数，并且截距和每个 B 样条参数的先验分布取为 \(N(0, 10)\)。MCMC 模拟过程先进行 1,000 次热身，然后进行 10,000 次模拟。
> bayespe<-stan_surv(Surv(days,status) ~ linkage+age+gender+housing,
basehaz="bs",basehaz_ops=list(degree=0,df=5),
chains=1,iter=11000,warmup=1000,seed=1234,
prior_intercept=normal(0,10),prior_aux=normal(0,10),
prior=normal(0,1),data=help)使用 print(bayespe) 得到的参数估计的中位数以及相应的风险比，如下。
Median MAD_SD exp(Median)
(Intercept) -7.6055 0.5038 NA
linkage 1.5339 0.1819 4.6362
age 0.0219 0.0102 1.0222
gender 0.2586 0.1926 1.2951
housing -0.1165 0.1542 0.8901
b-splines-coef1 0.1457 0.2733 NA
b-splines-coef2 1.6314 0.2695 NA
b-splines-coef3 -0.1793 0.2723 NA
b-splines-coef4 -0.8054 0.2731 NA
b-splines-coef5 -1.9390 0.2688 NA通过 MCMC 模拟得到的模型中解释变量的系数估计与拟合 Weibull 模型的贝叶斯版本时得到的系数估计非常相似。截距 Intercept 是 \(\mu\) 的估计，标记为 b-splines-coef 的系数是对样条基函数的系数 \(\alpha_j,j=2,3,4,5,6\) 的估计。基线风险函数在第一区间为 \(\exp(\mu)\)，在其余五个区间为 \(\exp{\mu + \alpha_j}\)。这些数值需要从每个单独的 10,000 次 MCMC 模拟中推导出来，而不是通过该输出的中位数得出。
对于具有给定特征的个体，posterior_survfit 可用于从 MCMC 模拟中模型参数的单个估计中获取后验预测生存函数和风险函数。具体来说，基线函数可以使用如下代码获得。
basevals<-data.frame(linkage=0,age=0,gender=0,housing=0)
basesurv<-posterior_survfit(bayespe,type="surv",newdata=basevals)
basehaz<-posterior_survfit(bayespe,type="haz",newdata=basevals)使用 plot 对它们绘图，例如 plot(basesurv)。
使用 waic 函数可以轻松得到使用 stan_surv 函数拟合的模型用于贝叶斯模型比较的 \(WAIC\) 统计量，如第 16 章 16.8 节所述。例如，waic(bayespe) 给出
其中拟合分段指数模型的 \(WAIC\) 统计量为 2237.4，相应的模型复杂度估计为 \(p_W = 9.6\)。标记为 elpd_waic 的行是预期对数逐点预测密度 (expected log-pointwise predictive density)，是一种预测准确性的度量，将其乘以 -2 得到 \(WAIC\) 统计量。
17.16 延伸阅读
本章只能介绍一小部分用于生存分析的 R 包，以及展示数据操作、向量和矩阵代数等功能的部分内容。此外，对于 R 软件在数据可视化和图形化表示方面的应用，本章所涉非常有限。然而，有很多书籍都介绍了 R 软件及其在统计分析中的应用。其中包括 Crawley (2013), Dalgaard (2008), Davies (2016), Field, Miles and Field (2012), Verzani (2014) 以及 Wickham (2019) 的著作。此外，还有一本名为 R for Dummies的书籍 (de Vries and Meys, 2015)，以及同一系列中介绍使用 R 进行统计分析的书籍 (Schmuller, 2017).
Moore (2016) 描述了生存分析的原则，并说明了使用 R 进行数据分析。Broström (2022) 介绍了生存分析和事件史分析，并提供了使用 R 包 eha 分析社会科学领域数据的实例。
所有 R 包都在 CRAN 网站上进行了描述，但其中一些包还提供了更详细的文档，称为 vignette. 这些 vignette 包括 survival 包的文档 (Therneau, 2022)、使用时依变量进行建模的文档 (Therneau, Crowson and Atkinson, 2022a)、parfm 包的文档 (Munda, Rotolo and Legrand, 2017)以及多状态模型和竞争风险模型的文档 (Therneau, Crowson and Atkinson, 2022b). 还可参考 Munda, Rotolo and Legrand (2012) 关于 parfm 的早期论文。Rizopoulos (2022) 描述了 JM 包，该包用于纵向和事件发生时间数据的联合建模。Brilleman et al. (2020) 描述了使用 rstanarm 包进行贝叶斯生存分析的方法。