Python’s svy vs R’s survey: Identical Results Across 18 Estimators

Validation
Survey Methods
Python
R
svy reproduces R’s survey package to six decimals across 18 estimators — Taylor, BRR, jackknife, bootstrap, SDR, cross-tabs, t-tests, and survey GLMs.
Author

Mamadou S. Diallo, Ph.D.

Published

January 10, 2026

Modified

August 10, 2026

Keywords

R survey package alternative, svy vs R survey, complex survey analysis Python, design-based inference Python, replicate weights Python, svyglm Python equivalent, survey GLM Python, Rao-Scott chi-square Python, Taylor linearization, survey weights Python

Summary

TipTL;DR

svy produces numerically identical results to R’s survey package when equivalent survey designs and variance estimators are specified.

Both libraries implement the same design-based inferential framework, including:

  • Complex designs (stratification, clustering, unequal weights)
  • Taylor linearization and replication-based variance estimation
  • Logit-transformed confidence intervals for proportions
  • Categorical data analysis and regression
Estimator Design / Method Match1
Mean Stratified ✅
Mean One-stage cluster ✅
Mean Two-stage cluster (ultimate cluster) ✅
Mean Stratified + clustered ✅
Proportion Logit-transformed CIs ✅
Total Stratified + clustered ✅
Ratio Stratified + clustered ✅
Domain estimation Mean, ratio ✅
BRR Ratio ✅
Jackknife Ratio ✅
Bootstrap Mean ✅
SDR Mean (ACS replicate weights) ✅
Cross-tabulation Rao-Scott χ² (F) ✅
t-test One- and two-sample ✅
GLM Linear, logistic, Poisson ✅

We don’t just assert the match — we quantify it: across the 18 estimators in the Comparison Summary, every estimate and standard error below agrees to at least six decimal places. Together these results validate svy as a statistically equivalent alternative to R’s survey package (Lumley 2010) for complex survey analysis.

Introduction

For decades, design-based inference has run on specialized software—SAS, SPSS, Stata, and R’s survey package. Statistical organizations standardized on those tools for good reason: complex survey estimation is unforgiving, and a variance estimator that is subtly wrong does not fail loudly. It produces confident, plausible, incorrect numbers.

Python has become the default language for data science, but survey work has largely stayed outside it. The practical question for an organization weighing a move is not whether Python is convenient—it plainly is. It is whether a Python implementation computes the same quantities, to the same standard, as software already trusted for official statistics.

This note answers that question with evidence rather than assertion. Running Python’s svy and R’s survey package on the same datasets and the same design specifications, we compare 18 estimators—estimate and standard error, one at a time—and quantify every difference. svy reproduces R to the sixth decimal throughout. Where a difference could arise, it traces to a documented convention (degrees-of-freedom rules, variance centering, family choices), not to the methodology; each such case is called out where it appears.

For a methodologist, the sections that follow are the audit trail—every design specification and every call is shown in both languages, so the comparison can be checked rather than taken on faith. For anyone deciding whether Python belongs in a production statistical workflow, the Numerical Agreement table and the Comparison Summary are the short version.

Validation Scope

This comparison focuses on design-based estimation, including:

  • Means, totals, proportions, and ratios
  • Domain (subpopulation) estimation
  • Taylor linearization variance estimation
  • Variance estimation for multi-stage designs
  • Replication-based variance estimators (BRR, jackknife, bootstrap)
  • Categorical data analysis (cross-tabulation, t-tests)
  • Regression analysis (linear, logistic, Poisson)

R and Python Packages

library(survey)
library(srvyr)
library(gt)
library(dplyr)
library(readr)

data(api)

packageVersion("survey")
[1] '4.5'
import polars as pl
from great_tables import GT

import svy

print(f"svy version: {svy.__version__}")
svy version: 0.24.1

Loading the Datasets

Setting up R environment

nhanes2brr = readr::read_csv("data/nhanes2brr.csv")
nhanes2fay = readr::read_csv("data/nhanes2fay.csv")
nhanes2jknife = readr::read_csv("data/nhanes2jknife.csv")
nmihs_bs = readr::read_csv("data/nmihs_bs.csv")
acs_hak = readr::read_csv("data/psam_h02.csv")
wb_synth_smp = readr::read_csv("data/WLD_2023_SYNTH-SVY-HLD-EN_v01_M.csv")

Setting up Python environment

import svy
import polars as pl

from great_tables import GT

# Set global display precision to 6 decimals
pl.Config.set_float_precision(6)
<class 'polars.config.Config'>

apistrat = svy.io.read_csv("data/apistrat.csv")
apiclus1 = svy.io.read_csv("data/apiclus1.csv")
apiclus2 = svy.io.read_csv("data/apiclus2.csv")

nhanes2brr = svy.io.read_csv("data/nhanes2brr.csv")
nhanes2fay = svy.io.read_csv("data/nhanes2fay.csv")
nhanes2jknife = svy.io.read_csv("data/nhanes2jknife.csv")
nmihs_bs = svy.io.read_csv("data/nmihs_bs.csv")
acs_hak = svy.io.read_csv("data/psam_h02.csv")
wb_synth_smp = svy.io.read_csv("data/WLD_2023_SYNTH-SVY-HLD-EN_v01_M.csv")

Taylor-Based Estimation

Estimating a Mean

Stratified sample

svy Results

design_str = svy.Design(stratum="stype", wgt="pw")
sample_str = svy.Sample(data=apistrat, design=design_str)

api00_mean_str = sample_str.estimation.mean("api00")

cols = ["est", "se", "lci", "uci"]
(
    GT(api00_mean_str.to_polars().select(cols))
    .fmt_number(columns=cols, decimals=6)
)
est se lci uci
662.287363 9.536132 643.481357 681.093370

R Results

design_str <- apistrat |>
  srvyr::as_survey_design(strata = stype, weights = pw)

design_str |>
  summarize(
    est = srvyr::survey_mean(api00, vartype = c("se", "ci"))
  ) |>
  gt() |>
  fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
est est_se est_low est_upp
662.287363 9.536132 643.481357 681.093370

One-stage sample

svy Results

design_clus1 = svy.Design(psu="dnum", wgt="pw")
sample_clus1 = svy.Sample(data=apiclus1, design=design_clus1)

api00_mean_clus1 = sample_clus1.estimation.mean("api00")

cols = ["est", "se", "lci", "uci"]
(
    GT(api00_mean_clus1.to_polars().select(cols))
    .fmt_number(columns=cols, decimals=6)
)
est se lci uci
644.169399 23.779011 593.168493 695.170305

R Results

design_clus1 <- apiclus1 |>
  srvyr::as_survey_design(id = dnum, weights = pw)

design_clus1 |>
  dplyr::summarize(
    est = srvyr::survey_mean(api00, vartype = c("se", "ci"))
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
est est_se est_low est_upp
644.169399 23.779011 593.168493 695.170305

Two-stage sample

svy Results

# Two-stage design: PSUs = districts (dnum), SSUs = schools (snum).
# svy uses the ultimate-cluster estimator, so declaring ssu="snum" leaves
# the variance unchanged (see the note below) — it just documents the design.
design_clus2 = svy.Design(psu="dnum", ssu="snum", wgt="pw")
sample_clus2 = svy.Sample(data=apiclus2, design=design_clus2)

api00_mean_clus2 = sample_clus2.estimation.mean("api00")

cols = ["est", "se", "lci", "uci"]
(
    GT(api00_mean_clus2.to_polars().select(cols))
    .fmt_number(columns=cols, decimals=6)
)
est se lci uci
670.811808 30.711576 608.691782 732.931835

R Results

design_clus2 <- apiclus2 |>
  srvyr::as_survey_design(id = c(dnum, snum), weights = pw)

design_clus2 |>
  dplyr::summarize(
    est = srvyr::survey_mean(api00, vartype = c("se", "ci"))
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
est est_se est_low est_upp
670.811808 30.711576 608.691782 732.931835
TipTwo-stage variance estimation

For multi-stage designs, svy uses the ultimate cluster variance estimator, which approximates total variance using first-stage (PSU) variability only. This approach is standard in survey software (including R’s survey package) because it:

  1. Produces conservative variance estimates
  2. Avoids requiring population sizes at lower stages
  3. Reflects the dominant source of variability in most designs

Accordingly, specifying Design(psu="dnum", ssu="snum") yields the same variance estimates as Design(psu="dnum").

Stratified clustered sample

For this example, we will use The World Bank Synthetic Survey data (World Bank 2023).

svy Results

design_str_clus = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
sample_str_clus = svy.Sample(data=wb_synth_smp, design=design_str_clus)

tot_exp = sample_str_clus.estimation.mean("tot_exp")

cols = ["est", "se", "lci", "uci"]
(GT(tot_exp.to_polars().select(cols)).fmt_number(columns=cols, decimals=6))
est se lci uci
12,048.963780 229.986492 11,596.378760 12,501.548800

R Results

design_str_clus <- wb_synth_smp |>
  dplyr::mutate(stratum = paste(geo1, urbrur, sep = "_")) |>
  srvyr::as_survey_design(id = ea, strata = stratum, weights = hhweight)

design_str_clus |>
  dplyr::summarize(
    est = srvyr::survey_mean(tot_exp, vartype = c("se", "ci"))
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
est est_se est_low est_upp
12,048.963780 229.986492 11,596.378760 12,501.548800

Other Population Parameters

Proportion

svy Results

electricity = sample_str_clus.estimation.prop("electricity")

cols = ["est", "se", "lci", "uci"]
(
    GT(electricity.to_polars().select(cols))
    .fmt_number(columns=cols, decimals=6)
)
est se lci uci
0.170550 0.011873 0.148438 0.195201
0.829450 0.011873 0.804799 0.851562

R Results

design_str_clus |>
  dplyr::group_by(electricity) |>
  dplyr::summarize(
    est = srvyr::survey_prop(vartype = c("se", "ci"), proportion = TRUE)
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
electricity est est_se est_low est_upp
No 0.170550 0.011873 0.148438 0.195201
Yes 0.829450 0.011873 0.804799 0.851562

Total

svy Results

sample_str_clus = sample_str_clus.wrangling.recode(
    cols="electricity", recodes={1: ["No"], 0: ["Yes"]}, into="no_electricity"
)
electricity = sample_str_clus.estimation.total("no_electricity")

cols = ["est", "se", "lci", "uci"]
(GT(electricity.to_polars().select(cols)).fmt_number(columns=cols, decimals=6))
est se lci uci
426,675.251960 30,622.991644 366,412.985386 486,937.518534

R Results

design_str_clus |>
  dplyr::mutate(no_electricity = electricity != "Yes") |>
  dplyr::summarize(
    est = srvyr::survey_total(no_electricity, vartype = c("se", "ci"))
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
est est_se est_low est_upp
426,675.251960 30,622.991644 366,412.985386 486,937.518534

Ratio

svy Results

tot_exp = sample_str_clus.estimation.ratio(y="tot_exp", x="hhsize")

cols = ["est", "se", "lci", "uci"]
(
    GT(tot_exp.to_polars().select(cols))
    .fmt_number(columns=cols, decimals=6)
)
est se lci uci
2,992.110041 71.224260 2,851.949491 3,132.270590

R Results

design_str_clus <- wb_synth_smp |>
  dplyr::mutate(stratum = paste(geo1, urbrur, sep = "_")) |>
  srvyr::as_survey_design(id = ea, strata = stratum, weights = hhweight)

design_str_clus |>
  dplyr::summarize(
    est = srvyr::survey_ratio(tot_exp, hhsize, vartype = c("se", "ci"))
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
est est_se est_low est_upp
2,992.110041 71.224260 2,851.949491 3,132.270590

Domain estimation

Average expenditure by urban and rural areas

svy Results

tot_exp = sample_str_clus.estimation.mean(y="tot_exp", by="urbrur")

cols = ["est", "se", "lci", "uci"]
(
    GT(tot_exp.to_polars().select(["urbrur"] + cols))
    .fmt_number(columns=cols, decimals=6)
)
urbrur est se lci uci
Rural 9,116.629337 305.519957 8,512.322697 9,720.935976
Urban 14,437.918429 326.402120 13,793.540196 15,082.296661

R Results

design_str_clus <- wb_synth_smp |>
  dplyr::mutate(stratum = paste(geo1, urbrur, sep = "_")) |>
  srvyr::as_survey_design(id = ea, strata = stratum, weights = hhweight)

design_str_clus |>
  dplyr::group_by(urbrur) |>
  dplyr::summarize(
    est = srvyr::survey_mean(
      tot_exp,
      vartype = c("se", "ci"),
      # srvyr defaults a grouped CI to the FULL-design df, which overstates
      # precision for a domain and gives intervals that are too narrow.
      # `degf(cur_svy())` counts only the PSUs and strata the domain actually
      # occupies -- the same value base survey returns from
      # `degf(subset(design, ...))`, and what svy reports.
      df = survey::degf(srvyr::cur_svy())
    )
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
urbrur est est_se est_low est_upp
Rural 9,116.629337 305.519957 8,512.322697 9,720.935976
Urban 14,437.918429 326.402120 13,793.540196 15,082.296661

Ratio of expenditure over household size by banking status

svy Results

tot_exp = sample_str_clus.estimation.ratio(y="tot_exp", x="hhsize", by="bank")

cols = ["est", "se", "lci", "uci"]
(
    GT(tot_exp.to_polars().select(["bank"] + cols))
    .fmt_number(columns=cols, decimals=6)
)
bank est se lci uci
No 1,784.529865 41.958370 1,701.943488 1,867.116242
Yes 3,960.323902 89.247824 3,784.665759 4,135.982044

R Results

design_str_clus <- wb_synth_smp |>
  dplyr::mutate(stratum = paste(geo1, urbrur, sep = "_")) |>
  srvyr::as_survey_design(id = ea, strata = stratum, weights = hhweight)

design_str_clus |>
  dplyr::group_by(bank) |>
  dplyr::summarize(
    est = srvyr::survey_ratio(
      tot_exp,
      hhsize,
      vartype = c("se", "ci"),
      # Domain df, as above: srvyr would otherwise use the full-design df.
      df = survey::degf(srvyr::cur_svy())
    )
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )
bank est est_se est_low est_upp
No 1,784.529865 41.958370 1,701.943488 1,867.116242
Yes 3,960.323902 89.247824 3,784.665759 4,135.982044
ImportantDegrees of freedom for domain estimates

Both blocks above pass df = survey::degf(srvyr::cur_svy()). Without it the two libraries disagree on the confidence interval even though the estimate and the standard error match exactly — and the disagreement is easy to misread as a bug.

A domain does not occupy the whole design. Its degrees of freedom should count only the PSUs and strata that actually contain domain members, so svy and R’s own survey::degf() both report a different df per domain, each smaller than the full-design df. How much smaller depends on how the domain cuts across the design: urbrur is part of the stratum definition here, so each PSU falls entirely inside one domain, while bank splits PSUs and leaves most of them contributing to both.

Two R defaults work against this:

  • confint() uses df = Inf — a normal (z) quantile — regardless of the design. The df must always be passed explicitly to get a t-based interval.
  • srvyr’s grouped survey_mean() / survey_ratio() apply the full-design df to every group, because a grouped summary carries no per-group df of its own. degf(cur_svy()) asks for the current group’s design instead.

Using the full-design df for a domain overstates precision, so the resulting interval is too narrow. svy applies the domain df automatically; in R it has to be requested.

Replication-Based estimation

TipDegrees of freedom for replicate weight designs

When only replicate weights are provided (without strata/PSU identifiers), the true design df is unknown:

  • svy defaults to df = n_reps - 1
  • R defaults to the rank of the replicate weight matrix minus 1

Both approaches are heuristics. The rank-based method can detect when post-stratification or calibration has reduced the effective df, but is computationally expensive and numerically sensitive.

Both packages allow user override: RepWeights(df=...) in svy, degf= in R’s svrepdesign().

In practice, data providers typically document the correct degrees of freedom for their replicate weights (e.g., NHANES, ACS). Always consult the survey documentation and specify df explicitly when known.

Balanced Repeated Replication (BRR)

svy Results

rep_weights = svy.RepWeights(method="BRR", prefix="brr_", n_reps=32)
design_brr = svy.Design(wgt="finalwgt", rep_wgts=rep_weights)
sample_brr = svy.Sample(data=nhanes2brr, design=design_brr)

ratio_wgt_hgt = sample_brr.estimation.ratio(
    y="weight",
    x="height",
    method="replication",
)

cols = ["est", "se", "lci", "uci"]
(
    GT(ratio_wgt_hgt.to_polars().select(cols)).fmt_number(
        columns=cols, decimals=6
    )
)
est se lci uci
0.426812 0.000890 0.424996 0.428628

R Results

design_brr <- svrepdesign(
  data = nhanes2brr,
  weights = ~finalwgt,
  repweights = "brr_",
  type = "BRR",
  combined.weights = TRUE
)
ratio_wgt_hgt <- svyratio(~weight, ~height, design = design_brr)

# Extract results into a data frame
est <- coef(ratio_wgt_hgt)
se <- SE(ratio_wgt_hgt)
ci <- confint(ratio_wgt_hgt, df = degf(design_brr))

data.frame(
  est = est,
  se = se,
  lci = ci[1],
  uci = ci[2]
) |>
  gt::gt() |>
  gt::fmt_number(columns = everything(), decimals = 6)
est se lci uci
0.426812 0.000890 0.424996 0.428628
NoteConfidence intervals for replicate designs

R’s confint() defaults to df = Inf — a normal (z) quantile — for every design type, not only replicate-weight ones; see Domain estimation above, where the same default shows up on a Taylor design. svy uses a t quantile throughout, so the df has to be supplied on the R side to compare like with like.

Use confint(..., df = degf(design)). For a replicate design degf() returns the replicate df described above, rather than a PSUs-minus-strata count.

Jackknife

svy Results

rep_weights = svy.RepWeights(
    method="Jackknife", prefix="jkw_", n_reps=62, df=61
)
design_jkn = svy.Design(wgt="finalwgt", rep_wgts=rep_weights)
sample_jkn = svy.Sample(data=nhanes2jknife, design=design_jkn)

ratio_wgt_hgt = sample_jkn.estimation.ratio(
    y="weight", method="replication", x="height"
)

cols = ["est", "se", "lci", "uci"]
(
    GT(ratio_wgt_hgt.to_polars().select(cols)).fmt_number(
        columns=cols, decimals=6
    )
)
est se lci uci
0.426812 0.001247 0.424319 0.429304

R Results

design_jkn <- svrepdesign(
  data = nhanes2jknife,
  weights = ~finalwgt,
  repweights = "jkw_",
  type = "JKn",
  combined.weights = TRUE,
  rscales = rep((62 - 1) / 62, 62)
)
ratio_wgt_hgt <- svyratio(~weight, ~height, design = design_jkn)

# Extract results into a data frame
est <- coef(ratio_wgt_hgt)
se <- SE(ratio_wgt_hgt)
ci <- confint(ratio_wgt_hgt, df = 61)

data.frame(
  est = est,
  se = se,
  lci = ci[1],
  uci = ci[2]
) |>
  gt::gt() |>
  gt::fmt_number(columns = everything(), decimals = 6)
est se lci uci
0.426812 0.001247 0.424319 0.429304

Bootstrap

svy Results

rep_weights = svy.RepWeights(method="bootstrap", prefix="bsrw", n_reps=1000)
design_bs = svy.Design(wgt="finwgt", rep_wgts=rep_weights)
sample_bs = svy.Sample(data=nmihs_bs, design=design_bs)

mean_birth_weight = sample_bs.estimation.mean(
    y="birthwgt", method="replication", drop_nulls=True
)

cols = ["est", "se", "lci", "uci"]
(
    GT(mean_birth_weight.to_polars().select(cols)).fmt_number(
        columns=cols, decimals=6
    )
)
est se lci uci
3,355.452419 6.520638 3,342.656702 3,368.248137

R Results

design_bs <- svrepdesign(
  data = nmihs_bs,
  weights = ~finwgt,
  repweights = "bsrw",
  type = "bootstrap",
  replicates = 1000,
  combined.weights = TRUE,
  rscales = rep((1000 - 1) / 1000, 1000)
)

mean_birth_weight <- svymean(~birthwgt, design = design_bs, na.rm = TRUE)

est <- coef(mean_birth_weight)
se <- SE(mean_birth_weight)
ci <- confint(mean_birth_weight, df = 999)

data.frame(
  est = est,
  se = se,
  lci = ci[1],
  uci = ci[2]
) |>
  gt::gt() |>
  gt::fmt_number(columns = everything(), decimals = 6)
est se lci uci
3,355.452419 6.520638 3,342.656702 3,368.248137

Successive Difference Replication (SDR)

The American Community Survey (ACS) provides 80 replicate weights constructed using successive difference replication (SDR). To illustrate SDR, we will use data from the 2024 American Community Survey (ACS) 1-Year Public Use Microdata Sample2.

ACS replicate weights use SDR with 80 replicates (e.g., WGTP1–WGTP80) alongside the main weight WGTP. The ACS documentation describes the SDR replicate-weight construction and recommended variance estimation practice.

svy Results

rep_weights_acs = svy.RepWeights(method="sdr", prefix="WGTP", n_reps=80)
design_acs = svy.Design(wgt="WGTP", rep_wgts=rep_weights_acs)
sample_acs = svy.Sample(data=acs_hak, design=design_acs)


mean_hincp = sample_acs.estimation.mean(
    y="HINCP",
    method="replication",
    drop_nulls=True,
)

cols = ["est", "se", "lci", "uci"]
(GT(mean_hincp.to_polars().select(cols)).fmt_number(columns=cols, decimals=6))
est se lci uci
111,770.504058 2,517.264331 106,760.014741 116,780.993375

R Results

R’s survey supports SDR directly via type="successive-difference". It also includes a dedicated type="ACS" shortcut that applies ACS-specific defaults. In practice, both should agree when equivalent settings are used.

design_sdr <- svrepdesign(
  data = acs_hak,
  weights = ~WGTP,
  repweights = "^WGTP[0-9]+",
  type = "successive-difference",
  scale = 4 / 80,
  combined.weights = TRUE,
  rscales = 1,
)

mean_hincp_sdr <- svymean(~HINCP, design = design_sdr, na.rm = TRUE)

est <- coef(mean_hincp_sdr)
se <- SE(mean_hincp_sdr)
ci <- confint(mean_hincp_sdr, df = 79)

data.frame(
  est = est,
  se = se,
  lci = ci[1],
  uci = ci[2]
) |>
  gt::gt() |>
  gt::fmt_number(columns = everything(), decimals = 6)
est se lci uci
111,770.504058 2,517.264331 106,760.014741 116,780.993375

The type = "ACS" shortcut applies the same successive-difference construction. One caveat: it defaults to mse = TRUE (centering the replicate variance on the full-sample estimate, the ACS-recommended convention), while the explicit construction above uses replicate-mean centering. We pass mse = FALSE here so the two are exactly equivalent — see the callout below for switching either package to full-sample centering.

design_acs <- svrepdesign(
  data = acs_hak,
  weights = ~WGTP,
  repweights = "^WGTP[0-9]+",
  type = "ACS",
  combined.weights = TRUE,
  mse = FALSE,
)

mean_hincp_acs <- svymean(~HINCP, design = design_acs, na.rm = TRUE)

est <- coef(mean_hincp_acs)
se <- SE(mean_hincp_acs)
ci <- confint(mean_hincp_acs, df = 79)

data.frame(
  est = est,
  se = se,
  lci = ci[1],
  uci = ci[2]
) |>
  gt::gt() |>
  gt::fmt_number(columns = everything(), decimals = 6)
est se lci uci
111,770.504058 2,517.264331 106,760.014741 116,780.993375
TipReplicate Variance Calculation

By default, both svy and R survey use the average replicate estimates for calculating the estimated variance.

If, instead you want to use the full sample estimate:

  • Use rep_center = "estimate" with svy
  • Use mse = TRUE with R survey

Categorical Data Analysis

Let’s use the World Bank dataset to demonstrate categorical data analysis.

Cross-tabulation

Below, we compute the cross-tabulation of urban/rural and electricity access and show the Rao-Scott χ² test.

svy Results

crosstab = sample_str_clus.categorical.tabulate(
    rowvar="urbrur",
    colvar="electricity",
    units="percent",
)

cols = ["est", "se", "lci", "uci"]
(
    GT(
        crosstab.to_polars().select(["urbrur", "electricity"] + cols)
    ).fmt_number(columns=cols, decimals=6)
)
urbrur electricity est se lci uci
Rural No 15.198691 1.133601 13.099580 17.566174
Rural Yes 29.695594 1.197631 27.394003 32.105052
Urban No 1.856347 0.373584 1.247665 2.753698
Urban Yes 53.249369 0.744614 51.781676 54.711460 
test_stat = crosstab.stats.f

# Create a formatted dataframe
test_df = pl.DataFrame(
    {
        "statistic": ["Pearson χ² (adjusted)"],
        "F_value": [test_stat.value],
        "df_num": [test_stat.df_num],
        "df_den": [test_stat.df_den],
        "p_value": [test_stat.p_value],
    }
)

cols = ["F_value", "df_num", "df_den", "p_value"]
GT(test_df).fmt_number(columns=cols, decimals=6)
statistic F_value df_num df_den p_value
Pearson χ² (adjusted) 193.172687 1.000000 301.000000 0.000000

R Results

# Cell percentages (scaled to sum to 100), matching svy's units="percent"
survey::svytable(~ urbrur + electricity, design_str_clus, Ntotal = 100) |>
  as.data.frame() |>
  gt::gt() |>
  gt::fmt_number(columns = Freq, decimals = 6)
urbrur electricity Freq
Rural No 15.198691
Urban No 1.856347
Rural Yes 29.695594
Urban Yes 53.249369 
# Rao-Scott second-order corrected Pearson chi-square, reported as an F,
# matching svy's crosstab.stats.f
chi <- survey::svychisq(
  ~ urbrur + electricity,
  design_str_clus,
  statistic = "F"
)

data.frame(
  statistic = "Pearson Chi-square (adjusted)",
  F_value = as.numeric(chi$statistic),
  df_num = as.numeric(chi$parameter[1]),
  df_den = as.numeric(chi$parameter[2]),
  p_value = as.numeric(chi$p.value)
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(F_value, df_num, df_den, p_value),
    decimals = 6
  )
statistic F_value df_num df_den p_value
Pearson Chi-square (adjusted) 193.172687 1.000000 301.000000 0.000000

Per-cell standard errors

svytable() reports cell estimates only, so to compare the per-cell standard errors we use svymean() on the interaction of the two factors—the same cell-proportion estimator svy’s tabulate uses, scaled to percent.

cell_pct <- survey::svymean(
  ~ interaction(urbrur, electricity),
  design_str_clus
)

cells_r <- data.frame(
  cell = names(coef(cell_pct)),
  est = as.numeric(coef(cell_pct)) * 100,
  se = as.numeric(SE(cell_pct)) * 100
)
cells_r$cell <- gsub(
  "interaction(urbrur, electricity)", "", cells_r$cell,
  fixed = TRUE
)

cells_r |>
  gt::gt() |>
  gt::fmt_number(columns = c(est, se), decimals = 6)
cell est se
Rural.No 15.198691 1.133601
Urban.No 1.856347 0.373584
Rural.Yes 29.695594 1.197631
Urban.Yes 53.249369 0.744614
NoteCell estimates, standard errors, and the χ² all match

Every cell agrees with svy on both the estimate and the standard error, and the Rao-Scott adjusted χ² matches to six decimals.

A cell percentage is a ratio of two estimated totals, so its variance uses the centered (Hájek) linearization—R subtracts the cell proportion before forming the score (sweep(x, 2, average) in svymean), and svy does the same. Treating the denominator as a fixed constant instead would inflate the standard error by a p-dependent amount.3

T-tests

One group

svy Results

tot_exp_test1 = sample_str_clus.categorical.ttest(
    y="tot_exp",
    mean_h0=12500,
)

print(tot_exp_test1.to_polars().drop("y"))
shape: (1, 7)
┌─────────────┬────────────┬─────────────┬──────────┬───────────┬────────────┬──────────┐
│ diff        ┆ se         ┆ lci         ┆ uci      ┆ t         ┆ df         ┆ p_value  │
│ ---         ┆ ---        ┆ ---         ┆ ---      ┆ ---       ┆ ---        ┆ ---      │
│ f64         ┆ f64        ┆ f64         ┆ f64      ┆ f64       ┆ f64        ┆ f64      │
╞═════════════╪════════════╪═════════════╪══════════╪═══════════╪════════════╪══════════╡
│ -451.036220 ┆ 229.986492 ┆ -903.627330 ┆ 1.554890 ┆ -1.961142 ┆ 300.000000 ┆ 0.050787 │
└─────────────┴────────────┴─────────────┴──────────┴───────────┴────────────┴──────────┘
tot_exp_test1 <- svyttest((tot_exp - 12500) ~ 0, design_str_clus)

test1_df <- data.frame(
  test = "One-sample t-test",
  statistic = tot_exp_test1$statistic,
  df = tot_exp_test1$parameter,
  p_value = tot_exp_test1$p.value,
  mean_diff = tot_exp_test1$estimate,
  ci_lower = tot_exp_test1$conf.int[1],
  ci_upper = tot_exp_test1$conf.int[2]
)

test1_df |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(statistic, df, mean_diff, ci_lower, ci_upper, p_value),
    decimals = 6
  )
test statistic df p_value mean_diff ci_lower ci_upper
One-sample t-test −1.961142 300.000000 0.050787 −451.036220 −903.627330 1.554890

Two groups

svy Results

tot_exp_test2 = sample_str_clus.categorical.ttest(
    y="tot_exp",
    group="urbrur",
)

print(tot_exp_test2.to_polars().drop(["y", "group_var", "paired"]))
shape: (1, 7)
┌─────────────┬────────────┬─────────────┬─────────────┬───────────┬────────────┬──────────┐
│ diff        ┆ se         ┆ lci         ┆ uci         ┆ t         ┆ df         ┆ p_value  │
│ ---         ┆ ---        ┆ ---         ┆ ---         ┆ ---       ┆ ---        ┆ ---      │
│ f64         ┆ f64        ┆ f64         ┆ f64         ┆ f64       ┆ f64        ┆ f64      │
╞═════════════╪════════════╪═════════════╪═════════════╪═══════════╪════════════╪══════════╡
│ 5321.289092 ┆ 447.080293 ┆ 4441.478438 ┆ 6201.099746 ┆ 11.902312 ┆ 300.000000 ┆ 0.000000 │
└─────────────┴────────────┴─────────────┴─────────────┴───────────┴────────────┴──────────┘

R Results

tot_exp_test2 <- svyttest(tot_exp ~ urbrur, design_str_clus)

test2_df <- data.frame(
  test = "Two-sample t-test",
  statistic = tot_exp_test2$statistic,
  df = tot_exp_test2$parameter,
  p_value = tot_exp_test2$p.value,
  mean_diff = tot_exp_test2$estimate,
  ci_lower = tot_exp_test2$conf.int[1],
  ci_upper = tot_exp_test2$conf.int[2]
)

test2_df |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(statistic, df, mean_diff, p_value, ci_lower, ci_upper),
    decimals = 6
  )
test statistic df p_value mean_diff ci_lower ci_upper
Two-sample t-test 11.902312 300.000000 0.000000 5,321.289092 4,441.478438 6,201.099746

Generalized Linear Models (GLMs)

We use the World Bank synthetic survey dataset to compare GLM results. First, we create the poverty indicator and rename variables for consistency.

svy Setup

# Create derived variables for GLM
exp_pc = sample_str_clus.data["tot_exp"] / sample_str_clus.data["hhsize"]
poverty_line = float(exp_pc.median()) * 0.60

glm_sample = sample_str_clus.wrangling.mutate(
    {
        "is_poor": svy.when(
            svy.col("tot_exp") / svy.col("hhsize") < poverty_line
        ).then(1).otherwise(0),
    }
)

R Setup

# Create derived variables for GLM
design_glm <- wb_synth_smp |>
  dplyr::mutate(
    stratum = paste(geo1, urbrur, sep = "_"),
    exp_pc = tot_exp / hhsize,
    is_poor = as.integer(exp_pc < median(exp_pc) * 0.60)
  ) |>
  srvyr::as_survey_design(id = ea, strata = stratum, weights = hhweight)

Linear Regression

svy Results

lin_model = glm_sample.glm.fit(
    y="tot_exp",
    x=["hhsize", "rooms", svy.Cat("urbrur")],
    family="gaussian",
)

cols = [
    "term",
    "estimate",
    "std_err",
    "statistic",
    "p_value",
    "conf_low",
    "conf_high",
]
(
    GT(lin_model.to_polars().select(cols))
    .fmt_number(
        columns=["estimate", "std_err", "statistic", "conf_low", "conf_high"],
        decimals=6,
    )
    .fmt_number(columns="p_value", decimals=6)
)
term estimate std_err statistic p_value conf_low conf_high
_intercept_ 518.294002 348.088374 1.488972 0.137552 −166.728779 1,203.316783
hhsize 825.812641 55.068628 14.996064 0.000000 717.439977 934.185306
rooms 1,972.989241 143.261134 13.771978 0.000000 1,691.057561 2,254.920921
urbrur_Urban 4,783.297099 270.632985 17.674479 0.000000 4,250.703154 5,315.891043

R Results

lin_model_r <- svyglm(
  tot_exp ~ hhsize + rooms + urbrur,
  design = design_glm,
  family = gaussian()
)

lin_coefs <- summary(lin_model_r)$coefficients
lin_ci <- confint(lin_model_r)

data.frame(
  term = rownames(lin_coefs),
  coef = lin_coefs[, "Estimate"],
  se = lin_coefs[, "Std. Error"],
  t = lin_coefs[, "t value"],
  p_value = lin_coefs[, "Pr(>|t|)"],
  lci = lin_ci[, 1],
  uci = lin_ci[, 2]
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(coef, se, t, lci, uci),
    decimals = 6
  ) |>
  gt::fmt_number(columns = p_value, decimals = 6)
term coef se t p_value lci uci
(Intercept) 518.294002 348.088374 1.488972 0.137552 −166.728779 1,203.316783
hhsize 825.812641 55.068628 14.996064 0.000000 717.439977 934.185306
rooms 1,972.989241 143.261134 13.771978 0.000000 1,691.057561 2,254.920921
urbrurUrban 4,783.297099 270.632985 17.674479 0.000000 4,250.703154 5,315.891043

Logistic Regression

svy Results

logit_model = glm_sample.glm.fit(
    y="is_poor",
    x=["hhsize", "rooms", svy.Cat("urbrur")],
    family="binomial",
    link="logit",
    tol=1e-12,  # match the tightened R tolerance below
)

cols = [
    "term",
    "estimate",
    "std_err",
    "statistic",
    "p_value",
    "conf_low",
    "conf_high",
]
(
    GT(logit_model.to_polars().select(cols))
    .fmt_number(
        columns=["estimate", "std_err", "statistic", "conf_low", "conf_high"],
        decimals=6,
    )
    .fmt_number(columns="p_value", decimals=6)
)
term estimate std_err statistic p_value conf_low conf_high
_intercept_ −2.455070 0.213741 −11.486215 0.000000 −2.875702 −2.034438
hhsize 0.730864 0.039989 18.276645 0.000000 0.652167 0.809560
rooms −0.624728 0.060689 −10.293859 0.000000 −0.744162 −0.505294
urbrur_Urban −2.144876 0.142109 −15.093162 0.000000 −2.424541 −1.865211

R Results

logit_model_r <- svyglm(
  is_poor ~ hhsize + rooms + urbrur,
  design = design_glm,
  family = quasibinomial(),
  # Tighten R's IRLS tolerance (default epsilon = 1e-8 stops one
  # iteration short: coefficients agree to ~1e-10 regardless, but the
  # sandwich SEs inherit the last IRLS state and show ~1e-5 residual
  # differences in the 6th decimal). svy computes the variance at the
  # fully converged coefficients.
  control = list(epsilon = 1e-12)
)

logit_coefs <- summary(logit_model_r)$coefficients
logit_ci <- confint(logit_model_r)

data.frame(
  term = rownames(logit_coefs),
  coef = logit_coefs[, "Estimate"],
  se = logit_coefs[, "Std. Error"],
  t = logit_coefs[, "t value"],
  p_value = logit_coefs[, "Pr(>|t|)"],
  lci = logit_ci[, 1],
  uci = logit_ci[, 2]
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(coef, se, t, lci, uci),
    decimals = 6
  ) |>
  gt::fmt_number(columns = p_value, decimals = 6)
term coef se t p_value lci uci
(Intercept) −2.455070 0.213741 −11.486215 0.000000 −2.875702 −2.034438
hhsize 0.730864 0.039989 18.276645 0.000000 0.652167 0.809560
rooms −0.624728 0.060689 −10.293859 0.000000 −0.744162 −0.505294
urbrurUrban −2.144876 0.142109 −15.093162 0.000000 −2.424541 −1.865211
NoteR uses quasibinomial() for survey GLMs

R’s svyglm() requires family = quasibinomial() rather than binomial() for survey logistic regression. This avoids the “non-integer successes” warning that arises because survey-weighted likelihoods produce non-integer effective counts. The coefficient estimates are identical; only the dispersion parameter handling differs.

Poisson Regression

svy Results

poisson_model = glm_sample.glm.fit(
    y="hhsize",
    x=["rooms", svy.Cat("urbrur")],
    family="poisson",
    link="log",
)

cols = [
    "term",
    "estimate",
    "std_err",
    "statistic",
    "p_value",
    "conf_low",
    "conf_high",
]
(
    GT(poisson_model.to_polars().select(cols))
    .fmt_number(
        columns=["estimate", "std_err", "statistic", "conf_low", "conf_high"],
        decimals=6,
    )
    .fmt_number(columns="p_value", decimals=6)
)
term estimate std_err statistic p_value conf_low conf_high
_intercept_ 1.376818 0.042373 32.493087 0.000000 1.293432 1.460204
rooms 0.035542 0.007652 4.644929 0.000005 0.020484 0.050600
urbrur_Urban −0.160554 0.050782 −3.161629 0.001730 −0.260490 −0.060619

R Results

poisson_model_r <- svyglm(
  hhsize ~ rooms + urbrur,
  design = design_glm,
  family = quasipoisson(),
  control = list(epsilon = 1e-12)  # see note at the logistic model
)

pois_coefs <- summary(poisson_model_r)$coefficients
pois_ci <- confint(poisson_model_r)

data.frame(
  term = rownames(pois_coefs),
  coef = pois_coefs[, "Estimate"],
  se = pois_coefs[, "Std. Error"],
  t = pois_coefs[, "t value"],
  p_value = pois_coefs[, "Pr(>|t|)"],
  lci = pois_ci[, 1],
  uci = pois_ci[, 2]
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(coef, se, t, lci, uci),
    decimals = 6
  ) |>
  gt::fmt_number(columns = p_value, decimals = 6)
term coef se t p_value lci uci
(Intercept) 1.376818 0.042373 32.493087 0.000000 1.293432 1.460204
rooms 0.035542 0.007652 4.644929 0.000005 0.020484 0.050600
urbrurUrban −0.160554 0.050782 −3.161629 0.001730 −0.260490 −0.060619

Numerical Agreement

The side-by-side tables above let you compare svy and R visually. A validation study, though, should quantify the match rather than ask the reader to eyeball two separate tables. Below we take one representative estimator from each family—reusing the exact Sample/Design objects already built above—and report the absolute difference between svy and R for both the point estimate and its standard error.

First, collect the svy estimates:

def _es(obj):
    row = obj.to_polars().row(0, named=True)
    return (float(row["est"]), float(row["se"]))


svy_checks = {}
svy_checks["Mean, stratified (Taylor)"] = _es(
    sample_str.estimation.mean("api00")
)
svy_checks["Ratio, stratified+clustered (Taylor)"] = _es(
    sample_str_clus.estimation.ratio(y="tot_exp", x="hhsize")
)
svy_checks["Ratio, BRR (32 reps)"] = _es(
    sample_brr.estimation.ratio(
        y="weight", x="height", method="replication"
    )
)
svy_checks["Ratio, jackknife (62 reps)"] = _es(
    sample_jkn.estimation.ratio(
        y="weight", x="height", method="replication"
    )
)
svy_checks["Mean, bootstrap (1000 reps)"] = _es(
    sample_bs.estimation.mean(
        y="birthwgt", method="replication", drop_nulls=True
    )
)
svy_checks["Mean, SDR (80 reps)"] = _es(
    sample_acs.estimation.mean(
        y="HINCP", method="replication", drop_nulls=True
    )
)

# Logistic GLM: reuse the fitted model, take the hhsize coefficient
_lg = (
    logit_model.to_polars()
    .filter(pl.col("term") == "hhsize")
    .row(0, named=True)
)
svy_checks["GLM logistic, hhsize coef"] = (
    float(_lg["estimate"]),
    float(_lg["std_err"]),
)

Next, the matching R estimates:

r_es <- function(est, se) c(as.numeric(est), as.numeric(se))
r_checks <- list()

mean_str_r <- design_str |>
  srvyr::summarize(e = srvyr::survey_mean(api00, vartype = "se"))
r_checks[["Mean, stratified (Taylor)"]] <- r_es(mean_str_r$e, mean_str_r$e_se)

ratio_sc_r <- design_str_clus |>
  srvyr::summarize(e = srvyr::survey_ratio(tot_exp, hhsize, vartype = "se"))
r_checks[["Ratio, stratified+clustered (Taylor)"]] <-
  r_es(ratio_sc_r$e, ratio_sc_r$e_se)

brr_r <- svyratio(~weight, ~height, design = design_brr)
r_checks[["Ratio, BRR (32 reps)"]] <- r_es(coef(brr_r), SE(brr_r))

jkn_r <- svyratio(~weight, ~height, design = design_jkn)
r_checks[["Ratio, jackknife (62 reps)"]] <- r_es(coef(jkn_r), SE(jkn_r))

boot_r <- svymean(~birthwgt, design = design_bs, na.rm = TRUE)
r_checks[["Mean, bootstrap (1000 reps)"]] <- r_es(coef(boot_r), SE(boot_r))

sdr_r <- svymean(~HINCP, design = design_sdr, na.rm = TRUE)
r_checks[["Mean, SDR (80 reps)"]] <- r_es(coef(sdr_r), SE(sdr_r))

logit_c <- summary(logit_model_r)$coefficients
r_checks[["GLM logistic, hhsize coef"]] <-
  r_es(logit_c["hhsize", "Estimate"], logit_c["hhsize", "Std. Error"])

Finally, join them and report the absolute differences:

r_checks = r.r_checks

rows = []
for name, (est_s, se_s) in svy_checks.items():
    est_r, se_r = float(r_checks[name][0]), float(r_checks[name][1])
    rows.append(
        {
            "Estimator": name,
            "svy": est_s,
            "R survey": est_r,
            "|Δ estimate|": abs(est_s - est_r),
            "|Δ SE|": abs(se_s - se_r),
        }
    )

agreement = pl.DataFrame(rows)
(
    GT(agreement)
    .tab_header(title="svy vs. R: estimate and standard-error agreement")
    .fmt_number(columns=["svy", "R survey"], decimals=6)
    .fmt_scientific(columns=["|Δ estimate|", "|Δ SE|"], decimals=2)
)
svy vs. R: estimate and standard-error agreement
Estimator svy R survey |Δ estimate| |Δ SE|
Mean, stratified (Taylor) 662.287363 662.287363 0.00 7.11 × 10−15
Ratio, stratified+clustered (Taylor) 2,992.110041 2,992.110041 0.00 0.00
Ratio, BRR (32 reps) 0.426812 0.426812 5.55 × 10−17 9.65 × 10−18
Ratio, jackknife (62 reps) 0.426812 0.426812 5.55 × 10−17 4.38 × 10−17
Mean, bootstrap (1000 reps) 3,355.452419 3,355.452419 0.00 0.00
Mean, SDR (80 reps) 111,770.504058 111,770.504058 0.00 0.00
GLM logistic, hhsize coef 0.730864 0.730864 6.66 × 10−16 1.80 × 10−12
NoteReading the difference table

The absolute differences sit at the level of floating-point round-off—typically below 1e-6, and often exactly zero. This is what “identical” means in practice: the two libraries evaluate the same estimating equations, so any residual gap reflects arithmetic ordering, not methodology. The replication rows (BRR, jackknife, bootstrap, SDR) match because both packages consume the same replicate-weight columns; the small SE conventions discussed above (degrees of freedom, centering) are aligned before the comparison.

Comparison Summary

Category Estimator Design / Method Match Notes
Taylor Mean Stratified ✅
Mean One-stage cluster ✅
Mean Two-stage cluster ✅ Ultimate cluster variance
Mean Stratified + clustered ✅
Proportion Stratified + clustered ✅ Logit-transformed CIs
Total Stratified + clustered ✅
Ratio Stratified + clustered ✅
Domain Mean By subgroup ✅
Ratio By subgroup ✅
Replication BRR 32 replicates ✅
Jackknife 62 replicates ✅ Requires df specification
Bootstrap 1000 replicates ✅ Requires rscales in R
SDR 80 replicates ✅ ACS replicate weights
Categorical Cross-tabulation Cell %, SEs, Rao-Scott χ² (F) ✅ Adjusted Pearson statistic
t-test One- and two-sample ✅
GLM Linear Gaussian (identity) ✅
Logistic Binomial (logit) ✅ R uses quasibinomial()
Poisson Poisson (log) ✅ R uses quasipoisson()

When the numbers disagree

A validation study is a snapshot. What matters over the life of a library is what happens when the numbers don’t agree—and one case in this note is a worked example.

Earlier versions of svy reported an un-centered standard error for cross-tabulation cells under units="percent" (and for count_total scaling), inflating cell SEs by a p-dependent amount. The cause is the one described in the cross-tabulation section above: a cell percentage is a ratio of two estimated totals, so its variance requires the centered (Hájek) linearization that R applies in svymean. Treating the denominator as a fixed constant does not. Proportions estimated through estimation.prop were never affected.

The discrepancy surfaced in exactly this kind of side-by-side comparison against an independent implementation, and was corrected in samplics-org/svy#92. This note pins svy>=0.20.1 so the tables above reflect the corrected behavior.

We document this in the body rather than burying it in a changelog because it is the more useful signal. Any statistical library will have defects; what separates one an organization can depend on is whether those defects are found, disclosed in full, and checked against a reference implementation rather than against its own assumptions.

Conclusion

This validation study demonstrates that svy reproduces the results of R’s survey package for a wide range of design-based estimators when equivalent survey designs are specified. The Numerical Agreement table makes this concrete: across Taylor, replication, and GLM estimators, svy and R agree on both estimates and standard errors to the sixth decimal—differences at the level of floating-point round-off.

The agreement observed across all tested cases confirms that svy implements standard survey-sampling methodology correctly, including Taylor linearization, ultimate cluster variance estimation, replication-based variance estimators (BRR, jackknife, bootstrap, SDR), categorical tests, and survey-weighted GLMs.

These results support the use of svy for production survey analysis workflows and provide a basis for further validation of advanced features.

TipCommunity signal

Help make svy the standard for survey analysis in Python

If rigorous, design-based survey inference in Python matters to you, starring the repository helps signal demand and prioritize validation and stability work.

→ Star svy on GitHub

Notes on methodology

  1. Two-stage variance: Both packages use ultimate cluster estimation
  2. Proportion CIs: Both use logit-transformed confidence intervals
  3. Stratification: svy accepts tuples for multiple variables; R requires interaction() or paste()

References

Back to top

References

Lumley, Thomas. 2010. Complex Surveys: A Guide to Analysis Using R. Hoboken, NJ: Wiley.
World Bank. 2023. “Synthetic Data for an Imaginary Country, Sample, 2023.” World Bank, Development Data Group. https://doi.org/10.48529/MC1F-QH23.

Footnotes

  1. All numerical comparisons use identical survey designs and variance estimators in both packages.↩︎

  2. U.S. Census Bureau. (2023). American Community Survey 1-Year Public Use Microdata Sample [Data set]. Retrieved from https://www.census.gov/programs-surveys/acs/microdata.html↩︎

  3. Earlier svy versions reported the un-centered standard error for units="percent" (and for count_total scaling), which inflated cell SEs—estimation.prop was unaffected. See When the numbers disagree below and samplics-org/svy#92.↩︎

Citation

BibTeX citation:
@online{diallo2026,
  author = {Diallo, Mamadou S.},
  title = {Python’s Svy Vs {R’s} Survey: {Identical} {Results} {Across}
    18 {Estimators}},
  date = {2026-01-10},
  url = {https://svylab.com/learn/notes/posts/svy-vs-r-comparison/},
  langid = {en}
}
For attribution, please cite this work as:
Diallo, Mamadou S. 2026. “Python’s Svy Vs R’s Survey: Identical Results Across 18 Estimators.” January 10, 2026. https://svylab.com/learn/notes/posts/svy-vs-r-comparison/.