Can Python Match R for Survey Statistics? A Validation Study

true
Validation
Survey Methods
Python
R
We tested Python’s svy against R’s survey package across means, totals, ratios, BRR, jackknife, bootstrap, and SDR. The results match to 6 decimal places.
Author

Mamadou S. Diallo, Ph.D.

Published

January 10, 2026

Modified

April 18, 2026

Keywords

Python survey analysis, svy Python library, R survey package alternative, complex survey statistics Python, weighted survey analysis Python, Taylor linearization, BRR variance estimation, jackknife variance estimation, bootstrap survey, design-based inference Python, survey sampling Python, replicate weights Python, stratified cluster sample Python, svy vs R survey, survey design Python, svy package

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 ✅

These results validate svy as a statistically equivalent alternative to R’s survey package for complex survey analysis by (Lumley (2010)).

Introduction

For decades, survey statisticians have relied on specialized software for design-based inference—tools like SAS, SPSS, Stata, and R’s survey package. These have long been the trusted workhorses for analyzing complex survey data with proper variance estimation. Now Python joins that group. With the svy library, Python users can perform rigorous, design-based analysis while staying within the modern data science ecosystem.

But can svy be trusted?

This note answers that question by comparing results from Python’s svy and R’s survey package. Using the same dataset and survey design specifications, we show that svy produces identical estimates and standard errors—validating its statistical rigor and establishing it as a reliable choice for complex survey analysis.

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.4.8'
import polars as pl
from great_tables import GT

import svy

print(f"svy version: {svy.__version__}")
svy version: 0.15.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: districts (dnum) → schools (snum)
# Note: ssu is specified for documentation, but variance uses ultimate cluster
design_clus2 = svy.Design(psu="dnum", 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

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,515.403784 9,717.854889
Urban 14,437.918429 326.402120 13,795.599356 15,080.237501

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"))
  ) |>
  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,515.403784 9,717.854889
Urban 14,437.918429 326.402120 13,795.599356 15,080.237501

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.960973 1,867.098757
Yes 3,960.323902 89.247824 3,784.695203 4,135.952601

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"))
  ) |>
  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.960973 1,867.098757
Yes 3,960.323902 89.247824 3,784.695203 4,135.952601

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() uses normal (z) approximation by default for replicate weight designs.

To match svy’s t-based CIs, use confint(..., df = degf(design)).

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

equivalently

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

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.165349 12.914649 17.482732
Rural Yes 29.695594 1.309476 27.129069 32.262119
Urban No 1.856347 0.374548 1.122246 2.590448
Urban Yes 53.249369 0.831273 51.620104 54.878634 
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

crosstab <- survey::svytable(
  ~ urbrur + electricity,
  design_str_clus,
  Ntotal = 100
)

test_result <- summary(crosstab)

test_result$table
       electricity
urbrur         No       Yes
  Rural 15.198691 29.695594
  Urban  1.856347 53.249369
test_result$statistic

    Pearson's X^2: Rao & Scott adjustment

data:  NextMethod()
F = 193.17, ndf = 1, ddf = 301, p-value < 2.2e-16

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.137542 −166.700941 1,203.288945
hhsize 825.812641 55.068628 14.996064 0.000000 717.444381 934.180902
rooms 1,972.989241 143.261134 13.771978 0.000000 1,691.069018 2,254.909464
urbrur_Urban 4,783.297099 270.632985 17.674479 0.000000 4,250.724798 5,315.869399

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",
)

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.875685 −2.034455
hhsize 0.730864 0.039989 18.276645 0.000000 0.652171 0.809557
rooms −0.624728 0.060689 −10.293860 0.000000 −0.744157 −0.505299
urbrur_Urban −2.144876 0.142109 −15.093162 0.000000 −2.424529 −1.865223

R Results

logit_model_r <- svyglm(
  is_poor ~ hhsize + rooms + urbrur,
  design = design_glm,
  family = quasibinomial()
)

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.486178 0.000000 −2.875704 −2.034437
hhsize 0.730864 0.039989 18.276552 0.000000 0.652167 0.809561
rooms −0.624728 0.060692 −10.293466 0.000000 −0.744166 −0.505289
urbrurUrban −2.144876 0.142109 −15.093180 0.000000 −2.424540 −1.865212
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.293434 1.460202
rooms 0.035542 0.007652 4.644929 0.000005 0.020484 0.050600
urbrur_Urban −0.160554 0.050782 −3.161629 0.001729 −0.260487 −0.060621

R Results

poisson_model_r <- svyglm(
  hhsize ~ rooms + urbrur,
  design = design_glm,
  family = quasipoisson()
)

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

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 ✅
GLM Linear Gaussian (identity) ✅
Logistic Binomial (logit) ✅ R uses quasibinomial()
Poisson Poisson (log) ✅ R uses quasipoisson()

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 agreement observed across all tested cases confirms that svy implements standard survey-sampling methodology correctly, including Taylor linearization, ultimate cluster variance estimation, and replication-based variance estimators.

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 were conducted using identical survey designs and variance estimators.↩︎

  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↩︎

Citation

BibTeX citation:
@online{s._diallo2026,
  author = {S. Diallo, Mamadou and S. Diallo, Mamadou},
  title = {Can {Python} {Match} {R} for {Survey} {Statistics?} {A}
    {Validation} {Study}},
  date = {2026-01-10},
  url = {https://svylab.com/learn/notes/posts/svy-vs-r-comparison/},
  langid = {en}
}
For attribution, please cite this work as:
S. Diallo, Mamadou, and Mamadou S. Diallo. 2026. “Can Python Match R for Survey Statistics? A Validation Study.” January 10, 2026. https://svylab.com/learn/notes/posts/svy-vs-r-comparison/.