The Fay–Herriot Model in Python with svy-sae

Area-level EBLUP with analytical and bootstrap MSE

Tutorials
Area-level
Fay-Herriot
Python
Fit the Fay-Herriot area-level model in Python with svy-sae: EBLUP estimates, shrinkage, out-of-sample areas, and Prasad-Rao or bootstrap MSE, checked against known true values.
Author

Mamadou S. Diallo, Ph.D.

Modified

November 5, 2026

Keywords

Fay-Herriot model Python, area-level small area estimation, EBLUP area level Python, Prasad-Rao MSE, parametric bootstrap MSE, out-of-sample small areas

WarningAbout these examples

This tutorial illustrates the svy-sae API; it is not a rigorous small area analysis. The models are kept simple so each step stays readable, and in many places, from the choice of covariates to the variance model, better choices could be made to improve them.

The Fay–Herriot model (Fay and Herriot 1979) combines two sources of information about each area. The direct estimate \(\hat\theta_i\) is unbiased but noisy, with known sampling variance \(D_i\); a regression on covariates \(x_i\), observed for every area, is stable but ignores what is specific to the area:

\[ \hat\theta_i = \theta_i + e_i, \qquad \theta_i = x_i'\beta + u_i, \qquad e_i \sim N(0, D_i), \quad u_i \sim N(0, \sigma^2_u). \]

The empirical best linear unbiased predictor (EBLUP) is a weighted average of the two:

\[ \tilde\theta_i = \gamma_i \hat\theta_i + (1 - \gamma_i)\, x_i'\hat\beta, \qquad \gamma_i = \frac{\hat\sigma^2_u}{\hat\sigma^2_u + D_i}. \]

An area with a precise direct estimate keeps most of it; an area with a noisy one leans on the regression. Areas outside the sample get the regression prediction alone, the synthetic estimate \(x_i'\hat\beta\).

This tutorial demonstrates how to fit the model with AreaLevel.fh(), read its results, predict districts the survey did not reach, and measure uncertainty with the analytical (Prasad–Rao) and bootstrap MSE. Because the data are synthetic, each step is checked against the true district poverty rates.

Setting up the inputs

We continue with the district poverty rates and census covariates from Preparing inputs:

Rebuild the inputs from Preparing inputs 
import polars as pl
import svy

import svy_sae as sae

hld_data = svy.datasets.load("hld_sample_wb_2023").data
hld_design = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
hld_sample = svy.Sample(data=hld_data, design=hld_design).wrangling.mutate(
    {
        "hhpovline": svy.col("hhsize") * 1800,
        "pov_status": svy.when(svy.col("tot_exp") < svy.col("hhpovline")).then(1).otherwise(0),
    }
)

direct = (
    hld_sample.estimation.prop(y="pov_status", by="geo2")
    .to_polars()
    .filter(pl.col("pov_status") == 1)
    .select("geo2", "est", "se", "cv", "lci", "uci", "df", "n")
    .with_columns((pl.col("se") ** 2).alias("var_direct"))
)

hld_census = svy.datasets.load("hld_pop_wb_2023").data
yes_no = ["electricity", "car", "refrigerator", "tv", "cell", "bank"]
district_x = (
    hld_census.group_by("geo1", "geo2")
    .agg(
        (pl.col("urbrur") == "Urban").mean().alias("urban"),
        (pl.col("piped_water") != "No piped water").mean().alias("piped_water"),
        (pl.col("flush_toilet") == "Flush toilet").mean().alias("flush_toilet"),
        *[(pl.col(c) == "Yes").mean().alias(c) for c in yes_no],
        pl.col("hhsize").mean().alias("hhsize"),
        pl.col("rooms").mean().alias("rooms"),
        pl.len().alias("n_hh"),
    )
    .sort("geo2")
)

truth = (
    hld_census.group_by("geo2")
    .agg((pl.col("tot_exp") < pl.col("hhsize") * 1800).mean().alias("true_rate"))
    .sort("geo2")
)

Smooth the sampling variances as in Smoothing variances:

direct = direct.with_columns(pl.col("n").log().alias("log_n"))
smoothed = sae.smooth_variance(direct, x="log_n", variance="var_direct", df="df")

A baseline without covariates

AreaLevel takes the frame of direct estimates, one row per area. Start with a model that has no covariates, x=[]: every district then borrows from the national average.

model = sae.AreaLevel(smoothed)

fit_base = model.fh(
    y="est",
    x=[],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
)

print(fit_base)
╭─────────────── SAE: Fay-Herriot (MSE: Prasad-Rao) - Area Level ───────────────╮
│ Method     REML  Log Likelihood   24.1529                                     │
│ No. Areas    59  AIC             -44.3057                                     │
│ No. Obs      59  BIC             -40.1506                                     │
│                  KIC             -42.3057                                     │
│                                                                               │
│                                                                               │
│  Term        Estimate   Std.Err       z   p-value   [0.025   0.975]           │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━           │
│  Intercept     0.2309    0.0205   11.25    0.0000   0.1906   0.2711           │
│                                                                               │
│                                                                               │
│ Variance components: area variance: 0.0124191 (se 0.00464)                    │
│                                                                               │
│ Area keys (geo2): data 59 rows, pop_data 61 rows, in both 59, out-of-sample 2 │
│                                                                               │
│ Findings:                                                                     │
│ - Residuals not normal                                                        │
│ - Area effects not normal                                                     │
│ - 2 influential areas                                                         │
│ - 41 of 61 areas with CV above 0.3                                            │
│ - 6 sampled areas without a precision gain                                    │
╰───────────────────────────────────────────────────────────────────────────────╯

pop_data is the census frame from Preparing inputs. It lists every district to predict, and supplies the covariates that the frame of direct estimates does not have. The summary reports how the districts matched: 59 in both frames, 2 in the census only.

The findings at the bottom are the model checks covered in Checking the model.

How much the model borrows

The per-district diagnostics include the shrinkage factor \(\gamma_i\), the weight on the direct estimate. Compare the districts with the smallest and largest samples:

borrowing = (
    fit_base.to_polars("diagnostics")
    .join(direct.select(pl.col("geo2").alias("area"), "n"), on="area")
    .select("area", "n", "direct", "direct_var", "synthetic", "gamma", "pred")
    .sort("n")
)

pl.concat([borrowing.head(3), borrowing.tail(3)])
shape: (6, 7)
area n direct direct_var synthetic gamma pred
str i64 f64 f64 f64 f64 f64
"geo_02_01" 25 0.12 0.044846 0.230852 0.216871 0.206811
"geo_09_01" 25 0.0 0.044846 0.230852 0.216871 0.180787
"geo_10_07" 25 0.44 0.044846 0.230852 0.216871 0.27621
"geo_03_01" 300 0.279496 0.005448 0.230852 0.69507 0.264663
"geo_10_08" 300 0.150211 0.005448 0.230852 0.69507 0.174801
"geo_04_02" 475 0.148295 0.003689 0.230852 0.770962 0.167204

Districts of 25 households keep less than a quarter of their direct estimate; the largest keep more than half. Each estimate pred lies between direct and synthetic, at the point set by gamma.

Adding a covariate

The share of households with piped water, from the census, tracks poverty closely. It is the covariate chosen in Choosing covariates:

fit = model.fh(
    y="est",
    x=["piped_water"],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
)

print(fit)
╭─────────────── SAE: Fay-Herriot (MSE: Prasad-Rao) - Area Level ───────────────╮
│ Method     REML  Log Likelihood   41.7020                                     │
│ No. Areas    59  AIC             -77.4040                                     │
│ No. Obs      59  BIC             -71.1714                                     │
│                  KIC             -74.4040                                     │
│                                                                               │
│                                                                               │
│  Term          Estimate   Std.Err       z   p-value    [0.025    0.975]       │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━       │
│  Intercept       0.6029    0.0527   11.43    0.0000    0.4995    0.7062       │
│  piped_water    -0.6602    0.0894   -7.38    0.0000   -0.8354   -0.4850       │
│                                                                               │
│                                                                               │
│ Variance components: area variance: 5.86419e-07 (se 4.607e-05)                │
│                                                                               │
│ Area keys (geo2): data 59 rows, pop_data 61 rows, in both 59, out-of-sample 2 │
│ Covariates from data: -; from pop_data: piped_water                           │
│                                                                               │
│ Findings:                                                                     │
│ - Area variance at the boundary                                               │
│ - 2 influential areas                                                         │
│ - 12 of 61 areas with CV above 0.3                                            │
╰───────────────────────────────────────────────────────────────────────────────╯

piped_water is not in smoothed, so it is read from pop_data for every district, as the summary reports.

An area variance at the boundary

The area variance has collapsed:

print(fit_base.component("area"))
print(fit.component("area"))
VarianceComponent(name='area', kind=variance, est=0.0124191, se=0.00463974)
VarianceComponent(name='area', kind=variance, est=5.86419e-07, se=4.60727e-05)

With \(\hat\sigma^2_u\) near zero, every \(\gamma_i\) is near zero and every estimate is essentially synthetic. The fit says so in its findings (SAE_SIGMA_U_BOUNDARY). It means the regression on piped water leaves no between-district variation that the model can tell apart from sampling error.

Here the truth can confirm it. Regress the true district rates on piped water and measure the variance left over:

import numpy as np

sampled = district_x.join(truth, on="geo2").join(direct.select("geo2"), on="geo2")
X = np.column_stack([np.ones(sampled.height), sampled["piped_water"].to_numpy()])
resid = sampled["true_rate"].to_numpy() - X @ np.linalg.lstsq(X, sampled["true_rate"].to_numpy())[0]

print(f"True area variance around the regression: {resid.var(ddof=2):.5f}")
print(f"Median smoothed sampling variance:         {smoothed['variance_smoothed'].median():.5f}")
True area variance around the regression: 0.00049
Median smoothed sampling variance:         0.01384

The true rates sit close to the line: what remains is small next to the sampling variances, too small for 59 districts to estimate. The synthetic estimate is the right answer for these data.

NoteMethod note

A boundary estimate is not always this benign. With weak covariates it can also mean the sampling variances are too large, which is one reason to smooth them and check them against the direct variances. The MSE also ignores the uncertainty of a boundary estimate.

Reading the results

Fixed effects, with standard errors and confidence intervals:

fit.to_polars("fixed")
shape: (2, 7)
term est se t_score p_value lci uci
str f64 f64 f64 f64 f64 f64
"Intercept" 0.60287 0.05273 11.433154 0.0 0.499519 0.70622
"piped_water" -0.660175 0.089397 -7.38478 1.5277e-13 -0.835393 -0.484957

Variance components, one row each (the Fay–Herriot model has one, area):

fit.to_polars("components")
shape: (1, 4)
name kind est se
str str f64 f64
"area" "variance" 5.8642e-7 0.000046

Model statistics. The information criteria are the ones used in Choosing covariates:

s = fit.stats
print(f"log-likelihood {s.log_likelihood:.2f}, AIC {s.aic:.2f}, BIC {s.bic:.2f}, KIC {s.kic:.2f}")
log-likelihood 41.70, AIC -77.40, BIC -71.17, KIC -74.40

Predictions

The estimates, with their MSE, coefficient of variation, and 95% interval:

pred = fit.to_polars("predictions")

pred.select("area", "pred", "mse", "cv", "lci", "uci", "in_sample").head()
shape: (5, 7)
area pred mse cv lci uci in_sample
str f64 f64 f64 f64 f64 bool
"geo_01_01" 0.088987 0.001532 0.439896 0.012263 0.165711 true
"geo_01_02" 0.075821 0.000858 0.38636 0.018404 0.133238 true
"geo_01_03" 0.086521 0.001712 0.478209 0.005426 0.167616 true
"geo_01_04" 0.093721 0.00087 0.314671 0.035918 0.151524 true
"geo_01_05" 0.0776 0.001053 0.418075 0.014012 0.141187 true

Districts outside the sample

in_sample is False for the districts that only pop_data lists:

pred.filter(~pl.col("in_sample")).select("area", "pred", "mse", "cv")
shape: (2, 4)
area pred mse cv
str f64 f64 f64
"geo_02_02" 0.385972 0.000651 0.066115
"geo_09_03" 0.313023 0.000323 0.057378

They get the synthetic estimate \(x_i'\hat\beta\), with an MSE of \(\hat\sigma^2_u\) plus the uncertainty of \(\hat\beta\). fit.area_match records how the keys of the two frames matched and where each covariate came from:

print(fit.area_match)
Area keys (geo2): data 59 rows, pop_data 61 rows, in both 59, out-of-sample 2
Covariates from data: -; from pop_data: piped_water

Measuring uncertainty

mse selects how the MSE is computed, and method how \(\sigma^2_u\) is estimated:

Table 1: MSE and estimation options for fh()
Argument Value Description
mse "analytical" (default) Prasad–Rao second-order approximation (Prasad and Rao 1990), with a bias term under ML
"bootstrap" Parametric bootstrap with n_reps replicates (default 1,000); rstate seeds it
method "reml" (default) Restricted maximum likelihood
"ml" Maximum likelihood

Refit with the bootstrap MSE. The estimates are unchanged; only the MSE differs:

fit_boot = model.fh(
    y="est",
    x=["piped_water"],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    mse="bootstrap",
    n_reps=1000,
    rstate=12345,
)

mse_compare = (
    pred.select("area", "pred", pl.col("mse").alias("analytical"))
    .join(fit_boot.to_polars("predictions").select("area", pl.col("mse").alias("bootstrap")), on="area")
    .join(truth, left_on="area", right_on="geo2")
    .with_columns(((pl.col("pred") - pl.col("true_rate")) ** 2).alias("squared error"))
)

mse_compare.select(pl.col("analytical", "bootstrap", "squared error").mean())
shape: (1, 3)
analytical bootstrap squared error
f64 f64 f64
0.000969 0.000385 0.001083

Here the two disagree: the bootstrap MSE is less than half the analytical one. The parametric bootstrap draws new data from the fitted model, and with \(\hat\sigma^2_u\) at zero those data have no area effects, so its MSE only reflects the uncertainty of \(\hat\beta\). The Prasad–Rao MSE also accounts for the uncertainty of \(\hat\sigma^2_u\), and its average matches the actual squared error against the truth.

WarningThe bootstrap at the boundary

When the fit reports SAE_SIGMA_U_BOUNDARY, prefer the analytical MSE where it exists. The bootstrap is still required where no analytical MSE exists: for most transformations, and for benchmarked estimates.

How close did we get?

Compare the direct estimates and both models with the true rates: error (RMSE), districts with a CV above 30%, and how many 95% intervals contain the truth.

def scores(name, frame):
    scored = frame.join(truth, left_on="area", right_on="geo2")
    covered = (pl.col("true_rate") >= pl.col("lci")) & (pl.col("true_rate") <= pl.col("uci"))
    return scored.select(
        pl.lit(name).alias("estimates"),
        pl.len().alias("districts"),
        ((pl.col("pred") - pl.col("true_rate")) ** 2).mean().sqrt().alias("rmse"),
        (pl.col("cv") > 0.3).sum().alias("cv above 30%"),
        covered.sum().alias("interval covers truth"),
    )


direct_pred = direct.select(
    pl.col("geo2").alias("area"), pl.col("est").alias("pred"), "cv", "lci", "uci"
)

pl.concat([
    scores("direct", direct_pred),
    scores("FH, no covariates", fit_base.to_polars("predictions")),
    scores("FH, piped water", pred),
])
shape: (3, 5)
estimates districts rmse cv above 30% interval covers truth
str u32 f64 u32 u32
"direct" 59 0.155022 46 52
"FH, no covariates" 61 0.106323 41 56
"FH, piped water" 61 0.032911 12 58

The baseline already cuts the error by a third by borrowing from the national average. Piped water does far better, for all 61 districts including the two the survey did not reach, and most of its intervals contain the true rate.

Summary

This tutorial fitted the Fay–Herriot model with AreaLevel.fh():

  1. A baseline with x=[], shrinking each district toward the national average, more for smaller samples
  2. A covariate read from the census through pop_data, which also predicts the districts outside the sample
  3. An area variance at the boundary, flagged by the fit and explained by the truth
  4. Results: fixed effects, variance components, statistics, and predictions with MSE, CV, and intervals
  5. Uncertainty from the Prasad–Rao approximation or the parametric bootstrap, which understates the MSE when \(\hat\sigma^2_u\) is at the boundary

Next steps

Piped water was given here. The next tutorial chooses it among the census covariates.

Choose the covariates
Continue to Choosing Covariates →

References

Fay, Robert E., and Roger A. Herriot. 1979. “Estimates of Income for Small Places: An Application of James-Stein Procedures to Census Data.” Journal of the American Statistical Association 74 (366): 269–77.
Prasad, N. G. N., and J. N. K. Rao. 1990. “The Estimation of the Mean Squared Error of Small-Area Estimators.” Journal of the American Statistical Association 85 (409): 163–71.