6. Choosing between statistical models

# GUSTO-I age effect modeling, August 2018
library(rms)
library(foreign)
# Import gusto; not publicly available; could use a smaller part, as is made avaliable
gusto <- as.data.frame(read.spss('gusto.sav'))

gusto$AGE <- gusto$AGE/10
gusto     <- gusto[,c("AGE","DAY30")]

## Age linear; add square; rcs
agegusto.linear <- lrm(DAY30~AGE,	data=gusto, x=T,y=T)
agegusto.square <- lrm(DAY30~pol(AGE,2),	data=gusto, x=T,y=T)
agegusto.rcs    <- lrm(DAY30~rcs(AGE,5),	data=gusto, x=T,y=T)

## linear spline, age 50 (equals value 5.0 for decades)
library(lspline)
agegusto.linearspline    <- lrm(DAY30~lspline(AGE,5.0),	data=gusto, x=T,y=T)

## dichotomize
agegusto.cat65    <- lrm(DAY30~ifelse(AGE<6.5,0,1),	data=gusto, x=T,y=T)


######################################################
## Plots; first at lp scale, then probability scale ##
par(mfrow=c(1,2))

newdata.age =data.frame("AGE"=seq(2,9.5, by=0.01))
pred.agegusto.linear  <- predict(agegusto.linear, newdata.age)
pred.agegusto.square  <- predict(agegusto.square,newdata.age)
pred.agegusto.rcs     <- predict(agegusto.rcs,newdata.age)
pred.agegusto.linearspline <- predict(agegusto.linearspline,newdata.age)
pred.agegusto.cat65   <- predict(agegusto.cat65,newdata.age)

plot(x=10*newdata.age[,1], y=pred.agegusto.linear, xlim=c(9,95), ylim=c(-6,0), las=1,
     xlab='Age in years', ylab='logit of 30-day mortality', cex.lab=1.2, type="l", lwd=2)
lines(x=10*newdata.age[,1], y=pred.agegusto.square, lty=2)
lines(x=10*newdata.age[,1], y=pred.agegusto.rcs, lty=3)
lines(x=10*newdata.age[,1], y=pred.agegusto.linearspline, lty=4)
lines(x=10*newdata.age[,1], y=pred.agegusto.cat65, lty=5)

legend("topleft", legend=c("linear","square","rcs","lspline","dichotomize at 65 years"),lty=1:5)

scat1d(x=10*gusto$AGE, side=1)

##################################
## Repeat for probability scale ##
newdata.age =data.frame("AGE"=seq(2,9.5, by=0.01))
pred.agegusto.linear  <- plogis(predict(agegusto.linear, newdata.age))
pred.agegusto.square  <- plogis(predict(agegusto.square,newdata.age))
pred.agegusto.rcs     <- plogis(predict(agegusto.rcs,newdata.age))
pred.agegusto.linearspline <- plogis(predict(agegusto.linearspline,newdata.age))
pred.agegusto.cat65   <- plogis(predict(agegusto.cat65,newdata.age))

plot(x=10*newdata.age[,1], y=pred.agegusto.linear, xlim=c(9,95), ylim=c(0,.5), las=1,
     xlab='Age in years', ylab='probability of 30-day mortality', cex.lab=1.2, type="l", lwd=2)
lines(x=10*newdata.age[,1], y=pred.agegusto.square, lty=2)
lines(x=10*newdata.age[,1], y=pred.agegusto.rcs, lty=3)
lines(x=10*newdata.age[,1], y=pred.agegusto.linearspline, lty=4)
lines(x=10*newdata.age[,1], y=pred.agegusto.cat65, lty=5)
scat1d(x=10*gusto$AGE, side=1)