由于这是一个很长的问题,我将其分为两个部分;第一个部分是基本问题,第二个部分提供了我迄今为止尝试过的细节。
问题 - 简短
如何在R中拟合个体脆弱性生存模型?特别是我正在尝试重新创建下表中从对该数据集链接进行半参数脆弱性模型拟合得到的系数估计和SE。该模型的形式为:
h_i(t) = z_i h_0(t) exp(\beta'X_i)
其中,z_i
是每个患者的未知脆弱性参数,X_i
是解释变量向量,\beta
是相应系数向量,h_0(t)
是基准风险函数,使用解释变量disease、gender、bmi和age(我在下面的代码中包含了清理因子参考水平的代码)。
问题 - 长篇
我试图按照医学研究中生存数据建模教材的例子来拟合脆弱性模型。特别是我关注半参数模型,该教材提供了正常 cox 模型、对数正态脆弱性和伽马脆弱性的参数和方差估计,如上表所示。
我能够重现没有脆弱性模型估计的结果。
library(dplyr)
library(survival)
dat <- read.table(
"./Survival of patients registered for a lung transplant.dat",
header = T
) %>%
as_data_frame %>%
mutate( disease = factor(disease, levels = c(3,1,2,4))) %>%
mutate( gender = factor(gender, levels = c(2,1)))
mod_cox <- coxph( Surv(time, status) ~ age + gender + bmi + disease ,data = dat)
mod_cox
然而,我真的很难找到一个能够可靠地重新创建第2列结果的软件包。在网上搜索时,我发现了这个表格,试图总结可用的软件包:
https://istack.dev59.com/2wkKb.webp https://cran.r-project.org/web/packages/frailtyEM/vignettes/frailtyEM_manual.pdf
下面我发布了我目前的发现以及我使用的代码,如果有人能够确定我是否只是错误地指定了函数,则可以帮助他们:
frailtyEM - 似乎对于gamma效果最好,但不提供对数正态模型
frailtyEM::emfrail(
Surv(time, status) ~ age + gender + bmi + disease + cluster(patient),
data = dat ,
distribution = frailtyEM::emfrail_dist(dist = "gamma")
)
survival - 给出伽马警告,从我阅读的所有内容来看,它的脆弱功能被视为已过时,建议使用coxme代替。
coxph(
Surv(time, status) ~ age + gender + bmi + disease + frailty.gamma(patient),
data = dat
)
coxph(
Surv(time, status) ~ age + gender + bmi + disease + frailty.gaussian(patient),
data = dat
)
coxme - 看起来可以工作,但提供的估计值与表格中的不同,并且不支持伽马分布。
coxme::coxme(
Surv(time, status) ~ age + gender + bmi + disease + (1|patient),
data = dat
)
frailtySurv - 我无法使其正常运行,似乎总是将方差参数以1的平稳值来适配,并且提供的系数估计好像是拟合了没有脆弱性模型。此外,文档未说明哪些字符串支持脆弱性参数,因此我无法确定如何使其适应对数正态分布。
frailtySurv::fitfrail(
Surv(time, status) ~ age + gender + bmi + disease + cluster(patient),
dat = dat,
frailty = "gamma"
)
frailtyHL - 会产生警告消息,显示“未收敛”,但仍会生成系数估计值,但它们与教科书上的不同。
mod_n <- frailtyHL::frailtyHL(
Surv(time, status) ~ age + gender + bmi + disease + (1|patient),
data = dat,
RandDist = "Normal"
)
mod_g <- frailtyHL::frailtyHL(
Surv(time, status) ~ age + gender + bmi + disease + (1|patient),
data = dat,
RandDist = "Gamma"
)
frailtypack - 我对其实现方式并不理解(或者说与教科书中的有很大不同)。该函数需要指定节点和平滑器,这似乎会极大地影响结果估计。
parfm - 仅适用于参数模型;话虽如此,但每次我尝试使用它来拟合Weibull比例风险模型时,它都会出错。
phmm - 尚未尝试
考虑到我已经不成功地尝试了大量软件包,很可能问题在于我没有正确理解实现方式并且误用了这些软件包。任何有关如何成功重新创建上述估计的帮助或示例都将不胜感激。
?survreg
中的生存率,其中提供了一个“分布”。 - Parfaitsurvreg
主要用于拟合参数模型,在生存分析包中半参数模型的等效物是coxph
,但杂项成分似乎被列为废弃,推荐使用coxme
包,但我无法使其正常工作。 - gowerc