{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(hierarchical_partial_pooling)=\n",
"# Hierarchical Partial Pooling\n",
":::{post} Jan 28, 2023\n",
":tags: hierarchical model \n",
":category: intermediate\n",
":author: Vladislavs Dovgalecs, Adrian Seybolt, Christian Luhmann\n",
":::"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Suppose you are tasked with estimating baseball batting skills for several players. One such performance metric is batting average. Since players play a different number of games and bat in different positions in the order, each player has a different number of at-bats. However, you want to estimate the skill of all players, including those with a relatively small number of batting opportunities.\n",
"\n",
"So, suppose a player came to bat only 4 times, and never hit the ball. Are they a bad player?\n",
"\n",
"As a disclaimer, the author of this notebook assumes little to non-existent knowledge about baseball and its rules. The number of times at bat in his entire life is around \"4\".\n",
"\n",
"\n",
"## Data\n",
"\n",
"We will use the baseball [data](http://www.swarthmore.edu/NatSci/peverso1/Sports%20Data/JamesSteinData/Efron-Morris%20Baseball/EfronMorrisBB.txt) {cite:p}`efron1975data`.\n",
"\n",
"\n",
"## Approach\n",
"\n",
"We will use PyMC to estimate the batting average for each player. Having estimated the averages across all players in the datasets, we can use this information to inform an estimate of an additional player, for which there is little data (*i.e.* 4 at-bats).\n",
"\n",
"In the absence of a Bayesian hierarchical model, there are two approaches for this problem:\n",
"\n",
"1. independently compute batting average for each player (no pooling)\n",
"2. compute an overall average, under the assumption that everyone has the same underlying average (complete pooling)\n",
"\n",
"Of course, neither approach is realistic. Clearly, all players aren't equally skilled hitters, so the global average is implausible. At the same time, professional baseball players are similar in many ways, so their averages aren't entirely independent either. \n",
"\n",
"It may be possible to cluster groups of \"similar\" players, and estimate group averages, but using a hierarchical modeling approach is a natural way of sharing information that does not involve identifying *ad hoc* clusters.\n",
"\n",
"The idea of hierarchical partial pooling is to model the global performance, and use that estimate to parameterize a population of players that accounts for differences among the players' performances. This tradeoff between global and individual performance will be automatically tuned by the model. Also, uncertainty due to different number of at bats for each player (*i.e.* information) will be automatically accounted for, by shrinking those estimates closer to the global mean.\n",
"\n",
"For far more in-depth discussion please refer to Stan [tutorial](http://mc-stan.org/documentation/case-studies/pool-binary-trials.html) {cite:p}`carpenter2016hierarchical` on the subject. The model and parameter values were taken from that example."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2022-01-09T12:13:28.055166Z",
"iopub.status.busy": "2022-01-09T12:13:28.054155Z",
"iopub.status.idle": "2022-01-09T12:13:31.399492Z",
"shell.execute_reply": "2022-01-09T12:13:31.398489Z"
}
},
"outputs": [],
"source": [
"import arviz as az\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import pymc as pm\n",
"import pytensor.tensor as pt\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2022-01-09T12:13:31.406489Z",
"iopub.status.busy": "2022-01-09T12:13:31.405489Z",
"iopub.status.idle": "2022-01-09T12:13:31.417488Z",
"shell.execute_reply": "2022-01-09T12:13:31.416489Z"
}
},
"outputs": [],
"source": [
"RANDOM_SEED = 8927\n",
"rng = np.random.default_rng(RANDOM_SEED)\n",
"az.style.use(\"arviz-darkgrid\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can load the dataset using pandas:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2022-01-09T12:13:31.423488Z",
"iopub.status.busy": "2022-01-09T12:13:31.423488Z",
"iopub.status.idle": "2022-01-09T12:13:31.663747Z",
"shell.execute_reply": "2022-01-09T12:13:31.663747Z"
}
},
"outputs": [],
"source": [
"data = pd.read_csv(pm.get_data(\"efron-morris-75-data.tsv\"), sep=\"\\t\")\n",
"at_bats, hits = data[[\"At-Bats\", \"Hits\"]].to_numpy().T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's develop a generative model for these data.\n",
"\n",
"We will assume that there exists a hidden factor (`phi`) related to the expected performance for all players (not limited to our 18). Since the population mean is an unknown value between 0 and 1, it must be bounded from below and above. Also, we assume that nothing is known about global average. Hence, a natural choice for a prior distribution is the uniform distribution.\n",
"\n",
"Next, we introduce a hyperparameter `kappa` to account for the variance in the population batting averages, for which we will use a bounded Pareto distribution. This will ensure that the estimated value falls within reasonable bounds. These hyperparameters will be, in turn, used to parameterize a beta distribution, which is ideal for modeling quantities on the unit interval. The beta distribution is typically parameterized via a scale and shape parameter, it may also be parametrized in terms of its mean $\\mu \\in [0,1]$ and sample size (a proxy for variance) $\\nu = \\alpha + \\beta (\\nu > 0)$.\n",
"\n",
"The final step is to specify a sampling distribution for the data (hit or miss) for every player, using a Binomial distribution. This is where the data are brought to bear on the model."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We could use `pm.Pareto('kappa', m=1.5)`, to define our prior on `kappa`, but the Pareto\n",
"distribution has very long tails. Exploring these properly\n",
"is difficult for the sampler, so we use an equivalent\n",
"but faster parametrization using the exponential distribution.\n",
"We use the fact that the log of a Pareto distributed\n",
"random variable follows an exponential distribution."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2022-01-09T12:13:31.675751Z",
"iopub.status.busy": "2022-01-09T12:13:31.675751Z",
"iopub.status.idle": "2022-01-09T12:13:32.615492Z",
"shell.execute_reply": "2022-01-09T12:13:32.614492Z"
}
},
"outputs": [],
"source": [
"N = len(hits)\n",
"player_names = data[\"FirstName\"] + \" \" + data[\"LastName\"]\n",
"coords = {\"player_names\": player_names.tolist()}\n",
"\n",
"with pm.Model(coords=coords) as baseball_model:\n",
" phi = pm.Uniform(\"phi\", lower=0.0, upper=1.0)\n",
"\n",
" kappa_log = pm.Exponential(\"kappa_log\", lam=1.5)\n",
" kappa = pm.Deterministic(\"kappa\", pt.exp(kappa_log))\n",
"\n",
" theta = pm.Beta(\"theta\", alpha=phi * kappa, beta=(1.0 - phi) * kappa, dims=\"player_names\")\n",
" y = pm.Binomial(\"y\", n=at_bats, p=theta, dims=\"player_names\", observed=hits)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recall our original question was with regard to the true batting average for a player with only 4 at bats and no hits. We can add this as an additional variable in the model"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2022-01-09T12:13:32.628501Z",
"iopub.status.busy": "2022-01-09T12:13:32.627493Z",
"iopub.status.idle": "2022-01-09T12:13:32.630495Z",
"shell.execute_reply": "2022-01-09T12:13:32.630495Z"
}
},
"outputs": [],
"source": [
"with baseball_model:\n",
" theta_new = pm.Beta(\"theta_new\", alpha=phi * kappa, beta=(1.0 - phi) * kappa)\n",
" y_new = pm.Binomial(\"y_new\", n=4, p=theta_new, observed=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The model can visualized like this"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2022-01-09T12:13:32.634494Z",
"iopub.status.busy": "2022-01-09T12:13:32.634494Z",
"iopub.status.idle": "2022-01-09T12:13:32.790964Z",
"shell.execute_reply": "2022-01-09T12:13:32.789964Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pm.model_to_graphviz(baseball_model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now fit the model using MCMC:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2022-01-09T12:13:32.795964Z",
"iopub.status.busy": "2022-01-09T12:13:32.794964Z",
"iopub.status.idle": "2022-01-09T12:14:57.857829Z",
"shell.execute_reply": "2022-01-09T12:14:57.856829Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [phi, kappa_log, theta, theta_new]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
"