# edps 589
# Fall 2018
# c.j.anderson
#
# baseline multinomial logit & equivalent loglinear for 
#   contingency table data
#
library(vcd)
library(nnet)     # for baseline multinomial logistic regression


var <- expand.grid(educ=c("school","working","inactive"),
                   father=c("<=12",">12"),
		           race=c("white","black") )
counts <- c(204, 195, 131,  
            100,  53,  67,
             78,  90,  28,
             12,   5,   9)
w <- cbind(var,counts)
	
#####################
# (RF,E) loglinear  #
#####################
summary(rf.e <- glm(counts ~ race + father + race*father
                             + educ, 
					 data=w,family=poisson))
anova(rf.e,test="LRT")	
###############################
# null multinomial logitistic #
###############################
summary(null <- multinom(educ ~ 1, weight=counts,data = w) )
# z=sqrt(wald) statistics for parameters
(z <- summary(null)$coefficients/summary(null)$standard.errors)
(p <- (1 - pnorm(abs(z), 0, 1)) * 2)

########################
# (RF,RE) loglinear    #
########################
summary(rf.re <- glm(counts ~ race + father + race*father
                             +educ + educ*race , 
					 data=w,family=poisson))
anova(rf.re,test="LRT")	

#############################
# (R) multinomial logistic  #
#############################
summary( mlogit.r <- multinom(educ ~ race , weight=counts, data = w) )
(z <- summary(mlogit.r)$coefficients/summary(mlogit.r)$standard.errors)
(p <- (1 - pnorm(abs(z), 0, 1)) * 2)


########################
# (RF,FE) loglinear    #
########################
summary(rf.fe <- glm(counts ~ race + father + race*father
                             +educ + educ*father, 
					 data=w,family=poisson))
anova(rf.fe,test="LRT")	

#############################
# (F) multinomial logistic  #
#############################
summary( mlogit.f <- multinom(educ ~ father , weight=counts, data = w) )
(z <- summary(mlogit.f)$coefficients/summary(mlogit.f)$standard.errors)
(p <- (1 - pnorm(abs(z), 0, 1)) * 2)


########################
# (RF,RE,FE) loglinear #
########################
summary(rf.re.fe <- glm(counts ~ race + father + race*father
                             +educ + educ*race + educ*father, 
					 data=w,family=poisson))
anova(rf.re.fe,test="LRT")	

#############################
# (R,F) multinomial logistic#
#############################
summary( mlogit.rf <- multinom(educ ~ race + father, weight=counts, data = w) )
(z <- summary(mlogit.rf)$coefficients/summary(mlogit.rf)$standard.errors)
(p <- (1 - pnorm(abs(z), 0, 1)) * 2)