倾家荡产相求：关于weibull regression和Cox PH regression

cgdModProject2013.txt (4.6 KB)

2013-5-11 16:17:04 上传

data file

分析一个Gamma Interferon的Trial
背景见这里
http://www.nejm.org/doi/full/10.1056/NEJM199102213240801
数据里有变量 ID，Z1-Z9，T1，D
其中Z1是，T1是时间，D是event/censoring indicator
参见下面的代码里有其他变量的解释

1. rm(list=ls())
   
2. setwd(dir="c:/bios622/final")
   
3. library(survival)
   
4. library(MASS)
   
5. 
   
6. #load data
   
7. cgd <- read.table(file='cgdModProject2013.csv', header=T, sep=',')
   
8. names(cgd)
   
9. 
   
10. #variable labels of Z1-Z9, T1, and D
    
11. #won't apply due to 'if (FALSE)'
    
12. if (FALSE) {
    
13. names(cgd)[c(2:12)] <- c('Treatment Code',
    
14. 'Pattern of inheritance',
    
15. 'Age (year)',
    
16. 'Height (cm)',
    
17. 'Weight (kg)',
    
18. 'Using corticosteroids at time of study entry',
    
19. 'Using prophylactic antibiotics at time of study entry',
    
20. 'Gender',
    
21. 'Hospital Category',
    
22. 'Elapsed time (in days) from randomization to the diagnosis of a serious infection or censoring',
    
23. 'Censoring indicator'
    
24. )
    
25. }
    
26. 
    
27. #label values
    
28. cgd$Z1 <- factor(cgd$Z1, levels=c(2,1), labels=c('Placebo','Gamma Interferon'))
    
29. cgd$Z2 <- factor(cgd$Z2, levels=c(1,2), labels=c('X-linked','Autosomal recessive'))
    
30. cgd$Z6 <- factor(cgd$Z6, levels=c(2,1), labels=c('No','Yes'))
    
31. cgd$Z7 <- factor(cgd$Z7, levels=c(1,2), labels=c('Yes','No'))
    
32. cgd$Z8 <- factor(cgd$Z8, levels=c(1,2), labels=c('Male','Female'))
    
33. cgd$Z9 <- factor(cgd$Z9, levels=c(1,2,3,4), labels=c('US-Other','US-NIH','Europe-Amsterdam','Europe-Other'))
    
34. cgd$D <- factor(cgd$D, levels=c(1,2), labels=c('Occurred', 'Censored'))
    
35. 
    
36. #time-to-event setup
    
37. survData <- Surv(time=cgd$T1, event=(cgd$D=='Occurred'))
    
38. 
    
39. #overall
    
40. summary(survData)
    
41. km1 <- survfit(survData ~ 1, error="greenwood", conf.type="log-log", data=cgd)
    
42. #png(file='01.overall.png', width=1200, height=960, res=144, bg='transparent')
    
43. plot(km1, xlab="Survival time (days)", ylim=c(.25,1),
    
44. ylab="Estimated survival probability",
    
45. main="Figure 1. Overall Kaplan-Meier Estimator with pointwise CI", conf=T, mark=3, cex=0.5)
    
46. #dev.off()

复制代码

我现在知道把COX PH MODEL估计的概率画到KAPLAN-MEIER曲线上

1. kmZ1 <- survfit(survData ~ Z1, data=cgd)
   
2. plot(kmZ1, xlab="Time (days)", ylab="Estimated probability of no infection", ylim=c(0,1),
   
3. main="Figure 2. Kaplan-Meier Estimator with pointwise CI by Treatment", lty=1:2, mark=3, cex=0.5, col=c("red","blue"))
   
4. cox1 <- coxph(survData ~ Z1+Z2, data=cgd)
   
5. lines(survexp(~Z1, ratetable=cox1, data=cgd), col=c('purple', 'orange'), lty=c(3,4))

复制代码

请问
（1）如何把WEIBULL　REGRESSION拟合的概率也画上去，我在网上找的代码不管用
SOLVED

1. kmZ1 <- survfit(survData ~ Z1, data=cgd)
   
2. plot(kmZ1, xlab="Time (days)", ylab="Estimated probability of no infection", ylim=c(0,1),
   
3. main="Figure 2. Kaplan-Meier Estimator with pointwise CI by Treatment", lty=1:2, mark=3, cex=0.5, col=c("red","blue"))
   
4. #overlay predicted survival rate by COX PH model
   
5. cox1 <- coxph(survData ~ Z1+Z2, data=cgd)
   
6. lines(survexp(~Z1, ratetable=cox1, data=cgd), col=c('purple', 'orange'), lty=c(3,4))
   
7. #overlay predicated survival rate by Weibull regression
   
8. sWei <- survreg(survData ~ Z1,dist='weibull',data=cgd)
   
9. time <- predict(sWei, newdata=list(Z1=as.factor('Gamma Interferon')),type="quantile",p=seq(.01,.99,by=.01))
   
10. lines(time, 1-seq(.01,.99,.01), lty=5)
    
11. time <- predict(sWei, newdata=list(Z1=as.factor('Placebo')),type="quantile",p=seq(.01,.99,by=.01))
    
12. lines(time, 1-seq(.01,.99,.01), lty=6)

复制代码

（2）怎样确定WEIBULL和COX　PH哪个更适合这个数据？（500点）
（3）WEIBULL　REGRESSION怎么检查ASSUMPTION和DIAGNOSTICS？（414点）

看起来好高深

想采访一下楼主，倾家荡产后问题还没解决，论坛币也拿不回来的感觉怎么样是不是很酸爽

楼主，想请教一个问题，一篇论文中对变量的定义和数据来源的说明中，发现有些变量是时间序列数据，而有些变量是截面数据，论文的目的是用Weibull hazard model做生存分析。

比如 variable1  它的数据为1970-2015的时间序列数据
     variable2  它的数据为2000-2010年的平均值
     实在想不通它的数据结构是什么样子的？

是不是以下这种形式呢？先谢谢啦。

year  variable1               varible2
1970        数值                 缺失
1971        数值                 缺失
1972        数值                 缺失
1973        数值                 缺失
1974        数值                 缺失
1975        数值                 缺失
1976        数值                 缺失
1977        数值                 缺失
1978        数值                 缺失
…         …                  …
2000        数值        2000-2015年变量2的平均值
2001        数值        2000-2015年变量3的平均值
2002        数值        2000-2015年变量4的平均值
2003        数值        2000-2015年变量5的平均值
2004        数值        2000-2015年变量6的平均值
2005        数值        2000-2015年变量7的平均值
2006        数值        2000-2015年变量8的平均值
2007        数值        2000-2015年变量9的平均值
2008        数值        2000-2015年变量10的平均值
2009        数值        2000-2015年变量11的平均值
2010        数值        2000-2015年变量12的平均值
2011        数值        2000-2015年变量13的平均值
2012        数值        2000-2015年变量14的平均值
2013        数值        2000-2015年变量15的平均值
2014        数值        2000-2015年变量16的平均值
2015        数值        2000-2015年变量17的平均值