Survival Analysis Part 3: Quartile Design Variable and Fractional Polynomials (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")
library(survival)
model1.coxph = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + bmi + gender + chf,
                      data = whas500, method = "breslow")
# Table 5.6
summary(model1.coxph)
summary(whas500$bmi)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 13.05   23.22   25.95   26.61   29.39   44.84
whas500$bmicat = cut(whas500$bmi, c(13, 23.22, 25.95, 29.39, 44.84))
summary(whas500$bmicat)
# (13,23.2] (23.2,25.9] (25.9,29.4] (29.4,44.8]
# 125         126         124         125
model2.coxph = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + bmicat + gender + chf,
                      data = whas500, method = "breslow")
summary(model2.coxph)
bmi.coeff = data.frame(model2.coxph$coefficients)[4:6,]
(summary(whas500$bmi)[[1]] + summary(whas500$bmi)[[2]])/2
midpt = c((summary(whas500$bmi)[[1]] + summary(whas500$bmi)[[2]])/2,
          (summary(whas500$bmi)[[2]] + summary(whas500$bmi)[[3]])/2,
          (summary(whas500$bmi)[[3]] + summary(whas500$bmi)[[5]])/2,
          (summary(whas500$bmi)[[5]] + summary(whas500$bmi)[[6]])/2
          )
bmi.plot = cbind(midpt, rbind(0,data.frame(bmi.coeff)))
bmi.plot
# Figure 5.1
plot(bmi.plot[,1], bmi.plot[,2], ylab = "Log Hazard", xlab = "BMI", type = "b")
# Table 5.8
library(mfp)
model.null = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf,
                    data = whas500, method = "breslow")
-2*model.null$loglik
model.lin = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf + bmi,
                   data = whas500, method = "breslow")
-2*model.lin$loglik
# J = 1, df = 2
model.J1 = mfp( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf + fp(bmi, df = 2),
                  data = whas500, method = "breslow", family = cox)
-2*model.J1$loglik
model.J1$powers
# J = 2, df = 4
model.J2 = mfp( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf + fp(bmi, df = 4),
                data = whas500, method = "breslow", family = cox)
-2*model.J2$loglik
model.J2$powers
# Figure 5.2
whas500$resid = residuals(model.null, type="martingale", data=whas500)
plot(whas500$bmi, whas500$resid, xlab="BMI", ylab="Martingale Residuals")
lines(lowess(whas500$bmi, whas500$resid))