# Serum estradiol and PHQ-9 score in women older than 20 who had both ovaries removed, NHANES 2013-2016.
# Chen et al. (2020), Ann Gen Psychiatry, doi:10.1186/s12991-020-00315-1.
# Headline: beta 0.014 (95% CI 0.001-0.028) PHQ-9 points per 1 pg/mL of estradiol, multivariate model
# (Table 2), P = 0.040.
#
# The paper's numbers show how it was run. Women with both ovaries removed (RHQ305), estradiol and
# testosterone (TST_H, TST_I), and a PHQ-9 score are 548, and with the MEC examination weight Table 1
# is reproduced exactly: the PHQ-9 mean (4.522 +- 0.316) as SAS's SUM() gives it (the answered
# items, refused and don't know left out; 544 answer all nine), estradiol 16.175 +- 1.981 pg/mL,
# testosterone 16.507 +- 1.328 ng/dL, and every count and weighted share. "History of estrogen use"
# is RHQ554 (ever took estrogen-only pills: 267 yes, 18 no) and its duration RHQ560Q, read as years
# whatever its unit (RHQ560U: 43 of 261 answered in months): that gives Table 1's 10.697 +- 1.108
# years and Table 2's univariate -0.104 (-0.148, -0.060). The multivariate model adjusts for that
# duration, which only women who took estrogen-only pills have, so complete-case analysis fitted it
# on those 261 (the paper reports 548): weighted, with strata and PSUs and t on the design's degrees
# of freedom, that reproduces Table 2's multivariate estimates and intervals to the digit (estradiol
# 0.0143, 0.0008-0.0278, P = 0.038 against 0.040; testosterone -0.033, -0.048 to -0.018; duration
# -0.110). The chosen version keeps that sample and converts months to years, as "duration (years)"
# means.

# PHQ-9 total as the paper's numbers show it was scored: the sum of the items answered 0-3 (refused
# and don't know left out), missing only when no item was answered.
row189_phq9 <- function(cycle) {
  d <- component("DPQ", cycle)
  items <- sapply(sprintf("DPQ%03d", seq(10, 90, 10)), function(item) ifelse(d[[item]] %in% 0:3, d[[item]], NA))
  data.frame(SEQN = d$SEQN, phq9 = ifelse(rowSums(!is.na(items)) > 0, rowSums(items, na.rm = TRUE), NA),
             phq9_all_items = rowSums(items))
}

# Taking estrogen-only pills (RHQ554) and for how long (RHQ560Q in RHQ560U's months or years; 77 and
# 99 are refused and don't know). The questions follow ever having used female hormones (RHQ540) and
# taking them as pills; women who never took estrogen-only pills get 0 years in duration_all, which
# a variant uses. 2021-2023 asks none of these.
row189_estrogen <- function(cycle, rhq) {
  if (!has(rhq, "RHQ554")) {
    return(data.frame(SEQN = rhq$SEQN, estrogen_pills = NA_real_, duration = NA_real_, duration_as_coded = NA_real_, duration_all = NA_real_))
  }
  q <- ifelse(rhq$RHQ560Q %in% 1:59, rhq$RHQ560Q, NA)
  years <- ifelse(rhq$RHQ560U %in% 1, q / 12, ifelse(rhq$RHQ560U %in% 2, q, NA))
  never <- rhq$RHQ540 %in% 2 | rhq$RHQ554 %in% 2 | (rhq$RHQ540 %in% 1 & is.na(rhq$RHQ554))
  data.frame(SEQN = rhq$SEQN, estrogen_pills = ifelse(rhq$RHQ554 %in% 1, 1, ifelse(rhq$RHQ554 %in% 2, 0, NA)),
             duration = years, duration_as_coded = q, duration_all = ifelse(never & is.na(years), 0, years))
}

row189_build <- function(cycle) {
  demo <- demographics(cycle)
  rhq <- component("RHQ", cycle)
  mcq <- component("MCQ", cycle)
  d <- merge_all(demo, row189_phq9(cycle), component("TST", cycle, c("LBXEST", "LBXTST")), rhq[, c("SEQN", "RHQ305")],
                 row189_estrogen(cycle, rhq), medical_conditions(cycle), component("DIQ", cycle, c("DIQ010", "DIQ050", "DIQ070")))
  m <- function(v) if (has(mcq, v)) mcq[[v]][match(d$SEQN, mcq$SEQN)] else rep(NA, nrow(d))
  d$estradiol <- d$LBXEST
  d$testosterone <- d$LBXTST
  # Below high school (DMDEDUC2 1-2) against high school or above (3-5): Table 1's 110 and 438.
  d$education <- factor(ifelse(d$education %in% 1:2, "below high school", ifelse(d$education %in% 3:5, "high school or above", NA)),
                        levels = c("below high school", "high school or above"))
  yes_no <- function(x) factor(ifelse(x %in% TRUE, "yes", "no"), levels = c("no", "yes"))
  # Ever told of heart failure, coronary heart disease, angina, heart attack, or stroke; anyone else
  # is "no", as Table 1 classifies all 548 (109 yes).
  d$cvd <- yes_no(d$heart_failure %in% 1 | d$chd %in% 1 | d$angina %in% 1 | d$heart_attack %in% 1 | d$stroke %in% 1)
  d$arthritis <- yes_no(d$arthritis %in% 1)
  # Chronic respiratory tract disease as Table 1 counts it (117, 23.60%): still has asthma (MCQ035),
  # emphysema (MCQ160G), or still has chronic bronchitis (MCQ170K). 2021-2023 asks emphysema and
  # chronic bronchitis only within one question, ever told of COPD, emphysema, or chronic bronchitis
  # (MCQ160P); respiratory_h, which every cycle can build, is still having asthma or ever any of
  # those (MCQ160G, MCQ160K, MCQ160O through 2018, the library's copd).
  d$respiratory <- yes_no(m("MCQ035") %in% 1 | (if (has(mcq, "MCQ160G")) m("MCQ160G") %in% 1 | m("MCQ170K") %in% 1 else d$copd %in% 1))
  d$respiratory_h <- yes_no(m("MCQ035") %in% 1 | d$copd %in% 1)
  # Diabetes (in the univariate analysis only) as Table 1 counts it (137, 21.18%): told, insulin, or pills.
  d$diabetes <- yes_no(d$DIQ010 %in% 1 | d$DIQ050 %in% 1 | d$DIQ070 %in% 1)
  d$in_population <- d$sex == 2 & d$age > 20 & d$RHQ305 %in% 1 & !is.na(d$estradiol) & !is.na(d$testosterone) & !is.na(d$phq9)
  d
}

ROW189_MODEL <- phq9 ~ estradiol + testosterone + education + duration + cvd + respiratory + arthritis

association <- list(
  id = "row189", row = 189, doi = "10.1186/s12991-020-00315-1",
  cycles = c("2013-2014", "2015-2016"),
  weight = "WTMEC2YR", blood_file = "TST",
  family = "linear", term = "estradiol",
  published = list(measure = "beta", estimate = 0.014, low = 0.001, high = 0.028, n = 548, events = NULL,
                   contrast = "PHQ-9 points per 1 pg/mL of serum estradiol"),
  left_out = c("duration of estrogen-only pill use (RHQ560Q/U): 2021-2023's reproductive health file asks no hormone-use questions; leaving it out also lifts the restriction it put on the paper's sample, which complete-case analysis limited to the 261 women who took estrogen-only pills and said for how long, so the harmonized version takes every woman with both ovaries removed (548 in 2013-2016)",
               "chronic respiratory disease as the paper counted it (still having asthma, emphysema, or still having chronic bronchitis, MCQ035, MCQ160G, MCQ170K): 2021-2023 asks emphysema and chronic bronchitis only within one question on ever having COPD, emphysema, or chronic bronchitis (MCQ160P), so the harmonized version counts still having asthma or ever any of those, in every cycle"),
  build = row189_build,
  # The paper's computation read RHQ560Q as years whatever its unit (43 answers given in months);
  # that reproduces Table 2 exactly, so the paper version follows it. The harmonized version has
  # no duration at all (2021-2023 asks no hormone-use questions).
  formula = update(ROW189_MODEL, . ~ . - duration + duration_as_coded),
  formula_harmonized = phq9 ~ estradiol + testosterone + education + cvd + respiratory_h + arthritis,
  variants = list(
    list(label = "duration with months converted to years", formula = ROW189_MODEL),
    list(label = "testosterone, duration as coded (published -0.033, -0.048 to -0.018)", term = "testosterone",
         formula = update(ROW189_MODEL, . ~ . - duration + duration_as_coded)),
    list(label = "unweighted", weighted = FALSE),
    list(label = "univariate, all 548 (published -0.003, -0.011 to 0.004)", formula = phq9 ~ estradiol, sample = "own"),
    list(label = "users of estrogen-only pills, no duration (Table 3: 0.014, 0.001-0.027)",
         formula = phq9 ~ estradiol + testosterone + education + cvd + respiratory + arthritis, sample = "own",
         derive = function(data, constants) { data$in_population <- data$in_population & data$estrogen_pills %in% 1; data }),
    list(label = "all women, 0 years for never taking estrogen-only pills", sample = "own",
         formula = update(ROW189_MODEL, . ~ . - duration + duration_all)),
    list(label = "PHQ-9 only when all nine items are answered", sample = "own",
         formula = update(ROW189_MODEL, phq9_all_items ~ .))
  ),
  # 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 = "sample", affects_headline = TRUE, followed = TRUE,
         detail = "The paper reports n = 548 and states the result for ovariectomized women, but its multivariate model adjusts for the duration of estrogen-only pill use, which only the 261 women who took them and gave a duration have, and only that sample reproduces Table 2's multivariate estimates and intervals to the digit (estradiol 0.0143, 0.0008 to 0.0278). With never-users given 0 years (528 women), the estimate is 0.006 (-0.006 to 0.018)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "Table 2 gives the duration of estrogen use in years, but Table 1's mean (10.697 years) and Table 2's estimates for it (univariate -0.104, multivariate -0.110) are reproduced only with RHQ560Q read as years whatever its unit, so the 43 durations given in months count as years. Converted to years, they give 0.0150 for the headline, against 0.0143.")
  ),
  choices = list(
    list(choice = "the multivariate model's sample", decision = "complete cases: the women with a duration of estrogen-only pill use (261 of 548)",
         reason = "the model adjusts for that duration, which only women who took estrogen-only pills have; Table 1 gives its mean over them alone (10.697 years), and Table 3's models by estrogen use leave it out. The paper reports 548 and states the result for ovariectomized women in general, but only this sample reproduces Table 2 (to the digit with the duration as coded); 0 years for never-users keeps 528 women and gives 0.006 (variant)"),
    list(choice = "duration of estrogen use", decision = "RHQ560Q read as years whatever its unit, as the paper's computation read it",
         reason = "the paper's numbers (Table 1's 10.697 +- 1.108 years, the univariate -0.104 and the multivariate -0.110) are RHQ560Q read as years whatever its unit, so 43 answers in months count as years; the published estimate comes from that coding, and months divided by 12 is a variant"),
    list(choice = "history of estrogen use", decision = "ever took estrogen-only pills (RHQ554)", reason = "Table 1's 267 yes and 18 no (55.42% and 4.92%), reproduced exactly; it enters only a variant"),
    list(choice = "PHQ-9 score", decision = "sum of the items answered 0-3, missing only if none was answered",
         reason = "unstated; this gives the paper's n of 548 and Table 1's 4.522 +- 0.316 exactly (all nine items required: 544 women, 4.491)"),
    list(choice = "weights", decision = "examination weight (WTMEC2YR over 2, with strata and PSUs; WTPH2YR from TST_L in 2021-2023)",
         reason = "'discharge weights' (sic); the MEC weight reproduces Table 1's shares and means and Table 2's estimates and CIs; unweighted fits do not (variant)"),
    list(choice = "age", decision = "older than 20, as written", reason = "RHQ305 is asked from 20, and every woman in 2013-2016 with both ovaries removed was 31 or older, so it changes nothing there"),
    list(choice = "CVD, chronic respiratory disease, arthritis", decision = "ever told of heart failure, coronary heart disease, angina, heart attack or stroke; still having asthma, emphysema, or still having chronic bronchitis; ever told of arthritis; everyone else no",
         reason = "unstated; these reproduce Table 1's counts and weighted shares exactly (109, 117, 342 yes; COPD (MCQ160O) is not part of the respiratory count)"),
    list(choice = "education", decision = "below high school (DMDEDUC2 1-2) against high school or above", reason = "Table 1's 110 and 438"),
    list(choice = "estradiol below the detection limit", decision = "NCHS's fill value (the limit over the square root of 2), as released",
         reason = "unstated; 135 of the 548 are below the limit. 2021-2023 measures estradiol by the same ID-LC-MS/MS method family in a wider steroid panel with a lower limit (1.72 pg/mL against 2.994 in 2013-2016), so fill values drop from 2.12 to 1.22 pg/mL: a measurement change recorded, not corrected")
  )
)