目录
一、生存分析
基本术语
二、R语言实现
0.加载包
1.示例数据，R包自带数据集"lung"
2.创建生存对象并拟合函数
3.可视化
3.1 默认设置作图
3.2 我的偏好作图
3.3 我的作图函数
三、其他
1.p值如何计算
一、生存分析
一些医学事件所经历的时间，从开始观察到事件发生的时间，不是短期内可以明确判断的。
e.g.乳腺癌病人术后生存时间
针对这类生存资料的分析方法：生存分析。
基本术语
event 事件
例如，患者死亡，患者复发……取决于研究要观察的事件
time 时间
一般指，从开始观察到事件发生的时间。
| xx期 | 定义 |
|---|---|
| OS（Overall survival） | 总体生存期：从初次治疗到患者死亡的时间。 |
| RFS（Recurrence free survival） | 无复发生存期：从初次治疗到出现复发或者随访截止的时间。 |
| PFS（Progression-free survival） | 从初次治疗到疾病进展或因任何原因死亡的时间。一般PFS越长，意味着患者有质量的生存时间越长，活得越好。 |
status 状态
| status | 含义 |
|---|---|
| 1 | 观察对象出现了终点事件，e.g.观察对象在手术后5年死亡 |
| 0 | 观察对象没有出现终点事件，e.g.观察对象在手术后第3年，与该患者失去联络，该患者的生存情况未知；或者该患者仍存活；或者…… |
二、R语言实现
0.加载包（如果没有，就安装）
library(survival)
library(survminer)
1. 示例数据，R包自带数据集"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
| 变量 | 含义 |
|---|---|
| inst | Institution code |
| time | Survival time in days |
| status | censoring status 1=censored, 2=dead |
| age | Age in years |
| sex | Male=1 Female=2 |
| ph.ecog | ECOG performance score as rated by the physician. 0=asymptomatic, 1= symptomatic but completely ambulatorph.karno: Karnofsky performance score (bad=0-good=100) rated by physician |
| pat.karno | Karnofsky performance score as rated by patient |
| meal.cal | Calories consumed at meals |
| wt.loss | Weight loss in last six months (pounds) |
2.创建生存对象并拟合函数
创建生存对象需要三个变量
| 变量 | 注 |
|---|---|
| 时间，即time | 例如生存的时间 |
| 状态 ，即status | 例如这个样本是否出现终点事件（0或1） |
| 分组标签 | 本例中分为男女两组 |
本文按照性别分组，男性是1，女性是2
table(lung$sex)
1 2
138 90
接下来是创建生存对象
fit <- survfit(Surv(time, status) ~ sex, data = lung)
3.可视化
3.1 默认设置作图
使用函数ggsurvplot()按照默认设置作图
ggsurvplot(fit)
Survival analysis
图中曲线下降说明某些观察对象出现了了终点事件；
图中的"+"表示是删失数据（Cencor），即未出现终点事件的观察对象，可能由于失去联络而未知该观察对象的生存情况。
3.2 我的作图偏好
其中HR(95%CI)和C-index是用函数annotate()添加，这部分要先计算cox模型。
注意点：
- R语言运行的HR的含义。男性是1，女性是2。按顺序排序，男性是第一组，女性是第二组，HR=0.59 < 1的含义是第二组（女性）相对于第一组（男性）是一个保护因素。为了提前标记好分组的顺序，我的习惯是提前将分组标记为0，1……。
拟合cox模型，及一些必要的准备。
cox <- coxph(Surv(time, status) ~ sex, data = lung)
res <- summary(cox)
x <- lung$time
作图代码如下：
p<-ggsurvplot(fit,
pval = TRUE, pval.coord = c(0.7*max(x), 1),
pval.method = TRUE, pval.method.coord=c(0.4*max(x),1),
risk.table = TRUE, risk.table.col = "strata",
ggtheme = theme_linedraw(),
palette = c("darkblue", "red"),
legend.title = "Gender",
legend.labs = c("Male", "Female"),
xlab = "Time(days)", xlim = c(0,1000), break.x.by = 200,
ylab = "Survival probability"
)
p$plot <- p$plot+
theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+
ggtitle("My Survival Plot")+
theme(plot.title = element_text(hjust = 0.5))+
annotate("text",x=0.7*max(x),y=0.91, label=paste0("HR(95%CI) = ", round(res$conf.int[1],2), "(",round(res$conf.int[3],2),"-", round(res$conf.int[4],2), ")" ))+
annotate("text",x=0.7*max(x),y=0.83, label=paste0("C-index = ", round(res$concordance[1],2)))
p$table <- p$table+
theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())
png("My Survival analysis.png", width = 4290, height = 4930, res = 1000)
p
dev.off()
| 参数 | 含义 |
|---|---|
| pval | 是否显示p-value |
| pval.coord | p-value的位置 |
| pval.method | 是否显示计算p-value的方法 |
| pval.method.coord | pval.method的位置，此处的代码意思是 |
| risk.table | 是否显示风险表 |
| risk.table.col | 设置风险表的颜色，此处设置为与分组（palette）颜色一致 |
| ggtheme | 设置主题，见ggplot2的主题设置 |
| palette | 设置分组颜色 |
| legend.title | 图例标题 |
| legend.labs | 两组的名称，此处设为"Male"和"Female" |
| xlab | 设置x轴名称 |
| xlim | 设置x轴的取值范围 |
| break.x.by | 设置x轴的刻度 |
| ylab | 设置y轴名称 |
My Survival analysis
3.3 我的作图函数
输入数据object0是如下格式的数据框：
| time | status | groups |
|---|---|---|
| xx | 1 | 1 |
| xx | 0 | 0 |
| …… | …… | …… |
file_name 是保存的文件名，我要输出的图片格式是tiff，所以此处应该是"file_name.tiff"这样形式，如果要保存为其他图片格式，需要修改函数中的代码。
其余的参数是我要修改的地方，也就是按我的偏好设置图的地方。
要注意的是这些值都需要设置，可以根据自己的偏好去修改函数，不要的删去，要的留下。
| ylab0 | xlab0 | break_x0 | xlim0 | legend_title0 | legend_labs0 |
|---|---|---|---|---|---|
| y轴名称 | x轴名称 | x轴刻度 | x轴显示范围 | 图例标题 | 分组标签 |
loop_survival_analysis <- function(object0, file_name, ylab0, xlab0, break_x0,
xlim0, legend_title0, legend_labs0){
library(survival)
library(survminer)
# survfit()-------------------------------------------------------------------
fit <- surv_fit(Surv(time, status) ~ groups, data = object0)
# cox-------------------------------------------------------------------------
cox <- coxph(Surv(time, status) ~ groups, data = object0)
res <- summary(cox)
x <- object0$time
# fig-------------------------------------------------------------------------
p<-ggsurvplot(fit,
pval = TRUE, pval.coord = c(0.7*max(x), 1),
pval.method = TRUE, pval.method.coord=c(0.4*max(x),1),
risk.table = TRUE, risk.table.col = "strata",
ggtheme = theme_linedraw(),
palette = c("darkblue", "red"),
legend.title = legend_title0,
legend.labs = legend_labs0,
xlab = xlab0, xlim = xlim0, break.x.by = break_x0,
ylab = ylab0
)
p$plot <- p$plot+
theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+
# ggtitle("My Survival Plot")+
theme(plot.title = element_text(hjust = 0.5))+
annotate("text",x=0.7*max(x),y=0.91, label=paste0("HR(95%CI) = ", round(res$conf.int[1],2), "(",round(res$conf.int[3],2),"-", round(res$conf.int[4],2), ")" ))+
annotate("text",x=0.7*max(x),y=0.83, label=paste0("C-index = ", round(res$concordance[1],2)))
p$table <- p$table+
theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())
tiff(file_name, width = 4290, height = 4930, res = 1000)
print(p)
dev.off()
}
三、其他
1. p值如何计算
Log-Rank test是无参数检验，近似于卡方检验
H0：两组生存过程相同，H1：两组生存过程不同
取出p值
surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)
p.value <- 1 - pchisq(surv_diff$chisq, length(surv_diff$n) -1)
surv_diffs 如下所示
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.001
