Selecting log-linear models for a contingency table

Bayesian variable selection, Part 4 — Albert (1996)

Variable selection for categorical data means choosing which associations belong in a log-linear model for a contingency table. 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 their marginal likelihoods — so the data, not the prior's tails, decide which associations are real.

The showcase is the Berkeley graduate-admissions table (UCBAdmissions): 4526 applicants cross-classified by admission (admitted/rejected), 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? This is Simpson's paradox, posed as model selection.

The model and the method

The 24 cell counts $y$ are modelled as independent Poisson with a log-linear mean. For factors $A$ (admit), $G$ (gender), $D$ (department): $$ \log \mathbb E[y] = \mu + A + G + D + \underbrace{(AG) + (AD) + (GD)}_{\text{pairwise associations}} + \underbrace{(AGD)}_{\text{three-way}}. $$ Main effects are always present (they only fix the marginal totals); 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 interesting contrast is $[AD][GD]$ — admission and gender conditionally independent given department — versus $[AG][AD][GD]$, which adds a direct admission-gender link.

Marginal likelihood by Laplace. Each model's evidence is $$ \log p(y\mid\text{model}) \approx \ell(\hat\theta) + \log \pi(\hat\theta) + \tfrac d2\log 2\pi - \tfrac12\log|H|,\qquad H = X'\,\mathrm{diag}(\hat\mu)\,X, $$ with $\hat\theta$ the Poisson-GLM fit and $H$ its Fisher information (the mode $\approx$ MLE for a table this well-populated).

The prior (Albert's robustness). Intercept and main effects get a diffuse normal ($\tau=10$). The interaction coefficients — the ones being tested — get a heavy-tailed Cauchy$(0,s)$ prior, $s=1$: robust, so a genuine association is not shrunk away and a spurious one is not kept on the strength of prior tails. A uniform prior is placed over the nine models; posterior model probabilities and association inclusion probabilities (summing model probabilities over models containing each term) follow. Engine: loglinsel.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import loglinsel as L                                               # from-scratch log-linear model selection

df = pd.read_csv("loglin_data.csv")                                 # UCBAdmissions as a long table
R, incl = L.loglin_select(df, "Freq", ["Admit","Gender","Dept"], s_cauchy=1.0)

print("log-linear models ranked by posterior probability:")
for _, r in R.head(3).iterrows():
    nm = " + ".join(t.replace(":", ":") for t in r["terms"]) or "independence"
    print("  %-40s k=%d   %.3f" % (nm if r["model"]!="saturated" else "saturated (adds the 3-way term)", r["k"], r["postprob"]))
print("\nassociation inclusion probabilities:")
for a, p in incl.items(): print("  %-18s %.3f" % (a, p))

# --- Fig: Simpson's paradox -- the aggregate gap vs the within-department picture ---
depts = ["A", "B", "C", "D", "E", "F"]
adm = df.pivot_table(index=["Gender", "Dept"], columns="Admit", values="Freq", aggfunc="sum")
rate = adm["Admitted"] / (adm["Admitted"] + adm["Rejected"])       # admission rate by gender x dept
agg = df.pivot_table(index="Gender", columns="Admit", values="Freq", aggfunc="sum")
agg_rate = agg["Admitted"] / (agg["Admitted"] + agg["Rejected"])   # aggregate admission rate by gender

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12.4, 4.4))
# left: the raw aggregate gap
ax0.bar(["Male", "Female"], [agg_rate["Male"], agg_rate["Female"]],
        color=["firebrick", "navy"], width=.6)
for i, g in enumerate(["Male", "Female"]):
    ax0.text(i, agg_rate[g] + .01, "%.0f%%" % (100 * agg_rate[g]), ha="center")
ax0.set_ylim(0, .5); ax0.set_ylabel("admission rate")
ax0.set_title("Aggregate: men admitted at a higher rate")
ax0.grid(alpha=.25, axis="y")
# right: within each department the gap disappears
ax1.plot(depts, [rate.loc[("Male", d)] for d in depts], "o-", color="firebrick", label="Male")
ax1.plot(depts, [rate.loc[("Female", d)] for d in depts], "s-", color="navy", label="Female")
ax1.set_ylim(0, 1); ax1.set_xlabel("department"); ax1.set_ylabel("admission rate")
ax1.set_title("Within department: rates track each other")
ax1.grid(alpha=.25); ax1.legend()
plt.tight_layout(); plt.show()

log-linear models ranked by posterior probability:
  Admit:Dept + Gender:Dept                 k=18   0.904
  Admit:Gender + Admit:Dept + Gender:Dept  k=19   0.094
  saturated (adds the 3-way term)          k=24   0.002

association inclusion probabilities:
  Admit:Gender       0.096
  Admit:Dept         1.000
  Gender:Dept        1.000
  Admit:Gender:Dept  0.002

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\regression\_tools.py:121: RuntimeWarning: divide by zero encountered in scalar divide
  scale = np.dot(wresid, wresid) / df_resid
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\genmod\generalized_linear_model.py:1342: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified
  warnings.warn(msg, category=PerfectSeparationWarning)

Reading the selection

• The conditional-independence model wins outright. $[AD][GD]$ — admission depends on department, gender-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 (the next model) actually lowers the evidence.
• The associations, quantified. Admit:Dept and Gender:Dept have inclusion probability 1.00 — departments differ sharply in how selective they are, and men and women apply to different departments. But Admit:Gender sits at 0.10: once department is in the model, the data give little support for a direct gender effect on admission. The three-way term is essentially absent (0.002).
• Simpson's paradox, resolved (figure). The left panel shows the raw aggregate gap; the right shows that within each department admission rates are similar across genders — women simply applied more to the competitive departments. The log-linear selection encodes exactly this: drop $AG$, keep $AD$ and $GD$.

import pymc as pm, patsy                                            # cross-check: group-SSVS on the interaction BLOCKS
y = df["Freq"].values.astype(float)
Xm = np.asarray(patsy.dmatrix("C(Admit)+C(Gender)+C(Dept)", df))     # main effects (always in)
full = patsy.dmatrix("C(Admit)*C(Gender)*C(Dept)", df, return_type="dataframe"); slc = full.design_info.term_slices
blk = lambda key: full.values[:, [s for t, s in slc.items() if t.name()==key][0]]
XAG, XAD, XGD = blk("C(Admit):C(Gender)"), blk("C(Admit):C(Dept)"), blk("C(Gender):C(Dept)")
with pm.Model() as mod:
    tau, c = 0.05, 30.0
    eta = pm.math.dot(Xm, pm.Normal("bm", 0, 10, shape=Xm.shape[1]))
    gg = {}
    for nm, XX in [("AG", XAG), ("AD", XAD), ("GD", XGD)]:          # one on/off switch per association block
        g = pm.Bernoulli("g_"+nm, 0.5); b = pm.Normal("b_"+nm, 0, 1, shape=XX.shape[1])
        eta = eta + pm.math.dot(XX, b*pm.math.switch(g, c*tau, tau)); gg[nm] = g
    pm.Poisson("y", pm.math.exp(eta), observed=y)
    idata = pm.sample(3000, tune=2000, chains=4, step=[pm.BinaryGibbsMetropolis(list(gg.values()))],
                      target_accept=0.9, random_seed=1, progressbar=False)
ssvs = {n: float(idata.posterior['g_'+n].mean()) for n in ["AG", "AD", "GD"]}
print("pymc group-SSVS block inclusion:  ", "   ".join("%s=%.3f" % (n, ssvs[n]) for n in ["AG","AD","GD"]))
print("-> both engines agree: the direct gender-admission association is NOT supported;\n   department alone explains the aggregate admission gap (Simpson's paradox).")

# --- Fig: association inclusion probabilities, enumeration vs group-SSVS ---
assoc = ["Admit:Gender", "Admit:Dept", "Gender:Dept"]
enum_p = [incl[a] for a in assoc]                                   # exact enumeration + Laplace (cell above)
ssvs_p = [ssvs[k] for k in ["AG", "AD", "GD"]]                      # PyMC group-SSVS
x = np.arange(len(assoc)); w = .38
fig, ax = plt.subplots(figsize=(8.0, 4.4))
ax.bar(x - w/2, enum_p, w, color="firebrick", label="enumeration (Laplace)")
ax.bar(x + w/2, ssvs_p, w, color="navy", label="group-SSVS (PyMC)")
ax.axhline(.5, color="0.4", ls="--", lw=1)                          # inclusion decision threshold
ax.set_xticks(x); ax.set_xticklabels(assoc)
ax.set_ylim(0, 1.08); ax.set_ylabel("inclusion probability")
ax.set_title("Association inclusion probabilities: two engines agree")
ax.grid(alpha=.25, axis="y"); ax.legend()
plt.tight_layout(); plt.show()

Multiprocess sampling (4 chains in 4 jobs)

CompoundStep

>BinaryGibbsMetropolis: [g_AG, g_AD, g_GD]

>NUTS: [bm, b_AG, b_AD, b_GD]

Sampling 4 chains for 2_000 tune and 3_000 draw iterations (8_000 + 12_000 draws total) took 11 seconds.

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\arviz_stats\base\diagnostics.py:90: RuntimeWarning: invalid value encountered in scalar divide
  (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)
There were 538 divergences after tuning. Increase `target_accept` or reparameterize.

pymc group-SSVS block inclusion:   AG=0.049   AD=1.000   GD=1.000
-> both engines agree: the direct gender-admission association is NOT supported;
   department alone explains the aggregate admission gap (Simpson's paradox).

Results — selection answers the causal-looking question

• Two engines, one verdict. Exact enumeration with Laplace marginal likelihoods and a PyMC group-SSVS (an on/off indicator per interaction block) agree: Admit:Dept and Gender:Dept are certain, Admit:Gender is not supported (~0.05–0.10). The apparent admissions bias is an artefact of aggregation across departments. (The discrete on/off indicators make the PyMC sampler mix poorly — divergences, and an undefined $\hat R$ for the indicators pinned at 1 — but the exact enumeration is the primary result and the SSVS agrees with it.)
• Robust by design. The Cauchy prior on the interaction terms means the near-zero inclusion of Admit:Gender is a statement about the data, not an artefact of a tight prior shrinking it — Albert's point. A genuine association of that size would still have been detected.
• Variable selection, categorical-data style. Here "which predictors matter" becomes "which associations are present in the table" — the same spike-and-slab / marginal-likelihood logic as Parts 1–3, now over the interaction terms of a log-linear model. It is the natural selection tool for the contingency-table models that fill Congdon's book.

References

• Albert, J. H. (1996). The Bayesian selection of log-linear models. Canadian Journal of Statistics 24, 327–347.
• Bickel, P. J., Hammel, E. A. & O'Connell, J. W. (1975). Sex bias in graduate admissions: data from Berkeley. Science 187, 398–404.
• Congdon, P. (2005). Bayesian Models for Categorical Data, Ch. 7 (log-linear models). Wiley.

Data: loglin_data.csv — the Berkeley UCBAdmissions table: admission × gender × department, 4526 applicants. Next: loglin_R.ipynb — the R cross-check (glm Poisson models over the hierarchy, ranked by BIC).