# Edps 589
# c.j.anderson
#
#  McNemar's test of homogeneity for a 2 x 2 table 
#
#
library(vcd)
library(vcdExtra)

def <- expand.grid(control=c("yes","no"),
                   hazard=c("yes","no"))
cell <- data.frame(def,count=c(13,157,24,505))
cell.tab <- xtabs(count ~ hazard + control,data=cell)
(with.marginals <- addmargins(cell.tab))

# Compare test statistic to chi-square distribution
(mh <- mcnemar.test(cell.tab,correct=FALSE))

# For test statistic as a z ~ N(0,1) as given in lecture notes
round(sqrt(mh$statistic),digits=2)

# Compute difference and CI for proprotions 
(props <- prop.table(cell.tab)) 

(p.row1 <- with.marginals[1,3]/with.marginals[3,3])
(p.col1 <- with.marginals[3,1]/with.marginals[3,3])
(diff <- p.row1 - p.col1)

# se of difference
(var.diff <- (p.row1*(1-p.row1) + p.col1*(1-p.col1) 
            - 2*((props[1,1]*props[2,2])- (props[1,2]*props[2,1])))/with.marginals[3,3])
(se.diff <- sqrt(var.diff))

# 96% CI for difference
(lower <- diff - 1.96*se.diff)
(upper <- diff + 1.96*se.diff)


# Risk of collision
157/24