R进行生存资料分析学习笔记

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

与SAS 的proc lifetest及proc phreg相比，R进行生存资料分析的相对优势体现在：

1、能够更加方便快捷绘制Kaplan-Meier曲线，
2、能够相对更加美观地绘制Kaplan-Meier曲线
3、相对容易给出中位生存时间及其置信区间
4、检测比例风险假设的方法更加灵活且完善
5、易于构建时间依赖模型

以下是R绘制的Kaplan-Meier曲线图形示例：


以下是学习笔记的编码，供参考，欢迎交流：

1. library(stats)
   
2. library(survival)
   
3. 
   
4. ## Information of data
   
5. 
   
6. data(package = "survival")  # List datasets in survival package
   
7. help(bladder1)              # Description of data
   
8. head(bladder1)              # Show first 6 rows
   
9. str(bladder1)               # Check type of variables
   
10. summary(bladder1)           # Statistical summary
    
11. 
    
12. ## Get the final data with nonzero follow-up
    
13. bladder1$time <-
    
14.   as.numeric(bladder1$stop -
    
15.                bladder1$start)
    
16. summary(bladder1$time)
    
17. bladder1 <- subset(bladder1,status<=1 & time>0)
    
18. 
    
19. 
    
20. ## Step1 Create Kaplan-Meier curve and estimate median survial/event time
    
21. ## The "log-log" confidence interval is preferred.
    
22. ## Create overval Kaplan-Meier curve
    
23. km.as.one <- survfit(Surv(time, status) ~ 1, data = bladder1,
    
24.                      conf.type = "log-log")
    
25. 
    
26. ## Create Kaplan-Meier curve stratified by treatment
    
27. km.by.trt <- survfit(Surv(time, status) ~ treatment, data = bladder1,
    
28.                      conf.type = "log-log")
    
29. 
    
30. ## Show simple statistics of Kaplan-Meier curve
    
31. km.as.one
    
32. km.by.trt
    
33. 
    
34. ## See survival estimates at given time (lots of outputs)
    
35. summary(km.as.one)
    
36. summary(km.by.trt)
    
37. 
    
38. ## Plot Kaplan-Meier curve without any specification
    
39. plot(km.as.one)
    
40. plot(km.by.trt)
    
41. 
    
42. ## Plot Kaplan-Meier curve Without confidence interval and mark of event
    
43. plot(km.as.one, conf = F, mark.time = F)
    
44. plot(km.by.trt, conf = F, mark.time = F)
    
45. 
    
46. 
    
47. ## step2 Create a simple cox regression and estimate HR:
    
48. model1 <-coxph(Surv(time, status) ~ treatment + number +size,
    
49.                data=bladder1)
    
50. 
    
51. ## Model output
    
52. summary(model1)                    # Output summary information
    
53. confint(model1)                    # Output 95% CI for the coefficients
    
54. exp(coef(model1))                  # Output HR (exponentiated coefficients)
    
55. exp(confint(model1))               # 95% CI for exponentiated coefficients
    
56. predict(model1, type="risk")       # predicted values
    
57. residuals(model1, type="deviance") # residuals
    
58. 
    
59. 
    
60. ## Step3 Check for violation of proportional hazard (constant HR over time)
    
61. model1.zph <- cox.zph(model1)
    
62. model1.zph                    
    
63. ## Note: p value of treatmentthiotepa <0.05
    
64. ## GLOBAL p value is more important than p values of treatmentthiotepa
    
65. ## Only GLOBAL p value <0.05 is neccesary for solution of violation
    
66. ## Following steps are just for providing examples
    
67. 
    
68. ## Displays a graph of the scaled Schoenfeld residuals, along with a smooth curve.
    
69. plot(model1.zph)
    
70. 
    
71. ## Step4 Solution of violating of proportional hazard
    
72. 
    
73. ## Solution 1 Stratify: create a stratified cox model,namely "contitional cox model "
    
74. model2 <-coxph(Surv(start, stop, status) ~ number +size + strata (treatment),
    
75.                data=bladder1)
    
76. summary(model2)
    
77. model2.zph <- cox.zph(model2)
    
78. model2.zph
    
79. 
    
80. ## Solution 1 Create a simple time-dependent model without time-dependent variable
    
81. model3 <-coxph(Surv(start, stop, status) ~ number +size + treatment,
    
82.                data=bladder1)
    
83. summary(model3)
    
84. 
    
85. ## Create a time-dependent model with time-dependent variable
    
86. bladder1$thiotepa <-
    
87.   as.numeric (bladder1$treatment=="thiotepa")
    
88. model4 <-coxph(Surv(start, stop, status) ~ number +size + treatment
    
89.                + thiotepa:time,     ## Note: "thiotepa:time" but not "thiotepa*time"
    
90.                data=bladder1)
    
91. summary(model4)
    
92. 
    
93. ## Step5 Performs stepwise model selection by AIC
    
94. model5 <-step(model4)  ## Default direction is "backward"
    
95. summary(model5)
    
96. 
    
97. 
    
98. ## rms::survplot() function
    
99. ## Load rms package
    
100. library(rms)
     
101. 
     
102. ## Plot Kaplan-Meier curve without any specification
     
103. fit1 <- npsurv(Surv(time, status) ~ 1, data = bladder1)
     
104. fit2 <- npsurv(Surv(time, status) ~ treatment, data = bladder1)
     
105. survplot(fit1)
     
106. survplot(fit2)
     
107. 
     
108. ## Plot Kaplan-Meier curve Without confidence interval
     
109. survplot(fit1, conf = "none")
     
110. survplot(fit2, conf = "none")
     
111. 
     
112. ## More options for Kaplan-Meier curve
     
113. survplot(fit2,
     
114.          xlab = "Time, month",                  # add x-axis label
     
115.          ylab = "Survival probability, %",      # add y-axis label
     
116.          ## xlim=c(0,60),                       # add x-axis limits
     
117.          ## ylim=c(-0.1,1),                     # add y-axis limits
     
118.          ## conf.int=.95,                       # show 95% CI,
     
119.          conf='none',                           # change type of CI
     
120.          label.curves = list(keys = "lines"),   # legend instead of direct label
     
121.          levels.only  = TRUE,                   # show only levels, no label
     
122.          col=c('red','green','blue'),           # change legend color
     
123.          ## fun = function(x) {1 - x},          # Cumulative probability plot
     
124.          ## loglog   = TRUE,                    # log(-log Survival) plot
     
125.          ## logt     = TRUE,                    # log time
     
126.          time.inc = 5,                          # time increment
     
127.          ## dots     = TRUE,                    # dot grid
     
128.          n.risk   = TRUE,                       # show number at risk
     
129.          y.n.risk = 0.01,                       # Change position of number at risk
     
130.          cex.n.risk = 0.6                       # change character size for number at risk
     
131.          )

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏13 回帖

关键词：资料分析 学习笔记 分析学 习笔记 information 资料

竟然可没有点赞么？这个小教程写得真心不错

太厉害了！

谢谢楼主 非常厉害 能否将上图的原始数据发一份 弄不清楚代码里面的哪些参数是指X 哪些是生存事件变量？想运用这些代码去做别的分析，谢谢！

moonstone 发表于 2015-12-11 21:32
与SAS 的proc lifetest及proc phreg相比，R进行生存资料分析的相对优势体现在：

1、能够更加方便快捷绘制 ...

谢谢你

支持楼主，互相学习！

谢谢楼主！想请教一下，我在编写程序时，用了conf.type="log-log",为什么出来的图像只是普通的曲线而不是log-log plot呢？

> plot(km.as.one, conf = F, mark.time = F)
Error in plot.survfit(km.as.one, conf = F, mark.time = F) :
  argument 2 matches multiple formal arguments

此步骤运行错误，问题出在哪儿？

survplot和plot在使用上哪个更好呢？
可信区间如何去掉的呢？
legend是怎么加进去的呢

lanhong1993 发表于 2017-4-18 10:42
> plot(km.as.one, conf = F, mark.time = F)
Error in plot.survfit(km.as.one, conf = F, mark.time = F ...

因为包升级了，conf选项变成其他的，改成conf.int = F就好啦