# Triglyceride-glucose index and depression (PHQ-9 of 10 or more) in adults with type 2 diabetes,
# NHANES 2005-2020. Ren et al. (2024), Medicine, doi:10.1097/MD.0000000000039258.
# Headline: OR 1.54 (95% CI 1.21-1.95) per unit of TyG, Model III (Table 3).

# The analysis frame for a cycle. `hba1c_at` is the HbA1c (%) that counts as diabetes,
# `borderline` whether a self-reported borderline diagnosis counts, and `lipids` whether
# participants missing LDL or HDL cholesterol (study variables in Table 1) are excluded.
row311_build <- function(cycle, hba1c_at = 6.5, borderline = TRUE, lipids = TRUE) {
  demo <- demographics(cycle)
  alq <- alcohol(cycle)
  mcq <- medical_conditions(cycle)
  d <- merge_all(demo, body_measures(cycle)[, c("SEQN", "bmi")], triglycerides(cycle), hdl_cholesterol(cycle),
                 fasting_glucose(cycle), hba1c(cycle), phq9(cycle), smoking(cycle), diabetes_questions(cycle))
  d$tyg <- tyg(d$tg, d$glucose)
  d$depression <- as.integer(d$phq9 >= 10)
  told <- d$told_diabetes %in% 1 | (borderline & d$borderline %in% 1)
  d$diabetes <- told | d$insulin %in% 1 | d$pills %in% 1 | (d$glucose >= 126) %in% TRUE | (d$hba1c >= hba1c_at) %in% TRUE
  d$sex <- factor(d$sex)
  d$race <- factor(d$race)
  # Married or living with a partner; widowed, divorced, or separated; never married.
  d$marital <- factor(d$marital)
  d$bmi3 <- cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("normal", "overweight", "obese"))
  d$smoking <- factor(d$smoking)
  # Average drinks a week over the past 12 months: drinking days a year times drinks on a
  # drinking day, over 52. Never drinkers and those who didn't drink in the past year drank none.
  a <- alq[match(d$SEQN, alq$SEQN), ]
  d$weekly_drinks <- ifelse(a$alcohol3 %in% 1:2, 0, ifelse(a$alcohol3 %in% 3, a$drinking_days * a$drinks_per_day / 52, NA))
  d$alcohol <- factor(ifelse(d$weekly_drinks %in% 0, "none", ifelse(d$weekly_drinks < 8, "moderate", ifelse(d$weekly_drinks >= 8, "heavy", NA))),
                      levels = c("none", "moderate", "heavy"))
  m <- mcq[match(d$SEQN, mcq$SEQN), ]
  d$chf <- factor(m$heart_failure)
  d$cad <- factor(m$chd)
  d$in_population <- d$age >= 18 & d$diabetes & !is.na(d$tyg) & !is.na(d$phq9) &
    (!lipids | (!is.na(d$ldl) & !is.na(d$hdl)))
  d
}

row311_model2 <- depression ~ tyg + age + sex + marital + race

association <- list(
  id = "row311", row = 311, doi = "10.1097/MD.0000000000039258",
  # The paper's starting count is the sum of the 2005-2018 files and the 2017-March 2020 files,
  # which repeat the 2017-2018 participants; the 2017-March 2020 files alone cover those years.
  cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2020"),
  weight = "WTSAF2YR",
  family = "logistic", term = "tyg",
  published = list(measure = "OR", estimate = 1.54, low = 1.21, high = 1.95, n = 3225, events = 364,
                   contrast = "per unit of the TyG index, ln(fasting triglycerides x fasting glucose / 2), both in mg/dL"),
  left_out = character(),
  build = function(cycle) row311_build(cycle),
  formula = depression ~ tyg + age + sex + marital + race + bmi3 + smoking + alcohol + chf + cad,
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "examination weight (WTMEC2YR), keeping those without a fasting weight", weight = "WTMEC2YR"),
    list(label = "2017-2018 and 2017-March 2020 files stacked, as the paper's counts show",
         cycles = c("2005-2006", "2007-2008", "2009-2010", "2011-2012", "2013-2014", "2015-2016", "2017-2018", "2017-2020")),
    list(label = "borderline diabetes not counted", build = function(cycle) row311_build(cycle, borderline = FALSE)),
    list(label = "HbA1c 5.29% or more counts as diabetes, as printed", build = function(cycle) row311_build(cycle, hba1c_at = 5.29)),
    list(label = "missing LDL or HDL cholesterol not excluded", build = function(cycle) row311_build(cycle, lipids = FALSE)),
    list(label = "Model I (published 1.60, 1.29-1.99)", formula = depression ~ tyg),
    list(label = "Model II (published 1.64, 1.30-2.07)", formula = row311_model2)
  ),
  # 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 = "model", affects_headline = TRUE, followed = TRUE,
         detail = "Section 2.3 says its confounders, among them blood glucose, triglycerides, HDL, and LDL, were all included in the analyses, but Model III's estimate is reproduced without them, with the Table 3 footnote's covariates (1.53 against 1.54), and not with glucose and triglycerides added (3.05, 1.17-7.91)."),
    list(kind = "reporting", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods count an HbA1c of 5.29% or more as diabetes, but that cutoff gives 11,598 participants against the published 3,225; with 6.5%, the sample gives the Methods' TyG quartile cutpoints (8.60, 9.03, 9.47 against 8.61, 9.03, 9.47) and Model III's estimate (1.53 against 1.54)."),
    list(kind = "sample", affects_headline = FALSE, followed = NA,
         detail = "Figure 1 starts from 85,750 participants, the 2005-2018 files plus the 2017-March 2020 files, which repeat the 2017-2018 participants, but the final sample and estimate fit the files without the repeats (weighted fits: 3,111 participants and 1.53 without them, 3,626 and 1.61 with them; published 3,225 and 1.54)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Tables 1 and 2 label triglycerides, HDL, and LDL in mmol/L, but the values are in mg/dL (triglycerides 142.66, against 142.6 mg/dL here)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Tables 1 and 2 swap the smoking labels Never and Now: their 533 'never' smokers are the current smokers (556 here) and their 1,664 'now' smokers the never smokers (1,671)."),
    list(kind = "reporting", affects_headline = FALSE, followed = NA,
         detail = "Table 1's alcohol counts (1,204, 1,162, 368) sum to 2,734 and their percentages to 84.9%, while Table 2 gives 1,204, 1,530, and 491 for the same 3,225 participants.")
  ),
  choices = list(
    list(choice = "cycles", decision = "2005-2016 and the 2017-March 2020 files, without the 2017-2018 files",
         reason = "the paper's starting count (85,750) is the 2005-2018 files plus the 2017-March 2020 files, which repeat the 2017-2018 participants; stacked (variant) the estimate rises to 1.61 on 3,624, farther from the published, so the paper's final sample probably lost the duplicates"),
    list(choice = "weights", decision = "fasting subsample weight (WTSAF2YR; WTSAFPRP in 2017-March 2020), pooled by years",
         reason = "the paper used 'the recommended weights', which NCHS gives as the fasting subsample weight for fasting triglycerides and glucose; the examination weight (variant) gives nearly the same estimate (1.53)"),
    list(choice = "HbA1c cutoff for diabetes", decision = "6.5%",
         reason = "printed as 5.29%, which gives an analytic sample of 11,598 against the paper's 3,225 (variant); 6.5% is the standard cutoff its cited definitions use"),
    list(choice = "borderline diabetes", decision = "a self-reported borderline diagnosis (DIQ010 = 3) counts as self-reported diabetes",
         reason = "the TyG quartile cutpoints in the Methods (8.61, 9.03, 9.47) and Table 1's mean TyG (9.07) and glucose (146.3) match only with it (8.60, 9.03, 9.47 over everyone with TyG values; 9.06; 147.3; without it 8.65, 9.06, 9.52; 9.10; 151.3), and Models I-III then match (1.60, 1.61, 1.53 against 1.60, 1.64, 1.54)"),
    list(choice = "LDL and HDL cholesterol", decision = "participants missing either are excluded",
         reason = "the paper excludes those missing 'other study variables' and lists LDL and HDL among them; Table 1's triglyceride mean (142.7 mg/dL) is that of participants with LDL, which NHANES computes only for triglycerides of 400 or less (160 without the step, whose estimate is 1.34)"),
    list(choice = "type 2 diabetes", decision = "every diabetic adult; type 1 not separated",
         reason = "the paper doesn't say how it separated type 1; leaving out those diagnosed before 30 who take only insulin (70) gives 1.49"),
    list(choice = "depression", decision = "PHQ-9 total of 10 or more, missing if any item is missing, refused, or don't know", reason = "as stated; item nonresponse unstated"),
    list(choice = "alcohol", decision = "average drinks a week over the past 12 months (drinking days a year times drinks a drinking day over 52): none 0, moderate under 8, heavy 8 or more; never drinkers and past-year abstainers drink none",
         reason = "the paper's categories (0, 1-8, 8 or more), which overlap at 8; Table 1's counts (1,204, 1,162, 368: 2,734) and Table 2's (1,204, 1,530, 491: 3,225) disagree, unresolved"),
    list(choice = "smoking", decision = "never, former, current from SMQ020 and SMQ040",
         reason = "convention; Table 1's 'Never' (533) and 'Now' (1,664) labels look swapped (here 1,671 never and 556 current)"),
    list(choice = "covariate coding", decision = "age continuous; marital status in three groups (married or with a partner; widowed, divorced, or separated; never married); race in five (RIDRETH1); BMI under 25, 25 to under 30, 30 or more; heart failure (MCQ160B) and coronary heart disease (MCQ160C) as told; complete case",
         reason = "Table 1's categories (its 'Divorced' holds widowed and separated, as its counts show); the paper excludes missing covariates"),
    list(choice = "participants with TyG values but no fasting weight", decision = "in the population, out of the weighted fit (206; they fasted too briefly)",
         reason = "the paper's counts include everyone with TyG values; a fasting-weighted analysis gives them no weight")
  )
)