基于R语言实现协变量校正的生存曲线绘制

近期有一位小伙伴遇到了一个问题：如何在R中实现控制协变量的生存曲线可视化。

如何在R中实现控制协变量的生存曲线可视化？

?

小伙伴

我之前单考虑我分组因素做出来没意义，后来用SPSS纳入多个因素后结果还可以

小编

生存曲线的协变量校正应该是ggadjustedcurves这个函数，但是只能输出生存曲线，也没有P值。。

在经过了一番折腾后，小编建议他在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,1
lung$status <- lung$status - 1
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  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 curves
adj <- 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（）中提取。该方法中亚组人群得到了很好的均衡，只不过方法与第三种有所不同。

现在，让我们看看第2-4方法对应的协变量校正生存曲线有何不同吧。

# 协变量均值替代法(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马上发表。