R Survey Skill

R's survey package by Thomas Lumley -- the canonical tool for design-based analysis of complex survey data. Covers survey design specification (strata, PSU, weights, FPC), variance estimation (Taylor linearization, BRR, jackknife, bootstrap), descriptive estimation (means, totals, proportions, ratios, quantiles), survey-weighted GLM regression (gaussian, binomial, Poisson, ordinal, Cox PH), domain/subpopulation analysis, calibration, and replicate weight handling. The Python svy package is modeled after this R implementation. Use when execution language is R and the data comes from a complex sample survey (NHANES, CPS, ACS PUMS, BRFSS, DHS, ECLS-K, MEPS-HC). For non-survey regression in R, use r-stats; for fixed effects use fixest; for panel RE/between use plm.

What is the survey Package?

The survey package is the definitive implementation of design-based inference for complex survey data. First released in 2003, it is the reference implementation that Python's svy and Stata's svy: prefix are modeled after.

• Design object: svydesign() -- specifies strata, PSU, weights, FPC
• Survey-aware estimation: svymean(), svytotal(), svyratio(), svyquantile() with design-based standard errors
• GLM regression: svyglm() for linear, logistic, Poisson, and more with design-adjusted inference
• Domain estimation: svyby() for correct subpopulation analysis without pre-filtering
• Replicate weights: svrepdesign() for BRR, jackknife, bootstrap, Fay
• Calibration: calibrate(), postStratify(), rake() for weight adjustment
• Model extensions: svyolr() for ordinal logistic, svycoxph() for survival analysis -- models not available in Python's svy

Version Notes

This skill targets survey 4.5 on R 4.5.3. Key features available:

• svydesign() and svrepdesign() for all design types
• svyglm() with quasi-family support (quasibinomial, quasipoisson)
• as.svrepdesign() to convert Taylor designs to replicate weight designs
• svycontrast() for custom linear/nonlinear contrasts
• Lonely PSU handling via options(survey.lonely.psu = ...)
• svyby() for domain estimation across arbitrary grouping variables

Note on emmeans: The emmeans package is NOT installed in this environment. For marginal means or pairwise comparisons from survey models, use svycontrast() for custom contrasts, or compute predicted values manually via predict(svyglm_fit, newdata = ...). See ./references/regression.md for patterns.

How to Use This Skill

Reference File Structure

File	Purpose	When to Read
quickstart.md	svydesign(), svymean(), svytotal(), confint(), basic workflow	Starting with survey or need design setup
estimation.md	Point estimates, ratios, quantiles, svyby() for domains, svycontrast()	Descriptive survey statistics
regression.md	svyglm() for linear/logit/Poisson, update(), summary output, marginal effects without emmeans	Survey regression models
replication.md	svrepdesign(), as.svrepdesign(), BRR, jackknife, bootstrap, replicate weight surveys	Working with replicate weights
domains.md	svyby() domain estimation, interaction of strata with subgroups, conditional analysis	Subpopulation analysis
gotchas.md	Lonely PSU handling, degrees of freedom, WLS-is-not-survey, subsetting (subset= not filter), calibration	Debugging common issues

Reading Order

1. New to survey in R? Start with quickstart.md then estimation.md
2. Need survey-weighted regression? Read quickstart.md then regression.md
3. Have replicate weights? Read replication.md then estimation.md or regression.md
4. Need subpopulation estimates? Read domains.md
5. Setting up a federal survey (NHANES, CPS)? Read quickstart.md (federal survey patterns table)
6. Coming from Python svy? Read quickstart.md -- the R API is the original that svy was modeled after

Related Skills

Skill	Relationship
svy	Python equivalent for complex survey analysis. Python's svy is modeled after R's survey. If execution language is Python, use svy instead
data-scientist	Provides methodology guidance (especially survey-analysis.md); survey-r provides R implementation. Load data-scientist for "when and why"
r-stats	Complement for non-survey regression (OLS, GLM, diagnostics). base R lm() with weights is NOT survey-weighted regression -- it does not account for stratification or clustering
fixest	Complement for fixed effects models and DiD. fixest does not handle complex survey designs; use survey-r for survey-weighted estimation
plm	Complement for panel models (RE, between, first difference). Does not handle survey designs
tidyverse	Use for data preparation before passing to survey design objects. survey uses base R data.frames

Quick Decision Trees

"I need to analyze survey data in R"

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
|   +-- Ordinal logistic --> ./references/regression.md (svyolr)
|   +-- Cox survival --> ./references/regression.md (svycoxph)
+-- Set up the survey design object
|   +-- ./references/quickstart.md
+-- Subpopulation / domain analysis
|   +-- ./references/domains.md
+-- Replicate weight design (BRR, jackknife, bootstrap)
|   +-- ./references/replication.md
+-- Something isn't working
    +-- ./references/gotchas.md

"I need survey-weighted regression in R"

What model?
+-- Linear regression (continuous Y)
|   +-- svyglm(family = gaussian()) --> ./references/regression.md
+-- Logistic regression (binary Y)
|   +-- svyglm(family = quasibinomial()) --> ./references/regression.md
+-- Poisson regression (count Y)
|   +-- svyglm(family = quasipoisson()) --> ./references/regression.md
+-- Ordinal logistic (ordered categories)
|   +-- svyolr() --> ./references/regression.md
+-- Cox proportional hazards (survival)
|   +-- svycoxph() --> ./references/regression.md
+-- Negative binomial (overdispersed counts)
|   +-- svyglm(family = quasipoisson()) as approximation
|   +-- Or: MASS::glm.nb() ignores survey design -- use with caution
+-- Fixed effects + survey weights
    +-- Methodologically complex -- consult data-scientist skill

"I need to set up variance estimation"

What do you have?
+-- Design variables (strata, PSU, weights)
|   +-- Taylor linearization via svydesign() --> ./references/quickstart.md
+-- Pre-computed replicate weights
|   +-- BRR --> ./references/replication.md
|   +-- Jackknife --> ./references/replication.md
|   +-- Bootstrap --> ./references/replication.md
+-- Need to create replicate weights from design
|   +-- as.svrepdesign() --> ./references/replication.md
+-- Not sure what I have
    +-- Read survey documentation --> ./references/quickstart.md (federal survey
        table)

Boundaries

survey covers:

• Design-based estimation (descriptive and regression) for complex surveys
• Taylor and replicate-weight variance estimation
• Domain/subpopulation analysis via svyby() and subset()
• Calibration (post-stratification, raking, GREG)
• Ordinal logistic (svyolr), Cox PH (svycoxph) -- models NOT in Python svy
• Interaction with broom::tidy() for tidy output

survey does NOT cover (use other tools):

• Fixed effects models -- use fixest (survey weights + FE is methodologically complex; consult data-scientist skill)
• Panel data models (RE, FD, between) -- use plm
• Difference-in-differences -- use fixest
• Machine learning -- use tidymodels or scikit-learn via reticulate
• Survey sampling design and sample size calculation -- use data-scientist skill for methodology
• Marginal means via emmeans -- emmeans is NOT installed; use svycontrast() or manual predict() instead

File-First Execution in Research Workflows

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

The pattern:

1. Write survey analysis code to scripts/stage8_analysis/{step}_{task-name}.R
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 Library Loading

# --- Config ---
library(survey)

Core Workflow

# 1. Load data
data <- arrow::read_parquet("data/raw/nhanes_demo.parquet")

# 2. Specify 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
des <- svydesign(
  ids = ~sdmvpsu,
  strata = ~sdmvstra,
  weights = ~wtmec2yr,
  data = data,
  nest = TRUE
)

# 3. Estimate
svymean(~bmxbmi, design = des, na.rm = TRUE)
svytotal(~bmxbmi, design = des, na.rm = TRUE)
confint(svymean(~bmxbmi, design = des, na.rm = TRUE))

# 4. Regression
fit <- svyglm(bmxbmi ~ ridageyr + factor(riagendr), design = des)
summary(fit)

Core Operations

Operation	Code
Design (Taylor)	svydesign(ids = ~psu, strata = ~strat, weights = ~wt, data = df, nest = TRUE)
Design (replicate)	svrepdesign(weights = ~wt, repweights = "wt[0-9]+", type = "BRR", data = df)
Mean	svymean(~var, design = des, na.rm = TRUE)
Total	svytotal(~var, design = des, na.rm = TRUE)
Proportion	svymean(~factor(var), design = des)
Ratio	svyratio(~num, ~denom, design = des)
Quantile	svyquantile(~var, design = des, quantiles = 0.5)
Domain estimation	svyby(~var, ~group, design = des, svymean)
Linear regression	svyglm(y ~ x, design = des, family = gaussian())
Logistic regression	svyglm(y ~ x, design = des, family = quasibinomial())
Poisson regression	svyglm(y ~ x, design = des, family = quasipoisson())
Ordinal logistic	svyolr(ordered(y) ~ x, design = des)
Cox PH	svycoxph(Surv(time, event) ~ x, design = des)
Subset design	subset(des, age >= 18)
Confidence intervals	confint(svymean(...))
Tidy output	broom::tidy(svyglm_fit, conf.int = TRUE)
Convert to replicate	as.svrepdesign(des, type = "JKn")
Calibrate	calibrate(des, formula = ~age_group, population = pop_totals)

Formula Syntax in survey

# Single variable
svymean(~income, design = des)

# Multiple variables
svymean(~income + age, design = des)

# Factor variable (for proportions)
svymean(~factor(education), design = des)

# Domain estimation
svyby(~income, ~gender, design = des, svymean)

# Regression with interaction
svyglm(y ~ x1 * factor(x2), design = des)

Topic Index

Topic	Reference File
svydesign() setup	./references/quickstart.md
Taylor linearization	./references/quickstart.md
svymean() / svytotal()	./references/quickstart.md
confint() for survey estimates	./references/quickstart.md
Federal survey design patterns	./references/quickstart.md
Point estimates	./references/estimation.md
Proportions	./references/estimation.md
Ratios	./references/estimation.md
Quantiles / medians	./references/estimation.md
svyby() for domain estimation	./references/estimation.md
svycontrast()	./references/estimation.md
Design effects (DEFF)	./references/estimation.md
Cross-tabulations (svytable)	./references/estimation.md
Chi-squared test (svychisq)	./references/estimation.md
svyglm() linear	./references/regression.md
svyglm() logistic	./references/regression.md
svyglm() Poisson	./references/regression.md
svyolr() ordinal logistic	./references/regression.md
svycoxph() Cox PH	./references/regression.md
Model comparison / update()	./references/regression.md
Marginal effects without emmeans	./references/regression.md
Odds ratios and IRR	./references/regression.md
svrepdesign()	./references/replication.md
as.svrepdesign()	./references/replication.md
BRR / Fay's method	./references/replication.md
Jackknife (JK1, JKn)	./references/replication.md
Bootstrap replicate weights	./references/replication.md
ACS PUMS successive difference	./references/replication.md
Domain estimation (svyby)	./references/domains.md
Subsetting designs (subset=)	./references/domains.md
Interaction domains	./references/domains.md
Lonely PSU handling	./references/gotchas.md
Degrees of freedom	./references/gotchas.md
WLS is not survey regression	./references/gotchas.md
subset() vs filter()	./references/gotchas.md
Calibration (calibrate, rake, postStratify)	./references/gotchas.md
Weight selection	./references/gotchas.md

Citation

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

Lumley, T. (2004). "Analysis of Complex Survey Samples." Journal of Statistical Software, 9(8), 1-19.

Lumley, T. (2023). survey: Analysis of Complex Survey Samples. R package version 4.5.

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