Hierarchical Logistic regression with Binomial family#

This notebook shows how to build a hierarchical logistic regression model with the Binomial family in Bambi.

This example is based on the Hierarchical baseball article in Bayesian Analysis Recipes, a collection of articles on how to do Bayesian data analysis with PyMC3 made by Eric Ma.

Problem description#

Extracted from the original work:

Baseball players have many metrics measured for them. Let’s say we are on a baseball team, and would like to quantify player performance, one metric being their batting average (defined by how many times a batter hit a pitched ball, divided by the number of times they were up for batting (“at bat”)). How would you go about this task?

Load libraries and data#

[1]:
import arviz as az
import bambi as bmb
import matplotlib.pyplot as plt
import numpy as np

from matplotlib.lines import Line2D
from matplotlib.patches import Patch
[2]:
az.style.use("arviz-darkgrid")
random_seed = 1234

We first need some measurements of batting data. Today we’re going to use data from the Baseball Databank. It is a compilation of historical baseball data in a convenient, tidy format, distributed under Open Data terms.

This repository contains several datasets in the form of .csv files. This example is going to use the Batting.csv file, which can be loaded directly with Bambi in a convenient way.

[3]:
df = bmb.load_data("batting")

# Then clean some of the data
df["AB"] = df["AB"].replace(0, np.nan)
df = df.dropna()
df["batting_avg"] = df["H"] / df["AB"]
df = df[df["yearID"] >= 2016]
df = df.iloc[0:15]
df.head(5)
[3]:
playerID yearID stint teamID lgID G AB R H 2B ... SB CS BB SO IBB HBP SH SF GIDP batting_avg
101348 abadfe01 2016 1 MIN AL 39 1.0 0 0 0 ... 0.0 0.0 0 1.0 0.0 0.0 0.0 0.0 0.0 0.000000
101350 abreujo02 2016 1 CHA AL 159 624.0 67 183 32 ... 0.0 2.0 47 125.0 7.0 15.0 0.0 9.0 21.0 0.293269
101352 ackledu01 2016 1 NYA AL 28 61.0 6 9 0 ... 0.0 0.0 8 9.0 0.0 0.0 0.0 1.0 0.0 0.147541
101353 adamecr01 2016 1 COL NL 121 225.0 25 49 7 ... 2.0 3.0 24 47.0 0.0 4.0 3.0 0.0 5.0 0.217778
101355 adamsma01 2016 1 SLN NL 118 297.0 37 74 18 ... 0.0 1.0 25 81.0 1.0 2.0 0.0 3.0 5.0 0.249158

5 rows × 23 columns

From all the columns above, we’re going to use the following:

  • playerID: Unique identification for the player.

  • AB: Number of times the player was up for batting.

  • H: Number of times the player hit the ball while batting.

  • batting_avg: Simply ratio between H and AB.

Explore our data#

It’s always good to explore the data before starting to write down our models. This is very useful to gain a good understanding of the distribution of the variables and their relationships, and even anticipate some problems that may occur during the sampling process.

The following graph summarizes the percentage of hits, as well as the number of times the players were up for batting and the number of times they hit the ball.

[4]:
BLUE = "#2a5674"
RED = "#b13f64"
[5]:
_, ax = plt.subplots(figsize=(10, 6))

# Customize x limits.
# This adds space on the left side to indicate percentage of hits.
ax.set_xlim(-120, 320)

# Add dots for the times at bat and the number of hits
ax.scatter(df["AB"], list(range(15)), s=140, color=BLUE, zorder=10)
ax.scatter(df["H"], list(range(15)), s=140, color=RED, zorder=10)

# Also a line connecting them
ax.hlines(list(range(15)), df["AB"], df["H"], color="#b3b3b3", lw=4)

ax.axvline(ls="--", lw=1.4, color="#a3a3a3")
ax.hlines(list(range(15)), -110, -50, lw=6, color="#b3b3b3", capstyle="round")
ax.scatter(60 * df["batting_avg"] - 110, list(range(15)), s=28, color=RED, zorder=10)

# Add the percentage of hits
for j in range(15):
    text = f"{round(df['batting_avg'].iloc[j] * 100)}%"
    ax.text(-12, j, text, ha="right", va="center", fontsize=14, color="#333")

# Customize tick positions and labels
ax.yaxis.set_ticks(list(range(15)))
ax.yaxis.set_ticklabels(df["playerID"])
ax.xaxis.set_ticks(range(0, 400, 100))

# Create handles for the legend (just dots and labels)
handles = [
    Line2D(
        [0], [0], label="At Bat", marker="o", color="None", markeredgewidth=0,
        markerfacecolor=RED, markersize=12
    ),
    Line2D(
        [0], [0], label="Hits", marker="o", color="None", markeredgewidth=0,
        markerfacecolor=BLUE, markersize=13
    )
]

# Add legend on top-right corner
legend = ax.legend(
    handles=handles,
    loc=1,
    fontsize=14,
    handletextpad=0.4,
    frameon=True
)

# Finally add labels and a title
ax.set_xlabel("Count", fontsize=14)
ax.set_ylabel("Player", fontsize=14)
ax.set_title("How often do batters hit the ball?", fontsize=20);
../_images/notebooks_hierarchical_binomial_bambi_10_0.png

The first thing one can see is that the number of times players were up for batting varies quite a lot. Some players have been there for very few times, while there are others who have been there hundreds of times. We can also note the percentage of hits is usually a number between 12% and 29%.

There are two players, alberma01 and abadfe01, who had only one chance to bat. The first one hit the ball, while the latter missed. That’s why alberma01 as a 100% hit percentage, while abadfe01 has 0%. There’s another player, aguilje01, who has a success record of 0% because he missed all the few opportunities he had to bat. These extreme situations, where the empirical estimation lives in the boundary of the parameter space, are associated with estimation problems when using a maximum-likelihood estimation approach. Nonetheless, they can also impact the sampling process, especially when using wide priors.

As a final note, abreujo02, has been there for batting 624 times, and thus the grey dot representing this number does not appear in the plot.

Non-hierarchical model#

Let’s get started with a simple cell-means logistic regression for \(p_i\), the probability of hitting the ball for the player \(i\)

\[\begin{array}{lr} \displaystyle \text{logit}(p_i) = \beta_i & \text{with } i = 0, \cdots, 14 \end{array}\]

Where

\[\beta_i \sim \text{Normal}(0, \ \sigma_{\beta}),\]

\(\sigma_{\beta}\) is a common constant for all the players, and \(\text{logit}(p_i) = \log\left(\frac{p_i}{1 - p_i}\right)\).

Specifying this model is quite simple in Bambi thanks to its formula interface.

First of all, note this is a Binomial family and the response involves both the number of hits (H) and the number of times at bat (AB). We use the p(x, n) function for the response term. This just tells Bambi we want to model the proportion resulting from dividing x over n.

The right-hand side of the formula is "0 + playerID". This means the model includes a coefficient for each player ID, but does not include a global intercept.

Finally, using the Binomial family is as easy as passing family="binomial". By default, the link function for this family is link="logit", so there’s nothing to change there.

[6]:
model_non_hierarchical = bmb.Model("p(H, AB) ~ 0 + playerID", df, family="binomial")
model_non_hierarchical
[6]:
Formula: p(H, AB) ~ 0 + playerID
Family name: Binomial
Link: logit
Observations: 15
Priors:
  Common-level effects
    playerID ~ Normal(mu: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.], sigma: [10.0223 10.0223 10.0223 10.0223 10.0223 10.0223 10.0223 10.0223 10.0223
 10.0223 10.0223 10.0223 10.0223 10.0223 10.0223])
[7]:
idata_non_hierarchical = model_non_hierarchical.fit(random_seed=random_seed)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [playerID]
100.00% [4000/4000 00:03<00:00 Sampling 2 chains, 0 divergences]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 4 seconds.

Next we observe the posterior of the coefficient for each player. The compact=False argument means we want separated panels for each player.

[8]:
az.plot_trace(idata_non_hierarchical, compact=False);
../_images/notebooks_hierarchical_binomial_bambi_17_0.png

So far so good! The traceplots indicate the sampler worked well.

Now, let’s keep this posterior aside for later use and let’s fit the hierarchical version.

Hierarchical model#

This model incorporates a group-specific intercept for each player:

\[\begin{array}{lr} \displaystyle \text{logit}(p_i) = \alpha + \gamma_i & \text{with } i = 0, \cdots, 14 \end{array}\]

where

\[\begin{split}\begin{array}{c} \alpha \sim \text{Normal}(0, \ \sigma_{\alpha}) \\ \gamma_i \sim \text{Normal}(0, \ \sigma_{\gamma}) \\ \sigma_{\gamma} \sim \text{HalfNormal}(\tau_{\gamma}) \end{array}\end{split}\]

The group-specific terms are indicated with the | operator in the formula. In this case, since there is an intercept for each player, we write 1|playerID.

[9]:
model_hierarchical = bmb.Model("p(H, AB) ~ 1 + (1|playerID)", df, family="binomial")
model_hierarchical
[9]:
Formula: p(H, AB) ~ 1 + (1|playerID)
Family name: Binomial
Link: logit
Observations: 15
Priors:
  Common-level effects
    Intercept ~ Normal(mu: 0, sigma: 2.5)

  Group-level effects
    1|playerID ~ Normal(mu: 0, sigma: HalfNormal(sigma: 2.5))
[10]:
idata_hierarchical = model_hierarchical.fit(random_seed=random_seed)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [1|playerID_offset, 1|playerID_sigma, Intercept]
100.00% [4000/4000 00:04<00:00 Sampling 2 chains, 23 divergences]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 5 seconds.
There were 21 divergences after tuning. Increase `target_accept` or reparameterize.
There were 2 divergences after tuning. Increase `target_accept` or reparameterize.
The number of effective samples is smaller than 25% for some parameters.

And there we got several divergences… What can we do?

One thing we could try is to increase target_accept as suggested in the message above, but there are so many divergences that instead we are going to first take a look at the prior predictive distribution to check whether our priors are too informative or too wide.

The Model instance has a method called prior_predictive() that generates samples from the prior predictive distribution. It returns an InferenceData object that contains the values of the prior predictive distribution.

[11]:
idata_prior = model_hierarchical.prior_predictive()
prior = idata_prior["prior_predictive"]["p(H, AB)"].values.squeeze()

If we inspect the array, we see there are 500 draws (rows) for each of the 15 players (columns)

[12]:
print(prior.shape)
print(prior)
(500, 15)
[[  0 262  39 ...   7 162  98]
 [  1 599  46 ...  11 127  44]
 [  0 205  60 ...   6 182 112]
 ...
 [  1 619  56 ...  15 182 109]
 [  0 124   0 ...   8   1  75]
 [  1 617  35 ...   6 172   0]]

Let’s plot these distributions together with the observed proportion of hits for every player here.

[13]:
# We define this function because this plot is going to be repeated below.
def plot_prior_predictive(df, prior):
    AB = df["AB"].values
    H = df["H"].values

    fig, axes = plt.subplots(5, 3, figsize=(10, 6), sharex="col")

    for idx, ax in enumerate(axes.ravel()):
        pps = prior[:, idx]
        ab = AB[idx]
        h = H[idx]
        hist = ax.hist(pps / ab, bins=25, color="#a3a3a3")
        ax.axvline(h / ab, color=RED, lw=2)
        ax.set_yticks([])
        ax.tick_params(labelsize=12)

    fig.subplots_adjust(left=0.025, right=0.975, hspace=0.05, wspace=0.05, bottom=0.125)
    fig.legend(
        handles=[Line2D([0], [0], label="Observed proportion", color=RED, linewidth=2)],
        handlelength=1.5,
        handletextpad=0.8,
        borderaxespad=0,
        frameon=True,
        fontsize=11,
        bbox_to_anchor=(0.975, 0.92),
        loc="right"

    )
    fig.text(0.5, 0.025, "Prior probability of hitting", fontsize=15, ha="center", va="baseline")
[14]:
plot_prior_predictive(df, prior)
/tmp/ipykernel_10811/4093018835.py:17: UserWarning: This figure was using constrained_layout, but that is incompatible with subplots_adjust and/or tight_layout; disabling constrained_layout.
  fig.subplots_adjust(left=0.025, right=0.975, hspace=0.05, wspace=0.05, bottom=0.125)
../_images/notebooks_hierarchical_binomial_bambi_30_1.png

Indeed, priors are too wide! Let’s use tighter priors and see what’s the result

[15]:
priors = {
    "Intercept": bmb.Prior("Normal", mu=0, sigma=1),
    "1|playerID": bmb.Prior("Normal", mu=0, sigma=bmb.Prior("HalfNormal", sigma=1))
}
model_hierarchical = bmb.Model("p(H, AB) ~ 1 + (1|playerID)", df, family="binomial", priors=priors)
model_hierarchical
[15]:
Formula: p(H, AB) ~ 1 + (1|playerID)
Family name: Binomial
Link: logit
Observations: 15
Priors:
  Common-level effects
    Intercept ~ Normal(mu: 0, sigma: 1)

  Group-level effects
    1|playerID ~ Normal(mu: 0, sigma: HalfNormal(sigma: 1))

Now let’s check the prior predictive distribution for these new priors.

[16]:
model_hierarchical.build()
idata_prior = model_hierarchical.prior_predictive()
prior = idata_prior["prior_predictive"]["p(H, AB)"].values.squeeze()
plot_prior_predictive(df, prior)
/tmp/ipykernel_10811/4093018835.py:17: UserWarning: This figure was using constrained_layout, but that is incompatible with subplots_adjust and/or tight_layout; disabling constrained_layout.
  fig.subplots_adjust(left=0.025, right=0.975, hspace=0.05, wspace=0.05, bottom=0.125)
../_images/notebooks_hierarchical_binomial_bambi_34_1.png

Definetely it looks much better. Now the priors tend to have a symmetric shape with a mode at 0.5, with substantial probability on the whole domain.

[17]:
idata_hierarchical = model_hierarchical.fit(random_seed=random_seed)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [1|playerID_offset, 1|playerID_sigma, Intercept]
100.00% [4000/4000 00:03<00:00 Sampling 2 chains, 38 divergences]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 4 seconds.
There were 3 divergences after tuning. Increase `target_accept` or reparameterize.
There were 35 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.7033138328272285, but should be close to 0.8. Try to increase the number of tuning steps.
The estimated number of effective samples is smaller than 200 for some parameters.

Let’s try with increasing target_accept and the number of tune samples.

[18]:
idata_hierarchical = model_hierarchical.fit(tune=2000, draws=2000, target_accept=0.95, random_seed=random_seed)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [1|playerID_offset, 1|playerID_sigma, Intercept]
100.00% [8000/8000 00:09<00:00 Sampling 2 chains, 0 divergences]
Sampling 2 chains for 2_000 tune and 2_000 draw iterations (4_000 + 4_000 draws total) took 10 seconds.
[19]:
var_names = ["Intercept", "1|playerID", "1|playerID_sigma"]
az.plot_trace(idata_hierarchical, var_names=var_names, compact=False);
../_images/notebooks_hierarchical_binomial_bambi_39_0.png

Let’s jump onto the next section where we plot and compare the probability of hit for the players using both models.

Compare predictions#

Now we’re going to plot the distribution of the probability of hit for each player, using both models.

But before doing that, we need to obtain the posterior in that scale. We could manually take the posterior of the coefficients, compute the linear predictor, and transform that to the probability scale. But that’s a lot of work!

Fortunately, Bambi models have a method called .predict() that we can use to predict in the probability scale. By default, it modifies in-place the InferenceData object we pass to it. Then, the posterior samples can be found in the variable p(H, AB)_mean.

[20]:
model_non_hierarchical.predict(idata_non_hierarchical)
model_hierarchical.predict(idata_hierarchical)

Let’s create a forestplot using the posteriors obtained with both models so we can compare them very easily .

[23]:
_, ax = plt.subplots(figsize = (8, 8))

# Add vertical line for the global probability of hitting
ax.axvline(x=(df["H"] / df["AB"]).mean(), ls="--", color="black", alpha=0.5)

# Create forestplot with ArviZ, only for the mean.
az.plot_forest(
    [idata_non_hierarchical, idata_hierarchical],
    var_names="p(H, AB)_mean",
    combined=True,
    colors=["#666666", RED],
    linewidth=2.6,
    markersize=8,
    ax=ax
)

# Create custom y axis tick labels
ylabels = [f"H: {round(h)}, AB: {round(ab)}" for h, ab in zip(df["H"].values, df["AB"].values)]
ylabels = list(reversed(ylabels))

# Put the labels for the y axis in the mid of the original location of the tick marks.
ax.set_yticklabels(ylabels, ha="right")

# Create legend
handles = [
    Patch(label="Non-hierarchical", facecolor="#666666"),
    Patch(label="Hierarchical", facecolor=RED),
    Line2D([0], [0], label="Mean probability", ls="--", color="black", alpha=0.5)
]

legend = ax.legend(handles=handles, loc=4, fontsize=14, frameon=True, framealpha=0.8);
../_images/notebooks_hierarchical_binomial_bambi_45_0.png

One of the first things one can see is that not only the center of the distributions varies but also their dispersion. Those posteriors that are very wide are associated with players who have batted only once or few times, while tighter posteriors correspond to players who batted several times.

Players who have extreme empirical proportions have similar extreme posteriors under the non-hierarchical model. However, under the hierarchical model, these distributions are now shrunk towards the global mean. Extreme values are very unlikely under the hierarchical model.

And finally, paraphrasing Eric, there’s nothing ineherently right or wrong about shrinkage and hierarchical models. Whether this is reasonable or not depends on our prior knowledge about the problem. And to me, after having seen the hit rates of the other players, it is much more reasonable to shrink extreme posteriors based on very few data points towards the global mean rather than just let them concentrate around 0 or 1.

[22]:
%load_ext watermark
%watermark -n -u -v -iv -w
Last updated: Sun Feb 13 2022

Python implementation: CPython
Python version       : 3.8.5
IPython version      : 7.29.0

matplotlib: 3.5.1
bambi     : 0.7.1
sys       : 3.8.5 (default, Sep  4 2020, 07:30:14)
[GCC 7.3.0]
numpy     : 1.20.0
arviz     : 0.11.4

Watermark: 2.1.0

Footnotes#

  1. By default, the .predict() method obtains the posterior for the mean of the likelihood distribution. This mean would be \(np\) for the Binomial family. However, since \(n\) varies from observation to observation, it returns the value of \(p\), as if it was a Bernoulli family.

  2. .predict()just appends _mean to the name of the response to indicate it is the posterior of the mean.