分析一个Gamma Interferon的Trial
背景见这里
http://www.nejm.org/doi/full/10.1056/NEJM199102213240801
数据里有变量 ID，Z1-Z9，T1，D
其中Z1是，T1是时间，D是event/censoring indicator
参见下面的代码里有其他变量的解释
- rm(list=ls())
- setwd(dir="c:/bios622/final")
- library(survival)
- library(MASS)
- #load data
- cgd <- read.table(file='cgdModProject2013.csv', header=T, sep=',')
- names(cgd)
- #variable labels of Z1-Z9, T1, and D
- #won't apply due to 'if (FALSE)'
- if (FALSE) {
- names(cgd)[c(2:12)] <- c('Treatment Code',
- 'Pattern of inheritance',
- 'Age (year)',
- 'Height (cm)',
- 'Weight (kg)',
- 'Using corticosteroids at time of study entry',
- 'Using prophylactic antibiotics at time of study entry',
- 'Gender',
- 'Hospital Category',
- 'Elapsed time (in days) from randomization to the diagnosis of a serious infection or censoring',
- 'Censoring indicator'
- )
- }
- #label values
- cgd$Z1 <- factor(cgd$Z1, levels=c(2,1), labels=c('Placebo','Gamma Interferon'))
- cgd$Z2 <- factor(cgd$Z2, levels=c(1,2), labels=c('X-linked','Autosomal recessive'))
- cgd$Z6 <- factor(cgd$Z6, levels=c(2,1), labels=c('No','Yes'))
- cgd$Z7 <- factor(cgd$Z7, levels=c(1,2), labels=c('Yes','No'))
- cgd$Z8 <- factor(cgd$Z8, levels=c(1,2), labels=c('Male','Female'))
- cgd$Z9 <- factor(cgd$Z9, levels=c(1,2,3,4), labels=c('US-Other','US-NIH','Europe-Amsterdam','Europe-Other'))
- cgd$D <- factor(cgd$D, levels=c(1,2), labels=c('Occurred', 'Censored'))
- #time-to-event setup
- survData <- Surv(time=cgd$T1, event=(cgd$D=='Occurred'))
- #overall
- summary(survData)
- km1 <- survfit(survData ~ 1, error="greenwood", conf.type="log-log", data=cgd)
- #png(file='01.overall.png', width=1200, height=960, res=144, bg='transparent')
- plot(km1, xlab="Survival time (days)", ylim=c(.25,1),
- ylab="Estimated survival probability",
- main="Figure 1. Overall Kaplan-Meier Estimator with pointwise CI", conf=T, mark=3, cex=0.5)
- #dev.off()
我现在知道把COX PH MODEL估计的概率画到KAPLAN-MEIER曲线上
- kmZ1 <- survfit(survData ~ Z1, data=cgd)
- plot(kmZ1, xlab="Time (days)", ylab="Estimated probability of no infection", ylim=c(0,1),
- main="Figure 2. Kaplan-Meier Estimator with pointwise CI by Treatment", lty=1:2, mark=3, cex=0.5, col=c("red","blue"))
- cox1 <- coxph(survData ~ Z1+Z2, data=cgd)
- lines(survexp(~Z1, ratetable=cox1, data=cgd), col=c('purple', 'orange'), lty=c(3,4))
请问
（1）如何把WEIBULLREGRESSION拟合的概率也画上去，我在网上找的代码不管用
SOLVED
- kmZ1 <- survfit(survData ~ Z1, data=cgd)
- plot(kmZ1, xlab="Time (days)", ylab="Estimated probability of no infection", ylim=c(0,1),
- main="Figure 2. Kaplan-Meier Estimator with pointwise CI by Treatment", lty=1:2, mark=3, cex=0.5, col=c("red","blue"))
- #overlay predicted survival rate by COX PH model
- cox1 <- coxph(survData ~ Z1+Z2, data=cgd)
- lines(survexp(~Z1, ratetable=cox1, data=cgd), col=c('purple', 'orange'), lty=c(3,4))
- #overlay predicated survival rate by Weibull regression
- sWei <- survreg(survData ~ Z1,dist='weibull',data=cgd)
- time <- predict(sWei, newdata=list(Z1=as.factor('Gamma Interferon')),type="quantile",p=seq(.01,.99,by=.01))
- lines(time, 1-seq(.01,.99,.01), lty=5)
- time <- predict(sWei, newdata=list(Z1=as.factor('Placebo')),type="quantile",p=seq(.01,.99,by=.01))
- lines(time, 1-seq(.01,.99,.01), lty=6)
（2）怎样确定WEIBULL和COXPH哪个更适合这个数据？（500点）
（3）WEIBULLREGRESSION怎么检查ASSUMPTION和DIAGNOSTICS？（414点）
|
熟悉论坛请点击新手指南
京公网安备 11010802022788号
论坛法律顾问：王进律师
知识产权保护声明
免责及隐私声明
