# A Body Shape Index (ABSI) and prostate cancer in men aged 40 and older, NHANES 2001-2018.
# Liu et al. (2024), Int Urol Nephrol, doi:10.1007/s11255-023-03917-2.
# Headline: OR 1.91 (95% CI 1.12-3.27), ABSI quartile 4 vs quartile 1, Model III (Table 2).

# The cycles that gave the oral glucose tolerance test (OGTT_D to OGTT_I).
ROW220_OGTT_CYCLES <- c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016")

# Diabetes as the paper defines it: told by a doctor, HbA1c 6.5% or more, fasting glucose 7.0
# mmol/L (126 mg/dL) or more, two-hour OGTT glucose 11.1 mmol/L (200 mg/dL) or more, or taking
# insulin or diabetes pills. 2005-2006 names the pills item DID070. A part a person lacks counts as
# not met; with every part unknown, diabetes is unknown. Returns diabetes with the OGTT where the
# cycle gave it, and diabetes_no_ogtt without it.
row220_diabetes <- function(cycle) {
  q <- component("DIQ", cycle)
  pills <- if (has(q, "DIQ070")) q$DIQ070 else if (has(q, "DID070")) q$DID070 else NA
  h <- hba1c(cycle)
  g <- fasting_glucose(cycle)
  parts <- cbind(ifelse(q$DIQ010 %in% 1, 1, ifelse(q$DIQ010 %in% 2:3, 0, NA)),
                 ifelse(q$DIQ050 %in% 1 | pills %in% 1, 1, ifelse(q$DIQ050 %in% 2 | pills %in% 2, 0, NA)),
                 as.integer(h$hba1c[match(q$SEQN, h$SEQN)] >= 6.5),
                 as.integer(g$glucose[match(q$SEQN, g$SEQN)] >= 126))
  status <- function(parts) ifelse(rowSums(parts == 1, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(parts)) > 0, 0, NA))
  out <- data.frame(SEQN = q$SEQN, diabetes_no_ogtt = status(parts))
  if (cycle %in% ROW220_OGTT_CYCLES) {
    o <- component("OGTT", cycle, "LBXGLT")
    parts <- cbind(parts, as.integer(o$LBXGLT[match(q$SEQN, o$SEQN)] >= 200))
  }
  out$diabetes <- status(parts)
  out
}

# Hypertension: told by a doctor (BPQ020), taking prescribed medicine for it (BPQ050A; BPQ150 in
# 2021-2023), or a mean systolic pressure of 140 mmHg or more or diastolic of `dbp_at` or more over
# the second and third readings. A part a person lacks counts as not met.
row220_hypertension <- function(cycle, dbp_at) {
  q <- bp_questions(cycle)
  bp <- blood_pressure(cycle, readings = 2:3)
  flags <- cbind(q$told_hypertension == 1, q$bp_medication == 1,
                 bp$sbp[match(q$SEQN, bp$SEQN)] >= 140, bp$dbp[match(q$SEQN, bp$SEQN)] >= dbp_at)
  data.frame(SEQN = q$SEQN, hypertension = ifelse(rowSums(flags, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(flags)) > 0, 0, NA)))
}

# Drinking status in the groups of Rattan et al. (2022), which the paper cites: never (fewer than 12
# drinks in life), heavy (on drinking days `heavy_at` drinks or more, which Rattan sets at 4 for men,
# or 5 or more drinks on one day on 5 or more days a month, 60 or more days a year), and mild
# (everyone else, former drinkers included). 2001-2016 ask about 12 drinks in a year (ALD100 in
# 2001-2002, ALQ101 after) and in life (ALQ110), how often one drank in the past year (ALQ120Q/U),
# drinks on a drinking day (ALQ130), and days with 5 or more drinks (ALQ140Q/U to 2010, ALQ141Q/U
# after; 5 or more for men). 2017 on drops the 12-drink items, so never is never having had a drink
# (ALQ111); ALQ121 and ALQ142 ask how often in categories (ALQ142 1-4, twice a week or more, is 5
# or more days a month; 5, once a week, is not). A heavy-drinking criterion a person lacks counts as
# not met. Also returns former (drank, but not in the past year).
row220_drinking <- function(cycle, heavy_at = 4) {
  d <- component("ALQ", cycle)
  answer <- function(x) ifelse(x %in% c(77, 99, 777, 999), NA, x)
  drinks <- answer(d$ALQ130)
  if (has(d, "ALQ111")) {
    never <- d$ALQ111 %in% 2
    ever <- d$ALQ111 %in% 1
    past_year_none <- d$ALQ121 %in% 0
    binge <- d$ALQ142 %in% 1:4
  } else {
    twelve_in_a_year <- if (has(d, "ALQ101")) d$ALQ101 else d$ALD100
    never <- d$ALQ110 %in% 2
    ever <- twelve_in_a_year %in% 1 | d$ALQ110 %in% 1
    past_year_none <- d$ALQ120Q %in% 0
    days <- answer(if (has(d, "ALQ141Q")) d$ALQ141Q else d$ALQ140Q)
    unit <- if (has(d, "ALQ141U")) d$ALQ141U else d$ALQ140U
    per_year <- days * ifelse(unit %in% 1, 52, ifelse(unit %in% 2, 12, ifelse(unit %in% 3, 1, NA)))
    binge <- (per_year >= 60) %in% TRUE
  }
  heavy <- (drinks >= heavy_at) %in% TRUE | binge
  data.frame(SEQN = d$SEQN, drinking = ifelse(never, "never", ifelse(ever & heavy, "heavy", ifelse(ever, "mild", NA))),
             former = as.integer(ever & past_year_none))
}

# Prostate cancer: ever told of cancer (MCQ220) with prostate (code 30) among the first three kinds
# (MCQ230A-C). Men with any other kind, more than three kinds, or a kind they could not name are
# left out (missing), as the paper keeps only men with no tumor history but prostate cancer; men
# never told of cancer are the controls.
row220_prostate <- function(cycle) {
  d <- component("MCQ", cycle)
  kinds <- as.matrix(d[, c("MCQ230A", "MCQ230B", "MCQ230C", "MCQ230D")])
  prostate <- rowSums(kinds[, 1:3] == 30, na.rm = TRUE) > 0
  other <- rowSums(!is.na(kinds) & kinds != 30) > 0
  data.frame(SEQN = d$SEQN, pca = ifelse(d$MCQ220 %in% 2, 0, ifelse(d$MCQ220 %in% 1 & prostate & !other, 1, NA)))
}

# Table 2's models.
ROW220_MODEL_II <- pca ~ absi_q + age3 + race + education2 + pir3 + living
ROW220_MODEL_III <- update(ROW220_MODEL_II, . ~ . + bmi3 + drinking + smoking + hypertension + diabetes)

association <- list(
  id = "row220", row = 220, doi = "10.1007/s11255-023-03917-2",
  cycles = c("2001-2002", "2003-2004", "2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  # Body measures come from the examination, and diabetes uses HbA1c.
  weight = "WTMEC2YR", blood_file = "GHB",
  family = "logistic", term = "absi_qQ4",
  published = list(measure = "OR", estimate = 1.91, low = 1.12, high = 3.27, n = 11013, events = 492,
                   contrast = "ABSI (x 1000) quartile 4 (86.42 or more) vs quartile 1 (below 80.94)"),
  left_out = c("the two-hour oral glucose tolerance test from the diabetes definition: 2021-2023 gave none, so diabetes there is told, HbA1c, fasting glucose, or medication (the harmonized model uses that definition in every cycle)",
               "the fewer-than-12-drinks-in-life item behind never drinkers (asked 2001-2016): 2021-2023, like 2017-2018, asks only whether one ever had a drink (ALQ111), so never drinkers there are those who never had one, as in the paper's own 2017-2018 data"),
  build = function(cycle) {
    d <- merge_all(demographics(cycle), body_measures(cycle), row220_prostate(cycle), row220_drinking(cycle),
                   smoking(cycle), row220_diabetes(cycle))
    alq3 <- row220_drinking(cycle, heavy_at = 3)
    d$drinking_3 <- alq3$drinking[match(d$SEQN, alq3$SEQN)]
    h90 <- row220_hypertension(cycle, dbp_at = 90)
    h80 <- row220_hypertension(cycle, dbp_at = 80)
    d$hypertension <- h90$hypertension[match(d$SEQN, h90$SEQN)]
    d$hypertension_80 <- h80$hypertension[match(d$SEQN, h80$SEQN)]
    # ABSI = waist (m) x height (m)^(5/6) x weight (kg)^(-2/3), times 1000.
    d$absi <- 1000 * (d$waist / 100) * (d$height / 100)^(5 / 6) * d$weight^(-2 / 3)
    d$absi_q <- cut(d$absi, c(-Inf, 80.94, 83.66, 86.42, Inf), right = FALSE, labels = c("Q1", "Q2", "Q3", "Q4"))
    d$age3 <- cut(d$age, c(-Inf, 60, 80, Inf), right = FALSE, labels = c("<60", "60-79", ">=80"))
    d$race <- factor(d$race)
    # Education as Table 1 splits it (the Methods split below and at high school instead).
    d$education2 <- factor(ifelse(d$education %in% 1:3, "high school or less", ifelse(d$education %in% 4:5, "above high school", NA)),
                           levels = c("high school or less", "above high school"))
    d$education2_methods <- factor(ifelse(d$education %in% 1:2, "less than high school", ifelse(d$education %in% 3:5, "high school or above", NA)),
                                   levels = c("less than high school", "high school or above"))
    # Income: low 1 or less, middle over 1 to 3, high over 3, as Table 1's counts show.
    d$pir3 <- cut(d$pir, c(-Inf, 1, 3, Inf), labels = c("low", "middle", "high"))
    d$living <- factor(ifelse(d$marital %in% 1, "with a partner", ifelse(d$marital %in% 2:3, "alone", NA)), levels = c("with a partner", "alone"))
    d$bmi3 <- cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("<25", "25-29.9", ">=30"))
    d$drinking <- factor(d$drinking, levels = c("never", "mild", "heavy"))
    d$drinking_3 <- factor(d$drinking_3, levels = c("never", "mild", "heavy"))
    d$smoking <- factor(d$smoking)
    d$hypertension <- factor(d$hypertension)
    d$hypertension_80 <- factor(d$hypertension_80)
    d$diabetes <- factor(d$diabetes)
    d$diabetes_no_ogtt <- factor(d$diabetes_no_ogtt)
    d$in_population <- d$sex == 1 & d$age >= 40 & !is.na(d$pca) & !is.na(d$absi)
    d
  },
  # The paper's computation dropped former drinkers in 2017-2018, whose drinks-a-day and binge items
  # (ALQ130, ALQ142) are skipped; that reproduces its n of 11,013. 2021-2023 asks the same items
  # with the same skips, so the same step applies there.
  derive = function(data, constants) {
    data$in_population <- data$in_population & !(data$cycle %in% c("2017-2018", REPLICATION_CYCLE) & data$former %in% 1)
    data
  },
  formula = ROW220_MODEL_III,
  formula_harmonized = update(ROW220_MODEL_III, . ~ . - diabetes + diabetes_no_ogtt),
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "Model I (published 6.43, 3.96-10.42)", formula = pca ~ absi_q),
    list(label = "Model II (published 1.92, 1.15-3.23)", formula = ROW220_MODEL_II),
    list(label = "2017-2018 former drinkers kept, as mild drinkers", derive = function(data, constants) data),
    list(label = "heavy drinking at 3 drinks a day, as printed", formula = update(ROW220_MODEL_III, . ~ . - drinking + drinking_3)),
    list(label = "hypertension at diastolic 80 mmHg, as printed", formula = update(ROW220_MODEL_III, . ~ . - hypertension + hypertension_80)),
    list(label = "education as the Methods split it", formula = update(ROW220_MODEL_III, . ~ . - education2 + education2_methods)),
    list(label = "age continuous", formula = update(ROW220_MODEL_III, . ~ . - age3 + age)),
    list(label = "ABSI continuous, per unit (published 1.05, 1.02-1.08)", formula = update(ROW220_MODEL_III, . ~ . - absi_q + absi), term = "absi")
  ),
  # 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 paper's definitions class men who drank before but not in the past year as mild drinkers, but its n of 11,013 (492 with prostate cancer) and Table 1's margins are matched (11,022 and 493 here, every margin within 20 men) only if the 333 such men of 2017-2018, whose drinks-a-day and binge items (ALQ130, ALQ142) are skipped, are dropped. Kept as mild drinkers, they give 11,355 men (509 with prostate cancer) and Model III 1.90 (1.11-3.27), against 1.92 (1.12-3.29) without them."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods define heavy drinking as 3 or more drinks a day, but Table 1's 2,048 heavy drinkers (31 with prostate cancer) are those of 4 or more for men (2,068 and 33 here), not of 3 (3,164 and 65); Model III is 1.92 either way."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods define hypertension by a diastolic pressure of 80 mmHg or more, but Table 1's 5,909 hypertensive men (358 with prostate cancer) are matched exactly with 90 or more over the second and third readings, while 80 gives 6,855 (363); Model III is 1.92 with 90 and 1.94 with 80."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods split education into less than high school and high school or above, but Table 1's counts (5,690 and 5,323) are those of its own split, high school or less and above high school (5,693 and 5,329 here), not the Methods' (3,141 and 7,881); Model III is 1.92 with Table 1's split and 1.93 with the Methods'."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods put an income-to-poverty ratio of exactly 3 in the high group, but Table 1's counts (1,898, 4,525, and 4,590) are matched within 3 men in each group with it in the middle group, and are off by 7 and 13 in the middle and high groups with it in the high group (10 men have exactly 3).")
  ),
  choices = list(
    list(choice = "weights", decision = "examination weight (WTMEC2YR) over nine cycles, each divided by nine, with strata and PSUs; WTPH2YR in 2021-2023",
         reason = "the paper weights as CDC directs without naming the weight; ABSI comes from the examination; NCHS uses 2-year weights for 2001-2002 when 1999-2000 is not pooled; weighted confidence intervals match the published widths, unweighted ones are far narrower"),
    list(choice = "quartile cutpoints", decision = "the published 80.94, 83.66, and 86.42, each opening the higher quartile",
         reason = "Table 2's row labels; boundary side unstated"),
    list(choice = "prostate cancer and other cancers", decision = "cases told of cancer with prostate among the first three kinds; men with any other kind, more than three kinds, or a kind not named are left out; controls never told of cancer",
         reason = "the paper keeps men with no tumor history but prostate cancer; Fig. 1's 91,351, 45,010, and 16,127 are reproduced exactly"),
    list(choice = "drinking status", decision = "Rattan et al.'s groups for men: heavy at 4 or more drinks on a drinking day or 5 or more drinks on 5 or more days a month; never fewer than 12 drinks in life (never a drink from 2017 on); mild otherwise, former drinkers included",
         reason = "the text prints 3 drinks, Rattan's threshold for women, but cites Rattan (ref. 20), and Table 1's 2,048 heavy drinkers (31 with prostate cancer) fit 4 (2,068 and 33 here) rather than 3 (3,164); its 801 never drinkers (47 with prostate cancer) match exactly; 3 drinks is a variant"),
    list(choice = "2017-2018 former drinkers", decision = "dropped, as the paper's computation dropped them (and in 2021-2023, whose drinking items skip them the same way)",
         reason = "its n (11,013, 492 with prostate cancer) and Table 1's margins are matched within 20 men (most within 5) only if the 333 men of 2017-2018 who drank before but not in the past year, whose ALQ130 and ALQ142 are skipped, are dropped, as complete-case coding of those items does; keeping them as mild drinkers, as the paper's definitions classify them, is a variant"),
    list(choice = "hypertension", decision = "told, taking prescribed medicine, or a mean systolic pressure of 140 mmHg or more or diastolic of 90 or more over the second and third readings",
         reason = "the text prints diastolic 80 but cites the 2018 ESC/ESH guideline (ref. 21), whose threshold is 90; on the sample the paper ran (2017-2018 former drinkers dropped), Table 1's 5,909 hypertensive men (358 with prostate cancer) are matched exactly with 90 over readings 2 and 3 (5,949 over all readings; 6,855 with 80); 80 is a variant"),
    list(choice = "diabetes", decision = "told, HbA1c 6.5% or more, fasting glucose (LBXGLU) 126 mg/dL or more, two-hour OGTT glucose 200 mg/dL or more where the cycle gave the test (2005-2016), or insulin or diabetes pills; no random glucose",
         reason = "the paper's list; NHANES has no random glucose measure (refrigerated serum glucose would add 3 men); Table 1's 2,752 with diabetes (149 with prostate cancer) compare with 2,747 (148) on the sample the paper ran; fasting glucose as recorded, as is usual in this literature"),
    list(choice = "income to poverty ratio", decision = "low 1 or less, middle over 1 to 3, high over 3",
         reason = "the text puts 3 itself in high, but Table 1's counts fit R's default intervals (10 men have exactly 3)"),
    list(choice = "education", decision = "Table 1's high school or less versus above high school",
         reason = "the Methods say less than high school versus high school or above, but Table 1's counts (5,690 and 5,323) fit its own split; the Methods' split is a variant"),
    list(choice = "age", decision = "three groups: under 60, 60 to 79, 80 or older", reason = "the Methods divide age into these groups; continuous age is a variant"),
    list(choice = "living status", decision = "married or living with a partner versus widowed, divorced, separated, or never married (DMDMARTL; DMDMARTZ in 2021-2023)", reason = "Methods"),
    list(choice = "race, BMI, smoking", decision = "RIDRETH1's five groups; BMI under 25, 25 to under 30, 30 or more; never, former, current smoking from SMQ020 and SMQ040", reason = "Methods and Table 1"),
    list(choice = "2021-2023 items", decision = "never drinker from ALQ111 and heavy drinking from ALQ130 and ALQ142, as in 2017-2018; oscillometric readings 2 and 3 (BPXO); medication from BPQ150",
         reason = "2021-2023 asks the 2017-2018 alcohol items and measures blood pressure only by the oscillometric device")
  )
)