import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az
import warnings
warnings.filterwarnings("ignore")
RANDOM_STATE = 42
weekly = pd.read_csv("../data/processed/m5_weekly.csv", parse_dates=["week_start"])
weekly["series"] = weekly.store_id + "_" + weekly.cat_id
weekly["hi_snap"] = (weekly.snap_days >= 5).astype(int)
df = weekly.dropna(subset=["units"]).copy()
df["ly"] = np.log(df.units.clip(lower=1))
print(f"{df.series.nunique()} series | {len(df):,} rows | "
f"{df.hi_snap.mean():.1%} of weeks are high-SNAP")13 Lab 11 Part 2: Partial Pooling — and the Question That Comes Before It
Run this one live: the executable lab.
13.1 What we are doing
A hierarchical model lets thirty series share information: each gets its own effect, drawn from a common distribution the data also estimates. Data-poor series borrow from data-rich ones, and nobody has to choose the weights.
That is genuinely powerful, and today you will watch it work. You will also watch it produce a confidently wrong answer — on the same data, with clean diagnostics, textbook shrinkage, and no warning of any kind.
The difference between the two runs is not the code. It is a decision made before any code was written: which series belong in the same hierarchy.
Suggested time: ~45 minutes · You will need: pandas, numpy, matplotlib, pymc, arviz, and data/processed/m5_weekly.csv. Six model fits, most under 30 seconds.
13.2 Setup
Our question: does a high-SNAP week lift units, and by how much in each store?
We will need one helper repeatedly — the SNAP effect actually observed in a set of weeks, with no model involved. It is just the difference in mean log units between high-SNAP weeks and the rest.
def empirical_effect(frame, series_name):
"""Observed high-SNAP effect on log units for one series."""
z = frame[frame.series == series_name]
if z.hi_snap.nunique() < 2:
return np.nan
return z.ly[z.hi_snap == 1].mean() - z.ly[z.hi_snap == 0].mean()
Walkthrough
30 series | 8,310 rows | 22.1% of weeks are high-SNAPThe question: does a high-SNAP week lift units, and by how much in each store?
And note this helper, because it is what makes the whole lab possible. It computes the SNAP effect actually observed in a set of weeks — the difference in mean log units between high-SNAP weeks and the rest — with no model involved. Later we will hold weeks back, fit a model without them, and then use this function to ask what really happened in those weeks.
Most of the time you cannot grade a model this way. Today you can, because we are going to create the missing data ourselves.
13.3 Step 1 — Build the imbalance
Real panels are lopsided: some series have years of history, others are new stores or newly re-categorized products with a handful of usable weeks. We create that imbalance deliberately so we can check the answer later.
Thin every third series down to 8 weeks. The weeks we remove are not thrown away — they become a held-out set that the model never sees, and that we will use in Step 3 to grade it.
codes, uniq = pd.factorize(df.series)
df["idx"] = codes
J = len(uniq)
rng = np.random.RandomState(RANDOM_STATE)
keep = []
for j in range(J):
rows = df.index[df.idx == j].tolist()
n = 8 if j % 3 == 0 else len(rows)
keep += list(rng.choice(rows, size=min(n, len(rows)), replace=False))
train = df.loc[sorted(keep)].copy()
heldout = df.drop(index=keep).copy()
thin = np.array([j % 3 == 0 for j in range(J)])
print(f"train {len(train):,} rows | held-out {len(heldout):,} rows")
print(f"{thin.sum()} thin series (8 weeks each), {(~thin).sum()} full series")Here is the model. It is the one from lecture, written non-centered, and it takes a flag so we can reuse it for the unpooled fit.
def fit(data, kind, n_series, centered=False, target_accept=0.9):
"""kind: 'unpooled' or 'hier'. Returns the InferenceData."""
i = data.idx.values
x = data.hi_snap.values.astype(float)
y = data.ly.values
with pm.Model():
a = pm.Normal("a", y.mean(), 1, shape=n_series)
s = pm.Exponential("s", 1.0)
if kind == "unpooled":
b = pm.Normal("b", 0, 0.5, shape=n_series)
else:
mu_b = pm.Normal("mu_b", 0, 0.5)
sd_b = pm.HalfNormal("sd_b", 0.25)
if centered:
b = pm.Normal("b", mu_b, sd_b, shape=n_series)
else:
z = pm.Normal("z", 0, 1, shape=n_series)
b = pm.Deterministic("b", mu_b + z * sd_b)
pm.Normal("obs", mu=a[i] + b[i] * x, sigma=s, observed=y)
return pm.sample(1000, tune=1500, chains=4, progressbar=False,
random_seed=RANDOM_STATE, target_accept=target_accept)
def b_means(idata, n_series):
return az.summary(idata, var_names=["b"], ci_kind="hdi",
ci_prob=0.94)["mean"].values
Walkthrough
train 5,620 rows | held-out 2,690 rows
10 thin series (8 weeks each), 20 full seriesReal panels are lopsided — some series have years of history, others are new stores with a handful of weeks. We manufacture that imbalance: every third series cut down to 8 weeks.
The weeks we remove are not discarded. They go into heldout, the model never sees them, and in Step 3 they become the answer key.
Remember the rule: j % 3 == 0. Every third series. It looks arbitrary and it is about to matter enormously.
13.4 Step 2 — Fit it, and admire the shrinkage
Fit both models on all thirty series and measure how far each estimate moved.
idata_u30 = fit(train, "unpooled", J)
idata_h30 = fit(train, "hier", J)
bu30, bh30 = b_means(idata_u30, J), b_means(idata_h30, J)
def diag(idata, label):
su = az.summary(idata, var_names=["b"], ci_kind="hdi", ci_prob=0.94)
print(f" {label:<14s} r_hat {su.r_hat.max():.3f} "
f"ess {su.ess_bulk.min():5.0f} "
f"divergences {int(idata.sample_stats.diverging.sum()):4d}")
print("Diagnostics")
diag(idata_u30, "unpooled")
diag(idata_h30, "hierarchical")move_thin = np.abs(bh30[thin] - bu30[thin]).mean()
move_full = np.abs(bh30[~thin] - bu30[~thin]).mean()
print("\n unpooled spread after pooling mean move")
print(f" thin (8w) {bu30[thin].std():.3f} "
f"{bh30[thin].std():.3f} {move_thin:.3f}")
print(f" full (277w) {bu30[~thin].std():.3f} "
f"{bh30[~thin].std():.3f} {move_full:.3f}")
print(f"\n ratio: thin series moved {move_thin / move_full:.1f}x further")
g = az.summary(idata_h30, var_names=["mu_b", "sd_b"], ci_kind="hdi", ci_prob=0.94)
print(f" mu_b {g.loc['mu_b','mean']:.3f} sd_b {g.loc['sd_b','mean']:.3f}")Everything a hierarchical model is supposed to do, it did. Diagnostics are clean. The data-poor series moved an order of magnitude further toward the group mean than the data-rich ones, and nobody told the model which was which.
The shrinkage plot makes it visual.
fig, ax = plt.subplots(figsize=(7, 5))
for j in range(J):
c = "crimson" if thin[j] else "steelblue"
ax.plot([0, 1], [bu30[j], bh30[j]], "-o", color=c, alpha=0.7, ms=4)
ax.axhline(g.loc["mu_b", "mean"], ls="--", c="gray", lw=1)
ax.set_xticks([0, 1]); ax.set_xticklabels(["unpooled", "hierarchical"])
ax.set_ylabel("estimated high-SNAP effect (log units)")
ax.set_title("Shrinkage: thin series (red) collapse onto the group mean")
plt.tight_layout(); plt.show()
Before you go on
Write down, in one sentence, what you would now tell the business. Keep it — Step 3 will ask whether it was right.
Walkthrough
Diagnostics
unpooled r_hat 1.004 ess 6637 divergences 0
hierarchical r_hat 1.001 ess 5203 divergences 0Both models sampled cleanly. R-hat near one, effective sample sizes in the thousands, not one divergence. If you were checking your work the way Lecture 10 taught you, you would stop here satisfied.
unpooled spread after pooling mean move
thin (8w) 0.104 0.000 0.110
full (277w) 0.012 0.001 0.010
ratio: thin series moved 11.5x further
mu_b 0.011 sd_b 0.008And this is textbook. The thin series — 8 weeks each — moved 11.5 times further toward the group mean than the full series did. Nobody told the model which series were starved. It worked that out from how much each one had to say.
The picture makes it obvious: the red lines, the thin series, collapse onto the dashed group mean. The blue lines barely move.
Try it
Pause the video and write one sentence: what would you now tell the business about SNAP weeks? Keep it where you can see it. Step 3 is going to ask whether it was right. Before we grade it, look at one number I have not commented on. sd_b is 0.008.
That is the model’s estimate of how much stores differ from one another. Eight thousandths of a log unit — which is to say, the model has concluded that all thirty of these series respond to SNAP essentially identically.
Hold that up against what you know about the world. Thirty series spanning three states, three product categories, stores of different sizes. Is “they are all the same” plausible?
It is not a diagnostic. Nothing flags it. But it is a claim about the world, and claims about the world can be sanity-checked in a way that R-hat cannot.
13.5 Step 3 — Grade it against the weeks it never saw
We thinned those ten series ourselves. The weeks we removed are sitting in heldout, and each thinned series has hundreds of them. So we can compute the effect that actually occurred and compare.
thin_names = [uniq[j] for j in range(J) if thin[j]]
truth = np.array([empirical_effect(heldout, s) for s in thin_names])
est_u = np.array([bu30[j] for j in range(J) if thin[j]])
est_h = np.array([bh30[j] for j in range(J) if thin[j]])
out = pd.DataFrame({"series": thin_names, "unpooled": est_u.round(3),
"hierarchical": est_h.round(3), "held_out_truth": truth.round(3)})
print(out.to_string(index=False))
rmse = lambda e: np.sqrt(np.nanmean((e - truth) ** 2))
print(f"\n RMSE vs truth: unpooled {rmse(est_u):.3f} "
f"hierarchical {rmse(est_h):.3f}")
print(f" mean held-out effect across these series: {np.nanmean(truth):.3f}")Read that table before reading on.
The hierarchical model reported an effect near mu_b for every one of these series. The held-out weeks say the true effects average close to 0.09 — several times larger. The model did not merely lose precision on these series; it moved them all to the wrong place, and reported tight intervals around it.
And notice what the RMSE column says: pooling barely improved on the unpooled fit at all. All that shrinkage bought us essentially nothing.
The diagnostics could never have caught this
\(\hat{R}\), ESS and divergences ask whether the sampler explored the posterior of the model you wrote. They cannot ask whether that model was the right one. Every diagnostic here passed.
Walkthrough
series unpooled hierarchical held_out_truth
CA_1_FOODS 0.294 0.011 0.087
CA_2_FOODS 0.017 0.011 0.011
CA_3_FOODS 0.035 0.011 0.105
CA_4_FOODS 0.132 0.011 0.039
TX_1_FOODS 0.089 0.011 0.082
TX_2_FOODS 0.171 0.011 0.067
TX_3_FOODS 0.132 0.011 0.074
WI_1_FOODS -0.015 0.011 0.033
WI_2_FOODS 0.289 0.011 0.209
WI_3_FOODS 0.015 0.011 0.177
RMSE vs truth: unpooled 0.102 hierarchical 0.097
mean held-out effect across these series: 0.089Read the middle column first. 0.011, ten times. The hierarchical model reported essentially the same effect for every one of these ten series.
Now read the last column. The held-out weeks — hundreds of them per series, weeks the model never saw — say the true effects average 0.089, and range from 0.011 up to 0.209.
The model was not merely imprecise. It moved all ten series to a number roughly eight times too small, and reported tight intervals around it.
And look at the RMSE line, because it takes away the last consolation. Unpooled 0.102, hierarchical 0.097. All of that shrinkage — the beautiful collapsing lines, the 11.5× ratio — bought a five percent improvement over just fitting each series separately.
Now the sentence that matters most in this lecture. The diagnostics could never have caught this. R-hat, ESS and divergences ask whether the sampler explored the posterior of the model you wrote. They cannot ask whether that was the right model to write. Every diagnostic passed, and every diagnostic was answering a different question from the one you cared about.
13.6 Step 4 — Find out what went wrong
Two things are wrong, and one of them is invisible until you look.
First, look at which series we thinned.
print("thinned series:", thin_names)print("\nEmpirical high-SNAP effect by category (all data, no model):")
for cat, grp in df.groupby("cat_id"):
effs = [empirical_effect(grp, s) for s in grp.series.unique()]
print(f" {cat:<10s} n={len(effs):2d} mean {np.mean(effs):+.3f} "
f"sd {np.std(effs):.3f}")Two findings, and they compound.
The groups were not exchangeable. SNAP is a food assistance program. FOODS responds; HOBBIES and HOUSEHOLD barely move. Putting all thirty in one hierarchy told the model these are thirty interchangeable draws from one distribution — so it estimated a group mean near zero, dominated by the twenty series with no effect, and pulled the FOODS series down to it.
And j % 3 == 0 selected exactly the FOODS series. The series are ordered store-major with three categories each, so every third one is a FOODS series. We starved precisely the ten series that had something to say, then pooled them toward twenty that did not.
Exchangeability is a modeling assumption, and it is yours to defend
Groups are exchangeable when, before seeing the data, you have no reason to expect one to differ from another in a particular direction. Ten FOODS series in ten stores: plausible. FOODS against HOBBIES: you knew the answer before you looked.
Nothing in PyMC checks this. No diagnostic reports it. It is subject-matter knowledge, applied before the model is written.
# ---- your turn: uncomment and answer in the strings ----
# whose_fault = "..." # <- would a bigger sample have fixed this? why or why not?
# print(whose_fault)
Walkthrough
Two things went wrong, and they compound.
thinned series: ['CA_1_FOODS', 'CA_2_FOODS', 'CA_3_FOODS', 'CA_4_FOODS', 'TX_1_FOODS',
'TX_2_FOODS', 'TX_3_FOODS', 'WI_1_FOODS', 'WI_2_FOODS', 'WI_3_FOODS']Every thinned series is a FOODS series. All ten. That is not a coincidence and it is not bad luck — the panel is ordered store-major with three categories each, so j % 3 == 0 selects the same category every time.
Empirical high-SNAP effect by category (all data, no model):
FOODS n=10 mean +0.089 sd 0.058
HOBBIES n=10 mean +0.006 sd 0.010
HOUSEHOLD n=10 mean +0.011 sd 0.014And there it is. SNAP is a food assistance program. FOODS responds at +0.089. HOBBIES and HOUSEHOLD are at +0.006 and +0.011 — indistinguishable from nothing. So put the two together, because neither alone is fatal.
We told the model that thirty series were thirty interchangeable draws from one distribution. They are not: twenty of them have no effect and ten do. The group mean therefore landed near zero, dominated by the twenty. And then — the second error — we starved precisely the ten series that had something to say, so those ten had almost no data with which to resist being pulled toward a mean built from series unlike them.
That is why sd_b came out at 0.008. Averaged across a group where two thirds genuinely have no effect, “stores are all the same” is the best summary of a badly-posed question.
Try it
Pause and answer the question in the code cell: would a bigger sample have fixed this?
No — and this is worth being clear about, because “get more data” is the reflex. More weeks per series would sharpen each estimate, and the hierarchy would still be pulling FOODS series toward a mean computed from HOBBIES and HOUSEHOLD. The error is in the structure we imposed, and no quantity of data corrects a structural mistake. It just makes you more confident about it.
Exchangeability is a modeling assumption and it is yours to defend. Groups are exchangeable when, before seeing the data, you have no reason to expect one to differ from another in a particular direction. Ten FOODS series across ten stores: plausible. FOODS against HOBBIES: you knew the answer before you looked.
Nothing in PyMC checks this. No diagnostic reports it. It is subject-matter knowledge, applied before the model is written.
13.7 Step 5 — Do it on the exchangeable set
Restrict to the ten FOODS series — the set we can actually defend pooling — and rebuild the same imbalance: thin every other one to 8 weeks.
foods = df[df.cat_id == "FOODS"].copy()
codes_f, uniq_f = pd.factorize(foods.series)
foods["idx"] = codes_f
Jf = len(uniq_f)
rng = np.random.RandomState(RANDOM_STATE)
keep_f = []
for j in range(Jf):
rows = foods.index[foods.idx == j].tolist()
n = 8 if j % 2 == 0 else len(rows)
keep_f += list(rng.choice(rows, size=min(n, len(rows)), replace=False))
train_f = foods.loc[sorted(keep_f)].copy()
heldout_f = foods.drop(index=keep_f).copy()
thin_f = np.array([j % 2 == 0 for j in range(Jf)])
print(f"{Jf} FOODS series | {thin_f.sum()} thinned to 8 weeks")
print(f"train {len(train_f):,} rows | held-out {len(heldout_f):,} rows")idata_uf = fit(train_f, "unpooled", Jf)
idata_hf = fit(train_f, "hier", Jf)
buf, bhf = b_means(idata_uf, Jf), b_means(idata_hf, Jf)
gf = az.summary(idata_hf, var_names=["mu_b", "sd_b"], ci_kind="hdi", ci_prob=0.94)
print(gf.round(3).to_string())mt = np.abs(bhf[thin_f] - buf[thin_f]).mean()
mf = np.abs(bhf[~thin_f] - buf[~thin_f]).mean()
print(" unpooled spread after pooling mean move")
print(f" thin (8w) {buf[thin_f].std():.3f} "
f"{bhf[thin_f].std():.3f} {mt:.3f}")
print(f" full (277w) {buf[~thin_f].std():.3f} "
f"{bhf[~thin_f].std():.3f} {mf:.3f}")
print(f"\n ratio: thin moved {mt / mf:.1f}x further")Compare that middle column to Step 2’s. There, everything collapsed onto the group mean. Here the thin series collapse but the full series keep most of their spread — the model discarded variation that 8 weeks could not support and kept variation that 277 weeks could. That is what partial pooling is supposed to look like.
Now grade it the same way.
tn = [uniq_f[j] for j in range(Jf) if thin_f[j]]
truth_f = np.array([empirical_effect(heldout_f, s) for s in tn])
eu = np.array([buf[j] for j in range(Jf) if thin_f[j]])
eh = np.array([bhf[j] for j in range(Jf) if thin_f[j]])
print(pd.DataFrame({"series": tn, "unpooled": eu.round(3),
"hierarchical": eh.round(3),
"held_out_truth": truth_f.round(3)}).to_string(index=False))
r_u = np.sqrt(np.nanmean((eu - truth_f) ** 2))
r_h = np.sqrt(np.nanmean((eh - truth_f) ** 2))
print(f"\n RMSE vs truth: unpooled {r_u:.3f} hierarchical {r_h:.3f}"
f" ({r_u / r_h:.1f}x better)")That is borrowing strength, and now it is real. The same code, the same priors, the same diagnostics — and the thin series are estimated well because the other five exist. The only thing that changed was which series were allowed into the hierarchy.
# ---- your turn: uncomment and fill in the comparison ----
# print("30-series pooling improved RMSE by ___%") # <- from Step 3's two numbers
# print("FOODS-only pooling improved RMSE by ___%") # <- from the two above
Walkthrough
10 FOODS series | 5 thinned to 8 weeks
train 1,425 rows | held-out 1,345 rowsSame code. Same priors. Same thinning trick — every other series cut to 8 weeks. The only change is which series are allowed into the hierarchy.
mean sd hdi94_lb hdi94_ub ess_bulk ess_tail r_hat
mu_b 0.077 0.035 0.016 0.15 1315 1619 1.00
sd_b 0.066 0.038 0.012 0.15 756 1016 1.00mu_b is 0.077 — a real SNAP effect, where before it was 0.011. And sd_b is 0.066, where before it was 0.008. Eight times larger.
That second number is the interesting one. The model now says stores genuinely differ from one another, because on a defensible group they do. The earlier 0.008 was not a finding about stores; it was an artifact of asking about the wrong group.
unpooled spread after pooling mean move
thin (8w) 0.108 0.010 0.120
full (277w) 0.055 0.038 0.015
ratio: thin moved 8.2x furtherCompare the middle column with Step 2’s. There, everything collapsed — thin spread went to 0.000, full to 0.001. Here the thin series collapse to 0.010 but the full series keep most of their spread: 0.055 down to 0.038, about seventy percent retained.
That is what partial pooling is supposed to look like. The model discarded variation that 8 weeks could not support, and kept variation that 277 weeks could. In Step 2 it discarded everything, which should have been the tell.
series unpooled hierarchical held_out_truth
CA_1_FOODS 0.292 0.097 0.087
CA_3_FOODS 0.103 0.080 0.102
TX_1_FOODS 0.052 0.074 0.082
TX_3_FOODS 0.211 0.090 0.072
WI_2_FOODS 0.338 0.101 0.207
RMSE vs truth: unpooled 0.126 hierarchical 0.049 (2.6x better)Now read the middle column against the last one. 0.097 against 0.087. 0.080 against 0.102. 0.090 against 0.072. These are close, and each one came from a series with eight weeks of data. RMSE against the truth: 0.126 unpooled, 0.049 hierarchical. Two and a half times better.
That is borrowing strength, and now it is real. Same code, same priors, same diagnostics — and the thin series are estimated well because the other five exist. The only thing that changed was which series were allowed into the hierarchy.
Try it
Fill in the comparison in the last cell: what percentage did pooling improve RMSE in Step 3, and what percentage here?
13.8 Step 6 — The trap in the code
One more failure, and this one is mechanical. The natural way to write the hierarchy reads exactly like the equation:
b = pm.Normal("b", mu_b, sd_b, shape=J) # the centered formFit it that way and watch the diagnostics.
idata_cent = fit(train_f, "hier", Jf, centered=True)
print("Diagnostics")
diag(idata_cent, "centered")
diag(idata_hf, "non-centered")Same model. Same data. Same answer, if you could sample it — but the centered form makes the sampler walk a funnel: when sd_b is small the b_j must crowd into a narrow neck, and a sampler tuned for the wide part cannot get in.
Which alarm fired
Look at the two columns. \(\hat{R}\) is near the threshold and ESS clears 400 — on those two numbers alone you might have shipped it. The divergence count is unambiguous.
On hierarchical models, divergences are the sensitive instrument. Do not wait for \(\hat{R}\).
Walkthrough
One more failure, and this one is mechanical rather than conceptual.
The natural way to write a hierarchy reads exactly like the equation: b = pm.Normal("b", mu_b, sd_b, shape=J). Everybody writes it that way first. It is called the centered form.
Diagnostics
centered r_hat 1.044 ess 413 divergences 127
non-centered r_hat 1.002 ess 1728 divergences 1Same model. Same data. The same posterior, if you could sample it — but the centered form makes the sampler walk a funnel: when sd_b is small the individual b_j must crowd into a narrow neck, and a sampler tuned for the wide part of the funnel cannot get into the neck.
Now read the three numbers as three different kinds of warning, because that is the lesson. ESS passes. 413, thirteen over the threshold. You would not look twice.
R-hat fails — 1.044 against 1.01. But be honest about what happens in practice when a number is a little over a line: people rerun it, or call it close enough, or decide the threshold is conservative.
127 divergences against a threshold of exactly zero. There is no version of “close to zero” that gets you to a hundred and twenty-seven. It admits no conversation.
So the rule: on hierarchical models, divergences fire first and least ambiguously. Do not wait for R-hat to decide for you.
13.9 Before you leave
Discuss (the failure): Step 2’s model had clean diagnostics, textbook shrinkage, and a defensible-sounding conclusion. Step 3 showed it was wrong. What kind of check caught it, and why could no amount of sampler diagnostics have done so?
Discuss (the assumption): You are asked to pool weekly demand across the ten stores in one chain. Name one variable that, if it differed across those stores, would break exchangeability the way category did here — and say how you would handle it without abandoning the hierarchy.
Discuss (the judgment): In Step 5,
sd_bcame out around 0.07 — stores genuinely differ. The unpooled fit put the range at roughly 5% to 33%. Given both numbers, what do you actually recommend the business do for the next SNAP week?Discuss (the thread): This is the sixth time this semester that a model has been improved by pulling estimates toward a common value. Name the other five, and say what is different about the shrinkage strength here.
Solutions posted to Canvas after class.
Walkthrough
Try it
Four questions for the board.
The failure. Step 2 had clean diagnostics, textbook shrinkage and a defensible conclusion. Step 3 showed it was wrong. What kind of check caught it, and why could no amount of sampler diagnostics have done so? model never saw can audit the model.*
The assumption. You are asked to pool weekly demand across the ten stores in one chain. Name one variable that, if it differed across those stores, would break exchangeability the way category did here — and say how you would handle it without abandoning the hierarchy.
as a predictor, or nest the hierarchy — stores within regions — rather than giving up on pooling.*
The judgment. In Step 5, sd_b came out around 0.066, so stores genuinely differ. The unpooled fit put the range at roughly 5% to 33%. Given both numbers, what do you actually recommend the business do for the next SNAP week?
The thread. This is the sixth time this semester a model improved by pulling estimates toward a common value. Name the other five, and say what is different about the shrinkage strength here.
and the difference is that here the strength is* estimated from the data rather than chosen by you and tuned on a validation set. sd_b is the shrinkage parameter, and the model infers it.
Solutions for the coded parts go up on Canvas after the deadline. Next week: Bayesian linear regression — and a preprocessing step you have used all semester that silently destroys the answer.
13.10 What your run should produce
Your numbers will differ where a model is stochastic; these are the values this lab produces on the pinned environment.
random_seed=42 throughout; six fits, whole lab ≈ 3 minutes.
| Quantity | Value |
|---|---|
| Panel | 30 series · 8,310 rows · 22.1% high-SNAP weeks |
| Split | train 5,620 · held-out 2,690 · 10 thin (8 wks) + 20 full |
| Step 2 diagnostics | unpooled 1.004 / 6,637 / 0 · hierarchical 1.001 / 5,203 / 0 — all clean |
| Step 2 shrinkage | thin 0.104 → 0.000 (move 0.110) · full 0.012 → 0.001 (move 0.010) · ratio 11.5× |
| Step 2 hierarchy | mu_b 0.011 · sd_b 0.008 ← implausibly small; the substantive tell |
| Step 3 grading | hierarchical reports 0.011 for all ten; held-out truth averages 0.089 |
| Step 3 RMSE | unpooled 0.102 · hierarchical 0.097 — pooling bought ~5% |
| Step 4 by category | FOODS +0.089 (sd 0.058) · HOBBIES +0.006 · HOUSEHOLD +0.011 |
| Step 4 thinning | j % 3 == 0 selects all ten FOODS series |
| Step 5 hierarchy | mu_b 0.077 [0.016, 0.150] · sd_b 0.066 — 8× Step 2’s |
| Step 5 shrinkage | thin 0.108 → 0.010 · full 0.055 → 0.038 (70% retained) · ratio 8.2× |
| Step 5 grading | unpooled 0.126 · hierarchical 0.049 — 2.6× better |
| Step 6 centered | R-hat 1.044 · ESS 413 · divergences 127 |
| Step 6 non-centered | R-hat 1.002 · ESS 1,728 · divergences 1 |