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.

Author

Mamadou S. Diallo, Ph.D.

Published

January 18, 2026

Modified

April 18, 2026

Keywords

replicate weights variance estimation Python, bootstrap weights survey Python, jackknife variance estimation Python, BRR balanced repeated replication Python, Fay BRR survey Python, successive difference replication Python, survey variance estimation Python, replicate weight adjustment Python

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 np
import polars as pl
import svy

rng = np.random.default_rng(12345)

hld_data = svy.datasets.load(name="hld_sample_wb_2023")

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

print(hld_sample)

╭────────────── Sample ───────────────╮
│ Survey Data                         │
│   Rows     : 8000                   │
│   Columns  : 52                     │
│   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  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 five replication methods. The API is consistent across all of them — only the statistical method and its requirements differ:

Method	Function	Requirements	Typical Use
Bootstrap	`create_bs_wgts()`	≥ 2 PSUs per stratum	Most flexible; complex designs; non-linear statistics
Jackknife (JKn)	`create_jk_wgts(paired=False)`	≥ 1 PSU per stratum	General purpose; # replicates = # PSUs
Jackknife (JK2)	`create_jk_wgts(paired=True)`	2–3 PSUs per stratum	Paired designs; fewer replicates
BRR	`create_brr_wgts()`	Exactly 2 PSUs per stratum	Balanced half-samples
Fay-BRR	`create_brr_wgts(fay_coef=...)`	Exactly 2 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)
`rstate`	Random state for reproducibility
`drop_nulls`	Drop rows with missing values in design columns before creation

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.

hld_sample = hld_sample.weighting.create_bs_wgts(
    n_reps=500,
    rstate=rng,
)

print(hld_sample)

╭────────────────── Sample ──────────────────╮
│ Survey Data                                │
│   Rows     : 8000                          │
│   Columns  : 552                           │
│   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   : hhweight                  │
│       N reps   : 500                       │
│       DF       : 499.0                     │
╰────────────────────────────────────────────╯

The replicate weight columns are named automatically from the active weight: hhweight1, hhweight2, …, hhweight500. You can override this with rep_prefix:

# Custom prefix
hld_sample = hld_sample.weighting.create_bs_wgts(
    n_reps=500,
    rep_prefix="bs_wgt",
    rstate=rng,
)
# Produces: bs_wgt1, bs_wgt2, ..., bs_wgt500

Examine a few replicate weights alongside the base weight:

print(
    hld_sample.show_data(
        columns=["ea", "geo1", "urbrur", "hhweight", "hhweight1", "hhweight2", "hhweight500"],
        order_type="random",
        rstate=42,
    )
)

shape: (5, 7)
┌───────┬────────┬────────┬────────────┬────────────┬────────────┬─────────────┐
│ ea    ┆ geo1   ┆ urbrur ┆ hhweight   ┆ hhweight1  ┆ hhweight2  ┆ hhweight500 │
│ ---   ┆ ---    ┆ ---    ┆ ---        ┆ ---        ┆ ---        ┆ ---         │
│ i64   ┆ str    ┆ str    ┆ f64        ┆ f64        ┆ f64        ┆ f64         │
╞═══════╪════════╪════════╪════════════╪════════════╪════════════╪═════════════╡
│ 61100 ┆ geo_06 ┆ Rural  ┆ 317.581823 ┆ 0.0        ┆ 635.163645 ┆ 317.581823  │
│ 94126 ┆ geo_09 ┆ Rural  ┆ 317.483899 ┆ 0.0        ┆ 317.483899 ┆ 317.483899  │
│ 43017 ┆ geo_04 ┆ Urban  ┆ 454.456638 ┆ 0.0        ┆ 908.913277 ┆ 908.913277  │
│ 17060 ┆ geo_01 ┆ Urban  ┆ 324.186382 ┆ 648.372764 ┆ 324.186382 ┆ 324.186382  │
│ 44140 ┆ geo_04 ┆ Urban  ┆ 374.628156 ┆ 374.628156 ┆ 749.256311 ┆ 0.0         │
└───────┴────────┴────────┴────────────┴────────────┴────────────┴─────────────┘

Choosing the number of replicates

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.

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.

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.

# JKn: one replicate per PSU
jkn_sample = hld_sample.weighting.create_jk_wgts(paired=False)

JK2 (paired jackknife) is designed for paired PSU designs (2–3 PSUs per variance stratum). It creates one replicate per stratum, where one PSU is deleted and the others upweighted. This produces far fewer replicates than JKn.

Variant	`paired`	# Replicates	Best for
JKn	`False`	Total # of PSUs	General purpose; moderate # of PSUs
JK2	`True`	# of variance strata	Paired designs; fewer replicates

JK2 requires 2–3 PSUs per stratum. Most survey designs have more than that, so you typically need to create variance strata first — new strata that pair PSUs together. The create_variance_strata() method handles this, including the common case where a stratum has an odd number of PSUs (it creates a triplet for the last three):

# Create variance strata suitable for JK2
jk2_sample = hld_sample.weighting.create_variance_strata(
    method="jk2",
    into="var_stratum",
)

print(f"Original strata: {hld_sample.strata.height}")

Original strata: 155

Now create the JK2 replicates:

jk2_sample = jk2_sample.weighting.create_jk_wgts(
    paired=True,
    rep_prefix="jk_wgt",
)

print(jk2_sample)

╭───────────────── Sample ──────────────────╮
│ Survey Data                               │
│   Rows     : 8000                         │
│   Columns  : 708                          │
│   Strata   : 155                          │
│   PSUs     : 320                          │
│                                           │
│ Survey Design                             │
│   Row index                 svy_row_index │
│   Stratum                   var_stratum   │
│   PSU                       ea            │
│   SSU                       None          │
│   Weight                    hhweight      │
│   With replacement          False         │
│   Prob                      None          │
│   Hit                       None          │
│   MOS                       None          │
│   Population size           None          │
│   Replicate weights                       │
│       Method   : Bootstrap                │
│       Prefix   : hhweight                 │
│       N reps   : 500                      │
│       DF       : 499.0                    │
╰───────────────────────────────────────────╯

Examine the replicate weight pattern — notice how each replicate drops one PSU from one variance stratum while leaving all other strata unchanged:

print(
    jk2_sample.show_data(
        columns=["var_stratum", "ea", "hhweight", "jk_wgt1", "jk_wgt2", "jk_wgt3"],
        n=15,
    )
)

shape: (15, 6)
┌─────────────┬───────┬────────────┬────────────┬────────────┬────────────┐
│ var_stratum ┆ ea    ┆ hhweight   ┆ jk_wgt1    ┆ jk_wgt2    ┆ jk_wgt3    │
│ ---         ┆ ---   ┆ ---        ┆ ---        ┆ ---        ┆ ---        │
│ i64         ┆ i64   ┆ f64        ┆ f64        ┆ f64        ┆ f64        │
╞═════════════╪═══════╪════════════╪════════════╪════════════╪════════════╡
│ 3           ┆ 11066 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 │
│ 3           ┆ 11066 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 │
│ 3           ┆ 11087 ┆ 256.192235 ┆ 256.192235 ┆ 256.192235 ┆ 256.192235 │
│ 4           ┆ 11094 ┆ 268.334047 ┆ 268.334047 ┆ 268.334047 ┆ 268.334047 │
│ 4           ┆ 11094 ┆ 268.334047 ┆ 268.334047 ┆ 268.334047 ┆ 268.334047 │
│ …           ┆ …     ┆ …          ┆ …          ┆ …          ┆ …          │
│ 0           ┆ 11024 ┆ 381.859989 ┆ 763.719979 ┆ 381.859989 ┆ 381.859989 │
│ 3           ┆ 11066 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 │
│ 4           ┆ 11112 ┆ 277.440406 ┆ 277.440406 ┆ 277.440406 ┆ 277.440406 │
│ 4           ┆ 11094 ┆ 268.334047 ┆ 268.334047 ┆ 268.334047 ┆ 268.334047 │
│ 3           ┆ 11066 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 ┆ 322.972201 │
└─────────────┴───────┴────────────┴────────────┴────────────┴────────────┘

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 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.

BRR requires exactly 2 PSUs per stratum. If your design has more, use create_variance_strata(method="brr") to pair them first. Unlike JK2, BRR cannot handle odd PSU counts — every stratum must have an even number of PSUs.

The number of replicates defaults to the smallest Hadamard matrix size ≥ the number of strata.

# Pair PSUs first (BRR requires exactly 2 per stratum)
brr_sample = hld_sample.weighting.create_variance_strata(
    method="brr",
    into="var_stratum",
)

# Then create BRR replicates
brr_sample = brr_sample.weighting.create_brr_wgts(rep_prefix="brr_wgt")

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.

# Fay-BRR with ρ = 0.5
# Selected PSUs get weight × 1.5, non-selected get weight × 0.5
fay_sample = brr_sample.weighting.create_brr_wgts(fay_coef=0.5)

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.

# SDR — specify the ordering column
sdr_sample = hld_sample.weighting.create_sdr_wgts(
    n_reps=4,
    rep_prefix="sdr_wgt",
    order_col="sort_order",
)

Automatic Propagation of Weight Adjustments

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.

Step 1: Simulate Response Status

resp_status = rng.choice(
    ("ineligible", "respondent", "non-respondent", "unknown"),
    p=(0.03, 0.82, 0.10, 0.05),
    size=hld_sample.n_records,
)

hld_sample = hld_sample.wrangling.mutate({"resp_status": resp_status})

Step 2: Nonresponse Adjustment

status_mapping = {
    "in": "ineligible",
    "rr": "respondent",
    "nr": "non-respondent",
    "uk": "unknown",
}

hld_sample = hld_sample.weighting.adjust(
    resp_status="resp_status",
    by=("geo1", "geo2"),
    resp_mapping=status_mapping,
    wgt_name="nr_wgt",
    respondents_only=True,
)

print(hld_sample)

╭───────────────── Sample ──────────────────╮
│ Survey Data                               │
│   Rows     : 6629                         │
│   Columns  : 1210                         │
│   Strata   : 155                          │
│   PSUs     : 320                          │
│                                           │
│ Survey Design                             │
│   Row index                 svy_row_index │
│   Stratum                   var_stratum   │
│   PSU                       ea            │
│   SSU                       None          │
│   Weight                    nr_wgt        │
│   With replacement          False         │
│   Prob                      None          │
│   Hit                       None          │
│   MOS                       None          │
│   Population size           None          │
│   Replicate weights                       │
│       Method   : Bootstrap                │
│       Prefix   : nr_wgt                   │
│       N reps   : 500                      │
│       DF       : 499.0                    │
╰───────────────────────────────────────────╯

Notice that the design now shows:

Main weight: nr_wgt
Replicate weights prefix: nr_wgt (i.e., nr_wgt1, …, nr_wgt500)

Both the main weight and all 500 replicates were adjusted in one call.

Step 3: Raking

raking_controls = {
    "statocc": {
        "Occupied for free": 250_000,
        "Owned": 1_932_000,
        "Rented": 325_000,
    },
    "electricity": {"No": 425_000, "Yes": 2_082_000},
}

hld_sample = hld_sample.weighting.rake(
    controls=raking_controls,
    wgt_name="final_wgt",
)

print(hld_sample)

╭───────────────── Sample ──────────────────╮
│ Survey Data                               │
│   Rows     : 6629                         │
│   Columns  : 1711                         │
│   Strata   : 155                          │
│   PSUs     : 320                          │
│                                           │
│ Survey Design                             │
│   Row index                 svy_row_index │
│   Stratum                   var_stratum   │
│   PSU                       ea            │
│   SSU                       None          │
│   Weight                    final_wgt     │
│   With replacement          False         │
│   Prob                      None          │
│   Hit                       None          │
│   MOS                       None          │
│   Population size           None          │
│   Replicate weights                       │
│       Method   : Bootstrap                │
│       Prefix   : final_wgt                │
│       N reps   : 500                      │
│       DF       : 499.0                    │
╰───────────────────────────────────────────╯

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 weight
for col in ["final_wgt", "final_wgt1", "final_wgt250", "final_wgt500"]:
    totals = hld_sample.data.group_by("statocc").agg(
        pl.col(col).sum().alias("total")
    ).sort("statocc")
    owned = totals.filter(pl.col("statocc") == "Owned")["total"][0]
    print(f"  {col}: Owned total = {owned:,.1f} (target: 1,932,000)")

  final_wgt: Owned total = 1,932,000.0 (target: 1,932,000)
  final_wgt1: Owned total = 1,932,000.0 (target: 1,932,000)
  final_wgt250: Owned total = 1,932,000.0 (target: 1,932,000)
  final_wgt500: Owned total = 1,932,000.0 (target: 1,932,000)

Skipping Replicate Adjustment

In some cases — preliminary analysis, memory constraints, or externally managed replicates — you may want to adjust only the main weight:

hld_sample.weighting.normalize(
    controls=1_000,
    wgt_name="norm_wgt",
    ignore_reps=True,  # Skip replicate adjustment
)

Warning

If you skip replicate adjustment (ignore_reps=True), variance estimates computed from those replicates will not account for the additional variability introduced by the weight adjustment. This typically leads to underestimated standard errors.

Design Metadata

When you create replicate weights, svy stores metadata in Sample.design.rep_wgts:

rw = hld_sample.design.rep_wgts
print(f"Method:     {rw.method}")
print(f"Prefix:     {rw.prefix}")
print(f"N reps:     {rw.n_reps}")
print(f"DF:         {rw.df}")
print(f"Fay coef:   {rw.fay_coef}")

Method:     Bootstrap
Prefix:     final_wgt
N reps:     500
DF:         499.0
Fay coef:   0.0

This metadata is used automatically by estimation functions to compute replicate-based variance estimates.

Tip

When replicate weights are provided by the data producer (e.g., in DHS or ACS public-use files), register them directly via the RepWeights configuration in svy.Design() instead of creating them. See the Design object documentation for details.

Choosing a Replication Method

If your design has…	Recommended method
Exactly 2 PSUs per stratum	BRR (most efficient)
2–3 PSUs per stratum	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	`create_variance_strata()` first

Next Steps

With your replicate weights created and adjusted, you’re ready to compute estimates with proper variance estimation. Continue to the estimation tutorial.

Ready to compute estimates?
Learn estimation methods in Survey Estimation →

References

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.