# Serum uric acid and creatine phosphokinase (CPK) in adults aged 30 and older, NHANES 2015-2018.
# Chen et al. (2023), BMC Cardiovasc Disord, doi:10.1186/s12872-023-03333-5.
# Headline: beta 0.12 (95% CI 0.09-0.15) IU/L of CPK per 1 umol/L of uric acid, Model 3 (Table 2).
#
# The flow chart is reproduced exactly: 9,486 aged 30 or older, 1,012 without uric acid or CPK
# (1,004 and 8), 43 with CPK over 1000 U/L, 8,431 analyzed. Model 3 keeps all 8,431 although LDL
# and fasting triglycerides exist only for the morning fasting subsample (3,864 and 3,907 of them):
# Table 1 shows drinking with a "Missing data" group, and its weighted means and SDs are those of
# covariates whose missing values were filled with the sample mean: LDL 2.9 +- 0.6 mmol/L (2.94 +-
# 0.63 filled, 2.96 +- 0.92 not), triglycerides 1.3 +- 0.7 (1.33 +- 0.66, against 0.97), income-
# poverty ratio 3.1 +- 1.6 (3.11 +- 1.57, against 3.17 +- 1.64), waist 102.4 +- 16.0 (102.3 +- 16.00,
# against 16.37). So the chosen version fills missing continuous covariates with the sample mean and
# keeps missing categorical ones as their own level. Table 1's shares also show which codings Model 3
# ("adjusted for covariates, as shown in Table 1") used: BMI cut at 24 and 28 (17.9, 25.5, 55.6%
# against the printed 17.9, 25.4, 55.7; the labels say 25 and 30, which give 23.4, 32.7, 42.8),
# hypertension as self-report alone (BPQ020: 38.1% against 38.0; the text's definition with measured
# pressure of 130/80 or more and treatment gives 57.9%), and diabetes as self-report alone (DIQ010:
# 13.7% yes, 83.6% no against 13.6 and 83.7; adding treatment gives 14.8%).
#
# Weighted fits put Model 3 inside the published interval (0.110) and come closest to Table 2's
# subgroup estimates (Black adults' Model 1 exactly), but Models 1 and 2 come out a little below the
# published (0.282 against 0.30, 0.133 against 0.14; unweighted 0.293 and 0.154), and Table 1's CPK
# (139.0 +- 120.7 IU/L) matches neither our weighted (136.6 +- 107.4) nor unweighted (146.3 +- 117.0)
# values; no population, weight, or coding tried explains the gap.

# Drinking frequency in the paper's three groups. The text: none or rarely, sometimes (1 to 3 times
# a month), often (once a week or more); past-year frequency from ALQ120Q/U (2015-2016) or ALQ121
# (2017 on), as days a year (the library's alcohol()). Lifetime abstainers, whom the frequency
# question skips, drink none. table1 = TRUE codes it as Table 1's shares show it was coded (none or
# rarely 25.3%, sometimes 18.6, often 29.3, missing 26.7; this gives 25.3, 18.5, 29.3, 27.0): by the
# answer's unit in 2015-2016 (week often, month sometimes, year or never in the past year none or
# rarely), ALQ121's once a week or more often and one to three times a month sometimes, never in the
# last year none or rarely, and its less-than-monthly answers (7 to 11, 3 to 6, 1 to 2 times in the
# last year) left missing, as are lifetime abstainers in every cycle.
row194_drinking <- function(cycle, table1 = FALSE) {
  if (!table1) {
    a <- alcohol(cycle)
    days <- a$drinking_days
    group <- ifelse((days >= 52) %in% TRUE, "often", ifelse((days >= 12) %in% TRUE, "sometimes",
                    ifelse((days < 12) %in% TRUE | a$alcohol3 %in% 1, "none or rarely", NA)))
    return(data.frame(SEQN = a$SEQN, drinking = group))
  }
  d <- component("ALQ", cycle)
  if (has(d, "ALQ121")) {
    group <- ifelse(d$ALQ121 %in% 1:5, "often", ifelse(d$ALQ121 %in% 6:7, "sometimes", ifelse(d$ALQ121 %in% 0, "none or rarely", NA)))
  } else {
    q <- ifelse(d$ALQ120Q %in% 0:365, d$ALQ120Q, NA)
    group <- ifelse(q %in% 0, "none or rarely", ifelse(!is.na(q) & d$ALQ120U %in% 1, "often",
                    ifelse(!is.na(q) & d$ALQ120U %in% 2, "sometimes", ifelse(!is.na(q) & d$ALQ120U %in% 3, "none or rarely", NA))))
  }
  data.frame(SEQN = d$SEQN, drinking = group)
}

# Moderate recreational activity, yes or no: "In a typical week do you do any moderate-intensity
# sports, fitness, or recreational activities ... for at least 10 minutes continuously" (PAQ665).
# 2021-2023 asks how often instead (PAD790Q, 0 = never); any is yes.
row194_activity <- function(cycle) {
  p <- component("PAQ", cycle)
  if (has(p, "PAQ665")) return(data.frame(SEQN = p$SEQN, active = ifelse(p$PAQ665 %in% 1, "yes", ifelse(p$PAQ665 %in% 2, "no", NA))))
  q <- p$PAD790Q
  data.frame(SEQN = p$SEQN, active = ifelse(q %in% 0, "no", ifelse(!is.na(q) & !q %in% c(7777, 9999), "yes", NA)))
}

row194_answer <- function(x) ifelse(x %in% 1, "yes", ifelse(x %in% 2, "no", NA))

row194_build <- function(cycle) {
  demo <- demographics(cycle)
  bio <- component("BIOPRO", cycle, c("LBDSUASI", "LBXSCK", "LBXSATSI", "LBXSASSI", "LBDSTPSI"))
  d <- merge_all(demo, bio, component("HDL", cycle, "LBDHDDSI"), component("TCHOL", cycle, "LBDTCSI"),
                 component("TRIGLY", cycle, c("LBDTRSI", "LBDLDLSI")), body_measures(cycle)[, c("SEQN", "bmi", "waist")],
                 row194_activity(cycle), component("SMQ", cycle, "SMQ020"), row194_drinking(cycle),
                 setNames(row194_drinking(cycle, table1 = TRUE), c("SEQN", "drinking_table1")),
                 component("BPQ", cycle, "BPQ020"), component("DIQ", cycle, "DIQ010"), component("KIQ_U", cycle, "KIQ022"),
                 setNames(hypertension_status(cycle, sbp_at = 130, dbp_at = 80), c("SEQN", "hypertension_text")),
                 setNames(diabetes_status(cycle, parts = c("told", "medication")), c("SEQN", "diabetes_text")))
  d$uric_acid <- d$LBDSUASI
  d$cpk <- d$LBXSCK
  d$alt <- d$LBXSATSI
  d$ast <- d$LBXSASSI
  d$total_protein <- d$LBDSTPSI
  d$hdl <- d$LBDHDDSI
  d$tc <- d$LBDTCSI
  d$tg <- d$LBDTRSI
  d$ldl <- d$LBDLDLSI
  d$race4 <- ifelse(d$race %in% 3, "white", ifelse(d$race %in% 4, "black", ifelse(d$race %in% 1, "mexican american", ifelse(d$race %in% c(2, 5), "other", NA))))
  # Less than high school, high school or GED, college graduate or above (some college included:
  # Table 1's 13.1, 22.8, 64.1% are DMDEDUC2 1-2, 3, 4-5, which give 13.3, 22.7, 63.9).
  d$education <- education3(d$education)
  d$bmi_group <- as.character(cut(d$bmi, c(-Inf, 24, 28, Inf), right = FALSE, labels = c("normal", "overweight", "obese")))
  d$bmi_group_text <- as.character(cut(d$bmi, c(-Inf, 25, 30, Inf), right = FALSE, labels = c("normal", "overweight", "obese")))
  d$smoker <- row194_answer(d$SMQ020)
  d$hypertension <- row194_answer(d$BPQ020)
  d$diabetes <- row194_answer(d$DIQ010)
  d$kidney <- row194_answer(d$KIQ022)
  d$hypertension_text <- ifelse(d$hypertension_text %in% 1, "yes", ifelse(d$hypertension_text %in% 0, "no", NA))
  d$diabetes_text <- ifelse(d$diabetes_text %in% 1, "yes", ifelse(d$diabetes_text %in% 0, "no", NA))
  d$in_population <- d$age >= 30 & !is.na(d$uric_acid) & !is.na(d$cpk) & d$cpk <= 1000
  d
}

ROW194_CONTINUOUS <- c("pir", "waist", "alt", "ast", "hdl", "tc", "tg", "ldl", "total_protein")
ROW194_CATEGORICAL <- list(sex = c(1, 2), race4 = c("white", "black", "mexican american", "other"), education = 1:3,
                           bmi_group = c("normal", "overweight", "obese"), bmi_group_text = c("normal", "overweight", "obese"),
                           active = c("no", "yes"), smoker = c("no", "yes"),
                           drinking = c("none or rarely", "sometimes", "often"), drinking_table1 = c("none or rarely", "sometimes", "often"),
                           hypertension = c("no", "yes"), diabetes = c("no", "yes"), kidney = c("no", "yes"),
                           hypertension_text = c("no", "yes"), diabetes_text = c("no", "yes"))

# Missing data as the paper's sample shows it was handled: a continuous covariate takes the mean of
# the analytic sample (of the data being fitted, so 2021-2023 takes its own), a categorical one keeps
# missing as its own level. The unfilled values stay in <name>_cc for the complete-case variant. A
# level nobody in the sample has is dropped, so no model gets an empty column.
row194_derive <- function(data, constants) {
  population <- data$in_population %in% TRUE
  for (v in ROW194_CONTINUOUS) {
    data[[paste0(v, "_cc")]] <- data[[v]]
    data[[v]] <- ifelse(is.na(data[[v]]), mean(data[[v]][population], na.rm = TRUE), data[[v]])
  }
  present <- function(f) factor(f, levels = levels(f)[levels(f) %in% as.character(f[population])])
  for (v in names(ROW194_CATEGORICAL)) {
    levels <- ROW194_CATEGORICAL[[v]]
    data[[paste0(v, "_cc")]] <- present(factor(data[[v]], levels = levels))
    data[[v]] <- present(with_unclear(data[[v]], levels))
  }
  data
}

ROW194_MODEL3 <- cpk ~ uric_acid + age + sex + race4 + education + pir + bmi_group + waist + active + smoker + drinking +
  alt + ast + hdl + tc + tg + ldl + total_protein + hypertension + diabetes + kidney
ROW194_MODEL2 <- cpk ~ uric_acid + age + sex + race4

# The fit among non-Hispanic Black adults (Table 2's subgroup, adjusted for nothing).
row194_black <- function(data, constants) {
  data <- row194_derive(data, constants)
  data$in_population <- data$in_population & data$race4 %in% "black"
  data
}

association <- list(
  id = "row194", row = 194, doi = "10.1186/s12872-023-03333-5",
  cycles = c("2015-2016", "2017-2018"),
  weight = "WTMEC2YR", blood_file = "BIOPRO",
  family = "linear", term = "uric_acid",
  published = list(measure = "beta", estimate = 0.12, low = 0.09, high = 0.15, n = 8431, events = NULL,
                   contrast = "IU/L of serum creatine phosphokinase per 1 umol/L of serum uric acid"),
  left_out = character(),
  build = row194_build,
  derive = row194_derive,
  # Drinking as the paper's computation coded it (Table 1's groups, with its "Missing data" group).
  formula = update(ROW194_MODEL3, . ~ . - drinking + drinking_table1),
  variants = list(
    list(label = "unweighted", weighted = FALSE),
    list(label = "Model 1, weighted (published 0.30, 0.27-0.33)", formula = cpk ~ uric_acid),
    list(label = "Model 1, unweighted", formula = cpk ~ uric_acid, weighted = FALSE),
    list(label = "Model 2, weighted (published 0.14, 0.11-0.17)", formula = ROW194_MODEL2),
    list(label = "Model 2, unweighted", formula = ROW194_MODEL2, weighted = FALSE),
    list(label = "Black adults, Model 1, weighted (published 0.42, 0.34-0.51)", formula = cpk ~ uric_acid, derive = row194_black),
    list(label = "Black adults, Model 1, unweighted", formula = cpk ~ uric_acid, derive = row194_black, weighted = FALSE),
    list(label = "complete case (LDL and TG only in the fasting subsample)", sample = "own",
         formula = cpk ~ uric_acid + age + sex_cc + race4_cc + education_cc + pir_cc + bmi_group_cc + waist_cc + active_cc + smoker_cc +
           drinking_cc + alt_cc + ast_cc + hdl_cc + tc_cc + tg_cc + ldl_cc + total_protein_cc + hypertension_cc + diabetes_cc + kidney_cc),
    list(label = "text's definitions: BMI at 25 and 30, BP 130/80 or treated, diabetes treated",
         formula = update(ROW194_MODEL3, . ~ . - bmi_group - hypertension - diabetes + bmi_group_text + hypertension_text + diabetes_text)),
    list(label = "drinking as the text defines it", formula = ROW194_MODEL3)
  ),
  # 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,
         detail = "The Methods and Table 1's labels cut BMI at 25 and 30 kg/m2, but Table 1's shares (17.9, 25.4, 55.7%) are those of cuts at 24 and 28 (17.9, 25.5, 55.6%; 25 and 30 give 23.4, 32.7, 42.8%), and Model 3 adjusts for the covariates as shown in Table 1."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods define hypertension as mean pressure of 130/80 mmHg or more, antihypertensive drugs, or a diagnosis, but Table 1's 38.0% is self-reported diagnosis alone (38.1%; the stated definition gives 57.9%)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods define diabetes as a diagnosis or use of insulin or diabetes medication, but Table 1's shares (13.6% yes, 83.7% no) are self-reported diagnosis alone (13.7% and 83.6%; adding treatment gives 14.8% yes)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods group drinking as none or rarely, sometimes (1 to 3 times a month), and often (once a week or more), but Table 1's shares (25.3, 18.6, 29.3%, and 26.7% missing) are those of 2015-2016's answers grouped by the unit they were given in, with 2017-2018's less-than-monthly answers and all lifetime abstainers left missing (25.3, 18.5, 29.3, 27.0%)."),
    list(kind = "coding", affects_headline = TRUE, followed = TRUE,
         detail = "The Methods and Table 1 name the top education group college graduate or above, but Table 1's shares (13.1, 22.8, 64.1%) are those with some college counted in it (13.3, 22.7, 63.9%).")
  ),
  choices = list(
    list(choice = "weights", decision = "examination weight (WTMEC2YR over 2, with strata and PSUs; WTPH2YR from BIOPRO_L in 2021-2023)",
         reason = "the paper says it used NHANES sample weights without naming one; weighted fits come closest to Table 2's subgroup estimates (Black adults' Model 1: 0.42, 0.34-0.51, exactly, against 0.36 unweighted; Mexican Americans 0.36 against the published 0.41, 0.26 unweighted; Whites 0.21 against 0.22, 0.18 unweighted)"),
    list(choice = "missing covariates", decision = "continuous covariates filled with the analytic sample's mean, categorical ones with missing as its own level; all 8,431 kept",
         reason = "unstated; the paper analyzed 8,431 although LDL and fasting triglycerides exist for only 3,864 and 3,907 of them, Table 1 shows drinking with a 'Missing data' group, and Table 1's means and SDs are the mean-filled values' (LDL 2.9 +- 0.6 mmol/L, 0.92 unfilled; triglycerides 1.3 +- 0.7, 0.97 unfilled; income-poverty ratio 3.1, 3.17 unfilled; waist SD 16.0, 16.37 unfilled). Complete-case analysis keeps 2,966 (variant)"),
    list(choice = "triglycerides and total cholesterol", decision = "fasting triglycerides and LDL (TRIGLY, LBDTRSI and LBDLDLSI) and total cholesterol from TCHOL (LBDTCSI)",
         reason = "unstated; Table 1's means fit these (TC 5.0 against 5.03 from TCHOL and 5.08 from the biochemistry profile; TG 1.3 against 1.33 fasting and 1.77 non-fasting)"),
    list(choice = "BMI", decision = "three groups cut at 24 and 28 kg/m2",
         reason = "the text and Table 1's labels cut at 25 and 30, but Table 1's shares are the cuts at 24 and 28 (17.9, 25.5, 55.6% against 17.9, 25.4, 55.7), so Model 3 used them; 25 and 30 is part of a variant"),
    list(choice = "hypertension and diabetes", decision = "self-reported diagnosis (BPQ020; DIQ010 yes or no, borderline with missing as its own level)",
         reason = "the text adds measured pressure of 130/80 or more, treatment, and diabetes treatment, but Table 1's shares are self-report alone (38.1% hypertension against 38.0; 13.7% and 83.6% for diabetes against 13.6 and 83.7); the text's definitions are a variant"),
    list(choice = "drinking", decision = "the text's groups (often once a week or more, sometimes 1 to 3 times a month, none or rarely less than monthly or never) as the paper's computation coded them: ALQ121's less-than-monthly answers (2017-2018) and lifetime abstainers missing, 2015-2016 answers grouped by their unit, missing as its own level",
         reason = "Table 1's shares show this coding; the published estimate comes from it, and the text's coding (lifetime abstainers none) is a variant"),
    list(choice = "race", decision = "non-Hispanic White, non-Hispanic Black, Mexican American, other (other Hispanic included)",
         reason = "the text's four groups; Table 1's 16.0% other is RIDRETH1 2 and 5 (16.1%)"),
    list(choice = "education", decision = "less than high school, high school or GED, college graduate or above with some college included",
         reason = "Table 1's shares (13.1, 22.8, 64.1%) are DMDEDUC2 1-2, 3, 4-5 (13.3, 22.7, 63.9%)"),
    list(choice = "covariate coding", decision = "age, income-poverty ratio, waist, ALT, AST, HDL, total cholesterol, triglycerides, LDL (mmol/L) and total protein (g/L) continuous; activity, smoking, kidney disease yes or no",
         reason = "Table 1 gives the continuous ones as means and the others as yes or no"),
    list(choice = "moderate recreational activity in 2021-2023", decision = "any moderate leisure-time activity (PAD790Q above 0)",
         reason = "2021-2023 does not ask PAQ665 (in a typical week, at least 10 minutes); it asks how often"),
    list(choice = "CPK and uric acid across analyzers", decision = "values as released in every cycle, no crosswalk",
         reason = "the paper pooled the DxC800 (2015-2016) and Cobas 6000 (2017-2018) values as released. 2021-2023 measured on a Cobas 8000; NCHS's bridging study recommends no adjustment for uric acid but gives Cobas 8000 = -0.5307 + 1.017 x Cobas 6000 for CPK (and Cobas 8000 = 1.477 + 0.9663 x Cobas 6000 for ALT, a covariate). These are recorded as measurement changes, not applied, as the plan reports measurement changes without correcting them")
  )
)