基于R语言实现协变量校正的生存曲线绘制
Editor's Note
推荐关注
The following article is from 劝人学医TDLP Author Roger 不言
如何在R中实现控制协变量的生存曲线可视化？
?
小伙伴
我之前单考虑我分组因素做出来没意义，后来用SPSS纳入多个因素后结果还可以
小编
生存曲线的协变量校正应该是ggadjustedcurves这个函数，但是只能输出生存曲线，也没有P值。。
在经过了一番折腾后，小编建议他在jamovi里面做协变量校正的生存曲线，尽管小编在平常的统计分析中很少用jamovi，但是清楚记得可以在它里面做协变量校正的生存曲线。
临床数据统计|小看你了，jamovi！！！
另外，凑巧下午看到了精鼎统计的一篇关于SPSS中如何实现协变量校正的生存曲线分析（SCI编辑让做控制协变量的生存曲线）。
小编也想还是不要不了了之，把这个话题好好探讨一下。因此，在查阅了一些资料后，我对R中实现协变量校正的生存曲线分析将进行总结。
首先，目前主要有两种方法实现R中协变量控制的生存曲线绘制，第一种是利用adjustedCurves包，而第二种则利用survminer包中的ggadjustedcurves函数。下面我们分别进行介绍。
依赖包及数据加载
library(survival)library(survminer)library(adjustedCurves)data(cancer, package="survival")数据预处理
我们使用survival包中的lung数据。首先我们对原始数据中的结局变量和一部分协变量做简单处理。
data(cancer, package="survival")str(lung)# 'data.frame': 228 obs. of 10 variables:# $ inst : num 3 3 3 5 1 12 7 11 1 7 ...# $ time : num 306 455 1010 210 883 ...# $ status : num 2 2 1 2 2 1 2 2 2 2 ...# $ age : num 74 68 56 57 60 74 68 71 53 61 ...# $ sex : num 1 1 1 1 1 1 2 2 1 1 ...# $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...# $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...# $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...# $ meal.cal : num 1175 1225 NA 1150 NA ...# $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...na.omit(lung)#去除缺失值lung[,c(5,6)]<-lapply(lung[,c(5,6)],factor)#协变量因子化# 结局变量为0,1lung$status <- lung$status - 1str(lung)# 'data.frame': 228 obs. of 10 variables:# $ inst : num 3 3 3 5 1 12 7 11 1 7 ...# $ time : num 306 455 1010 210 883 ...# $ status : num 1 1 0 1 1 0 1 1 1 1 ...# $ age : num 74 68 56 57 60 74 68 71 53 61 ...# $ sex : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 2 2 1 1 ...# $ ph.ecog : Factor w/ 4 levels "0","1","2","3": 2 1 1 2 1 2 3 3 2 3 ...# $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...# $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...# $ meal.cal : num 1175 1225 NA 1150 NA ...# $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...不控制协变量的生存曲线
关于协变量控制的意义，大家可以参考精鼎的这篇推文，这里不再过多赘述。首先，我们看看控制协变量前后生存曲线有什么差异。这里，我们首先不控制协变量，仅仅绘制ECOG评分的KM曲线。
fit1<-survfit(Surv(time,status)~ph.ecog,data = lung)ggsurvplot(fit1,data = lung, conf.int=F, palette = "aaas")未控制协变量ECOG的KM曲线
控制协变量的生存曲线
首先，我们使用adjustedCurves包来实现这一过程。adjustedCurves包由Robin Denz开发，可通过R镜像或者gith安装。
install.packages("adjustedCurves")该包主要使用adjustedsurv函数来实现协变量校正的生存曲线分析，用法及参数如下：
adjustedsurv(data, variable, ev_time, event, method, conf_int=FALSE, conf_level=0.95, times=NULL, bootstrap=FALSE, n_boot=500, n_cores=1, na.action=options()$na.action, clean_data=TRUE, ...)adjustedsurv函数通过method参数基于多种方法实现生存曲线的协变量校正，这里我们使用direct方法。
Currently the following methods can be used:"direct": Direct Standardization based on a previously fit model (Cox-Regression, ...)."direct_pseudo": Direct Standardization based on Pseudo-Values."iptw_km": A weighted Kaplan-Meier estimator."iptw_cox": A weighted estimator based on a stratified weighted Cox-Regression model."iptw_pseudo": A weighted estimator based on Pseudo-Values."matching": Using Propensity Score Matching to estimate the adjusted survival curves."emp_lik": An Empirical Likelihood based estimator."aiptw": An Augmented Inverse Probability of Treatment Weighting estimator."aiptw_pseudo": An Augmented Inverse Probability of Treatment Weighting estimator using Pseudo-Values."strat_amato": A method based on stratification and weighting by Amato (1988)."strat_nieto": A method based on stratification and weighting by Gregory (1988) and Nieto & Coresh (1996)."strat_cupples": A method based on stratification and weighting by Cupples et al. (1995)."km": A simple stratified Kaplan-Meier estimator without any form of adjustment.首先，我们拟合一个协变量校正（多因素）的Cox回归模型。我们在模型中纳入ECOG评分、年龄和性别，通过adjustedsurv函数看校正了病人年龄和性别后，不同ECOG评分患者之间的生存差异。
fit2 <- coxph(Surv(time, status) ~ph.ecog + age + sex, data=lung)协变量校正生存曲线实现及可视化
# calculate and plot curvesadj <- adjustedsurv(data=lung, variable="ph.ecog", ev_time="time", event="status", method="direct", outcome_model=fit2)plot(adj)adjustedsurv函数控制协变量后ECOG的KM曲线
可以看到，协变量控制前后不同ECOG评分组的生存曲线还是有点儿差别的。
下面我们使用survminer包里面自带的函数ggadjustedcurves实现上述过程，首先让我们了解一下ggadjustedcurves的使用方法。其中fit为拟合的cox等比例风险假定模型，variable为我们指定的分组变量。如果在函数中不指定variable，那么则会以fit中的strata（）参数作为分组变量。
ggadjustedcurves( fit, variable = NULL, data = NULL, reference = NULL, method = "conditional", fun = NULL, palette = "hue", ylab = "Survival rate", size = 1, ggtheme = theme_survminer(), ...)当然，这里面最重要的还是method变量，即通过何种方法估计生存曲线。这里主要有四种方法，即single（单一法）、average（均值替代法）、marginal（边际法）和conditional（条件法）。
对于 method = "single"，计算并绘制一条生存曲线。该曲线仅显示总体人群的生存曲线。 对于 method = "average"，为vaviable变量中的每个层级单独绘制生存曲线。如果未指定此参数，则将从 fit 参数的 strata（）中提取。每条曲线基于 Cox 模型拟合的数据估计对应亚组的生存时间。值得注意的是，在此方法中，亚组人群并不均衡。 对于 method = "marginal"，为变量参数选择的分组变量的每个级别绘制生存曲线。如果未指定此参数，则将从 fit 对象的strata（）中提取。 值得注意的是，该方法中亚组人群得到了很好的均衡，以保持与参考群体中的分布相似。如果未指定参考人群，则将整个数据集用作参考人群。 对于 method = "conditional"，为变量参数选择的分组变量的每个级别绘制单独的生存曲线。如果未指定此参数，则将从 fit 对象的strata（）中提取。该方法中亚组人群得到了很好的均衡，只不过方法与第三种有所不同。
# 协变量均值替代法(average in groups)adj1<-ggadjustedcurves(fit2, data = lung, onf.int=T, variable = "ph.ecog", method = "average", palette = "aaas")
# 边际方法(marginal balancing in groups)adj2<-ggadjustedcurves(fit2, data =lung, variable = "ph.ecog", method = "marginal", reference =lung, palette = "aaas")
# 条件方法(conditional balancing in groups)adj3<-ggadjustedcurves(fit2, data = lung, variable = "ph.ecog", method = "conditional", palette = "aaas")
ggadjustedcurves函数控制协变量后ECOG的KM曲线-均值替代法
ggadjustedcurves函数控制协变量后ECOG的KM曲线-边际法
ggadjustedcurves函数控制协变量后ECOG的KM曲线-条件法
控制协变量前后生存曲线确实略有差异；
对于不同的协变量校正方法结果大同小异。
希望这次推文解决了小伙伴的困扰！
文章来源于劝人学医TDLP
关注下方公众号，分享更多更好玩的R语言知识。
点个在看，SCI马上发表。