# Workday sleep duration and systolic blood pressure in people with BMI under 25, NHANES 2015-2018.
# Su et al. (2022), Sci Rep, doi:10.1038/s41598-022-05124-y.
# Headline: beta 3.58 (95% CI 1.60-5.56) mmHg for workday sleep under 6 hours against 6 to under
# 8 hours, Model II (Table 3).

association <- list(
  id = "row068", row = 68, doi = "10.1038/s41598-022-05124-y",
  cycles = c("2015-2016", "2017-2018"),
  weight = "WTMEC2YR", blood_file = "BIOPRO",
  family = "linear", term = "sleep3short",
  published = list(measure = "beta", estimate = 3.58, low = 1.60, high = 5.56, n = 2887, events = NA,
                   contrast = "workday sleep under 6 hours vs 6 to under 8 hours (systolic blood pressure, mmHg)"),
  left_out = c("snort, gasp, or stop breathing while asleep (SLQ040): not asked in 2021-2023"),
  build = function(cycle) {
    demo <- demographics(cycle)
    slq <- component("SLQ", cycle)
    s <- data.frame(SEQN = slq$SEQN, sleep_hours = slq$SLD012,
                    snort_item = if (has(slq, "SLQ040")) slq$SLQ040 else NA)
    # Taking prescribed medicine for high blood pressure: BPQ050A in 2015-2018, BPQ150 in
    # 2021-2023; both are asked only of those told they had high blood pressure.
    bpq <- component("BPQ", cycle)
    b <- data.frame(SEQN = bpq$SEQN, BPQ020 = bpq$BPQ020, bp_medication = if (has(bpq, "BPQ050A")) bpq$BPQ050A else bpq$BPQ150)
    # Drinking in the past 12 months: how often (ALQ120Q, with its unit) in 2015-2016, ALQ121
    # from 2017 on; those who never drank (or had under 12 drinks in life, in 2015-2016) skip it.
    alq <- component("ALQ", cycle)
    drink <- if (has(alq, "ALQ121")) ifelse(alq$ALQ121 %in% 0, "no", ifelse(alq$ALQ121 %in% 1:10, "drinking", NA)) else
      ifelse(alq$ALQ120Q %in% 0, "no", ifelse(alq$ALQ120Q %in% 1:365, "drinking", NA))
    d <- merge_all(demo, s, blood_pressure(cycle), b, body_measures(cycle), data.frame(SEQN = alq$SEQN, drink = drink),
                   component("BIOPRO", cycle, c("LBDSALSI", "LBDSCRSI", "LBXSASSI")), component("CBC", cycle, "LBXHGB"),
                   component("TCHOL", cycle, "LBDTCSI"), component("HDL", cycle, "LBDHDDSI"),
                   component("DIQ", cycle, "DIQ010"), component("SMQ", cycle, "SMQ040"))
    d$sleep3 <- factor(ifelse(d$sleep_hours < 6, "short", ifelse(d$sleep_hours < 8, "mid", "long")), levels = c("mid", "short", "long"))
    d$sex <- factor(d$sex)
    d$race4 <- factor(ifelse(d$race %in% 1, "mexican", ifelse(d$race %in% 3, "white", ifelse(d$race %in% 4, "black", "other"))),
                      levels = c("mexican", "white", "black", "other"))
    d$alcohol <- with_unclear(d$drink, c("no", "drinking"))
    d$diabetes <- with_unclear(ifelse(d$DIQ010 %in% 1, "yes", ifelse(d$DIQ010 %in% 2, "no", ifelse(d$DIQ010 %in% 3, "borderline", NA))),
                               c("yes", "no", "borderline"))
    d$hypertension <- with_unclear(ifelse(d$BPQ020 %in% 1, "yes", ifelse(d$BPQ020 %in% 2, "no", NA)), c("yes", "no"))
    d$snort <- with_unclear(ifelse(d$snort_item %in% 1:3, "yes", ifelse(d$snort_item %in% 0, "no", NA)), c("no", "yes"))
    # The paper's Table 1 shows told hypertension and snorting recorded only for 2017-2018: its
    # "not recorded" shares are 2015-2016's weighted share plus a few nonresponses. These copies
    # reproduce that, for the paper version; the harmonized version records them in every cycle.
    lost <- cycle == "2015-2016"
    d$hypertension_as_published <- if (lost) factor(rep("unclear", nrow(d)), levels = levels(d$hypertension)) else d$hypertension
    d$snort_as_published <- if (lost) factor(rep("unclear", nrow(d)), levels = levels(d$snort)) else d$snort
    d$smoke <- with_unclear(ifelse(d$SMQ040 %in% 1:2, "smoking", ifelse(d$SMQ040 %in% 3, "not smoking", NA)), c("smoking", "not smoking"))
    # Missing laboratory covariates: an indicator for each, with the missing value set to 0 (any
    # constant gives the same fit once the indicator is in the model).
    lab <- function(x) ifelse(is.na(x), 0, x)
    flag <- function(x) as.integer(is.na(x))
    d$albumin <- lab(d$LBDSALSI); d$albumin_missing <- flag(d$LBDSALSI)
    d$creatinine <- lab(d$LBDSCRSI); d$creatinine_missing <- flag(d$LBDSCRSI)
    d$ast <- lab(d$LBXSASSI); d$ast_missing <- flag(d$LBXSASSI)
    d$hemoglobin <- lab(d$LBXHGB); d$hemoglobin_missing <- flag(d$LBXHGB)
    d$tc <- lab(d$LBDTCSI); d$tc_missing <- flag(d$LBDTCSI)
    d$hdl <- lab(d$LBDHDDSI); d$hdl_missing <- flag(d$LBDHDDSI)
    # Everyone with workday sleep hours (asked from age 16) and a blood pressure reading, not
    # taking blood pressure medicine, with BMI under 25.
    d$in_population <- !is.na(d$sleep_hours) & !is.na(d$sbp) & !(d$bp_medication %in% 1) & !is.na(d$bmi) & d$bmi < 25
    d
  },
  variants = list(
    list(label = "told hypertension and snorting recorded in both cycles, as the text defines them", formula =
           sbp ~ sleep3 + sex + age + race4 + alcohol + albumin + albumin_missing + creatinine + creatinine_missing + hemoglobin +
           hemoglobin_missing + diabetes + hypertension + snort + smoke + tc + tc_missing + bmi + ast + ast_missing + hdl + hdl_missing),
    list(label = "unweighted", weighted = FALSE),
    list(label = "crude, weighted (published 6.15, 3.88-8.42)", formula = sbp ~ sleep3),
    list(label = "crude, unweighted", formula = sbp ~ sleep3, weighted = FALSE),
    list(label = "Model I, weighted (published 4.17, 2.19-6.15)", formula = sbp ~ sleep3 + sex + age + race4),
    list(label = "Model I, unweighted", formula = sbp ~ sleep3 + sex + age + race4, weighted = FALSE)
  ),
  # 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 = "Model II adjusts for told hypertension and snort or stop breathing, but Table 1's 'not recorded' shares for them (50.04% and 51.37%) equal 2015-2016's weighted share plus a few nonresponses, so the analysis had neither item for 2015-2016. With both unrecorded there, Model II gives 3.581 (1.70 for 8 hours or more), as published, and the paper version codes them so; recorded in both cycles, as the text defines them and the harmonized version records them, it gives 3.678."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's snort or stop breathing rows swap yes and no in the columns by sex: 6.72% of all participants are yes (as the Results say) and 41.91% no, which are the averages of the men's and women's 'No' shares (8.56%, 5.22%) and 'Yes' shares (39.08%, 44.22%).")
  ),
  choices = list(
    list(choice = "population", decision = "everyone with workday sleep hours (SLD012, asked from age 16), a mean systolic reading, not taking prescribed blood pressure medicine (BPQ050A), and BMI under 25; no age limit",
         reason = "the paper states no age limit; these steps reproduce its flow exactly (6,818, 1,055, 2,944, and 5,521 excluded; 2,887 left; 1,378 men, 1,509 women)"),
    list(choice = "weights", decision = "examination weight, WTMEC2YR, each cycle's halved; in 2021-2023 the phlebotomy weight (WTPH2YR from BIOPRO_L), because the model holds blood analytes",
         reason = "the paper names no weight; Table 1's weighted means and percentages (age 38.54, SBP 115.33, BMI 21.97, sleep groups 40.53/6.26/53.21) are reproduced exactly with it, and so are the crude and Model I estimates and CIs (as weighted least squares, which the paper's CIs are); unweighted fits do not reproduce them. NCHS directs the phlebotomy weight for blood analytes in 2021-2023; it leaves out those with no blood sample (about 6% there), whom the paper's missing indicators would keep"),
    list(choice = "sleep groups", decision = "under 6 hours, 6 to under 8 (reference), 8 or more",
         reason = "the paper's labels (< 6 h, 6-8 h, >= 8 h); they reproduce Table 1's weighted shares exactly"),
    list(choice = "systolic blood pressure", decision = "mean of all available auscultatory readings (BPXSY1-4); 2021-2023 measured only with the oscillometric device (BPXO, three readings)",
         reason = "the paper's definition; the library reads BPX in 2015-2018 and BPXO in 2021-2023"),
    list(choice = "race", decision = "Mexican American, non-Hispanic White, non-Hispanic Black, and other (other Hispanic and other or multiracial, RIDRETH1 2 and 5)",
         reason = "reproduces Table 1's weighted shares and Table 4's subgroup counts (302, 968, 550, 1,067)"),
    list(choice = "alcohol, smoking, diabetes", decision = "drinking in the past 12 months (ALQ120Q above 0 in 2015-2016, ALQ121 1-10 after) or not (0), else not recorded; smoking now (SMQ040 every day or some days), not at all, else not recorded (never smokers and under 18s skip it); DIQ010 yes, no, borderline, else not recorded",
         reason = "the paper's definitions; each reproduces Table 1's weighted percentages exactly"),
    list(choice = "snort or stop breathing", decision = "SLQ040 rarely, occasionally, or frequently as yes, never as no, don't know as not recorded",
         reason = "unstated; this split reproduces Table 1's percentages exactly (for 2017-2018, the only cycle the paper's data recorded)"),
    list(choice = "told hypertension and snorting in 2015-2016", decision = "unrecorded in 2015-2016 in the paper version, as the paper's data had them; recorded from BPQ020 in every cycle in the harmonized version (2021-2023 has no snorting item)",
         reason = "Table 1's 'not recorded' shares (50.04% and 51.37%) equal 2015-2016's weighted share plus a few nonresponses: the paper's data lost both items for 2015-2016, a data-handling failure of its own files, which the plan follows in the paper version only. That version reproduces Model II to the second decimal (3.58, 1.60-5.56; 1.70 for 8 hours or more); recorded in both cycles (variant) it gives 3.68"),
    list(choice = "missing laboratory covariates", decision = "albumin, creatinine, AST, hemoglobin, total cholesterol, and HDL each with a missing indicator and the missing value set to 0",
         reason = "the paper's missing-indicator method; the fill value does not change the fit")
  ),
  formula = sbp ~ sleep3 + sex + age + race4 + alcohol + albumin + albumin_missing + creatinine + creatinine_missing + hemoglobin +
    hemoglobin_missing + diabetes + hypertension_as_published + snort_as_published + smoke + tc + tc_missing + bmi + ast +
    ast_missing + hdl + hdl_missing,
  formula_harmonized = sbp ~ sleep3 + sex + age + race4 + alcohol + albumin + albumin_missing + creatinine + creatinine_missing +
    hemoglobin + hemoglobin_missing + diabetes + hypertension + smoke + tc + tc_missing + bmi + ast + ast_missing + hdl + hdl_missing
)