### Berkeley admissions data for log-linear models

# dataset already exists in R library
UCBAdmissions 
berk.data = as.data.frame(UCBAdmissions)

# setting references to last category (same as SAS)
berk.data$Gender = relevel(berk.data$Gender, ref='Female')
berk.data$Dept = relevel(berk.data$Dept, ref='F')

### To test the odds-ratios in the marginal table and each of the subtables
library(vcd)

# saturated log-linear model 
berk.sat = glm(Freq~Admit*Gender*Dept, family=poisson(), data=berk.data)
summary(berk.sat)
fitted(berk.sat)

# complete independence model
berk.ind = glm(Freq~Admit+Gender+Dept, family=poisson(), data=berk.data)
summary(berk.ind)
fits = fitted(berk.ind)
resids = residuals(berk.ind,type="pearson")
h = lm.influence(berk.ind)$hat
adjresids = resids/sqrt(1-h)
round(cbind(berk.data$Freq,fits,adjresids),2)

# joint independence model
berk.jnt = glm(Freq~Admit+Gender+Dept+Gender*Dept, family=poisson(), data=berk.data)
summary(berk.jnt)
fits = fitted(berk.jnt)
resids = residuals(berk.jnt,type="pearson")
h = lm.influence(berk.jnt)$hat
adjresids = resids/sqrt(1-h)
round(cbind(berk.data$Freq,fits,adjresids),2)

# conditional independence model
berk.cnd = glm(Freq~Admit+Gender+Dept+Admit*Gender+Dept*Gender, family=poisson(), data=berk.data)
summary(berk.cnd)
fits = fitted(berk.cnd)
resids = residuals(berk.cnd,type="pearson")
h = lm.influence(berk.cnd)$hat
adjresids = resids/sqrt(1-h)
round(cbind(berk.data$Freq,fits,adjresids),2)
anova(berk.cnd, test="LR")

# model of homogeneous association
berk.ha = glm(Freq~(Admit+Gender+Dept)^2, family=poisson(), data=berk.data)
summary(berk.ha)
fits = fitted(berk.ha)
resids = residuals(berk.ha,type="pearson")
h = lm.influence(berk.ha)$hat
adjresids = resids/sqrt(1-h)
round(cbind(berk.data$Freq,fits,adjresids),2)