# Triglyceride glucose-body mass index (TyG-BMI) and stress urinary incontinence (KIQ042), adults,
# NHANES 2001-2018. Li, Xie, Tian, Wang, Mo, Yang, and Guo (2024), Lipids Health Dis,
# doi:10.1186/s12944-024-02306-7.
# Headline: OR 2.36 (95% CI 2.03-2.78), TyG-BMI quartile 4 (287.23-679.46) vs quartile 1
# (112.58-204.87), Model 3 (Table 2, SUI).
#
# The population below gives the paper's flow exactly (50,201 adults; 6,365, 24,593, and 492
# excluded; 18,751 left), and the covariates below give its Table 1 counts exactly or within a few
# people, hypertension aside (see choices).

# Alcohol use as the paper codes it: at least 12 drinks in any one year, yes or no. 2001-2002 names
# that item ALD100 and 2003-2016 ALQ101. 2017-2018 and 2021-2023 no longer ask it; there, ever
# having had a drink (ALQ111) is yes and never is no, which gives Table 1's 13,793 yes, 4,907 no,
# and 51 missing exactly (any drinking in the past year, ALQ121, would give 2 to 3 points fewer yes).
row236_alcohol <- function(cycle) {
  d <- component("ALQ", cycle)
  item <- if (has(d, "ALQ101")) d$ALQ101 else if (has(d, "ALD100")) d$ALD100 else d$ALQ111
  data.frame(SEQN = d$SEQN, drinker = ifelse(item %in% 1, 1, ifelse(item %in% 2, 0, NA)))
}

# Vigorous and moderate leisure-time activity, yes or no. 2001-2006 ask about the past 30 days
# (PAD200, PAD320; "unable to do" counts as no), 2007-2018 about a typical week (PAQ650, PAQ665),
# and 2021-2023 how often (PAD810Q, PAD790Q), where any is yes. Work activity is not counted:
# leisure alone gives Table 1's counts exactly, and with work the shares would be 34% and 57%.
row236_activity <- function(cycle) {
  d <- component("PAQ", cycle)
  answer <- function(x) ifelse(x %in% 1, 1, ifelse(x %in% 2:3, 0, NA))
  any <- function(q) ifelse(q %in% 0, 0, ifelse(!is.na(q) & !q %in% c(7777, 9999), 1, NA))
  if (has(d, "PAD200")) return(data.frame(SEQN = d$SEQN, vigorous = answer(d$PAD200), moderate = answer(d$PAD320)))
  if (has(d, "PAQ650")) return(data.frame(SEQN = d$SEQN, vigorous = answer(d$PAQ650), moderate = answer(d$PAQ665)))
  data.frame(SEQN = d$SEQN, vigorous = any(d$PAD810Q), moderate = any(d$PAD790Q))
}

# Covered by health insurance: HID010 in 2001-2004, HIQ011 from 2005 on.
row236_insurance <- function(cycle) {
  d <- component("HIQ", cycle)
  item <- if (has(d, "HIQ011")) d$HIQ011 else d$HID010
  data.frame(SEQN = d$SEQN, insured = ifelse(item %in% 1, 1, ifelse(item %in% 2, 0, NA)))
}

# Yes when any part is met, no when none is and at least one is known.
row236_any <- function(...) {
  parts <- cbind(...)
  ifelse(rowSums(parts, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(parts)) > 0, 0, NA))
}

# Diabetes as the paper defines it: told by a doctor (DIQ010), taking insulin (DIQ050) or diabetes
# pills (DIQ070; DID070 in 2005-2008), or HbA1c over 6.5% (`at_least` counts 6.5% itself). A part a
# person lacks counts as not met.
row236_diabetes <- function(cycle, at_least = FALSE) {
  q <- component("DIQ", cycle)
  pills <- if (has(q, "DIQ070")) q$DIQ070 else q$DID070
  h <- hba1c(cycle)
  a1c <- h$hba1c[match(q$SEQN, h$SEQN)]
  data.frame(SEQN = q$SEQN, diabetes = row236_any(q$DIQ010 == 1, q$DIQ050 == 1, pills == 1, if (at_least) a1c >= 6.5 else a1c > 6.5))
}

# Hypertension with the thresholds the paper prints: told (BPQ020), taking medication (BPQ050A;
# BPQ150 in 2021-2023), or a mean pressure over 140/90 mmHg. hypertension_status counts 140/90 itself.
row236_hypertension_over <- function(cycle) {
  q <- bp_questions(cycle)
  bp <- blood_pressure(cycle)
  i <- match(q$SEQN, bp$SEQN)
  data.frame(SEQN = q$SEQN, hypertension_printed = row236_any(q$told_hypertension == 1, q$bp_medication == 1, bp$sbp[i] > 140, bp$dbp[i] > 90))
}

# Ever told to take medicine for high cholesterol (BPQ090D, 2001-2018; not asked in 2021-2023).
row236_told_cholesterol_medicine <- function(cycle) {
  d <- component("BPQ", cycle)
  data.frame(SEQN = d$SEQN, told_cholesterol_medicine = if (has(d, "BPQ090D")) yes(d$BPQ090D) else NA)
}

ROW236_COVARIATES <- c("education", "marital", "pir3", "smoking", "alcohol", "vigorous", "moderate",
                       "diabetes", "hypertension", "high_cholesterol", "insurance")

row236_build <- function(cycle, cuts = c(204.87, 242.54, 287.30)) {
  kiq <- component("KIQ_U", cycle)
  kiq <- data.frame(SEQN = kiq$SEQN, KIQ042 = kiq$KIQ042, KIQ044 = kiq$KIQ044, KIQ046 = if (has(kiq, "KIQ046")) kiq$KIQ046 else NA)
  diabetes_conventional <- row236_diabetes(cycle, at_least = TRUE)
  d <- merge_all(demographics(cycle), kiq, body_measures(cycle)[, c("SEQN", "bmi", "waist")],
                 triglycerides(cycle)[, c("SEQN", "tg")], fasting_glucose(cycle), smoking(cycle),
                 row236_alcohol(cycle), row236_activity(cycle), row236_insurance(cycle), row236_diabetes(cycle),
                 hypertension_status(cycle), row236_hypertension_over(cycle), bp_questions(cycle),
                 row236_told_cholesterol_medicine(cycle), total_cholesterol(cycle))
  d$diabetes_conventional <- diabetes_conventional$diabetes[match(d$SEQN, diabetes_conventional$SEQN)]
  # Stress UI: leakage with coughing, lifting, or exercise in the past 12 months (KIQ042), with or
  # without urge leakage (KIQ044), as Table 1's overlapping SUI, UUI, and MUI counts show.
  d$sui <- ifelse(d$KIQ042 %in% 1, 1, ifelse(d$KIQ042 %in% 2, 0, NA))
  d$uui <- ifelse(d$KIQ044 %in% 1, 1, ifelse(d$KIQ044 %in% 2, 0, NA))
  # TyG-BMI = ln(triglycerides (mg/dL) x fasting glucose (mg/dL) / 2) x BMI (kg/m2), in quartiles.
  d$tyg_bmi <- tyg(d$tg, d$glucose) * d$bmi
  d$tyg_bmi_q <- cut(d$tyg_bmi, c(-Inf, cuts, Inf), labels = c("Q1", "Q2", "Q3", "Q4"))
  d$sex <- factor(d$sex)
  d$race <- factor(d$race, levels = 1:5)
  d$education <- factor(education3(d$education), levels = 1:3)
  d$marital <- factor(d$marital, 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"))
  d$smoking <- factor(d$smoking, levels = 1:3)
  d$alcohol <- factor(d$drinker, levels = 0:1)
  d$vigorous <- factor(d$vigorous, levels = 0:1)
  d$moderate <- factor(d$moderate, levels = 0:1)
  d$diabetes <- factor(d$diabetes, levels = 0:1)
  d$diabetes_conventional <- factor(d$diabetes_conventional, levels = 0:1)
  d$hypertension <- factor(d$hypertension, levels = 0:1)
  d$hypertension_printed <- factor(d$hypertension_printed, levels = 0:1)
  # High cholesterol: told by a doctor (BPQ080), told to take medicine for it (BPQ090D), or total
  # cholesterol 240 mg/dL or more, which gives Table 1's 7,659 and 11,092 exactly. 2021-2023 asks
  # only whether one now takes such medicine (BPQ101D), so the harmonized version counts taking it
  # (BPQ100D before 2021) in place of being told to.
  d$high_cholesterol <- factor(row236_any(d$told_cholesterol == 1, d$told_cholesterol_medicine == 1, d$tc >= 240), levels = 0:1)
  d$high_cholesterol_taking <- factor(row236_any(d$told_cholesterol == 1, d$cholesterol_medication == 1, d$tc >= 240), levels = 0:1)
  d$insurance <- factor(d$insured, levels = 0:1)
  # Missing covariates kept as their own level, for the variant standing in for the paper's
  # multiple imputation.
  for (v in ROW236_COVARIATES) d[[paste0(v, "_u")]] <- with_unclear(d[[v]], levels(d[[v]]))
  # Under 20 (41,150 in the paper); not answering yes or no to each of the three leakage items
  # (6,365): with physical activity (KIQ042), urge (KIQ044), and nonphysical activity (KIQ046);
  # no BMI, triglycerides, or fasting glucose, and also no waist circumference, which the paper's
  # counts show it required (24,593); pregnant (492).
  d$in_population_harmonized <- d$age >= 20 & !is.na(d$sui) & !is.na(d$uui) & !is.na(d$tyg_bmi) & !is.na(d$waist) &
    !(d$pregnant %in% 1)
  d$in_population <- d$in_population_harmonized & d$KIQ046 %in% 1:2
  d
}

ROW236_MODEL3 <- sui ~ tyg_bmi_q + age + sex + race + education + marital + pir3 + smoking + alcohol + vigorous + moderate +
  diabetes + hypertension + high_cholesterol + insurance

association <- list(
  id = "row236", row = 236, doi = "10.1186/s12944-024-02306-7",
  cycles = c("2001-2002", "2003-2004", "2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  # The exposure needs the fasting triglycerides and glucose, so the fasting subsample weight (in
  # GLU) applies; it stays WTSAF2YR in 2021-2023.
  weight = "WTSAF2YR", blood_file = "GLU",
  family = "logistic", term = "tyg_bmi_qQ4",
  published = list(measure = "OR", estimate = 2.36, low = 2.03, high = 2.78, n = 18751, events = 4352,
                   contrast = "TyG-BMI quartile 4 (287.23-679.46) vs quartile 1 (112.58-204.87)"),
  left_out = c("the requirement that the question on leakage during nonphysical activities (KIQ046) be answered: 2021-2023 does not ask it (on the paper's cycles it removes 20 people)",
               "being told to take cholesterol medicine (BPQ090D) as a criterion of high cholesterol: 2021-2023 asks only whether one takes it (BPQ101D), which the harmonized version uses in every cycle (BPQ100D before 2021)"),
  build = function(cycle) row236_build(cycle),
  formula = ROW236_MODEL3,
  formula_harmonized = update(ROW236_MODEL3, . ~ . - high_cholesterol + high_cholesterol_taking),
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "Model 1, crude (published 1.86, 1.63-2.11)", formula = sui ~ tyg_bmi_q),
    list(label = "Model 1, crude, unweighted", formula = sui ~ tyg_bmi_q, weighted = FALSE),
    list(label = "Model 2, age, sex, race (published 2.46, 2.14-2.84)", formula = sui ~ tyg_bmi_q + age + sex + race),
    list(label = "Model 2, unweighted", formula = sui ~ tyg_bmi_q + age + sex + race, weighted = FALSE),
    list(label = "examination weight (WTMEC2YR)", weight = "WTMEC2YR"),
    list(label = "Q3-Q4 cutpoint as printed (287.23)", build = function(cycle) row236_build(cycle, cuts = c(204.87, 242.54, 287.23))),
    list(label = "diabetes counts HbA1c 6.5% itself", formula = update(ROW236_MODEL3, . ~ . - diabetes + diabetes_conventional)),
    list(label = "hypertension over 140/90 mmHg, as printed", formula = update(ROW236_MODEL3, . ~ . - hypertension + hypertension_printed)),
    list(label = "missing covariates as their own level, for imputation", sample = "own",
         formula = sui ~ tyg_bmi_q + age + sex + race + education_u + marital_u + pir3_u + smoking_u + alcohol_u + vigorous_u +
           moderate_u + diabetes_u + hypertension_u + high_cholesterol_u + insurance_u)
  ),
  # 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 those with incomplete UI data and those lacking BMI, triglycerides, or fasting glucose, but its exclusion counts (6,365, 24,593, and 492 pregnant, leaving 18,751) and Table 1's counts are matched only if waist circumference is also required and the third leakage item (KIQ046, during nonphysical activities), which none of its outcomes uses, is also answered; without these two requirements 19,138 remain."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE,
         detail = "The Results and Table 1 give Q4 as 287.23 to 679.46, but Table 1's equal quartile sizes (4,688, 4,687, 4,688, 4,688) put the Q3-Q4 boundary at 287.30 on the reproduced population, where its SUI, UUI, and MUI counts by quartile are matched exactly; at 287.23 six people move into Q4 (4,682 and 4,693), and Model 3 is 2.42 either way."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "The Results give the sex split (49.45% women) as weighted proportions, but it is unweighted (9,273 of 18,751), as are Table 1's other shares (SUI 4,352 of 18,751, 23.21%).")
  ),
  choices = list(
    list(choice = "cycles", decision = "the nine cycles 2001-2002 to 2017-2018, without the 2017-March 2020 files",
         reason = "the paper's 91,351 participants are the nine cycles' DEMO files exactly, and its 41,150 under 20 leave 50,201 adults, as here"),
    list(choice = "population", decision = "adults answering yes or no to all three leakage items (KIQ042, KIQ044, KIQ046), with BMI, triglycerides, fasting glucose, and waist circumference, not pregnant (RIDEXPRG)",
         reason = "the paper names only the UI items and BMI, triglycerides, and glucose, but only with KIQ046 and waist circumference also required are its steps (6,365, 24,593, 492) and its 18,751 matched, with Table 1's counts (sex, race, education, marital status, income, smoking, alcohol, activity, insurance, UI types, missing counts) exact and its TyG-BMI range (112.58-679.46) too; without them there are 19,138"),
    list(choice = "weights", decision = "fasting subsample weight (WTSAF2YR) over 9 cycles, strata and PSUs",
         reason = "the paper says it weighted but not with which weight; the exposure needs fasting labs, so NCHS's guidelines call for the fasting weight; with it the abstract's UUI and MUI prevalences (19.42%, 9.32%) are matched exactly (SUI 23.44% against 23.59%), and the crude (1.89 against 1.86) and Model 2 (2.56 against 2.46) estimates are reproduced, while the unweighted crude estimate (2.14) is not"),
    list(choice = "missing covariates", decision = "complete cases",
         reason = "the paper imputed them (multiple imputation, no details); keeping everyone, with missing values as their own level, gives 2.36 (2.02-2.77), the published estimate (variant)"),
    list(choice = "stress UI", decision = "KIQ042 yes against no, with or without urge leakage",
         reason = "Table 1's SUI (4352), UUI (4264), and MUI (1970) counts are reproduced exactly with overlapping groups"),
    list(choice = "quartiles", decision = "the paper's cutpoints 204.87 and 242.54, and 287.30 between Q3 and Q4 where it prints 287.23",
         reason = "the paper's quartiles are equal groups of its sample (4688, 4687, 4688, 4688), whose third quartile is 287.30 on the reproduced sample; with 287.23 six people move to Q4 (4682 and 4693), while with 287.30 Table 1's SUI, UUI, and MUI counts by quartile are matched exactly (variant: as printed)"),
    list(choice = "alcohol use", decision = "at least 12 drinks in any one year (ALD100 in 2001-2002, ALQ101 in 2003-2016); ever a drink (ALQ111) in 2017-2018 and 2021-2023",
         reason = "the paper's item was not asked from 2017 on; ever a drink gives Table 1's 13,793 yes, 4,907 no, and 51 missing exactly"),
    list(choice = "vigorous and moderate activity", decision = "leisure-time only: PAD200 and PAD320 (2001-2006), PAQ650 and PAQ665 (2007-2018), any PAD810Q or PAD790Q (2021-2023)",
         reason = "unstated; leisure alone gives Table 1's counts exactly (4,530 and 8,131 yes; 2 and 7 missing)"),
    list(choice = "diabetes", decision = "told, insulin or pills (DID070 in 2005-2008), or HbA1c over 6.5% as printed",
         reason = "gives Table 1's 2,900 within 2; counting 6.5% itself gives 2,978 (variant)"),
    list(choice = "hypertension", decision = "told, taking medication, or a mean of all readings of 140/90 mmHg or more",
         reason = "the paper prints over 140/90; at or over gives 8,092 and over 8,015 against Table 1's 8,182 (neither exact; which readings is unstated), so the conventional threshold, nearer Table 1, is used (variant: as printed)"),
    list(choice = "high cholesterol", decision = "told (BPQ080), told to take medicine (BPQ090D), or total cholesterol 240 mg/dL or more",
         reason = "the paper's sentence is garbled; this reading gives Table 1's 7,659 and 11,092 exactly, while taking medicine (BPQ100D) gives 7,605"),
    list(choice = "covariate coding", decision = "age continuous; RIDRETH1 five groups; DMDEDUC2 1-2, 3, 4-5; marital status in three groups; income-to-poverty ratio under 1.3, 1.3 to 3.5, 3.5 or more; smoking never, former, current; insurance HID010 (2001-2004) or HIQ011",
         reason = "Table 1's categories, whose counts these codings give exactly")
  )
)