########## R script: ICUstepGAM ##########

library(gam)

ICU <- read.table("ICU.txt",header=TRUE)
raceBlack <- as.numeric(ICU$race==2)
raceOther <- as.numeric(ICU$race==3)
deepStupor <- as.numeric(ICU$LOC==2)
coma       <- as.numeric(ICU$LOC==3)
myICU <- ICU[,-c(4,20)]
myICU <- data.frame(myICU,raceBlack,raceOther,deepStupor,coma)

gam.object <-  gam(died~ age+female+surgical+cancer+chRenFail
                         +infection+CPR+SBP+heartRate+prevAdmin
                         +emergency+fracture+hiPO2+hiPH+hiPCO2
                         +hiBIC+hiCRE+raceBlack+raceOther+deepStupor
                         +coma,data=myICU,family=binomial)

stepFit <- step.gam(gam.object, 
               scope=list("age"=~1+age+s(age,3)+s(age,6)+s(age,9),
                          "female"=~1+female,
                          "surgical"=~1+surgical,
                          "cancer"=~1+cancer,
                          "chRenFail"=~1+chRenFail,
                          "infection"=~1+infection,
                          "CPR"=~1+CPR,
                          "SBP"=~1+SBP+s(SBP,3)+s(SBP,6)+s(SBP,9),
                          "heartRate"=~1+heartRate+s(heartRate,3)
                                       +s(heartRate,6)+s(heartRate,9),
                          "prevAdmin"=~1+prevAdmin,
                          "emergency"=~1+emergency,
                          "fracture"=~1+fracture,
                          "hiPO2"=~1+hiPO2,
                          "hiPH"=~1+hiPH,
                          "hiPCO2"=~1+hiPCO2,
                          "hiBIC"=~1+hiBIC,
                          "hiCRE"=~1+hiCRE,
                          "raceBlack"=~1+raceBlack,
                          "raceOther"=~1+raceOther,
                          "deepStupor"=~1+deepStupor,
                          "coma"=~1+coma))

print(summary(stepFit))

par(mfrow=c(2,5))
plot(stepFit,se=TRUE)

########## End of ICUstepGAM ##########