Preparing Unit-Level Inputs in Python with svy-sae

Household records from the survey, the same covariates from the census

Tutorials
Unit-level
Python
Prepare the inputs of a unit-level small area model in Python: household records and covariates from a survey, the same covariates for every census household, and district means for EBLUP with svy-sae.
Author

Mamadou S. Diallo, Ph.D.

Modified

November 5, 2026

Keywords

unit-level small area estimation Python, census auxiliary data, household survey covariates, EBLUP EBP inputs, 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.

Area-level models start from one direct estimate per area. Unit-level models start from the survey records themselves: they fit a regression of the outcome on covariates across the sampled units, with an effect for each area, and use it to predict every unit of the population that was not sampled. The covariates must therefore be known for every unit of the population, typically from a census, and measured the same way in the survey and the census.

This tutorial prepares the inputs used by every unit-level tutorial that follows: household records from the survey, the same household covariates for every census household, and their district means. The survey, districts, and census are those of the area-level tutorials, so the two approaches can be compared on the same districts.

Setting up the sample

We use the synthetic household survey and census of an imaginary country (World Bank 2023b, 2023a):

import polars as pl
import svy

import svy_sae as sae

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

print(f"Survey households: {hld_data.height:,}")
print(f"Census households: {hld_census.height:,}")
Survey households: 8,000
Census households: 2,501,755

The outcome is total household expenditure per person. A household is poor when it falls below 1,800, the poverty line of the area-level tutorials and of svy’s estimation tutorial.

NoteMethod note

The data also have a pc_exp column, defined differently. We use tot_exp / hhsize so that the poverty rates match those of the area-level tutorials.

Household covariates

The covariates must be built the same way in both files. Define them once and apply the same function to the survey and the census:

yes_no = ["electricity", "car", "refrigerator", "tv", "cell", "bank"]
covariates = [
    "urban", "piped_water", "flush_toilet", *yes_no, "hhsize", "rooms",
]


def household_covariates(frame):
    return frame.with_columns(
        (pl.col("tot_exp") / pl.col("hhsize")).alias("exp_pp"),
        (pl.col("urbrur") == "Urban").cast(pl.Float64).alias("urban"),
        (pl.col("piped_water") != "No piped water").cast(pl.Float64).alias("piped_water"),
        (pl.col("flush_toilet") == "Flush toilet").cast(pl.Float64).alias("flush_toilet"),
        *[(pl.col(c) == "Yes").cast(pl.Float64).alias(c) for c in yes_no],
        pl.col("hhsize").cast(pl.Float64),
        pl.col("rooms").cast(pl.Float64),
    )


hld_data = household_covariates(hld_data)
hld_census = household_covariates(hld_census)

hld_data.select("geo2", "exp_pp", *covariates).head()
shape: (5, 13)
geo2 exp_pp urban piped_water flush_toilet electricity car refrigerator tv cell bank hhsize rooms
str f64 f64 f64 f64 f64 f64 f64 f64 f64 f64 f64 f64
"geo_01_01" 5044.0 1.0 0.0 0.0 1.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0
"geo_01_01" 4652.0 1.0 1.0 0.0 1.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
"geo_01_01" 4978.5 1.0 1.0 1.0 1.0 0.0 1.0 1.0 1.0 1.0 2.0 4.0
"geo_01_01" 5766.5 1.0 1.0 0.0 1.0 0.0 1.0 1.0 1.0 1.0 2.0 1.0
"geo_01_01" 14549.0 1.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 3.0

These are the household-level versions of the district shares used as area-level covariates: a district’s share of households with piped water is the mean of piped_water over its households.

Are they measured alike?

A covariate that means something different in the survey and the census biases every prediction. One check: the survey’s weighted estimate of each covariate should agree with the census value, within sampling error.

hld_design = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
hld_sample = svy.Sample(data=hld_data, design=hld_design)

survey_means = hld_sample.estimation.mean(y=covariates)

pl.DataFrame(
    {
        "covariate": c,
        "survey": est.to_polars()["est"].item(),
        "se": est.to_polars()["se"].item(),
        "census": hld_census[c].mean(),
    }
    for c, est in zip(covariates, survey_means)
).with_columns(((pl.col("survey") - pl.col("census")) / pl.col("se")).alias("z"))
shape: (11, 5)
covariate survey se census z
str f64 f64 f64 f64
"urban" 0.551057 0.006773 0.551057 -1.5080e-12
"piped_water" 0.590111 0.010401 0.570747 1.86172
"flush_toilet" 0.467043 0.013231 0.425793 3.117665
"electricity" 0.82945 0.011873 0.814937 1.222264
"car" 0.247051 0.010566 0.223529 2.226198
… … … … …
"tv" 0.695723 0.015122 0.657454 2.53069
"cell" 0.629976 0.00962 0.607843 2.300723
"bank" 0.572641 0.013986 0.544033 2.045569
"hhsize" 4.026912 0.103417 3.998749 0.272321
"rooms" 2.822785 0.051621 2.678295 2.799049

z is the difference in standard errors. The share of urban households matches exactly, but the survey reports more of every asset than the census, by up to about three standard errors. Differences that large, all in the same direction, are a pattern rather than sampling noise. With real data they call for a look at the questionnaires and reference dates before these covariates enter the model; here we keep them, as the tutorials illustrate the API rather than a full analysis.

Population data

The models need the census in one of two forms. The household records, one row per household:

census_units = hld_census.select("geo1", "geo2", *covariates)

print(f"{census_units.height:,} rows")
2,501,755 rows

Or the district means of the covariates, with the number of households:

census_means = hld_census.group_by("geo2").agg(
    *[pl.col(c).mean() for c in covariates], pl.len().alias("n_hh")
)

census_means.head()
shape: (5, 13)
geo2 urban piped_water flush_toilet electricity car refrigerator tv cell bank hhsize rooms n_hh
str f64 f64 f64 f64 f64 f64 f64 f64 f64 f64 f64 u32
"geo_02_01" 0.0 0.26893 0.176675 0.55396 0.141427 0.295909 0.373368 0.471715 0.259356 5.176675 2.400348 2298
"geo_02_02" 0.039719 0.328546 0.179873 0.64065 0.125495 0.364402 0.451961 0.475733 0.296355 4.332706 2.493364 10096
"geo_02_03" 0.102981 0.357212 0.202843 0.666462 0.129775 0.391624 0.488879 0.518689 0.333487 4.461881 2.411668 19557
"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_02_04" 0.152849 0.365568 0.234068 0.665498 0.140358 0.384591 0.475529 0.499809 0.332877 4.273356 2.441403 31384
Table 1: Population data by model
Model Population data Why
EBLUP (BHF) census_means with pop_size="n_hh", or census_units The district mean is linear in the covariates, so their means suffice
EBP census_units Nonlinear indicators such as poverty rates are simulated for every household

Districts outside the sample

As at the area level, the census covers two districts the survey did not reach:

census_means.join(hld_data.select("geo2").unique(), on="geo2", how="anti").select("geo2", "n_hh")
shape: (2, 2)
geo2 n_hh
str u32
"geo_02_02" 10096
"geo_09_03" 13912

Unit-level models predict them from their census households alone.

The true values

The census also records every household’s expenditure, so the true district mean expenditure and poverty rate are known:

truth = hld_census.group_by("geo2").agg(
    pl.col("exp_pp").mean().alias("true_mean"),
    (pl.col("exp_pp") < 1800).mean().alias("true_rate"),
)

The district mean is the average over households of expenditure per person.

WarningOnly with synthetic data

In practice the truth is unknown. Later tutorials use it to show how close the estimates are, which the model’s own MSE cannot.

Summary

This tutorial prepared the inputs of the unit-level models:

  1. Household covariates built by one function applied to both the survey and the census
  2. A consistency check of the survey’s covariate estimates against the census
  3. Population data as census records for EBP, or district means with sizes for EBLUP
  4. True values of district mean expenditure and poverty rates

Next steps

Fit the first unit-level model
Continue to The BHF Model →

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.