# Natural log of the non-HDL to HDL cholesterol ratio (NHHR) and gallstones in adults aged 50 or
# younger, NHANES 2017-March 2020. Cheng et al. (2024), Lipids Health Dis, doi:10.1186/s12944-024-02262-2.
# Headline: OR 1.77 (95% CI 1.11-2.83) per unit of ln NHHR, Model 3 (Table 2).

association <- list(
  id = "row040", row = 40, doi = "10.1186/s12944-024-02262-2",
  cycles = "2017-2020",
  weight = "WTMEC2YR", blood_file = "HDL",
  # The paper names NHANES's sample weights only for Table 1's tests. Its crude quartile odds
  # ratios equal the unweighted ones from Table 1's counts, and its crude and Model 2 confidence
  # intervals are as wide as unweighted fits give for Table 1's 2,117 people with 157 cases
  # (weighted intervals would be wider), so the models were unweighted.
  weighted = FALSE,
  family = "logistic", term = "ln_nhhr",
  published = list(measure = "OR", estimate = 1.77, low = 1.11, high = 2.83, n = 3772, events = NULL,
                   contrast = "per unit of ln NHHR ((total cholesterol - HDL cholesterol) / HDL cholesterol)"),
  left_out = c("physical activity at work (PAQ620, whether work involves moderate-intensity activity): 2021-2023 asks only about leisure-time activity"),
  build = function(cycle) {
    demo <- demographics(cycle)
    conditions <- medical_conditions(cycle)[, c("SEQN", "gallstones", "chd", "heart_attack", "copd")]
    paq <- component("PAQ", cycle)
    d <- merge_all(demo, total_cholesterol(cycle), hdl_cholesterol(cycle), conditions,
                   body_measures(cycle)[, c("SEQN", "bmi")], component("SMQ", cycle, "SMQ020"),
                   component("ALQ", cycle, c("ALQ111", "ALQ151")), component("BPQ", cycle, "BPQ020"),
                   component("DIQ", cycle, "DIQ010"), leisure_activity(cycle)[, c("SEQN", "sedentary_minutes")])
    d$work_moderate <- if (has(paq, "PAQ620")) yes(paq$PAQ620[match(d$SEQN, paq$SEQN)]) else NA
    d$nhhr <- (d$tc - d$hdl) / d$hdl
    d$ln_nhhr <- log(d$nhhr)
    # The paper's quartiles of NHHR (Table 1), for the quartile variant.
    d$nhhr_q <- cut(d$nhhr, c(-Inf, 1.79, 2.51, 3.46, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))
    d$tc_mmol <- cholesterol_mmol(d$tc)
    d$hdl_mmol <- cholesterol_mmol(d$hdl)
    d$sex <- factor(d$sex)
    d$race_code <- d$race
    d$race <- factor(d$race)
    d$pir3 <- factor(ifelse(d$pir <= 1.3, "<=1.3", ifelse(d$pir < 3.5, "1.3-3.5", ">=3.5")), levels = c("<=1.3", "1.3-3.5", ">=3.5"))
    d$education3 <- factor(education3(d$education))
    d$bmi3 <- factor(ifelse(d$bmi <= 25, "<=25", ifelse(d$bmi < 30, "25-30", ">=30")), levels = c("<=25", "25-30", ">=30"))
    d$smoking <- yes(d$SMQ020)
    # Fig. 2: "Yes: ever have 4/5 or more drinks every day" (ALQ151); those who never drank skip
    # the item and have not.
    d$alcohol <- ifelse(d$ALQ151 %in% 1, 1, ifelse(d$ALQ151 %in% 2 | d$ALQ111 %in% 2, 0, NA))
    d$alcohol_asked <- yes(d$ALQ151)
    d$hypertension <- yes(d$BPQ020)
    d$diabetes <- ifelse(d$DIQ010 %in% 1, 1, ifelse(d$DIQ010 %in% 2:3, 0, NA))
    d$sedentary <- d$sedentary_minutes
    # Fig. 1: aged 50 or younger, with total and HDL cholesterol, and with an answer to the
    # gallstone question (asked from age 20).
    d$in_population <- d$age <= 50 & !is.na(d$tc) & !is.na(d$hdl) & !is.na(d$gallstones)
    d
  },
  formula = gallstones ~ ln_nhhr + sex + age + race + pir3 + education3 + bmi3 + smoking + alcohol + tc_mmol + hdl_mmol +
    sedentary + hypertension + diabetes + chd + heart_attack + copd + work_moderate,
  formula_harmonized = gallstones ~ ln_nhhr + sex + age + race + pir3 + education3 + bmi3 + smoking + alcohol + tc_mmol + hdl_mmol +
    sedentary + hypertension + diabetes + chd + heart_attack + copd,
  variants = list(
    list(label = "weighted (WTMECPRP)", weighted = TRUE),
    list(label = "crude, unweighted (published 1.42, 1.00-2.00)", formula = gallstones ~ ln_nhhr),
    list(label = "crude, weighted", formula = gallstones ~ ln_nhhr, weighted = TRUE),
    list(label = "Model 2, unweighted (published 1.80, 1.22-2.66)", formula = gallstones ~ ln_nhhr + sex + age + race),
    list(label = "Model 2, weighted", formula = gallstones ~ ln_nhhr + sex + age + race, weighted = TRUE),
    list(label = "Model 3 without total and HDL cholesterol", formula = gallstones ~ ln_nhhr + sex + age + race + pir3 + education3 + bmi3 +
      smoking + alcohol + sedentary + hypertension + diabetes + chd + heart_attack + copd + work_moderate),
    list(label = "Model 3 without BMI, total and HDL cholesterol", formula = gallstones ~ ln_nhhr + sex + age + race + pir3 + education3 +
      smoking + alcohol + sedentary + hypertension + diabetes + chd + heart_attack + copd + work_moderate),
    list(label = "Model 3, BMI and PIR continuous, race as one 1-5 code", formula = gallstones ~ ln_nhhr + sex + age + race_code + pir +
      education3 + bmi + smoking + alcohol + tc_mmol + hdl_mmol + sedentary + hypertension + diabetes + chd + heart_attack + copd + work_moderate),
    list(label = "Model 3, never-drinkers missing for alcohol", sample = "own", formula = gallstones ~ ln_nhhr + sex + age + race + pir3 +
      education3 + bmi3 + smoking + alcohol_asked + tc_mmol + hdl_mmol + sedentary + hypertension + diabetes + chd + heart_attack + copd + work_moderate),
    list(label = "Model 3, NHHR Q4 vs Q1 (published 1.86, 1.04-3.32)", term = "nhhr_qQ4", formula = gallstones ~ nhhr_q + sex + age + race +
      pir3 + education3 + bmi3 + smoking + alcohol + tc_mmol + hdl_mmol + sedentary + hypertension + diabetes + chd + heart_attack + copd + work_moderate)
  ),
  # 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 = FALSE,
         detail = "The text analyzes 3,772 participants with missing covariates multiply imputed, but Table 1 covers 2,117, Table 2's crude quartile estimates are those of Table 1's counts (quartile 4: 1.51, 0.95 to 2.41), and its crude and Model 2 intervals for ln NHHR are as wide as unweighted fits on 2,117 people with 157 cases give (standard errors 0.177 and 0.199 here, scaled to that size, against the published 0.177 and 0.199)."),
    list(kind = "model", affects_headline = TRUE, followed = FALSE,
         detail = "Table 2's note puts BMI and continuous total and HDL cholesterol in Model 3, but its interval (1.11 to 2.83) is about a quarter as wide on the log scale as a model with both cholesterol terms gives (0.26 to 8.68). Its estimate is also close to Model 2's (1.77 against 1.80), while adjusting for BMI lowers it here (1.09 with BMI, 1.66 without, both without the cholesterol terms) and the paper's own BMI subgroups give 0.94, 1.35, and 1.37."),
    list(kind = "coding", affects_headline = TRUE, followed = FALSE,
         detail = "Fig. 2 defines alcohol consumption as ever having 4/5 or more drinks every day (ALQ151), but Table 1's 18.7% drinkers (396 of 2,117) matches ALQ151 yes together with those who never drank (704 of the 3,772 here, 18.7%), not ALQ151 yes alone (422, 11.2%)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's 'Male' row (70.89%, 61.81%, 50.85%, 32.26% by NHHR quartile) gives the share of women (70.70%, 60.31%, 50.79%, 30.93% of the study population here)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's PIR rows hold more people than quartiles 2 and 3 have (629 and 644 of 529, their shares summing to 118.85% and 121.65%)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's note and the Methods give continuous variables as mean and standard error, but the values are standard deviations (age 32.67 +- 8.96 in a quartile of 529)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "In Fig. 4 the small COPD and total cholesterol groups carry the narrow intervals and the large ones the wide intervals, the reverse of what their sizes give: COPD yes 1.58 (1.04 to 2.41) and no 1.72 (0.42 to 7.16), with 157 and 3,610 people here; total cholesterol 6.2 mmol/L or more 1.74 (1.11 to 2.73) and under 6.2 1.43 (0.49 to 4.18), with 287 and 3,485.")
  ),
  choices = list(
    list(choice = "weights", decision = "unweighted",
         reason = "the paper does not say the regressions were weighted; Table 2's crude quartile ORs equal the unweighted ORs from Table 1's counts, and its crude and Model 2 intervals are as wide as unweighted fits on 2,117 people with 157 cases give (our unweighted standard errors on all 3,772, scaled to that size: 0.177 and 0.199; published 0.177 and 0.199)"),
    list(choice = "analytic sample", decision = "Fig. 1's 3,772 (aged 50 or younger, total and HDL cholesterol, a yes or no to MCQ550), complete cases for Model 3's covariates",
         reason = "Fig. 1 is reproduced exactly (15,560, 10,719, 6,683, 3,772); the paper says it multiply imputed missing covariates, but its models were fitted on Table 1's 2,117 participants, a subset the paper does not define (not the fasting subsample, people who fasted 8 or 9 hours, complete cases, or those under 50), so complete cases of the stated population are used"),
    list(choice = "total and HDL cholesterol in Model 3", decision = "continuous, mmol/L, as Table 2's note and the covariate text say",
         reason = "stated; with both in the model ln NHHR is close to a function of them and its interval is about four times as wide as the published one, so the published Model 3 cannot have held both as continuous terms (variant without them)"),
    list(choice = "BMI", decision = "Fig. 2's categories: 25 or less, 25 to 30, 30 or more", reason = "Fig. 2 defines the covariate by these categories; Table 1 uses them"),
    list(choice = "income to poverty ratio", decision = "Table 1's categories: 1.3 or less, 1.3 to 3.5, 3.5 or more", reason = "Table 1"),
    list(choice = "race and education", decision = "RIDRETH1's five groups; education as less than high school (DMDEDUC2 1-2), high school (3), college or above (4-5)",
         reason = "Table 1's categories; its 61% 'college graduate or above' matches DMDEDUC2 4 and 5 together (61.1%)"),
    list(choice = "smoking, alcohol, physical activity", decision = "SMQ020 yes; ALQ151 yes, with never-drinkers (ALQ111 no) as no; PAQ620 yes",
         reason = "Fig. 2's definitions; never-drinkers have not had 4/5 drinks every day"),
    list(choice = "hypertension and diabetes", decision = "self-reported diagnosis: BPQ020 yes; DIQ010 yes, borderline as no",
         reason = "unstated; Table 1's shares (19.0% and 5.6%) match BPQ020 (18.9%) and DIQ010 yes (5.6%) in the population, not borderline counted as diabetes (7.4%)"),
    list(choice = "CHD, heart attack, COPD", decision = "MCQ160C, MCQ160E, MCQ160P yes", reason = "unstated; the 2017-March 2020 questions for these conditions"),
    list(choice = "sedentary time", decision = "PAD680 minutes, continuous", reason = "the covariate text says minutes")
  )
)