Hibbs regression: interval coverage

Source: ElectionsEconomy/hibbs_coverage.Rmd

The coverage simulation repeatedly fits the same Gaussian model with lapylace and checks whether posterior intervals cover the generating slope.

Setup

from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lapylace as lp

root = Path("../../ROS-Examples")

def coef_medians(fit):
    return pd.Series(
        [np.median(fit.alpha_draws()), *np.median(fit.beta_draws(), axis=0)],
        index=["Intercept", *fit.columns],
    )

def fit_glm(formula, data, family=None, seed=1, prior_scale=2.5, intercept_scale=5, aux_scale=10, **kwargs):
    return lp.stan_glm(
        formula,
        data=data,
        family=family or lp.gaussian(),
        prior=lp.normal(0, prior_scale),
        prior_intercept=lp.normal(0, intercept_scale),
        prior_aux=lp.exponential(aux_scale),
        chains=2,
        parallel_chains=2,
        iter_warmup=300,
        iter_sampling=500,
        seed=seed,
        refresh=100,
        **kwargs,
    )

hibbs = pd.read_csv(root / "ElectionsEconomy/data/hibbs.dat", sep=r"\s+")
base_fit = fit_glm("vote ~ growth", hibbs, seed=17801, prior_scale=10, intercept_scale=100, aux_scale=10)
coef_medians(base_fit).round(2)

Intercept    46.19
growth        3.09
dtype: float64

rng = np.random.default_rng(17802)
x = hibbs.growth.to_numpy()
true_alpha, true_beta, sigma = 46.0, 3.0, 3.5
rows = []
for s in range(8):
    fake = pd.DataFrame({"x": x, "y": true_alpha + true_beta*x + rng.normal(0, sigma, len(x))})
    fit = fit_glm("y ~ x", fake, seed=17810+s, prior_scale=10, intercept_scale=100, aux_scale=10)
    beta = fit.beta_draws()[:, 0]
    lo, hi = np.quantile(beta, [.1, .9])
    rows.append({"sim": s+1, "median": np.median(beta), "lo80": lo, "hi80": hi, "covers": lo <= true_beta <= hi})
pd.DataFrame(rows)

	sim	median	lo80	hi80	covers
0	1	3.804800	3.008717	4.617024	False
1	2	3.839445	2.797728	4.673550	True
2	3	2.482030	1.887453	3.176729	True
3	4	2.821425	1.820033	3.964031	True
4	5	2.786345	1.549725	3.924712	True
5	6	3.276315	2.682132	3.885255	True
6	7	3.147620	2.264992	4.052663	True
7	8	2.000405	0.949227	2.936976	False