################## ## Challenger e.g. ################## shuttle <- read.table("Data Sets/shuttle.txt", header = TRUE) attach(shuttle) ## Teleconference plot plot(jitter(damage, 1) ~ tempF, data = subset(shuttle, damage == 1), ylim = c(0, 1.2)) ## Better plot! plot(jitter(damage, 0.1) ~ tempF) ## With sectioned proportions plot(jitter(damage, 0.1) ~ tempF) sections <- cut(tempF, c(50, 55, 60, 65, 70, 75, 80, 85)) meanD <- tapply(damage, sections, mean) low <- c(50, 55, 60, 65, 70, 75, 80) high <- c(55, 60, 65, 70, 75, 80,85) segments(low, meanD, high, meanD) ## Logistic model shu.glm <- glm(damage ~ tempF, data = shuttle, family = binomial) temps <- seq(50, 85, 0.5) lines(temps, predict(shu.glm, type = "response", newdata = data.frame(tempF = temps))) ## Group data launch <- gl(138/6, 6) s <- tapply(damage, launch, sum) f <- 6 - s temp <- tapply(tempF, launch, min) grp.glm <- glm(cbind(s, f) ~ temp, data = shuttle, family = binomial) plot(jitter(c(s/f), 0.5) ~ temp, ylim = c(0, 1.2)) temps <- seq(50, 85, 0.5) lines(temps, predict(shu.glm, type = "response", newdata = data.frame(tempF = temps))) ########################### ### Car purchase e.g. ########################## dat <- read.table("car_income.txt", header = TRUE) attach(dat) boxplot(income ~ purchase) plot(income ~ age) glm1 <- glm(purchase ~ income + age, family="binomial") glm2 <- glm(purchase ~ income + as.factor(age), family="binomial") old <- as.numeric(age > 3) glm3 <- glm(purchase ~ income + old, family="binomial") summary(glm3) anova(glm3, test = "Chisq")