如何在R中拟合虚弱生存模型

3

由于这是一个很长的问题,我将其分为两个部分;第一个部分是基本问题,第二个部分提供了我迄今为止尝试过的细节。

问题 - 简短

如何在R中拟合个体脆弱性生存模型?特别是我正在尝试重新创建下表中从对该数据集链接进行半参数脆弱性模型拟合得到的系数估计和SE。该模型的形式为:

h_i(t) = z_i h_0(t) exp(\beta'X_i)

其中,z_i 是每个患者的未知脆弱性参数,X_i 是解释变量向量,\beta 是相应系数向量,h_0(t) 是基准风险函数,使用解释变量diseasegenderbmiage(我在下面的代码中包含了清理因子参考水平的代码)。

问题 - 长篇

我试图按照医学研究中生存数据建模教材的例子来拟合脆弱性模型。特别是我关注半参数模型,该教材提供了正常 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 - 尚未尝试

考虑到我已经不成功地尝试了大量软件包,很可能问题在于我没有正确理解实现方式并且误用了这些软件包。任何有关如何成功重新创建上述估计的帮助或示例都将不胜感激。


如果你的问题只是关于如何使用R、哪个包、什么代码等,那么它在这里是不适合讨论的。这是一个大问题,在[SO]上可能得不到太多回应。你最好去r-help寻求帮助。如果你想要我把这个问题迁移到[SO],请告诉我。 - gung - Reinstate Monica
嗨 @gung ,事实上,我认为在SO上得到回复的可能性很小,所以我在这里发布了帖子,因为我假设这里更多统计方面的用户会更有可能知道答案或能够提供帮助。如果你仍然坚持认为这是离题的话,请将其转移到SO;如果我在那里没有得到任何回复,我会尝试按照你建议的联系r-help。 - gowerc
你是否有“肺移植登记患者的生存率”数据可用?请查看 ?survreg 中的生存率,其中提供了一个“分布”。 - Parfait
嗨 @Parfait,我已更新问题并附上了数据集的直接下载链接以尝试澄清问题。感谢您的建议;survreg主要用于拟合参数模型,在生存分析包中半参数模型的等效物是 coxph,但杂项成分似乎被列为废弃,推荐使用 coxme 包,但我无法使其正常工作。 - gowerc
1个回答

1

关于

我真的很难找到一个可以可靠地重新创建第二列和第三列结果的软件包。

请参阅 Survival Analysis CRAN task view 中的 Random Effect Models,或在 R Site Search 上搜索例如 "survival frailty"。


感谢您的建议;根据这个资源,frailtypack 是最高评分的生存衰弱建模软件包。我曾尝试使用这个软件包,但很难理解它或使其正常工作,因为其实现似乎与教科书中的非常不同。特别是,您似乎必须指定许多结点和平滑值,这些值似乎极大地影响结果(我不知道如何使用这些值来重新创建上面的数字)。 - gowerc
这里有两篇来自《统计软件杂志》的论文[1][2],可能会对您有所帮助。 [1] https://www.jstatsoft.org/article/view/v047i04 [2] https://www.jstatsoft.org/article/view/v081i03 - Benjamin Christoffersen
谢谢,我会看一下并尝试使用这些来解决我的问题。 - gowerc
第一篇论文可能是与您最相关的。尽管我记得打包是关于_相关生存数据_,而您的问题只涉及个体脆弱性,如果我正确理解了您的问题?我不确定您是否可以使用 frailtypack 中的任何函数来适应这样的模型。 - Benjamin Christoffersen
嵌套脆弱性模型是您想要的,但第3页上没有共享$v_i$效应。 - Benjamin Christoffersen

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接