Skip to content

Commit 45aace2

Browse files
authored
update GLM binomial (#868)
1 parent c121d94 commit 45aace2

2 files changed

Lines changed: 191 additions & 135 deletions

File tree

examples/generalized_linear_models/GLM-binomial-regression.ipynb

Lines changed: 170 additions & 85 deletions
Large diffs are not rendered by default.

examples/generalized_linear_models/GLM-binomial-regression.myst.md

Lines changed: 21 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -73,7 +73,7 @@ from scipy.special import expit
7373

7474
```{code-cell} ipython3
7575
%config InlineBackend.figure_format = 'retina'
76-
az.style.use("arviz-darkgrid")
76+
az.style.use("arviz-variat")
7777
rng = np.random.default_rng(1234)
7878
```
7979

@@ -119,8 +119,8 @@ freq.set_ylabel("number of successes")
119119
freq.scatter(x, y, color="k")
120120
# plot proportion related stuff on ax[1]
121121
ax[1].plot(x, p_true, label=r"$g^{-1}(β_0 + β_1 \cdot x_i)$")
122-
ax[1].set_ylabel("proportion successes", color="b")
123-
ax[1].tick_params(axis="y", labelcolor="b")
122+
ax[1].set_ylabel("proportion successes", color="C0")
123+
ax[1].tick_params(axis="y", labelcolor="C0")
124124
ax[1].set(xlabel="$x$", title="Binomial regression")
125125
ax[1].legend()
126126
# get y-axes to line up
@@ -142,7 +142,7 @@ Technically, we don't need to supply `coords`, but providing this (a list of obs
142142
coords = {"observation": data.index.values}
143143
144144
with pm.Model(coords=coords) as binomial_regression_model:
145-
x = pm.ConstantData("x", data["x"], dims="observation")
145+
x = pm.Data("x", data["x"], dims="observation")
146146
# priors
147147
beta0 = pm.Normal("beta0", mu=0, sigma=1)
148148
beta1 = pm.Normal("beta1", mu=0, sigma=1)
@@ -165,7 +165,7 @@ with binomial_regression_model:
165165
Confirm no inference issues by visual inspection of chain. We've got no warnings about divergences, $\hat{R}$, or effective sample size. Everything looks good.
166166

167167
```{code-cell} ipython3
168-
az.plot_trace(idata, var_names=["beta0", "beta1"]);
168+
az.plot_rank_dist(idata, var_names=["beta0", "beta1"]);
169169
```
170170

171171
## Examine results
@@ -174,55 +174,25 @@ The code below plots out model predictions in data space, and our posterior beli
174174
```{code-cell} ipython3
175175
:tags: [hide-input]
176176
177-
fig, ax = plt.subplots(1, 2, figsize=(9, 4), gridspec_kw={"width_ratios": [2, 1]})
178-
179-
# Data space plot ========================================================
180-
az.plot_hdi(
181-
data["x"],
182-
idata.posterior.p,
183-
hdi_prob=0.95,
184-
fill_kwargs={"alpha": 0.25, "linewidth": 0},
185-
ax=ax[0],
186-
color="C1",
177+
pc = az.plot_lm(
178+
idata,
179+
y="p",
180+
y_obs="y",
181+
group="posterior",
182+
ci_prob=0.95,
183+
visuals={"observed_scatter": False},
187184
)
188-
# posterior mean
189-
post_mean = idata.posterior.p.mean(("chain", "draw"))
190-
ax[0].plot(data["x"], post_mean, label="posterior mean", color="C1")
191-
# plot truth
192-
ax[0].plot(data["x"], p_true, "--", label="true", color="C2")
193-
# formatting
194-
ax[0].set(xlabel="x", title="Data space")
195-
ax[0].set_ylabel("proportion successes", color="C1")
196-
ax[0].tick_params(axis="y", labelcolor="C1")
197-
ax[0].legend()
198-
# instantiate a second axes that shares the same x-axis
199-
freq = ax[0].twinx()
200-
freq.set_ylabel("number of successes")
201-
freq.scatter(data["x"], data["y"], color="k", label="data")
202-
# get y-axes to line up
203-
y_buffer = 1
204-
freq.set(ylim=[-y_buffer, n + y_buffer])
205-
ax[0].set(ylim=[-(y_buffer / n), 1 + (y_buffer / n)])
206-
freq.grid(None)
207-
# set both y-axis to have 5 ticks
208-
ax[0].set(yticks=np.linspace(0, 20, 5) / n)
209-
freq.set(yticks=np.linspace(0, 20, 5))
210-
211-
# Parameter space plot ===================================================
212-
az.plot_kde(
213-
az.extract(idata, var_names="beta0"),
214-
az.extract(idata, var_names="beta1"),
215-
contourf_kwargs={"cmap": "Blues"},
216-
ax=ax[1],
217-
)
218-
ax[1].plot(beta0_true, beta1_true, "C2o", label="true")
219-
ax[1].set(xlabel=r"$\beta_0$", ylabel=r"$\beta_1$", title="Parameter space")
220-
ax[1].legend(facecolor="white", frameon=True);
185+
pc.map(az.visuals.line_xy, "true", x=data["x"], y=p_true, color="C1", label="true")
186+
pc.map(az.visuals.scatter_xy, "observed", x=data["x"], y=data["y"] / n, color="k")
187+
plt.legend()
188+
189+
pc = az.plot_dist(idata, var_names=["beta0", "beta1"])
190+
az.add_lines(pc, {"beta0": beta0_true, "beta1": beta1_true});
221191
```
222192

223-
The left panel shows the posterior mean (solid line) and 95% credible intervals (shaded region). Because we are working with simulated data, we know what the true model is, so we can see that the posterior mean compares favourably with the true data generating model.
193+
The top panel shows the posterior mean (solid line) and 95% credible intervals (shaded region). Because we are working with simulated data, we know what the true model is, so we can see that the posterior mean compares favourably with the true data generating model.
224194

225-
This is also shown by the posterior distribution over parameter space (right panel), which does well when comparing to the true data generating parameters.
195+
This is also shown by the posterior distribution over parameter space (bottom panel), which does well when comparing to the true data generating parameters.
226196

227197
Using binomial regression in real data analysis situations would probably involve more predictor variables, and correspondingly more model parameters, but hopefully this example has demonstrated the logic behind binomial regression.
228198

@@ -235,6 +205,7 @@ A good introduction to generalized linear models is provided by {cite:t}`roback2
235205
- Updated by [Benjamin T. Vincent](https://github.com/drbenvincent) in February 2022
236206
- Updated by Benjamin T. Vincent in February 2023 to run on PyMC v5
237207
- Updated to use `az.extract` by [Benjamin T. Vincent](https://github.com/drbenvincent) in February 2023, ([pymc-examples#522](https://github.com/pymc-devs/pymc-examples/pull/522))
208+
- Updated by Osvaldo Martin in April 2026
238209

239210
+++
240211

0 commit comments

Comments
 (0)