Benchmarking Small Area Estimates in Python with svy-sae

Ratio and raking benchmarks to known or estimated targets

Benchmark Fay-Herriot small area estimates in Python with svy-sae so they agree with regional or national figures: ratio and raking adjustments, targets estimated from the same survey, and bootstrap MSE of the benchmarked estimates.

Keywords

benchmarking small area estimates, Fay-Herriot benchmarking Python, ratio benchmarking raking, consistent small area estimates, bootstrap MSE benchmarked estimates

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.

District estimates are rarely published alone. The survey that produced them also publishes regional and national figures, and readers expect the districts of a region to add up to the region. Model-based estimates do not do so on their own: Checking the model found regions where the districts’ aggregate and the regional estimate are far apart.

Benchmarking adjusts the estimates so each group’s size-weighted mean equals a target. It buys consistency, and its cost depends on the targets: adjusting to an accurate figure can improve the estimates, while adjusting to a noisy one passes its error on to every district of the group.

This tutorial demonstrates how to benchmark the district poverty rates to regional and national targets with Benchmark, choose between ratio and raking adjustments, account for targets that are themselves estimates, and measure what benchmarking does to accuracy.

NoteNot a search for the best model

This tutorial does not try to find the best model for these data. It shows the options svy-sae offers to evaluate a model and the ways to improve it; the choices made below are examples, not recommendations.

Setting up the inputs

We continue with the inputs and the model of the previous tutorials:

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")
)

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

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

Targets and sizes

The targets are the regional poverty rates, estimated directly from the survey:

regional = (
    hld_sample.estimation.prop(y="pov_status", by="geo1")
    .to_polars()
    .filter(pl.col("pov_status") == 1)
    .select("geo1", pl.col("est").alias("target"))
)

A region’s rate is the mean of its districts’ rates weighted by their number of households. The sizes frame gives every district its region and its size, n_hh from the census. It also holds each district’s weight in the regional direct estimate, the sum of its survey weights, used below:

weight_total = hld_data.group_by("geo2").agg(pl.col("hhweight").sum().alias("weight_total"))

sizes = (
    district_x.select("geo1", "geo2", "n_hh")
    .join(weight_total, on="geo2", how="left")
    .with_columns(pl.col("weight_total").fill_null(0))
)

The two districts outside the sample count toward their region’s size but have no survey weight, so their weight_total is 0.

Benchmarking to regional estimates

Pass a Benchmark to fh(). Benchmarked estimates need the bootstrap MSE: the analytical one is for the estimates before benchmarking.

fit_bench = model.fh(
    y="est",
    x=["piped_water"],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    benchmark=sae.Benchmark(
        method="ratio",
        targets=regional,
        sizes=sizes,
        area="geo2",
        size="n_hh",
        group="geo1",
        direct_weight="weight_total",
    ),
    mse="bootstrap",
    n_reps=1000,
    rstate=12345,
)

The predictions are the benchmarked estimates. The report shows each region’s target, the aggregate before benchmarking, and the relative change:

fit_bench.to_polars("benchmark")

shape: (10, 4)

group	target	before	change
str	f64	f64	f64
"geo_01"	0.085854	0.084013	0.021922
"geo_02"	0.341601	0.327366	0.043484
"geo_03"	0.218635	0.193967	0.127176
"geo_04"	0.171531	0.17307	-0.008893
"geo_05"	0.323229	0.307004	0.05285
"geo_06"	0.218746	0.26723	-0.181434
"geo_07"	0.227652	0.235576	-0.033638
"geo_08"	0.252795	0.348402	-0.274415
"geo_09"	0.320075	0.268993	0.189902
"geo_10"	0.227131	0.23195	-0.020776

fit_bench.plots.benchmark();

Change made to each region

Four regions move by more than 10%, which the fit flags as SAE_BENCHMARK_LARGE. The estimates before benchmarking stay in the diagnostics as pred_unbenchmarked, and the calibration check of Checking the model now reports no gap:

fit_bench.diagnostics.calibration(
    targets=regional, sizes=sizes, area="geo2", size="n_hh", group="geo1"
)

shape: (10, 4)

geo1	target	aggregate	difference
str	f64	f64	f64
"geo_01"	0.085854	0.085854	-3.3307e-16
"geo_02"	0.341601	0.341601	2.2204e-16
"geo_03"	0.218635	0.218635	2.2204e-16
"geo_04"	0.171531	0.171531	0.0
"geo_05"	0.323229	0.323229	0.0
"geo_06"	0.218746	0.218746	0.0
"geo_07"	0.227652	0.227652	0.0
"geo_08"	0.252795	0.252795	2.2204e-16
"geo_09"	0.320075	0.320075	-2.2204e-16
"geo_10"	0.227131	0.227131	0.0

Ratio or raking

Table 1: Options of Benchmark

Benchmark argument	Values	Description
method	"ratio"	Multiply a group’s estimates by target / aggregate
	"raking"	Add target − aggregate to a group’s estimates
targets	number	One target for all areas, with no group
	frame	One row per group, with group and a target column
sizes, area, size, group	frame and columns	Each area’s key, size, and group; sizes are normalised within each group
direct_weight	column of sizes	Each area’s weight in its group’s direct estimate, when the targets are direct estimates from the same survey

Ratio benchmarking keeps rates positive and changes each district in proportion to its level; raking shifts every district of a region by the same amount, which can push a low rate below zero.

Targets that are estimates

The regional targets come from the same survey as the district direct estimates, so they have sampling error too, correlated with the districts’. direct_weight tells the bootstrap how each regional estimate is built from the district ones, so it re-estimates the targets in every replicate. Without it, the targets are treated as known exactly, and the MSE leaves out their error. Compare the two against the true rates:

def scores(name, fit):
    scored = fit.to_polars("predictions").join(truth, left_on="area", right_on="geo2")
    error2 = (pl.col("pred") - pl.col("true_rate")) ** 2
    covered = (pl.col("true_rate") >= pl.col("lci")) & (pl.col("true_rate") <= pl.col("uci"))
    return scored.select(
        pl.lit(name).alias("estimates"),
        error2.mean().sqrt().alias("rmse"),
        pl.col("mse").mean().alias("mean mse"),
        error2.mean().alias("mean squared error"),
        covered.sum().alias("interval covers truth"),
    )


fit_fixed = model.fh(
    y="est",
    x=["piped_water"],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    benchmark=sae.Benchmark(
        method="ratio", targets=regional, sizes=sizes, area="geo2", size="n_hh", group="geo1"
    ),
    mse="bootstrap",
    n_reps=1000,
    rstate=12345,
)

pl.concat([
    scores("targets as known", fit_fixed),
    scores("targets as estimates", fit_bench),
])

shape: (2, 5)

estimates	rmse	mean mse	mean squared error	interval covers truth
str	f64	f64	f64	u32
"targets as known"	0.043754	0.001391	0.001914	41
"targets as estimates"	0.043754	0.002338	0.001914	60

The estimates are the same; only the MSE differs. Treated as known, the targets’ error is missing from the MSE, and the intervals miss the truth for a third of the districts. With direct_weight, the MSE covers it.

WarningEstimated targets

When the targets are direct estimates from the same survey, pass direct_weight. Leave it out only for targets known without sampling error, such as census counts or administrative totals.

What benchmarking does to accuracy

Benchmarking to the regional direct estimates moves the districts toward those estimates, errors included. In geo_08, the regional estimate is about one standard error below the true rate, and benchmarking pulls all its districts down with it.

The census gives the true regional and national rates. Here they stand in for figures known from another source, such as a census or an administrative register:

is_poor = pl.col("tot_exp") < pl.col("hhsize") * 1800
national_rate = hld_census.select(is_poor.mean()).item()
regional_rates = hld_census.group_by("geo1").agg(is_poor.mean().alias("target"))


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


pl.concat([
    scores("not benchmarked", fit),
    scores("regional direct estimates", fit_bench),
    scores("known national rate", benchmarked(sae.Benchmark(
        method="ratio", targets=national_rate, sizes=sizes, area="geo2", size="n_hh"
    ))),
    scores("known regional rates", benchmarked(sae.Benchmark(
        method="ratio", targets=regional_rates, sizes=sizes, area="geo2", size="n_hh",
        group="geo1",
    ))),
]).select("estimates", "rmse", "interval covers truth")

shape: (4, 3)

estimates	rmse	interval covers truth
str	f64	u32
"not benchmarked"	0.032911	58
"regional direct estimates"	0.043754	60
"known national rate"	0.022151	60
"known regional rates"	0.018418	56

Benchmarking to the regional direct estimates makes the district estimates less accurate: it trades accuracy for consistency with the published regional figures. Benchmarking to accurate figures improves them, because it corrects the model’s own errors at the regional level.

TipWhether to benchmark

Benchmark when consistency with published figures is required, or when the targets are more accurate than the model’s aggregates. Check calibration first: a large gap may come from the targets rather than the model.

Summary

This tutorial benchmarked the district poverty rates:

1. Benchmark with regional targets, district sizes, and groups, passed to fh() with the bootstrap MSE
2. The report, to_polars("benchmark") and plots.benchmark(), and calibration with no gap afterwards
3. Ratio or raking: scale or shift each group’s estimates
4. direct_weight for targets estimated from the same survey, so the MSE includes their error
5. Accuracy: benchmarking passes the targets’ errors to the districts; accurate targets improve the estimates

Next steps

The area-level workflow is complete, from direct estimates to benchmarked district estimates. With unit records and a census, unit-level models use more of the data.

Go to the unit level
Continue to The BHF Model →