## Donner
donner <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/donner.txt")
fit1 <- glm(Status~Age*Sex, donner, family=binomial)
Female <- 1*(donner$Sex=="Female")
fit2 <- glm(Status~Age*Female, donner, family=binomial)
db1 <- dfbeta(fit1)[,"Age"]
db2 <- dfbeta(fit2)[,"Age"]
db <- ifelse(donner$Sex=="Female", db1, db2)
df <- data.frame(donner, Leverage=hatvalues(fit1), Cook=cooks.distance(fit1), db=db, pi=fit1$fitted.values, d=rstudent(fit1))
df$Age <- jitter(donner$Age, amount=0.5)

## Simple Cook's D plot
plot(fit1, 4)

## VIF
require(car)
vif(fit1)
vif(lm(Died~Age*Sex, donner)) ## Not exactly the same
vif(lm(Died~Age*Sex, donner, weight=fit1$weights)) ## Same as vif(fit1)

## Slide 13
xyplot(Leverage~Age|Sex,df,type="h", ylim=c(-0.02,max(df$Leverage)+.02))

## Slide 14
xyplot(Cook~Age|Sex,df,type="h",lwd=2,group=Status,auto.key=list(columns=2,points=FALSE,lines=TRUE),ylab="Cook's distance")

## Slide 15
xyplot(db~Age|Sex, df, type="h", group=Status, auto.key=list(columns=2, points=FALSE, lines=TRUE), ylab=expression(Delta[beta]))

## Slide 16
col <- c("#FF4E37FF", "#008DFFFF")[donner$Status]
plot(fit1$fitted.values, rstudent(fit1), pch=19,cex=2*df$Leverage/max(df$Leverage),xlab=expression(pi),ylab=expression(d*"(Studentized deleted)"),col=col)
plot(fit1$fitted.values, rstudent(fit1)^2, pch=19,cex=2*df$Leverage/max(df$Leverage),xlab=expression(pi),ylab=expression(d^2*"(Studentized deleted)"),col=col)