## AIDS 
aids <- read.delim("../../data/aids.txt")
attach(aids)

plot(Year,Cases,pch=19,ylim=c(0,300))

fit <- glm(Cases~Year,data=aids,family=poisson)
summary(fit)
confint(fit)
fit0 <- glm(Cases~1,data=aids,family=poisson)
anova(fit0,fit,test="Chisq")

xx <- seq(min(Year),max(Year),length=51)
pred <- predict(fit,newdata=data.frame(Year=xx),se.fit=TRUE)
L <- pred$fit-qnorm(.975)*pred$se.fit
U <- pred$fit+qnorm(.975)*pred$se.fit
par(mfrow=c(1,2))
plot(range(Year),range(c(L,U)),type="n",xlab="Year",ylab="Linear predictor")
polygon(c(xx,rev(xx)),c(L,rev(U)),border=FALSE,col="gray")
lines(xx,pred$fit,lwd=2)

plot(Year,Cases,type="n",ylim=c(0,350))
polygon(c(xx,rev(xx)),exp(c(L,rev(U))),border=FALSE,col="gray")
lines(xx,exp(pred$fit),lwd=2)
points(Year,Cases,pch=19)

## Diagnostics
df <- data.frame(aids,Leverage=hatvalues(fit),Cook=cooks.distance(fit),db=dfbeta(fit)[,"Year"],pi=fit$fitted.values,d=rstudent(fit))
xyplot(Leverage~Year,df,type="h",lwd=2,ylim=c(-0.02,max(df$Leverage)+.02))
xyplot(Cook~Year,df,type="h",lwd=2,ylab="Cook's distance")
xyplot(db~Year,df,type="h",lwd=2,ylab=expression(Delta[beta]))
plot(fit$fitted.values,rstudent(fit),pch=19,cex=3*df$Leverage/max(df$Leverage),xlab=expression(lambda),ylab=expression(d*"(Studentized deleted)"))
plot(fit$fitted.values,rstudent(fit)^2,pch=19,cex=3*df$Leverage/max(df$Leverage),xlab=expression(lambda),ylab=expression(d^2*"(Studentized deleted)"))
abline(h=qchisq(.975,1),col="red")

## Predictive power
1-crossprod(fit$y-fit$fitted)/crossprod(fit$y-mean(fit$y))
1-crossprod(fit$y-fit$fitted)/fit$df.residual/(crossprod(fit$y-mean(fit$y))/fit$df.null)
1-fit$deviance/fit$null.deviance
1-fit$deviance/fit$df.residual/(fit$null.deviance/fit$df.null)

fit <- glm(Cases~Year,data=aids,family=poisson)
fit2 <- glm(Cases~Year,data=aids,family=poisson,subset=(Year <= 1991))
fit3 <- glm(Cases~Year+I(Year^2),data=aids,family=poisson)

## British doctors
britdoc <- read.delim("../../data/britdoc.txt")
PY <- britdoc$PersonYears/1000

fit <- glm(Deaths~Age+Smoking,britdoc,offset=log(PY),family=poisson)
df <- data.frame(britdoc,PY=1)
rates <- cbind(df[,1:2],rate=round(predict(fit,df,type="response"),2))
xtable(R <- cbind(rates[1:5,"rate"],rates[6:10,"rate"]))

exp(coef(fit)["Age65 to 74"] - coef(fit)["Age45 to 54"])