# Composite Dietary Antioxidant Index (CDAI) and hyperlipidemia (NCEP ATP III), adults, NHANES
# 2005-2020. Zhao, Zhang, Zhang, Lin, and Cao (2024), Sci Rep, doi:10.1038/s41598-024-66922-0.
# Headline: OR 0.98 (95% CI 0.97-0.99) per unit of CDAI, Model 3 (Table 2).

# The intakes the CDAI sums: vitamin A (RAE), vitamin C, vitamin E (alpha-tocopherol), zinc,
# selenium, and five carotenoids, as the library's cdai_components() reads them, and energy.
ROW256_NUTRIENTS <- c("VARA", "VC", "ATOC", "ZINC", "SELE", "ACAR", "BCAR", "CRYP", "LYCO", "LZ", "KCAL")

# Total physical activity in MET-minutes a week in 2005-2006, which asked no GPAQ: leisure-time
# activities from the individual activity file (PAQIAF: each activity's MET score times minutes
# each time times times in the past 30 days), walking or cycling to get places (PAD020; times a
# day, week, or month times minutes a day, at 4 METs, the GPAQ's value for travel), and tasks in or
# around the home or yard needing moderate or greater effort (PAQ100; times in the past 30 days
# times minutes each time, at 4 METs, the GPAQ's moderate value), all per week. Missing where an
# answer needed is.
row256_met_2005 <- function(cycle) {
  d <- component("PAQ", cycle)
  iaf <- component("PAQIAF", cycle)
  valid <- function(x, top) ifelse(x >= 0 & x <= top, x, NA)
  per_week <- 7 / 30.4375
  leisure_rows <- valid(iaf$PADMETS, 20) * valid(iaf$PADDURAT, 1440) * valid(iaf$PADTIMES, 900) * per_week
  leisure <- tapply(leisure_rows, iaf$SEQN, sum)
  any_leisure <- d$PAD200 %in% 1 | d$PAD320 %in% 1
  no_leisure <- d$PAD200 %in% 2:3 & d$PAD320 %in% 2:3
  leisure_met <- ifelse(no_leisure, 0, ifelse(any_leisure, unname(leisure[as.character(d$SEQN)]), NA))
  times <- valid(d$PAQ050Q, 1000) * ifelse(d$PAQ050U %in% 1, 7, ifelse(d$PAQ050U %in% 2, 1, ifelse(d$PAQ050U %in% 3, per_week, NA)))
  travel_met <- ifelse(d$PAD020 %in% 2:3, 0, ifelse(d$PAD020 %in% 1, 4 * times * valid(d$PAD080, 1440), NA))
  home_met <- ifelse(d$PAQ100 %in% 2:3, 0, ifelse(d$PAQ100 %in% 1, 4 * valid(d$PAD120, 900) * valid(d$PAD160, 1440) * per_week, NA))
  data.frame(SEQN = d$SEQN, met = leisure_met + travel_met + home_met)
}

# Taking prescribed medicine to lower cholesterol: 1 yes, 0 no, missing where it can't be told.
# 2005-2020 ask it (BPQ100D) only of those told their cholesterol was high (BPQ080) and told to
# take medicine for it (BPQ090D), and 2005-2006 ask BPQ080 only of those whose cholesterol was ever
# checked; 2021-2023 ask everyone (BPQ101D).
row256_lipid_medicine <- function(cycle) {
  d <- component("BPQ", cycle)
  get <- function(v) if (has(d, v)) d[[v]] else rep(NA, nrow(d))
  taking <- if (has(d, "BPQ101D")) d$BPQ101D else get("BPQ100D")
  data.frame(SEQN = d$SEQN,
             lipid_medicine = ifelse(taking %in% 1, 1, ifelse(get("BPQ080") %in% 2 | get("BPQ090D") %in% 2 | taking %in% 2, 0, NA)))
}

# Drinking as the paper defines it: heavy is 3 or more drinks a day for women or 4 or more for
# men, or binge drinking on 5 or more days a month; moderate is 2 or more a day for women or 3 or
# more for men, or binge drinking on 2 or more days a month; mild is every other answer, never and
# former drinkers included (Table 1 has no other group). The library's alcohol5() also counts as
# heavy anyone who ever drank 4 or 5 drinks every day (ALQ151), which the paper does not.
row256_alcohol <- function(alq, sex) {
  binge_month <- alq$binge_days / 12
  heavy <- (sex == 2 & alq$drinks_per_day >= 3) | (sex == 1 & alq$drinks_per_day >= 4) | binge_month >= 5
  moderate <- (sex == 2 & alq$drinks_per_day >= 2) | (sex == 1 & alq$drinks_per_day >= 3) | binge_month >= 2
  group <- ifelse(alq$alcohol3 %in% 1:2, "mild", ifelse(alq$alcohol3 %in% 3 & heavy %in% TRUE, "heavy",
                  ifelse(alq$alcohol3 %in% 3 & moderate %in% TRUE, "moderate", ifelse(alq$alcohol3 %in% 3, "mild", NA))))
  factor(group, levels = c("mild", "moderate", "heavy"))
}

# One cycle's frame. In 2021-2023 the total physical activity of every domain can't be built.
row256_frame <- function(cycle) {
  demo <- demographics(cycle)
  diet <- dietary_totals(cycle, ROW256_NUTRIENTS, days = "mean_or_one")
  diet <- cbind(diet[, c("SEQN", "recall_day1", "recall_day2", "KCAL")], cdai_components(diet))
  alq <- alcohol(cycle)
  met <- if (cycle == REPLICATION_CYCLE) data.frame(SEQN = demo$SEQN, met = NA_real_) else if (cycle == "2005-2006") row256_met_2005(cycle) else met_minutes(cycle)
  d <- merge_all(demo, diet, component("DR1TOT", cycle, "DRQSDIET"), body_measures(cycle)[, c("SEQN", "bmi")],
                 total_cholesterol(cycle), hdl_cholesterol(cycle), triglycerides(cycle), row256_lipid_medicine(cycle),
                 smoking(cycle), diabetes_status(cycle, parts = c("told", "medication")), hypertension_status(cycle), met)
  d$alcohol <- row256_alcohol(alq[match(d$SEQN, alq$SEQN), ], d$sex)
  # Hyperlipidemia as the paper states it: any of TG 150 mg/dL or more, TC 200 or more, LDL 130 or
  # more, HDL 40 or less (men) or 50 or less (women), or cholesterol medicine. TG and LDL exist
  # only for the morning fasting subsample; elsewhere the other criteria decide.
  low_hdl <- (d$sex == 1 & d$hdl <= 40) | (d$sex == 2 & d$hdl <= 50)
  flags <- cbind(d$tg >= 150, d$tc >= 200, d$ldl >= 130, low_hdl, d$lipid_medicine == 1)
  d$hyperlipidemia <- ifelse(rowSums(flags, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(flags)) > 0, 0, NA))
  # The same, but undetermined where a criterion not measured could have been met.
  d$hyperlipidemia_measured <- ifelse(rowSums(flags, na.rm = TRUE) > 0, 1, ifelse(rowSums(is.na(flags)) == 0, 0, NA))
  d$sex <- factor(d$sex, levels = 1:2)
  d$race <- factor(d$race, levels = 1:5)
  d$education <- factor(education3(d$education), levels = 1:3)
  d$pir3 <- factor(ifelse(d$pir < 1.5, "low", ifelse(d$pir <= 3.5, "middle", ifelse(d$pir > 3.5, "high", NA))), levels = c("low", "middle", "high"))
  # Married (living with a partner too, as DMDMARTZ groups them from 2017) against single.
  d$married <- factor(ifelse(d$marital %in% 1, 1, ifelse(d$marital %in% 2:3, 0, NA)), levels = 0:1)
  d$smoking <- factor(d$smoking, levels = 1:3)
  d$diabetes <- factor(d$diabetes, levels = 0:1)
  d$hypertension <- factor(d$hypertension, levels = 0:1)
  d$energy <- d$KCAL
  d$stack_copy <- FALSE
  # Adults (18 or older) with a reliable first-day recall, TC and HDL measured and cholesterol
  # medicine known (the flow chart's "hyperlipidemia data"), not on a special diet (DRQSDIET), and
  # mean energy intake of 500 to under 5000 kcal a day.
  base <- d$age >= 18 & d$recall_day1 %in% 1 & !is.na(d$tc) & !is.na(d$hdl) & !(d$DRQSDIET %in% 1) & (d$KCAL >= 500 & d$KCAL < 5000) %in% TRUE
  d$in_population <- base & !is.na(d$lipid_medicine)
  d$in_population_any_medicine <- base
  d
}

# The cycle's frame; with stack = TRUE, the 2017-March 2020 frame also carries a second copy of
# every 2017-2018 participant, from the 2017-2018 files, as the paper's stacked files did. The
# 2017-March 2020 files number their participants anew (SEQN 109263 on), so the copies keep their
# 2017-2018 SEQNs without colliding, and carry their 2017-2018 day-one dietary weight, which no
# 2017-March 2020 weight file holds.
row256_build <- function(cycle, stack = FALSE) {
  d <- row256_frame(cycle)
  d$copy_weight <- NA_real_
  if (stack && cycle == "2017-2020") {
    copy <- row256_frame("2017-2018")
    copy$stack_copy <- TRUE
    w <- weights_for("WTDRD1", "2017-2018")
    copy$copy_weight <- w$weight_2yr[match(copy$SEQN, w$SEQN)]
    d <- rbind(d, copy[, names(d)])
  }
  d
}

# Each CDAI component's mean and SD over the paper's sample (unweighted, as the paper reports no
# weights), used unchanged in 2021-2023.
row256_constants <- function(data) {
  s <- data[data$in_population %in% TRUE, ]
  list(cdai = lapply(setNames(CDAI_COMPONENTS, CDAI_COMPONENTS), function(k) list(mean = mean(s[[k]], na.rm = TRUE), sd = stats::sd(s[[k]], na.rm = TRUE))))
}

row256_derive <- function(data, constants) {
  if (is.null(constants)) constants <- row256_constants(data)
  data$cdai <- cdai(data, constants$cdai)
  data
}

# The stacked variant: the copies' weight is their 2017-2018 weight with a 2-year cycle's share of
# the pooled years (the variant is unweighted, as the chosen version is, so this only keeps them
# in the sample), and the CDAI is standardized over the stacked sample, as the paper's was.
row256_derive_stacked <- function(data, constants) {
  copy <- data$stack_copy %in% TRUE
  years <- sum(vapply(unique(data$cycle), function(cycle) cycle_info(cycle)$years, 0))
  data$w[copy] <- data$copy_weight[copy] * 2 / years
  row256_derive(data, row256_constants(data))
}

row256_model3 <- hyperlipidemia ~ cdai + age + sex + bmi + race + education + pir3 + married + alcohol + smoking + diabetes +
  hypertension + met + energy

association <- list(
  id = "row256", row = 256, doi = "10.1038/s41598-024-66922-0",
  # The flow chart starts from 85,750 participants, the sum of the 2005-2018 files and the
  # 2017-March 2020 files, which hold the 2017-2018 participants a second time. The 2017-March
  # 2020 files cover those years once; a variant stacks them as the paper did.
  cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2020"),
  weight = "WTDRD1", blood_file = "TCHOL",
  # The paper never mentions weights, and its Table 1 counts are unweighted.
  weighted = FALSE,
  family = "logistic", term = "cdai",
  published = list(measure = "OR", estimate = 0.98, low = 0.97, high = 0.99, n = 30788, events = 25525,
                   contrast = "per unit of the CDAI (sum of the six intakes, each minus its mean over its SD, means and SDs from the paper's sample)"),
  left_out = c("physical activity as total MET-minutes a week over work, travel, and leisure (2021-2023 asks only about leisure-time activity)"),
  build = function(cycle) row256_build(cycle),
  constants = row256_constants,
  derive = row256_derive,
  # Hyperlipidemia as the paper computed it: undetermined (and so out of the complete-case model,
  # where the paper imputed it) where a criterion not measured outside the fasting subsample could
  # have been met. Its 82.91% is reproduced so (83.4%); the stated definition gives 70.7%.
  formula = update(row256_model3, hyperlipidemia_measured ~ .),
  formula_harmonized = hyperlipidemia_measured ~ cdai + age + sex + bmi + race + education + pir3 + married + alcohol + smoking + diabetes +
    hypertension + energy,
  variants = list(
    list(label = "stacked 2017-2018 and 2017-March 2020 files, as the paper did", build = function(cycle) row256_build(cycle, stack = TRUE),
         derive = row256_derive_stacked),
    list(label = "Model 1, crude (published 0.97, 0.96-0.98)", formula = hyperlipidemia_measured ~ cdai),
    list(label = "Model 1, crude, own sample", formula = hyperlipidemia_measured ~ cdai, sample = "own"),
    list(label = "Model 2, age and sex (published 0.98, 0.97-0.98)", formula = hyperlipidemia_measured ~ cdai + age + sex),
    list(label = "weighted (dietary day-one weight)", weighted = TRUE),
    list(label = "hyperlipidemia as stated, criteria not measured counted as not met", sample = "own", formula = row256_model3),
    list(label = "stacked, and hyperlipidemia as stated", sample = "own", formula = row256_model3,
         build = function(cycle) row256_build(cycle, stack = TRUE), derive = row256_derive_stacked),
    list(label = "unknown cholesterol medicine as none (no medicine step)", sample = "own",
         derive = function(data, constants) { data <- row256_derive(data, constants); data$in_population <- data$in_population_any_medicine; data }),
    list(label = "two reliable recalls required", sample = "own",
         derive = function(data, constants) { data <- row256_derive(data, constants); data$in_population <- data$in_population & data$recall_day2 %in% 1; data })
  ),
  # 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,
         detail = "The paper analyzes NHANES 2005-2020, but the 85,750 participants of its flow chart (Fig. 1; 33,914 under 18, 51,836 adults) are exactly the 2005-2018 files together with the 2017-March 2020 files, which hold the 2017-2018 participants a second time. Stacked that way the population here is 30,755 (published 30,788), 3,868 of them second copies, and Model 3 gives 0.9876 (0.9782-0.9971), against 0.9883 (0.9783-0.9984) with each participant once."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE,
         detail = "The text gathers participants who had two dietary recalls and averages the two days, but its 6,261 excluded for missing dietary data are exactly those without a reliable first-day recall (requiring both recalls would exclude 12,119), so its population includes those with one; requiring two gives 19,895 people and Model 3 0.9866 (0.9756-0.9977)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The text defines hyperlipidemia as any of five criteria, two of which (TG and LDL) NHANES measures only in its fasting subsample, and its population keeps those outside it (its 8,497 without hyperlipidemia data are matched by those without TC, HDL, or medicine status, 8,492 here). That definition gives 70.7% with hyperlipidemia here, against the paper's 82.91% (25,525 of 30,788); leaving undetermined those whom an unmeasured criterion could decide gives 83.4%, so the outcome was coded that way, with the undetermined filled in by the paper's imputation."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1 gives physical activity as a median of 2,160 MET-minutes a week (720 to 6,000), but the stated scoring gives 1,100 (0 to 3,840) over the population here, a quarter of whom report none, and 2,100 (724 to 5,700) among those with any activity.")
  ),
  choices = list(
    list(choice = "cycles", decision = "2005-2006 through 2015-2016 and the 2017-March 2020 files, each participant once",
         reason = "the flow chart's 85,750, 33,914 under 18, and 51,836 adults are exactly the 2005-2018 files stacked with the 2017-March 2020 files, which repeat the 2017-2018 participants; a variant stacks them"),
    list(choice = "missing dietary data", decision = "no reliable first-day recall (DR1DRSTZ not 1)", reason = "gives the flow chart's 6261 exactly; requiring both days would exclude 12,119"),
    list(choice = "daily intakes", decision = "mean of the two recalls, the first day's where the second is missing or unreliable, for the CDAI and energy",
         reason = "the paper averages two days yet keeps everyone with a first-day recall; this gives Table 1's CDAI median and quartiles (-0.63, -2.68 to 1.94; published -0.61, -2.68 to 1.94) and its energy quartiles closely; a variant requires two recalls"),
    list(choice = "carotenoids and vitamins", decision = "carotenoids as alpha- plus beta-carotene, beta-cryptoxanthin, lycopene, and lutein with zeaxanthin; vitamin A as RAE; vitamin E as alpha-tocopherol",
         reason = "unstated beyond 'total carotenoids'; the library's CDAI components, which reproduce Table 1's CDAI quartiles"),
    list(choice = "standardizing means and SDs", decision = "unweighted, over the paper's sample (in_population)", reason = "the paper standardizes over its participants and reports no weights"),
    list(choice = "missing hyperlipidemia data", decision = "TC or HDL unmeasured, or cholesterol medicine unknown (BPQ080 unanswered, or told of high cholesterol with the medicine questions unanswered)",
         reason = "gives 8492 against the flow chart's 8497; TC and HDL alone give 2505. A variant drops the medicine step"),
    list(choice = "special diet and energy", decision = "DRQSDIET = 1; mean energy under 500 or 5000 kcal or more", reason = "gives 6061 and 351 against the flow chart's 5914 and 376"),
    list(choice = "hyperlipidemia", decision = "any criterion met (HDL 40 or less for men, 50 or less for women, as stated; medicine BPQ100D, BPQ101D in 2021-2023), and undetermined where a criterion not measured (TG and LDL outside the fasting subsample) could have been met, as the paper's computation coded it; the complete-case model leaves the undetermined out, where the paper imputed them",
         reason = "the stated definition, with criteria not measured counted as not met, gives 70.7% against the paper's 82.91% (25,525 of 30,788); leaving undetermined those whom an unmeasured criterion could decide gives 83.4%, so the published estimate comes from that coding; the stated definition is a variant"),
    list(choice = "weights", decision = "unweighted", reason = "the paper never mentions weights or design and its Table 1 counts are unweighted; a variant uses the dietary day-one weight"),
    list(choice = "missing covariates", decision = "complete cases", reason = "the paper used unspecified multiple imputation; complete cases also drop every 18-19 year old, whom NHANES asks no adult education or marital status"),
    list(choice = "marital status", decision = "married or living with a partner against single", reason = "2017-March 2020 groups living with a partner with married; that grouping gives Table 1's 59.86% married (60.5%), leaving partners out gives 54.8%"),
    list(choice = "alcohol", decision = "the paper's three groups, never and former drinkers as mild, without the library's daily binge (ALQ151) item",
         reason = "the stated definition; it gives 65.9, 15.4, 18.7% against Table 1's 63.6, 16.3, 20.1 (the library's alcohol5 gives 62.7, 13.6, 23.6)"),
    list(choice = "hypertension", decision = "mean blood pressure 140/90 or more, told, or taking medicine (library hypertension_status)", reason = "the stated definition; 'exceeding' read as at or above; gives 42.3% against Table 1's 41.85%"),
    list(choice = "diabetes", decision = "told, insulin, or diabetes pills", reason = "the stated definition, with no glucose or HbA1c criterion"),
    list(choice = "physical activity", decision = "total GPAQ MET-minutes a week (library met_minutes) in 2007-2020; in 2005-2006 leisure activities (PAQIAF), travel, and home or yard tasks",
         reason = "the stated definition; 2005-2006 asked no GPAQ, and leaving it missing would drop the cycle from a complete-case analysis. Table 1's median (2160) is about twice the GPAQ's (1100), which no stated scoring explains; leaving activity out changes the estimate by under 0.001"),
    list(choice = "covariate forms", decision = "age, BMI, MET-minutes, and energy continuous; income-to-poverty ratio in the paper's three groups (under 1.5, 1.5-3.5, over 3.5); education in three", reason = "Table 1 reports age, BMI, activity, and energy as means or medians and the others as groups")
  )
)