Logistic regression priors – Applied Bayesian Computing in Python

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,
    )

from scipy.special import expit
rng = np.random.default_rng(18501)
data = pd.DataFrame({"x": rng.normal(size=120)})
data["p"] = expit(-0.2 + 1.8 * data.x)
data["y"] = rng.binomial(1, data.p)
data.head()

	x	p	y
0	-0.823600	0.156767	0
1	-0.303863	0.321485	1
2	0.172216	0.527470	0
3	-1.658365	0.039732	0
4	0.202274	0.540931	0

weak = fit_glm("y ~ x", data, family=lp.bernoulli(), seed=18502, prior_scale=5, intercept_scale=5)
regularized = fit_glm("y ~ x", data, family=lp.bernoulli(), seed=18503, prior_scale=1, intercept_scale=2)
pd.concat({"weak": coef_medians(weak), "regularized": coef_medians(regularized)}, axis=1).round(3)

	weak	regularized
Intercept	-0.129	-0.121
x	1.187	1.116

grid = pd.DataFrame({"x": np.linspace(-3, 3, 160)})
fig, ax = plt.subplots()
for label, fit in {"weak": weak, "regularized": regularized}.items():
    ax.plot(grid.x, fit.posterior_epred(grid).mean(axis=0), label=label)
ax.scatter(data.x, data.y, alpha=.25, color="0.4")
ax.legend(frameon=False)

Simulated data