Chapter 4 Analysis of Contingency Tables and Compact Reporting
This chapter is about categorical data — variables measured in categories rather than numbers. You’ll learn to build contingency tables, test for association and trend, measure agreement between raters, handle paired categorical data, control for confounding variables, and produce clean summary tables for reports.
Most of this chapter uses course_dataset (180 simulated patients across three hospital sites) because it was specifically built to contain every kind of categorical relationship this chapter covers: a treatment/response 2×2 table, a site-stratified design, an ordinal exposure variable, two raters’ severity gradings, and paired pre/post symptom data. We’ll bring in InsectSprays (built into R) for one contrasting example.
#> 'data.frame': 180 obs. of 14 variables:
#> $ id : chr "P001" "P002" "P003" "P004" ...
#> $ site : chr "Site C" "Site B" "Site A" "Site C" ...
#> $ treatment : chr "New" "Standard" "New" "Standard" ...
#> $ age : int 59 68 33 53 58 49 42 40 67 49 ...
#> $ sex : chr "Male" "Male" "Male" "Male" ...
#> $ activity_level : chr "Moderate" "High" "Moderate" "Low" ...
#> $ baseline_pain : num 8.7 6.7 7.6 5.5 6 4.1 6.8 6 5.8 6.6 ...
#> $ week4_pain : num 5.8 6.5 3.9 2.7 5.9 4.6 6.6 2.5 4.4 2.9 ...
#> $ week8_pain : num 5.1 5.6 0.1 2.4 5.1 4.6 4.2 2.4 3.5 1.3 ...
#> $ response : chr "Responder" "Non-Responder" "Responder" "Responder" ...
#> $ symptom_pre : chr "Present" "Absent" "Present" "Present" ...
#> $ symptom_post : chr "Absent" "Present" "Present" "Absent" ...
#> $ rater1_severity: chr "Severe" "Mild" "Moderate" "Severe" ...
#> $ rater2_severity: chr "Severe" "Mild" "Moderate" "Severe" ...4.1 Contingency Table Structure
#>
#> Non-Responder Responder
#> New 4 92
#> Standard 45 39#>
#> Mild Moderate Severe
#> Mild 53 10 2
#> Moderate 8 62 6
#> Severe 3 4 32Adding margins (row/column totals):
#>
#> Non-Responder Responder Sum
#> New 4 92 96
#> Standard 45 39 84
#> Sum 49 131 180Using xtabs() with formula notation (handy for more than two variables):
#> response
#> treatment Non-Responder Responder
#> New 4 92
#> Standard 45 39Manual Table Generation Using Loops
Before using R’s built-in table() or xtabs() functions, it is highly instructive to see how R builds a contingency table algorithmically. We can initialize a blank matrix and loop through the dataset row by row, incrementing the counts:
# 1. Initialize an empty 2x2 matrix with zeros
row_labels <- c("New", "Standard")
col_labels <- c("Non-Responder", "Responder")
manual_tab <- matrix(0, nrow = 2, ncol = 2,
dimnames = list(row_labels, col_labels))
# 2. Loop through each row of the dataset and increment counts
for (i in 1:nrow(course_df)) {
trt <- course_df$treatment[i]
resp <- course_df$response[i]
if (trt %in% row_labels && resp %in% col_labels) {
manual_tab[trt, resp] <- manual_tab[trt, resp] + 1
}
}
manual_tab#> Non-Responder Responder
#> New 4 92
#> Standard 45 39This manual loop matches the output of table(course_df$treatment, course_df$response) exactly, showing the underlying logic of cross-tabulation.
4.2 Chi-Square and Fisher’s Exact Test
The Chi-Square (\(\chi^2\)) Test of Independence is used to determine whether there is a statistically significant association between two categorical variables.
- Null Hypothesis (\(H_0\)): The two variables are independent (no association).
- Alternative Hypothesis (\(H_1\)): The two variables are dependent (associated).
4.2.1 Step 1: Running the Chi-Square Test
Let’s run the test on our tab_2x2 (Treatment vs. Pain Relief):
#>
#> Pearson's Chi-squared test with Yates' continuity correction
#>
#> data: tab_2x2
#> X-squared = 52.729, df = 1, p-value = 3.83e-13Check Before You Run This (Chi-Square & Fisher’s Diagnostics) Before running a Chi-Square test of independence: 1. Categorical Variables: Both variables must be nominal or ordinal. 2. Independence of Observations: Each subject must contribute to exactly one cell in the table (no repeated measures, which require McNemar’s test). 3. Adequate Expected Counts: The expected count in each cell must be \(\ge 5\) (or at least 80% of cells \(\ge 5\) and no cells \(< 1\)). If this is violated, use Fisher’s Exact Test instead.
4.2.2 Step 2: Understanding the Output
The output lists three key values:
1. X-squared: The test statistic (\(\chi^2 = 52.73\)). A larger value indicates a greater difference between what we observed and what we would expect if the variables were independent.
2. df: Degrees of freedom. For a table of \(r\) rows and \(c\) columns, \(\text{df} = (r-1) \times (c-1)\). In our \(2 \times 2\) table, \(\text{df} = (2-1) \times (2-1) = 1\).
3. p-value: The probability of observing such an extreme association by chance. Here, \(p = 3.83 \times 10^{-13}\), which is extremely small (far below the standard \(0.05\) threshold), indicating a very strong, highly significant association.
The Math Behind It: Manual Chi-Square Calculation & Yates’ Correction The mathematical formula for the standard Pearson Chi-Square statistic is: \[\chi^2 = \sum \frac{(O - E)^2}{E}\] where \(O\) represents the observed cell counts and \(E\) represents the expected cell counts.
However, for \(2 \times 2\) tables, R’s chisq.test() defaults to applying Yates’ continuity correction to prevent overestimating significance in small tables. The corrected formula is:
\[\chi^2_{\text{Yates}} = \sum \frac{(|O - E| - 0.5)^2}{E}\]
Let’s compute both manually in R to see how they match R’s output:
# 1. Extract observed and expected values as vectors
obs <- as.vector(tab_2x2)
exp <- as.vector(chisq_result$expected)
# 2. Standard Pearson Chi-Square (without correction)
chi_standard <- sum((obs - exp)^2 / exp)
chi_standard#> [1] 55.19418#> [1] 1.091917e-13# 3. Yates' Corrected Chi-Square (matching R's default chisq.test output)
chi_yates <- sum((abs(obs - exp) - 0.5)^2 / exp)
chi_yates#> [1] 52.72863#> [1] 3.829673e-13Notice that chi_yates matches R’s default output (\(52.729\)) exactly. If you disable the correction in R using chisq.test(tab_2x2, correct = FALSE), it will output chi_standard (\(55.194\)).
How to Report This in a Paper: > “A Pearson’s chi-square test with Yates’ continuity correction revealed a statistically significant association between treatment and pain response, \(\chi^2(1) = 52.73\), \(p < .001\). The association was strong, with a Phi coefficient (effect size) of \(\phi = \sqrt{52.73 / 180} = 0.54\).”
4.2.3 Step 3: Checking Assumptions (Expected Cell Counts)
The Chi-Square test is an approximation that assumes a large sample size. It becomes unreliable if the expected cell counts are too small. * The Rule of Thumb: At least 80% of expected cell counts must be \(\ge 5\), and no cell should have an expected count \(< 1\).
Let’s check the expected counts calculated by R:
#>
#> Non-Responder Responder
#> New 26.13333 69.86667
#> Standard 22.86667 61.13333Pedagogical note: R calculates expected values using the formula: \[\text{Expected Count} = \frac{\text{Row Total} \times \text{Column Total}}{\text{Grand Total}}\] For example, the expected count for (New, Non-Responder) is \((96 \times 49) / 180 = 26.13\). Since all our expected values are \(\ge 22.87\) (which is well above the threshold of \(5\)), the Chi-Square test assumptions are perfectly satisfied.
4.2.4 Step 4: Small Samples — Fisher’s Exact Test
If any expected cell count is \(< 5\), the Chi-Square test is invalid. In such cases, we use Fisher’s Exact Test, which calculates the exact hypergeometric probability of the table.
Let’s run Fisher’s test:
#>
#> Fisher's Exact Test for Count Data
#>
#> data: tab_2x2
#> p-value = 1.341e-14
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#> 0.009392081 0.115710445
#> sample estimates:
#> odds ratio
#> 0.03841056Notice that Fisher’s test also gives you: * Odds Ratio (OR): Under independence, \(\text{OR} = 1\). Here, the odds ratio is \(0.038\), which represents the odds of being a Non-Responder in the New treatment group compared to the Standard treatment group. Since the odds ratio is much less than \(1\), patients on the New treatment have a massive \(96.2\%\) reduction in the odds of being a Non-Responder (meaning they are dramatically more likely to be Responders!). * 95% Confidence Interval for the OR: Since this interval (\(0.009\) to \(0.116\)) does not contain \(1\) and is far below it, the association is highly statistically significant, matching the extremely small p-value (\(p = 1.34 \times 10^{-14}\)).
4.2.5 Step 5: Testing Larger Tables (\(r \times c\))
Let’s see what happens when we test a larger \(3 \times 3\) table (tab_3x3 representing Rater 1 vs Rater 2 Severity classifications):
#>
#> Pearson's Chi-squared test
#>
#> data: tab_3x3
#> X-squared = 189.53, df = 4, p-value < 2.2e-16#>
#> Mild Moderate Severe
#> Mild 23.11111 27.44444 14.444444
#> Moderate 27.02222 32.08889 16.888889
#> Severe 13.86667 16.46667 8.666667Note: In this case, all expected cell counts are above 5 (ranging from 8.67 to 32.09), so the Chi-Square approximation is technically valid. We show Fisher’s exact test here as a demonstration of how it can be applied to larger tables, or as a robustness check.
Because we want to demonstrate the method, we can run Fisher’s Exact Test:
#>
#> Fisher's Exact Test for Count Data
#>
#> data: tab_3x3
#> p-value < 2.2e-16
#> alternative hypothesis: two.sidedExercise: Test on InsectSprays
Using R’s built-in InsectSprays dataset:
1. Re-categorize the insect counts into a binary column: InsectSprays$high_count <- ifelse(InsectSprays$count > 10, "High", "Low").
2. Generate a table of spray (type of spray) × high_count.
3. Check the expected counts for this table. Can you use a Chi-Square test, or should you run Fisher’s exact test? Execute the appropriate test.
4.3 Trend Tests for Ordinal Data
When your exposure variable is ordinal (has a natural order), a trend test uses that order and is more powerful than a generic chi-square test.
#>
#> Non-Responder Responder
#> High 6 31
#> Low 19 48
#> Moderate 24 52Using base R’s prop.trend.test() (a Cochran-Armitage-style trend test built into R):
# prop.trend.test needs: number of "successes",
# total n, and a numeric score per group
successes <- trend_tab[, "Responder"]
totals <- rowSums(trend_tab)
# Low=1, Moderate=2, High=3
prop.trend.test(successes, totals, score = c(1, 2, 3))#>
#> Chi-squared Test for Trend in Proportions
#>
#> data: successes out of totals ,
#> using scores: 1 2 3
#> X-squared = 2.6318, df = 1, p-value = 0.1047Using the coin package for a linear-by-linear association test:
library(coin)
course_df$activity_level <- factor(
course_df$activity_level,
levels = c("Low", "Moderate", "High"),
ordered = TRUE
)
lbl_test(table(course_df$activity_level, course_df$response))#>
#> Asymptotic Linear-by-Linear Association Test
#>
#> data: Var2 by Var1 (Low < Moderate < High)
#> Z = -1.0883, p-value = 0.2765
#> alternative hypothesis: two.sidedBoth trend tests above assume the categories are genuinely ordered and roughly evenly spaced in the underlying construct they represent. If your “ordinal” categories are really just unordered labels with no natural progression, a trend test doesn’t make sense — use a regular chi-square test instead.
Compare the p-value from prop.trend.test() above to a plain chisq.test(trend_tab). Which is smaller? Given that activity_level genuinely has an order, which result should you trust more?
4.4 Agreement Measures: Kappa
When two raters independently classify the same subjects, percent agreement overstates how good the agreement really is, because some agreement happens by chance alone. Cohen’s kappa corrects for this.
Check Before You Run This (Cohen’s Kappa Diagnostics) Before running a Kappa test: 1. Identical Categories: Both raters must use the exact same set of mutually exclusive, categorical rating levels (e.g., both use Mild/Moderate/Severe). 2. Paired Observations: The subjects being rated must be identical, and each subject must be rated by both raters. 3. Independent Raters: The two raters must perform their ratings independently of each other.
library(psych)
kappa_result <- cohen.kappa(cbind(
course_df$rater1_severity,
course_df$rater2_severity
))
kappa_result#> Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels,
#> w.exp = w.exp)
#>
#> Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries
#> lower estimate upper
#> unweighted kappa 0.63 0.72 0.80
#> weighted kappa 0.67 0.76 0.86
#>
#> Number of subjects = 180Since severity is ordinal, the weighted kappa (which penalizes a “Mild vs. Severe” disagreement more than a “Mild vs. Moderate” disagreement) is more appropriate.
Under the hood:
* Unweighted Kappa: Treats all off-diagonal disagreements equally.
* Weighted Kappa: Assigns weights to off-diagonal cells depending on their distance from the diagonal (representing agreement).
- Linear weights (the default or weight = "linear"): Penalty increases linearly with the number of categories separating the ratings.
- Quadratic weights (often preferred in medical literature): Penalty increases quadratically with the distance, placing much heavier penalties on extreme disagreements.
Let’s inspect the weighted kappa output:
#> [1] 0.7612732The Math Behind It: Manual Kappa Calculation Cohen’s Kappa (\(\kappa\)) is calculated as: \[\kappa = \frac{p_o - p_e}{1 - p_e}\] where \(p_o\) represents the observed proportion of agreement, and \(p_e\) represents the expected proportion of agreement by chance alone.
Let’s calculate this manually in R to verify it matches cohen.kappa():
# 1. Create a table of the two raters
rater_a <- factor(c("Mild","Mild","Moderate","Severe","Severe","Moderate"))
rater_b_bad <- factor(c("Severe","Mild","Mild","Moderate","Mild","Severe"))
tab <- table(rater_a, rater_b_bad)
# 2. Observed agreement (sum of the diagonal cells / grand total)
po <- sum(diag(tab)) / sum(tab)
po#> [1] 0.1666667# 3. Expected agreement (sum of row probability * col probability)
row_prob <- rowSums(tab) / sum(tab)
col_prob <- colSums(tab) / sum(tab)
pe <- sum(row_prob * col_prob)
pe#> [1] 0.3333333#> [1] -0.25This matches R’s cohen.kappa(cbind(rater_a, rater_b_bad))$kappa output (\(-0.25\)) exactly!
Landis & Koch benchmarks for interpreting kappa: < 0 Poor, 0.00–0.20 Slight, 0.21–0.40 Fair, 0.41–0.60 Moderate, 0.61–0.80 Substantial, 0.81–1.00 Almost Perfect. These are widely used rules of thumb, not universal laws — always report the actual value, not just the label.
How to Report This in a Paper: > “Inter-rater reliability between the two clinical raters on patient symptom severity was substantial (unweighted Cohen’s \(\kappa = 0.72\), quadratic weighted \(\kappa = 0.76\)).”
Second example — perfect agreement vs. poor agreement:
# identical ratings:
rater_b_good <- factor(c(
"Mild","Mild","Moderate","Severe","Severe","Moderate"
))
cohen.kappa(cbind(rater_a, rater_b_good))$kappa#> [1] 1#> [1] -0.25Simulate your own pair of raters with sample(), deliberately making them agree about 90% of the time (hint: use ifelse(runif(n) < 0.9, rater_a, sample(...)) as in the dataset generator). Compute kappa and confirm it lands in the “Substantial” or “Almost Perfect” range.
4.5 Paired Categorical Tests: McNemar and Cochran’s Q
4.5.1 McNemar’s test — two paired binary measurements
#>
#> Absent Present
#> Absent 24 20
#> Present 79 57#>
#> McNemar's Chi-squared test with continuity correction
#>
#> data: mcnemar_tab
#> McNemar's chi-squared = 33.98, df = 1, p-value = 5.569e-09McNemar’s test ignores the patients whose status didn’t change (Present→Present or Absent→Absent) and focuses entirely on the discordant pairs — patients who changed status in one direction or the other. It tests whether those changes are symmetric (equally likely in both directions) or not.
4.5.2 Cochran’s Q test — generalizing McNemar to 3+ repeated measurements
Cochran’s Q isn’t available in a commonly pre-packaged apt-installed library, but the formula is simple enough to implement directly — which is also a great way to really understand what it’s doing:
cochran_q <- function(x) {
# x: matrix, rows = subjects,
# columns = repeated binary measurements (0/1)
x <- as.matrix(x)
k <- ncol(x)
Ri <- rowSums(x)
Cj <- colSums(x)
Q <- (k - 1) * (k * sum(Cj^2) - sum(Cj)^2) /
(k * sum(Ri) - sum(Ri^2))
list(
statistic = Q,
df = k - 1,
p.value = pchisq(Q, k - 1, lower.tail = FALSE)
)
}
# Simulate three repeated binary measurements per patient
# (e.g., symptom present at 3 visits)
set.seed(42)
n_pat <- 40
visit1 <- rbinom(n_pat, 1, 0.6)
visit2 <- rbinom(n_pat, 1, 0.45)
visit3 <- rbinom(n_pat, 1, 0.30)
repeated_symptoms <- cbind(visit1, visit2, visit3)
cochran_q(repeated_symptoms)#> $statistic
#> [1] 2.066667
#>
#> $df
#> [1] 2
#>
#> $p.value
#> [1] 0.3558189If you’d rather use a maintained package for Cochran’s Q, RVAideMemoire::cochran.qtest() provides one — install it with install.packages("RVAideMemoire"). We implemented it manually above so this book keeps working even without every optional package installed, and so you can see precisely how the test statistic is built.
4.6 Controlling for a Stratifying Variable: The CMH Test
A treatment-response association might look different — or vanish, or even reverse — once you account for which hospital site a patient came from. The Cochran-Mantel-Haenszel test checks the association while controlling for a stratifying variable.
#>
#> Pearson's Chi-squared test with Yates' continuity correction
#>
#> data: table(course_df$treatment, course_df$response)
#> X-squared = 52.729, df = 1, p-value = 3.83e-13# Stratified by site
cmh_array <- table(course_df$treatment, course_df$response, course_df$site)
mantelhaen.test(cmh_array)#>
#> Mantel-Haenszel chi-squared test with continuity correction
#>
#> data: cmh_array
#> Mantel-Haenszel X-squared = 51.393, df = 1, p-value = 7.559e-13
#> alternative hypothesis: true common odds ratio is not equal to 1
#> 95 percent confidence interval:
#> 0.01378706 0.12237804
#> sample estimates:
#> common odds ratio
#> 0.04107595If the stratified (CMH) result tells a meaningfully different story than the unstratified chi-square result, site is likely a confounder — a variable associated with both the exposure and outcome that was distorting the raw, unstratified relationship. This is one of the most important ideas in observational data analysis.
Look at the individual 2×2 tables inside cmh_array (try cmh_array[,, "Site A"], and so on for the other sites). Does the treatment-response relationship look consistent across all three sites, or does it vary a lot? What does that suggest about whether CMH pooling is appropriate here?
4.7 Compact Reporting: Tables That Publish
4.7.1 A fully-manual baseline table (always works, no extra packages needed)
build_baseline_table <- function(data, vars, strata) {
groups <- unique(data[[strata]])
result <- data.frame(Variable = vars)
for (g in groups) {
subset_data <- data[data[[strata]] == g, ]
result[[g]] <- sapply(vars, function(v) {
x <- subset_data[[v]]
if (is.numeric(x)) {
sprintf("%.1f (%.1f)", mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
} else {
paste0(names(sort(table(x), decreasing = TRUE))[1], " most common")
}
})
}
result
}
build_baseline_table(course_df, vars = c("age", "sex"), strata = "treatment")This is exactly what packages like tableone and gtsummary automate for you — computing mean (SD) for numeric variables and counts/percentages for categorical ones, split by group. Understanding the manual version means you’ll never be stuck if a package isn’t available, and you’ll better understand what those packages’ output actually means.
4.7.2 Using tableone and gtsummary for Professional Reports
In practice, medical and statistical researchers do not compile these tables manually. They use packages like tableone or gtsummary which automate this process and generate beautifully formatted tables with a single line of code.
Here is how you generate a baseline table using the tableone package:
library(tableone)
# Generate a Table 1 object
t1 <- CreateTableOne(vars = c("age", "sex", "baseline_pain"),
strata = "treatment", data = course_df)
# Print it (renders in standard markdown/text format)
print(t1)#> Stratified by treatment
#> New Standard p test
#> n 96 84
#> age (mean (SD)) 53.74 (13.29) 50.98 (11.72) 0.143
#> sex = Male (%) 46 (47.9) 39 (46.4) 0.960
#> baseline_pain (mean (SD)) 6.19 (2.09) 6.37 (2.13) 0.564Here is how you generate a highly polished HTML baseline table using the gtsummary package:
library(gtsummary)
course_df |>
tbl_summary(by = treatment, include = c(age, sex, baseline_pain)) |>
add_p()| Characteristic | New N = 961 |
Standard N = 841 |
p-value |
|---|---|---|---|
| age | 56 (44, 65) | 52 (43, 58) | |
| sex | |||
| Female | 50 (52%) | 45 (54%) | |
| Male | 46 (48%) | 39 (46%) | |
| baseline_pain | 6.10 (4.40, 7.70) | 6.45 (4.55, 7.80) | |
| 1 Median (Q1, Q3); n (%) | |||
Why do the p-values differ between tableone and gtsummary?
A student comparing the two tables above side-by-side will notice that for age, the tableone package reports a p-value of \(0.143\), whereas gtsummary reports \(0.083\). This is not an error!
* tableone defaults to continuous variables being summarized by mean (SD) and compared using a parametric t-test (or ANOVA).
* gtsummary defaults to continuous variables being summarized by median (IQR) and compared using a non-parametric Wilcoxon rank-sum test (or Kruskal-Wallis).
Always check the defaults of the reporting packages you use!
Extend build_baseline_table() above to also report a percentage breakdown for categorical variables (currently it just reports the most common category, which loses information). This is a genuinely useful exercise — you’ll end up with something close to what tableone does internally.
table()andxtabs()build contingency tables;addmargins()adds totals.- Fisher’s exact test is safer than chi-square when expected counts are small; check
chisq.test(...)$expectedto decide. - Trend tests (
prop.trend.test(),coin::lbl_test()) use the natural ordering of an ordinal exposure for more statistical power. - Cohen’s kappa corrects raw percent agreement for chance — always prefer weighted kappa for ordinal categories.
- McNemar’s test handles paired binary data; Cochran’s Q generalizes it to 3+ repeated binary measurements.
- The CMH test lets you check an association while controlling for a stratifying variable, revealing potential confounding.
- Building a manual summary table once is worth doing — it demystifies what packages like
tableoneandgtsummaryautomate.
Next, in Chapter 5, we turn to hypothesis testing for continuous or ordinal outcomes when parametric assumptions (like normality) don’t hold.