# Systemic immune-inflammation index (SII) and MASLD (fatty liver index 60 or more, or US fatty
# liver index 30 or more) in adults aged 20 and older, NHANES 2007-2018. Wang et al. (2024),
# Front Nutr, doi:10.3389/fnut.2024.1415484.
# Headline: OR 1.47 (95% CI 1.24-1.74), SII quartile 4 (625.74-28,397.28) vs quartile 1
# (1.53-313.50), Model 3 (Figure 3; Supplementary Table S4(A)).

# Fasting insulin (uU/mL): in the glucose file (GLU) through 2011-2012, in its own file (INS) from
# 2013-2014 on (2021-2023 included).
fasting_insulin <- function(cycle) {
  g <- component("GLU", cycle)
  d <- if (has(g, "LBXIN")) g else component("INS", cycle)
  data.frame(SEQN = d$SEQN, insulin = d$LBXIN)
}

# Viral hepatitis markers, 1 if positive: hepatitis B surface antigen (LBDHBG), hepatitis C
# antibody (confirmed), and HCV RNA (LBXHCR). The confirmed antibody is LBDHCV through 2011-2012
# and LBDHCI from 2017-2018 on; 2013-2016 released only the RNA result (tested when the antibody
# screen was reactive). The library's hepatitis() reads only LBDHCI, so it has no antibody result
# for 2007-2012.
viral_markers <- function(cycle) {
  b <- component("HEPBD", cycle)
  c <- component("HEPC", cycle)
  antibody <- if (has(c, "LBDHCI")) c$LBDHCI else if (has(c, "LBDHCV")) c$LBDHCV else rep(NA, nrow(c))
  out <- data.frame(SEQN = union(b$SEQN, c$SEQN))
  out$hbsag <- b$LBDHBG[match(out$SEQN, b$SEQN)]
  out$hcv_antibody <- antibody[match(out$SEQN, c$SEQN)]
  out$hcv_rna <- c$LBXHCR[match(out$SEQN, c$SEQN)]
  out
}

# Prescription medicines that can cause steatosis, taken in the past 30 days (RXQ_RX), by Multum
# drug ID: amiodarone, methotrexate, tamoxifen, valproate (valproic acid, divalproex), and systemic
# glucocorticoids (prednisone, prednisolone, methylprednisolone, dexamethasone, hydrocortisone,
# cortisone; topical, eye, and ear forms have IDs of their own and are not counted). 1 if any was
# taken, 0 if the prescription question was answered. 2021-2023 releases no drug names.
STEATOGENIC_DRUGS <- c("d00002", "d00060", "d00381", "d00083", "d03833",
                       "d00350", "d00084", "d00293", "d00206", "d00254", "d00609")
steatogenic_medication <- function(cycle) {
  rx <- component("RXQ_RX", cycle)
  people <- data.frame(SEQN = unique(rx$SEQN))
  users <- rx$SEQN[rx$RXDDRGID %in% STEATOGENIC_DRUGS]
  answered <- rx$SEQN[rx$RXDUSE %in% 1:2]
  people$steatogenic <- ifelse(people$SEQN %in% users, 1, ifelse(people$SEQN %in% answered, 0, NA))
  people
}

# Moderate-to-vigorous activity in minutes a week by domain: minutes a day times days a week, with
# vigorous minutes counted twice. 2007-2018 asked the GPAQ about work, walking or cycling to get
# places, and leisure; 2021-2023 asks only about leisure time (leisure_activity() reads both).
# The "_printed" sums count vigorous minutes once, as the paper's formula is printed.
activity_minutes <- function(cycle) {
  leisure <- leisure_activity(cycle)
  out <- data.frame(SEQN = leisure$SEQN, leisure = leisure$mvpa_equivalent,
                    leisure_printed = leisure$leisure_moderate + leisure$leisure_vigorous,
                    work = NA_real_, work_printed = NA_real_, transport = NA_real_)
  d <- component("PAQ", cycle)
  if (has(d, "PAQ605")) {
    part <- function(answer, days, minutes) {
      m <- ifelse(minutes %in% c(7777, 9999), NA, minutes)
      n <- ifelse(days %in% 1:7, days, NA)
      ifelse(answer %in% 2, 0, ifelse(answer %in% 1, n * m, NA))
    }
    rows <- match(out$SEQN, d$SEQN)
    vigorous <- part(d$PAQ605, d$PAQ610, d$PAD615)
    moderate <- part(d$PAQ620, d$PAQ625, d$PAD630)
    out$work <- (2 * vigorous + moderate)[rows]
    out$work_printed <- (vigorous + moderate)[rows]
    out$transport <- part(d$PAQ635, d$PAQ640, d$PAD645)[rows]
  }
  out
}

# Fatty liver index (Bedogni et al. 2006): triglycerides mg/dL, GGT U/L, waist cm.
fatty_liver_index <- function(tg, bmi, ggt, waist) {
  100 * stats::plogis(0.953 * log(tg) + 0.139 * bmi + 0.718 * log(ggt) + 0.053 * waist - 15.745)
}

# US fatty liver index (Ruhl and Everhart 2015): GGT U/L, waist cm, insulin uU/mL, glucose mg/dL.
# The published index takes ln(insulin); the paper prints insulin untransformed (log_insulin = FALSE).
us_fatty_liver_index <- function(race, age, ggt, waist, insulin, glucose, log_insulin = TRUE) {
  100 * stats::plogis(0.3458 * (race == 1) - 0.8073 * (race == 4) + 0.0093 * age + 0.6151 * log(ggt) + 0.0249 * waist +
                        1.1792 * (if (log_insulin) log(insulin) else insulin) + 0.8242 * log(glucose) - 14.7812)
}

# The paper's models (Statistical analysis; Supplementary Tables S2-S4).
MODEL_2 <- masld ~ sii_q + sex + age + race + pir_poor + education + marital + insurance
MODEL_3 <- update(MODEL_2, . ~ . + tobacco + alcohol + hypertension + diabetes + cvd + waist_high + active + bmi3 +
                    tg_high + hdl_high + alt + ast + ggt)

# The MASLD definition's exclusions as the text states them, with those excluded dropped.
stated_exclusions <- function(data, constants) {
  data$sii_q <- sii_quartile(data$sii)
  data$in_population <- data$in_population & data$other_liver_disease %in% 0
  data
}

# SII quartiles at the paper's cutpoints (Figure 3 legend): each boundary lies between the printed
# end of one quartile and the printed start of the next.
sii_quartile <- function(x) cut(x, c(-Inf, 313.505, 440.005, 625.735, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))

association <- list(
  id = "row327", row = 327, doi = "10.3389/fnut.2024.1415484",
  cycles = c("2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  # USFLI needs fasting insulin and glucose and FLI fasting triglycerides, so a weighted analysis
  # takes the fasting subsample weight (in GLU, GLU_L in 2021-2023).
  weight = "WTSAF2YR", blood_file = "CBC",
  weighted = FALSE,
  family = "logistic", term = "sii_qQ4",
  published = list(measure = "OR", estimate = 1.47, low = 1.24, high = 1.74, n = 14413, events = 6518,
                   contrast = "SII quartile 4 (625.74-28,397.28) vs quartile 1 (1.53-313.50); SII = platelets x neutrophils / lymphocytes (10^3 cells/uL)"),
  left_out = c(
    "physical activity at work and walking or cycling to get places (GPAQ, PAQ605-PAD645): 2021-2023 asks only about leisure time; replaced by leisure-time moderate-to-vigorous activity of 150 minutes a week or more, vigorous minutes counted twice"
  ),
  build = function(cycle) {
    demo <- demographics(cycle)
    alq <- alcohol(cycle)
    bio <- biochemistry(cycle)
    d <- merge_all(demo, blood_count(cycle)[, c("SEQN", "platelets", "neutrophils", "lymphocytes")],
                   body_measures(cycle)[, c("SEQN", "bmi", "waist")], triglycerides(cycle)[, c("SEQN", "tg")],
                   fasting_glucose(cycle), fasting_insulin(cycle), hdl_cholesterol(cycle),
                   component("HIQ", cycle, "HIQ011"), component("SMQ", cycle, "SMQ020"),
                   hypertension_status(cycle), diabetes_status(cycle), activity_minutes(cycle),
                   viral_markers(cycle), liver_cancer(cycle))
    rows <- match(d$SEQN, bio$SEQN)
    d$alt <- bio$alt[rows]
    d$ast <- bio$ast[rows]
    d$ggt <- bio$ggt[rows]
    conditions <- medical_conditions(cycle)
    d$cvd <- cvd(conditions)[match(d$SEQN, conditions$SEQN)]
    # Autoimmune hepatitis (MCQ510E, asked from 2017-2018 on).
    mcq <- component("MCQ", cycle)
    d$autoimmune_hepatitis <- if (has(mcq, "MCQ510E")) as.integer(mcq$MCQ510E[match(d$SEQN, mcq$SEQN)] %in% 5) else 0
    # Two-hour OGTT glucose: given 2005-2006 through 2015-2016 only.
    has_ogtt <- cycle %in% c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016")
    ogtt <- if (has_ogtt) component("OGTT", cycle, "LBXGLT") else data.frame(SEQN = d$SEQN, LBXGLT = NA_real_)
    d$ogtt <- ogtt$LBXGLT[match(d$SEQN, ogtt$SEQN)]
    if (cycle != REPLICATION_CYCLE) {
      d <- merge_all(d, steatogenic_medication(cycle))
    } else {
      d$steatogenic <- NA
    }
    a <- alq[match(d$SEQN, alq$SEQN), ]
    # Drinks in the past 12 months: drinks on a drinking day times drinking days; none for those who
    # never drank or did not drink in the past year.
    drinks_year <- ifelse(a$alcohol3 %in% 1:2, 0, a$drinks_per_day * a$drinking_days)

    d$sii <- sii(d$platelets, d$neutrophils, d$lymphocytes)
    d$fli <- fatty_liver_index(d$tg, d$bmi, d$ggt, d$waist)
    d$usfli <- us_fatty_liver_index(d$race, d$age, d$ggt, d$waist, d$insulin, d$glucose)
    printed <- us_fatty_liver_index(d$race, d$age, d$ggt, d$waist, d$insulin, d$glucose, log_insulin = FALSE)
    # Steatosis by either index; someone with only one index is classified by it.
    d$masld <- as.integer((d$fli >= 60) %in% TRUE | (d$usfli >= 30) %in% TRUE)
    d$masld_printed <- as.integer((d$fli >= 60) %in% TRUE | (printed >= 30) %in% TRUE)
    # The definition's exclusions: hepatitis B or C, more than 2 drinks a day (men) or 1 (women)
    # averaged over the past year, steatogenic medicines, liver cancer, autoimmune hepatitis.
    heavy <- (d$sex == 1 & drinks_year / 365 > 2) | (d$sex == 2 & drinks_year / 365 > 1)
    d$other_liver_disease <- as.integer(d$hbsag %in% 1 | d$hcv_antibody %in% 1 | d$hcv_rna %in% 1 | heavy %in% TRUE |
                                          d$steatogenic %in% 1 | d$liver_cancer %in% 1 | d$autoimmune_hepatitis %in% 1)
    d$masld_stated <- as.integer(d$masld == 1 & d$other_liver_disease == 0)

    d$sex <- factor(d$sex)
    d$race <- factor(d$race)
    d$pir_poor <- factor(ifelse(d$pir < 1.3, "poor", ifelse(d$pir >= 1.3, "not poor", NA)), levels = c("not poor", "poor"))
    d$education <- factor(ifelse(d$education %in% 1:2, "below high school", ifelse(d$education %in% 3:5, "high school or above", NA)),
                          levels = c("high school or above", "below high school"))
    d$marital <- factor(d$marital)
    d$insurance <- factor(ifelse(d$HIQ011 %in% 1, "yes", ifelse(d$HIQ011 %in% 2, "no", NA)), levels = c("no", "yes"))
    d$tobacco <- factor(ifelse(d$SMQ020 %in% 1, "yes", ifelse(d$SMQ020 %in% 2, "no", NA)), levels = c("no", "yes"))
    d$alcohol <- factor(ifelse(drinks_year >= 12, "yes", ifelse(drinks_year < 12, "no", NA)), levels = c("no", "yes"))
    d$hypertension <- factor(d$hypertension)
    d$diabetes_ogtt <- factor(ifelse(d$diabetes %in% 1 | (d$ogtt >= 200) %in% TRUE, 1, d$diabetes))
    d$diabetes <- factor(d$diabetes)
    d$cvd <- factor(d$cvd)
    d$waist_high <- factor(as.integer((d$sex == 1 & d$waist >= 102) | (d$sex == 2 & d$waist >= 88)))
    # 150 minutes a week or more; no answer counts as none.
    d$active <- factor(as.integer(rowSums(cbind(d$work, d$transport, d$leisure), na.rm = TRUE) >= 150))
    d$active_printed <- factor(as.integer(rowSums(cbind(d$work_printed, d$transport, d$leisure_printed), na.rm = TRUE) >= 150))
    d$active_leisure <- factor(as.integer((d$leisure >= 150) %in% TRUE))
    d$bmi3 <- cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-29.9", ">=30"))
    # Fasting triglycerides and HDL cholesterol, each split where Table 1's counts put it; no
    # fasting triglyceride counts as not high, as Table 1 has none missing.
    d$tg_high <- factor(as.integer((d$tg >= 160) %in% TRUE))
    d$tg_150 <- factor(as.integer((d$tg >= 150) %in% TRUE))
    d$hdl_high <- factor(as.integer(d$hdl >= 90))
    d$alcohol_unclear <- with_unclear(as.character(d$alcohol), c("no", "yes"))
    d$pir_unclear <- with_unclear(as.character(d$pir_poor), c("not poor", "poor"))
    # Adults 20 and older with an SII and at least one of the two indices.
    d$in_population <- d$age >= 20 & !is.na(d$sii) & (!is.na(d$fli) | !is.na(d$usfli))
    d
  },
  derive = function(data, constants) {
    data$sii_q <- sii_quartile(data$sii)
    data
  },
  formula = MODEL_3,
  formula_harmonized = update(MODEL_3, . ~ . - active + active_leisure),
  variants = list(
    list(label = "crude, all 14,414 (published 1.62, 1.48-1.78)", formula = masld ~ sii_q, sample = "own"),
    list(label = "Model 2, own sample (published 1.63, 1.48-1.80)", formula = MODEL_2, sample = "own"),
    list(label = "weighted (fasting subsample weight)", weighted = TRUE),
    list(label = "MASLD exclusions as stated, excluded dropped", derive = stated_exclusions, sample = "own"),
    list(label = "MASLD exclusions as stated, counted as non-MASLD", formula = update(MODEL_3, masld_stated ~ .)),
    list(label = "USFLI as printed (insulin not logged)", formula = update(MODEL_3, masld_printed ~ .)),
    list(label = "diabetes with 2-h OGTT glucose (text)", formula = update(MODEL_3, . ~ . - diabetes + diabetes_ogtt)),
    list(label = "activity as printed (vigorous counted once)", formula = update(MODEL_3, . ~ . - active + active_printed)),
    list(label = "high triglycerides at 150 mg/dL", formula = update(MODEL_3, . ~ . - tg_high + tg_150)),
    list(label = "alcohol, income missing as a level (n near MI)", formula = update(MODEL_3, . ~ . - alcohol - pir_poor + alcohol_unclear + pir_unclear),
         sample = "own")
  ),
  # 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 = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The text defines MASLD as steatosis (FLI 60 or more, or USFLI 30 or more) without hepatitis B or C, excessive alcohol or drugs, liver cancer, or autoimmune liver disease, but the counts fit steatosis with none of these exclusions: the 14,413 participants are everyone with an index and a blood count (14,414 here; 12,858 with those excluded dropped), and the 6,518 cases exceed the 6,243 with steatosis here (5,512 with those excluded counted as non-cases)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "unresolved",
         detail = "The Methods print USFLI with a 1.1792 x insulin term, but the paper's share with steatosis (6,518 of 14,413, 45.2%) fits the published index's ln(insulin) (43.3%), not the printed term (85.1%)."),
    list(kind = "coding", affects_headline = TRUE, followed = FALSE, evidence = "unresolved",
         detail = "With USFLI as published, the paper's 6,518 cases still exceed the 6,243 that FLI and USFLI give for the same population (14,414 here, with Table 1's ALT, AST, GGT, and SII means matched), the extra cases falling where USFLI decides (182 cases with BMI under 25 in Table 1 against 97 here); no coding tried gives Table 1's race and BMI splits."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The text's diabetes definition includes a 2-hour OGTT glucose of 200 mg/dL or more, but Table 1's count with diabetes (2,806, under a swapped label) is that of the definition without it (2,807 here; 3,118 with it)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The printed MVPA formula adds moderate and vigorous minutes, but Table 1's 8,732 sufficiently active participants are reproduced with vigorous minutes counted twice (8,725; 8,505 counted once)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1 swaps the labels of the diabetes, hypertension, and income rows (its 2,806 without diabetes, 6,007 without hypertension, and 67.38% poor fit those with diabetes, 2,807 here, those with hypertension, 6,109, and those not poor, 67.9%), and its cardiovascular disease row repeats the health insurance row's figures for those without MASLD."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The counts in the text and Figure 1 do not add up: 34,770 less the 20,000 and 57 excluded is 14,713, not 14,413, and 14,413 less 6,518 is 7,895 (Table 1), not 7,985; the text also gives the age exclusion as 2,572 where Figure 1 has 25,072."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract's quartile 4 odds ratios for NLR (1.25, 1.04-1.49) and PLR (1.29, 1.09-1.53) differ from Figure 3's Model 3 estimates (1.29, 1.09-1.52; 1.13, 0.96-1.33)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The text reports SII quartile 4 as 1.62 (1.48, 1.78) in the propensity-matched sample, but that is Figure 3's unadjusted estimate in the full sample; Table S5(A) gives the matched estimate as 1.556 (1.277, 1.897).")
  ),
  choices = list(
    list(choice = "weights", decision = "unweighted",
         reason = "the paper mentions no weights, strata, or PSUs; Table 1 gives unweighted counts and the supplement's models are SPSS's unweighted logistic regression output. Weighting with the fasting subsample weight (WTSAF2YR) is a variant"),
    list(choice = "study population", decision = "adults 20 and older with an SII (platelet, neutrophil, and lymphocyte counts) and at least one of FLI or USFLI",
         reason = "gives 14,414 against the paper's 14,413, with Table 1's ALT, AST, and GGT means (25.02, 25.35, 29.67 against 25.01, 25.35, 29.67) and SII mean and SD (514.04 and 394.52 against 513.93 and 394.28), and 56 lacking a blood count against the paper's 57; requiring both indices gives 14,051. The flow chart's 'missing liver data (n = 20,000)' does not add up (34,770 - 20,000 - 57 = 14,713; 34,770 - 14,413 - 57 = 20,300)"),
    list(choice = "MASLD exclusions", decision = "none applied: the outcome is steatosis by FLI or USFLI",
         reason = "the paper's numbers show the definition's exclusions (hepatitis B or C, heavy drinking, drugs, liver cancer, autoimmune liver disease) were not applied: its n is everyone with index and blood count data (dropping the 1,556 with an exclusion would leave 12,858), and its 6,518 cases exceed the 6,243 with steatosis (counting the excluded as non-cases would leave 5,512); Table 1's drinking shares (55.1% of non-cases, 51.1% of cases) fit no reclassification (ours 54.1% and 49.6%; 55.8% and 46.2% with heavy drinkers made non-cases). The stated definition is a variant both ways; its drug step (our list, since the paper names no drugs: amiodarone, methotrexate, tamoxifen, valproate, systemic glucocorticoids; 339 users) can't be built in 2021-2023 but is not part of the chosen version"),
    list(choice = "USFLI insulin term", decision = "ln(insulin), as Ruhl and Everhart published the index",
         reason = "the paper prints 1.1792 x insulin; computed that way 85.1% of the sample has steatosis (variant), against the paper's 45.2%; with ln(insulin) 43.3% do"),
    list(choice = "index inputs", decision = "fasting triglycerides (LBXTR; LBXTLG in 2021-2023) in mg/dL, GGT in U/L, waist in cm, fasting insulin in uU/mL (GLU file through 2011-2012, INS after), fasting glucose in mg/dL; steatosis if either available index meets its cutoff",
         reason = "the published indices' units; fasting triglycerides give the paper's sample size (serum triglycerides would make FLI available for 29,639 adults)"),
    list(choice = "case count", decision = "kept as the published formulas give it (6,243 cases against 6,518)",
         reason = "the paper's 275 extra cases fall mostly where USFLI decides (Table 1: BMI under 25, 182 cases against our 97; Mexican Americans 1,260 against 1,168; Black participants 1,282 against 1,269), which points to larger USFLI values in part of the data. Codings tried that change the count: FLI from serum triglycerides (6,406), insulin in pmol/L (7,033; in 2009-2012 alone exactly 6,518, but other pairs of cycles come as close and none gives Table 1's race and BMI splits), lower cutoffs; none is used"),
    list(choice = "SII quartiles", decision = "the published cutpoints, boundaries between the printed ranges (313.505, 440.005, 625.735)",
         reason = "Figure 3 legend; quartiles of 3,604, 3,605, 3,601, and 3,604, with the legend's minimum (1.53) and maximum (28,397.28) exactly"),
    list(choice = "diabetes", decision = "told by a doctor, insulin or diabetes pills, HbA1c 6.5% or more, or fasting glucose 126 mg/dL or more; no OGTT",
         reason = "the text also lists 2-hour OGTT glucose of 200 or more, but Table 1's 2,806 with diabetes (under its swapped labels) match the definition without it (2,807) and not with it (3,118); the text's version is a variant"),
    list(choice = "physical activity", decision = "work, walking or cycling, and leisure minutes a week (minutes a day times days), vigorous minutes counted twice, 150 or more as sufficiently active; no answer counts as none",
         reason = "the paper prints moderate plus vigorous minutes without the doubling and names no domains; with all domains and the doubling 8,725 are active against Table 1's 8,732 (8,505 without the doubling, variant; 4,675 with leisure only)"),
    list(choice = "triglycerides and HDL covariates", decision = "fasting triglycerides 160 mg/dL or more; HDL cholesterol 90 mg/dL or more",
         reason = "Table 1 gives only yes and no counts (3,018 and 424, no missing): 160 gives 3,042 (150 gives 3,555, variant) and 90 gives 414 (89 gives 446); Table S4(A)'s HDL OR (0.325) fits HDL high, not low"),
    list(choice = "alcohol use", decision = "12 or more drinks in the past 12 months (drinking days times drinks a day)",
         reason = "the paper's 'more than 12 drinks in the past year'; 52.1% of those with data drink so, against Table 1's 53.26% (more than 12: 47.8%)"),
    list(choice = "hypertension", decision = "told by a doctor, taking medicine for it, or mean pressure 140/90 or more (auscultatory readings; oscillometric in 2021-2023)",
         reason = "the paper's definition; 6,109 against the 6,007 Table 1 shows under its swapped labels"),
    list(choice = "other covariates", decision = "sex; age continuous; RIDRETH1's five groups; income to poverty ratio under 1.3; education below high school or not; marital status in three groups (DMDMARTZ's); health insurance (HIQ011); smoked 100 cigarettes (SMQ020); CVD as told of heart failure, coronary heart disease, angina, heart attack, or stroke; waist 102 cm or more (men) or 88 (women); BMI under 25, 25-29.9, 30 or more; ALT, AST, GGT continuous",
         reason = "the paper's definitions; Table 1's counts are matched (education 3,590 against 3,593, marital groups within 2, insurance 11,257 against 11,273, tobacco 6,397 against 6,399, waist 8,204 exactly)"),
    list(choice = "missing covariates", decision = "complete case (n = 12,089)",
         reason = "the paper imputed income, education, marital status, insurance, tobacco, alcohol, and BMI by chained equations (m = 5), which this pipeline can't; keeping missing alcohol and income as their own levels (n = 14,314) gives the same estimate (variant)"),
    list(choice = "harmonized physical activity", decision = "leisure-time moderate-to-vigorous minutes a week, vigorous counted twice, 150 or more",
         reason = "2021-2023 asks only about leisure-time activity")
  )
)