How to Smooth Unstable Area Rates with Empirical Bayes in Python
Problem statement
A rate map of small areas tends to rank areas by noise. The Census Bureau's Vintage 2024 county estimates give deaths and population for all 3,144 US counties in 2023, and the highest crude death rate belongs to Loving County, Texas: 3 deaths among 46 residents, 65.2 per 1,000 — seven times the national rate of 9.4. One death more or fewer moves that rate by 21.7 per 1,000.
Empirical Bayes (EB) smoothing pulls each rate towards a prior mean. How far it pulls depends on how little information the area carries. Loving County ends up near the national figure, and a county of a million people barely moves.
It is also easy to get wrong. Checked against deaths from years the smoother never saw, a plain global EB made these county rates less accurate overall, not more: its error rose from 0.92 to 1.30 deaths per 1,000. This guide covers both the smoothing and the check that tells you whether it helped.
Quick answer
import pandas as pd
from esda.smoothing import Empirical_Bayes
est = pd.read_csv("co-est2024-alldata.csv", encoding="latin-1",
dtype={"SUMLEV": str, "STATE": str, "COUNTY": str})
counties = est[est["SUMLEV"] == "050"].copy()
deaths = counties["DEATHS2023"].to_numpy(float)
people = counties["POPESTIMATE2023"].to_numpy(float)
counties["raw"] = deaths / people * 1000
counties["eb"] = Empirical_Bayes(deaths, people).r.ravel() * 1000Measured on the 3,144 counties: Loving County went from 65.2 to 11.3 per 1,000. The range of county rates narrowed from 0.0–65.2 to 2.7–30.1. Counties of 100,000 people or more moved by a median of 0.08%.
Before you map the smoothed column, test it against held-out years (step 6). On this data the answer was not the one you would expect.
Step-by-step solution
1. Start from counts and populations, not from rates
EB needs the numerator and the denominator separately, because the denominator decides how far each rate is pulled. A column of percentages cannot be smoothed.
est = pd.read_csv("co-est2024-alldata.csv", encoding="latin-1",
dtype={"SUMLEV": str, "STATE": str, "COUNTY": str})
counties = est[est["SUMLEV"] == "050"].copy() # 040 rows are states
counties["GEOID"] = counties["STATE"] + counties["COUNTY"]Read the codes as strings. SUMLEV read as an integer becomes 50, and STATE loses its leading zero, which breaks every later join to boundaries.
2. Look at how the spread depends on population
Group the raw rates by population class before smoothing anything:
population counties sd of raw rate
under 1,000 36 11.10
1,000–2,500 111 4.81
2,500–5,000 178 3.65
25,000–100,000 1,000 2.67
100,000 and over 620 2.34
The spread in the smallest class is almost five times the spread in the largest. That is the signature of sampling noise. Real differences in mortality do not grow as the population shrinks.
3. Apply global Empirical Bayes
from esda.smoothing import Empirical_Bayes
eb = Empirical_Bayes(deaths, people)
eb.r.shape # (3144, 1) — a column vector, not a Series
counties["eb"] = eb.r.ravel() * 1000The estimator is short enough to read in the esda source. The prior mean is the pooled rate, total deaths over total population (9.365 per 1,000 here). The between-county variance is estimated by method of moments. Each county then keeps a share of its own rate:
kept = variance / (variance + prior_mean / population)
smoothed = kept * raw_rate + (1 - kept) * prior_mean
Measured, Loving County kept 3.4% of its own rate, a county of about 5,000 people kept 79%, and one of about 100,000 kept 98.7%.
4. Measure the shrinkage by size class
population median kept median shift sd raw sd EB
under 1,000 0.344 25.7% 11.10 2.20
1,000–2,500 0.588 15.2% 4.81 2.71
5,000–10,000 0.852 4.6% 3.41 2.89
25,000–100,000 0.971 0.7% 2.67 2.58
100,000 and over 0.994 0.08% 2.34 2.33
After smoothing, the spread is roughly flat across classes, which is what the method is designed to do. The top 50 counties by raw rate had a median population of 3,948; the top 50 by EB rate had 13,706, and only 27 counties appeared in both lists.
5. Add neighbours with Spatial Empirical Bayes
The spatial version uses each county's neighbours as its prior instead of the national mean:
import geopandas as gpd
import libpysal
from esda.smoothing import Spatial_Empirical_Bayes
shapes = gpd.read_file("cb_2023_us_county_500k.zip")[["GEOID", "geometry"]]
g = shapes.merge(counties, on="GEOID").to_crs(5070).reset_index(drop=True)
w = libpysal.weights.Queen.from_dataframe(g, use_index=True)
g["seb"] = Spatial_Empirical_Bayes(g["DEATHS2023"].to_numpy(float),
g["POPESTIMATE2023"].to_numpy(float), w).r.ravel() * 1000Build the weights with use_index=False and esda (2.10) refuses to run:
ValueError: The `id_order` of `w` must be set to align with the order of `e` and `b`.
The check is deliberate. Weights whose order is not declared could silently pair each county with another county's neighbours. With use_index=True, or with KNN.from_dataframe, the order is set and the smoother runs in about 0.03 s. Loving County comes out at 10.1.
6. Validate against years the smoother never saw
This is the step most tutorials skip. Build a reference rate for each county from 2021, 2022 and 2024, rescale it to the 2023 national level, and score each 2023 estimate against it:
RMSE, deaths per 1,000 raw global EB spatial EB EB by size class
under 2,500 2.52 4.49 4.54 3.74
2,500–5,000 1.60 2.02 1.97 1.45
5,000–10,000 1.12 1.32 1.38 1.08
100,000 and over 0.33 0.33 0.33 0.33
all counties 0.92 1.30 1.32 1.08
Smoothing made the small counties worse. The reason is in the bias. Global EB under-estimated counties under 2,500 people by 2.03 deaths per 1,000 on average, because their rate in the other years was 13.6, not the national 9.4. Small rural counties are older, and older populations have higher crude death rates. The smoother assumed every county was average and pulled real differences away.
7. Give the prior a better mean
Stratify the prior so that each county is shrunk towards counties like it. Running EB separately within population classes removed the bias (+0.07 in the smallest class, against −2.03 before). It beat the raw rate in the 2,500–5,000 class (1.45 against 1.60) and the 5,000–10,000 class (1.08 against 1.12).
It still lost in the smallest class. Loving County recorded 3, 3, 3 and 2 deaths in 2021 to 2024, so its other-year rate is 49.9. That is a persistent pattern, not a one-off. A smoother that pulls it to 16 is removing information.
8. Publish the smoothed and the raw rate together
Keep both columns and label the map with the method: "Empirical Bayes rate, prior stratified by county population". Somebody looking up one county will want its observed rate; somebody reading the map needs the stable one.
Code examples
Example 1 — EB with an optional stratified prior and the weight each area kept
import numpy as np
import pandas as pd
def smooth_rates(frame, events, population, per=1000, strata=None):
"""Raw and Empirical Bayes rates, plus the weight each area kept."""
out = frame.copy()
e = out[events].to_numpy(float)
b = out[population].to_numpy(float)
out["raw"] = e / b * per
out["eb"] = np.nan
out["kept"] = np.nan
groups = out.groupby(strata).groups if strata else {"all": out.index}
for _, index in groups.items():
i = out.index.get_indexer(index)
ei, bi = e[i], b[i]
mean = ei.sum() / bi.sum()
rate = ei / bi
variance = (bi * (rate - mean) ** 2).sum() / bi.sum() - mean / bi.mean()
variance = max(variance, 0.0)
kept = variance / (variance + mean / bi)
out.iloc[i, out.columns.get_loc("kept")] = kept
out.iloc[i, out.columns.get_loc("eb")] = (kept * rate + (1 - kept) * mean) * per
return outWithout strata it reproduces esda.smoothing.Empirical_Bayes exactly: the largest difference across the 3,144 counties was 7.1e-15. The extra kept column is what makes the result explainable. It also clips a negative variance estimate to zero, which esda's global class does not do.
Example 2 — a shrinkage report by size class
def shrinkage_report(frame, population,
bins=(0, 1000, 2500, 5000, 10000, 25000, 100000, np.inf)):
"""How far smoothing moved each size class."""
classes = pd.cut(frame[population], bins, right=False)
shift = (frame["raw"] - frame["eb"]).abs()
grouped = frame.groupby(classes, observed=True)
report = pd.DataFrame({
"areas": grouped.size(),
"raw_sd": grouped["raw"].std(),
"eb_sd": grouped["eb"].std(),
"median_kept": grouped["kept"].median(),
"median_shift_pct": (shift / frame["raw"].where(frame["raw"] > 0))
.groupby(classes, observed=True).median() * 100,
})
return report.round(3) areas raw_sd eb_sd median_kept median_shift_pct
[0.0, 1000.0) 36 11.098 2.197 0.344 25.713
[1000.0, 2500.0) 111 4.813 2.705 0.588 15.201
[2500.0, 5000.0) 178 3.652 2.687 0.739 9.171
[5000.0, 10000.0) 417 3.406 2.892 0.852 4.571
[10000.0, 25000.0) 782 2.902 2.672 0.926 2.248
[25000.0, 100000.0) 1000 2.668 2.581 0.971 0.701
[100000.0, inf) 620 2.341 2.325 0.994 0.079
Example 3 — scoring estimates against other years
def holdout_rmse(frame, predictions, reference, population,
bins=(0, 2500, 5000, 10000, 25000, 100000, np.inf)):
"""Score each estimate against a rate built from other years."""
classes = pd.cut(frame[population], bins, right=False)
rows = {}
for name in predictions:
error = (frame[name] - frame[reference]) ** 2
rows[name] = np.sqrt(error.groupby(classes, observed=True).mean())
rows[name]["all"] = np.sqrt(error.mean())
return pd.DataFrame(rows).round(2)
def national(frame, year):
return frame[f"DEATHS{year}"].sum() / frame[f"POPESTIMATE{year}"].sum()
years = (2021, 2022, 2024)
expected = sum(counties[f"POPESTIMATE{y}"] * national(counties, y) for y in years)
observed = counties[[f"DEATHS{y}" for y in years]].sum(axis=1)
counties["other_years"] = observed / expected * national(counties, 2023) * 1000Rescaling matters. Raw pooled rates for 2021–2024 run higher than 2023 because of the pandemic years, and that level shift alone makes every estimate look biased.
Explanation
Why small denominators produce extreme rates
Deaths behave roughly like a Poisson count, whose standard error is the square root of the count. For Loving County that is √3 deaths on 46 people: a standard error of 37.7 per 1,000 on a rate of 65.2. Across all counties, 306 have a 95% interval wider than a quarter of their own rate, and their median population is 2,582.
A choropleth ignores all of that. The county with the widest interval gets the darkest colour.
Why the weight depends on both population and variance
The shrinkage weight compares two variances. prior_mean / population is the noise expected from counting alone. variance is the real spread between counties once that noise is removed.
When real differences are large compared with the noise, areas keep their own rates. When noise dominates, as in a county of 46 people, they are pulled to the prior. Nothing about this uses geography until the spatial version replaces the national prior with a local one.
Why the global prior made these rates worse
EB assumes each area's true rate is drawn from one distribution with one mean. County crude death rates do not fit that assumption. Age structure pushes rural counties up systematically, so the counties that get shrunk hardest are exactly the ones whose true rate sits furthest above the prior mean.
The held-out check shows it directly. The bias in the smallest class was −2.03 per 1,000 under the global prior and +0.07 under a prior stratified by size. The variance went down in both cases; only the stratified version stopped trading that for bias.
Why the neighbours did not fix it
Spatial EB with six nearest neighbours cut the small-county bias from −2.03 to −1.17. It did not improve the RMSE (4.53 against 4.49). A small county's neighbours include larger towns, and the local prior is population-weighted, so the towns dominate it.
A prior based on something that actually drives the rate works better than proximity alone. Population class did that here; an age-standardised expected count would do it better still.
When smoothing earns its keep
EB helps most when differences between areas are mainly noise: rare events, similar underlying risk, and many areas with few cases. Crude county mortality has large real differences, so plain EB can easily over-smooth it. Run the held-out check. It takes seconds, and it is the only way to know which situation you are in.
Edge cases or notes
.ris a column vector. Use.ravel()before assigning it to a DataFrame column, or pandas raises on the shape.- Zero events are fine. Kalawao County, Hawaii (0 deaths, 81 people) becomes 8.8 per 1,000 instead of 0.
- esda's global class does not clip a negative variance estimate. With very homogeneous data the weights can fall outside 0–1; the function in Example 1 clips to zero.
- Islands have no neighbours. Queen contiguity left 8 US counties without neighbours; spatial EB then uses only the county's own values.
- Weights must declare their order.
use_index=Falsefails theid_ordercheck;use_index=Trueor a KNN weight passes it. - EB does not age-standardise. A smoothed crude rate still compares old counties with young ones.
- Do not smooth survey estimates as if they were counts. ACS estimates carry margins of error that the Poisson assumption does not describe.
- Keep the raw rate in the output. A smoothed value published without its method will be quoted as an observation.
Internal links
- Small numbers explained: why rates in small areas jump around — the problem this guide solves
- Fixing a rate map dominated by tiny, nearly empty areas — the map-level fixes, including pooling years
- How to normalise census counts into rates and densities — building the raw rates
- Counts, rates and densities: what a demographic map should show — choosing the measure first
- Margins of error explained — the survey equivalent of the Poisson noise here
- Fixing spatial weights that warn about islands — the neighbourless counties in step 5
- How to calculate Moran's I in Python — the same libpysal weights, used to test for pattern
- Choropleth classification explained — how the classes then hide or show the smoothing
FAQ
What does Empirical Bayes smoothing do to a rate?
It replaces each area's rate with a weighted average of that rate and a prior mean. Areas with small populations get most of their weight from the prior; large areas keep their own rate almost unchanged.
Which Python library implements it?
esda, in esda.smoothing: Empirical_Bayes(e, b) for a global prior and Spatial_Empirical_Bayes(e, b, w) for a neighbour-based one. Both return an object whose .r attribute is an n × 1 array of rates.
Why does Spatial Empirical Bayes raise an id_order error?
The weights object does not declare the order of its observations, so esda cannot be sure they line up with the counts. Build the weights with use_index=True or with a KNN constructor and the check passes.
Does smoothing always make rates more accurate?
No. On 2023 county death rates, global EB raised the error against other years from 0.92 to 1.30 deaths per 1,000. It pulled small rural counties towards a national mean that did not describe them.
How do I choose the prior?
Group areas by whatever drives the rate. A prior stratified by county population removed the bias here and beat the raw rate for counties of 2,500–10,000 people; age-standardised expected counts would be better where age is the driver.
Should I map the smoothed rate or the raw one?
Map the smoothed one if the held-out check says it is closer to other years. Publish the raw rate alongside it, and name the method in the legend.