#  R code to fit the Cox model to   
#  enter data on disease free survival 
#  for bone marrow transplant patients

 
# Enter the data into a data frame.

mtb <- read.table("c:/stat565/data/BMT2.txt",
      header=F, col.names=c("patient", "group", "time", 
       "status", "agep", "aged", "wait", "fab", "mtx")) 
mtb$z1 <- mtb$agep-28
mtb$z2 <- mtb$aged-28
mtb$z3 <- mtb$z1*mtb$z2    
mtb$z4 <- 0
mtb$z4[mtb$group>1] <- 1
mtb$z5 <- 0
mtb$z5[mtb$group>2] <- 1
mtb$z4 <- mtb$z4-mtb$z5
mtb$z6 <- mtb$fab
mtb$z7 <- 12*(mtb$wait/365)
mtb$z8 <- mtb$mtx

mtb

library(survival)

mtbfit <- coxph(Surv(time, status) ~ z1+z2+z3+z4+z5+z6+z7+z8, 
                data = mtb, x=T, y=T, robust=T,
               method="efron", singular.ok=T)

summary(mtb) 

# Martingale residuals

mart.res <- resid(mtb)
par( lwd=2, mex=1, fin=c(6,6))
plot(mtb$z1, mart.res, xlab="Z1", ylab="Residuals",
   main="Martingale Residuals", cex=1.0, lwd=2)
lines(lowess(mtb$z1, mart.res, iter=0), lty=1 )


# Testing PH assumption 

mtbfit3 <- coxph(Surv(time,status) ~ z1+z2+z3+z4+z5+z6+z7+z8,
               , data=mtb)
zph.vet <- cox.zph(mtbfit3, transform="log")
# output from cox.zph gives 
# 1) pearson correlation between scaled schoenfeld residuals
#and g(t) for each covariate
# 2) test statistic for test
# 3) global test statistic over all covariates

# Plot Schoenfeld residuals against the 
# transformed time scale.

 plot(zph.vet[1])
 


# Compute dfbeta (score) residuals

dfbeta <- resid(mtbfit3, type="dfbeta")
n <- length(mtb$z1)

plot(mtb$patient, dfbeta[,1], ylab="Influence for Z1",
  xlab="subject")