# 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)