# Edps 589
# Fall 2021
# c.j.anderson
#
# Contents (all using vglm):
#  o Proportional odds model
#  o Adjacent-categories model for ordinal responses
#  o Multinomial (baseline)
#
# Happiness Data from Agresti (2013) http://www.stat.ufl.edu/~aa/cda/data.html
#    which is from GSS.
# Table 8.5
#  Varnames:    race trauma happy 
#               where  race is 0=white, 1=black
#
library(MASS)
library(VGAM)

setwd("D:\\Dropbox\\edps 589\\8 Multicategory_logit")

gss <- read.table("happiness_data.txt",header=TRUE)

table(gss$happy)

table(gss$race)
table(gss$trauma)


###########################
# Proportional Odds model #
###########################

summary(po1 <- vglm(happy ~ race + trauma, data=gss, family=cumulative(parallel=TRUE)))

summary(po2 <- vglm(happy ~ race + trauma, data=gss, family=cumulative(parallel=FALSE)))

# Quick check on proportional odds assumption...
# deviance po1= 148.407  on 190 degrees of freedom
# deviance po2= 146.9951 on 188 
(lr <- 148.407-146.9951)
1-pchisq(lr,2)

###########################
# Adjacent categories                    #
###########################

summary(adj.cat1 <- vglm(happy ~ race + trauma, data=gss, family=acat(parallel=TRUE, reverse=TRUE)))

#odds ratio for race (see notes for interpretation)
exp(1.8423)
#odds ratio for number of traumtic events
exp(.3570)

# Relax assumption on equality of slopes
summary(adj.cat2 <- vglm(happy ~ race + trauma, data=gss, family=acat(parallel=FALSE, reverse=TRUE)))

# Warning: Hauck-Donner effect detected in the following estimate(s): '(Intercept):2', 'race:1'
# This means that the Wald statistics does not increase monotonically with larger differences between
#  null and estimated parameter.  Can result in a large effect be said to be non-statistically significant.
#  Likelihood ratio test does not suffer from this.


# Difference in deviances
lr <- 148.1996 -  146.8737
df <- 190-188
1-pchisq(lr,df)

###########################
# Baseline multinomial    #
###########################

summary(multinom <- vglm(happy ~ race + trauma, data=gss, family=multinomial))

#### Compare results of adj.cat2 and baseline multinomial...they are same...
#### Same deviance, df, and a simple transformation for parameters of one to the
#### the other.