Preparing Area-Level Inputs in Python with svy-sae

Direct estimates from the survey, covariates from the census

Prepare the inputs of an area-level small area model in Python: design-based direct estimates and sampling variances by district with svy, and district covariates from a census.

Keywords

direct estimation Python, survey direct estimates by domain, sampling variance by area Python, Fay-Herriot inputs, census covariates small area estimation, poverty mapping Python

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.

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, so the quality of the whole workflow rests on computing them with the weights, strata, and clusters of the actual design. The model then borrows strength through covariates observed for every area, typically from a census or administrative registers.

svy-sae does not compute direct estimates. That is survey estimation, and it belongs to svy, which carries the design. The division of labor is:

1. svy: from microdata and a survey design to direct estimates and variances per area,
2. svy-sae: from those per-area inputs to model-based small area estimates.

This tutorial prepares the two inputs used by every area-level tutorial that follows: district poverty rates from a household survey, and district covariates from a census. Later tutorials smooth the variances, choose covariates, fit and check the model, and benchmark the estimates, all on these same districts.

Setting up the sample

We use the synthetic household survey of an imaginary country from World Bank (2023b), the same data as the svy tutorials: 8,000 households, with strata from region and urban/rural, enumeration areas (EAs) as PSUs, and household weights.

import polars as pl
import svy

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)

As in svy’s estimation tutorial, a household is poor when its total expenditure falls below 1,800 per person:

hld_sample = hld_sample.wrangling.mutate(
    {
        "hhpovline": svy.col("hhsize") * 1800,
        "pov_status": svy.when(svy.col("tot_exp") < svy.col("hhpovline")).then(1).otherwise(0),
    }
)

NoteMethod note

Scaling the poverty line by household size keeps the example simple. Applied work typically uses adult-equivalence scales and regional price indices.

Direct estimates by district

The small areas are the districts (geo2), nested in ten regions (geo1). Estimate the poverty rate of each district:

pov_geo2 = hld_sample.estimation.prop(y="pov_status", by="geo2")

pov_geo2.to_polars().head(4)

shape: (4, 13)

y_var	geo2	pov_status	est	se	lci	uci	cv	df	n	param	variance	alpha
str	str	i64	f64	f64	f64	f64	f64	f64	i64	str	str	f64
"pov_status"	"geo_01_01"	0	0.928575	0.043272	0.748208	0.982723	0.046601	9.0	250	"Proportion"	"Taylor"	0.05
"pov_status"	"geo_01_01"	1	0.071425	0.043272	0.017277	0.251792	0.60584	9.0	250	"Proportion"	"Taylor"	0.05
"pov_status"	"geo_01_02"	0	0.874377	0.056925	0.009521	0.999802	0.065104	1.0	50	"Proportion"	"Taylor"	0.05
"pov_status"	"geo_01_02"	1	0.125623	0.056925	0.000198	0.990479	0.453142	1.0	50	"Proportion"	"Taylor"	0.05

prop() returns a row per category. Keep the share of poor households, and add the sampling variance, the squared standard error, which the Fay–Herriot model takes as input:

direct = (
    pov_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"))
)

print(f"Districts in the sample: {direct.height}")
direct.head()

Districts in the sample: 59

shape: (5, 9)

geo2	est	se	cv	lci	uci	df	n	var_direct
str	f64	f64	f64	f64	f64	f64	i64	f64
"geo_01_01"	0.071425	0.043272	0.60584	0.017277	0.251792	9.0	250	0.001872
"geo_01_02"	0.125623	0.056925	0.453142	0.000198	0.990479	1.0	50	0.00324
"geo_01_03"	0.048261	0.014688	0.304357	0.02447	0.092976	11.0	300	0.000216
"geo_01_04"	0.232343	0.075106	0.323254	0.047113	0.649463	2.0	75	0.005641
"geo_01_05"	0.092234	0.083132	0.901315	0.004293	0.705405	3.0	100	0.006911

Keep n (households sampled) and df (design degrees of freedom): the variance smoothing tutorial uses both.

How precise are they?

A common publication rule withholds estimates with a coefficient of variation above 30%:

n_unreliable = direct.filter(pl.col("cv") > 0.3).height

print(f"Districts with a CV above 30%: {n_unreliable} of {direct.height}")

Districts with a CV above 30%: 46 of 59

Most districts fail it. With 25 to 475 households per district, the survey supports regional estimates, not district ones. Raising that precision is what small area estimation is for.

Districts that need care

A few districts have a sampling variance of zero:

direct.filter(pl.col("var_direct") < 1e-12)

shape: (4, 9)

geo2	est	se	cv	lci	uci	df	n	var_direct
str	f64	f64	f64	f64	f64	f64	i64	f64
"geo_02_01"	0.12	3.4694e-18	2.8912e-17	NaN	NaN	0.0	25	1.2037e-35
"geo_09_01"	0.0	0.0	inf	NaN	NaN	0.0	25	0.0
"geo_10_01"	0.0	0.0	inf	NaN	NaN	1.0	50	0.0
"geo_10_07"	0.44	9.3675e-17	2.1290e-16	NaN	NaN	0.0	25	8.7750e-33

They are zero for two different reasons:

• geo_02_01 and geo_10_07 were sampled in a single EA. With one PSU, the design gives no information on the variance (df is 0), and the computed variance is zero up to rounding.
• geo_09_01 and geo_10_01 have no poor household in the sample. The estimate is 0, and so is its estimated variance.

None of these variances is a real zero. Fay–Herriot would trust these direct estimates completely. The next tutorial replaces them with variances predicted from the other districts.

District covariates from the census

The model needs covariates for every district, measured the same way everywhere. Here they come from the census of the same country (World Bank 2023a), 2.5 million households:

hld_census = svy.datasets.load("hld_pop_wb_2023").data

print(f"Census households: {hld_census.height:,}")

Census households: 2,501,755

Summarize each district as shares of households and averages:

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

district_x.head()

shape: (5, 14)

geo1	geo2	urban	piped_water	flush_toilet	electricity	car	refrigerator	tv	cell	bank	hhsize	rooms	n_hh
str	str	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	u32
"geo_01"	"geo_01_01"	1.0	0.778402	0.611179	0.940061	0.288855	0.710622	0.79275	0.691367	0.718783	3.609577	2.750666	59310
"geo_01"	"geo_01_02"	1.0	0.798348	0.667126	0.969459	0.343759	0.791465	0.855172	0.724749	0.781221	3.828793	3.091558	31237
"geo_01"	"geo_01_03"	1.0	0.782133	0.618622	0.950656	0.295008	0.732066	0.81118	0.694966	0.731907	3.709907	2.7907	81591
"geo_01"	"geo_01_04"	1.0	0.77124	0.605938	0.949469	0.299552	0.726602	0.803569	0.693654	0.727941	3.86503	2.898997	50048
"geo_01"	"geo_01_05"	1.0	0.795654	0.643752	0.952594	0.308929	0.741894	0.821767	0.717052	0.757229	3.768314	2.88201	22824

n_hh, the number of households in each district, is the size used to aggregate district estimates to regions in the benchmarking tutorial.

Districts outside the sample

The census covers districts the survey did not reach:

district_x.join(direct, on="geo2", how="anti").select("geo1", "geo2", "n_hh")

shape: (2, 3)

geo1	geo2	n_hh
str	str	u32
"geo_02"	"geo_02_02"	10096
"geo_09"	"geo_09_03"	13912

These two districts have no direct estimate, but they have covariates, so the model can still predict them. Keep the census frame separate rather than joining it to direct: AreaLevel.fh() takes it as pop_data and matches the districts itself.

The true poverty rates

Because the data are synthetic, the census also records every household’s expenditure, so the true poverty rate of each district is known:

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

Measure how far the direct estimates are from the truth:

check = direct.join(truth, on="geo2").with_columns(
    (pl.col("est") - pl.col("true_rate")).alias("error")
)

rmse_direct = check.select((pl.col("error") ** 2).mean().sqrt()).item()
print(f"Direct estimates, RMSE against the truth: {rmse_direct:.3f}")

check.select("geo2", "n", "est", "true_rate", "error").sort(
    pl.col("error").abs(), descending=True
).head()

Direct estimates, RMSE against the truth: 0.155

shape: (5, 5)

geo2	n	est	true_rate	error
str	i64	f64	f64	f64
"geo_02_01"	25	0.12	0.572237	-0.452237
"geo_10_01"	50	0.0	0.423086	-0.423086
"geo_09_01"	25	0.0	0.412918	-0.412918
"geo_05_01"	50	0.055238	0.432254	-0.377016
"geo_05_03"	100	0.69343	0.367367	0.326063

The largest errors are in the districts with the fewest sampled households. Later tutorials report the same RMSE for the model-based estimates.

WarningOnly with synthetic data

In practice the truth is unknown, which is why the model’s own measures of uncertainty (MSE and CV) matter. The comparison here shows what those measures cannot: how close the estimates actually are.

Summary

This tutorial prepared the inputs of an area-level model:

1. Direct estimates with svy’s prop(by="geo2"): one poverty rate per sampled district, with its sampling variance, n, and df
2. Covariates from the census: shares and averages for every district, including the two outside the sample
3. True rates from the census, to check each estimate against

Next steps

The variances of four districts cannot be used as they are. Smooth them before fitting the model.

Fix the variances
Continue to Smoothing Sampling Variances →

References

World Bank. 2023a. “Synthetic Data for an Imaginary Country, Full Population, 2023.” World Bank, Development Data Group. https://doi.org/10.48529/78M1-AE09.

———. 2023b. “Synthetic Data for an Imaginary Country, Sample, 2023.” World Bank, Development Data Group. https://doi.org/10.48529/MC1F-QH23.