# Ordinal Association:  school of psychiatric though x Orgina of Schizopherniauare 
# Edps 589
# Fall 2018
# C.J.Anderson
#
#  library MASS
#          vcd
#          vcdextra

# Set working directory
setwd("D:\\Dropbox\\edps 589\\2 Chi-square")

# set libraries 
library(MASS)       # Needed for loglm
library(vcd)
library(vcdExtra)   # Needed for CMHtest
library(gnm)        # Fits Generalized Non-linear Models (is somewhat limited but
                    #  will fit RC association model

# Make data.frame
var.levels <- expand.grid(school=c("eclectic", "medical","psychoa"),
                          origins=c("biogenic","environ","combo"))
( schizo <- data.frame(var.levels, count=c(90,13,19, 12,1,13, 78,6,50))		 )				  


# Test of independence in whole table
loglm(formula = count ~ school + origins, data = schizo)

# Table
schizo.tab <- xtabs(count ~ origins + school, data=schizo)
plot(schizo.tab,main="Schools of Psychiatric Thought",col="cyan")


n <- sum(schizo.tab)

CMHtest(schizo.tab,
        strata=NULL,
		rscores=1:3,  
		cscores=1:3,
		types=c("cor", "cmeans", "rmeans", "general")
		)
( r.bad.scores <- sqrt(10.735/(n-1))	)	

# Better scores		
CMHtest(schizo.tab,
        strata=NULL,
		rscores=c(1,3,2),  
		cscores=c(2,1,3),
		types=c("cor", "cmeans", "rmeans", "general")
		)
( r.better.scores <- sqrt(20.260/(n-1))	)	
		
# near best possible:  RC-association model
#  -- independence
( rc.0 <- gnm(count ~ school + origins,data=schizo,family=poisson) )
#  -- RC(1) association model
( rc.1 <- gnm(count ~ school + origins + Mult(school,origins),data=schizo,family=poisson) )

# This is not working but I don't see why not:
#    CMHtest can handle non-integer values
#                       negative & positive numbers
CMHtest(schizo.tab,
        strata=NULL,
		rscores=c(-0.4643, 0.4889, 0.2032),
		cscores=c( 0.1285,-0,4181,-1.2660),
		types=c("cor", "cmeans", "rmeans", "general")
		)
( r.c <- sqrt(22.029/(n-1))	)	# using 22.029 from SAS

		
# thet best possible:  Correspondence analysis