### script for MT05 week 6 practical

library(MASS)
options(contrasts = c("contr.treatment", "contr.poly"))


## 1 Low Birth Weights

?birthwt
attach(birthwt)
race <- factor(race, labels = c("white", "black", "other"))
table(ptl)
ptd <- factor(ptl > 0)
ftv <- factor(ftv)
table(ftv)
levels(ftv)[-(1:2)] <- "2+"
table(ftv)  # as a check
bwt <- data.frame(low = factor(low), age, lwt, race,
    smoke = (smoke > 0), ptd, ht = (ht > 0), ui = (ui > 0), ftv)
detach(); rm(race, ptd, ftv)

birthwt.glm <- glm(low ~ ., family = binomial, data = bwt)
summary(birthwt.glm, cor = F)

birthwt.step <- stepAIC(birthwt.glm)
birthwt.step$anova
birthwt.step2 <- stepAIC(birthwt.glm, ~ .^2 + I(scale(age)^2)
                         + I(scale(lwt)^2), trace = F)
birthwt.step2$anova
summary(birthwt.step2, cor = F)$coef
table(bwt$low, predict(birthwt.step2) > 0)
plot(birthwt.step2)


## 2 Speed Limits in Sweden

try(rm(y)) # to be careful
attach(Traffic)
tr <- data.frame(y, day = factor(day), year = factor(year), limit)
detach()
fm <- glm(y ~ day + year + limit, data = tr, family = poisson)
fm2 <- stepAIC(fm, scope = list(upper = . ~ . + limit:year))
plot(fm2)
exp(coef(fm2)["limit"])


## 3 Cancer Deaths of Atomic Bomb Survivors

options(digits=5)
fm <- glm(Deaths ~ tyears + dose + offset(log(AtRisk)), family = poisson,
          data = Atomic)

dropterm(fm, test="Chisq")

Atomic$ot <- c(4,10,14,18,22,26,30)[Atomic$tyears]
fm2 <- glm(Deaths ~ ot + dose + offset(log(AtRisk)), family = poisson,
        data = Atomic)
anova(fm, fm2, test="Chisq")

preddata <- Atomic
preddata$AtRisk <- rep(1,42)
Atomic$lrates <- predict(fm2, preddata)
tab <- matrix(exp(-Atomic$lrates),6,7, byrow=T)
dimnames(tab) <- list(levels(Atomic$dose), levels(Atomic$tyears))
round(tab,1)

trellis.device()
xyplot(exp(-lrates) ~ ot | dose, data=Atomic, type="b")
Atomic$selrates <- predict(fm2, preddata, se=T)$se
Atomic$up <- Atomic$lrates + 2*Atomic$selrates
Atomic$low <- Atomic$lrates - 2*Atomic$selrates

xyplot(exp(lrates) ~ ot | dose, data=Atomic, type="b", subscripts=T,
  panel = function(x,y,subscripts, ...) {
     panel.xyplot(x,y, ...)
     error.bar(x, y, exp(Atomic$low[subscripts]),
               exp(Atomic$up[subscripts]), add=T, incr=F)
  }
)

xyplot(lrates ~ as.numeric(as.character(dose)) | tyears, data=Atomic,
       type="b", subscripts=T,
  panel = function(x,y,subscripts, ...) {
     panel.xyplot(x,y, ...)
     error.bar(x, y, Atomic$low[subscripts],
               Atomic$up[subscripts], add=T, incr=F)
  }
)