# Systemic immune-inflammation index (SII, per 100 units) and self-reported diagnosed diabetes,
# adults aged 20 and older, NHANES 2017-March 2020. Nie, Zhou, Wang, and Kan (2023), Front
# Endocrinol, doi:10.3389/fendo.2023.1245199.
# Headline: OR 1.04 (95% CI 1.02-1.06), p = 0.0006, per 100 units of SII, Model 3 (Table 3).

# A yes/no answer with refused and don't know counted as no: Table 1's percentages for smoking,
# high blood pressure, and vigorous work activity have those answers in their denominators.
row303_yes <- function(x) factor(ifelse(x %in% 1, "yes", ifelse(is.na(x), NA, "no")), levels = c("no", "yes"))

row303_build <- function(cycle) {
  demo <- demographics(cycle)
  diq <- component("DIQ", cycle, "DIQ010")
  bio <- component("BIOPRO", cycle, c("LBDSBUSI", "LBXSCLSI", "LBDSALSI"))
  d <- merge_all(demo, blood_count(cycle)[, c("SEQN", "wbc", "platelets", "neutrophils", "lymphocytes")], diq, bio,
                 component("ALB_CR", cycle, "URXUMA"), dietary_totals(cycle, c("PROT", "CARB", "SUGR", "TFAT", "CHOL"), days = 1),
                 component("ALQ", cycle, "ALQ111"), component("SMQ", cycle, "SMQ020"), component("BPQ", cycle, "BPQ020"),
                 body_measures(cycle)[, c("SEQN", "bmi")])
  # Vigorous work activity (PAQ605) is not asked in 2021-2023, which asks only about leisure time.
  if (cycle == REPLICATION_CYCLE) d$PAQ605 <- NA else d <- merge_all(d, component("PAQ", cycle, "PAQ605"))
  d$sii <- sii(d$platelets, d$neutrophils, d$lymphocytes)
  d$sii100 <- d$sii / 100
  # Diagnosed diabetes (DIQ010) yes against no; the flow chart's 244 with missing diabetes data are
  # exactly the 239 borderline and 5 don't know answers.
  d$diabetes <- ifelse(d$DIQ010 %in% 1, 1, ifelse(d$DIQ010 %in% 2, 0, NA))
  d$diabetes_borderline_no <- ifelse(d$DIQ010 %in% 1, 1, ifelse(d$DIQ010 %in% 2:3, 0, NA))
  d$sex <- factor(d$sex)
  d$race <- factor(d$race)
  d$marital3 <- factor(d$marital, levels = 1:3, labels = c("married or partner", "widowed, divorced, separated", "never married"))
  d$pir3 <- cut(d$pir, c(-Inf, 1.5, 3.5, Inf), labels = c("0-1.5", "1.5-3.5", ">3.5"))
  d$education3 <- factor(education3(d$education), levels = 1:3, labels = c("less than high school", "high school", "more than high school"))
  d$bmi3 <- cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-<30", ">=30"))
  d$smoke <- row303_yes(d$SMQ020)
  d$high_bp <- row303_yes(d$BPQ020)
  d$vigorous_work <- row303_yes(d$PAQ605)
  d$drink <- factor(ifelse(d$ALQ111 %in% 1, "yes", ifelse(d$ALQ111 %in% 2, "no", NA)), levels = c("no", "yes"))
  d$bun <- d$LBDSBUSI
  d$chloride <- d$LBXSCLSI
  d$albumin <- d$LBDSALSI
  d$urine_albumin <- d$URXUMA
  d$in_population <- d$age >= 20 & !is.na(d$sii) & !is.na(d$diabetes)
  d
}

# Quartiles of SII/100 for a variant (Table 3's quartile rows): unweighted quartiles of the study
# population, which Table 2's group sizes (1,969, 1,969, 1,969, 1,970) show.
row303_constants <- function(data) {
  p <- data$in_population %in% TRUE
  list(sii100_quartiles = unname(stats::quantile(data$sii100[p], c(0.25, 0.5, 0.75))))
}

row303_derive <- function(data, constants) {
  if (is.null(constants)) constants <- row303_constants(data)
  data$sii100_q <- cut(data$sii100, c(-Inf, constants$sii100_quartiles, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))
  data
}

row303_model3 <- diabetes ~ sii100 + age + sex + race + marital3 + pir3 + education3 + drink + smoke + bmi3 + high_bp + bun + wbc +
  platelets + chloride + PROT + CARB + SUGR + TFAT + CHOL + vigorous_work

association <- list(
  id = "row303", row = 303, doi = "10.3389/fendo.2023.1245199",
  # The flow chart's 15,560 participants are P_DEMO's; the paper calls them 2017-2018 and 2019-2020.
  cycles = "2017-2020",
  weight = "WTMEC2YR", blood_file = "CBC",
  # The paper says it used a weighting approach, and Table 1's means and percentages are those
  # with the examination weight (WTMECPRP), but the crude and Model 2 estimates and their CIs, by
  # SII/100 and by quartile, are reproduced only unweighted (variants), as EmpowerStats fits them.
  weighted = FALSE,
  family = "logistic", term = "sii100",
  published = list(measure = "OR", estimate = 1.04, low = 1.02, high = 1.06, n = 7877, events = 1266,
                   contrast = "per 100 units of SII (platelets x neutrophils / lymphocytes, counts in 10^3 cells/uL)"),
  left_out = c("vigorous work activity (PAQ605): 2021-2023 asks no work-activity questions"),
  build = row303_build,
  constants = row303_constants,
  derive = row303_derive,
  formula = row303_model3,
  formula_harmonized = update(row303_model3, . ~ . - vigorous_work),
  variants = list(
    list(label = "weighted (examination weight), as the paper says", weighted = TRUE),
    list(label = "crude, unweighted (published 1.04, 1.02-1.05)", formula = diabetes ~ sii100, sample = "own"),
    list(label = "crude, weighted", formula = diabetes ~ sii100, sample = "own", weighted = TRUE),
    list(label = "Model 2, unweighted (published 1.04, 1.02-1.05)", formula = diabetes ~ sii100 + age + sex + race, sample = "own"),
    list(label = "Model 2, weighted", formula = diabetes ~ sii100 + age + sex + race, sample = "own", weighted = TRUE),
    list(label = "Model 3 without drink", formula = update(row303_model3, . ~ . - drink), sample = "own"),
    list(label = "Model 3 as Table 1 and the Covariates list: serum albumin, no drink",
         formula = update(row303_model3, . ~ . - drink + albumin), sample = "own"),
    list(label = "Model 3 as the Methods list: albumins, no drink, WBC, platelets",
         formula = update(row303_model3, . ~ . - drink - wbc - platelets + albumin + urine_albumin), sample = "own"),
    list(label = "quartile 4 vs 1, Model 3 (published 1.31, 1.05-1.63)", formula = update(row303_model3, . ~ . - sii100 + sii100_q), term = "sii100_qQ4"),
    list(label = "borderline counted as no diabetes", formula = update(row303_model3, diabetes_borderline_no ~ .), sample = "own",
         build = function(cycle) { d <- row303_build(cycle); d$in_population <- d$age >= 20 & !is.na(d$sii) & !is.na(d$diabetes_borderline_no); d })
  ),
  # 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 abstract reports weighted multivariate regression, and Table 1's means are those with the examination weight, but Table 3's Model 1 and Model 2 estimates and intervals are reproduced only unweighted (1.037, 1.022-1.053 and 1.036, 1.020-1.053 against 1.04, 1.02-1.05 each; weighted 1.044, 1.012-1.077 and 1.031, 0.997-1.066), and Model 3's interval has the unweighted width (ours 1.002-1.052 unweighted, 0.955-1.070 weighted)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The Methods compute SII by dividing the lymphocyte count by the platelet count and multiplying by the neutrophil count, but the published mean SII (524.91, SD 358.90), reproduced here, is that of platelets x neutrophils / lymphocytes, the formula the Introduction gives."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The abstract reads the Model 3 odds ratio (1.04, 1.02-1.06) as the change for each additional unit of SII, but Table 3 and the Results give it per 100 units (SII/100), the scale on which Models 1 and 2 are reproduced."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's income-to-poverty rows 1.5-3.5 and >3.5 repeat its BMI rows 25-30 and >30 in both columns (27.08 and 62.77 with diabetes, 32.94 and 38.99 without), so the income percentages with diabetes sum to 114.5.")
  ),
  choices = list(
    list(choice = "cycles", decision = "2017-March 2020 pre-pandemic files only", reason = "the flow chart starts from P_DEMO's 15,560 participants"),
    list(choice = "diabetes", decision = "told by a doctor (DIQ010 = 1) against not (2); borderline and don't know excluded",
         reason = "the flow chart's 244 with missing diabetes data are exactly the 239 borderline and 5 don't know answers; this gives the published 7,877 and 1,266 cases, and the published mean SII (524.91, SD 358.90). Counting borderline as no changes the estimate by under 0.001 (variant)"),
    list(choice = "weights", decision = "unweighted, with the examination weight (WTMECPRP) as the paper's weight",
         reason = "Table 1's weighted means match WTMECPRP exactly (age 61.09 and 46.55, SII 597.58 and 532.32), but only unweighted fits reproduce the crude and Model 2 estimates with their CIs (1.037, 1.022-1.053 and 1.036, 1.020-1.053 against 1.04, 1.02-1.05 each; weighted 1.044, 1.012-1.077 and 1.031, 0.997-1.066) and, within 0.01, the quartile rows of Models 1 and 2, as EmpowerStats fits them"),
    list(choice = "Model 3 covariates", decision = "the Table 3 footnote's list: age, sex, race, marital status, income, education, drink, smoke, BMI, high blood pressure, BUN, WBC, platelets, chloride, day-1 protein, carbohydrate, total sugars, total fat, cholesterol, vigorous work activity",
         reason = "it describes the table that holds the headline, and its names match Table 1's rows. The Covariates paragraph and Table 1 have serum albumin in place of drink, and the Methods' Model 3 sentence has serum and urine albumin and no WBC, platelets, or drink (variants: 1.019 and 1.020). The footnote's list is also the closest to the published estimate"),
    list(choice = "drink", decision = "ever had a drink of alcohol (ALQ111), yes or no", reason = "unstated; it is the yes/no question parallel to smoke (SMQ020, 100 cigarettes in life); leaving drink out changes the estimate by under 0.001"),
    list(choice = "smoke, high blood pressure, vigorous work activity", decision = "SMQ020, BPQ020, PAQ605, each yes against any other answer (refused and don't know as no)",
         reason = "as the paper names them; coded so, Table 1's weighted percentages are reproduced exactly (51.14, 68.85, 20.53 with diabetes; 41.60, 27.14, 28.05 without)"),
    list(choice = "categorical covariates", decision = "race in RIDRETH1's five groups, marital status in DMDMARTZ's three, education in three (DMDEDUC2 1-2, 3, 4-5), income to poverty 1.5 or less, over 1.5 to 3.5, over 3.5, BMI under 25, 25 to under 30, 30 or more",
         reason = "the groups Table 1 shows; its weighted percentages for BMI are reproduced exactly and for education and marital status within 0.05 points"),
    list(choice = "continuous covariates", decision = "age, BUN (mmol/L), chloride, WBC, platelets, and day-1 dietary intakes (reliable recalls) as linear terms; BUN and chloride as measured", reason = "Table 1 reports them as means, which match ours for BUN, chloride, and the intakes"),
    list(choice = "missing covariates", decision = "complete cases (5,903 of 7,877)",
         reason = "unstated, and the flow chart excludes no one for covariates; filling them in (missing groups for categorical covariates, medians for continuous ones, all 7,877) gives 1.025. No covariate set or missing-data handling tried reaches the published 1.04 (p = 0.0006); Models 1 and 2 are reproduced exactly")
  )
)