How Fast Is svy? Performance Benchmarks Against R’s survey

Timings for means, totals, ratios, proportions, domain estimates and replicate-weight variance in Python’s svy 0.22.1, compared with svy 0.18.2 and R’s survey 4.5, at 50,000 to 1,000,000 rows.

Keywords

survey analysis performance Python, R survey package speed, complex survey analysis Python, design-based inference Python, Taylor linearization performance, replicate weights Python, survey benchmark Python, svy vs R survey, large survey data Python, polars survey analysis

Summary

TipAt a glance

On a stratified two-stage sample of 1,000,000 rows, here is how long svy 0.22.1 takes to produce an estimate, next to svy 0.18.2 and to R’s survey package.

	Seconds		Times faster
	svy	R survey	vs R survey	vs svy 0.18.2
Mean	0.027	0.359	13.5×	9.3×
Ratio	0.040	0.790	19.6×	6.1×
Proportion	0.087	1.743	20.0×	4.3×
Mean by domain (12 domains)	0.124	1.955	15.7×	10.4×
8 means, single batched call	0.047	0.565	12.0×	41.0×
svy 0.22.1 at 1,000,000 rows. The batched row is compared against a loop of single-variable calls, which is the fastest route to the same result in svy 0.18.2.

Estimating eight means in a single call is about 41 times quicker than the loop svy 0.18.2 required and domain estimation runs roughly 16 times faster than R’s survey package.¹

What these numbers measure

svy and R’s survey package are not doing identical work when you ask each for a mean. Both compute the same design-based estimate and the same standard error, but around that arithmetic each library performs its own input validation, type checking, copying, and object construction, and these differ substantially. The numbers below are therefore not a like-for-like measurement of estimation arithmetic.

What they do measure is the thing an analyst actually experiences: you ask for an estimate, and this is how long you wait. That includes everything each library does to answer the question, which is the honest unit of comparison for someone deciding what to run their work in.

Whether the two libraries agree on the answer is a separate question, and it has its own note. The validation study shows that svy reproduces R’s survey package to six decimal places across 37 estimators. Every timing below was also checked to confirm both libraries returned the same estimate and the same standard error before the timing was recorded, so nothing here is fast because it computed something different.

How the benchmark was run

Data. A synthetic stratified two-stage sample with 50 strata, 2,000 primary sampling units, unequal weights, unbalanced cluster sizes and non-zero intra-cluster correlation. Sizes of 50,000, 250,000, 1,000,000 rows.

Machine. Apple M1 Max, 10 cores, 34 GB RAM. Python 3.13.12, polars 1.43.2, R survey 4.5.

Protocol. Each estimate is run once to warm up, then 7 times; the median is reported. Data is already in memory, so no file reading is included. Building the survey design is timed separately from estimating, since you normally build a design once and estimate from it many times.

NoteThe timing loop

Every case below is a zero-argument callable handed to the same harness: one discarded warm-up call, then seven timed calls, median reported.

def time_it(fn):
    fn()  # one discarded warm-up call
    times = []
    for _ in range(7):
        t0 = time.perf_counter()
        fn()
        times.append(time.perf_counter() - t0)
    return statistics.median(times)

The R script mirrors this loop exactly, timing with proc.time(). The full harness is published with this note as bench.zip: the data generator, both benchmark runners, and the run_all.sh driver that installs each svy version in its own environment. With uv and an R installation carrying survey, jsonlite and readr, unzipping it and running ./run_all.sh regenerates every number on this page — on your hardware, with your timings. The synthetic data is rebuilt from a fixed seed, so the inputs are identical to the ones used here; only the machine changes.

Results come from one machine and one design shape. Designs with very many small strata, deeply nested domains, or wide categorical tabulations will behave differently, so treat these as an indication of scale rather than a guarantee for every workload. Core count matters too, and not only for the absolute timings: svy parallelizes across cores through polars while R’s survey package is essentially single-threaded, so the ratios between the two shrink on machines with fewer cores and grow on machines with more. To see what your own hardware does, bench.zip reruns the whole comparison.

The data and the design

The generated sample has the design columns stratum, psu, and wgt, eight numeric outcomes y0–y7, a ratio denominator x, a binary outcome b, and a 12-level domain variable grp. Both libraries read the same generated file and declare the same design.

svy

import polars as pl
import svy

data = pl.read_parquet("bench_1000000.parquet")

design = svy.Design(stratum="stratum", psu="psu", wgt="wgt")
sample = svy.Sample(data, design)

R survey

library(survey)

data <- readr::read_csv("bench_1000000.csv")

design <- svydesign(ids = ~psu, strata = ~stratum, weights = ~wgt, data = data)

Every estimation benchmark reuses the sample and design objects above, so the timings reflect the estimator alone. Construction is timed as its own case, next.

Estimation timings

Building the design

The first case times the full path from raw data to a first estimate — design construction included.

svy

svy.Sample(data, design).estimation.mean("y0")

R survey

svymean(~y0, svydesign(ids = ~psu, strata = ~stratum, weights = ~wgt, data = data))

	50,000 rows			250,000 rows			1,000,000 rows
	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey
Build design and estimate a mean	0.017	0.010	0.139	0.070	0.036	0.709	0.273	0.133	3.150
Median seconds; lower is better.

Constructing the design object is a fixed cost you pay once. It is kept out of the remaining rows so they reflect what a repeated call actually costs.

One estimate at a time

The bread-and-butter estimators, each on a single variable against the design built above.

svy

sample.estimation.mean("y0")
sample.estimation.total("y0")
sample.estimation.ratio("y0", "x")
sample.estimation.prop("b", ci_method="logit")

R survey

svymean(~y0, design)
svytotal(~y0, design)
svyratio(~y0, ~x, design)
svyciprop(~I(b == 1), design, method = "logit")

The proportion row uses svyciprop() rather than svymean(~b) because svy’s prop builds a logit-transformed confidence interval; svyciprop(method = "logit") is the call that computes the same quantity.

	50,000 rows			250,000 rows			1,000,000 rows
	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey
Mean	0.013	0.002	0.023	0.060	0.008	0.094	0.248	0.027	0.359
Total	0.013	0.002	0.022	0.059	0.008	0.095	0.235	0.029	0.361
Ratio	0.014	0.003	0.043	0.062	0.010	0.188	0.247	0.040	0.790
Proportion	0.020	0.006	0.094	0.092	0.022	0.477	0.375	0.087	1.743
Median seconds; lower is better.

Domain estimation

Estimating by subgroup is common in survey work and often the slowest part of a tabulation run. Here each library estimates the mean of y0 within the 12 levels of grp.

svy

sample.estimation.mean("y0", by="grp")

R survey

svyby(~y0, ~grp, design, svymean)

	50,000 rows			250,000 rows			1,000,000 rows
	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey
Mean by domain (12 domains)	0.070	0.009	0.127	0.324	0.031	0.463	1.294	0.124	1.955
Median seconds; lower is better.

This is where the current release improves most on its predecessor, and where the gap to R’s survey package is widest.

Eight means: loop versus one batched call

svy accepts a list of columns and evaluates them together, sharing one pass over the design. R’s svymean takes a multi-variable formula natively, which is the fair counterpart.

svy

y_vars = [f"y{j}" for j in range(8)]  # y0 ... y7

# One call per variable
for v in y_vars:
    sample.estimation.mean(v)

# Single batched call
sample.estimation.mean(y_vars)

R survey

# One call per variable
for (v in paste0("y", 0:7)) svymean(as.formula(paste0("~", v)), design)

# Single batched call
svymean(~ y0 + y1 + y2 + y3 + y4 + y5 + y6 + y7, design)

	50,000 rows			250,000 rows			1,000,000 rows
	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey	svy 0.18.2	svy 0.22.1	R survey
8 means, one call each	0.106	0.016	0.172	0.476	0.059	0.722	1.926	0.187	2.822
8 means, single batched call	—	0.005	0.031	—	0.014	0.131	—	0.047	0.565
Median seconds; lower is better. An em dash marks an estimate the version cannot express: passing several variables in one call was introduced after 0.18.2.

The batched call is several times faster than calling the estimator once per variable, and it is the form worth reaching for whenever you need a block of indicators from the same sample.

Replicate-weight variance

Replication is where survey computation gets expensive, because the work grows with the number of replicate weights, and published microdata routinely ships several hundred of them. The harness generates bootstrap replicate weight columns (bsrw1, bsrw2, …) once, and both libraries consume the same columns — so neither is spending time generating weights, only computing a replicate variance from weights it was given.

The code below is shown for 400 replicates; the 100-replicate runs differ only in the count.

svy

rep_wgts = svy.RepWeights(method="bootstrap", prefix="bsrw", n_reps=400)
design_rep = svy.Design(wgt="wgt", rep_wgts=rep_wgts)
sample_rep = svy.Sample(data, design_rep)

sample_rep.estimation.mean("y0", method="replication")
sample_rep.estimation.total("y0", method="replication")
sample_rep.estimation.ratio("y0", "x", method="replication")

R survey

design_rep <- svrepdesign(
  data = data,
  weights = ~wgt,
  repweights = "^bsrw[0-9]+$",
  type = "bootstrap",
  replicates = 400,
  combined.weights = TRUE,
  rscales = rep(399 / 400, 400),
  # degf is supplied so both runtimes use B - 1 degrees of freedom;
  # deriving it from the replicate-weight matrix is not part of what is
  # timed. On this data the derived value is also B - 1.
  degf = 399
)

svymean(~y0, design_rep)
svytotal(~y0, design_rep)
svyratio(~y0, ~x, design_rep)

	100 replicate weights			400 replicate weights
	svy 0.22.1	R survey	times faster	svy 0.22.1	R survey	times faster
Build design and estimate a mean	0.010	0.125	12.7×	0.081	0.305	3.8×
Mean	0.003	0.022	6.9×	0.008	0.080	9.8×
Total	0.003	0.011	3.7×	0.008	0.041	4.9×
Ratio	0.003	0.047	14.9×	0.009	0.156	17.9×
8 means, one call each	0.024	0.194	8.2×	0.068	0.739	10.9×
Median seconds at 50,000 rows, bootstrap replicate weights supplied to both libraries.

svy 0.18.2 is not shown here. Several replication estimators were corrected between that release and this one, so the two versions no longer compute the same variance, and comparing their timings would not mean anything.

In short

For everyday survey estimation, svy 0.22.1 returns results in a fraction of a second on samples of a million rows, and it is consistently quicker than R’s survey package on the same designs.² The gap is widest for domain estimation and for estimating many variables in one call, which together cover most of what a tabulation run actually does.

If you are still on svy 0.18.2, upgrading is the single easiest speed improvement available, and it requires no change to your code unless you want to adopt the batched call.

---

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

References

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

Stay in the loop

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

Email Organization (optional)

No spam. Unsubscribe anytime.

Footnotes

1. For the estimators timed here: means, totals, ratios, proportions, domain means and a batched multi-variable call. Other estimators were not benchmarked. Measured on Apple M1 Max, 10 cores, 34 GB RAM, with polars using every core; the ratios shrink on machines with fewer. See How the benchmark was run.

2. The caveat on the summary applies here too: a handful of estimators, timed on one machine.