# Blood cadmium and depression (PHQ-9 of 10 or more) in women aged 20 and older, NHANES 2005-2016.
# Ji and Wang (2024), Environ Health Prev Med, doi:10.1265/ehpm.24-00050.
# Headline: OR 1.33 (95% CI 1.21-1.45) per unit of ln blood cadmium, Model 3 (Table 2).

# NCHS's weights for blood metals in these cycles: the examination weight in 2005-2012, when every
# examined participant had them measured, and the blood metal subsample weight (WTSH2YR) in
# 2013-2016, when only a random half of those aged 12 and older did. A variant's derive puts them in
# place of the examination weight, each cycle's weight times its share of the pooled years.
subsample_weights <- function(data, constants) {
  cycles <- unique(data$cycle)
  years <- stats::setNames(vapply(cycles, function(cycle) cycle_info(cycle)$years, 0), cycles)
  rows <- !is.na(data$metal_weight)
  data$w[rows] <- data$metal_weight[rows] * years[data$cycle[rows]] / sum(years)
  data
}

association <- list(
  id = "row106", row = 106, doi = "10.1265/ehpm.24-00050",
  cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016"),
  weight = "WTMEC2YR", blood_file = "PBCD",
  # The paper names no weight and says it accounted for the survey design, but its crude and
  # age-race estimates are reproduced exactly only without weights (variants below).
  weighted = FALSE,
  family = "logistic", term = "ln_cadmium",
  published = list(measure = "OR", estimate = 1.33, low = 1.21, high = 1.45, n = 10868, events = 1173,
                   contrast = "per unit of ln blood cadmium (ug/L)"),
  left_out = c("marital status as married against everything else: 2021-2023's DMDMARTZ puts living with a partner together with married, so the harmonized version counts married or living with a partner as yes"),
  build = function(cycle) {
    demo <- demographics(cycle)
    dm <- component("DEMO", cycle)
    pbcd <- component("PBCD", cycle)
    d <- merge_all(demo, phq9(cycle), pbcd[, c("SEQN", "LBXBCD", "LBXBPB", "LBXTHG")], body_measures(cycle)[, c("SEQN", "bmi")],
                   component("SMQ", cycle, "SMQ020"), component("BPQ", cycle, "BPQ020"), component("DIQ", cycle, "DIQ010"))
    d$metal_weight <- if (has(pbcd, "WTSH2YR")) {
      x <- pbcd$WTSH2YR[match(d$SEQN, pbcd$SEQN)]
      ifelse(is.na(x), 0, x)
    } else NA_real_
    d$depression <- as.integer(d$phq9 >= 10)
    d$ln_cadmium <- log(d$LBXBCD)
    d$race <- factor(d$race)
    d$education <- factor(education3(d$education))
    d$smoking <- factor(yes(d$SMQ020))
    # Married (DMDMARTL 1) against every other answer, as Table 1's shares show.
    d$married <- if (has(dm, "DMDMARTL")) {
      m <- dm$DMDMARTL[match(d$SEQN, dm$SEQN)]
      factor(ifelse(m %in% 1, 1, ifelse(m %in% 2:6, 0, NA)))
    } else NA
    d$partnered <- factor(ifelse(d$marital %in% 1, 1, ifelse(d$marital %in% 2:3, 0, NA)))
    d$hypertension <- factor(yes(d$BPQ020))
    d$diabetes <- factor(ifelse(d$DIQ010 %in% 1:3, d$DIQ010, NA))
    d$diabetes2 <- factor(ifelse(d$DIQ010 %in% 1, 1, ifelse(d$DIQ010 %in% 2:3, 0, NA)))
    d$bmi3 <- cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-29.9", ">=30"))
    d$in_population <- d$sex == 2 & d$age >= 20 & !is.na(d$phq9) &
      !is.na(d$LBXBCD) & !is.na(d$LBXBPB) & !is.na(d$LBXTHG)
    d
  },
  variants = list(
    list(label = "weighted (WTMEC2YR in every cycle)", weighted = TRUE),
    list(label = "weighted as NCHS directs (WTSH2YR in 2013-2016)", weighted = TRUE, derive = subsample_weights),
    list(label = "crude (published 1.56, 1.45-1.69)", formula = depression ~ ln_cadmium),
    list(label = "crude, weighted (WTMEC2YR)", formula = depression ~ ln_cadmium, weighted = TRUE),
    list(label = "age, race (published 1.67, 1.55-1.81)", formula = depression ~ ln_cadmium + age + race),
    list(label = "age, race, weighted (WTMEC2YR)", formula = depression ~ ln_cadmium + age + race, weighted = TRUE),
    list(label = "BMI continuous, as Table 1 summarizes it", formula = depression ~ ln_cadmium + age + race + bmi + pir + education + smoking + married + hypertension + diabetes),
    list(label = "diabetes yes or no (borderline as no)", formula = depression ~ ln_cadmium + age + race + bmi3 + pir + education + smoking + married + hypertension + diabetes2)
  ),
  # 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 = "weighting", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods say the analyses adjusted for the complex multistage cluster survey design, but Table 2's Model 1 (1.56, 1.45-1.69) and Model 2 (1.67, 1.55-1.81) are reproduced to every printed digit only unweighted (Model 1: 1.565, 1.451-1.687; with the examination weight 1.618, 1.465-1.787)."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The abstract reports OR 1.33 (1.21-1.45) for each incremental unit of blood cadmium, but Table 2 and the Results give it per unit of ln-transformed blood cadmium, the scale on which it is reproduced (1.318, 1.202-1.446)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The text defines diabetes as a physician's diagnosis, yes or no, but Tables 1 and 3 show it in three levels (yes, no, borderline); the two codings give nearly the same estimate here (1.318 and 1.317).")
  ),
  choices = list(
    list(choice = "weights", decision = "unweighted", reason = "the paper names no weight and says it accounted for the design, but only unweighted fits reproduce its crude (1.56, 1.45-1.69) and age-race (1.67, 1.55-1.81) estimates to the second decimal, CIs included; Table 1's shares are unweighted; weighted fits are kept as variants, with the examination weight and with the blood metal subsample weight NCHS directs for 2013-2016"),
    list(choice = "missing depression data", decision = "PHQ-9 total missing if any of the nine items is missing, refused, or don't know", reason = "reproduces the flow chart exactly: 29,745 excluded for depression, then 6,265 for metals, 12,309 men, and 723 under 20"),
    list(choice = "blood metals", decision = "cadmium, lead, and mercury all measured (LBXBCD, LBXBPB, LBXTHG); values below detection as NHANES releases them (detection limit over the square root of 2)", reason = "the flow chart excludes those missing any of the three; detection limits are not discussed"),
    list(choice = "exposure scale", decision = "natural log of blood cadmium (ug/L), per unit", reason = "Table 2's title and the Results; the abstract's 'each incremental unit of blood cadmium' is per unit of ln cadmium"),
    list(choice = "pregnancy", decision = "pregnant women kept", reason = "no pregnancy exclusion is listed, and the flow chart and n = 10,868 are reproduced without one"),
    list(choice = "missing covariates", decision = "complete case", reason = "1,026 excluded for missing covariates, reproduced exactly (most for income to poverty ratio)"),
    list(choice = "BMI", decision = "three groups (under 25, 25 to 29.9, 30 or more), as the covariate text says it was stratified", reason = "stated; Table 1 summarizes BMI by its mean, which says nothing of how the model held it, and nothing in the paper's numbers shows the model held it otherwise (continuous BMI is a variant)"),
    list(choice = "income to poverty ratio", decision = "continuous", reason = "Table 1 reports it as a mean; no categories are given"),
    list(choice = "education", decision = "DMDEDUC2 in three levels: less than high school (1-2), high school or GED (3), above high school (4-5)", reason = "Table 1's categories; its shares (23.77, 21.83, 54.41%) are reproduced exactly"),
    list(choice = "smoking", decision = "smoked 100 cigarettes in life (SMQ020) yes or no", reason = "the covariate text; Table 1's 37.15% yes is reproduced exactly"),
    list(choice = "marital status", decision = "married (DMDMARTL 1) against widowed, divorced, separated, never married, or living with a partner", reason = "unstated mapping; Table 1's 47.52% yes is reproduced exactly by married alone (married or living with a partner gives 54.91%)"),
    list(choice = "hypertension", decision = "told by a doctor (BPQ020) yes or no", reason = "the covariate text: self-reported physician diagnosis; Table 1's 35.48% is reproduced exactly"),
    list(choice = "diabetes", decision = "DIQ010 in three levels: yes, no, borderline", reason = "Table 1 and Table 3 show three levels (text says yes or no); the two-level coding is a variant"),
    list(choice = "race", decision = "RIDRETH1 in five groups", reason = "the covariate text; Table 1's shares are reproduced exactly")
  ),
  formula = depression ~ ln_cadmium + age + race + bmi3 + pir + education + smoking + married + hypertension + diabetes,
  formula_harmonized = depression ~ ln_cadmium + age + race + bmi3 + pir + education + smoking + partnered + hypertension + diabetes
)