# 3-way table:   marignal dependence, conditional independence (blue colar worker)
# Edps 589
# Fall 2018
# C.J.Anderson
#

library(vcd)
library(vcdExtra)
library(MASS)
library(DescTools)
library(lawstat)

var.values <- expand.grid(worker=c("low","high"),
                          superv=c("low","high"),
			              manager=c("bad","good"))
counts <- c(103, 87, 32, 42, 59, 109, 78, 205)	
bcolar <- cbind(var.values,counts)
bcolar


# 2-way marginal table
worker.by.superv <- xtabs(counts ~ worker+superv,data=bcolar) 
addmargins(worker.by.superv)

# 2-way marginal analysis:  
assocstats(worker.by.superv)             # Compute G2 and X2
( OR <- oddsratio(worker.by.superv)  )   # Compute odds ratio and view
summary(OR)                              # signficant test for odds ratio -- don't really 
                                         #   need because have X2 and G2 from assocstats()
exp(0.6183896)                           # Convert log odds ratio to odds ratio
exp(confint(OR) )                        # Confidnece interval of odds ratio



# 3-way  Table of data
#  -- partial or conditional tables
bcolar.tab <- xtabs(counts ~ worker+superv+manager,data=bcolar)
addmargins(bcolar.tab)

# 3-way contitional/partial table analysis 
assocstats(bcolar.tab)                    # Compute X2 and G2
( OR <- oddsratio(bcolar.tab)  )              # Compute odds ratio & view
( sum.or <-summary(OR)    )                            # Signficant test for odds ratios
exp(sum.or[1:2])
exp( confint(OR) )                        # Confidnece interval of odds ratios

# The set of tests:
# First: Test for homogenity (shouldn't be significant...p should = 1
WoolfTest(bcolar.tab)

# Brewsow-Day
BreslowDayTest(bcolor.tab, OR = NA, correct = FALSE)


# If reatain, then Test for conditional independence
CMHtest(bcolar.tab)

# If data are conditionally independent, estimate their Common odds ratio
mantelhaen.test(bcolar.tab,alternative = c("two.sided"),
                correct = FALSE, exact = FALSE, conf.level = 0.95)

# From lawstat package:  test and common odds ratio
cmh.test(bcolar.tab)