> ################################# > #### 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") trying URL 'http://lib.stat.cmu.edu/R/CRAN/bin/macosx/leopard/contrib/2.14/vcd_1.2-13.tgz' Content type 'application/x-gzip' length 1349845 bytes (1.3 Mb) opened URL ================================================== downloaded 1.3 Mb The downloaded packages are in /var/folders/rl/g9pzw2wj2ql6mwznzb_z6h680000gn/T//Rtmpd8EDBP/downloaded_packages > 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 [1] 19 132 11 52 0 9 6 97 > > #### we can represent this table also in 3 dimensions > deathp <- array(deathp, dim=c(2,2,2)) > deathp , , 1 [,1] [,2] [1,] 19 11 [2,] 132 52 , , 2 [,1] [,2] [1,] 0 6 [2,] 9 97 > dimnames(deathp) <- list(DeathPen=c("yes","no"), + Defendant=c("white","black"), + Victim=c("white","black")) > deathp , , Victim = white Defendant DeathPen white black yes 19 11 no 132 52 , , Victim = black Defendant DeathPen white black yes 0 6 no 9 97 > #### create a flat contingency table > ftable(deathp, row.vars=c("Defendant","Victim"), + col.vars="DeathPen") DeathPen yes no Defendant Victim white white 19 132 black 0 9 black white 11 52 black 6 97 > > ### you can get a same as above > ftable(deathp, row.vars=c(1,3), col.vars=2) Defendant white black DeathPen Victim yes white 19 11 black 0 6 no white 132 52 black 9 97 > > > ##### 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 , , DeathPen = yes Victim Defendant white black white 11.59848 6.070232 black 12.03342 6.297866 , , DeathPen = no Victim Defendant white black white 93.4322 48.89909 black 96.9359 50.73281 > > ### chi-squared statistic > chisqr<-sum(((deathexp-deathp)^2)/deathexp) > 1-pchisq(chisqr,4) [1] 0 > > > ####### 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) Pearson's Chi-squared test with Yates' continuity correction data: AB X-squared = 112.5201, df = 1, p-value < 2.2e-16 > chisq.test(AC) Pearson's Chi-squared test with Yates' continuity correction data: AC X-squared = 0.0863, df = 1, p-value = 0.7689 > chisq.test(BC) Pearson's Chi-squared test with Yates' continuity correction data: BC X-squared = 4.7678, df = 1, p-value = 0.029 > > ### compute the odds ratios for AC table for example > ## you can use any functions we have seen previously for two-way tables > assocstats(AC) X^2 df P(> X^2) Likelihood Ratio 0.22145 1 0.63794 Pearson 0.22145 1 0.63794 Phi-Coefficient : 0.026 Contingency Coeff.: 0.026 Cramer's V : 0.026 > oddsratio(AC, log=FALSE) [1] 1.18106 > exp(confint(oddsratio(AC))) lwr upr [1,] 0.5953049 2.343172 > > ### test of conditional independence > ### Cochran-Mantel-Haenszel test) > mantelhaen.test(deathp) Mantel-Haenszel chi-squared test with continuity correction data: deathp Mantel-Haenszel X-squared = 0.7963, df = 1, p-value = 0.3722 alternative hypothesis: true common odds ratio is not equal to 1 95 percent confidence interval: 0.2864136 1.4091583 sample estimates: common odds ratio 0.6352968 > mantelhaen.test(deathp,correct=FALSE) Mantel-Haenszel chi-squared test without continuity correction data: deathp Mantel-Haenszel X-squared = 1.2097, df = 1, p-value = 0.2714 alternative hypothesis: true common odds ratio is not equal to 1 95 percent confidence interval: 0.2864136 1.4091583 sample estimates: common odds ratio 0.6352968 > > ### via odds-ratios > oddsratio(deathp, 3, log=FALSE) white black 0.6804408 0.7894737 > ##log odds ratio for a 2x2 table given the levels of the 3rd variable > lor=oddsratio(deathp,3) > exp(confint(lor)) ## CI lwr upr white 0.30704729 1.50791 black 0.04121472 15.12248 > summary(lor) Log Odds Ratio Std. Error z value Pr(>|z|) white -0.38501 0.40600 -0.9483 0.1715 black -0.23639 1.50644 -0.1569 0.4377 > plot(lor, xlab="Victim") > > > ## getting the the AC table for B=Victim=white > deathp[,,1] ## getting the the AC table for B=Victim=white Defendant DeathPen white black yes 19 11 no 132 52 > oddsratio(deathp[,,1], log=FALSE) ## estimate of a conditional OR [1] 0.6804408 > exp(confint(oddsratio(deathp[,,1]))) ## CI of a conditional OR lwr upr [1,] 0.3070473 1.50791 > chisq.test(deathp[,,1], correct=FALSE) Pearson's Chi-squared test data: deathp[, , 1] X-squared = 0.8774, df = 1, p-value = 0.3489 > > > deathp[,,2] ## getting the the AC table for B=Victim=black Defendant DeathPen white black yes 0 6 no 9 97 > ##In case of zero entries, 0.5 will be added to the table. > oddsratio(deathp[,,2], log=FALSE) ## estimate of a conditional OR [1] 0.7894737 > exp(confint(oddsratio(deathp[,,2]))) ## CI of a conditional OR lwr upr [1,] 0.04121472 15.12248 > chisq.test(deathp[,,2]) Pearson's Chi-squared test with Yates' continuity correction data: deathp[, , 2] X-squared = 8e-04, df = 1, p-value = 0.978 Warning message: In chisq.test(deathp[, , 2]) : Chi-squared approximation may be incorrect > > ######################################### > ### 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 \$Frequency penalty defendant no yes black 149 17 white 141 19 \$RowPct penalty defendant no yes black 0.8975904 0.1024096 white 0.8812500 0.1187500 > result=chisq.test(data) > result Pearson's Chi-squared test with Yates' continuity correction data: data X-squared = 0.0863, df = 1, p-value = 0.7689 > > ### 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 \$Frequency victim defendant black white black 103 63 white 9 151 \$Percent victim defendant black white black 0.31595092 0.19325153 white 0.02760736 0.46319018 \$RowPct victim defendant black white black 0.6204819 0.3795181 white 0.0562500 0.9437500 \$ColPct victim defendant black white black 0.91964286 0.29439252 white 0.08035714 0.70560748 > result=chisq.test(data) > summary(result) Length Class Mode statistic 1 -none- numeric parameter 1 -none- numeric p.value 1 -none- numeric method 1 -none- character data.name 1 -none- character observed 4 xtabs numeric expected 4 -none- numeric residuals 4 xtabs numeric stdres 4 xtabs numeric > result Pearson's Chi-squared test with Yates' continuity correction data: data X-squared = 112.5201, df = 1, p-value < 2.2e-16 > > ### 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 \$Frequency penalty victim no yes black 106 6 white 184 30 \$Percent penalty victim no yes black 0.32515337 0.01840491 white 0.56441718 0.09202454 \$RowPct penalty victim no yes black 0.94642857 0.05357143 white 0.85981308 0.14018692 \$ColPct penalty victim no yes black 0.3655172 0.1666667 white 0.6344828 0.8333333 > result=chisq.test(data) > summary(result) Length Class Mode statistic 1 -none- numeric parameter 1 -none- numeric p.value 1 -none- numeric method 1 -none- character data.name 1 -none- character observed 4 xtabs numeric expected 4 -none- numeric residuals 4 xtabs numeric stdres 4 xtabs numeric > result Pearson's Chi-squared test with Yates' continuity correction data: data X-squared = 4.7678, df = 1, p-value = 0.029 > > ### 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 \$Frequency defendant penalty black white no 97 9 yes 6 0 \$RowPct defendant penalty black white no 0.9150943 0.08490566 yes 1.0000000 0.00000000 > result=chisq.test(data) Warning message: In chisq.test(data) : Chi-squared approximation may be incorrect > result Pearson's Chi-squared test with Yates' continuity correction data: data X-squared = 8e-04, df = 1, p-value = 0.978 > > oddsratio(table) Error in lor(x) : (list) object cannot be coerced to type 'double' > > > ### 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)) Error in list(Frequency = data, , RowPct = prop.table(data, 1)) : argument 2 is empty > table \$Frequency defendant penalty black white no 97 9 yes 6 0 \$RowPct defendant penalty black white no 0.9150943 0.08490566 yes 1.0000000 0.00000000 > result=chisq.test(data) > result Pearson's Chi-squared test with Yates' continuity correction data: data X-squared = 0.5194, df = 1, p-value = 0.4711 > > ### Mantel-Haenszel Test on victim*defendant*penalty > mantelhaen.test(data3) Mantel-Haenszel chi-squared test with continuity correction data: data3 Mantel-Haenszel X-squared = 0.7963, df = 1, p-value = 0.3722 alternative hypothesis: true common odds ratio is not equal to 1 95 percent confidence interval: 0.2864136 1.4091583 sample estimates: common odds ratio 0.6352968 > > ### 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) Breslow and Day test (with Tarone correction): Breslow-Day X-squared = 0.3805608 Breslow-Day-Tarone X-squared = 0.3566714 Test for test of a common OR: p-value = 0.5503607