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

How to fit frailty survival models in R

提问人:gowerc 提问时间:3/10/2018 最后编辑:gowerc 更新时间:10/8/2019 访问量:4284

问:

因为这是一个很长的问题,所以我把它分成了两部分;第一个只是基本问题,第二个提供了我到目前为止尝试的细节。

问题 - 短

如何在 R 中拟合个体脆弱生存模型?特别是,我正在尝试重新创建下表中的系数估计值和 SE,这些估计值是通过将半参数脆弱模型拟合到此数据集链接中找到的。该模型采用以下形式:

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

其中,是每个患者的未知衰弱参数,是解释变量的向量,是系数的相应向量,是使用解释变量疾病性别BMI年龄的基线风险函数(我在下面包含代码以清理因素参考水平)。z_iX_i\betah_0(t)enter image description here

问题 - 长

我正在尝试遵循并重新创建医学研究教科书中的建模生存数据示例,以拟合虚弱的 mdoels。我特别关注半参数模型,教科书为该模型提供了正态 cox 模型、对数正态衰弱性和 Gamma 衰弱的参数和方差估计,如上表所示

我能够使用

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 列结果的包。在网上搜索,我找到了这个表格,它试图总结可用的软件包:

enter image description here

下面我发布了我当前的发现以及我使用的代码,它可以帮助某人识别我是否只是错误地指定了函数:

frailtyEM - 似乎最适合伽玛,但不提供对数正态模型

frailtyEM::emfrail(
    Surv(time, status) ~ age + gender + bmi + disease + cluster(patient), 
    data = dat ,
    distribution = frailtyEM::emfrail_dist(dist = "gamma")
)

生存 - 对伽马发出警告,从我读到的所有内容来看,它的脆弱功能似乎被归类为贬值,建议改用 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 - 仅拟合参数化模型;话虽如此,每次我试图用它来拟合威布尔比例风险模型时,它都会出错。

phmm - 还没试过

鉴于我获得的大量软件包没有成功,我完全理解问题很可能是我自己没有正确理解实现并错过了使用这些软件包。不过,有关如何成功重新创建上述估计的任何帮助或示例将不胜感激。

R 生存

评论

0赞 gung - Reinstate Monica 3/10/2018
如果你的问题只是关于如何使用 R、哪个包、什么代码等,那么这里就偏离了主题。这是一个很大的问题,它可能在 Stack Overflow 上没有得到太多回应。你可能在 r-help 上做得最好。如果您希望我将其迁移到 SO,请告诉我。
0赞 gowerc 3/10/2018
嗨@gung,事实上,我认为我不太可能得到关于 SO 的回复,这就是我在这里发帖的原因,因为我认为这里更注重统计的用户将有更高的机会知道答案或能够提供帮助。如果您仍然坚持认为这是题外话,那么是的,请您将其转移到 SO 吗?如果我在那里没有得到任何回复,我会按照您的建议尝试联系 r-help。
0赞 Parfait 3/10/2018
您是否有肺移植登记患者的生存数据?查看您提供分配的生存。?survreg
0赞 gowerc 3/10/2018
嗨,@Parfait,我已经更新了问题以尝试清除它,并附上了下载数据集的直接链接。谢谢你的建议; 虽然专注于拟合参数模型,但半参数模型生存包中的等价物是,脆弱的组件似乎被归类为折旧,以支持我无法正常工作的包survregcoxphcoxme

答:

1赞 Benjamin Christoffersen 3/10/2018 #1

关于

我真的很难找到一个可以可靠地重新创建后 2 列结果的包。

请参阅“随机效应模型”下的“生存分析 CRAN”任务视图,或在 R 站点搜索中搜索“生存衰弱”。

评论

0赞 gowerc 3/10/2018
谢谢你的建议;从这个资源来看,它似乎是生存衰弱建模的最高分包。我曾尝试使用这个包,但很难理解它或让它正常工作,因为实现似乎与教科书所做的非常不同。特别是,您似乎必须指定一些结点和平滑值,这似乎对结果有很大影响(我不确定如何使用这些值来重新创建上述数字)。frailtypack
0赞 Benjamin Christoffersen 3/10/2018
这篇论文和这篇来自《统计软件杂志》的论文可能会对你有所帮助。
0赞 gowerc 3/10/2018
谢谢;会看看我是否可以用这些来解决我的问题
0赞 Benjamin Christoffersen 3/10/2018
第一篇论文可能与你最相关。不过,据我所知,包装是关于相关的生存数据,而如果我正确理解您的问题,您的问题只涉及个人虚弱?我不确定您是否可以将这样的模型与中的任何功能相匹配。frailtypack
0赞 Benjamin Christoffersen 3/10/2018
例如,嵌套的脆弱模型是你想要的,但没有第 3 页上共享的 $v_i$ 效应。