P4 — Design-aware main estimate: step-by-step reproducible code

SPPH 604 · demonstration paper (Kueh et al., BMJ Open 2026;16:e113719)

Replicates the paper’s unweighted Cox model, then re-estimates it with the NHANES fasting-subsample survey design (weights + strata + clusters).

The rule this file exists to demonstrate: build the survey design on the whole sampled frame, then reach your analysis population with subset(design, ...). Never filter the data frame and then call svydesign() on the remains.

source("../_shared/paths.R")   # find_repro(), emit_results(), ci(), fmtn(), endash()
repro <- find_repro()
suppressPackageStartupMessages({library(dplyr); library(survival); library(survey)})
options(survey.lonely.psu = "adjust")

frame <- readRDS(file.path(repro, "data/derived/merged_all.rds"))   # 59,842 records
d     <- readRDS(file.path(repro, "data/derived/masld_analytic.rds"))
frame$in_analysis <- frame$SEQN %in% d$SEQN
frame$group <- factor(frame$group, levels = c("I","II","III","IV"))

hr <- function(m,g){s<-summary(m); sprintf("%.2f (%.2f-%.2f)",
        s$conf.int[g,1], s$conf.int[g,"lower .95"], s$conf.int[g,"upper .95"])}
psu <- function(df) length(unique(paste(df$SDMVSTRA, df$SDMVPSU)))

Step 1 — Replicate the paper (unweighted Cox)

Two row sets, labelled, because they are not interchangeable. The full analytic file (n = 6,371) is what reproduces the paper. The locked domain (n = 6,048, those with a valid fasting weight) is what P1 locks and what P2, P3 and P5 analyse — and the only set on which a design-aware estimate is possible.

dom0 <- d[!is.na(d$WTSAF2YR) & d$WTSAF2YR > 0, ]
mu_full <- coxph(Surv(time_yr,dead) ~ group, data=d)
mu_dom  <- coxph(Surv(time_yr,dead) ~ group, data=dom0)
cat("Unweighted, full analytic file  n =", nrow(d),
    ": IV", hr(mu_full,"groupIV"), " <- reproduces the paper
")

Unweighted, full analytic file  n = 6371 : IV 15.14 (6.19-37.04)  <- reproduces the paper

cat("Unweighted, locked domain       n =", nrow(dom0),
    ": IV", hr(mu_dom,"groupIV"),  " <- the rows P2/P3/P5 use
")

Unweighted, locked domain       n = 6048 : IV 16.96 (6.24-46.05)  <- the rows P2/P3/P5 use

The difference between them is the records without a valid fasting weight, not a modelling choice. Everything below uses the locked domain.

Step 2 — Improve: the survey design, built before restricting

FLI needs fasting triglycerides, so the binding constraint is the fasting subsample and the weight is the pooled fasting weight WTSAF2YR / 6 (6 cycles).

The design is constructed on everyone with a valid fasting weight — not on the MASLD cohort — and the analysis population is then reached with subset(). This is a domain (subpopulation) analysis: subset.survey.design retains the original strata and PSU structure so the variance is estimated against the design that actually produced the sample.

fast <- frame %>% filter(!is.na(WTSAF2YR) & WTSAF2YR > 0) %>% mutate(wt = WTSAF2YR/6)
des_full <- svydesign(ids=~SDMVPSU, strata=~SDMVSTRA, weights=~wt, data=fast, nest=TRUE)

des <- subset(des_full, in_analysis)      # <- the domain. NOT a new svydesign().

dom <- fast %>% filter(in_analysis)
cat("design frame: n =", nrow(fast), " PSUs =", psu(fast),
    " strata =", length(unique(fast$SDMVSTRA)), "\n")

design frame: n = 17208  PSUs = 184  strata = 90 

cat("analysis domain: n =", nrow(dom),
    " weighted population =", format(round(sum(dom$wt)), big.mark=","), "\n")

analysis domain: n = 6048  weighted population = 94,334,495 

mu <- coxph(Surv(time_yr,dead) ~ group, data=dom)          # unweighted, same rows
mw <- svycoxph(Surv(time_yr,dead) ~ group, design=des)
cat("Unweighted (domain): IV", hr(mu,"groupIV"), "\n")

Unweighted (domain): IV 16.96 (6.24-46.05) 

cat("Weighted:  II", hr(mw,"groupII"), " III", hr(mw,"groupIII"),
    " IV", hr(mw,"groupIV"), "\n")

Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(des_full, in_analysis)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(des_full, in_analysis)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(des_full, in_analysis)
Weighted:  II 5.38 (1.19-24.39)  III 4.04 (0.89-18.34)  IV 12.40 (2.69-57.15) 

Step 3 — Direct contrasts the paper did not report (weighted)

Change the reference level inside the existing design with update(). Do not rebuild svydesign() — that would discard the frame a second time.

des3 <- update(des, group = relevel(group, "III"))
des2 <- update(des, group = relevel(group, "II"))
mw3 <- svycoxph(Surv(time_yr,dead) ~ group, design=des3)
mw2 <- svycoxph(Surv(time_yr,dead) ~ group, design=des2)
cat("Weighted IV vs III:", hr(mw3,"groupIV"), "   IV vs II:", hr(mw2,"groupIV"), "\n")

Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
update(des, group = relevel(group, "III"))
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
update(des, group = relevel(group, "II"))
Weighted IV vs III: 3.07 (2.11-4.46)    IV vs II: 2.30 (1.76-3.02) 

Step 3b — Separating adjustment from design on one row set

The ledger 15.1 -> 2.9 -> 12.4 mixes crude and adjusted models on different rows, so it is not a sequence of successive corrections. Hold the rows fixed and the two effects separate cleanly.

COV <- c("age","female","dyslip","htn","t2dm","ckd","mi","cancer",
         "on_htn_med","on_dm_med","on_lipid_med")
fast$cc <- fast$in_analysis & complete.cases(fast[,COV])
dcc <- fast[fast$cc, ]
f0 <- Surv(time_yr,dead) ~ group
f1 <- as.formula(paste("Surv(time_yr,dead) ~ group +", paste(COV, collapse=" + ")))
gg <- function(f){s<-summary(f); sprintf('%.2f (%.2f-%.2f)',
        s$conf.int['groupIV',1], s$conf.int['groupIV',3], s$conf.int['groupIV',4])}
rl <- function(x,r){x$group <- relevel(x$group, r); x}
dR <- function(r) subset(svydesign(ids=~SDMVPSU, strata=~SDMVSTRA, weights=~wt,
                    data=rl(fast,r), nest=TRUE), cc)
mat <- data.frame(
  cell = c('unweighted crude','unweighted Model 1','weighted crude','weighted Model 1'),
  IV_vs_I = c(gg(coxph(f0,data=dcc)), gg(coxph(f1,data=dcc)),
              gg(svycoxph(f0,design=dR('I'))), gg(svycoxph(f1,design=dR('I')))),
  IV_vs_III = c(gg(coxph(f0,data=rl(dcc,'III'))), gg(coxph(f1,data=rl(dcc,'III'))),
                gg(svycoxph(f0,design=dR('III'))), gg(svycoxph(f1,design=dR('III')))))

Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt, 
    data = rl(fast, r), nest = TRUE), cc)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt, 
    data = rl(fast, r), nest = TRUE), cc)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt, 
    data = rl(fast, r), nest = TRUE), cc)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt, 
    data = rl(fast, r), nest = TRUE), cc)

mat

cell	IV_vs_I	IV_vs_III
unweighted crude	16.05 (5.91-43.58)	2.34 (1.72-3.18)
unweighted Model 1	3.81 (1.39-10.43)	1.05 (0.77-1.45)
weighted crude	11.93 (2.60-54.84)	3.07 (2.11-4.48)
weighted Model 1	2.98 (0.69-12.82)	1.29 (0.89-1.86)

Read down a column for adjustment, across a row for design. On the locked IV-vs-III estimand adjustment dominates and design changes no conclusion; on the replication IV-vs-I estimand design changes whether the adjusted interval excludes 1.

Step 4 — Why the rule matters (and exactly when it bites)

For this main estimate both approaches give the same answer, because the MASLD domain happens to contain every PSU and stratum in the fasting frame. Pre-filtering loses nothing here. That is luck, not method.

show <- function(lbl, df) cat(sprintf("%-32s n=%5d  PSUs=%4d  strata=%3d
", lbl, nrow(df), psu(df), length(unique(df$SDMVSTRA))))
show("full fasting frame", fast)

full fasting frame               n=17208  PSUs= 184  strata= 90

show("MASLD analysis domain", dom)

MASLD analysis domain            n= 6048  PSUs= 184  strata= 90

show("  Group I only", dom %>% filter(group=="I"))

  Group I only                   n=  375  PSUs= 155  strata= 89

show("  Group I, women only", dom %>% filter(group=="I", female==1))

  Group I, women only            n=   58  PSUs=  51  strata= 44

The domain keeps all 184 PSUs. A subgroup does not: Group I women sit in ~51 PSUs spread over 44 strata, which leaves strata containing a single cluster. That is where pre-filtering stops being harmless.

d2 <- dom %>% filter(group=="I", female==1)
for (opt in c("fail","adjust")) {
  options(survey.lonely.psu = opt)
  dfull <- svydesign(ids=~SDMVPSU, strata=~SDMVSTRA, weights=~wt, data=fast, nest=TRUE)
  pre <- tryCatch(sprintf("%.4f", SE(svymean(~age, svydesign(ids=~SDMVPSU,
           strata=~SDMVSTRA, weights=~wt, data=d2, nest=TRUE)))),
         error=function(e) paste0("ERROR: ", conditionMessage(e)))
  sub <- sprintf("%.4f", SE(svymean(~age,
           subset(dfull, in_analysis & group=="I" & female==1))))
  cat(sprintf("lonely.psu=%-7s  pre-filtered: %-45s  subset(): %s
", opt, pre, sub))
}

lonely.psu=fail     pre-filtered: ERROR: Stratum (59) has only one PSU at stage 1  subset(): 2.3852
lonely.psu=adjust   pre-filtered: 2.3852                                         subset(): 2.3852

options(survey.lonely.psu = "adjust")

Read that table carefully, because the lesson is subtler than “pre-filtering gives the wrong number”:

• Under R’s default (lonely.psu = "fail"), pre-filtering does not give a wrong answer — it gives no answer. The design you handed svydesign() has strata with one cluster and no variance is estimable.
• Setting lonely.psu = "adjust" makes it run, and here it returns the same SE. But you have silenced a symptom of having thrown the design away, and you are now relying on an approximation instead of the stratum structure you actually had.
• subset() needs no rescue. It works under either setting, because the strata and PSUs were never discarded.

So the rule is not “pre-filtering biases your estimate.” It is: pre-filtering can create lonely-PSU strata that would not arise under proper domain analysis, making the variance estimate depend on the chosen lonely-PSU fallback. Restricting the 6,048-person MASLD cohort did nothing at all, because it retained every PSU and stratum. The failure appears once restriction starts dropping them. Sometimes the fallback agrees; you will not know which time you are in.

Note what this separates: a correct procedure is not the same as a procedure that changes today’s answer. The wrong implementation returns an identical HR and CI in the main analysis and is still wrong, because it is not guaranteed to preserve the variance structure for narrower domains.

Your own paper’s subgroups will be sparser than this one. Build the design once, then subset.

Interpretation. On these 6,048 rows (crude, IV vs I), weighting moves the estimate from about 17 to about 12 and widens the interval. Weighting has no fixed direction: in the 2×2 on the matched 5,939 rows, the crude IV-vs-III estimate rises when weighted. The wider interval is the price of estimating for the US population rather than for this sample; the unweighted interval is not wrong, it answers the unweighted question. Crude and weighted, Group IV still has the highest hazard against the other high-risk groups (IV vs III ≈ 3.1; IV vs II ≈ 2.3).

On which weight. CDC’s rule is the least common denominator: when an analysis draws on variables measured in different components, use the weight belonging to the smallest applicable subsample. FLI requires fasting triglycerides, so the fasting weight governs, even though GGT is measured on the larger MEC sample. This is a rule, not a judgment call. Both weights target the same US population, but only the fasting weight accounts for selection into the fasting subsample and its nonresponse, which is where FLI comes from.

Results for M1

cell <- function(row, col) endash(mat[[col]][mat$cell == row])
emit_results("P4", list(
  uw_crude_i    = cell("unweighted crude",   "IV_vs_I"),
  uw_crude_iii  = cell("unweighted crude",   "IV_vs_III"),
  uw_m1_i       = cell("unweighted Model 1", "IV_vs_I"),
  uw_m1_iii     = cell("unweighted Model 1", "IV_vs_III"),
  w_crude_i     = cell("weighted crude",     "IV_vs_I"),
  w_crude_iii   = cell("weighted crude",     "IV_vs_III"),
  w_m1_i        = cell("weighted Model 1",   "IV_vs_I"),
  w_m1_iii      = cell("weighted Model 1",   "IV_vs_III"),
  n_matched     = fmtn(nrow(dcc)),
  psu_full      = fmtn(psu(fast)),
  psu_domain    = fmtn(psu(dom))
), list(matrix = mat))

wrote ../_results/P4.json  (11 facts, 1 tables)