# Metabolic score for visceral fat (METS-VF) and osteoarthritis, NHANES 2009-2018.
# Xue et al. (2024), BMC Public Health, doi:10.1186/s12889-024-19722-0.
# Headline: OR 2.33 (95% CI 1.65-3.28), METS-VF quartile 4 vs quartile 1, Model 3 (Table 2, the main
# text's only regression table; the per-unit estimate is in a supplement).

# METS-VF (Bello-Chavolla et al.), as the paper prints it, with glucose, triglycerides and HDL in mg/dL.
mets_vf_published <- function(glucose, tg, hdl, bmi, waist, height, male, age) {
  mets_ir <- log(2 * glucose + tg) * bmi / log(hdl)
  4.466 + 0.011 * log(mets_ir)^3 + 3.239 * log(waist / height)^3 + 0.319 * male + 0.594 * log(age)
}

# What the paper's own numbers show it computed (variants): METS-IR without its outer log and with
# glucose in mmol/L (LBDGLUSI). On our analytic sample this gives Table 1's mean METS-VF (9.54
# without osteoarthritis and 10.14 with, against the published 9.54 and 10.21; the printed formula
# gives 6.73 and 7.15), Supplementary Table 1's profiles by quartile (mean age 38, 48, 51, 53 against
# 38, 48, 51, 52; osteoarthritis 7, 13, 17, 20% against 6, 12, 16, 20%), and Supplementary Table 2's
# per-unit odds ratios (1.31, 1.21, 1.22 against 1.35, 1.23, 1.23).
mets_vf_as_computed <- function(glucose, tg, hdl, bmi, waist, height, male, age) {
  mets_ir <- (2 * glucose * 0.0555 + tg) * bmi / log(hdl)
  4.466 + 0.011 * log(mets_ir)^3 + 3.239 * log(waist / height)^3 + 0.319 * male + 0.594 * log(age)
}

row295_model3 <- osteoarthritis ~ mets_vf_q + age + sex + race + education3 + pir3 + alcohol12 + smoking3 + work +
  recreation + calcium + vitamin_d + hypertension + diabetes + chd + energy

row295_build <- function(cycle) {
  demo <- demographics(cycle)
  d <- merge_all(demo, body_measures(cycle)[, c("SEQN", "bmi", "waist", "height")], fasting_glucose(cycle),
                 triglycerides(cycle)[, c("SEQN", "tg")], hdl_cholesterol(cycle), arthritis(cycle),
                 dietary_totals(cycle, "KCAL", days = 2), smoking(cycle), alcohol(cycle),
                 biochemistry(cycle)[, c("SEQN", "calcium")], component("VID", cycle, "LBXVIDMS"),
                 component("BPQ", cycle, "BPQ020"), component("DIQ", cycle, "DIQ010"), component("MCQ", cycle, "MCQ160C"))
  male <- d$sex == 1
  d$mets_vf <- mets_vf_published(d$glucose, d$tg, d$hdl, d$bmi, d$waist, d$height, male, d$age)
  d$mets_vf_computed <- mets_vf_as_computed(d$glucose, d$tg, d$hdl, d$bmi, d$waist, d$height, male, d$age)
  # Variant: the 2009-2010 type question (MCQ191: 1 rheumatoid, 2 osteoarthritis) read with the
  # later codes (MCQ195: 1 osteoarthritis), so rheumatoid arthritis counts as osteoarthritis.
  later_codes <- has(component("MCQ", cycle), "MCQ195")
  d$osteoarthritis_misread <- if (later_codes) d$osteoarthritis else d$rheumatoid
  # Energy: the mean of the two 24-hour recalls ("the average consumption of the two recalls").
  d$energy <- d$KCAL
  implausible <- (male & (d$energy < 500 | d$energy > 8000)) | (!male & (d$energy < 500 | d$energy > 5000))
  one_or_two <- dietary_totals(cycle, "KCAL", days = "mean_or_one")
  d$energy_one_or_two <- one_or_two$KCAL[match(d$SEQN, one_or_two$SEQN)]
  implausible_one_or_two <- (male & (d$energy_one_or_two < 500 | d$energy_one_or_two > 8000)) |
    (!male & (d$energy_one_or_two < 500 | d$energy_one_or_two > 5000))
  d$education3 <- factor(education3(d$education))
  d$pir3 <- cut(d$pir, c(-Inf, 1.3, 3.5, Inf), labels = c("<1.3", "1.3-3.5", ">3.5"))
  # At least 12 drinks in a year: ALQ101 through 2015-2016. 2017 on dropped that question; there
  # it is built from past-year drinking days (ALQ121) times drinks a drinking day (ALQ130).
  alq <- component("ALQ", cycle)
  if (has(alq, "ALQ101")) {
    d$alcohol12 <- yes(alq$ALQ101[match(d$SEQN, alq$SEQN)])
  } else {
    yearly <- d$drinking_days * d$drinks_per_day
    twelve <- (d$drinking_days >= 12) %in% TRUE | (yearly >= 12) %in% TRUE
    fewer <- d$alcohol3 %in% 1:2 | (yearly < 12) %in% TRUE
    d$alcohol12 <- ifelse(twelve, 1, ifelse(fewer, 0, NA))
  }
  d$smoking3 <- factor(d$smoking)
  paq <- component("PAQ", cycle)
  at <- function(v) if (has(paq, v)) paq[[v]][match(d$SEQN, paq$SEQN)] else NA
  activity <- function(vigorous, moderate) factor(ifelse(vigorous %in% TRUE, "vigorous", ifelse(moderate %in% TRUE, "moderate", "other")),
                                                  levels = c("other", "moderate", "vigorous"))
  # Work and leisure activity as Table 1 groups them (vigorous, moderate, other), from the yes/no
  # items; 2021-2023 asks how often leisure activity is done instead (more than 0 times is yes).
  d$work <- if (has(paq, "PAQ605")) activity(at("PAQ605") == 1, at("PAQ620") == 1) else NA
  times <- function(v) ifelse(at(v) %in% c(7777, 9999), NA, at(v))
  d$recreation <- if (has(paq, "PAQ650")) activity(at("PAQ650") == 1, at("PAQ665") == 1) else
    activity(times("PAD810Q") > 0, times("PAD790Q") > 0)
  d$vitamin_d <- d$LBXVIDMS
  d$hypertension <- yes(d$BPQ020)
  # Borderline diabetes (DIQ010 = 3) is neither yes nor no: the flow chart's 224 exclusions for
  # missing diabetes fit it being missing.
  d$diabetes <- yes(d$DIQ010)
  d$chd <- yes(d$MCQ160C)
  d$sex <- factor(d$sex)
  d$race <- factor(d$race)
  base <- !is.na(d$mets_vf) & d$age >= 20
  d$in_population <- base & !is.na(d$osteoarthritis) & !is.na(d$energy) & !(implausible %in% TRUE)
  d$in_population_one_or_two <- base & !is.na(d$osteoarthritis) & !is.na(d$energy_one_or_two) & !(implausible_one_or_two %in% TRUE)
  d
}

# Variant: energy from the first recall where there is no second (the flow chart's 558 exclusions
# for a "total energy deficit" are close to the number with no first-day recall).
row295_build_one_or_two <- function(cycle) {
  d <- row295_build(cycle)
  d$energy <- d$energy_one_or_two
  d$in_population <- d$in_population_one_or_two
  d
}

association <- list(
  id = "row295", row = 295, doi = "10.1186/s12889-024-19722-0",
  cycles = c("2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTSAF2YR",
  # The paper's own numbers show it computed METS-VF another way than it printed (see choices), so
  # the chosen version uses the index the published estimate came from.
  family = "logistic", term = "mets_vf_computed_qQ4",
  published = list(measure = "OR", estimate = 2.33, low = 1.65, high = 3.28, n = 7639, events = 937,
                   contrast = "METS-VF quartile 4 vs quartile 1 (quartiles of the analytic sample)"),
  left_out = c("alcohol intake as at least 12 drinks a year (ALQ101, asked through 2015-2016; 2021-2023 asks only whether one ever drank, how often, and how much)",
               "work activity (PAQ605, PAQ620; 2021-2023 asks only about leisure-time activity)"),
  build = row295_build,
  # Quartile cutpoints (unpublished): unweighted quartiles of the analytic sample of Model 3.
  constants = function(data) {
    variables <- setdiff(all.vars(row295_model3), "mets_vf_q")
    p <- data$in_population %in% TRUE & stats::complete.cases(data[, variables])
    list(mets_vf_quartiles = unname(stats::quantile(data$mets_vf[p], c(0.25, 0.5, 0.75))),
         mets_vf_computed_quartiles = unname(stats::quantile(data$mets_vf_computed[p], c(0.25, 0.5, 0.75))))
  },
  derive = function(data, constants) {
    quartile <- function(x, at) cut(x, c(-Inf, at, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))
    if (!is.null(constants)) {
      data$mets_vf_q <- quartile(data$mets_vf, constants$mets_vf_quartiles)
      data$mets_vf_computed_q <- quartile(data$mets_vf_computed, constants$mets_vf_computed_quartiles)
    }
    data
  },
  formula = update(row295_model3, . ~ . - mets_vf_q + mets_vf_computed_q),
  formula_harmonized = osteoarthritis ~ mets_vf_computed_q + age + sex + race + education3 + pir3 + smoking3 + recreation +
    calcium + vitamin_d + hypertension + diabetes + chd + energy,
  variants = list(
    list(label = "METS-VF as printed (mg/dL, with the outer log), Model 3", formula = row295_model3, term = "mets_vf_qQ4"),
    list(label = "METS-VF as printed, crude", formula = osteoarthritis ~ mets_vf_q, term = "mets_vf_qQ4"),
    list(label = "METS-VF as printed, age and sex", formula = osteoarthritis ~ mets_vf_q + age + sex, term = "mets_vf_qQ4"),
    list(label = "crude (published 3.90, 2.94-5.15)", formula = osteoarthritis ~ mets_vf_computed_q),
    list(label = "Model 2, age and sex (published 2.27, 1.67-3.08)", formula = osteoarthritis ~ mets_vf_computed_q + age + sex),
    list(label = "METS-VF as the paper computed it, per unit, Model 3 (supplement 1.23, 1.14-1.33)",
         formula = update(row295_model3, . ~ . - mets_vf_q + mets_vf_computed), term = "mets_vf_computed"),
    list(label = "unweighted", weighted = FALSE),
    list(label = "examination weight", weight = "WTMEC2YR"),
    list(label = "2009-2010 arthritis type read with the later codes", formula = update(row295_model3, osteoarthritis_misread ~ . - mets_vf_q + mets_vf_computed_q)),
    list(label = "energy from one recall where there is no second", build = row295_build_one_or_two)
  ),
  # 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 Methods print METS-IR as Ln[(2 x fasting glucose + fasting triglycerides) x BMI] / Ln(HDL-C), but Table 1's mean METS-VF (9.54 without osteoarthritis, 10.21 with) and Table 2's odds ratios are reproduced by METS-IR with no logarithm in its numerator and glucose in mmol/L (9.54 and 10.14; Model 3, quartile 4 vs 1, 2.24 against 2.33), not by the printed formula in mg/dL (6.73 and 7.15; 3.18)."),
    list(kind = "weighting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's note calls its percentages survey-weighted, but they are the unweighted ones: non-Hispanic White 42.42% without osteoarthritis (41.4% unweighted here, 66.8% weighted) and diabetes 11.77% and 21.83% (11.7% and 22.1% unweighted, 8.6% and 16.3% weighted)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's Former and Never smoking rows are swapped: its 56.02% former and 23.89% never without osteoarthritis are 58.6% never and 22.8% former here, and its 46.64% former and 33.34% never with osteoarthritis are 48.5% never and 33.3% former."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Results give Model 3's quartile 2 odds ratio as 1.14 (1.02-2.02), where Table 2 has 1.43 (1.02-2.01)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The text excludes 233 participants under 20 after those without osteoarthritis data, but the arthritis questions (MCQ160A) are asked only from age 20, so no one under 20 remains at that step.")
  ),
  choices = list(
    list(choice = "headline", decision = "quartile 4 vs quartile 1 in Model 3 (Table 2)", reason = "the abstract gives no estimate; Table 2, the main text's only regression table, has quartiles only; the per-unit estimate (1.23, 1.14-1.33) is in Supplementary Table 2"),
    list(choice = "METS-VF formula", decision = "as the paper's numbers show it computed: METS-IR without its outer log and with glucose in mmol/L, then the printed METS-VF formula", reason = "the formula as printed (mg/dL) gives Table 1 means of 6.73 and 7.15 against the published 9.54 and 10.21, and a crude OR of 6.23 against 3.90; the computed index gives 9.54 and 10.14, Supplementary Table 1's age by quartile, and the crude, Model 3, and per-unit estimates (3.47, 2.24, 1.224 against 3.90, 2.33, 1.23). The published estimate comes from that index, so the replication uses it; the printed formula is fitted as a variant"),
    list(choice = "quartile cutpoints", decision = "unweighted quartiles of the index in the Model 3 analytic sample of the paper's cycles, kept for 2021-2023", reason = "the paper prints none"),
    list(choice = "osteoarthritis", decision = "told of arthritis (MCQ160A) and type osteoarthritis (MCQ191 = 2 in 2009-2010, MCQ195 = 1 from 2011); other types count as no; a type answered don't know or refused is missing", reason = "the paper's question wording; with the computed METS-VF, osteoarthritis by quartile is 7, 13, 17, 20% against Supplementary Table 1's 6, 12, 16, 20%, and reading 2009-2010 with the later codes fits the published estimates less well (variant)"),
    list(choice = "dietary energy", decision = "mean of two reliable recalls, required; under 500 or over 5,000 kcal (women) or 8,000 (men) excluded", reason = "the paper says all participants had two recalls and their average was used; requiring both gives 7,623 participants (7,379 with a positive fasting weight) against the published 7,639; allowing one recall gives 8,581 (variant)"),
    list(choice = "weights", decision = "fasting subsample weight (WTSAF2YR), each cycle's weight over five", reason = "unstated; NCHS's guidelines use the smallest subsample's weight, and METS-VF needs fasting glucose and triglycerides. The paper divides the 2-year weights by 2, a constant that changes no estimate or interval"),
    list(choice = "alcohol intake", decision = "at least 12 drinks in any one year (ALQ101) through 2015-2016; in 2017-2018, past-year drinking days (ALQ121) times drinks a drinking day (ALQ130) of 12 or more", reason = "the paper's definition; 2017-2018 no longer asks ALQ101. With this coding the share who drink (69.9% without osteoarthritis, 65.2% with) matches Table 1's (69.8%, 66.1%); ever having had a drink (ALQ111) would give about 77%, and the same estimate"),
    list(choice = "work and leisure activity", decision = "vigorous if any vigorous activity (PAQ605, PAQ650), moderate if any moderate (PAQ620, PAQ665), other otherwise, with missing answers as other", reason = "Table 1's groups; the flow chart excludes no one for missing activity, and the groups' shares match Table 1. In 2021-2023 leisure activity is vigorous or moderate if done more than 0 times (PAD810Q, PAD790Q)"),
    list(choice = "hypertension, diabetes, coronary heart disease", decision = "self-reported diagnosis (BPQ020, DIQ010, MCQ160C), yes or no, with borderline diabetes missing", reason = "unstated; the flow chart's counts of missing values (13, 224, 17) fit these items, and Table 1's diabetes share (11.8% and 21.8%) matches ours (11.7% and 22.1%)"),
    list(choice = "other covariates", decision = "education in three levels (DMDEDUC2 1-2, 3, 4-5); income to poverty ratio under 1.3, 1.3-3.5, over 3.5; smoking never, former, current; age, serum calcium, vitamin D and energy continuous", reason = "Table 1's categories; it reports the continuous ones as means")
  )
)