R数据分析：cox模型如何做预测，高分文章复现

今天要给大家分享的文章是

Cone EB, Marchese M, Paciotti M, Nguyen DD, Nabi J, Cole AP, Molina G, Molina RL, Minami CA, Mucci LA, Kibel AS, Trinh QD. Assessment of Time-to-Treatment Initiation and Survival in a Cohort of Patients With Common Cancers. JAMA Netw Open. 2020 Dec 1;3(12):e2030072. doi: 10.1001/jamanetworkopen.2020.30072. PMID: 33315115; PMCID: PMC7737088.

作者做了癌症结局与延迟治疗时间的关系。作者希望能给在疫情背景下怎么样更好地分流癌症患者这一实际问题提供实证依据。作者纳入了4个癌症，发现了基本上Time-to-treatment initiation (TTI)约长，癌症患者的5年和10年死亡率越高。这个结果和目前的部分指南其实是矛盾的，所以整篇文章还是有一定实际意义的。

本文依然是关注文章如何在做法上在统计上进行复现，启发大家如何用自己的数据做一个同样设计的研究。

文章回顾

作者将关心的主变量TTI进行了分类处理，将延迟治疗时间分了4类：

TTI was a categorical variable with 4 levels: 8 to 60, 61 to 120, 121 to 180, and greater than 180 days

为了准确地回答TTI如何影响了全因死亡率，作者还纳入了一系列的协变量包括：clinical demographic variables，clinical demographic variables，Modality of any treatment。构建了一个cox模型。

作者呈现出来的结果包括两个表，如下图，分别是纳入的不同癌症人群的人群特征和按照延迟治疗时间分组后的癌症特征：

 

这两个表都是属于非常简单的描述统计，大家用用tableone或者gtsummary都可以很方便的做出来（请参考之前文章）。

具体到主分析，作者以all-cause mortality为结局变量，用训练好的cox模型以图的形式展示了不同TTI组随着时间变化predicted 5-year and 10-year all-cause mortality的变化情况：

 

上图就是文章中截取的，展示的是通过cox模型预测出来的不同TTI类别组的随着时间变化的生存概率。整篇文章的主要统计工作就是这些，那么我们今天要做的就是复现出来整个预测概率并出图。

实例操练

我手上有数据如下：

 

我希望训练一个cox模型然后像文章一样通过模型预测情况看一看，不同用药组随着时间变化其生存概率的变化是如何的。

首先我用我的数据训练一个简单的cox模型，需要用到coxph函数，整体代码如下：

fit <- coxph(Surv(Time, death) ~ CD4 + drug + gender, data = data)

在上面的代码中我只是随意选择了3个变量：CD4 + drug + gender，或者说我在CD4 + gender存在的情况下我关心不同的drug（对应文章中的TTI）对患者随着时间变化的生存概率的影响情况。给大家展示下模型结果

 

可以看到drug的效应并不显著，但是不影响我们示例作图。

到目前我们的cox模型训练好了，接下来我们需要用这个模型预测各组（此时是Drug变量的不同水平）从而得到模型的预测值，这个预测值一定是用控制了一系列的协变量之后的cox模型得到的预测值，作者的处理办法是将协变量都固定在均值水平：

Deriving predicted outcomes, while setting the covariate values at their mean, allowed us to estimate the mean difference in outcomes for each patient under different TTI while conserving their clinical and demographic characteristics

对于我们的cox模型来讲，我们关心的是drug这个变量，其余的变量同样的均需要固定在均值水平，生成新数据的代码如下：

ND1 <- with(aids.id, expand.grid(CD4 = 3, drug = levels(drug), gender = factor("male", levels = levels(gender))))

在上面的代码中，我们将协变量1--CD4固定在了3，将协变量性别固定为了男性，这个时候我们将数据带入训练好的cox模型中可以得到预测值surv_probs_Cox，直接用plot进行作图：

plot(surv_probs_Cox, col = c("black", "red"), xlab = "Follow-Up Time (months)", ylab = "Survival Probabilities", bty='l', lwd=2 ) grid(ny=5,nx=NA,col='gray',lty = 1)

运行代码后得到下图，下面的图就展示了随着时间变化不同用药的病人死亡概率的差异。

 

当然啦，为了和原文更像一点，我们可以用同样的数据在ggplot2中重新画下，并且像原文章一样将曲线变为smoothed curve，效果如下：

 

基本上和原文的展示形式一样了。每种癌症做完，形成大图，整篇文章的主分析也就完成了，好了，本文毕。

小结

今天给大家写了一篇文章的结果呈现方法，希望传递的信息是我们做cox模型，不一定都是在正文中给出一系列的HR，当我们能收到的临床数据很少的时候我们用图来展示结果有可能更加的easily interpreted by clinicians and policymakers。