Residuals and fitted values

Source: Residuals/residuals.Rmd

The original examples use fitted values and residuals from stan_glm; here those come from posterior expected values under 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,
    )

rng = np.random.default_rng(18701)
fake = pd.DataFrame({"x": rng.normal(size=150), "z": rng.binomial(1, .45, size=150)})
fake["y"] = 1 + 2*fake.x + 1.5*fake.z + rng.normal(0, 1, size=len(fake))
fit = fit_glm("y ~ x + z", fake, seed=18702, prior_scale=5, intercept_scale=5, aux_scale=2)
fake["fitted"] = fit.posterior_epred(fake).mean(axis=0)
fake["resid"] = fake["y"] - fake["fitted"]
fake.head()

	x	z	y	fitted	resid
0	-2.264170	0	-2.156759	-3.810696	1.653937
1	1.134729	0	5.796919	3.185905	2.611014
2	0.183398	1	2.225688	2.826868	-0.601180
3	0.473300	0	1.628876	1.824360	-0.195484
4	-0.435070	0	-0.243358	-0.045511	-0.197846

fig, axes = plt.subplots(1, 2, figsize=(8,3))
axes[0].scatter(fake.fitted, fake.y, alpha=.5)
axes[0].set_xlabel("Fitted")
axes[0].set_ylabel("Observed")
axes[1].scatter(fake.fitted, fake.resid, alpha=.5)
axes[1].axhline(0, color="black", lw=1)
axes[1].set_xlabel("Fitted")
axes[1].set_ylabel("Residual")

Text(0, 0.5, 'Residual')