import numpy as np
import svy
svy.Estimate.PRINT_WIDTH = 95
from rich import print as rprint
rng = np.random.default_rng(12345)
hld_data = svy.datasets.load(name="hld_sample_wb_2023", source="bundled")
hld_design = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
hld_sample = svy.Sample(data=hld_data, design=hld_design)Survey Parameter Estimation in Python with svy
Means, totals, proportions, ratios, and medians with proper variance estimation
Tutorials
Estimation
Variance Estimation
Python
Estimate means, totals, proportions, ratios, and medians from complex survey data in Python. Learn Taylor linearization and replicate-based variance estimation with design effects using the svy library.
Keywords
survey estimation Python, survey mean estimation Python, survey proportion estimation Python, survey total estimation Python, survey ratio estimation Python, Taylor linearization variance Python, design effect DEFF Python, complex survey variance estimation Python, survey confidence interval Python, weighted poverty rate estimation Python, bootstrap variance estimation survey Python, replicate weights estimation Python, BRR jackknife survey Python
The purpose of a survey is to draw inferences about the target population from a sampled subset. Because most surveys employ some combination of stratification, clustering, unequal probabilities of selection, and post-collection adjustments, treating observations as independent and identically distributed (i.i.d.) can misstate uncertainty.
Proper variance estimation uses design information (weights, strata, PSUs) via Taylor linearization or replication methods (BRR, jackknife, bootstrap) to produce accurate standard errors and confidence intervals.
This tutorial demonstrates how to use the svy library to produce point estimates for means, totals, proportions, ratios, and medians, along with corresponding Taylor-based and replication-based measures of uncertainty.
Setting Up the Sample
We’ll use the imaginary country household dataset from World Bank (2023):
Taylor-Based Estimation
Taylor linearization is the standard approach for variance estimation in complex surveys. It accounts for stratification, clustering, and unequal weights without requiring replicate weight columns.
Estimating Means
Estimate the average total household expenditure:
tot_exp_mean = hld_sample.estimation.mean(y="tot_exp")
print(tot_exp_mean)╭─────────────────── Estimate: MEAN (TAYLOR) ───────────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 12,373.0892 744.3261 10,843.1051 13,903.0733 6.02 │ │ │ ╰───────────────────────────────────────────────────────────────╯
Export the result to a Polars DataFrame for further analysis:
tot_exp_mean.to_polars()
shape: (1, 6)
| est | se | lci | uci | cv | df |
|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | i64 |
| 12373.08919 | 744.326051 | 10843.105081 | 13903.073299 | 0.060157 | 26 |
Estimating Several Outcomes at Once
Every estimation method accepts a list of variables for y, computing all of them in one call that indexes the survey design a single time and fans the outcomes out. It returns one Estimate per variable, in order:
indicators = ["tot_exp", "pc_exp", "hhsize"]
estimates = hld_sample.estimation.mean(y=indicators)
for name, est in zip(indicators, estimates):
print(f"{name}:")
print(est)tot_exp:╭─────────────────── Estimate: MEAN (TAYLOR) ───────────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 12,373.0892 744.3261 10,843.1051 13,903.0733 6.02 │ │ │ ╰───────────────────────────────────────────────────────────────╯
pc_exp:╭───────────────── Estimate: MEAN (TAYLOR) ──────────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 3,810.0716 216.5374 3,364.9725 4,255.1706 5.68 │ │ │ ╰────────────────────────────────────────────────────────────╯
hhsize:╭────────── Estimate: MEAN (TAYLOR) ───────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 3.8679 0.2755 3.3016 4.4341 7.12 │ │ │ ╰──────────────────────────────────────────────╯
Collect them into a single tidy table for reporting:
import polars as pl
summary = pl.concat(
est.to_polars().select(pl.lit(name).alias("indicator"), "est", "se", "cv")
for name, est in zip(indicators, estimates)
)
print(summary)shape: (3, 4)
┌───────────┬─────────────┬────────────┬──────────┐
│ indicator ┆ est ┆ se ┆ cv │
│ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ f64 ┆ f64 ┆ f64 │
╞═══════════╪═════════════╪════════════╪══════════╡
│ tot_exp ┆ 12373.08919 ┆ 744.326051 ┆ 0.060157 │
│ pc_exp ┆ 3810.071566 ┆ 216.537414 ┆ 0.056833 │
│ hhsize ┆ 3.867853 ┆ 0.275491 ┆ 0.071226 │
└───────────┴─────────────┴────────────┴──────────┘The same list form works for total(), prop(), ratio(), and median(), and composes with by and where.
Estimating Proportions
We’ll estimate the share of households living below the poverty line.
In this imaginary country, assume the person-level poverty threshold is 1,800 local currency units. We convert it to a household threshold by multiplying by household size—a household is classified as poor if its welfare measure falls below this threshold.
NoteMethod note
This household scaling is for illustration only. Applied analyses commonly use adult-equivalence scales, economies-of-scale adjustments, and regional price indices.
First, create the poverty status variable using mutate():
# Create household poverty line and binary poverty status
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),
}
)
rprint(
hld_sample.show_data(
columns=[
"hid",
"geo1",
"urbrur",
"hhsize",
"tot_exp",
"hhpovline",
"pov_status",
"hhweight",
],
order_type="random",
n=9,
rstate=rng,
),
)shape: (9, 8) ┌─────────────┬────────┬────────┬────────┬─────────┬───────────┬────────────┬──────────┐ │ hid ┆ geo1 ┆ urbrur ┆ hhsize ┆ tot_exp ┆ hhpovline ┆ pov_status ┆ hhweight │ │ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │ │ str ┆ str ┆ str ┆ i64 ┆ i64 ┆ i64 ┆ i32 ┆ f64 │ ╞═════════════╪════════╪════════╪════════╪═════════╪═══════════╪════════════╪══════════╡ │ 2e2b2dab80d ┆ geo_04 ┆ Urban ┆ 4 ┆ 11102 ┆ 7200 ┆ 0 ┆ 84.432 │ │ fa338a1a01a ┆ geo_04 ┆ Urban ┆ 1 ┆ 5551 ┆ 1800 ┆ 0 ┆ 84.432 │ │ e79f5c19794 ┆ geo_04 ┆ Urban ┆ 3 ┆ 16076 ┆ 5400 ┆ 0 ┆ 84.432 │ │ 965e12a3207 ┆ geo_02 ┆ Rural ┆ 3 ┆ 4949 ┆ 5400 ┆ 1 ┆ 51.816 │ │ 82e2017a716 ┆ geo_04 ┆ Rural ┆ 3 ┆ 9870 ┆ 5400 ┆ 0 ┆ 58.992 │ │ e2c284b84bb ┆ geo_04 ┆ Urban ┆ 2 ┆ 10847 ┆ 3600 ┆ 0 ┆ 84.432 │ │ e0e4799d198 ┆ geo_02 ┆ Rural ┆ 5 ┆ 5552 ┆ 9000 ┆ 1 ┆ 51.816 │ │ 396d979d814 ┆ geo_04 ┆ Rural ┆ 3 ┆ 9547 ┆ 5400 ┆ 0 ┆ 58.992 │ │ f038415917c ┆ geo_02 ┆ Urban ┆ 4 ┆ 18610 ┆ 7200 ┆ 0 ┆ 90.25 │ └─────────────┴────────┴────────┴────────┴─────────┴───────────┴────────────┴──────────┘
Estimate the overall poverty rate with the design effect (deff). The argument names which simple-random-sample reference the comparison uses — see Design Effects below for which to choose:
hld_pov_ratio = hld_sample.estimation.prop(y="pov_status", deff="wor", drop_nulls=True)
print(hld_pov_ratio)╭──────────────── Estimate: PROP (TAYLOR, deff=wor) ─────────────────╮ │ │ │ pov_status est se lci uci cv (%) deff │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 0 0.7936 0.0389 0.7022 0.8624 4.91 7.7365 │ │ 1 0.2064 0.0389 0.1376 0.2978 18.86 7.7365 │ │ │ ╰────────────────────────────────────────────────────────────────────╯
hld_pov_ratio.to_polars()
shape: (2, 8)
| pov_status | est | se | lci | uci | cv | deff | df |
|---|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 | i64 |
| "0" | 0.793567 | 0.038938 | 0.70224 | 0.862373 | 0.049067 | 7.736536 | 26 |
| "1" | 0.206433 | 0.038938 | 0.137627 | 0.29776 | 0.188623 | 7.736536 | 26 |
The ci_method parameter controls how confidence intervals for proportions are constructed. The choice matters most when proportions are near 0 or 1 or when sample sizes are small — for moderate proportions with adequate sample sizes, all methods produce similar results. All methods use a \(t\)-quantile with survey degrees of freedom and produce intervals within \([0, 1]\).
ci_method |
Description |
|---|---|
"logit" (default) |
Wald interval on the logit scale, back-transformed to \([0, 1]\) |
"wilson" |
Score-test inversion; coverage closest to nominal across a wide range of scenarios (Franco et al. 2019) |
"beta" |
Korn–Graubard CI using the incomplete Beta function with effective sample size; wider and more conservative |
"korn-graubard" |
Extends "beta" by truncating effective sample size at \(n\) and handling \(p = 0\) and \(p = 1\) |
Estimating Totals
Estimate the total count of poor households:
hld_pov_count = hld_sample.estimation.total(y="pov_status", deff="wor", drop_nulls=True)
print(hld_pov_count)╭────────────────── Estimate: TOTAL (TAYLOR, deff=wor) ───────────────────╮ │ │ │ est se lci uci cv (%) deff │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 11,945.6420 2,253.2203 7,314.0813 16,577.2027 18.86 7.7365 │ │ │ ╰─────────────────────────────────────────────────────────────────────────╯
hld_pov_count.to_polars()
shape: (1, 7)
| est | se | lci | uci | cv | deff | df |
|---|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | f64 | i64 |
| 11945.642 | 2253.220313 | 7314.081315 | 16577.202685 | 0.188623 | 7.736536 | 26 |
Estimating Ratios
Ratio estimation divides the weighted total of one variable by the weighted total of another. This is useful for per-capita measures, expenditure shares, and similar derived quantities. Here we estimate per-capita expenditure as the ratio of total expenditure to household size:
hld_ratio = hld_sample.estimation.ratio(y="tot_exp", x="hhsize")
print(hld_ratio)╭───────────────── Estimate: RATIO (TAYLOR) ─────────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 3,198.9557 192.9835 2,802.2724 3,595.6390 6.03 │ │ │ ╰────────────────────────────────────────────────────────────╯
hld_ratio.to_polars()
shape: (1, 6)
| est | se | lci | uci | cv | df |
|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | i64 |
| 3198.955695 | 192.983516 | 2802.272397 | 3595.638992 | 0.060327 | 26 |
Estimating Medians
Estimate the median household total expenditure:
hld_median_exp = hld_sample.estimation.median(y="tot_exp")
print(hld_median_exp)╭───────────────── Estimate: MEDIAN (TAYLOR) ──────────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 10,608.0000 630.4945 9,239.0000 11,831.0000 5.94 │ │ │ ╰──────────────────────────────────────────────────────────────╯
hld_median_exp.to_polars()
shape: (1, 6)
| est | se | lci | uci | cv | df |
|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | i64 |
| 10608.0 | 630.494497 | 9239.0 | 11831.0 | 0.059436 | 26 |
Design Effects
A design effect is a ratio: the variance your design achieved, over the variance a simple random sample of the same size would have achieved. deff selects which SRS goes in that denominator — it says nothing about how your own sample was drawn.
deff="wor"compares against SRS without replacement. This is Kish’s design effect, and the conventional choice.deff="wr"compares against SRS with replacement — the square of Kish’s “deft”.
The two differ by exactly the finite-population correction \(1 - n/N\). For a national household survey sampling a fraction of a percent, they agree to within a percent and the choice does not matter. As the sampling fraction grows it matters a great deal: at \(n/N = 0.5\) they differ two-fold, and "wor" grows without bound as the sample approaches a census — correctly, since a census has no sampling variance to compare against, but startling if you are not expecting it. Evaluation studies that enrol most of a beneficiary frame live in exactly this regime.
"wor" assumes your weights are reciprocals of selection probabilities. It infers \(N\) from their sum, so once weights have been normalized — and to a lesser extent raked or calibrated — that sum is no longer a population count and the design effect is quietly wrong. svy cannot detect this in general: weights normalized to twice the sample size look perfectly ordinary and produce a plausible, wrong answer. If your weights have been rescaled, use deff="wr", which has no \(N\) in it at all.
The one case svy can prove is when the weights sum to no more than the sample size — \(1 - n/N\) is then zero or negative, and no design effect exists. That raises rather than returning a blank, and the message explains both possible causes.
The reference is recorded on the result, so a design effect is never ambiguous later: it appears in the printed header, as Estimate.deff_ref, and in serialized payloads.
The two references on the same estimate, so the difference is concrete rather than theoretical:
for reference in ("wor", "wr"):
est = hld_sample.estimation.mean(y="tot_exp", deff=reference)
row = est.to_dicts()[0]
print(f"{reference:>4}: deff = {row['deff']:.4f} (header says: {est.deff_ref})") wor: deff = 7.6658 (header says: wor)
wr: deff = 7.5565 (header says: wr)The ratio between them is exactly 1 / (1 - n/N) — here the sampling fraction is small, so they nearly agree. Asking for a design effect on weights that have been rescaled raises rather than returning a blank:
import rich
normalized = hld_sample.weighting.normalize(hld_sample.data.height)
try:
normalized.estimation.mean(y="tot_exp", deff="wor")
except Exception as err:
rich.print(err)╭─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮ │ │ │ ✗ Design effect is not computable for this design [DEFF_NOT_COMPUTABLE] │ │ │ │ The without-replacement reference divides by 1 - n/N, with N taken from the sum of the weights. Here that sum │ │ is no greater than the sample size, so the correction is zero or negative and no design effect exists. Either │ │ the weights have been rescaled -- normalize() makes them sum to the sample size -- and no longer count │ │ population units, or this is a census, in which case there is no sampling variance to compare against. │ │ │ │ where estimation.deff │ │ param deff │ │ expected weights summing to more than the sample size │ │ got sum(weights) <= n │ │ hint Use deff='wr', which compares against a with-replacement reference, needs no population size and is │ │ unaffected by rescaled weights. │ │ │ ╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯
Domain Estimation: by and where
In survey analysis, you often need estimates for subgroups (domains) of the population. svy provides two parameters for this, available on all estimation methods — mean(), prop(), total(), ratio(), and median() — as well as on categorical analysis methods like ttest() and ranktest():
by: produces separate estimates for each level (or combination of levels) of the grouping variable(s)where: restricts estimation to a subpopulation while preserving the full sample for variance calculation
Domain estimation with by
Estimate mean expenditure by wealth quintile — a variable that cuts across the sampling strata:
tot_exp_mean_quint = hld_sample.estimation.mean(y="tot_exp", by="quint_nat")
print(tot_exp_mean_quint)╭────────────────────────── Estimate: MEAN (TAYLOR) ──────────────────────────╮ │ │ │ quint_nat est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ average 9,790.2350 561.3281 8,631.7108 10,948.7593 5.73 │ │ poor 8,146.8227 609.4174 6,875.6003 9,418.0451 7.48 │ │ poorest 6,706.3287 626.9333 5,326.4577 8,086.1997 9.35 │ │ rich 11,722.3861 871.2422 9,928.0291 13,516.7431 7.43 │ │ richest 17,545.6057 1,585.9602 14,213.6271 20,877.5844 9.04 │ │ │ ╰─────────────────────────────────────────────────────────────────────────────╯
Domain estimates export to a tidy DataFrame — one row per domain level:
tot_exp_mean_quint.to_polars()
shape: (5, 7)
| quint_nat | est | se | lci | uci | cv | df |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | i64 |
| "average" | 9790.235028 | 561.328102 | 8631.710766 | 10948.759291 | 0.057336 | 24 |
| "poor" | 8146.822719 | 609.417381 | 6875.600339 | 9418.045099 | 0.074804 | 20 |
| "poorest" | 6706.328714 | 626.93334 | 5326.457737 | 8086.199691 | 0.093484 | 11 |
| "rich" | 11722.386106 | 871.242231 | 9928.029143 | 13516.743069 | 0.074323 | 25 |
| "richest" | 17545.605731 | 1585.960162 | 14213.627071 | 20877.58439 | 0.090391 | 18 |
Multiple grouping variables produce estimates for each combination:
tot_exp_mean_quint_elec = hld_sample.estimation.mean(
y="tot_exp", by=("quint_nat", "electricity")
)
print(tot_exp_mean_quint_elec)╭───────────────────────────────── Estimate: MEAN (TAYLOR) ─────────────────────────────────╮ │ │ │ quint_nat electricity est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ average No 7,912.6937 972.8433 4,816.6722 11,008.7152 12.29 │ │ average Yes 9,969.9220 564.7388 8,804.3584 11,135.4856 5.66 │ │ poor No 5,998.8381 705.0824 4,331.5832 7,666.0930 11.75 │ │ poor Yes 8,738.9583 655.5131 7,361.7764 10,116.1403 7.50 │ │ poorest No 5,924.0546 832.4232 3,955.6866 7,892.4226 14.05 │ │ poorest Yes 7,745.8252 683.5590 6,199.5074 9,292.1431 8.82 │ │ rich No 5,023.4178 973.7919 833.5292 9,213.3063 19.39 │ │ rich Yes 12,072.4782 852.4220 10,316.8822 13,828.0742 7.06 │ │ richest No 5,782.0000 0.0000 nan nan 0.00 │ │ richest Yes 17,580.5853 1,583.7851 14,253.1763 20,907.9943 9.01 │ │ │ ╰───────────────────────────────────────────────────────────────────────────────────────────╯
Poverty rates by wealth quintile:
hld_pov_ratio_quint = hld_sample.estimation.prop(
y="pov_status", by="quint_nat", deff="wor", drop_nulls=True
)
print(hld_pov_ratio_quint)╭────────────────────── Estimate: PROP (TAYLOR, deff=wor) ───────────────────────╮ │ │ │ quint_nat pov_status est se lci uci cv (%) deff │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ average 0 1.0000 0.0000 1.0000 1.0000 0.00 nan │ │ average 1 0.0000 0.0000 0.0000 0.0000 inf nan │ │ poor 0 0.2583 0.0503 0.1677 0.3758 19.46 nan │ │ poor 1 0.7417 0.0503 0.6242 0.8323 6.78 nan │ │ poorest 0 0.0000 0.0000 0.0000 0.0000 inf nan │ │ poorest 1 1.0000 0.0000 1.0000 1.0000 0.00 nan │ │ rich 0 1.0000 0.0000 1.0000 1.0000 0.00 0.0000 │ │ rich 1 0.0000 0.0000 0.0000 0.0000 inf nan │ │ richest 0 1.0000 0.0000 1.0000 1.0000 0.00 0.0000 │ │ richest 1 0.0000 0.0000 0.0000 0.0000 inf nan │ │ │ ╰────────────────────────────────────────────────────────────────────────────────╯
hld_pov_ratio_quint.to_polars()
shape: (10, 9)
| quint_nat | pov_status | est | se | lci | uci | cv | deff | df |
|---|---|---|---|---|---|---|---|---|
| str | str | f64 | f64 | f64 | f64 | f64 | f64 | i64 |
| "average" | "0" | 1.0 | 0.0 | 1.0 | 1.0 | 0.0 | NaN | 24 |
| "average" | "1" | 0.0 | 0.0 | 0.0 | 0.0 | inf | NaN | 24 |
| "poor" | "0" | 0.258331 | 0.050283 | 0.167686 | 0.375847 | 0.194647 | NaN | 20 |
| "poor" | "1" | 0.741669 | 0.050283 | 0.624153 | 0.832314 | 0.067797 | NaN | 20 |
| "poorest" | "0" | 0.0 | 0.0 | 0.0 | 0.0 | inf | NaN | 11 |
| "poorest" | "1" | 1.0 | 0.0 | 1.0 | 1.0 | 0.0 | NaN | 11 |
| "rich" | "0" | 1.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 25 |
| "rich" | "1" | 0.0 | 0.0 | 0.0 | 0.0 | inf | NaN | 25 |
| "richest" | "0" | 1.0 | 0.0 | 1.0 | 1.0 | 0.0 | 0.0 | 18 |
| "richest" | "1" | 0.0 | 0.0 | 0.0 | 0.0 | inf | NaN | 18 |
Subpopulation analysis with where
Unlike filtering the data before analysis (which leads to incorrect standard errors), where preserves the full design structure. It sets the weights of excluded observations to zero but keeps them in the data, matching the behavior of R’s subset() applied to a survey design object.
Conditions are expressed using svy.col():
# Mean expenditure among urban households only
mean_exp_urban = hld_sample.estimation.mean(
y="tot_exp",
where=svy.col("urbrur") == "Urban",
)
print(mean_exp_urban)╭──────────────────── Estimate: MEAN (TAYLOR) ────────────────────╮ │ where: urbrur == "Urban" │ │ │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 14,573.2217 1,011.0039 12,418.3178 16,728.1256 6.94 │ │ │ ╰─────────────────────────────────────────────────────────────────╯
# Mean expenditure among large households (5+ members)
mean_exp_large_hh = hld_sample.estimation.mean(
y="tot_exp",
where=svy.col("hhsize") >= 5,
)
print(mean_exp_large_hh)╭──────────────────── Estimate: MEAN (TAYLOR) ────────────────────╮ │ where: hhsize >= 5 │ │ │ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 15,084.2829 1,329.8489 12,290.3739 17,878.1918 8.82 │ │ │ ╰─────────────────────────────────────────────────────────────────╯
where and by compose naturally. Here we estimate mean expenditure by wealth quintile, restricted to urban households:
mean_exp_quint_urban = hld_sample.estimation.mean(
y="tot_exp",
by="quint_nat",
where=svy.col("urbrur") == "Urban",
)
print(mean_exp_quint_urban)╭────────────────────────── Estimate: MEAN (TAYLOR) ──────────────────────────╮ │ where: urbrur == "Urban" │ │ │ │ │ │ quint_nat est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ average 10,300.0423 636.5568 8,924.8449 11,675.2398 6.18 │ │ poor 10,213.8726 649.0476 8,767.7045 11,660.0407 6.35 │ │ poorest 8,509.0964 876.8918 5,718.4354 11,299.7574 10.31 │ │ rich 12,987.8830 864.0032 11,146.3038 14,829.4621 6.65 │ │ richest 18,150.7339 1,676.3669 14,555.2844 21,746.1833 9.24 │ │ │ ╰─────────────────────────────────────────────────────────────────────────────╯
Complex conditions can be composed using boolean operators:
# Poverty rate among large households without electricity
pov_offgrid_large = hld_sample.estimation.prop(
y="pov_status",
drop_nulls=True,
where=(svy.col("electricity") == "No") & (svy.col("hhsize") >= 5),
)
print(pov_offgrid_large)╭───────────────── Estimate: PROP (TAYLOR) ─────────────────╮ │ where: ([electricity == "No"]) & ([hhsize >= 5]) │ │ │ │ │ │ pov_status est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 0 0.0524 0.0332 0.0112 0.2122 63.47 │ │ 1 0.9476 0.0332 0.7878 0.9888 3.51 │ │ │ ╰───────────────────────────────────────────────────────────╯
Replicate-Based Estimation
Replicate-based variance estimation uses replicate weights (bootstrap, BRR, or jackknife) instead of Taylor linearization. This approach is especially useful for non-linear statistics where linearization may be inaccurate.
TipWhen to prefer replication over Taylor
Replication methods shine when estimating medians, percentiles, or other non-smooth statistics; when the number of PSUs per stratum is very small; or when sharing data with analysts who may not have access to the full design specification. See the Replicate Weights tutorial for a comprehensive guide to creating and adjusting replicate weights.
Setting Up a Replicate Design
There are two ways to set up a replicate design: create replicate weights from an existing design, or declare a design around replicate weight columns that already exist in your data.
Creating Replicate Weights from a Design
The hld_sample does not ship with replicate weight columns, so we create bootstrap weights from the existing design. Bootstrap replication is the most flexible method—it works with any number of PSUs per stratum and handles non-linear statistics well.
hld_rep_sample = hld_sample.weighting.create_bs_wgts(
n_reps=500,
rep_prefix="bs_wgt",
rstate=42,
)
print(hld_rep_sample)╭────────────────── Sample ──────────────────╮ │ Survey Data │ │ Rows : 825 │ │ Columns : 522 │ │ Strata : 7 │ │ PSUs : 33 │ │ │ │ Survey Design │ │ Row index svy_row_index │ │ Stratum (geo1, urbrur) │ │ PSU ea │ │ SSU None │ │ Weight norm_wgt │ │ With replacement False │ │ Prob None │ │ Hit None │ │ MOS None │ │ Population size None │ │ Replicate weights │ │ Method : Bootstrap │ │ Prefix : bs_wgt │ │ N reps : 500 │ │ DF : 499.0 │ ╰────────────────────────────────────────────╯
The design now carries replicate metadata that the estimation methods pick up automatically:
rep_info = hld_rep_sample.design.rep_wgts
print(f"Method: {rep_info.method}")
print(f"Number of replicates: {rep_info.n_reps}")
print(f"Degrees of freedom: {rep_info.df}")Method: Bootstrap
Number of replicates: 500
Degrees of freedom: 499.0Using Pre-Existing Replicate Weights
Many public-use survey files ship with replicate weight columns already computed by the data producer (e.g., columns named repwgt1, repwgt2, …). In that case, you declare the replicate design directly using svy.RepWeights and pass it to svy.Design:
# Example: declaring a design from pre-existing BRR replicate weights
pre_existing_design = svy.Design(
wgt="hhweight",
rep_wgts=svy.RepWeights(
method=svy.EstimationMethod.BRR,
prefix="repwgt",
n_reps=80,
df=39, # Documented by the data producer
fay_coef=0.5, # Fay coefficient, if applicable
),
)
print(pre_existing_design)╭──────────── Design ─────────────╮ │ Row index None │ │ Stratum None │ │ PSU None │ │ SSU None │ │ Weight hhweight │ │ With replacement False │ │ Prob None │ │ Hit None │ │ MOS None │ │ Population size None │ │ Replicate weights │ │ Method : BRR │ │ Prefix : repwgt │ │ N reps : 80 │ │ DF : 39 │ │ Fay coef : 0.5 │ ╰─────────────────────────────────╯
ImportantAlways check the documentation
Data producers should document the replication method, the number of replicates, the Fay coefficient (for BRR), and especially the degrees of freedom. Do not assume a default—specifying the wrong df will produce incorrect confidence intervals and p-values.
Degrees of Freedom
The degrees of freedom (df) control the width of confidence intervals and the reference distribution for hypothesis tests. When you create replicate weights with svy, the default is n_reps - 1 for bootstrap and n_strata for jackknife, which is appropriate in most cases.
However, many survey programs (e.g., the U.S. Census Bureau’s ACS, CPS) document a specific df that reflects the original design—often n_PSUs - n_strata or a similar design-based formula. When the data producer specifies a value, always use it:
# Override df when creating weights
hld_rep_custom_df = hld_sample.weighting.create_bs_wgts(
n_reps=500,
rep_prefix="bs_custom",
rstate=42,
)
# Check the default df
print(f"Default df: {hld_rep_custom_df.design.rep_wgts.df}")Default df: 499.0To override the degrees of freedom after the fact, use design.update():
# Suppose the data producer documents df = 301
updated_design = hld_rep_custom_df.design.update(
rep_wgts=svy.RepWeights(
method=svy.EstimationMethod.BOOTSTRAP,
prefix="bs_custom",
n_reps=500,
df=301,
),
)
print(f"Updated df: {updated_design.rep_wgts.df}")Updated df: 301Variance Center (variance_center)
Replicate-based variance is computed as a weighted sum of squared deviations of replicate estimates from a center value. The variance_center parameter controls what that center is:
variance_center |
Center value | Use case |
|---|---|---|
"rep_mean" (default) |
Mean of the replicate estimates | Standard choice for bootstrap and BRR |
"estimate" |
Full-sample point estimate | Required by some methods (SDR); also called the “conservative” or MSE estimator |
For most applications the default ("rep_mean") is appropriate. The "estimate" option computes a mean-squared-error (MSE) style variance that includes any bias of the replicate distribution relative to the full-sample estimate. SDR replication requires "estimate", and some agencies mandate it for their specific methods.
Estimation with Replicate Weights
The default variance estimation method is always Taylor linearization. To use replication-based variance, pass method="replication" explicitly. This makes the choice clear and reproducible — there’s no implicit switching based on whether replicate weights exist.
Means
rep_mean = hld_rep_sample.estimation.mean(y="tot_exp", method="replication")
print(rep_mean)╭───────────────── Estimate: MEAN (BOOTSTRAP) ──────────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 12,373.0892 737.0626 10,924.9606 13,821.2177 5.96 │ │ │ ╰───────────────────────────────────────────────────────────────╯
rep_mean.to_polars()
shape: (1, 6)
| est | se | lci | uci | cv | df |
|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | i64 |
| 12373.08919 | 737.062597 | 10924.960643 | 13821.217737 | 0.05957 | 499 |
Domain estimation works the same way:
rep_mean_admin1 = hld_rep_sample.estimation.mean(
y="tot_exp", by="geo1", method="replication"
)
print(rep_mean_admin1)╭─────────────────────── Estimate: MEAN (BOOTSTRAP) ───────────────────────╮ │ │ │ geo1 est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ geo_01 14,314.0160 709.5101 12,920.0206 15,708.0114 4.96 │ │ geo_02 13,912.7257 2,202.8153 9,584.7898 18,240.6616 15.83 │ │ geo_03 12,417.6946 1,382.9522 9,700.5677 15,134.8215 11.14 │ │ geo_04 9,934.2525 795.4253 8,371.4571 11,497.0479 8.01 │ │ │ ╰──────────────────────────────────────────────────────────────────────────╯
rep_mean_admin1.to_polars()
shape: (4, 7)
| geo1 | est | se | lci | uci | cv | df |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | i64 |
| "geo_01" | 14314.016 | 709.510132 | 12920.020591 | 15708.011409 | 0.049568 | 499 |
| "geo_02" | 13912.725721 | 2202.815268 | 9584.789821 | 18240.66162 | 0.158331 | 499 |
| "geo_03" | 12417.694581 | 1382.952219 | 9700.56771 | 15134.821451 | 0.111369 | 499 |
| "geo_04" | 9934.252498 | 795.425288 | 8371.457056 | 11497.04794 | 0.080069 | 499 |
Proportions
rep_pov_ratio = hld_rep_sample.estimation.prop(
y="pov_status", drop_nulls=True, method="replication"
)
print(rep_pov_ratio)╭─────────────── Estimate: PROP (BOOTSTRAP) ────────────────╮ │ │ │ pov_status est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 0 0.7936 0.0388 0.7071 0.8596 4.89 │ │ 1 0.2064 0.0388 0.1404 0.2929 18.79 │ │ │ ╰───────────────────────────────────────────────────────────╯
rep_pov_ratio.to_polars()
shape: (2, 7)
| pov_status | est | se | lci | uci | cv | df |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | i64 |
| "0" | 0.793567 | 0.038779 | 0.707127 | 0.859562 | 0.048867 | 499 |
| "1" | 0.206433 | 0.038779 | 0.140438 | 0.292873 | 0.187854 | 499 |
The same ci_method options described in Table 1 apply to replicate-based proportion estimates. The CI method is a post-processing step applied to the point estimate and standard error, so the choice is independent of how the variance was computed.
Proportions by domain:
rep_pov_ratio_admin1 = hld_rep_sample.estimation.prop(
y="pov_status", by="geo1", drop_nulls=True, method="replication"
)
print(rep_pov_ratio_admin1)╭──────────────────── Estimate: PROP (BOOTSTRAP) ────────────────────╮ │ │ │ geo1 pov_status est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ geo_01 0 0.8480 0.0959 0.5639 0.9601 11.31 │ │ geo_01 1 0.1520 0.0959 0.0399 0.4361 63.11 │ │ geo_02 0 0.8271 0.0634 0.6668 0.9195 7.67 │ │ geo_02 1 0.1729 0.0634 0.0805 0.3332 36.68 │ │ geo_03 0 0.7376 0.1012 0.5014 0.8871 13.73 │ │ geo_03 1 0.2624 0.1012 0.1129 0.4986 38.58 │ │ geo_04 0 0.7800 0.0612 0.6376 0.8772 7.85 │ │ geo_04 1 0.2200 0.0612 0.1228 0.3624 27.82 │ │ │ ╰────────────────────────────────────────────────────────────────────╯
rep_pov_ratio_admin1.to_polars()
shape: (8, 8)
| geo1 | pov_status | est | se | lci | uci | cv | df |
|---|---|---|---|---|---|---|---|
| str | str | f64 | f64 | f64 | f64 | f64 | i64 |
| "geo_01" | "0" | 0.848 | 0.095924 | 0.563862 | 0.960119 | 0.113118 | 499 |
| "geo_01" | "1" | 0.152 | 0.095924 | 0.039881 | 0.436138 | 0.631081 | 499 |
| "geo_02" | "0" | 0.827051 | 0.063429 | 0.666772 | 0.919541 | 0.076693 | 499 |
| "geo_02" | "1" | 0.172949 | 0.063429 | 0.080459 | 0.333228 | 0.366752 | 499 |
| "geo_03" | "0" | 0.737566 | 0.101239 | 0.501435 | 0.887052 | 0.137262 | 499 |
| "geo_03" | "1" | 0.262434 | 0.101239 | 0.112948 | 0.498565 | 0.385772 | 499 |
| "geo_04" | "0" | 0.779989 | 0.061197 | 0.637598 | 0.877208 | 0.078459 | 499 |
| "geo_04" | "1" | 0.220011 | 0.061197 | 0.122792 | 0.362402 | 0.278156 | 499 |
Totals
rep_pov_count = hld_rep_sample.estimation.total(
y="pov_status", drop_nulls=True, method="replication"
)
print(rep_pov_count)╭──────────── Estimate: TOTAL (BOOTSTRAP) ────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 170.3070 31.9928 107.4498 233.1642 18.79 │ │ │ ╰─────────────────────────────────────────────────────╯
rep_pov_count.to_polars()
shape: (1, 6)
| est | se | lci | uci | cv | df |
|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | i64 |
| 170.306991 | 31.992791 | 107.449814 | 233.164168 | 0.187854 | 499 |
Ratios
Ratio estimation divides the weighted total of one variable by the weighted total of another:
rep_ratio = hld_rep_sample.estimation.ratio(
y="tot_exp", x="hhsize", method="replication"
)
print(rep_ratio)╭─────────────── Estimate: RATIO (BOOTSTRAP) ────────────────╮ │ │ │ est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 3,198.9557 190.8525 2,823.9822 3,573.9292 5.97 │ │ │ ╰────────────────────────────────────────────────────────────╯
rep_ratio.to_polars()
shape: (1, 6)
| est | se | lci | uci | cv | df |
|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | i64 |
| 3198.955695 | 190.852504 | 2823.982171 | 3573.929219 | 0.059661 | 499 |
Ratios by domain:
rep_ratio_admin1 = hld_rep_sample.estimation.ratio(
y="tot_exp", x="hhsize", by="geo1", method="replication"
)
print(rep_ratio_admin1)╭──────────────────── Estimate: RATIO (BOOTSTRAP) ────────────────────╮ │ │ │ geo1 est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ geo_01 3,434.2649 587.7246 2,279.5451 4,588.9846 17.11 │ │ geo_02 3,439.2063 302.5232 2,844.8301 4,033.5824 8.80 │ │ geo_03 2,675.0031 374.0250 1,940.1453 3,409.8610 13.98 │ │ geo_04 3,401.3023 379.6461 2,655.4005 4,147.2042 11.16 │ │ │ ╰─────────────────────────────────────────────────────────────────────╯
rep_ratio_admin1.to_polars()
shape: (4, 7)
| geo1 | est | se | lci | uci | cv | df |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | i64 |
| "geo_01" | 3434.264875 | 587.724568 | 2279.545143 | 4588.984607 | 0.171135 | 499 |
| "geo_02" | 3439.206273 | 302.523153 | 2844.830144 | 4033.582403 | 0.087963 | 499 |
| "geo_03" | 2675.003135 | 374.024987 | 1940.145251 | 3409.861019 | 0.139822 | 499 |
| "geo_04" | 3401.30234 | 379.646091 | 2655.400506 | 4147.204174 | 0.111618 | 499 |
Comparing Taylor and Replicate Standard Errors
For smooth statistics like means and proportions, Taylor linearization and replication typically produce very similar standard errors. Larger differences may appear for non-linear statistics or when the number of PSUs per stratum is small. Here we compare the two approaches for the overall mean of total expenditure:
import polars as pl
taylor_row = tot_exp_mean.to_polars().select("est", "se", "lci", "uci").with_columns(
pl.lit("Taylor").alias("method")
)
rep_row = rep_mean.to_polars().select("est", "se", "lci", "uci").with_columns(
pl.lit("Bootstrap (500 reps)").alias("method")
)
comparison = pl.concat([taylor_row, rep_row]).select("method", "est", "se", "lci", "uci")
print(comparison)shape: (2, 5)
┌──────────────────────┬─────────────┬────────────┬──────────────┬──────────────┐
│ method ┆ est ┆ se ┆ lci ┆ uci │
│ --- ┆ --- ┆ --- ┆ --- ┆ --- │
│ str ┆ f64 ┆ f64 ┆ f64 ┆ f64 │
╞══════════════════════╪═════════════╪════════════╪══════════════╪══════════════╡
│ Taylor ┆ 12373.08919 ┆ 744.326051 ┆ 10843.105081 ┆ 13903.073299 │
│ Bootstrap (500 reps) ┆ 12373.08919 ┆ 737.062597 ┆ 10924.960643 ┆ 13821.217737 │
└──────────────────────┴─────────────┴────────────┴──────────────┴──────────────┘The point estimates are identical (both use the same base weights), while the standard errors differ slightly due to the different variance estimation strategies.
Next Steps
Real-world designs sometimes have a stratum with a single PSU, which blocks variance estimation. Continue to Singleton PSUs to learn how to detect and resolve them.
Hit a singleton?
Detect and resolve them in Singleton PSUs →
References
Franco, Carolina, Roderick J. A. Little, Thomas A. Louis, and Eric V. Slud. 2019. “Comparative Study of Confidence Intervals for Proportions in Complex Sample Surveys.” Journal of Survey Statistics and Methodology 7 (3): 334–64. https://doi.org/10.1093/jssam/smy019.
World Bank. 2023. “Synthetic Data for an Imaginary Country, Sample, 2023.” World Bank, Development Data Group. https://doi.org/10.48529/MC1F-QH23.