#####################
## British doctors ##
#####################
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/britdoc.txt")
PY <- Data$PersonYears/1000

## No interaction
fit <- glm(Deaths~Age+Smoking, Data, offset=log(PY), family=poisson)
New <- data.frame(Data, PY=1)
p <- predict(fit, New, type="response")
matrix(p, 5, 2, dimnames=list(levels(Data$Age), levels(Data$Smoking)))

## Interaction
fit2 <- glm(Deaths~Age*Smoking, Data, offset=log(PY), family=poisson)
anova(fit2, fit, test="Chisq")
AIC(fit)
AIC(fit2)
p <- predict(fit2, New, type="response")
R <- matrix(p, 5, 2, dimnames=list(levels(Data$Age), levels(Data$Smoking)))
R <- cbind(R, RR=R[,1]/R[,2])
R

####################
## Overdispersion ##
####################
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/aids.txt")
attach(Data)
require(visreg)
require(MASS)

## Quasipoisson
fit <- glm(Cases~Year, Data, family=poisson)
visreg(fit, scale="response", ylim=c(0, 350), partial=FALSE, ylab="Cases")
points(Year,Cases,pch=19, cex=0.7)

fit.qp <- glm(Cases~Year, Data, family=quasipoisson)
visreg(fit.qp, scale="response", ylim=c(0, 350), partial=FALSE, ylab="Cases")
points(Year,Cases,pch=19, cex=0.7)

## Negative binomial: Mean-variance relationship
fit.nb <- glm.nb(Cases~Year, Data)
l <- seq(0,300,len=101)
plot(l,l,type="l",xlab=expression(lambda),ylab="Variance",lwd=3,col="gray",ylim=range(l+l^2/fit.nb$theta), las=1)
lines(l,l+l^2/fit.nb$theta,lwd=3,col="red")
legend("topleft", legend=c("Poisson","Negative Binomial"),lwd=3,col=c("gray","red"), ncol=2)

fit <- glm(Cases~Year, Data, family=poisson)
visreg(fit, scale="response", ylim=c(0, 350), partial=FALSE, ylab="Cases", main="Poisson")
points(Year,Cases,pch=19, cex=0.7)

visreg(fit.qp, scale="response", ylim=c(0, 350), partial=FALSE, ylab="Cases", main="Quasipoisson")
points(Year,Cases,pch=19, cex=0.7)

visreg(fit.nb, scale="response", ylim=c(0, 350), partial=FALSE, ylab="Cases", main="Negative binomial")
points(Year,Cases,pch=19, cex=0.7)