Log-linear model selection — R cross-check

Bayesian variable selection, Part 4 (R)

Validates loglin_python.ipynb. R's native tool for contingency tables is a Poisson glm, and the classical model-selection currency is BIC — which is itself a Laplace/unit-information approximation to (minus twice) the log marginal likelihood. So fitting the competing hierarchical log-linear models and ranking them by BIC should reproduce the Bayesian result: the conditional-independence model [Admit·Dept][Gender·Dept] wins, and adding a direct Admit:Gender term is not supported. We confirm the latter with an explicit likelihood-ratio test.

Data: loglin_data.csv — the Berkeley UCBAdmissions table (admission × gender × department, 4526 applicants).

.libPaths('C:/Users/user/R/win-library/4.6')
d <- read.csv('loglin_data.csv'); N <- sum(d$Freq)                                        # N = total count = 4526

# competing hierarchical log-linear models, fit as Poisson GLMs on the cell counts
mods <- list('[A][G][D]'    = glm(Freq ~ Admit + Gender + Dept, poisson, d),              # mutual independence
             '[AD][GD]'     = glm(Freq ~ Admit + Gender + Dept + Admit:Dept + Gender:Dept, poisson, d),
             '[AG][AD][GD]' = glm(Freq ~ (Admit + Gender + Dept)^2, poisson, d),          # + Admit:Gender
             'saturated'    = glm(Freq ~ Admit * Gender * Dept, poisson, d))

# BIC as a marginal-likelihood approximation.  For a contingency table the effective
# sample size is the TOTAL COUNT N, not the number of cells (R's default nobs=24) --
# the information about each parameter grows with N.  This matches the Bayesian marginal
# likelihood; the default cell-based penalty is too weak and over-rewards the saturated model.
bicN <- sapply(mods, function(m) -2*as.numeric(logLik(m)) + length(coef(m))*log(N))
pp   <- exp(-0.5 * (bicN - min(bicN))); pp <- pp / sum(pp)                                # -> posterior model prob
tab  <- data.frame(model=names(mods), k=sapply(mods, function(m) length(coef(m))),
                   BIC_N=round(bicN,1), postprob=round(pp,3))
print(tab, row.names=FALSE)
m_adgd <- mods[['[AD][GD]']]; m_all2 <- mods[['[AG][AD][GD]']]

cat('\nLikelihood-ratio test: does adding the Admit:Gender association help?\n')
print(anova(m_adgd, m_all2, test='Chisq'))
cat('\ncross-check -> Python: [AD][GD] wins; Admit:Gender not supported (inclusion ~0.10). Agree.\n')

options(repr.plot.width=7, repr.plot.height=4.2)
barplot(pp, names.arg=tab$model, col=c('grey70','firebrick','grey70','grey70'), ylim=c(0,1),
        ylab='posterior model probability (from BIC)', main='Berkeley admissions: conditional independence [AD][GD] wins', las=1)

        model  k  BIC_N postprob
    [A][G][D]  8 2324.1    0.000
     [AD][GD] 18  332.3    0.969
 [AG][AD][GD] 19  339.2    0.031
    saturated 24  361.1    0.000

Likelihood-ratio test: does adding the Admit:Gender association help?

Analysis of Deviance Table

Model 1: Freq ~ Admit + Gender + Dept + Admit:Dept + Gender:Dept
Model 2: Freq ~ (Admit + Gender + Dept)^2
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1         6     21.735                     
2         5     20.204  1   1.5312   0.2159

cross-check -> Python: [AD][GD] wins; Admit:Gender not supported (inclusion ~0.10). Agree.

Same conclusion, classical route

• BIC agrees with the Bayesian enumeration. The [AD][GD] model has the lowest BIC and almost all of the BIC-based posterior mass — matching the Cauchy-prior Laplace result in Python. Adding the Admit:Gender term raises BIC (a worse model).
• The likelihood-ratio test confirms it. Adding Admit:Gender to [AD][GD] gives a small, non-significant change in deviance — no evidence for a direct gender-admission association once department is controlled. The aggregate gap is Simpson's paradox, exactly as the inclusion probabilities said.

References

• Albert, J. H. (1996). The Bayesian selection of log-linear models. Canadian J. Statistics 24, 327–347.
• Bickel, P. J., Hammel, E. A. & O'Connell, J. W. (1975). Sex bias in graduate admissions. Science 187, 398–404.