# Serum albumin and depressive symptoms (PHQ-9 of 10 or more) in adults aged 20 and older, NHANES 2005-2018.
# Zhang et al. (2023), BMC Psychiatry, doi:10.1186/s12888-023-04935-1.
# Headline: OR 0.77 (95% CI 0.60-0.99) for the highest against the lowest albumin quartile, Model III
# (abstract; Fig. 2 and the Results print 0.60-0.97).

# Model III (Fig. 2 legend).
model3 <- depressed ~ albumin_q + age + race + sex + education + bmi + drinking + smoking_now + heart_failure + chd +
  liver_condition + cancer + diabetes + thyroid

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

# Albumin quartiles (g/L): unweighted quartiles of the analytic sample, intervals closed on the right.
quartile_cutpoints <- function(x) unname(stats::quantile(x, c(0.25, 0.5, 0.75)))
albumin_quartile <- function(x, cutpoints) cut(x, c(-Inf, cutpoints, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))

# The sample the paper analyzed, with quartiles of its own (see choices).
published_sample <- function(data, constants) {
  data$albumin_q <- albumin_quartile(data$albumin, constants$cutpoints_published_sample)
  data$in_population <- data$in_population & data$published_sample
  data
}

# All adults, as the paper's text describes its sample.
all_adults <- function(data, constants) {
  data$albumin_q <- albumin_quartile(data$albumin, constants$cutpoints)
  data
}

association <- list(
  id = "row191", row = 191, doi = "10.1186/s12888-023-04935-1",
  cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTMEC2YR", blood_file = "BIOPRO",
  family = "logistic", term = "albumin_qQ4",
  published = list(measure = "OR", estimate = 0.77, low = 0.60, high = 0.99, n = 13681, events = 1551,
                   contrast = "highest vs lowest quartile of serum albumin"),
  left_out = character(),
  build = function(cycle) {
    demo <- demographics(cycle)
    alq <- component("ALQ", cycle)
    d <- merge_all(demo, phq9(cycle), component("BIOPRO", cycle, "LBDSALSI"), body_measures(cycle)[, c("SEQN", "bmi")],
                   smoking(cycle), alcohol(cycle)[, c("SEQN", "alcohol3")], medical_conditions(cycle),
                   component("MCQ", cycle, "MCQ160M"), component("DIQ", cycle, "DIQ010"))
    # Ever drank 4 (women) or 5 (men) drinks almost every day: ALQ151 from 2011, ALQ150 (5 drinks)
    # before. It is not asked of those who never had 12 drinks in their life (from 2017, any drink).
    daily <- if (has(alq, "ALQ151")) alq$ALQ151 else alq$ALQ150
    d$daily_answer <- daily[match(d$SEQN, alq$SEQN)]
    d$depressed <- as.integer(d$phq9 >= 10)
    d$albumin <- d$LBDSALSI
    d$sex <- factor(d$sex)
    d$race <- factor(d$race)
    d$education <- factor(d$education)
    # "Do you smoke now?" (SMQ040, every day or some days), never smokers counted as nonsmokers.
    d$smoking_now <- factor(ifelse(d$smoking %in% 3, 1, ifelse(d$smoking %in% 1:2, 0, NA)))
    # The daily-drinking item Table 1's counts show was used, lifetime abstainers counted as no.
    d$drinking <- factor(ifelse(d$daily_answer %in% 1, 1, ifelse(d$daily_answer %in% 2 | d$alcohol3 %in% 1, 0, NA)))
    d$drinking_now <- factor(ifelse(d$alcohol3 %in% 3, 1, ifelse(d$alcohol3 %in% 1:2, 0, NA)))
    for (v in c("heart_failure", "chd", "liver_condition", "cancer")) d[[v]] <- factor(d[[v]])
    d$thyroid <- factor(yes(d$MCQ160M))
    d$diabetes <- factor(ifelse(d$DIQ010 %in% 1:3, d$DIQ010, NA))
    # The published sample: complete-case analysis on SMQ040 (asked only of those who smoked 100
    # cigarettes) and on the daily-drinking item (not asked of lifetime abstainers) kept only ever
    # smokers who had drunk, and PHQ-9 scores of 27 were lost.
    d$published_sample <- d$smoking %in% 2:3 & d$daily_answer %in% 1:2 & (d$phq9 < 27) %in% TRUE
    d$in_population <- d$age >= 20 & !is.na(d$phq9) & !is.na(d$albumin)
    d
  },
  constants = function(data) {
    sample <- analytic(data)
    list(cutpoints = quartile_cutpoints(data$albumin[sample]),
         cutpoints_published_sample = quartile_cutpoints(data$albumin[sample & data$published_sample]))
  },
  # The published estimate comes from the sample the paper's complete-case coding left (ever
  # smokers who had drunk, PHQ-9 under 27), so the replication follows it; all adults, as the text
  # describes, is a variant. In 2021-2023 the same questions skip the same people.
  derive = published_sample,
  variants = list(
    list(label = "Model I (published 0.52, 0.42-0.65)", formula = depressed ~ albumin_q),
    list(label = "Model I, unweighted", formula = depressed ~ albumin_q, weighted = FALSE),
    list(label = "Model II (published 0.60, 0.48-0.76)", formula = depressed ~ albumin_q + age + sex + race),
    list(label = "per g/L (published 0.98, 0.95-1.00)", term = "albumin",
         formula = depressed ~ albumin + age + race + sex + education + bmi + drinking + smoking_now + heart_failure + chd + liver_condition + cancer + diabetes + thyroid),
    list(label = "PHQ-9 score, linear (Fig. 3: -0.39, -0.67 to -0.11)", family = "linear",
         formula = phq9 ~ albumin_q + age + race + sex + education + bmi + drinking + smoking_now + heart_failure + chd + liver_condition + cancer + diabetes + thyroid),
    list(label = "all adults, as the text describes the sample", derive = all_adults),
    list(label = "all adults, Model I", derive = all_adults, formula = depressed ~ albumin_q),
    list(label = "all adults, Model II", derive = all_adults, formula = depressed ~ albumin_q + age + sex + race),
    list(label = "all adults, unweighted", derive = all_adults, weighted = FALSE),
    list(label = "all adults, drinking in the past 12 months (the text's question)", derive = all_adults,
         formula = depressed ~ albumin_q + age + race + sex + education + bmi + drinking_now + smoking_now + heart_failure + chd + liver_condition + cancer + diabetes + thyroid)
  ),
  # 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 text describes adults aged 20 and older with complete data, but n = 13,681, its 1,551 cases, every count in Table 1, and Models I and II are reproduced only in ever smokers who had drunk alcohol (the smoking question used, SMQ040, skips never smokers, and the drinking item skips lifetime abstainers), less 8 with a PHQ-9 score of 27. All adults give n = 31,716 and an OR of 0.89 (0.75-1.06), against 0.82 (0.64-1.06) in that sample."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods code drinking status from the question 'Do you drink alcohol now?', which NHANES doesn't ask, and Table 1's 3,588 drinkers (25.23% weighted) are reproduced only by ever having had 4 or 5 or more drinks almost every day (ALQ151; ALQ150, 5 drinks, through 2009-2010)."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE,
         detail = "The abstract prints the headline's interval as 0.60-0.99, and the Results and Fig. 2 print it as 0.60-0.97 (P = 0.044); the file keeps the abstract's."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's header gives quartile sizes of 2,745, 2,744, 3,292, and 4,900, those of the same cutpoints with groups closed on the left, while its body (3,971, 3,159, 3,167, 3,384) and Fig. 2's cases per quartile (577, 369, 316, 289) are those of groups closed on the right.")
  ),
  choices = list(
    list(choice = "study population", decision = "the sample the paper's computation analyzed: adults aged 20 and older with a complete PHQ-9 under 27, serum albumin, and Model III covariates, which complete-case coding limited to ever smokers (SMQ040 is asked only of them) who had drunk (the daily-drinking item is not asked of lifetime abstainers)", reason = "that sample reproduces n = 13,681 and 1,551 cases, every count in Table 1 (8,226 men, 6,263 current smokers, 3,588 drinkers, quartiles of 3,971, 3,159, 3,167, 3,384), its weighted shares (56.74% men, 44.52% current smokers), and Models I and II; the published estimate comes from it. The text describes adults aged 20 and older; that version is a variant (n = 31,716)"),
    list(choice = "confidence interval of the headline", decision = "the abstract's 0.60-0.99 kept as published", reason = "the abstract prints 0.77 (0.60-0.99); the Results text and Fig. 2 print 0.77 (0.60-0.97), P = 0.044"),
    list(choice = "albumin quartiles", decision = "unweighted quartiles of albumin in g/L (LBDSALSI) over the analytic sample, intervals closed on the right: 40 or less, 41-42, 43-44, 45 or more", reason = "cutpoints and unit are not reported. Unweighted quartiles of the published sample (40, 42, 44 g/L) reproduce Table 1's quartile sizes and Fig. 2's cases per quartile (577, 369, 316, 289) exactly; weighted ones (41, 43, 45) do not. All adults give the same unweighted cutpoints (weighted: 40, 43, 45). Table 1's header row (2,745, 2,744, 3,292, 4,900) is the same cutpoints with intervals closed on the left, which its body and Fig. 2 do not use. The continuous estimate is per g/L (variant)"),
    list(choice = "albumin across analyzers", decision = "values as released in every cycle, no crosswalk", reason = "the paper mentions none and Table 1 is reproduced without one. NCHS's bridging study says 2017-2018's Roche Cobas 6000 reads 4.4% below the earlier Beckman analyzer and recommends its equation when pooling; 2021-2023's Cobas 8000 differs from the Cobas 6000 by 1.5%, with no adjustment recommended. So the fixed cutpoints, set mostly on Beckman values, will put more of 2021-2023 in the lowest quartile (about half of 2017-2018's adults were in it)"),
    list(choice = "weights", decision = "examination weight (WTMEC2YR) over 7, with strata and PSUs (WTPH2YR from BIOPRO_L in 2021-2023)", reason = "the paper combined MEC examination weights over seven cycles as NCHS directs; on the published sample, weighted fits reproduce Model I (Q2-Q4 0.70, 0.62, 0.52) and Model II (0.74, 0.70, 0.60), CIs included, and unweighted ones do not"),
    list(choice = "smoking status", decision = "current smoker (SMQ040 every day or some days) against former or never smoker", reason = "the text's 'Do you smoke now?'; Table 1's 6,263 smokers are SMQ040 1-2"),
    list(choice = "drinking status", decision = "ever drank 4 (women) or 5 (men) drinks almost every day (ALQ151; ALQ150, 5 drinks for everyone, in 2005-2010), lifetime abstainers as no", reason = "no NHANES item reads 'Do you drink alcohol now?'; this item alone reproduces Table 1's 3,588 drinkers (25.23% weighted); drinking in the past 12 months is a variant"),
    list(choice = "PHQ-9", decision = "sum of the nine items, missing if any is missing, refused, or don't know; 10 or more", reason = "the paper's definition; with the published sample's steps it reproduces n and cases exactly"),
    list(choice = "BMI", decision = "continuous", reason = "Table 1 reports 'BMI status' as a mean, and with continuous BMI the same covariates reproduce Fig. 3's linear Model III exactly (variant); BMI in three groups does not"),
    list(choice = "education", decision = "DMDEDUC2 in its five levels", reason = "Table 1's five categories, reproduced exactly"),
    list(choice = "diabetes", decision = "DIQ010 yes, no, borderline", reason = "Table 1's three levels"),
    list(choice = "chronic conditions", decision = "ever told of congestive heart failure (MCQ160B), coronary heart disease (MCQ160C), a liver condition (MCQ160L), cancer (MCQ220), a thyroid problem (MCQ160M); yes or no", reason = "the Methods' self-reported history; Table 1's counts are reproduced exactly"),
    list(choice = "pregnancy", decision = "pregnant women kept", reason = "no pregnancy exclusion is described, and the published sample is reproduced without one"),
    list(choice = "Model III quartile estimates", decision = "the paper's covariates and codings, unchanged", reason = "on the published sample, Fig. 3's linear Model III (Q2-Q4 -0.31, -0.52, -0.39, CIs included) and Fig. 2's continuous Model III (0.98, 0.95-1.00) are reproduced exactly, but Fig. 2's logistic Model III quartiles are not: we get 0.85, 0.86, 0.82 (0.64-1.06) against 0.81, 0.80, 0.77. Codings tried without a match: education, race, or diabetes as numbers; three or four BMI groups; SMQ040 in three levels")
  ),
  formula = model3
)