# Non-HDL cholesterol and depression (PHQ-9 total of 10 or more) in adults, NHANES 2005-2018.
# Zhu et al. (2023), Front Psychiatry, doi:10.3389/fpsyt.2023.1274648.
# Headline: OR 1.22 (95% CI 1.03-1.45), highest vs lowest quartile of non-HDL cholesterol,
# Model 3 (Table 2).

# The paper's models (Statistical analyses), with HDL in Model 3 as Table 2's note has it.
MODEL_1 <- depression ~ nonhdl_q + age + sex + race
MODEL_2 <- update(MODEL_1, . ~ . + marital + pir + education + bmi + smoking + alcohol)
MODEL_3 <- update(MODEL_2, . ~ . + cvd + diabetes + hdl + hypertension)

# The analytic sample: the study population with complete data for Model 3.
analytic <- function(data) {
  variables <- c("nonhdl", setdiff(all.vars(MODEL_3), "nonhdl_q"))
  data$in_population %in% TRUE & stats::complete.cases(data[, variables]) & !is.na(data$w) & data$w > 0
}

# Quartiles of non-HDL cholesterol (mmol/L), intervals closed on the right.
quartile <- function(x, cutpoints) cut(x, c(-Inf, cutpoints, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))

# A weighted quantile: the smallest value at which the cumulative share of the weight reaches p.
weighted_quantile <- function(x, w, probs) {
  o <- order(x)
  share <- cumsum(w[o]) / sum(w)
  vapply(probs, function(p) x[o][which(share >= p)[1]], 0)
}

association <- list(
  id = "row058", row = 58, doi = "10.3389/fpsyt.2023.1274648",
  cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTMEC2YR", blood_file = "TCHOL",
  family = "logistic", term = "nonhdl_qQ4",
  published = list(measure = "OR", estimate = 1.22, low = 1.03, high = 1.45, n = 28041, events = 2419,
                   contrast = "highest vs lowest quartile of non-HDL cholesterol (total minus HDL cholesterol, mmol/L)"),
  left_out = c(
    "marital status in four groups (widowed apart from divorced or separated): 2021-2023 releases only DMDMARTZ's three groups; the harmonized version uses married or living with a partner; widowed, divorced, or separated; never married",
    "two-hour OGTT glucose in the diabetes definition: 2021-2023 (like 2017-2018) has no OGTT; the harmonized version counts diabetes as told by a doctor, insulin or pills, HbA1c 6.5% or more, or fasting glucose 126 mg/dL or more"
  ),
  build = function(cycle) {
    demo <- demographics(cycle)
    raw <- component("DEMO", cycle)
    d <- merge_all(demo, phq9(cycle), total_cholesterol(cycle), hdl_cholesterol(cycle),
                   body_measures(cycle)[, c("SEQN", "bmi")], smoking(cycle), alcohol(cycle)[, c("SEQN", "alcohol3")],
                   medical_conditions(cycle), diabetes_status(cycle), hypertension_status(cycle))
    # Two-hour OGTT glucose: given 2005-2006 through 2015-2016 only.
    has_ogtt <- cycle %in% c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016")
    ogtt <- if (has_ogtt) component("OGTT", cycle, "LBXGLT") else data.frame(SEQN = d$SEQN, LBXGLT = NA_real_)
    d$ogtt <- ogtt$LBXGLT[match(d$SEQN, ogtt$SEQN)]
    d$depression <- as.integer(d$phq9 >= 10)
    d$nonhdl <- cholesterol_mmol(d$tc - d$hdl)
    d$hdl <- cholesterol_mmol(d$hdl)
    d$sex <- factor(d$sex)
    d$race <- factor(ifelse(d$race %in% 3, "white", ifelse(d$race %in% 1, "mexican", ifelse(d$race %in% 4, "black", "other"))),
                     levels = c("white", "mexican", "black", "other"))
    # Marital status: DMDMARTZ's three groups (demographics()) for the harmonized version, and the
    # paper's four from DMDMARTL (2005-2018; 2021-2023 has only DMDMARTZ).
    d$marital3 <- factor(d$marital, levels = 1:3, labels = c("married", "widowed or divorced", "never"))
    m <- if (has(raw, "DMDMARTL")) raw$DMDMARTL[match(d$SEQN, raw$SEQN)] else rep(NA, nrow(d))
    d$marital <- factor(ifelse(m %in% c(1, 6), "married", ifelse(m %in% 5, "never", ifelse(m %in% 2, "widowed", ifelse(m %in% 3:4, "divorced", NA)))),
                        levels = c("married", "never", "widowed", "divorced"))
    d$education <- factor(ifelse(d$education %in% 1:2, "less than high school", ifelse(d$education %in% 3, "high school",
                                 ifelse(d$education %in% 4, "some college", ifelse(d$education %in% 5, "college", NA)))),
                          levels = c("less than high school", "high school", "some college", "college"))
    d$smoking <- factor(d$smoking)
    d$alcohol <- factor(d$alcohol3)
    # Heart failure, angina, heart attack, or stroke, as the paper lists CVD.
    parts <- cbind(d$heart_failure, d$angina, d$heart_attack, d$stroke)
    d$cvd <- factor(ifelse(rowSums(parts == 1, na.rm = TRUE) > 0, 1, ifelse(rowSums(is.na(parts)) == 4, NA, 0)))
    d$diabetes_harmonized <- factor(d$diabetes)
    d$diabetes <- factor(ifelse(d$diabetes %in% 1 | (d$ogtt >= 200) %in% TRUE, 1, d$diabetes))
    d$hypertension <- factor(d$hypertension)
    d$in_population <- d$age >= 18 & !is.na(d$phq9) & !is.na(d$nonhdl)
    d
  },
  constants = function(data) {
    sample <- analytic(data)
    list(cutpoints = unname(stats::quantile(data$nonhdl[sample], c(0.25, 0.5, 0.75))),
         cutpoints_weighted = weighted_quantile(data$nonhdl[sample], data$w[sample], c(0.25, 0.5, 0.75)))
  },
  derive = function(data, constants) {
    data$nonhdl_q <- quartile(data$nonhdl, constants$cutpoints)
    data
  },
  formula = MODEL_3,
  formula_harmonized = update(MODEL_3, . ~ . - marital - diabetes + marital3 + diabetes_harmonized),
  variants = list(
    list(label = "crude (published 1.30, 1.14-1.50)", formula = depression ~ nonhdl_q),
    list(label = "Model 1 (published 1.40, 1.22-1.62)", formula = MODEL_1),
    list(label = "Model 2 (published 1.13, 0.96-1.33)", formula = MODEL_2),
    list(label = "unweighted", weighted = FALSE),
    list(label = "crude, unweighted", formula = depression ~ nonhdl_q, weighted = FALSE),
    list(label = "Model 1, unweighted", formula = MODEL_1, weighted = FALSE),
    list(label = "Model 3 without HDL (the Methods' list)", formula = update(MODEL_3, . ~ . - hdl)),
    list(label = "fasting subsample and weight (Methods' 'fasting')", weight = "WTSAF2YR"),
    list(label = "weighted quartile cutpoints", derive = function(data, constants) {
      data$nonhdl_q <- quartile(data$nonhdl, constants$cutpoints_weighted)
      data
    })
  ),
  # 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 computes non-HDL cholesterol from fasting subjects' lipids, but Figure 1 keeps 34,585 adults with a PHQ-9 and non-HDL cholesterol (34,459 here with no fasting restriction), about twice the fasting subsample (13,393 of the 29,467 analyzed here have a fasting weight)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Figure 1 excludes 5,575 under 18, though it goes from 70,190 participants to 42,143 adults, 28,047 fewer; the abstract gives the 42,143 adults as the study's participants, of whom 28,041 were analyzed."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract reports a positive association in participants of normal BMI with OR 0.93 (0.66 to 1.32), where Table 3 gives 1.39 (0.92 to 2.11), and 1.29 (1.01 to 1.66) in those without hypertension, where Table 3 gives 1.30 (1.01 to 1.68).")
  ),
  choices = list(
    list(choice = "weights", decision = "examination weight (WTMEC2YR) over 7, with strata and PSUs; in 2021-2023 the phlebotomy weight (WTPH2YR from TCHOL_L)",
         reason = "the paper built 14-year weights as NHANES recommends, and the PHQ-9 and lipids come from the examination; weighted fits reproduce the crude and Model 1 estimates (1.30 and 1.40), unweighted ones do not (1.25 and 1.32). NCHS directs the phlebotomy weight for blood analytes in 2021-2023"),
    list(choice = "fasting", decision = "no fasting restriction: everyone with total and HDL cholesterol",
         reason = "the Methods say non-HDL cholesterol was computed from fasting subjects' lipids, but Figure 1 keeps 34,585 of 42,143 adults with a PHQ-9 and non-HDL cholesterol (ours 34,459), about twice the morning fasting subsample; the text's version (fasting subsample, fasting weight) is a variant"),
    list(choice = "quartiles", decision = "unweighted quartiles of non-HDL cholesterol (mmol/L) in the analytic sample, intervals closed on the right (2.84, 3.52, 4.27)",
         reason = "unreported; they reproduce Table 2's crude and Model 1 Q4 estimates to the second decimal (Q2 and Q3 within 0.03); weighted quartiles (2.87, 3.54, 4.29) put the crude Q4 at 1.27 (variant)"),
    list(choice = "HDL in Model 3", decision = "included (continuous, mmol/L)",
         reason = "Table 2's note lists HDL in Model 3 and the Methods do not; the published estimate is closer with it (variant without it)"),
    list(choice = "study population", decision = "adults 18 and older with a PHQ-9 and non-HDL cholesterol, complete case for Model 3's covariates; 18- and 19-year-olds drop out because education and the medical conditions are asked from 20",
         reason = "the paper's Figure 1 (42,143 adults, as ours). We keep 29,467 (2,564 depressed) against the paper's 28,041 (2,419): the paper excluded 6,544 for missing covariates against our 4,992, and no coding we tried accounts for the gap (requiring a measured blood pressure drops 636, counting borderline diabetes as missing 411)"),
    list(choice = "PHQ-9", decision = "sum of the nine items, missing if any is missing, refused, or don't know; depression at 10 or more",
         reason = "the paper's cutoff and its exclusion of incomplete responses; 36,259 adults have a PHQ-9 against the paper's 36,397 (allowing one missing item would give 36,401)"),
    list(choice = "age, income to poverty ratio, BMI", decision = "continuous",
         reason = "Table 1 reports them as means (weighted 47.3, 3.04, 29.1 against 47.26, 3.06, 29.08); BMI or income in groups changes Model 3 by under 0.02"),
    list(choice = "race, marital status, education", decision = "non-Hispanic White, Mexican American, non-Hispanic Black, other; married or living with a partner, never married, widowed, divorced or separated (DMDMARTL); less than high school, high school or GED, some college, college graduate",
         reason = "the paper's groups; weighted shares within 0.3 points of Table 1's"),
    list(choice = "smoking and drinking", decision = "never (fewer than 100 cigarettes), former, current (SMQ020, SMQ040); never (fewer than 12 drinks in life; never a drink from 2017-2018), former (none in the past 12 months), current",
         reason = "unstated beyond never, former, current; smoking matches Table 1 (54.7, 25.2, 20.2% against 54.8, 25.0, 20.2%), drinking less closely (10.2, 15.4, 74.4% against 10.4, 13.2, 76.5%)"),
    list(choice = "CVD", decision = "ever told of heart failure, angina, heart attack, or stroke (coronary heart disease not included)",
         reason = "the paper's list; 7.36% weighted against Table 1's 7.12%"),
    list(choice = "diabetes", decision = "told by a doctor, insulin or diabetes pills, HbA1c 6.5% or more, fasting glucose 126 mg/dL or more, or two-hour OGTT glucose 200 mg/dL or more (2005-2016)",
         reason = "unstated; this matches Table 1 (13.97% weighted against 13.89%), while told only gives 9.42% and the same without the OGTT 13.07%"),
    list(choice = "hypertension", decision = "told by a doctor, taking medicine for it, or mean blood pressure 140/90 mmHg or more",
         reason = "unstated; this matches Table 1 (37.62% weighted against 37.53%), while told only gives 31.80%"),
    list(choice = "harmonized marital status and diabetes", decision = "DMDMARTZ's three marital groups; diabetes without the OGTT",
         reason = "2021-2023 releases only three marital groups and gave no OGTT")
  )
)