donner <- read.delim("donner.txt")

## Slide 9
fit <- glm(Status~Age*Sex,donner,family="binomial")
summary(fit)

## Fit the model directly (Slides 11-12)
X <- model.matrix(~Age*Sex,donner)
y <- 1*(donner$Status=="Survived")
b <- rep(0,4)
converged <- function(new,old,eps=.0001)
  {
    return(max(abs(new-old)) < eps)
  }
repeat
  {
    old <- b
    ## Calculate z and W
    eta <- X%*%b
    pi <- exp(eta)/(1+exp(eta))
    W <- diag(as.numeric(pi*(1-pi)))
    z <- eta + solve(W)%*%(y-pi)
    ## Update b
    b <- solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%z
    print(b)
    ## Check for convergence
    if (converged(b,old)) break
  }
VarB <- solve(t(X)%*%W%*%X)
SE <- sqrt(diag(VarB))
z <- b/SE
p <- 2*pnorm(-abs(z))

## Estimate probabilities (Slides 14-17)
XX <- data.frame(Age=rep(c(20,40,60),2),Sex=c(rep("Male",3),rep("Female",3)))
predict(fit,newdata=XX)
predict(fit,newdata=XX,type="response")
predict(fit,newdata=XX,type="response",se.fit=TRUE)

## Slide 18
nn <- 51
xx <- seq(15,65,len=nn)
XX <- data.frame(Age=rep(x,2),Sex=c(rep("Male",nn),rep("Female",nn)))
eta <- predict(fit,newdata=XX)
pi <- predict(fit,newdata=XX,type="response")
df <- data.frame(pi=pi,eta=eta,XX)
xyplot(eta~Age,df,group=Sex,type="l",ylab="Linear predictor",auto.key=list(points=FALSE,lines=TRUE,columns=2))
xyplot(pi~Age,df,group=Sex,type="l",ylab="Pr(Surviving)",auto.key=list(points=FALSE,lines=TRUE,columns=2))

## Slide 20
fit <- lm(y~Age*Sex,donner)
eta <- predict(fit,newdata=XX)
pi <- predict(fit,newdata=XX,type="response")
df1 <- cbind(df,Method="Logistic")
df2 <- data.frame(pi=pi,eta=eta,XX,Method="Linear")
df <- rbind(df1,df2)
xyplot(pi~Age|Sex,df,group=Method,type="l",auto.key=list(points=FALSE,lines=TRUE,columns=2),ylab="Pr(Surviving)")

## Slide 22
fit <- glm(Status~Age*Sex + Sex*I(Age^2),donner,family="binomial")
eta <- predict(fit,newdata=XX)
pi <- predict(fit,newdata=XX,type="response")
df1$Method <- "Logistic (linear)"
df3 <- data.frame(pi=pi,eta=eta,XX,Method="Logistic (quadratic)")
df <- rbind(df1,df3)
xyplot(pi~Age|Sex,df,group=Method,type="l",auto.key=list(points=FALSE,lines=TRUE,columns=2),ylab="Pr(Surviving)")