# Model building:  gss data
# Edps 589
# Fall 2018
# C.J.Anderson
#
#  -- Use gss data
#  -- CMH test
#  -- chi-squares tests
#  -- Linear by linear
#  -- RC(1) association

#  Packages that I will probably use
library(MASS)
library(vcd)
library(vcdExtra)

# Where things will live on laptop
setwd("D:\\Dropbox\\edps 589\\7Model_building")

# Use GSS data
gss <- read.table("gss_data.txt",header=TRUE)
head(gss)

# Make these factors 
gss$fechld <- as.factor(gss$fechld)
gss$mapaid <- as.factor(gss$mapaid)

# Fit of independence
summary(independence <- glm(count~ mapaid + fechld,data=gss,family=poisson))

gss.tab <- xtabs(count ~ fechld + mapaid, data=gss)

# Cochran-Mantel-Haenszel test of association
#   (vcdExtra package)
# Note:	Can get more than just M:  
#           types=c("cor","rmeans","cmeans","general")

CMHtest(gss.tab,
        strata=NULL,
		rscores=1:4,  
		cscores=1:5,
		 types=c("cor","general")
		)
		
# using the fact that M= (n-1)r^2

n <- sum(gss.tab)
( r <- sqrt( 36.26132 /(n-1)) )


# Linear by linear (uniform scores)

# Create scores -- for use later
gss$row.scores <- as.numeric(gss$fechld)
gss$col.scores <- as.numeric(gss$mapaid)

summary(uniform<-glm(count~mapaid + fechld + row.scores*col.scores, data=gss,family=poisson))
anova(independence,uniform)

# RC(1) 
summary( rc1 <- rc(gss.tab,nd=1,
   weighting=c("none"),
   rowsup = NULL, colsup = NULL,
   se = c("jackknife"),
   nreplicates = 100, 
   family = poisson, 
   )
)   
1-pchisq(rc1$deviance,rc1$df.residuals)
plot(rc1, conf.int=0.95,main="RC(1) Association Model + 95% CI on Scale Values")
s


#################################################
# Nominal x Ordinal                             #
################################################# 
# Cochran-Mantel-Haenszel 

CMHtest(gss.tab,
        strata=NULL,
		rscores=1:4,  
		cscores=1:5,
		 types=c("cor","general","rmeans")
		)