R语言实现非等比例风险生存资料分析（1）

1. #非等比例风险的生存资料分析
2. ###1 生成模拟数据###
3. library(flexsurv)
4. set.seed(123)
5. # 生成样本数量
6. n <- 100
7. # 生成时间数据
8. time <- sample(1:1000,n,replace=F)  # 调整shape和scale参数以控制生存曲线形状
9. # 生成事件数据（假设按比例风险模型）
10. status <- rbinom(n, size = 1, prob = plogis(-0.5 + 0.01 * time))
11. # 生成一些协变量数据
12. covariate1 <- sample(c(0,1,100),n,replace=T)
13. covariate2 <- rbinom(n, size = 3, prob = 0.5)
14. group<-sample(c(0,1),n,replace=T)
15. # 创建生存数据框
16. survival_data <-data.frame(time,status,covariate1,covariate2,group)
18. ###2 绘制KM曲线###
19. library(survminer)
20. library(survival)
21. library(ggplot2)
22. fit<-survfit(Surv(time,status) ~group,data=survival_data)
23. ggsurvplot(fit,title="图1")
25. ###3 检验等比例风险###
26. #观察图片
27. plot(fit,fun='cloglog',col=c("red","green"),xlab="生存时间对数",ylab="二次对数生存率")
28. #残差法
29. coxfit<-coxph(Surv(time,status) ~group,data=survival_data)
30. zph<-cox.zph(coxfit,transform = "km")
31. ggcoxzph(zph,title="图2")#不满足等比例风险
32. cox.zph(coxfit,transform = "rank")
33. cox.zph(coxfit,transform = "identity")
35. ###4 two-stage 组间差异比较###
36. library(TSHRC)
37. twostage(time,
38.          status,group,
39.          nboot=1000)#仅可以说明组间存在差异
40. #LRPV:log-rank检验；MTPV第二阶段检验；TSPV：两阶段检验结果
41. #install.packages("ComparisonSurv")
42. library(ComparisonSurv)
43. help(package="ComparisonSurv")
44. Overall.test(time,status,group)#包含two-stage在内的9种方法
46. ###5 比较大小###
47. # 找出交叉点对应的时间
48. crosspoint(time,status,group)#768/774/775/999
49. #使用ComparisonSurv中的长期和短期比较
50. Long.test(time,status,group,t0=775)
51. Short.test(time,status,group,t0=775)
52. #RMST法:限制平均生存时间
53. #使用surv2sampleComp
54. library(surv2sampleComp)
55. help(package="surv2sampleComp")
56. sur1<-surv2sample(time,status,group,tau_start=0,tau=775,procedure="KM")
57. max(time);min(time)
58. sur2<-surv2sample(time,status,group,tau_start=775,tau=938,procedure="KM")#tau为两组中最长观测值的最小值
59. #使用survRM2
60. library(survRM2)
61. help(package="survRM2")
62. sur3<-rmst2(survival_data$time,
63.       survival_data$status, arm=survival_data$group,
64.       tau = 755, alpha = 0.05)
65. plot(sur3)
66. #landmark方法
67. dat1<-survival_data[survival_data$time<775,]
68. pzg<-coxph(Surv(time,status) ~group,data=dat1) %>% summary
69. pz<-pzg[["coefficients"]][1,5]
70. hr<-pzg[["coefficients"]][1,1]
71. dat2<-survival_data[survival_data$time>=775,]
72. pzg2<-coxph(Surv(time,status) ~group,data=dat2) %>% summary
73. pz2<-pzg2[["coefficients"]][1,5]
74. hr2<-pzg2[["coefficients"]][1,1]
75. duan_1<-survfit(Surv(time,status)~group,dat1)
76. duan_2<-survfit(Surv(time,status)~group,dat2)
77. #绘图
78. plot(duan_1, lty = c('solid', 'dashed'), col = c('red', 'blue'), xlim = c(0, 1000), xlab = "生存时间(天)", ylab = '累计生存率', axes = FALSE, lwd = 2)
79. text(695, 1.0, paste('P=',round(pz,3)))
80. text(695, 0.8, paste('HR=',round(hr,3)))
81. axis(1)
82. axis(2)
83. axis(4)
84. duan2a <- survfit(Surv(dat2$time, dat2$status == "1") ~ dat2$group, subset = (dat2$group == 1))
85. duan2b <- survfit(Surv(dat2$time, dat2$status == "1") ~ dat2$group, subset = (dat2$group == 0))
86. lines(c(775, duan2b$time), c(1, duan2b$surv), type = "s", lty = 'solid', lwd = 2, col = 'red')
87. lines(c(775, duan2a$time), c(1, duan2a$surv), type = "s", lty = 'dashed', lwd = 2, col = 'blue')
88. text(940, 1.0, paste('P<0.01'))
89. text(940, 0.8, paste('HR=',round(hr2,3)))
90. abline(v =775, lty = 2, col = "darkgreen", lwd = 3)

复制代码

 
 
 
 

来源：https://blog.csdn.net/weixin_49320263/article/details/132323890
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！