# Systemic immune-inflammation index (SII) tertiles and self-reported stroke, adults, NHANES
# 1999-2020. Liu, Qian, Wang, and Wu (2024), Rev Cardiovasc Med, doi:10.31083/j.rcm2504130.
# Headline: OR 1.18 (95% CI 1.01-1.42), highest SII tertile (597 or more) against the lowest
# (under 384), Model 3 (Table 2).

# Leisure-time and transport activity in 1999-2006, which asked other questions than the GPAQ of
# 2007 on. Leisure activities come from the individual activity file (PAQIAF: times in the past 30
# days, minutes each time, and the activity's MET value), walking or bicycling to get places from
# PAQ050Q/U and PAD080 at 4 METs, as the GPAQ scores transport. `met` is MET-minutes a week over
# leisure and transport; `met_leisure` counts leisure only, at 4 METs for moderate and 8 for
# vigorous activities as the GPAQ and the 2021-2023 questions are scored.
row028_activity_early <- function(cycle) {
  p <- component("PAQ", cycle)
  iaf <- component("PAQIAF", cycle)
  weekly_minutes <- iaf$PADTIMES * iaf$PADDURAT * 7 / 30
  per_person <- function(x) {
    total <- tapply(x, iaf$SEQN, sum)
    unname(total[as.character(p$SEQN)])
  }
  # No leisure activity: answered no (or unable) to both the vigorous and the moderate question.
  none <- p$PAD200 %in% 2:3 & p$PAD320 %in% 2:3
  leisure <- function(x) ifelse(p$SEQN %in% iaf$SEQN, per_person(x), ifelse(none, 0, NA))
  level_met <- ifelse(iaf$PADLEVEL %in% 2, 8, ifelse(iaf$PADLEVEL %in% 1, 4, NA))
  codes <- function(x) ifelse(x %in% c(77777, 99999), NA, x)
  times_a_week <- codes(p$PAQ050Q) * ifelse(p$PAQ050U %in% 1, 7, ifelse(p$PAQ050U %in% 2, 1, ifelse(p$PAQ050U %in% 3, 7 / 30, NA)))
  transport <- ifelse(p$PAD020 %in% 2:3, 0, ifelse(p$PAD020 %in% 1, 4 * times_a_week * codes(p$PAD080), NA))
  data.frame(SEQN = p$SEQN, met = leisure(iaf$PADMETS * weekly_minutes) + transport,
             met_leisure = leisure(level_met * weekly_minutes))
}

# MET-minutes a week in any cycle: the GPAQ's work, transport, and leisure domains from 2007 on
# (none in 2021-2023, which asks only about leisure time), and leisure time alone in every cycle.
row028_activity <- function(cycle) {
  if (cycle %in% c("1999-2000", "2001-2002", "2003-2004", "2005-2006")) return(row028_activity_early(cycle))
  leisure <- leisure_activity(cycle)
  total <- met_minutes(cycle)
  data.frame(SEQN = leisure$SEQN, met = total$met[match(leisure$SEQN, total$SEQN)],
             met_leisure = 4 * leisure$leisure_moderate + 8 * leisure$leisure_vigorous)
}

# One cycle's files, as the paper builds each variable.
row028_frame <- function(cycle) {
  demo <- demographics(cycle)
  bpq <- component("BPQ", cycle)
  # Triglycerides and LDL count only for those NCHS counts as fasted (a positive fasting-subsample
  # weight); the files also hold values for some who weren't.
  lipids <- triglycerides(cycle)
  lipids$tg[!(lipids$fasted %in% 1)] <- NA
  lipids$ldl[!(lipids$fasted %in% 1)] <- NA
  d <- merge_all(demo, blood_count(cycle)[, c("SEQN", "platelets", "neutrophils", "lymphocytes")],
                 component("MCQ", cycle, c("MCQ160F", "MCQ220")), smoking(cycle), body_measures(cycle)[, c("SEQN", "bmi")],
                 diabetes_status(cycle), hypertension_status(cycle, 140, 90), total_cholesterol(cycle), hdl_cholesterol(cycle),
                 lipids[, c("SEQN", "tg", "ldl")], bp_questions(cycle)[, c("SEQN", "cholesterol_medication")], row028_activity(cycle))
  d$sii <- sii(d$platelets, d$neutrophils, d$lymphocytes)
  d$stroke <- ifelse(d$MCQ160F %in% 1, 1, ifelse(d$MCQ160F %in% 2, 0, NA))
  d$cancer <- ifelse(d$MCQ220 %in% 1, 1, ifelse(d$MCQ220 %in% 2, 0, NA))
  # Hyperlipidemia as the paper defines it (HDL at or under 40 mg/dL in men, 50 in women). Someone
  # not asked about cholesterol medicine (not told of high cholesterol) is not taking it, so the
  # status is missing only without any lipid value and without an answer about medicine.
  low_hdl <- (d$sex == 1 & d$hdl <= 40) | (d$sex == 2 & d$hdl <= 50)
  criterion <- (d$tg >= 150) %in% TRUE | (d$tc >= 200) %in% TRUE | (d$ldl >= 130) %in% TRUE | low_hdl %in% TRUE |
    d$cholesterol_medication %in% 1
  medicine <- if (has(bpq, "BPQ101D")) bpq$BPQ101D else bpq$BPQ100D
  medicine_unknown <- medicine %in% c(7, 9) | (if (has(bpq, "BPQ090D")) bpq$BPQ090D %in% c(7, 9) else FALSE)
  medicine_known <- d$SEQN %in% bpq$SEQN[!medicine_unknown]
  labs_known <- !is.na(d$tg) | !is.na(d$tc) | !is.na(d$ldl) | !is.na(d$hdl)
  d$hyperlipidemia <- ifelse(criterion, 1, ifelse(labs_known | medicine_known, 0, NA))
  d
}

# The analysis frame for a cycle. `sources` lets a variant build the 2017-March 2020 slot from
# other files: the 2017-2018 files instead, or both, stacked as the paper did. The pre-pandemic
# files renumbered everyone (SEQN 109263 on), so the 2017-2018 participants they contain don't
# collide with DEMO_J's SEQNs; rows from another cycle's files carry that cycle's own examination
# weight in source_weight, which derive puts on the pooled scale.
row028_build <- function(cycle, sources = cycle) {
  frames <- lapply(sources, function(source) {
    d <- row028_frame(source)
    d$source_cycle <- source
    d$source_weight <- NA_real_
    if (source != cycle) {
      w <- weights_for("WTMEC2YR", source, blood_file = "CBC")
      d$source_weight <- w$weight_2yr[match(d$SEQN, w$SEQN)]
    }
    d
  })
  d <- do.call(rbind, frames)
  # The paper's flow chart starts from 66,568 participants aged over 18 (19 and older: the count
  # matches exactly); the stroke question is asked from age 20, so the analytic sample is 20 and
  # older either way.
  d$in_population <- d$age >= 18 & !is.na(d$sii) & !is.na(d$stroke) & !is.na(d$smoking) &
    !is.na(d$hyperlipidemia) & !is.na(d$hypertension)
  d
}

# Medians for imputing covariates, over the paper cycles' study population: the paper imputes every
# covariate's missing values with its median, except income, whose missing values are a category.
row028_constants <- function(data) {
  p <- data$in_population %in% TRUE
  middle <- function(x) unname(stats::quantile(x[p], 0.5, type = 1, na.rm = TRUE))
  list(bmi = middle(data$bmi), met = middle(data$met), met_leisure = middle(data$met_leisure),
       diabetes = middle(data$diabetes), cancer = middle(data$cancer))
}

row028_activity_groups <- function(met) {
  factor(ifelse(met <= 0, "sedentary", ifelse(met < 500, "insufficient", ifelse(met <= 1000, "moderate", "high"))),
         levels = c("sedentary", "insufficient", "moderate", "high"))
}

row028_derive <- function(data, constants) {
  if (is.null(constants)) constants <- row028_constants(data)
  # Rows a variant took from another cycle's files: their own 2-year weight, on the scale
  # pooled_weight put the other cycles on (each weight times its years over the pooled years).
  borrowed <- !is.na(data$source_weight)
  if (any(borrowed)) {
    total <- sum(vapply(unique(data$cycle), function(cycle) cycle_info(cycle)$years, 0))
    years <- vapply(data$source_cycle[borrowed], function(cycle) cycle_info(cycle)$years, 0)
    data$w[borrowed] <- data$source_weight[borrowed] * years / total
  }
  impute <- function(x, value) ifelse(is.na(x), value, x)
  data$sii_group <- cut(data$sii, c(-Inf, 384, 597, Inf), right = FALSE, labels = c("low", "median", "high"))
  data$log10_sii <- log10(data$sii)
  data$sex <- factor(data$sex)
  data$race4 <- factor(ifelse(data$race == 3, "white", ifelse(data$race == 4, "black", ifelse(data$race == 1, "mexican", "other"))),
                       levels = c("white", "black", "mexican", "other"))
  data$smoking3 <- factor(data$smoking, levels = 1:3, labels = c("never", "former", "current"))
  # Table 1's counts put the few with no education answer under high school.
  data$education3 <- factor(impute(education3(data$education), 1), levels = 1:3,
                            labels = c("under high school", "high school", "college or higher"))
  data$pir4 <- factor(ifelse(is.na(data$pir), "unknown", ifelse(data$pir <= 1, "<=1", ifelse(data$pir <= 3, "1-3", ">3"))),
                      levels = c("<=1", "1-3", ">3", "unknown"))
  data$bmi3 <- cut(impute(data$bmi, constants$bmi), c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-29.9", ">=30"))
  data$activity4 <- row028_activity_groups(impute(data$met, constants$met))
  data$activity4_leisure <- row028_activity_groups(impute(data$met_leisure, constants$met_leisure))
  data$diabetes2 <- factor(impute(data$diabetes, constants$diabetes))
  data$cancer2 <- factor(impute(data$cancer, constants$cancer))
  data$hyperlipidemia2 <- factor(data$hyperlipidemia)
  data$hypertension2 <- factor(data$hypertension)
  # Complete-case versions (no imputation) for a variant.
  data$bmi3_cc <- cut(data$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-29.9", ">=30"))
  data$activity4_cc <- row028_activity_groups(data$met)
  data$education3_cc <- factor(education3(data$education))
  data$diabetes2_cc <- factor(data$diabetes)
  data$cancer2_cc <- factor(data$cancer)
  data
}

row028_model1 <- stroke ~ sii_group + age + sex + race4
row028_model2 <- stroke ~ sii_group + age + sex + race4 + smoking3 + activity4 + education3 + pir4 + bmi3

association <- list(
  id = "row028", row = 28, doi = "10.31083/j.rcm2504130",
  # The paper says 1999 to 2020. Its 66,568 participants aged over 18 are exactly the count over
  # the ten 2-year files from 1999-2000 to 2017-2018 plus the 2017-March 2020 pre-pandemic files,
  # so it stacked the 2017-2018 participants twice (the pre-pandemic files renumber them). The
  # primary analysis uses 1999-2016 and 2017-March 2020, which cover those years once; variants
  # stack them as the paper did, or use 2017-2018 alone.
  cycles = c("1999-2000", "2001-2002", "2003-2004", "2005-2006", "2007-2008", "2009-2010", "2011-2012",
             "2013-2014", "2015-2016", "2017-2020"),
  weight = "WTMEC2YR", blood_file = "CBC",
  family = "logistic", term = "sii_grouphigh",
  published = list(measure = "OR", estimate = 1.18, low = 1.01, high = 1.42, n = 57600, events = 2368,
                   contrast = "highest SII tertile (597 or more) against the lowest (under 384); SII = platelets x neutrophils / lymphocytes, counts in 10^3 cells/uL"),
  left_out = c("physical activity over work, transport, and leisure: 2021-2023 asks only about leisure time, so the harmonized version groups leisure-time MET-minutes (4 METs moderate, 8 vigorous) into the paper's four groups in every cycle"),
  build = function(cycle) row028_build(cycle),
  constants = row028_constants,
  derive = row028_derive,
  formula = stroke ~ sii_group + age + sex + race4 + smoking3 + activity4 + education3 + pir4 + bmi3 +
    diabetes2 + hyperlipidemia2 + cancer2 + hypertension2,
  formula_harmonized = stroke ~ sii_group + age + sex + race4 + smoking3 + activity4_leisure + education3 + pir4 + bmi3 +
    diabetes2 + hyperlipidemia2 + cancer2 + hypertension2,
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "crude (published 1.46, 1.23-1.74)", formula = stroke ~ sii_group),
    list(label = "Model 1 (published 1.32, 1.11-1.58)", formula = row028_model1),
    list(label = "Model 2 (published 1.22, 1.02-1.47)", formula = row028_model2),
    list(label = "2017-2018 stacked with 2017-March 2020, as the paper did",
         build = function(cycle) row028_build(cycle, if (cycle == "2017-2020") c("2017-2018", "2017-2020") else cycle)),
    list(label = "2017-2018 files in place of 2017-March 2020",
         build = function(cycle) row028_build(cycle, if (cycle == "2017-2020") "2017-2018" else cycle)),
    list(label = "Model 3 without cancer (the Methods' list)",
         formula = stroke ~ sii_group + age + sex + race4 + smoking3 + activity4 + education3 + pir4 + bmi3 + diabetes2 + hyperlipidemia2 + hypertension2),
    list(label = "complete cases, no imputation", sample = "own",
         formula = stroke ~ sii_group + age + sex + race4 + smoking3 + activity4_cc + education3_cc + pir4 + bmi3_cc +
           diabetes2_cc + hyperlipidemia2 + cancer2_cc + hypertension2),
    list(label = "per unit of log10 SII (published 1.30, 0.99-1.70)", term = "log10_sii",
         formula = stroke ~ log10_sii + age + sex + race4 + smoking3 + activity4 + education3 + pir4 + bmi3 +
           diabetes2 + hyperlipidemia2 + cancer2 + hypertension2)
  ),
  # 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, evidence = "data",
         detail = "The text describes adults in NHANES 1999 to 2020, but its 66,568 adults and 59,732 with SII are exactly the counts over the ten 2-year files to 2017-2018 plus the 2017-March 2020 files, which hold the 2017-2018 participants a second time under new SEQNs. Stacked that way, the sample here is 57,606 (published 57,600) and Model 3 gives 1.18 (0.99 to 1.42)."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "Model 3's interval for the highest tertile, 1.18 (1.01 to 1.42), is asymmetric on the log scale beyond rounding (0.16 below the estimate, 0.19 above), unlike Table 2's other intervals: its estimate and upper bound imply a lower bound of 0.97 to 0.99, and the stacked sample here gives 1.18 (0.99 to 1.42)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The text imputes each covariate's missing values with its median, which for education is college or higher, but Table 1's education counts (14,960, 13,379, 29,261) put those with no answer under high school, as coding them so does in the stacked sample here (14,963, 13,379, 29,264)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The text and Table 1 give 2,368 strokes among the 57,600 analyzed, more than the 2,310 adults with SII who answered yes to the stroke question in the stacked files; the stacked sample here has 2,305."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1 labels the sexes the other way round: its 29,870 males (51.86%, repeated in the Results) are the stacked sample's 29,873 women."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's note calls its percentages weighted, but the Total column's are unweighted shares: its 4.11% with stroke is 2,368 of 57,600, above every tertile's weighted share (2.58% to 3.72%), and the stacked sample's weighted share is 2.94%."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's race and ethnicity rows count 2,986 participants in the Total column (1,707, 668, 456, and 155), not 57,600."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 3 gives 1,757 strokes among the 37,767 under 60 and 611 among the 19,833 aged 60 or more, where the stacked sample has 587 and 1,718, and 4,871 strokes among the 52,200 without cancer, more than the 2,368 in all (1,832 in the stacked sample)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The text's 66,568 adults 'aged 18 or above' are exactly those aged 19 or older in the stacked files; the analytic sample is 20 and older either way, as the stroke question is asked from age 20."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract gives the sample as 53,600 people and the Discussion as 53,111, against the 57,600 of the Methods and Table 1."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract says the risk of stroke 'decreased by 34%' for every unit of log SII and cites Model 3's OR of 1.30 (0.99 to 1.70), which is an increase; the 34% is Model 2's (1.34).")
  ),
  choices = list(
    list(choice = "cycles", decision = "1999-2000 to 2015-2016 and 2017-March 2020 (each participant once); variants stack 2017-2018 with 2017-March 2020 as the paper did, or use 2017-2018 alone",
         reason = "the paper says only 1999 to 2020; its 66,568 participants over 18 and 59,732 with SII are exactly the counts over the ten 2-year files to 2017-2018 plus the pre-pandemic files, which hold the 2017-2018 participants again under new SEQNs. Stacked, our sample is 57,606 (published 57,600) and Model 3 gives 1.18"),
    list(choice = "weights", decision = "examination weight (WTMEC2YR, the 4-year WTMEC4YR for 1999-2002, WTMECPRP for 2017-March 2020; WTPH2YR in 2021-2023), pooled by years as NCHS directs, with strata and PSUs",
         reason = "the paper says it used appropriate weights with the R survey package; the weighted crude OR implied by Table 1's weighted stroke percentages (1.46) is the published crude OR, and weighted fits come closest (unweighted gives 1.29)"),
    list(choice = "age", decision = "18 and older, as stated", reason = "the flow chart's 66,568 counts those over 18; the stroke question is asked from age 20, so the analytic sample is 20 and older either way"),
    list(choice = "stroke", decision = "MCQ160F yes or no; refused and don't know are missing and excluded", reason = "this excludes 2,062 for missing stroke (flow chart 2,059)"),
    list(choice = "SII tertiles", decision = "the published cutpoints: under 384, 384 to under 597, 597 or more", reason = "stated; they are the unweighted tertiles of the stacked sample (384.0 and 596.7)"),
    list(choice = "income to poverty ratio", decision = "1 or less, over 1 to 3, over 3, and missing as unknown", reason = "the Methods' groups; on the stacked sample they reproduce Table 1's counts (10,721, 22,022, 19,444, 5,419 against 10,721, 22,019, 19,443, 5,417), which Table 1 labels <1, 1 to <3, 3 or more"),
    list(choice = "education", decision = "under high school (DMDEDUC2 1-2), high school (3), college or higher (4-5); those with no answer (61) go under high school", reason = "on the stacked sample this reproduces Table 1's counts (14,963, 13,379, 29,264 against 14,960, 13,379, 29,261); the median would put them in college or higher"),
    list(choice = "missing covariates", decision = "BMI, MET-minutes, diabetes, and cancer imputed with the study population's median (BMI 27.9, so the 25-29.9 group), income missing as its own group", reason = "the paper imputes each covariate's median and keeps unknown income as a category; Table 1's BMI counts match"),
    list(choice = "physical activity", decision = "MET-minutes a week in the paper's four groups: 2007-March 2020 GPAQ work, transport, and leisure (8 METs vigorous, 4 moderate, 4 walking or cycling); 1999-2006 leisure activities from PAQIAF at NCHS's MET values plus walking or cycling to get places at 4 METs",
         reason = "the paper gives the groups but not the items; this uses the GPAQ's domains that 1999-2006 asked in minutes. No construction tried matches Table 1's groups (published 27.7%, 20.6%, 11.4%, 40.2%; ours 29.7%, 15.0%, 11.2%, 44.1%); dropping activity moves the estimate from 1.16 to 1.19"),
    list(choice = "diabetes", decision = "told by a doctor, insulin or diabetes pills, fasting glucose 126 mg/dL or more, or HbA1c 6.5% or more", reason = "stated; on the stacked sample Table 1 has 10,335 against our 9,939 (counting borderline gives 10,716); either changes the estimate by under 0.01"),
    list(choice = "hyperlipidemia", decision = "triglycerides 150 mg/dL or more, total cholesterol 200 or more, LDL 130 or more, HDL 40 or less (men) or 50 or less (women), or now taking cholesterol medicine (BPQ100D; BPQ101D in 2021-2023); not asked counts as not taking",
         reason = "stated; the flow chart excludes only 2 for missing hyperlipidemia, which counting unasked as not taking reproduces (1), while treating it as missing would exclude 687. Triglycerides and LDL only from those who fasted, as NCHS directs: on the stacked sample this gives 41,012 with hyperlipidemia (Table 1: 40,996; 41,101 with every value in the files)"),
    list(choice = "hypertension", decision = "told by a doctor, taking medicine for it, or mean measured pressure 140/90 or more", reason = "unstated; the flow chart's 24 missing fits a combined definition (13 missing) rather than self-report alone (205); the estimate changes by under 0.01 with self-report alone"),
    list(choice = "cancer", decision = "in Model 3 (MCQ220)", reason = "Table 2's note lists it; the Methods' Model 3 sentence omits it (variant)"),
    list(choice = "log base for the continuous SII", decision = "log10 (variant only)", reason = "unstated; per unit of log10 SII, Model 3 gives 1.29 (1.29 stacked, 1.30 published, 0.99-1.70), per unit of the natural log 1.12")
  )
)