# Farmer et al. data
# Edps 589
# Fall 2018
# C.J.Anderson

library(MASS)
library(vcd)
library(vcdExtra)

setwd("D:\\Dropbox\\edps 589\\2 Chi-square")
farmer <- read.table("farmer_et_al_data.txt",header=TRUE)
farmer

farmer.tab <- xtabs(count ~ gender + prestige, data=farmer)
addmargins(farmer.tab)


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

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

n <- sum(farmer.tab)
( r <- sqrt(19.840 /(n-1)) )


# Other column scores
CMHtest(farmer.tab,
        strata=NULL,
		rscores=1:2,  
		cscores=c(-4,-2,-1,1,2,4),
		types=c("cor", "cmeans", "rmeans", "general")
		)

# Even another column scores
CMHtest(farmer.tab,
        strata=NULL,
		rscores=1:2,  
		cscores=c(-4,-3,-0.5,0.5,3,4),
		types=c("cor", "cmeans", "rmeans", "general")
		)