R语言医学数据分析实战（第九章文中源码解析）

Chapter 9

9.1 生存对象

library(survival)
data(ovarian)
str(ovarian)
# 从书中可以看出resid.ds、rx、ecog.ps这几个变量都是
# 名义型或者有序型变量，常转化因子进行分析，书中（16页）
ovarian$resid.ds <- factor(ovarian$resid.ds,
                           levels = c(1, 2),
                           labels = c("no", "yes"))
ovarian$rx <- factor(ovarian$rx,
                     levels = c(1, 2),
                     labels = c("A", "B"))
ovarian$ecog.ps <- factor(ovarian$ecog.ps,
                          levels = c("1", "2"),
                          labels = c("good", "bad"))
#  查看连续型变量年龄的分布，从直方图中可以看出非对称分布。
# 为了便于分析，将age转化为因子，cut函数将age以50为分界线进行划分
hist(ovarian$age) 
ovarian$agegr <- cut(ovarian$age, 
                     breaks = c(0, 50, 75),
                     labels = c("<=50", ">50"))
# 查看因子（age）的分布情况
table(ovarian$agegr)
# Surv函数在survical库中，创建一个生存对象。
# time是直到事件发生的跟踪时间，event指示预期事件的发生的状态。
surv.obj <- Surv(time = ovarian$futime, event = ovarian$fustat)
surv.obj

9.2 生存率的估计与生存曲线

# 函数对Surv函数计算后的数据进行生存率估计。~后面的1表示没有自变量
# summary函数中的censored = TRUE表示将删失的数据也显示出来
survfit(surv.obj ~ 1)
surv.all <- survfit(surv.obj ~ 1)
summary(surv.all)
summary(surv.all, censored = TRUE)
# plot是绘图函数，mark.time = TRUE这里是指对删失的数据进行“+”标记，conf.int = FALSE表示不显示置信区间
plot(surv.all, mark.time = TRUE,conf.int = FALSE)
# abline常和plot函数一起使用，具体参考https://www.cnblogs.com/xudongliang/p/6755727.html
abline(h = .5, lty = 2, col = "red")

9.3 生存率的比较（这里加入了自变量）

# survfit中对生存对象加入了治疗方法rx进行对比
surv.treat <- survfit(surv.obj ~ rx, data = ovarian)
summary(surv.treat)
# 将两种质量方法在同一幅图显示出来，lty表示使用不同线型，col表示不同颜色
plot(surv.treat, mark.time = T, conf.int = TRUE, 
     lty = c(1, 2), col = c("blue", "red"))
# 添加图例，分表代表x,y,图例，线型，颜色。参考：https://blog.csdn.net/glodon_mr_chen/article/details/79496403
legend(60, .3, legend = c("A", "B"), 
       lty = c(1, 2), col = c("blue", "red"))

# 使用survminer库中ggsurvplot函数进行生存率曲线绘制
install.packages("survminer")
library(survminer)
ggsurvplot(surv.treat, data = ovarian, pval = TRUE)
# survdiff进行统计学检验，最常用的是时序检验
survdiff(surv.obj ~ rx, data = ovarian)

9.4 Cox回归

9.4.1 建立Cox回归模型

# 加入rx、resid.ds 、 agegr 、 ecog.ps将所有协变量添加进去
# coxph函数用于建立Cox回归模型
cox1 <- coxph(surv.obj ~ rx + resid.ds + agegr + ecog.ps, 
              data = ovarian)
summary(cox1)
# 剔除效果不显著的变量，AIC值较小的
drop1(cox1, test = "Chisq")
cox2 <- coxph(surv.obj ~ rx + resid.ds + agegr, data = ovarian)
# 使用step函数自动完成剔除过程
step.cox <- step(cox1)

9.4.2 比例风险假定的检验

# 使用cox.zph函数执行PH检验、检验模型是否满足等比例风险
cox.zph(cox2)
summary(cox2)

9.4.3 生存的预测

# data.frame建立新的数据框
newdata <- data.frame(rx = c("A" , "B"), 
                      resid.ds = c("no", "no"),
                      agegr = c(">50", ">50"))
newdata

# 对新数据使用predict进行预测
hr <- predict(cox2, newdata = newdata, type = "risk")
hr
hr[1]/hr[2]

cox.fit <- survfit(cox2, newdata = newdata, type = "kaplan-meier")
plot(cox.fit, lty = c(1, 2), col = c(2, 4))
title(main="Cox survival curves by treatment
      for age > 50, no residual disease patients",
      xlab="Duration in days", 
      ylab="Survival probability",
      las = 1)
legend(5, 0.3, c("Treatment A", "Treatment B"), 
       lty = c(1, 2), col = c(2, 4))