【R语言-生存分析之观察生存率计算】_cumulative mortality r语言

生存分析

全人群肿瘤登记资料常用的统计分析指标包括发病率、死亡率、现患率以及生存率等，其中肿瘤生存率是评估肿瘤治疗效果和肿瘤负担的必要指标，其计算涉及肿瘤患者的发病、 死亡和随访三个方面的资料， 数据整理和计算过程均较为复杂， 如何及时、 准确地计算肿瘤生存率， 并使其可以在不同地区、 不同人群、不同时期间被客观公正地比较和评价，仍然是统计学家们一直在研究的主题。人群肿瘤登记生存率的分析指标包括观察生存率（observed survival rate，OSR） 和相对生存率（relativesurvival rate，RSR）等，本文首先利用肿瘤随访数据来展示计算观察生存率和绘制生存曲线的方法。（相对生存率计算方法见下篇）

1.生存分析的相关概念

死亡率 mortality rate：某单位时间里的平均死亡强度频率
死亡概率mortality proportion:为q，往后一个时段内死亡的可能性大小
生存概率survival probability：死亡概率的对应面，往后活满1个时间段的可能性大小
生存率：S(t)病人经历给定时间t个单位时间后仍存活的概率，实质为累积生存概率(cumulative probability of survival)

2.生存分析的特点

观察截止时间：无法等到每个人都死亡
删失（截尾）censoring:我们知道某个体在某时间点存活，但是无法知道其确切的生存时间
截尾数据原因：
随访终止时间结束，终点事件没有发生
随访终止之前，患者失访
随访终止之前，患者退出研究

3.人群生存研究的常用指标

观察生存率 observed survival
净生存率 net survival
相对生存率 relative survival ratio
期望生存率 expect survival proportion
粗死亡概率crude probability of death

4.生存率概念

观察生存率
-终点事件：全部死亡（死于癌症和其他）
-所有随访对象中，生存期大于等于该时间的研究对象的比例，如五年生存率
-包含普通人群的死亡概率和癌症患者预后两部分信息
反映肿瘤患者的整体死亡状况。
净生存率（用相对生存率估计）
-终点事件：死于癌症
-若肿瘤患者死于肿瘤之外的其他原因，将与存活状态同等处理。
-理论上的数值，是癌症生存率比较的关键指标
比较不同年龄、性别、社会经济学状况下癌症患者的生存率时，使用净生存率显得尤为重要，因为肿瘤外其它死因会影响癌症患者的生存状况。

5.实际操作代码实现

以xx市肿瘤随访数据库为例计算观察生存率，其原始数据包括性别、年龄、出生日期、肿瘤诊断日期、ICD10编码、生存状态、死亡原因、死亡日期等变量（见下图）。

5.1 原始资料转换为生存分析数据

#####生存分析前期准备####
rm(list = ls())
setwd("E://surviVal")
if(!requireNamespace("stringr","lubridate","purrr","dplyr","ggplot2","readxl","survival","relsurv","survminer",quietly=TRUE))
 install.packages(c("stringr","lubridate","purrr","dplyr","ggplot2","readxl","survival","relsurv","survminer"))
##载入软件包##
library(stringr)
library(purrr)
library(dplyr)
library(ggplot2)
library(readxl)
library(readxl)
library(survival)
library(survminer)
library(relsurv)
library(lubridate)
##读取原始数据##
read_xlsx("数据整理//原始数据.xlsx")->sv_data###此处可替换各区数据###
##整理转换为生存数据##
sv<-sv_data%>%data.frame()
lapply(sv[,c(1,2,9,10)],as.numeric)->sv[,c(1,2,9,10)]
lapply(sv[,c(3,4,8,11)],as.Date)->sv[,c(3,4,8,11)]
sv<-mutate(sv,inciden_year=year(Inciden),start=Inciden,outcome=Stat*1)
sv$deathcause<-ifelse(is.na(sv$Caus)==T,sv$deathcause<-0,sv$deathcause<-sv$Caus*1)
sv$failure_crude<-ifelse(sv$outcome==2,sv$failure_crude<-1,sv$failure_crude<-0)
sv$failure_adjust<-ifelse(sv$outcome==2&sv$deathcause==1,sv$failure_adjust<-1,sv$failure_adjust<-0)
sv$terminal<-rep(as.Date("2021-12-31"),length(sv$Sex))##可在此处修改观察截止日期“YY-MM-DD"
sv$end<-ifelse(sv$outcome==2,sv$end<-sv$Deathda,sv$end<-sv$terminal)%>%as.Date(origin="1970-1-1")
sv<-sv%>%mutate(surv_year=((end-start)/365.25)%>%as.numeric(),surv_month=((end-start)/30)%>%as.numeric(),surv_day=difftime(end,start,units="days")%>%as.numeric())
sv$area<-ifelse(sv$cDepID%in%c("大港区","汉沽区","塘沽区"),sv$area<-"滨海新区",sv$area<-sv$cDepID)
sv$ICD<-str_sub(sv$Icd10,1,3)
ccdd<-read.table("初始数据//ccdd.txt")
sv<-left_join(sv,ccdd,by=c("ICD"="ICD"))
source("初始数据//function_trans.R")
sv$agegroup<-lapply(sv$Age,age_trans)%>%as.character()
sv$cx<-lapply(sv$area,city_rural)%>%as.character()
sv$agenumber<-lapply(sv$agegroup,agegtrans)%>%as.numeric()
sv$cxnumber<-lapply(sv$cx,cxtrans)%>%as.numeric()
gc()##清除缓存##
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

5.2 转换后数据如下

其中failure_adjust和surv_year是进行生存分析的主要变量

5.3 构建生存分析模型并保存输出分析结果

#####构建生存分析模型######
###全人群生存分析模型###
dir.create("输出结果")
sink(file = "输出结果//surv_result.txt", append = T, type = c("output", "message"),split = F)##建立输出结果存储文本##
fit.all<-survfit(Surv(surv_year,failure_crude)~1,data=sv)
cat("全人群生存分析结果\n")
cat("平均生存时间，中位生存时间\n")
print(fit.all,print.rmean = T)##查看中位生存时间、受限的平均生存时间##
cat("1/3/5年生存率\n")
summary(fit.all,times = c(1,3,5))##此处可定义查看期间年度生存率times=c(n)##
###分年份生存分析模型###
fit.year<-survfit(Surv(surv_year,failure_crude)~inciden_year,data=sv)
cat("分年份生存分析结果\n")
cat("平均生存时间，中位生存时间\n")
print(fit.year,print.rmean = T)
cat("1/3/5年生存率\n")
summary(fit.year,times = c(1,3,5))
cat("logrank模型结果\n")
survdiff(Surv(surv_year,failure_crude)~inciden_year,data=sv)##logrank模型比较不同组的生存曲线（生存率）差异##
###分性别生存分析模型##
fit.sex<-survfit(Surv(surv_year,failure_crude)~Sex,data=sv)
cat("分性别生存分析结果\n")
cat("平均生存时间，中位生存时间\n")
print(fit.sex,print.rmean = T)
cat("1/3/5年生存率\n")
summary(fit.sex,times = c(1,3,5))
cat("logrank模型结果\n")
survdiff(Surv(surv_year,failure_crude)~Sex,data=sv)
###分城乡生存分析模型##
fit.cx<-survfit(Surv(surv_year,failure_crude)~cx,data=sv)
cat("分城乡生存分析结果\n")
cat("平均生存时间，中位生存时间\n")
print(fit.cx,print.rmean = T)
cat("1/3/5年生存率\n")
summary(fit.cx,times = c(1,3,5))
cat("logrank模型结果\n")
survdiff(Surv(surv_year,failure_crude)~cx,data=sv)
###分年龄组生存分析模型##
fit.agegroup<-survfit(Surv(surv_year,failure_crude)~agegroup,data=sv)
cat("分年龄组生存分析结果\n")
cat("平均生存时间，中位生存时间\n")
print(fit.agegroup,print.rmean = T)
cat("1/3/5年生存率\n")
summary(fit.agegroup,times = c(1,3,5))
cat("logrank模型结果\n")
survdiff(Surv(surv_year,failure_crude)~agegroup,data=sv)
###分癌种生存分析模型##
fit.class<-survfit(Surv(surv_year,failure_crude)~CCDD,data=sv)
cat("分癌种生存分析结果\n")
cat("平均生存时间，中位生存时间\n")
print(fit.class,print.rmean = T)
cat("1/3/5年生存率\n")
summary(fit.class,times = c(1,3,5))
cat("logrank模型结果\n")
survdiff(Surv(surv_year,failure_crude)~CCDD,data=sv)
##cox风险比例模型探索某个变量对生存率的影响显著程度##
cat("cox风险比例模型输出结果\n")
coxph(Surv(surv_year,failure_crude)~Sex+Age+inciden_year+cx+CCDD+agenumber,data = sv)##p<0.05具有显著影响##
sink()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

5.4 输出结果如下（举例）

可见分性别平均生存时间（remean）、中位生存时间（median）、1/3/5年生存率、logrank生存曲线差异检验、cox风险比例等分析结果【可直接摘出相应指标分析撰写论文】

5.5 绘制生存曲线批量输出

#####绘制生存曲线图#####
ggsurvplot(fit.sex,data=sv,pval = T,pval.size=4,pval.coord=c(1,0.3),conf.int=T,surv.median.line = "hv",
           xlim=c(1,13),break.time.by=1,xlab="Time in years",
           palette = c("#E7B800","#2E9FDF"),
           legend.labs=c("male","female"),legend.title="Sex",title="Sex Survival Curve",
           ggtheme = theme(axis.text.x=element_text(angle = 0)))
###批量出图###
dir.create("输出图片")
nm<-c("分性别生存曲线图","分年龄组生存曲线图","分城乡生存曲线图","分年份生存曲线图","分癌症生存曲线图")
ff<-c("fit.sex","fit.agegroup","fit.cx","fit.year","fit.class")
for(i in 1:length(nm)){
  png(filename = paste0("输出图片//",nm[i],".jpg"))
  print(ggsurvplot(eval(parse(text=ff[i])),data=sv,pval = T,pval.size=4,pval.coord=c(1,0.3),conf.int=T,surv.median.line = "hv"))##循环中ggplot2绘图结果不能直接显示出来，需要用print()函数将结果显现
  dev.off()
1
2
3
4
5
6
7
8
9
10
11
12
13
14

5.5 输出结果如下