# Edps 589 - Homework 3
# Fall 2018
# C.J.Anderson
#

library(vcd)            # for xtabs and structable
library(vcdExtra)       # for collapse.table
library(lawstat)        # for mantel-haenszeltest

###########################################
# Problem 1:   2.30 from Agresti (2007)   #
###########################################

rad.def <- expand.grid(Rx = c("surgery","radiation"),
                   cancer=c("controled","not controled"))
				   
(rad <- data.frame(rad.def,count=c(21,15,2,3)))	

rad.tab <- xtabs(count ~ Rx + cancer,data=rad)

(addmargins(rad.tab) )

fisher.test(rad.tab,alternative="greater",conf.int=TRUE,conf.level=.99)
			   

##########################################
# Problem 2:  2.33 from Agresti (2007)   #
##########################################

s.def <- expand.grid(
            victim <- c("white","black"),
            defendent <- c("white","black"),
            death <- c("yes","no")	)

(Dsent <- data.frame(s.def,
            count=c(19,0,11,6, 132,9,52,97)))
# --- names got changed so I change to what I wanted			
names(Dsent) <- c("victim","defendent","death","count")

# part (a)			
(d.tab <- xtabs(count ~   victim  + death+ defendent , data=Dsent))	

# ... and even better
structable(d.tab)

# part (b)
d.plus <- d.tab + 0.5
(sub.white <- d.plus[1, ,] )
# --- take exp of the coefficients, because the comand give log(odds ratio)
(sub.white <- d.plus[1, ,] )
exp(oddsratio(sub.white)$coefficients)

# --- or all in one step --- 
exp(oddsratio(d.plus[1, ,] )$coefficients)

(sub.black <- d.plus[2, , ])
exp(oddsratio(sub.black)$coefficients)


# part (c)
(d.c <- collapse.table(d.tab, victim=c("sum","sum"))[1, , ] )

exp(oddsratio(d.c)$coefficients)


####################################################
# Problem 3:  No computations                      #
####################################################

####################################################
# Problem 4:  No computations                      #
####################################################

####################################################
# Problem 5: 3.8 and 3.9 from Agresti (1996)       #
####################################################

countries <- read.table("D:\\Dropbox\\edps 589\\homework\\R_a3_8_n_9.txt",header=TRUE)
c.tab <- xtabs(count ~  spouse + group + country,data=countries)
structable(c.tab)


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

# Test of conditional independence & common odds ratio
mantelhaen.test(c.tab,alternative = c("two.sided"),
                correct = FALSE, exact = FALSE, conf.level = 0.95)