# Edps 589
# Fall 2018
# c.j.anderson
#
# Log-linear models and gss presidential choice data
#
library(vcd)
library(MASS)   # to use loglm
library(rms)    # see what this gives when using "Glm")   

setwd("C:\\Users\\cja\\Documents\\Dropbox\\edps 589\\6 loglinear")


###########################################################
# GSS 1989 Presidential Race between Bush, Clinton & Perot#
###########################################################

###########################################################
# 2 - way table                                           #
###########################################################
# create cross-classification
var.levels <- expand.grid(view=c("liberal","moderate","conservative"),
                          choice=c("Bush","Clinton","Perot"))
count <- c(70,195,382,342,332,199,56,101,117)
(gss <- cbind(var.levels,count))

# as a data frame
gss.data <- as.data.frame(gss)


# Display table
view_by_choice <- xtabs(count ~  view + choice, data=gss)
addmargins(view_by_choice)

# Table of proportions
pgss0 <- prop.table(view_by_choice)
round(pgss0,2)

# Row wise table of proportions
pgss <- prop.table(view_by_choice,1)
round(pgss,2)

# Fit independence model
summary( indep.mod <- glm(count ~ view + choice, data=gss.data, family=poisson) )

# --- if you also want Pearson's Xsq
X2 <- sum(residuals(indep.mod,type=c("pearson"))**2)

# --- table of fitted values
gss.data$fit <- indep.mod$fitted
xtabs(fit ~ view + choice,data=gss.data)

#--- Another way to get both G2 and X2
loglm(count ~view + choice, data=gss.data)

#--- What rms package gives
I.test <- Glm(count~ view + choice, data=gss.data, family=poisson)

###########################################################
# 3-way table                                             #
###########################################################
# create cross-classification
var.levels <- expand.grid(view=c("liberal","moderate","conservative"),
                          gender = c("male","female"),
                          choice=c("Bush","Clinton","Perot"))
count <- c(26, 82,202, 44,113,180,
          121,128, 75,221,204,124,
		   24, 52, 74, 32, 49, 43)
(gss3 <- cbind(var.levels,count))
gss3.tab <- xtabs(count ~ gender + view + choice,data=gss3)

# as a data frame
gss3.data <- as.data.frame(gss3)

# Display 3-way table
structable(choice ~ gender + view, data=gss3.tab)