# Serum 25-hydroxyvitamin D and self-reported osteoarthritis in adults 40 and older, NHANES
# 2001-2018. Yu, Lin, Dai, Xu, and Liu (2023), Front Nutr, doi:10.3389/fnut.2023.1016809.
# Headline: OR 1.91 (95% CI 1.59-2.30) for 25(OH)D of 100 nmol/L or more against under 50 nmol/L,
# Model 4 (Table 3).

# Total 25(OH)D in nmol/L: NCHS's LC-MS/MS-equivalent values for 2001-2006 (LBDVIDMS, converted
# from the RIA measurements and released in October 2015) and the LC-MS/MS values from 2007 on
# (LBXVIDMS).
row014_vitamin_d <- function(cycle) {
  d <- component("VID", cycle)
  data.frame(SEQN = d$SEQN, vitd = if (has(d, "LBXVIDMS")) d$LBXVIDMS else d$LBDVIDMS)
}

# Alcohol as "in any one year, had at least 12 drinks" (ALD100 in 2001-2002, ALQ101 in 2003-2016).
# 2017-2018 dropped the question; there the past year's drinks are counted from how often
# (ALQ121, the library's days a year) and how many on a drinking day (ALQ130), with never
# drinkers (ALQ111) and no drinking in the past year as fewer than 12. That reading gives Table 1's
# 72.1% with 12 or more (72.0%; ever drinking, ALQ111, gives 77.1%). Not built in 2021-2023, whose
# analysis leaves it out.
row014_alcohol <- function(cycle) {
  d <- component("ALQ", cycle)
  if (has(d, "ALD100") || has(d, "ALQ101")) {
    return(data.frame(SEQN = d$SEQN, alcohol12 = yes(if (has(d, "ALD100")) d$ALD100 else d$ALQ101), alcohol_ever = NA_real_))
  }
  a <- alcohol(cycle)
  drinks <- ifelse(a$drinking_days %in% 0, 0, a$drinking_days * a$drinks_per_day)
  drinks[a$alcohol3 %in% 1] <- 0
  data.frame(SEQN = a$SEQN, alcohol12 = ifelse(drinks >= 12, 1, ifelse(drinks < 12, 0, NA)), alcohol_ever = yes(d$ALQ111))
}

# Any moderate or vigorous recreational activity: over the past 30 days in 2001-2006 (PAD200,
# PAD320: 1 yes, 2 no, 3 unable to do), in a typical week in 2007-2018 (PAQ650, PAQ665), and any
# reported frequency of moderate or vigorous leisure-time activity in 2021-2023 (PAD790Q, PAD810Q,
# 0 for never).
row014_activity <- function(cycle) {
  d <- component("PAQ", cycle)
  if (has(d, "PAD790Q")) {
    count <- function(q) ifelse(q %in% c(7777, 9999), NA, q)
    moderate <- count(d$PAD790Q); vigorous <- count(d$PAD810Q)
    return(data.frame(SEQN = d$SEQN, active = ifelse(moderate > 0 | vigorous > 0, 1, ifelse(moderate %in% 0 & vigorous %in% 0, 0, NA))))
  }
  vigorous <- if (has(d, "PAD200")) d$PAD200 else d$PAQ650
  moderate <- if (has(d, "PAD320")) d$PAD320 else d$PAQ665
  data.frame(SEQN = d$SEQN, active = ifelse(vigorous %in% 1 | moderate %in% 1, 1, ifelse(vigorous %in% 2:3 & moderate %in% 2:3, 0, NA)))
}

# Use of a dietary supplement containing vitamin D in the past 30 days, as NCHS recommends finding
# it: each product a participant reported (DSQ2 in 2001-2006, DSQIDS from 2007) is looked up in
# the NHANES Dietary Supplement Database (DSPI maps the old product IDs of 2001-2016 to the current
# ones; DSII lists each product's ingredients, vitamin D being ingredient 364). Users have at least
# one such product; non-users said they took no supplement (DSD010, in DSQ1 or DSQTOT) or reported
# only products whose ingredients are known and lack vitamin D; anyone else (a product with no
# information, or no answer) is missing. 2021-2023 names products by their current ID (DSDPID in
# DSQIDS_L) and asked these questions by telephone after the first dietary recall.
row014_supplement <- function(cycle) {
  early <- cycle %in% c("1999-2000", "2001-2002", "2003-2004", "2005-2006")
  answers <- component(if (early) "DSQ1" else "DSQTOT", cycle)
  products <- component(if (early) "DSQ2" else "DSQIDS", cycle)
  catalog <- component("DSPI", "1999-2000")
  ingredients <- component("DSII", "1999-2000")
  if (!all(ingredients$DSDINGR[ingredients$DSDIID %in% 364] == "VITAMIN D")) stop("DSII ingredient 364 is not vitamin D")
  id <- if (has(products, "DSDPID")) products$DSDPID else catalog$DSDPID[match(as.character(products$DSDSUPID), catalog$DSDSUPID)]
  users <- unique(products$SEQN[id %in% ingredients$DSDPID[ingredients$DSDIID %in% 364]])
  unknown <- unique(products$SEQN[!id %in% ingredients$DSDPID])
  took <- answers$DSD010
  data.frame(SEQN = answers$SEQN,
             vitd_supplement = ifelse(answers$SEQN %in% users, 1, ifelse(took %in% 2 | (took %in% 1 & !answers$SEQN %in% unknown), 0, NA)),
             any_supplement = ifelse(took %in% 1, 1, ifelse(took %in% 2, 0, NA)))
}

association <- list(
  id = "row014", row = 14, doi = "10.3389/fnut.2023.1016809",
  cycles = c("2001-2002", "2003-2004", "2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTMEC2YR", blood_file = "VID",
  family = "logistic", term = "vitd_group100+",
  published = list(measure = "OR", estimate = 1.91, low = 1.59, high = 2.30, n = 21334, events = NA,
                   contrast = "serum 25(OH)D of 100 nmol/L or more vs under 50 nmol/L"),
  left_out = c(
    "alcohol as at least 12 drinks in any one year (ALD100/ALQ101, asked through 2015-2016; 2017-2018 and 2021-2023 ask only about ever drinking and the past year)"
  ),
  build = function(cycle) {
    demo <- demographics(cycle)
    d <- merge_all(demo, component("DEMO", cycle, "RIDEXMON"), row014_vitamin_d(cycle), arthritis(cycle),
                   body_measures(cycle)[, c("SEQN", "bmi")], smoking(cycle), row014_activity(cycle),
                   component("HUQ", cycle, "HUQ010"), row014_supplement(cycle))
    if (cycle == REPLICATION_CYCLE) {
      d$alcohol12 <- NA_real_; d$alcohol_ever <- NA_real_
    } else {
      d <- merge_all(d, row014_alcohol(cycle))
    }
    # Clinical cut-offs, the lowest the reference. The paper prints "<50" and ">=100"; 50 and 75
    # are taken to open their groups, as 100 does.
    groups <- c("<50", "50-75", "75-100", "100+")
    d$vitd_group <- cut(d$vitd, c(-Inf, 50, 75, 100, Inf), right = FALSE, labels = groups)
    d$vitd_group_upper <- cut(d$vitd, c(-Inf, 50, 75, 100, Inf), right = TRUE, labels = c("<=50", "50-75", "75-100", "over100"))
    d$vitd10 <- d$vitd / 10
    d$sex <- factor(d$sex)
    d$race4 <- factor(ifelse(d$race == 1, "mexican american", ifelse(d$race == 3, "white", ifelse(d$race == 4, "black", ifelse(d$race %in% c(2, 5), "other", NA)))),
                      levels = c("white", "black", "mexican american", "other"))
    d$education3 <- factor(education3(d$education))
    d$pir3 <- cut(d$pir, c(-Inf, 1.3, 3.5, Inf), right = FALSE, labels = c("<1.3", "1.3-3.5", "3.5+"))
    d$bmi3 <- cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-30", "30+"))
    d$season <- factor(ifelse(d$RIDEXMON %in% 1:2, d$RIDEXMON, NA), levels = 1:2, labels = c("november-april", "may-october"))
    d$smoking <- factor(d$smoking, levels = 1:3)
    d$health3 <- factor(ifelse(d$HUQ010 %in% 1:2, "excellent or very good", ifelse(d$HUQ010 %in% 3, "good", ifelse(d$HUQ010 %in% 4:5, "fair or poor", NA))),
                        levels = c("excellent or very good", "good", "fair or poor"))
    d$in_population <- d$age >= 40 & !is.na(d$vitd) & !is.na(d$osteoarthritis)
    d
  },
  formula = osteoarthritis ~ vitd_group + age + sex + race4 + education3 + pir3 + bmi + season + alcohol12 + smoking + active + vitd_supplement + health3,
  formula_harmonized = osteoarthritis ~ vitd_group + age + sex + race4 + education3 + pir3 + bmi + season + smoking + active + vitd_supplement + health3,
  variants = list(
    list(label = "BMI in the text's groups (<25, 25-30, 30+)",
         formula = osteoarthritis ~ vitd_group + age + sex + race4 + education3 + pir3 + bmi3 + season + alcohol12 + smoking + active + vitd_supplement + health3),
    list(label = "Model 1, crude (published 2.56, 2.16-3.04)", formula = osteoarthritis ~ vitd_group),
    list(label = "Model 2 (published 1.93, 1.61-2.32)", formula = osteoarthritis ~ vitd_group + age + sex + race4 + education3 + pir3 + bmi + season),
    list(label = "Model 3 (published 1.86, 1.55-2.25)",
         formula = osteoarthritis ~ vitd_group + age + sex + race4 + education3 + pir3 + bmi + season + alcohol12 + smoking + active + vitd_supplement),
    list(label = "per 10 nmol/L (published 1.07, 1.05-1.09)", term = "vitd10",
         formula = osteoarthritis ~ vitd10 + age + sex + race4 + education3 + pir3 + bmi + season + alcohol12 + smoking + active + vitd_supplement + health3),
    list(label = "50, 75, 100 closing their groups (>100 vs <=50)", term = "vitd_group_upperover100",
         formula = osteoarthritis ~ vitd_group_upper + age + sex + race4 + education3 + pir3 + bmi + season + alcohol12 + smoking + active + vitd_supplement + health3),
    list(label = "2017-2018 alcohol as ever drinking (ALQ111)", sample = "own",
         derive = function(data, constants) { late <- data$cycle == "2017-2018"; data$alcohol12[late] <- data$alcohol_ever[late]; data }),
    list(label = "any supplement (DSD010) for vitamin D supplements", sample = "own",
         formula = osteoarthritis ~ vitd_group + age + sex + race4 + education3 + pir3 + bmi + season + alcohol12 + smoking + active + any_supplement + health3),
    list(label = "without vitamin D supplements",
         formula = osteoarthritis ~ vitd_group + age + sex + race4 + education3 + pir3 + bmi + season + alcohol12 + smoking + active + health3),
    list(label = "unweighted", weighted = FALSE)
  ),
  # 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 covariate text groups BMI (under 25, 25 to 30, 30 or more), but Table 3's adjusted estimates are reproduced with BMI continuous (Model 4's three upper groups 1.26, 1.62, and 1.90 here against 1.25, 1.62, and 1.91; Model 3's top group 1.85 against 1.86), while the text's groups give 1.82, 1.73, and 1.80 for the top group in Models 2, 3, and 4 (published 1.93, 1.86, 1.91)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The covariate text defines vitamin D supplement use as taking more than one supplement containing vitamin D in the past 30 days, but Table 1's 46.6% users is the weighted share taking one or more (47.4% here); two or more gives 12.5%."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract gives the subgroup estimate for BMI under 25 as 1.01 (1.04 to 1.08), outside its own interval; the Results give 1.04 (1.00 to 1.08).")
  ),
  choices = list(
    list(choice = "BMI in the models", decision = "continuous", reason = "the covariate text groups BMI (<25, 25-30, 30+), as Table 1 does, but Table 3's numbers show BMI entered continuously: so entered, every adjusted estimate is reproduced (Model 2 1.25, 1.63, 1.93 against 1.23, 1.62, 1.93; Model 3 1.85 against 1.86; Model 4 1.26, 1.62, 1.90 against 1.25, 1.62, 1.91; per 10 nmol/L 1.07), while the groups give 1.82, 1.73, and 1.80 for the top group (a variant; both reproduce the headline)"),
    list(choice = "osteoarthritis", decision = "told of arthritis (MCQ160A) and its type osteoarthritis, with the type codes of each cycle (library arthritis()); other types are not osteoarthritis, an unknown type is missing",
         reason = "the paper's definition; this gives its flow exactly (91,351 participants, 58,284 under 40, 3,702 without 25(OH)D, 3,641 without osteoarthritis information) and its 14.3% prevalence in 2001-2006, which reading MCQ190 with the 2011 codes would make 6.5%"),
    list(choice = "25(OH)D", decision = "LBDVIDMS (2001-2006) and LBXVIDMS (2007 on), groups <50, 50 to <75, 75 to <100, 100+ nmol/L", reason = "the LC-MS/MS-equivalent values the paper describes; it prints <50 and >=100 and leaves 50 and 75 unstated, so all groups are closed below (a variant closes them above)"),
    list(choice = "alcohol in 2017-2018", decision = "12 or more drinks in the past year from ALQ121 and ALQ130", reason = "the 12-drinks question was dropped; this matches Table 1's 72.1% (72.0%) and Table 3, and keeps the cycle (requiring the old question would drop it, leaving 18,858)"),
    list(choice = "vitamin D supplements", decision = "any reported product containing vitamin D in the Dietary Supplement Database, missing when a product has no information",
         reason = "the paper's 'supplement containing vitamin D' (its 'more than one' read as one or more: Table 1 has 46.6% users, this 47.4%); NCHS recommends the database over nutrient totals; with products of unknown content missing, the analytic sample is 21,365 against the paper's 21,334 (counting them as non-users gives 21,579)"),
    list(choice = "vitamin D supplements in 2021-2023", decision = "kept, built the same way from DSQTOT_L and DSQIDS_L, under the dietary day-one weight (WTDRD1) as the plan's weight rule directs",
         reason = "they can be built, and leaving them out moves the paper-cycle estimate from 1.90 to 2.05 (a variant); 2021-2023 asked them by telephone after the first dietary recall (6,754 participants have the data, by the codebook), so its complete-case sample is smaller, and NCHS directs WTDRD1 for any analysis using them, with examination data too; the sensitivity analysis uses the examination weight"),
    list(choice = "recreational activity", decision = "any moderate or vigorous recreational activity (PAD200/PAD320 in 2001-2006, 'unable' as no; PAQ650/PAQ665 from 2007)", reason = "the paper's definition; gives Table 1's 54.8% active"),
    list(choice = "covariate coding", decision = "race: Mexican American, non-Hispanic White (reference), non-Hispanic Black, other (other Hispanic included); education under high school, high school or GED, more; PIR <1.3, 1.3-3.5, 3.5+; season of examination; smoking never, former, current; self-rated health excellent or very good, good, fair or poor",
         reason = "the paper's definitions; Table 1's weighted shares are reproduced to within 0.4 points"),
    list(choice = "weights", decision = "examination weight over nine cycles (WTMEC2YR/9; the phlebotomy weight in 2021-2023)", reason = "the paper weighted 'under the NHANES analytical guidelines' without naming the weight; 25(OH)D is measured at the examination")
  )
)