# Dietary energy intake per kg of body weight and asthma in adults, NHANES 2009-2018.
# Cao et al. (2024), BMC Nutrition, doi:10.1186/s40795-024-00938-7.
# Headline: OR 0.61 (95% CI 0.53-0.69), energy intake quartile 4 (> 33.29 kcal/kg/day) vs quartile 1
# (< 17.73 kcal/kg/day), Model 3 (Table 2).

# Leisure-time activity in the paper's three groups (vigorous, moderate, sedentary): 2009-2018 ask
# whether one does vigorous (PAQ650) or moderate (PAQ665) recreational activity for at least 10
# minutes; 2021-2023 asks how often (PAD810Q vigorous, PAD790Q moderate, 0 = never).
row249_activity <- function(cycle) {
  d <- component("PAQ", cycle)
  if (has(d, "PAQ650")) {
    vigorous <- yes(d$PAQ650)
    moderate <- yes(d$PAQ665)
  } else {
    often <- function(q) ifelse(q %in% c(7777, 9999) | is.na(q), NA, as.integer(q > 0))
    vigorous <- often(d$PAD810Q)
    moderate <- often(d$PAD790Q)
  }
  level <- ifelse(vigorous %in% 1, "vigorous", ifelse(moderate %in% 1, "moderate",
                  ifelse(vigorous %in% 0 & moderate %in% 0, "sedentary", NA)))
  data.frame(SEQN = d$SEQN, activity = level)
}

# The paper's quartile cutpoints of energy intake (kcal/kg/day; Methods and Table 2).
ROW249_CUTS <- c(17.73, 24.52, 33.29)
row249_quartile <- function(x) cut(x, c(-Inf, ROW249_CUTS, Inf), right = FALSE, labels = c("Q1", "Q2", "Q3", "Q4"))

# Table 2's models, as its note lists them (the Methods also put activity and smoking in Models 2
# and 3).
ROW249_MODEL_1 <- asthma ~ energy_q + sex + age + race + family_asthma + activity + smoking + pir3
ROW249_MODEL_2 <- asthma ~ energy_q + sex + age + race + family_asthma + pir3 + wbc + eosinophil_pct + hemoglobin + vitamin_d
ROW249_MODEL_3 <- update(ROW249_MODEL_2, . ~ . + n3 + n6 + fiber + zinc)

association <- list(
  id = "row249", row = 249, doi = "10.1186/s40795-024-00938-7",
  cycles = c("2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  # The exposure is the first day's recall, so a weighted fit takes the day-one dietary weight, in
  # 2021-2023 too (blood_file changes only an examination weight).
  weight = "WTDRD1", blood_file = "CBC",
  # The paper mentions no weights, and only unweighted fits reproduce Table 2.
  weighted = FALSE,
  family = "logistic", term = "energy_qQ4",
  published = list(measure = "OR", estimate = 0.61, low = 0.53, high = 0.69, n = 21354, events = 3154,
                   contrast = "energy intake quartile 4 (> 33.29 kcal/kg/day) vs quartile 1 (< 17.73 kcal/kg/day)"),
  left_out = c("familial asthma (MCQ300B, a close relative ever told they had asthma): 2021-2023 does not ask it; it leaves the model and the complete-case step"),
  build = function(cycle) {
    demo <- demographics(cycle)
    nutrients <- c("KCAL", "P183", "P184", "P205", "P225", "P226", "P182", "P204", "FIBE", "ZINC")
    day1 <- dietary_totals(cycle, nutrients, days = 1)
    two_day <- dietary_totals(cycle, "KCAL", days = 2)
    mcq <- component("MCQ", cycle)
    d <- merge_all(demo, day1, body_measures(cycle)[, c("SEQN", "weight", "bmi")],
                   component("CBC", cycle, c("LBXWBCSI", "LBXEOPCT", "LBXHGB")),
                   component("VID", cycle, "LBXVIDMS"), smoking(cycle), row249_activity(cycle))
    d$KCAL2 <- two_day$KCAL[match(d$SEQN, two_day$SEQN)]
    d$asthma <- yes(mcq$MCQ010)[match(d$SEQN, mcq$SEQN)]
    # 2021-2023 does not ask about family history of asthma.
    d$family_asthma <- if (has(mcq, "MCQ300B")) factor(yes(mcq$MCQ300B)[match(d$SEQN, mcq$SEQN)]) else NA
    d$energy_kg <- d$KCAL / d$weight
    d$energy_q <- row249_quartile(d$energy_kg)
    d$energy_q_two_day <- row249_quartile(d$KCAL2 / d$weight)
    d$sex <- factor(d$sex)
    d$race <- factor(d$race)
    # Income in the paper's groups: low (PIR 1.3 or less), medium (over 1.3 to 3.5), high (over 3.5).
    d$pir3 <- cut(d$pir, c(-Inf, 1.3, 3.5, Inf), labels = c("low", "medium", "high"))
    d$wbc <- d$LBXWBCSI
    d$eosinophil_pct <- d$LBXEOPCT
    d$hemoglobin <- d$LBXHGB
    d$vitamin_d <- d$LBXVIDMS
    # Fatty acids in g/day: n-3 is ALA, SDA, EPA, DPA, and DHA; n-6 is linoleic and arachidonic acid.
    d$n3 <- d$P183 + d$P184 + d$P205 + d$P225 + d$P226
    d$n6 <- d$P182 + d$P204
    d$n3_kg <- 1000 * d$n3 / d$weight
    d$n6_kg <- 1000 * d$n6 / d$weight
    d$fiber <- d$FIBE
    d$zinc <- d$ZINC
    d$smoking <- factor(d$smoking)
    d$activity <- factor(d$activity, levels = c("sedentary", "moderate", "vigorous"))
    # Adults with an asthma answer and the first day's energy intake, and, as the paper's
    # complete-case step excludes anyone missing a covariate it describes, BMI, smoking status, and
    # physical activity known (with the Model 3 covariates, 21,354 people, as published).
    d$in_population_model_covariates <- d$age >= 20 & !is.na(d$asthma) & !is.na(d$energy_kg) & !is.na(d$bmi)
    d$in_population <- d$in_population_model_covariates & !is.na(d$smoking) & !is.na(d$activity)
    d
  },
  formula = ROW249_MODEL_3,
  formula_harmonized = update(ROW249_MODEL_3, . ~ . - family_asthma),
  variants = list(
    list(label = "weighted (day-one dietary weight)", weighted = TRUE),
    list(label = "unadjusted (published 0.74, 0.67-0.83)", formula = asthma ~ energy_q),
    list(label = "Model 1 (published 0.77, 0.69-0.86)", formula = ROW249_MODEL_1),
    list(label = "Model 2 (published 0.76, 0.68-0.85)", formula = ROW249_MODEL_2),
    list(label = "Model 3 with activity and smoking, as the Methods list it", formula = update(ROW249_MODEL_3, . ~ . + activity + smoking)),
    list(label = "energy as the two-day mean the Methods describe", formula = update(ROW249_MODEL_3, . ~ . - energy_q + energy_q_two_day),
         term = "energy_q_two_dayQ4", sample = "own"),
    list(label = "n-3 and n-6 PUFAs in mg/kg/day (Table 1's unit)", formula = update(ROW249_MODEL_3, . ~ . - n3 - n6 + n3_kg + n6_kg)),
    list(label = "smoking and activity not required (n 21,366)",
         derive = function(data, constants) { data$in_population <- data$in_population_model_covariates; data }),
    list(label = "unadjusted, weighted", formula = asthma ~ energy_q, weighted = TRUE)
  ),
  # 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,
         detail = "The Methods average energy intake over the two 24-hour recalls, but the authors' dataset (Supplementary Material 2) holds only the first day's energy, which matches this build person by person, and the first day's energy gives the paper's n of 21,354 and Table 2 (Model 3 0.612, 0.533-0.702); the two-day mean leaves 18,816 people and gives 0.653 (0.567-0.753)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "Table 1 and Fig. 2's list of adjustments give n-3 and n-6 PUFAs in mg/kg/day, but Table 1's medians (1.5 and 14.4) and the n3 and n6 columns of the authors' dataset are in g/day, and Model 3 is reproduced with them in g/day (0.612; in mg/kg/day, 0.632)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Methods, Table 2's note, and Fig. 3 give Q3 as 21.52 to 33.29 kcal/kg/day, but Q2 ends at 24.52, which is also the second quartile boundary of the authors' dataset.")
  ),
  choices = list(
    list(choice = "weights", decision = "unweighted",
         reason = "the paper never mentions weights, strata, or PSUs; unweighted fits reproduce Table 2's unadjusted, Model 1, Model 2, and Model 3 quartile 4 estimates within 0.01, weighted ones do not (unadjusted 0.82 against 0.74)"),
    list(choice = "energy intake", decision = "the first day's recall (DR1TKCAL, reliable recalls) divided by measured weight (BMXWT)",
         reason = "the Methods describe the mean of two recalls, but the authors' dataset holds only the first day's energy, and it matches this build person by person; the two-day mean is a variant"),
    list(choice = "quartile cutpoints", decision = "the published 17.73, 24.52, and 33.29 kcal/kg/day, each cutpoint opening the higher quartile",
         reason = "Methods and Table 2; the Methods, Table 2's note, and Fig. 3 misprint Q3's lower bound as 21.52"),
    list(choice = "study population", decision = "aged 20 or older, asthma answered yes or no, first-day energy and weight known, and BMI, smoking status, and physical activity known besides the Model 3 covariates",
         reason = "the paper drops anyone missing a covariate it describes; requiring these gives the published 21,354 people and 3,154 with asthma exactly"),
    list(choice = "familial asthma", decision = "MCQ300B yes or no, refused or don't know as missing",
         reason = "unstated; complete-case analysis"),
    list(choice = "income to poverty ratio", decision = "the paper's three groups: 1.3 or less, over 1.3 to 3.5, over 3.5", reason = "Methods"),
    list(choice = "race and ethnicity", decision = "RIDRETH1's five groups", reason = "Table 1 shows five; the authors' dataset codes RIDRETH1"),
    list(choice = "n-3 and n-6 PUFAs", decision = "g/day; n-3 is 18:3, 18:4, 20:5, 22:5, and 22:6, n-6 is 18:2 and 20:4",
         reason = "these sums equal the authors' n3 and n6 columns for every participant; Table 1 labels them mg/kg/day, but its medians are g/day (mg/kg/day is a variant)"),
    list(choice = "Model 3 covariates", decision = "Table 2's note: no physical activity or smoking",
         reason = "the Methods build Models 2 and 3 on Model 1, which has both; the note's list gives Model 3's printed quartile 4 estimate (0.612 against 0.603 with them), though with them Model 2's (0.762 against 0.766 for a printed 0.76) fits better; both reproduce the headline, and the Methods' version is a variant"),
    list(choice = "physical activity", decision = "vigorous or moderate recreational activity (PAQ650, PAQ665) for 2009-2018; how often one does vigorous or moderate leisure-time activity (PAD810Q, PAD790Q) for 2021-2023",
         reason = "the paper's three groups (sedentary, moderate, vigorous) without saying which items; it enters only Model 1 and the complete-case step"),
    list(choice = "pregnancy", decision = "pregnant women kept", reason = "the paper names no such exclusion, and its sample size is reproduced without one"),
    list(choice = "2021-2023 recall", decision = "the first day's recall as released",
         reason = "2021-2023 took both recalls by telephone, where 2009-2018 took the first in person at the examination; nothing else in the exposure changes")
  )
)