2. Example 1 — Chi-Square Test of Independence
Research question: Is the presence of Batrachochytrium dendrobatidis (Bd), the chytrid fungus driving global amphibian declines, associated with elevation zone in the Sierra Nevada?
Data: 240 swab samples from three elevation zones (low, mid, high) scored as Bd-positive or Bd-negative.
The hypotheses
The chi-square statistic
The test compares observed counts in each cell of a contingency table to expected counts — the counts we would predict if the null were true (i.e., if infection rate were uniform across elevations).
\[\chi^2 = \sum_{\text{all cells}} \frac{(O_i - E_i)^2}{E_i}\]
Expected count for cell (i, j): \(E_{ij} = \dfrac{(\text{row total}_i)(\text{col total}_j)}{n}\)
The statistic follows a χ² distribution with \((r-1)(c-1)\) degrees of freedom, where \(r\) = number of rows and \(c\) = number of columns.
⚠️ The expected-count rule
The χ² approximation is reliable only when every expected cell count is ≥ 5. If any expected cell is < 5, use Fisher’s exact test (see Example 3) or collapse categories if scientifically defensible.
R walkthrough
set.seed(101)
# Simulate Bd survey data
elevation <- rep(c("Low", "Mid", "High"), times = c(80, 80, 80))
p_infect <- c(Low = 0.55, Mid = 0.35, High = 0.20)
infected <- rbinom(240, 1,
prob = p_infect[elevation]) == 1
bd <- data.frame(
elevation = factor(elevation, levels = c("Low", "Mid", "High")),
infected = factor(ifelse(infected, "Positive", "Negative"))
)
# Observed contingency table
obs_table <- table(bd$elevation, bd$infected)
addmargins(obs_table)##
## Negative Positive Sum
## Low 42 38 80
## Mid 51 29 80
## High 61 19 80
## Sum 154 86 240# Expected counts under H₀
chisq_result <- chisq.test(obs_table)
round(chisq_result$expected, 1)##
## Negative Positive
## Low 51.3 28.7
## Mid 51.3 28.7
## High 51.3 28.7chisq_result##
## Pearson's Chi-squared test
##
## data: obs_table
## X-squared = 9.8218, df = 2, p-value = 0.007366# Summarise infection rates per zone
bd_summary <- bd |>
group_by(elevation) |>
summarise(
n = n(),
pos = sum(infected == "Positive"),
rate = pos / n,
se = sqrt(rate * (1 - rate) / n),
lo = pmax(0, rate - 1.96 * se),
hi = pmin(1, rate + 1.96 * se),
.groups = "drop"
)
par(mar = c(4, 5, 3, 2), las = 1)
bp <- barplot(
bd_summary$rate,
names.arg = levels(bd_summary$elevation),
col = c("#c8832a", "#4a7c5e", "#1e3a2f"),
border = "white",
ylim = c(0, 0.8),
ylab = "Proportion Bd-positive",
main = expression(italic("Bd") ~ "infection rate by elevation zone"),
cex.main = 0.92
)
arrows(bp, bd_summary$lo, bp, bd_summary$hi,
angle = 90, code = 3, length = 0.06, col = "#3a3a38", lwd = 1.4)
abline(h = mean(bd_summary$rate), lty = 2, col = "#7a7a72", lwd = 1.2)
text(bp[3] + 0.5, mean(bd_summary$rate) + 0.03,
"Overall mean", col = "#7a7a72", cex = 0.76, adj = 0)
Bd infection rate by elevation zone. Error bars show 95% binomial confidence intervals. The declining gradient from low to high elevation is the pattern driving the significant chi-square result.
Interpreting the output
The χ² statistic summarises how much the observed counts deviate from what perfect independence would predict, scaled by the expected counts. A large χ² means the data are inconsistent with independence.
Effect size — Cramér’s V:
\[V = \sqrt{\frac{\chi^2}{n \cdot (\min(r,c) - 1)}}\]
Cramér’s V ranges from 0 (no association) to 1 (perfect association). Benchmarks: V ≈ 0.10 = small, V ≈ 0.30 = medium, V ≈ 0.50 = large.
# Cramér's V
n <- sum(obs_table)
k <- min(nrow(obs_table), ncol(obs_table))
V <- sqrt(chisq_result$statistic / (n * (k - 1)))
cat(sprintf("Cramér's V = %.3f\n", V))## Cramér's V = 0.202💡 Reporting a chi-square result
“Bd infection status was significantly associated with elevation zone, χ²(2) = , p = , Cramér’s V = ___. Infection rates declined from 48% at low elevation to 24% at high elevation.”
Always state: test statistic, degrees of freedom, p-value, effect size, and a plain-language description of the direction of the association.
