# NES: logistic regression
Source: `NES/nes_logistic.Rmd`
The original page fits `stan_glm(..., family=binomial(link="logit"))`. This translation uses the same Bayesian logistic-regression path with `lapylace`.
## Setup
```{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,
)
```
```{python}
nes = pd.read_csv(root / "NES/data/nes.txt", sep=r"\s+", index_col=0, na_values=["NA"])
nes92 = nes.loc[nes.year == 1992, ["rvote", "income", "female", "black"]].dropna().copy()
nes92["rvote"] = nes92["rvote"].astype(int)
nes92.head()
```
```{python}
fit_1 = fit_glm("rvote ~ income", nes92, family=lp.bernoulli(), seed=18001, prior_scale=2.5, intercept_scale=5)
fit_2 = fit_glm("rvote ~ female + black + income", nes92, family=lp.bernoulli(), seed=18002, prior_scale=2.5, intercept_scale=5)
pd.concat({"income": coef_medians(fit_1), "demographics": coef_medians(fit_2)}, axis=1).round(3)
```
```{python}
grid = pd.DataFrame({"income": np.arange(1, 6), "female": 0, "black": 0})
prob = fit_2.posterior_epred(grid)
pd.DataFrame({"income": grid.income, "Pr(Republican vote)": prob.mean(axis=0), "lo80": np.quantile(prob,.1,axis=0), "hi80": np.quantile(prob,.9,axis=0)})
```