svy on NHANES, August 2021–August 2023

Health

United States

Validation

Chronic-condition prevalence and multimorbidity from NHANES public microdata with svy: reproducing two NCHS data briefs at published precision including their Korn-Graubard intervals, agreeing with R’s survey package to machine precision, and working through the weight, skip-pattern and age-standardization decisions the cycle demands.

Author

Mamadou S. Diallo, Ph.D.

Published

August 30, 2026

Modified

August 31, 2026

Keywords

NHANES, NCHS, complex survey analysis Python, survey weighted logistic regression, Korn-Graubard confidence interval, XPT SAS transport, comorbidity multimorbidity, design-based estimation, chronic disease prevalence

Summary

TL;DR

NHANES is the reference dataset for American chronic-disease prevalence, and almost all of the difficulty in using it is upstream of the estimator. The design is a plain stratified two-stage sample. What is not plain is everything about how the eight condition indicators get built.

Working through it on the August 2021–August 2023 cycle, svy reproduces the published NCHS estimates at published precision, including the Korn–Graubard confidence intervals:

NCHS Data Brief 508 / 511	svy	Published
Obesity, adults 20+	40.26 (se 1.84)	40.3 (se 1.8)
Obesity, men	39.20 (se 1.93)	39.2 (se 1.9)
Obesity, women	41.28 (se 2.20)	41.3 (se 2.2)
Severe obesity, adults 20+	9.41 (se 0.66)	9.4 (se 0.7)
Hypertension, adults 18+	47.70 [45.4, 50.0]	47.7 [45.4, 49.9]
Hypertension, men	50.87 [48.5, 53.2]	50.8 [48.4, 53.2]
Awareness, among hypertensive	59.17 [56.4, 61.9]	59.2 [56.4, 61.9]
Treatment, among hypertensive	51.19 [48.6, 53.7]	51.2 [48.6, 53.7]

Notethat the svy estimators have been validated separately and in depth against R’s survey package — see svy vs R’s survey.

Background

NHANES combines a household interview with a standardized physical examination in a mobile examination centre (MEC). This cycle interviewed 11,933 people and examined 8,860 of them.

The design is the straightforward part:

View code

import svy

DESIGN = svy.Design(stratum="SDMVSTRA", psu="SDMVPSU", wgt="WTMEC2YR")
print(DESIGN)

╭────────── Design ───────────╮
│ Row index          None     │
│ Stratum            SDMVSTRA │
│ PSU                SDMVPSU  │
│ SSU                None     │
│ Weight             WTMEC2YR │
│ With replacement   False    │
│ Prob               None     │
│ Hit                None     │
│ MOS                None     │
│ Population size    None     │
│ Replicate weights  None     │
╰─────────────────────────────╯

Fifteen masked variance strata, two masked PSUs in each, so 30 PSUs and 15 degrees of freedom; regardless of the six thousand people in the analysis.

What the files actually carry

svy.read_sas dispatches .xpt to a native Rust reader, so there is no conversion step. create_from_sas goes one further and returns a Sample with the file’s metadata already imported:

View code

import svy

demo = svy.read_sas("data/DEMO_L.xpt")
print(type(demo).__name__, demo.shape)

demo_sample = svy.create_from_sas("data/DEMO_L.xpt")
print(demo_sample.n_records, "records, labels imported:")
print(" ", demo_sample.meta.resolve_labels("RIDAGEYR").var_label)
print(" ", demo_sample.meta.resolve_labels("WTMEC2YR").var_label)

DataFrame (11933, 27)
11933 records, labels imported:
  Age in years at screening
  Full sample 2-year MEC exam weight

Reading the metadata from all eight modules gives a uniform and consequential result:

Module	Variables	With a label
DEMO_L	27	27
BMX_L	22	22
BPXO_L	12	12
BPQ_L	6	6
DIQ_L	9	9
GHB_L	3	3
KIQ_U_L	9	9
MCQ_L	35	35

Variable labels survive. All 123 variables across the eight modules carry one, as the two above show.

Value labels do not, because they were never there. A SAS transport file stores a format name per variable, never the format’s contents; the mapping itself lives in a SAS catalogue that NHANES does not ship. And NHANES sets a format name on exactly two of the 123 variables — WTINT2YR and WTMEC2YR, both 20.6, a numeric width-and-decimals format for printing the weights. Not one variable carries a format that would map 1 to “Yes” or 9 to “Don’t know”. That mapping exists only in the HTML codebooks on CDC’s website.

The user_missing and tagged_missings blocks are empty too, and for a separate reason: NHANES encodes refusals and don’t-knows as ordinary numbers rather than as SAS special missing values. There is nothing in the file to import.

So MCQ160C == 9 — a respondent who does not know whether they have coronary heart disease — arrives as the number 9, and there are 34 of them. Write svy.col("MCQ160C") == 1 and every one is counted as not having the disease.

svy will not rescue you here, and says so. Its test suite asserts the absence of a missing-value model:

svy has no model of missing values. A 99 labelled “Refusal” is the integer 99 with the label “Refusal”. svy reads it, prints it, and forms no opinion about it. Absence is a polars null, which needs no metadata.

Labelling a code changes what it is called and nothing else. There is no set_missing to reach for, and MissingKind in svy.core.enumerations is vestigial. The recode is the analyst’s job, which is what the next section is about.

One thing the reader does handle for you

The 449 people who provided no blood sample have WTPH2YR set to zero. In the XPT file that zero is encoded as an IBM-format denormal rather than a true zero, and svy normalises it on read:

View code

w = svy.read_sas("data/GHB_L.xpt")["WTPH2YR"]
print("min:", w.min(), "| == 0:", (w == 0).sum(), "| > 0:", (w > 0).sum())

min: 0.0 | == 0: 449 | > 0: 6750

So the obvious filter works: WTPH2YR > 0 selects exactly the people with a usable phlebotomy weight, and == 0 finds all 449 without one. Worth knowing because that is not guaranteed — a reader that passes the denormal through would make == 0 match nothing at all, and silently admit 449 zero-weight records into the design.

The recode, and what each step is worth

The design goes on first, and everything after it runs through Sample.wrangling, so the design travels with the data rather than being bolted on at the end:

View code

sample = svy.Sample(data=df, design=DESIGN)

Two operations then stand between the raw file and a usable indicator.

Sentinels. Every questionnaire item uses 7 = Refused and 9 = Don’t know. They have to become null before any indicator is built. recode maps new values to the old ones they replace, so this is one call over all sixteen items:

View code

sample = sample.wrangling.recode(
    ITEMS,  # all sixteen questionnaire items
    {None: [7, 9]},  # Refused and Don't know -> null
    replace=True,  # overwrite in place
)

Skip patterns. BPQ150 (“now taking blood pressure medication”) is asked only of people who answered yes to BPQ020 (“ever told you had high blood pressure”). For everyone else it is null — and that null means not on medication, not unknown. This one is conditional on another column, so it is a mutate:

View code

sample = sample.wrangling.mutate({
    "BPQ150": (
        svy.when(svy.col("BPQ150").is_not_null()).then(svy.col("BPQ150"))
        .when(svy.col("BPQ020") == 2).then(svy.lit(2.0))  # asked, said no -> No
        .otherwise(None)                                   # gate unknown -> unknown
    )
})

The two together are the whole cleaning step, and both return a new Sample, so the pipeline chains.

Running it three times — correct, then with each step disabled — prices them separately. Adults 20+, MEC weight:

Condition	Correct	Skips not closed	Sentinels not zapped
Hypertension	49.12 (se 1.03, n=5,821)	95.13 (se 0.34, n=3,367)	49.10 (se 1.04, n=5,822)
Diabetes (diagnosed)	11.81 (se 0.70, n=6,023)	11.81 (se 0.70, n=6,023)	11.81 (se 0.70, n=6,023)
Obesity	40.26 (se 1.84, n=5,929)	40.26 (se 1.84, n=5,929)	40.26 (se 1.84, n=5,929)
High cholesterol	38.61 (se 0.91, n=5,985)	38.61 (se 0.91, n=5,985)	38.37 (se 0.90, n=6,023)
Cardiovascular disease	9.73 (se 0.51, n=5,995)	9.73 (se 0.51, n=5,995)	9.70 (se 0.51, n=6,021)
Chronic lung disease	5.74 (se 0.60, n=6,013)	5.74 (se 0.60, n=6,013)	5.73 (se 0.60, n=6,023)
Kidney disease	2.95 (se 0.33, n=6,012)	2.95 (se 0.33, n=6,012)	2.95 (se 0.33, n=6,023)
Cancer	11.70 (se 0.53, n=6,017)	11.70 (se 0.53, n=6,017)	11.69 (se 0.53, n=6,023)

The two failures are nothing alike.

The skip-pattern error is catastrophic and it disguises itself as precision. Leaving the structural nulls as missing, then dropping them, restricts the analysis to the 3,367 people who had already been told they have high blood pressure — among whom the prevalence of high blood pressure is, unsurprisingly, 95%. The estimate nearly doubles, and the standard error falls from 1.03 to 0.34. Judged on its confidence interval alone, the wrong answer looks like the better one.

The sentinel error is small here, and it is a bias rather than noise. Across 6,023 adults there are only 173 refusals and don’t-knows, so the effect is a quarter of a point at most. But it is always in the same direction — a don’t-know becomes a “no”, so prevalence is always understated — and it silently inflates the denominator, since those records are no longer dropped. In a cycle or a subgroup with more nonresponse, the same mechanism scales up. It is worth doing correctly regardless of what it costs today.

Do not drop records to subset

CDC’s own guidance is blunt about this:

When working with complex survey data such as NHANES, you should never delete or drop records from your analysis dataset before executing your analysis procedures.

Filtering to adults treats the number of adults in each cluster as fixed, when it is itself a random quantity, and understates the variance. In svy the subpopulation is a where= argument and the design stays whole:

View code

adults = (svy.col("RIDAGEYR") >= 20) & ~svy.col("pregnant")
sample.estimation.prop("obesity", where=adults, ci_method="korn-graubard")

Building the conditions

Each chronic condition is conventionally built as an OR across up to three sources:

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

Self-report alone misses undiagnosed disease; a lab threshold alone misses controlled disease. The union is the standard operationalisation.

It is also the part most often ported unexamined from one cycle to the next, and this cycle is where that breaks. Whether a source contributes anything at all depends on its skip pattern, so the recipe has to be tested rather than assumed. Counting the cases each medication source adds beyond self-report alone:

Source	New cases	Why
Insulin (`DIQ050`)	0	Asked only of `DIQ010 == 1`. It is nested inside the source it is meant to supplement, so it can never add anyone.
Oral agents (`DIQ070`)	120	Asked of `DIQ010 ∈ {1,3}` or `DIQ160 == 1`. The extra 120 are 60 borderline and 60 prediabetic respondents on metformin — not people with diabetes.
Statins (`BPQ101D`)	272	Ungated: asked of everyone. These are real cases, 12.8% of all medication users, that self-report misses.

So for diabetes the medication source is either vacuous or actively wrong, and diagnosed diabetes here is DIQ010 == 1 alone. For cholesterol the medication source is essential.

BPQ101D is the trap for anyone porting a 2017–2018 pipeline

In 2017–2018, cholesterol medication (BPQ100D) sat behind a gate, itself behind BPQ080. In 2021–2023 the variable was renamed to BPQ101D and the gate was removed — CDC states that “the skip pattern for cholesterol medication use changed.”

Code that reproduces the old logic as BPQ080 == 1 & BPQ101D == 1 still runs, raises nothing, and drops 272 treated people. The same rename happened to blood-pressure medication: BPQ050A became BPQ150, where the gate was kept.

Two more renames matter for anyone trending against earlier cycles. MCQ160P (“ever told you had COPD, emphysema, ChB”) is one question that replaced three separately-asked items — MCQ160G, MCQ160K and MCQ160O. A union of three items and a single item are not the same measurement and should not be plotted as one series.

Once the definitions do line up, svy.combine_samples() stacks cycles into a single stratified design, nesting wave inside stratum inside PSU so that Taylor variance treats the waves as independent. It picks up an existing wave column such as NHANES’s own SDDSRVYR, and adjust="average" divides the weights by the number of cycles — which matters for totals and not for proportions. Nothing here pools cycles.

The eight conditions used here:

Condition	Definition
Hypertension	mean measured SBP ≥ 130 or mean DBP ≥ 80 or `BPQ150 == 1`
Diabetes (diagnosed)	`DIQ010 == 1`
Obesity	`BMXBMI ≥ 30`
High cholesterol	`BPQ080 == 1` or `BPQ101D == 1`
Cardiovascular disease	any of `MCQ160B`–`MCQ160F` (heart failure, CHD, angina, MI, stroke)
Chronic lung disease	`MCQ160P`
Kidney disease	`KIQ022 == 1`
Cancer	`MCQ220 == 1`

In code, the OR needs to distinguish “no” from “not known”, so a plain any_horizontal will not do — a confirmed Yes has to win over a missing source, while nothing-says-yes-and-something-is-unknown has to stay null:

View code

def any_yes(*sources):
    """OR across sources, propagating nulls only when undecided."""
    return (
        svy.when(svy.any_horizontal(*[s.fill_null(False) for s in sources]))
        .then(True)
        .when(svy.any_horizontal(*[s.is_null() for s in sources]))
        .then(None)
        .otherwise(False)
    )


sample = sample.wrangling.mutate(
    {
        "high_chol": any_yes(yes("BPQ080"), yes("BPQ101D")),
        "cvd": any_yes(
            yes("MCQ160B"),
            yes("MCQ160C"),
            yes("MCQ160D"),
            yes("MCQ160E"),
            yes("MCQ160F"),
        ),
        "obesity": at_least("BMXBMI", 30),
        "kidney": yes("KIQ022"),
    }
)

Counting conditions per person needs the same care. A respondent with two confirmed conditions and one unknown does not have “exactly two”, so the count is null unless every flag is resolved:

View code

CONDITIONS = ["hypertension", "diabetes", "obesity", "high_chol",
              "cvd", "lung", "kidney", "cancer"]

flags = [svy.col(c) for c in CONDITIONS]
n = svy.sum_horizontal(*flags)                 # booleans sum without a cast
complete = svy.all_horizontal(*[f.is_not_null() for f in flags])

sample = sample.wrangling.mutate({
    "n_conditions": svy.when(complete).then(n).otherwise(None),
    "multimorbid": svy.when(complete).then(n >= 2).otherwise(None),
})

That costs 5.7% of adults, which is worth reporting rather than hiding. It has to be a second mutate — like with_columns, one call evaluates every expression against the frame as it was on entry, so a column created in the same call is not visible to its siblings.

The demographic groupings are categorize and recode rather than hand-written when chains:

View code

sample = (
    sample.wrangling.categorize(
        "RIDAGEYR",
        bins=[-float("inf"), 39, 59, float("inf")],
        labels=["20-39", "40-59", "60+"],
        into="age_group",
    )
    .wrangling.recode("RIAGENDR", {"Men": [1], "Women": [2]}, into="sex")
    .wrangling.recode(
        "RIDRETH3",
        {
            "Mexican American": [1],
            "Other Hispanic": [2],
            "White, non-Hispanic": [3],
            "Black, non-Hispanic": [4],
            "Asian, non-Hispanic": [6],
            "Other/multiracial": [7],
        },
        into="race_eth",
    )
)

Note RIDRETH3 has no code 5 — it is reserved in RIDRETH1 and unused here. Writing the map out explicitly rather than looping over a range is what makes that visible.

Three judgement calls are worth stating rather than burying.

The blood-pressure threshold is a choice. 130/80 is the 2017 ACC/AHA definition and is what NCHS uses for this cycle; the older 140/90 gives a substantially lower prevalence.

A measured blood pressure is required, not merely used when present. Without one, “not on medication” cannot be turned into “not hypertensive”. This also keeps interview-only records out of the denominator, where the medication item alone would otherwise admit 514 people at zero weight — contributing nothing to the estimate while inflating the reported sample size.

DIQ010 == 3 is “Borderline”, a real third category and not a missing code. It is resolved to “no” here, which is what NCHS does for diagnosed diabetes, but folding it into “yes” instead is a defensible choice that changes the answer.

The weight must match the most restrictive source

CDC’s rule is explicit:

You must use the weight of the smallest subpopulation that includes all the variables you want to include in your analysis.

This cycle introduced a new weight that makes the rule bite. HbA1c (LBXGH) is attempted on all MEC examinees aged 12+, but roughly 5% of examined adults gave no blood, and nonresponse to the blood draw differs sharply by age and race. So NCHS shipped WTPH2YR, a phlebotomy weight, in every lab file:

The phlebotomy weight should be used for analyses that use variables derived from blood analytes.

The arithmetic shows exactly what it is for. Adults 20+:

	Population represented
`WTMEC2YR` over the MEC sample	245,061,169
`WTPH2YR` over the phlebotomy subset	245,061,169
`WTMEC2YR` over the phlebotomy subset	233,934,098

The phlebotomy weight re-inflates a smaller set of people to the same national total. The MEC weight, applied to those same people, represents 11.1 million fewer Americans — and reports no error while doing it.

On the point estimate the effect is modest: total diabetes (diagnosed or HbA1c ≥ 6.5%) comes to 14.09% under WTPH2YR and 13.98% under WTMEC2YR on identical records, well inside a standard error of 0.75. The damage is not usually in the first estimate. It is that the analysis silently stops representing the population it claims to.

Switching weight is switching Design, and the frame has to move with it — the phlebotomy subsample is the records with a positive WTPH2YR:

View code

phleb = sample.wrangling.filter_records(svy.col("WTPH2YR") > 0).update_design(wgt="WTPH2YR")
print(
    phleb.estimation.prop(
        "diabetes_total",
        where=adults,
        ci_method="korn-graubard",
        drop_nulls=True,
    ),
)

╭─────────────────── Estimate: PROP (TAYLOR) ───────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())                  │
│                                                               │
│                                                               │
│  diabetes_total      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  false            0.8591   0.0075   0.8424   0.8747     0.87  │
│  true             0.1409   0.0075   0.1253   0.1576     5.30  │
│                                                               │
╰───────────────────────────────────────────────────────────────╯

The rule cuts the other way too. Diagnosed diabetes rests on DIQ010 alone, a household-interview item, so the interview weight is the correct one:

Diagnosed diabetes, adults 20+	Estimate
`WTINT2YR`, all interviewed adults	11.23 (se 0.67)
`WTMEC2YR`, MEC examinees only	11.81 (se 0.70)
NCHS Data Brief 516, published	11.3 [9.3, 13.5]

Reaching for WTMEC2YR because the rest of the analysis is MEC-based moves the estimate by 0.6 points, away from the published figure. (Brief 516 estimates this on the fasting subsample under yet another weight, so this is corroboration rather than a replication.)

Prevalence

The sample is already built and cleaned, so estimation is one call:

View code

adults = (svy.col("RIDAGEYR") >= 20) & ~svy.col("pregnant")

print(
    sample.estimation.prop(
        "obesity",
        where=adults,
        ci_method="korn-graubard",
        drop_nulls=True,
    )
)

╭─────────────── Estimate: PROP (TAYLOR) ────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())           │
│                                                        │
│                                                        │
│  obesity      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  false     0.5974   0.0184   0.5568   0.6369     3.09  │
│  true      0.4026   0.0184   0.3631   0.4432     4.58  │
│                                                        │
╰────────────────────────────────────────────────────────╯

y also takes a sequence, so all eight conditions come back from one call as an EstimateList:

View code

est = sample.estimation.prop(
    CONDITIONS,
    where=adults,
    ci_method="korn-graubard",
    drop_nulls=True,
)
print(est)

╭────────────────────── Estimate: PROP (TAYLOR) ──────────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())                        │
│                                                                     │
│                                                                     │
│  y              level      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  hypertension   false   0.5088   0.0103   0.4866   0.5310     2.03  │
│  hypertension   true    0.4912   0.0103   0.4690   0.5134     2.10  │
│  diabetes       false   0.8819   0.0070   0.8662   0.8964     0.79  │
│  diabetes       true    0.1181   0.0070   0.1036   0.1338     5.90  │
│  obesity        false   0.5974   0.0184   0.5568   0.6369     3.09  │
│  obesity        true    0.4026   0.0184   0.3631   0.4432     4.58  │
│  high_chol      false   0.6139   0.0091   0.5942   0.6333     1.48  │
│  high_chol      true    0.3861   0.0091   0.3667   0.4058     2.36  │
│  cvd            false   0.9027   0.0051   0.8913   0.9133     0.56  │
│  cvd            true    0.0973   0.0051   0.0867   0.1087     5.22  │
│  lung           false   0.9426   0.0060   0.9283   0.9548     0.64  │
│  lung           true    0.0574   0.0060   0.0452   0.0717    10.48  │
│  kidney         false   0.9705   0.0033   0.9625   0.9771     0.34  │
│  kidney         true    0.0295   0.0033   0.0229   0.0375    11.21  │
│  cancer         false   0.8830   0.0053   0.8712   0.8941     0.60  │
│  cancer         true    0.1170   0.0053   0.1059   0.1288     4.53  │
│                                                                     │
╰─────────────────────────────────────────────────────────────────────╯

Adults 20+, not pregnant, MEC weight, Korn–Graubard intervals. Every interval carries 15 degrees of freedom.

by= estimates every subgroup in one pass, and to_polars() turns any result into a frame:

View code

print(
    sample.estimation.prop(
        "obesity",
        by="race_eth",
        where=adults,
        ci_method="korn-graubard",
        drop_nulls=True,
    )
)

╭────────────────────────── Estimate: PROP (TAYLOR) ───────────────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())                                 │
│                                                                              │
│                                                                              │
│  race_eth              obesity      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  Asian, non-Hispanic   false     0.8663   0.0331   0.7774   0.9292     3.82  │
│  Asian, non-Hispanic   true      0.1337   0.0331   0.0708   0.2226    24.75  │
│  Black, non-Hispanic   false     0.4831   0.0246   0.4293   0.5371     5.09  │
│  Black, non-Hispanic   true      0.5169   0.0246   0.4629   0.5707     4.76  │
│  Mexican American      false     0.5243   0.0494   0.4134   0.6334     9.43  │
│  Mexican American      true      0.4757   0.0494   0.3666   0.5866    10.39  │
│  Other Hispanic        false     0.5909   0.0308   0.5221   0.6572     5.22  │
│  Other Hispanic        true      0.4091   0.0308   0.3428   0.4779     7.53  │
│  Other/multiracial     false     0.5531   0.0451   0.4520   0.6512     8.16  │
│  Other/multiracial     true      0.4469   0.0451   0.3488   0.5480    10.10  │
│  White, non-Hispanic   false     0.6044   0.0174   0.5660   0.6418     2.89  │
│  White, non-Hispanic   true      0.3956   0.0174   0.3582   0.4340     4.41  │
│                                                                              │
╰──────────────────────────────────────────────────────────────────────────────╯

Note the df column: three domains report 14 rather than 15, because one stratum contributes no observations to them. svy adjusts the reference distribution per domain rather than reusing the full-sample degrees of freedom.

Which interval method

ci_method takes logit (the default, matching Stata’s svy: prop), beta (matching R’s svyciprop(method = "beta")), korn-graubard (matching the NCHS SAS macro) and wilson.

beta and korn-graubard are both Korn–Graubard and differ in two places, both of which only bite in small domains and on rare conditions. korn-graubard truncates the effective sample size at the number of people actually observed; beta leaves it uncapped, because R does. In 12 of the 72 condition-by-domain cells here the uncapped n_eff exceeds the number of people observed — for chronic lung disease among Asian adults it reaches 854 from 335 respondents, and the interval comes out 26% narrower than the capped one.

The second difference is sharper. Among Asian adults aged 20–39, no one out of 102 reported kidney disease:

Method	Estimate	95% CI
logit	0.00	[0.00, 0.0000]
beta	0.00	[0.00, 0.0000]
wilson	0.00	[0.00, 0.0000]
korn-graubard	0.00	[0.00, 0.0655]

Three of the four report a zero-width 95% interval — the claim that no Asian adult under 40 in the United States has kidney disease, from 102 observations. Only korn-graubard returns a proper one-sided interval.

For beta this is correct behaviour: it is faithfully reproducing R’s svyciprop(method = "beta"), which degenerates at the boundary. korn-graubard follows the NCHS SAS macro, which adds explicit p = 0 and p = 1 handling. The wilson case is a genuine limitation — the Wilson interval is normally well-behaved at zero, but the design-based effective sample size is p(1-p)/se², which is 0/0 there, and svy short-circuits to a point.

For rare conditions in small domains, which is most subgroup work on this data, use korn-graubard. It is also what NCHS uses, which is the other reason the validation below lands.

Multimorbidity

tabulate gives the whole distribution of the condition count with design-based standard errors. A second positional argument would make it a cross-tab against another variable:

View code

print(sample.categorical.tabulate("n_conditions", units="percent",
                                  drop_nulls=True))

╭────────────────── Table: n_conditions ──────────────────╮
│                                                         │
│  Row   Estimate   Std Err       CV     Lower     Upper  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  0      25.7296    1.4773   0.0574   22.7084   29.0019  │
│  1      25.1980    0.8040   0.0319   23.5231   26.9500  │
│  2      22.8175    0.8716   0.0382   21.0131   24.7283  │
│  3      14.2637    0.6027   0.0423   13.0266   15.5972  │
│  4       7.1989    0.3268   0.0454    6.5327    7.9273  │
│  5       3.3623    0.3157   0.0939    2.7506    4.1042  │
│  6       1.1661    0.1248   0.1070    0.9279    1.4645  │
│  7       0.2592    0.0748   0.2888    0.1400    0.4794  │
│  8       0.0048    0.0048   0.9965    0.0006    0.0401  │
│                                                         │
╰─────────────────────────────────────────────────────────╯

One row variable at a time

estimation.mean and estimation.prop take y as either a name or a sequence of names, so eight conditions come back in one call. categorical.tabulate takes rowvar as a single string, so tabulating each of the eight means eight calls and stacking the results yourself. Not wrong, just an asymmetry worth knowing before you write the loop.

multimorbid — the two-or-more flag built during the recode — is an ordinary binary outcome, so it goes through prop like any other. The overall figure and each breakdown are separate calls:

View code

print(
    sample.estimation.prop(
        "multimorbid",
        where=adults,
        ci_method="korn-graubard",
        drop_nulls=True,
    )
)

for domain in ("sex", "age_group", "race_eth"):
    print(
        sample.estimation.prop(
            "multimorbid",
            by=domain,
            where=adults,
            ci_method="korn-graubard",
            drop_nulls=True,
        )
    )

╭───────────────── Estimate: PROP (TAYLOR) ──────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())               │
│                                                            │
│                                                            │
│  multimorbid      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  false         0.5052   0.0142   0.4745   0.5358     2.81  │
│  true          0.4948   0.0142   0.4642   0.5255     2.87  │
│                                                            │
╰────────────────────────────────────────────────────────────╯

╭───────────────────── Estimate: PROP (TAYLOR) ──────────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())                       │
│                                                                    │
│                                                                    │
│  sex     multimorbid      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  Men     false         0.4960   0.0169   0.4594   0.5326     3.41  │
│  Men     true          0.5040   0.0169   0.4674   0.5406     3.35  │
│  Women   false         0.5141   0.0168   0.4776   0.5506     3.27  │
│  Women   true          0.4859   0.0168   0.4494   0.5224     3.46  │
│                                                                    │
╰────────────────────────────────────────────────────────────────────╯

╭─────────────────────── Estimate: PROP (TAYLOR) ────────────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())                           │
│                                                                        │
│                                                                        │
│  age_group   multimorbid      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  20-39       false         0.7722   0.0174   0.7325   0.8085     2.25  │
│  20-39       true          0.2278   0.0174   0.1915   0.2675     7.64  │
│  40-59       false         0.4730   0.0157   0.4391   0.5071     3.32  │
│  40-59       true          0.5270   0.0157   0.4929   0.5609     2.98  │
│  60+         false         0.2430   0.0160   0.2093   0.2791     6.60  │
│  60+         true          0.7570   0.0160   0.7209   0.7907     2.12  │
│                                                                        │
╰────────────────────────────────────────────────────────────────────────╯

╭──────────────────────────── Estimate: PROP (TAYLOR) ─────────────────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())                                     │
│                                                                                  │
│                                                                                  │
│  race_eth              multimorbid      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  Asian, non-Hispanic   false         0.6913   0.0340   0.6123   0.7628     4.92  │
│  Asian, non-Hispanic   true          0.3087   0.0340   0.2372   0.3877    11.01  │
│  Black, non-Hispanic   false         0.4319   0.0165   0.3943   0.4702     3.81  │
│  Black, non-Hispanic   true          0.5681   0.0165   0.5298   0.6057     2.90  │
│  Mexican American      false         0.5995   0.0388   0.5115   0.6831     6.47  │
│  Mexican American      true          0.4005   0.0388   0.3169   0.4885     9.68  │
│  Other Hispanic        false         0.5707   0.0377   0.4863   0.6521     6.61  │
│  Other Hispanic        true          0.4293   0.0377   0.3479   0.5137     8.79  │
│  Other/multiracial     false         0.4996   0.0283   0.4379   0.5613     5.66  │
│  Other/multiracial     true          0.5004   0.0283   0.4387   0.5621     5.65  │
│  White, non-Hispanic   false         0.4787   0.0156   0.4450   0.5126     3.26  │
│  White, non-Hispanic   true          0.5213   0.0156   0.4874   0.5550     3.00  │
│                                                                                  │
╰──────────────────────────────────────────────────────────────────────────────────╯

Stacking those four results gives the table below — two or more conditions, adults 20+:

Domain	%	SE	95% CI	df
All adults	49.48	1.42	[46.42, 52.55]	15
Men	50.40	1.69	[46.74, 54.06]	15
Women	48.59	1.68	[44.94, 52.24]	15
Age 20-39	22.78	1.74	[19.15, 26.75]	15
Age 40-59	52.70	1.57	[49.29, 56.09]	15
Age 60+	75.70	1.60	[72.09, 79.07]	15
Asian, non-Hispanic	30.87	3.40	[23.72, 38.77]	14
Black, non-Hispanic	56.81	1.65	[52.98, 60.57]	14
Mexican American	40.05	3.88	[31.69, 48.85]	14
Other Hispanic	42.93	3.77	[34.79, 51.37]	15
Other/multiracial	50.04	2.83	[43.87, 56.21]	15
White, non-Hispanic	52.13	1.56	[48.74, 55.50]	15

Half of American adults carry at least two of these eight conditions, and three quarters of those over 60 do.

Note the degrees of freedom. Three domains report 14 rather than 15, because one stratum contributes no observations to them; svy adjusts the reference distribution per domain rather than reusing the full-sample df. The width of the Asian and Mexican American intervals — more than twice the White interval — is the reduced-oversampling caveat from the top of this page, showing up exactly where CDC said it would.

Regression

View code

X_DIABETES = [
    "RIDAGEYR",
    svy.Cat("sex", ref="Men"),
    svy.Cat("race_eth", ref="White, non-Hispanic"),
    "BMXBMI",
]

m = sample.glm.fit("diabetes", x=X_DIABETES, family="binomial", where=adults)
print(m)

╭───────────────────────────────────── GLM: Binomial (logit) ─────────────────────────────────────╮
│ Modeling: diabetes                                                                              │
│                                                                                                 │
│ Observations       5929  AIC           3667.6193                                                │
│ DF Residuals          7  BIC                   -                                                │
│ Deviance      3648.5231  Scale            1.0000                                                │
│ R-squared       0.14498  R-sq (adj)      0.14383                                                │
│                          Iterations            6                                                │
│ F-stat (adj)   74.38227  Prob (F-adj)     <0.001                                                │
│                                                                                                 │
│                                                                                                 │
│  Term                              Coef.   Std.Err.           t    P>|t|     [0.025     0.975]  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  _intercept_                    -7.69223    0.37550   -20.48529   <0.001   -8.58014   -6.80431  │
│  RIDAGEYR                        0.05632    0.00383    14.71298   <0.001    0.04727    0.06537  │
│  sex_Women                      -0.30559    0.08475    -3.60597   0.0087   -0.50599   -0.10520  │
│  race_eth_Asian, non-Hispanic    0.63179    0.21924     2.88169   0.0236    0.11336    1.15021  │
│  race_eth_Black, non-Hispanic    0.77659    0.15543     4.99648   0.0016    0.40906    1.14412  │
│  race_eth_Mexican American       0.68245    0.16193     4.21441   0.0040    0.29954    1.06537  │
│  race_eth_Other Hispanic         0.45821    0.20336     2.25316   0.0589   -0.02267    0.93909  │
│  race_eth_Other/multiracial      0.67436    0.17393     3.87710   0.0061    0.26307    1.08565  │
│  BMXBMI                          0.07925    0.00567    13.98860   <0.001    0.06585    0.09265  │
│                                                                                                 │
╰─────────────────────────────────────────────────────────────────────────────────────────────────╯

That gives Taylor-linearised standard errors, t-based intervals on the design df, and an adjusted Wald F, all on the log-odds scale. exponentiate=True reports the same fit as odds ratios, with White non-Hispanic and Men as reference:

View code

m.fitted.show(exponentiate=True)

╭───────────────────────────────────── GLM: Binomial (logit) ─────────────────────────────────────╮
│ Modeling: diabetes                                                                              │
│                                                                                                 │
│ Observations       5929  AIC           3667.6193                                                │
│ DF Residuals          7  BIC                   -                                                │
│ Deviance      3648.5231  Scale            1.0000                                                │
│ R-squared       0.14498  R-sq (adj)      0.14383                                                │
│                          Iterations            6                                                │
│ F-stat (adj)   74.38227  Prob (F-adj)     <0.001                                                │
│                                                                                                 │
│                                                                                                 │
│  Term                           Odds ratio   Std.Err.           t    P>|t|    [0.025    0.975]  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  _intercept_                       0.00046    0.37550   -20.48529   <0.001   0.00019   0.00111  │
│  RIDAGEYR                          1.05793    0.00383    14.71298   <0.001   1.04840   1.06755  │
│  sex_Women                         0.73669    0.08475    -3.60597   0.0087   0.60291   0.90014  │
│  race_eth_Asian, non-Hispanic      1.88097    0.21924     2.88169   0.0236   1.12004   3.15885  │
│  race_eth_Black, non-Hispanic      2.17404    0.15543     4.99648   0.0016   1.50540   3.13966  │
│  race_eth_Mexican American         1.97873    0.16193     4.21441   0.0040   1.34924   2.90190  │
│  race_eth_Other Hispanic           1.58125    0.20336     2.25316   0.0589   0.97759   2.55766  │
│  race_eth_Other/multiracial        1.96278    0.17393     3.87710   0.0061   1.30092   2.96136  │
│  BMXBMI                            1.08248    0.00567    13.98860   <0.001   1.06807   1.09708  │
│                                                                                                 │
│ Std.Err., t and P>|t| are on the link scale; the interval is exp() of the link-scale bounds.    │
╰─────────────────────────────────────────────────────────────────────────────────────────────────╯

Note the footer: the standard error, t and p stay on the link scale, and the interval is the exponentiated link-scale bound rather than a delta-method interval on the ratio scale. The t statistic and p-value are unchanged either way, since \(\beta = 0 \iff \mathrm{OR} = 1\).

The model has 5,929 observations and 7 residual degrees of freedom

DF Residuals: 7. The design supplies 15 degrees of freedom and the model spends 8 of them on parameters. The reference distribution for every t statistic above is \(t_7\), not the normal — which is why the Asian coefficient, at t = 2.88, gives p = 0.024 rather than the 0.004 a normal approximation would report.

Adding one more categorical predictor with a few levels would exhaust what is left. Six thousand respondents do not buy model complexity on a 30-PSU design; this is the single most common way NHANES regressions overclaim.

For a prevalence paper, predicted probabilities usually communicate better than odds ratios, and margins() gives them with design-correct standard errors:

View code

print(m.margins(at={"BMXBMI": [20.0, 25.0, 30.0, 35.0, 40.0]}))

╭────── GLM Margins: BMXBMI (predictive, 95% CI) ──────╮
│                                                      │
│  Value     Margin         SE                 95% CI  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  20.00   0.054976   0.003682   [0.046269, 0.063683]  │
│  25.00   0.077945   0.003927   [0.068658, 0.087232]  │
│  30.00   0.108704   0.004659   [0.097689, 0.119720]  │
│  35.00   0.148596   0.006714   [0.132719, 0.164473]  │
│  40.00   0.198457   0.010438   [0.173775, 0.223138]  │
│                                                      │
╰──────────────────────────────────────────────────────╯

variables= gives the average marginal effect instead — the average change in predicted probability across the sample, which for a categorical predictor is a contrast against its reference level:

View code

for r in m.margins(variables=["sex", "BMXBMI"]):
    print(r)

╭─────────────── GLM Margins: sex (ame, 95% CI) ────────────────╮
│                                                               │
│     Contrast      Margin         SE                   95% CI  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  Women - Men   -0.027832   0.007545   [-0.045673, -0.009990]  │
│                                                               │
╰───────────────────────────────────────────────────────────────╯

╭───── GLM Margins: BMXBMI (ame, 95% CI) ──────╮
│                                              │
│    Margin         SE                 95% CI  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  0.007208   0.000560   [0.005884, 0.008532]  │
│                                              │
╰──────────────────────────────────────────────╯

Women are about 2.8 percentage points less likely to have diagnosed diabetes than men, adjusting for age, BMI, and race and ethnicity. That is a sentence a reader can act on; “odds ratio 0.737” is not.

The link matters less than it looks

link= takes probit and cloglog alongside the default logit. Their coefficients are not comparable — the BMI coefficient is 0.076 under logit, 0.042 under probit and 0.065 under cloglog, because each lives on its own scale. The marginal effects, which are on the probability scale, very nearly agree:

link	AME, women vs men	AME, per BMI unit
`logit`	-0.02783	+0.00721
`probit`	-0.02803	+0.00733
`cloglog`	-0.02673	+0.00678

The spread across links is smaller than the standard error on any one of them, which is the practical argument for reporting margins rather than coefficients: the marginal effect is a property of the fitted probabilities, so it barely notices which link produced them.

Validation

Four NCHS data briefs cover this cycle. They are the stronger test: they exercise the whole pipeline — definitions, weight choice and estimator together — not just the variance formula. All exclude pregnant women, and the widely-quoted headline figures are the crude ones reproduced here.

Obesity, adults 20+ (Brief 508)

	svy	Published
Total	40.26 (se 1.84)	40.3 (se 1.8)
Men	39.20 (se 1.93)	39.2 (se 1.9)
Women	41.28 (se 2.20)	41.3 (se 2.2)
Age 20-39	35.53	35.5
Age 40-59	46.38	46.4
Age 60+	38.90	38.9
Severe obesity (BMI ≥ 40)	9.41 (se 0.66)	9.4 (se 0.7)

Every published figure at published precision, standard errors included, on an analytic sample of 5,929 — the brief’s stated n exactly.

Hypertension, adults 18+ (Brief 511)

	svy	Published
Total	47.70 [45.4, 50.0]	47.7 [45.4, 49.9]
Men	50.87 [48.5, 53.2]	50.8 [48.4, 53.2]
Women	44.61 [41.5, 47.8]	44.6 [41.4, 47.8]
Aware, among hypertensive	59.17 [56.4, 61.9]	59.2 [56.4, 61.9]
Treated, among hypertensive	51.19 [48.6, 53.7]	51.2 [48.6, 53.7]

The confidence intervals match to the last published digit. NCHS computes these with its own Korn–Graubard SAS macro; ci_method="korn-graubard" reproduces it, effective-sample-size truncation included, on a 15-df design.

The estimators are separately validated against R’s survey package — see svy vs R’s survey.

Age standardization

Every NCHS brief leads with age-adjusted rates, so this is not optional for public-health work. svy has it as weighting.standardize.

A directly standardized proportion is a fixed-weight combination of the age-specific ones,

\[p_{\text{adj}} = \sum_i w_i \, p_i, \qquad w_i = \frac{N_i}{\sum_j N_j}\]

where \(N_i\) is the year-2000 U.S. standard population in age group \(i\). NCHS collapses that standard to three groups for these briefs — 20–39, 40–59 and 60 and over.

The part that is easy to get wrong is the variance. Because the age-specific estimates come from the same sample, they are correlated, so the variance needs the full covariance matrix and not just the diagonal:

\[\operatorname{Var}(p_{\text{adj}}) = \sum_i \sum_j w_i w_j \operatorname{Cov}(p_i, p_j)\]

weighting.standardize takes the reweighting route: reshape the sample so its age composition matches the standard, then estimate as usual.

View code

STD_2000 = {"20-39": 77_670_618, "40-59": 72_816_615, "60+": 45_363_752}

std = sample.weighting.standardize(
    "age_group",  # the composition axis
    shares=STD_2000,  # the standard population
    where=adults & svy.col("obesity").is_not_null(),
)

print(
    std.estimation.prop(
        "obesity",
        where=adults,
        ci_method="korn-graubard",
        drop_nulls=True,
    ),
)

╭─────────────── Estimate: PROP (TAYLOR) ────────────────╮
│ where: ([RIDAGEYR >= 20]) & (pregnant.not())           │
│                                                        │
│                                                        │
│  obesity      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  false     0.5966   0.0188   0.5552   0.6370     3.15  │
│  true      0.4034   0.0188   0.3630   0.4448     4.66  │
│                                                        │
╰────────────────────────────────────────────────────────╯

std_pop <- data.frame(age_group = names(STD_2000), Freq = STD_2000)

ds <- svystandardize(subset(adults, !is.na(obesity)),
                     by = ~age_group,             # the composition axis
                     over = ~1,                   # domains to standardize within
                     population = std_pop)

svymean(~as.numeric(obesity), ds, na.rm = TRUE)

R	svy
`by =`	`cells`	the composition axis
`population =`	`shares=`	the standard population
`over =`	`by=`	domains to standardize within
`excluding.missing =`	`where=`	the scope

If you are coming from survey, the R panel above is the same call; the arguments map one to one, and the two agree to round-off on all nine indicators here.

Standardized weights belong to one analysis

where bakes one variable’s missingness into the weights, and by bakes in the domain structure. Reusing a standardized sample to estimate a different variable, or a different breakdown, is silently wrong rather than an error. Call standardize once per estimate — which is why the code above passes svy.col("obesity").is_not_null() rather than standardizing once and reusing.

Crude against age-standardized, adults 20+:

Condition	Crude	SE	Age-standardized	SE	Shift
Hypertension	49.12	1.03	46.05	1.13	-3.06
Obesity	40.26	1.84	40.34	1.88	+0.08
High cholesterol	38.61	0.91	35.38	0.87	-3.24
Diabetes (diagnosed)	11.81	0.70	10.62	0.58	-1.19
Cancer	11.70	0.53	9.79	0.34	-1.91
Cardiovascular disease	9.73	0.51	8.23	0.42	-1.50
Severe obesity (BMI ≥ 40)	9.41	0.66	9.73	0.70	+0.32
Chronic lung disease	5.74	0.60	5.12	0.51	-0.62
Kidney disease	2.95	0.33	2.59	0.30	-0.36

Almost everything falls, because this cycle’s sample is older than the 2000 standard — the age-based oversampling of 60+ that replaced the race and income oversampling.

Three independent checks confirm the standard population and the method, since Brief 508 publishes both columns:

	Crude	Age-adjusted
Obesity, adults 20+	40.26 vs 40.3	40.34 vs 40.3
Severe obesity	9.41 vs 9.4	9.73 vs 9.7
Obesity, men	39.20 vs 39.2	39.29 vs 39.3
Obesity, women	41.28 vs 41.3	41.39 vs 41.4

Reproducing the crude-to-adjusted shift — +0.3 for severe obesity, +0.1 for each sex — is the discriminating test, and it lands in all four cases. Adding Brief 511’s hypertension figure on the 18+ base makes five.

Brief 511 uses an 18+ base instead, which needs the standard population’s 15–19 group split to get 18–19. That number is published — 8,001,203, in the master list behind the collapsed table, and confirmed to the person by the SEER single-year standard populations (3,936,904 + 4,064,299). It is worth chasing down rather than approximating: splitting 15–19 pro rata instead moves the answer by a tenth of a point, enough to change the rounded figure.

With the published value, Brief 511’s age-adjusted hypertension comes out at 44.54 against a published 44.5, with a 95% interval of [42.13, 46.96] against a published [42.1, 46.9]:

View code

STD_18PLUS = {"18-39": 85_671_821, "40-59": 72_816_615, "60+": 45_363_752}

adults18 = (svy.col("RIDAGEYR") >= 18) & ~svy.col("pregnant")

std18 = sample.weighting.standardize(
    "ag18",
    shares=STD_18PLUS,
    where=adults18 & svy.col("hypertension").is_not_null(),
)
print(
    std18.estimation.prop(
        "hypertension",
        where=adults18,
        ci_method="korn-graubard",
        drop_nulls=True,
    )
)

╭────────────────── Estimate: PROP (TAYLOR) ──────────────────╮
│ where: ([RIDAGEYR >= 18]) & (pregnant.not())                │
│                                                             │
│                                                             │
│  hypertension      est       se      lci      uci   cv (%)  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  false          0.5546   0.0112   0.5304   0.5787     2.02  │
│  true           0.4454   0.0112   0.4213   0.4696     2.51  │
│                                                             │
╰─────────────────────────────────────────────────────────────╯

The data

Everything above is computed from eight public NHANES modules for the August 2021–August 2023 cycle. They are U.S. government public-domain files, downloadable without registration from https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2021/DataFiles/:

File	Component	Respondents	Variables used here
`DEMO_L`	Demographics and sample weights	11,933	`SDMVSTRA`, `SDMVPSU`, `WTINT2YR`, `WTMEC2YR`, `RIDSTATR`, `RIDAGEYR`, `RIAGENDR`, `RIDRETH3`, `RIDEXPRG`
`BMX_L`	Body measures	8,860	`BMXBMI`, `BMDSTATS`
`BPXO_L`	Blood pressure, oscillometric	7,801	`BPXOSY1`–`BPXOSY3`, `BPXODI1`–`BPXODI3`
`BPQ_L`	Blood pressure and cholesterol questionnaire	8,501	`BPQ020`, `BPQ080`, `BPQ150`, `BPQ101D`
`DIQ_L`	Diabetes questionnaire	11,744	`DIQ010`, `DIQ050`, `DIQ070`
`GHB_L`	Glycohemoglobin	7,199	`LBXGH`, `WTPH2YR`
`KIQ_U_L`	Kidney conditions and urology	7,809	`KIQ022`
`MCQ_L`	Medical conditions	11,744	`MCQ010`, `MCQ160B`–`MCQ160F`, `MCQ160P`, `MCQ220`

The linking key is SEQN throughout, and every module is left-joined onto DEMO_L so that skipping a component does not drop a respondent.

The R cross-check uses survey 4.5 and haven 2.5.5 against the same eight files.

References

Korn, E. L. and Graubard, B. I. (1998). Confidence intervals for proportions with small expected number of positive counts estimated from survey data. Survey Methodology, 24(2), 193–201.

Parker, J. D. et al. (2017). National Center for Health Statistics Data Presentation Standards for Proportions. Vital and Health Statistics 2(175).

Emmerich, S. D., Fryar, C. D., Stierman, B. and Ogden, C. L. (2024). Obesity and severe obesity prevalence in adults: United States, August 2021–August 2023. NCHS Data Brief 508.

Fryar, C. D., Kit, B., Carroll, M. D. and Afful, J. (2024). Hypertension prevalence, awareness, treatment, and control in adults age 18 and older: United States, August 2021–August 2023. NCHS Data Brief 511.

Analyses use public-use NHANES microdata from the National Center for Health Statistics. This does not constitute an endorsement by NCHS or CDC of this product.

Community 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