# Weekday sleep under 7 hours (against 7 to 9 hours) and the visceral adiposity index (VAI),
# adults, NHANES 2007-2018. Liu et al. (2024), BMJ Open, doi:10.1136/bmjopen-2023-082601.
# Headline: beta 0.15 (95% CI 0.01-0.28) for short against middle sleep, Model III (Table 2).

# Usual sleep on weekdays or workdays, in hours. 2007-2014 asked it directly in whole hours
# (SLD010H, 1-12 with 12 meaning 12 or more; 77 and 99 refused or don't know). 2015-2016 asked it
# directly in half hours (SLD012). 2017-2018 and 2021-2023 derive SLD012 from the usual sleep and
# wake times (SLQ300, SLQ310), in half hours from 3 to 13.5, with 2 for under 3 hours and 14 for 14
# or more. The question and unit are the same throughout, so the hours are used as recorded.
row133_sleep <- function(cycle) {
  s <- component("SLQ", cycle)
  hours <- if (has(s, "SLD012")) s$SLD012 else ifelse(s$SLD010H %in% 1:12, s$SLD010H, NA)
  data.frame(SEQN = s$SEQN, sleep_hours = hours)
}

# Physical activity in minutes a week over the Global Physical Activity Questionnaire's five
# domains (vigorous and moderate work, walking or cycling to get places, vigorous and moderate
# leisure): days a week times minutes a day. Missing if any domain's answer is missing, refused,
# or "don't know". 2021-2023 asks about leisure time only, so it is NA there.
row133_activity <- function(cycle) {
  d <- component("PAQ", cycle)
  if (!has(d, "PAQ605")) return(data.frame(SEQN = d$SEQN, vigorous_minutes = NA_real_, moderate_minutes = NA_real_))
  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))
  }
  data.frame(SEQN = d$SEQN,
             vigorous_minutes = part(d$PAQ605, d$PAQ610, d$PAD615) + part(d$PAQ650, d$PAQ655, d$PAD660),
             moderate_minutes = part(d$PAQ620, d$PAQ625, d$PAD630) + part(d$PAQ635, d$PAQ640, d$PAD645) +
               part(d$PAQ665, d$PAQ670, d$PAD675))
}

# Model III with the seven self-reported conditions as separate terms (a variant).
ROW133_SEPARATE <- vai ~ sleep3 + age + sex + race4 + married + education4 + pir + smoking + alcohol + activity + energy +
  hypertension + diabetes + stroke + heart_attack + heart_failure + chd + cancer

association <- list(
  id = "row133", row = 133, doi = "10.1136/bmjopen-2023-082601",
  cycles = c("2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  # HDL cholesterol (every examinee with blood drawn) and fasting triglycerides make up VAI; HDL_L
  # carries 2021-2023's phlebotomy weight (TRIGLY_L carries only the fasting weight).
  weight = "WTMEC2YR", blood_file = "HDL",
  family = "linear", term = "sleep3short",
  published = list(measure = "beta", estimate = 0.15, low = 0.01, high = 0.28, n = 11252, events = NA,
                   contrast = "weekday sleep under 7 hours vs 7 to 9 hours (difference in mean VAI)"),
  left_out = c(
    "physical activity (inactive, insufficient, sufficient, from minutes a week over the GPAQ's work, transport, and leisure domains): 2021-2023 asks only about leisure-time activity, so the harmonized model leaves it out",
    "married vs unmarried: 2021-2023 records marital status only as married or living with a partner, widowed, divorced, or separated, and never married (DMDMARTZ), so the harmonized model codes married or living with a partner against the rest",
    "Healthy Eating Index-2015 (a completeness requirement only, not in the model): computing it needs USDA food pattern equivalents from outside NHANES, so neither version computes it; both require the reliable day-1 recall it is computed from, which energy intake requires anyway"
  ),
  build = function(cycle) {
    demo <- demographics(cycle)
    dm <- component("DEMO", cycle)
    alq <- alcohol(cycle)
    alq$alcohol5 <- alcohol5(alq, demo$sex[match(alq$SEQN, demo$SEQN)])
    mcq <- medical_conditions(cycle)
    d <- merge_all(demo, row133_sleep(cycle), body_measures(cycle)[, c("SEQN", "waist", "bmi")],
                   triglycerides(cycle)[, c("SEQN", "tg", "fasted")], hdl_cholesterol(cycle), smoking(cycle),
                   alq[, c("SEQN", "alcohol5")], row133_activity(cycle), dietary_totals(cycle, "KCAL", days = 1),
                   mcq[, c("SEQN", "heart_failure", "chd", "heart_attack", "stroke", "cancer")],
                   bp_questions(cycle)[, c("SEQN", "told_hypertension")], component("DIQ", cycle, "DIQ010"))
    # VAI (Amato et al. 2010), with the fasting subsample's triglycerides and HDL cholesterol in
    # mmol/L; every triglyceride value in the files, as the paper's counts show, and a copy for
    # those NCHS counts as fasted (positive fasting weight) only (a variant).
    d$vai <- vai(d$sex, d$waist, d$bmi, tg_mmol(d$tg), cholesterol_mmol(d$hdl))
    d$vai_fasted <- ifelse(d$fasted %in% 1, d$vai, NA)
    d$sleep3 <- factor(ifelse(d$sleep_hours < 7, "short", ifelse(d$sleep_hours <= 9, "middle", "long")),
                       levels = c("middle", "short", "long"))
    d$sleep3_9long <- factor(ifelse(d$sleep_hours < 7, "short", ifelse(d$sleep_hours < 9, "middle", "long")),
                             levels = c("middle", "short", "long"))
    d$sex <- factor(d$sex)
    d$race4 <- factor(ifelse(d$race %in% 3, "white", ifelse(d$race %in% 1, "mexican", ifelse(d$race %in% 4, "black", "other"))),
                      levels = c("white", "mexican", "black", "other"))
    # Married (DMDMARTL 1) against widowed, divorced, separated, never married, or living with a
    # partner, 2007-2018; 2021-2023 can't separate married from living with a partner.
    d$married <- if (has(dm, "DMDMARTL")) {
      code <- dm$DMDMARTL[match(d$SEQN, dm$SEQN)]
      factor(ifelse(code %in% 1, "married", ifelse(code %in% 2:6, "unmarried", NA)), levels = c("unmarried", "married"))
    } else factor(rep(NA, nrow(d)), levels = c("unmarried", "married"))
    d$married_or_partner <- factor(ifelse(d$marital %in% 1, "married", ifelse(d$marital %in% 2:3, "unmarried", NA)),
                                   levels = c("unmarried", "married"))
    d$education4 <- factor(ifelse(d$education %in% 1:2, "grade or less", ifelse(d$education %in% 3, "high school",
                             ifelse(d$education %in% 4, "some college", ifelse(d$education %in% 5, "college or more", NA)))),
                           levels = c("grade or less", "high school", "some college", "college or more"))
    d$smoking <- factor(d$smoking)
    d$alcohol <- factor(d$alcohol5)
    # Moderate-equivalent minutes (vigorous minutes count twice): none is inactive, under 150 a
    # week insufficient, 150 or more sufficient.
    equivalent <- d$moderate_minutes + 2 * d$vigorous_minutes
    d$activity <- factor(ifelse(equivalent == 0, "inactive", ifelse(equivalent < 150, "insufficient", "sufficient")),
                         levels = c("inactive", "insufficient", "sufficient"))
    d$energy <- d$KCAL
    d$hypertension <- factor(d$told_hypertension)
    # Told by a doctor they have diabetes (DIQ010 1) or not (2); borderline (3), refused, and
    # don't know are missing.
    told_diabetes <- ifelse(d$DIQ010 %in% 1, 1, ifelse(d$DIQ010 %in% 2, 0, NA))
    d$diabetes <- factor(told_diabetes)
    for (v in c("stroke", "heart_attack", "heart_failure", "chd", "cancer")) d[[v]] <- factor(d[[v]])
    # Self-reported chronic disease: told of any of the seven conditions (yes or no to each),
    # missing if any answer is.
    conditions <- cbind(d$told_hypertension, told_diabetes, as.numeric(as.character(d$stroke)),
                        as.numeric(as.character(d$heart_attack)), as.numeric(as.character(d$heart_failure)),
                        as.numeric(as.character(d$chd)), as.numeric(as.character(d$cancer)))
    d$chronic_disease <- factor(as.integer(rowSums(conditions) > 0))
    conditions[, 2] <- ifelse(d$DIQ010 %in% 1, 1, ifelse(d$DIQ010 %in% 2:3, 0, NA))
    d$chronic_disease_borderline_no <- factor(as.integer(rowSums(conditions) > 0))
    # Participants 16 and older (the sleep questionnaire's age, which the flow chart's counts use;
    # the text says 18) with VAI, weekday sleep hours, and a reliable day-1 dietary recall (which
    # HEI-2015 and energy intake need); the model's covariates are then complete cases. Marital
    # status and education are asked from age 20, so the analytic sample is 20 and older either way.
    d$in_population <- d$age >= 16 & !is.na(d$vai) & !is.na(d$sleep_hours) & d$recall_day1 %in% 1
    d
  },
  formula = vai ~ sleep3 + age + sex + race4 + married + education4 + pir + smoking + alcohol + activity + energy + chronic_disease,
  formula_harmonized = vai ~ sleep3 + age + sex + race4 + married_or_partner + education4 + pir + smoking + alcohol + energy + chronic_disease,
  variants = list(
    list(label = "seven conditions as separate terms", formula = ROW133_SEPARATE),
    list(label = "unweighted", weighted = FALSE),
    list(label = "only those NCHS counts as fasted", sample = "own",
         formula = vai_fasted ~ sleep3 + age + sex + race4 + married + education4 + pir + smoking + alcohol + activity + energy + chronic_disease),
    list(label = "fasting subsample weight (WTSAF2YR)", weight = "WTSAF2YR", weight_file = "TRIGLY"),
    list(label = "Model I, crude (published 0.20, 0.06-0.33)", formula = vai ~ sleep3),
    list(label = "Model I, crude, unweighted", formula = vai ~ sleep3, weighted = FALSE),
    list(label = "Model II, age (published 0.21, 0.07-0.34)", formula = vai ~ sleep3 + age),
    list(label = "Model II, age, unweighted", formula = vai ~ sleep3 + age, weighted = FALSE),
    list(label = "borderline diabetes counted as no diabetes", sample = "own",
         formula = vai ~ sleep3 + age + sex + race4 + married + education4 + pir + smoking + alcohol + activity + energy + chronic_disease_borderline_no),
    list(label = "middle 7 to under 9 h, long 9 h or more", term = "sleep3_9longshort",
         formula = vai ~ sleep3_9long + age + sex + race4 + married + education4 + pir + smoking + alcohol + activity + energy + chronic_disease)
  ),
  # 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 = FALSE, followed = NA,
         detail = "Figure 1 labels its 38,562 starting participants as aged 18 or older, but that count and its 15,786 with VAI and sleep data are exactly those aged 16 or older (15,012 at 18 or older). The analytic sample is 20 and older either way (the Results give ages 20 to 80), since marital status and education are asked from 20."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's mean sleep in the middle and long groups (7.54 and 9.51 hours) fits 9 hours counted as long, but its group sizes (3,885, 6,823, 544) and mean VAI fit the stated groups, with 9 hours in the middle (3,865, 6,868, 585 here); with 9 hours counted as long, the long group would be 2.3 times as large."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2 prints Model II's long-sleep estimate as 0.20 (0.09 to 0.49) with P = 0.17, which needs a negative lower limit (Model I's is -0.08)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's 'Never' drinking row repeats its 'Heavy' row (15.27, 11.09, 15.17%)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The abstract and Results give the participants' mean sleep and VAI as 7.05 hours and 2.03, which match the analytic sample's unweighted means (7.07 and 2.02 here; weighted 7.15 and 1.98), though the Methods say continuous variables are presented as weighted means.")
  ),
  choices = list(
    list(choice = "weights", decision = "examination weight, WTMEC2YR divided by 6 (WTPH2YR from HDL_L in 2021-2023)",
         reason = "as the paper states; the weighted crude and age-adjusted fits reproduce Models I and II (0.195 and 0.203 against 0.20 and 0.21, with matching intervals) and unweighted ones do not (0.13). NCHS would weight fasting triglycerides with the fasting subsample weight (variant)"),
    list(choice = "age", decision = "16 and older (the text says 18)",
         reason = "the flow chart's 38,562 'aged 18 or older' and 15,786 with VAI and sleep are exactly the counts at 16 and older (the sleep questionnaire's age) in these cycles (15,012 at 18 and older); 16 and 18 give the same analytic sample, because marital status and education are asked from 20 (Results: ages 20 to 80)"),
    list(choice = "triglycerides for VAI", decision = "the fasting subsample's values (TRIGLY files), in mmol/L, including the 5% NCHS doesn't count as fasted (fasting weight 0)",
         reason = "unstated; every value in the files gives the flow chart's 15,786 with VAI and sleep exactly (14,928 counted as fasted); those counted as fasted only is a variant"),
    list(choice = "sleep hours across cycles", decision = "SLD010H (whole hours, 2007-2014) and SLD012 (half hours, asked in 2015-2016 and derived from sleep and wake times in 2017-2018, as in 2021-2023), used as recorded",
         reason = "the same question (usual sleep at night on weekdays or workdays) in the same unit throughout"),
    list(choice = "sleep groups", decision = "under 7 hours short, 7 to 9 inclusive middle (reference), over 9 long, so 6.5 is short and 9.5 long",
         reason = "the paper's labels (<7, 7-9, >9); Table 1's group sizes (3,885, 6,823, 544) fit them (3,865, 6,868, 585), and so do the group VAI means (2.11, 1.92, 2.12 against 2.13, 1.92, 2.12). Table 1's group sleep means (7.54, 9.51) fit 9 hours in the long group instead, which would make it 2.3 times larger (variant)"),
    list(choice = "self-reported chronic diseases", decision = "one covariate: told of any of hypertension, diabetes, stroke, heart attack, congestive heart failure, coronary heart disease, or cancer",
         reason = "the Methods and Table 2's note adjust for 'self-reported chronic diseases, including' the seven, which reads as one covariate or seven; one term reproduces Model III's short-sleep estimate (0.140 against 0.15) and its long-sleep estimate (0.015, -0.25 to 0.28, against 0.01, -0.26 to 0.29), and seven terms give 0.119 and -0.040 (variant)"),
    list(choice = "diabetes", decision = "told by a doctor (DIQ010 1) or not (2); borderline missing",
         reason = "the paper codes yes or no with no borderline group and excludes missing self-reported diseases; this gives 11,318 participants (published 11,252) against 11,602 with borderline counted as no (variant, 0.140 either way)"),
    list(choice = "marital status", decision = "married (DMDMARTL 1) against all others, living with a partner included",
         reason = "Table 1's married shares (53.1, 58.8, 41.0) match married alone (52.8, 58.6, 40.5), not married or living with a partner (62.0, 66.2, 54.5)"),
    list(choice = "education, race, smoking", decision = "DMDEDUC2 1-2, 3, 4, 5; RIDRETH1 white, Mexican American, black, other (other Hispanic and other or multiracial); never, former, current from SMQ020 and SMQ040",
         reason = "the paper's groups; Table 1's weighted shares match each within half a point in the short and middle groups and within three points in the small long group"),
    list(choice = "energy intake", decision = "day-1 recall (DR1TKCAL, reliable recall), kcal, continuous",
         reason = "Table 1's means (2,248, 2,176, 1,995) match day 1 (2,252, 2,181, 1,942), not the two-day mean (about 2,140, 2,090)"),
    list(choice = "physical activity", decision = "moderate-equivalent minutes a week (vigorous counted twice) over the GPAQ's five domains: none inactive, 1-149 insufficient, 150 or more sufficient; missing if any domain is",
         reason = "Table 1's minutes a week (1,144, 895, 752) match these (1,154, 904, 790); the 150-minute cutpoint is the physical activity guidelines', the convention for these labels"),
    list(choice = "drinking status", decision = "never, former, mild, moderate, heavy by Rattan and colleagues' definitions (library alcohol and alcohol5, which read 2007-2010's binge items for 5 or more drinks, ALQ140Q/U and ALQ150, in place of the later ones)",
         reason = "unstated; the five labels are this convention's. Table 1 can't check it (its 'Never' row repeats 'Heavy'), and its former (24.0, 19.6, 23.0) and heavy (15.3, 11.1, 15.2) shares match no definition tried (these give about 17, 13, 24 and 28, 24, 30); alcohol coding moves the estimate by under 0.01"),
    list(choice = "covariate form and missing data", decision = "age, poverty income ratio (INDFMPIR), and energy continuous; complete-case analysis",
         reason = "Table 1 reports them as means; the paper excludes participants with missing covariates")
  )
)