Logistic regression graphs

Source: LogitGraphs/logitgraphs.Rmd

The R page uses logistic regression fits for plotting. Here the fits are Bayesian Bernoulli GLMs from lapylace.

Simulated data

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(18401)
fake = pd.DataFrame({"x": rng.normal(size=250)})
fake["p"] = expit(-0.3 + 1.4 * fake.x)
fake["y"] = rng.binomial(1, fake.p)
fake.head()

	x	p	y
0	0.535821	0.610675	1
1	0.236657	0.507829	1
2	-0.679870	0.222386	0
3	-1.496108	0.083589	0
4	1.403239	0.840846	1

fit = fit_glm("y ~ x", fake, family=lp.bernoulli(), seed=18402, prior_scale=2.5, intercept_scale=5)
coef_medians(fit).round(3)

Intercept   -0.194
x            1.437
dtype: float64

grid = pd.DataFrame({"x": np.linspace(fake.x.min(), fake.x.max(), 200)})
prob = fit.posterior_epred(grid)
fig, ax = plt.subplots()
ax.scatter(fake.x, fake.y, alpha=.25)
ax.plot(grid.x, prob.mean(axis=0), color="black")
ax.fill_between(grid.x, np.quantile(prob,.1,axis=0), np.quantile(prob,.9,axis=0), color="0.85")
ax.set_xlabel("x")
ax.set_ylabel("Pr(y=1)")

Text(0, 0.5, 'Pr(y=1)')