Skip to contents

NHANES Mortality Linkage: A Complete Workflow

Source: vignettes/nhanes-mortality-workflow.Rmd
nhanes-mortality-workflow.Rmd

Overview

This vignette walks through the complete nhanesR workflow using a concrete example: the association between serum total cholesterol and all-cause mortality across ten NHANES cycles (1999–2018), adjusting for HDL cholesterol, prior myocardial infarction, and cholesterol-lowering medication use.

None of the code chunks run automatically — copy and paste each one into your R console and run it interactively. Downloaded files are cached locally, so re-running any step is fast after the first time.

The workflow has nine steps:

  1. Browse available cycles and files
  2. Discover variables with nhanes_search_variables() and nhanes_variable_map()
  3. Download laboratory data with nhanes_download_analyte()
  4. Download questionnaire data (MI history, cholesterol medication)
  5. Harmonize variable names and units across cycles with nhanes_harmonize()
  6. Recode questionnaire variables
  7. Merge all components with demographics
  8. Link mortality follow-up and prepare the survival dataset
  9. Fit a survey-weighted Cox proportional hazards model

Package options

Three options control nhanesR behavior. The package sets defaults at load time, but any option defined in your .Rprofile before loading takes precedence.

Option Default Purpose
nhanesR.cache_dir file.path(tempdir(), "nhanesR") Root path for all cached RDS files
nhanesR.verbose TRUE Print progress messages during downloads
nhanesR.timeout 120L HTTP timeout in seconds

By default, nhanesR caches files inside R’s session-temporary directory (tempdir()). No files are written to your home directory without your explicit consent. The trade-off is that downloads are repeated in each new R session. To keep a persistent cache, set nhanesR.cache_dir in your ~/.Rprofile:

To make changes permanent, add lines like these to your ~/.Rprofile:

options(
  nhanesR.cache_dir = "/data/nhanes_cache",  # e.g. a shared server path
  nhanesR.verbose   = FALSE,
  nhanesR.timeout   = 300L
)

To check or change settings interactively during a session:

# View current cache location
nhanes_cache_dir()

# Opt in to a persistent home-directory cache for this session
nhanes_cache_dir("~/my_nhanes_cache")

# Suppress download messages for this session
options(nhanesR.verbose = FALSE)

Background: NHANES structure

NHANES (National Health and Nutrition Examination Survey) is conducted in two-year cycles (e.g. 1999–2000, 2001–2002, …, 2017–2018). Within each cycle, data are organized into five components:

Component What it contains
Demographics Age, sex, race/ethnicity, income, survey weights and design variables
Laboratory Blood and urine measurements
Examination Physical exam, blood pressure, anthropometry
Questionnaire Self-reported health history, medications, behaviors
Dietary 24-hour dietary recall interviews

Each participant has a unique identifier, SEQN, that links files within a cycle. SEQNs are not reused across cycles, so always include "cycle" in merge keys when pooling multiple cycles.

1. What cycles and files are available? 

# All continuous NHANES cycles known to nhanesR
nhanes_cycles()

# Just the cycle labels for the first ten continuous cycles (1999-2018)
cycles <- nhanes_cycles()[1:10, "cycle"]
cycles

To see what files are available for a specific cycle and component, use nhanes_manifest():

nhanes_manifest("2015-2016", "Laboratory")
nhanes_manifest("2013-2014", "Questionnaire")

2. Discover variables

NHANES analytes are often stored under different variable names in different cycles. nhanes_search_variables() searches the CDC variable catalog by keyword. nhanes_variable_map() returns a one-row-per-cycle lookup table showing the exact file name and variable name to use.

# Find total cholesterol across all cycles (summarized by default)
nhanes_search_variables("total cholesterol", component = "Laboratory")

# Raw one-row-per-cycle output
nhanes_search_variables("total cholesterol", component = "Laboratory",
                         summarize = FALSE)
# Per-cycle lookup: which file and variable name holds total cholesterol?
nhanes_variable_map("total cholesterol")

# HDL changed variable name three times across cycles
nhanes_variable_map("HDL")

# Questionnaire: history of MI (keep_vars filters out false positives)
nhanes_variable_map("heart attack", component = "Questionnaire",
                     keep_vars = c("MCQ160E", "MCQ160e"))

The nhanes_variable_map() output directly informs the keep_vars argument used in download and harmonization below.

3. Download laboratory data

nhanes_download_analyte() uses the variable catalog to look up the correct CDC file name for each cycle, then downloads it. This resolves cross-cycle file renames automatically — for example, total cholesterol was in LAB13 (1999–2000), L13_B (2001–2002), L13_C (2003–2004), and TCHOL_D onward.

cycles <- nhanes_cycles()[1:10, "cycle"]  # 1999-2018

# Demographics — file name has always been DEMO; nhanes_download() works fine
demo_list <- nhanes_download("DEMO", cycles)

# Total cholesterol — file renamed across early cycles; use download_analyte()
tchol_list <- nhanes_download_analyte("total cholesterol", cycles)

# HDL cholesterol
hdl_list <- nhanes_download_analyte("HDL", cycles)

Files are downloaded in SAS transport (XPT) format, parsed, and cached locally. Subsequent calls load from cache.

4. Download questionnaire data

The same nhanes_download_analyte() function works for any component. Use keep_vars when a search term would otherwise match false positives.

# History of myocardial infarction (MCQ file)
# MCQ160E (1999-2010) and MCQ160e (2011-2018) are the same question;
# keep_vars filters out RXQ510 which also mentions "heart attack"
mi_list <- nhanes_download_analyte(
  "heart attack", cycles,
  component = "Questionnaire",
  keep_vars = c("MCQ160E", "MCQ160e")
)

# Cholesterol-lowering medication (BPQ file)
# "Ever told to take prescribed medicine to lower blood cholesterol?"
chol_med_list <- nhanes_download_analyte(
  "cholesterol", cycles,
  component = "Questionnaire",
  keep_vars = c("BPQ090D", "BPQ101D")
)

5. Harmonize across cycles

nhanes_harmonize() renames per-cycle variables to a single common name and optionally stacks the cycles into one data frame.

Unit-based harmonization (laboratory data): specify unit and name; the function finds the right column in each cycle by matching its label attribute, no variable codes needed. prefer_mgdl = TRUE (default) drops mmol/L duplicates automatically. trim = TRUE (default) returns only SEQN, cycle, and the target column — ready for merging.

# Total cholesterol — LBXTC throughout, but label_pattern narrows the match
# in 1999-2004 when TC and HDL were bundled in the same file
TC <- nhanes_harmonize(
  tchol_list,
  unit          = "mg/dL",
  name          = "TC_mgdl",
  label_pattern = "total cholesterol"
)

# HDL — three different variable names across cycles; unit approach handles all
HDL <- nhanes_harmonize(
  hdl_list,
  unit          = "mg/dL",
  name          = "HDL_mgdl",
  label_pattern = "HDL"
)

str(TC)   # SEQN (chr), cycle (chr), TC_mgdl (num)
str(HDL)  # SEQN (chr), cycle (chr), HDL_mgdl (num)

Mapping-based harmonization (questionnaire data): use mapping when there is no unit to match. The same trim = TRUE default applies.

MI <- nhanes_harmonize(
  mi_list,
  mapping = c(MCQ160E = "MI_history", MCQ160e = "MI_history")
)

chol_med <- nhanes_harmonize(
  chol_med_list,
  mapping = c(BPQ090D = "chol_med", BPQ101D = "chol_med")
)

# Each result is a trim 3-column data frame ready for merging
str(MI)        # SEQN, cycle, MI_history
str(chol_med)  # SEQN, cycle, chol_med

6. Recode questionnaire variables

NHANES questionnaire responses use a numeric coding convention:

Code Meaning
1 Yes
2 No
7 Refused
9 Don’t know

For analysis, recode to 0/1 and treat 7 and 9 as NA:

nhanes_recode_yn <- function(x) {
  out        <- rep(NA_integer_, length(x))
  out[x == 1] <- 1L
  out[x == 2] <- 0L
  out
}

MI$MI_history      <- nhanes_recode_yn(MI$MI_history)
chol_med$chol_med  <- nhanes_recode_yn(chol_med$chol_med)

# Verify: should see 0, 1, and NA only
table(MI$MI_history,      useNA = "always")
table(chol_med$chol_med,  useNA = "always")

7. Stack demographics and merge all components

Stack the per-cycle demographics list, then merge all components by SEQN and cycle. Use all.x = TRUE (left join) from the demographics outward so that participants without lab values are retained with NA.

demo <- nhanes_stack(demo_list)

# Inner join lab data (keeps only participants who attended the exam)
analytic <- Reduce(
  function(a, b) merge(a, b, by = c("SEQN", "cycle")),
  list(demo, TC, HDL)
)

# Left join questionnaire data (all interviewed participants have these)
analytic <- merge(analytic, MI,       by = c("SEQN", "cycle"), all.x = TRUE)
analytic <- merge(analytic, chol_med, by = c("SEQN", "cycle"), all.x = TRUE)

nrow(analytic)
names(analytic)

# Check key variables arrived
c("TC_mgdl", "HDL_mgdl", "MI_history", "chol_med",
  "RIDAGEYR", "RIAGENDR", "WTMEC2YR", "SDMVPSU", "SDMVSTRA") %in%
  names(analytic)

nhanes_mortality_link() downloads the NCHS Public-Use Linked Mortality Files and left-joins them by SEQN. Follow-up runs through December 31, 2019.

analytic_mort <- nhanes_mortality_link(analytic)

# Key variables added:
#   ELIGSTAT     1=eligible, 2=under 18, 3=insufficient data for linkage
#   MORTSTAT     0=assumed alive 31-Dec-2019, 1=assumed deceased
#   UCOD_LEADING Underlying cause of death (11-category ICD-10 recode)
#   PERMTH_EXM   Months from examination date to death or Dec 31 2019
#   PERMTH_INT   Same, from interview date

table(analytic_mort$MORTSTAT, useNA = "always")

nhanes_survival_prep() removes ineligible participants, creates time and event columns, and warns about asymmetric follow-up across cycles. Use origin = "exam" when laboratory measurements are the exposure — they were collected at the exam visit.

surv_data <- nhanes_survival_prep(
  analytic_mort,
  origin     = "exam",
  time_unit  = "years",
  weight_var = "WTMEC2YR"
)

# Follow-up by cycle — note shrinking maximum as cycles approach 2019
nhanes_followup_summary(surv_data)

For cause-specific mortality:

nhanes_ucod_labels()   # see available cause-of-death codes

surv_cvd <- nhanes_survival_prep(
  analytic_mort,
  origin     = "exam",
  time_unit  = "years",
  cause      = "001",   # Diseases of heart
  weight_var = "WTMEC2YR"
)

table(event = surv_cvd$event, cvd_death = surv_cvd$event_cause)

9. Survey-weighted Cox model

NHANES uses a complex multi-stage probability sample. Standard errors must account for the sampling design or they will be anti-conservative.

Choosing the correct survey weight

NHANES provides three families of survey weight. Using the wrong one produces biased population estimates and incorrect standard errors.

Weight Use when
WTINT2YR Interview-only data (questionnaires, no lab or exam)
WTMEC2YR Any examination or laboratory component
WTSAF2YR Analytes from the fasting subsample

The fasting subsample weight (WTSAF2YR) is a statistical probability weight — not a body-weight measurement — that accounts for an additional random subsampling step: only a subset of MEC attendees are asked to fast before their blood draw. Analytes that require fasting include triglycerides, glucose, insulin, and Friedewald-calculated LDL. Using WTMEC2YR for these analytes ignores the fasting subsampling and will over- or under-represent the population.

For total cholesterol and HDL — which do not require fasting — WTMEC2YR is the correct weight.

Pooling across cycles

When combining data from multiple two-year cycles, the 2-year weight must be adjusted. The simplest approach is to divide by the number of cycles pooled:

surv_data$wt_pooled <- surv_data$survey_weight / n_cycles

Some NHANES files include pre-computed 4-year weights (WTMEC4YR, WTSAF4YR). Use these when available rather than dividing manually.

Weight adjustment for pooled cycles: divide the two-year exam weight WTMEC2YR by the number of cycles pooled (here, 10).

nest = TRUE is the correct specification for NHANES — PSU labels may repeat across strata.

surv_data$wt_pooled <- surv_data$survey_weight / 10

# Scale continuous predictors to per-SD units for interpretable hazard ratios
surv_data$TC_sd  <- scale(surv_data$TC_mgdl)[,  1]
surv_data$HDL_sd <- scale(surv_data$HDL_mgdl)[, 1]

design <- svydesign(
  id      = ~SDMVPSU,
  strata  = ~SDMVSTRA,
  weights = ~wt_pooled,
  nest    = TRUE,
  data    = surv_data
)

Fit a Cox model for all-cause mortality adjusting for age, sex, HDL, prior MI, and cholesterol-lowering medication:

fit <- svycoxph(
  Surv(time, event) ~ TC_sd + HDL_sd + RIDAGEYR + RIAGENDR +
                      MI_history + chol_med,
  design = design
)

summary(fit)
round(exp(cbind(HR = coef(fit), confint(fit))), 3)

Interpreting the output:

  • TC_sd: hazard ratio per one-SD higher total cholesterol, adjusted for all other covariates. The direction often reverses after adjusting for HDL and statin use — an important confounding structure in lipid epidemiology.
  • HDL_sd: higher HDL is typically protective (HR < 1).
  • RIAGENDR: coded 1 = male, 2 = female; HR compares females to males.
  • MI_history: coded 1 = prior MI, 0 = none; HR estimates excess mortality risk in those with a history of heart attack.
  • chol_med: coded 1 = currently on cholesterol-lowering medication, 0 = no. Two distinct biases apply simultaneously:
    • Confounding by indication: participants prescribed a statin typically had higher pre-treatment cholesterol and greater cardiovascular risk, so statin users are sicker on average than their measured TC suggests.
    • Exposure mismeasurement: statins lower TC by approximately 30–40 mg/dL, so the measured TC in treated individuals systematically underestimates their underlying lipid burden. Including chol_med as a covariate adjusts for the group difference but does not recover the pre-treatment value. Common analytic responses include restricting the analysis to untreated participants, imputing pre-treatment TC by adding back an estimated treatment effect, or stratifying by medication status and reporting separate associations.

Notes on data management

Cross-cycle variable harmonization

Many NHANES analytes changed variable names or file names across cycles. Use nhanes_search_variables() to discover what exists, nhanes_variable_map() to get the per-cycle file names, and nhanes_download_analyte() + nhanes_harmonize() to download and rename consistently:

# General pattern for any analyte
analyte_list <- nhanes_download_analyte("search term", cycles,
                                         component = "Laboratory")
analyte      <- nhanes_harmonize(analyte_list,
                                  unit  = "mg/dL",
                                  name  = "my_variable",
                                  label_pattern = "search term")

# For questionnaire variables (no unit to match), use mapping instead
quest_list <- nhanes_download_analyte("keyword", cycles,
               component = "Questionnaire",
               keep_vars = c("VAR_OLD", "VAR_NEW"))
quest <- nhanes_harmonize(quest_list,
           mapping = c(VAR_OLD = "my_flag", VAR_NEW = "my_flag"))

SEQN is a character identifier

nhanesR automatically converts SEQN to character on download. Never use it in arithmetic. Always include "cycle" as a second join key when merging pooled multi-cycle data — SEQNs are unique only within a cycle.

NHANES questionnaire coding

Most Yes/No questionnaire items use: 1 = Yes, 2 = No, 7 = Refused, 9 = Don't know. Always recode to 0/1 before analysis and set 7 and 9 to NA.

Notes on the public-use LMF

  • Asymmetric follow-up: all public-use LMF files censor at December 31, 2019 regardless of cycle. Participants from 2017–2018 have at most ~2 years of follow-up; those from 1999–2000 have up to ~20 years.
  • Data perturbation: PERMTH_EXM, PERMTH_INT, and UCOD_LEADING contain synthetic values for select records to reduce re-identification risk. MORTSTAT and ELIGSTAT are not perturbed.
  • Restricted-use files: the 2022-linked files extend follow-up to December 31, 2022. They require an approved project and RDC access.

Further reading: