whas.table <- read.table('http://www.stanford.edu/class/stats262/data/whas.csv', sep=',', header=T)
whas.table <- whas.table[(whas.table$yrgrp == '1975 & 1978'),]
attach(whas.table)

library(survival)

# Using AIC 

for (var in c('age', 'sex', 'cpk', 'sho', 'chf', 'lenstay')) {
  print(eval(parse(text=paste("coxph(Surv(lenfol, (fstat == 'Dead')) ~", var, ')'))))
}

# Variables 'age', 'sho', 'chf', 'lenstay' make the "cut"

whas.initial.coxph <- coxph(Surv(lenfol, (fstat == 'Dead')) ~ age + sho + chf + lenstay)
print(whas.initial.coxph)

postscript('lenstay_smooth.ps', height=4, width=5, paper='special', horizontal=F)
c.smooth <- lowess(lenstay, (fstat == 'Dead'))$y
h.smooth <- lowess(lenstay, (fstat == 'Dead') - residuals.coxph(whas.initial.coxph, type='mart'))$y

plot(lenstay, log(c.smooth / h.smooth), xlab='Length of stay', ylab='Smoothed log(c/H)')
lines(lowess(lenstay, log(c.smooth / h.smooth)))
dev.off()

postscript('lenstay_mart.ps', height=4, width=5, paper='special', horizontal=F)
mart.resid <- residuals.coxph(coxph(Surv(lenfol, (fstat == 'Dead')) ~ 1), type='mart')
plot(lenstay, mart.resid, xlab='Length of stay', ylab='Martingale residual')
lines(lowess(lenstay, mart.resid))
dev.off()

# Suggests quadratic effect for lenstay

postscript('age_smooth.ps', height=4, width=5, paper='special', horizontal=F)
c.smooth <- lowess(age, (fstat == 'Dead'))$y
h.smooth <- lowess(age, (fstat == 'Dead') - residuals.coxph(whas.initial.coxph, type='mart'))$y

plot(age, log(c.smooth / h.smooth), xlab='Age', ylab='Smoothed log(c/H)')
lines(lowess(age, log(c.smooth / h.smooth)))
dev.off()

postscript('age_mart.ps', height=4, width=5, paper='special', horizontal=F)
mart.resid <- residuals.coxph(whas.initial.coxph, type='mart')
plot(age, mart.resid, xlab='Length of stay', ylab='Martingale residual')
lines(lowess(age, mart.resid))
dev.off()

# Look at plot for CPK just to make sure we didn't miss anything

postscript('cpk_smooth.ps', height=4, width=5, paper='special', horizontal=F)
c.smooth <- lowess(cpk, (fstat == 'Dead'))$y
h.smooth <- lowess(cpk, (fstat == 'Dead') - residuals.coxph(whas.initial.coxph, type='mart'))$y

plot(cpk, log(c.smooth / h.smooth), xlab='CPK', ylab='Smoothed log(c/H)')
lines(lowess(cpk, log(c.smooth / h.smooth)))
dev.off()

postscript('cpk_mart.ps', height=4, width=5, paper='special', horizontal=F)
mart.resid <- residuals.coxph(whas.initial.coxph, type='mart')
plot(cpk, mart.resid, xlab='Length of stay', ylab='Martingale residual')
lines(lowess(cpk, mart.resid))
dev.off()

# Suggests linear effect for lenstay is adequate

whas.final.coxph <- coxph(Surv(lenfol, (fstat == 'Dead')) ~ sho + age + poly(lenstay, 2) + chf)
print(whas.final.coxph)