热门标签
热门文章
- 1ssh ServerAliveInterval 设置保持连接
- 2CWE-191: Integer Underflow (Wrap or Wraparound)整数下溢
- 3linux下的socket在哪个头文件,socket函数在哪个文件定义，解决方案
- 4springboot整合Redis详解完整篇
- 5使用esp32+micropython+microdot搭建web(http+websocket)服务器（超详细）第三部分_micropython web
- 6踩内存（内存重叠）的处理_uclibc 踩内存问题
- 7【Verifying transaction: failed】
- 87-liunx服务器规范
- 9基于uniapp的仿书旗小说的小程序_uniapp小说阅读器
- 10C# wpf 使用GDI+实现截屏_wpf截图
当前位置: article > 正文
【R语言】生存分析_r语言生存分析
作者：菜鸟追梦旅行 | 2024-02-23 00:06:55
赞
踩
r语言生存分析
生存分析的基本内容
生成分析的目的：研究某对象某一时间某一事件的发生的概率，以及影响对象时间发生的因素。
生存分析通常被定义为一组用于分析数据的方法，其中结果变量是一个时间点到任何感兴趣事件发生的时间。
这个事件可能是死亡、疾病发生、婚姻、离婚等，或者任何与时间相关的事件。
使用生存分析的原因是它具备处理删失数据的条件（测量或观察的数据仅部分已知的条件），而其他技术（包括线性回归）不能够很好地解决这类问题。
假设： 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)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
> 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
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
● 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
- 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
> ggsurvplot(km.by.one)
> ggsurvplot(km.by.sex)
- 1
- 2
半参数模型生存分析方法
半参数模型同样不会对生存函数或者危险函数的形式作出任何假设，但是其对于协变量存在一个很强的假设。
使用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个观察量被删除了)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
cox.zph( )检验Cox模型的比例风险假设
> res <- cox.zph(cox)
> plot(res)
> plot(res)
- 1
- 2
- 3
autoplot( )函数 进行可视化
> autoplot(survfit(cox))
- 1
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)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
如果有分类变量不满足比例风险假定，我们可以使用分层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)
- 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
参数模型
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
- 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
> 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)
- 1
- 2
- 3
这里比较特殊，需要使用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
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
随机生成森林模型
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"))
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/菜鸟追梦旅行/article/detail/132203
推荐阅读
相关标签
