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")
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,048.9638 229.9865 11,596.3788 12,501.5488 1.91 │ │ │ ╰───────────────────────────────────────────────────────────────╯
Export the result to a Polars DataFrame for further analysis:
tot_exp_mean.to_polars()
shape: (1, 5)
| est | se | lci | uci | cv |
|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 |
| 12048.96378 | 229.986492 | 11596.37876 | 12501.5488 | 0.019088 |
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 │ ╞═════════════╪════════╪════════╪════════╪═════════╪═══════════╪════════════╪════════════╡ │ 0d2b023b3d2 ┆ geo_04 ┆ Rural ┆ 7 ┆ 8402 ┆ 12600 ┆ 1 ┆ 156.667553 │ │ fe9f7ae0188 ┆ geo_01 ┆ Urban ┆ 4 ┆ 47047 ┆ 7200 ┆ 0 ┆ 329.043107 │ │ b99c71c17ea ┆ geo_07 ┆ Urban ┆ 5 ┆ 18850 ┆ 9000 ┆ 0 ┆ 277.005333 │ │ c56212914f4 ┆ geo_08 ┆ Rural ┆ 4 ┆ 10692 ┆ 7200 ┆ 0 ┆ 263.047094 │ │ 6efd4139054 ┆ geo_01 ┆ Urban ┆ 3 ┆ 12814 ┆ 5400 ┆ 0 ┆ 301.116939 │ │ 93b4cc76330 ┆ geo_09 ┆ Rural ┆ 5 ┆ 9416 ┆ 9000 ┆ 0 ┆ 210.268025 │ │ 43d7c99a694 ┆ geo_04 ┆ Rural ┆ 5 ┆ 7118 ┆ 9000 ┆ 1 ┆ 156.667553 │ │ 37041d9504f ┆ geo_04 ┆ Urban ┆ 7 ┆ 24284 ┆ 12600 ┆ 0 ┆ 324.971068 │ │ 6d57a16d9d6 ┆ geo_04 ┆ Urban ┆ 3 ┆ 15840 ┆ 5400 ┆ 0 ┆ 374.628156 │ └─────────────┴────────┴────────┴────────┴─────────┴───────────┴────────────┴────────────┘
Estimate the overall poverty rate with the design effect (deff):
hld_pov_ratio = hld_sample.estimation.prop(y="pov_status", deff=True, drop_nulls=True)
print(hld_pov_ratio)╭───────────────────── Estimate: PROP (TAYLOR) ──────────────────────╮ │ │ │ pov_status est se lci uci cv (%) deff │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 0 0.7712 0.0140 0.7425 0.7977 1.82 8.9479 │ │ 1 0.2288 0.0140 0.2023 0.2575 6.13 8.9479 │ │ │ ╰────────────────────────────────────────────────────────────────────╯
hld_pov_ratio.to_polars()
shape: (2, 7)
| pov_status | est | se | lci | uci | cv | deff |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 |
| "0" | 0.771237 | 0.014026 | 0.742474 | 0.797663 | 0.018186 | 8.947936 |
| "1" | 0.228763 | 0.014026 | 0.202337 | 0.257526 | 0.061313 | 8.947936 |
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. The default is "logit", a Wald interval on the logit scale that is back-transformed to \([0, 1]\). Three alternatives are available: "wilson", a score-test inversion with coverage closest to nominal across a wide range of scenarios (Franco et al. 2019); "beta", the Korn-Graubard CI using the incomplete Beta function with effective sample size, which produces wider, more conservative intervals; and "korn-graubard", which extends "beta" by truncating effective sample size at \(n\) and handling \(p = 0\) and \(p = 1\). All methods use a \(t\)-quantile with survey degrees of freedom and produce intervals within \([0, 1]\).
Estimating Totals
Estimate the total count of poor households:
hld_pov_count = hld_sample.estimation.total(y="pov_status", deff=True, drop_nulls=True)
print(hld_pov_count)╭────────────────────────── Estimate: TOTAL (TAYLOR) ──────────────────────────╮ │ │ │ est se lci uci cv (%) deff │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 572,308.8749 36,062.1842 501,342.9490 643,274.8008 6.30 9.4508 │ │ │ ╰──────────────────────────────────────────────────────────────────────────────╯
hld_pov_count.to_polars()
shape: (1, 6)
| est | se | lci | uci | cv | deff |
|---|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 | f64 |
| 572308.874886 | 36062.184185 | 501342.94896 | 643274.800811 | 0.063012 | 9.450779 |
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 (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 2,992.1100 71.2243 2,851.9495 3,132.2706 2.38 │ │ │ ╰───────────────────────────────────────────────────────────╯
hld_ratio.to_polars()
shape: (1, 5)
| est | se | lci | uci | cv |
|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 |
| 2992.110041 | 71.22426 | 2851.949491 | 3132.27059 | 0.023804 |
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,268.0000 220.5454 9,844.0000 10,712.0123 2.15 │ │ │ ╰──────────────────────────────────────────────────────────────╯
hld_median_exp.to_polars()
shape: (1, 5)
| est | se | lci | uci | cv |
|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 |
| 10268.0 | 220.545427 | 9844.0 | 10712.012338 | 0.021479 |
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 tenure status — a variable that cuts across the sampling strata:
tot_exp_mean_tenure = hld_sample.estimation.mean(y="tot_exp", by="statocc")
print(tot_exp_mean_tenure)╭───────────────────────────── Estimate: MEAN (TAYLOR) ─────────────────────────────╮ │ │ │ statocc est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ Occupied for free 10,642.8433 336.8647 9,979.9352 11,305.7514 3.17 │ │ Owned 12,291.7846 251.7569 11,796.3582 12,787.2110 2.05 │ │ Rented 11,683.1705 314.6398 11,063.9982 12,302.3427 2.69 │ │ │ ╰───────────────────────────────────────────────────────────────────────────────────╯
Domain estimates export to a tidy DataFrame — one row per domain level:
tot_exp_mean_tenure.to_polars()
shape: (3, 6)
| statocc | est | se | lci | uci | cv |
|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 |
| "Occupied for free" | 10642.8433 | 336.864689 | 9979.935184 | 11305.751417 | 0.031652 |
| "Owned" | 12291.784592 | 251.756872 | 11796.358152 | 12787.211032 | 0.020482 |
| "Rented" | 11683.17046 | 314.639773 | 11063.998235 | 12302.342686 | 0.026931 |
Multiple grouping variables produce estimates for each combination:
tot_exp_mean_tenure_elec = hld_sample.estimation.mean(
y="tot_exp", by=("statocc", "electricity")
)
print(tot_exp_mean_tenure_elec)╭────────────────────────────────── Estimate: MEAN (TAYLOR) ──────────────────────────────────╮ │ │ │ statocc electricity est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ Occupied for No 5,939.2262 276.2697 5,395.5615 6,482.8908 4.65 │ │ free │ │ Occupied for Yes 11,489.1496 334.3824 10,831.1264 12,147.1729 2.91 │ │ free │ │ Owned No 6,869.8860 190.3772 6,495.2471 7,244.5249 2.77 │ │ Owned Yes 13,554.2039 261.1314 13,040.3296 14,068.0782 1.93 │ │ Rented No 6,837.7979 598.0019 5,661.0040 8,014.5918 8.75 │ │ Rented Yes 12,065.4141 318.9201 11,437.8188 12,693.0095 2.64 │ │ │ ╰─────────────────────────────────────────────────────────────────────────────────────────────╯
Poverty rates by tenure status:
hld_pov_ratio_tenure = hld_sample.estimation.prop(
y="pov_status", by="statocc", deff=True, drop_nulls=True
)
print(hld_pov_ratio_tenure)╭─────────────────────────────── Estimate: PROP (TAYLOR) ────────────────────────────────╮ │ │ │ statocc pov_status est se lci uci cv (%) deff │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ Occupied for free 0 0.8128 0.0201 0.7700 0.8491 2.47 2.1008 │ │ Occupied for free 1 0.1872 0.0201 0.1509 0.2300 10.72 2.1008 │ │ Owned 0 0.7453 0.0161 0.7123 0.7757 2.16 8.4940 │ │ Owned 1 0.2547 0.0161 0.2243 0.2877 6.33 8.4940 │ │ Rented 0 0.8965 0.0125 0.8692 0.9186 1.39 1.7199 │ │ Rented 1 0.1035 0.0125 0.0814 0.1308 12.06 1.7199 │ │ │ ╰────────────────────────────────────────────────────────────────────────────────────────╯
hld_pov_ratio_tenure.to_polars()
shape: (6, 8)
| statocc | pov_status | est | se | lci | uci | cv | deff |
|---|---|---|---|---|---|---|---|
| str | str | f64 | f64 | f64 | f64 | f64 | f64 |
| "Occupied for free" | "0" | 0.812767 | 0.020072 | 0.770038 | 0.849112 | 0.024696 | 2.100841 |
| "Occupied for free" | "1" | 0.187233 | 0.020072 | 0.150888 | 0.229962 | 0.107205 | 2.100841 |
| "Owned" | "0" | 0.745317 | 0.016118 | 0.712322 | 0.775719 | 0.021626 | 8.493961 |
| "Owned" | "1" | 0.254683 | 0.016118 | 0.224281 | 0.287678 | 0.063289 | 8.493961 |
| "Rented" | "0" | 0.896467 | 0.01249 | 0.869183 | 0.918593 | 0.013932 | 1.719937 |
| "Rented" | "1" | 0.103533 | 0.01249 | 0.081407 | 0.130817 | 0.120638 | 1.719937 |
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,437.9184 326.4021 13,793.5402 15,082.2967 2.26 │ │ │ ╰───────────────────────────────────────────────────────────────╯
# 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 (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 13,699.7382 343.4833 13,022.9148 14,376.5616 2.51 │ │ │ ╰───────────────────────────────────────────────────────────────╯
where and by compose naturally. Here we estimate mean expenditure by tenure status, restricted to urban households:
mean_exp_tenure_urban = hld_sample.estimation.mean(
y="tot_exp",
by="statocc",
where=svy.col("urbrur") == "Urban",
)
print(mean_exp_tenure_urban)╭───────────────────────────── Estimate: MEAN (TAYLOR) ─────────────────────────────╮ │ where: urbrur == "Urban" │ │ │ │ │ │ statocc est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ Occupied for free 12,306.6092 438.5758 11,440.7793 13,172.4391 3.56 │ │ Owned 15,415.0163 336.6602 14,750.3868 16,079.6458 2.18 │ │ Rented 12,143.8980 354.5300 11,443.9901 12,843.8059 2.92 │ │ │ ╰───────────────────────────────────────────────────────────────────────────────────╯
Complex conditions can be composed using boolean operators:
# Poverty rate among renters in large households
pov_renter_large = hld_sample.estimation.prop(
y="pov_status",
drop_nulls=True,
where=(svy.col("statocc") == "Rented") & (svy.col("hhsize") >= 5),
)
print(pov_renter_large)╭───────────────── Estimate: PROP (TAYLOR) ─────────────────╮ │ where: ([statocc == "Rented"]) & ([hhsize >= 5]) │ │ │ │ │ │ pov_status est se lci uci cv (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 0 0.7132 0.0340 0.6414 0.7757 4.77 │ │ 1 0.2868 0.0340 0.2243 0.3586 11.85 │ │ │ ╰───────────────────────────────────────────────────────────╯
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 : 8000 │ │ Columns : 554 │ │ Strata : 19 │ │ PSUs : 320 │ │ │ │ Survey Design │ │ Row index svy_row_index │ │ Stratum (geo1, urbrur) │ │ PSU ea │ │ SSU None │ │ Weight hhweight │ │ 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,048.9638 219.7485 11,617.2175 12,480.7100 1.82 │ │ │ ╰───────────────────────────────────────────────────────────────╯
rep_mean.to_polars()
shape: (1, 5)
| est | se | lci | uci | cv |
|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 |
| 12048.96378 | 219.748452 | 11617.217537 | 12480.710023 | 0.018238 |
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,288.0552 632.8557 13,044.6650 15,531.4454 4.43 │ │ geo_02 11,332.3640 632.9519 10,088.7849 12,575.9431 5.59 │ │ geo_03 12,892.2285 701.7075 11,513.5632 14,270.8937 5.44 │ │ geo_04 12,508.1484 567.2538 11,393.6482 13,622.6487 4.54 │ │ geo_05 10,288.6744 624.4159 9,061.8662 11,515.4826 6.07 │ │ geo_06 11,866.6574 677.4616 10,535.6286 13,197.6861 5.71 │ │ geo_07 11,255.0069 893.3712 9,499.7742 13,010.2396 7.94 │ │ geo_08 11,633.2305 1,147.3705 9,378.9578 13,887.5031 9.86 │ │ geo_09 12,265.4656 681.1446 10,927.2008 13,603.7303 5.55 │ │ geo_10 10,460.2557 638.3975 9,205.9774 11,714.5340 6.10 │ │ │ ╰──────────────────────────────────────────────────────────────────────────╯
rep_mean_admin1.to_polars()
shape: (10, 6)
| geo1 | est | se | lci | uci | cv |
|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 |
| "geo_01" | 14288.055196 | 632.855704 | 13044.664994 | 15531.445398 | 0.044293 |
| "geo_02" | 11332.363998 | 632.951851 | 10088.784893 | 12575.943103 | 0.055853 |
| "geo_03" | 12892.228463 | 701.707458 | 11513.563196 | 14270.893729 | 0.054429 |
| "geo_04" | 12508.148448 | 567.253822 | 11393.648192 | 13622.648703 | 0.045351 |
| "geo_05" | 10288.674404 | 624.415856 | 9061.86622 | 11515.482589 | 0.06069 |
| "geo_06" | 11866.657351 | 677.461603 | 10535.628629 | 13197.686074 | 0.05709 |
| "geo_07" | 11255.006899 | 893.371246 | 9499.774155 | 13010.239643 | 0.079375 |
| "geo_08" | 11633.230467 | 1147.370547 | 9378.957831 | 13887.503102 | 0.098629 |
| "geo_09" | 12265.465563 | 681.144559 | 10927.200827 | 13603.730298 | 0.055534 |
| "geo_10" | 10460.255709 | 638.397465 | 9205.977446 | 11714.533972 | 0.061031 |
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.7712 0.0139 0.7427 0.7975 1.81 │ │ 1 0.2288 0.0139 0.2025 0.2573 6.09 │ │ │ ╰───────────────────────────────────────────────────────────╯
rep_pov_ratio.to_polars()
shape: (2, 6)
| pov_status | est | se | lci | uci | cv |
|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 |
| "0" | 0.771237 | 0.013936 | 0.742714 | 0.797461 | 0.018069 |
| "1" | 0.228763 | 0.013936 | 0.202539 | 0.257286 | 0.060918 |
The same ci_method options described in ?@tbl-ci-methods 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.9141 0.0170 0.8743 0.9422 1.86 │ │ geo_01 1 0.0859 0.0170 0.0578 0.1257 19.83 │ │ geo_02 0 0.6584 0.0531 0.5479 0.7540 8.06 │ │ geo_02 1 0.3416 0.0531 0.2460 0.4521 15.54 │ │ geo_03 0 0.7814 0.0357 0.7034 0.8434 4.57 │ │ geo_03 1 0.2186 0.0357 0.1566 0.2966 16.32 │ │ geo_04 0 0.8285 0.0330 0.7538 0.8840 3.98 │ │ geo_04 1 0.1715 0.0330 0.1160 0.2462 19.22 │ │ geo_05 0 0.6768 0.0602 0.5493 0.7824 8.90 │ │ geo_05 1 0.3232 0.0602 0.2176 0.4507 18.63 │ │ geo_06 0 0.7813 0.0490 0.6704 0.8625 6.27 │ │ geo_06 1 0.2187 0.0490 0.1375 0.3296 22.39 │ │ geo_07 0 0.7723 0.0483 0.6642 0.8533 6.25 │ │ geo_07 1 0.2277 0.0483 0.1467 0.3358 21.20 │ │ geo_08 0 0.7472 0.0734 0.5793 0.8639 9.83 │ │ geo_08 1 0.2528 0.0734 0.1361 0.4207 29.05 │ │ geo_09 0 0.6799 0.0458 0.5842 0.7626 6.73 │ │ geo_09 1 0.3201 0.0458 0.2374 0.4158 14.30 │ │ geo_10 0 0.7729 0.0357 0.6953 0.8353 4.62 │ │ geo_10 1 0.2271 0.0357 0.1647 0.3047 15.71 │ │ │ ╰────────────────────────────────────────────────────────────────────╯
rep_pov_ratio_admin1.to_polars()
shape: (20, 7)
| geo1 | pov_status | est | se | lci | uci | cv |
|---|---|---|---|---|---|---|
| str | str | f64 | f64 | f64 | f64 | f64 |
| "geo_01" | "0" | 0.914146 | 0.017021 | 0.874268 | 0.942212 | 0.01862 |
| "geo_01" | "1" | 0.085854 | 0.017021 | 0.057788 | 0.125732 | 0.198258 |
| "geo_02" | "0" | 0.658399 | 0.053094 | 0.547943 | 0.753984 | 0.08064 |
| "geo_02" | "1" | 0.341601 | 0.053094 | 0.246016 | 0.452057 | 0.155426 |
| "geo_03" | "0" | 0.781365 | 0.035678 | 0.703351 | 0.843429 | 0.045662 |
| … | … | … | … | … | … | … |
| "geo_08" | "1" | 0.252795 | 0.073446 | 0.136142 | 0.420721 | 0.290536 |
| "geo_09" | "0" | 0.679925 | 0.045786 | 0.584208 | 0.762562 | 0.067339 |
| "geo_09" | "1" | 0.320075 | 0.045786 | 0.237438 | 0.415792 | 0.143047 |
| "geo_10" | "0" | 0.772869 | 0.03568 | 0.695347 | 0.835337 | 0.046166 |
| "geo_10" | "1" | 0.227131 | 0.03568 | 0.164663 | 0.304653 | 0.157092 |
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 (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 572,308.8749 35,702.0622 502,163.9839 642,453.7659 6.24 │ │ │ ╰─────────────────────────────────────────────────────────────────────╯
rep_pov_count.to_polars()
shape: (1, 5)
| est | se | lci | uci | cv |
|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 |
| 572308.874886 | 35702.062236 | 502163.983892 | 642453.765879 | 0.062383 |
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 (%) │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ 2,992.1100 69.6747 2,855.2182 3,129.0019 2.33 │ │ │ ╰───────────────────────────────────────────────────────────╯
rep_ratio.to_polars()
shape: (1, 5)
| est | se | lci | uci | cv |
|---|---|---|---|---|
| f64 | f64 | f64 | f64 | f64 |
| 2992.110041 | 69.674651 | 2855.218206 | 3129.001876 | 0.023286 |
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,901.5668 226.6067 3,456.3460 4,346.7876 5.81 │ │ geo_02 2,541.5788 219.9765 2,109.3845 2,973.7731 8.66 │ │ geo_03 3,243.7571 197.8127 2,855.1086 3,632.4056 6.10 │ │ geo_04 3,245.3405 194.8500 2,862.5130 3,628.1680 6.00 │ │ geo_05 2,425.0481 177.9375 2,075.4491 2,774.6471 7.34 │ │ geo_06 3,066.4355 262.8349 2,550.0361 3,582.8350 8.57 │ │ geo_07 2,786.6999 200.4976 2,392.7763 3,180.6235 7.19 │ │ geo_08 2,713.7455 263.9658 2,195.1241 3,232.3668 9.73 │ │ geo_09 2,555.7910 173.4632 2,214.9828 2,896.5993 6.79 │ │ geo_10 3,015.4543 209.7698 2,603.3135 3,427.5951 6.96 │ │ │ ╰─────────────────────────────────────────────────────────────────────╯
rep_ratio_admin1.to_polars()
shape: (10, 6)
| geo1 | est | se | lci | uci | cv |
|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 |
| "geo_01" | 3901.566785 | 226.606673 | 3456.345995 | 4346.787575 | 0.058081 |
| "geo_02" | 2541.578792 | 219.976519 | 2109.384459 | 2973.773124 | 0.086551 |
| "geo_03" | 3243.757116 | 197.812725 | 2855.108641 | 3632.405592 | 0.060983 |
| "geo_04" | 3245.340486 | 194.849982 | 2862.512999 | 3628.167972 | 0.06004 |
| "geo_05" | 2425.048094 | 177.937471 | 2075.449114 | 2774.647074 | 0.073375 |
| "geo_06" | 3066.435534 | 262.834915 | 2550.036051 | 3582.835016 | 0.085713 |
| "geo_07" | 2786.699877 | 200.497643 | 2392.776265 | 3180.623489 | 0.071948 |
| "geo_08" | 2713.745458 | 263.965808 | 2195.124076 | 3232.36684 | 0.09727 |
| "geo_09" | 2555.791035 | 173.463202 | 2214.982785 | 2896.599285 | 0.067871 |
| "geo_10" | 3015.454306 | 209.769759 | 2603.313494 | 3427.595117 | 0.069565 |
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 ┆ 12048.96378 ┆ 229.986492 ┆ 11596.37876 ┆ 12501.5488 │
│ Bootstrap (500 reps) ┆ 12048.96378 ┆ 219.748452 ┆ 11617.217537 ┆ 12480.710023 │
└──────────────────────┴─────────────┴────────────┴──────────────┴──────────────┘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
Continue to Categorical Data Analysis to learn how to create weighted cross-tabulations and perform chi-square tests.
Ready for more analysis?
Learn categorical analysis in Categorical Data Analysis →
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.