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

svy reproduces R’s survey package to six decimals across 37 estimators — Taylor, BRR, jackknife, bootstrap, SDR, quantiles, singleton PSUs, correlation, covariance, cross-tabs, t-tests, contrasts, marginal effects, and survey GLMs.

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, design-based correlation Python, svyvar Python equivalent, svyquantile Python equivalent, survey quantiles Python, lonely PSU Python, survey.lonely.psu equivalent, svycontrast Python equivalent, survey odds ratios Python, marginal effects survey Python, marginaleffects Python equivalent, 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
• Quantiles with Woodruff intervals
• Singleton (lonely) PSU handling
• Correlation and covariance
• Categorical data analysis and regression
• Marginal effects and post-estimation linear contrasts

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	✅
Correlation	Fisher-z CIs	✅
Covariance	Stratified + clustered	✅
Quantile	p = 0.01, 0.25, 0.5, 0.75, 0.99	✅
Domain estimation	Mean, proportion, total, ratio, quantile	✅
Singleton PSUs	lonely.psu remove, certainty, average, adjust	✅
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, probit, Poisson	✅
GLM	Gamma, inverse Gaussian, rate model	✅
Odds ratios	exp(β) with exp(CI)	✅
Margins	Average marginal effects, predictive	✅
Contrast	Between domains, between coefficients	✅

We don’t just assert the match — we quantify it: across the 37 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 37 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
• Quantiles, from the extreme tails to the median
• Correlation and covariance
• Domain (subpopulation) estimation for each of the above
• Singleton PSU handling, against every survey.lonely.psu option
• 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 across the exponential families, with offsets and non-canonical links
• Odds and rate ratios, marginal effects, and predictive margins
• Linear contrasts over domain estimates and model coefficients

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

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

Correlation

The R blocks below need the design as a base survey object rather than an srvyr one, because neither statistic has an srvyr verb.

wb_str_clus <- wb_synth_smp |>
  dplyr::mutate(stratum = paste(geo1, urbrur, sep = "_"))

design_str_clus_base <- svydesign(
  id = ~ea, strata = ~stratum, weights = ~hhweight,
  data = wb_str_clus, nest = TRUE
)

df_str_clus <- survey::degf(design_str_clus_base)
t_str_clus <- qt(0.975, df_str_clus)

svy Results

exp_size_corr = sample_str_clus.estimation.corr(("tot_exp", "hhsize"))

cols = ["est", "se", "lci", "uci"]
(
    GT(exp_size_corr.to_polars().select(cols))
    .fmt_number(columns=cols, decimals=6)
)

est	se	lci	uci
0.238605	0.029082	0.180606	0.294950

R Results

# A design-based correlation is a smooth function of five estimated means, so
# svycontrast() applies its own delta method over them. R has no correlation
# interval, so the Fisher-z transform svy uses by default is built by hand.
moments_r <- svymean(
  ~ tot_exp + hhsize + I(tot_exp^2) + I(hhsize^2) + I(tot_exp * hhsize),
  design_str_clus_base
)

corr_r <- svycontrast(moments_r, quote(
  (`I(tot_exp * hhsize)` - tot_exp * hhsize) /
    sqrt((`I(tot_exp^2)` - tot_exp^2) * (`I(hhsize^2)` - hhsize^2))
))

corr_est <- as.numeric(coef(corr_r))
corr_se <- as.numeric(SE(corr_r))
corr_z <- atanh(corr_est)
corr_z_se <- corr_se / (1 - corr_est^2)

dplyr::tibble(
  est = corr_est,
  se = corr_se,
  lci = tanh(corr_z - t_str_clus * corr_z_se),
  uci = tanh(corr_z + t_str_clus * corr_z_se)
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )

est	se	lci	uci
0.238605	0.029082	0.180606	0.294950

Covariance

svy Results

exp_size_cov = sample_str_clus.estimation.cov(("tot_exp", "hhsize"))

cols = ["est", "se", "lci", "uci"]
(
    GT(exp_size_cov.to_polars().select(cols))
    .fmt_number(columns=cols, decimals=6)
)

est	se	lci	uci
3,723.494515	561.238547	2,619.046350	4,827.942680

R Results

# svyvar() returns the full 2x2 matrix flattened; element 2 is the covariance.
var_r <- svyvar(~ tot_exp + hhsize, design_str_clus_base)

cov_est <- as.numeric(coef(var_r))[2]
cov_se <- as.numeric(SE(var_r))[2]

dplyr::tibble(
  est = cov_est,
  se = cov_se,
  lci = cov_est - t_str_clus * cov_se,
  uci = cov_est + t_str_clus * cov_se
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )

est	se	lci	uci
3,723.494515	561.238547	2,619.046350	4,827.942680

NoteR reports neither interval, and has no correlation estimator at all

svyvar() gives the covariance and its standard error but no confidence interval — confint() on it would default to df = Inf. Correlation has no survey function whatsoever: the estimate above comes from svycontrast() running its own delta method over five estimated means, so the agreement is an independent check of svy’s linearization rather than a restatement of the same formula.

Both intervals are therefore constructed explicitly in the R blocks, matching svy’s defaults: Wald on the design df for the covariance, which is unbounded, and Fisher’s z for the correlation, which is confined to [-1, 1] and would otherwise spill past the bound. corr(..., ci_method="wald") gives the symmetric interval if you want it.

svy’s covariance carries the same n / (n - 1) factor as svyvar(), so the two are on one scale. Naming a self-pair — cov(("tot_exp", "tot_exp")) — returns the variance on svyvar’s diagonal.

Quantiles

A quantile is not a smooth function of estimated totals, so the linearization behind every estimator above does not apply. Both libraries instead follow Woodruff (1952): estimate the weighted CDF, take the design-based variance of the proportion \(P(Y \le \hat{q})\) on the probability scale, and invert the CDF at \(F(\hat{q}) \pm t \cdot se_p\) to get the interval. The probabilities below run from the sparse tails, where implementations usually diverge, through the quartiles to the median.

svy Results

probs = [0.01, 0.25, 0.5, 0.75, 0.99]
exp_quant = sample_str_clus.estimation.quantile("tot_exp", p=probs)

cols = ["est", "se", "lci", "uci"]
(
    GT(exp_quant.to_polars().select(["prob"] + cols))
    .fmt_number(columns=cols, decimals=6)
)

prob	est	se	lci	uci
0.01	2,472.000000	180.397497	2,159.000000	2,869.000000
0.25	6,989.000000	172.520986	6,654.000000	7,333.000000
0.5	10,268.000000	222.066778	9,851.000000	10,725.000000
0.75	14,851.000000	315.568579	14,277.000000	15,519.000000
0.99	36,951.000000	1,063.582987	35,563.000000	39,749.000000

R Results

# Base survey rather than srvyr, so the defaults are visible: qrule = "math"
# is svy's q_method="higher", and interval.type = "mean" is the Woodruff
# interval on the design df.
probs <- c(0.01, 0.25, 0.5, 0.75, 0.99)
quant_r <- svyquantile(
  ~tot_exp, design_str_clus_base, quantiles = probs, ci = TRUE
)

dplyr::tibble(
  prob = probs,
  est = as.numeric(coef(quant_r)),
  se = as.numeric(SE(quant_r)),
  lci = confint(quant_r)[, 1],
  uci = confint(quant_r)[, 2]
) |>
  gt::gt() |>
  gt::fmt_number(columns = c(est, se, lci, uci), decimals = 6)

prob	est	se	lci	uci
0.01	2,472.000000	180.397497	2,159.000000	2,869.000000
0.25	6,989.000000	172.520986	6,654.000000	7,333.000000
0.50	10,268.000000	222.066778	9,851.000000	10,725.000000
0.75	14,851.000000	315.568579	14,277.000000	15,519.000000
0.99	36,951.000000	1,063.582987	35,563.000000	39,749.000000

Interpolation rules

A weighted CDF is a step function, so a rule is needed for where the quantile sits when \(p\) falls between two steps. The default above takes the next observation up. The other rule the two libraries share interpolates linearly between the neighbouring steps, which places the quantile and its limits between data values rather than on them.

svy Results

exp_quant_lin = sample_str_clus.estimation.quantile(
    "tot_exp", p=probs, q_method="linear"
)

(
    GT(exp_quant_lin.to_polars().select(["prob"] + cols))
    .fmt_number(columns=cols, decimals=6)
)

prob	est	se	lci	uci
0.01	2,464.139365	176.806789	2,151.578513	2,847.446373
0.25	6,988.302752	172.279860	6,651.949013	7,330.000000
0.5	10,265.851696	220.542292	9,844.000000	10,712.000000
0.75	14,850.237444	312.113324	14,275.000000	15,503.400970
0.99	36,920.027245	1,033.584426	35,525.907606	39,593.840677

R Results

quant_lin_r <- svyquantile(
  ~tot_exp, design_str_clus_base, quantiles = probs, ci = TRUE, qrule = "hf4"
)

dplyr::tibble(
  prob = probs,
  est = as.numeric(coef(quant_lin_r)),
  se = as.numeric(SE(quant_lin_r)),
  lci = confint(quant_lin_r)[, 1],
  uci = confint(quant_lin_r)[, 2]
) |>
  gt::gt() |>
  gt::fmt_number(columns = c(est, se, lci, uci), decimals = 6)

prob	est	se	lci	uci
0.01	2,464.139365	176.806789	2,151.578513	2,847.446373
0.25	6,988.302752	172.279860	6,651.949013	7,330.000000
0.50	10,265.851696	220.542292	9,844.000000	10,712.000000
0.75	14,850.237444	312.113324	14,275.000000	15,503.400970
0.99	36,920.027245	1,033.584426	35,525.907606	39,593.840677

NoteWhich q_method is which qrule

Under the default rule a quantile and its Woodruff limits are observed data values, so there is no rounding to compare: svy and R return the same numbers bit for bit, at every probability, including 0.01 and 0.99 where the CDF is inverted over a handful of observations. The standard error is the back-solved half-width \((uci - lci) / 2t\), so it agrees to the precision of the t quantile itself. Under the linear rule the values are interpolated, and the two agree to floating-point round-off.

svy q_method	R qrule	Rule
"higher" (default)	"math" (default), "hf1"	Next observation at or above \(p\)
"linear"	"hf4"	Linear interpolation of the weighted CDF
"lower", "nearest", "middle"	none	Previous observation, closest step, or their midpoint

R’s "school", "hf2", and "hf3" differ from "hf1" only when the CDF lands exactly on \(p\), which never happens with continuous weights, so they reproduce the default table here too. "hf5" to "hf9" and "shahvaish" are other interpolations with no svy counterpart.

The interval is centered at \(F(\hat{q})\) rather than at \(p\); on a step CDF those differ, and R’s retired oldsvyquantile(interval.type = "Wald") gives other limits for that reason.

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

Electricity access by urban and rural areas

svy Results

electricity_dom = sample_str_clus.estimation.prop("electricity", by="urbrur")

cols = ["est", "se", "lci", "uci"]
(
    GT(electricity_dom.to_polars().select(["urbrur", "electricity"] + cols))
    .fmt_number(columns=cols, decimals=6)
)

urbrur	electricity	est	se	lci	uci
Rural	No	0.338544	0.024732	0.291472	0.389044
Rural	Yes	0.661456	0.024732	0.610956	0.708528
Urban	No	0.033687	0.006756	0.022618	0.049896
Urban	Yes	0.966313	0.006756	0.950104	0.977382

R Results

# srvyr's grouped survey_prop() builds the logit interval on the full-design df
# even when `df` is passed, so go through base survey instead: svyciprop() on
# each domain subset, with that subset's own df.
prop_domain_r <- dplyr::bind_rows(lapply(c("Rural", "Urban"), function(area) {
  sub <- subset(design_str_clus_base, urbrur == area)
  dplyr::bind_rows(lapply(c("No", "Yes"), function(level) {
    ci <- svyciprop(
      as.formula(sprintf('~I(electricity == "%s")', level)),
      sub,
      method = "logit",
      df = survey::degf(sub)
    )
    dplyr::tibble(
      urbrur = area,
      electricity = level,
      est = as.numeric(coef(ci)),
      se = as.numeric(SE(ci)),
      lci = confint(ci)[1],
      uci = confint(ci)[2]
    )
  }))
}))

prop_domain_r |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )

urbrur	electricity	est	se	lci	uci
Rural	No	0.338544	0.024732	0.291472	0.389044
Rural	Yes	0.661456	0.024732	0.610956	0.708528
Urban	No	0.033687	0.006756	0.022618	0.049896
Urban	Yes	0.966313	0.006756	0.950104	0.977382

Households without electricity by urban and rural areas

svy Results

no_electricity_dom = sample_str_clus.estimation.total(
    "no_electricity", by="urbrur"
)

cols = ["est", "se", "lci", "uci"]
(
    GT(no_electricity_dom.to_polars().select(["urbrur"] + cols))
    .fmt_number(columns=cols, decimals=6)
)

urbrur	est	se	lci	uci
Rural	380,234.000686	29,154.166533	322,568.188595	437,899.812778
Urban	46,441.251273	9,370.282331	27,942.578658	64,939.923889

R Results

design_str_clus |>
  dplyr::mutate(no_electricity = electricity != "Yes") |>
  dplyr::group_by(urbrur) |>
  dplyr::summarize(
    est = srvyr::survey_total(
      no_electricity,
      vartype = c("se", "ci"),
      # Domain df, as above.
      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	380,234.000686	29,154.166533	322,568.188594	437,899.812778
Urban	46,441.251273	9,370.282331	27,942.578658	64,939.923889

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

Expenditure quantiles by urban and rural areas

svy Results

exp_quant_dom = sample_str_clus.estimation.quantile(
    "tot_exp", p=probs, by="urbrur"
)

cols = ["est", "se", "lci", "uci"]
(
    GT(
        exp_quant_dom.to_polars()
        .sort(["urbrur", "prob"])
        .select(["urbrur", "prob"] + cols)
    )
    .fmt_number(columns=cols, decimals=6)
)

urbrur	prob	est	se	lci	uci
Rural	0.01	2,159.000000	190.347509	1,730.000000	2,483.000000
Rural	0.25	5,865.000000	175.938734	5,484.000000	6,180.000000
Rural	0.5	8,027.000000	256.324534	7,532.000000	8,546.000000
Rural	0.75	10,982.000000	367.550171	10,309.000000	11,763.000000
Rural	0.99	28,207.000000	1,923.950721	25,360.000000	32,971.000000
Urban	0.01	3,042.000000	180.074292	2,675.000000	3,386.000000
Urban	0.25	8,985.000000	271.757686	8,414.000000	9,487.000000
Urban	0.5	12,672.000000	279.102489	12,108.000000	13,210.000000
Urban	0.75	17,891.000000	477.158883	17,079.000000	18,963.000000
Urban	0.99	40,692.000000	664.071439	39,312.000000	41,934.000000

R Results

# A subset() design carries the domain's own df, so svyquantile() needs no
# `df` argument here, unlike the grouped srvyr verbs above.
quant_domain_r <- dplyr::bind_rows(lapply(c("Rural", "Urban"), function(area) {
  sub <- subset(design_str_clus_base, urbrur == area)
  q <- svyquantile(~tot_exp, sub, quantiles = probs, ci = TRUE)
  dplyr::tibble(
    urbrur = area,
    prob = probs,
    est = as.numeric(coef(q)),
    se = as.numeric(SE(q)),
    lci = confint(q)[, 1],
    uci = confint(q)[, 2]
  )
}))

quant_domain_r |>
  gt::gt() |>
  gt::fmt_number(columns = c(est, se, lci, uci), decimals = 6)

urbrur	prob	est	se	lci	uci
Rural	0.01	2,159.000000	190.347509	1,730.000000	2,483.000000
Rural	0.25	5,865.000000	175.938734	5,484.000000	6,180.000000
Rural	0.50	8,027.000000	256.324534	7,532.000000	8,546.000000
Rural	0.75	10,982.000000	367.550171	10,309.000000	11,763.000000
Rural	0.99	28,207.000000	1,923.950721	25,360.000000	32,971.000000
Urban	0.01	3,042.000000	180.074292	2,675.000000	3,386.000000
Urban	0.25	8,985.000000	271.757686	8,414.000000	9,487.000000
Urban	0.50	12,672.000000	279.102489	12,108.000000	13,210.000000
Urban	0.75	17,891.000000	477.158883	17,079.000000	18,963.000000
Urban	0.99	40,692.000000	664.071439	39,312.000000	41,934.000000

ImportantDegrees of freedom for domain estimates

Every R block above asks for the domain’s own degrees of freedom. Without that 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_total(), and 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.
• Grouped survey_prop() goes further: it ignores df when building the logit interval and uses the full-design df regardless, which is why the proportion block drops to base survey and calls svyciprop() on each domain subset.
• svyquantile() on a subset() design is the exception: it takes that subset’s df on its own, so the quantile block passes nothing.

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.

Singleton PSUs

A stratum with a single PSU offers no between-PSU variability to measure, so its variance contribution is undefined. Both libraries refuse by default: R’s survey.lonely.psu = "fail" and svy’s SingletonError stop the estimation until a treatment is chosen. To compare the treatments, the World Bank sample is thinned with one deterministic rule in both languages: the three strata with the fewest enumeration areas keep only their lowest-numbered one.

svy Results

few_ea = (
    wb_synth_smp.group_by(["geo1", "urbrur"])
    .agg(pl.col("ea").n_unique().alias("n_ea"), pl.col("ea").min().alias("keep"))
    .sort(["n_ea", "geo1", "urbrur"])
    .head(3)
)
wb_single = (
    wb_synth_smp.join(
        few_ea.select(["geo1", "urbrur", "keep"]), on=["geo1", "urbrur"], how="left"
    )
    .filter(pl.col("keep").is_null() | (pl.col("ea") == pl.col("keep")))
    .drop("keep")
)

sample_single = svy.Sample(data=wb_single, design=design_str_clus)
sample_single.singleton.show()

shape: (3, 5)

singleton_key	n_obs	psu	geo1	urbrur
str	i64	str	str	str
"geo_02__by__Urban"	25	"23010"	"geo_02"	"Urban"
"geo_05__by__Urban"	25	"52035"	"geo_05"	"Urban"
"geo_08__by__Urban"	25	"83094"	"geo_08"	"Urban"

try:
    sample_single.estimation.mean("tot_exp")
except Exception as e:
    print(type(e).__name__)

SingletonError

Each treatment returns a new Sample; the estimator call does not change.

treatments = {
    "skip()": sample_single.singleton.skip(),
    "scale()": sample_single.singleton.scale(),
    "center()": sample_single.singleton.center(),
}

rows = []
for name, smp in treatments.items():
    est_row = smp.estimation.mean("tot_exp").to_polars().row(0, named=True)
    rows.append({"svy": name} | {k: est_row[k] for k in ["est", "se", "lci", "uci", "df"]})
single_svy = pl.DataFrame(rows)

GT(single_svy).fmt_number(columns=cols, decimals=6)

svy	est	se	lci	uci	df
skip()	11,876.948066	232.289128	11,419.741712	12,334.154419	287
scale()	11,876.948066	253.131208	11,378.718993	12,375.177139	287
center()	11,876.948066	234.420703	11,415.546208	12,338.349924	287

R Results

few_ea <- wb_synth_smp |>
  dplyr::group_by(geo1, urbrur) |>
  dplyr::summarize(n_ea = dplyr::n_distinct(ea), keep = min(ea), .groups = "drop") |>
  dplyr::arrange(n_ea, geo1, urbrur) |>
  head(3)

wb_single <- wb_synth_smp |>
  dplyr::left_join(
    dplyr::select(few_ea, geo1, urbrur, keep), by = c("geo1", "urbrur")
  ) |>
  dplyr::filter(is.na(keep) | ea == keep) |>
  dplyr::select(-keep) |>
  dplyr::mutate(stratum = paste(geo1, urbrur, sep = "_"))

design_single <- svydesign(
  id = ~ea, strata = ~stratum, weights = ~hhweight, data = wb_single, nest = TRUE
)

options(survey.lonely.psu = "fail")
tryCatch(svymean(~tot_exp, design_single), error = function(e) conditionMessage(e))

[1] "Stratum (geo_02_Urban) has only one PSU at stage 1"

lonely_mean <- function(option) {
  options(survey.lonely.psu = option)
  m <- svymean(~tot_exp, design_single)
  ci <- confint(m, df = degf(design_single))
  dplyr::tibble(
    est = as.numeric(coef(m)), se = as.numeric(SE(m)),
    lci = ci[1], uci = ci[2], df = degf(design_single)
  )
}

single_r <- dplyr::bind_rows(
  remove = lonely_mean("remove"),
  certainty = lonely_mean("certainty"),
  average = lonely_mean("average"),
  adjust = lonely_mean("adjust"),
  .id = "lonely.psu"
)
options(survey.lonely.psu = "fail")

single_r |>
  gt::gt() |>
  gt::fmt_number(columns = c(est, se, lci, uci), decimals = 6)

lonely.psu	est	se	lci	uci	df
remove	11,876.948066	232.289128	11,419.741712	12,334.154419	287
certainty	11,876.948066	232.289128	11,419.741712	12,334.154419	287
average	11,876.948066	253.131208	11,378.718993	12,375.177139	287
adjust	11,876.948066	234.420703	11,415.546208	12,338.349924	287

The point estimate is the same on every row of both tables: no treatment drops a record, they only decide what the lonely stratum contributes to the variance. Pairing each R option with its svy counterpart:

r_single = {
    name: (float(est), float(se))
    for name, est, se in zip(r.single_r["lonely.psu"], r.single_r["est"], r.single_r["se"])
}
svy_single = {row["svy"]: row for row in single_svy.iter_rows(named=True)}

pairs = [
    ("remove", "skip()"),
    ("certainty", "skip()"),
    ("average", "scale()"),
    ("adjust", "center()"),
]
(
    GT(
        pl.DataFrame(
            [
                {
                    "R lonely.psu": r_opt,
                    "svy": s_opt,
                    "svy se": svy_single[s_opt]["se"],
                    "R se": r_single[r_opt][1],
                    "|Δ estimate|": abs(svy_single[s_opt]["est"] - r_single[r_opt][0]),
                    "|Δ SE|": abs(svy_single[s_opt]["se"] - r_single[r_opt][1]),
                }
                for r_opt, s_opt in pairs
            ]
        )
    )
    .fmt_number(columns=["svy se", "R se"], decimals=6)
    .fmt_scientific(columns=["|Δ estimate|", "|Δ SE|"], decimals=2)
)

R lonely.psu	svy	svy se	R se	|Δ estimate|	|Δ SE|
remove	skip()	232.289128	232.289128	0.00	0.00
certainty	skip()	232.289128	232.289128	0.00	0.00
average	scale()	253.131208	253.131208	0.00	2.84 × 10−14
adjust	center()	234.420703	234.420703	0.00	0.00

NoteThe same word, two different treatments

R’s "certainty" and "remove" give the same numbers: in survey:::onestrat both let the lonely stratum contribute nothing to a Taylor variance. Their svy counterpart is skip().

svy’s certainty() is a different operation: it makes every unit in the lonely stratum its own PSU, so the stratum keeps a variance contribution. R has no option for that, nor for pool() and collapse(), which merge the lonely strata into a pseudo-stratum or a neighbour. R reproduces each of those only by recoding the id or stratum variable before svydesign(), and then agrees exactly, because the two designs are the same design.

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

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

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),
        # Log household size, used as the exposure offset of the rate model
        "log_hhsize": svy.col("hhsize").log(),
        # Expenditure in thousands: the inverse-Gaussian variance goes as
        # mu^3, so the raw scale overflows before IRLS can converge
        "exp_k": svy.col("tot_exp") / 1000.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),
    log_hhsize = log(hhsize),
    exp_k = tot_exp / 1000
  ) |>
  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, and quasibinomial(link = "probit") for the probit model below. 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. Gamma() and inverse.gaussian() need no quasi-variant — the warning is specific to the binomial and Poisson likelihoods.

Probit Regression

svy Results

probit_model = glm_sample.glm.fit(
    y="is_poor",
    x=["hhsize", "rooms", svy.Cat("urbrur")],
    family="binomial",
    link="probit",
    tol=1e-12,
)

cols = [
    "term",
    "estimate",
    "std_err",
    "statistic",
    "p_value",
    "conf_low",
    "conf_high",
]
(
    GT(probit_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.453682	0.137932	−10.539085	0.000000	−1.725127	−1.182237
hhsize	0.394777	0.025427	15.525738	0.000000	0.344737	0.444817
rooms	−0.296267	0.062206	−4.762685	0.000003	−0.418685	−0.173848
urbrur_Urban	−1.178814	0.077744	−15.162731	0.000000	−1.331811	−1.025817

R Results

probit_model_r <- svyglm(
  is_poor ~ hhsize + rooms + urbrur,
  design = design_glm,
  family = quasibinomial(link = "probit"),
  # Probit needs a tighter stop than the logit above: at epsilon = 1e-12 R's
  # IRLS is still moving in the 6th decimal of the sandwich SE. At 1e-16 it
  # settles on the same point svy converges to.
  control = list(epsilon = 1e-16, maxit = 200)
)

probit_coefs <- summary(probit_model_r)$coefficients
probit_ci <- confint(probit_model_r)

data.frame(
  term = rownames(probit_coefs),
  coef = probit_coefs[, "Estimate"],
  se = probit_coefs[, "Std. Error"],
  t = probit_coefs[, "t value"],
  p_value = probit_coefs[, "Pr(>|t|)"],
  lci = probit_ci[, 1],
  uci = probit_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.453682	0.137932	−10.539083	0.000000	−1.725127	−1.182237
hhsize	0.394777	0.025427	15.525735	0.000000	0.344737	0.444817
rooms	−0.296267	0.062206	−4.762685	0.000003	−0.418685	−0.173848
urbrurUrban	−1.178814	0.077744	−15.162732	0.000000	−1.331811	−1.025817

Odds Ratios

to_polars(exponentiate=True) reports exp(β) with the exponentiated link-scale interval — an odds ratio under logit, a rate ratio under log, a hazard ratio under cloglog — and names the column for the link. The same result in R is exp(coef(f)) with exp(confint(f)).

svy Results

cols = ["term", "odds_ratio", "std_err", "conf_low", "conf_high"]
(
    GT(logit_model.fitted.to_polars(exponentiate=True).select(cols))
    .fmt_number(
        columns=["odds_ratio", "std_err", "conf_low", "conf_high"],
        decimals=6,
    )
)

term	odds_ratio	std_err	conf_low	conf_high
_intercept_	0.085857	0.213741	0.056377	0.130754
hhsize	2.076874	0.039989	1.919697	2.246920
rooms	0.535407	0.060689	0.475132	0.603328
urbrur_Urban	0.117083	0.142109	0.088519	0.154863

R Results

data.frame(
  term = rownames(logit_coefs),
  odds_ratio = exp(logit_coefs[, "Estimate"]),
  se = logit_coefs[, "Std. Error"],
  lci = exp(logit_ci[, 1]),
  uci = exp(logit_ci[, 2])
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(odds_ratio, se, lci, uci),
    decimals = 6
  )

term	odds_ratio	se	lci	uci
(Intercept)	0.085857	0.213741	0.056377	0.130754
hhsize	2.076874	0.039989	1.919697	2.246920
rooms	0.535407	0.060689	0.475132	0.603328
urbrurUrban	0.117083	0.142109	0.088519	0.154863

NoteThe standard error stays on the link scale

Both libraries exponentiate the estimate and both interval bounds, but leave std_err — and the statistic and p-value with it — on the link scale, where the Wald test is computed. The resulting interval is therefore not symmetric about the ratio, and there is deliberately no “standard error of the odds ratio” column: a symmetric error bar around a ratio is the mistake this layout is meant to prevent.

svy refuses exponentiate=True on the links where exp(β) is not a ratio (identity, probit, inverse, inverse_squared) rather than printing a number with no interpretation.

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

Poisson Rate Model

An offset enters the linear predictor with its coefficient fixed at 1. With log exposure as the offset, a Poisson model’s coefficients are log rate ratios rather than log counts — here, deaths per person-year rather than deaths per household.

svy Results

rate_model = glm_sample.glm.fit(
    y="deaths_12m",
    x=["rooms", svy.Cat("urbrur")],
    family="poisson",
    offset="log_hhsize",
    tol=1e-12,
)

cols = [
    "term",
    "estimate",
    "std_err",
    "statistic",
    "p_value",
    "conf_low",
    "conf_high",
]
(
    GT(rate_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_	−4.626185	0.120970	−38.242277	0.000000	−4.864250	−4.388121
rooms	0.055654	0.022458	2.478137	0.013761	0.011458	0.099851
urbrur_Urban	0.932823	0.141609	6.587304	0.000000	0.654142	1.211504

R Results

rate_model_r <- svyglm(
  deaths_12m ~ rooms + urbrur + offset(log_hhsize),
  design = design_glm,
  family = quasipoisson(),
  control = list(epsilon = 1e-12)
)

rate_coefs <- summary(rate_model_r)$coefficients
rate_ci <- confint(rate_model_r)

data.frame(
  term = rownames(rate_coefs),
  coef = rate_coefs[, "Estimate"],
  se = rate_coefs[, "Std. Error"],
  t = rate_coefs[, "t value"],
  p_value = rate_coefs[, "Pr(>|t|)"],
  lci = rate_ci[, 1],
  uci = rate_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)	−4.626185	0.120970	−38.242277	0.000000	−4.864250	−4.388121
rooms	0.055654	0.022458	2.478137	0.013761	0.011458	0.099851
urbrurUrban	0.932823	0.141609	6.587304	0.000000	0.654142	1.211504

Notepredict() and the offset

The two libraries part company on prediction, not on fitting. R’s own two answers disagree with each other: on the same rows fitted(f) returns exp(Xβ + offset), while predict(f, newdata = ...) drops the offset and returns exp(Xβ) — a rate, not a count. svy follows fitted(), which is also what Stata’s predict after glm, exposure() does, and raises if newdata arrives without the offset column rather than silently predicting a rate.

Gamma Regression

Expenditure is positive and right-skewed, which is the shape a Gamma model with a log link is for. Note that R takes Gamma() here, not a quasi-family — the non-integer-successes warning is specific to the binomial and Poisson likelihoods.

svy Results

gamma_model = glm_sample.glm.fit(
    y="tot_exp",
    x=["hhsize", "rooms", svy.Cat("urbrur")],
    family="gamma",
    link="log",
    tol=1e-12,
)

cols = [
    "term",
    "estimate",
    "std_err",
    "statistic",
    "p_value",
    "conf_low",
    "conf_high",
]
(
    GT(gamma_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_	8.337323	0.034906	238.852278	0.000000	8.268630	8.406016
hhsize	0.074083	0.005082	14.578301	0.000000	0.064082	0.084084
rooms	0.165183	0.005784	28.557342	0.000000	0.153800	0.176566
urbrur_Urban	0.398991	0.024234	16.463802	0.000000	0.351298	0.446683

R Results

gamma_model_r <- svyglm(
  tot_exp ~ hhsize + rooms + urbrur,
  design = design_glm,
  family = Gamma(link = "log"),
  control = list(epsilon = 1e-12)
)

gamma_coefs <- summary(gamma_model_r)$coefficients
gamma_ci <- confint(gamma_model_r)

data.frame(
  term = rownames(gamma_coefs),
  coef = gamma_coefs[, "Estimate"],
  se = gamma_coefs[, "Std. Error"],
  t = gamma_coefs[, "t value"],
  p_value = gamma_coefs[, "Pr(>|t|)"],
  lci = gamma_ci[, 1],
  uci = gamma_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)	8.337323	0.034906	238.852278	0.000000	8.268630	8.406016
hhsize	0.074083	0.005082	14.578301	0.000000	0.064082	0.084084
rooms	0.165183	0.005784	28.557342	0.000000	0.153800	0.176566
urbrurUrban	0.398991	0.024234	16.463802	0.000000	0.351298	0.446683

Inverse Gaussian Regression

The last of the exponential families, fitted here on expenditure in thousands: the inverse-Gaussian variance goes as mu^3, so the raw scale overflows before IRLS can converge — in both libraries. The canonical inverse_squared link is unusable at this scale for the same reason, so the fit uses a log link.

svy Results

invgauss_model = glm_sample.glm.fit(
    y="exp_k",
    x=["hhsize", "rooms", svy.Cat("urbrur")],
    family="inverse_gaussian",
    link="log",
    tol=1e-12,
)

cols = [
    "term",
    "estimate",
    "std_err",
    "statistic",
    "p_value",
    "conf_low",
    "conf_high",
]
(
    GT(invgauss_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.309056	0.039373	33.247552	0.000000	1.231571	1.386540
hhsize	0.093496	0.006274	14.903200	0.000000	0.081150	0.105842
rooms	0.183885	0.006259	29.378393	0.000000	0.171567	0.196203
urbrur_Urban	0.393878	0.025014	15.746008	0.000000	0.344651	0.443105

R Results

invgauss_model_r <- svyglm(
  exp_k ~ hhsize + rooms + urbrur,
  design = design_glm,
  family = inverse.gaussian(link = "log"),
  control = list(epsilon = 1e-12)
)

ig_coefs <- summary(invgauss_model_r)$coefficients
ig_ci <- confint(invgauss_model_r)

data.frame(
  term = rownames(ig_coefs),
  coef = ig_coefs[, "Estimate"],
  se = ig_coefs[, "Std. Error"],
  t = ig_coefs[, "t value"],
  p_value = ig_coefs[, "Pr(>|t|)"],
  lci = ig_ci[, 1],
  uci = ig_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.309056	0.039373	33.247550	0.000000	1.231571	1.386540
hhsize	0.093496	0.006274	14.903199	0.000000	0.081150	0.105842
rooms	0.183885	0.006259	29.378391	0.000000	0.171567	0.196203
urbrurUrban	0.393878	0.025014	15.746007	0.000000	0.344651	0.443105

Marginal Effects

Coefficients on a non-identity link are not effects on the response scale. margins() moves them there: an average marginal effect is the mean over the sample of dmu/dx, and a predictive margin is the mean fitted response with one predictor fixed at a chosen value. R has no equivalent in survey; the reference here is the marginaleffects package, whose avg_slopes() and avg_predictions() accept a svyglm fit directly.

Average marginal effects

A continuous predictor gets a derivative. A categorical one gets a discrete contrast instead — the mean change in fitted response when every household is moved from the reference level to level k — reported one row per non-reference level.

svy Results

ame = pl.concat(
    [m.to_polars() for m in logit_model.margins()], how="diagonal_relaxed"
)

cols = ["margin", "se", "lci", "uci"]
(
    GT(ame.select(["term", "value"] + cols))
    .fmt_number(columns=cols, decimals=6)
    .sub_missing(missing_text="—")
)

term	value	margin	se	lci	uci
hhsize	—	0.066884	0.002588	0.061791	0.071977
rooms	—	−0.057171	0.005184	−0.067373	−0.046969
urbrur	Urban - Rural	−0.203761	0.014832	−0.232950	−0.174573

R Results

library(marginaleffects)

# marginaleffects defaults to a normal (z) interval; pass the model's residual
# df for the t interval svy reports.
avg_slopes(
  logit_model_r,
  wts = weights(design_glm, "sampling"),
  df = logit_model_r$df.residual
) |>
  as.data.frame() |>
  dplyr::transmute(
    term,
    value = contrast,
    margin = estimate,
    se = std.error,
    lci = conf.low,
    uci = conf.high
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(margin, se, lci, uci),
    decimals = 6
  )

term	value	margin	se	lci	uci
hhsize	dY/dX	0.066884	0.002588	0.061791	0.071977
rooms	dY/dX	−0.057171	0.005184	−0.067373	−0.046969
urbrur	Urban - Rural	−0.203761	0.014832	−0.232950	−0.174573

Predictive margins

svy Results

pred_margins = logit_model.margins(at={"rooms": [1, 3, 5]})

cols = ["margin", "se", "lci", "uci"]
(
    GT(pred_margins.to_polars().select(["term", "value"] + cols))
    .fmt_number(columns=cols, decimals=6)
)

term	value	margin	se	lci	uci
rooms	1	0.296965	0.012542	0.272282	0.321648
rooms	3	0.164715	0.007718	0.149526	0.179903
rooms	5	0.081562	0.009326	0.063209	0.099916

R Results

avg_predictions(
  logit_model_r,
  variables = list(rooms = c(1, 3, 5)),
  wts = weights(design_glm, "sampling"),
  df = logit_model_r$df.residual
) |>
  as.data.frame() |>
  dplyr::transmute(
    term = "rooms",
    value = rooms,
    margin = estimate,
    se = std.error,
    lci = conf.low,
    uci = conf.high
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(margin, se, lci, uci),
    decimals = 6
  )

term	value	margin	se	lci	uci
rooms	1	0.296965	0.012542	0.272282	0.321648
rooms	3	0.164715	0.007718	0.149526	0.179903
rooms	5	0.081562	0.009326	0.063209	0.099916

NoteDegrees of freedom, again

marginaleffects builds a normal (z) interval by default, so its bounds are narrower than svy’s until it is passed df = f$df.residual. This is the same default that shows up in confint() and in srvyr’s grouped summaries: the estimate and the standard error agree before the df is set, and only the interval moves.

Linear Contrasts

Both libraries expose the same post-estimation step: given a set of estimates and their covariance matrix, a linear combination Lθ̂ has variance LVLᵀ and t-based inference on the design degrees of freedom. R does it with svycontrast(); svy exposes contrast() on both estimation results and fitted models, with the combination written as an expression over svy.estd(...) references.

Between domains

The difference in average expenditure between urban and rural households.

svy Results

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

urban_rural_gap = tot_exp_dom.contrast(
    {"Urban - Rural": svy.estd("Urban") - svy.estd("Rural")}
)

cols = ["est", "se", "lci", "uci", "t", "p_value"]
(
    GT(urban_rural_gap.to_polars().select(["contrast"] + cols))
    .fmt_number(columns=cols, decimals=6)
)

contrast	est	se	lci	uci	t	p_value
Urban - Rural	5,321.289092	447.080293	4,441.490276	6,201.087908	11.902312	0.000000

R Results

# svyby() must be asked for covmat = TRUE: without the between-domain
# covariance there is nothing for svycontrast() to combine.
tot_exp_dom_r <- svyby(
  ~tot_exp, ~urbrur, design_str_clus_base, svymean, covmat = TRUE
)

gap_r <- svycontrast(
  tot_exp_dom_r, list(`Urban - Rural` = c(Rural = -1, Urban = 1))
)

# confint() defaults to df = Inf, so pass the design df for a t interval.
gap_ci <- confint(gap_r, df = df_str_clus)
gap_t <- as.numeric(coef(gap_r)) / as.numeric(SE(gap_r))

dplyr::tibble(
  contrast = "Urban - Rural",
  est = as.numeric(coef(gap_r)),
  se = as.numeric(SE(gap_r)),
  lci = gap_ci[1],
  uci = gap_ci[2],
  t = gap_t,
  p_value = 2 * pt(-abs(gap_t), df_str_clus)
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )

contrast	est	se	lci	uci	t	p_value
Urban - Rural	5,321.289092	447.080293	4,441.490276	6,201.087908	11.902312	0.000000

Between model coefficients

Whether an extra room and an extra household member carry the same expenditure association, from the linear model fitted above.

svy Results

rooms_vs_hhsize = lin_model.contrast(
    {"rooms - hhsize": svy.estd("rooms") - svy.estd("hhsize")}
)

cols = ["est", "se", "lci", "uci", "t", "p_value"]
(
    GT(rooms_vs_hhsize.to_polars().select(["contrast"] + cols))
    .fmt_number(columns=cols, decimals=6)
)

contrast	est	se	lci	uci	t	p_value
rooms - hhsize	1,147.176600	167.053398	818.422778	1,475.930421	6.867125	0.000000

R Results

rooms_vs_hhsize_r <- svycontrast(
  lin_model_r, list(`rooms - hhsize` = c(hhsize = -1, rooms = 1))
)

# A survey GLM spends a df per estimated coefficient, so the residual df is
# degf(design) - p + 1 -- what svy reports on the contrast, and what confint()
# needs to be told.
glm_df <- lin_model_r$df.residual
ct_ci <- confint(rooms_vs_hhsize_r, df = glm_df)
ct_t <- as.numeric(coef(rooms_vs_hhsize_r)) / as.numeric(SE(rooms_vs_hhsize_r))

dplyr::tibble(
  contrast = "rooms - hhsize",
  est = as.numeric(coef(rooms_vs_hhsize_r)),
  se = as.numeric(SE(rooms_vs_hhsize_r)),
  lci = ct_ci[1],
  uci = ct_ci[2],
  t = ct_t,
  p_value = 2 * pt(-abs(ct_t), glm_df)
) |>
  gt::gt() |>
  gt::fmt_number(
    columns = where(is.numeric),
    decimals = 6
  )

contrast	est	se	lci	uci	t	p_value
rooms - hhsize	1,147.176600	167.053398	818.422778	1,475.930421	6.867125	0.000000

NoteDegrees of freedom for contrasts

A contrast draws on every estimate it touches, so it takes the full design df — not the domain-aware df each individual row carries. That is R’s degf convention, and svy follows it: the domain contrast above reports df = 301 even though the rural and urban rows report 133 and 168. On a fitted model the df drops to the residual degf(design) - p + 1.

svy applies both automatically. In R, confint() on a svycontrast result still defaults to df = Inf, so the df has to be passed each time.

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["Correlation, stratified+clustered (Taylor)"] = _es(
    sample_str_clus.estimation.corr(("tot_exp", "hhsize"))
)
svy_checks["Covariance, stratified+clustered (Taylor)"] = _es(
    sample_str_clus.estimation.cov(("tot_exp", "hhsize"))
)

# Domain estimate: pull the rural row out of the by-group result
_rural = (
    sample_str_clus.estimation.mean(y="tot_exp", by="urbrur")
    .to_polars()
    .filter(pl.col("urbrur") == "Rural")
    .row(0, named=True)
)
svy_checks["Mean, rural domain (Taylor)"] = (
    float(_rural["est"]),
    float(_rural["se"]),
)
svy_checks["Median, stratified+clustered (Taylor)"] = _es(
    sample_str_clus.estimation.median("tot_exp")
)
svy_checks["Quantile p = 0.01, linear rule (Taylor)"] = _es(
    sample_str_clus.estimation.quantile("tot_exp", p=0.01, q_method="linear")
)
svy_checks["Mean, lonely strata centered (Taylor)"] = _es(
    treatments["center()"].estimation.mean("tot_exp")
)

# Domain quantile: the sparse upper tail of the rural distribution
_q99_rural = (
    sample_str_clus.estimation.quantile("tot_exp", p=0.99, by="urbrur")
    .to_polars()
    .filter(pl.col("urbrur") == "Rural")
    .row(0, named=True)
)
svy_checks["Quantile p = 0.99, rural domain (Taylor)"] = (
    float(_q99_rural["est"]),
    float(_q99_rural["se"]),
)
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
    )
)

svy_checks["GLM gamma, hhsize coef"] = (
    float(
        gamma_model.to_polars()
        .filter(pl.col("term") == "hhsize")
        .row(0, named=True)["estimate"]
    ),
    float(
        gamma_model.to_polars()
        .filter(pl.col("term") == "hhsize")
        .row(0, named=True)["std_err"]
    ),
)

# Margins: the categorical AME row, urban vs rural
_ame_cat = (
    ame.filter(pl.col("term") == "urbrur").row(0, named=True)
)
svy_checks["AME, urbrur contrast"] = (
    float(_ame_cat["margin"]),
    float(_ame_cat["se"]),
)

svy_checks["Contrast, urban - rural domain mean"] = _es(urban_rural_gap)
svy_checks["Contrast, GLM rooms - hhsize"] = _es(rooms_vs_hhsize)

# 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)

r_checks[["Correlation, stratified+clustered (Taylor)"]] <-
  r_es(coef(corr_r), SE(corr_r))

r_checks[["Covariance, stratified+clustered (Taylor)"]] <-
  r_es(as.numeric(coef(var_r))[2], as.numeric(SE(var_r))[2])

mean_rural_r <- svymean(
  ~tot_exp, subset(design_str_clus_base, urbrur == "Rural")
)
r_checks[["Mean, rural domain (Taylor)"]] <-
  r_es(coef(mean_rural_r), SE(mean_rural_r))

r_checks[["Median, stratified+clustered (Taylor)"]] <-
  r_es(coef(quant_r)[3], SE(quant_r)[3])

r_checks[["Quantile p = 0.01, linear rule (Taylor)"]] <-
  r_es(coef(quant_lin_r)[1], SE(quant_lin_r)[1])

adjust_r <- single_r[single_r$lonely.psu == "adjust", ]
r_checks[["Mean, lonely strata centered (Taylor)"]] <- r_es(adjust_r$est, adjust_r$se)

q99_rural_r <- svyquantile(
  ~tot_exp, subset(design_str_clus_base, urbrur == "Rural"),
  quantiles = 0.99, ci = TRUE
)
r_checks[["Quantile p = 0.99, rural domain (Taylor)"]] <-
  r_es(coef(q99_rural_r), SE(q99_rural_r))

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

gamma_c <- summary(gamma_model_r)$coefficients
r_checks[["GLM gamma, hhsize coef"]] <-
  r_es(gamma_c["hhsize", "Estimate"], gamma_c["hhsize", "Std. Error"])

ame_r <- as.data.frame(avg_slopes(
  logit_model_r,
  wts = weights(design_glm, "sampling"),
  df = logit_model_r$df.residual
))
ame_cat <- ame_r[ame_r$term == "urbrur", ]
r_checks[["AME, urbrur contrast"]] <-
  r_es(ame_cat$estimate, ame_cat$std.error)

r_checks[["Contrast, urban - rural domain mean"]] <-
  r_es(coef(gap_r), SE(gap_r))

r_checks[["Contrast, GLM rooms - hhsize"]] <-
  r_es(coef(rooms_vs_hhsize_r), SE(rooms_vs_hhsize_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
Correlation, stratified+clustered (Taylor)	0.238605	0.238605	5.31 × 10−14	2.34 × 10−15
Covariance, stratified+clustered (Taylor)	3,723.494515	3,723.494515	1.36 × 10−12	5.68 × 10−13
Mean, rural domain (Taylor)	9,116.629337	9,116.629337	1.64 × 10−11	0.00
Median, stratified+clustered (Taylor)	10,268.000000	10,268.000000	0.00	3.35 × 10−11
Quantile p = 0.01, linear rule (Taylor)	2,464.139365	2,464.139365	6.37 × 10−12	1.18 × 10−11
Mean, lonely strata centered (Taylor)	11,876.948066	11,876.948066	0.00	0.00
Quantile p = 0.99, rural domain (Taylor)	28,207.000000	28,207.000000	0.00	1.72 × 10−8
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 gamma, hhsize coef	0.074083	0.074083	5.45 × 10−15	1.08 × 10−15
AME, urbrur contrast	−0.203761	−0.203761	2.53 × 10−15	3.93 × 10−9
Contrast, urban - rural domain mean	5,321.289092	5,321.289092	8.55 × 10−11	0.00
Contrast, GLM rooms - hhsize	1,147.176600	1,147.176600	1.84 × 10−11	2.84 × 10−14
GLM logistic, hhsize coef	0.730864	0.730864	1.11 × 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	✅
	Correlation	Stratified + clustered	✅	R has no correlation estimator; svycontrast delta method
	Covariance	Stratified + clustered	✅	R’s svyvar scale, n/(n-1) factor
	Quantile	p = 0.01, 0.25, 0.5, 0.75, 0.99	✅	Woodruff interval; higher is qrule = "math", linear is "hf4"
Domain	Mean	By subgroup	✅	Requires domain df in R
	Proportion	By subgroup	✅	Requires base survey::svyciprop on the subset
	Total	By subgroup	✅	Requires domain df in R
	Ratio	By subgroup	✅	Requires domain df in R
	Quantile	By subgroup	✅	subset() design carries the domain df
Singletons	skip()	lonely.psu = "remove", "certainty"	✅	Lonely stratum contributes nothing
	scale()	lonely.psu = "average"	✅	Full sample; domain fraction differs, see note
	center()	lonely.psu = "adjust"	✅	Grand-mean centering
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()
	Probit	Binomial (probit)	✅	R needs epsilon = 1e-16
	Poisson	Poisson (log)	✅	R uses quasipoisson()
	Rate model	Poisson + log-exposure offset	✅	Coefficients are log rate ratios
	Gamma	Gamma (log)	✅	No quasi-family needed
	Inverse Gaussian	Inverse Gaussian (log)	✅	Response rescaled; mu^3 variance
	Odds ratios	exp(β), exp(CI)	✅	SE stays on the link scale
Margins	Average marginal effect	Continuous predictor	✅	vs marginaleffects::avg_slopes()
	Average marginal effect	Categorical contrast	✅	One row per non-reference level
	Predictive margin	Response at fixed values	✅	vs marginaleffects::avg_predictions()
Contrast	Between domains	svyby(covmat=TRUE)	✅	Full-design df
	Between coefficients	Linear GLM	✅	Residual df degf - p + 1

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, quantile, singleton, replication, association, contrast, margin, 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, Woodruff quantile intervals, domain estimation, lonely-PSU treatments, replication-based variance estimators (BRR, jackknife, bootstrap, SDR), design-based correlation and covariance, categorical tests, survey-weighted GLMs across the exponential families, marginal effects, and post-estimation contrasts.

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()
4. Domain df: svy applies the domain’s own degrees of freedom automatically; in R it must be requested per group
5. Quantiles: both invert the weighted CDF at \(F(\hat{q}) \pm t \cdot se_p\) (Woodruff); svy’s higher and linear rules are R’s qrule = "math" and "hf4"; the other three svy rules have no R counterpart
6. Association: R has no correlation estimator, so the reference values come from svycontrast()’s independent delta method
7. Contrasts: both take the full-design df (the residual df on a fitted model); R’s svyby() needs covmat = TRUE for the between-domain covariance
8. Margins: R’s survey has none, so the reference is marginaleffects, which needs wts = for the sampling weights and df = for a t interval
9. Singletons: svy’s skip(), scale(), and center() are R’s remove/certainty, average, and adjust; svy’s certainty() is a PSU recode with no R option

References

• svy Documentation
• R survey package
• Stata Survey Data Reference Manual

Stay in the loop

New releases, tutorials, and svyLab updates — a few emails a year.

Email Organization (optional)

No spam. Unsubscribe anytime.

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.