# Red cell distribution width (RDW) and coronary heart disease in adults with rheumatoid arthritis,
# NHANES 2011-March 2020. Zhang et al. (2024), Medicine, doi:10.1097/MD.0000000000037315.
# Headline: OR 1.145 (95% CI 1.036-1.266) per 1% of RDW, Model 1 (Table 3, unadjusted, as is Table
# 2's univariate estimate; the abstract calls it adjusted). Model 3 (1.187, 1.065-1.322) is a variant.

# Model 2: age, gender, and race. Model 3: "all covariants in Table 1".
row072_model2 <- chd ~ rdw + age + sex + race4
row072_model3 <- chd ~ rdw + age + pir + bmi + sex + race4 + education3 + hypertension + smoking + drinking + diabetes +
  hypercholesterolemia

association <- list(
  id = "row072", row = 72, doi = "10.1097/MD.0000000000037315",
  cycles = c("2011-2012", "2013-2014", "2015-2016", "2017-2020"),
  weight = "WTMEC2YR", blood_file = "CBC",
  family = "logistic", term = "rdw",
  published = list(measure = "OR", estimate = 1.145, low = 1.036, high = 1.266, n = 1236, events = 102,
                   contrast = "per 1% of red cell distribution width"),
  left_out = character(),
  build = function(cycle) {
    paq <- component("PAQ", cycle)
    # Moderate recreational activity (a variant's covariate): PAQ665 through March 2020; 2021-2023
    # asks how often instead (PAD790Q, 0 for none).
    active <- if (has(paq, "PAQ665")) paq$PAQ665 %in% 1 else paq$PAD790Q %in% 1:180
    d <- merge_all(demographics(cycle), blood_count(cycle)[, c("SEQN", "rdw")], arthritis(cycle),
                   component("MCQ", cycle, "MCQ160C"), body_measures(cycle)[, c("SEQN", "bmi")],
                   component("BPQ", cycle, c("BPQ020", "BPQ080")), component("SMQ", cycle, "SMQ020"),
                   component("DIQ", cycle, "DIQ010"), component("ALQ", cycle, "ALQ151"),
                   data.frame(SEQN = paq$SEQN, active = active))
    d$chd <- ifelse(d$MCQ160C %in% 1, 1, ifelse(d$MCQ160C %in% 2, 0, NA))
    d$sex <- factor(d$sex, levels = 1:2, labels = c("male", "female"))
    # Mexican American (reference), non-Hispanic White, non-Hispanic Black, and other (other
    # Hispanic with other or multiracial).
    d$race4 <- factor(ifelse(d$race %in% 1, "mexican", ifelse(d$race %in% 3, "white", ifelse(d$race %in% 4, "black", "other"))),
                      levels = c("mexican", "white", "black", "other"))
    # Table 1's groups, with every answer but the named ones in the reference or "no" group.
    d$education3 <- factor(ifelse(d$education %in% 2:3, "high school", ifelse(d$education %in% 4:5, "university", "less than high school")),
                           levels = c("less than high school", "high school", "university"))
    d$hypertension <- factor(ifelse(d$BPQ020 %in% 1, "yes", "no"), levels = c("yes", "no"))
    d$hypercholesterolemia <- factor(ifelse(d$BPQ080 %in% 1, "yes", "no"), levels = c("yes", "no"))
    d$smoking <- factor(ifelse(d$SMQ020 %in% 1, "yes", "no"), levels = c("yes", "no"))
    d$diabetes <- factor(ifelse(d$DIQ010 %in% c(1, 3), "yes", "no"), levels = c("yes", "no"))
    # "Regular alcohol consumption": ever had 4 (women) or 5 (men) or more drinks almost every day
    # (ALQ151), with those not asked or not answering as "other".
    d$drinking <- factor(ifelse(d$ALQ151 %in% 1, "yes", ifelse(d$ALQ151 %in% 2, "no", "other")), levels = c("yes", "no", "other"))
    d$activity <- factor(ifelse(d$active %in% TRUE, "yes", "no"), levels = c("yes", "no"))
    d$in_population <- !is.na(d$rdw) & !is.na(d$chd) & d$rheumatoid %in% 1
    d
  },
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "each cycle's weight over 4", weight_rule = "equal"),
    list(label = "interview weight (WTINT2YR, WTINTPRP)", weight = "WTINT2YR"),
    list(label = "interview weight, each cycle's over 4", weight = "WTINT2YR", weight_rule = "equal"),
    list(label = "Model 2, weighted (published 1.168, 1.060-1.287)", formula = row072_model2),
    list(label = "Model 2, unweighted", formula = row072_model2, weighted = FALSE),
    list(label = "Model 3, weighted (published 1.187, 1.065-1.322)", formula = row072_model3),
    list(label = "Model 3, unweighted", formula = row072_model3, weighted = FALSE),
    list(label = "Model 3 with moderate recreational activity, weighted", formula = update(row072_model3, . ~ . + activity))
  ),
  # 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 = "reporting", affects_headline = TRUE, followed = TRUE, evidence = "paper",
         detail = "The abstract presents OR 1.145 (1.036-1.266) as holding after adjusting for age, gender, race, education level, smoking, and drinking, but it is Table 3's Model 1, which the table note says adjusts for no covariates, and Table 2's univariate estimate; RDW alone reproduces it (1.152, 1.045-1.269)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract, the text, and Tables 1, 2, and 4 name two race groups 'Hispanic White' and 'Hispanic Black', but Table 1's counts for them (434 and 398) are those of non-Hispanic White and non-Hispanic Black participants."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The text and Figure 1 give 1,236 participants, but Table 1's tertile columns sum to 1,235 (402, 408, 425), the number the paper's exclusion steps give when reproduced (with its 36,770 with RDW and 23,659 with heart disease answered exactly)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2's univariate odds ratios for the covariates are not reproduced under any weighting by data that reproduce Table 1's counts exactly (women: 0.235 published; 0.66 with the examination weight, 0.37 unweighted)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Results give the participants' mean age as 61.107 years, citing Table 1, but Table 1's tertile means (57.858, 61.027, and 62.915 years for 402, 408, and 425 participants) average to 60.645."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract names the age subgroup 55-66 years, which the Results and Table 4 give as 56 to 66, and the Results give the male subgroup's P as .084, which its interval (1.034-1.442) rules out and Table 4 gives as .0184.")
  ),
  choices = list(
    list(choice = "headline model", decision = "Model 1: RDW alone",
         reason = "the abstract's estimate is Table 3's Model 1 and Table 2's univariate estimate, which the table note says adjusts for nothing, although the abstract says it holds after adjustment"),
    list(choice = "cycles", decision = "2011-2012, 2013-2014, 2015-2016, and the 2017-March 2020 files, without the 2017-2018 files",
         reason = "the paper's 45,462 participants are exactly these four cycles' DEMO files; with the 2017-2018 files instead of the 2017-March 2020 ones the total would be 39,156"),
    list(choice = "population", decision = "RDW measured (LBXRDW), coronary heart disease answered yes or no (MCQ160C 1 or 2), and told of arthritis (MCQ160A = 1) of the rheumatoid type (MCQ195 = 2); no age step, since the arthritis and heart disease items are asked from age 20",
         reason = "the paper's steps; they reproduce its flow exactly (36,770 with RDW, 23,659 with heart disease answered) and give 1,235 with rheumatoid arthritis and 102 cases, Table 1's totals (its tertiles sum to 1,235, its cases to 102; the text says 1,236). Counting refused or don't know as answered would give 23,739 and 1,243"),
    list(choice = "weights", decision = "examination weight (WTMEC2YR; WTMECPRP for 2017-March 2020), each cycle's times its share of the 9.2 pooled years as NCHS directs, with strata and PSUs; WTPH2YR from CBC_L in 2021-2023",
         reason = "the paper says only 'a weighted approach'; RDW is measured at the examination, so NCHS's guidelines call for the examination weight. Weighted fits come closer than unweighted ones to Model 1 (1.152 against 1.145; unweighted 1.110) and Model 2 (1.158 against 1.168; unweighted 1.134). No weighting reproduces every model: the interview weight with each cycle's weight over 4 (weights stacked unscaled) matches Model 1 most closely (1.148, 1.037-1.271, P .0098 as published) but gives 1.149 (P .04) for Model 2, against the published 1.168 (P .0026). The interview weight and dividing each cycle's weight by 4 are variants"),
    list(choice = "RDW scale", decision = "continuous, per percentage point (LBXRDW, %)",
         reason = "the paper reports the OR per unit of RDW and describes RDW as a percentage"),
    list(choice = "coronary heart disease", decision = "ever told of coronary heart disease (MCQ160C = 1) against not (2)",
         reason = "the paper's 'CHD values'; refused or don't know answers are missing, which reproduces the flow's 13,111 excluded"),
    list(choice = "Model 2 and Model 3 covariates (variants only)",
         decision = "age, income-to-poverty ratio, and BMI continuous; gender; race in four groups (other Hispanic with other); education less than high school (DMDEDUC2 1, or unanswered), high school (2-3), university or above (4-5); hypertension (BPQ020), hypercholesterolemia (BPQ080), smoking (SMQ020), each yes against anything else; diabetes (DIQ010) yes or borderline against anything else; drinking as ever 4 or 5 or more drinks almost every day (ALQ151) yes, no, or other; complete cases for income and BMI",
         reason = "unstated definitions; each grouping reproduces Table 1's counts exactly (race 161/434/398/242, education 169/502/564, hypertension 773/462, hypercholesterolemia 632/603, smoking 641/594, diabetes 389/846, drinking 242/750/243). Income is missing for 133 and BMI for 24, so Model 3 has 1,082 participants and 92 cases. Model 3 is not reproduced closely (1.127 weighted, 1.170 unweighted, against 1.187), and neither are Table 2's univariate estimates for the covariates under any weighting (women 0.66 weighted and 0.37 unweighted against 0.235; poverty 0.97 and 1.07 against 1.205), so those rows come from some other computation"),
    list(choice = "physical activity in Model 3", decision = "left out",
         reason = "the Methods list moderate recreational activity, but Table 1, which defines Model 3, does not; adding it (PAQ665) is a variant")
  ),
  formula = chd ~ rdw
)