# Caffeine intake and obesity (BMI-for-age overweight or obese) in children and adolescents aged
# 2-19, NHANES 2011-2016 and 2017-March 2020.
# Liu and Cui (2024), PLOS ONE, doi:10.1371/journal.pone.0300566.
# Headline: OR 1.0020 (95% CI 1.0005-1.0036), labeled caffeine quartile 4 vs quartile 1, Model 3 (Table 2).
#
# Table 2's estimates are not quartile contrasts. Its "Quartile k" cells fit the odds ratio per
# milligram of caffeine among the participants in quartile k of unadjusted caffeine (quartile 1 is
# no caffeine, so it has no slope and prints as the reference), from unweighted models. The chosen
# version computes the headline that way (per milligram within quartile 4, unweighted); the
# analysis the text describes (quartile 4 against quartile 1 of energy-adjusted caffeine, weighted)
# is a variant, the contrast the paper's subgroup estimates appear to be (Table 4, Model 3: 1.60 in
# boys, 1.44 in girls).

# Two-day mean intakes, named after the DR1T/DR2T prefix as dietary_totals() takes them.
ROW066_DIET <- c(caffeine = "CAFF", energy = "KCAL", protein = "PROT", sugar = "SUGR", fat = "TFAT",
                 cholesterol = "CHOL", vitamin_b6 = "VB6", vitamin_c = "VC", vitamin_d = "VD", calcium = "CALC",
                 phosphorus = "PHOS", magnesium = "MAGN", zinc = "ZINC", copper = "COPP", sodium = "SODI",
                 potassium = "POTA", selenium = "SELE")

# The dietary covariates of Model 3 as the Table 2 note lists them.
ROW066_NUTRIENTS <- c("energy", "protein", "sugar", "fat", "cholesterol", "vitamin_b6", "vitamin_c", "calcium",
                      "phosphorus", "zinc", "copper", "sodium", "selenium")

ROW066_MODEL3 <- overweight_obese ~ caffeine_q + sex + age + race + pir + energy + protein + sugar + fat + cholesterol +
  vitamin_b6 + vitamin_c + calcium + phosphorus + zinc + copper + sodium + selenium

# Per milligram of caffeine in place of the quartiles, with Model 3's covariates. With energy in
# the model, raw and residual-adjusted caffeine give the same coefficient.
ROW066_PER_MG <- update(ROW066_MODEL3, . ~ . - caffeine_q + caffeine)

# Two-day mean nutrient intake from the 24-hour supplement recalls (DS1TOT, DS2TOT). A nutrient no
# supplement supplied is missing there and counts as zero, as does a day whose supplement count
# is missing. 2021-2023 has no 24-hour supplement recall (only the 30-day DSQTOT).
row066_supplements <- function(cycle, variables) {
  s1 <- component("DS1TOT", cycle)
  s2 <- component("DS2TOT", cycle)
  out <- data.frame(SEQN = s1$SEQN)
  rows <- match(out$SEQN, s2$SEQN)
  zero <- function(x) ifelse(is.na(x), 0, x)
  for (v in variables) out[[v]] <- (zero(s1[[paste0("DS1T", v)]]) + zero(s2[[paste0("DS2T", v)]][rows])) / 2
  out
}

# The residual method (caffeine regressed on energy by least squares in the analytic population,
# residual plus the mean caffeine) and the quartile cutpoints, all from the paper's own cycles.
row066_constants <- function(data) {
  p <- data$in_population %in% TRUE
  fit <- stats::lm(caffeine ~ energy, data = data[p, ])
  mean_caffeine <- mean(data$caffeine[p])
  adjusted <- stats::residuals(fit) + mean_caffeine
  list(intercept = unname(stats::coef(fit)[1]), slope = unname(stats::coef(fit)[2]), mean = mean_caffeine,
       adjusted_cuts = unname(stats::quantile(adjusted, c(0.25, 0.5, 0.75))),
       raw_cuts = unname(stats::quantile(data$caffeine[p], c(0.25, 0.5, 0.75))))
}

# Quartile 1 below the first cutpoint, quartile 4 at or above the third. Unadjusted caffeine ties at
# zero (a fifth of the sample) and its first cutpoint is 0.5 mg, so its quartile 1 is those with no
# caffeine, as Table 1's quartile 1 (mean 0.00 mg/day) is.
row066_quartile <- function(x, cuts) {
  labels <- c("Q1", "Q2", "Q3", "Q4")
  factor(labels[findInterval(x, cuts) + 1], levels = labels)
}

# Single regression imputation of a missing income-to-poverty ratio (the only covariate with
# missing values) within the study population, from every other Model 3 variable, the exposure, and
# the outcome; predictions are kept within the variable's range (0 to 5). It stands in for the
# paper's "8 multiple imputations using fully conditional specifications", which with one
# incomplete variable regress it on the others; one conditional-mean imputation gives about the
# same point estimate.
row066_impute_pir <- function(data) {
  predictors <- c("age", "sex", "race", ROW066_NUTRIENTS, "caffeine_adjusted", "overweight_obese")
  ok <- data$in_population %in% TRUE & stats::complete.cases(data[, predictors])
  missing <- ok & is.na(data$pir)
  if (!any(missing)) return(data$pir)
  fit <- stats::lm(stats::reformulate(predictors, "pir"), data = data[ok & !is.na(data$pir), ])
  pir <- data$pir
  pir[missing] <- pmin(pmax(stats::predict(fit, data[missing, ]), 0), 5)
  pir
}

# `within` keeps only the given quartiles of unadjusted caffeine, as Table 2's cells are
# per-milligram odds ratios within a quartile. `impute` fills a missing income-to-poverty ratio
# (a variant; the chosen version uses complete cases).
row066_derive <- function(data, constants, impute = FALSE, within = NULL) {
  if (is.null(constants)) constants <- row066_constants(data)
  data$caffeine_adjusted <- data$caffeine - (constants$intercept + constants$slope * data$energy) + constants$mean
  data$caffeine_q <- row066_quartile(data$caffeine_adjusted, constants$adjusted_cuts)
  data$caffeine_raw_q <- row066_quartile(data$caffeine, constants$raw_cuts)
  if (impute) data$pir <- row066_impute_pir(data)
  if (!is.null(within)) data$in_population <- data$in_population & data$caffeine_raw_q %in% within
  data
}

association <- list(
  id = "row066", row = 66, doi = "10.1371/journal.pone.0300566",
  # The paper's 16,225 participants aged 2-19 are exactly those of DEMO_G, DEMO_H, DEMO_I and
  # P_DEMO, so 2017-2018 entered only through the 2017-March 2020 files (no overlap).
  cycles = c("2011-2012", "2013-2014", "2015-2016", "2017-2020"),
  weight = "WTDR2D",
  # Table 2's "Quartile 4" cell is not a contrast of quartiles: each cell is the per-mg odds ratio of
  # unadjusted caffeine among participants in that quartile of unadjusted caffeine, unweighted (its
  # crude and minimally adjusted cells are reproduced to three decimals; variants). The published
  # estimate comes from that computation, so the replication follows it; the quartile contrast the
  # text describes is a variant.
  weighted = FALSE,
  family = "logistic", term = "caffeine",
  published = list(measure = "OR", estimate = 1.0020, low = 1.0005, high = 1.0036, n = 10001, events = NULL,
                   contrast = "per mg of caffeine a day among those in the highest quartile of caffeine intake (the cell Table 2 labels quartile 4 against quartile 1)"),
  left_out = character(),
  build = function(cycle) {
    demo <- demographics(cycle)
    # Two-day means for those with two reliable recalls (status 1 both days). Through March 2020 the
    # first recall was in person and the second by telephone; in 2021-2023 both were by telephone:
    # same method, foods database, and units, a change of mode listed, not corrected.
    diet <- dietary_totals(cycle, unname(ROW066_DIET), days = 2)
    d <- merge_all(demo, body_measures(cycle)[, c("SEQN", "bmi_child")], diet)
    for (v in names(ROW066_DIET)) d[[v]] <- d[[ROW066_DIET[[v]]]]
    # Dietary covariates with 24-hour supplement intakes added, for the variant that reads the
    # Table 2 note's "total intake (from dietary and supplements)" that way.
    if (cycle == REPLICATION_CYCLE) {
      for (v in ROW066_NUTRIENTS) d[[paste0(v, "_total")]] <- NA_real_
    } else {
      supplements <- row066_supplements(cycle, unname(ROW066_DIET[ROW066_NUTRIENTS]))
      rows <- match(d$SEQN, supplements$SEQN)
      for (v in ROW066_NUTRIENTS) {
        extra <- supplements[[ROW066_DIET[[v]]]][rows]
        d[[paste0(v, "_total")]] <- d[[v]] + ifelse(is.na(extra), 0, extra)
      }
    }
    # BMDBMIC: CDC BMI-for-age category (1 underweight, 2 normal, 3 overweight, 4 obese), ages 2-19.
    d$overweight_obese <- ifelse(d$bmi_child %in% 3:4, 1, ifelse(d$bmi_child %in% 1:2, 0, NA))
    d$obese <- ifelse(d$bmi_child %in% 4, 1, ifelse(d$bmi_child %in% 1:3, 0, NA))
    d$sex <- factor(d$sex, levels = 1:2, labels = c("male", "female"))
    d$race <- factor(d$race, levels = 1:5, labels = c("mexican american", "other hispanic", "white", "black", "other"))
    d$in_population <- d$age >= 2 & d$age <= 19 & !is.na(d$caffeine) & !is.na(d$energy) & !is.na(d$bmi_child)
    d
  },
  constants = function(data) row066_constants(data),
  derive = function(data, constants) row066_derive(data, constants, within = "Q4"),
  formula = ROW066_PER_MG,
  variants = list(
    list(label = "crude (published 1.0008, 0.9996-1.0019)", formula = overweight_obese ~ caffeine),
    list(label = "sex, age, race (published 1.0009, 0.9997-1.0022)", formula = overweight_obese ~ caffeine + sex + age + race),
    list(label = "weighted", weighted = TRUE),
    list(label = "within quartile 3, crude (published 1.0171, 1.0005-1.0340)", formula = overweight_obese ~ caffeine,
         derive = function(data, constants) row066_derive(data, constants, within = "Q3")),
    list(label = "within quartile 3, sex, age, race (published 1.0088, 0.9918-1.0262)", formula = overweight_obese ~ caffeine + sex + age + race,
         derive = function(data, constants) row066_derive(data, constants, within = "Q3")),
    list(label = "as the text describes: quartile 4 vs 1 of energy-adjusted caffeine, weighted", formula = ROW066_MODEL3, term = "caffeine_qQ4",
         weighted = TRUE, derive = function(data, constants) row066_derive(data, constants)),
    list(label = "as the text describes, crude", formula = overweight_obese ~ caffeine_q, term = "caffeine_qQ4", weighted = TRUE,
         derive = function(data, constants) row066_derive(data, constants)),
    list(label = "as the text describes, quartiles of unadjusted caffeine (Table 1's grouping)",
         formula = update(ROW066_MODEL3, . ~ . - caffeine_q + caffeine_raw_q), term = "caffeine_raw_qQ4", weighted = TRUE,
         derive = function(data, constants) row066_derive(data, constants)),
    list(label = "per mg of caffeine, all participants", derive = function(data, constants) row066_derive(data, constants)),
    list(label = "obese only (BMDBMIC 4) as the outcome", formula = update(ROW066_PER_MG, obese ~ .)),
    list(label = "Methods' Model 3 list (+K, Mg, vitamin D, -selenium)",
         formula = update(ROW066_PER_MG, . ~ . - selenium + potassium + magnesium + vitamin_d)),
    list(label = "income-to-poverty ratio imputed by regression, for the paper's multiple imputation", sample = "own",
         derive = function(data, constants) row066_derive(data, constants, impute = TRUE, within = "Q4"))
  ),
  # 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 = "reporting", affects_headline = TRUE, followed = TRUE,
         detail = "Table 2 and the Results give 1.0020 (1.0005 to 1.0036) as Model 3's odds ratio for quartile 4 against quartile 1, but Table 1's weighted obesity rates (33.80% in quartile 1, 40.77% in quartile 4) give a crude contrast near 1.35. Table 2's cells are per-milligram odds ratios of caffeine among those in each quartile: so computed, Model 3 gives 1.0015 (1.0000 to 1.0029) for quartile 4 here, and the crude cells for quartiles 3 and 4 give 1.0185 and 1.0011 (published 1.0171 and 1.0008)."),
    list(kind = "weighting", affects_headline = TRUE, followed = TRUE,
         detail = "The abstract says all data were survey-weighted and the Methods name combined dietary sample weights, but Table 2's cells and their interval widths are reproduced only unweighted (Model 3, quartile 4: 1.0015, 1.0000 to 1.0029 unweighted, and 1.0017, 0.9993 to 1.0041 weighted, against the published 1.0020, 1.0005 to 1.0036)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods adjust caffeine for energy by the residual method, and Tables 1 and 2 print quartile cutpoints of 236.25, 472.5, and 708.05, but Table 1's quartiles are of unadjusted caffeine, with mean intakes of 0.00, 1.79, 9.86, and 70.09 mg a day and mean energy (1,602 to 2,098 kcal) matching unadjusted quartiles here (1,594 to 2,077), not energy-adjusted ones (2,351, 1,647, 1,418, 1,884); Table 2's cells are reproduced within unadjusted quartiles."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's Overall column is the mean of its four quartile columns (caffeine 20.43, SD 18.86; 507 Mexican American participants, the mean of 434, 466, 598, and 530), not the whole sample's, and its 214 other-race participants in quartile 2 should be 428 for that quartile's counts to sum to its 2,497."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's weighted 'rate of obesity' (34.87%) is the share overweight or obese (BMI-for-age category 3 or 4: 35.4% here), not the share obese (19.3%)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Fig 1 starts from 45,932 participants and excludes 29,707 outside ages 2 to 19, but the four cycles' files hold 45,462, among them exactly the paper's 16,225 aged 2 to 19."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Results give the inflection point of caffeine as 479.15 where Table 3 prints 26.99, and Table 3's interval below that point (1.0009 to 1.0013) excludes its own estimate (1.0022)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract and Results give per-quartile effects that match no printed estimate (a 0.05% increase in Model 3 and a 0.03% decrease in Model 1, where Table 2's per-quarter odds ratios are 1.0022 and 1.0009), and the abstract gives Table 4's quartile 4 against quartile 1 estimates (1.5961 for boys, 1.4418 for girls) as per-quartile increments.")
  ),
  choices = list(
    list(choice = "outcome", decision = "overweight or obese (BMDBMIC 3 or 4) against underweight or normal weight (1 or 2)",
         reason = "the Methods compute ORs 'for overweight and obese individuals', and Table 1's weighted 'rate of obesity' of 34.87% fits BMDBMIC 3 or 4 (35.4% in this sample) and not 4 alone (19.3%)"),
    list(choice = "study population", decision = "aged 2-19 with two reliable 24-hour recalls and a BMI category (n 10,960), and of them, for the headline, those in the highest quartile of caffeine intake",
         reason = "the paper's flow (16,225 aged 2-19, less those without caffeine and without a BMI category); the 16,225 match, the 192 without a BMI category are close to the 198 here, but the 10,193 with caffeine (10,001 in the end) match no definition tried (two recalls 11,158, the first recall 13,144)"),
    list(choice = "caffeine", decision = "mean of the two recalls' totals (DR1TCAFF, DR2TCAFF)", reason = "the Methods: 'the average caffeine intake from two total nutrient intake recalls'"),
    list(choice = "energy adjustment", decision = "for the variants that follow the text's description: residual method, caffeine regressed on energy (two-day means) by least squares in the analytic population, the residual plus the mean caffeine; the paper's cycles fix the coefficients",
         reason = "the Methods say caffeine was adjusted for total energy 'with a residual model'; this is the usual (Willett) form"),
    list(choice = "quartiles", decision = "unweighted quartiles of unadjusted caffeine in the analytic population of the paper's cycles (quartile 1: no caffeine), reused in 2021-2023; the text's quartiles of energy-adjusted caffeine are used in its variants",
         reason = "Table 1's quartiles are of unadjusted caffeine (quartile 1 has no caffeine; unadjusted quartiles here have mean energy 1,594, 1,727, 1,878, 2,077 kcal against Table 1's 1,602, 1,741, 1,890, 2,098, adjusted ones 2,351, 1,647, 1,418, 1,884), and Table 2's cells are reproduced within them; the printed cutpoints (236.25, 472.5, 708.05: steps of 236.25) have no units and fit no caffeine scale"),
    list(choice = "weights", decision = "unweighted, as the paper's computation was; the two-day dietary weight (WTDR2D, WTDR2DPP for 2017-March 2020, pooled 2/9.2 and 3.2/9.2) in a variant and in the 'as the text describes' variants",
         reason = "the paper names 'combined dietary sample weights', but Table 2's cells and their interval widths are reproduced only by unweighted models"),
    list(choice = "Model 3 covariates", decision = "the Table 2 note's list: sex, age (continuous), race, income-to-poverty ratio (continuous), and two-day mean energy, protein, total sugar, total fat, cholesterol, vitamins B6 and C, calcium, phosphorus, zinc, copper, sodium, selenium",
         reason = "the note of the headline's table; the Methods' list (with potassium, magnesium and vitamin D, without selenium) is a variant"),
    list(choice = "'total intake (from dietary and supplements)'", decision = "not a covariate of its own; the dietary covariates come from the recalls' food and beverage totals",
         reason = "it names no variable, and the Methods say the dietary factors are 'the average of two total nutrient intakes' (DR1TOT, DR2TOT); adding the 24-hour supplement intakes to each dietary covariate is a variant (2021-2023 has no 24-hour supplement recall)"),
    list(choice = "energy adjustment of the dietary covariates", decision = "raw two-day means, with energy in the model",
         reason = "the notes say dietary confounders were also residual-adjusted; with energy in the model a linear residual is a reparametrization that leaves the caffeine estimate unchanged"),
    list(choice = "race", decision = "RIDRETH1's five groups as indicators", reason = "Table 1's groups; 'characteristics with three or more categories as indicator variables'"),
    list(choice = "missing income-to-poverty ratio", decision = "complete cases",
         reason = "the paper reports 8 multiple imputations by fully conditional specification, which kept those missing the ratio (Table 1's 2,028 Mexican American participants exceed the 1,889 with a ratio in this larger sample); multiple imputation's random draws can't be reproduced, so, as for every paper that used it, the re-implementation uses complete cases, and a single regression imputation of the ratio is a variant"),
    list(choice = "what Table 2 reports", decision = "the per-milligram odds ratio of unadjusted caffeine among those in its highest quartile, unweighted, as the paper computed the cell it labels quartile 4 against quartile 1; the quartile contrast the Methods and Results describe is a variant",
         reason = "Table 1's weighted prevalences (33.80% in Q1, 40.77% in Q4) give a crude Q4 vs Q1 OR near 1.35, not Table 2's 1.0008. Table 2's cells match per-milligram ORs of unadjusted caffeine within quartiles of unadjusted caffeine from unweighted models (crude Q3 1.0181, 1.0017-1.0348 here against 1.0171, 1.0005-1.0340; crude Q4 1.0010 against 1.0008; Model 2 Q3 1.0080 against 1.0088; Model 3 Q4 1.0015, 1.0001-1.0029 against 1.0020, 1.0005-1.0036), with interval widths that scale with each quartile's spread of caffeine; weighted fits and energy-adjusted caffeine match less well")
  )
)