# METS-IR (metabolic score for insulin resistance) and serum ferritin in women aged 20-49, NHANES
# 2005-2010 and 2015-2018. Hao et al. (2022), Front Med, doi:10.3389/fmed.2022.925344.
# Headline: beta 0.29 (95% CI 0.14-0.44) ng/mL of ferritin per unit of METS-IR, Model 3 (Table 2).
#
# The paper's own numbers show how its analysis was run. Its flow chart (17,704 with ferritin, 6,193
# with METS-IR, 1,990 under 20, 4,182 in the end), Table 1's counts (all but the sugar row's) and its
# means and SDs of age, BMI, HDL, glucose, triglycerides, CRP and METS-IR are reproduced exactly by
# women aged 20-49 from 2005-2010 alone, with METS-IR computed from the biochemistry profile's
# non-fasting serum glucose and triglycerides (LBXSGL, LBXSTR) where the paper says fasting glucose
# and triglycerides; no participant from 2015-2018, the cycles it says it used, got a METS-IR. On that sample the published
# Model 1 and Model 2 estimates are reproduced to the second decimal unweighted, and Model 3 is
# reproduced when it leaves out METS-IR's own components (BMI, HDL, glucose, triglycerides), which
# the text says it adjusted for. The chosen version runs the analysis the paper describes (fasting
# glucose and triglycerides, all five cycles, ferritin converted across analyzers as it says) with
# the model and codings its numbers show; the first variant is the analysis as the authors ran it.

# Table 1 splits each dietary intake at its mean in the analytic sample, with missing as "Unclear".
# Energy, sugar and iron as printed. Fat at 68.69 g: Table 1's counts for fat (1,961 below) are the
# split at 68.69 g, the sample's mean fat, which Table 1 prints in its moisture row, while the fat
# row prints 73.01. The moisture row repeats the fat row's counts, so the moisture cut is not
# printed; 2,607.39 g is the mean two-day moisture intake in the paper's sample as reconstructed.
ROW096_DIET_CUTS <- c(KCAL = 1863.94, SUGR = 112.20, TFAT = 68.69, IRON = 13.95, MOIS = 2607.39)

ROW096_EARLY <- c("2005-2006", "2007-2008", "2009-2010")

# METS-IR (Bello-Chavolla et al. 2018): ln(2 x glucose + triglycerides) x BMI / ln(HDL), mg/dL.
row096_mets_ir <- function(glucose, tg, bmi, hdl) log(2 * glucose + tg) * bmi / log(hdl)

# Activity intensity from work and recreation, as Table 1's counts show it was built: vigorous if
# any vigorous work or leisure activity, else moderate if any moderate, else mild. 2005-2006 asked
# vigorous and moderate activity in the past 30 days (PAD200, PAD320) in place of the GPAQ's work
# and leisure items (PAQ605, PAQ620, PAQ650, PAQ665). 2021-2023 asks no work activity: NA there.
row096_activity <- function(cycle) {
  p <- component("PAQ", cycle)
  levels <- c("mild", "moderate", "vigorous")
  if (has(p, "PAQ605")) {
    vigorous <- p$PAQ605 %in% 1 | p$PAQ650 %in% 1
    moderate <- p$PAQ620 %in% 1 | p$PAQ665 %in% 1
  } else if (has(p, "PAD200")) {
    vigorous <- p$PAD200 %in% 1
    moderate <- p$PAD320 %in% 1
  } else {
    return(data.frame(SEQN = p$SEQN, activity = factor(rep(NA, nrow(p)), levels = levels)))
  }
  data.frame(SEQN = p$SEQN, activity = factor(ifelse(vigorous, "vigorous", ifelse(moderate, "moderate", "mild")), levels = levels))
}

# One cycle's frame. mets_from = "fasting" uses the fasting subsample's plasma glucose and
# triglycerides (GLU, TRIGLY), as the paper describes; "serum" uses the biochemistry profile's
# non-fasting values, as its numbers show it did. `keep_cycles` limits the population to some
# cycles. `deming` converts 2005-2008 ferritin (Hitachi 912) to the Elecsys 170 scale of 2009 on
# with NCHS's equation, log10(E170) = 0.989 log10(H912) + 0.049, as the paper says it did.
row096_build <- function(cycle, mets_from = "fasting", keep_cycles = NULL, deming = TRUE) {
  demo <- demographics(cycle)
  mcq <- component("MCQ", cycle)
  alq <- component("ALQ", cycle)
  d <- merge_all(demo, ferritin(cycle), fasting_glucose(cycle), triglycerides(cycle), hdl_cholesterol(cycle),
                 body_measures(cycle), crp(cycle),
                 dietary_totals(cycle, names(ROW096_DIET_CUTS), days = 2),
                 component("SMQ", cycle, "SMQ020"), component("DIQ", cycle, "DIQ010"), component("BPQ", cycle, "BPQ020"),
                 mcq[, c("SEQN", intersect(c("MCQ160L", "MCQ220", "MCQ092"), names(mcq))), drop = FALSE],
                 alq[, c("SEQN", intersect("ALQ101", names(alq))), drop = FALSE],
                 alcohol(cycle)[, c("SEQN", "alcohol3", "drinking_days")],
                 row096_activity(cycle))
  if (deming && cycle %in% c("2005-2006", "2007-2008")) d$ferritin <- 10^(0.989 * log10(d$ferritin) + 0.049)
  if (mets_from == "serum") {
    bio <- biochemistry(cycle)
    d$glucose_mets <- bio$glucose_serum[match(d$SEQN, bio$SEQN)]
    d$tg_mets <- bio$tg_serum[match(d$SEQN, bio$SEQN)]
  } else {
    # Fasting values count only for those who met the 8 to 24 hour fasting criteria; NCHS gives the
    # others a fasting-subsample weight of 0 while keeping their values in the file.
    glu <- component("GLU", cycle)
    fasted <- (glu[[weight_variable("WTSAF2YR", cycle)]][match(d$SEQN, glu$SEQN)] > 0) %in% TRUE
    d$glucose_mets <- ifelse(fasted, d$glucose, NA)
    d$tg_mets <- ifelse(fasted, d$tg, NA)
  }
  d$mets_ir <- row096_mets_ir(d$glucose_mets, d$tg_mets, d$bmi, d$hdl)
  # Table 1's "Mexican American" group (1,322) is RIDRETH1 1 and 2: Mexican American and other Hispanic.
  d$race4 <- factor(ifelse(d$race %in% 1:2, "hispanic", ifelse(d$race %in% 3, "white", ifelse(d$race %in% 4, "black", ifelse(d$race %in% 5, "other", NA)))),
                    levels = c("hispanic", "white", "black", "other"))
  d$education3 <- factor(education3(d$education))
  d$education5 <- factor(d$education)
  answer <- function(x, levels = c("no", "yes")) factor(ifelse(x %in% 1, "yes", ifelse(x %in% 2, "no", NA)), levels = levels)
  d$smoker <- answer(d$SMQ020)
  d$hypertension <- answer(d$BPQ020)
  d$liver <- answer(d$MCQ160L)
  d$cancer <- answer(d$MCQ220)
  d$diabetes <- factor(ifelse(d$DIQ010 %in% 1, "yes", ifelse(d$DIQ010 %in% 2, "no", ifelse(d$DIQ010 %in% 3, "borderline", NA))),
                       levels = c("no", "yes", "borderline"))
  d$transfusion <- with_unclear(if (has(d, "MCQ092")) ifelse(d$MCQ092 %in% 1, "yes", ifelse(d$MCQ092 %in% 2, "no", NA)) else rep(NA, nrow(d)), c("no", "yes"))
  # Drinking as the text defines it: alcohol at least 12 times in the past year (past-year frequency,
  # ALQ120Q/U through 2015-2016, ALQ121 at least once a month from 2017); never drinkers are "no".
  past_year <- ifelse((d$drinking_days >= 12) %in% TRUE, "yes", ifelse((d$drinking_days < 12) %in% TRUE | d$alcohol3 %in% 1, "no", NA))
  d$drinking_past_year <- with_unclear(past_year, c("no", "yes"))
  # Drinking as Table 1 counts it (2,478 yes, 1,270 no, 434 unclear): ALQ101, at least 12 drinks in
  # any one year. 2017-2018 no longer asks it, so the text's definition stands in there.
  d$drinking <- with_unclear(if (has(d, "ALQ101")) ifelse(d$ALQ101 %in% 1, "yes", ifelse(d$ALQ101 %in% 2, "no", NA)) else past_year, c("no", "yes"))
  split <- function(x, at) with_unclear(ifelse(x < at, "low", ifelse(x >= at, "high", NA)), c("low", "high"))
  d$energy <- split(d$KCAL, ROW096_DIET_CUTS[["KCAL"]])
  d$sugar <- split(d$SUGR, ROW096_DIET_CUTS[["SUGR"]])
  d$fat <- split(d$TFAT, ROW096_DIET_CUTS[["TFAT"]])
  d$iron <- split(d$IRON, ROW096_DIET_CUTS[["IRON"]])
  d$moisture <- split(d$MOIS, ROW096_DIET_CUTS[["MOIS"]])
  # Women 20-49 (the paper's stated range; 2017-2018 measured ferritin at every age and in men), with
  # ferritin and METS-IR. Categorical covariates missing for 10 or fewer were dropped, as Figure 1
  # lists (education, diabetes, smoking, liver condition, malignancy, hypertension).
  d$in_population <- d$sex == 2 & d$age >= 20 & d$age <= 49 & !is.na(d$ferritin) & !is.na(d$mets_ir) &
    !is.na(d$education) & d$DIQ010 %in% 1:3 & d$SMQ020 %in% 1:2 & d$MCQ160L %in% 1:2 & d$MCQ220 %in% 1:2 & d$BPQ020 %in% 1:2
  if (!is.null(keep_cycles)) d$in_population <- d$in_population & cycle %in% keep_cycles
  d
}

ROW096_MODEL3 <- ferritin ~ mets_ir + age + race4 + pir_cov + education3 + smoker + drinking + activity + hypertension +
  diabetes + liver + cancer + energy + sugar + moisture + fat + iron + transfusion + crp_cov
# Model 3 as the text states it: every Table 1 covariate, METS-IR's components included.
ROW096_MODEL3_TEXT <- update(ROW096_MODEL3, . ~ . + bmi + hdl + glucose_mets + tg_mets)
# The analysis as the text describes it: fasting values, all five cycles, ferritin converted.
ROW096_TEXT <- list(build = function(cycle) row096_build(cycle), weight = "WTSAF2YR",
                    cycles = c("2005-2006", "2007-2008", "2009-2010", "2015-2016", "2017-2018"))

association <- list(
  id = "row096", row = 96, doi = "10.3389/fmed.2022.925344",
  # The published estimate comes from the analysis as the paper ran it (its flow chart, Table 1, and
  # Models 1-3 are reproduced exactly that way): METS-IR from the biochemistry profile's non-fasting
  # serum glucose and triglycerides, 2005-2010 only (no 2015-2018 participant got a METS-IR in the
  # paper's data), ferritin not converted between analyzers. The text's version is a variant.
  cycles = ROW096_EARLY,
  weight = "WTMEC2YR", blood_file = "FERTIN",
  # The text says it used 2-year sample weights, but Models 1 and 2 are reproduced to the second
  # decimal only unweighted (variants), as is common in EmpowerStats papers.
  weighted = FALSE,
  family = "linear", term = "mets_ir",
  published = list(measure = "beta", estimate = 0.29, low = 0.14, high = 0.44, n = 4182, events = NULL,
                   contrast = "ng/mL of serum ferritin per unit of METS-IR"),
  left_out = c("activity intensity (2021-2023 asks no work activity, and the paper's intensity combines work and recreation)",
               "blood transfusion (MCQ092 is not asked in 2021-2023)",
               "drinking as Table 1 codes it, from ALQ101 (at least 12 drinks in any one year), which 2021-2023 does not ask; the harmonized version uses the text's definition, drinking at least 12 times in the past year, from the past-year frequency questions every cycle asks"),
  build = function(cycle) row096_build(cycle, mets_from = "serum", deming = FALSE),
  # The paper's missing-data rule: a continuous covariate missing for 10% or less of the sample takes
  # the sample mean; otherwise it is split (at its mean, as the paper split its diet variables) with
  # missing as its own group. PIR (6.4% missing in the paper's sample) and CRP are imputed.
  derive = function(data, constants) {
    population <- data$in_population %in% TRUE
    covariate <- function(x) {
      m <- mean(x[population], na.rm = TRUE)
      if (mean(is.na(x[population])) <= 0.10) return(ifelse(is.na(x), m, x))
      with_unclear(ifelse(x < m, "low", ifelse(x >= m, "high", NA)), c("low", "high"))
    }
    data$pir_cov <- covariate(data$pir)
    data$crp_cov <- covariate(data$crp)
    data
  },
  formula = ROW096_MODEL3,
  formula_harmonized = ferritin ~ mets_ir + age + race4 + pir_cov + education3 + smoker + drinking_past_year + hypertension +
    diabetes + liver + cancer + energy + sugar + moisture + fat + iron + crp_cov,
  variants = list(
    list(label = "Model 1 (published 0.54, 0.41-0.68)", formula = ferritin ~ mets_ir),
    list(label = "Model 2 (published 0.49, 0.36-0.63)", formula = ferritin ~ mets_ir + age + race4),
    list(label = "Model 1, weighted (MEC weight)", formula = ferritin ~ mets_ir, weighted = TRUE),
    list(label = "Model 2, weighted (MEC weight)", formula = ferritin ~ mets_ir + age + race4, weighted = TRUE),
    list(label = "Model 3 with METS-IR's components, as the text says", formula = ROW096_MODEL3_TEXT),
    list(label = "weighted (MEC weight)", weighted = TRUE),
    list(label = "education in DMDEDUC2's five levels", formula = update(ROW096_MODEL3, . ~ . - education3 + education5)),
    c(list(label = "as the text describes it: fasting values, all five cycles, ferritin converted"), ROW096_TEXT),
    c(list(label = "as the text describes it, weighted (fasting subsample weight)", weighted = TRUE), ROW096_TEXT),
    c(list(label = "as the text describes it, with METS-IR's components", formula = ROW096_MODEL3_TEXT), ROW096_TEXT),
    list(label = "fasting values, 2005-2010 only, no ferritin conversion", weight = "WTSAF2YR",
         build = function(cycle) row096_build(cycle, deming = FALSE)),
    list(label = "serum values, all five cycles, no ferritin conversion", cycles = ROW096_TEXT$cycles,
         build = function(cycle) row096_build(cycle, mets_from = "serum", deming = FALSE))
  ),
  # What the coding check found this paper's analysis to have computed otherwise than its text
  # says (the plan's section 6.6).
  departures = list(
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The paper computes METS-IR from fasting glucose and triglycerides, but its flow chart (6,193 with METS-IR) and Table 1's FPG and fasting TG means (89.92 and 95.00, 123.47 and 136.72 mg/dL) are those of the biochemistry profile's non-fasting serum values (LBXSGL, LBXSTR), which alone reproduce Models 1, 2, and 3 (0.544, 0.493, 0.298 against 0.54, 0.49, 0.29). Fasting values give 0.172 on the same cycles."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The paper says it used NHANES 2005-2010 and 2015-2018, but no 2015-2018 participant has a METS-IR in its analysis: the flow chart's 6,193 with METS-IR and 1,990 under 20 are exactly 2005-2010's. Serum values on all five cycles give 6,699 participants and 0.228."),
    list(kind = "model", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "Table 2's note says Model 3 adjusts for all the covariates in Table 1, which include BMI, HDL, FPG, and fasting TG, the components of METS-IR, but the published 0.29 (0.14-0.44) is reproduced only without them (0.298, 0.146-0.450); with them Model 3 gives 2.40."),
    list(kind = "weighting", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods say 2-year sample weights were used, but Models 1 and 2 are reproduced to every printed digit, intervals included, only unweighted (0.544, 0.412-0.677 and 0.493, 0.360-0.625 against 0.54, 0.41-0.68 and 0.49, 0.36-0.63)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The paper says 2005-2008 ferritin was converted to the Elecsys 170 scale with the Deming equation, but Table 1's mean ferritin (54.04 ng/ml) and the groups split at it (2,766 and 1,416) match the unconverted values."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods give the race groups as Mexican, White, Black, and other races, but Table 1's 'Mexican American' count (1,322) is Mexican American and other Hispanic participants together, and Model 3 adjusts for all the covariates in Table 1."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's fat row prints its cut as 73.01 g, but its counts (1,961 below) are the split at 68.69 g, the sample's mean fat intake; the moisture row prints 68.69 g with the fat row's counts, so the moisture cut is not printed."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The flow chart starts from 84,724 participants and excludes 67,020 without ferritin, but the five cycles it names hold 50,259 participants; the 17,704 with ferritin it leaves are reproduced exactly.")
  ),
  choices = list(
    list(choice = "glucose and triglycerides in METS-IR", decision = "the biochemistry profile's non-fasting serum glucose and triglycerides (LBXSGL, LBXSTR), as the paper's computation used them",
         reason = "the paper defines METS-IR with fasting values, but its flow chart and Table 1 (FPG 89.92 and 95.00, TG 123.47 and 136.72 mg/dL) match the serum values, and only they reproduce its Models 1-3; the published estimate comes from them, so the replication follows them (the text's version is a variant)"),
    list(choice = "cycles", decision = "2005-2010",
         reason = "the paper names 2005-2010 and 2015-2018, but no 2015-2018 participant got a METS-IR in its data (its 6,193 with METS-IR and 1,990 under 20 are exactly 2005-2010's)"),
    list(choice = "Model 3 covariates", decision = "every Table 1 covariate but METS-IR's components (BMI, HDL, glucose, triglycerides)",
         reason = "the text says all Table 1 covariates, but only leaving the components out reproduces the published 0.29 (0.14-0.44): 0.298 (0.146-0.450) on the paper's sample; with them it is 2.40"),
    list(choice = "weights", decision = "unweighted", reason = "Models 1 and 2 are reproduced to the second decimal only unweighted; the text names 2-year weights"),
    list(choice = "ferritin across analyzers", decision = "2005-2008 values converted to the Elecsys scale with NCHS's Deming equation",
         reason = "the paper says it converted them; its mean ferritin (54.04) and group sizes match unconverted values, and the conversion barely moves the estimate (variant)"),
    list(choice = "pregnancy", decision = "not excluded", reason = "no exclusion is listed, and the paper's 4,182 include 420 pregnant women"),
    list(choice = "age range", decision = "20-49", reason = "the paper reports ages 20-49 (Discussion; Table 1 means); ferritin was measured at every age in 2017-2018 and only to 49 otherwise"),
    list(choice = "education", decision = "less than high school, high school, more than high school",
         reason = "the text and Table 1; DMDEDUC2's five levels give the printed Model 3 to the digit on the paper's sample, 0.291 (0.138-0.444) against 0.298 (0.146-0.450), so its model may have used them (variant)"),
    list(choice = "race", decision = "Hispanic (RIDRETH1 1-2), White, Black, other", reason = "Table 1's 'Mexican American' count (1,322) is Mexican American and other Hispanic together"),
    list(choice = "drinking", decision = "ALQ101 (at least 12 drinks in any one year) with missing as its own level; in 2017-2018, which no longer asks it, the text's definition (at least 12 times in the past year)",
         reason = "Table 1's counts are ALQ101's"),
    list(choice = "activity intensity", decision = "vigorous if any vigorous work or leisure activity, else moderate if any moderate, else mild (PAD200 and PAD320 in 2005-2006)",
         reason = "reproduces Table 1's counts (1,462, 1,369, 1,351) exactly"),
    list(choice = "diet covariates", decision = "two-day means split at the sample mean (Table 1's cuts; fat at 68.69 g and moisture at 2,607.39 g) with missing as its own level",
         reason = "Table 1 shows means as cuts and 621 unclear, the count missing either day; its fat row prints 73.01 but its counts are the 68.69 split, and its moisture row repeats the fat row; its sugar row has 907 unclear where we find 621"),
    list(choice = "missing PIR and CRP", decision = "the sample mean when 10% or less is missing, else split at the mean with missing as its own level", reason = "the paper's stated rule")
  )
)