################################# #### Death Penalty Example: a 2x2x2 table #### let A=defendant, B=vicitim, C=penalty #### A simple line by line analysis #### Nice R code that corresponds to SAS output ############################################### install.packages("vcd") library(vcd) #### reading the data into a table #### first step is just a vector of values deathp <- c(19,132, 11,52,0,9, 6,97) deathp #### we can represent this table also in 3 dimensions deathp <- array(deathp, dim=c(2,2,2)) deathp dimnames(deathp) <- list(DeathPen=c("yes","no"), Defendant=c("white","black"), Victim=c("white","black")) deathp #### create a flat contingency table ftable(deathp, row.vars=c("Defendant","Victim"), col.vars="DeathPen") ### you can get a same as above ftable(deathp, row.vars=c(1,3), col.vars=2) ##### test of mutual independence ### these are the marginal counts penalty<-margin.table(deathp,1) defand<-margin.table(deathp, 2) victim<-margin.table(deathp,3) ### expected value under the mutual independence model deathexp<-(defand%o%victim%o%penalty)/(sum(deathp)^2) deathexp ### chi-squared statistic chisqr<-sum(((deathexp-deathp)^2)/deathexp) 1-pchisq(chisqr,4) ####### test of mutual independence via marginal tables ### recall that this means that we need to test that ALL odds ratios are equal to 1 ### since this is a 2x2x2 table and all marginal tables are 2x2, then we can just ### do the chi-square test of independence ### let A=defendant, B=vicitim, C=penalty AB<-margin.table(deathp, 2:3) AC<-margin.table(deathp, c(2,1)) BC<-margin.table(deathp, c(3,1)) chisq.test(AB) chisq.test(AC) chisq.test(BC) ### compute the odds ratios for AC table for example ## you can use any functions we have seen previously for two-way tables assocstats(AC) oddsratio(AC, log=FALSE) exp(confint(oddsratio(AC))) ### test of conditional independence ### Cochran-Mantel-Haenszel test) mantelhaen.test(deathp) mantelhaen.test(deathp,correct=FALSE) ### via odds-ratios oddsratio(deathp, 3, log=FALSE) ##log odds ratio for a 2x2 table given the levels of the 3rd variable lor=oddsratio(deathp,3) exp(confint(lor)) ## CI summary(lor) plot(lor, xlab="Victim") ## getting the the AC table for B=Victim=white deathp[,,1] ## getting the the AC table for B=Victim=white oddsratio(deathp[,,1], log=FALSE) ## estimate of a conditional OR exp(confint(oddsratio(deathp[,,1]))) ## CI of a conditional OR chisq.test(deathp[,,1], correct=FALSE) deathp[,,2] ## getting the the AC table for B=Victim=black ##In case of zero entries, 0.5 will be added to the table. oddsratio(deathp[,,2], log=FALSE) ## estimate of a conditional OR exp(confint(oddsratio(deathp[,,2]))) ## CI of a conditional OR chisq.test(deathp[,,2]) ######################################### ### Here as another way to handle this same data defendant=rep(c("white","black"),c(4,4)) victim=rep(rep(c("white","black"),c(2,2)),2) penalty=rep(c("yes","no"),4) count=c(19,132,0,9,11,52,6,97) ### Table of defendant by penalty data=xtabs(count~defendant+penalty) table=list(Frequency=data,RowPct=prop.table(data,1)) table result=chisq.test(data) result ### Table of defendant by victim data=xtabs(count~defendant+victim) table=list(Frequency=data,Percent=prop.table(data),RowPct=prop.table(data,1),ColPct=prop.table(data,2)) table result=chisq.test(data) summary(result) result ### Table of victim by penalty data=xtabs(count~victim+penalty) table=list(Frequency=data,Percent=prop.table(data),RowPct=prop.table(data,1),ColPct=prop.table(data,2)) table result=chisq.test(data) summary(result) result ### Table1 of defendant by penalty given victim data=xtabs(count~penalty+defendant+victim) data=addmargins(data,3)[,,1] table=list(Frequency=data,RowPct=prop.table(data,1)) table result=chisq.test(data) result oddsratio(table) ### Table1 of defendant by penalty data3=xtabs(count~penalty+defendant+victim) data=addmargins(data3,3)[,,2] table=list(Frequency=data,,RowPct=prop.table(data,1)) table result=chisq.test(data) result ### Mantel-Haenszel Test on victim*defendant*penalty mantelhaen.test(data3) ### Test for homogenous associations: ### This is not a built in R function. Make sure you compile this function (i.e., breslowday.test(), copy it into R) from the course site first breslowday.test(deathp)