# Sitting 7 hours a day or more (PAD680) and urgency urinary incontinence, adults, NHANES
# 2007-2018. Di, Yuan, Xiang, Wang, and Liao (2024), Heliyon, doi:10.1016/j.heliyon.2024.e27764.
# Headline: OR 1.2 (95% CI 1.1-1.3), p = 0.001, sitting 7 hours or more against less, Model 1
# (age, race, education, income-to-poverty ratio, marital status), whole population (Table 2).

# Drinking days a week. The paper's "drinks per week" groups (Table 1) are those of how often
# someone drank in the past year (2007-2016: ALQ120Q per week, per month of 4 weeks, or per year
# of 52; 2017 on: ALQ121's categories as the library's days a year, over 52), not of drinks.
# Past-year abstainers drink on 0 days. Those items skip lifetime abstainers (fewer than 12
# drinks in life through 2015-2016, ALQ101 and ALQ110 no; never a drink from 2017, ALQ111 no),
# who drink on 0 days too; `as_coded` leaves them missing, as the paper's counts show it did.
row221_drinking_days <- function(cycle, as_coded = FALSE) {
  d <- component("ALQ", cycle)
  if (has(d, "ALQ121")) {
    a <- alcohol(cycle)
    days <- a$drinking_days[match(d$SEQN, a$SEQN)] / 52
    lifetime <- d$ALQ111 %in% 2
  } else {
    q <- ifelse(d$ALQ120Q %in% 0:365, d$ALQ120Q, NA)
    per_week <- ifelse(d$ALQ120U %in% 1, 1, ifelse(d$ALQ120U %in% 2, 1 / 4, ifelse(d$ALQ120U %in% 3, 1 / 52, NA)))
    days <- ifelse(q %in% 0, 0, q * per_week)
    lifetime <- d$ALQ101 %in% 2 & d$ALQ110 %in% 2
  }
  days[lifetime] <- if (as_coded) NA else 0
  data.frame(SEQN = d$SEQN, drinking_days = days)
}

# Any moderate or vigorous recreational activity. 2007-2018 ask yes or no (PAQ665, PAQ650: in a
# typical week, for at least 10 minutes at a time); 2021-2023 ask how often (PAD790Q, PAD810Q),
# and any is yes.
row221_recreation <- function(cycle) {
  d <- component("PAQ", cycle)
  if (has(d, "PAQ665")) return(data.frame(SEQN = d$SEQN, moderate = yes(d$PAQ665), vigorous = yes(d$PAQ650)))
  any <- function(q) ifelse(q %in% 0, 0, ifelse(!is.na(q) & !q %in% c(7777, 9999), 1, NA))
  data.frame(SEQN = d$SEQN, moderate = any(d$PAD790Q), vigorous = any(d$PAD810Q))
}

row221_build <- function(cycle, as_coded = FALSE, keep_pregnant = FALSE) {
  demo <- demographics(cycle)
  dm <- component("DEMO", cycle)
  d <- merge_all(demo, component("KIQ_U", cycle, c("KIQ042", "KIQ044")),
                 leisure_activity(cycle)[, c("SEQN", "sedentary_minutes")],
                 body_measures(cycle)[, c("SEQN", "bmi")], component("SMQ", cycle, "SMQ020"),
                 row221_drinking_days(cycle, as_coded), row221_recreation(cycle),
                 component("MCQ", cycle, "MCQ160C"), diabetes_status(cycle), hypertension_status(cycle))
  stress <- ifelse(d$KIQ042 %in% 1, 1, ifelse(d$KIQ042 %in% 2, 0, NA))
  urge <- ifelse(d$KIQ044 %in% 1, 1, ifelse(d$KIQ044 %in% 2, 0, NA))
  known <- !is.na(stress) & !is.na(urge)
  # Urgency UI in the models is any urge leakage (KIQ044), with or without stress leakage: so
  # coded, with stress UI as any stress leakage and mixed UI as both, the crude, Model 1, and
  # Model 2 estimates of all four outcomes in Table 2's whole-population column are reproduced.
  # Table 1's groups are exclusive instead; `urge_only` is its urgency UI group.
  d$urge <- ifelse(known, urge, NA)
  d$urge_only <- ifelse(known, as.integer(urge == 1 & stress == 0), NA)
  d$sit7 <- as.integer(d$sedentary_minutes >= 420)
  d$race <- factor(d$race, levels = 1:5)
  d$education <- factor(education3(d$education), levels = 1:3)
  d$pir3 <- cut(d$pir, c(-Inf, 1.3, 3.5, Inf), right = FALSE, labels = c("<1.3", "1.3-3.5", ">=3.5"))
  # Six groups as Table 1 shows them (DMDMARTL); 2021-2023 has only DMDMARTZ's three, which the
  # library's `marital` gives in every cycle, missing exactly where the six groups are.
  marital6 <- if (has(dm, "DMDMARTL")) dm$DMDMARTL[match(d$SEQN, dm$SEQN)] else NA
  d$marital6 <- factor(ifelse(marital6 %in% 1:6, marital6, NA), levels = 1:6)
  d$marital3 <- factor(d$marital, levels = 1:3)
  # Model 2's covariates, coded as Table 1 shows them.
  d$bmi4 <- cut(d$bmi, c(-Inf, 20, 25, 30, Inf), labels = c("<=20", "20-25", "25-30", ">30"))
  d$smoker <- factor(yes(d$SMQ020), levels = 0:1)
  d$drinks3 <- cut(d$drinking_days, c(-Inf, 1, 4, Inf), right = FALSE, labels = c("<1", "1-3", ">=4"))
  d$moderate <- factor(d$moderate, levels = 0:1)
  d$vigorous <- factor(d$vigorous, levels = 0:1)
  d$chd <- factor(yes(d$MCQ160C), levels = 0:1)
  d$diabetes <- factor(d$diabetes, levels = 0:1)
  d$hypertension <- factor(d$hypertension, levels = 0:1)
  # Under 20, no yes or no to either leakage question or no sitting time (n = 4863), and any
  # Model 2 covariate unknown (n = 6991). Pregnant women are not in the paper's counts (Table 1's
  # 10,768 women, by race, marital status, and UI group, are those left without them).
  covariates <- !is.na(d$education) & !is.na(d$pir) & !is.na(d$marital3) & !is.na(d$bmi) & !is.na(d$smoker) &
    !is.na(d$drinks3) & !is.na(d$moderate) & !is.na(d$vigorous) & !is.na(d$chd) &
    !is.na(d$diabetes) & !is.na(d$hypertension)
  d$in_population <- d$age >= 20 & (keep_pregnant | !(d$pregnant %in% 1)) & known & !is.na(d$sit7) & covariates
  d
}

row221_model1 <- urge ~ sit7 + age + race + education + pir3 + marital6
row221_model2 <- urge ~ sit7 + age + race + education + pir3 + marital6 + bmi4 + smoker + drinks3 + moderate + vigorous + diabetes + hypertension + chd

association <- list(
  id = "row221", row = 221, doi = "10.1016/j.heliyon.2024.e27764",
  cycles = c("2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTMEC2YR",
  family = "logistic", term = "sit7",
  # Events: any urge leakage, Table 1's urgency (1612 + 1339) and mixed (344 + 1927) UI groups.
  published = list(measure = "OR", estimate = 1.2, low = 1.1, high = 1.3, n = 22916, events = 5222,
                   contrast = "sitting 7 hours a day or more against less than 7 hours"),
  left_out = c("marital status in six groups (2021-2023's DMDMARTZ has three: married or living with a partner; widowed, divorced, or separated; never married); the harmonized version uses those three, which leaves the estimate on the paper's cycles unchanged (1.144 either way)"),
  # The paper's population: its "unknown covariates" step dropped lifetime abstainers, whom the
  # drinking items skip. That population reproduces n = 22,916 and every Table 2 estimate, so the
  # published estimate comes from it; the same items skip the same people in 2021-2023.
  build = function(cycle) row221_build(cycle, as_coded = TRUE),
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "crude (published 1.1, 1.0-1.2)", formula = urge ~ sit7),
    list(label = "Model 2 (published 1.1, 1.0-1.2)", formula = row221_model2),
    list(label = "lifetime abstainers kept, as drinking on 0 days a week", build = function(cycle) row221_build(cycle)),
    list(label = "lifetime abstainers kept, crude", formula = urge ~ sit7, build = function(cycle) row221_build(cycle)),
    list(label = "lifetime abstainers kept, Model 2", formula = row221_model2, build = function(cycle) row221_build(cycle)),
    list(label = "pregnant women kept", build = function(cycle) row221_build(cycle, as_coded = TRUE, keep_pregnant = TRUE)),
    list(label = "urge without stress leakage (Table 1's UUI group)", formula = urge_only ~ sit7 + age + race + education + pir3 + marital6),
    list(label = "marital status left out", formula = urge ~ sit7 + age + race + education + pir3)
  ),
  # 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 = "sample", affects_headline = TRUE, followed = TRUE,
         detail = "The text excludes people with unknown covariates, but its population also lacks the 3,562 lifetime abstainers, whom the alcohol items skip and whose other covariates are known. Without them the population is the paper's 22,916 (12,148 men, 10,768 women) exactly and every Table 2 estimate is reproduced (Model 1 1.18, 1.07-1.30); kept, as drinking on 0 days a week, they give 26,478 people and 1.14 (1.05-1.25)."),
    list(kind = "sample", affects_headline = TRUE, followed = TRUE,
         detail = "The text lists three exclusions (under 20, no answer on UI or sitting time, unknown covariates), but its population also leaves out pregnant women: Table 1's 10,768 women and its counts by race and UI group are matched exactly without them, and keeping them adds 228 (Model 1 1.18 either way)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The text groups alcohol use by drinks a week (under 1, 1 to 3, 4 or more), but Table 1's counts in those groups are those of drinking days a week (men exactly, women within one). Alcohol is in Model 2 and reaches Model 1 only through the complete-case step, where drinking days keep 11 drinkers who gave no number of drinks a day, whom drinks a week would leave out."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The text and Fig. 1 exclude 25,372 people under 20, but the 59,842 participants less the 34,770 adults Fig. 1 next shows, which the data reproduce, leave 25,072."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Results give the headline estimate (OR 1.2, 1.1 to 1.3, p = 0.001) as that of men, but Table 2 prints it in the whole-population column, where it is reproduced here (1.18, 1.07-1.30); Table 2's Model 1 estimate for men is 1.2 (0.99 to 1.4), p = 0.1."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's title calls it weighted, but its percentages are unweighted shares (1,798 of 12,148 men, 14.8%, are Mexican American), as are the weighted proportions with UI that the Results give (2,202 of 12,148 men, 18.1%; 5,860 of 10,768 women, 54.4%).")
  ),
  choices = list(
    list(choice = "outcome", decision = "any urge leakage (KIQ044 yes against no), mixed UI included", reason = "unstated (Table 1's groups are exclusive); on the paper's population this coding reproduces Table 2's urgency UI estimates (crude 1.11, p 0.02; Model 1 1.18, 1.07-1.30, p 0.001; Model 2 1.12, p 0.03), and coding stress UI as any stress leakage and mixed UI as both reproduces theirs; urge without stress leakage gives 1.13 (p 0.02)"),
    list(choice = "weights", decision = "examination weight (WTMEC2YR), each cycle's over 6", reason = "the leakage questions are asked in the examination center; weighted fits reproduce Table 2, unweighted ones do not (crude 1.14, p < 0.001, against the published p 0.02)"),
    list(choice = "lifetime abstainers", decision = "dropped, as the paper's computation dropped them", reason = "the paper's drinking variable left them missing (the alcohol items skip them), so its 'unknown covariates' step dropped 3562 people whose other covariates were known; with them dropped the population is the paper's 22,916 (12,148 men, 10,768 women) exactly and every Table 2 estimate is reproduced (Model 1 1.18, 1.07-1.30). Keeping them, as drinking on 0 days a week, is a variant"),
    list(choice = "pregnant women", decision = "excluded (RIDEXPRG 1)", reason = "unstated, but the paper's counts leave them out: without them its 10,768 women and Table 1's counts by race and UI group are matched exactly; with them there would be 228 more"),
    list(choice = "missing UI or sitting time", decision = "KIQ042 and KIQ044 each yes or no, and PAD680 0-1380 minutes", reason = "gives the paper's n = 4863 and 29,907 exactly"),
    list(choice = "unknown covariates", decision = "complete data on every Model 2 covariate", reason = "with the paper's alcohol coding this gives its 22,916 (12,148 men, 10,768 women) exactly"),
    list(choice = "covariate coding", decision = "age continuous; RIDRETH1 five groups; DMDEDUC2 1-2, 3, 4-5; income-to-poverty ratio under 1.3, 1.3 to 3.5, 3.5 or more; DMDMARTL six groups", reason = "Table 1's categories, whose counts these codings match"),
    list(choice = "drinking", decision = "drinking days a week (under 1, 1 to 3, 4 or more)", reason = "the paper says drinks a week, but Table 1's counts are those of drinking days (men exactly, women within one), not of drinks; only Model 2 and the population step use it"),
    list(choice = "recreational activity in 2021-2023", decision = "any moderate or vigorous leisure-time activity (PAD790Q, PAD810Q above 0)", reason = "2021-2023 asks how often, not yes or no for a typical week of 10-minute bouts; it enters only the complete-data step"),
    list(choice = "Model 2's diabetes and hypertension", decision = "diabetes: told, medication, HbA1c 6.5% or more, or fasting glucose 126 mg/dL or more; hypertension: told, medication, or a mean pressure of 140/90 mmHg or more", reason = "unstated; Table 1's counts (diabetes 2427 men, 1806 women; hypertension 5370, 4497) are far above those told alone and near these definitions' (2309, 1672; 5348, 4482)")
  ),
  formula = row221_model1,
  formula_harmonized = urge ~ sit7 + age + race + education + pir3 + marital3
)