R语言进阶之生存分析
在医学研究中，生存分析是一类非常重要的统计方法，它主要的目的是对生存率和时间进行建模，计算患者在特定时间段内生存的概率，主要用于评估治疗的效果和疾病的危险程度。由于患者可能在研究结束时或存活或死亡，还有一部分可能失联（可能活着），导致出现一定比例的删失值。因此，生存分析的数据分析也有其自身的特点。
R语言中进行生存分析主要使用“survival”这个包，一看名字就感觉这个包专业而靠谱，实际上确实如此。“survial”包可以针对单样本或者多样本进行生存分析，可以使用的模型有参数加速失效模型（parametric accelerated failure models）和Cox比例风险模型（Cox proportional hazards model）。
生存分析的数据一般包括三个部分：起始时间、结束时间和结束时患者的状态（通常1代表目标事件发生，0代表删失数据，即目标事件未发生或者失访等），这里的目标事件一般是指疾病或者死亡。另外，数据也可以是到结束时所经历的时间段和结束时患者的状态。通常，我们使用Surv()函数来将数据进行格式转化，便于进行后续分析。
在生存分析中，我们最常用的是如下三个函数：
survfit( ) # 主要用于计算单个或多个组的生存分布
survdiff( ) # 主要用于检验不同组的生存分布差异
coxph( ) # 主要用于拟合Cox比例风险模型这里我们使用Mayo诊所的临床数据（lung）进行生存分析，这个数据集收集了228肺癌患者的信息，包括年龄、性别、生存天数、能量摄入和Karnofsky得分等，整个分析共分为5个步骤。
第一步 加载R包并导入数据
# 加载R包
library(survival)
# 查询数据集的相关信息
help(lung)
# 创建Surv对象
survobj <- with(lung, Surv(time,status)) # 这里主要是指定时间和生存状态第二步 评估整体生存率
# 绘制总样本的生存分布曲线
# 使用Kaplan-Meier估计
fit0 <- survfit(survobj~1, data=lung) # 构建总体生存分析模型
summary(fit0) # 查看结果
plot(fit0, xlab="Survival Time in Days",
ylab="% Surviving", yscale=100,
main="Survival Distribution (Overall)") # 绘图
上图中的横坐标代表天数，纵坐标代表生存率；实线是估计的生存率，虚线代表的是估计生存率的95%置信区间。
第三步 比较不同性别之间生存率的差异
# 比较男性和女性的生存差异
fit1 <- survfit(survobj~sex,data=lung) # 在公式里指定“sex”就可以比较性别差异了
# 依据性别绘制生存分布曲线
plot(fit1, xlab="Survival Time in Days",
ylab="% Surviving", yscale=100,col=c("red","blue"),
main="Survival Distributions by Gender")
legend("topright", title="Gender",c("Male", "Female"),
fill=c("red", "blue")) # 绘图
从上图我们可以看出，女性的整体生存率是高于男性的，但是这种差异是否显著还并不确定，这就需要我们进行下面的统计检验。
# 对性别的生存差异进行统计检验
survdiff(survobj~sex, data=lung) # 在公式里指定“sex“
结果的P值小于0.05，我们可以认为男女之间的生存率是有差异的。
第四步 基于Cox比例风险模型的生存分析
# 基于Cox比例风险模型从年龄和医学评分来预测男性的生存情况
MaleMod <- coxph(survobj~age+ph.ecog+ph.karno+pat.karno, data=lung, subset=sex==1)
MaleMod # 输出结果
从这上面的结果可以看出，ph.ecog和ph.karno对生存率的影响是显著的，而其他变量则并不显著。
第五步 评估Cox比例风险模型的假设检验条件
# 评估cox比例风险的假设条件
cox.zph(MaleMod)
除ph.karno外，其余P值都很大，因此ph.karno可能会违背Cox比例风险模型的假设，需要谨慎对待。
关于绘图的相关内容请参见R语言入门系列教程。
生存分析的内容就先分享到这里，欢迎大家持续关注【生信与临床】。