# Workday sleep duration and serum albumin in adults, NHANES 2015-2018.
# Li and Guo (2022), BMC Public Health, doi:10.1186/s12889-022-13524-y.
# Headline: beta -1.00 (95% CI -1.26 to -0.74) g/L for workday sleep of 5 hours or less against
# more than 7 up to 8 hours, Model 3 (Table 2 and the abstract).
# The text says it used the examination weight, but its Table 1 and all three models of Table 2 are
# reproduced to every printed digit only with the interview weight, so that is the chosen weight
# (see choices); the examination weight is a variant.

# Continuous covariates whose missing values the paper filled with the median of its 9,973
# participants (Table 1's means and SDs are reproduced exactly only with that).
row284_imputed <- c("tp", "alt", "ast", "creatinine", "uacr", "crp", "glucose", "bmi")

row284_constants <- function(data) {
  p <- data$in_population %in% TRUE
  sapply(row284_imputed, function(v) stats::median(data[[v]][p], na.rm = TRUE), simplify = FALSE)
}

row284_derive <- function(data, constants) {
  if (is.null(constants)) constants <- row284_constants(data)
  for (v in row284_imputed) data[[paste0(v, "_filled")]] <- ifelse(is.na(data[[v]]), constants[[v]], data[[v]])
  data$creatinine_log2 <- log2(data$creatinine_filled)
  data$uacr_log2 <- log2(data$uacr_filled)
  data
}

row284_model2 <- albumin ~ sleep6 + sex + age + race4 + marital2 + work_moderate
row284_model3 <- albumin ~ sleep6 + sex + age + race4 + marital2 + work_moderate + tp_filled + alt_filled + ast_filled +
  creatinine_filled + uacr_filled + crp_filled + glucose_filled + bmi_filled + hypertension + high_cholesterol + cancer

association <- list(
  id = "row284", row = 284, doi = "10.1186/s12889-022-13524-y",
  cycles = c("2015-2016", "2017-2018"),
  weight = "WTINT2YR", blood_file = "BIOPRO",
  family = "linear", term = "sleep6short",
  published = list(measure = "beta", estimate = -1.00, low = -1.26, high = -0.74, n = 9973, events = NA,
                   contrast = "workday sleep of 5 hours or less vs more than 7 up to 8 hours (serum albumin, g/L)"),
  left_out = c("moderate work activity (PAQ620): 2021-2023 asks only about leisure-time activity"),
  build = function(cycle) {
    slq <- component("SLQ", cycle)
    bio <- component("BIOPRO", cycle, c("LBDSALSI", "LBDSTPSI", "LBXSATSI", "LBXSASSI", "LBDSCRSI", "LBDSGLSI"))
    paq <- component("PAQ", cycle)
    d <- merge_all(demographics(cycle),
                   data.frame(SEQN = slq$SEQN, sleep_hours = slq$SLD012, onset = as.character(slq$SLQ300)),
                   bio, urine_albumin_creatinine(cycle), crp(cycle), body_measures(cycle)[, c("SEQN", "bmi")],
                   component("BPQ", cycle, c("BPQ020", "BPQ080")), component("MCQ", cycle, "MCQ220"),
                   data.frame(SEQN = paq$SEQN, PAQ620 = if (has(paq, "PAQ620")) paq$PAQ620 else NA))
    d$albumin <- d$LBDSALSI
    # Weekday or workday sleep hours (SLD012): asked directly in 2015-2016, derived from usual
    # sleep and wake times (SLQ300, SLQ310) in 2017-2018 and 2021-2023.
    d$sleep6 <- stats::relevel(cut(d$sleep_hours, c(-Inf, 5, 6, 7, 8, 9, Inf), labels = c("short", "5-6", "6-7", "7-8", "8-9", "long")),
                               ref = "7-8")
    d$sex <- factor(d$sex, levels = 1:2, labels = c("male", "female"))
    # Table 1's "Mexican American" holds other Hispanic participants too; "other race" is RIDRETH1 5.
    d$race4 <- factor(ifelse(d$race %in% 1:2, "mexican", ifelse(d$race %in% 3, "white", ifelse(d$race %in% 4, "black", "other"))),
                      levels = c("mexican", "white", "black", "other"))
    d$race4_text <- 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"))
    # Living alone: widowed, divorced, separated, or never married; anyone else (no answer
    # included) as married or living with a partner.
    d$marital2 <- factor(ifelse(d$marital %in% 2:3, "alone", "married or partner"), levels = c("married or partner", "alone"))
    # Moderate work activity, "don't know" as no; not asked in 2021-2023.
    d$work_moderate <- if (cycle == REPLICATION_CYCLE) NA else factor(ifelse(d$PAQ620 %in% 1, "yes", "no"), levels = c("yes", "no"))
    told <- function(x) factor(ifelse(x %in% 1, "yes", ifelse(x %in% 2, "no", "unknown")), levels = c("yes", "no", "unknown"))
    d$hypertension <- told(d$BPQ020)
    d$high_cholesterol <- told(d$BPQ080)
    d$cancer <- told(d$MCQ220)
    d$tp <- d$LBDSTPSI
    d$alt <- d$LBXSATSI
    d$ast <- d$LBXSASSI
    d$creatinine <- d$LBDSCRSI
    d$uacr <- d$acr
    d$glucose <- d$LBDSGLSI
    # Adults with workday sleep hours, a usual sleep time ('HH:MM'; NCHS blanks the times of those
    # whose derived hours are top- or bottom-coded), and serum albumin.
    timed <- grepl("^[0-9]{2}:[0-9]{2}$", d$onset)
    d$in_population <- d$age >= 20 & !is.na(d$sleep_hours) & timed & !is.na(d$albumin)
    d
  },
  constants = row284_constants,
  derive = row284_derive,
  variants = list(
    list(label = "examination weight, as the text says (WTMEC2YR)", weight = "WTMEC2YR"),
    list(label = "unweighted", weighted = FALSE),
    list(label = "Model 1 (published -1.09, -1.43 to -0.76)", formula = albumin ~ sleep6),
    list(label = "Model 1, examination weight", formula = albumin ~ sleep6, weight = "WTMEC2YR"),
    list(label = "Model 2 (published -0.93, -1.24 to -0.61)", formula = row284_model2),
    list(label = "Model 2, examination weight", formula = row284_model2, weight = "WTMEC2YR"),
    list(label = "creatinine and UACR log2-transformed, as Table 1 shows them",
         formula = update(row284_model3, . ~ . - creatinine_filled - uacr_filled + creatinine_log2 + uacr_log2)),
    list(label = "complete cases, no imputation", sample = "own",
         formula = albumin ~ sleep6 + sex + age + race4 + marital2 + work_moderate + tp + alt + ast + creatinine + uacr + crp + glucose +
           bmi + hypertension + high_cholesterol + cancer),
    list(label = "other Hispanic with other race, as the text's labels read", formula = update(row284_model3, . ~ . - race4 + race4_text)),
    list(label = "with a cycle term (2017-2018's analyzer reads albumin lower)", formula = update(row284_model3, . ~ . + cycle))
  ),
  # 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 text says it used the combined examination weight (WTMEC2YR/2), but Table 1's weighted means and shares and all three Table 2 models in all five sleep groups, intervals included, are reproduced to the printed digit only with the interview weight (WTINT2YR). With the examination weight Models 1 and 2 give -1.10 and -0.94 against the printed -1.09 and -0.93; the headline is -1.00 with either."),
    list(kind = "weighting", affects_headline = TRUE, followed = FALSE, evidence = "data",
         detail = "The text says its weighted data were calculated according to the NHANES analytical guidelines, but its intervals are those of weighted least squares without the survey design: so fitted, Model 3 gives the published -1.00 (-1.26 to -0.74), while with the strata and PSUs the interval is -1.42 to -0.58."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The text's race groups are Mexican American, non-Hispanic white, non-Hispanic black, and other race, but Table 1's weighted shares are reproduced exactly only with other Hispanic participants in the Mexican American group; with them in other race, as the labels read, Model 3 is -1.004 against -1.001."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The text and Table 1 give continuous variables as means +- SEs, but the values after +- are weighted SDs (age 49.01 +- 16.10 years in the 5 hours or less group).")
  ),
  choices = list(
    list(choice = "population", decision = "adults 20 and older with workday sleep hours (SLD012), a usual workday sleep time (SLQ300 as 'HH:MM'), and serum albumin (LBDSALSI)",
         reason = "the paper's steps; they reproduce its flow exactly: 11,288 adults, 79 without sleep hours, 35 without a sleep time (all 2017-2018 participants whose derived hours NCHS coded as under 3 or 14 or more and whose times it blanked; 2021-2023 blanks them too), 1,201 without albumin, 9,973 left, and Table 1's group sizes (694, 1,051, 2,248, 2,926, 1,890, 1,164)"),
    list(choice = "weights", decision = "interview weight, WTINT2YR over 2, with strata and PSUs (WTINT2YR in 2021-2023)",
         reason = "the text says it combined the examination weight (WTMEC2YR/2), but the interview weight reproduces every printed number checked: Table 1's weighted means (sleep hours, age, albumin, laboratory values, BMI) and shares (sex, race, marital status, work activity, the three conditions) in all six groups, and Table 2's Models 1, 2, and 3 in all five sleep groups, intervals included (as weighted least squares, which the paper's intervals are; the survey design gives the same estimates with wider intervals). The examination weight misses Models 1 and 2 by 0.01 (-1.10, -0.94) and Table 1's shares by up to half a point; it is a variant, and it too gives -1.00 for the headline. The plan follows the computation the paper's numbers reveal"),
    list(choice = "sleep groups", decision = "5 hours or less, over 5 to 6, over 6 to 7, over 7 to 8 (reference), over 8 to 9, over 9; SLD012's codes for under 3 hours (2) and 14 or more (14) as values",
         reason = "the paper's groups; they reproduce Table 1's sizes and mean hours exactly"),
    list(choice = "race", decision = "Mexican American with other Hispanic (RIDRETH1 1-2), non-Hispanic White, non-Hispanic Black, other race (5)",
         reason = "the text names four groups without placing other Hispanic participants; Table 1's weighted shares are reproduced exactly only with them under 'Mexican American' (the text's labels, other Hispanic with other race, are a variant)"),
    list(choice = "marital status and moderate work activity", decision = "living alone if widowed, divorced, separated, or never married, otherwise married or living with a partner (no answer included); moderate work activity (PAQ620) yes, otherwise no (don't know included)",
         reason = "unstated handling of missing answers; Table 1's weighted shares are reproduced exactly this way"),
    list(choice = "hypertension, high cholesterol, cancer", decision = "told yes, no, or unknown (refused, don't know), unknown as its own level",
         reason = "the paper's three levels; Table 1's weighted shares are reproduced"),
    list(choice = "laboratory covariates", decision = "total protein (g/L), ALT, AST, creatinine (umol/L), UACR (mg/g), hs-CRP (mg/L), and serum glucose from the biochemistry profile (LBDSGLSI, mmol/L), continuous and untransformed",
         reason = "Table 1's glucose means are serum glucose, not fasting plasma glucose; whether creatinine and UACR entered log2-transformed is unstated (Table 1 shows their logs), and untransformed values reproduce Model 3 in all five groups where log2 values give -0.97 (variant)"),
    list(choice = "missing covariates", decision = "each continuous covariate's missing values (BMI 127, UACR 158, hs-CRP 34, AST 16, a few others) filled with its unweighted median over the 9,973 participants on the paper's cycles, reused in 2021-2023; all 9,973 analyzed",
         reason = "unstated; Table 1's weighted means and SDs of BMI, UACR, hs-CRP, and AST are reproduced exactly only with median filling, and so is Model 3. Complete cases (9,675) are a variant"),
    list(choice = "albumin across analyzers", decision = "values as released in both cycles, no crosswalk",
         reason = "a measurement change, as the paper pooled them: NCHS's bridging study found 2017-2018's Roche Cobas 6000 reads albumin about 4.4% lower than 2015-2016's Beckman analyzer, while 2021-2023's Cobas 8000 differs from the Cobas 6000 by 1.5% with no adjustment recommended. Sleep hours also changed form between the two cycles (asked directly, then derived from times), so the pooled contrast mixes cycles; a variant adds a cycle term")
  ),
  formula = row284_model3,
  formula_harmonized = update(row284_model3, . ~ . - work_moderate)
)