# Dietary fiber intake (highest versus lowest tertile) and self-reported nonfatal myocardial
# infarction and/or stroke in adults, NHANES 2011-2018. Dong and Yang (2022), Front Nutr,
# doi:10.3389/fnut.2022.936926.
# Headline: OR 0.64 (95% CI 0.46-0.91), tertile 3 versus tertile 1 of fiber intake, Model 3 (Table 3).

# Yes, no, or unknown (anyone not asked, refused, or not sure), as the paper codes its medicine
# and smoking covariates.
row131_three <- function(yes, no) factor(ifelse(yes, "yes", ifelse(no, "no", "unknown")), levels = c("yes", "no", "unknown"))

row131_build <- function(cycle) {
  replication <- cycle == REPLICATION_CYCLE
  demo <- demographics(cycle)
  diet <- dietary_totals(cycle, c("FIBE", "KCAL"), days = 1)
  bp <- blood_pressure(cycle, readings = 1:3, zero_diastolic = "keep")
  bio <- biochemistry(cycle)
  smq <- component("SMQ", cycle, c("SMQ020", "SMQ040"))
  diq <- component("DIQ", cycle, c("DIQ010", "DIQ070"))
  bpq <- component("BPQ", cycle)
  d <- merge_all(demo, body_measures(cycle)[, c("SEQN", "bmi")], diet, bp, total_cholesterol(cycle), hdl_cholesterol(cycle),
                 hba1c(cycle), component("MCQ", cycle, c("MCQ160E", "MCQ160F")), smq, diq,
                 bpq[, c("SEQN", intersect(c("BPQ020", "BPQ040A", "BPQ050A", "BPQ150", "BPQ080", "BPQ090D", "BPQ100D", "BPQ101D"), names(bpq))), drop = FALSE],
                 component("RXQASA", cycle, "RXQ520"))
  at <- function(frame, column) frame[[column]][match(d$SEQN, frame$SEQN)]
  # Fiber from foods on the first recall day plus fiber from supplements: the same day's
  # supplement recall (DS1TFIBE) through 2018; 2021-2023 dropped that recall, so there the
  # supplement part is the 30-day supplement questionnaire's average daily fiber (DSQTFIBE).
  zero <- function(x) ifelse(is.na(x), 0, x)
  dsq <- zero(at(component("DSQTOT", cycle), "DSQTFIBE"))
  ds1 <- if (replication) dsq else zero(at(component("DS1TOT", cycle), "DS1TFIBE"))
  d$fiber_food <- d$FIBE
  d$fiber <- d$FIBE + ds1
  d$fiber_dsq <- d$FIBE + dsq
  d$kcal <- d$KCAL
  d$event <- ifelse(d$MCQ160E %in% 1 | d$MCQ160F %in% 1, 1, ifelse(d$MCQ160E %in% 2 & d$MCQ160F %in% 2, 0, NA))
  d$sex <- factor(d$sex, levels = 1:2, labels = c("men", "women"))
  d$race4 <- factor(ifelse(d$race %in% c(2, 5), 2, d$race), levels = 1:4, labels = c("Mexican American", "other", "NH White", "NH Black"))
  d$marital3 <- factor(d$marital, levels = 1:3)
  d$education3 <- factor(education3(d$education), levels = 1:3)
  d$pir2 <- factor(ifelse(d$pir < 1.2, "<1.2", ">=1.2"), levels = c("<1.2", ">=1.2"))
  d$bmi3 <- cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-30", ">=30"))
  # Mean of the three readings; Figure 1 drops anyone missing one (2008 people).
  complete_bp <- d$sbp_n %in% 3 & d$dbp_n %in% 3
  d$sbp <- ifelse(complete_bp, d$sbp, NA)
  d$dbp <- ifelse(complete_bp, d$dbp, NA)
  d$glucose <- at(bio, "glucose_serum")
  d$tg <- at(bio, "tg_serum")
  # Current smoker from SMQ040 (every day or some days, yes; not at all, no), unknown where it
  # was not asked (never smokers) or not answered, as the paper codes its other items. Table 2's
  # counts (546 yes, 583 no, 7743 unknown) and odds ratios are those of SMQ040 in 2011-2012 alone.
  d$current_smoker <- row131_three(d$SMQ040 %in% 1:2, d$SMQ040 %in% 3)
  d$current_smoker_table2 <- if (cycle == "2011-2012") d$current_smoker else row131_three(FALSE, FALSE)
  d$current_smoker_never_no <- row131_three(d$SMQ040 %in% 1:2, d$SMQ020 %in% 2 | d$SMQ040 %in% 3)
  d$diabetes3 <- factor(ifelse(d$DIQ010 %in% 1:3, d$DIQ010, NA), levels = c(2, 1, 3), labels = c("no", "yes", "borderline"))
  d$hypertension <- factor(ifelse(d$BPQ020 %in% 1:2, d$BPQ020, NA), levels = 2:1, labels = c("no", "yes"))
  d$hypercholesterolemia <- factor(ifelse(d$BPQ080 %in% 1:2, d$BPQ080, NA), levels = 2:1, labels = c("no", "yes"))
  d$hypoglycemic <- row131_three(d$DIQ070 %in% 1, d$DIQ070 %in% 2)
  d$aspirin <- row131_three(d$RXQ520 %in% 1, d$RXQ520 %in% 2)
  # Blood pressure and cholesterol medicines as the paper's cycles ask them: through 2018 only
  # those told to take a medicine (BPQ040A, BPQ090D) are asked whether they take it.
  d$antihypertensive <- if (has(d, "BPQ050A")) row131_three(d$BPQ050A %in% 1, d$BPQ050A %in% 2) else factor(NA, levels = c("yes", "no", "unknown"))
  d$lipid_lowering <- if (has(d, "BPQ100D")) row131_three(d$BPQ100D %in% 1, d$BPQ100D %in% 2) else factor(NA, levels = c("yes", "no", "unknown"))
  # The same medicines as 2021-2023 asks them: BPQ150 asks everyone told of high blood pressure,
  # BPQ101D asks everyone. Through 2018, those never told to take a medicine are not taking one.
  if (replication) {
    d$antihypertensive_h <- row131_three(d$BPQ150 %in% 1, d$BPQ020 %in% 1 & d$BPQ150 %in% 2)
    d$lipid_lowering_h <- row131_three(d$BPQ101D %in% 1, d$BPQ101D %in% 2)
  } else {
    d$antihypertensive_h <- row131_three(d$BPQ050A %in% 1, d$BPQ020 %in% 1 & (d$BPQ040A %in% 2 | d$BPQ050A %in% 2))
    unsure <- !(d$BPQ080 %in% 1:2) | d$BPQ090D %in% c(7, 9) | d$BPQ100D %in% c(7, 9)
    d$lipid_lowering_h <- row131_three(d$BPQ100D %in% 1, !(d$BPQ100D %in% 1) & !unsure)
  }
  d$vigorous <- d$sleep <- factor(NA, levels = c("no", "yes"))
  if (!replication) {
    paq <- component("PAQ", cycle, "PAQ605")
    slq <- component("SLQ", cycle, "SLQ050")
    d$vigorous <- factor(ifelse(at(paq, "PAQ605") %in% 1:2, at(paq, "PAQ605"), NA), levels = 2:1, labels = c("no", "yes"))
    d$sleep <- factor(ifelse(at(slq, "SLQ050") %in% 1:2, at(slq, "SLQ050"), NA), levels = 2:1, labels = c("no", "yes"))
  }
  # Figure 1: 18 and older (education, asked from 20, makes it 20 and older) with every key
  # variable known.
  known <- d$age >= 18 & !is.na(d$pir) & !is.na(d$education3) & !is.na(d$marital3) & !is.na(d$bmi) & !is.na(d$sbp) & !is.na(d$dbp) &
    !is.na(d$glucose) & !is.na(d$tc) & !is.na(d$tg) & !is.na(d$hdl) & !is.na(d$hba1c) &
    !is.na(d$diabetes3) & !is.na(d$hypertension) & !is.na(d$hypercholesterolemia) & !is.na(d$event) & !is.na(d$fiber)
  d$in_population <- known & !is.na(d$vigorous) & !is.na(d$sleep)
  d$in_population_harmonized <- known
  d
}

row131_formula <- event ~ fiber_t + age + sex + race4 + marital3 + education3 + pir2 + bmi3 + current_smoker + sbp + dbp + glucose + tc + tg + hdl +
  hba1c + kcal + vigorous + diabetes3 + hypertension + hypercholesterolemia + sleep + hypoglycemic + antihypertensive + lipid_lowering + aspirin

# Tertiles of fiber intake on the paper's analytic sample, unweighted: Table 2's counts (2929,
# 2959, 2984) fall the way left-closed groups split values tied at a cutpoint.
row131_tertile <- function(x, at) cut(x, c(-Inf, at, Inf), right = FALSE, labels = c("T1", "T2", "T3"))

association <- list(
  id = "row131", row = 131, doi = "10.3389/fnut.2022.936926",
  cycles = c("2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTDRD1",
  family = "logistic", term = "fiber_tT3",
  published = list(measure = "OR", estimate = 0.64, low = 0.46, high = 0.91, n = 8872, events = NULL,
                   contrast = "highest versus lowest tertile of fiber intake from foods and supplements, first 24-hour recall (tertile means 7.51 and 28.79 g/day)"),
  left_out = c(
    "vigorous work activity (PAQ605; 2021-2023 asks only leisure-time activity), as a covariate and as Figure 1's exclusion of those missing it",
    "sleep disorder (told a doctor of trouble sleeping, SLQ050, not asked in 2021-2023), as a covariate and as an exclusion",
    "fiber from supplements on the recall day (DS1TFIBE): 2021-2023 dropped the 24-hour supplement recall, so there the supplement part is the 30-day supplement questionnaire's average daily fiber (DSQTFIBE)",
    "antihypertensive and lipid-lowering medicine coded as the paper's cycles ask them: 2021-2023 asks everyone told of high blood pressure (BPQ150) and everyone (BPQ101D) whether they take one, so the harmonized model codes both cycles' answers that way"
  ),
  build = function(cycle) row131_build(cycle),
  constants = function(data) {
    keep <- data$in_population %in% TRUE & stats::complete.cases(data[, setdiff(all.vars(row131_formula), "fiber_t")]) & !is.na(data$w) & data$w > 0
    third <- function(x) unname(stats::quantile(x[keep], c(1 / 3, 2 / 3)))
    list(fiber = third(data$fiber), fiber_food = third(data$fiber_food), fiber_dsq = third(data$fiber_dsq))
  },
  derive = function(data, constants) {
    data$fiber_t <- row131_tertile(data$fiber, constants$fiber)
    data$fiber_food_t <- row131_tertile(data$fiber_food, constants$fiber_food)
    data$fiber_dsq_t <- row131_tertile(data$fiber_dsq, constants$fiber_dsq)
    data
  },
  # The published estimate's smoking covariate had current smoking only for 2011-2012 (Table 2's
  # counts and ORs fit SMQ040 in that cycle alone), so the paper version reproduces that; the
  # harmonized version, which also runs on 2021-2023, codes current smoking in every cycle.
  formula = update(row131_formula, . ~ . - current_smoker + current_smoker_table2),
  formula_harmonized = event ~ fiber_t + age + sex + race4 + marital3 + education3 + pir2 + bmi3 + current_smoker + sbp + dbp + glucose + tc + tg + hdl +
    hba1c + kcal + diabetes3 + hypertension + hypercholesterolemia + hypoglycemic + antihypertensive_h + lipid_lowering_h + aspirin,
  variants = list(
    list(label = "crude (published 0.64, 0.48-0.84)", formula = event ~ fiber_t),
    list(label = "age, sex, race (published 0.56, 0.42-0.76)", formula = event ~ fiber_t + age + sex + race4),
    list(label = "per g/day, Model 3 (published 0.98, 0.96-1.00)", formula = update(row131_formula, . ~ . - fiber_t + fiber), term = "fiber"),
    list(label = "examination weight", weight = "WTMEC2YR"),
    list(label = "unweighted", weighted = FALSE),
    list(label = "current smoking coded in every cycle", formula = row131_formula),
    list(label = "current smoker yes or no, never smokers no", formula = update(row131_formula, . ~ . - current_smoker + current_smoker_never_no)),
    list(label = "fiber from foods only", formula = update(row131_formula, . ~ . - fiber_t + fiber_food_t), term = "fiber_food_tT3"),
    list(label = "supplement fiber from the 30-day questionnaire", formula = update(row131_formula, . ~ . - fiber_t + fiber_dsq_t), term = "fiber_dsq_tT3")
  ),
  # 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 Methods adjust for current smoking (yes, no, or unknown), but Table 2's counts (546 yes, 583 no, 7,743 unknown) and odds ratios (no 1.32, unknown 0.86; ours 1.36, 0.89) are those of SMQ040 in 2011-2012 alone, with everyone in 2013-2018 unknown. Coded so, Model 3 gives 0.69 for the headline, and coded in every cycle 0.78."),
    list(kind = "sample", affects_headline = TRUE, followed = FALSE, evidence = "unresolved",
         detail = "Figure 1 drops 6,590 participants for missing vigorous activity, but the vigorous activity item (PAQ605) is missing for fewer than 10 of the 16,223 at that step, and no variable tried has missing values that fit the per-cycle losses Supplementary Figure S1 implies (71%, 53%, 58%, and 54% kept). The file keeps them (n = 15,065 against 8,872), and its weighted descriptives still match Table 1."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods exclude only those younger than 18, but the education question is asked from age 20, so Figure 1's 1,071 excluded for missing education include the 18- and 19-year-olds (its 20,220 is reproduced exactly so), and the sample is 20 and older."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Discussion says higher dietary fiber intake was independently associated with an increased prevalence of stroke and MI, but every odds ratio for higher intake in Table 3 is below 1.")
  ),
  choices = list(
    list(choice = "headline", decision = "tertile 3 versus tertile 1, Model 3 (Table 3)", reason = "the study's rule: the main table's most-adjusted highest-versus-reference estimate, since the per-gram estimate's CI (0.96-1.00) does not exclude 1 and no p-value is given"),
    list(choice = "weights", decision = "dietary day-one weight (WTDRD1) over the four cycles", reason = "the paper names only 'an appropriate NHANES sample weight'; NCHS directs the weight of the smallest subsample, the day-one recall; weighted univariate odds ratios match Table 2 (race, BMI, HbA1c) where unweighted ones do not"),
    list(choice = "population", decision = "18 and older with every key variable known, in Figure 1's steps; education, asked from 20, makes it 20 and older", reason = "reproduces Figure 1's 20220 and 16223 exactly"),
    list(choice = "vigorous activity exclusion", decision = "drop only those missing PAQ605 (fewer than 10)", reason = "Figure 1 drops 6590 for missing vigorous activity, which PAQ605 cannot produce; no NHANES variable's missingness fits the per-cycle losses Supplementary Figure S1 implies (keeping 71%, 53%, 58%, 54%), so the sample is 15065 against 8872; weighted descriptives still match Table 1 (tertile means, prevalence by tertile 6.34/5.50/4.39 against 6.50/5.45/4.25, sex, age, vigorous 23%, sleep 30%)"),
    list(choice = "blood pressure", decision = "mean of the first three readings, anyone missing one excluded, diastolic zeros kept as recorded", reason = "the paper averages three readings; requiring all three reproduces Figure 1's 2008"),
    list(choice = "glucose, triglycerides, cholesterol", decision = "serum glucose and triglycerides from the biochemistry panel (LBXSGL, LBXSTR), total cholesterol LBXTC, HDL LBDHDD, HbA1c LBXGH", reason = "unstated; the fasting values exist for half the sample, while these reproduce Figure 1's 996 missing and Table 1's means"),
    list(choice = "fiber", decision = "day-one food fiber (DR1TFIBE, reliable recalls) plus that day's supplement fiber (DS1TFIBE, none counted as 0)", reason = "the paper sums food and supplements from the first 24-hour recall without naming the supplement file; the day-one supplement recall is that recall's"),
    list(choice = "tertiles", decision = "cutpoints at the unweighted thirds of the analytic sample (11.3 and 19.3 g/day), groups closed on the left", reason = "cutpoints are unpublished; Table 1's unequal weighted shares (30.88/34.82/34.29) show unweighted tertiles, Table 2's sizes (2929/2959/2984) fall as left-closed groups split values tied at a cutpoint, and the weighted tertile means come out 7.51/15.07/29.23 against 7.51/14.90/28.79"),
    list(choice = "outcome", decision = "yes to heart attack or stroke (MCQ160E, MCQ160F); no to both; missing otherwise", reason = "the paper's composite"),
    list(choice = "current smoker", decision = "yes if SMQ040 every day or some days, no if not at all, unknown where not asked (never smokers) or not answered; in the paper version for 2011-2012 only (unknown in the other cycles), as the paper's computation merged it, and in the harmonized version for every cycle", reason = "Table 2's counts (546/583/7743) and odds ratios (no 1.32, unknown 0.86 against our 1.36, 0.89) are SMQ040 coded this way in 2011-2012 alone, so the paper merged smoking for one cycle only; a data-handling failure of the paper's own files, followed in the paper version and not in the harmonized one, as the plan sets; the paper version with smoking coded in every cycle is a variant"),
    list(choice = "medicines and aspirin", decision = "DIQ070, BPQ050A, BPQ100D, and RXQ520 (taking low-dose aspirin on one's own) as yes, no, or unknown where not asked", reason = "Table 2's counts fit these items, aspirin's 190/3437 fitting RXQ520 and not RXQ515"),
    list(choice = "conditions", decision = "diabetes DIQ010 in three levels; hypertension BPQ020; hypercholesterolemia BPQ080; sleep disorder SLQ050 (told a doctor of trouble sleeping); vigorous activity PAQ605 (work)", reason = "the paper's categories; SLQ060 was asked only through 2014, and Table 1's 30% fits SLQ050; the paper's definition of vigorous activity is PAQ605's wording"),
    list(choice = "other covariates", decision = "race from RIDRETH1 with other Hispanic among other races; marital status in three groups; education in three; PIR under 1.2 or not; BMI under 25, 25-30, 30 or more; age, blood pressure, labs, and day-one energy continuous", reason = "the paper's categories"),
    list(choice = "2021-2023 blood pressure", decision = "the oscillometric readings (BPXO)", reason = "2021-2023 measures blood pressure only by the oscillometric device; the library maps it")
  )
)