####### Model Checking for Cox PH Models #######

library(survival)
library(MASS)

data(VA)
attach(VA)
# VA data 
VA <- data.frame(stime,status,treat=factor(treat),age,Karn,diag.time,
                 cell=factor(cell), prior=factor(prior))
va.cox1 <- coxph(Surv(stime, status) ~ Karn + cell, data=VA)
# Martingale Residual Plots (with a smoothed fit) 
par(mfrow=c(2,1))
scatter.smooth(VA$Karn, residuals(va.cox1), main="Residual Plot", 
               ylab="martingale residuals", xlab="Karnofsky score") 
scatter.smooth(VA$age, residuals(va.cox1), main="Residual Plot", 
               ylab="martingale residuals", xlab="age")
# Check for the constancy of the coefficients: Schoenfeld residual against 
# a nonlinear scale of time, with spline smoothers. 
par(mfrow=c(2,2))
plot(cox.zph(va.cox1))

# Motors data 
data(motors)
attach(motors)
motor.cox <- coxph(Surv(time, cens) ~ temp, motors)
scatter.smooth(motors$temp, residuals(motor.cox), main="Residual Plot", 
               ylab="martingale residuals", xlab="temperature")
plot(cox.zph(motor.cox))