Svy

svy Skill

svy: design-based analysis of complex survey data in Python. Covers survey design specification (strata, PSU, weights, FPC), variance estimation (Taylor linearization, BRR, jackknife, bootstrap), descriptive estimation (means, totals, proportions, ratios, medians), survey-weighted GLM regression (gaussian, binomial, Poisson), domain/subpopulation analysis, calibration, and survey data I/O (SAS, SPSS, Stata). Uses Polars DataFrames natively. Use when analyzing data from complex sample surveys (NHANES, CPS, ACS PUMS, MEPS, ECLS-K, BRFSS, DHS). For non-survey regression, use statsmodels; for fixed effects, use pyfixest; for panel/IV models, use linearmodels.

Comprehensive skill for complex survey data analysis with svy. Use decision trees below to find the right guidance, then load detailed references.

What is svy?

svy is the Python package for design-based analysis of complex survey data:

• Survey-aware estimation: Means, totals, proportions, ratios, medians with proper design-based standard errors
• GLM regression: Survey-weighted linear, logistic, and Poisson regression with design-adjusted inference
• Flexible variance estimation: Taylor linearization (default), bootstrap, BRR (including Fay's modification), and jackknife (JK1, JKn) replicate methods
• Domain estimation: Correct subpopulation analysis without pre-filtering (preserves design structure)
• Native Polars: Built on Polars DataFrames, not pandas
• Survey data I/O: Read SAS (.sas7bdat), SPSS (.sav), Stata (.dta), and CSV with metadata
• Calibration: Post-stratification, raking, and GREG calibration for weight adjustment
• Validated: Results numerically equivalent to R's survey package across all methods

Version Notes

This skill targets svy 0.13.0 (released 2026-03-25). svy supersedes samplics (archived 2026-03-10), an earlier library by the same author (Mamadou S. Diallo, Ph.D.). Key differences from samplics:

• Unified Sample object replaces separate TaylorEstimator / ReplicateEstimator classes
• Polars-native (samplics used numpy arrays)
• Expanded GLM support and data I/O module (svy.io)
• The API is substantially different from samplics — do not assume samplics patterns carry over

How to Use This Skill

Reference File Structure

Reading Order

1. New to svy? Start with design-weights.md then estimation.md 2. Need survey-weighted regression? Read design-weights.md then regression.md 3. Have replicate weights already? Read design-weights.md (replicate design section) then estimation.md or regression.md 4. Setting up a federal survey (NHANES, CPS, etc.)? Read design-weights.md (federal survey patterns table) 5. Coming from samplics? Read design-weights.md for the new API; the Sample object replaces TaylorEstimator/ReplicateEstimator

Related Skills

Quick Decision Trees

"I need to analyze survey data"

What task?
├─ Descriptive statistics (mean, total, proportion)
│   └─ ./references/estimation.md
├─ Regression model
│   ├─ Linear (continuous outcome) → ./references/regression.md
│   ├─ Logistic (binary outcome) → ./references/regression.md
│   └─ Poisson (count outcome) → ./references/regression.md
├─ Set up the survey design object
│   └─ ./references/design-weights.md
├─ Read survey data from SAS/SPSS/Stata
│   └─ ./references/design-weights.md
├─ Subpopulation / domain analysis
│   └─ ./references/estimation.md
└─ Cross-tabulation
    └─ ./references/estimation.md

"I need survey-weighted regression"

What model?
├─ Linear regression (continuous Y)
│   └─ family="gaussian" → ./references/regression.md
├─ Logistic regression (binary Y)
│   └─ family="binomial" → ./references/regression.md
├─ Poisson regression (count Y)
│   └─ family="poisson" → ./references/regression.md
├─ Ordinal logistic / Cox survival / IV
│   └─ Not in svy — use rpy2 + R survey package (see rpy2 bridge below)
└─ Fixed effects + survey weights
    └─ Not directly supported — see Boundaries below

"I need to set up variance estimation"

What do you have?
├─ Design variables (strata, PSU, weights)
│   └─ Taylor linearization → ./references/design-weights.md
├─ Pre-computed replicate weights
│   ├─ BRR weights → ./references/design-weights.md
│   ├─ Jackknife weights → ./references/design-weights.md
│   └─ Bootstrap weights → ./references/design-weights.md
├─ Need to create replicate weights from design
│   └─ ./references/design-weights.md
└─ Not sure what I have
    └─ Read survey documentation first → ./references/design-weights.md (federal survey table)

"I need descriptive statistics from a survey"

What statistic?
├─ Population mean → ./references/estimation.md
├─ Population total → ./references/estimation.md
├─ Proportion → ./references/estimation.md
├─ Ratio (Y/X) → ./references/estimation.md
├─ Median / quantile → ./references/estimation.md
├─ Cross-tabulation → ./references/estimation.md
├─ By subgroup (domain estimation) → ./references/estimation.md
└─ Hypothesis test (t-test) → ./references/estimation.md

Boundaries

svy covers:

• Design-based estimation (descriptive and regression) for complex surveys
• Taylor and replicate-weight variance estimation
• Domain/subpopulation analysis
• Calibration and weight adjustment
• Survey data I/O

svy does NOT cover (use other tools):

• Fixed effects models — use pyfixest (survey weights + FE is methodologically complex; consult data-scientist skill)
• Panel data models (RE, FD, between) — use linearmodels
• Difference-in-differences — use pyfixest
• Causal inference methods (IV, RD, synthetic control) — use pyfixest/linearmodels/statsmodels
• Time series analysis — use statsmodels
• Machine learning — use scikit-learn
• Ordinal logistic, Cox proportional hazards, negative binomial — use rpy2 + R survey package
• Survey sampling design and sample size calculation — use data-scientist skill for methodology

The rpy2 Bridge

For models svy does not support (ordinal logistic, survival models, negative binomial GLM, cumulative link models), fall back to R's survey package via rpy2:

Decision rule: If the model family is not "gaussian", "binomial", or "poisson", use rpy2.

The R survey package (survey::svyglm, survey::svyolr, survey::svycoxph) covers the full range of survey-weighted models. Set up the survey design in R using the same design variables you would pass to svy.Design. See R survey package documentation at r-survey.r-forge.r-project.org for API details.

Legacy: samplics

samplics (2020-2026) is archived. svy supersedes it with a cleaner API, Polars integration, and expanded methods. If working with legacy code that uses samplics:

• The API is substantially different — TaylorEstimator/ReplicateEstimator classes are replaced by svy.Sample
• samplics used numpy arrays; svy uses Polars DataFrames
• Consult samplics documentation at samplics-org.github.io/samplics/ for legacy reference
• Migration requires rewriting, not find-and-replace

File-First Execution in Research Workflows

Important: In data research pipelines (see CLAUDE.md), svy analyses are executed through script files, not interactively. This ensures auditability and reproducibility.

The pattern: 1. Write estimation/regression code to scripts/stage8_analysis/{step}_{task-name}.py 2. Execute via Bash with automatic output capture wrapper script 3. Validation results get automatically embedded in scripts as comments 4. If failed, create versioned copy for fixes

Closely read agent_reference/SCRIPT_EXECUTION_REFERENCE.md for the mandatory file-first execution protocol. All survey analysis scripts must follow the Inline Audit Trail (IAT) standard — document design specification choices (why these strata/PSU/weights, what variance method, domain definitions) with # INTENT:, # REASONING:, and # ASSUMES: comments.

---

Quick Reference

Essential Import

import svy

Core Workflow

# 1. Load data
data = svy.io.read_stata("nhanes.dta")

# 2. Specify design
design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec2yr")

# 3. Create sample object
sample = svy.Sample(data=data, design=design)

# 4. Estimate
mean_bmi = sample.estimation.mean("bmxbmi")
model = sample.glm.fit(y="bmxbmi", x=["ridageyr", svy.Cat("riagendr")], family="gaussian")

Core Operations

Topic Index

Citation

When this library is used as a primary analytical tool, include in the report's Software & Tools references:

Diallo, M.S. svy: Python package for complex survey sampling and analysis [Computer software]. (Formerly samplics.)

Cite when: svy is used for survey-weighted estimation with complex survey designs (strata, PSU, replicate weights). Do not cite when: Only imported but no survey estimation performed.

For method-specific citations (e.g., variance estimation techniques), consult the reference files in this skill and agent_reference/CITATION_REFERENCE.md.

svy Design and Weights Reference

svy v0.13.0 — syntax and library guidance only.

---

Contents

1. Core Concepts 2. Creating a Taylor Linearization Design 3. Creating a Replicate Weight Design 4. Finite Population Correction (FPC) 5. Weight Types and Selection 6. Calibration and Post-Stratification 7. Singleton PSU Handling 8. Reading Survey Data Files 9. Federal Survey Design Quick-Reference 10. Combining Survey Cycles 11. Polars Integration Notes 12. Migration from samplics

---

Core Concepts

Every svy analysis starts with two objects:

1. `svy.Design` — describes the sampling structure (how units were selected) 2. `svy.Sample` — binds data to a design, enabling estimation

The design determines how variance is estimated. Two approaches:

Both produce valid design-based inference. Taylor linearization requires design variables (strata, PSU); replicate weights encode the design information within the weight columns themselves.

---

Creating a Taylor Linearization Design

Taylor linearization (also called the "ultimate cluster" method) is the default and most common approach. It requires knowing the stratification and primary sampling unit (PSU) variables.

Minimal Design (Weight Only)

import svy

# Simple random sample with unequal weights
design = svy.Design(wgt="weight")
sample = svy.Sample(data=data, design=design)

This assumes no stratification and no clustering — only unequal selection probabilities. Variance is computed using a with-replacement approximation.

Stratified Design (No Clustering)

# Stratified random sample
design = svy.Design(stratum="region", wgt="weight")
sample = svy.Sample(data=data, design=design)

Clustered Design (No Stratification)

# Cluster sample (e.g., schools as PSUs)
design = svy.Design(psu="school_id", wgt="weight")
sample = svy.Sample(data=data, design=design)

Full Complex Design (Stratified Clustered)

# Complex multi-stage design — most federal surveys
design = svy.Design(
    stratum="sdmvstra",
    psu="sdmvpsu",
    wgt="wtmec2yr"
)
sample = svy.Sample(data=data, design=design)

Multiple Stratification Variables

Some designs have nested stratification (e.g., region within urban/rural):

# Multiple strata specified as a tuple
design = svy.Design(
    stratum=("region_id", "urban_rural"),
    psu="psu_id",
    wgt="final_weight"
)
sample = svy.Sample(data=data, design=design)

svy.Design Parameters

---

Creating a Replicate Weight Design

When the data provider supplies replicate weights (common for public-use files from NCES, Census, BLS), use them instead of Taylor linearization. Replicate weights encode the design information and often provide more accurate variance estimates for complex statistics.

Bootstrap Replicate Weights

# Create a RepWeights object specifying the column prefix, count, and method
rep_wgts = svy.RepWeights(
    prefix="pwgtp",                       # Column name prefix (matches pwgtp1, pwgtp2, ...)
    n_reps=80,                            # Number of replicate weight columns
    method=svy.EstimationMethod.BOOTSTRAP # Replication method
)
design = svy.Design(
    wgt="pwgtp",                          # Main analysis weight
    rep_wgts=rep_wgts                     # RepWeights object
)
sample = svy.Sample(data=data, design=design)

BRR (Balanced Repeated Replication) Weights

rep_wgts = svy.RepWeights(
    prefix="brr_wt",                      # Column name prefix (matches brr_wt1, brr_wt2, ...)
    n_reps=64,                            # Number of replicate weight columns
    method=svy.EstimationMethod.BRR       # BRR method
)
design = svy.Design(wgt="finalwgt", rep_wgts=rep_wgts)
sample = svy.Sample(data=data, design=design)

Fay's BRR Modification

Fay's method is a variant of BRR that perturbs rather than deletes half-samples, reducing instability for small domains. The Fay coefficient (rho) is typically between 0.3 and 0.5.

rep_wgts = svy.RepWeights(
    prefix="brr_wt",
    n_reps=64,
    method=svy.EstimationMethod.BRR,
    fay_coef=0.5                          # Fay's rho (parameter name is fay_coef)
)
design = svy.Design(wgt="finalwgt", rep_wgts=rep_wgts)
sample = svy.Sample(data=data, design=design)

Jackknife Replicate Weights

rep_wgts = svy.RepWeights(
    prefix="jk_wt",
    n_reps=50,
    method=svy.EstimationMethod.JACKKNIFE # JKn (delete-one-PSU jackknife)
)
design = svy.Design(wgt="finalwgt", rep_wgts=rep_wgts)
sample = svy.Sample(data=data, design=design)

Jackknife methods:

• JK1: Delete-one jackknife (unstratified designs)
• JKn: Delete-one-PSU jackknife (stratified designs — most common)

Replicate Weight Column Specification

Replicate weight columns are specified through the svy.RepWeights object using a prefix and n_reps count. The prefix must match the column naming pattern in the data (e.g., prefix "pwgtp" matches columns pwgtp1, pwgtp2, ..., pwgtp80).

# Specify columns by prefix and count
rep_wgts = svy.RepWeights(prefix="pwgtp", n_reps=80, method=svy.EstimationMethod.BOOTSTRAP)
design = svy.Design(wgt="pwgtp", rep_wgts=rep_wgts)

Available `EstimationMethod` values: TAYLOR, BRR, BOOTSTRAP, JACKKNIFE, SDR

When to Use Replicate Weights vs. Taylor

---

Finite Population Correction (FPC)

The FPC adjusts variance estimates when a substantial fraction of the population is sampled (typically > 5-10%). Without FPC, variance estimates are conservative (too large).

# The Design parameter for FPC is `pop_size` (not `fpc`):
design = svy.Design(
    stratum="stratum",
    psu="psu_id",
    wgt="weight",
    pop_size="pop_size"                   # Column name containing population size
)
sample = svy.Sample(data=data, design=design)

Known issue (v0.13.0): The pop_size parameter is accepted by svy.Design as a column name string (type str | None), but Sample._calculate_fpc() performs isinstance(value, Number) on the raw string rather than resolving the column from the DataFrame, causing a TypeError. FPC is non-functional in 0.13.0. Workaround: Omit pop_size and note in your analysis that variance estimates are conservative (not FPC-adjusted). If FPC is critical, apply the finite population correction factor manually to variance estimates. Monitor future releases for a fix.

When to apply FPC:

• Sampling fraction > 5% of the population within strata
• Self-representing (certainty) strata where all units are selected
• Small population surveys (e.g., all schools in a small state)

When to skip FPC:

• Large national surveys sampling << 1% of the population (NHANES, CPS, etc.)
• When population size is unknown
• When a conservative variance estimate is acceptable

---

Weight Types and Selection

Common Weight Types in Federal Surveys

Choosing the Right Weight

1. Identify the analysis population: Which respondents have non-missing data for your variables? 2. Match the weight to the population: Use the weight designed for that subset 3. Consult the survey documentation: Weight variable names and usage instructions are survey-specific

Example — NHANES weight selection:

• Interview data only: use wtint2yr (interview weight)
• Examination data: use wtmec2yr (MEC exam weight)
• Fasting subsample: use wtsaf2yr (fasting subsample weight)
• Diet recall (day 1): use wtdrd1 (dietary day 1 weight)

Using the wrong weight produces biased estimates. The weight must correspond to the most restrictive component of data collection used in your analysis.

---

Calibration and Post-Stratification

svy supports weight adjustment methods to improve estimates by incorporating known population totals.

Post-Stratification

Post-stratification adjusts weights so that weighted sample totals match known population totals (e.g., Census counts by age/sex/race).

# Post-stratification via the weighting accessor
sample = sample.weighting.poststratify(
    controls={"18-34": 50_000_000, "35-64": 60_000_000, "65+": 40_000_000},
    by="age_group"
)

Raking (Iterative Proportional Fitting)

Raking adjusts weights to match marginal distributions of multiple variables simultaneously.

# Raking via the weighting accessor — controls is a dict of {variable: {level: target_total}}
sample = sample.weighting.rake(
    controls={
        "gender": {"Male": 160_000_000, "Female": 165_000_000},
        "age_group": {"18-34": 50_000_000, "35-64": 60_000_000, "65+": 40_000_000},
        "region": {"NE": 55_000_000, "MW": 65_000_000, "S": 75_000_000, "W": 80_000_000},
    }
)

GREG (Generalized Regression Estimator)

GREG calibration uses a regression model to adjust weights, incorporating both categorical and continuous auxiliary variables.

# GREG calibration via the weighting accessor
sample = sample.weighting.calibrate(
    controls={svy.Cat("gender"): {"Male": 160_000_000, "Female": 165_000_000}}
)

Note: Calibration is typically performed by the data producer before public release. Analysts working with public-use files usually do not need to calibrate — the provided weights already incorporate these adjustments. Only calibrate if you are working with raw sampling weights or need to adjust for a specific target population.

---

Singleton PSU Handling

A "singleton PSU" (or "lonely PSU") occurs when a stratum contains only one primary sampling unit. This makes within-stratum variance undefined, because variance estimation requires at least two PSUs per stratum.

Common Causes

• Rare subpopulations where domain estimation leaves some strata with only one PSU
• Data subsetting that removes PSUs from strata
• Design strata that genuinely have only one PSU (certainty selections)

Handling Options

Singleton handling is on the sample.singleton accessor, not on Design:

# Detect singleton strata
if sample.singleton.exists():
    print(sample.singleton.summary())

    # Handle singletons — choose one method:
    sample = sample.singleton.certainty()    # zero variance contribution (for true certainty strata)
    sample = sample.singleton.center()       # center at grand mean (equivalent to R's "adjust")
    sample = sample.singleton.combine()      # combine singleton strata
    sample = sample.singleton.collapse()     # collapse into other strata
    sample = sample.singleton.pool()         # pool with neighboring strata
    sample = sample.singleton.scale()        # scale variance contribution
    sample = sample.singleton.skip()         # skip singleton strata

`SingletonHandling` enum values: ERROR, CERTAINTY, SKIP, COMBINE, COLLAPSE, POOL, SCALE, CENTER

Best practice: The center() or pool() approach is generally safest. The certainty() approach (zero variance contribution) is appropriate only when the stratum truly is a certainty selection. Always report how singleton PSUs were handled.

Prevention

• Avoid domain estimation on very small subgroups
• If subsetting is necessary, consider collapsing strata before analysis
• Use replicate weights when available (the replication method handles singletons implicitly)

---

Reading Survey Data Files

svy provides svy.io methods for reading common survey data formats. These return Polars DataFrames.

Stata (.dta)

data = svy.io.read_stata("nhanes_2017_2020.dta")

Most federal surveys distribute data in Stata format. Value labels and variable labels may be preserved as metadata.

SAS (.sas7bdat)

data = svy.io.read_sas("meps_h233.sas7bdat")

MEPS and some NCES surveys use SAS format. Ensure the SAS format catalog (.sas7bcat) is in the same directory if needed for value labels.

SPSS (.sav)

data = svy.io.read_spss("ecls_k_2011.sav")

NCES surveys (ECLS-K, ELS, HSLS) often distribute data in SPSS format.

CSV with Metadata

data = svy.io.read_csv("survey_data.csv", metadata="codebook.json")

Parquet (Via Polars Directly)

svy does not have a dedicated parquet reader — use Polars directly:

import polars as pl
data = pl.read_parquet("data/raw/survey_data.parquet")

In DAAF research pipelines, data is typically stored as parquet after initial conversion. Use svy.io for the initial read from the original format, then save as parquet for subsequent use.

---

Federal Survey Design Quick-Reference

A quick-reference table for setting up svy designs for commonly used federal surveys. Always verify against the current survey documentation — design variable names can change across survey cycles.

NHANES Example

import svy
import polars as pl

data = pl.read_parquet("data/raw/nhanes_demo.parquet")

# INTENT: Standard NHANES complex design for MEC-examined participants
# REASONING: sdmvstra and sdmvpsu are masked design variables;
#   wtmec2yr is the 2-year MEC exam weight for exam-based analyses
# ASSUMES: All analysis variables were collected during the MEC examination
design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec2yr")
sample = svy.Sample(data=data, design=design)

ACS PUMS Example (Replicate Weights)

import svy
import polars as pl

data = pl.read_parquet("data/raw/acs_pums_2022.parquet")

# INTENT: ACS PUMS uses successive-difference replication (SDR)
# REASONING: No design variables in public-use file; must use replicate weights
# ASSUMES: Person-level analysis using person weight and person replicate weights
rep_wgts = svy.RepWeights(prefix="pwgtp", n_reps=80, method=svy.EstimationMethod.SDR)
design = svy.Design(wgt="pwgtp", rep_wgts=rep_wgts)
sample = svy.Sample(data=data, design=design)

MEPS Example

import svy

data = svy.io.read_sas("h233.sas7bdat")

# INTENT: MEPS household component, full-year consolidated file
# REASONING: varstr/varpsu are the standard MEPS design variables
# ASSUMES: Person-level analysis for FY 2021
design = svy.Design(stratum="varstr", psu="varpsu", wgt="perwt21f")
sample = svy.Sample(data=data, design=design)

---

Combining Survey Cycles

Some analyses require combining multiple cycles of a survey (e.g., NHANES 2017-2018 + 2019-2020) to increase sample size for rare subpopulations.

Weight Adjustment

When combining N two-year cycles, divide the survey weight by N:

import polars as pl

# Combine two NHANES cycles
cycle1 = pl.read_parquet("data/raw/nhanes_2017_2018.parquet")
cycle2 = pl.read_parquet("data/raw/nhanes_2019_2020.parquet")
combined = pl.concat([cycle1, cycle2])

# INTENT: Adjust weights for combined 4-year analysis
# REASONING: NHANES analytic guidelines require dividing 2-year weights by
#   the number of cycles combined to produce correct population estimates
# ASSUMES: Both cycles use the same design structure and weight definitions
combined = combined.with_columns(
    (pl.col("wtmec2yr") / 2).alias("wtmec4yr")
)

design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec4yr")
sample = svy.Sample(data=combined, design=design)

Important Caveats

• Only combine cycles with the same design structure
• Consult survey-specific guidelines for the correct weight adjustment
• For NHANES: the analytic guidelines at wwwn.cdc.gov/nchs/nhanes/tutorials/ provide detailed instructions
• For ACS: combining 1-year and 5-year estimates requires separate methodology
• The variance structure may change across cycles — check for redesign years

---

Polars Integration Notes

svy Expects Polars DataFrames

svy.Sample expects a Polars DataFrame as the data argument. Data loaded via svy.io methods returns Polars DataFrames automatically.

Converting from Other Formats

import polars as pl

# From pandas
import pandas as pd
pd_df = pd.read_csv("survey.csv")
pl_df = pl.from_pandas(pd_df)

# From numpy arrays
import numpy as np
arr = np.load("survey_data.npz")
pl_df = pl.DataFrame({"var1": arr["var1"], "var2": arr["var2"]})

# From parquet (native Polars)
pl_df = pl.read_parquet("data.parquet")

Column Type Requirements

• Weight column: Must be numeric (Float64 or Int64)
• Strata column: Can be string or integer; treated as categorical
• PSU column: Can be string or integer; treated as categorical
• Analysis variables: Numeric for estimation; string/categorical for proportions and grouping

Accessing the Underlying DataFrame

# Access the data within a Sample object
sample.data  # Returns the Polars DataFrame

---

Migration from samplics

svy replaces samplics with a fundamentally different API. Key migration points:

Design Specification

# samplics (OLD — archived)
from samplics.estimation import TaylorEstimator
estimator = TaylorEstimator("mean")
estimator.estimate(
    y=data["income"].to_numpy(),
    samp_weight=data["weight"].to_numpy(),
    stratum=data["stratum"].to_numpy(),
    psu=data["psu"].to_numpy()
)

# svy (NEW)
import svy
design = svy.Design(stratum="stratum", psu="psu", wgt="weight")
sample = svy.Sample(data=data, design=design)
result = sample.estimation.mean("income")

Key Differences

Why the Change Matters

samplics required passing design variables on every estimation call. svy's Sample object binds data and design once, then all estimation and regression methods automatically use the correct design. This reduces errors from inconsistent design specification across calls.

svy Estimation Reference

svy v0.13.0 — syntax and library guidance only.

---

Contents

1. Prerequisites: Design and Sample Setup 2. Population Means 3. Population Totals 4. Proportions 5. Ratios 6. Medians and Quantiles 7. Domain / Subpopulation Estimation 8. Cross-Tabulations 9. Hypothesis Testing (Survey-Weighted t-Tests) 10. Design Effects (DEFF) 11. Working with Polars DataFrames 12. Common Patterns and Pitfalls

---

Prerequisites: Design and Sample Setup

All estimation requires a svy.Sample object combining data with a design specification. See design-weights.md for full design setup. Brief recap:

import svy

# --- Taylor linearization design (most common) ---
design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec2yr")
sample = svy.Sample(data=data, design=design)

For replicate weight designs, see design-weights.md. Once the Sample is created, estimation methods are identical regardless of the variance estimation method — the design object determines how SEs are computed.

---

Population Means

Basic Mean

# Population mean of BMI with design-based SE
result = sample.estimation.mean("bmxbmi")
print(result)

The result includes: point estimate, standard error (SE), 95% confidence interval, and design effect (DEFF).

Multiple Variables

# Call mean() once per variable — no multi-variable shorthand
result_income = sample.estimation.mean("income")
result_age = sample.estimation.mean("age")

Handling Missing Data

# Drop nulls (default behavior)
result = sample.estimation.mean("bmxbmi", drop_nulls=True)

svy uses Polars null handling. Rows with null values in the analysis variable are excluded from the estimate by default. The effective sample size after dropping nulls is reported.

Critical: Dropping nulls changes the effective domain of estimation. If missingness is non-random (which it usually is in surveys), acknowledge this limitation in your analysis. Document it with an # ASSUMES: comment in research scripts.

---

Population Totals

Basic Total

# Estimated population total
result = sample.estimation.total("income")
print(result)

Totals estimate the sum of a variable across the entire target population, not just the sample. The SE reflects the uncertainty of this population-level estimate.

When to Use Totals vs. Means

• Totals for aggregate quantities: total enrollment, total expenditure, total population count
• Means for per-unit averages: mean income, mean BMI, mean test score
• Proportions for binary/categorical shares: percent employed, percent below poverty

---

Proportions

Basic Proportion

# Proportion of a binary or categorical variable
result = sample.estimation.prop("employed")
print(result)

The variable should contain categorical or binary values. svy computes the proportion in each category with design-based SEs.

Multi-Category Proportions

# Proportions across all categories of a variable
result = sample.estimation.prop("education_level")

This returns the estimated population proportion for each level of the variable (e.g., "High School": 0.28, "College": 0.35, "Graduate": 0.12, etc.) with SEs and CIs for each.

---

Ratios

Basic Ratio

# Ratio estimation: y is numerator, x is denominator
result = sample.estimation.ratio(y="total_expenditure", x="household_size")
print(result)

Ratio estimation is used when the quantity of interest is a ratio of two survey variables (e.g., per-capita expenditure = total expenditure / household size). The SE accounts for the covariance between numerator and denominator.

Ratio vs. Mean of a Derived Variable

Do not compute expenditure / household_size as a new column and then estimate its mean. This gives incorrect SEs because it ignores the covariance structure. Use estimation.ratio() for proper variance estimation of ratios.

# WRONG: pre-computing the ratio then estimating the mean
# data = data.with_columns((pl.col("expenditure") / pl.col("hh_size")).alias("per_capita"))
# sample.estimation.mean("per_capita")  # <-- incorrect SEs

# CORRECT: use ratio estimation
sample.estimation.ratio(y="expenditure", x="hh_size")  # <-- correct SEs

---

Medians and Quantiles

Median

# Population median with design-based SE
result = sample.estimation.median("income")
print(result)

Median estimation for survey data uses weighted quantile computation with linearization-based or replicate-weight-based variance estimation. SEs for medians are typically larger than for means.

---

Domain / Subpopulation Estimation

The by Parameter

Domain estimation computes statistics for subgroups of the population while preserving the full survey design structure.

# Mean BMI by gender
result = sample.estimation.mean("bmxbmi", by="riagendr")
print(result)

# Mean income by education level
result = sample.estimation.mean("income", by="education")

# Proportions by region
result = sample.estimation.prop("employed", by="region")

Why Not Pre-Filter?

Never pre-filter the data for domain estimation. Pre-filtering removes observations needed for correct variance estimation.

# WRONG: filtering before estimation
# females_only = data.filter(pl.col("gender") == "Female")
# female_sample = svy.Sample(data=females_only, design=design)
# female_sample.estimation.mean("income")  # <-- WRONG SEs

# CORRECT: use domain estimation
sample.estimation.mean("income", by="gender")  # <-- correct SEs for each gender

Pre-filtering discards PSUs and strata from the design, which can: 1. Produce incorrect variance estimates (too small or too large) 2. Create singleton PSU problems (strata with only one PSU after filtering) 3. Change the degrees of freedom for inference

The by parameter handles domain estimation correctly by keeping the full design structure and computing conditional estimates.

Multiple Grouping Variables

# Multiple grouping variables passed as a tuple
result = sample.estimation.mean("income", by=("gender", "education"))

---

Cross-Tabulations

Survey-Weighted Contingency Tables

# Cross-tabulation via prop() with by=
result = sample.estimation.prop("employment_status", by="education_level")

This produces a survey-weighted cross-tabulation showing the estimated population proportion in each cell, with design-based SEs. Equivalent to R's svytable() or svyby(~var, ~by_var, design, svymean).

For a full contingency table with chi-square test, use the categorical.tabulate() method:

# Formal cross-tabulation with test statistics
table = sample.categorical.tabulate(rowvar="employment_status", colvar="education_level")

---

Hypothesis Testing (Survey-Weighted t-Tests)

Comparing Domain Means

# Domain means via the by= parameter
result = sample.estimation.mean("income", by="gender")

# Formal two-group t-test via the categorical accessor
ttest_result = sample.categorical.ttest(y="income", group="gender")
print(ttest_result)

For formal hypothesis testing of differences between domains, svy computes design-adjusted t-statistics that account for the complex sampling structure via sample.categorical.ttest(). The group parameter specifies the binary grouping variable. The degrees of freedom are based on the number of PSUs minus the number of strata (not the sample size), which can substantially affect p-values for small numbers of clusters.

Key Difference from Unweighted Tests

Standard t-tests assume simple random sampling with known, equal variance. Survey-weighted tests:

• Use the survey weights in the point estimate
• Use the design-based variance (accounting for stratification and clustering)
• Use design-based degrees of freedom (typically much smaller than n - 1)
• Produce wider confidence intervals when there is substantial clustering

---

Design Effects (DEFF)

The design effect (DEFF) measures how much the variance of an estimate is inflated (or deflated) by the complex design compared to a simple random sample of the same size.

DEFF = Var_complex / Var_SRS

• DEFF = 1.0: The complex design is as efficient as SRS
• DEFF > 1.0: The complex design increases variance (common with clustered designs)
• DEFF < 1.0: The complex design decreases variance (common with stratified designs)
• Typical range: 1.5 to 5.0 for clustered household surveys

DEFF values are included in svy estimation output. Report them alongside estimates — they communicate how much the survey design affects precision.

Effective Sample Size

The effective sample size is:

n_eff = n / DEFF

A survey of 10,000 respondents with DEFF = 4.0 has the statistical precision of an SRS of only 2,500. Always consider the effective sample size when evaluating whether a survey has adequate power for a particular analysis.

---

Working with Polars DataFrames

svy uses Polars DataFrames natively. Data loaded via svy.io methods returns Polars DataFrames. If you have data in other formats:

From Parquet (Common in DAAF Pipelines)

import polars as pl
import svy

# Load data as Polars DataFrame
data = pl.read_parquet("data/raw/nhanes_demo.parquet")

# Proceed with design specification
design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec2yr")
sample = svy.Sample(data=data, design=design)

From Pandas

import pandas as pd
import polars as pl
import svy

# Convert pandas to Polars
pd_data = pd.read_csv("survey_data.csv")
data = pl.from_pandas(pd_data)

# Proceed with svy
design = svy.Design(stratum="stratum", psu="psu", wgt="weight")
sample = svy.Sample(data=data, design=design)

Data Wrangling Within svy

svy includes a wrangling module for common survey data preparation tasks. These operate on the Sample object directly, preserving the design linkage.

from svy import CaseStyle, LetterCase

# Clean column names
sample = sample.wrangling.clean_names(
    case_style=CaseStyle.SNAKE,
    letter_case=LetterCase.LOWER
)

# Recode categorical variables
sample = sample.wrangling.recode(
    "education",
    {"High School": ["HS", "high_school"],
     "College": ["BA", "BS", "college"]}
)

# Bin continuous variables into categories
sample = sample.wrangling.categorize(
    "age",
    bins=[0, 18, 35, 65, 100],
    labels=["0-17", "18-34", "35-64", "65+"]
)

# Cap extreme values (winsorize)
sample = sample.wrangling.bottom_and_top_code(
    {"income": (0, 200000)}
)

# Create derived variables
from svy.core.expr import col
sample = sample.wrangling.mutate({
    "income_thousands": col("income") / 1000,
    "age_squared": col("age") ** 2
})

Note: For complex data preparation (joins, reshaping, filtering), use the polars skill directly, then pass the prepared Polars DataFrame to svy.Sample. The wrangling module is for convenience on simple transformations.

---

Common Patterns and Pitfalls

Pattern: Complete Estimation Workflow

import svy
import polars as pl

# --- Config ---
DATA_PATH = "data/raw/nhanes_demo.parquet"

# --- Load ---
data = pl.read_parquet(DATA_PATH)

# --- Design ---
# INTENT: NHANES uses a complex multi-stage stratified cluster design
# REASONING: sdmvstra = pseudo-strata, sdmvpsu = pseudo-PSU, wtmec2yr = 2-year MEC exam weight
# ASSUMES: Analysis population is the MEC-examined subsample
design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec2yr")
sample = svy.Sample(data=data, design=design)

# --- Estimate ---
mean_bmi = sample.estimation.mean("bmxbmi")
print(mean_bmi)

mean_bmi_by_gender = sample.estimation.mean("bmxbmi", by="riagendr")
print(mean_bmi_by_gender)

prop_obese = sample.estimation.prop("obese_flag")
print(prop_obese)

# --- Validate ---
print(f"Sample size: {data.shape[0]}")
assert data.shape[0] > 0, "No data loaded"

Pitfall: Using Unweighted Statistics

Never use pl.col("var").mean() or pandas .mean() on survey data. Unweighted statistics are biased for the target population and do not have correct standard errors.

Pitfall: Ignoring Weight Variable Selection

Surveys often provide multiple weight variables for different analysis populations (e.g., NHANES has wtint2yr for interview data and wtmec2yr for examination data). Using the wrong weight produces biased estimates. Always consult the survey documentation to select the appropriate weight.

Pitfall: Treating Survey SEs as Cluster-Robust SEs

Survey-weighted SEs and cluster-robust SEs (e.g., from pyfixest or statsmodels with cov_type="cluster") are not the same thing:

• Survey SEs account for stratification, clustering, and unequal probability of selection
• Cluster-robust SEs only account for within-cluster correlation
• Survey SEs use design-based degrees of freedom; cluster-robust SEs use large-sample approximations

Use svy when you have a complex survey with known design variables. Use cluster-robust SEs when you have non-survey data with clustered observations.

svy Regression Reference

svy v0.13.0 — syntax and library guidance only.

---

Contents

1. Overview: Survey-Weighted GLM 2. Survey-Weighted Linear Regression (Gaussian) 3. Survey-Weighted Logistic Regression (Binomial) 4. Survey-Weighted Poisson Regression 5. Specifying Predictors 6. Extracting Results 7. Survey Regression vs. WLS vs. Cluster-Robust SEs 8. Diagnostics and Model Fit 9. Domain-Specific Regression 10. The rpy2 Bridge for Unsupported Models 11. Complete Regression Workflow Example

---

Overview: Survey-Weighted GLM

svy fits generalized linear models (GLMs) that account for the complex survey design in both point estimation and variance estimation. The interface is through sample.glm.fit():

model = sample.glm.fit(
    y="outcome_variable",
    x=["predictor1", "predictor2", svy.Cat("categorical_var")],
    family="gaussian"  # or "binomial" or "poisson"
)

Supported families:

Not supported in svy (use rpy2 + R survey package):

• Negative binomial
• Ordinal logistic (svyolr in R)
• Cox proportional hazards (svycoxph in R)
• Multinomial logistic
• Quasi-families

The design-based SEs produced by sample.glm.fit() are model-robust "sandwich" estimators that account for stratification, clustering, and unequal probability of selection. These are analogous to R's survey::svyglm().

---

Survey-Weighted Linear Regression (Gaussian)

Basic Linear Model

import svy
import polars as pl

# Load and set up design (see design-weights.md)
data = pl.read_parquet("data/raw/nhanes.parquet")
design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec2yr")
sample = svy.Sample(data=data, design=design)

# INTENT: Estimate association between age and BMI, controlling for gender
# REASONING: Linear model appropriate for continuous outcome (BMI)
# ASSUMES: Linear relationship between age and BMI within the modeled range
model = sample.glm.fit(
    y="bmxbmi",
    x=["ridageyr", svy.Cat("riagendr")],
    family="gaussian"
)
print(model)

Multiple Continuous Predictors

model = sample.glm.fit(
    y="bmxbmi",
    x=["ridageyr", "indfmpir", "lbxtc"],
    family="gaussian"
)

Interpreting Gaussian GLM Output

The coefficients from a gaussian-family GLM are interpreted identically to OLS:

• Each coefficient represents the expected change in Y for a one-unit change in X, holding other variables constant
• The intercept is the expected value of Y when all predictors are zero
• SEs, t-statistics, and p-values reflect the complex survey design, not simple random sampling assumptions

---

Survey-Weighted Logistic Regression (Binomial)

Basic Logistic Model

# INTENT: Model probability of obesity as a function of demographics
# REASONING: Binary outcome (obese yes/no) requires logistic regression
# ASSUMES: The outcome variable is coded 0/1
model = sample.glm.fit(
    y="obese_flag",
    x=["ridageyr", "indfmpir", svy.Cat("riagendr"), svy.Cat("ridreth1")],
    family="binomial"
)
print(model)

Interpreting Logistic GLM Output

• Coefficients are on the log-odds scale
• To get odds ratios: exponentiate the coefficients (exp(coef))
• A positive coefficient means higher odds of the outcome for a one-unit increase in the predictor
• SEs are on the log-odds scale; exponentiate the confidence interval bounds for OR CIs

import numpy as np

# sample.glm.fit() returns the GLM accessor; the GLMFit result is at glm.fitted
glm = sample.glm.fit(y="obese_flag", x=["ridageyr", svy.Cat("riagendr")], family="binomial")
fit = glm.fitted                          # GLMFit object

# Convert to odds ratios via the GLMFit.to_polars() DataFrame
coef_df = fit.to_polars()                 # columns: term, estimate, std_err, conf_low, conf_high, statistic, p_value, df
odds_ratios = np.exp(coef_df["estimate"])
or_ci_lower = np.exp(coef_df["conf_low"])
or_ci_upper = np.exp(coef_df["conf_high"])

# Or access individual coefficients via fit.coefs (list of GLMCoef objects)
for c in fit.coefs:
    print(f"{c.term}: OR={np.exp(c.est):.3f}  [{np.exp(c.lci):.3f}, {np.exp(c.uci):.3f}]")

Important: Quasi-Binomial in R vs. svy

In R's survey::svyglm(), family=quasibinomial() is preferred over family=binomial() to avoid warnings about non-integer successes (which arise from weighted data). svy handles this internally — use family="binomial" directly.

---

Survey-Weighted Poisson Regression

Basic Poisson Model

# INTENT: Model count of doctor visits as a function of health indicators
# REASONING: Count outcome with no upper bound suits Poisson regression
# ASSUMES: Conditional mean equals conditional variance (Poisson assumption)
model = sample.glm.fit(
    y="doctor_visits",
    x=["ridageyr", svy.Cat("health_status"), svy.Cat("insurance_status")],
    family="poisson"
)
print(model)

Interpreting Poisson GLM Output

• Coefficients are on the log scale
• To get incidence rate ratios (IRRs): exponentiate the coefficients
• A coefficient of 0.3 means exp(0.3) = 1.35, i.e., a 35% higher rate for a one-unit increase
• SEs are on the log scale; exponentiate CI bounds for IRR CIs

Overdispersion Caveat

Poisson regression assumes the conditional variance equals the conditional mean. Survey data often exhibits overdispersion (variance > mean). The design-based SEs from svy partially accommodate this because they are sandwich estimators, but severe overdispersion may still produce misleading inference. Consider: 1. Checking whether the mean-variance assumption is reasonable 2. For negative binomial (overdispersed Poisson), use rpy2 + R's survey::svyglm(family=quasipoisson()) or MASS::glm.nb with survey design

---

Specifying Predictors

Continuous Predictors

Pass variable names as strings in the x list:

x=["age", "income", "bmi"]

Categorical Predictors with svy.Cat

Wrap categorical variable names in svy.Cat() to get proper dummy coding:

x=["age", svy.Cat("education"), svy.Cat("region")]

svy.Cat() creates indicator (dummy) variables for each level, with one level omitted as the reference category. The reference category is typically the first level alphabetically or numerically.

Mixing Continuous and Categorical

model = sample.glm.fit(
    y="income",
    x=[
        "age",                          # continuous
        "years_education",              # continuous
        svy.Cat("gender"),              # categorical (2 levels)
        svy.Cat("race_ethnicity"),      # categorical (5 levels)
        svy.Cat("marital_status"),      # categorical (3 levels)
    ],
    family="gaussian"
)

Interactions

# Pre-compute interaction columns in Polars before passing to svy
data = data.with_columns(
    (pl.col("age") * pl.col("income")).alias("age_x_income")
)
sample = svy.Sample(data=data, design=design)
model = sample.glm.fit(
    y="bmi",
    x=["age", "income", "age_x_income"],
    family="gaussian"
)

---

Extracting Results

Model Output Structure

The model object returned by sample.glm.fit() contains point estimates, standard errors, test statistics, p-values, and confidence intervals for each coefficient.

model = sample.glm.fit(y="bmxbmi", x=["ridageyr", svy.Cat("riagendr")], family="gaussian")

# Print full summary
print(model)

Accessing Individual Components

sample.glm.fit() returns the GLM accessor. The fitted result (GLMFit) is stored at glm.fitted:

glm = sample.glm.fit(y="bmxbmi", x=["ridageyr", svy.Cat("riagendr")], family="gaussian")
fit = glm.fitted                          # GLMFit object

# --- Tabular output (recommended) ---
coef_df = fit.to_polars()
# columns: term, estimate, std_err, conf_low, conf_high, statistic, p_value, df

# --- Individual coefficient objects ---
for c in fit.coefs:                       # list of GLMCoef
    print(c.term, c.est, c.se, c.lci, c.uci)
    # c.wald  — TDist with .value (t-statistic) and .p_value

# --- Model statistics ---
stats = fit.stats                         # GLMStats
print(stats.n, stats.r_squared, stats.aic, stats.bic, stats.deviance)

Important: sample.glm creates a new GLM object each time it is accessed (it is a property). You must hold a reference to the return value of .fit() to access glm.fitted and glm.predict() later.

Reporting Results

When reporting survey regression results: 1. State the survey design (strata, PSU, weight variable) 2. Report the variance estimation method (Taylor linearization or replicate type) 3. Report design-based degrees of freedom (not sample size minus parameters) 4. Report coefficients with design-based SEs and CIs 5. Note the effective sample size if substantially smaller than n

---

Survey Regression vs. WLS vs. Cluster-Robust SEs

This distinction is critical. Three approaches look superficially similar but are fundamentally different:

1. Survey-Weighted Regression (svy)

# svy: accounts for strata, PSU, and unequal selection probabilities
design = svy.Design(stratum="strata", psu="psu", wgt="weight")
sample = svy.Sample(data=data, design=design)
model = sample.glm.fit(y="y", x=["x1", "x2"], family="gaussian")

What it does:

• Uses weights in point estimation (weighted least squares for coefficients)
• Uses the full design (strata + PSU + weights) for variance estimation
• Produces design-based degrees of freedom (# PSUs - # strata)
• Accounts for stratification (reduces variance), clustering (increases variance), and unequal selection

2. Weighted Least Squares (statsmodels WLS)

# statsmodels WLS: uses weights for point estimates only
import statsmodels.formula.api as smf
results = smf.wls("y ~ x1 + x2", data=df, weights=df["weight"]).fit()

What it does:

• Uses weights in point estimation
• Does NOT account for stratification or clustering in variance
• Assumes independent observations with heterogeneous variance
• Produces incorrect SEs for survey data — typically too small

3. Cluster-Robust Standard Errors (statsmodels / pyfixest)

# Cluster-robust SEs: accounts for clustering only
import statsmodels.formula.api as smf
results = smf.ols("y ~ x1 + x2", data=df).fit(cov_type="cluster", cov_kwds={"groups": df["psu"]})

What it does:

• Does NOT use survey weights in point estimation (unless combined with WLS)
• Accounts for within-cluster correlation in variance
• Does NOT account for stratification or unequal selection
• Approximation that works for some settings but is not a substitute for design-based inference

Summary Comparison

Rule of thumb: If you have data from a complex survey with known design variables, use svy. If you have non-survey data with clustered observations, cluster-robust SEs are appropriate.

---

Diagnostics and Model Fit

R-Squared in Survey Context

Traditional R-squared is not well-defined for survey-weighted regression because the weights change the effective sample. Some implementations report a pseudo-R-squared. Interpret with caution.

Residual Analysis

# Residuals and fitted values are obtained via glm.predict()
glm = sample.glm.fit(y="bmxbmi", x=["ridageyr", svy.Cat("riagendr")], family="gaussian")

# Pass y_col to get residuals (without it, pred.residuals is None)
pred = glm.predict(new_data=data, y_col="bmxbmi")
residuals = pred.residuals               # numpy array
fitted_values = pred.yhat                 # numpy array

# As a Polars DataFrame:
pred_df = pred.to_polars()                # columns: yhat, se, lci, uci, residuals

Residual plots for survey regressions should use weighted residuals. Unweighted residual plots can be misleading because they treat all observations equally regardless of their representation of the population.

Specification Concerns

Standard diagnostic tests (Breusch-Pagan, RESET, VIF) from statsmodels are designed for simple random samples. For survey data:

• Multicollinearity: VIF computed on the unweighted design matrix is still informative for detecting collinearity, though not for assessing its impact on survey SEs
• Functional form: Plot weighted residuals against predictors visually; formal tests require survey-adjusted versions
• Influential observations: Observations with large weights are inherently influential by design — they represent more of the population. Do not simply remove them without substantive justification

Model Comparison

glm = sample.glm.fit(y="bmxbmi", x=["ridageyr", svy.Cat("riagendr")], family="gaussian")
fit = glm.fitted

# AIC/BIC are on the GLMStats object (fit.stats), not the model directly
print(f"AIC: {fit.stats.aic}")
print(f"BIC: {fit.stats.bic}")

For comparing nested survey regression models, use Wald tests based on the design-based covariance matrix rather than likelihood ratio tests (which assume independent observations).

---

Domain-Specific Regression

To fit a regression model for a subpopulation (domain), do NOT pre-filter the data. Instead, use domain estimation:

# WRONG: pre-filtering breaks the design
# females = data.filter(pl.col("gender") == 2)
# female_sample = svy.Sample(data=females, design=design)
# model = female_sample.glm.fit(...)  # <-- WRONG SEs

Limitation (v0.13.0): Domain-restricted regression (fitting a model within a subpopulation while preserving the full design) is not supported in sample.glm.fit(). The fit() method has no by=, subset=, or where= parameter. For domain estimation of means/totals/proportions, use sample.estimation.mean(..., by="group"), which does support correct domain analysis. For regression within a subpopulation, the only current option is to pre-filter the data, with the caveat that SEs will not fully account for the original design structure. Document this limitation in your analysis with an # ASSUMES: comment.

Domain regression preserves the full design structure for variance estimation while restricting the model to the subpopulation of interest. This is methodologically equivalent to R's svyglm(..., design=subset(design, gender == 2)) where the subsetting is done within the survey design framework. Monitor future svy releases for this feature.

---

The rpy2 Bridge for Unsupported Models

When svy does not support the needed model (ordinal logistic, Cox survival, negative binomial, etc.), use R's survey package via rpy2.

Setup

# rpy2 bridge pattern — this is standard rpy2, not svy-specific
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr

pandas2ri.activate()
survey_r = importr("survey")

# Convert Polars to pandas for rpy2 (rpy2 uses pandas bridge)
df_pandas = data.to_pandas()

# Create R survey design
r_design = survey_r.svydesign(
    ids=ro.Formula("~psu_id"),
    strata=ro.Formula("~stratum"),
    weights=ro.Formula("~weight"),
    data=df_pandas,
    nest=True
)

Ordinal Logistic

# Ordinal logistic regression in R via rpy2
r_model = survey_r.svyolr(
    ro.Formula("ordered_outcome ~ age + education"),
    design=r_design
)
print(ro.r.summary(r_model))

Cox Proportional Hazards

# Cox PH model in R via rpy2
r_model = survey_r.svycoxph(
    ro.Formula("Surv(time, event) ~ age + treatment"),
    design=r_design
)
print(ro.r.summary(r_model))

Negative Binomial (via Quasi-Poisson)

# Quasi-Poisson as overdispersion-robust alternative
r_model = survey_r.svyglm(
    ro.Formula("count ~ age + region"),
    design=r_design,
    family=ro.r("quasipoisson()")
)
print(ro.r.summary(r_model))

Decision rule: If family is not "gaussian", "binomial", or "poisson", use the rpy2 bridge. Document the bridge usage with # REASONING: comments explaining why svy was insufficient.

---

Complete Regression Workflow Example

import svy
import polars as pl
import numpy as np

# --- Config ---
DATA_PATH = "data/raw/2026-03-27_nhanes_demo_exam.parquet"
OUTPUT_PATH = "output/analysis/2026-03-27_bmi_regression.parquet"

# --- Load ---
data = pl.read_parquet(DATA_PATH)
print(f"Loaded {data.shape[0]} rows, {data.shape[1]} columns")

# --- Design ---
# INTENT: NHANES 2017-2020 uses a complex multi-stage probability design
# REASONING: sdmvstra/sdmvpsu are masked design variables; wtmec2yr is the
#   2-year MEC exam weight appropriate for variables collected during the exam
# ASSUMES: Analysis restricted to MEC-examined participants aged 20+
design = svy.Design(stratum="sdmvstra", psu="sdmvpsu", wgt="wtmec2yr")
sample = svy.Sample(data=data, design=design)

# --- Transform ---
# INTENT: Create age-squared term for nonlinear age effect on BMI
from svy.core.expr import col
sample = sample.wrangling.mutate({
    "age_sq": col("ridageyr") ** 2
})

# --- Analysis ---
# INTENT: Estimate association between demographics and BMI
# REASONING: Gaussian family for continuous outcome; design-based SEs
# ASSUMES: Linear in parameters; additive effects; no unmeasured confounders
model = sample.glm.fit(
    y="bmxbmi",
    x=[
        "ridageyr",
        "age_sq",
        "indfmpir",
        svy.Cat("riagendr"),
        svy.Cat("ridreth1"),
    ],
    family="gaussian"
)
print(model)

# --- Validate ---
# Verify model ran on expected sample size
print(f"Input data rows: {data.shape[0]}")
fit = model.fitted
print(f"Model N: {fit.stats.n}")
print(f"R-squared: {fit.stats.r_squared:.4f}")
print(f"AIC: {fit.stats.aic:.2f}")