Author

Mamadou S. Diallo, Ph.D.

Modified

August 23, 2026

NHANES Case Study Outline — svy

Companion to the DHS Ethiopia case study. Demonstrates svy end to end on real NHANES data, from XPT import through survey-weighted logistic regression, with a comorbidity-construction section as the centrepiece.

All variable names below were verified against the 2017–2018 (cycle J) codebooks on 2026-08-23. All svy API calls were verified against svy 0.24.1.

0. Scope decisions (state these up front in the case study)

Decision Choice Why
Cycles 2017–2018 only (cycle J) Avoids multi-cycle weight rescaling entirely. Pooling is a separate, harder topic — keep it out.
Weight WTMEC2YR Analysis uses MEC exam + lab data, so the MEC weight, not WTINT2YR.
Fasting labs Excluded LBXGLU requires WTSAF2YR and would force the whole analysis onto the fasting subsample. Use LBXGH (HbA1c) instead — full MEC sample.
Population Adults 20+ The MCQ160 series is only asked of 20+. Matches the denominator to the questions.
Subsetting where= everywhere, never filter CDC’s strongest analytic warning: filtering treats within-cluster subgroup size as fixed and underestimates variance.

Callout worth making explicitly in the narrative: svy nests PSU within stratum automatically. NHANES reuses SDMVPSU = 1/2 inside every stratum, so R users must remember nest = TRUE or get wrong variances. svy follows Stata here — the footgun does not exist.

1. Data acquisition

Files needed (all .XPT, 2017–2018 = _J suffix):

File Pulls
DEMO_J SEQN, SDMVSTRA, SDMVPSU, WTMEC2YR, RIDAGEYR, RIAGENDR, RIDRETH3, DMDEDUC2, INDFMPIR
BMX_J BMXBMI
BPX_J BPXSY1–BPXSY4, BPXDI1–BPXDI4
BPQ_J BPQ020, BPQ050A, BPQ080, BPQ100D
DIQ_J DIQ010, DIQ050, DIQ070
GHB_J LBXGH (HbA1c)
MCQ_J MCQ010, MCQ160b–MCQ160g, MCQ160k–MCQ160o, MCQ220
KIQ_U_J KIQ022

svy.read_sas() auto-dispatches .xpt / .xport to the native Rust XPT reader, so no conversion step and no haven equivalent needed.

import svy
demo = svy.read_sas("DEMO_J.XPT")

Show in the case study: value labels and variable labels survive the XPT read. This is a genuine differentiator over reading NHANES with plain pandas.

2. Merging modules

svy has no join — merge with polars on SEQN before constructing the Sample. Keep this short in the narrative; it is not the interesting part.

df = demo.join(bmx, on="SEQN", how="left").join(bpq, on="SEQN", how="left")  # etc.

Left joins from DEMO_J, not inner — an inner join silently drops anyone who skipped a module and quietly changes the denominator.

3. Comorbidities — the centrepiece

3.1 The three-source pattern

Each chronic condition in NHANES is constructed as an OR across up to three sources. Make this the organising idea of the section:

condition = self-reported diagnosis  OR  measured/lab threshold  OR  on medication

Using only self-report undercounts undiagnosed disease. Using only the lab misses controlled disease. The OR is the standard operationalisation.

3.2 Condition definitions (all codes verified, cycle J)

Condition Definition
Hypertension mean of valid BPXSY1–BPXSY4 ≥ 130 or mean BPXDI1–BPXDI4 ≥ 80 or BPQ050A == 1
Diabetes DIQ010 == 1 or LBXGH ≥ 6.5 or DIQ050 == 1 or DIQ070 == 1
Obesity BMXBMI ≥ 30
High cholesterol BPQ080 == 1 or BPQ100D == 1
Cardiovascular disease any of MCQ160b (CHF), MCQ160c (CHD), MCQ160d (angina), MCQ160e (MI), MCQ160f (stroke) == 1
Chronic lung disease any of MCQ160g (emphysema), MCQ160k (chronic bronchitis), MCQ160o (COPD) == 1; MCQ010 (asthma) kept separate
Kidney disease KIQ022 == 1
Cancer MCQ220 == 1

Then: n_conditions = sum of flags, and multimorbidity = n_conditions >= 2.

Note the BP threshold is a choice: 130/80 is the 2017 ACC/AHA definition, 140/90 the older one. State which you used — prevalence differs substantially.

3.3 The four traps (this is the real content)

(a) Sentinel codes. Every item above uses 7 = Refused, 9 = Don't know. These must become null before any indicator is built. Otherwise MCQ160c == 9 evaluates as “not 1” and a don’t-know is silently counted as no disease. This is the single largest source of wrong published NHANES prevalence. XPT files frequently do not carry these as user-missing, so zap_missing will not catch them — declare them per variable.

(b) Skip patterns — null is often a legitimate “no”. BPQ050A (now taking BP medicine) is only asked of respondents who said yes to BPQ020. A null there means not asked, i.e. not on medication — not missing. Same for BPQ100D. Treating structural nulls as item nonresponse corrupts the denominator. svy’s MissingKind.STRUCTURAL vs DONT_KNOW / REFUSED is exactly the right vocabulary for this distinction.

(c) Age eligibility differs per item. MCQ160 series and KIQ022 are 20+; BPQ020/BPQ050A/BPQ080/BPQ100D are 16+; DIQ010/DIQ050/DIQ070 are 1+. Restricting to 20+ makes them consistent — say so rather than letting it happen silently.

(d) The weight must match the most restrictive source used. Every source above is MEC or interview, so WTMEC2YR is correct. Adding one fasting lab would force WTSAF2YR on the entire analysis and shrink the sample. This is why fasting glucose was excluded in §0.

3.4 svy surfaces to demonstrate

sample.wrangling: recode, categorize, mutate, cast, fill_null, apply_labels, top_code / bottom_code.

4. Design object

design = svy.Design(stratum="SDMVSTRA", psu="SDMVPSU", wgt="WTMEC2YR")
sample = svy.Sample(df, design)

If any domain produces a single-PSU stratum, svy raises rather than silently under-reporting variance. Show sample.singleton.detected() and sample.singleton.handle("center") — the latter is Stata’s singleunit(centered). Worth demonstrating: svy is stricter than R’s default here, and that is the correct behaviour.

5. Prevalence estimation

# Overall and by subgroup
sample.estimation.prop("diabetes", ci_method="korn-graubard")
sample.estimation.prop("diabetes", by="RIAGENDR", ci_method="korn-graubard")

# Subpopulation — correct domain SEs, not a filter
sample.estimation.prop("diabetes", where=svy.col("RIDAGEYR") >= 40,
                       ci_method="korn-graubard")

# Cross-tab
sample.categorical.tabulate("multimorbidity", "RIDRETH3", units="percent")

Say which CI method and why. ci_method options are logit (default, matches Stata svy: prop), beta (matches R svyciprop(method="beta")), korn-graubard (matches the NCHS SAS reference macro), wilson.

beta and korn-graubard are both Korn–Graubard and differ only in two places, both of which bite in small domains and rare conditions:

  • beta leaves n_eff uncapped; korn-graubard truncates at n. Uncapped n_eff can exceed the number of people observed — one synthetic case gave n = 240, n_eff = 12,879, producing an interval 7× too narrow.
  • At p = 0 or p = 1, beta returns a degenerate zero-width interval; korn-graubard returns a proper one-sided interval.

For a rare comorbidity in a small domain — which is most of this section — use korn-graubard.

6. Survey-weighted logistic regression

m = sample.glm.fit(
    "diabetes",
    x=["RIDAGEYR", svy.Cat("RIAGENDR"), svy.Cat("RIDRETH3"), "BMXBMI"],
    family="binomial",
)

Gives Taylor-linearised SEs, design df, t-based CIs, and an adjusted Wald F. fit(..., where=...) does correct domain estimation matching R’s svyglm(..., design = subset(d, ...)).

Odds ratios: coefficients are on the log-odds scale. As of 0.24.1 there is no exponentiate= flag, so exponentiate manually:

import math
for c in m.coefs:
    print(c.term, math.exp(c.est), math.exp(c.lci), math.exp(c.uci))

Exponentiate the CI endpoints — not a delta-method SE on the ratio scale. The t statistic and p-value are unchanged (β = 0 ⟺ OR = 1).

Also demonstrate m.margins() — predicted probabilities and marginal effects, which base R’s svyglm does not provide natively. For a prevalence paper these are often the better thing to report than ORs.

7. Validation

Pin the numbers against published NCHS estimates rather than against another library. CDC’s Sample Code module (https://wwwn.cdc.gov/nchs/nhanes/tutorials/samplecode.aspx) replicates seven NCHS publications in SAS, SUDAAN, Stata and R — including hypertension and prescription medication use, both of which carry published ORs and Korn–Graubard proportion tables.

Suggested: two-tier tolerance. Published tables are rounded, so compare at published precision for correctness, and separately pin a full-precision reference generated once from CDC’s own code for regression detection.

8. Gaps to flag honestly in the write-up

Gap Status
Odds ratios Manual exp(); no exponentiate= flag yet
Age standardization / age-adjusted prevalence Not implemented — many NCHS publications report age-adjusted rates
Multi-cycle pooling Deliberately out of scope; no helper exists
Probit / cloglog links Not implemented
Joins Use polars; svy has no join

Contrast with the DHS Ethiopia case study

Worth a short closing section — the two exercise genuinely different paths:

NHANES DHS
Design Stratified multistage Two-stage cluster
Weights Inflation (WTMEC2YR) Normalized to n (v005/1e6)
deff reference wor valid wor raises by construction — must use wr
Format XPT Stata .dta / SPSS .sav
Validation oracle Seven NCHS publications Appendix B tables in every final report
Back to top