########################################################################
### File loglinear.r                                                 ###
### Tom A.B. Snijders                                                ###
### Further Statistical Methods: Loglinear Models                    ###
### Version February 22, 2011                                        ###
########################################################################

########################################################################
### loglinear models                                                 ###
########################################################################

# First data entry.
# From Andersen (1985) see http://www.ed.uiuc.edu/courses/EdPsy490AT
# worker satisfaction
worker_satisf <- factor(rep(c("Low", "High"), 4))
# supervisor satisfaction
superv_satisf <- factor(rep(c("Low", "Low", "High", "High"), 2))
# quality of management
management <- factor(rep(c("Bad", "Good"), c(4,4)))
# frequencies
freq <- c(103, 87, 32, 42,   59, 109, 78, 205) 

BlueCollar.data <- data.frame(management, superv_satisf, worker_satisf, freq)

library(MASS)

# Various loglinear models
loglm(freq ~ management * superv_satisf * worker_satisf, data = BlueCollar.data)
(lm_2 <- loglm(freq ~ management * superv_satisf + management* worker_satisf
      + superv_satisf * worker_satisf, data = BlueCollar.data))
(lm_3 <- loglm(freq ~ management * superv_satisf + management* worker_satisf,
      data = BlueCollar.data))
loglm(freq ~ management * superv_satisf 
      + superv_satisf * worker_satisf, data = BlueCollar.data)
loglm(freq ~ management* worker_satisf
      + superv_satisf * worker_satisf, data = BlueCollar.data)
loglm(freq ~ management * superv_satisf +worker_satisf,
      data = BlueCollar.data)
loglm(freq ~  management* worker_satisf + superv_satisf,
      data = BlueCollar.data)

anova(lm_2, lm_3)

# Chi-squared tests conditional on management quality;
# directly use the frequency numbers.
chisq.test(matrix(c(103, 87, 32, 42), 2, 2))
chisq.test(matrix(c(59, 109, 78, 205), 2, 2))

########################################################################
### proportional odds model                                          ###
### See Venables & Ripley (2002), second part of Section 7.3         ###
########################################################################

# For a description of the data, see
?housing

# First, as a baseline comparison, we fit the multinomial logit model.

library(nnet)
(house.mult = multinom(Sat ~ Infl + Type + Cont, weights = Freq, data = housing))


# Instead of using nnet::multinom, we now utilize the correspondence
# between loglinear models and multinomial regression,
# and fit the equivalent loglinear model.
# See Venables & Ripley (2002) p. 200-205.
(house.glm <- glm(Freq ~ Infl*Type*Cont + Sat*(Infl+Type+Cont), 
                 family = poisson, data = housing))

# Are the log-odds ratios approximately the same?

# We calculate expected cell counts
house.pm <- predict(house.glm, type="response")
# and arrange these into a matrix
house.pm <- matrix(house.pm, ncol=3, byrow=TRUE)
# and compute probabilities
house.pr <- house.pm/rowSums(house.pm)
# and their cumulative sums
house.cpr <- apply(house.pr, 1, cumsum)
# To compute the log-odds we use
logit <- function(x) log(x/(1-x))
# and apply this to the fitted cumulative probabilities
(house.ld <- logit(house.cpr[2, ]) - logit(house.cpr[1, ]))
# These numbers do not appear to be strongly variable.
# Their mean is
mean(house.ld)
# Thus, the proportional odds assumption seems to be reasonable,
# when compared to this fitted multinomial logit model.

# The proportional odds model can be estimated by the function polr
?polr
summary(house.plr <- 
         polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing),
         digits=3 )
# The residual deviance of house.mult was  3470.084,
# now it is 3479.149
# Loss of 9.1, for 6 parameters less.
# AIC goes down from 3498.08 to 3495.15. 

# Fitted probabilities are obtained from
house.plr.pr <- predict(house.plr, housing, type = "probs")
# but this gives 3 duplicates, so we are more interested in
house.plr.pr <- house.plr.pr[3*(1:24)-2,]
# The difference w.r.t to the earlier predictions is not great: 
house.plr.pr - house.pr
# which all are less than
max(abs(house.plr.pr - house.pr))

# Further model improvements can be investigated from
addterm(house.plr, ~.^2, test = "Chisq")
house.plr2 <- stepAIC(house.plr, ~.^2)
house.plr2$anova
# This indeed leads to considerable improvement:
anova(house.plr, house.plr2)

# The simpler, and less well fitting model, is easier to interpret.
# Let us try to make some interpretation:
summary(house.plr)
# On the logit scale, the influence of householders on management
# of the property has positive effects on satisfaction: 0, 0.57, 1.29;
# Tower blocks are most attractive, the other three types less,
# with terrace houses least attractive 
# (difference w.r.t. tower blocks of 1.09 on the logit scale);
# having more contact adds 0.36 on the logit scale.