# Visceral adiposity index (natural log) and chronic kidney disease in adults aged 60 or older,
# NHANES 2011-2018. Peng et al. (2023), Prev Med Rep, doi:10.1016/j.pmedr.2023.102306.
# Headline: OR 1.23 (95% CI 1.02-1.48) per unit of lnVAI, Model 3 (Table 2).

# Model 1 and Model 2 of Table 2, run as variants on their own (full) analytic sample.
row121_model2 <- ckd ~ ln_vai + age + sex + race + education + activity + smoking

association <- list(
  id = "row121", row = 121, doi = "10.1016/j.pmedr.2023.102306",
  cycles = c("2011-2012", "2013-2014", "2015-2016", "2017-2018"),
  weight = "WTMEC2YR", blood_file = "BIOPRO",
  family = "logistic", term = "ln_vai",
  published = list(measure = "OR", estimate = 1.23, low = 1.02, high = 1.48, n = 6085, events = NULL,
                   contrast = "per unit of ln(VAI); VAI from waist, BMI, and triglycerides and HDL cholesterol in mmol/L"),
  left_out = c("physical activity (high or low at 600 MET-minutes a week, from the GPAQ's work, transport, and leisure items; 2021-2023 asks only about leisure-time activity, so the score can't be built)"),
  build = function(cycle) {
    demo <- demographics(cycle)
    bio <- biochemistry(cycle)
    bpq <- component("BPQ", cycle)
    diq <- diabetes_questions(cycle)
    mcq <- medical_conditions(cycle)
    glu <- fasting_glucose(cycle)
    d <- merge_all(demo, body_measures(cycle)[, c("SEQN", "waist", "bmi")], hdl_cholesterol(cycle), urine_albumin_creatinine(cycle),
                   total_cholesterol(cycle), triglycerides(cycle), hba1c(cycle), smoking(cycle),
                   blood_pressure(cycle)[, c("SEQN", "sbp", "dbp")], met_minutes(cycle))
    d$creatinine <- bio$creatinine[match(d$SEQN, bio$SEQN)]
    d$tg_serum <- bio$tg_serum[match(d$SEQN, bio$SEQN)]
    d$glucose <- glu$glucose[match(d$SEQN, glu$SEQN)]
    # VAI by Amato et al. (2010), with the standard biochemistry profile's (non-fasting)
    # triglycerides, which nearly every examinee has, and triglycerides and HDL in mmol/L.
    d$vai <- vai(d$sex, d$waist, d$bmi, tg_mmol(d$tg_serum), cholesterol_mmol(d$hdl))
    d$ln_vai <- log(d$vai)
    # CKD: eGFR (CKD-EPI 2009) under 60 or urine albumin-to-creatinine ratio over 30 mg/g.
    d$egfr <- egfr(d$creatinine, d$age, d$sex, d$race)
    d$ckd <- as.integer(d$egfr < 60 | d$acr > 30)
    d$sex <- factor(d$sex)
    d$race <- factor(d$race)
    # Categorical covariates keep missing values as their own level, as the paper says.
    d$education <- with_unclear(education3(d$education), 1:3)
    d$smoking <- with_unclear(d$smoking, 1:3)
    d$activity <- with_unclear(ifelse(d$met >= 600, "high", ifelse(d$met < 600, "low", NA)), c("low", "high"))
    m <- mcq[match(d$SEQN, mcq$SEQN), ]
    parts <- cbind(m$chd, m$angina, m$heart_attack, m$stroke)
    d$cvd <- with_unclear(ifelse(rowSums(parts == 1, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(parts)) > 0, 0, NA)), 0:1)
    d$cancer <- with_unclear(m$cancer, 0:1)
    q <- diq[match(d$SEQN, diq$SEQN), ]
    dm_text <- q$told_diabetes %in% 1 | q$insulin %in% 1 | q$pills %in% 1 | (d$hba1c >= 6.5) %in% TRUE
    dm_known <- !is.na(q$told_diabetes) | !is.na(d$hba1c)
    d$diabetes <- with_unclear(ifelse(dm_text | (d$glucose >= 126) %in% TRUE, 1, ifelse(dm_known, 0, NA)), 0:1)
    d$diabetes_text <- with_unclear(ifelse(dm_text, 1, ifelse(dm_known, 0, NA)), 0:1)
    b <- bpq[match(d$SEQN, bpq$SEQN), ]
    told_bp <- yes(b$BPQ020)
    # Taking blood pressure medicine: BPQ050A (1999-2018) became BPQ150 in 2021-2023.
    bp_med <- yes(if (has(b, "BPQ050A")) b$BPQ050A else b$BPQ150)
    high_bp <- (d$sbp >= 140) | (d$dbp >= 90)
    bp_known <- !is.na(told_bp) | !is.na(high_bp)
    d$hypertension <- with_unclear(ifelse(bp_med %in% 1 | high_bp %in% TRUE, 1, ifelse(bp_known, 0, NA)), 0:1)
    d$hypertension_text <- with_unclear(ifelse(told_bp %in% 1 | bp_med %in% 1 | high_bp %in% TRUE, 1, ifelse(bp_known, 0, NA)), 0:1)
    # Cholesterol-lowering medicine: 1999-2018 ask it (BPQ100D) only of those told to take a
    # prescription (BPQ090D), itself asked only of those whose cholesterol was checked (BPQ060);
    # 2021-2023 asks everyone (BPQ101D).
    if (has(b, "BPQ101D")) {
      taking <- yes(b$BPQ101D)
    } else {
      taking <- ifelse(b$BPQ100D %in% 1, 1, ifelse(b$BPQ100D %in% 2 | b$BPQ090D %in% 2 | b$BPQ060 %in% 2 |
                                                     (b$BPQ080 %in% 2 & is.na(b$BPQ090D)), 0, NA))
    }
    d$lipid_medication <- with_unclear(taking, 0:1)
    d$tc_mmol <- cholesterol_mmol(d$tc)
    # LDL cholesterol as NHANES reports it (Friedewald, fasting subsample only), and the same
    # equation on the non-fasting serum lipids that everyone has, for a variant.
    d$ldl_mmol <- cholesterol_mmol(d$ldl)
    d$ldl_nonfasting_mmol <- cholesterol_mmol(d$tc - d$hdl - d$tg_serum / 5)
    d$in_population <- d$age >= 60 & !is.na(d$creatinine) & !is.na(d$acr) & is.finite(d$ln_vai)
    d
  },
  formula = ckd ~ ln_vai + age + sex + race + education + activity + smoking + cvd + cancer + diabetes + hypertension +
    sbp + dbp + hba1c + tc_mmol + ldl_mmol + lipid_medication,
  formula_harmonized = ckd ~ ln_vai + age + sex + race + education + smoking + cvd + cancer + diabetes + hypertension +
    sbp + dbp + hba1c + tc_mmol + ldl_mmol + lipid_medication,
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "fasting subsample weight (WTSAF2YR), as NCHS directs for LDL-C", weight = "WTSAF2YR"),
    list(label = "Model 3 without physical activity, as Table 2's note lists it",
         formula = ckd ~ ln_vai + age + sex + race + education + smoking + cvd + cancer + diabetes + hypertension +
           sbp + dbp + hba1c + tc_mmol + ldl_mmol + lipid_medication),
    list(label = "diabetes and hypertension as the Methods text defines them",
         formula = ckd ~ ln_vai + age + sex + race + education + activity + smoking + cvd + cancer + diabetes_text +
           hypertension_text + sbp + dbp + hba1c + tc_mmol + ldl_mmol + lipid_medication),
    list(label = "LDL-C from non-fasting lipids (Friedewald), keeping the full sample", sample = "own",
         formula = ckd ~ ln_vai + age + sex + race + education + activity + smoking + cvd + cancer + diabetes +
           hypertension + sbp + dbp + hba1c + tc_mmol + ldl_nonfasting_mmol + lipid_medication),
    list(label = "Model 1, full sample (published 1.35, 1.20-1.52)", formula = ckd ~ ln_vai, sample = "own"),
    list(label = "Model 1, full sample, unweighted", formula = ckd ~ ln_vai, sample = "own", weighted = FALSE),
    list(label = "Model 2, full sample (published 1.42, 1.26-1.61)", formula = row121_model2, sample = "own"),
    list(label = "Model 2, full sample, unweighted", formula = row121_model2, sample = "own", weighted = FALSE)
  ),
  # 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, evidence = "data",
         detail = "Table 2's caption gives every model n = 6,085, but Model 3 adjusts for LDL-C, which NHANES measures only in the fasting subsample, and the paper's numbers fit that subsample: Table 1's LDL-C mean (2.8 mmol/L) is the fasting value's (2.83; 2.73 by Friedewald from non-fasting lipids), and Model 3's intervals are about 1.5 times as wide as Models 1 and 2's. Our Model 3, on that subsample, has n = 2,906."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods define diabetes as HbA1c of 6.5% or more, a self-reported diagnosis, or diabetes medication, but Table 1's 26.5% is reproduced exactly only with fasting glucose of 126 mg/dL or more added (24.8% as defined)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE, evidence = "data",
         detail = "The Methods define hypertension as mean pressure of 140/90 mmHg or more, antihypertensive medication, or a self-reported diagnosis, but Table 1's 63.0% is reproduced without the self-reported diagnosis (63.1%; 66.8% as defined)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 2 prints Model 2's Q3 odds ratio as 1.04 (0.81, 1.33), the same as its Q2, with P = 0.017, which that interval rules out.")
  ),
  choices = list(
    list(choice = "weights", decision = "examination weight (WTMEC2YR pooled over four cycles; WTPH2YR in 2021-2023)",
         reason = "Table 1 is survey-weighted over the examined sample and the text doesn't say the models were unweighted; the published Model 1 and 2 intervals are as wide as design-based ones (unweighted ones are about 40% narrower), and only the weighted Model 3 shows the published Q3 above Q4"),
    list(choice = "LDL cholesterol and Model 3's sample", decision = "NHANES's LDL-C (fasting subsample) with complete-case analysis, so Model 3 is fitted on the fasting subsample (2,906) while Models 1 and 2 (variants) use the full sample",
         reason = "Table 1's LDL-C mean (2.8 mmol/L) is the fasting value's (2.83; 2.73 by Friedewald from non-fasting lipids); Model 3's published intervals are about 1.5 times as wide as Models 1 and 2's, as halving the sample would make them; Table 2's caption repeats the cohort's 6,085"),
    list(choice = "triglycerides for VAI", decision = "standard biochemistry profile (LBXSTR, non-fasting), in mmol/L",
         reason = "only 696 of 6,781 lost for VAI (fasting values exist for about half); Table 1's triglyceride mean (1.7 mmol/L) is the serum value's (1.72, fasting 1.35); VAI's sample quartiles (1.07, 1.82, 3.06) and mean (2.49) match the paper's (1.1, 1.8, 3.0; 2.5)"),
    list(choice = "renal function exclusion", decision = "both serum creatinine and urine ACR required",
         reason = "the paper excludes 'missing data on renal function, including serum creatinine and urine ACR'; this gives 6,159 against the paper's 6,085 (its intermediate 6,781 matches no definition tried)"),
    list(choice = "CKD", decision = "CKD-EPI 2009 creatinine equation with its race coefficient; ACR over 30 mg/g", reason = "as stated (Levey et al., 2009, and 'ACR >30 mg/g')"),
    list(choice = "diabetes", decision = "HbA1c 6.5% or more, fasting glucose 126 mg/dL or more, told by a doctor, or taking insulin or pills",
         reason = "Table 1's 26.5% matches this exactly; the Methods text's definition, without fasting glucose, gives 24.8% (variant)"),
    list(choice = "hypertension", decision = "taking blood pressure medicine or mean pressure 140/90 mmHg or more",
         reason = "Table 1's 63.0% matches this (63.1%); the Methods text adds a self-reported diagnosis, which gives 66.8% (variant)"),
    list(choice = "physical activity in Model 3", decision = "included", reason = "the Methods, the Results sentence reporting the headline, and the figure captions include it; only the Table 2 and 3 notes omit it (variant)"),
    list(choice = "CVD history", decision = "told of coronary heart disease, angina, heart attack, or stroke", reason = "as stated; gives 19.7% against Table 1's 22.5% (21.3% with heart failure added), unresolved"),
    list(choice = "cholesterol-lowering medicine", decision = "taking prescribed cholesterol medicine (BPQ100D; BPQ101D in 2021-2023); those never asked count as not taking", reason = "unstated item; gives Table 1's share (43.6% against 43.7%)"),
    list(choice = "blood pressure", decision = "mean of up to four auscultatory readings (oscillometric in 2021-2023), a diastolic 0 as missing", reason = "number of readings unstated; SBP and DBP means match Table 1 (132.0 and 68.8 against 131.9 and 68.2)"),
    list(choice = "covariate coding", decision = "age, SBP, DBP, HbA1c, total and LDL cholesterol continuous (complete case); sex, race (five groups), education (three), smoking (never, former, current), activity, and the yes/no conditions as factors with missing as its own level",
         reason = "the paper groups missing categorical values; Table 1's categories")
  )
)