生存分析的基本内容
生成分析的目的：研究某对象某一时间某一事件的发生的概率，以及影响对象时间发生的因素。
生存分析通常被定义为一组用于分析数据的方法，其中结果变量是一个时间点到任何感兴趣事件发生的时间。
这个事件可能是死亡、疾病发生、婚姻、离婚等，或者任何与时间相关的事件。
使用生存分析的原因是它具备处理删失数据的条件（测量或观察的数据仅部分已知的条件），而其他技术（包括线性回归）不能够很好地解决这类问题。
假设： T T T为非负数随机变量，代表直到时间发生时的等待时间
- 【Def-1】时间发生时间小于 t t t的概率。 F ( t ) = P r ( T < t ) F(t)=Pr(T<t) F(t)=Pr(T<t)
- 【Def-2】生存函数：为患者、设备或其他感兴趣对象在任何给定时间内存活的概率的函数。（在 t t t时间内事物没有发生某件事的概率） S ( t ) = P r ( T ≥ t ) = 1 − F ( t ) = ∫ t ∞ ( x ) d x S(t)=Pr(T\geq t)=1-F(t)=\int_{t}^{\infty} (x) dx S(t)=Pr(T≥t)=1−F(t)=∫t∞(x)dx
- 【Def-3】风险函数：事物能够在某个时间点存活的概率。 λ ( t ) lim d t → 0 ( P r ( t ≤ T < t ∣ T ≥ t ) d t ) \lambda(t)\lim_{dt \to 0}\left(\frac{Pr(t\leq T<t|T\geq t)}{dt}\right) λ(t)dt→0lim(dtPr(t≤T<t∣T≥t))
删失数据：在研究某事物的观察过程中，该对象生存时间没有被完全观测到，造成生存数据不完整的现象。
删失数据一般分为三种：即左删失、右删失和区间删失。
●左删失（Left Censored）：指的是事件的发生时间只能确定在某一时间点之前。
●右删失（Right Censored）：指的是事件的发生时间只能确定在某一时间点之后。
●区间删失（Interval Censored）：指的是事件的发生时间只能确定在某一时间区间内。
生存分析方法通常有3种，包括非参数生存分析方法，半参数生存分析方法及参数生存分析方法，不同的方法有不同的使用条件。
使用R语言进行生存分析
使用的包
使用的数据集是lung，其已经包含在survival包中
install.packages("survival")
install.packages("survminer")
install.packages("flexsurv")
install.packages("ranger")
install.packages("ggfortify")
install.packages("knitr")
install.packages("ggpubr")
library(ggplot2)
library(ggpubr)
library(magrittr)
library(survival)
library(survminer)
library(knitr)
library(flexsurv)
library(ranger)
library(ggfortify)
> attach(lung)
> head(lung,5)
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
● inst：机构的代码。
● time：生存的时间，单位是天。
● status：1表示没有死亡，2表示死亡。
● age：年龄。● sex：1表示男性，2表示女性。
● ph.ecog：一个表现分数。0表示最好，5表示最糟。
● ph.karno：医生评价的Karnofsky评价。0表示最早，100表示最好。
● pat.karno：患者评价的Karnofsky评价。0表示最早，100表示最好。
● meal.cal：用餐时消耗的卡路里。
● wt.loss：最近六个月下降的体重数。
这一份数据集不需要进行进一步的处理，可以直接使用。
非参数模型
Surv( )函数创建生存对象
survfit( )函数后见Kaplan-Meier模型
> lung$SurvObj <- with(lung,Surv(time = time,event = status))
> head(lung,5)
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss SurvObj
1 3 306 2 74 1 1 90 100 1175 NA 306
2 3 455 2 68 1 0 90 90 1225 15 455
3 3 1010 1 56 1 0 90 90 NA 15 1010+
4 5 210 2 57 1 1 90 60 1150 11 210
5 1 883 2 60 1 0 100 90 NA 0 883
> km.by.sex <- survfit(Surv(time,status)~sex,data=lung)
> km.by.one <- survfit(Surv(time,status)~1,data=lung)
> km.by.one
Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
n events median 0.95LCL 0.95UCL
[1,] 228 165 310 285 363
> summary(km.by.one,times = c(1,50,100,200,300*(1:3)))
Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 228 0 1.0000 0.0000 1.0000 1.000
50 217 11 0.9518 0.0142 0.9243 0.980
100 196 20 0.8640 0.0227 0.8206 0.910
200 144 41 0.6803 0.0311 0.6219 0.744
300 92 29 0.5306 0.0346 0.4669 0.603
600 24 47 0.2136 0.0335 0.1571 0.290
900 3 17 0.0503 0.0228 0.0207 0.123
> km.by.sex
Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
n events median 0.95LCL 0.95UCL
sex=1 138 112 270 212 310
sex=2 90 53 426 348 550
> summary(km.by.sex,times = c(1,50,100,200,300*(1:3)))
Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
sex=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 138 0 1.0000 0.0000 1.0000 1.000
50 128 10 0.9275 0.0221 0.8853 0.972
100 114 14 0.8261 0.0323 0.7652 0.892
200 78 30 0.6073 0.0417 0.5309 0.695
300 49 20 0.4411 0.0439 0.3629 0.536
600 13 29 0.1451 0.0353 0.0900 0.234
900 2 9 0.0357 0.0216 0.0109 0.117
sex=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 90 0 1.0000 0.0000 1.0000 1.000
50 89 1 0.9889 0.0110 0.9675 1.000
100 82 6 0.9221 0.0283 0.8683 0.979
200 66 11 0.7946 0.0432 0.7142 0.884
300 43 9 0.6742 0.0523 0.5791 0.785
600 11 18 0.3433 0.0634 0.2390 0.493
900 1 8 0.0832 0.0499 0.0257 0.270
> ggsurvplot(km.by.one)
> ggsurvplot(km.by.sex)
半参数模型生存分析方法
半参数模型同样不会对生存函数或者危险函数的形式作出任何假设，但是其对于协变量存在一个很强的假设。
使用Cox模型来解决生存分析问题的时候，这里有两个强假设需要满足，第一个是对数线性假定（模型中的协变量应与对数风险比呈线性关系），第二个是比例风险假定（各危险因素的作用不随时间的变化而变化）。
coxph( )构建Cox模型
> cox <- coxph(Surv(time,status)~age+sex+ph.karno+wt.loss,data = lung)
> cox
Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.karno + wt.loss,
data = lung)
coef exp(coef) se(coef) z p
age 0.015140 1.015255 0.009837 1.539 0.12379
sex -0.513955 0.598125 0.174410 -2.947 0.00321
ph.karno -0.012871 0.987211 0.006184 -2.081 0.03741
wt.loss -0.002246 0.997757 0.006357 -0.353 0.72389
Likelihood ratio test=18.84 on 4 df, p=0.0008436
n= 214, number of events= 152
(因为不存在，14个观察量被删除了)
cox.zph( )检验Cox模型的比例风险假设
> res <- cox.zph(cox)
> plot(res)
> plot(res)
autoplot( )函数 进行可视化
> autoplot(survfit(cox))
survfit( )进行预测，第一个参数是构建好的模型，第二个参数是需要预测的新数据
> pc <- survfit(cox,newdata = lung[2,])
> summary(pc,times = c(1,50,100,200,300*(1:3)))
Call: survfit(formula = cox, newdata = lung[2, ])
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 214 0 1.000 0.0000 1.0000 1.000
50 204 10 0.948 0.0170 0.9147 0.982
100 186 17 0.857 0.0291 0.8021 0.916
200 140 35 0.667 0.0446 0.5850 0.760
300 89 29 0.491 0.0523 0.3984 0.605
600 24 44 0.176 0.0462 0.1053 0.294
900 3 17 0.039 0.0229 0.0123 0.123
> autoplot(pc)
如果有分类变量不满足比例风险假定，我们可以使用分层Cox回归模型。strata( )函数
结果和书上不一样，暂时没找到原因，图也不对
> cox <- coxph(Surv(time=time,time2=status)~age+sex+strata(ph.karno)+wt.loss,data = lung)
> summary(cox)
Call:
coxph(formula = Surv(time = time, time2 = status) ~ age + sex +
strata(ph.karno) + wt.loss, data = lung)
n= 214, number of events= 152
(因为不存在，14个观察量被删除了)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.018253 1.018421 0.010021 1.821 0.068536 .
sex -0.603327 0.546989 0.180373 -3.345 0.000823 ***
wt.loss -0.005364 0.994651 0.006693 -0.801 0.422881
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0184 0.9819 0.9986 1.039
sex 0.5470 1.8282 0.3841 0.779
wt.loss 0.9947 1.0054 0.9817 1.008
Concordance= 0.615 (se = 0.028 )
Likelihood ratio test= 15.38 on 3 df, p=0.002
Wald test = 14.58 on 3 df, p=0.002
Score (logrank) test = 14.94 on 3 df, p=0.002
> pre <- survfit(cox,newdata = lung[2,])
>
> pre
Call: survfit(formula = cox, newdata = lung[2, ])
n events median 0.95LCL 0.95UCL
[1,] 68 43 345 268 457
> autoplot(pre)
参数模型
AIC值越小，模型拟合得越好。
flexsurv包中flexsurvreg( )
AIC函数
> attach(lung)
> S <- Surv(time,status)
> Dist <- c("exp","weibull","lnorm","gamma","gompertz",
+ "gengamma","gengamma.orig","genf","genf.orig","llogis")
> AIC <- matrix(ncol = 2,nrow = 10)
> for(i in 1:10){
+ AIC[i,1] <- Dist[i]
+ model <- flexsurvreg(S~1,dist = Dist[i])
+ AIC[i,2] <- AIC(model)
+ }
>
> colnames(AIC) <- c("Distribution","AIC")
> AIC <- base::transform(AIC,Distribution=as.character(Distribution),
+ AIC=as.factor(AIC))
> AIC$AIC <- as.numeric(levels(AIC$AIC)[AIC$AIC])
> AIC[order(AIC$AIC),]
Distribution AIC
2 weibull 2311.702
7 gengamma.orig 2313.380
6 gengamma 2313.380
4 gamma 2313.469
5 gompertz 2314.711
8 genf 2315.153
9 genf.orig 2315.153
10 llogis 2325.862
1 exp 2326.676
3 lnorm 2342.538
> weibmodel <- flexsurvreg(S~1,dist="weibull")
> plot(weibmodel,ylab="Survival probability",xlab="Time",main="Weibull Survival Plot")
> legend("topright",legend = c("KM Plot","Fitted"),lty = c(1,1),col = c("black","red"),cex=0.75)
这里比较特殊，需要使用summary( )函数进行预测，第一个参数是构建好的模型；第二个参数是newdata，是需要预测的数据集；第三个参数是t，表示预测第多少天的生存概率。
> pre <- summary(weibmodel,newdata = lung[2,],t=c(100,200,300,400,500,600))
> pre
time est lcl ucl
1 100 0.8588388 0.8225402 0.8949835
2 200 0.6844801 0.6334088 0.7376424
3 300 0.5238260 0.4679668 0.5818579
4 400 0.3889118 0.3333643 0.4491276
5 500 0.2816790 0.2264905 0.3396391
6 600 0.1997281 0.1467537 0.2556226
随机生成森林模型
ranger包ran在这里插入图片描述
ger( )函数
> r_fit <- ranger(Surv(time,status)~age+sex+ph.karno+wt.loss,data = na.omit(lung),mtry = 4,importance = "permutation",splitrule = "extratrees",verbose = TRUE)
>
> death_times <- r_fit$unique.death.times
> surv_prob <- data.frame(r_fit$survival)
> avg_prob <- sapply(surv_prob, mean)
>
> plot(r_fit$unique.death.times,r_fit$survival[1,],
+ type="l",
+ ylim = c(0,1),
+ col="red",
+ xlab="days",
+ ylab="survival",
+ main="Patient Survival Curves")
>
> cols <- colors()
> for (n in sample(c(2:dim(veteran)[1]),20)) {
+ lines(r_fit$unique.death.times,r_fit$survival[n,],type = "l",col=cols[n])
+ }
> lines(death_times,avg_prob,lwd=2)
> legend(500,0.7,legend = c("Average=black"))
1980
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
