rm(list=ls())

##########################
##### LAB. SESSION 2: ####
### FITTING LINEAR AND ###
### GENERALIZED LINEAR ###
######## MODELS ##########
####### Exercises ########
##########################

data <- read.table("admissions.txt",header=TRUE)

tab <- table(data$admit,data$rank)
prop.table(tab,margin=2)

gpa.q <- quantile(data$gpa,seq(0,1,by=0.1))
gpa.class <- cut(data$gpa,
             breaks=gpa.q)
gpa.class.mean <- tapply(data$gpa,gpa.class,mean)
admit.gpa.class <- tapply(data$admit,gpa.class,mean)
plot(gpa.class.mean,
     log(admit.gpa.class/(1-admit.gpa.class)))
     
# Do the same for GRE

data$rank <- factor(data$rank)
contrasts(data$rank)[lower.tri(contrasts(data$rank))] <- 1

fit.glm <- glm(admit ~ gpa + gre + rank,
           family=binomial, data=data)
summary(fit.glm)

beta.hat <- coef(fit.glm)
logit.p.hat <- beta.hat[1]+beta.hat[2]*mean(data$gpa)+
               beta.hat[3]*mean(data$gre)+
               beta.hat[4]+beta.hat[5]+beta.hat[6]
p.hat <- exp(logit.p.hat)/(1+exp(logit.p.hat))
p.hat