1. 数据准备
install.packages(c("survival","survminer"))
> 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
数据字段说明
inst：机构代码
time：以天为单位的生存时间
status：删失状态1 = 删失，2 = 出现失效事件
age：岁
sex：性别，男= 1女= 2
ph.ecog：ECOG评分（0 =好，5 =死）
ph.karno：医师进行的Karnofsky评分（0 = 差，100 = 好）
pat.karno：患者自行进行的Karnofsky评分（0 = 差，100 = 好）
meal.cal：用餐时消耗的卡路里
wt.loss：最近六个月的体重减轻
2. 拟合生存曲线
lung$sex <- factor(lung$sex,
levels = c(1,2),
labels = c("male", "female"))
dd=datadist(lung)
options(datadist="dd")
fit <- survfit(Surv(time,status) ~ sex, data = lung)
print(fit)
Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
n events median 0.95LCL 0.95UCL
sex=male 138 112 270 212 310
sex=female 90 53 426 348 550
默认情况下，函数print()显示生存曲线的简短摘要。它打印观察次数、事件次数、中位数存活率和中位数的置信度。
如果要显示更完整的生存曲线摘要，请键入以下内容
# Summary of survival curves
summary(fit)
# Access to the sort summary table
summary(fit)$table
访问survfit()返回的值
函数survfit()返回变量列表，包括以下组件：
n: 每条曲线上的受试者总数。
time: 曲线上的时间点。
n.risk: 在时间t处于危险状态的受试者数量
n.event: 在时间t发生的事件数。
n.censor: 在时间t无事件情况下退出风险集的被审查对象的数量。
lower,upper: 曲线的下限和上限。
strata: 表示曲线估计的分层。如果层不为空，则结果中有多条曲线。地层级别(因子)是曲线的标签。
可以按如下方式访问组件：
d <- data.frame(time = fit$time, n.risk = fit$n.risk, n.event = fit$n.event, n.censor = fit$n.censor, surv = fit$surv, upper = fit$upper, lower = fit$lower )
head(d)
3. 绘制基础曲线
> library(survival)
> library(survminer)
#可选调色板有 "grey","npg","aaas","lancet","jco", "ucscgb","uchicago","simpsons"和"rickandmorty".
ggsurvplot(fit, # 创建的拟合对象
data = lung, # 指定变量数据来源
conf.int = TRUE, # 显示置信区间
pval = TRUE, # 添加P值
risk.table = TRUE, # 绘制累计风险曲线
surv.median.line = "hv", # 添加中位生存时间线
add.all = TRUE, # 添加总患者生存曲线
palette = "hue") # 自定义调色板
可以使用以下参数进一步自定义绘图：
- conf.int.style = “step” 更改置信区间波段的样式
- xlab 更改x轴标签
- break.time.by = 200 将x轴在时间间隔中中断by200
- risk.table = “abs_pct”显示处于危险状态的个人的绝对数量和百分比
- risk.table.y.text.col = TRUE 和 risk.table.y.text = FALSE 在风险表图例的文本注释中提供条形图而不是名称。
- ncensor.plot = TRUE to 画出在时间t处被审查的对象的数量。正如马尔辛·科辛斯基所建议的那样, 这是对生存曲线的一个很好的补充反馈，这样人们就可以意识到：生存曲线是什么样子的，风险集的数量是多少，风险集变小的原因是什么：是由事件造成的，还是由审查事件造成的？
- legend.labs 若要更改图例标签。
Kaplan-Meier Plot可以解释如下：
横轴(x轴)表示时间(以天为单位)，纵轴(y轴)表示存活的概率或存活人数的比例。这两条线代表了两组人的生存曲线。曲线上的垂直下落表示事件。曲线上的垂直刻度线意味着患者在这个时候被审查了。
在时间零时，存活概率为1.0(或100%的参与者还活着)。
在时间250，性别=1的生存概率约为0.55(或55%)，性别=2的生存概率约为0.75(或75%)。
性别=1的中位生存期约为270天，性别=2的中位生存期约为426天，这表明与性别=1相比，性别=2的中位生存期较好
使用下面的代码可以获得每组的中位生存时间:
summary(fit)$table
比较生存曲线的对数秩检验：Survdiff()
对数秩检验是比较两条或两条以上生存曲线最常用的方法。零假设是两组患者的存活率没有差异。对数秩检验是一种非参数检验，对生存分布没有任何假设。本质上，对数秩检验将每组中观察到的事件数量与如果零假设为真(即，如果生存曲线相同)的预期数量进行比较。对数秩统计量近似分布为卡方检验统计量。
[生存包]中的函数survdiff()可用于计算比较两条或多条生存曲线的对数秩检验。
Survdiff()的用法如下：
surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)
surv_diff
Call:
survdiff(formula = Surv(time, status) ~ sex, data = lung)
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 138 112 91.6 4.55 10.3
sex=2 90 53 73.4 5.68 10.3
Chisq= 10.3 on 1 degrees of freedom, p= 0.00131
该函数返回组件列表，包括:
n: 每组中受试者的数量。
obs: 每组中观察到的加权事件数。
exp: 每组事件的加权预期数量。
chisq: 检验平等的chisquare统计量。平等性检验的平方统计量。
strata: 可选地，每个层中包含的受试者的数量。
生存差异的对数秩检验p值为p=0.0013，表明不同性别群体的生存有显著差异。
还可以拟合复杂生存曲线
4. 构建Cox模型
#用 cph()函数来拟合Cox比例风险回归模
res.cox <- cph(Surv(time, status) ~ age + sex, data = lung)
#构建生存函数
med <- Quantile(res.cox) # 计算中位生存时间
surv <- Survival(res.cox) # 构建生存概率函数
#绘图：预测中位生存时间
nom <- nomogram(res.cox, fun=function(x) med(lp=x),
funlabel="Median Survival Time")
plot(nom)
绘图：预测生存概率
nom <- nomogram(res.cox, fun=list(function(x) surv(365, x),
function(x) surv(730, x)),
funlabel=c("1-year Survival Probability",
"2-year Survival Probability"))
plot(nom, xfrac=.6)
评价COX回归的预测效果
我们主要计算“ C-index ”即C-指数和绘制矫正曲线，来对Cox回归的预测效果进行评价。
1. 计算C-指数
> rcorrcens(Surv(time,status) ~ predict(res.cox), data =lung)
Somers' Rank Correlation for Censored Data Response variable:Surv(time, status)
C Dxy aDxy SD Z P n
predict(res.cox) 0.397 -0.206 0.206 0.051 4.03 1e-04 228
2. 绘制校正曲线
这里对参数做一些说明：
- 绘制校正曲线前需要在模型函数中添加参数x=T, y=T，详细参考帮助
- u需要与之前模型中定义好的time.inc一致，即365或730；
- m要根据样本量来确定，由于标准曲线一般将所有样本分为3组（在图中显示3个点）
- m代表每组的样本量数，因此m*3应该等于或近似等于样本量；
- b代表最大再抽样的样本量
重新调整模型函数res.cox2
这里是加上了x=T, y=T，time.inc = 365三个参数：
res.cox2 <- cph(Surv(time,status) ~ age+sex, data =lung, x=T, y=T, surv=TRUE, time.inc = 365)
#构建校正曲线
cal1 <- calibrate(res.cox2, cmethod='KM', method="boot", u=365, m=76, B=228)
#绘制校正曲线
par(mar=c(8,5,3,2),cex = 1.0)
plot(cal1,lwd=2,lty=1,
errbar.col=c(rgb(0,118,192,maxColorValue=255)),
xlim=c(0.25,0.6),ylim=c(0.15,0.70),
xlab="Nomogram-Predicted Probability of 1-Year DFS",
ylab="Actual 1-Year DFS (proportion)",
col=c(rgb(192,98,83,maxColorValue=255)))
参考文献：
