Reproduction Report

Paper: Kueh MTW, Intaran MAU, Goh R, et al. Body weight categories and fat distribution in relation to all-cause mortality among adults with metabolic dysfunction-associated steatotic liver disease: a population-based analysis of NHANES, 2007–2018. BMJ Open 2026;16:e113719.

Data: NHANES 2007–2018 (6 cycles) + NCHS Public-Use Linked Mortality Files (2019 release), all public. Software: R 4.5.1.

Scope note. These materials document a teaching reproduction of published analyses using the specified public data and course implementation. Differences from published results are reproducibility findings, not allegations of research misconduct, and may reflect analytic, implementation, or data-version differences. The findings below stand as recorded; this note frames how to read them, not whether they hold.

Verdict: Table 1, Table 2 and Figure 1 reproduce closely, though not exactly, from the public data. The analytic cohort is 1.1% larger than the paper’s (6,371 against about 6,300; the gap is unresolved, §6), group sizes differ by at most 2.3% (Group III, 787 against 769), there is one extra death (586 against 585), Group IV’s numbers at risk are within about 1% (§4), and the hazard ratios compared in §3 agree to within a few percent. Figure 2 reproduces for BMI (panel a) but not for WHtR (panel b), where our curve turns down above WHtR ≈ 0.75 while the published curve keeps rising; the paper does not report the spline specification, so the difference cannot be resolved from what is published. Six reproducibility findings surfaced along the way (alcohol exclusion, BMI cut-point labelling, CV-death definition, FLI circularity, survey design, the WHtR spline) and are documented below — these are the pedagogically useful part.

1. Sample funnel — typed, and two row sets kept apart

Every step carries a type label. A funnel shows N falling; it does not show why, and the reasons are not interchangeable.

After eligibility, the two largest drops are Step 2 (design) and Step 4 (target population), and neither is missing data. Step 2 is the NHANES fasting subsample: it was selected deliberately, and WTSAF2YR exists to account for that selection stage and its nonresponse, so the correct response is the right weight and a domain analysis (P4), not imputation. Step 4 is people without hepatic steatosis, who are not missing — they are not who the paper is about. Genuine missingness here is 792 item-nonresponse losses at Step 3 and, separately, 2 outcome-ascertainment losses at Step 7: different problems, not added together, and not 30,532, which is what reading every post-age-eligibility reduction as attrition would give.

Two row sets, kept apart.

The two differ by 323 records and neither is an error. Dropping Step 2 lets 819 adults through Step 3 who have all four FLI inputs but no valid fasting weight — that is why the full analytic file’s own FLI step is 14,989 rather than 14,170 — and 323 of those 819 survive to the end. The paper is unweighted, so 6,371 is what reproduces it; 6,048 is a different (correct) denominator, not a correction to the first. Quote 6,371 when you mean what reproduces the paper, and 6,048 when you mean what a design-aware estimate can be computed on. logs/sample_funnel.csv logs both chains, typed, and closes with the Alcohol exclusion applied = FALSE switch, which is why the full analytic file stays at 6,371 rather than dropping further.

Fatty Liver Index (Bedogni 2006): FLI = e^L / (1+e^L) × 100, L = 0.953·ln(TG) + 0.139·BMI + 0.718·ln(GGT) + 0.053·WC − 15.745, steatosis if FLI ≥ 60.

Phenotype groups: obesity = BMI ≥ 30 (≥ 25 for Asians); high central adiposity = WHtR ≥ 0.6.

2. Table 1 — baseline characteristics (Overall column)

All group-stratified columns and p-values are in output/tables/table1_reproduced.{csv,html}.

Comorbidity definitions (calibrated to the paper’s prevalences — the paper uses composites, not pure self-report): - Hypertension = self-report or measured BP ≥ 140/90 - T2DM = self-report or HbA1c ≥ 6.5% - Dyslipidaemia = self-report or lipid-lowering meds or LDL ≥ 130 - MI = prior heart attack or coronary heart disease (MCQ160E ∪ MCQ160C) - CKD = eGFR < 60 (CKD-EPI 2021) or ACR ≥ 30 - Current smoking is expressed as a % of ever-smokers (matches the paper’s 40.6%)

3. Table 2 — Cox regression (HR, 95% CI)

Every reproduced estimate in this section is unweighted and drawn from the full analytic file (6,371 rows, §1), with each group compared against Group I. The unadjusted fits use all 6,371 rows; the Model 1 and Model 2 fits are complete case: they keep only the rows with every covariate in that model, a little under 6,371 (04_table2_cox.R does not print the exact N).

All-cause mortality

Cardiovascular mortality (CV death = underlying cause = diseases of heart)

Event counts match exactly except all-cause Group II (389 against 388, so 586 deaths against 585); the cardiovascular counts match exactly (158). The unadjusted HRs agree to the 2nd–3rd significant figure; the Model 2 HRs differ by up to about 3% (Group II, 2.37 against 2.30). Full three-model table in output/tables/table2_reproduced.csv.

• Model 1: age, sex, dyslipidaemia, hypertension, T2DM, CKD, MI, cancer, antihypertensive / antidiabetic / lipid-lowering medication use.
• Model 2: Model 1 + current smoking + sedentary minutes/day.

4. Figures

• output/figures/figure1_km.png — Kaplan–Meier by group. Group IV falls to ≈ 0.80 survival at ~7 y; Group I stays highest; log-rank p < 0.0001. Number-at-risk matches (e.g. Group IV 608 / 549 / 431 / 330 vs paper 604 / 546 / 433 / 334).
• output/figures/figure2_rcs.png — restricted cubic splines (4 knots, unadjusted, probability scale, matching the published caption). Panel (a) BMI reproduces: a steep fall from ≈ 0.22 at BMI 22 to ≈ 0.08 by BMI 33, then a shallow tail — the paper’s “J-shaped trend towards reduced all-cause mortality as BMI increased” (nonlinear p = 0.003). Panel (b) WHtR does not reproduce at the upper tail. The published curve rises and then flattens, staying highest at WHtR 1.0 (≈ 0.045 → ≈ 0.12). Ours rises to a peak of ≈ 0.11 near WHtR 0.75 and then declines to ≈ 0.086 by 1.0 — an inverted U, not a monotonic increase. The overall association is significant (p = 0.004) but the non-linearity is not (p = 0.073), and the confidence band above WHtR ≈ 0.85 is very wide. See finding #6.

5. Key reproduction findings (the interesting part)

1. The Table 1/2 column headers and the Methods text state different BMI cut-points. The headers read “BMI ≥ 25 / < 25”, while the reported BMI medians in the non-obese groups are 27.6 and 28.8, which a 25 cut cannot produce. The Methods text gives obesity = BMI ≥ 30, and coding to 30 reproduces the published group sizes; coding to 25 does not. We therefore read the analysis as using 30 and the printed headers as a labelling inconsistency.

2. The stated alcohol exclusion does not match the reported numbers. The paper says it excluded “significant alcohol” (> 2 drinks/day men, > 1/day women), yet its N = 6,300 and 585 deaths match a cohort with no material alcohol exclusion (our 6,371, 586 deaths). A faithful average-daily exclusion removes ~430 people (N → 5,944). We keep the paper-matching cohort as primary and expose the exclusion as a sensitivity switch (APPLY_ALCOHOL_EXCL in 02_build_analytic.R). We could not reconcile the stated exclusion with the reported N under any coding of the public alcohol variables we implemented; a different implementation of the same sentence may reconcile them.

3. CV mortality = “diseases of heart” only. Using heart + cerebrovascular gives 185 CV deaths; heart-only gives exactly the paper’s 158. Table 2 uses the heart-only definition, built inside 04_table2_cox.R. The cv_death variable saved in masld_analytic.rds is still the broader heart + cerebrovascular one (185), so it does not reproduce Table 2 if used directly.

4. FLI circularity (Reviewer 1’s point). FLI uses BMI and waist circumference, which also define the BMI/WHtR exposure groups — exposure and steatosis ascertainment share inputs. The paper’s TyG-index sensitivity analysis addresses this; it is one of the main methodological threats to interrogate, and one of the candidates for the dominant threat discussed in §7.

5. Complex survey design not applied. The paper (and this reproduction) treat NHANES as a simple cohort — no svydesign, weights, strata, or PSUs. Defensible for internal HR estimation but affects national representativeness of the descriptives.

6. The WHtR spline does not reproduce, and the specification needed to check it is not reported. The published Figure 2b rises across the whole WHtR range and stays highest at the top; ours peaks near WHtR 0.75 and turns down (see §4). The Methods say only that “restricted cubic splines were constructed” — no knot count, no knot placement, no statement of whether the model was adjusted, and no trimming rule for the sparse upper tail. Those are exactly the choices that determine the shape where the data run out. Our panels use rcs(x, 4) with the rms default knot placement, so a different knot count, knot placement or trimming rule is a plausible explanation for the divergence where the data thin out. This is the clearest case in the paper of a published result that cannot be checked from what is reported, and it is a better worked example of under-specification than anything above.

6. Assumptions / deviations from an unknown ground truth

• Medication classes for the Cox models use NHANES questionnaire self-report (BPQ/DIQ) as a proxy; GLP-1 / SGLT2i-specific flags would need the RXQ_RX drug file parsed against the Multum lexicon (the lexicon file 404’d at download and is not required for the primary results).
• Residual ±1–3 pt differences in HTN/T2DM/dyslipidaemia/CKD reflect the paper’s under-specified composite definitions (“as previously described”).
• N differs by +71 (1.1%), and the gap is unresolved. The alcohol exclusion does not explain it: applying the exclusion overshoots, taking N to 5,944, well below the paper’s approximately 6,300 (§5, finding #2). Reading the paper’s age criterion strictly as > 18 gives 6,301 with 585 deaths, close to the published totals though not to every group cell, but the paper’s Table 1 mean age, 51.0 (16.6), matches our ≥ 18 cohort instead (§2). The evidence conflicts, and the +71 cannot be attributed to either reading from what is published.
• Blood pressure follows the documented NHANES averaging protocol, not a plain mean. The Blood Pressure Procedures Manual (s3.1.7.3) excludes the first reading whenever more than one was obtained, treats an all-zero diastolic set as zero, and omits a zero that sits alongside a positive diastolic reading. The published paper does not specify its BP averaging implementation in enough detail to establish that it used the same algorithm, so this is a documentation-informed analytic choice, not a claim of exact computational reproduction - the same status as the alcohol-exclusion switch above. Correcting it changed the derived diastolic BP for 5,523 of 6,371 participants (and the systolic BP for a similar number) and moved hypertension from 3,444 to 3,448. Cohort N, phenotype cell counts and deaths did not move, and the Model 1 and Model 2 hazard ratios shifted only in the third significant figure (Group IV Model 2: 2.91 (1.18-7.20) -> 2.91 (1.18-7.21)). A coding error can be real and worth fixing even when the headline estimate turns out to be robust to it.

7. Map to the course project increments (P1–P5)

This reproduction is a worked example of the term project spine — and at two of the five increments, a worked example of what the increment has to repair:

Naming the single dominant threat is not a weekly increment; it is the analytic-judgment core of M1. The candidates this reproduction surfaces are FLI circularity (#4), the survey-design omission (#5), and the under-specified WHtR spline (#6).

What this reproduction does and does not establish. The pipeline reproduces the paper’s analysis, and the discrepancies above are documented rather than tuned away. Read crudely, Group IV has the poorest observed survival, and the paper reports exactly that — that descriptive claim does follow. Read with the paper’s Model 2 adjustment, the picture is different: the Model 2 estimates for Groups III and IV are near enough to each other to be indistinguishable (repro 3.00 vs 2.91, §3; paper 3.03 vs 2.89; all unweighted and against Group I), so the Model 2 results do not provide clear evidence that Group IV has higher mortality than Group III. Estimated directly as Group IV against Group III, that contrast is compatible with no difference as well as with a meaningful increase — 1.05 (0.77–1.45) unweighted and 1.29 (0.89–1.86) design-weighted (WTSAF2YR/6, strata and PSU; course increments P2 and P4). These two are adjusted for Model 1, not Model 2, and both are fitted on the 5,939 locked-domain rows complete on the Model 1 covariates, so they differ from each other only in the weighting. They belong to the course’s paper replication (the paper’s four groups and covariates); the Week 4–7 lectures adjust for a different covariate set, so their estimates for related contrasts differ.

So the replication succeeds, and the disagreement is not with the authors’ arithmetic but with which contrast supports the stated conclusion. Do not report this reproduction as confirming the stronger, comparative reading — that the Group IV phenotype carries higher mortality than Group III, and does so independently of overall body weight. Neither the published table nor this pipeline establishes that.