OK, now we’re ready to do some analyses. This vignette focuses on
relatively simple non-parametric tests and measures of association.
CrossTable
For tabular displays, the CrossTable() function in the
gmodels package produces cross-tabulations modeled after
PROC FREQ in SAS or CROSSTABS in SPSS. It has
a wealth of options for the quantities that can be shown in each
cell.
Recall the GSS data used earlier.
# Agresti (2002), table 3.11, p. 106
GSS <- data.frame(
expand.grid(sex = c("female", "male"),
party = c("dem", "indep", "rep")),
count = c(279,165,73,47,225,191))
(GSStab <- xtabs(count ~ sex + party, data=GSS))
## party
## sex dem indep rep
## female 279 73 225
## male 165 47 191
Generate a cross-table showing cell frequency and the cell
contribution to \(\chi^2\).
# 2-Way Cross Tabulation
library(gmodels)
CrossTable(GSStab, prop.t=FALSE, prop.r=FALSE, prop.c=FALSE)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## |-------------------------|
##
##
## Total Observations in Table: 980
##
##
## | party
## sex | dem | indep | rep | Row Total |
## -------------|-----------|-----------|-----------|-----------|
## female | 279 | 73 | 225 | 577 |
## | 1.183 | 0.078 | 1.622 | |
## -------------|-----------|-----------|-----------|-----------|
## male | 165 | 47 | 191 | 403 |
## | 1.693 | 0.112 | 2.322 | |
## -------------|-----------|-----------|-----------|-----------|
## Column Total | 444 | 120 | 416 | 980 |
## -------------|-----------|-----------|-----------|-----------|
##
##
There are options to report percentages (row, column, cell), specify
decimal places, produce Chi-square, Fisher, and McNemar tests of
independence, report expected and residual values (pearson,
standardized, adjusted standardized), include missing values as valid,
annotate with row and column titles, and format as SAS or SPSS style
output! See help(CrossTable) for details.
Chi-square test
For 2-way tables you can use chisq.test() to test
independence of the row and column variable. By default, the \(p\)-value is calculated from the asymptotic
chi-squared distribution of the test statistic. Optionally, the \(p\)-value can be derived via Monte Carlo
simulation.
(HairEye <- margin.table(HairEyeColor, c(1, 2)))
## Eye
## Hair Brown Blue Hazel Green
## Black 68 20 15 5
## Brown 119 84 54 29
## Red 26 17 14 14
## Blond 7 94 10 16
chisq.test(HairEye)
##
## Pearson's Chi-squared test
##
## data: HairEye
## X-squared = 138.29, df = 9, p-value < 2.2e-16
chisq.test(HairEye, simulate.p.value = TRUE)
##
## Pearson's Chi-squared test with simulated p-value (based on 2000
## replicates)
##
## data: HairEye
## X-squared = 138.29, df = NA, p-value = 0.0004998
Fisher Exact Test
fisher.test(X) provides an exact test
of independence. X must be a two-way contingency table in
table form. Another form, fisher.test(X, Y) takes two
categorical vectors of the same length.
For tables larger than \(2 \times 2\)
the method can be computationally intensive (or can fail) if the
frequencies are not small.
fisher.test(GSStab)
##
## Fisher's Exact Test for Count Data
##
## data: GSStab
## p-value = 0.03115
## alternative hypothesis: two.sided
Fisher’s test is meant for tables with small total sample size. It
generates an error for the HairEye data with \(n\)=592 total frequency.
fisher.test(HairEye)
## Error in fisher.test(HairEye): FEXACT error 6 (f5xact). LDKEY=618 is too small for this problem: kval=238045028.
## Try increasing the size of the workspace.
Mantel-Haenszel test and conditional association
Use the mantelhaen.test(X) function to perform a
Cochran-Mantel-Haenszel \(\chi^2\) chi
test of the null hypothesis that two nominal variables are
conditionally independent, \(A \perp
B \; | \; C\), in each stratum, assuming that there is no
three-way interaction. X is a 3 dimensional contingency
table, where the last dimension refers to the strata.
The UCBAdmissions serves as an example of a \(2 \times 2 \times 6\) table, with
Dept as the stratifying variable.
# UC Berkeley Student Admissions
mantelhaen.test(UCBAdmissions)
##
## Mantel-Haenszel chi-squared test with continuity correction
##
## data: UCBAdmissions
## Mantel-Haenszel X-squared = 1.4269, df = 1, p-value = 0.2323
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.7719074 1.0603298
## sample estimates:
## common odds ratio
## 0.9046968
The results show no evidence for association between admission and
gender when adjusted for department. However, we can easily see that the
assumption of equal association across the strata (no 3-way association)
is probably violated. For \(2 \times 2 \times
k\) tables, this can be examined from the odds ratios for each
\(2 \times 2\) table
(oddsratio()), and tested by using
woolf_test() in vcd.
oddsratio(UCBAdmissions, log=FALSE)
## odds ratios for Admit and Gender by Dept
##
## A B C D E F
## 0.3492120 0.8025007 1.1330596 0.9212838 1.2216312 0.8278727
lor <- oddsratio(UCBAdmissions) # capture log odds ratios
summary(lor)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## A -1.052076 0.262708 -4.0047 6.209e-05 ***
## B -0.220023 0.437593 -0.5028 0.6151
## C 0.124922 0.143942 0.8679 0.3855
## D -0.081987 0.150208 -0.5458 0.5852
## E 0.200187 0.200243 0.9997 0.3174
## F -0.188896 0.305164 -0.6190 0.5359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
woolf_test(UCBAdmissions)
##
## Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)
##
## data: UCBAdmissions
## X-squared = 17.902, df = 5, p-value = 0.003072
Some plot methods
Fourfold displays
We can visualize the odds ratios of Admission for each department
with fourfold displays using fourfold(). The cell
frequencies \(n_{ij}\) of each \(2 \times 2\) table are shown as a quarter
circle whose radius is proportional to \(\sqrt{n_{ij}}\), so that its area is
proportional to the cell frequency.
UCB <- aperm(UCBAdmissions, c(2, 1, 3))
dimnames(UCB)[[2]] <- c("Yes", "No")
names(dimnames(UCB)) <- c("Sex", "Admit?", "Department")
Confidence rings for the odds ratio allow a visual test of the null
of no association; the rings for adjacent quadrants overlap iff
the observed counts are consistent with the null hypothesis. In the
extended version (the default), brighter colors are used where the odds
ratio is significantly different from 1. The following lines produce .
col <- c("#99CCFF", "#6699CC", "#F9AFAF", "#6666A0", "#FF0000", "#000080")
fourfold(UCB, mfrow=c(2,3), color=col)
Another vcd function, cotabplot(), provides
a more general approach to visualizing conditional associations in
contingency tables, similar to trellis-like plots produced by
coplot() and lattice graphics. The panel
argument supplies a function used to render each conditional subtable.
The following gives a display (not shown) similar to .
cotabplot(UCB, panel = cotab_fourfold)
Doubledecker plots
When we want to view the conditional probabilities of a response
variable (e.g., Admit) in relation to several factors, an
alternative visualization is a doubledecker() plot. This
plot is a specialized version of a mosaic plot, which highlights the
levels of a response variable (plotted vertically) in relation to the
factors (shown horizontally). The following call produces , where we use
indexing on the first factor (Admit) to make
Admitted the highlighted level.
In this plot, the association between Admit and
Gender is shown where the heights of the highlighted
conditional probabilities do not align. The excess of females admitted
in Dept A stands out here.
doubledecker(Admit ~ Dept + Gender, data=UCBAdmissions[2:1,,])
Odds ratio plots
Finally, the there is a plot() method for
oddsratio objects. By default, it shows the 95% confidence
interval for the log odds ratio. is produced by:
plot(lor,
xlab="Department",
ylab="Log Odds Ratio (Admit | Gender)")
