## Read in data; fit model
donner <- read.delim("../../data/donner.txt")
Female <- 1*(donner$Sex=="Female")

## Fit model
fit <- glm(Status~Age*Female,donner,family="binomial")

## Inference for beta
summary(fit)
confint(fit)
cbind(coef(fit)-qnorm(.975)*sqrt(diag(vcov(fit))),
      coef(fit)+qnorm(.975)*sqrt(diag(vcov(fit))))

## LRTs
fit0 <- glm(Status~Female+Age,donner,family="binomial")
anova(fit0,fit,test="Chisq") ## Test of interaction
fit0 <- glm(Status~Female+Age:Female,donner,family="binomial")
anova(fit0,fit,test="Chisq") ## Test of age effect in males
fit0 <- glm(Status~Age+Age:Female,donner,family="binomial")
anova(fit0,fit,test="Chisq") ## Test of sex effect in 0-year olds
Male <- 1-Female
fit1 <- glm(Status~Male*Age,donner,family="binomial")
fit0 <- glm(Status~Male+Age:Male,donner,family="binomial")
anova(fit0,fit1,test="Chisq") ## Test of age effect in females
Age30 <- donner$Age - 30
fit1 <- glm(Status~Age30*Female,donner,family="binomial")
fit0 <- glm(Status~Age30+Age30:Female,donner,family="binomial")
anova(fit0,fit1,test="Chisq") ## Test of sex effect in 30-year olds
Age45 <- donner$Age - 45
fit1 <- glm(Status~Age45*Female,donner,family="binomial")
fit0 <- glm(Status~Age45+Age45:Female,donner,family="binomial")
anova(fit0,fit1,test="Chisq") ## Test of sex effect in 45-year olds

## Inference for probabilities
XX <- data.frame(Age=rep(c(20,40,60),2),Female=c(rep(0,3),rep(1,3)))
predict(fit,newdata=XX,type="response")
pred <- predict(fit,newdata=XX,se.fit=TRUE)
HW <- qnorm(.975)*pred$se.fit
CI <- cbind(pred$fit-HW,pred$fit+HW)
exp(CI)/(1+exp(CI))

## Inference for OR
exp(coef(fit))
exp(confint(fit))
exp(10*c(coef(fit)["Age"],confint(fit,"Age"))) ## 10-yr diff for males
fit1 <- glm(Status~Age*Male,donner,family="binomial")
exp(10*c(coef(fit1)["Age"],confint(fit1,"Age"))) ## 10-yr diff for females
fit1 <- glm(Status~Age30*Female,donner,family="binomial")
exp(c(coef(fit1)["Female"],confint(fit1,"Female"))) ## Sex diff at age 30
lambda <- c(0,-40,1,20) ## 20-yo female vs. 60-yo male
tau <- crossprod(lambda,coef(fit))
SE <- sqrt(lambda%*%vcov(fit)%*%lambda)
exp(c(tau,tau+c(-1,1)*qnorm(.975)*SE))