13 Bayesian Statistics II — Hierarchical Models
Lecture 11, Part 2 of ECON 8310.
13.1 Lecture Outline
- Which Series May Be Pooled
- Partial Pooling
- Fitting It, and the Trap
- Reading the Result
- Key Takeaways
Five parts. Which series may be pooled — the question before the model. Then partial pooling itself, which is one extra line of code. Then fitting it, and a trap that catches nearly everyone the first time. Then reading the result, checked against data the model never saw. And takeaways.
13.2 Which Series May Be Pooled
Before you pool anything, ask what the groups have in common.
13.2.1 The Question
A concrete one, and you have met it before. Does demand rise in weeks when SNAP benefits are disbursed, and by how much — store by store? Homework 3’s memo had you stock up for SNAP weeks store by store, which assumed those differences are real. Today we check.
But one question comes before the model. SNAP is a food assistance program. Measure the raw SNAP-week effect on log units in all thirty series, grouped by category:
| Category | Mean effect | Spread across stores |
|---|---|---|
| FOODS (10 series) | \(+0.089\) | 0.058 |
| HOBBIES (10 series) | \(+0.006\) | 0.010 |
| HOUSEHOLD (10 series) | \(+0.011\) | 0.014 |
These are not thirty versions of the same thing. FOODS responds; the other two barely move. Pooling all thirty would tell the model they are interchangeable.
Here’s a concrete business question, and you’ve met it before. Does demand rise in weeks when SNAP benefits are disbursed, and by how much — store by store?
Homework 3’s memo had you decide which store-categories to stock up for a SNAP week rather than applying one national uplift. That question assumed the per-series differences were real. Today we check. But one question comes before the model, and it costs nothing to ask.
SNAP is a food assistance program. So let’s measure the raw SNAP-week effect — no model, just the difference in mean log units — in all thirty series, grouped by category. FOODS, ten series: plus point-oh-eight-nine, with a spread across stores of point-oh-five-eight. HOBBIES, ten series: plus point-oh-oh-six. HOUSEHOLD, ten series: plus point-oh-one-one.
Look at that first column. FOODS responds. The other two barely move — and that is not a surprise once you say it out loud, because SNAP is food assistance. You could have predicted the sign of that table before computing it.
These are not thirty versions of the same thing. Pooling all thirty would be telling the model they are interchangeable, and the model would believe you.
And I want to stress how cheap that check was. Three group-bys and a subtraction. No model, no sampler, no priors. Two minutes of work.
That is characteristic of the most valuable checks in this course — the leakage checks, the benchmark, the planted noise control in Homework 3. They are nearly free, they happen before the modeling, and they are the ones people skip because they don’t feel like real analysis.
13.2.2 Exchangeability: the Assumption Under Every Hierarchy
Partial pooling rests on one assumption, worth naming before we use it.
Exchangeable
Groups are exchangeable when, before seeing the data, you have no reason to expect one to differ from another in a particular direction. Swapping their labels would not change what you believe.
The ten FOODS series qualify: same category, different stores, and you could not say in advance which store’s SNAP effect should be largest. FOODS against HOBBIES does not — you know before looking that a food-assistance program moves food.
The setup for today. We work with the ten FOODS series, and deliberately thin five of them down to 8 weeks. Five series are data-rich, five are data-poor. That imbalance is the whole point.
The lab pools all thirty anyway, and measures what it costs.
So let’s name the assumption properly, because it sits under every hierarchical model you will ever build.
Groups are exchangeable when, before seeing the data, you have no reason to expect one to differ from another in a particular direction. Swapping their labels would not change what you believe. That’s the test. Not “are they similar?” — but “would relabelling them change my expectations?”
That distinction matters, so let me push on it. Exchangeable does not mean identical. The ten FOODS stores certainly differ — different neighborhoods, different sizes, different customers. Exchangeability doesn’t deny that.
It says something weaker and more useful: that before seeing data, you have no basis for ordering them. You cannot say in advance that the Texas store’s SNAP effect should exceed California’s.
That’s a claim about your knowledge, not about the world. Which is why it’s a judgment call you have to make and defend — and why no diagnostic can make it for you.
Apply it. The ten FOODS series qualify: same category, ten different stores, and you genuinely could not say in advance which store’s SNAP effect should be largest. Swap the store labels and you’d believe exactly the same things.
FOODS against HOBBIES does not qualify. You know before looking that a food-assistance program moves food. There’s a direction, and you knew it in advance.
So for today we work with the ten FOODS series — the defensible set — and we deliberately thin five of them down to eight weeks of data. Five data-rich, five data-poor. That imbalance is the whole point of the exercise.
And in the lab, we pool all thirty anyway — on purpose — and measure exactly what it costs. I’d encourage you to do that one carefully. It is the most instructive forty minutes of the semester, because the wrong model looks completely fine until you check it against held-out data.
13.2.3 Complete Pooling and No Pooling
Within that exchangeable set there are two obvious approaches, and each fails in a different direction.
| What it does | How it fails | |
|---|---|---|
| Complete pooling | One SNAP effect for all ten series. b = pm.Normal("b", 0, 1) |
Denies that stores differ at all. A genuinely different store is averaged away. |
| No pooling | A separate effect per series, estimated independently. shape=10 |
A series with 8 weeks gets an estimate built from 8 weeks. Mostly noise, reported as a finding. |
You have used both already without naming them. Homework 3’s random forest pooled all thirty series into one model with store dummies — close to complete pooling. Homework 1 fitted ARIMA to one series at a time — no pooling. The choice has always been made by convenience; partial pooling makes it a modeling decision and lets the data settle it.
Within that exchangeable set, there are two obvious things you could do, and each fails in a different direction.
Complete pooling: one SNAP effect for all ten series. One parameter. This denies that stores differ at all, so a genuinely different store gets averaged away into the crowd.
No pooling: a separate effect per series, estimated independently — shape=10. Now a series with eight weeks of data gets an estimate built from eight weeks of data. That’s mostly noise, reported as a finding. And here’s the thing — you have used both of these already, without naming them.
Homework 3’s random forest pooled all thirty series into one model with store dummies. That’s close to complete pooling. Homework 1 fitted ARIMA to one series at a time, completely independently. That’s no pooling.
You made that choice twice, and both times you made it by convenience — by what the library made easy. Partial pooling turns it into a modeling decision, with a parameter, and lets the data settle where between the two extremes you should sit.
It’s worth noticing that this is a recurring shape in the course. Something you had been deciding implicitly, by tooling default, gets promoted into an explicit parameter with a posterior.
Lecture 6 did it for which predictors survive — you were choosing by intuition, then LASSO made it a penalty. Lecture 11 Part 1 did it for how fast the level may move. Today it’s how much groups share.
Each time, the gain is the same: the decision becomes visible, arguable, and estimable rather than buried in a default nobody wrote down.
13.3 Partial Pooling
Let each series have its own effect — drawn from a shared distribution the data also estimates.
13.3.1 The Hierarchy Is One Extra Line
Give each series its own effect \(b_j\), but say where those effects come from:
\[b_j \sim \mathcal{N}(\mu_b,\ \sigma_b), \qquad j = 1, \ldots, 10\]
\(\mu_b\) is the typical SNAP effect across FOODS, and \(\sigma_b\) is how much stores genuinely differ. Crucially, both are estimated from the data rather than assumed.
| \(\sigma_b \rightarrow 0\) | The model concludes stores do not differ, and reproduces complete pooling. |
| \(\sigma_b\) large | The model concludes they differ a lot, and approaches no pooling. |
| In between | Each \(b_j\) is pulled toward \(\mu_b\) — how far depending on how much data that series has. |
That last row is the whole idea. A series with 277 weeks barely moves; one with 8 weeks moves a long way. Nobody chose those weights — they fall out of the model.
Here is the whole idea, and it really is one line.
Give each series its own effect, b-j. But then say where those effects come from: b-j is drawn from a Normal with mean mu-b and standard deviation sigma-b. mu-b is the typical SNAP effect across FOODS. sigma-b is how much stores genuinely differ. And critically — both of those are estimated from the data, not assumed by you. Think about the two extremes.
If sigma-b goes to zero, the model has concluded stores don’t differ, and every b-j collapses onto mu-b. That reproduces complete pooling.
If sigma-b is large, the model has concluded they differ a lot, and each b-j is free to sit wherever its own data says. That approaches no pooling.
And in between — which is where you almost always end up — each b-j is pulled toward mu-b, and how far depends on how much data that series has.
That last row is the whole idea, so let me say it plainly. A series with two hundred seventy-seven weeks barely moves, because its own estimate is already trustworthy. A series with eight weeks moves a long way, because its own estimate carries almost no information.
Nobody chose those weights. There is no threshold, no rule, no “if n is less than thirty.” They fall out of estimating sigma-b.
13.3.2 Why This Is the Right Answer, Not a Compromise
Partial pooling sounds like splitting the difference between two bad options. It is better than that, and the reason is worth stating.
When a series has 8 observations, its unpooled estimate is mostly sampling noise. The hierarchy’s pull toward \(\mu_b\) is not a fudge — it is the correct response to a noisy measurement. You already believe this in another form: it is why you would not conclude a coin is biased after three flips.
This is regularization again, for the sixth time. A hierarchical prior shrinks per-group estimates toward a common value exactly as Ridge shrinks coefficients toward zero — and here the shrinkage strength \(\sigma_b\) is estimated rather than tuned by cross-validation.
GAM smoothing penalty (L3) \(\rightarrow\) tree pruning (L4) \(\rightarrow\) reg_lambda (L5) \(\rightarrow\) Ridge and LASSO (L6) \(\rightarrow\) weight decay (L7) \(\rightarrow\) hierarchical priors (today).
I want to address an objection that’s probably forming, because it’s a reasonable one.
Partial pooling sounds like splitting the difference between two bad options. A fudge. Meeting in the middle because you can’t decide. It is better than that, and the reason is worth stating.
When a series has eight observations, its unpooled estimate is mostly sampling noise. Pulling it toward the group mean is not a compromise — it is the correct response to a noisy measurement.
You already believe this in another form. You would not conclude a coin is biased after three flips, even if all three came up heads. You’d pull your estimate toward one-half, because three flips is not much evidence. That’s shrinkage, and you do it instinctively. And this is regularization, for the sixth time this semester.
A hierarchical prior shrinks per-group estimates toward a common value exactly as Ridge shrinks coefficients toward zero. Same operation, different name.
The smoothing penalty in Lecture 3. Cost-complexity pruning in Lecture 4. reg_lambda in Lecture 5. Ridge and LASSO in Lecture 6. Weight decay in Lecture 7. And now the hierarchical prior.
But this version has something none of the other five had. In every previous case, you chose the shrinkage strength by cross-validation — you tried values and picked one. Here, sigma-b is estimated, as part of the model, with its own posterior and its own uncertainty. That makes it the cleanest version of the idea in the course.
13.4 Fitting It, and the Trap
The obvious code does not work, and the reason is geometry rather than statistics.
13.4.1 The Model in PyMC
Non-centered hierarchical model
with pm.Model() as model:
mu_b = pm.Normal("mu_b", 0, 0.5) # typical SNAP effect
sd_b = pm.HalfNormal("sd_b", 0.25) # how much stores differ
z = pm.Normal("z", 0, 1, shape=J) # non-centered
b = pm.Deterministic("b", mu_b + z * sd_b)
a = pm.Normal("a", y.mean(), 1, shape=J) # per-series level
s = pm.Exponential("s", 1.0)
pm.Normal("obs", mu=a[idx] + b[idx]*snap, sigma=s, observed=y)idx maps each row to its series, so b[idx] picks out the right effect for every observation. That indexing pattern is how every hierarchical model in PyMC is written.
Note what is not here. No decision about how much to pool, no weighting scheme, no threshold for “enough data.” Those all come out of estimating \(\sigma_b\).
Here it is in PyMC. Nine lines.
mu-b and sd-b are the group-level parameters — the typical effect and how much stores differ. Then z, a vector of standard Normals. Then b is constructed as mu-b plus z times sd-b. Then a per-series intercept, observation noise, and the likelihood.
The pattern to internalize is b[idx]. idx maps each row to its series number, so b[idx] picks out the right effect for every observation. Every hierarchical model in PyMC is written that way — once you can read that line, you can read most of the applied Bayesian literature. And notice what is not in that code.
There is no decision about how much to pool. No weighting scheme. No threshold for “enough data.” No rule that says series with fewer than twenty observations get shrunk by half. All of that comes out of estimating one parameter, sigma-b. That’s the elegance of it.
One practical note on priors, since we’re looking at the code. sd_b gets a HalfNormal with scale point-two-five — positive-only, because a standard deviation must be, and concentrated on small values because we’re working in log units where an effect of point-two-five is already very large.
That prior is doing real work. With only ten groups, there is not much information about how much groups differ, so the data cannot pin sigma-b down on its own. A vague prior here lets it wander large, which weakens the pooling — you’d drift toward the unpooled answer.
This is the general situation with group-level scale parameters: they are the hardest thing in the model to estimate, and the place where your prior matters most.
13.4.2 The Trap: Write It the Obvious Way and It Breaks
The natural way to write the hierarchy is b = pm.Normal("b", mu_b, sd_b, shape=J) — the centered form, which reads exactly like the equation. On our data it produces:
| Centered | Non-centered | Want | |
|---|---|---|---|
| Divergences | 127 | 1 | 0 |
| \(\hat{R}\) | 1.044 | 1.002 | \(< 1.01\) |
| ESS (bulk) | 413 | 1,728 | \(> 400\) |
Why: when \(\sigma_b\) is small, the \(b_j\) must crowd into a narrow band, and the posterior takes on a funnel shape — wide where \(\sigma_b\) is large, pinched where it is small. A sampler tuned for the wide part cannot navigate the neck, and reports divergences.
Note how the alarms differ in kind. ESS passes — 413, just over the line, and you would not look twice. \(\hat{R}\) at 1.044 fails, but it is the sort of near-miss people talk themselves out of. 127 divergences against a threshold of exactly zero admits no such conversation. On hierarchical models divergences fire first and least ambiguously; do not wait for \(\hat{R}\) to decide for you.
BMCP §4.6.1, Posterior Geometry Matters — the single most common failure in applied hierarchical modeling.
Now, why did I write b as mu_b + z * sd_b instead of the obvious thing?
The obvious thing is pm.Normal("b", mu_b, sd_b, shape=J) — which reads exactly like the equation on slide 8. It is the natural translation. Everybody writes it that way first. And on our data it produces this. Divergences: one hundred twenty-seven, where we want zero. R-hat: one-point-oh-four-four. Fails, but only just enough to argue about. ESS: four hundred thirteen. Passes — barely, thirteen over the line.
Look at those three numbers again, because they do not fail in the same way.
ESS clears its threshold. You would not look twice at four hundred thirteen. R-hat is over the line at one-point-oh-four-four — and be honest about what happens in practice when a number is a little over a threshold: people run it again, or say it is close enough, or quietly widen the threshold.
One hundred twenty-seven divergences against a threshold of exactly zero admits none of that. There is no version of “close to zero” that gets you to a hundred and twenty-seven.
So the rule to carry forward: on hierarchical models, divergences fire first and least ambiguously. Do not wait for R-hat to decide for you. Do not wait for R-hat to tell you something is wrong. It may never say so.
Why does it happen? When sigma-b is small, all the b-j have to crowd into a very narrow band. So the posterior takes on a funnel shape — wide and open where sigma-b is large, pinched to a point where sigma-b is small.
A sampler tunes its step size for the wide part of the funnel. Then it tries to get into the neck with steps that are far too big, fails, and reports a divergence.
The non-centered form samples z from a plain standard Normal — a nice round shape with no funnel — and constructs b afterwards. Same model. Same posterior. Geometry the sampler can actually walk. This is BMCP section four-point-six-point-one, and it is the single most common failure in applied hierarchical modeling.
The instinct when you see divergences is to raise target_accept — force smaller steps. Sometimes that helps a little. Here it mostly doesn’t: at point-nine-five you still get divergences, you just get fewer, and the sampler runs much slower for the privilege.
Reparameterizing is the actual fix, and it costs nothing. Same answer, no divergences, and it runs faster than the centered version did.
So the habit: when a hierarchical model diverges, your first move is to check whether you wrote it centered — not to tune the sampler.
13.5 Reading the Result
Shrinkage, measured — and checked against data the model never saw.
13.5.1 Shrinkage, Measured
Fit both models. Compare each series’ estimated SNAP effect before and after pooling.
| Series group | Unpooled spread | After pooling | Mean move |
|---|---|---|---|
| Thin — 8 weeks (5 series) | 0.108 | 0.010 | 0.120 |
| Full — 277 weeks (5 series) | 0.055 | 0.037 | 0.015 |
The data-poor series moved 8.2 times further toward the group mean than the data-rich ones. Nobody told the model which series were thin. It inferred that from how much each series had to say.
Read the two right-hand columns together. The thin series’ apparent spread of \(0.108\) nearly vanishes, while the full series keep most of theirs (\(0.055 \rightarrow 0.037\)). The hierarchy did not flatten everything: it discarded spread that 8 weeks could not support and kept spread that 277 weeks could.
That is the claim to be skeptical of, so the next slide tests it.
Now let’s see what the pooling actually did. Fit both models — unpooled and hierarchical — and compare each series’ estimate before and after.
The thin series, eight weeks each: unpooled spread of point-one-oh-eight, collapsing to point-oh-one-oh after pooling. Mean movement: point-one-two-zero. The full series, two hundred seventy-seven weeks: spread of point-oh-five-five, going to point-oh-three-seven. Mean movement: point-oh-one-five.
So the data-poor series moved eight times further toward the group mean than the data-rich ones. And nobody told the model which series were thin. It inferred that entirely from how much each series had to say.
Now read those two columns together, because this is what separates partial pooling from just averaging everything.
The thin series’ apparent spread of point-one-oh-eight nearly vanishes — down to point-oh-one. But the full series keep most of theirs — point-oh-five-five down to point-oh-three-seven.
The hierarchy did not flatten everything. It discarded spread that eight weeks could not support, and kept spread that two hundred seventy-seven weeks could.
That’s a strong claim, and you should be skeptical of it. So the next slide tests it against data the model never saw.
13.5.2 Checked Against Data the Model Never Saw
We thinned those five series ourselves, so the weeks we removed are still there. Compare each estimate against the SNAP effect actually observed in the held-out weeks.
| Thin series | Unpooled | Hierarchical | Held-out truth |
|---|---|---|---|
| CA_1_FOODS | 0.286 | 0.097 | 0.087 |
| CA_3_FOODS | 0.105 | 0.082 | 0.102 |
| TX_1_FOODS | 0.052 | 0.075 | 0.082 |
| TX_3_FOODS | 0.219 | 0.091 | 0.072 |
| WI_2_FOODS | 0.338 | 0.103 | 0.207 |
| RMSE against truth | 0.126 | 0.049 | — |
Partial pooling was two and a half times more accurate, beating the unpooled fit on four series out of five. Those estimates were not just noisy — they were confidently wrong, and the confidence came from 8 weeks of data.
Borrowing strength is not a figure of speech: the thin series were estimated better because the other five existed.
Here’s the nice thing about having thinned those series ourselves: the weeks we removed still exist. So we can compute the SNAP effect that actually occurred in them, and grade both models against it. CA_1 FOODS. Unpooled says point-two-eight-six. Hierarchical says point-oh-nine-seven. The truth in the held-out weeks: point-oh-eight-seven. CA_3. Unpooled point-one-oh-five, hierarchical point-oh-eight-two, truth point-one-oh-two.
TX_1. Unpooled point-oh-five-two, hierarchical point-oh-seven-five, truth point-oh-eight-two — note the hierarchical estimate is higher here, and closer. Shrinkage moves estimates both ways. TX_3. Two-one-nine, oh-nine-one, truth oh-seven-two.
WI_2. Three-three-eight, one-oh-three, truth two-oh-seven. This is the one where the unpooled fit was directionally right and the hierarchical was too conservative. RMSE against the truth: unpooled point-one-two-six. Hierarchical point-oh-four-eight.
Partial pooling was two and a half times more accurate, and it beat the unpooled fit on four of the five series.
And notice the character of the unpooled errors. They were not merely noisy — they were confidently wrong. Point-three-three-eight is a specific, precise-looking claim that a store’s SNAP uplift is thirty-four percent. It came from eight weeks of data, and the truth was about twenty-one percent.
Borrowing strength is not a figure of speech. The thin series were estimated better because the other five existed.
13.5.3 What the Model Actually Concluded
The estimated hierarchy: \(\mu_b = 0.077\) with a 94% HDI of \([0.016,\ 0.150]\), and \(\sigma_b = 0.066\).
In words: a SNAP week lifts FOODS units by roughly 8%, and stores genuinely differ around that by about 7 percentage points. Both numbers are answers, and the second one is the one no point estimate would have given you.
It does not license the store-level policy the unpooled fit suggested. An analyst reading the unpooled numbers would report uplifts from \(5\%\) to \(40\%\) and stock accordingly. The held-out weeks say that range was mostly noise — but \(\sigma_b = 0.066\) says it is not all noise either.
The defensible recommendation: plan on an 8% FOODS uplift everywhere, allow real but modest store variation around it, and do not act on a store’s own estimate until that store has the weeks to support one.
So what’s the business answer?
The estimated hierarchy: mu-b of point-oh-seven-eight, with a ninety-four percent HDI from point-oh-one-seven to point-one-five-zero. And sigma-b of point-oh-six-five.
In words: a SNAP week lifts FOODS units by roughly eight percent, and stores genuinely differ around that by about seven percentage points.
Both of those are answers, and the second one is the one no point estimate would ever have given you. “How much do stores differ?” is a question with a number attached now. But be careful what you conclude from it.
This does not license the store-level policy the unpooled fit suggested. An analyst reading those unpooled numbers would report uplifts from five percent to forty percent and build a differentiated stocking plan on it. The held-out weeks say most of that range was noise.
But — and this is the part that requires care — sigma-b of point-oh-six-five says it is not all noise either. Stores do differ. Just far less than eight weeks of data made it look.
So the defensible recommendation is a three-part sentence: plan on an eight percent FOODS uplift everywhere; allow for real but modest store variation around it; and do not act on any individual store’s own estimate until that store has the weeks to support one.
That’s a more nuanced answer than either extreme, it’s backed by held-out data, and it is exactly what the hierarchy was for.
13.6 Key Takeaways
What hierarchy buys, and what it costs.
13.6.1 Lecture 11 Part 2: Key Takeaways
Ask what may be pooled before you pool it. Groups must be exchangeable. FOODS and HOBBIES are not, and no amount of good sampling fixes a hierarchy built over the wrong set.
Complete pooling denies groups differ; no pooling lets a group with 8 observations report noise as a finding. You have used both without naming them.
Partial pooling gives each group its own parameter drawn from a shared distribution, and estimates how much groups genuinely differ.
Shrinkage is automatic and proportionate. Our thin series moved 8.1\(\times\) further toward the group mean than data-rich ones, with nobody specifying that — and on held-out weeks it was 2.6\(\times\) more accurate.
This is regularization for the sixth time — and the only version where the shrinkage strength is estimated rather than tuned.
Write it non-centered. The obvious form gave 127 divergences; the reparameterized form gave 1. Same model, different geometry — and divergences objected loudest, against a threshold that allows no argument.
Next: Bayesian linear regression — coefficients as distributions, and what you can ask them.
One, and it’s first for a reason. Ask what may be pooled before you pool it. Groups must be exchangeable. FOODS and HOBBIES are not, and no amount of good sampling fixes a hierarchy built over the wrong set.
Two. Complete pooling denies groups differ; no pooling lets a group with eight observations report noise as a finding. You’ve used both without naming them.
Three. Partial pooling gives each group its own parameter drawn from a shared distribution, and estimates how much groups genuinely differ.
Four. Shrinkage is automatic and proportionate. Our thin series moved eight point one times further toward the group mean than data-rich ones, with nobody specifying that — and on held-out weeks it was two point six times more accurate.
Five. This is regularization for the sixth time, and the only version where the shrinkage strength is estimated rather than tuned.
And six. Write it non-centered. The obvious form gave a hundred and twenty-seven divergences; the reparameterized form gave one. Same model, different geometry — and divergences objected loudest, against a threshold that allows no argument.
Next time: Bayesian linear regression. Coefficients as distributions, what you can ask them — and a preprocessing step you have used all semester that will silently destroy the answer.