Weighted EBP for Informative Sampling in Python with svy-sae

The pseudo-EBP with survey weights

Tutorials
Unit-level
EBP
Survey Weights
Python
Use survey weights in the EBP for small area poverty estimation in Python with svy-sae: the pseudo-EBP for informative sampling, what the weights change, and its limits.
Author

Mamadou S. Diallo, Ph.D.

Modified

November 5, 2026

Keywords

pseudo EBP survey weights, informative sampling small area estimation, weighted EBP Python, poverty mapping survey weights, complex sampling design SAE

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 EBP assumes that the model fitted to the sample holds for the population: within a district, which households were selected says nothing about their expenditure once the covariates are known. When selection depends on the outcome in a way the covariates do not capture, the sampling is informative, and a model fitted without the weights describes the sample rather than the population.

The pseudo-EBP of Guadarrama, Molina, and Rao (Guadarrama, Molina, and Rao 2018) uses the survey weights to guard against it. It keeps the variance components of the unweighted fit, then estimates the regression coefficients and area effects with the weights, and uses those in the simulation and in every bootstrap replicate.

This tutorial demonstrates how to fit the weighted EBP with weights=, what the weights change, and its limits.

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 of the previous unit-level tutorials:

Rebuild the inputs from Preparing unit-level inputs 
import polars as pl
import svy

import svy_sae as sae

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(svy.datasets.load("hld_sample_wb_2023").data)
hld_census = household_covariates(svy.datasets.load("hld_pop_wb_2023").data)

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

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 household weights vary across strata, from about 160 to 550:

hld_data.select(
    pl.col("hhweight").min().alias("smallest"),
    pl.col("hhweight").max().alias("largest"),
    pl.corr("hhweight", "exp_pp").alias("correlation with expenditure"),
)
shape: (1, 3)
smallest largest correlation with expenditure
f64 f64 f64
156.667553 552.230851 0.01796

Fitting the weighted EBP

Pass the weight column as weights. The weighted EBP takes the log or no transformation:

model = sae.UnitLevel(hld_data)

settings = dict(
    y="exp_pp",
    x=covariates,
    area="geo2",
    pop_data=census_units,
    transformation="log",
    indicators=["headcount", "mean"],
    threshold=1800.0,
    n_mc=100,
    n_reps=50,
    rstate=12345,
)

fit_weighted = model.ebp(weights="hhweight", **settings)
fit_unweighted = model.ebp(**settings)

print(fit_weighted)
╭─────────────────── SAE: EBP (log, weighted) - Unit Level ────────────────────╮
│ Method     REML  Log Likelihood  -1440.6473                                  │
│ No. Areas    61  AIC              2909.2947                                  │
│ No. Obs    8000  BIC              3007.1154                                  │
│                  KIC              2923.2947                                  │
│                                                                              │
│                                                                              │
│  Term           Estimate   Std.Err     z   p-value   [0.025   0.975]         │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━         │
│  Intercept        7.7750       nan   nan       nan      nan      nan         │
│  urban            0.1699       nan   nan       nan      nan      nan         │
│  piped_water     -0.0611       nan   nan       nan      nan      nan         │
│  flush_toilet     0.1529       nan   nan       nan      nan      nan         │
│  electricity      0.0892       nan   nan       nan      nan      nan         │
│  car              0.3571       nan   nan       nan      nan      nan         │
│  refrigerator     0.1675       nan   nan       nan      nan      nan         │
│  tv               0.1411       nan   nan       nan      nan      nan         │
│  cell             0.2444       nan   nan       nan      nan      nan         │
│  bank             0.0283       nan   nan       nan      nan      nan         │
│  hhsize          -0.1567       nan   nan       nan      nan      nan         │
│  rooms            0.0644       nan   nan       nan      nan      nan         │
│                                                                              │
│                                                                              │
│ Variance components: area variance: 0.00258136; residual variance: 0.0820895 │
│                                                                              │
│ Findings:                                                                    │
│ - 17 influential areas                                                       │
╰──────────────────────────────────────────────────────────────────────────────╯

What the weights change

The fixed effects are the weighted coefficients. Compare them with the unweighted ones:

fit_unweighted.to_polars("fixed").select("term", pl.col("est").alias("unweighted")).join(
    fit_weighted.to_polars("fixed").select("term", pl.col("est").alias("weighted")), on="term"
)
shape: (12, 3)
term unweighted weighted
str f64 f64
"Intercept" 7.767783 7.775044
"urban" 0.157851 0.169861
"piped_water" -0.059935 -0.061071
"flush_toilet" 0.152042 0.152877
"electricity" 0.089006 0.089216
… … …
"tv" 0.141204 0.141087
"cell" 0.243403 0.244403
"bank" 0.031311 0.028346
"hhsize" -0.15631 -0.156749
"rooms" 0.067804 0.064402

They differ little. The weighted coefficients are reported without standard errors, which the method does not provide. The variance components, and the diagnostics, are those of the unweighted model:

pl.concat([
    fit_unweighted.to_polars("components").select(pl.lit("unweighted").alias("fit"), "name", "est"),
    fit_weighted.to_polars("components").select(pl.lit("weighted").alias("fit"), "name", "est"),
])
shape: (4, 3)
fit name est
str str f64
"unweighted" "area" 0.002581
"unweighted" "residual" 0.08209
"weighted" "area" 0.002581
"weighted" "residual" 0.08209

Limits

Box–Cox has no weighted estimator for its parameter, so the weighted EBP refuses it:

try:
    model.ebp(weights="hhweight", **(settings | {"transformation": "boxcox", "n_reps": 0}))
except Exception as error:
    print(error)

  ❌ No Box-Cox with survey weights [USAGE_ERROR]
  The weighted EBP is defined for transformation='none' or 'log'; the Box-Cox parameter has no weighted estimator here.
  - param: transformation
  - expected: ['none', 'log']
  - got: boxcox
  Docs: https://svylab.com/docs/svy/reference/codes.html#usage-error

Constant weights give back the unweighted EBP:

point = settings | {"n_reps": 0}
constant = sae.UnitLevel(hld_data.with_columns(pl.lit(1.0).alias("one"))).ebp(
    weights="one", **point
)

constant.to_polars("predictions").join(
    model.ebp(**point).to_polars("predictions"), on=["area", "indicator"], suffix="_unweighted"
).select((pl.col("pred") - pl.col("pred_unweighted")).abs().max().alias("largest difference"))
shape: (1, 1)
largest difference
f64
0.0

How close did we get?

Compare both fits with the true poverty rates and means:

true_values = hld_census.group_by("geo2").agg(
    (pl.col("exp_pp") < 1800).mean().alias("Poverty Headcount"),
    pl.col("exp_pp").mean().alias("Mean"),
).unpivot(index="geo2", variable_name="indicator", value_name="truth")


def scores(name, fit):
    scored = fit.to_polars("predictions").join(
        true_values, left_on=["area", "indicator"], right_on=["geo2", "indicator"]
    )
    covered = (pl.col("truth") >= pl.col("lci")) & (pl.col("truth") <= pl.col("uci"))
    return scored.group_by("indicator").agg(
        pl.lit(name).alias("fit"),
        ((pl.col("pred") - pl.col("truth")) ** 2).mean().sqrt().alias("rmse"),
        covered.sum().alias("interval covers truth"),
    )


pl.concat([scores("unweighted", fit_unweighted), scores("weighted", fit_weighted)]).sort(
    "indicator", "fit"
)
shape: (4, 4)
indicator fit rmse interval covers truth
str str f64 u32
"Mean" "unweighted" 176.308549 35
"Mean" "weighted" 168.115579 37
"Poverty Headcount" "unweighted" 0.024349 43
"Poverty Headcount" "weighted" 0.023684 48

The weights improve the estimates slightly. This survey’s weights barely relate to expenditure, so its sampling is close to non-informative and the two fits agree; the weights matter when the design selects households in a way tied to the outcome. On these data, the choice of transformation matters far more: Box–Cox, which only the unweighted EBP allows, does better still (see the EBP tutorial).

Summary

This tutorial fitted the weighted EBP with ebp(weights=...):

  1. The pseudo-EBP: unweighted variance components, weighted coefficients and area effects
  2. What changes: the fixed effects, reported without standard errors; the diagnostics describe the unweighted model
  3. Limits: the log or no transformation only; constant weights give back the unweighted EBP
  4. Accuracy: a small gain here, as the sampling is close to non-informative

Next steps

This completes the unit-level tutorials.

Back to the overview
Return to the tutorials →

References

Guadarrama, María, Isabel Molina, and J. N. K. Rao. 2018. “Small Area Estimation of General Parameters Under Complex Sampling Designs.” Computational Statistics & Data Analysis 121: 20–40.