Source: Roaches/roaches.Rmd
The core count model can be represented as a Poisson log-link GLM. lapylace supplies the formula interface and posterior predictive checks.
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,
)
roaches = pd.read_csv(root / "Roaches/data/roaches.csv", index_col=0)
roaches["log_roach1"] = np.log1p(roaches["roach1"])
roaches["log_exposure2"] = np.log(roaches["exposure2"])
roaches.head()
| 1 |
153 |
308.00 |
1 |
0 |
0.800000 |
5.733341 |
-0.223144 |
| 2 |
127 |
331.25 |
1 |
0 |
0.600000 |
5.805888 |
-0.510826 |
| 3 |
7 |
1.67 |
1 |
0 |
1.000000 |
0.982078 |
0.000000 |
| 4 |
7 |
3.00 |
1 |
0 |
1.000000 |
1.386294 |
0.000000 |
| 5 |
0 |
2.00 |
1 |
0 |
1.142857 |
1.098612 |
0.133531 |
fit = fit_glm("y ~ log_roach1 + treatment + senior + log_exposure2", roaches, family=lp.poisson(), seed=18201, prior_scale=2.5, intercept_scale=5)
fit.summary(["alpha", "beta"])
| alpha |
1.702370 |
0.002044 |
0.039898 |
1.637050 |
1.701730 |
1.769040 |
380.945 |
2842.87 |
1.004810 |
| beta[1] |
0.587532 |
0.000412 |
0.008703 |
0.573138 |
0.587481 |
0.601825 |
447.154 |
3336.97 |
1.003180 |
| beta[2] |
-0.586398 |
0.000793 |
0.024583 |
-0.627744 |
-0.586724 |
-0.547011 |
960.436 |
7167.43 |
0.999164 |
| beta[3] |
-0.241843 |
0.001135 |
0.032736 |
-0.294662 |
-0.240005 |
-0.187897 |
832.077 |
6209.53 |
0.999140 |
| beta[4] |
0.324706 |
0.002008 |
0.051172 |
0.239293 |
0.323782 |
0.407600 |
649.244 |
4845.10 |
1.002240 |
mu = fit.posterior_epred(roaches).mean(axis=0)
pd.DataFrame({"observed_mean": [roaches.y.mean()], "predicted_mean": [mu.mean()], "observed_zero_rate": [(roaches.y == 0).mean()]})
| 0 |
25.648855 |
25.643976 |
0.358779 |