Survival Analysis Part 4: The Proportional Hazards Assumption and Leverage (Based on the book by Hosmer, Lemeshow, May)

rm(list = ls())
setwd("C:\\ASA")
# import data
whas500 = read.table("C:\\ASA\\whas500.DAT", sep = "", header = F)
names(whas500) = c("id",
                   "age",
                   "gender",
                   "hr",
                   "sysbp",
                   "diasbp",
                   "bmi",
                   "cvd",
                   "afb",
                   "sho",
                   "chf",
                   "av3",
                   "miord",
                   "mitype",
                   "year",
                   "admitdate",
                   "disdate",
                   "fdate",
                   "los",
                   "dstat",
                   "lenfol",
                   "fstat")
whas500$bmifp1 = (whas500$bmi/10)^2
whas500$bmifp2 = (whas500$bmi/10)^3
library(survival)
model.coxph = coxph( Surv(lenfol, fstat)~ bmifp1
                      + bmifp2
                      + age
                      + hr
                      + diasbp
                      + gender
                      + chf
                      + age*gender,
                      data = whas500,
                      method = "breslow")
summary(model.coxph)
# Table 6.1
zph.1 = cox.zph(model.coxph, transform = "identity")
zph.1
zph.log = cox.zph(model.coxph, transform = "log")
zph.log
zph.km = cox.zph(model.coxph, transform = "km")
zph.km
zph.rank = cox.zph(model.coxph, transform = "rank")
zph.rank
# Figure 6.1
plot(zph.log[4])
abline(h=0, lty=3)
plot(zph.1[4])
abline(h=0, lty=3)
plot(zph.km[4])
abline(h=0, lty=3)
plot(zph.rank[4])
abline(h=0, lty=3)

# Leverage: Figure 6.4 and Figure 6.5
score = residuals(model.coxph, type="score")
plot(whas500$bmi, score[,1], ylab="Score Residuals", xlab="bmifp1")
plot(whas500$hr, score[,4], ylab="Score Residuals", xlab="heart rate")
plot(whas500$age, score[,3], ylab="Score Residuals", xlab="age")