# Paritioning chi-square:  school of psychiatric though x Orgina of Schizopherniauare 
# Edps 589
# Fall 2018
# C.J.Anderson
#
#  library MASS
#          vcd

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

# set libraries
library(MASS)
library(vcd)

# 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 independence in full table
(fit.whole <- glm(count~school+origins,family=poisson,data=schizo))
1-pchisq(fit.whole$deviance,4)


# The partioning:

# Table 1:  sub-table from original table:  (eclectic, medical) x (biogenic, medical)
(fit1 <- update(fit.whole,subset=(school=="eclectic" | school=="medical") &
             (origins=="biogenic" | origins=="environ"))$deviance)
	
1-pchisq(fit1,1)	

	
# Table 2 (biogenic or environ vs combo) x (eclectic, medical)
# --  new variable added to schizo data
schizo$select<- c(rep("bio+env",6), rep("combo",3))
# -- sum over selected for each school
schizo.sub2 <- aggregate(schizo$count, by=list(select=schizo$select,school=schizo$school),sum)
# -- fit model with further sub-setting and ask for deviance
(fit2<- glm(x~select+school, family=poisson, data=schizo.sub2, subset=(school=="eclectic" | school=="medical"))$deviance)
1-pchisq(fit2,1)

# Table 3:   (bio vs env) x (ecl+med vs psychoa) --------------- not working, find error
schizo$select3 <- rep(c("ecl+med","ecl+med","psychoa"),3)
schizo.sub3 <- aggregate(schizo$count, by=list(school=schizo$select3,origin=schizo$origins),sum)
(fit3 <- glm(x~school+origin,family=poisson,data=schizo.sub3, subest=(origin=="biogenic" | origin=="environ")))
1-pchisq(fit3,1)


glm(x ~ school + origin,family=poisson,data=schizo.sub3, subset=(origin=="biogenic" | origin=="environ"))

# Table 4:  (biogen+envi vs combo) x (ecletic+medical vs psyhco)
schizo$select4 <- rep(c("eclectic+medical","eclectic+medical","psychoanal"),3)
(schizo.sub4 <- aggregate(schizo$count,by=list(orgins=schizo$select,sch=schizo$select4),sum))
(fit4<- glm(x~orgins+sch,family=poisson,data=schizo.sub4)$deviance)

1-pchisq(fit4,1)