# Waist circumference and the square root of serum total testosterone in men aged 20 and older,
# NHANES 2013-2016. Zhu et al. (2024), Int J Endocrinol, doi:10.1155/2024/4306797.
# Headline: beta -0.117 (95% CI -0.136 to -0.098) per cm of waist circumference, fully adjusted
# model (Table 2).

# Session of the examination: 0 morning, 1 afternoon, 2 evening (PHDSESN) through 2017-2018;
# 2021-2023 codes 0 morning, 1 afternoon or evening (PHDSESNZ).
row296_session <- function(cycle) {
  d <- component("FASTQX", cycle)
  if (has(d, "PHDSESN")) {
    data.frame(SEQN = d$SEQN, session3 = ifelse(d$PHDSESN %in% 0:2, d$PHDSESN, NA),
               session2 = ifelse(d$PHDSESN %in% 0, 0, ifelse(d$PHDSESN %in% 1:2, 1, NA)))
  } else {
    data.frame(SEQN = d$SEQN, session3 = NA, session2 = ifelse(d$PHDSESNZ %in% 0:1, d$PHDSESNZ, NA))
  }
}

# Hours a week in Table 1's five groups: none, 0.1-0.9, 1.0-3.4, 3.5-5.9, 6 or more.
row296_hours <- function(minutes) {
  cut(minutes / 60, c(-Inf, 0, 0.95, 3.45, 5.95, Inf), labels = c("none", "0.1-0.9", "1.0-3.4", "3.5-5.9", "6+"))
}

# Minutes a week in one GPAQ domain (2007-2018): answer (1 yes, 2 no), days a week, minutes a day.
row296_domain <- function(answer, days, minutes) {
  m <- ifelse(minutes %in% c(7777, 9999), NA, minutes)
  ifelse(answer %in% 2, 0, ifelse(answer %in% 1 & days %in% 1:7, days * m, NA))
}

# One row per SEQN. `steroids` are the Multum categories counted as glucocorticoid use.
row296_build <- function(cycle, steroids = c(301, 296)) {
  replication <- cycle == REPLICATION_CYCLE
  demo <- demographics(cycle)
  alq <- component("ALQ", cycle)
  mcq <- component("MCQ", cycle)
  heq <- component("HEQ", cycle)
  d <- merge_all(demo, component("TST", cycle, c("LBXTST", "LBXEST", "LBXSHBG")), body_measures(cycle), smoking(cycle),
                 alcohol(cycle), alq[, c("SEQN", intersect("ALQ101", names(alq))), drop = FALSE], diabetes_status(cycle),
                 row296_session(cycle), component("COT", cycle, "LBXCOT"),
                 mcq[, c("SEQN", "MCQ230A", "MCQ230B", "MCQ230C", "MCQ230D", "MCQ160M", "MCQ160L")],
                 heq[, c("SEQN", intersect(c("HEQ010", "HEQ030"), names(heq))), drop = FALSE])
  at <- function(frame, column) frame[[column]][match(d$SEQN, frame$SEQN)]
  d$sqrt_tt <- sqrt(d$LBXTST)
  d$race6 <- factor(d$race3, levels = c(1, 2, 3, 4, 6, 7))
  d$education5 <- factor(d$education, levels = 1:5)
  d$diabetes <- factor(d$diabetes, levels = 0:1)
  d$smoking3 <- factor(d$smoking, levels = 1:3)
  d$session <- factor(d$session3, levels = 0:2, labels = c("morning", "afternoon", "evening"))
  d$session_2 <- factor(d$session2, levels = 0:1, labels = c("morning", "afternoon or evening"))
  d$cotinine <- d$LBXCOT
  # Alcohol as Table 1 counts it: nondrinkers never had 12 drinks in any one year (ALQ101, asked
  # through 2016); everyone else falls by drinking days a month in the past year (ALQ120Q/U), those
  # who did not drink in the past year in the lowest group.
  levels4 <- c("none", "1-5", "5-10", "10+")
  days_month <- d$drinking_days / 12
  frequency <- ifelse(is.na(days_month) | days_month < 5, "1-5", ifelse(days_month < 10, "5-10", "10+"))
  d$alcohol <- if (has(d, "ALQ101")) factor(ifelse(d$ALQ101 %in% 2, "none", ifelse(is.na(d$ALQ101), NA, frequency)), levels = levels4) else factor(NA, levels = levels4)
  # The same groups from what every cycle asks: nondrinkers had no drink in the past year.
  d$alcohol_h <- factor(ifelse(d$alcohol3 %in% 1:2, "none", ifelse(d$alcohol3 %in% 3 & !is.na(days_month), frequency, NA)), levels = levels4)
  # Physical activity time: work activity, vigorous plus moderate (GPAQ, 2007-2018), which
  # 2021-2023 does not ask; leisure-time and moderate work activity are for variants.
  d$pa_work <- d$pa_work_moderate <- d$pa_leisure <- row296_hours(NA)
  if (!replication) {
    paq <- component("PAQ", cycle)
    work_v <- row296_domain(paq$PAQ605, paq$PAQ610, paq$PAD615)
    work_m <- row296_domain(paq$PAQ620, paq$PAQ625, paq$PAD630)
    leisure <- row296_domain(paq$PAQ650, paq$PAQ655, paq$PAD660) + row296_domain(paq$PAQ665, paq$PAQ670, paq$PAD675)
    d$pa_work <- row296_hours(work_v + work_m)[match(d$SEQN, paq$SEQN)]
    d$pa_work_moderate <- row296_hours(work_m)[match(d$SEQN, paq$SEQN)]
    d$pa_leisure <- row296_hours(leisure)[match(d$SEQN, paq$SEQN)]
  }
  # Figure 1's exclusions, in its order.
  cancers <- as.matrix(d[, c("MCQ230A", "MCQ230B", "MCQ230C", "MCQ230D")])
  hormone_cancer <- rowSums(matrix(cancers %in% c(14, 22, 30, 36, 37), nrow = nrow(d))) > 0
  thyroid <- d$MCQ160M %in% 1
  hepatitis_b <- d$HEQ010 %in% 1
  hepatitis_c <- if (has(d, "HEQ030")) d$HEQ030 %in% 1 else FALSE
  liver <- d$MCQ160L %in% 1
  d$glucocorticoid <- d$infection <- d$hiv <- NA
  if (!replication) {
    d$glucocorticoid <- at(medication_use(cycle, list(glucocorticoid = steroids)), "glucocorticoid")
    d$infection <- at(component("HSQ", cycle, "HSQ520"), "HSQ520")
    hiv <- component("HIV", cycle)
    d$hiv <- at(hiv, if (has(hiv, "LBDHI")) "LBDHI" else "LBXHIVC")
  }
  base <- d$sex == 1 & d$age >= 20 & !is.na(d$LBXTST) & !is.na(d$LBXEST) & !is.na(d$LBXSHBG) & !is.na(d$waist) & !is.na(d$bmi) &
    !is.na(d$smoking) & !hormone_cancer & !thyroid & !liver & !hepatitis_b
  d$in_population <- base & !is.na(d$alcohol) & !hepatitis_c & !(d$glucocorticoid %in% 1) & !(d$infection %in% 1) & !(d$hiv %in% 1)
  d$in_population_harmonized <- base & !is.na(d$alcohol_h)
  d
}

row296_formula <- sqrt_tt ~ waist + age + race6 + education5 + pir + diabetes + session + cotinine + alcohol + smoking3 + pa_work

association <- list(
  id = "row296", row = 296, doi = "10.1155/2024/4306797",
  cycles = c("2013-2014", "2015-2016"),
  weight = "WTMEC2YR", blood_file = "TST",
  family = "linear", term = "waist",
  published = list(measure = "beta", estimate = -0.117, low = -0.136, high = -0.098, n = 3359, events = NULL,
                   contrast = "per 1 cm of waist circumference; outcome the square root of total testosterone (ng/dL)"),
  left_out = c(
    "exclusion of glucocorticoid users (2021-2023 prescription data name no drugs)",
    "exclusion of recent infection (flu, pneumonia, or ear infection in the past 30 days, HSQ520, not asked in 2021-2023)",
    "exclusion of HIV infection (no HIV test in 2021-2023)",
    "hepatitis C told by a doctor (HEQ030) in the liver disease exclusion (dropped from 2021-2023; any liver condition and hepatitis B told remain)",
    "physical activity (the paper's is time in work activity; 2021-2023 asks only leisure-time activity)",
    "session of the blood draw in three levels: 2021-2023 codes morning versus afternoon or evening, so the harmonized model uses those two levels",
    "nondrinkers as those who never had 12 drinks in any one year (ALQ101, not asked from 2017): the harmonized model counts those with no drink in the past year as nondrinkers, with the same drinking-day groups"
  ),
  build = function(cycle) row296_build(cycle),
  formula = row296_formula,
  formula_harmonized = sqrt_tt ~ waist + age + race6 + education5 + pir + diabetes + session_2 + cotinine + alcohol_h + smoking3,
  variants = list(
    list(label = "nonadjusted (published -0.117, -0.127 to -0.107)", formula = sqrt_tt ~ waist),
    list(label = "age and race (published -0.120, -0.131 to -0.108)", formula = sqrt_tt ~ waist + age + race6),
    list(label = "nonadjusted, everyone in the study population", formula = sqrt_tt ~ waist, sample = "own"),
    list(label = "unweighted", weighted = FALSE),
    list(label = "activity: moderate work only", formula = update(row296_formula, . ~ . - pa_work + pa_work_moderate)),
    list(label = "activity: leisure time instead of work", formula = update(row296_formula, . ~ . - pa_work + pa_leisure)),
    list(label = "harmonized covariates, paper's population", formula = sqrt_tt ~ waist + age + race6 + education5 + pir + diabetes + session_2 + cotinine + alcohol_h + smoking3),
    list(label = "glucocorticoids: systemic only (Multum 301)", build = function(cycle) row296_build(cycle, steroids = 301))
  ),
  # 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 Covariates section and Table 1 group alcohol intake as nondrinker and 1-5, 5-10, and 10+ drinks a month (a drink being 12 oz of beer, a 5 oz glass of wine, or 1.5 oz of liquor), but Table 1's counts (592, 1,816, 331, 620) are reproduced only by drinking days a month over the past year, with nondrinkers those who never had 12 drinks in any one year and past-year abstainers in the 1-5 group (590, 1,816, 330, 621)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2 gives the fully adjusted T/E2 estimate for waist quartile 2 as -0.932 (-0.150, -0.354), an interval that excludes its estimate.")
  ),
  choices = list(
    list(choice = "outcome scale", decision = "square root of total testosterone in ng/dL", reason = "the paper does not say which unit was square-rooted; Table 1 reports ng/dL and Table 2's quartile contrasts fit the ng/dL scale (the extraction's check)"),
    list(choice = "weights", decision = "examination weight (WTMEC2YR) halved for the two cycles; WTPH2YR from TST_L in 2021-2023", reason = "the paper says only that it weighted; testosterone was measured on the examined sample"),
    list(choice = "starting sample", decision = "men 20 and older with total testosterone, estradiol, SHBG, waist, BMI, the 12-drinks-a-year item (ALQ101), and smoking status", reason = "reproduces Figure 1's 5290 men and its missing counts (751, 210, 14, 244) exactly; the 3 men left out of 5290 are those who refused or did not know whether they smoked 100 cigarettes"),
    list(choice = "cancer affecting hormones", decision = "any reported cancer of the prostate, testis, thyroid, liver, or breast (MCQ230A-D)", reason = "prostate, testis, thyroid, and liver reproduce Figure 1's 133 and the following 171 with thyroid problems exactly; no man in the paper's sample reported breast cancer, which is hormone dependent"),
    list(choice = "thyroid disorder", decision = "ever told of a thyroid problem (MCQ160M)", reason = "reproduces Figure 1's 171"),
    list(choice = "liver disease", decision = "ever told of any liver condition (MCQ160L) or of hepatitis B or C (HEQ010, HEQ030)", reason = "reproduces Figure 1's 211, 46 of them hepatitis B or C alone, and the 3556 remaining"),
    list(choice = "glucocorticoid use", decision = "any prescription medicine in the Multum glucocorticoid (301) or inhaled corticosteroid (296) category in the past month", reason = "unstated; the paper counts 43, systemic plus inhaled gives 46 and systemic alone 31; the final sample is 3357 against the paper's 3359"),
    list(choice = "recent infection", decision = "flu, pneumonia, or ear infection in the past 30 days (HSQ520)", reason = "unstated; gives 136 against the paper's 137"),
    list(choice = "HIV/AIDS", decision = "reactive HIV antibody test (LBDHI 2013-2014, LBXHIVC 2015-2016)", reason = "unstated; gives the paper's 17"),
    list(choice = "alcohol intake", decision = "nondrinker if never 12 drinks in any one year (ALQ101); otherwise drinking days a month in the past year under 5, 5 to under 10, or 10 and more, past-year abstainers in the lowest group", reason = "the paper's 'drinks/month' groups are reproduced only this way: 590/1816/330/621 against Table 1's 592/1816/331/620"),
    list(choice = "physical activity time", decision = "hours a week of vigorous plus moderate work activity (GPAQ PAQ605-PAD630) in Table 1's five groups", reason = "unstated; of the GPAQ domains and their sums, work activity's univariate associations with testosterone match Supplementary Table 2 (-0.86, -0.33, -0.85, 0.31 against -0.88, -0.36, -0.96, 0.24) and its small groups match Table 1 (103, 308, 140 against 102, 297, 125), though it has 1718 inactive against 1907"),
    list(choice = "diabetes", decision = "told by a doctor, insulin or diabetes pills, HbA1c 6.5% or more, or fasting glucose 126 mg/dL or more; borderline counts as no", reason = "the paper's definition; 574 against Table 1's 579"),
    list(choice = "session, race, education, smoking", decision = "PHDSESN in three levels; RIDRETH3 in six; DMDEDUC2 in five; never, former, current from SMQ020 and SMQ040", reason = "the paper's categories; their counts match Table 1"),
    list(choice = "missing PIR, cotinine, or work activity", decision = "complete case", reason = "unstated; Figure 1 excludes no one for them; the fully adjusted model keeps 3064 of the 3357 men, nearly all of the 293 dropped for missing PIR")
  )
)