# Composite Dietary Antioxidant Index (CDAI) and atherosclerotic cardiovascular disease (ASCVD:
# coronary heart disease, angina, heart attack, or stroke) in postmenopausal women, NHANES
# 2013-2018. Liu, Lai, Zhao, Zhang, and Hu (2023), Antioxidants, doi:10.3390/antiox12091740.
# Headline: OR 0.67 (95% CI 0.51-0.88) per SD of CDAI, Model C (Table 2).
# The population is reproduced exactly (3109 women, 453 cases, Table 1's age and race counts), but
# not the strength of the association: the crude OR per SD is 0.78 here against the published
# 0.65, under every CDAI construction tried, while Table 1's CDAI medians by ASCVD status (-1.14
# and -0.05) fit this data (-1.04 and -0.09).

# 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.
ROW319_NUTRIENTS <- c("VARA", "VC", "ATOC", "ZINC", "SELE", "ACAR", "BCAR", "CRYP", "LYCO", "LZ")

# "Alcohol use", at least 12 drinks of any alcoholic beverage in any one year. 2013-2016 ask it
# (ALQ101); 2017-2018 does not, and there it is read as at least 12 drinks in the past year
# (drinking days a year from ALQ121 times drinks a drinking day, ALQ130), never and former drinkers
# none (alcohol_12). `alcohol_monthly` is the other reading of 2017-2018: drinking at least once a
# month in the past year (ALQ121 codes 1-7), never drinkers missing.
row319_alcohol <- function(cycle) {
  d <- component("ALQ", cycle)
  if (has(d, "ALQ101")) return(data.frame(SEQN = d$SEQN, alcohol_12 = yes(d$ALQ101), alcohol_monthly = yes(d$ALQ101)))
  a <- alcohol(cycle)
  drinks <- ifelse(a$alcohol3 %in% 1:2, 0, a$drinking_days * a$drinks_per_day)
  data.frame(SEQN = a$SEQN, alcohol_12 = ifelse(drinks >= 12, 1, ifelse(drinks < 12, 0, NA)),
             alcohol_monthly = ifelse(d$ALQ121 %in% 1:7, 1, ifelse(d$ALQ121 %in% c(0, 8:10), 0, NA))[match(a$SEQN, d$SEQN)])
}

# Any moderate or vigorous recreational activity: 2013-2018 ask whether one does any in a typical
# week (PAQ650 vigorous, PAQ665 moderate); 2021-2023 ask how often (PAD790Q moderate, PAD810Q
# vigorous, 0 for never).
row319_recreation <- function(cycle) {
  d <- component("PAQ", cycle)
  if (has(d, "PAQ650")) {
    any <- ifelse(d$PAQ650 %in% 1 | d$PAQ665 %in% 1, 1, ifelse(d$PAQ650 %in% 2 & d$PAQ665 %in% 2, 0, NA))
  } else {
    often <- function(q) ifelse(q %in% 0, 0, ifelse(!is.na(q) & q > 0 & q < 7777, 1, NA))
    moderate <- often(d$PAD790Q)
    vigorous <- often(d$PAD810Q)
    any <- ifelse(moderate %in% 1 | vigorous %in% 1, 1, ifelse(moderate %in% 0 & vigorous %in% 0, 0, NA))
  }
  data.frame(SEQN = d$SEQN, recreation = any)
}

# The CDAI's standardizing constants: each component's mean and SD over `rows`, and the SD of the
# CDAI they give there (unweighted, and weighted by w for a variant).
row319_constants <- function(data, rows = data$in_population %in% TRUE) {
  s <- data[rows, ]
  stats <- lapply(setNames(CDAI_COMPONENTS, CDAI_COMPONENTS), function(k) list(mean = mean(s[[k]], na.rm = TRUE), sd = stats::sd(s[[k]], na.rm = TRUE)))
  score <- cdai(s, stats)
  weighted_mean <- sum(s$w * score) / sum(s$w)
  list(cdai = stats, cdai_sd = stats::sd(score), cdai_sd_weighted = sqrt(sum(s$w * (score - weighted_mean)^2) / sum(s$w)))
}

row319_derive <- function(data, constants, sd = "cdai_sd") {
  if (is.null(constants)) constants <- row319_constants(data)
  data$cdai <- cdai(data, constants$cdai)
  data$cdai_sd <- data$cdai / constants[[sd]]
  data
}

row319_build <- function(cycle) {
  replication <- cycle == REPLICATION_CYCLE
  demo <- demographics(cycle)
  diet <- dietary_totals(cycle, ROW319_NUTRIENTS, days = "mean_or_one")
  diet <- cbind(diet[, c("SEQN", "recall_day1")], cdai_components(diet))
  day1 <- dietary_totals(cycle, c("KCAL", "PFAT"), days = 1)[, c("SEQN", "KCAL", "PFAT")]
  mcq <- component("MCQ", cycle)
  d <- merge_all(demo, component("RHQ", cycle, "RHD043"), mcq[, c("SEQN", "MCQ160C", "MCQ160D", "MCQ160E", "MCQ160F")],
                 diet, day1, body_measures(cycle)[, c("SEQN", "bmi", "waist")], blood_count(cycle)[, c("SEQN", "neutrophils", "lymphocytes")],
                 hdl_cholesterol(cycle), total_cholesterol(cycle), component("SMQ", cycle, "SMQ020"), component("BPQ", cycle, "BPQ020"),
                 diabetes_status(cycle, parts = c("told", "medication", "hba1c")), row319_recreation(cycle))
  # 2021-2023 asks none of these: six-group marital status (DMDMARTZ has three), trouble sleeping
  # told to a doctor (SLQ050), a relative's heart attack or angina before 50 (MCQ300A), and
  # 12 drinks in any one year.
  if (replication) {
    d$DMDMARTL <- NA
    d$SLQ050 <- NA
    d$MCQ300A <- NA
    d$alcohol_12 <- NA
    d$alcohol_monthly <- NA
  } else {
    d <- merge_all(d, component("DEMO", cycle, "DMDMARTL"), component("SLQ", cycle, "SLQ050"), mcq[, c("SEQN", "MCQ300A")], row319_alcohol(cycle))
  }
  items <- d[, c("MCQ160C", "MCQ160D", "MCQ160E", "MCQ160F")]
  d$ascvd <- as.integer(rowSums(items == 1, na.rm = TRUE) > 0)
  answered <- rowSums(items == 1 | items == 2, na.rm = TRUE) == 4
  # Table S1's groups, Table 1's energy groups, and HDL-C as Table 1 shows it (0 or 1: 50 mg/dL or more).
  d$age5 <- cut(d$age, c(-Inf, 50, 60, 70, 80, Inf), right = FALSE, labels = c("<50", "50-59", "60-69", "70-79", "80"))
  d$race <- factor(d$race, levels = 1:5)
  d$education <- factor(education3(d$education), levels = 1:3)
  d$marital6 <- factor(ifelse(d$DMDMARTL %in% 1:6, d$DMDMARTL, NA), levels = 1:6)
  d$pir3 <- cut(d$pir, c(-Inf, 1, 3, Inf), labels = c("<=1.00", "1.01-3.00", ">3.00"))
  d$bmi4 <- cut(d$bmi, c(-Inf, 18.5, 25, 30, Inf), right = FALSE, labels = c("<18.5", "18.5-24.9", "25.0-29.9", "30+"))
  d$alcohol <- factor(d$alcohol_12, levels = 0:1)
  d$alcohol_m <- factor(d$alcohol_monthly, levels = 0:1)
  d$smoked_100 <- factor(yes(d$SMQ020), levels = 0:1)
  d$recreation <- factor(d$recreation, levels = 0:1)
  d$sleep_disorder <- factor(yes(d$SLQ050), levels = 0:1)
  d$hypertension <- factor(yes(d$BPQ020), levels = 0:1)
  d$diabetes <- factor(d$diabetes, levels = 0:1)
  d$family_mi <- factor(yes(d$MCQ300A), levels = 0:1)
  d$nlr <- d$neutrophils / d$lymphocytes
  d$hdl50 <- factor(ifelse(d$hdl >= 50, 1, ifelse(d$hdl < 50, 0, NA)), levels = 0:1)
  d$energy4 <- cut(d$KCAL, c(-Inf, 1550, 1973, 2555, Inf), right = FALSE, labels = c("<1550", "1550-1972", "1973-2554", "2555+"))
  d$pufa <- d$PFAT
  # Postmenopausal: "Menopause/Change of life" as the reason for no periods (RHD043 = 7), or any
  # woman 55 or older (4057); a reliable first-day dietary recall (3513); all four ASCVD questions
  # answered yes or no (3479); and an income-to-poverty ratio (3109, with 453 cases).
  d$in_population <- d$sex == 2 & (d$RHD043 %in% 7 | d$age >= 55) & d$recall_day1 %in% 1 & answered & !is.na(d$pir)
  d
}

row319_model_b <- ascvd ~ cdai_sd + age5 + race + education + marital6 + pir3 + bmi4

association <- list(
  id = "row319", row = 319, doi = "10.3390/antiox12091740",
  cycles = c("2013-2014", "2015-2016", "2017-2018"),
  # Table 1's weighted sample size (44,737,249, with 5,814,546 for ASCVD) is exactly the sum of
  # WTMEC2YR / 3 over the 3109 women.
  weight = "WTMEC2YR", blood_file = "CBC",
  family = "logistic", term = "cdai_sd",
  published = list(measure = "OR", estimate = 0.67, low = 0.51, high = 0.88, n = 3109, events = 453,
                   contrast = "per SD 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("marital status in six groups (2021-2023 releases only DMDMARTZ's three groups)",
               "sleep disorders (2021-2023 does not ask SLQ050, whether one told a doctor about trouble sleeping)",
               "family history of heart attack or angina before 50 (2021-2023 does not ask MCQ300A)",
               "alcohol use as at least 12 drinks in any one year (asked as ALQ101 only through 2016; 2021-2023 does not ask it)"),
  build = row319_build,
  constants = function(data) row319_constants(data),
  derive = function(data, constants) row319_derive(data, constants),
  formula = ascvd ~ cdai_sd + age5 + race + education + marital6 + pir3 + bmi4 + waist + alcohol + smoked_100 + recreation +
    sleep_disorder + hypertension + diabetes + family_mi + nlr + hdl50 + tc + energy4 + pufa,
  formula_harmonized = ascvd ~ cdai_sd + age5 + race + education + pir3 + bmi4 + waist + smoked_100 + recreation +
    hypertension + diabetes + nlr + hdl50 + tc + energy4 + pufa,
  variants = list(
    list(label = "Model A, crude (published 0.65, 0.53-0.80)", formula = ascvd ~ cdai_sd),
    list(label = "Model A, crude, own sample", formula = ascvd ~ cdai_sd, sample = "own"),
    list(label = "Model B (published 0.71, 0.58-0.87)", formula = row319_model_b),
    list(label = "Model B, own sample", formula = row319_model_b, sample = "own"),
    list(label = "unweighted", weighted = FALSE),
    list(label = "age and BMI continuous", formula = ascvd ~ cdai_sd + age + race + education + marital6 + pir3 + bmi + waist + alcohol +
           smoked_100 + recreation + sleep_disorder + hypertension + diabetes + family_mi + nlr + hdl50 + tc + energy4 + pufa),
    list(label = "HDL-C continuous", formula = ascvd ~ cdai_sd + age5 + race + education + marital6 + pir3 + bmi4 + waist + alcohol +
           smoked_100 + recreation + sleep_disorder + hypertension + diabetes + family_mi + nlr + hdl + tc + energy4 + pufa),
    list(label = "2017-2018 alcohol: monthly drinking, never drinkers missing", sample = "own",
         formula = ascvd ~ cdai_sd + age5 + race + education + marital6 + pir3 + bmi4 + waist + alcohol_m + smoked_100 + recreation +
           sleep_disorder + hypertension + diabetes + family_mi + nlr + hdl50 + tc + energy4 + pufa),
    list(label = "per weighted SD of the CDAI", derive = function(data, constants) row319_derive(data, constants, sd = "cdai_sd_weighted"))
  ),
  # 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 = FALSE, evidence = "paper",
         detail = "On the paper's population, reproduced exactly, the CDAI as described gives Table 1's weighted medians (-1.04 with ASCVD and -0.09 without, against -1.14 and -0.05), but not Table 2: its quartile cutpoints (-1.04, 1.11, 3.72) lie above Table 1's quartiles (-2.27, -0.18, 2.46), and its crude odds ratio per SD (0.65, 0.53-0.80) is 0.78 (0.65-0.95) here under every CDAI construction tried."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The text adds women 55 or older with no menstrual information to those giving menopause as the reason for no periods, but its 4,057 postmenopausal women are reproduced only by adding every woman 55 or older whatever her menstrual answers (2,415 plus 1,642)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The text computes intakes as the average of the two recalls, but the 544 excluded for missing dietary data are exactly those without a reliable first-day recall (requiring both recalls would exclude 891), so women with only a first-day recall were kept; averaging the recalls each woman has fits Table 1's CDAI medians more closely than first-day intakes."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's counts for BMI 18.5-24.9 and 25.0-29.9 (1,297 and 356) do not fit their weighted percentages (25.4% and 29.7%), which the reproduced population gives from 715 and 923 women."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1 labels HDL-C in mg/dL but reports a 0/1 variable (median 1.00, interquartile range 0.00 to 1.00)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Results give Model B's per-SD interval as 0.58-0.81, where Table 2 has 0.71 (0.58, 0.87).")
  ),
  choices = list(
    list(choice = "postmenopausal", decision = "RHD043 = 7 (menopause/change of life), or any woman aged 55 or older whatever her menstrual answers",
         reason = "the only reading that gives the flow chart's 4057 exactly (2415 by RHD043 plus 1642 women 55 or older who did not answer menopause); 2021-2023 asks RHD043 by self-interview (ACASI) rather than face to face"),
    list(choice = "missing dietary data", decision = "no reliable first-day recall (DR1DRSTZ not 1)", reason = "gives the flow chart's 544 exactly; requiring both days would exclude 891"),
    list(choice = "missing ASCVD data", decision = "any of MCQ160C, D, E, F not answered yes or no", reason = "gives the flow chart's 34 exactly, and the final 3109 with 453 cases, whose age and race counts match Table 1 exactly"),
    list(choice = "daily intakes for the CDAI", decision = "mean of the two recalls, the first day's where the second is missing or unreliable",
         reason = "the paper averages two days yet keeps everyone with a first-day recall; this gives Table 1's weighted CDAI medians (-0.26 overall, -1.04 with ASCVD, -0.09 without; published -0.18, -1.14, -0.05) more closely than first-day intakes (-0.41, -1.08, -0.27)"),
    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; the library's CDAI components, and the closest fit to Table 1's CDAI medians"),
    list(choice = "standardizing means and SDs", decision = "unweighted, over the paper's 3109 women (in_population)", reason = "the paper standardizes over its subjects; unweighted means and SDs fit Table 1's medians, weighted ones shift them by about 0.25"),
    list(choice = "SD for the per-SD estimate", decision = "unweighted SD of the CDAI over the same 3109 women", reason = "unstated; a variant uses the weighted SD"),
    list(choice = "weights", decision = "WTMEC2YR / 3 with strata and PSUs", reason = "Table 1's weighted sample size equals the sum of WTMEC2YR / 3 over the analytic sample, to the person"),
    list(choice = "missing covariates", decision = "complete cases", reason = "the paper imputed covariates with under 10% missing by random forest; complete cases drop women missing waist, lipids, NLR, BMI, family history, alcohol, or another covariate"),
    list(choice = "age, BMI, income", decision = "Table S1's groups (age 40-49 to 80, four BMI groups, income-to-poverty ratio up to 1.00, 1.01-3.00, above 3.00)", reason = "Table S1 says each was categorized; a variant enters age and BMI continuous"),
    list(choice = "HDL-C", decision = "two levels, 50 mg/dL or more against less", reason = "Table 1 shows HDL-C as 0 or 1 (median 1.00, IQR 0.00-1.00) and the stratified analysis splits it at 50 mg/dL; a variant enters it in mg/dL"),
    list(choice = "energy and PUFA intake", decision = "first-day recall: energy in Table 1's four groups, PUFA in grams", reason = "Table 1's energy counts (1452, 719, 580, 358) and PUFA quartiles (15.02, 9.64-22.02) are the first day's exactly; the two-day mean gives other counts"),
    list(choice = "alcohol use in 2017-2018", decision = "at least 12 drinks in the past year (ALQ121 days times ALQ130 drinks), never and former drinkers none",
         reason = "2017-2018 does not ask ALQ101 and the paper does not say what it used; Table 1's 1457 users (56.0% weighted) fit several readings once its imputed values are allowed for; a variant uses monthly drinking with never drinkers missing"),
    list(choice = "moderate to vigorous recreational activity", decision = "yes to PAQ650 or PAQ665 (in 2021-2023 any moderate or vigorous leisure-time activity, PAD790Q or PAD810Q above 0)", reason = "Table S1; gives 1232 against Table 1's 1179, which no stated reading matches"),
    list(choice = "hypertension", decision = "told of high blood pressure (BPQ020)", reason = "the advice to take medicine (BPQ040A) is asked only of those told, so the stated definition reduces to BPQ020"),
    list(choice = "diabetes", decision = "told (DIQ010 = 1), insulin or pills, or HbA1c 6.5% or more, a missing part counting as not met", reason = "Table S1; gives Table 1's 833 exactly")
  )
)