# Dietary Inflammatory Index (DII), tertile 3 vs tertile 1, and hyperuricemia, NHANES 2005-2018.
# Wang, Qin, Li, Zhang, and Zeng (2023), Front Nutr, doi:10.3389/fnut.2023.1218166.
# Headline: OR 1.31 (95% CI 1.19-1.44) for the highest against the lowest DII tertile, Model 3
# (abstract and Table 3).

# Day-one recall intakes of every nutrient the DII draws on, for those whose recall is reliable
# (status 1). dietary_totals() stops on a cycle whose file lacks a nutrient it is asked for
# (vitamin D before 2007), 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 DII_NUTRIENTS) {
    column <- paste0(prefix, v)
    out[[v]] <- if (has(d, column)) ifelse(status %in% 1, d[[column]], NA) else NA_real_
  }
  out
}

# The DII from parameter values named as DII_PARAMETERS names them (as dii_parameters() returns
# them), so that a paper's own grouping of a parameter can replace the library's. A parameter
# missing for a person is left out of that person's sum, as dii() leaves it out.
dii_from_values <- function(values, parameters = DII_PARAMETERS$parameter) {
  score <- rep(0, nrow(values))
  for (p in parameters) {
    k <- DII_PARAMETERS[DII_PARAMETERS$parameter == p, ]
    contribution <- (2 * stats::pnorm((values[[p]] - k$mean) / k$sd) - 1) * k$effect
    score <- score + ifelse(is.na(contribution), 0, contribution)
  }
  score
}

# 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))


# The analysis frame for a cycle. `kidney` is "ckd-epi" (CKD-EPI 2009, the chosen version) or
# "cockcroft" (the Cockcroft-Gault formula the Methods print), for both the eGFR exclusion and the
# eGFR covariate.
row062_build <- function(cycle, kidney = "ckd-epi") {
  demo <- demographics(cycle)
  dm <- component("DEMO", cycle)
  bio <- biochemistry(cycle)
  d <- merge_all(demo, dii_intakes(cycle), bio[, c("SEQN", "uric_acid", "creatinine", "glucose_serum")],
                 body_measures(cycle)[, c("SEQN", "bmi", "weight")], smoking(cycle), met_minutes(cycle),
                 diabetes_questions(cycle)[, c("SEQN", "told_diabetes", "insulin", "pills")], hba1c(cycle), fasting_glucose(cycle),
                 bp_questions(cycle)[, c("SEQN", "told_hypertension", "bp_medication")], blood_pressure(cycle)[, c("SEQN", "sbp", "dbp")],
                 hyperlipidemia_status(cycle, demo))
  names(d)[names(d) == "weight"] <- "body_weight"
  alq <- alcohol(cycle)
  d$drinking5 <- alcohol5(alq[match(d$SEQN, alq$SEQN), ], d$sex)
  d$ogtt <- if (has_file("OGTT", cycle)) { o <- component("OGTT", cycle, "LBXGLT"); o$LBXGLT[match(d$SEQN, o$SEQN)] } else NA_real_

  # The DII from 28 parameters with Shivappa and colleagues' (2014) global means, SDs, and effect
  # scores. The paper's Table 2 (each component's score, median and IQR, by sex) fixes how it
  # grouped the fatty acids: n-3 from EPA, DPA, and DHA (20:5, 22:5, 22:6) and n-6 from 18:2,
  # 18:3, 18:4, and 20:4. The library's grouping, with alpha-linolenic acid (18:3) in n-3, is a
  # variant.
  values <- dii_parameters(d)
  d$dii_library <- dii_from_values(values)
  values$N3FAT <- d$P205 + d$P225 + d$P226
  values$N6FAT <- d$P182 + d$P183 + d$P184 + d$P204
  d$dii <- dii_from_values(values)

  d$hyperuricemia <- as.integer(ifelse(d$sex == 1, d$uric_acid >= 7, d$uric_acid >= 6))
  creatinine <- standard_creatinine(d$creatinine, cycle)
  d$egfr <- if (kidney == "ckd-epi") egfr(creatinine, d$age, d$sex, d$race, "2009") else
    (140 - d$age) * d$body_weight * ifelse(d$sex == 1, 1.23, 1.03) / (creatinine * 88.42)  # creatinine in umol/L
  # Education: DMDEDUC2 for adults 20 and over; 18- and 19-year-olds by their grade (DMDEDUC3:
  # 0-12, 55, 66 under high school, 13-14 high school or GED, 15 more). 2021-2023 has no DMDEDUC3,
  # so the harmonized version uses DMDEDUC2 alone.
  grade <- if (has(dm, "DMDEDUC3")) dm$DMDEDUC3[match(d$SEQN, dm$SEQN)] else NA
  young <- ifelse(grade %in% c(0:12, 55, 66), 1, ifelse(grade %in% 13:14, 2, ifelse(grade %in% 15, 3, NA)))
  levels3 <- c("under high school", "high school", "over high school")
  d$education3 <- factor(ifelse(!is.na(d$education), education3(d$education), young), levels = 1:3, labels = levels3)
  d$education3_adults <- factor(education3(d$education), levels = 1:3, labels = levels3)
  d$sex <- factor(d$sex, levels = c(1, 2), labels = c("male", "female"))
  d$race4 <- factor(ifelse(d$race == 3, "white", ifelse(d$race == 4, "black", ifelse(d$race == 1, "mexican", ifelse(d$race %in% c(2, 5), "other", NA)))),
                    levels = c("white", "black", "mexican", "other"))
  d$smoking <- factor(d$smoking, levels = 1:3, labels = c("never", "former", "current"))
  d$drinking <- factor(d$drinking5, levels = 1:5, labels = c("never", "former", "mild", "moderate", "heavy"))
  d$met3 <- cut(d$met, c(-Inf, 600, 1200, Inf), right = FALSE, labels = c("low", "moderate", "vigorous"))
  # Diabetes: told by a doctor, HbA1c of 6.5% or more, fasting glucose of 126 mg/dL or more,
  # random (serum) glucose of 200 mg/dL or more, two-hour OGTT glucose of 200 mg/dL or more, or
  # insulin or diabetes pills. A test a person lacks counts as not met.
  base <- d$told_diabetes %in% 1 | (d$hba1c >= 6.5) %in% TRUE | (d$glucose >= 126) %in% TRUE | (d$glucose_serum >= 200) %in% TRUE |
    d$insulin %in% 1 | d$pills %in% 1
  d$diabetes_harmonized <- as.integer(base)
  d$diabetes <- as.integer(base | (d$ogtt >= 200) %in% TRUE)
  # Hypertension: mean measured pressure of 140/90 mmHg or more, told, or taking prescribed 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))
  pregnant <- if (has(dm, "RIDEXPRG")) dm$RIDEXPRG[match(d$SEQN, dm$SEQN)] %in% 1 else FALSE
  d$in_population <- d$recall_day1 %in% 1 & !is.na(d$uric_acid) & d$age >= 18 & !pregnant & (d$egfr >= 60) %in% TRUE
  d
}

association <- list(
  id = "row062", row = 62, doi = "10.3389/fnut.2023.1218166",
  cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTDRD1", blood_file = "BIOPRO",
  # The paper never mentions weights; it used EmpowerStats, its Table 1 is unweighted (42% non-
  # Hispanic White), and its crude and age-sex-BMI estimates are reproduced unweighted.
  weighted = FALSE,
  family = "logistic", term = "dii_tT3",
  published = list(measure = "OR", estimate = 1.31, low = 1.19, high = 1.44, n = 31781, events = 5491,
                   contrast = "DII tertile 3 (2.60 or more) vs tertile 1 (under 0.79)"),
  left_out = c("physical activity (MET-minutes a week over work, transport, and leisure): 2021-2023 asks about leisure-time activity only",
               "the two-hour OGTT glucose part of the diabetes definition: 2021-2023 has no OGTT",
               "education of 18- and 19-year-olds (DMDEDUC3): 2021-2023 records education from age 20 only, so the harmonized version, like a complete-case analysis there, leaves 18- and 19-year-olds out"),
  build = function(cycle) row062_build(cycle),
  derive = function(data, constants) {
    tertile <- function(x) cut(x, c(-Inf, 0.79, 2.60, Inf), right = FALSE, labels = c("T1", "T2", "T3"))
    data$dii_t <- tertile(data$dii)
    data$dii_library_t <- tertile(data$dii_library)
    u <- function(x) with_unclear(x, levels(x))
    data$education3_u <- u(data$education3)
    data$smoking_u <- u(data$smoking)
    data$drinking_u <- u(data$drinking)
    data$met3_u <- u(data$met3)
    data
  },
  formula = hyperuricemia ~ dii_t + age + sex + bmi + race4 + education3 + smoking + drinking + met3 + egfr + diabetes + hypertension + hyperlipidemia,
  formula_harmonized = hyperuricemia ~ dii_t + age + sex + bmi + race4 + education3_adults + smoking + drinking + egfr + diabetes_harmonized +
    hypertension + hyperlipidemia,
  variants = list(
    list(label = "weighted (WTDRD1)", weighted = TRUE),
    list(label = "T2 vs T1 (published 1.17, 1.07-1.29)", term = "dii_tT2"),
    list(label = "per unit of DII (published 1.06, 1.04-1.09)", formula = hyperuricemia ~ dii + age + sex + bmi + race4 + education3 + smoking +
           drinking + met3 + egfr + diabetes + hypertension + hyperlipidemia, term = "dii"),
    list(label = "Model 1, crude (published 1.21, 1.12-1.30)", formula = hyperuricemia ~ dii_t, sample = "own"),
    list(label = "Model 2, age sex BMI (published 1.26, 1.17-1.36)", formula = hyperuricemia ~ dii_t + age + sex + bmi, sample = "own"),
    list(label = "Model 1, crude, weighted", formula = hyperuricemia ~ dii_t, sample = "own", weighted = TRUE),
    list(label = "Model 2, weighted", formula = hyperuricemia ~ dii_t + age + sex + bmi, sample = "own", weighted = TRUE),
    list(label = "missing covariates as their own level", formula = hyperuricemia ~ dii_t + age + sex + bmi + race4 + education3_u + smoking_u +
           drinking_u + met3_u + egfr + diabetes + hypertension + hyperlipidemia, sample = "own"),
    list(label = "library DII grouping (alpha-linolenic acid in n-3)", formula = hyperuricemia ~ dii_library_t + age + sex + bmi + race4 + education3 +
           smoking + drinking + met3 + egfr + diabetes + hypertension + hyperlipidemia, term = "dii_library_tT3"),
    list(label = "Cockcroft-Gault as printed (exclusion and covariate)", build = function(cycle) row062_build(cycle, kidney = "cockcroft"))
  ),
  # 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 = "paper",
         detail = "The Methods name the CKD-EPI creatinine equation for eGFR, an exclusion step and a Model 3 covariate, but print a Cockcroft-Gault formula with creatinine in mmol/L; Table 1's eGFR means and SDs (98.23, SD 19.26 in men; 101.22, SD 20.45 in women) and the 2,956 excluded are those of CKD-EPI 2009 (2,954 excluded here), while the printed formula gives means near 123 with SDs near 44."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The DII's n-3 and n-6 fatty acids were summed as EPA, DPA, and DHA (n-3) and 18:2, 18:3, 18:4, and 20:4 (n-6), counting alpha-linolenic acid (18:3, an n-3 fat) as n-6: so grouped, the DII reproduces Table 2's component medians, Table 1's means by sex (1.15 and 1.87), and the tertile cutpoints (0.79 and 2.60; 0.786 and 2.599 here), while 18:3 counted as n-3 gives means of 0.74 and 1.52."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The text defines diabetes by a doctor's diagnosis, HbA1c greater than 6.5%, or fasting, random, or two-hour OGTT glucose, which gives 16.06% of men here, but Table 1's 16.53% is matched only with HbA1c of 6.5% or more and insulin or diabetes pills added (16.53% here)."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE, evidence = "unresolved",
         detail = "Table 3 gives n = 31,781 for every model, but Model 3's interval for tertile 3 is too wide for that many: it is 27% wider on the log scale than Model 2's (1.19 to 1.44 against 1.17 to 1.36), where keeping everyone (missing covariates as their own level) widens it by 5% here (1.14 to 1.34 against 1.17 to 1.36).")
  ),
  choices = list(
    list(choice = "weights", decision = "unweighted", reason = "the paper mentions no weights and used EmpowerStats; Table 1 is unweighted, and only unweighted fits reproduce Models 1 and 2 (1.21 and 1.26, with their intervals)"),
    list(choice = "dietary recall", decision = "day one, reliable recalls (DR1DRSTZ = 1)", reason = "unstated; the 9,549 excluded for lacking dietary data are exactly those without a reliable day-one recall"),
    list(choice = "DII parameters", decision = "the library's 28 (folic acid as DR1TFA) with n-3 as EPA + DPA + DHA and n-6 as 18:2 + 18:3 + 18:4 + 20:4", reason = "the paper names the parameters but not their nutrients; so built, the DII matches Table 2's component medians and IQRs, Table 1's mean and SD by sex (1.15, SD 1.88; 1.87, SD 1.80), the range (-5.28 to 5.79), and the tertile cutpoints (0.786 and 2.599 here); folic acid fits Table 2's folate scores and total folate does not, and with alpha-linolenic acid in n-3 the means are 0.74 and 1.52"),
    list(choice = "vitamin D in 2005-2006", decision = "left out of those participants' DII", reason = "NHANES dietary files have vitamin D only from 2007"),
    list(choice = "tertile cutpoints", decision = "the published 0.79 and 2.60, T2 from 0.79 and T3 from 2.60", reason = "Methods and Table 3; the sample's own tertiles of this DII are 0.786 and 2.599 (10,617, 10,579, and 10,588 people at the published cutpoints against 10,594, 10,593, and 10,594)"),
    list(choice = "eGFR", decision = "CKD-EPI 2009, with 2005-2006 creatinine recalibrated as NCHS directs, for both the exclusion (under 60) and the covariate", reason = "the Methods name CKD-EPI but print a Cockcroft-Gault formula (in mmol/L); Table 1's eGFR means and SDs by sex (98.2, SD 19.3; 101.2, SD 20.5) and the 2,956 excluded match CKD-EPI 2009 with the recalibration (2,954 here; 3,031 without it); Cockcroft-Gault gives means near 123 with SDs near 44"),
    list(choice = "pregnancy", decision = "RIDEXPRG = 1", reason = "unstated; excludes 641 (published 642)"),
    list(choice = "drinking", decision = "Rattan's five groups as alcohol5() builds them (never, former, mild, moderate, heavy)", reason = "the paper's definition; Table 1's shares are close (men: heavy 25.4% here against 25.9%, moderate 12.6% against 12.4%, mild 36.4% against 38.1%, former 17.4% against 15.2%); counting ever drinking 4 or 5 drinks almost every day (ALQ150 or ALQ151) as heavy would make 33.0% heavy"),
    list(choice = "physical activity", decision = "GPAQ MET-minutes a week (8 METs vigorous, 4 moderate and transport) over work, transport, and leisure, in the paper's three groups; missing in 2005-2006, which asked no GPAQ", reason = "the paper's groups; no definition tried reproduces Table 1's shares, and leaving MET out moves the estimate by under 1%"),
    list(choice = "education", decision = "under high school (DMDEDUC2 1-2), high school or GED (3), more (4-5); 18- and 19-year-olds by DMDEDUC3", reason = "Table 1's three groups; under high school matches (26.1% of men)"),
    list(choice = "race", decision = "non-Hispanic White, non-Hispanic Black, Mexican American, other (other Hispanic with other races)", reason = "Table 1; shares match"),
    list(choice = "diabetes", decision = "told, HbA1c 6.5% or more, fasting glucose 126 mg/dL or more, random glucose 200 or more, OGTT 200 or more, or insulin or diabetes pills", reason = "the Methods list the first five (HbA1c 'greater than' 6.5%); Table 1's share in men (16.53%) is matched only with 6.5% or more and the medicines (16.53%; 16.06% as printed)"),
    list(choice = "hypertension", decision = "mean of the available readings 140/90 or more, told, or taking prescribed medicine (BPQ050A)", reason = "unstated readings; Table 1's shares match (38.0% of men, 36.4% of women)"),
    list(choice = "hyperlipidemia", decision = "NCEP ATP III thresholds with fasting triglycerides and LDL, or cholesterol medicine (hyperlipidemia_status)", reason = "unstated sources; Table 1's shares are close (65.8% vs 67.0% of men); serum triglycerides overshoot (70.3%)"),
    list(choice = "missing covariates", decision = "complete case (25,142 of 31,784)", reason = "unstated; Table 3 heads every model n = 31,781, but its Model 3 interval is wider than Model 2's by more than the covariates alone explain, as a smaller complete-case sample would make it; keeping missing values as their own level changes the estimate little (variant)")
  )
)