######################################################################
# NAME:  Tom Loughin                                                 #
# DATE:  2012-01-23                                                  #
# Purpose: Analyze political ideology 3-way table using Poisson      #
#     regression/Loglinear Models, treating all factors as nominal   #
# NOTES:                                                             #
######################################################################

# Enter the data
#alldata <- read.table("\\\\ais-fs1.sfu.ca\\home\\users\\tloughin\\documents\\Dropbox\\master\\chapter4\\R\\Pol_ideol_data.csv", sep = ",", header = TRUE)
alldata <- read.table("C:\\Users\\Tom Loughin\\Dropbox\\master\\chapter4\\R\\Pol_ideol_data.csv", sep = ",", header = TRUE)
head(alldata)

# Saturated Model: GPI
mod.sat <- glm(formula = count ~ gender*party*ideol, family = poisson(link = "log"), data = alldata)

# Homogeneous association model in all 3 associations: GP,GI,PI
mod.homo <- glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), data = alldata)
anova(mod.homo, mod.sat, test = "Chisq")

# Model assuming only PI association G,PI
mod.homo.PI <- glm(formula = count ~ gender + party*ideol, family = poisson(link = "log"), data = alldata)
anova(mod.homo.PI, mod.homo, test = "Chisq")

# Model assuming pairwise independence: G,P,I
mod.indep <- glm(formula = count ~ gender + party + ideol, family = poisson(link = "log"), data = alldata)
anova(mod.indep, mod.homo.PI, test = "Chisq")

# Summary and LR tests of homogeneous association model 
summary(mod.homo)$coefficients
library(car)
Anova(mod.homo)

# Additional sub-models of homogeneous association
mod.homo.PIGI <- glm(formula = count ~ gender*ideol + party*ideol, family = poisson(link = "log"), data = alldata)
Anova(mod.homo.PIGI)
mod.homo.PIGP <- glm(formula = count ~ gender*party + party*ideol, family = poisson(link = "log"), data = alldata)
Anova(mod.homo.PIGP)

# Begin calculations for various ORs and  confidence intervals
# Already did math to determine which coefficients contribute to ORs.
contr.mat <- rbind(c(rep(0,7),1,0,0,0,0,0,0,0,0), 
                   c(rep(0,7),0,-1,0,1,0,0,0,0,0),
                   c(rep(0,7),0,0,0,1,0,0,0,0,0),
                   c(rep(0,7),0,0,-1,1,0,0,0,0,0),
                   c(rep(0,7),0,0,0,1,-1,0,0,0,0),
                   c(rep(0,7),0,1,0,0,0,0,0,0,0),
                   c(rep(0,7),0,1,-1,0,0,0,0,0,0),
                   c(rep(0,7),0,1,0,0,-1,0,0,0,0),
                   c(rep(0,7),0,0,-1,0,0,0,0,0,0),
                   c(rep(0,7),0,0,0,0,-1,0,0,0,0),
                   c(rep(0,7),0,0,1,0,-1,0,0,0,0),
                   c(rep(0,7),0,0,0,0,0,-1,0,1,0),
                   c(rep(0,7),0,0,0,0,0,0,0,1,0),
                   c(rep(0,7),0,0,0,0,0,0,-1,1,0),
                   c(rep(0,7),0,0,0,0,0,0,0,1,-1),
                   c(rep(0,7),0,0,0,0,0,1,0,0,0),
                   c(rep(0,7),0,0,0,0,0,1,-1,0,0),
                   c(rep(0,7),0,0,0,0,0,1,0,0,-1),
                   c(rep(0,7),0,0,0,0,0,0,-1,0,0),
                   c(rep(0,7),0,0,0,0,0,0,0,0,-1),
                   c(rep(0,7),0,0,0,0,0,0,1,0,-1)
                   )
row.names(contr.mat) <- c("GP Rep | M:F", "GI VC:SC | M:F", "GI VC:M | M:F", "GI VC:SL | M:F", "GI VC:VL | M:F", "GI SC:M | M:F", "GI SC:SL | M:F", "GI SC:VL | M:F", "GI M:SL | M:F", "GI M:VL | M:F", "GI SL:VL | M:F", "PI VC:SC | R:D", "PI VC:M | R:D", "PI VC:SL | R:D", "PI VC:VL | R:D", "PI SC:M | R:D", "PI SC:SL | R:D", "PI SC:VL | R:D", "PI M:SL | R:D", "PI M:VL | R:D", "PI SL:VL | R:D")

#library(mcprofile)
#LRCI <- mcpcalc(mod.homo, CM = contr.mat)
#exp(confint(LRCI)) #Note errors and NA's.  Procedure not reliable here.

# Wald inferences using multcomp package

library(multcomp)
wald <- glht(mod.homo, linfct = contr.mat)

# Defaults use multiplicity adjustment for simultaneous confidence level
summary(wald)
exp(confint(wald)$confint)
# Options to get unadjusted (univariate) tests and CIs
summary(wald, test = univariate())
exp(confint(wald, calpha = qnorm(0.975))$confint)


# Wald CIs
# Get out coefficients and variances
beta <- matrix(coef(mod.homo), ncol = 1)
v.beta <- vcov(mod.homo)
# Estimate Lin Combos and standard errors as matrix computations
log.contrasts <- contr.mat %*% beta
SElog.contrasts <- matrix(sqrt(diag(contr.mat %*% v.beta %*% t(contr.mat))), ncol = 1)
# Compute confidence intervals in linear predictor scale
alpha = 0.05
lower.log <- log.contrasts + qnorm(alpha/2)*SElog.contrasts
upper.log <- log.contrasts + qnorm(1-alpha/2)*SElog.contrasts
# Combine Lin Combo coefficients, estimates of contrasts, and confidence intervals in mean scale
wald.ci <- round(data.frame(exp(log.contrasts), exp(lower.log), exp(upper.log)), digits = 2)
# Attach contrast names to and columns.                  
colnames(wald.ci) <- c("Estimate", "Lower CI", "Upper CI")
wald.ci