Replicate Weights for Variance Estimation in Python
BRR, Jackknife, Bootstrap, and SDR methods for complex survey designs
Tutorials
Survey Weighting
Variance Estimation
Python
Learn how to create and adjust replicate weights for survey variance estimation using the svy library. Covers BRR, Jackknife (JKn and JK2), Bootstrap, and SDR methods with automatic propagation of weight adjustments.
Replicate weights provide a flexible approach to variance estimation for complex survey designs. Rather than relying on analytical formulas (Taylor linearization), replication methods estimate variance by repeatedly perturbing the sample weights and observing the resulting variation in estimates.
This tutorial assumes you’ve completed the Sample Selection and Weighting tutorials. We use the same World Bank synthetic sample data throughout.
When to Use Replicate Weights
Replicate weights are especially useful when:
Estimating non-linear statistics (medians, percentiles, ratios) where Taylor linearization may be inaccurate
The number of PSUs per stratum is small, making linearization-based variance estimates unstable
Sharing data with secondary analysts who may not have access to design details
Approach
Strengths
Limitations
Taylor linearization
Computationally efficient; no additional columns needed
Requires correct design specification; may be inaccurate for non-linear statistics
Replication
Flexible; works well for non-linear statistics; easy to share
Increases file size; computationally heavier for many replicates
Setting Up
import numpy as npimport polars as plimport svyrng = np.random.default_rng(12345)hld_data = svy.datasets.load(name="hld_sample_wb_2023", source="bundled")hld_sample = svy.Sample( data=hld_data, design=svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight"),)print(hld_sample)
╭────────────── Sample ───────────────╮│Survey Data││ Rows : 825 ││ Columns : 16 ││ Strata : 7 ││ PSUs : 33 ││││Survey Design││ Case id None ││ Wave None ││ Stratum (geo1, urbrur) ││ PSU ea ││ SSU None ││ Weight hhweight ││ With replacement False ││ Prob None ││ Hit None ││ MOS None ││ Population size None ││ Replicate weights None │╰─────────────────────────────────────╯
The World Bank sample includes a household weight (hhweight) computed from the selection probabilities. Although labeled as a survey weight in the dataset, no nonresponse adjustment or calibration has been applied — it is effectively a base weight. This makes it ideal for demonstrating the recommended workflow: create replicates from the base weight first, then apply adjustments that propagate automatically to both the main weight and the replicates.
Replication Methods Overview
svy supports the following replication methods. The API is consistent across all of them — only the statistical method and its requirements differ:
Method
Function
Requirements
Typical Use
Bootstrap (Rao–Wu)
create_bs_wgts()
≥ 2 PSUs per stratum
Most flexible; complex designs; non-linear statistics
Bootstrap (Poisson)
create_bs_wgts(kind="poisson")
A weight column — nothing else
Files that ship no PSU identifier
Jackknife (JKn)
create_jk_wgts()
≥ 1 PSU per stratum
General purpose; # replicates = # PSUs
Jackknife (JK2)
create_jk_wgts(paired=True)
≥ 2 PSUs per stratum
Paired designs; fewer replicates
BRR
create_brr_wgts()
Even # of PSUs per stratum
Balanced half-samples
Fay-BRR
create_brr_wgts(fay_coef=...)
Even # of PSUs per stratum
Damped BRR for stability
SDR
create_sdr_wgts()
Ordered/systematic samples
Systematic samples; time series
All methods share the same core parameters:
Parameter
Description
n_reps
Number of replicates (bootstrap and BRR; JKn/JK2 determine this automatically)
rep_prefix
Prefix for replicate weight column names (defaults to the active weight name)
stratum, psu
The columns to build the replicates over, when they are not the ones the Design names. Defaults to the design’s.
rstate
Random state for reproducibility
drop_nulls
Drop rows with missing values in design columns before creation
BRR and the paired jackknife pair PSUs into variance strata before replicating, and take three further parameters for it:
Parameter
Description
stratum_name
Name of the variance-stratum column the pairing creates (default svy_var_stratum)
order_by
Pair adjacent PSUs in this order — what a systematically sampled frame wants
shuffle
Pair PSUs at random instead (use rstate for reproducibility)
stratum and psu here name the units the replicates are drawn over, which is a separate question from Design.stratum / Design.psu — the analysis design, and what Taylor linearizes over. They coincide when you build replicates yourself from your own design, and often do not on a public-use file. See Declaring producer-supplied replicate weights.
Bootstrap
Bootstrap replication draws PSUs with replacement within each stratum. The Rao-Wu rescaled bootstrap adjusts weights to maintain unbiasedness under the sampling design. It is the most general method — it works with any number of PSUs per stratum (≥ 2) and handles non-linear statistics well.
The replicate weight columns are named automatically from the active weight: hhweight1, hhweight2, …, hhweight500. You can override this with rep_prefix:
For simple statistics (means, totals), 200–500 replicates usually suffice. For percentiles or other non-linear statistics, consider 1,000+. More replicates reduce Monte Carlo error but increase computation time and file size.
Poisson Bootstrap
The Rao–Wu bootstrap resamples PSUs within strata, so it needs a PSU column. Many public-use files do not have one — the producer suppressed it for disclosure control — which leaves the analyst with a weight and nothing to resample. The Beaumont–Patak generalized (Poisson) bootstrap exists for exactly that case: it draws an independent random factor per unit, so a weight column is the only requirement.
poisson_sample = svy.Sample( data=hld_data, design=svy.Design(wgt="hhweight"), # no stratum, no psu).weighting.create_bs_wgts( n_reps=200, kind="poisson", rep_prefix="pb_wgt", rstate=42,)print(repr(poisson_sample.design.rep_wgts))
Both kinds use the same \(1/R\) per-replicate variance coefficient — they differ in how the replicates are drawn, not in how the variance is scaled. Because the Poisson bootstrap draws independent per-unit factors, it has no sampling units by construction and records neither stratum nor psu.
Beaumont, J.-F. and Patak, Z. (2012). On the generalized bootstrap for sample surveys with special attention to Poisson sampling. International Statistical Review, 80(1), 127–148.
Jackknife Methods (JKn and JK2)
The jackknife estimates variance by systematically leaving out one PSU (or one group of PSUs) at a time and observing how the estimate changes. svy supports two variants controlled by the paired parameter.
The examples below each start from a fresh svy.Sample rather than from hld_sample, so every method is built from the original design weights instead of layering on the bootstrap replicates created above.
JKn (delete-one-PSU) creates one replicate per PSU across the entire sample. In each replicate, one PSU is dropped and the remaining PSUs within that stratum are upweighted to compensate. The number of replicates equals the total number of PSUs — with many PSUs this can produce a large number of columns, but the method is very general.
Three things in that output are worth reading carefully:
kind=jkn — the jackknife family, using R’s svrepdesign(type=) names: jk1 is the unstratified delete-one-PSU jackknife, jkn its stratified form, and jk2 the paired one-replicate-per-stratum scheme. They carry different variance coefficients, so this is not just a label.
stratum=('geo1', 'urbrur'), psu='ea' — the units the replicates were actually drawn over, recorded on the weights so a later reader does not have to assume they match the Design.
rep_coefs=0.8 x25, 0.75 x8 (derived) — the per-replicate variance coefficients. JKn uses \((n_h-1)/n_h\), which varies by stratum: 25 replicates come from five-PSU strata (\(4/5\)) and 8 from four-PSU strata (\(3/4\)). svy worked these out from the declared units, which is what (derived) means.
Sample.rep_coefs shows the same thing keyed by the actual column name — the form you can check rather than trust:
coefs = jkn_sample.rep_coefsprint({k: coefs[k] for k inlist(coefs)[:3]})print(f"distinct coefficients: {sorted(set(coefs.values()))}")
JK2 (paired jackknife) is designed for paired PSU designs. It creates one replicate per variance stratum — new strata that pair PSUs together — deleting one PSU and upweighting the others. This produces far fewer replicates than JKn.
Variant
paired
# Replicates
Best for
JKn
False (default)
Total # of PSUs
General purpose; moderate # of PSUs
JK2
True
# of variance strata
Paired designs; fewer replicates
Pairing is done for you. create_jk_wgts(paired=True) builds the variance strata itself, handles the common case of an odd number of PSUs in a stratum (the last three become a triplet), writes the pairing to the column named by stratum_name, and records it on the replicate weights. Your Design.stratum is left exactly as declared, so Taylor estimates keep linearizing over the true strata.
Thirty-three PSUs across 7 strata collapse to 14 variance strata, so JK2 produces 14 replicates against JKn’s 33. The coefficient is 1.0: JK2 has one delete-one replicate per stratum, and the global \((R-1)/R\) would understate the variance by exactly that factor.
Examine the replicate weight pattern — one row per PSU makes it legible. Each replicate drops one PSU from one variance stratum and leaves every other stratum untouched:
Variance stratum 0 is a pair — replicate 1 drops PSU 13033 and doubles 12022. Variance stratum 1 is the triplet that absorbed an odd PSU count — replicate 2 drops 17013 and multiplies the two survivors by \(3/2\) rather than \(2\).
Note
You can control how PSUs are paired with order_by (sort before pairing) or shuffle=True (random pairing, with rstate for reproducibility). Pairing similar-sized PSUs together (e.g., order_by="pop_size") can improve the efficiency of the variance estimate.
Balanced Repeated Replication (BRR)
BRR constructs balanced half-samples using a Hadamard matrix. In each replicate, one PSU is selected from each variance stratum: the selected PSU’s weight is doubled while the other PSU receives zero weight. This creates a systematic set of perturbations that efficiently captures between-PSU variability.
Like the paired jackknife, create_brr_wgts() pairs PSUs into variance strata itself — no pre-step, and no stratum required. Unlike JK2, though, BRR cannot absorb an odd PSU count into a triplet: every stratum must hold an even number of PSUs. Our frame has strata of four and five, so we keep an even number of EAs per stratum for this demonstration:
Twenty-eight PSUs pair into 14 variance strata, and the smallest Hadamard order at least that large is 16 — hence 16 replicates against 14 degrees of freedom.
Fay-BRR
Standard BRR can produce unstable estimates when PSU contributions vary greatly, because it zeros out half the sample in each replicate. Fay-BRR dampens the perturbation using a coefficient \(\rho \in (0, 1)\):
Selected PSU weight: multiplied by \((2 - \rho)\)
Non-selected PSU weight: multiplied by \(\rho\)
Common choices are \(\rho = 0.3\) to \(0.5\). As \(\rho \to 0\), Fay-BRR approaches standard BRR; as \(\rho \to 1\), the perturbation vanishes. The variance coefficient absorbs the damping — \(1 / (R(1-\rho)^2)\) rather than \(1/R\) — so the estimate stays consistent:
fay_sample = svy.Sample( data=brr_data, design=svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight"),).weighting.create_brr_wgts( fay_coef=0.5, stratum_name="var_stratum", rep_prefix="fay_wgt",)rw = fay_sample.design.rep_wgtsprint(repr(rw))print(f"per-replicate coefficient: {rw.coefficients()[0]} (BRR would be {1/ rw.n_reps})")
RepWeights(method=BRR, prefix='fay_wgt', n_reps=16, df=14.0, fay=0.5, wgt='hhweight', stratum='var_stratum', psu='ea')
per-replicate coefficient: 0.25 (BRR would be 0.0625)
Successive Difference Replication (SDR)
SDR is designed for systematic samples where units are ordered (e.g., by geography or time). It creates replicates based on successive differences between adjacent units, which better captures the correlation structure in ordered samples compared to methods that assume independent PSU selection.
The key feature of creating replicates early is that every subsequent weight adjustment automatically propagates to the replicate weights. This ensures that replication-based variance estimates reflect all sources of variability — not just sampling variability, but also uncertainty from nonresponse adjustment, calibration, and trimming.
When you call adjust(), rake(), calibrate(), poststratify(), trim(), or normalize(), svy applies the same adjustment logic to each replicate weight in the same pass as the main weight. The adjusted replicates are renamed to match the new main weight (e.g., nr_wgt → nr_wgt1, nr_wgt2, …).
Example: Full Adjustment Pipeline
Let’s walk through a complete pipeline using the bootstrap sample from above. We simulate nonresponse, then apply adjustment and raking — the replicates follow along automatically.
The design now shows final_wgt as the active weight with final_wgt1, …, final_wgt500 as the replicates — each one individually raked to the same control margins.
Verify: Replicates Match Controls
# Check that a few replicate weights hit the same margins as the main weightfor col in ["final_wgt", "final_wgt1", "final_wgt250", "final_wgt500"]: totals = hld_sample.data.group_by("urbrur").agg( pl.col(col).sum().alias("total") ).sort("urbrur") urban = totals.filter(pl.col("urbrur") =="Urban")["total"][0]print(f" {col}: Urban total = {urban:,.1f} (target: 38,322)")
final_wgt: Urban total = 38,321.9 (target: 38,322)
final_wgt1: Urban total = 38,321.6 (target: 38,322)
final_wgt250: Urban total = 38,322.0 (target: 38,322)
final_wgt500: Urban total = 38,322.0 (target: 38,322)
Skipping Replicate Adjustment
In some cases — preliminary analysis, memory constraints, or externally managed replicates — you may want to adjust only the main weight:
If you skip replicate adjustment (ignore_reps=True), the new weight has no replicate weights: the unadjusted replicates do not go with it, so the design drops them and variance on the new weight is Taylor only. Their columns stay in the data, and update_design(wgt=<previous weight>) brings them back with the weight they belong to.
Design Metadata
When you create replicate weights, svy records what it did on Sample.design.rep_wgts. There is one type per replication method — BootstrapWgts, JackknifeWgts, BrrWgts, SdrWgts — each carrying exactly the parameters its algorithm has, so a bootstrap cannot hold a Fay coefficient and a BRR cannot hold a jackknife kind:
coef_source names the provenance of the variance coefficients:
Value
Meaning
"default"
The method’s standard value, closed-form in n_reps
"derived"
svy computed them (JKn’s per-stratum \((n_h-1)/n_h\)) and cannot recompute them later
"scale"
You asserted them, via scale=
This metadata is what the estimation functions read to compute replicate-based variance. Note that using it is opt-in: the default variance estimator is always Taylor linearization, and replication requires method="replication" at the call site, in estimation as in tabulate(), ttest(), ranktest() and glm.fit().
╭───────────────── Estimate: MEAN (BOOTSTRAP) ──────────────────╮│y: tot_exp ││││ est se lci uci cv (%)││ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ ││ 12,049.1919 701.9875 10,669.9765 13,428.4074 5.83 │││╰───────────────────────────────────────────────────────────────╯
rw.wgt names the full-sample weight the replicates go with. To give a replicate set a new prefix, use wrangling.rename_rep_wgts(): each column keeps its number and padding, and the current design and earlier ones in the design history follow the rename:
Paired with: final_wgt
Old prefix: final_wgt
New prefix: final_rep
Columns: ['final_rep1', 'final_rep2', 'final_rep3'] ...
Declaring Producer-Supplied Replicate Weights
Public-use files usually ship replicate weights already computed — REP1…REP200, BRR_1…BRR_80, and so on. You declare those rather than create them, using the type for the method the producer documents.
To make this concrete without leaving the tutorial, we build a stand-in for such a file: take the JKn replicates from earlier, keep the replicate columns, and rename the design columns to the VARSTRAT / VARUNIT spelling producers commonly publish.
stratum and psu on the replicate weights name the columns the replicates were drawn over. That is a separate question from Design.stratum / Design.psu, which describe the analysis design and are what Taylor linearizes over. The two coincide when you build replicates from your own design; on a public-use file they routinely do not, because producers collapse strata and suppress PSUs for disclosure and publish a distinct pair alongside — or instead of — the design variables.
Declaring them is worth the keystrokes, because they are what svy counts \((n_h-1)/n_h\) from:
No kind was passed, and svy read jkn off the counts — one replicate per PSU across several strata is JKn, one per PSU in a single stratum is JK1, one per stratum is JK2 — then derived the per-stratum coefficients from the same units. The inference is recorded in sample.warnings at INFO level (JACKKNIFE_KIND_UNSPECIFIED), so it is auditable rather than silent. Naming the columns is ordinary survey vocabulary; jk1/jkn/jk2 is not, and requiring both would ask for the same fact twice in the harder of the two languages.
You can still state the family explicitly when the producer documents it:
svy.JackknifeWgts( prefix="REP", n_reps=33, kind="jkn", # "jk1", "jkn", "jk2"; "paired" is an alias for "jk2" stratum="VARSTRAT", psu="VARUNIT",)
A kind the units contradict is refused rather than quietly honoured. Declaring the unstratifiedjk1 against units that name seven strata is a claim svy can check, and getting it wrong overstates standard errors by \(\sqrt{R/n_h}\):
try: svy.Sample( data=pub_file, design=svy.Design( wgt="hhweight", rep_wgts=svy.JackknifeWgts( prefix="REP", n_reps=33, kind="jk1", # says "unstratified"… stratum=("VARSTRAT_A", "VARSTRAT_B"), # …but these say otherwise psu="VARUNIT", ), ), )except svy.MethodError as err:print(err)
❌ Invalid option [INVALID_CHOICE]
Parameter 'rep_wgts.kind' must be one of ['jkn'].
- where: Sample
- param: rep_wgts.kind
- expected: ['jkn']
- got: jk1
Hint: kind='jk1' is the *unstratified* delete-one-PSU jackknife, but the units declared on these weights have 7 strata. Use kind='jkn' (or drop 'kind' and let svy read it off the units); if the replicates really were drawn without regard to strata, say so by naming only psu on the weights and leaving stratum unset.
The estimates then run as usual, with method="replication":
╭───────────────── Estimate: MEAN (JACKKNIFE) ──────────────────╮│y: tot_exp ││││ est se lci uci cv (%)││ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ ││ 12,373.0892 744.3261 10,856.9466 13,889.2317 6.02 │││╰───────────────────────────────────────────────────────────────╯
When the producer publishes the coefficients
Replicate weights often ship with a documented per-replicate variance coefficient that is not the method’s default — typically because it bakes in a finite-population correction. Declare it with scale=, and it replaces the method default:
A scalar is broadcast to every replicate; a sequence must have n_reps entries and is checked at construction. Coming from R, svy’s scale is svrepdesign’s scale × rscales folded into one field: svrepdesign(type = "bootstrap", scale = c, rscales = rep(1, R)) is scale=c here, with no conversion factor.
TipWhich spelling to use
svy.RepWeights(method="jackknife", …) is the door for a method name that arrives as a string — from a codebook, a config file, or a producer’s documentation. It is a factory function that returns the variant for the name.
Code that knows the method when it is written should construct the variant directly — svy.JackknifeWgts(…) — which is typed and autocompletes. And because RepWeights is a function rather than a class, use svy.RepWgts (the union of the four variants) for isinstance checks and type annotations:
rw: svy.RepWgts = design.rep_wgts # ✅ annotationisinstance(rw, svy.RepWgts) # ✅ checkisinstance(rw, svy.RepWeights) # ❌ TypeError — it is a function
Choosing a Replication Method
If your design has…
Recommended method
An even number of PSUs per stratum
BRR (most efficient)
Any number of PSUs, few replicates wanted
JK2 (paired jackknife)
Varying PSUs per stratum
JKn or Bootstrap
Many strata, few PSUs each
Fay-BRR (dampened)
Systematic/ordered sample
SDR
Complex non-linear statistics
Bootstrap (≥ 500 reps)
Need to minimize file size
JK2 (fewest columns)
Multi-PSU strata needing BRR/JK2
Nothing extra — create_brr_wgts() and create_jk_wgts(paired=True) pair PSUs themselves
No PSU identifier in the file
Poisson bootstrap: create_bs_wgts(kind="poisson")
Next Steps
With your replicate weights created and adjusted, you’re ready to compute estimates with proper variance estimation. Continue to the estimation tutorial.
Fay, R. E. (1989). Theory and application of replicate weighting for variance calculations. Proceedings of the Survey Research Methods Section, American Statistical Association, 212–217.
Judkins, D. R. (1990). Fay’s method for variance estimation. Journal of Official Statistics, 6(3), 223–239.
Rao, J. N. K., & Wu, C. F. J. (1988). Resampling inference with complex survey data. Journal of the American Statistical Association, 83(401), 231–241.
Valliant, R., Dever, J. A., & Kreuter, F. (2018). Practical Tools for Designing and Weighting Survey Samples (2nd ed.). Springer.