data(ICU)
# remove redundant variables (race, coma)
ICU1 <- ICU[,-c(4,20)]
# fit full model
icu.full <- glm(died ~ ., data=ICU1, family=binomial)
summary(icu.full)
# simpler model (found from a "best" subsets procedure)
icu.mod1 <- glm(died ~ age + sex + cancer + systolic + admit + uncons,
data=ICU1,
family=binomial)
summary(icu.mod1)
# even simpler model
icu.mod2 <- glm(died ~ age + cancer + admit + uncons,
data=ICU1,
family=binomial)
summary(icu.mod2)
anova(icu.mod2, icu.mod1, icu.full, test="Chisq")
## Reproduce Fig 6.12 from VCD
icu.fit <- data.frame(ICU, prob=predict(icu.mod2, type="response"))
# combine categorical risk factors to a single string
risks <- ICU[, c("cancer", "admit", "uncons")]
risks[,1] <- ifelse(risks[,1]=="Yes", "Cancer", "")
risks[,2] <- ifelse(risks[,2]=="Emergency", "Emerg", "")
risks[,3] <- ifelse(risks[,3]=="Yes", "Uncons", "")
risks <- apply(risks, 1, paste, collapse="")
risks[risks==""] <- "(none)"
icu.fit$risks <- risks
library(ggplot2)
ggplot(icu.fit, aes(x=age, y=prob, color=risks)) +
geom_point(size=2) +
geom_line(size=1.25, alpha=0.5) +
theme_bw() + ylab("Probability of death")
Run the code above in your browser using DataLab