用survminer优雅的画生存曲线

生存曲线绘制

#Survival Curves
#The ggsurvplot() function creates ggplot2 plots from survfit objects.
data(lung)
head(lung)
  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

table(lung$sex)
  1   2 
138  90 

fit <- survfit(Surv(time,status)~ sex, data = lung)

class(fit)
[1] "survfit"

library(survminer)
ggsurvplot(fit,data = lung)

#Use the fun argument to set the transformation of the survival curve
#"event" for cumulative events, "cumhaz" for the cumulative hazard function or "pct" for survival probability in percentage.
ggsurvplot(fit, data = lung, fun = "event")
ggsurvplot(fit, data = lung, fun = "cumhaz")
ggsurvplot(fit, data = lung, fun = "pct")#默认pct

#很多的参数可以调，position and content of the legend; additional annotations like p-value, title, subtitle.
ggsurvplot(fit, data = lung,
           conf.int = TRUE,#添加置信区间
           pval = TRUE,#添加P值
           fun = "pct",
           risk.table = TRUE,
           size = 1,
           linetype = "strata",
           palette = c("#E7B800","#2E9FDF"),
           legend = "bottom",
           legend.title = "Sex",
           legend.labs = c("Male","Female"))

Diagnostics of Cox Model

Cox模型的诊断一般包括三方面的内容：

• 比例风险假定；
• 模型影响点（异常值）识别；

• 比例风险的对数值与协变量之间的非线性关系识别；
  对上述三方面的诊断，常见的方法为残差法。

• Schoenfeld残差用于检验比例风险假定；

• Deviance残差用于影响点（异常值）识别；
• Martingale残差用于非线性检验；

The function cox.zph() from survival package may be used to test the proportional hazards assumption for a Cox regression model fit. The graphical verification of this assumption may be performed with the function ggcoxzph() from the survminer package. For each covariate it produces plots with scaled Schoenfeld residuals against the time.

模型诊断——PH假定。PH假定可通过假设检验和残差图检验。正常情况下，Schoenfeld残差应该与时间无关，如果残差与时间有相关趋势，则违反PH假设的证据。残差图上，横轴代表时间，如果残差均匀的分布则，表示残差与时间相互独立。

library(survival)
fit <- coxph(Surv(time, status) ~ sex + age, data = lung)
ftest <- cox.zph(fit)
ftest
           rho chisq     p
sex     0.1236 2.452 0.117
age    -0.0275 0.129 0.719
GLOBAL      NA 2.651 0.266
#绘图
library(survminer)
ggcoxzph(ftest)

从上面的结果可以看出，三个变量的P值都大于0.05，说明每个变量均满足PH检验，而模型的整体检验P值0.266，模型整体满足PH检验。上图中实线是拟合的样条平滑曲线，虚线表示拟合曲线上下2个单位的标准差。如果曲线偏离2个单位的标准差则表示不满足比例风险假定。从上图中可见，各协变量满足PH风险假设。

模型诊断——模型影响点（异常值）识别

我们可以通过绘制Deviance残差图或者dfbeta值实现上述诊断。在R语言survminer中ggcoxdiagnostics()函数可以画出Deviance残差图。
The function ggcoxdiagnostics() plots different types of residuals as a function of time, linear predictor or observation id. The type of residual is selected with type argument. Possible values are “martingale”, “deviance”, “score”, “schoenfeld”, “dfbeta”’, “dfbetas”, and “scaledsch”. The ox.scale argument defines what shall be plotted on the OX axis. Possible values are “linear.predictions”, “observation.id”, “time”. Logical arguments hline and sline may be used to add horizontal line or
smooth line to the plot.

library("survival")
library("survminer")
fit <- coxph(Surv(time, status) ~
               sex + age, data = lung)
ggcoxdiagnostics(fit,
                 type = "deviance",
                 ox.scale = "linear.predictions")
ggcoxdiagnostics(fit,
                 type = "schoenfeld",
                 ox.scale = "time")

残差值均匀的分布在0上下，表明满足上述假定。

Summary of Cox Model

The function ggforest() from the survminer package creates a forest plot for a Cox regression model fit. Hazard ratio estimates along with confidence intervals and p-values are plotter for each variable.

> library("survival")
> library("survminer")
> lung$age <- ifelse(lung$age > 70, ">70","<= 70")
> fit <- coxph( Surv(time, status) ~ sex + ph.ecog + age, data = lung)
> fit
Call:
coxph(formula = Surv(time, status) ~ sex + ph.ecog + age, data = lung)

          coef exp(coef) se(coef)     z       p
sex     -0.567     0.567    0.168 -3.37 0.00075
ph.ecog  0.470     1.600    0.113  4.16 3.1e-05
age>70   0.307     1.359    0.187  1.64 0.10175

Likelihood ratio test=31.59  on 3 df, p=6e-07
n= 227, number of events= 164 
   (1 observation deleted due to missingness)

校正生存曲线绘制

The function ggcoxadjustedcurves() from the survminer package plots Adjusted Survival Curves for Cox Proportional Hazards Model. Adjusted Survival Curves show how a selected factor influences survival estimated from a Cox model. Note that these curves differ from Kaplan Meier estimates since they present expected ssurvival based on given Cox model.

library("survival")
library("survminer")
lung$sex <- ifelse(lung$sex == 1,"Male", "Female")
fit <- coxph(Surv(time, status) ~sex + ph.ecog + age,data = lung)
ggcoxadjustedcurves(fit, data=lung,variable=lung$sex)

Note that it is not necessary to include the grouping factor in the Cox model. Survival curves are estimated from Cox model for each group defined by the factor independently.

lung$age3 <- cut(lung$age, c(35,55,65,85))
ggcoxadjustedcurves(fit, data=lung, variable=lung$age3)