# Systemic immune-inflammation index (SII, per 100 units) and hyperlipidemia, NHANES 2015-2016 and
# 2017-March 2020. Mahemuti et al. (2023), Nutrients, doi:10.3390/nu15051177.
# Headline: OR 1.03 (95% CI 1.01-1.05) per 100 units of SII, Model 2 (age, sex, race; Table 3).

association <- list(
  id = "row306", row = 306, doi = "10.3390/nu15051177",
  # 25,531 participants at the start of the flow chart are DEMO_I (9,971) plus P_DEMO (15,560).
  cycles = c("2015-2016", "2017-2020"),
  # Table 1's weighted means and percentages are reproduced with the examination weight (not the
  # fasting subsample's), but the regressions are reproduced only unweighted (variants below).
  weight = "WTMEC2YR", blood_file = "CBC",
  weighted = FALSE,
  family = "logistic", term = "sii100",
  published = list(measure = "OR", estimate = 1.03, low = 1.01, high = 1.05, n = 6117, events = 4265,
                   contrast = "per 100 units of SII (platelets x neutrophils / lymphocytes, counts in 10^3 cells/uL)"),
  left_out = character(),
  build = function(cycle) {
    demo <- demographics(cycle)
    alq <- component("ALQ", cycle)
    d <- merge_all(demo, blood_count(cycle)[, c("SEQN", "platelets", "neutrophils", "lymphocytes")],
                   total_cholesterol(cycle), hdl_cholesterol(cycle), triglycerides(cycle),
                   bp_questions(cycle)[, c("SEQN", "cholesterol_medication")],
                   body_measures(cycle)[, c("SEQN", "bmi")], smoking(cycle), hypertension_status(cycle, 140, 90),
                   diabetes_status(cycle, parts = c("told", "medication")))
    # The paper's mean SII, 459.54 (SD 317.28), is below ours, 507.7 (SD 318.4), on the same 6,117
    # participants, whose mean age (50.70, SD 17.43) and share of men (48.08%) match its text exactly;
    # its quartile estimates differ from ours too (variant), though the per-100 ones agree.
    d$sii <- sii(d$platelets, d$neutrophils, d$lymphocytes)
    d$sii100 <- d$sii / 100
    # The flow chart's 7,298 participants with SII and "hyperlipidemia data" (before the age step)
    # are exactly those with all four lipid values, which only the fasting subsample has.
    lipids_known <- !is.na(d$tc) & !is.na(d$tg) & !is.na(d$ldl) & !is.na(d$hdl)
    low_hdl <- (d$sex == 1 & d$hdl < 40) | (d$sex == 2 & d$hdl < 50)
    d$hyperlipidemia <- ifelse(lipids_known, as.integer(d$tc >= 200 | d$tg >= 150 | d$ldl >= 130 | low_hdl |
                                                          d$cholesterol_medication %in% 1), NA)
    # Variant: also counting those told to take cholesterol medicine (BPQ090D, not asked in 2021-2023).
    bpq <- component("BPQ", cycle)
    told_rx <- if (has(bpq, "BPQ090D")) yes(bpq$BPQ090D[match(d$SEQN, bpq$SEQN)]) else NA
    d$hyperlipidemia_told_rx <- ifelse(lipids_known, as.integer(d$hyperlipidemia %in% 1 | told_rx %in% 1), NA)
    d$sex <- factor(d$sex)
    d$race <- factor(d$race)
    # Model 3's covariates (a variant), coded as Table 1's percentages show them: education splits
    # DMDEDUC2 into 1, 2, and 3-5 (5.4%, 8.4%, 86.2%, under the labels less than high school, high
    # school, more than high school); missing income falls in the middle income group; drinking
    # groups by drinks a day, with non-drinkers and unknowns as light; other missing values are
    # their own level, since the flow chart excludes no one for covariates.
    d$education3 <- with_unclear(ifelse(d$education %in% 1, "1", ifelse(d$education %in% 2, "2", ifelse(d$education %in% 3:5, "3-5", NA))), c("1", "2", "3-5"))
    d$marital3 <- with_unclear(d$marital, 1:3)
    d$pir3 <- factor(ifelse(is.na(d$pir), "1.5-3.5", ifelse(d$pir <= 1.5, "<=1.5", ifelse(d$pir <= 3.5, "1.5-3.5", ">3.5"))),
                     levels = c("<=1.5", "1.5-3.5", ">3.5"))
    drinks <- alq$ALQ130[match(d$SEQN, alq$SEQN)]
    drinks <- ifelse(drinks %in% 1:30, drinks, NA)
    male <- d$sex == "1"
    d$drinking <- factor(ifelse(((!male & drinks >= 3) | (male & drinks >= 4)) %in% TRUE, "excessive",
                                ifelse(((!male & drinks >= 2) | (male & drinks >= 3)) %in% TRUE, "moderate", "light")),
                         levels = c("light", "moderate", "excessive"))
    d$smoking3 <- with_unclear(d$smoking, 1:3)
    d$bmi3 <- with_unclear(cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-30", ">=30")), c("<25", "25-30", ">=30"))
    d$hypertension2 <- with_unclear(d$hypertension, 0:1)
    d$diabetes2 <- with_unclear(d$diabetes, 0:1)
    d$in_population <- !is.na(d$sii) & lipids_known & d$age >= 20
    d
  },
  # SII quartiles for a variant (Table 3's quartile rows): unweighted quartiles of the analytic
  # sample, which Table 2's group sizes (1,529, 1,529, 1,529, 1,530) show.
  constants = function(data) {
    p <- data$in_population %in% TRUE
    list(sii_quartiles = unname(stats::quantile(data$sii[p], c(0.25, 0.5, 0.75))))
  },
  derive = function(data, constants) {
    if (!is.null(constants)) data$sii_q <- cut(data$sii, c(-Inf, constants$sii_quartiles, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))
    data
  },
  formula = hyperlipidemia ~ sii100 + age + sex + race,
  variants = list(
    list(label = "weighted, examination weight (Table 1's weight)", weighted = TRUE),
    list(label = "weighted, fasting subsample weight", weight = "WTSAF2YR", weighted = TRUE),
    list(label = "unweighted, positive fasting subsample weight only", weight = "WTSAF2YR",
         derive = function(data, constants) { data$in_population <- data$in_population & !is.na(data$w) & data$w > 0; data }),
    list(label = "crude, unweighted (published 1.04, 1.02-1.06)", formula = hyperlipidemia ~ sii100),
    list(label = "crude, weighted", formula = hyperlipidemia ~ sii100, weighted = TRUE),
    list(label = "Model 3, unweighted (published 1.02, 1.00-1.04)",
         formula = hyperlipidemia ~ sii100 + age + sex + race + marital3 + pir3 + education3 + drinking + smoking3 + bmi3 + hypertension2 + diabetes2),
    list(label = "Model 3, weighted",
         formula = hyperlipidemia ~ sii100 + age + sex + race + marital3 + pir3 + education3 + drinking + smoking3 + bmi3 + hypertension2 + diabetes2,
         weighted = TRUE),
    list(label = "outcome also counts told to take cholesterol medicine (BPQ090D)", formula = hyperlipidemia_told_rx ~ sii100 + age + sex + race),
    list(label = "quartile 4 vs 1, Model 2, unweighted (published 1.27, 1.08-1.50)",
         formula = hyperlipidemia ~ sii_q + age + sex + race, term = "sii_qQ4")
  ),
  # 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,
         detail = "The Statistical Analysis says a weighting strategy was used, and Table 1's weighted means and percentages are those with the examination weight, but Table 3's estimates and interval widths are reproduced only unweighted (Model 2: 1.030, 1.009-1.051 unweighted and 1.061, 1.020-1.103 weighted, against 1.03, 1.01-1.05)."),
    list(kind = "coding", affects_headline = TRUE, followed = FALSE,
         detail = "The Methods compute SII as platelets x neutrophils / lymphocytes, but the Results' mean SII (459.54, SD 317.28) is below what that formula gives for the same 6,117 participants (507.7, SD 318.4), whose mean age (50.70) and share of men (48.08%) match the text exactly; Table 3's quartile estimates differ too (Model 2, quartile 4 vs 1: 1.27 published, 1.44 here), though the per-100 estimate agrees."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE,
         detail = "The abstract attributes its estimate, 1.03 (1.01, 1.05), to a multivariate linear regression analysis, but it is Table 3's Model 2 odds ratio from logistic regression, as the Methods describe and as reproduced (1.030)."),
    list(kind = "coding", affects_headline = FALSE, followed = NA,
         detail = "The Covariates section and Table 1 give education as less than high school, high school, and more than high school, but Table 1's shares (5.37%, 8.42%, 86.20% with hyperlipidemia) are those of less than 9th grade, 9th to 11th grade, and high school graduate or more (DMDEDUC2 1, 2, and 3-5)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 4 gives the estimate for women with SII under 958.14 as 1.0006 (1.0002, 1.1010), an upper bound that the estimate and its lower bound rule out (on the log scale they put it near 1.0010)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2's light alcohol consumption row repeats the never-smoker row's percentages (55.62, 58.51, 55.16, 52.74), so its drinking rows sum to less than 100% in every quartile.")
  ),
  choices = list(
    list(choice = "cycles", decision = "2015-2016 and 2017-March 2020", reason = "the paper says 2015-2020; its 25,531 participants are the DEMO_I and P_DEMO totals"),
    list(choice = "hyperlipidemia data required", decision = "total cholesterol, triglycerides, LDL (Friedewald, LBDLDL) and HDL all measured", reason = "this gives exactly the flow chart's 7,298 participants with SII and hyperlipidemia data, and 6,117 aged 20 or older"),
    list(choice = "weights", decision = "unweighted, with the examination weight (WTMEC2YR, WTMECPRP) as the paper's weight", reason = "Table 1's weighted means and percentages match the examination weight (not the fasting weight); the crude and Model 2 estimates are reproduced only unweighted; the text says only that a weighting strategy was used. The fasting weight would also drop 197 participants with lipid values but a zero fasting weight, all in 2017-March 2020"),
    list(choice = "lipid thresholds", decision = "total cholesterol 200 mg/dL or more, triglycerides 150 or more, LDL 130 or more, HDL under 40 (men) or 50 (women)", reason = "the paper prints no inequality signs; these are the NCEP ATP III thresholds it cites"),
    list(choice = "cholesterol-lowering drugs", decision = "now taking prescribed cholesterol medicine (BPQ100D; BPQ101D in 2021-2023), missing counted as no", reason = "the paper says persons who reported using cholesterol-lowering drugs. This gives 4,182 events against the published 4,265; adding those told to take medicine (BPQ090D) gives 4,272 and the same estimate (variant)"),
    list(choice = "age and race in Model 2", decision = "age continuous; race in RIDRETH1's five groups", reason = "Table 1 reports age as a mean and race in those five groups"),
    list(choice = "Model 3 covariates (variant only)", decision = "coded as Table 1's percentages show (see build)", reason = "Model 3 is not the headline")
  )
)