Count Data Tests (Chi-Square)

Classical Test	bossanova Equivalent	Model Comparison
Chi-square goodness of fit	`compare(null, full)`	Intercept-only vs category model
Chi-square independence	`compare(main_effects, interaction)`	Main effects vs interaction
Multi-way tables	`model("count ~ a * b * c", df, family="poisson")`	Hierarchical model comparisons

Poisson GLM with log link provides a unified framework for all chi-square tests via likelihood ratio comparisons using compare().

import numpy as np
from scipy import stats
import polars as pl
from bossanova.model import model
from bossanova import load_dataset, compare

np.random.seed(42)

# Load penguins dataset
penguins = load_dataset("penguins").drop_nulls().filter(pl.col("sex") != "NA")

# Count penguins per species
species_counts = penguins.group_by("species").agg(pl.len().alias("count")).sort("species")
observed = species_counts["count"].to_numpy()
categories = species_counts["species"].to_list()

# Cross-tabulation: species × sex (complete 3×2 table)
contingency = penguins.group_by("species", "sex").agg(pl.len().alias("count"))
contingency_wide = contingency.pivot(on="sex", index="species", values="count").fill_null(0)
observed_table = contingency_wide.select(pl.exclude("species")).to_numpy()

Chi-Square Goodness of Fit¶

Classical:

\chi^2 = \sum_{i=1}^{k} \frac{(O_i - E_i)^2}{E_i}, \quad \chi^2 \sim \chi^2(k-1) \text{ under } H_0 \text{ (equal proportions)}

(1)

As GLM:

y_i \sim \text{Poisson}(\mu_i), \quad \log(\mu_i) = \beta_0 + \beta_1 x_{1i} + \cdots

(2)

G^2 = 2\sum_{i=1}^{k} O_i \log\!\left(\frac{O_i}{\hat{E}_i}\right) \dot{\sim} \chi^2(k-1) \text{ under } H_0: \beta_1 = \cdots = 0

(3)

The Pearson $\chi^2$ and the likelihood ratio $G^2$ are asymptotically equivalent. bossanova reports $G^2$ via compare().

scipy¶

from scipy.stats import chisquare

scipy_chi2 = chisquare(observed)
scipy_chi2

Power_divergenceResult(statistic=np.float64(28.270270270270274), pvalue=np.float64(7.264217011785267e-07))

bossanova¶

df = pl.DataFrame({"count": observed, "species": categories})

m_null = model("count ~ 1", df, family="poisson").fit()
m_full = model("count ~ species", df, family="poisson").fit()

compare(m_null, m_full)

result = compare(m_null, m_full)
lr_stat = float(result.filter(pl.col("chi2").is_not_null())["chi2"][0])
lr_pvalue = float(result.filter(pl.col("p_value").is_not_null())["p_value"][0])
assert np.isclose(lr_stat, scipy_chi2.statistic, rtol=0.5), f"stat mismatch: {lr_stat} vs {scipy_chi2.statistic}"
# LR and Pearson chi-square p-values converge asymptotically; for small counts they can differ
# but both should agree on significance
assert (lr_pvalue < 0.05) == (scipy_chi2.pvalue < 0.05), f"Significance mismatch: {lr_pvalue} vs {scipy_chi2.pvalue}"

Chi-Square Test of Independence¶

Classical:

\chi^2 = \sum_{i=1}^{r} \sum_{j=1}^{c} \frac{(O_{ij} - E_{ij})^2}{E_{ij}}, \quad \chi^2 \sim \chi^2((r-1)(c-1)) \text{ under } H_0 \text{ (independence)}

(4)

As GLM:

y_{ij} \sim \text{Poisson}(\mu_{ij}), \quad \log(\mu_{ij}) = \beta_0 + \alpha_i + \gamma_j + (\alpha\gamma)_{ij}

(5)

Independence means no interaction term. The LRT comparing the main-effects model ( $\log\mu_{ij} = \beta_0 + \alpha_i + \gamma_j$ ) to the saturated model tests $H_0: (\alpha\gamma)_{ij} = 0$ for all $i,j$ .

scipy¶

from scipy.stats import chi2_contingency

scipy_indep = chi2_contingency(observed_table)
pl.DataFrame({
    "statistic": [scipy_indep.statistic],
    "p_value": [scipy_indep.pvalue],
    "df": [scipy_indep.dof]
})

bossanova¶

# Null model (independence = main effects only)
m_null = model("count ~ species + sex", contingency, family="poisson").fit()
# Full model (saturated, with interaction)
m_full = model("count ~ species * sex", contingency, family="poisson").fit()

compare(m_null, m_full)

result = compare(m_null, m_full)
lr_stat = float(result.filter(pl.col("chi2").is_not_null())["chi2"][0])
lr_pvalue = float(result.filter(pl.col("p_value").is_not_null())["p_value"][0])
assert np.isclose(lr_stat, scipy_indep.statistic, rtol=0.5), f"stat mismatch: {lr_stat} vs {scipy_indep.statistic}"
assert np.isclose(lr_pvalue, scipy_indep.pvalue, rtol=0.5), f"p mismatch: {lr_pvalue} vs {scipy_indep.pvalue}"

Multi-Factor Log-Linear Models¶

For tables with three or more factors, the GLM framework extends naturally—there is no simple classical test equivalent.

bossanova¶

# Three-way table: species × island × sex
counts_3way = penguins.group_by("species", "island", "sex").agg(pl.len().alias("count"))

m = model("count ~ species + island + sex", counts_3way, family="poisson").fit().infer()

m.params.select("term", "estimate", "p_value")