Lab6
Chidinma Emenike
11/21/2014
library(survival)## Loading required package: splinesBurns=read.csv("Burns.csv")
head(Burns)## Time Status Treatment Gender Race PercBurn Head Buttock Trunk UpperLeg
## 1 12 0 0 0 0 15 0 0 1 1
## 2 9 0 0 0 1 20 0 0 1 0
## 3 7 1 0 0 1 15 0 0 0 1
## 4 29 0 0 0 0 20 1 0 1 0
## 5 4 1 0 0 1 70 1 1 1 1
## 6 8 1 0 0 1 20 1 0 1 0
## LowerLeg RespTract BurnType ExcisionTime ExcisionStatus AntibioticTime
## 1 0 0 2 12 0 12
## 2 0 0 4 9 0 9
## 3 1 0 2 13 0 13
## 4 0 0 2 11 1 29
## 5 0 0 2 28 1 31
## 6 0 0 4 11 0 11
## AntibioticStatus
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0coxph(Surv(Time, Status)~ExcisionStatus, data=Burns)## Call:
## coxph(formula = Surv(Time, Status) ~ ExcisionStatus, data = Burns)
##
##
## coef exp(coef) se(coef) z p
## ExcisionStatus -0.77 0.463 0.316 -2.44 0.015
##
## Likelihood ratio test=5.7 on 1 df, p=0.0169 n= 154, number of events= 48The p-value of ExcisionStatus is less than .05 (.015), indicating an association between removal of the burn wound and lengthened time to infection.
- To make the c-log-log plot for the model in (1):
km=survfit(Surv(Time, Status)~ExcisionStatus, data=Burns) #use KM, not Cox PH
plot(km,fun='cloglog',mark.time=F,col=1:2)
The curves in the graph are not linear, but since the data is not necessarily Exponential, and we’re testing for appropriateness of a Cox PH model, this doesn’t matter. All we’re looking for is a constant gap between the curves.
The gap between the curves appears to be about the same value across the x-axis, so the PH assumption is adequate.
- To perform a formal test of th PH assumption and obtain a p-value:
mod=coxph(Surv(Time, Status)~ExcisionStatus, data=Burns)
cox.zph(mod)## rho chisq p
## ExcisionStatus 0.138 0.852 0.356The p-value from this test is greater than .05 (.356), indicating that a Cox PH assumption is acceptable for this data (i.e., agreeing with our answer from 2).
ExcisionStatus can change over time (i.e. depends on time), so it should be treated as a time-varying covariate. We could do this by creating a new data set using the information in Burns. Our new data would be broken up by day and would have a binary variable representing whether nor not the wound was removed (0=no, 1=yes) for each patient. All the patients would have an individual identifying variable.
To see all relevant data points associated with one particular patient:
BurnsTVC=read.csv("BurnsTVC.csv")
BurnsTVC[29:57, c(1,2,3,21)]## start stop Status.time Excision
## 29 0 1 0 0
## 30 1 2 0 0
## 31 2 3 0 0
## 32 3 4 0 0
## 33 4 5 0 0
## 34 5 6 0 0
## 35 6 7 0 0
## 36 7 8 0 0
## 37 8 9 0 0
## 38 9 10 0 0
## 39 10 11 0 1
## 40 11 12 0 1
## 41 12 13 0 1
## 42 13 14 0 1
## 43 14 15 0 1
## 44 15 16 0 1
## 45 16 17 0 1
## 46 17 18 0 1
## 47 18 19 0 1
## 48 19 20 0 1
## 49 20 21 0 1
## 50 21 22 0 1
## 51 22 23 0 1
## 52 23 24 0 1
## 53 24 25 0 1
## 54 25 26 0 1
## 55 26 27 0 1
## 56 27 28 0 1
## 57 28 29 0 1Time: 29 (how many days patient was studied)
Status: 0 (status at end of study, this patient was censored)
ExcisionTime: 11 (at what point excision occured)
ExcisionStatus: 1 (did excision ever occur)
- To fit a Cox model which treats excision as a time-varying covariate:
coxph(Surv(start, stop, Status.time)~Excision+Excision:stop, data=BurnsTVC)## Call:
## coxph(formula = Surv(start, stop, Status.time) ~ Excision + Excision:stop,
## data = BurnsTVC)
##
##
## coef exp(coef) se(coef) z p
## Excision -1.9252 0.146 1.3015 -1.479 0.14
## Excision:stop 0.0843 1.088 0.0954 0.884 0.38
##
## Likelihood ratio test=5.23 on 2 df, p=0.0733 n= 3357, number of events= 48This p-values found here (.14 for Excision and .38 for Excision:stop) are both insignificant
- To fit a Cox model for Time by Treatment using both Burns and BurnsTVC:
coxph(Surv(Time, Status)~Treatment, data=Burns)## Call:
## coxph(formula = Surv(Time, Status) ~ Treatment, data = Burns)
##
##
## coef exp(coef) se(coef) z p
## Treatment -0.561 0.57 0.293 -1.91 0.056
##
## Likelihood ratio test=3.73 on 1 df, p=0.0535 n= 154, number of events= 48coxph(Surv(start, stop, Status.time)~Treatment, data=BurnsTVC)## Call:
## coxph(formula = Surv(start, stop, Status.time) ~ Treatment, data = BurnsTVC)
##
##
## coef exp(coef) se(coef) z p
## Treatment -0.561 0.57 0.293 -1.91 0.056
##
## Likelihood ratio test=3.73 on 1 df, p=0.0535 n= 3357, number of events= 48Treatment is a fixed (i.e. not a time-varying) covariate, so using BurnsTVC makes no difference in the output. The only difference is the amount of data. Since BurnsTVC breaks the data up by day, there are more values to consider (3357 total days compared to 154 people).
- To check whether the PH assumption is valid using the interaction between Treatment and Time in BurnsTVC:
coxph(Surv(start, stop, Status.time)~Treatment+Treatment:stop, data=BurnsTVC)## Call:
## coxph(formula = Surv(start, stop, Status.time) ~ Treatment +
## Treatment:stop, data = BurnsTVC)
##
##
## coef exp(coef) se(coef) z p
## Treatment -0.3791 0.684 0.4192 -0.904 0.37
## Treatment:stop -0.0157 0.984 0.0264 -0.594 0.55
##
## Likelihood ratio test=4.1 on 2 df, p=0.129 n= 3357, number of events= 48The coefficient of the Treatment term (-.3791) indicates that the Treatment=1 group has a lower time to infection. The coefficient of the Treatment:stop term (-.0157) indicates that when stop=k, log of the hazard ratio is -.3791-.0157k.
The estimated HR is 0.585 at stop = 10 and 0.500 at stop = 20.
How do we know?
exp(-0.3791-0.0157*10) #HR at stop=10## [1] 0.585exp(-0.3791-0.0157*20) #HR at stop=20## [1] 0.5\(HR_{stop10}\) \(\neq\) \(HR_{stop20}\), so the PH assumption is not appropriate.
- To add the coefficients from 8 to the Schoenfeld residuals:
m=coxph(Surv(Time, Status)~Treatment, data=Burns)
plot(cox.zph(m, transform=identity))
abline(h=-0.77, col='red') #coeffiecient from original PH model
abline(-.3791, -.0157, col='blue') #coefficients from the model in (8)