Lab5

library(survival)

## Loading required package: splines

Myeloma=read.csv("Myeloma.csv")
head(Myeloma)

##   Time Status Age Sex Bun Ca   Hb Pcells Protein
## 1   13      1  66   1  25 10 14.6     18       1
## 2   52      0  66   1  13 11 12.0    100       0
## 3    6      1  53   0  15 13 11.4     33       1
## 4   40      1  69   1  10 10 10.2     30       1
## 5   10      1  65   1  20 10 13.2     66       0
## 6    7      0  57   0  12  8  9.9     45       0

1. To fit the Exponential AFT model:

survreg(Surv(Time, Status)~Bun, dist='exponential', data=Myeloma)

## Call:
## survreg(formula = Surv(Time, Status) ~ Bun, data = Myeloma, dist = "exponential")
## 
## Coefficients:
## (Intercept)         Bun 
##     3.90088    -0.01637 
## 
## Scale fixed at 1 
## 
## Loglik(model)= -155.8   Loglik(intercept only)= -159.8
##  Chisq= 7.95 on 1 degrees of freedom, p= 0.0048 
## n= 48

The time ratio for a 1-unit increase in Bun is -.0164.

To determine whether the model can also be interpreted as a proportional hazards model, create and evaluate the hazard function for each model.

m1=function(x) dexp(x,1/exp(3.90087890))/(1-pexp(x,1/exp(3.90087890)))
m2=function(x) dexp(x,1/exp(3.90087890-0.01636856))/(1-pexp(x,1/exp(3.90087890-0.01636856)))
m2(10)/m1(10)

## [1] 1.017

m2(25)/m1(25)

## [1] 1.017

exp(0.01636856) #should equal m2/m1 (when evaluated at the same point) for a PH model

## [1] 1.017

This can be evaluated as a proportional hazards model, with a hazard ratio of 1.017.

2. To fit the Weibull AFT model:

survreg(Surv(Time, Status)~Bun, dist='weibull', data=Myeloma)

## Call:
## survreg(formula = Surv(Time, Status) ~ Bun, data = Myeloma, dist = "weibull")
## 
## Coefficients:
## (Intercept)         Bun 
##      3.9095     -0.0164 
## 
## Scale= 0.9306 
## 
## Loglik(model)= -155.7   Loglik(intercept only)= -159.8
##  Chisq= 8.19 on 1 degrees of freedom, p= 0.0042 
## n= 48

To determine whether the Weibull model can be interpreted as a proportional hazards model:

m3 = function(x) dweibull(x,1/0.9306184,exp(3.90947497))/(1-pweibull(x,1/0.9306184,exp(3.90947497)))
m4 = function(x) dweibull(x,1/0.9306184,exp(3.90947497-0.01639564))/(1-pweibull(x,1/0.9306184,exp(3.90947497-0.01639564)))
m4(10)/m3(10)

## [1] 1.018

m4(25)/m3(25)

## [1] 1.018

These ratios are equal, so this can also be evaluated as a proportional hazards model with a hazard ratio of 1.018.

3. As the value of a covariate changes, the shape parameter of the Weibull changes.

4. The hazard ratio is \(b_{1}^{a}\)/\(b_{0}^{a}\), or \((b_{1}/b_{0})^{a}\).

5. a is the Scale output of the model, \(b_{0}\) is the Intercept term, and \(b_{1}\) is the intercept plus the coefficient for Bun.

6. The hazard ratio found in (2) equals the ratio of the value of the intercept plus the Bun coefficient to the intercept, all to the power of the scale value. This model uses the scale value (a) to create the ratio, whereas the Exponential AFT model only uses \(b_{0}\) and \(b_{1}\).

7. To fit the model and plot the hazard functions:

survreg(Surv(Time, Status)~Bun, dist='lognormal', data=Myeloma)

## Call:
## survreg(formula = Surv(Time, Status) ~ Bun, data = Myeloma, dist = "lognormal")
## 
## Coefficients:
## (Intercept)         Bun 
##     3.36851    -0.01481 
## 
## Scale= 1.146 
## 
## Loglik(model)= -155.4   Loglik(intercept only)= -159.6
##  Chisq= 8.57 on 1 degrees of freedom, p= 0.0034 
## n= 48

ln0 = function(x) dlnorm(x,3.3685146, 1.146415 )/(1-plnorm(x,3.3685146,1.146415))
ln1=function(x) dlnorm(x,3.3685146-0.0148148, 1.263288)/(1-plnorm(x,3.3685146-0.469791536,1.263288))
curve(ln1,0,100)
curve(ln0,add=T,lty=2)

After a point, both plots appear to be decreasing at a parallel rate, meaning that it is impossible for one line to always be proportional to the other.

Hazard ratio at time = 5 months and time = 20 months:

HR = function(x) ln1(x)/ln0(x)
HR(5)

## [1] 1.258

HR(20)

## [1] 1.229

8. To fit the Cox Proportional-Hazards model:

mod=coxph(Surv(Time, Status)~Bun, data=Myeloma)

The coefficient on the Bun term tells us that for this model, a 1-unit increase in Bun means an increase in survival time by exp(.0202) units.

9. The sign of the Bun coefficient in the Cox PH model is negative, as it was for the other models. This is okay because the parametric models are for survival time, so negative coefficients mean decreased hazard of death. This is also true for Cox PH models.

10. To plot the estimated survival curve from a Cox PH model:

plot(survfit(mod), mark.time=F) #draws mean survival curve by default

11. To graph the estimated survival curves for Bun=20 and Bun=50:

new = data.frame( Bun = c( 20 , 50 ) ) #create comparison of survival for two specific values of Bun
plot( survfit( mod , newdata = new ) , mark.time=F, lty = 1:2 ) #solid line (lty=1) for Bun=20, and dashed line (lty=2) for Bun=50

12. To compare the curves in (11) to the curves found using the log-normal model:

plot( survfit( mod , newdata = new ) , mark.time=F, lty = 1:2 )
curve(ln1, add=T, col='red')
curve(ln0,add=T, col='blue')

13. To create a Cox PH model for Time using only Age and Sex:

m=coxph(Surv(Time, Status)~Age+Sex, data=Myeloma)
new=data.frame(Sex=c(0,1), Age=rep(mean(Myeloma$Age), 2))
plot(survfit(m, newdata=new), mark.time=F, col=1:2)

The male-to-female hazard ratio is exp(-.085), or .919 (meaning that males have shorter survival times than females).

The hazard of a female at age x is the same as the hazard of a male at age (x+6.07).

How to calculate this:

\(e^{.014a}\)=\(e^{-.085+.014b}\) #since \(h_{0}\) is the same for both, these must be equal

.085=.014(b-a) #rearrange

6.071=b-a #simplify

14. A Cox PH model with no covariates reduces to a Kaplan-Meier curve because this type a model is only the baseline hazard, which we find nonparametrically.