# Blood manganese and MAFLD below the turning point log10 Mn = 1.10, NHANES 2015-2020 (in fact
# 2017-March 2020). Tang et al. (2024), Front Public Health, doi:10.3389/fpubh.2024.1280163.
# Headline: OR 3.936 (95% CI 2.631-5.887) per unit of log10 whole-blood Mn below log10 Mn = 1.10
# (Table 3, threshold-effect model). The published n (8,542) and events (4,370) are the whole
# sample's; in the authors' data 7,131 (3,624 with MAFLD) are below the knot.
#
# The authors' supplementary dataset shows how the estimate was made: an unweighted logistic model
# of MAFLD on log10 Mn with no covariates, fitted on the participants below the knot (each segment
# fitted separately), on the 2015-2016, 2017-2018 and 2017-March 2020 files stacked. The 2017-March
# 2020 files already hold the 2017-2018 participants under new SEQNs, so about 3,100 people appear
# twice; 2015-2016 has no elastography, so no one from it is in the sample. Requiring chromium and
# cobalt, measured only at ages 40 and over, made the sample 40-80. The chosen version is the
# analysis on 2017-March 2020 without the duplicates; a variant stacks the files as the authors did.

ROW087_KNOT <- 1.10

# The ten blood metals the paper requires: lead, total mercury, cadmium, manganese and selenium
# (PBCD), inorganic, methyl and ethyl mercury (IHGEM), chromium and cobalt (CRCO, measured at ages
# 40 and over; not measured in 2021-2023).
row087_metals <- function(cycle) {
  pbcd <- component("PBCD", cycle, c("LBXBPB", "LBXTHG", "LBXBCD", "LBXBMN", "LBXBSE"))
  ihgem <- component("IHGEM", cycle, c("LBXIHG", "LBXBGM", "LBXBGE"))
  d <- data.frame(SEQN = pbcd$SEQN, pb = pbcd$LBXBPB, hg = pbcd$LBXTHG, cd = pbcd$LBXBCD, mn = pbcd$LBXBMN, se = pbcd$LBXBSE)
  at <- match(d$SEQN, ihgem$SEQN)
  d$inhg <- ihgem$LBXIHG[at]
  d$mehg <- ihgem$LBXBGM[at]
  d$ethg <- ihgem$LBXBGE[at]
  if (cycle == REPLICATION_CYCLE) {
    d$cr <- NA_real_
    d$co <- NA_real_
  } else {
    crco <- component("CRCO", cycle, c("LBXBCR", "LBXBCO"))
    d$cr <- crco$LBXBCR[match(d$SEQN, crco$SEQN)]
    d$co <- crco$LBXBCO[match(d$SEQN, crco$SEQN)]
  }
  d
}

# MAFLD as the paper defines it (Eslam et al. 2020): steatosis of any grade and metabolic
# dysfunction, which is BMI 25 or more, or type 2 diabetes (antidiabetic drugs, fasting glucose
# 7.0 mmol/L or more, HbA1c over 6.4%), or at least two of: waist over 102 cm (men) or 88 cm
# (women); blood pressure 130/85 or more or antihypertensive drugs; triglycerides 1.70 mmol/L or
# more or lipid-lowering drugs; HDL under 1.0 (men) or 1.3 mmol/L (women) or lipid-lowering drugs;
# prediabetes (fasting glucose 5.6-6.9 mmol/L or HbA1c 5.7-6.4%); HOMA-IR 2.5 or more; hs-CRP over
# 2 mg/L. Steatosis is CAP over 268 dB/m, with partial exams kept, as the authors' data show.
# MAFLD is known wherever CAP is: a criterion whose measurement is missing counts as not met, as in
# the authors' data (fasting glucose, triglycerides and insulin exist only for the fasting subsample).
row087_mafld <- function(cycle, demo) {
  d <- merge_all(demo[, c("SEQN", "sex")], elastography(cycle), body_measures(cycle), fasting_glucose(cycle), hba1c(cycle),
                 triglycerides(cycle), hdl_cholesterol(cycle), crp(cycle), blood_pressure(cycle), bp_questions(cycle),
                 diabetes_questions(cycle), component("INS", cycle, "LBXIN"))
  met <- function(x) x %in% TRUE
  male <- d$sex == 1
  glucose <- d$glucose * 0.05551  # mmol/L
  hdl <- cholesterol_mmol(d$hdl)
  lipid_drugs <- met(d$cholesterol_medication == 1)
  diabetes <- met(d$insulin == 1) | met(d$pills == 1) | met(glucose >= 7.0) | met(d$hba1c > 6.4)
  abnormalities <- met(ifelse(male, d$waist > 102, d$waist > 88)) +
    (met(d$sbp >= 130) | met(d$dbp >= 85) | met(d$bp_medication == 1)) +
    (met(tg_mmol(d$tg) >= 1.70) | lipid_drugs) +
    (met(ifelse(male, hdl < 1.0, hdl < 1.3)) | lipid_drugs) +
    (met(glucose >= 5.6 & glucose < 7.0) | met(d$hba1c >= 5.7 & d$hba1c <= 6.4)) +
    met(d$LBXIN * glucose / 22.5 >= 2.5) +
    met(d$crp > 2)
  dysfunction <- met(d$bmi >= 25) | diabetes | abnormalities >= 2
  data.frame(SEQN = d$SEQN, cap = d$cap, mafld = ifelse(is.na(d$cap), NA, as.integer(d$cap > 268 & dysfunction)))
}

# One cycle's frame. `segment` picks the participants the model is fitted on: below the knot (the
# headline), at or above it, or both (for the hinge model).
row087_build <- function(cycle, segment = "below") {
  demo <- demographics(cycle)
  d <- merge_all(demo, row087_metals(cycle), row087_mafld(cycle, demo))
  d$log_mn <- log10(d$mn)
  d$log_mn_above <- pmax(d$log_mn - ROW087_KNOT, 0)
  measured <- stats::complete.cases(d[, c("pb", "hg", "cd", "mn", "se", "inhg", "mehg", "ethg")])
  in_segment <- switch(segment, below = d$log_mn < ROW087_KNOT, above = d$log_mn >= ROW087_KNOT, both = !is.na(d$log_mn)) %in% TRUE
  # Adults with MAFLD status and all ten metals. 2021-2023 measures no chromium or cobalt, which
  # NHANES measured only at ages 40 and over, so there an age of 40 or more stands in for them.
  d$in_population <- d$age >= 20 & !is.na(d$mafld) & measured & !is.na(d$cr) & !is.na(d$co) & in_segment
  d$in_population_harmonized <- d$age >= 40 & !is.na(d$mafld) & measured & in_segment
  d
}

# The authors' stacking: the 2017-2018 files' participants added again, under their 2017-2018 SEQNs,
# to the 2017-March 2020 files' (which already hold them under new SEQNs). Variants are fitted on the
# association's cycles, so this build returns both sets for 2017-2020; the fit is unweighted, so the
# 2017-2018 rows need no weight.
row087_stacked <- function(cycle, segment = "below") {
  rows <- row087_build(cycle, segment)
  if (cycle != "2017-2020") return(rows)
  rbind(rows, row087_build("2017-2018", segment))
}

association <- list(
  id = "row087", row = 87, doi = "10.3389/fpubh.2024.1280163",
  cycles = "2017-2020",
  weight = "WTMEC2YR", blood_file = "PBCD",
  # The paper mentions no weights, strata or PSUs, and an unweighted fit of the authors' own data
  # reproduces Table 3 exactly.
  weighted = FALSE,
  family = "logistic", term = "log_mn",
  published = list(measure = "OR", estimate = 3.936, low = 2.631, high = 5.887, n = 8542, events = 4370,
                   contrast = "per unit of log10 whole-blood manganese (ug/L), among participants below log10 Mn = 1.10"),
  left_out = c("the requirement that chromium and cobalt were measured (2021-2023 does not measure them): the harmonized version requires an age of 40 or more instead, since NHANES measured them only at those ages and requiring them limited the paper's sample to them"),
  build = function(cycle) row087_build(cycle),
  formula = mafld ~ log_mn,
  variants = list(
    list(label = "2017-2018 and 2017-2020 files stacked, as the authors did", build = function(cycle) row087_stacked(cycle)),
    list(label = "stacked, at or above the knot (published 0.458, 0.121-1.726)", build = function(cycle) row087_stacked(cycle, "above")),
    list(label = "at or above the knot", build = function(cycle) row087_build(cycle, "above")),
    list(label = "continuous hinge at the knot, both segments", build = function(cycle) row087_build(cycle, "both"), formula = mafld ~ log_mn + log_mn_above),
    list(label = "weighted (MEC weight)", weighted = TRUE)
  ),
  # 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 = "sample", affects_headline = TRUE, followed = FALSE, evidence = "paper",
         detail = "The 2017-2018 files are stacked with the 2017-March 2020 files, which already hold their participants: in the authors' supplementary dataset 3,141 of the 3,269 analytic rows labeled 2017-2018 have an exact twin among those labeled 2019-2020, so the 8,542 rows hold about 5,400 people. Stacked, the model gives 4.055 on 7,131 rows below the knot (the authors' count); counted once, 4.160 on 4,431."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The paper says it analyzed adults aged 20 and over, but it required blood chromium and cobalt, which NHANES measured only at ages 40 and over, so its sample is aged 40 to 80 (Table 1's medians are 60 and 60; the rows with chromium in the authors' dataset are aged 40 to 80)."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The paper says it used the 2015-2020 files, but 2015-2016 has no elastography, so no 2015-2016 participant has MAFLD status and none is in the analysis, in the authors' dataset as here."),
    list(kind = "model", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods describe a piecewise linear regression, but Table 3's estimates are reproduced exactly from the authors' dataset by separate unadjusted, unweighted logistic models below and above log10 Mn = 1.10 (3.936, 2.631-5.887; 0.458, 0.121-1.726); a continuous hinge model with a shared intercept gives 3.546."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The abstract gives the turning point as log10 Mn = 1.10 and as Mn = 12.61 ug/L, but 10^1.10 is 12.59 ug/L; Table 3 is reproduced exactly from the authors' dataset with the split at log10 Mn = 1.10, which the file uses.")
  ),
  choices = list(
    list(choice = "cycles", decision = "2017-March 2020 files only",
         reason = "the authors stacked the 2017-2018 and 2017-March 2020 files, which counts about 3,100 people twice (variant); 2015-2016 has no elastography, so it adds no one"),
    list(choice = "model", decision = "logistic, MAFLD on log10 Mn with no covariates, fitted on those below log10 Mn = 1.10",
         reason = "Table 3 names no covariates; a refit of the authors' data this way reproduces its four estimates exactly, while a hinge model with a shared intercept does not (variant)"),
    list(choice = "weights", decision = "unweighted", reason = "the paper mentions no weights; the authors' data reproduce Table 3 unweighted"),
    list(choice = "steatosis", decision = "CAP over 268 dB/m, partial exams kept, no reliability filter",
         reason = "unstated; in the authors' data every MAFLD case has CAP 269 or more and MAFLD is known for everyone with a CAP value"),
    list(choice = "metabolic criteria", decision = "the paper's thresholds in mmol/L; drugs from the questionnaire (DIQ050, DIQ070, BPQ050A, BPQ100D); BP the mean of the readings; a missing measurement counts as not met",
         reason = "unstated; this agrees with the authors' MAFLD for 9,612 of the 9,698 participants of 2017-March 2020 with a CAP value (and 5,901 of 5,948 in their 2017-2018 rows); no single change to a criterion explains the rest"),
    list(choice = "metals required", decision = "all ten (lead, cadmium, total, inorganic, methyl and ethyl mercury, manganese, selenium, chromium, cobalt)",
         reason = "the paper excludes those missing heavy metal data; chromium and cobalt limit the sample to ages 40 and over")
  )
)