import polars as pl
import svy
import svy_sae as sae
hld_data = svy.datasets.load("hld_sample_wb_2023")
hld_design = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
hld_sample = svy.Sample(data=hld_data, design=hld_design)Direct Estimates for Small Areas
Tutorials
Area-level
Direct Estimation
Python
Compute design-based direct estimates and sampling variances by area with svy, then feed them into a Fay-Herriot small area model in Python with svy-sae.
Keywords
direct estimation Python, survey direct estimates by domain, sampling variance by area Python, Fay-Herriot inputs, design-based estimation small areas, small area estimation workflow Python
Why direct estimates come first
Every area-level small area model starts from direct estimates: design-based estimates computed separately for each area, with their sampling variances. The Fay–Herriot model treats those variances as known inputs, so the quality of an SAE workflow rests on computing them properly — with the survey weights, stratification, and clustering of the actual design.
svy-sae deliberately does not compute direct estimates. That is survey estimation, and it belongs to svy, which handles design specification, weighting, and variance estimation. The division of labor is:
- svy — from microdata and a survey design to direct estimates and variances per area,
- svy-sae — from those per-area inputs to model-based small area estimates.
This page walks the first step end to end and hands the result to AreaLevel.fh(). For everything svy can do — totals, proportions, ratios, quantiles, replicate-weight variance — see the svy estimation tutorial.
Set up the sample
We use the synthetic household survey from World Bank (2023): 8,000 households with a stratified two-stage design (strata from region and urban/rural, EAs as PSUs, household weights).
Direct estimates by area
The small areas are the geo2 districts. One call produces the design-based mean of per-capita expenditure per district, with standard errors that respect the stratification and clustering:
direct = hld_sample.estimation.mean(y="pc_exp", by="geo2")
direct_df = direct.to_polars()
direct_df.head()
shape: (5, 7)
| geo2 | est | se | lci | uci | cv | df |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | i64 |
| "geo_01_01" | 5552.510026 | 510.352944 | 4398.011457 | 6707.008594 | 0.091914 | 9 |
| "geo_01_02" | 3243.737136 | 582.766639 | -4161.015097 | 10648.489368 | 0.179659 | 1 |
| "geo_01_03" | 4800.331744 | 358.215217 | 4011.905369 | 5588.75812 | 0.074623 | 11 |
| "geo_01_04" | 2820.689312 | 161.300991 | 2126.667162 | 3514.711462 | 0.057185 | 2 |
| "geo_01_05" | 5084.202403 | 1400.069362 | 628.556834 | 9539.847971 | 0.275376 | 3 |
The Fay–Herriot model needs the sampling variance, which is the squared standard error:
direct_df = direct_df.with_columns((pl.col("se") ** 2).alias("var_direct"))
direct_df.select(["geo2", "est", "se", "var_direct"]).head()
shape: (5, 4)
| geo2 | est | se | var_direct |
|---|---|---|---|
| str | f64 | f64 | f64 |
| "geo_01_01" | 5552.510026 | 510.352944 | 260460.127666 |
| "geo_01_02" | 3243.737136 | 582.766639 | 339616.955918 |
| "geo_01_03" | 4800.331744 | 358.215217 | 128318.141339 |
| "geo_01_04" | 2820.689312 | 161.300991 | 26018.009749 |
| "geo_01_05" | 5084.202403 | 1400.069362 | 1.9602e6 |
NoteSmoothing the variances
Fay–Herriot treats the sampling variances as known, but direct variance estimates from small samples are themselves noisy. In production workflows, practitioners often smooth them — for example with a Generalized Variance Function — before fitting. The intermediate frame above is in your hands for exactly that reason.
Add area-level covariates
The model borrows strength through covariates observed for every area. In real applications these come from a census or administrative registers; here we aggregate two household characteristics from the survey to keep the example self-contained:
covariates = hld_data.group_by("geo2").agg(
pl.col("hhsize").mean().alias("avg_hhsize"),
pl.col("rooms").mean().alias("avg_rooms"),
)
fh_input = direct_df.join(covariates, on="geo2", how="inner")
fh_input.head()
shape: (5, 10)
| geo2 | est | se | lci | uci | cv | df | var_direct | avg_hhsize | avg_rooms |
|---|---|---|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | i64 | f64 | f64 | f64 |
| "geo_01_05" | 5084.202403 | 1400.069362 | 628.556834 | 9539.847971 | 0.275376 | 3 | 1.9602e6 | 3.26 | 2.16 |
| "geo_03_04" | 4546.23416 | 630.008049 | 2541.267371 | 6551.200949 | 0.138578 | 3 | 396910.142115 | 3.83 | 3.14 |
| "geo_10_03" | 4721.753013 | 517.455436 | -1853.141693 | 11296.64772 | 0.10959 | 1 | 267760.127854 | 2.266667 | 2.506667 |
| "geo_02_04" | 2140.028199 | 185.839408 | -221.285364 | 4501.341762 | 0.08684 | 1 | 34536.285482 | 4.3 | 2.3 |
| "geo_06_03" | 4232.510167 | 576.128087 | 2751.525772 | 5713.494562 | 0.13612 | 5 | 331923.572379 | 4.16 | 3.508571 |
Fit the Fay–Herriot model
The frame now has one row per area with the three inputs the model needs — a direct estimate, its variance, and covariates:
model = sae.AreaLevel(fh_input)
res = model.fh(
y="est",
x=["avg_hhsize", "avg_rooms"],
variance="var_direct",
area="geo2",
method="reml",
)
pred = res.to_polars("predictions")
pred.head()
shape: (5, 7)
| area | pred | mse | cv | lci | uci | indicator |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | str |
| "geo_01_01" | 4801.573276 | 92125.45622 | 0.063213 | 4206.670635 | 5396.475917 | "mean" |
| "geo_01_02" | 3398.233821 | 100293.551357 | 0.093193 | 2777.518341 | 4018.949301 | "mean" |
| "geo_01_03" | 4585.581801 | 68315.531286 | 0.056999 | 4073.291898 | 5097.871704 | "mean" |
| "geo_01_04" | 2742.648625 | 23067.309035 | 0.055377 | 2444.965279 | 3040.33197 | "mean" |
| "geo_01_05" | 3161.856973 | 118377.888888 | 0.108816 | 2487.497648 | 3836.216297 | "mean" |
The FH estimates shrink each direct estimate toward the regression fit, with more shrinkage where the direct variance is large — precisely the areas where the survey alone is least reliable.
Compare direct and model-based estimates
comparison = (
fh_input.select(["geo2", pl.col("est").alias("direct"), "var_direct"])
.join(
pred.select([pl.col("area").alias("geo2"), pl.col("pred").alias("fh"), "mse"]),
on="geo2",
)
.with_columns((pl.col("mse") / pl.col("var_direct")).alias("mse_ratio"))
.sort("mse_ratio")
)
comparison.head(10)
shape: (10, 6)
| geo2 | direct | var_direct | fh | mse | mse_ratio |
|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 |
| "geo_09_07" | 4345.475799 | 2.1311e6 | 3499.806613 | 114913.141498 | 0.053922 |
| "geo_01_05" | 5084.202403 | 1.9602e6 | 3161.856973 | 118377.888888 | 0.060391 |
| "geo_10_04" | 5297.408579 | 2.1117e6 | 4375.762105 | 132798.440378 | 0.062887 |
| "geo_06_04" | 5584.606996 | 1.4511e6 | 4052.981934 | 112117.562029 | 0.077262 |
| "geo_02_05" | 3414.748731 | 1.4190e6 | 3040.412408 | 115272.328827 | 0.081232 |
| "geo_09_04" | 3703.240068 | 1.3415e6 | 2816.436193 | 121767.490898 | 0.090769 |
| "geo_10_09" | 5315.682314 | 1.2911e6 | 4341.779676 | 118443.2526 | 0.091737 |
| "geo_09_02" | 3281.342534 | 1.2672e6 | 2109.466515 | 129852.395442 | 0.102469 |
| "geo_10_05" | 4619.225316 | 1.0897e6 | 4515.007804 | 124434.788347 | 0.114195 |
| "geo_07_04" | 4482.657607 | 924187.53301 | 4087.091603 | 109102.175851 | 0.118052 |
An mse_ratio well below 1 means the model-based estimate is more precise than the direct one — the gain that small area estimation exists to deliver, largest in the areas with the smallest samples.
What’s next
- Fay–Herriot in depth — analytical and bootstrap MSE, categorical covariates, and diagnostics.
- svy documentation — the full direct estimation toolkit: totals, proportions, ratios, quantiles, and replicate-weight variance methods.
References
World Bank. 2023. “Synthetic Data for an Imaginary Country, Sample, 2023.” World Bank, Development Data Group. https://doi.org/10.48529/MC1F-QH23.