## Data
donner <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/donner.txt")

## Basic usage
fit <- glm(Status~Age*Sex,donner,family=binomial)
summary(fit)

## With female as reference group
Female <- 1*(donner$Sex=="Female")
fit <- glm(Status~Age*Female,donner,family=binomial)
summary(fit)

## Construct by hand (slides 11-12)
X <- model.matrix(~Age*Sex,donner)
y <- 1*(donner$Status=="Survived")
b <- rep(0,4)
converged <- function(new, old, eps=.0001) {
  max(abs((new-old)/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 (15-18)
fit <- glm(Status~Age*Sex, donner, family=binomial)
New <- data.frame(Age=rep(c(20,40,60),2), Sex=c(rep("Male",3),rep("Female",3)))
predict(fit, newdata=New, type="response")
predict(fit, newdata=New, type="response", se.fit=TRUE)
pred <- predict(fit, newdata=New, se.fit=TRUE)
CI <- cbind(pred$fit - qnorm(.975)*pred$se, pred$fit + qnorm(.975)*pred$se)
exp(CI)/(1+exp(CI))

## Visreg
visreg(fit, "Age", by="Sex")
visreg(fit, "Age", by="Sex", scale="response")

## Slide 19
fit <- glm(Status~Age*Sex, donner, family=binomial)
xx <- seq(15, 65, len=99)
XX <- rbind(data.frame(Age=xx, Sex="Male"), data.frame(Age=xx, Sex="Female"))
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))

## Linear regression vs. logistic (slide 21)
fit <- lm(y~Age*Sex, donner)
eta <- predict(fit, newdata=XX)
pi <- predict(fit, newdata=XX, type="response")
df1 <- data.frame(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)")

## Quadratic
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)")