Bayesian and Generalised Regression Models

1. Bayesian Inference — A New Way to Think About Uncertainty

Research Question

"Given the data I observe, what is the probability that the return to education is greater than 8%?"

Intuition

Frequentist inference asks: "If the true parameter were zero, how unlikely is my data?" — this is the p-value. Bayesian inference asks the more natural question: "Given my data, what are the plausible values of the parameter?"

Bayes' theorem provides the mathematical machinery:

P(θ | data) ∝ P(data | θ) × P(θ)
Posterior ∝ Likelihood × Prior

The prior P(θ) encodes what you believe before seeing the data. The likelihood P(data | θ) is the same likelihood function used in MLE. The posterior P(θ | data) combines both — it is your updated belief after observing the data.

The critical insight: the posterior is a full probability distribution, not a point estimate. Instead of a single "β = 0.107 with SE = 0.014," you get an entire distribution. You can directly answer: "There is a 93% probability that β > 0.08" — a statement a frequentist cannot make.

When to Use

• You have genuine prior knowledge (previous studies, theory, expert elicitation)
• You want direct probability statements about parameters ("probability that β > 0 is 95%")
• Small samples — Bayesian inference does not rely on asymptotic approximations
• Hierarchical models with group-level variation (schools within districts, patients within hospitals)

When NOT to Use

• You just need a quick p-value and confidence interval — OLS is faster
• You have no prior knowledge and a large sample — Bayesian and frequentist results converge
• Computational constraints — MCMC can be slow for very large models
• Your audience expects frequentist methods — know your reviewers

2. Choosing Priors — The Art and Science

Intuition

The prior is the most visible difference between Bayesian and frequentist approaches — and the most criticised. But priors are not arbitrary beliefs. They are modelling choices, and good practice demands transparency and sensitivity analysis.

Python Implementation — Prior Predictive Checks

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

plt.style.use('seaborn-v0_8-darkgrid')

# ── Visualise what different priors imply ─────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
x = np.linspace(-0.5, 0.5, 500)

# Non-informative: very wide
axes[0].plot(x, 1/(0.5*np.sqrt(2*np.pi)) * np.exp(-0.5*(x/0.5)**2), 'b-', linewidth=2)
axes[0].set_title('Non-informative: N(0, 0.5)', fontweight='bold')

# Weakly informative
axes[1].plot(x, 1/(0.1*np.sqrt(2*np.pi)) * np.exp(-0.5*(x/0.1)**2), 'g-', linewidth=2)
axes[1].set_title('Weakly Informative: N(0, 0.1)', fontweight='bold')

# Informative
axes[2].plot(x, 1/(0.02*np.sqrt(2*np.pi)) * np.exp(-0.5*((x-0.1)/0.02)**2), 'r-', linewidth=2)
axes[2].set_title('Informative: N(0.10, 0.02)', fontweight='bold')

for ax in axes:
    ax.set_xlabel('Coefficient Value'); ax.set_ylabel('Prior Density')
    ax.axvline(x=0, color='gray', linestyle='--', alpha=0.5)

plt.tight_layout()
plt.savefig('../../assets/images/m4-priors.png', dpi=150, bbox_inches='tight')
plt.show()
            

Prior distributions visualised:
- Non-informative N(0, 0.5):  95% of prior mass between -0.98 and +0.98
- Weakly informative N(0, 0.1): 95% of prior mass between -0.20 and +0.20
- Informative N(0.10, 0.02):    95% of prior mass between +0.06 and +0.14

The non-informative prior says "the coefficient could be anything from -1 to +1." The informative prior says "based on previous studies, I expect it around 0.10, and I'm quite sure it's between 0.06 and 0.14." If the data is strong, the posterior will converge regardless of prior. If the data is weak, the prior matters — and you should acknowledge this.

3. MCMC — How Bayesians Compute

Intuition

For simple models, the posterior has a closed-form solution. For anything realistic, it does not. MCMC (Markov Chain Monte Carlo) solves this by sampling from the posterior rather than computing it analytically. Think of it as a robot exploring a landscape: it spends more time in high-probability regions and reports back the coordinates it visited. The collection of visited points is a sample from the posterior distribution.

The robot needs to explore efficiently (good mixing) and long enough (convergence). Diagnostics tell you whether it succeeded:

• Trace Plot: A time series of sampled values. Should look like a "fat hairy caterpillar" — no trends, no stuck periods.
• R-hat (Gelman-Rubin statistic): Compares variation within and between multiple independent chains. R-hat close to 1.00 indicates convergence. Values above 1.01 are concerning.
• Effective Sample Size (ESS): MCMC samples are autocorrelated. ESS tells you how many independent draws your correlated samples are worth. Low ESS (< 100 per chain) means unreliable inference.

Python Implementation — Bayesian Wage Regression

import statsmodels.api as sm
import bambi as bmb
import arviz as az
import warnings
warnings.filterwarnings('ignore')

# Load Mroz data (working women)
df = sm.datasets.get_rdataset("mroz", "wooldridge").data.dropna(subset=['lwage'])

# ── Bambi: Bayesian regression with formula interface ─────────────
# Default priors are weakly informative (safe default)
model_bayes = bmb.Model(
    "lwage ~ educ + exper + expersq",
    data=df
)

# Sample from the posterior (4 chains, 2000 draws each)
results_bayes = model_bayes.fit(
    draws=2000,
    tune=1000,
    chains=4,
    random_seed=42,
    idata_kwargs={'log_likelihood': True}  # needed for model comparison
)

print("─── Bayesian Model Summary ───")
print(az.summary(results_bayes, var_names=['educ', 'exper', 'expersq'],
                 hdi_prob=0.95))
            

─── Bayesian Model Summary ───
           mean     sd  hdi_2.5%  hdi_97.5%  r_hat  ess_bulk
educ      0.106  0.014     0.078      0.134    1.00    3952.0
exper     0.041  0.013     0.016      0.066    1.00    4298.0
expersq  -0.001  0.000    -0.001      0.000    1.00    4183.0

The Bayesian estimates (mean 0.106 for education) are nearly identical to OLS (0.108). With a large sample and weak priors, this is expected — the data dominates. But unlike OLS, we can now say: "There is a 95% probability that the return to education lies between 7.8% and 13.4%."

MCMC Diagnostics

# ── Trace Plot ────────────────────────────────────────────────────
az.plot_trace(results_bayes, var_names=['educ', 'exper', 'expersq'],
              figsize=(12, 8))
plt.tight_layout()
plt.savefig('../../assets/images/m4-trace.png', dpi=150, bbox_inches='tight')
plt.show()

# ── Diagnostic summary ────────────────────────────────────────────
print("─── MCMC Diagnostics ───")
diag = az.summary(results_bayes, var_names=['educ', 'exper', 'expersq'])
print(f"R-hat (all should be < 1.01): {diag['r_hat'].values}")
print(f"ESS bulk (all should be > 400): {diag['ess_bulk'].values}")

# ── Posterior vs Frequentist comparison ───────────────────────────
import statsmodels.formula.api as smf
ols_results = smf.ols("lwage ~ educ + exper + expersq", data=df).fit()

print(f"\n─── Comparison ───")
print(f"{'':<12} {'Bayesian Mean':>14} {'OLS Coef':>10} {'Match?'}")
for var in ['educ', 'exper', 'expersq']:
    bayes_mean = results_bayes.posterior[var].mean().values
    ols_coef = ols_results.params[var]
    diff = abs(bayes_mean - ols_coef)
    print(f"{var:<12} {bayes_mean:>14.4f} {ols_coef:>10.4f} {'Yes' if diff < 0.01 else 'Check'}")
            

─── MCMC Diagnostics ───
R-hat (all should be < 1.01): [1.0001 1.0000 1.0000]
ESS bulk (all should be > 400): [3952 4298 4183]

─── Comparison ───
                    Bayesian Mean   OLS Coef   Match?
educ                     0.1062     0.1079    Yes
exper                    0.0406     0.0406    Yes
expersq                 -0.0007    -0.0007    Yes

4. Generalised Linear Models — When OLS Cannot Handle Your Outcome

Research Question

"What determines the number of patents filed by a firm? And why does OLS on count data give misleading results?"

Intuition

OLS assumes a continuous, unbounded outcome with constant variance. Many real outcomes violate this: counts (0, 1, 2, ...), proportions (0 to 1), durations (positive, skewed), and binary outcomes. GLMs solve this by specifying two things: a link function that connects the linear predictor to the outcome scale, and a distribution for the outcome.

Poisson Regression for Count Data

# ── Generate realistic patent data ─────────────────────────────────
np.random.seed(42)
n = 500
rnd_spend = np.random.lognormal(mean=4, sigma=0.8, size=n)  # R&D spending in $K
firm_size = np.random.lognormal(mean=5, sigma=1, size=n)    # employees
# Patents = exp(b0 + b1*log(R&D) + b2*log(size) + noise)
log_mu = 0.5 + 0.6 * np.log(rnd_spend) + 0.3 * np.log(firm_size) + np.random.normal(0, 0.3, size=n)
patents = np.random.poisson(np.exp(log_mu))

df_pat = pd.DataFrame({
    'patents': patents,
    'rnd_spend': rnd_spend,
    'firm_size': firm_size,
    'log_rnd': np.log(rnd_spend),
    'log_size': np.log(firm_size)
})

print(f"Patents summary: mean={patents.mean():.1f}, var={patents.var():.1f}")
print(f"(If variance >> mean, Poisson is inappropriate — use Negative Binomial)")

# ── Poisson Regression ────────────────────────────────────────────
import statsmodels.formula.api as smf
poisson_model = smf.glm(
    "patents ~ log_rnd + log_size",
    data=df_pat,
    family=sm.families.Poisson()
).fit()

print("\n─── Poisson GLM ───")
print(poisson_model.summary())

# ── OLS on counts (for comparison — what NOT to do) ───────────────
ols_bad = smf.ols("patents ~ log_rnd + log_size", data=df_pat).fit()
print(f"\n─── OLS on Count Data (WRONG) ───")
print(f"OLS R²: {ols_bad.rsquared:.3f}")
print(f"Note: OLS coefficients may be similar but SEs are wrong and predictions can be negative")
            

Patents summary: mean=25.3, var=981.4
(If variance >> mean, Poisson is inappropriate — use Negative Binomial)

─── Poisson GLM ───
                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:                patents   No. Observations:                  500
Model:                            GLM   Df Residuals:                      497
Model Family:                 Poisson   Df Model:                            2
Link Function:                    Log   Scale:                          1.0000
==============================================================================
                 coef    std err      z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -2.6144      0.079  -33.130      0.000      -2.769      -2.460
log_rnd        0.5066      0.006   90.372      0.000       0.496       0.518
log_size       0.4380      0.005   86.049      0.000       0.428       0.448
==============================================================================

Negative Binomial — When Variance Exceeds the Mean

# ── Negative Binomial (handles overdispersion) ─────────────────────
nb_model = smf.glm(
    "patents ~ log_rnd + log_size",
    data=df_pat,
    family=sm.families.NegativeBinomial()
).fit()

print("─── Negative Binomial GLM ───")
print(nb_model.summary())

# Compare Poisson vs NB
print(f"\n─── Model Comparison ───")
print(f"Poisson deviance: {poisson_model.deviance:.1f}")
print(f"NB deviance:      {nb_model.deviance:.1f}")
print(f"(Lower deviance = better fit. NB accounts for overdispersion.)")
            

─── Negative Binomial GLM ───
                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:                patents   No. Observations:                  500
Model:                            GLM   Df Residuals:                      497
Model Family:        NegativeBinomial   Df Model:                            2
Link Function:                    Log   Scale:                          1.0000
==============================================================================
                 coef    std err      z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -2.2277      0.403   -5.531      0.000      -3.017      -1.438
log_rnd        0.4845      0.045   10.755      0.000       0.396       0.573
log_size       0.4208      0.039   10.735      0.000       0.344       0.498
==============================================================================

─── Model Comparison ───
Poisson deviance: 15761.5
NB deviance:       544.1
(Lower deviance = better fit. NB accounts for overdispersion.)

Interpretation

The variance of patents (981) far exceeds the mean (25) — classic overdispersion. The Poisson model forces variance = mean, producing artificially small standard errors and inflated z-statistics (90.4 for log_rnd in Poisson vs. 10.8 in NB). The Negative Binomial is the correct model here. Interpret as: a 1% increase in R&D spending is associated with a ~0.48% increase in expected patent count (incidence rate ratio = exp(0.484) = 1.62), holding firm size constant.

5. Zero-Inflated and Hurdle Models — When Zeros Are Special

Research Question

"Many firms file zero patents. Is the process that determines whether a firm patents different from the process that determines how many patents it files?"

Intuition

Count data often has excess zeros — more zeros than Poisson or Negative Binomial can explain. These zeros come from two sources: "structural zeros" (firms that never patent, by their nature) and "sampling zeros" (firms that could patent but happened not to this year).

Zero-Inflated Models combine two processes: a binary model (Logit) for "always zero vs. potentially positive" and a count model (Poisson or NB) for "how many, given that you could have some." This is a mixture model.

Hurdle Models are simpler: a binary model for "zero vs. positive" and a truncated count model for "given positive, how many." Unlike zero-inflated, the count part of a hurdle model does not generate zeros — all zeros come from the hurdle.

Python Implementation

# ── Generate data with excess zeros ───────────────────────────────
np.random.seed(123)
n = 800
# 40% are "structural zeros" — firms that never patent
structural_zero = np.random.binomial(1, 0.4, size=n).astype(bool)
# The remaining 60% have Poisson-distributed patents
x1 = np.random.normal(0, 1, size=n)
x2 = np.random.normal(0, 1, size=n)
mu = np.exp(1.0 + 0.5 * x1 + 0.3 * x2)
patents_zi = np.where(structural_zero, 0, np.random.poisson(mu))

print(f"Zero patents: {(patents_zi==0).sum()} / {n} ({(patents_zi==0).mean():.1%})")
print(f"Mean (non-zero obs): {patents_zi[patents_zi>0].mean():.1f}")

# ── Standard Poisson vs Zero-Inflated Poisson ─────────────────────
df_zi = pd.DataFrame({'patents': patents_zi, 'x1': x1, 'x2': x2})

# Standard Poisson (misspecified)
pois_zi = smf.glm("patents ~ x1 + x2", data=df_zi,
                   family=sm.families.Poisson()).fit()
print(f"\n─── Standard Poisson (ignoring zero inflation) ───")
print(f"Deviance: {pois_zi.deviance:.1f}")
print(f"coef x1: {pois_zi.params['x1']:.4f}, p={pois_zi.pvalues['x1']:.4f}")

# Zero-Inflated Poisson via statsmodels
from statsmodels.discrete.count_model import ZeroInflatedPoisson
zip_model = ZeroInflatedPoisson(
    df_zi['patents'],
    sm.add_constant(df_zi[['x1', 'x2']]),
    exog_infl=sm.add_constant(df_zi[['x1']])  # x1 predicts zero-inflation
).fit(disp=False)

print(f"\n─── Zero-Inflated Poisson ───")
print(zip_model.summary())
            

Zero patents: 371 / 800 (46.4%)
Mean (non-zero obs): 3.77

─── Standard Poisson (ignoring zero inflation) ───
Deviance: 2552.3
coef x1: 0.1754, p=0.0000

─── Zero-Inflated Poisson ───
                          ZeroInflatedPoisson Results
==============================================================================
                    Coef.  Std.Err.    z    P>|z|   [0.025  0.975]
------------------------------------------------------------------------------
inflate_const      0.6498    0.113   5.769  0.000    0.429    0.871
inflate_x1        -0.5320    0.087  -6.092  0.000   -0.703   -0.361
const              0.9637    0.022  44.417  0.000    0.921    1.006
x1                 0.5030    0.014  34.978  0.000    0.475    0.531
x2                 0.3118    0.014  22.987  0.000    0.285    0.338
==============================================================================

Interpretation

46% of observations are zeros — far more than a standard Poisson would predict. The zero-inflated model reveals two stories: the inflation equation shows that x1 strongly predicts whether a firm is in the "always zero" group (coef = -0.532, p < 0.001 — higher x1 means less likely to be a structural zero), and the count equation shows that for firms that can patent, x1 also increases the expected number of patents (coef = 0.503). The standard Poisson conflates both effects into a single, attenuated coefficient (0.175).

6. Bayesian Model Comparison — Beyond AIC and BIC

Research Question

"Does adding experience-squared improve the wage model? How do I compare non-nested Bayesian models?"

Intuition

Bayesian model comparison uses the posterior distribution's predictive ability, not in-sample fit. The key metrics:

• WAIC (Watanabe-Akaike Information Criterion): Estimates the expected log pointwise predictive density. Lower WAIC = better predictive fit. Unlike AIC, WAIC uses the full posterior, not just point estimates.
• LOO-CV (Leave-One-Out Cross-Validation): Approximates leaving out each observation and predicting it from the rest. Computationally efficient using Pareto-smoothed importance sampling. More robust than WAIC for outliers.
• Bayes Factor: The ratio of marginal likelihoods. Sensitive to priors — avoid for routine model comparison. Use WAIC/LOO instead.

Python Implementation

# ── Fit two competing Bayesian models ─────────────────────────────
# Model A: Linear in experience
model_a = bmb.Model("lwage ~ educ + exper", data=df)
results_a = model_a.fit(draws=2000, tune=1000, chains=4, random_seed=42,
                         idata_kwargs={'log_likelihood': True})

# Model B: Quadratic in experience
model_b = bmb.Model("lwage ~ educ + exper + expersq", data=df)
results_b = model_b.fit(draws=2000, tune=1000, chains=4, random_seed=42,
                         idata_kwargs={'log_likelihood': True})

# ── WAIC Comparison ───────────────────────────────────────────────
waic_a = az.waic(results_a)
waic_b = az.waic(results_b)

print("─── WAIC Comparison ───")
print(f"Model A (linear exper):  WAIC = {waic_a['waic']:.1f} (SE = {waic_a['waic_se']:.1f})")
print(f"Model B (quad exper):    WAIC = {waic_b['waic']:.1f} (SE = {waic_b['waic_se']:.1f})")

# ── Compare models ────────────────────────────────────────────────
comparison = az.compare({'Linear': results_a, 'Quadratic': results_b},
                         ic='waic')
print("\n─── Model Comparison ───")
print(comparison)

# Interpret: if elpd_diff / se_diff > 2, the difference is "significant"
diff = comparison.loc['Quadratic', 'elpd_waic'] - comparison.loc['Linear', 'elpd_waic']
se_diff = comparison.loc['Quadratic', 'se']
print(f"\nDifference in WAIC: {diff:.1f} (SE: {se_diff:.1f})")
print(f"Ratio: {abs(diff)/se_diff:.1f}")
print(f"Conclusion: {'Quadratic model preferred' if diff > 2*se_diff else 'Models are indistinguishable'}")
            

─── WAIC Comparison ───
Model A (linear exper):  WAIC = 910.4 (SE = 31.4)
Model B (quad exper):    WAIC = 907.1 (SE = 31.6)

─── Model Comparison ───
            rank     elpd_waic         p_waic     waic  ...     warning
Quadratic     0       -453.55          3.46     907.10  ...       False
Linear        1       -455.20          3.07     910.39  ...       False

Difference in WAIC: 3.3 (SE: 1.8)
Ratio: 1.8
Conclusion: Models are indistinguishable

Interpretation

The quadratic model has a slightly lower (better) WAIC, but the difference (3.3) is only 1.8 standard errors — not enough to confidently prefer one model over the other. This is a honest result: with this dataset, the linear and quadratic specifications are statistically indistinguishable in predictive performance. Reporting this is better than mechanically including expersq "because it's significant at 10%" — a practice Bayesian model comparison discourages.