# Dietary Inflammatory Index (DII), quartile 4 vs quartile 1, and self-reported stroke, adults
# 18-79, NHANES 1999-2018. Mao et al. (2024), BMC Public Health, doi:10.1186/s12889-023-17556-w.
# Headline: OR 1.87 (95% CI 1.53-2.29) for the highest against the lowest DII quartile, Model II,
# weighted (abstract and Table 3).

# Day-one recall intakes of every nutrient the DII draws on, and total folate, for those whose
# recall is reliable (status 1). dietary_totals() reads the status as DR1DRSTZ or DRDDRSTZ, but
# 1999-2000 names it DRDDRSTS, and it stops on a cycle whose file lacks a nutrient it is asked
# for (vitamin D before 2007; beta-carotene, vitamin A as RAE, alpha-tocopherol, and folic acid
# in 1999-2000), so a nutrient a file lacks is left missing here instead, and the DII leaves it
# out of the sum.
dii_intakes <- function(cycle) {
  d <- component("DR1TOT", cycle)
  prefix <- if (has(d, "DR1TKCAL")) "DR1T" else "DRXT"
  status <- d[[intersect(c("DR1DRSTZ", "DRDDRSTZ", "DRDDRSTS"), names(d))[1]]]
  out <- data.frame(SEQN = d$SEQN, recall_day1 = as.integer(status %in% 1))
  for (v in c(DII_NUTRIENTS, "FOLA")) {
    column <- paste0(prefix, v)
    out[[v]] <- if (has(d, column)) ifelse(status %in% 1, d[[column]], NA) else NA_real_
  }
  out
}

# Whether the study carries a component's file for a cycle (NHANES ran the OGTT only in 2005-2016).
has_file <- function(name, cycle) any(SOURCES$cycle == cycle & SOURCES$file == file_name(name, cycle))

# Drinking as a yes/no. The paper counts as drinkers "those consuming at least 12 drinks during the
# year preceding the survey". `lifetime`: had 12 or more drinks in any one year or in life (ALQ100
# in 1999-2000, ALD100 in 2001-2002, ALQ101 in 2003-2016, then ALQ110); 2017 on asks neither, so
# there it is ever having had a drink (ALQ111). `past_year`: 12 or more drinks in the past 12
# months, drinking days (ALQ120Q/U, or ALQ121's categories from 2017) times drinks on a drinking
# day (ALQ130). alcohol() reads ALQ101 or ALD100 but not ALQ100, so it stops on 1999-2000.
drinking12 <- function(cycle) {
  d <- component("ALQ", cycle)
  if (has(d, "ALQ111")) {
    lifetime <- yes(d$ALQ111)
    days <- c(`0` = 0, `1` = 365, `2` = 300, `3` = 182, `4` = 104, `5` = 52, `6` = 30, `7` = 12, `8` = 9, `9` = 4.5, `10` = 1.5)
    per_year <- unname(days[as.character(d$ALQ121)])
  } else {
    any_year <- yes(d[[intersect(c("ALQ101", "ALD100", "ALQ100"), names(d))[1]]])
    lifetime <- ifelse(any_year %in% 1 | yes(d$ALQ110) %in% 1, 1, ifelse(yes(d$ALQ110) %in% 0 | (any_year %in% 0 & is.na(d$ALQ110)), 0, NA))
    q <- ifelse(d$ALQ120Q %in% 0:365, d$ALQ120Q, NA)
    unit <- ifelse(d$ALQ120U %in% 1, 52, ifelse(d$ALQ120U %in% 2, 12, ifelse(d$ALQ120U %in% 3, 1, NA)))
    per_year <- ifelse(q %in% 0, 0, q * unit)
  }
  per_year[lifetime %in% 0] <- 0
  drinks <- ifelse(d$ALQ130 >= 1 & d$ALQ130 < 77, d$ALQ130, NA)  # 77/99 and 777/999: refused, don't know
  past_year <- ifelse(per_year %in% 0, 0, ifelse(per_year * drinks >= 12, 1, ifelse(!is.na(per_year * drinks), 0, NA)))
  data.frame(SEQN = d$SEQN, drinker_lifetime = lifetime, drinker_past_year = past_year)
}

# The paper's 26 DII parameters: the library's 28 without n-3 and n-6 fat.
DII_26 <- setdiff(DII_PARAMETERS$parameter, c("N3FAT", "N6FAT"))
QUARTILE_CUTS <- c(0.23, 1.76, 2.95)

association <- list(
  id = "row034", row = 34, doi = "10.1186/s12889-023-17556-w",
  cycles = c("1999-2000", "2001-2002", "2003-2004", "2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014",
             "2015-2016", "2017-2018"),
  weight = "WTDRD1",
  family = "logistic", term = "dii_qQ4",
  published = list(measure = "OR", estimate = 1.87, low = 1.53, high = 2.29, n = 44019, events = 1486,
                   contrast = "DII quartile 4 (2.95 or more) vs quartile 1 (under 0.23)"),
  left_out = c("the two-hour OGTT glucose part of the diabetes definition: 2021-2023 has no OGTT (NHANES ran it in 2005-2016 only)"),
  build = function(cycle) {
    demo <- demographics(cycle)
    intake <- dii_intakes(cycle)
    d <- merge_all(demo, intake, component("MCQ", cycle, "MCQ160F"), body_measures(cycle)[, c("SEQN", "bmi")],
                   component("SMQ", cycle, c("SMQ020", "SMQ040")), drinking12(cycle),
                   diabetes_status(cycle, parts = c("told", "hba1c", "glucose")),
                   bp_questions(cycle)[, c("SEQN", "told_hypertension", "bp_medication")],
                   blood_pressure(cycle, readings = 1:3)[, c("SEQN", "sbp", "dbp")])
    d$ogtt <- if (has_file("OGTT", cycle)) { o <- component("OGTT", cycle, "LBXGLT"); o$LBXGLT[match(d$SEQN, o$SEQN)] } else NA_real_
    dm <- component("DEMO", cycle)
    pregnant <- if (has(dm, "RIDEXPRG")) dm$RIDEXPRG[match(d$SEQN, dm$SEQN)] %in% 1 else FALSE

    # The DII over the paper's 26 parameters with Shivappa and colleagues' (2014) global means,
    # SDs, and effect scores, folic acid as DR1TFA (the library's). The parameter list names total
    # folate, but the DII's weighted medians in Table 2 and the Results (1.39 overall, 1.99 with
    # stroke, 1.37 without) fit folic acid (1.36, 1.99, 1.34) and not total folate (1.19, 1.81,
    # 1.16), which is a variant.
    folate <- intake
    folate$FA <- intake$FOLA
    d$dii <- dii(intake, DII_26)[match(d$SEQN, intake$SEQN)]
    d$dii_total_folate <- dii(folate, DII_26)[match(d$SEQN, intake$SEQN)]

    d$stroke <- ifelse(d$MCQ160F %in% 1, 1, ifelse(d$MCQ160F %in% 2, 0, NA))
    d$sex <- factor(d$sex, levels = c(1, 2), labels = c("male", "female"))
    d$race <- factor(d$race, levels = c(3, 4, 1, 2, 5), labels = c("white", "black", "mexican", "other hispanic", "other"))
    # Table 1's education shares (5.4%, 35.2%, 59.3%) are those of under 9th grade, 9th grade to
    # high school graduate, and more.
    d$education3 <- factor(ifelse(d$education %in% 1, 1, ifelse(d$education %in% 2:3, 2, ifelse(d$education %in% 4:5, 3, NA))),
                           levels = 1:3, labels = c("under 9th grade", "9th grade to high school", "over high school"))
    ever <- yes(d$SMQ020)
    d$smoker_ever <- factor(ever, levels = c(0, 1), labels = c("no", "yes"))
    d$smoker_current <- factor(ifelse(ever %in% 0, 0, ifelse(ever %in% 1 & d$SMQ040 %in% 1:2, 1, ifelse(ever %in% 1 & d$SMQ040 %in% 3, 0, NA))),
                               levels = c(0, 1), labels = c("no", "yes"))
    d$drinker <- factor(d$drinker_lifetime, levels = c(0, 1), labels = c("no", "yes"))
    d$drinker_past_year <- factor(d$drinker_past_year, levels = c(0, 1), labels = c("no", "yes"))
    # Diabetes: told by a doctor, HbA1c 6.5% or more, fasting glucose 126 mg/dL or more, or a
    # two-hour OGTT glucose of 200 mg/dL or more; a test a person lacks counts as not met.
    d$diabetes_harmonized <- d$diabetes
    d$diabetes <- ifelse(d$diabetes %in% 1 | (d$ogtt >= 200) %in% TRUE, 1, d$diabetes)
    # Hypertension: mean of the first three readings 140/90 mmHg or more, told, or taking medicine.
    flags <- cbind(d$told_hypertension == 1, d$bp_medication == 1, d$sbp >= 140, d$dbp >= 90)
    d$hypertension <- ifelse(rowSums(flags, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(flags)) > 0, 0, NA))
    d$in_population <- d$age >= 18 & d$age < 80 & !pregnant & d$recall_day1 %in% 1 & !is.na(d$stroke)
    d
  },
  derive = function(data, constants) {
    quartile <- function(x, cuts) cut(x, c(-Inf, cuts, Inf), right = FALSE, labels = c("Q1", "Q2", "Q3", "Q4"))
    data$dii_q <- quartile(data$dii, QUARTILE_CUTS)
    data$dii_total_folate_q <- quartile(data$dii_total_folate, QUARTILE_CUTS)
    if (!is.null(constants)) data$dii_sample_q <- quartile(data$dii, constants$sample_cuts)
    u <- function(x) with_unclear(x, levels(x))
    data$education3_u <- u(data$education3)
    data$smoker_current_u <- u(data$smoker_current)
    data$drinker_u <- u(data$drinker)
    data$diabetes_u <- with_unclear(data$diabetes, c(0, 1))
    data$hypertension_u <- with_unclear(data$hypertension, c(0, 1))
    data
  },
  # The sample's own quartiles of this DII, for the variant that splits it into quarters.
  constants = function(data) {
    keep <- data$in_population %in% TRUE
    list(sample_cuts = unname(stats::quantile(data$dii[keep], c(0.25, 0.5, 0.75))))
  },
  formula = stroke ~ dii_q + age + sex + race + education3 + smoker_current + drinker + bmi + diabetes + hypertension,
  formula_harmonized = stroke ~ dii_q + age + sex + race + education3 + smoker_current + drinker + bmi + diabetes_harmonized + hypertension,
  variants = list(
    list(label = "unweighted (published 1.89, 1.60-2.23)", weighted = FALSE),
    list(label = "crude (published 2.47, 2.08-2.94)", formula = stroke ~ dii_q),
    list(label = "Model I (published 2.43, 2.02-2.93)", formula = stroke ~ dii_q + age + sex + race),
    list(label = "crude, unweighted (published 2.37, 2.04-2.77)", formula = stroke ~ dii_q, weighted = FALSE),
    list(label = "Model I, unweighted (published 2.23, 1.91-2.60)", formula = stroke ~ dii_q + age + sex + race, weighted = FALSE),
    list(label = "MEC exam weight", weight = "WTMEC2YR"),
    list(label = "ever smoker and past-year drinker, as the text says", formula = stroke ~ dii_q + age + sex + race + education3 + smoker_ever +
           drinker_past_year + bmi + diabetes + hypertension, sample = "own"),
    list(label = "DII with total folate, as the parameter list says", formula = stroke ~ dii_total_folate_q + age + sex + race + education3 +
           smoker_current + drinker + bmi + diabetes + hypertension, term = "dii_total_folate_qQ4"),
    list(label = "missing covariates as their own level", formula = stroke ~ dii_q + age + sex + race + education3_u + smoker_current_u + drinker_u +
           bmi + diabetes_u + hypertension_u, sample = "own"),
    list(label = "sample quartiles of this DII", formula = stroke ~ dii_sample_q + age + sex + race + education3 + smoker_current + drinker + bmi +
           diabetes + hypertension, term = "dii_sample_qQ4")
  ),
  # 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,
         detail = "The text counts as smokers everyone who smoked 100 cigarettes, quitters included, but Table 1's smoking shares (22.3% overall, 29.6% with stroke, 22.1% without) are those of current smoking (22.6%, 30.3%, 22.4% weighted here), not of ever smoking (46.8%, 59.9%, 46.5%)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The text counts as drinkers those with at least 12 drinks in the year before the survey, but Table 1's drinking shares (89.31% without stroke, 84.3% with) are matched by 12 or more drinks in any one year or in life (89.4% and 83.3% here), not by 12 or more in the past year (59.5% and 33.2%)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The text and Table 1 group education as below high school, high school, and above high school, but Table 1's shares (5.4%, 35.2%, 59.3%) are those of under 9th grade, 9th grade to high school graduate, and more (5.3%, 35.3%, 59.4% here); below high school is 16.5%."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The DII's parameter list names total folate, but the DII values in Table 2 and the Results (1.39 overall, 1.99 with stroke, 1.37 without) fit a DII with folic acid (1.36, 1.99, 1.34 here) and not one with total folate (1.19, 1.81, 1.16)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The exclusion counts (46,469 aged under 18 or 80 and over, 1,516 pregnant, 5,747 without dietary data, 3,665 without stroke status) leave 43,919 of the 101,316, not the 44,019 analyzed; 46,369 are under 18 or 80 and over in the data, which leaves 44,019."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's laboratory rows don't fit their labels: the WBC, neutrophil, lymphocyte, platelet, and hemoglobin rows repeat the age, male sex, non-Hispanic White, non-Hispanic Black, and Mexican American rows, the TC row repeats the monocyte row, and the FBG, HbA1c, eGFR, TG, HDL-C, and RBC values don't fit their labels and units."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's note calls its categorical shares unweighted, but they are weighted (68.3% non-Hispanic White; its smoking shares match weighted ones here), and its overall drinking share (83.4%) lies outside both groups' (89.31% and 84.3%)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2 gives a mean vitamin E intake of 66.6 mg (66.8 without stroke, 56.0 with); the same recalls give 8.4 mg here."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 3 gives Model I's estimate for quartile 2 as 1.26 (1.01 to 1.58) with P 0.002, which its interval rules out (a lower bound of 1.01 puts P near 0.04).")
  ),
  choices = list(
    list(choice = "weights", decision = "day-one dietary weight (WTDRD1; the four-year WTDR4YR for 1999-2002), pooled over the ten cycles as NCHS directs",
         reason = "the paper says only that it used 'the sample weights corresponding to different research periods'; NCHS directs the day-one dietary weight for day-one intakes; the examination weight is a variant"),
    list(choice = "dietary recall", decision = "day one (the MEC recall), reliable recalls only (status 1; 1999-2000 names it DRDDRSTS)",
         reason = "the Methods: intake 'recorded in the mobile examination center'"),
    list(choice = "DII parameters", decision = "the 26 the paper lists, with folic acid as DR1TFA (dii() as the library builds it)",
         reason = "the parameter list names total folate, but the weighted medians the paper reports (1.39 overall, 1.99 with stroke, 1.37 without) fit folic acid (1.36, 1.99, 1.34) and not total folate (1.19, 1.81, 1.16); total folate is a variant"),
    list(choice = "nutrients a cycle lacks", decision = "left out of those participants' DII: vitamin D in 1999-2006; beta-carotene, vitamin A as RAE, alpha-tocopherol, and folic acid in 1999-2000",
         reason = "NHANES files lack them; the paper kept those cycles and does not say how it scored them; dii() leaves a missing parameter out"),
    list(choice = "quartile cutpoints", decision = "the published ones: Q1 under 0.23, Q4 2.95 or more",
         reason = "Statistical analysis; this DII's own quartiles are lower (0.05, 1.49, 2.67), so the published cutpoints put 28%, 28%, 25%, and 19% of the sample in Q1 to Q4, yet they reproduce the crude and Model I estimates (2.47 and 2.39 against 2.47 and 2.43), while the sample's own quartiles give lower ones (variant)"),
    list(choice = "population", decision = "ages 18-79, not pregnant (RIDEXPRG = 1), a reliable day-one recall, and a yes or no to the stroke question (MCQ160F, asked from age 20)",
         reason = "the paper's exclusions; 44,013 here against 44,019 published, whose age step prints 46,469 where its other counts imply 46,369 (the count here); 1,421 strokes against 1,486"),
    list(choice = "smoking", decision = "current smoker (100 cigarettes in life and smoking every day or some days) against not",
         reason = "the text calls everyone who smoked 100 cigarettes a smoker, but Table 1's smoking shares (22.3% overall, 29.6% with stroke, 22.1% without) are those of current smoking (22.6%, 30.3%, 22.4% weighted here), not of ever smoking (46.8%, 59.9%, 46.5%); the text's version is a variant"),
    list(choice = "drinking", decision = "12 or more drinks in any one year or in life (ALQ100, ALD100, or ALQ101, then ALQ110); from 2017, which asks neither, ever having had a drink (ALQ111)",
         reason = "the text says 12 drinks in the year before the survey, but Table 1's drinker shares (89.3% without stroke, 84.3% with) are matched by this (89.4%, 83.3%) and not by 12 drinks in the past year (59.5%, 33.2%); the text's version is a variant"),
    list(choice = "education", decision = "under 9th grade, 9th grade to high school graduate or GED, more than high school",
         reason = "Table 1's shares (5.4%, 35.2%, 59.3%) are those of these groups (5.3%, 35.3%, 59.4%), not of under high school (16.5%)"),
    list(choice = "race", decision = "RIDRETH1's five groups, non-Hispanic White as reference", reason = "Covariates; the reference is unstated"),
    list(choice = "age and BMI", decision = "continuous", reason = "unstated; Table 1 reports age as a mean"),
    list(choice = "diabetes", decision = "told by a doctor, HbA1c 6.5% or more, fasting glucose 126 mg/dL or more, or two-hour OGTT glucose 200 mg/dL or more (2005-2016 only); a test a person lacks counts as not met",
         reason = "the paper's definition; 12.0% weighted here against Table 1's 12.3%"),
    list(choice = "hypertension", decision = "mean of the first three readings 140/90 mmHg or more, told (BPQ020), or taking prescribed medicine (BPQ050A, BPQ150 in 2021-2023)",
         reason = "the paper's definition; 2021-2023 measured blood pressure with an oscillometric device"),
    list(choice = "missing covariates", decision = "complete case (41,097 of 44,013; drinking is missing for 2,423, BMI for 536)",
         reason = "unstated; keeping missing values as their own level changes the estimate little (variant)")
  )
)