Many experiments can have binary outcomes. If we manipulate one or more experimental factors, we can aggregate these results by subcategories. This leads to contingency tables, which contain the total count per factor combination. These data can be understood in the same way as other ANOVA models, but with the aggregated responses as counts rather than averages. The following section is meant as both an introduction to the topic, showcasing examples of tests and their application.
With count data, we typically model observations as arising from a Poisson distribution, which takes values in \(0, 1, \ldots\). The Poisson distribution has mean and variance \(\mu>0\) and, to ensure the fitted average are positive, regression models typically consider the natural logarithm.
To make things concrete, consider a two-way between-subject design with crossed factors having \(I\) and \(J\) categories, respectively. The mean equation would be
\[\begin{align*} \ln \mu_{ij} = \alpha + \beta_i + \gamma_j + (\beta\gamma)_{ij}, \qquad i=1,\ldots, I; j =1, \ldots, J \end{align*}\] where \(\beta\), \(\gamma\), \((\beta\gamma)\)’s are subject to sum-to-zero categories. One major difference with regular ANOVA is that there will be at most \(IJ\) counts (one for each cell). If we fit a model with an interaction, we will overfit the data and predict back the entries. This is possible for the Poisson distribution because the variance is fully determined by the mean.
There are direct analog to the tests we normally consider: we can compare the full model with all interactions (termed saturated model) with simpler alternatives. The deviance statistic assesses whether the simpler model (null hypothesis) provides adequate fit, relative to the full model (alternative hypothesis).
We can compare different models in the same way than for multi-way between-subject designs. We fit two competing models (one simpler null model, say without interaction) and the alternative which includes the same specification, plus additional terms. We then compare the log-likelihood, a measure of fit, and form a statistic. As for ANOVA, the alternative models fit better (no matter what), but we can assess whether this improvement, due to including \(\nu\) additional parameters, could be due to chance alone. The likelihood ratio statistic assess this change, and our large-sample benchmark. we compare it to a chi-square random variable with \(\nu\) degrees of freedom, denoted \(\chi^2_\nu\).
In contingency tables, we can use Pearson’s \(\chi^2\) test to compare observed counts to postulated expected counts, comparing the statistic to a chi-square distribution.
The statistic is of the form \[P =\sum_{i=1}^I \sum_{j=1}^J \frac{(O_{ij}-E_{ij})^2}{E_{ij}}\] where \(E_{ij}\) is the expected count under the null hypothesis, and \(O_{ij}\) are the observed entries. The degrees of freedom are \(N\) (here \(IJ=N\), the number of cells) minus the number of parameters under the null hypothesis. The large-sample approximation is adequate provided that the expected counts in each subgroup, \(E_{ij}\)’s, are larger than 5.
For contingency tables, the usual effect size is Cramér’s \(V\), which is a transformation of the chi-square statistic to remove the dependence on sample size and varies from \(0\) (no association) to \(1\) (fully determined data). The estimator is typically biased upward, so slightly different estimators (recipe) are available, including many with small-sample corrections.
Example 15.1 (The importance of selling format) We revisit data from Duke and Amir (2023) on selling formats, which compare integrated purchases (e.g., proposing to choose first the item, then choosing the combination) or having a menu at checkout for the quantity (sequential). We look at whether the person purchased the item and aggregate counts in a \(2 \times 2\) table. We are again testing for an interaction versus only main effects.
In the paper, all participants were included for this test, but other were excluded at a latter stage. The database DA23_E2 only includes the 325 participants who were present at all stages.
data(DA23_E2, package = "hecedsm")
tabs <- with(DA23_E2, table(purchased, format))
tabs format
purchased quantity-integrated quantity-sequential
0 100 133
1 66 26# Chi-square test for interaction
chisq.test(tabs)
Pearson's Chi-squared test with Yates' continuity correction
data: tabs
X-squared = 20.786, df = 1, p-value = 5.135e-06The integrated quantity lead to a higher proportion of sales, and this difference is statistically significant at level 1%. The effect size (estimated with a small-sample correction of \(V=0.25\) indicates a moderate effect, even if the \(p\)-value is small since we have 39.8% who bought for quantity-integrated, versus 16.4% for sequential.
Example 15.2 (Spontaneous verbal rehersal of memory tasks) We consider Elliott et al. (2021) multi-lab replication of Flavell, Beach, and Chinsky (1966) study on spontaneous verbalization of children when asked to identify pictures of objects. We pool data from all labs and study the counts as a function of age and frequency. We are interested in assessing whether there is an interaction between the two. Using Pearson \(\chi^2\) test, we can fit the model in which counts in each cell \(E_{ij}\), which corresponds to the average predicted with a model with main effect only (overall average + row \(i\) and column \(j\) estimated deviation, from the main model).
data(MULTI21_D1, package = "hecedsm")
contingency <- xtabs(
count ~ age + frequency,
data = MULTI21_D1)
# Use chi-square directly - the correction is not applied
# to get the same result as with the Poisson regression model
chisqtest <- chisq.test(contingency, correct = FALSE)
chisqtest
Pearson's Chi-squared test
data: contingency
X-squared = 87.467, df = 6, p-value < 2.2e-16# Effect size (with adjustment for small sample)
effectsize::cramers_v(chisqtest)$Cramers_v_adjusted[1] 0.2043891# Aggregate data into long format
MULTI21_D1_long <- MULTI21_D1 |>
dplyr::group_by(age, frequency) |>
dplyr::summarize(total = sum(count))
mod_main <- glm(total ~ age + frequency,
family = poisson,
data = MULTI21_D1_long)
PearsonX2 <- sum(residuals(mod_main, type = "pearson")^2)
# p-value for Pearson chi-square test
pval_score <- pchisq(PearsonX2, df = mod_main$df.residual, lower.tail = FALSE)
# p-value for likelihood ratio test
pval_lrt <- pchisq(deviance(mod_main), df = mod_main$df.residual, lower.tail = FALSE)The two statistics, likelihood ratio test and Pearson’s chi-square tests, give similar answers (they are two ways of assessing the same hypothesis, and both have \(\chi^2\) distributions under the null hypothesis of no interaction.
Given our large sample size, it is unsurprising to find differences. Perhaps more interesting is looking only at 5 years old versus 7 years old, as this is the age range where most changes occur. We also report effect sizes.
# Take only a subset of the levels for age
counts_5vs7 <- xtabs(
count ~ age + frequency,
data = MULTI21_D1,
subset = age %in% c("5yo", "7yo"),
drop.unused.levels = TRUE)
counts_5vs7 frequency
age never sometimes usually
5yo 53 80 87
7yo 19 73 177# Compute chi-square test with 2x3 sub-table
chisq.test(counts_5vs7)
Pearson's Chi-squared test
data: counts_5vs7
X-squared = 42.575, df = 2, p-value = 5.688e-10# Even stronger evidence of more verbalization
# compute effect size - Cramer's V without small sample adjustment
effectsize::cramers_v(chisq.test(counts_5vs7), adjust = FALSE)$Cramers_v[1] 0.2950689Example 15.3 (Racial discrimination in hiring) We consider a study from Bertrand and Mullainathan (2004), who study racial discrimination in hiring based on the consonance of applicants names; a similar example was recently featured in selection at the M.Sc. level from Ondes. The authors created curriculum vitae for four applicants and randomly allocated them a name, either one typical of a white person or a black person. The response is a count indicating how many of the applicants were called back (out of two black and two white) depending on their origin.
If there was no racial discrimination (null hypothesis), we would expect the average number of times a white applicant was called back (but no black applicant) to be the same as a single black applicant (but no white). Only the entries for different numbers of call-back (either 0 vs 2, 0 vs 1 or 1 vs 2 for either race) are instructive about our question of interest.
| no white | 1 white | 2 white | |
|---|---|---|---|
| no black | 1103 | 74 | 19 |
| 1 black | 33 | 46 | 18 |
| 2 black | 6 | 7 | 17 |
Under the null hypothesis, the model would have equal off-diagonal entries in Table 15.1. We can test this null hypothesis of symmetry by creating a factor that has similar levels for off-diagonal entries.
data(BM04_T2, package = "hecedsm")
# Fit the Poisson regression models
# Saturated model - one average per cell, 9 parameters
mod_alt <- glm(count ~ white*black,
data = BM04_T2,
family = poisson) # default log-link
# Symmetric model with 6 parameters (3 diag + 3 upper triangular)
mod_null <- glm(count ~ gnm::Symm(black, white),
data = BM04_T2,
family = poisson)
# Compare the two nested models using a likelihood ratio test
anova(mod_null, mod_alt, test = "LRT")Analysis of Deviance Table
Model 1: count ~ gnm::Symm(black, white)
Model 2: count ~ white * black
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 3 28.232
2 0 0.000 3 28.232 3.246e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# The saturated model need not be fitted - use deviance and residuals directly
pval <- pchisq(deviance(mod_null), df = mod_null$df.residual, lower.tail = FALSE)We could use Pearson’s \(\chi^2\) test too. In our example, the entries \(E_{ij}\) for off-diagonal average counts entries \(E_{01}\), \(E_{12}\), \(E_{02}\), etc., are obtained as the average of counts \(E_{ij} = (O_{ij} + O_{ji})/2\) \((i\neq j)\), etc. The numeric value of the test statistic is 0, yielding a \(p\)-value less than \(10^{-5}\). There is thus strong evidence of racial discrimination.
## Alternative test (Pearson)
# Compute observed counts minus expected counts under model sum_i (O_i-E_i)^2/E_i
# E_i is the expected count, estimated via average of non-diagonal entries 0.5*mu(0,2) + 0.5*mu(2,0)
PearsonX2 <- sum(residuals(mod_null, type = "pearson")^2)
pchisq(PearsonX2, df = 3, lower.tail = FALSE)We can also compute contrasts. Because the model is multiplicative, it makes sense to report departures of symmetry by computing ratios \(\widehat{\mu}_{ji}/\widehat{\mu}_{ij}\) for \(i, j = 0, 1, 2\) and \(i \neq j\). The grid specification is as usual. Wald tests (computed on the log scale) are also reported.
library(emmeans)
# Compute contrasts as ratios counts01/counts10
emmeans(mod_alt, specs = c("black", "white")) |>
# Compute custom contrasts - USE THIS ORDER FOR CONTRASTS
contrast(method = # Compute joint test of
list("0vs1" = c(0,1,0,-1,0,0,0,0,0),
#(1 white, 0 black) vs (0 white, 1 black)
"0vs2" = c(0,0,1,0,0,0,-1,0,0),
"1vs2" = c(0,0,0,0,0,1,0,-1,0)),
type = "response") contrast ratio SE df null z.ratio p.value
0vs1 0.446 0.0933 Inf 1 -3.858 0.0001
0vs2 0.316 0.1479 Inf 1 -2.461 0.0138
1vs2 0.389 0.1732 Inf 1 -2.120 0.0340
Tests are performed on the log scale In more general regression models, we build a design matrix with explanatory variables \(\mathrm{X}_1, \ldots, \mathrm{X}_p\). The mean model then reads \[\begin{align*} \mu = \exp(\beta_0 + \beta_1 \mathrm{X}_{1} + \cdots + \beta_p \mathrm{X}_{p}), \end{align*}\] so the mean is multiplied by \(\exp(\beta_j)\) for an increase of one unit of \(\mathrm{X}_{j}\), ceteris paribus. If \(\beta_j < 0\), \(\exp(\beta_j) < 1\) and so we have a decrease of \(100\cdot(1-\exp(\beta_j))\)% of the mean response. Likewise, if \(\beta_j>0\), the mean number increases by \(100\cdot(\exp(\beta_j)-1)\)%.
Example 15.4 (Road accidents and speed limits on the motorway in Sweden) Sweden is a worldwide leader in road safety and has a long history of countermeasures to increase road traffic safety, including the Vision Zero program. Back in the 1960s, a study was conducted by the authorities to investigate the potential of speed limits on motorways to reduce the number of accidents. The sweden data contains the number of accidents on 92 matching days in both 1961 and 1962 (Svensson 1981); speed limits were in place on selected days in either year.
To study the impact of the restrictions we can fit a Poisson model. Let \(Y_{i1}\) (respectively \(Y_{i2}\)) denote the number of accidents in 1961 (1962) on day \(i\) and let \(\texttt{limit}_{ij}\) denote a binary indicator equal to one if speed limits were enforced on day \(i\) of year \(j\). We set \[\begin{align*} Y_{i1} &\sim \mathsf{Po}\{\exp(\delta_i + \alpha \texttt{limit}_{i1})\}, \\ Y_{i2} &\sim\mathsf{Po}\{\exp(\delta_i + \gamma + \alpha \texttt{limit}_{i2})\}, \qquad i=1, \ldots, 92. \end{align*}\] The nuisance parameters \(\delta_1, \ldots, \delta_{92}\) control for changes in background number of accidents and are of no practical interest, while \(\gamma\) denotes the change from 1961 to 1962. We are interested here in assessing changes in the number of accidents due to the policy, \(\alpha\); of secondary interest is \(\gamma\), which determines whether there has been a change in the number of accident in 1962 relative to 1961.
data(sweden, package = "hecedsm")
modswed <- glm(accidents ~ -1 + day + limit + year,
family = poisson("log"),
data = sweden)
tab <- car::Anova(modswed, type = 3)| variable | stat | df | p-value |
|---|---|---|---|
| day | 9395.22 | 92 | \(<10^{-5}\) |
| limit | 46.29 | 1 | \(<10^{-5}\) |
| year | 0.70 | 1 | \(0.401\) |
The residual deviance is 107.95 for 90 degrees of freedom, suggests the overall fit is good, despite the large number of nuisance parameters \(\delta_1, \ldots, \delta_{92}\). The coefficient associated to limit is strongly significant: the estimated coefficient is \(\widehat{\alpha}=-0.29\), indicates that speed limits reduce the mean number of accidents by \(25.3\)% on average. In contrast, the likelihood ratio test reported in Table 15.2 shows that the change in the yearly number of accident from 1961 to 1962, \(\gamma\), is not significantly different from zero.