source("../_shared/paths.R") # find_repro(), emit_results(), ci(), fmtn(), endash()
repro <- find_repro()
suppressPackageStartupMessages({library(dplyr); library(survey); library(survival)
library(mice); library(mitools)})
options(survey.lonely.psu = "adjust")
m <- readRDS(file.path(repro,"data/derived/merged_all.rds"))
a <- readRDS(file.path(repro,"data/derived/masld_analytic.rds"))
m$in_analysis <- m$SEQN %in% a$SEQN
m$group <- factor(m$group, levels=c("I","II","III","IV"))
fast <- m %>% filter(!is.na(WTSAF2YR), WTSAF2YR > 0) %>% mutate(wt = WTSAF2YR/6)
dom <- fast %>% filter(in_analysis)
## The paper's Model 1. Nothing else.
COV <- c("age","female","dyslip","htn","t2dm","ckd","mi","cancer",
"on_htn_med","on_dm_med","on_lipid_med")
f1 <- as.formula(paste("Surv(time_yr,dead) ~ group +", paste(COV, collapse=" + ")))
gv <- function(f,t="groupIV"){s<-summary(f); c(s$conf.int[t,1],s$conf.int[t,3],s$conf.int[t,4])}
fm <- function(x) sprintf("%.2f (%.2f-%.2f)", x[1],x[2],x[3])P5 — Missing data: step-by-step reproducible code
SPPH 604 · demonstration paper (Kueh et al., BMJ Open 2026;16:e113719)
Starts from exactly the target population and survey design used for the corrected main estimate in P4, and uses exactly the paper’s Model 1 covariates. No variable is added to the model in order to create a missing-data problem.
Step 0 — Three different kinds of downward step
A cohort funnel shows N falling. It does not show why, and the reasons have completely different inferential meanings. Sort every step into one of three bins before you call anything “missing data”.
w <- function(lbl, df, kind) cat(sprintf("%-38s %6d %13s %s\n", lbl, nrow(df),
format(round(sum(df$wt)), big.mark=","), kind))
ad <- fast %>% filter(!is.na(age), age>=18)
fli <- ad %>% filter(!is.na(fli))
st <- fli %>% filter(fli>=60)
cat(sprintf("%-38s %6s %13s %s\n","step","n","weighted pop","kind"))step n weighted pop kindw("fasting subsample frame", fast, "DESIGN - the weight handles this")fasting subsample frame 17208 260,290,473 DESIGN - the weight handles thisw(" adults 18+", ad, "ELIGIBILITY- definitional") adults 18+ 14962 234,756,163 ELIGIBILITY- definitionalw(" FLI computable", fli, "MISSINGNESS- item nonresponse") FLI computable 14170 224,125,729 MISSINGNESS- item nonresponsew(" steatosis, FLI >= 60", st, "DISEASE - target population") steatosis, FLI >= 60 6052 94,418,570 DISEASE - target populationw(" MASLD, linked, group known", dom, "DISEASE + small missingness") MASLD, linked, group known 6048 94,334,495 DISEASE + small missingness- Design. Fasting measures come from a designed NHANES subsample. Its absence elsewhere is not missing data; it is why
WTSAF2YRexists. - Disease. Most of the drop is people who do not have steatosis. They are not missing — they are not in the target population.
- Missingness. Only the item-nonresponse step is a missing-data problem, and it is 5.3% of adults in the fasting frame.
The single most common error on a funnel is to call every step “attrition”.
Step 1 — The missingness pattern, before choosing any method
miss <- sapply(COV, function(v) sum(is.na(dom[[v]])))
for (v in names(miss)[miss>0])
cat(sprintf(" %-14s missing %4d %.2f%% unweighted %.2f%% weighted\n", v, miss[v],
100*miss[v]/nrow(dom), 100*sum(dom$wt[is.na(dom[[v]])])/sum(dom$wt))) cancer missing 109 1.80% unweighted 1.18% weightedcc <- complete.cases(dom[,COV])
cat(sprintf("\n variables with any missingness: %d of %d\n", sum(miss>0), length(COV)))
variables with any missingness: 1 of 11cat(sprintf(" complete cases: %d of %d (loses %d; %.2f%% unweighted, %.2f%% weighted)\n",
sum(cc), nrow(dom), sum(!cc), 100*mean(!cc), 100*sum(dom$wt[!cc])/sum(dom$wt))) complete cases: 5939 of 6048 (loses 109; 1.80% unweighted, 1.18% weighted)One Model-1 variable carries any missingness at all. There is no overlap pattern to examine, because there is nothing to overlap with.
Step 2 — State the estimand before comparing anything
The weighted, Model-1-adjusted hazard ratio for Group IV versus Group I, in the MASLD domain of the NHANES fasting subsample, 2007–2018.
Fixing this first matters: complete-case analysis and multiple imputation can silently target different populations, and then a difference between them is not evidence about missing data.
Step 3 — Complete case, on the corrected design
des <- subset(svydesign(ids=~SDMVPSU, strata=~SDMVSTRA, weights=~wt,
data=fast, nest=TRUE), in_analysis)
fit_w <- svycoxph(f1, design=des)
fit_u <- coxph(f1, data=dom)
cat(sprintf("weighted Model-1, complete case: IV = %s (n = %d)\n", fm(gv(fit_w)), fit_w$n))Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = fast, nest = TRUE), in_analysis)
weighted Model-1, complete case: IV = 2.98 (0.69-12.82) (n = 5939)cat(sprintf("unweighted Model-1, complete case: IV = %s (n = %d)\n", fm(gv(fit_u)), fit_u$n))unweighted Model-1, complete case: IV = 3.81 (1.39-10.43) (n = 5939)Step 3b — A stress test that needs no model
The only incomplete Model-1 variable is binary. That permits a sensitivity analysis a reader can check by eye: assign every one of the 109 missing values to one extreme, then the other, and refit. This does not identify the missing values and it is not a formal bounds argument — it is a transparent stress test spanning the two extreme assignments the binary indicator permits.
fitv <- function(val, lbl){
fk <- fast
if (!is.na(val)) fk$cancer[fk$in_analysis & is.na(fk$cancer)] <- val
f <- svycoxph(f1, design = subset(svydesign(ids=~SDMVPSU, strata=~SDMVSTRA,
weights=~wt, data=fk, nest=TRUE), in_analysis))
g <- gv(f)
data.frame(scenario = lbl, HR = round(g[1],2),
lo = round(g[2],2), hi = round(g[3],2), n = f$n)
}
stress <- rbind(fitv(NA, "complete case (109 dropped)"),
fitv(0, "all 109 assigned NO cancer"),
fitv(1, "all 109 assigned CANCER"))Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = fk, nest = TRUE), in_analysis)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = fk, nest = TRUE), in_analysis)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = fk, nest = TRUE), in_analysis)stress| scenario | HR | lo | hi | n |
|---|---|---|---|---|
| complete case (109 dropped) | 2.98 | 0.69 | 12.82 | 5939 |
| all 109 assigned NO cancer | 2.99 | 0.69 | 12.89 | 6048 |
| all 109 assigned CANCER | 2.99 | 0.69 | 12.90 | 6048 |
The Group IV estimate spans 2.98 to 2.99 across both extremes — a total range of 0.01 on the hazard ratio. These are two global extreme scenarios, not bounds over all assignments. A pattern of missingness correlated with phenotype and covariates could move the estimate further than either uniform assignment does, so we do not claim the true value must lie inside this range — only that two transparent extremes leave it essentially unchanged.
Step 4 — Multiple imputation, survey-aware
Assumptions that would make this informative: the missingness is MAR given the Model-1 covariates and the outcome, and the imputation model is correctly specified. The outcome enters as the event indicator plus a Nelson–Aalen cumulative-hazard term — omitting it biases covariate–outcome associations toward the null.
This is not fully survey-modelled multiple imputation, and we do not claim it is. The imputation model is a standard mice() fit with the sampling weight and stratum entered as predictors; the design is reintroduced at the analysis stage, when each completed dataset is refitted with svycoxph() on a properly subsetted design. A fully design-aware imputation would model the sampling structure inside the imputation itself. We report this as a pragmatic sensitivity analysis and check below whether the shortcut matters here.
idat <- dom[, c(COV,"group","dead","time_yr","wt","SDMVSTRA")]
idat$na_cumhaz <- mice::nelsonaalen(idat, time_yr, dead)
ini <- mice(idat, maxit=0, printFlag=FALSE); pr <- ini$predictorMatrix
pr[,"time_yr"] <- 0 # raw time: not a predictor
pr[c("time_yr","dead","na_cumhaz","wt","SDMVSTRA"),] <- 0 # complete: never imputed
imp <- mice(idat, m=10, predictorMatrix=pr, printFlag=FALSE, seed=604)
fits <- lapply(1:10, function(i){
fk <- fast; fk[fk$in_analysis, COV] <- complete(imp,i)[,COV]
svycoxph(f1, design=subset(svydesign(ids=~SDMVPSU,strata=~SDMVSTRA,weights=~wt,
data=fk,nest=TRUE), in_analysis))
})
po <- summary(MIcombine(fits))["groupIV",]Multiple imputation results:
MIcombine.default(fits)
results se (lower upper) missInfo
groupII 0.82193871 0.741614630 -0.63159925 2.27547668 0 %
groupIII 0.84556803 0.741224004 -0.60720432 2.29834039 0 %
groupIV 1.09509245 0.745671891 -0.36639760 2.55658250 0 %
age 0.06764454 0.005422335 0.05701696 0.07827212 0 %
female -0.44538653 0.129727191 -0.69964716 -0.19112591 0 %
dyslip 0.09585611 0.135392950 -0.16950919 0.36122142 0 %
htn 0.11525235 0.229042886 -0.33366346 0.56416815 0 %
t2dm 0.29255642 0.174595221 -0.04964393 0.63475676 0 %
ckd 0.88356024 0.120471499 0.64744044 1.11968004 0 %
mi 0.53979621 0.151724807 0.24242105 0.83717137 0 %
cancer 0.11064838 0.136939861 -0.15774881 0.37904558 0 %
on_htn_med -0.12608022 0.185987246 -0.49060852 0.23844809 0 %
on_dm_med 0.12184895 0.164464946 -0.20049642 0.44419432 0 %
on_lipid_med -0.21377285 0.166595894 -0.54029480 0.11274910 0 %cat(sprintf("weighted Model-1, MI (m=10): IV = %.2f (%.2f-%.2f)\n",
exp(po[["results"]]), exp(po[["(lower"]]), exp(po[["upper)"]])))weighted Model-1, MI (m=10): IV = 2.99 (0.69-12.89)cat(sprintf("fraction of missing information: %s\n", po[["missInfo"]]))fraction of missing information: 0 %Step 4b — The same four approaches, on the locked estimand
Everything above is reported against Group I, the paper’s reference, because that is the contrast the published analysis reports. The term’s locked question is IV vs III. Re-levelling inside the same model gives it, and the missingness has to be checked there too rather than assumed to behave the same way.
rl3 <- function(x){ x$group <- relevel(x$group, "III"); x }
desL <- function(fk) subset(svydesign(ids=~SDMVPSU, strata=~SDMVSTRA, weights=~wt,
data = rl3(fk), nest = TRUE), in_analysis)
fitL <- function(val, lbl){
fk <- fast
if (!is.na(val)) fk$cancer[fk$in_analysis & is.na(fk$cancer)] <- val
f <- svycoxph(f1, design = desL(fk))
data.frame(scenario = lbl, HR = fm(gv(f)), n = f$n)
}
fitsL <- lapply(1:10, function(i){
fk <- fast; fk[fk$in_analysis, COV] <- complete(imp, i)[, COV]
svycoxph(f1, design = desL(fk))
})
poL <- summary(MIcombine(fitsL))["groupIV", ]Multiple imputation results:
MIcombine.default(fitsL)
results se (lower upper) missInfo
groupI -0.84556803 0.741224004 -2.29834039 0.60720432 0 %
groupII -0.02362932 0.173237813 -0.36316920 0.31591055 0 %
groupIV 0.24952442 0.186700778 -0.11640238 0.61545122 0 %
age 0.06764454 0.005422335 0.05701696 0.07827212 0 %
female -0.44538653 0.129727191 -0.69964716 -0.19112591 0 %
dyslip 0.09585611 0.135392950 -0.16950919 0.36122142 0 %
htn 0.11525235 0.229042886 -0.33366346 0.56416815 0 %
t2dm 0.29255642 0.174595221 -0.04964393 0.63475676 0 %
ckd 0.88356024 0.120471499 0.64744044 1.11968004 0 %
mi 0.53979621 0.151724807 0.24242105 0.83717137 0 %
cancer 0.11064838 0.136939861 -0.15774881 0.37904558 0 %
on_htn_med -0.12608022 0.185987246 -0.49060852 0.23844809 0 %
on_dm_med 0.12184895 0.164464946 -0.20049642 0.44419432 0 %
on_lipid_med -0.21377285 0.166595894 -0.54029480 0.11274910 0 %locked <- rbind(
fitL(NA, "complete case (109 dropped)"),
fitL(0, "all 109 assigned NO cancer"),
fitL(1, "all 109 assigned CANCER"),
data.frame(scenario = "MI (m = 10)", n = fitsL[[1]]$n,
HR = sprintf("%.2f (%.2f-%.2f)", exp(poL[["results"]]),
exp(poL[["(lower"]]), exp(poL[["upper)"]]))))Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = rl3(fk), nest = TRUE), in_analysis)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = rl3(fk), nest = TRUE), in_analysis)
Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = rl3(fk), nest = TRUE), in_analysis)locked[, c("scenario","HR","n")]| scenario | HR | n |
|---|---|---|
| complete case (109 dropped) | 1.29 (0.89-1.86) | 5939 |
| all 109 assigned NO cancer | 1.28 (0.89-1.85) | 6048 |
| all 109 assigned CANCER | 1.28 (0.89-1.85) | 6048 |
| MI (m = 10) | 1.28 (0.89-1.85) | 6048 |
The locked contrast behaves as the replication contrast does: the four approaches return estimates that differ by at most 0.01, far less than the width of any of the intervals. Note that this is a statement about this missingness, not a general result — it had to be computed on the locked contrast to be claimed about the locked contrast.
Step 5 — What this actually shows
cc_w <- gv(fit_w); mi_w <- c(exp(po[["results"]]), exp(po[["(lower"]]), exp(po[["upper)"]]))Stratified 1 - level Cluster Sampling design (with replacement)
With (184) clusters.
subset(svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wt,
data = fast, nest = TRUE), in_analysis)cat(sprintf(" complete case (weighted): %s\n imputed (weighted): %s\n difference in HR: %.3f\n",
fm(cc_w), fm(mi_w), mi_w[1]-cc_w[1])) complete case (weighted): 2.98 (0.69-12.82)
imputed (weighted): 2.99 (0.69-12.89)
difference in HR: 0.012co <- sapply(fits, function(f) coef(f)[["groupIV"]])
cat(sprintf(" between-imputation variance of the Group IV log-HR: %.3e
", var(co))) between-imputation variance of the Group IV log-HR: 3.331e-11cat(sprintf(" distinct fitted values across m=10: %d
", length(unique(round(co,10))))) distinct fitted values across m=10: 7The missingness is real and correcting it changes essentially nothing. One covariate, 1.80% of the sample, and a hazard ratio that moves in the third significant figure.
Be exact about the fraction of missing information: it prints as 0 %, but that is approximately zero at the displayed precision, not literally zero. The between-imputation variance is 3.3e-11 — the imputed cancer values really do differ between completions, they simply have almost no leverage on the Group IV coefficient. Reporting “0%” without that qualification would claim more than the object supports.
It also means the shortcut flagged above is immaterial here: entering the weight and stratum as predictors rather than modelling the design inside the imputation changes the estimate not at all. That is a statement about this dataset with 1.8% missingness on one covariate, not a general licence.
The conclusion, stated as the evidence supports it
Within the locked domain — the 6,048 rows P5 analyses — only 1.8% of records are incomplete on the paper’s Model-1 covariates, entirely because of missing cancer status. Complete-case analysis, extreme-value stress tests, and pragmatic multiple imputation yield materially similar estimates. In this dataset, handling this residual missingness does not alter the substantive inference; accounting for the complex survey design has substantially greater impact.
Note the two things that sentence is careful not to say. It does not say missingness is generally unimportant — it says it is numerically negligible here. And it does not say survey design is the paper’s dominant threat; it says design mattered more than missingness among the analyses run in P4 and P5. Which threat most endangers the paper’s conclusion is a judgment that needs all five increments, and it is made at M1, not here.
Methodological importance is measured, not presumed. A sensitivity analysis that returns “no material change” has done its job.
What does move the conclusion is the design, not the missing data. Compare the two complete-case rows in Step 3: unweighted, the adjusted interval excludes 1; correctly weighted, it does not. Respecting the sampling design changes what you can claim; imputing 1.8% of one covariate does not.
A sensitivity analysis we deliberately did not run. Adding income (PIR, ~9% missing) and systolic BP would raise complete-case loss to about 13% and make the missing-data exercise look substantial. Neither is in the paper’s Model 1. Adding variables to manufacture a missing-data problem changes the adjustment set and answers a different question. If you want them, introduce them as a separately labelled sensitivity model with its own estimand — not inside this comparison.
Results for M1
lk <- function(s) endash(locked$HR[locked$scenario == s])
emit_results("P5", list(
n_missing = fmtn(sum(is.na(dom$cancer))),
pct_uw = sprintf("%.2f", 100 * mean(!cc)),
pct_w = sprintf("%.2f", 100 * sum(dom$wt[!cc]) / sum(dom$wt)),
n_vars_miss = fmtn(sum(miss > 0)),
n_vars = fmtn(length(COV)),
locked_cc = lk("complete case (109 dropped)"),
locked_no = lk("all 109 assigned NO cancer"),
locked_yes = lk("all 109 assigned CANCER"),
locked_mi = lk("MI (m = 10)"),
n_cc = fmtn(locked$n[locked$scenario == "complete case (109 dropped)"]),
n_full = fmtn(locked$n[locked$scenario == "MI (m = 10)"])
), list(locked_estimand = locked[, c("scenario","HR","n")],
replication_estimand = data.frame(
scenario = stress$scenario,
HR = endash(sprintf("%.2f (%.2f-%.2f)", stress$HR, stress$lo, stress$hi)),
n = stress$n)))wrote ../_results/P5.json (11 facts, 2 tables)