library(MASS)
?Insurance
options(contrasts=c("contr.treatment", "contr.poly"))
fm <- glm(Claims ~ District*Group*Age + offset(log(Holders)),
          data = Insurance, family = poisson)
residuals(fm)
fm <- glm(Claims ~ District*Group*Age + offset(log(Holders)),
          data = Insurance, family = poisson, maxit=50, trace=T)
fm <- glm(Claims ~ (District+Group+Age)^2 + offset(log(Holders)),
          data = Insurance, family = poisson)
summary(fm, cor=F)

fm2 <- stepAIC(fm, direction="both")

dropterm(fm2, test="Chisq")

summary(fm2, cor=F)

fm3 <- update(fm2, . ~ District + as.integer(Group) + as.integer(Age) + offset(log(Holders)))
summary(fm3, cor=F)

fm4 <- update(fm2, . ~ District + C(Group, poly, 1) + C(Age, poly, 1) + offset(log(Holders)))
summary(fm4, cor=F)

plot(profile(fm4))

plot(fm4) # diagnostic plots
(dfm <- cbind(Insurance, fitted = round(fitted(fm4)),
              res=round(residuals(fm4, type="pearson"), 1)))

tab <- array(dfm$Claims/dfm$Holders, dim=c(4,4,4))
dimnames(tab) <- list(Age=levels(dfm$Age), Group=levels(dfm$Group), District=levels(dfm$district))
mosaicplot(aperm(tab, 3:1), color=T)

tab2 <- array(fitted(fm4)/dfm$Holders, dim=c(4,4,4))
dimnames(tab2) <- list(Age=levels(dfm$Age), Group=levels(dfm$Group), District=levels(dfm$district))

par(mfrow=c(1,2), pty="s")
mosaicplot(aperm(tab, 3:1), color=T, main="Observed")
mosaicplot(aperm(tab2, 3:1), color=T, main="Expected")

par(mfrow=c(1,3), pty="m")
nd <- Insurance[1:4,]
nd$Holders <- rep(1000,4)
pr <- predict(fm4, nd, se=T, type="response")
error.bar(1:4, pr$fit, -2*pr$se, +2*pr$se, xaxt="n", xlab="Age",
          ylab="claims/1000", ylim=c(0,350))
axis(1, at=1:4, labels=levels(Insurance$Age))

nd <- Insurance[seq(1, by=4, len=4),]
nd$Holders <- rep(1000,4)
pr <- predict(fm4, nd, se=T, type="response")
error.bar(1:4, pr$fit, -2*pr$se, +2*pr$se, xaxt="n", xlab="Group",
          ylab="claims/1000", ylim=c(0,350))
axis(1, at=1:4, labels=levels(Insurance$Group))

nd <- Insurance[seq(1, by=16, len=4),]
nd$Holders <- rep(1000,4)
pr <- predict(fm4, nd, se=T, type="response")
error.bar(1:4, pr$fit, -2*pr$se, +2*pr$se, xaxt="n", xlab="District",
          ylab="claims/1000", ylim=c(0,350))
axis(1, at=1:4, labels=levels(Insurance$District))

par(mfrow=c(1,1), pty="m")

trellis.device()
nd <- Insurance
nd$Holders <- rep(1000,64)
pr <- predict(fm4, nd, se=T, type="response")
barchart(Age ~ pr$fit | Group*District, data=nd, xlab="rate", xlim=c(0,350))
# or better with a log axis:
barchart(Age ~ pr$fit | Group*District, data=nd, xlab="rate",
         scale=list(x=list(log=10)))