# Dietary selenium intake (24-hour recall, day 1) and chronic kidney disease, adults aged 30-80,
# NHANES 2015-2018. Pi, Liao, Song, Cao, Tang, Lin, and Zhong (2024), Front Nutr,
# doi:10.3389/fnut.2024.1396470.
# Headline: OR 0.77 (95% CI 0.63-0.95), p = 0.016, selenium quartile 4 (over 0.144 mg/day) vs
# quartile 1 (0.072 mg/day or less), Model 2 (Table 2).

# eGFR (mL/min/1.73 m2) by the 4-variable MDRD Study equation for standardized creatinine (Levey
# 2006, the paper's reference 26): 175 x creatinine^-1.154 x age^-0.203, x 0.742 for women, and x
# 1.212 for Black people when `race_term` is TRUE. 2015-2023 creatinine is already standardized.
row315_mdrd <- function(creatinine, age, sex, black, race_term = FALSE) {
  175 * creatinine^-1.154 * age^-0.203 * ifelse(sex == 2, 0.742, 1) * ifelse(race_term & black, 1.212, 1)
}

# CKD: eGFR under 60 or urine albumin-to-creatinine ratio of 30 mg/g or more; unknown without both.
row315_ckd <- function(egfr, acr) ifelse(is.na(egfr) | is.na(acr), NA, as.integer(egfr < 60 | acr >= 30))

# Selenium from dietary supplements: any in the past 30 days (DSQTSELE; the file leaves it missing
# for those who took no supplement), or in the 24 hours the dietary recall covers (DS1TSELE, not
# released for 2021-2023).
row315_supplements <- function(cycle) {
  q <- component("DSQTOT", cycle)
  out <- data.frame(SEQN = q$SEQN, supplement_30d = as.integer(q$DSQTSELE > 0 & !is.na(q$DSQTSELE)))
  if (cycle != REPLICATION_CYCLE) {
    s <- component("DS1TOT", cycle)
    out$supplement_24h <- as.integer(s$DS1TSELE > 0 & !is.na(s$DS1TSELE))[match(out$SEQN, s$SEQN)]
  } else {
    out$supplement_24h <- NA
  }
  out
}

# Self-reported BMI from current height (inches, WHD010) and weight (pounds, WHD020).
row315_self_bmi <- function(cycle) {
  w <- component("WHQ", cycle, c("WHD010", "WHD020"))
  height <- ifelse(w$WHD010 %in% c(7777, 9999), NA, w$WHD010)
  weight <- ifelse(w$WHD020 %in% c(7777, 9999), NA, w$WHD020)
  data.frame(SEQN = w$SEQN, bmi_self = 703.0696 * weight / height^2)
}

row315_build <- function(cycle) {
  d <- merge_all(demographics(cycle), dietary_totals(cycle, c("SELE", "CARB"), days = 1), row315_supplements(cycle),
                 biochemistry(cycle)[, c("SEQN", "creatinine", "uric_acid", "tg_serum")], urine_albumin_creatinine(cycle),
                 crp(cycle), total_cholesterol(cycle), body_measures(cycle)[, c("SEQN", "bmi")], row315_self_bmi(cycle),
                 bp_questions(cycle), blood_pressure(cycle)[, c("SEQN", "sbp", "dbp")],
                 diabetes_questions(cycle)[, c("SEQN", "told_diabetes")], component("SMQ", cycle, "SMQ020"))
  # CKD by MDRD without its race term, which the paper's counts show it used (23.8% with CKD, and
  # Table 1's CKD shares by race group); ckd_race_term applies the 1.212 for Black people.
  d$ckd <- row315_ckd(row315_mdrd(d$creatinine, d$age, d$sex, d$race == 4), d$acr)
  d$ckd_race_term <- row315_ckd(row315_mdrd(d$creatinine, d$age, d$sex, d$race == 4, race_term = TRUE), d$acr)
  # Day-one selenium (mcg) in the paper's quartiles of mg/day: up to 0.072, 0.103, 0.144, and over.
  d$selenium_q <- cut(d$SELE, c(-Inf, 72, 103, 144, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))
  d$sex <- factor(d$sex)
  # Four race groups, as Table 1 has: Mexican American, non-Hispanic White, non-Hispanic Black, and
  # everyone else (other Hispanic and other races). Table 1's counts match these groups under other
  # labels (see choices).
  d$race4 <- factor(ifelse(d$race == 1, "Mexican American", ifelse(d$race == 3, "NH White", ifelse(d$race == 4, "NH Black",
                    ifelse(d$race %in% c(2, 5), "other", NA)))), levels = c("NH White", "Mexican American", "NH Black", "other"))
  d$living <- factor(ifelse(d$marital %in% 1, "with a partner", ifelse(d$marital %in% 2:3, "alone", NA)), levels = c("with a partner", "alone"))
  d$education <- factor(education3(d$education), levels = 1:3)
  d$hypertension <- factor(d$told_hypertension, levels = 0:1)
  high_bp <- cbind(d$told_hypertension == 1, d$sbp >= 140, d$dbp >= 90)
  d$hypertension_measured <- factor(ifelse(rowSums(high_bp, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(high_bp)) > 0, 0, NA)), levels = 0:1)
  d$diabetes <- factor(d$told_diabetes, levels = 0:1)
  d$supplement <- factor(d$supplement_30d, levels = 0:1)
  d$supplement_24h <- factor(d$supplement_24h, levels = 0:1)
  d$smoker <- factor(ifelse(d$SMQ020 %in% 1, 1, ifelse(d$SMQ020 %in% 2, 0, NA)), levels = 0:1)
  # Aged 30-80 (80 is 80 and over), not pregnant, with a reliable day-one recall, CKD known, and
  # every covariate the paper lists (Model 2's, and smoking and the income-to-poverty ratio), as
  # Table 1 has none missing.
  covariates <- c("sex", "race4", "living", "education", "bmi", "hypertension", "diabetes", "CARB", "supplement",
                  "crp", "tc", "tg_serum", "uric_acid", "pir", "smoker")
  d$in_population <- d$age >= 30 & d$age <= 80 & !(d$pregnant %in% 1) & !is.na(d$SELE) & !is.na(d$ckd) &
    stats::complete.cases(d[, covariates])
  d
}

ROW315_MODEL1 <- ckd ~ selenium_q + age + sex + race4 + living + education + bmi + hypertension + diabetes
ROW315_MODEL2 <- update(ROW315_MODEL1, . ~ . + CARB + supplement + crp + tc + tg_serum + uric_acid)
# Model 2 as the paper computed it: its adjusted estimates are reproduced only without age.
ROW315_CHOSEN <- update(ROW315_MODEL2, . ~ . - age)

association <- list(
  id = "row315", row = 315, doi = "10.3389/fnut.2024.1396470",
  cycles = c("2015-2016", "2017-2018"),
  # Unweighted, as the paper's estimates are (see choices); the dietary day-one weight is the one
  # NCHS directs for analyses of the day-one recall, used by the weighted variant and the weighted
  # replication.
  weight = "WTDRD1", blood_file = "BIOPRO", weighted = FALSE,
  family = "logistic", term = "selenium_qQ4",
  published = list(measure = "OR", estimate = 0.77, low = 0.63, high = 0.95, n = 6390, events = 1523,
                   contrast = "dietary selenium quartile 4 (over 0.144 mg/day) vs quartile 1 (0.072 mg/day or less)"),
  left_out = character(),
  build = function(cycle) row315_build(cycle),
  formula = ROW315_CHOSEN,
  variants = list(
    list(label = "weighted (WTDRD1)", weighted = TRUE),
    list(label = "crude (published 0.66, 0.56-0.78)", formula = ckd ~ selenium_q),
    list(label = "Model 1 without age (published 0.77, 0.64-0.92)", formula = update(ROW315_MODEL1, . ~ . - age)),
    list(label = "Model 2 with age, as the text describes", formula = ROW315_MODEL2),
    list(label = "Model 1 with age, as the text describes", formula = ROW315_MODEL1),
    list(label = "CKD by MDRD with its race term", formula = update(ROW315_CHOSEN, ckd_race_term ~ .)),
    list(label = "self-reported BMI (WHD010, WHD020)", formula = update(ROW315_CHOSEN, . ~ . - bmi + bmi_self)),
    list(label = "supplement selenium in the recall's 24 hours (DS1TSELE)", formula = update(ROW315_CHOSEN, . ~ . - supplement + supplement_24h)),
    list(label = "adding smoking and income (Figure 2's adjustment)", formula = update(ROW315_CHOSEN, . ~ . + smoker + pir)),
    list(label = "hypertension also by measured pressure 140/90 or more", formula = update(ROW315_CHOSEN, . ~ . - hypertension + hypertension_measured))
  ),
  # 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 = "model", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods and Table 2's note adjust Models 1 and 2 for age, but their estimates are reproduced only without age (quartile 4 vs 1: 0.76 and 0.76 against 0.77 and 0.77; with age 0.92 and 0.86), as are the per mg/day estimates (0.22 and 0.20 against 0.23 and 0.23; with age 0.73 and 0.43)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods compute eGFR with the 2006 MDRD Study equation (reference 26), which multiplies by 1.212 for Black people, but the paper's CKD share (23.8%) and Table 1's CKD shares by race are reproduced without that term (23.8%; non-Hispanic Black, the row Table 1 labels Others, 32.0% against 32.1%), not with it (21.5%; 21.1%)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods define hypertension by a doctor's diagnosis and systolic pressure of 140 mmHg or more or diastolic of 90 or more, but Table 1's shares (35.1% without CKD, 65.5% with) are reproduced by the diagnosis alone (35.5%, 65.3%), not with measured pressure added (44%, 74%)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's race labels are rotated: by their counts and CKD shares, its Non-Hispanic white row is other Hispanic and other races, its Non-Hispanic black row non-Hispanic White, and its Others row non-Hispanic Black."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract says the risk of CKD falls 7.7% for every additional 0.1 mg of selenium, but Table 2's odds ratios per mg/day (0.23 in Models 1 and 2, reproduced per mg/day) correspond to 0.86 per 0.1 mg, 14% lower odds."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2's per-SD row gives Model 1 as 0.93 (0.81-1.05) with p = 0.025, though its interval includes 1, and Model 2 as 0.93 (0.83-1.09) with p = 0.059, an interval not centered on its estimate on the log scale.")
  ),
  choices = list(
    list(choice = "cycles", decision = "2015-2016 and 2017-2018", reason = "the paper's 19,225 participants are DEMO_I and DEMO_J exactly"),
    list(choice = "weights", decision = "unweighted",
         reason = "the paper (SPSS) names no weights or design; its crude ORs are the unweighted ones of Table 2's counts (0.65 here against 0.66) and its CIs have unweighted widths; weighted variant"),
    list(choice = "age in the models", decision = "left out, as the paper's computation left it out",
         reason = "the text says the models adjust for age, but its adjusted estimates are reproduced only without age: Model 1 0.76 (0.63-0.90) against 0.77 (0.64-0.92), Model 2 0.76 (0.62-0.93) against 0.77 (0.63-0.95), and per mg/day 0.22 and 0.20 against 0.23 and 0.23 (0.73 and 0.43 with age); the published estimate comes from the models without age, and the text's models are variants"),
    list(choice = "CKD", decision = "MDRD eGFR (175, standardized creatinine) under 60 without the equation's race term, or ACR 30 mg/g or more",
         reason = "the paper cites the 2006 MDRD equation without its terms; without the race term 23.8% have CKD, as in the paper (21.5% with it), and the CKD shares of the four race groups (NH White 26.2%, Mexican American 19.6%, NH Black 32.0%, other 16.4%) match Table 1's (25.9, 20.2, 32.1, 16.4%), while with the term NH Black's is 21.1%; no creatinine correction, since 2015-2018 creatinine is standardized and the NHANES III correction the paper cites is for older surveys"),
    list(choice = "exposure", decision = "day-one recall selenium (DR1TSELE) in the paper's quartiles (up to 72, 103, 144 mcg, and over)",
         reason = "the Methods describe one day's recall, and day-one quartile shares (26.2, 25.3, 25.0, 23.5%) match Table 2's (25.8, 25.4, 24.9, 23.9%) while the two-day mean's (22.6, 27.6, 28.8, 21.0%) do not"),
    list(choice = "population", decision = "aged 30 to 80 (80 is 80 and over), not pregnant (RIDEXPRG), reliable day-one recall, serum creatinine and urine ACR, and every listed covariate known",
         reason = "the paper's exclusions for basic, laboratory, and survey information are unspecified; this gives 6,738 against its 6,390, with the same CKD share; leaving out age 80 gives 21.0% with CKD"),
    list(choice = "Model 2 covariates", decision = "Table 2's, without age (see above): sex, race, marital status, education, BMI, hypertension, diabetes, carbohydrate, selenium supplements, CRP, total cholesterol, triglycerides, uric acid",
         reason = "the headline is Table 2's; Figure 2's model adds smoking and the income-to-poverty ratio (variant)"),
    list(choice = "race", decision = "four groups: Mexican American, NH White, NH Black, everyone else (other Hispanic and other races)",
         reason = "Table 1 has four groups whose counts and CKD shares match these, though under rotated labels (its 'Non-Hispanic white' is other Hispanic and other races, its 'Non-Hispanic black' NH White, its 'Others' NH Black); the text's CKD patients 'more likely non-Hispanic whites' fits this reading"),
    list(choice = "hypertension", decision = "told by a doctor (BPQ020)", reason = "Table 1's 35.1% and 65.5% fit told alone (35.5%, 65.3%), not told or measured pressure 140/90 or more (44%, 74%; variant)"),
    list(choice = "diabetes", decision = "told by a doctor (DIQ010 = 1; borderline counts as no)", reason = "the paper's definition; Table 1's 13.2% and 34.7% against 13.0% and 33.4%"),
    list(choice = "selenium supplements", decision = "any selenium from supplements in the past 30 days (DSQTSELE above 0)",
         reason = "unstated; Table 1's 10.1% and 12.0% fit neither this (22.5%, 25.4%) nor the recall's 24 hours (16.4%, 21.5%), which gives the same estimate (variant); 2021-2023 releases only the 30-day file"),
    list(choice = "BMI", decision = "measured (BMXBMI)",
         reason = "Table 1's means (28.74, 29.59) are below measured BMI here (29.8, 31.0) and nearer self-reported BMI (29.0, 29.9), but neither matches; measured BMI is the examination's and gives the estimate nearer the published one (variant: self-reported, 0.89)"),
    list(choice = "covariate coding", decision = "age, BMI, carbohydrate (g), CRP (mg/L), total cholesterol, refrigerated-serum triglycerides, and uric acid continuous; living with a partner or alone; education in three levels",
         reason = "Table 1's categories; its triglycerides (1.75, 1.90 mmol/L) are the refrigerated-serum ones (1.75, 1.91) on a sample too large for the fasting subsample; units do not change the selenium estimate")
  )
)