# Diastolic blood pressure and depression (PHQ-9 total of 10 or more), NHANES 2005-2018.
# Zhang, Xu, and Xu (2024), Front Psychiatry, doi:10.3389/fpsyt.2024.1433990.
# Headline: OR 1.06 (95% CI 1.00-1.12) per 10 mmHg of diastolic pressure, Model 4 (Table 2).

# Mean pressure by the NHANES protocol the paper quotes: one reading is the average; with more
# than one, the first is left out and the rest averaged. 2005-2018 measured by auscultation (BPX,
# up to four readings, BPX_J in 2017-2018); 2021-2023 by an oscillometric device (BPXO, three
# readings). A diastolic reading of 0 is not a reading: NHANES leaves it out of the average when
# another reading is above zero, and someone whose diastolic readings are all 0 has none.
protocol_bp <- function(cycle) {
  oscillometric <- cycle %in% c("2017-2020", REPLICATION_CYCLE)
  d <- component(if (oscillometric) "BPXO" else "BPX", cycle)
  average <- function(prefix, diastolic) {
    m <- as.matrix(d[, intersect(paste0(prefix, 1:4), names(d)), drop = FALSE])
    if (diastolic) m[m %in% 0] <- NA
    apply(m, 1, function(r) {
      r <- r[!is.na(r)]
      if (length(r) == 0) NA_real_ else if (length(r) == 1) r else mean(r[-1])
    })
  }
  data.frame(SEQN = d$SEQN, sbp = average(if (oscillometric) "BPXOSY" else "BPXSY", FALSE),
             dbp = average(if (oscillometric) "BPXODI" else "BPXDI", TRUE))
}

# Drinking status (1 never, 2 former, 3 current) as the library's alcohol() defines it: never is
# fewer than 12 drinks in life (2005-2016) or never a drink (2017 on), former is no drinking in
# the past 12 months, current is any. alcohol() itself stops on 2005-2010, whose ALQ files have
# no ALQ151, so the status is built here.
drinking <- function(cycle) {
  d <- component("ALQ", cycle)
  if (has(d, "ALQ121")) {
    ever <- yes(d$ALQ111)
    past_year <- ifelse(d$ALQ121 %in% 1:10, 1, ifelse(d$ALQ121 %in% 0, 0, NA))
  } else {
    ever <- ifelse(yes(d$ALQ101) %in% 1 | yes(d$ALQ110) %in% 1, 1,
                   ifelse(yes(d$ALQ110) %in% 0 | (yes(d$ALQ101) %in% 0 & is.na(d$ALQ110)), 0, NA))
    q <- ifelse(d$ALQ120Q %in% 0:365, d$ALQ120Q, NA)
    per_year <- ifelse(d$ALQ120U %in% 1, 52, ifelse(d$ALQ120U %in% 2, 12, ifelse(d$ALQ120U %in% 3, 1, NA)))
    days <- ifelse(q %in% 0, 0, q * per_year)
    past_year <- ifelse(days > 0, 1, ifelse(days %in% 0, 0, NA))
    past_year[ever %in% 0] <- 0
  }
  data.frame(SEQN = d$SEQN, alcohol3 = ifelse(ever %in% 0, 1, ifelse(ever %in% 1 & past_year %in% 0, 2, ifelse(past_year %in% 1, 3, NA))))
}

# Multum Lexicon categories (RXQ_DRUG) of the prescription medicines taken in the past 30 days:
# antidepressants, and the classes of blood-pressure-lowering drugs (ACE inhibitors, peripherally
# and centrally acting antiadrenergic agents, beta blockers, calcium channel blockers, diuretics,
# vasodilators, antihypertensive combinations, angiotensin II inhibitors, aldosterone receptor
# antagonists, renin inhibitors).
ANTIDEPRESSANTS <- 249
ANTIHYPERTENSIVES <- c(42, 43, 44, 47, 48, 49, 53, 55, 56, 340, 342)

# The paper's models (Table 2 and its note). Age and BMI enter in Table 1's groups.
MODEL_1 <- depression ~ dbp10 + age60 + sex + education2 + race
MODEL_2 <- update(MODEL_1, . ~ . + marital + pir3 + smoking + alcohol + bmi25 + activity + energy)
MODEL_3 <- update(MODEL_2, . ~ . + arthritis + thyroid + cancer + diabetes + liver)
MODEL_4 <- update(MODEL_3, . ~ . + chd + heart_failure + ckd + stroke + hyperlipidemia + antihypertensive + antidepressant)

# The model variables that use constants fixed on the paper's cycles.
derive_row050 <- function(data, constants) {
  split <- function(x, at) factor(ifelse(x >= at, "high", "low"), levels = c("low", "high"))
  data$energy <- split(data$KCAL, constants$kcal_median)
  data$energy_two_day <- split(data$KCAL2, constants$kcal2_median)
  # Activity in Table 1's groups: low or high at the median of those who reported any, and
  # unknown for those who reported none or were not asked (2005-2006 had no GPAQ).
  level <- function(x, at) factor(ifelse(is.na(x) | x %in% 0, "unknown", ifelse(x >= at, "high", "low")),
                                  levels = c("low", "high", "unknown"))
  data$activity <- level(data$met, constants$met_median)
  data$activity_leisure <- level(data$met_leisure, constants$leisure_median)
  data$activity_missing <- factor(ifelse(is.na(data$met), "unknown", ifelse(data$met >= constants$met_known_median, "high", "low")),
                                  levels = c("low", "high", "unknown"))
  data
}

association <- list(
  id = "row050", row = 50, doi = "10.3389/fpsyt.2024.1433990",
  cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  # CKD (serum creatinine), diabetes (HbA1c), and hyperlipidemia (lipids) use blood analytes.
  weight = "WTMEC2YR", blood_file = "BIOPRO",
  family = "logistic", term = "dbp10",
  published = list(measure = "OR", estimate = 1.06, low = 1.00, high = 1.12, n = 26581, events = 2261,
                   contrast = "per 10 mmHg of diastolic blood pressure"),
  left_out = c(
    "antidepressant use (Multum category 249 in the prescription file): 2021-2023 releases no drug names, and no question asks about antidepressants",
    "antihypertensive drug use from the prescription file (Multum blood-pressure-lowering classes): 2021-2023 releases no drug names; replaced by self-reported current treatment for high blood pressure (BPQ050A; BPQ150 in 2021-2023)",
    "physical activity over work, transport, and leisure (GPAQ MET-minutes): 2021-2023 asks only about leisure time; replaced by leisure-time moderate (4 METs) and vigorous (8 METs) activity in the paper's groups, split at the paper cycles' median"
  ),
  build = function(cycle) {
    demo <- demographics(cycle)
    bpq <- component("BPQ", cycle)
    paq <- component("PAQ", cycle)
    # 2005-2006 asked another activity questionnaire, without the GPAQ items (leisure_activity
    # stops on it, met_minutes gives NA).
    leisure <- if (has(paq, "PAQ650") || has(paq, "PAD790Q")) leisure_activity(cycle)[, c("SEQN", "leisure_moderate", "leisure_vigorous")] else
      data.frame(SEQN = paq$SEQN, leisure_moderate = NA_real_, leisure_vigorous = NA_real_)
    two_day <- dietary_totals(cycle, "KCAL", days = 2)
    ckd <- ckd_status(cycle, demo)
    d <- merge_all(demo, phq9(cycle), protocol_bp(cycle), smoking(cycle), drinking(cycle),
                   body_measures(cycle)[, c("SEQN", "bmi")], met_minutes(cycle), leisure,
                   dietary_totals(cycle, "KCAL", days = 1), medical_conditions(cycle),
                   component("MCQ", cycle, "MCQ160M"), diabetes_status(cycle), ckd,
                   urine_albumin_creatinine(cycle), hyperlipidemia_status(cycle, demo))
    d$KCAL2 <- two_day$KCAL[match(d$SEQN, two_day$SEQN)]
    if (cycle != REPLICATION_CYCLE) {
      d <- merge_all(d, medication_use(cycle, list(antihypertensive = ANTIHYPERTENSIVES, antidepressant = ANTIDEPRESSANTS)))
    } else {
      d$antihypertensive <- NA
      d$antidepressant <- NA
    }
    # Self-reported current treatment for high blood pressure. BPQ050A (2005-2018) is asked of
    # those told they had it and told to take medicine (BPQ040A); BPQ150 (2021-2023) of those told
    # they had it. Those not told, or not told to take medicine, are not treated.
    treated <- if (has(bpq, "BPQ150")) bpq$BPQ150 else bpq$BPQ050A
    not_advised <- if (has(bpq, "BPQ040A")) bpq$BPQ040A %in% 2 else rep(FALSE, nrow(bpq))
    bp_med <- ifelse(treated %in% 1, 1, ifelse(treated %in% 2 | bpq$BPQ020 %in% 2 | not_advised, 0, NA))
    d$bp_medication <- bp_med[match(d$SEQN, bpq$SEQN)]
    d$depression <- as.integer(d$phq9 >= 10)
    d$dbp10 <- d$dbp / 10
    d$age60 <- factor(ifelse(d$age >= 60, ">=60", "<60"), levels = c("<60", ">=60"))
    d$age3 <- cut(d$age, c(-Inf, 40, 60, Inf), right = FALSE, labels = c("<40", "40-59", ">=60"))
    d$sex <- factor(d$sex)
    d$race <- factor(d$race, levels = c(3, 1, 2, 4, 5))
    d$education2 <- factor(ifelse(d$education %in% 1:3, "less than college", ifelse(d$education %in% 4:5, "college or higher", NA)),
                           levels = c("less than college", "college or higher"))
    d$marital <- factor(d$marital)
    d$pir3 <- cut(d$pir, c(-Inf, 1.3, 3.5, Inf), right = FALSE, labels = c("<1.3", "1.3-3.5", ">=3.5"))
    d$smoking <- factor(d$smoking)
    d$alcohol <- factor(d$alcohol3)
    d$bmi25 <- factor(ifelse(d$bmi >= 25, ">=25", "<25"), levels = c("<25", ">=25"))
    d$met_leisure <- 4 * d$leisure_moderate + 8 * d$leisure_vigorous
    d$thyroid <- yes(d$MCQ160M)
    d$liver <- d$liver_condition
    d$ckd_both <- ifelse(is.na(d$egfr) | is.na(d$acr), NA, d$ckd)
    d$in_population <- !is.na(d$phq9) & !is.na(d$dbp)
    d
  },
  constants = function(data) {
    # Table 1's low and high groups of energy intake and of activity are nearly equal in size, so
    # each is split at its median in the analytic sample (unweighted), high at the median or above.
    raw <- setdiff(all.vars(MODEL_4), c("energy", "activity"))
    sample <- data$in_population & stats::complete.cases(data[, c(raw, "KCAL")])
    sample_two_day <- data$in_population & stats::complete.cases(data[, c(raw, "KCAL2")])
    positive <- function(x) x[!is.na(x) & x > 0]
    list(kcal_median = stats::median(data$KCAL[sample]),
         kcal2_median = stats::median(data$KCAL2[sample_two_day]),
         met_median = stats::median(positive(data$met[sample])),
         met_known_median = stats::median(data$met[sample & !is.na(data$met)]),
         leisure_median = stats::median(positive(data$met_leisure[sample])))
  },
  derive = derive_row050,
  formula = MODEL_4,
  formula_harmonized = update(MODEL_4, . ~ . - activity + activity_leisure - antihypertensive + bp_medication - antidepressant),
  variants = list(
    list(label = "unweighted (published 1.05, 1.01-1.09)", weighted = FALSE),
    list(label = "crude (published 1.05, 1.00-1.10)", formula = depression ~ dbp10),
    list(label = "Model 1 (published 1.07, 1.01-1.13)", formula = MODEL_1),
    list(label = "Model 2 (published 1.08, 1.02-1.14)", formula = MODEL_2),
    list(label = "Model 3 (published 1.07, 1.01-1.13)", formula = MODEL_3),
    list(label = "age continuous", formula = update(MODEL_4, . ~ . - age60 + age)),
    list(label = "age in three groups (<40, 40-59, 60+)", formula = update(MODEL_4, . ~ . - age60 + age3)),
    list(label = "BMI continuous", formula = update(MODEL_4, . ~ . - bmi25 + bmi)),
    list(label = "activity unknown only if missing (no activity is low)", formula = update(MODEL_4, . ~ . - activity + activity_missing)),
    list(label = "energy intake as the two-day mean", formula = update(MODEL_4, . ~ . - energy + energy_two_day), sample = "own"),
    list(label = "CKD only where eGFR and ACR are both known", formula = update(MODEL_4, . ~ . - ckd + ckd_both)),
    list(label = "pregnant women left out", derive = function(data, constants) {
      data <- derive_row050(data, constants)
      data$in_population <- data$in_population & !(data$pregnant %in% 1)
      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, evidence = "paper",
         detail = "The covariate text gives age as '<40 years, >=60 years', but Table 1 groups it as under 60 and 60 or more, and those two groups reproduce Table 2's estimates from the crude model to Model 4 (1.049, 1.075, 1.085, 1.068, 1.059 against 1.05, 1.07, 1.08, 1.07, 1.06), where three groups (under 40, 40 to 59, 60 or more) give 1.072 for Model 4."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The text splits physical activity into high and low at the median MET and excludes participants missing covariates, but Table 1 has a third activity group, 'Unknown' (6,534 of the 26,581 analyzed; 32.38% of the depressed against 19.12% of the others)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2's row for diastolic pressure of 90 or more gives Model 2 the same estimate, interval, and P as Model 1 (1.34, 1.02 to 1.76, 0.04), while every other row changes between the two models, as this row does in Table 3's unweighted and imputed analyses (1.28 to 1.25, 1.46 to 1.40).")
  ),
  choices = list(
    list(choice = "weights", decision = "examination weight WTMEC2YR over 7 (NCHS's pooled weight for seven 2-year cycles)",
         reason = "the paper names 1/7 * WTMEC2YR; weighted fits reproduce its crude to Model 4 sequence, and the unweighted fit its sensitivity analysis 1"),
    list(choice = "blood pressure", decision = "auscultatory readings (BPX, BPX_J in 2017-2018) averaged by the protocol the paper quotes: the only reading, or the mean of those after the first; a diastolic reading of 0 is no reading",
         reason = "the paper quotes the protocol and names mercury sphygmomanometers; with all-zero diastolic readings as missing, 948 with a PHQ-9 lack a pressure against the paper's 954 (825 if zeros counted)"),
    list(choice = "depression", decision = "PHQ-9 total of 10 or more, missing if any item is missing, refused, or don't know",
         reason = "the paper's cutoff; the convention for incomplete answers (36,259 with a PHQ-9 against the paper's 36,397)"),
    list(choice = "study population", decision = "everyone with a PHQ-9 and a diastolic pressure and complete covariates; no age bound, which leaves adults 20 and older because education and the medical conditions are asked from 20",
         reason = "the paper's flow has no age or pregnancy step"),
    list(choice = "age", decision = "two groups, under 60 and 60 or older",
         reason = "Table 1's groups; the Methods list 'age (<40 years, >=60 years)', a typo for these groups or a list of three; two groups reproduce Table 2's sequence (1.049, 1.075, 1.085, 1.068, 1.059 against 1.05, 1.07, 1.08, 1.07, 1.06), three groups less closely, and continuous age puts Models 1-4 0.02-0.04 too high (variants)"),
    list(choice = "BMI", decision = "two groups, under 25 and 25 or more",
         reason = "Table 1's groups; the Methods give none, but with BMI continuous Model 2 (1.065) falls below Model 1 (1.075) where the paper's rises (1.07 to 1.08), and Model 4 is 1.051 (variant)"),
    list(choice = "education", decision = "less than college (DMDEDUC2 1-3) or college or higher (4-5)",
         reason = "the paper's groups; 62.3% college or higher (weighted) against Table 1's 62.56%"),
    list(choice = "income to poverty ratio", decision = "under 1.3, 1.3 to under 3.5, 3.5 or more", reason = "Table 1"),
    list(choice = "race", decision = "RIDRETH1's five groups", reason = "the paper's five groups"),
    list(choice = "drinking", decision = "never (fewer than 12 drinks in life; never a drink from 2017), former (none in the past 12 months), current",
         reason = "unstated; the library's definition, built in this file because alcohol() stops on 2005-2010; Table 1 has fewer former drinkers (13.1% against our 15.3%)"),
    list(choice = "smoking", decision = "never (fewer than 100 cigarettes), former, current (SMQ020, SMQ040)", reason = "unstated; convention"),
    list(choice = "energy intake", decision = "first-day recall split at the analytic sample's unweighted median, high at the median or above",
         reason = "unstated; Table 1's groups are near equal in size (13,284 and 13,297), so the median is the sample's; the two-day mean would drop another 3,106 people (variant)"),
    list(choice = "physical activity", decision = "GPAQ MET-minutes a week (8 METs vigorous work and leisure, 4 moderate and walking or cycling), low or high at the median of those who reported any, unknown for those who reported none or were not asked (2005-2006 had no GPAQ)",
         reason = "unstated; Table 1 keeps an 'Unknown' group that is far more common among the depressed (32.4% against 19.1%), as no activity is and missingness is not; counting no activity as low changes the estimate little (variant)"),
    list(choice = "medical conditions", decision = "ever told by a doctor: arthritis (MCQ160A), thyroid problem (MCQ160M), cancer (MCQ220), liver condition (MCQ160L), CHD (MCQ160C), CHF (MCQ160B), stroke (MCQ160F)",
         reason = "unstated; ever-told prevalences match Table 1 (thyroid 11.0% against 10.8%, liver 3.5% against 3.4%)"),
    list(choice = "diabetes", decision = "told by a doctor, insulin or pills, HbA1c 6.5% or more, or fasting glucose 126 mg/dL or more",
         reason = "unstated; common definition (12.9% against Table 1's 13.6%)"),
    list(choice = "CKD", decision = "eGFR (CKD-EPI 2009) under 60 or urine ACR 30 mg/g or more, from whichever was measured",
         reason = "the paper's A2 or G3a and above; requiring both measures drops 535 people and changes the estimate little (variant)"),
    list(choice = "hyperlipidemia", decision = "NCEP ATP III thresholds or cholesterol medication", reason = "unstated; 69.8% against Table 1's 70.3%"),
    list(choice = "antihypertensive and antidepressant drugs", decision = "prescription medicines taken in the past 30 days (RXQ_RX) in Multum's blood-pressure-lowering classes (42, 43, 44, 47, 48, 49, 53, 55, 56, 340, 342) and antidepressants (249)",
         reason = "the Results speak of use in the past month, the prescription file's window; 27.9% and 13.1% against Table 1's 26.8% and 12.9%"),
    list(choice = "missing covariates", decision = "complete case", reason = "the paper's flow excludes those missing covariates"),
    list(choice = "harmonized physical activity", decision = "leisure-time moderate (4 METs) and vigorous (8 METs) activity in the paper's groups: unknown if none or missing, low or high at the paper cycles' median of those with any (960 MET-minutes a week)",
         reason = "2021-2023 asks only about leisure time; the plan's rule for activity"),
    list(choice = "harmonized antihypertensive drug", decision = "self-reported current treatment for high blood pressure (BPQ050A; BPQ150 in 2021-2023); untreated if never told of high blood pressure or not told to take medicine",
         reason = "the closest 2021-2023 measure of the same covariate, since its prescription file has no drug names"),
    list(choice = "2021-2023 weight", decision = "phlebotomy weight (WTPH2YR, from BIOPRO_L)",
         reason = "CKD (serum creatinine), diabetes (HbA1c, glucose), and hyperlipidemia (lipids) use blood analytes, for which NCHS directs the phlebotomy weight")
  )
)