SPPH 604 · demonstration paper (Kueh et al., BMJ Open 2026;16:e113719)
What this file does. It re-derives the eligibility funnel and the four phenotype groups from the merged NHANES data and shows that we reach N = 6,371 and groups 391 / 4585 / 787 / 608. The full multi-cycle download + merge is in
../../reproduction/R/01_download.Rand02_build_analytic.R; here we read that merged frame and expose the eligibility logic itself, because deciding who is in the sample is the P1 skill.
source("../_shared/paths.R") # find_repro(), emit_results(), ci(), fmtn(), endash()
repro <- find_repro()
suppressPackageStartupMessages({library(dplyr)})
# Merged, all-participant frame produced by 02_build_analytic.R (raw + derived cols)
d <- readRDS(file.path(repro, "data/derived/merged_all.rds"))
yn <- function(x) as.integer(x == 1)
nrow(d) # all NHANES 2007-2018 records[1] 59842FLI (Bedogni 2006) needs fasting triglycerides, GGT, BMI and waist. Steatosis = FLI ≥ 60.
d <- d %>% mutate(
fli = { L <- 0.953*log(LBXTR) + 0.139*BMXBMI + 0.718*log(LBXSGTSI) +
0.053*BMXWAIST - 15.745; exp(L)/(1+exp(L))*100 },
steatosis = as.integer(fli >= 60))d <- d %>% mutate(
crit_adip = as.integer(bmi >= 25 | (female==0 & waist > 94) | (female==1 & waist > 80)),
crit_glu = as.integer(LBXGLU >= 100 | LBXGH >= 5.7 | yn(DIQ010)==1 |
yn(DIQ070)==1 | yn(DIQ050)==1),
crit_bp = as.integer(sbp >= 130 | dbp >= 85 | yn(BPQ050A)==1),
crit_tg = as.integer(LBXTR >= 150 | yn(BPQ100D)==1),
crit_hdl = as.integer((female==0 & LBDHDD <= 40) | (female==1 & LBDHDD < 50)),
n_crit = rowSums(cbind(crit_adip,crit_glu,crit_bp,crit_tg,crit_hdl), na.rm = TRUE),
masld = as.integer(steatosis == 1 & n_crit >= 1))The Table headers say “BMI < 25”, but that is a mislabel — the non-obese groups have BMI medians of 27.6 / 28.8. The Methods-text cut (30) is what reproduces the group sizes.
d <- d %>% mutate(
obese = as.integer(bmi >= 30 | (asian == 1 & bmi >= 25)),
central = as.integer(whtr >= 0.6),
group = factor(case_when(obese==1 & central==0 ~ "I",
obese==1 & central==1 ~ "II",
obese==0 & central==0 ~ "III",
obese==0 & central==1 ~ "IV"),
levels = c("I","II","III","IV")))Every step carries a type label. A funnel shows N falling; it does not show why, and the reasons are not interchangeable. Sorting them is the point of this step.
f <- function(cond) sum(cond, na.rm = TRUE)
inFast <- !is.na(d$WTSAF2YR) & d$WTSAF2YR > 0
funnel <- tibble::tribble(
~Step, ~Criterion, ~N, ~Type,
0, "NHANES 2007-2018 records", nrow(d), "-",
1, "Age >= 18", f(d$age>=18), "ELIGIBILITY",
2, "In fasting subsample", f(d$age>=18 & inFast), "DESIGN",
3, "FLI computable (TG+GGT+BMI+waist)", f(d$age>=18 & inFast & !is.na(d$fli)),
"MISSINGNESS",
4, "Steatosis (FLI >= 60)", f(d$age>=18 & inFast & d$steatosis==1),
"TARGET POPULATION",
5, "MASLD (>= 1 criterion)", f(d$age>=18 & inFast & d$masld==1),
"TARGET POPULATION",
6, "Phenotype group assignable", f(d$age>=18 & inFast & d$masld==1 & !is.na(d$group)),
"MISSINGNESS",
7, "Mortality-eligible and linked", f(d$age>=18 & inFast & d$masld==1 & !is.na(d$group) &
d$eligstat==1 & !is.na(d$mortstat) & !is.na(d$time_yr)),
"OUTCOME ASCERTAINMENT")
funnel$Dropped <- c(NA, -diff(funnel$N))
funnel| Step | Criterion | N | Type | Dropped |
|---|---|---|---|---|
| 0 | NHANES 2007-2018 records | 59842 | - | NA |
| 1 | Age >= 18 | 36580 | ELIGIBILITY | 23262 |
| 2 | In fasting subsample | 14962 | DESIGN | 21618 |
| 3 | FLI computable (TG+GGT+BMI+waist) | 14170 | MISSINGNESS | 792 |
| 4 | Steatosis (FLI >= 60) | 6052 | TARGET POPULATION | 8118 |
| 5 | MASLD (>= 1 criterion) | 6050 | TARGET POPULATION | 2 |
| 6 | Phenotype group assignable | 6050 | MISSINGNESS | 0 |
| 7 | Mortality-eligible and linked | 6048 | OUTCOME ASCERTAINMENT | 2 |
The five categories are not five flavours of the same thing. They have different inferential consequences, and a funnel that hides the difference invites the wrong fix.
WTSAF2YR to account for the extra selection stage and its nonresponse. The correct response is the right weight and a domain analysis, not imputation. Inference stays conditional on the weighting assumptions.The single most common error on a cohort funnel is to call every downward step “attrition” and then reach for imputation. Three of these five categories cannot be imputed, and one of them should not be.
coh <- d %>% filter(age>=18, masld==1, !is.na(group), eligstat==1,
!is.na(mortstat), !is.na(time_yr))
data.frame(
Group = c("I","II","III","IV"),
Ours_N = as.integer(table(coh$group)),
Paper_N= c(386,4541,769,604),
Ours_deaths = as.integer(tapply(coh$dead, coh$group, sum, na.rm=TRUE)),
Paper_deaths= c(5,388,74,118))| Group | Ours_N | Paper_N | Ours_deaths | Paper_deaths |
|---|---|---|---|---|
| I | 391 | 386 | 5 | 5 |
| II | 4585 | 4541 | 389 | 388 |
| III | 787 | 769 | 74 | 74 |
| IV | 608 | 604 | 118 | 118 |
cat("Total N:", nrow(coh), "(paper 6300) Total deaths:",
sum(coh$dead,na.rm=TRUE), "(paper 585)\n")Total N: 6371 (paper 6300) Total deaths: 586 (paper 585)One sentence in the Methods, coded three defensible ways, spans > 2,000 people. We do not tune to hit N; we keep the no-exclusion cohort (which matches the paper) and flag the discrepancy.
base <- coh
## NHANES codes "refused" as 777 and "don't know" as 999. Left raw, a rule like
## `ALQ130 > 2` silently reclassifies every refusal as a heavy drinker: here that
## is 9 people, which moves the excluded count from 2,063 to 2,072. Reserved
## codes are not measurements, and every increment must strip them before use.
base$alq <- ifelse(base$ALQ130 %in% c(777, 999), NA, base$ALQ130)
crude <- with(base, (female==0 & alq>2) | (female==1 & alq>1)) # ALQ130 alone
avg <- with(base, (female==0 & avg_drinks_day>2) | (female==1 & avg_drinks_day>1))
data.frame(
coding = c("Crude ALQ130 rule", "Frequency-weighted average", "No effective exclusion"),
excluded = c(sum(crude,na.rm=TRUE), sum(avg,na.rm=TRUE), 0),
N_remaining = c(nrow(base)-sum(crude,na.rm=TRUE), nrow(base)-sum(avg,na.rm=TRUE), nrow(base)))| coding | excluded | N_remaining |
|---|---|---|
| Crude ALQ130 rule | 2063 | 4308 |
| Frequency-weighted average | 427 | 5944 |
| No effective exclusion | 0 | 6371 |
cc <- coh %>% filter(!is.na(cancer), !is.na(stroke), !is.na(hf), !is.na(LBXSCR))
cat("Full cohort:", nrow(coh), " | complete-case on covariates:", nrow(cc),
" | missing serum creatinine:", sum(is.na(coh$LBXSCR)), "\n")Full cohort: 6371 | complete-case on covariates: 6246 | missing serum creatinine: 1 These counts do not account for the residual gap. Our N is above the paper’s ~6,300, and dropping rows with incomplete covariates would take it below. We report the +71 as unresolved rather than attribute it.
M1 cites these by name; it does not re-run this file. Every number is formatted once, here, so the deck cannot reformat or mistype it.
fn <- funnel$N; names(fn) <- paste0("s", funnel$Step)
emit_results("P1", list(
full_n = fmtn(nrow(coh)),
groups = paste(fmtn(as.integer(table(coh$group))), collapse = " / "),
deaths = fmtn(sum(coh$dead, na.rm = TRUE)),
locked_n = fmtn(sum(!is.na(coh$WTSAF2YR) & coh$WTSAF2YR > 0)),
funnel_records = fmtn(fn[["s0"]]),
funnel_adults = fmtn(fn[["s1"]]),
funnel_fasting = fmtn(fn[["s2"]]),
funnel_fli = fmtn(fn[["s3"]]),
funnel_steatosis = fmtn(fn[["s4"]]),
funnel_masld = fmtn(fn[["s5"]]),
funnel_linked = fmtn(fn[["s7"]]),
item_nonresponse = fmtn(funnel$Dropped[funnel$Step == 3]),
outcome_loss = fmtn(funnel$Dropped[funnel$Step == 7]),
# what "read every downward step as attrition" would actually total
all_downward_steps = fmtn(fn[["s1"]] - fn[["s7"]]),
alq_excluded = fmtn(sum(crude, na.rm = TRUE)),
alq_n = fmtn(nrow(base) - sum(crude, na.rm = TRUE)),
avg_excluded = fmtn(sum(avg, na.rm = TRUE)),
avg_n = fmtn(nrow(base) - sum(avg, na.rm = TRUE)),
alq_reserved = fmtn(sum(base$ALQ130 %in% c(777, 999)))
), list(funnel = as.data.frame(funnel)))wrote ../_results/P1.json (19 facts, 1 tables)