Robit and logit links

Source: Robit/robit.Rmd

lapylace currently supplies the logit GLM path. For comparison with heavier-tailed robit intuition, this page shows the logit fit and a Student-t CDF curve with matched slope.

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
from scipy.stats import t
rng = np.random.default_rng(19101)
data = pd.DataFrame({"x": rng.normal(size=180)})
data["p"] = expit(-0.4 + 1.6*data.x)
data["y"] = rng.binomial(1, data.p)
fit = fit_glm("y ~ x", data, family=lp.bernoulli(), seed=19102, prior_scale=2.5, intercept_scale=5)
coef_medians(fit).round(3)

Intercept   -0.431
x            1.701
dtype: float64

grid = pd.DataFrame({"x": np.linspace(data.x.min(), data.x.max(), 200)})
eta = fit.linear_predictor(grid).mean(axis=0)
fig, ax = plt.subplots()
ax.scatter(data.x, data.y, alpha=.25)
ax.plot(grid.x, 1/(1+np.exp(-eta)), label="logit")
ax.plot(grid.x, t.cdf(eta, df=4), label="t CDF, df=4")
ax.legend(frameon=False)