The paper states two aims, not one

• Aim 1 — the clinical correlates of four groups defined by BMI and waist-to-height ratio
• Aim 2 — the prognostic value of that classification for mortality

A single paper can carry more than one research question, and they have different clocks.

Five goals the same variables can serve

The variables can be identical. The target quantity is not.

Same variables, different analyses

A Cox model appears in several rows. The model family does not determine the goal.

Same variables, different evidence of success

HR + CI: what it is, and is not, evidence for

HR + CI is enough to report an adjusted association, or a prognostic-factor effect.

It is not evidence that a prediction model is good, and not by itself evidence that a causal effect has been identified.

The paper names “prognostic value” and supplies adjusted association.

What does this paper most directly support?

It directly supports: description of how the groups differ; association between phenotype and subsequent mortality; an adjusted prognostic association.

It does not, by itself, establish: that central adiposity causes the difference; that reducing WHtR would reduce mortality; that adding WHtR improves individual prediction.

First: how NHANES is organised

• Six two-year cycles, 2007–2018, each with its own file suffix: E F G H I J
• Files are .XPT (SAS transport) — haven::read_xpt(), or nhanesA::nhanes()
• Each file is one row per participant, for one topic — not one big table
• Three collection settings: Interview (home) · MEC exam (measured body size, BP) · Laboratory (triglycerides, GGT, glucose)

Mortality is not an NHANES file. It is a separate linked, fixed-width file from NCHS.

Where the measurements and labs come from

Anthropometry, blood pressure and laboratory values are measured or assayed. Age, sex and race/ethnicity come from the demographic file.

Where the questionnaire, design and outcome come from

These variables alone reach across seven examination and laboratory components, plus demographics, plus questionnaire components, plus a separate linked mortality file. That is what L1 assembles.

SEQN is the join key

• Every respondent has one SEQN — the identifier CDC uses to sort and merge files
• Analysis = pick the topic files you need, merge on SEQN, one row per person
• SEQN does not repeat across the continuous cycles: all 59,842 values across 2007–2018 are distinct, and the per-cycle ranges do not overlap — Analysis set: all NHANES 2007–2018 records, N = 59,842
• The ranges are not contiguous, though — a gap of 1,640 sits between 2011–2012 and 2013–2014, so do not infer a cycle from an ID
• Keep a cycle variable anyway: provenance, definitions that change between cycles, and weights built per cycle. Not to disambiguate SEQN

L1 in one screen: two files, one key

pacman::p_load(tidyverse, nhanesA, tableone, naniar)

demo <- nhanes("DEMO_H", translated = TRUE)   # 2013-2014 demographics
bmx  <- nhanes("BMX_H",  translated = TRUE)   # 2013-2014 body measures

merged <- full_join(demo, bmx, by = "SEQN")

• DEMO_H and BMX_H are the same DEMO and BMX from the “Where the measurements and labs come from” table, cycle H
• translated = TRUE turns stored codes into labels — convenient, and it changes what a value looks like, so recode rules written against raw codes will silently miss
• full_join keeps people measured in one file but not the other. That is a row-set decision made by a single word: inner_join would be a different study

Then do it for every cycle

cycles <- c("H", "I", "J")                     # 2013-14, 2015-16, 2017-18

all_cycles <- lapply(cycles, function(cy) {
  d <- nhanes(paste0("DEMO_", cy), translated = TRUE)
  b <- nhanes(paste0("BMX_",  cy), translated = TRUE)
  full_join(d, b, by = "SEQN") %>% mutate(cycle = cy)
})

stacked <- bind_rows(all_cycles)

• paste0("DEMO_", cy) is why the suffix table was worth memorising
• bind_rows() stacks cycles; the joins happen within a cycle first
• Our reproduction does the same over all six cycles, E through J

Raw names in, analysis names out

analysis <- stacked %>%
  rename(Sex = RIAGENDR, Age = RIDAGEYR,
         RaceEthnicity = RIDRETH3, BMI = BMXBMI) %>%
  filter(Age >= 20) %>%
  mutate(BMICat = cut(BMI, c(0, 18.5, 25, 30, Inf), right = FALSE,
                      labels = c("Underweight", "Normal weight",
                                 "Overweight", "Obese")))

• Every rename is a decision you have to be able to defend, not cosmetics
• right = FALSE is not optional. Without it cut() closes the intervals the other way and BMI exactly 30.0 lands in Overweight — one argument, and the categories are the mirror image of the ones you meant
• BMICat here uses the standard WHO cut-points. The paper does not: obesity at BMI ≥ 30, and ≥ 25 for Asian participants
• The lab keeps age ≥ 20; the paper says age > 18; our reproduction uses ≥ 18

Same variable, three eligibility rules, three different studies.

Look at what is missing before you model it

colSums(is.na(analysis))
naniar::gg_miss_var(analysis, show_pct = TRUE)

tab1_vars <- c("Age", "RaceEthnicity", "BMI", "BMICat")
CreateTableOne(vars = tab1_vars, data = analysis, strata = "Sex",
               addOverall = TRUE, test = FALSE)

• Missingness is inspected first, then decided about — never discovered inside a model
• CreateTableOne is the Table 1 you will rebuild for the paper’s cohort in M1
• test = FALSE on purpose: a p-value comparing baseline groups answers no question anyone asked

Population: who is being studied?

Adults with MASLD in NHANES, 2007–2018

MASLD = metabolic dysfunction-associated steatotic liver disease: excess liver fat together with metabolic dysfunction.

But nobody looked at a liver. Steatosis was estimated by a formula.

The study population is not the general population. Everyone in the main analysis already meets the study’s own definition of MASLD.

How MASLD was actually ascertained

Hepatic steatosis was estimated with the Fatty Liver Index (Bedogni 2006), thresholded at FLI ≥ 60:

The four enter a published logistic formula — not a product of the four. What matters here is that the index is computed, from four variables in three component files.

Plus at least one of five cardiometabolic criteria: adiposity; raised glucose/HbA1c or diabetes; raised blood pressure or treatment; raised triglycerides or treatment; low HDL.

So “MASLD” here is a derived variable, not a diagnosis — no imaging, no biopsy.

The funnel, with a type on every row

After eligibility, the two largest drops are design and disease definition — and neither is missing data.

Analysis set: locked mortality-linked cohort, N = 6,048.

What the funnel reveals

Reproducibility note. The published paper reports ~6,300 participants with MASLD; this pipeline’s full reconstructed analytic file has 6,371. The funnel’s 6,048 is a further subset, those with a valid fasting weight. Analysis set: full reconstructed analytic file, N = 6,371; locked cohort, N = 6,048.

The paper’s ~6,300 and our 6,371 should not be presented as if they described the same cohort.

The paper also states eligibility as age > 18, while this reproduction uses age ≥ 18. Small implementation differences like this are recorded, not silently harmonised.

The exposure is constructed, not measured

The exposure is not one measurement. It is built in two steps:

BMXWAIST ÷ BMXHT → waist-to-height ratio · BMXBMI → obesity category

→ BMI category × WHtR category → four phenotype groups

• All three inputs come from one component, BMX, in the MEC exam — so a person missing the exam is missing the exposure entirely
• Contrast that with eligibility, which needs BMX and the laboratory (LBXTR, LBXSGTSI) — eligibility spans two settings, the exposure does not

Exposure component 1: BMI

\[ \text{BMI} = \text{weight (kg)} / \text{height}^2 \text{ (m}^2\text{)} \]

• Weight and height were measured in the exam, not self-reported (BMXBMI)
• BMI marks overall body size — not body fat, and not where fat is stored
• Obesity defined as BMI ≥ 30; for Asian participants BMI ≥ 25 (RIDRETH3 == 6)

Paper reporting note. The Methods define obesity as BMI ≥ 30 (≥ 25 for Asian participants), but published Tables 1–2 label the groups as BMI ≥ 25 versus < 25. Those labels conflict with the Methods and with the reported group BMI distributions, so they appear to be a reporting error.

Exposure component 2: waist-to-height ratio

\[ \text{WHtR} = \text{waist (cm)} / \text{height (cm)} = \texttt{BMXWAIST} / \texttt{BMXHT} \]

A person 170 cm tall with a 102 cm waist has WHtR = 102/170 = 0.60.

• High central adiposity is defined as WHtR ≥ 0.60
• WHtR is a simple marker of abdominal adiposity — not CT or MRI visceral fat

The exposure is the combination

Group IV — not obese by BMI, but a large waist for their height. The study asks whether BMI alone misses this higher-risk phenotype.

• There is no treated arm and control arm: the four-category phenotype is the index factor, and mortality is compared across the categories
• Every hazard ratio is measured against Group I. A ratio without its referent means nothing — and here the referent is an obese group

Eligibility and exposure share their inputs

Being selected into the study and being placed in a phenotype group are not independent events. That is structural dependence.

For a causal question, conditioning on this selected population could open selection or collider paths under some structures; for a prognostic question, it primarily defines a selected target population. It is not automatically circularity, and not automatically bias.

Never trust a label you have not recomputed.

Outcome and the clock

• All-cause mortality (MORTSTAT) via the NCHS linked mortality file — not an NHANES file, a separate linked fixed-width one
• Time zero is the MEC examination; PERMTH_EXM is person-months from that exam
• Follow-up through 31 December 2019, stated in the NCHS 2019 Linked Mortality File linkage methodology document — not inferred from the filename
• Participants entered between 2007 and 2018, so follow-up lengths and censoring differ — time-to-event methods are appropriate. A risk difference at a fixed horizon can still be estimated, using methods that account for censoring

Time zero is a question the checklist never asks

• “Exposure, outcome, covariates” can all be specified while time zero stays ambiguous
• An ambiguous time zero silently changes who is at risk, and for how long
• Here exposure is measured at time zero — convenient, and rarer than you would think

Ask of your M0 candidates: when does follow-up begin, and does everyone start together?

“In the model” is not the same as “confounding”

• The paper adjusts for age, sex, comorbidities, medications
• “We adjusted for X” says where X sits in the syntax. It does not say what X does.
• The same variable can be a confounder, a mediator, a competing prognostic factor, or a collider — and which it is depends on the research goal from Q1

What a crude-to-adjusted change establishes

• Group IV vs I all-cause HR ≈ 15.1 crude → ≈ 2.9 with the paper’s Model 2
• Conditioning on the adjustment set moves the estimate enormously
• What that movement means — confounding controlled, mediation removed, a collider opened — depends on the goal and the causal role of those variables

Analysis set: full reconstructed analytic file, N = 6,371 (Model 2 on its complete cases); unweighted, as in the paper.

Five NHANES traps — each has bitten this project

• ALQ130 is drinks per drinking day, not drinks per day. A “>2 a day” rule applied to it directly removes 2,063 people here, against 427 for average intake (frequency × quantity): about 1,600 over-excluded. Week 3 owns this one
• 777 = refused, 999 = don’t know — codes, not counts. Left inside a threshold rule they reclassified nine such records as heavy drinkers. Two traps on one variable, a factor of about two hundred apart
• RIDRETH3 does not exist before 2011. It is what identifies Asian participants for the lower BMI threshold — so the first two cycles cannot apply the paper’s own rule
• Fasting labs carry their own weight. LBXTR and LBXGLU come from the fasting subsample, weighted by WTSAF2YR, not WTMEC2YR
• Blood pressure is not a row mean — next slide

Analysis set: full reconstructed analytic file, N = 6,371 (the alcohol counts).

The blood pressure trap, in full

The NHANES Blood Pressure Procedures Manual (s3.1.7.3) specifies the averaging:

• The first reading is excluded whenever more than one was obtained
• An all-zero diastolic set averages to zero — a diastolic 0 is a measurement, and NHANES codes a failure to obtain a reading separately, as missing
• A zero sitting beside a positive reading is omitted from the average

Applying it changed the derived diastolic pressure for 5,523 of 6,371 participants and moved hypertension from 3,444 to 3,448 — while the Model 1 and Model 2 hazard ratios shifted only in the third significant figure. Analysis set: full reconstructed analytic file, N = 6,371.

A coding error can be real, and worth fixing, even when the headline estimate survives it.

NHANES is not a simple random sample

• A complex multistage probability design: weights, strata (SDMVSTRA), PSUs (SDMVPSU)
• We use them in Week 6; today, note they exist and that ignoring them is a choice
• Fasting-subsample analyses carry their own weight — WTSAF2YR, not WTMEC2YR
• Build the design on the full frame, then subset() — never on the subsample alone

Designed subsampling is not missingness

• Sampling into the fasting subsample is statistically designed, with a weight built to keep it nationally representative; compliance and lab nonresponse also play a part
• Check rather than assume: adults in and out of it differ by 0.2 years of age, 0.02 kg/m² BMI, 0.07 cm waist
• It is balanced on the very measurements this paper’s exposure is built from

Analysis set: all NHANES adults 2007–2018, N = 36,580.

Designed subsampling ≠ ordinary missingness.

Survey design: paper versus reproduction

• The published paper’s Statistical analysis section does not report using NHANES survey weights, strata or PSUs
• Our pipeline does include design variables and the fasting-subsample weight
• Therefore a survey-weighted rebuild is a design-informed reanalysis, not a literal reproduction of every published analytic choice

Without the design information, estimates and standard errors should not be described as design-based nationally representative estimates.

PECOT, and the question rebuilt

Among adults with MASLD, are BMI/WHtR phenotypes associated with different subsequent risks of all-cause mortality?

M0 — ask the six of your own papers

• Shortlist two or three candidate NHANES papers; answer all six, each with a pointer to where in the paper you found it
• Reconstructable is necessary, not sufficient — a perfectly specified paper can still be a poor M0 choice
• M0 wants: enough observations and events, a design feature its analysis must respect, and a missing variable the analysis actually uses
• An unanswered question is not a disqualification — it is a candidate for the judgement you defend at M1

Bridge to L1 and M0

• Thursday (L1): the code you saw before the break — several cycles downloaded, topic files read, merged on SEQN, renamed, categorised, missingness inspected, Table 1 built
• Before you go: you should be able to name the component file each variable on these slides comes from
• Toward M0 (Wk3): bring six answers per candidate paper, not a summary

Reading a Methods section is the first act of reproduction.