Bayesian Selection of Log-Linear Models

Python · R · Part 4 of the variable-selection arc — categorical data

Overview

Part 4 of the variable-selection arc: selection for categorical data. Variable selection for a contingency table means choosing which associations belong in a log-linear model. Albert (1996) brought the Bayesian selection machinery to this setting — model the cell counts as Poisson with a log-linear mean, place robust heavy-tailed priors on the interaction terms, and compare models by marginal likelihood, so the data (not the prior's tails) decide which associations are real. The showcase is the Berkeley graduate-admissions table (4,526 applicants cross-classified by admission, gender, and department A–F): aggregated, men are admitted at a higher rate than women — apparent bias. The log-linear selection asks the precise question — is there a direct admission–gender association, once department is accounted for? — posing Simpson's paradox as model selection.

\log \mathbb E[y] = \mu + A + G + D + \underbrace{(AG) + (AD) + (GD)}_{\text{pairwise associations}} + \underbrace{(AGD)}_{\text{three-way}}

Poisson log-linear models, robust priors, Laplace evidence

The 24 cell counts are independent Poisson with a log-linear mean over the factors admit ( $A$ ), gender ( $G$ ), department ( $D$ ): main effects plus the pairwise associations $(AG),(AD),(GD)$ and the three-way $(AGD)$ . Main effects are always present — they only fix the marginal totals — so selection is over the interaction terms, which are the associations. Keeping models hierarchical (a term enters only with all its lower-order relatives) leaves nine candidates, from mutual independence $[A][G][D]$ to the saturated $[AGD]$ ; the key contrast is $[AD][GD]$ — admission and gender conditionally independent given department — versus $[AG][AD][GD]$ , which adds a direct admission–gender link. Each model's evidence is computed by a Laplace approximation around the Poisson-GLM fit ( $\log|H|$ from the Fisher information $H=X'\,\mathrm{diag}(\hat\mu)\,X$ ). Albert's robustness comes from the prior: diffuse normals on the intercept and main effects, but a heavy-tailed Cauchy(0, 1) on the interaction coefficients being tested — so a genuine association is not shrunk away, and a spurious one is not kept on the strength of prior tails. A uniform prior over the nine models gives posterior model probabilities and association inclusion probabilities. Engine: a self-contained loglinsel.py.

Results

The conditional-independence model wins outright. $[AD][GD]$ — admission depends on department, application patterns depend on department, but admission and gender are conditionally independent given department — carries 90% of the posterior; adding the direct admission–gender term actually lowers the evidence. Quantified: Admit:Dept and Gender:Dept have inclusion probability 1.00 (departments differ sharply in selectivity, and men and women apply to different departments), but Admit:Gender sits at 0.10 and the three-way term is essentially absent (0.002). Simpson's paradox, resolved: the raw aggregate gap disappears within department — women simply applied more to the competitive departments — and the selection encodes exactly this (drop $AG$ , keep $AD$ and $GD$ ). The near-zero inclusion of Admit:Gender is robust by design: the Cauchy prior means it is a statement about the data, not a tight prior shrinking it — a genuine association of that size would still have been detected. Exact enumeration with Laplace marginal likelihoods and a PyMC group-SSVS (an on/off indicator per interaction block) agree. Here "which predictors matter" becomes "which associations are present in the table" — the same spike-and-slab / marginal-likelihood logic as Parts 1–3, over the interactions of a log-linear model.

Notebooks

Python Log-linear model selection from scratch — Laplace marginal likelihoods with robust Cauchy interaction priors over the model hierarchy, plus a PyMC group-SSVS cross-check, on the Berkeley admissions table (Python) Open → R R cross-check — the competing hierarchical log-linear models as Poisson GLMs ranked by BIC, recovering the same conditional-independence conclusion (R) Open →

Downloads

loglinsel.py The from-scratch log-linear model-selection engine — Poisson-GLM fits, Laplace marginal likelihoods with Cauchy interaction priors, and inclusion probabilities over the model hierarchy (new self-contained module) loglin_data.csv The Berkeley UCBAdmissions table — admission × gender × department, 4,526 applicants (Bickel–Hammel–O'Connell 1975)

References

Albert, J. H. (1996). Bayesian selection of log-linear models. Canadian Journal of Statistics 24(3), 327–347. — the robust-prior log-linear model selection reproduced here
Bickel, P. J., Hammel, E. A. & O'Connell, J. W. (1975). Sex bias in graduate admissions: data from Berkeley. Science 187(4175), 398–404. — the Berkeley admissions table and the Simpson's-paradox case
Bishop, Y. M. M., Fienberg, S. E. & Holland, P. W. (1975). Discrete Multivariate Analysis: Theory and Practice. MIT Press. — hierarchical log-linear models for contingency tables
Kass, R. E. & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association 90(430), 773–795. — Laplace evidence, Bayes factors, and the BIC connection
Raftery, A. E. (1986). A note on Bayes factors for log-linear contingency table models with vague prior information. Journal of the Royal Statistical Society B 48(2), 249–250. — Bayes factors for log-linear models
Dellaportas, P. & Forster, J. J. (1999). Markov chain Monte Carlo model determination for hierarchical and graphical log-linear models. Biometrika 86(3), 615–633. — stochastic model search over the log-linear hierarchy