# Variables that many analyses share, built the same way in every cycle. Each function returns
# one row per SEQN of the file it reads; where NHANES renamed or changed a variable between
# cycles, the function says how the cycles are made to match.

# NHANES codes 7, 9, 77, 99, 777, 999, 7777, 9999 as refused or don't know.
refused <- function(x, codes = c(7, 9)) ifelse(x %in% codes, NA, x)

yes <- function(x) ifelse(x %in% 1, 1, ifelse(x %in% 2, 0, NA))

has <- function(data, variable) variable %in% names(data)

demographics <- function(cycle) {
  d <- component("DEMO", cycle)
  marital <- if (has(d, "DMDMARTL")) {
    # 1 married or living with a partner, 2 widowed, divorced, or separated, 3 never married:
    # the grouping DMDMARTZ uses from 2017-March 2020 on.
    ifelse(d$DMDMARTL %in% c(1, 6), 1, ifelse(d$DMDMARTL %in% 2:4, 2, ifelse(d$DMDMARTL %in% 5, 3, NA)))
  } else {
    ifelse(d$DMDMARTZ %in% 1:3, d$DMDMARTZ, NA)
  }
  data.frame(
    SEQN = d$SEQN, SDMVSTRA = d$SDMVSTRA, SDMVPSU = d$SDMVPSU,
    age = d$RIDAGEYR,
    sex = d$RIAGENDR,  # 1 male, 2 female
    race = d$RIDRETH1,  # 1 Mexican American, 2 other Hispanic, 3 NH White, 4 NH Black, 5 other
    race3 = if (has(d, "RIDRETH3")) d$RIDRETH3 else NA,  # adds 6 NH Asian (2011 on)
    education = ifelse(d$DMDEDUC2 %in% 1:5, d$DMDEDUC2, NA),  # adults 20+
    pir = d$INDFMPIR,
    marital = marital,
    pregnant = if (has(d, "RIDEXPRG")) ifelse(d$RIDEXPRG %in% 1, 1, 0) else 0,
    exam = ifelse(d$RIDSTATR %in% 2, 1, 0)
  )
}

# Education in three levels: less than high school, high school or GED, more than high school.
education3 <- function(education) ifelse(education %in% 1:2, 1, ifelse(education %in% 3, 2, ifelse(education %in% 4:5, 3, NA)))

smoking <- function(cycle) {
  d <- component("SMQ", cycle)
  ever <- yes(d$SMQ020)
  now <- ifelse(d$SMQ040 %in% 1:2, 1, ifelse(d$SMQ040 %in% 3, 0, NA))
  # 1 never (fewer than 100 cigarettes in life), 2 former, 3 current
  status <- ifelse(ever %in% 0, 1, ifelse(ever %in% 1 & now %in% 0, 2, ifelse(ever %in% 1 & now %in% 1, 3, NA)))
  data.frame(SEQN = d$SEQN, smoking = status)
}

body_measures <- function(cycle) {
  d <- component("BMX", cycle)
  data.frame(SEQN = d$SEQN, weight = d$BMXWT, height = d$BMXHT, bmi = d$BMXBMI, waist = d$BMXWAIST,
             hip = if (has(d, "BMXHIP")) d$BMXHIP else NA,
             bmi_child = if (has(d, "BMDBMIC")) d$BMDBMIC else NA)
}

# PHQ-9 total score (0-27), from the nine items scored 0-3; missing if any item is missing,
# refused, or answered "don't know".
phq9 <- function(cycle) {
  d <- component("DPQ", cycle)
  items <- sprintf("DPQ%03d", seq(10, 90, 10))
  scores <- sapply(items, function(item) ifelse(d[[item]] %in% 0:3, d[[item]], NA))
  data.frame(SEQN = d$SEQN, phq9 = rowSums(scores), phq9_item9 = scores[, "DPQ090"])
}

# Mean systolic and diastolic blood pressure over a cycle's readings. 1999-2018 measured blood
# pressure by auscultation (BPX, up to four readings); 2017-March 2020 and 2021-2023 by an
# oscillometric device (BPXO, three readings). method = "auscultatory" or "oscillometric" picks
# the protocol where a cycle has both (2017-2018). A diastolic reading of 0 is treated as missing.
blood_pressure <- function(cycle, method = "auto", readings = 1:4, zero_diastolic = "missing") {
  oscillometric <- if (method == "auto") cycle %in% c("2017-2020", REPLICATION_CYCLE) else method == "oscillometric"
  if (oscillometric) {
    d <- component("BPXO", cycle)
    sbp <- d[, intersect(sprintf("BPXOSY%d", readings), names(d)), drop = FALSE]
    dbp <- d[, intersect(sprintf("BPXODI%d", readings), names(d)), drop = FALSE]
  } else {
    d <- component("BPX", cycle)
    sbp <- d[, intersect(sprintf("BPXSY%d", readings), names(d)), drop = FALSE]
    dbp <- d[, intersect(sprintf("BPXDI%d", readings), names(d)), drop = FALSE]
  }
  if (zero_diastolic == "missing") dbp[dbp == 0] <- NA
  data.frame(SEQN = d$SEQN, sbp = rowMeans(sbp, na.rm = TRUE), dbp = rowMeans(dbp, na.rm = TRUE),
             sbp_n = rowSums(!is.na(sbp)), dbp_n = rowSums(!is.na(dbp)))
}

# Blood pressure questionnaire: told high blood pressure, taking medication for it, told high
# cholesterol, taking cholesterol medication. BPQ050A (1999-2018) became BPQ150 in 2021-2023,
# and BPQ100D became BPQ101D; both are asked of those told of the condition.
bp_questions <- function(cycle) {
  d <- component("BPQ", cycle)
  pick <- function(...) { for (v in c(...)) if (has(d, v)) return(d[[v]]); rep(NA, nrow(d)) }
  data.frame(SEQN = d$SEQN,
             told_hypertension = yes(d$BPQ020),
             bp_medication = yes(pick("BPQ050A", "BPQ150")),
             told_cholesterol = yes(pick("BPQ080")),
             cholesterol_medication = yes(pick("BPQ100D", "BPQ101D")))
}

diabetes_questions <- function(cycle) {
  d <- component("DIQ", cycle)
  data.frame(SEQN = d$SEQN,
             told_diabetes = ifelse(d$DIQ010 %in% 1, 1, ifelse(d$DIQ010 %in% 2:3, 0, NA)),
             borderline = ifelse(d$DIQ010 %in% 3, 1, 0),
             told_prediabetes = if (has(d, "DIQ160")) yes(d$DIQ160) else NA,
             insulin = yes(d$DIQ050),
             # Taking diabetic pills: DIQ070, named DID070 in 2005-2006.
             pills = if (has(d, "DIQ070")) yes(d$DIQ070) else if (has(d, "DID070")) yes(d$DID070) else NA)
}

hba1c <- function(cycle) {
  d <- component("GHB", cycle)
  data.frame(SEQN = d$SEQN, hba1c = d$LBXGH)
}

# Fasting subsample values. `fasted` is 1 for those with a positive fasting-subsample weight, whom
# NCHS counts as fasted (8 to 24 hours); values of others are in the files too.
fasting_weight <- function(d) {
  w <- if (has(d, "WTSAF2YR")) d$WTSAF2YR else if (has(d, "WTSAFPRP")) d$WTSAFPRP else if (has(d, "WTSAF4YR")) d$WTSAF4YR else NA
  as.integer(!is.na(w) & w > 0)
}

fasting_glucose <- function(cycle) {
  d <- component("GLU", cycle)
  data.frame(SEQN = d$SEQN, glucose = d$LBXGLU, fasted = fasting_weight(d))
}

# Triglycerides and LDL from the fasting subsample. 2021-2023 renamed LBXTR to LBXTLG.
triglycerides <- function(cycle) {
  d <- component("TRIGLY", cycle)
  tg <- if (has(d, "LBXTR")) d$LBXTR else d$LBXTLG
  data.frame(SEQN = d$SEQN, tg = tg, ldl = if (has(d, "LBDLDL")) d$LBDLDL else NA, fasted = fasting_weight(d))
}

total_cholesterol <- function(cycle) {
  d <- component("TCHOL", cycle)
  data.frame(SEQN = d$SEQN, tc = d$LBXTC)
}

# HDL cholesterol: LBDHDL in 1999-2002, LBXHDD in 2003-2004, LBDHDD from 2005 on.
hdl_cholesterol <- function(cycle) {
  d <- component("HDL", cycle)
  hdl <- if (has(d, "LBDHDD")) d$LBDHDD else if (has(d, "LBXHDD")) d$LBXHDD else d$LBDHDL
  data.frame(SEQN = d$SEQN, hdl = hdl)
}

blood_count <- function(cycle) {
  d <- component("CBC", cycle)
  data.frame(SEQN = d$SEQN, wbc = d$LBXWBCSI, lymphocytes = d$LBDLYMNO, monocytes = d$LBDMONO,
             neutrophils = d$LBDNENO, eosinophils = d$LBDEONO, basophils = d$LBDBANO,
             platelets = d$LBXPLTSI, rdw = d$LBXRDW, hemoglobin = d$LBXHGB, mcv = d$LBXMCVSI)
}

biochemistry <- function(cycle) {
  d <- component("BIOPRO", cycle)
  get <- function(v) if (has(d, v)) d[[v]] else NA
  data.frame(SEQN = d$SEQN, uric_acid = get("LBXSUA"), creatinine = get("LBXSCR"), albumin = get("LBXSAL"),
             alt = get("LBXSATSI"), ast = get("LBXSASSI"), ggt = get("LBXSGTSI"), alp = get("LBXSAPSI"),
             cpk = get("LBXSCK"), iron = get("LBXSIR"), bun = get("LBXSBU"), glucose_serum = get("LBXSGL"),
             tg_serum = get("LBXSTR"), cholesterol_serum = get("LBXSCH"), total_protein = get("LBXSTP"),
             bilirubin = get("LBXSTB"), calcium = get("LBXSCA"), phosphorus = get("LBXSPH"))
}

urine_albumin_creatinine <- function(cycle) {
  d <- component("ALB_CR", cycle)
  acr <- if (has(d, "URDACT")) d$URDACT else 100 * d$URXUMA / d$URXUCR
  data.frame(SEQN = d$SEQN, acr = acr)
}

# Serum creatinine (mg/dL) on the standardized scale the eGFR equations assume. NCHS's
# documentation directs correcting two cycles' values (Selvin et al., 2007): 1999-2000 by
# 1.013 x + 0.147 (LAB18) and 2005-2006 by -0.016 + 0.978 x (BIOPRO_D); the others need none.
standard_creatinine <- function(creatinine, cycle) {
  if (cycle == "1999-2000") return(1.013 * creatinine + 0.147)
  if (cycle == "2005-2006") return(-0.016 + 0.978 * creatinine)
  creatinine
}

# eGFR (mL/min/1.73 m2) by the CKD-EPI creatinine equation of 2009 (with its race term) or 2021,
# from standardized creatinine (standard_creatinine()).
egfr <- function(creatinine, age, sex, race, equation = "2009") {
  female <- sex == 2
  if (equation == "2009") {
    kappa <- ifelse(female, 0.7, 0.9)
    alpha <- ifelse(female, -0.329, -0.411)
    141 * pmin(creatinine / kappa, 1)^alpha * pmax(creatinine / kappa, 1)^-1.209 * 0.993^age *
      ifelse(female, 1.018, 1) * ifelse(race == 4, 1.159, 1)
  } else {
    kappa <- ifelse(female, 0.7, 0.9)
    alpha <- ifelse(female, -0.241, -0.302)
    142 * pmin(creatinine / kappa, 1)^alpha * pmax(creatinine / kappa, 1)^-1.200 * 0.9938^age * ifelse(female, 1.012, 1)
  }
}

medical_conditions <- function(cycle) {
  d <- component("MCQ", cycle)
  get <- function(v) if (has(d, v)) yes(d[[v]]) else NA
  copd <- if (has(d, "MCQ160P")) yes(d$MCQ160P) else {
    # 1999-2018 asked emphysema, chronic bronchitis, and (from 2013) COPD separately; 2017-March
    # 2020 and 2021-2023 ask them as one question (MCQ160P).
    parts <- cbind(get("MCQ160G"), get("MCQ160K"), get("MCQ160O"))
    ifelse(rowSums(parts == 1, na.rm = TRUE) > 0, 1, ifelse(rowSums(!is.na(parts)) > 0, 0, NA))
  }
  data.frame(SEQN = d$SEQN, asthma = get("MCQ010"), heart_failure = get("MCQ160B"), chd = get("MCQ160C"),
             angina = get("MCQ160D"), heart_attack = get("MCQ160E"), stroke = get("MCQ160F"),
             arthritis = get("MCQ160A"), copd = copd, cancer = get("MCQ220"),
             liver_condition = get("MCQ160L"), gallstones = get("MCQ550"))
}

# Cardiovascular disease: ever told of heart failure, coronary heart disease, angina, heart
# attack, or stroke.
cvd <- function(conditions) {
  parts <- conditions[, c("heart_failure", "chd", "angina", "heart_attack", "stroke")]
  ifelse(rowSums(parts == 1, na.rm = TRUE) > 0, 1, ifelse(rowSums(is.na(parts)) == 5, NA, 0))
}

# Leisure-time moderate and vigorous activity in minutes per week. 2007-2018 asked the GPAQ
# (PAQ650-PAQ675: days per week and minutes per day), 2017-March 2020 the same, and 2021-2023
# asks how often (per day, week, month, or year) and for how long (PAD790Q/U, PAD800, PAD810Q/U,
# PAD820). Vigorous minutes count twice in "mvpa_equivalent", as guidelines count them.
leisure_activity <- function(cycle) {
  d <- component("PAQ", cycle)
  if (!has(d, "PAD790Q") && !has(d, "PAQ650")) {
    # 1999-2006 asked other activity questions; only sedentary minutes, where asked, carry over.
    return(data.frame(SEQN = d$SEQN, leisure_moderate = NA_real_, leisure_vigorous = NA_real_, mvpa_equivalent = NA_real_,
                      sedentary_minutes = if (has(d, "PAD680")) ifelse(d$PAD680 %in% c(7777, 9999), NA, d$PAD680) else NA_real_))
  }
  if (has(d, "PAD790Q")) {
    per_week <- function(q, u) {
      q <- ifelse(q %in% c(7777, 9999), NA, q)
      factor <- ifelse(u %in% "D", 7, ifelse(u %in% "W", 1, ifelse(u %in% "M", 7 / 30.4375, ifelse(u %in% "Y", 7 / 365.25, NA))))
      ifelse(q %in% 0, 0, q * factor)
    }
    minutes <- function(m) ifelse(m %in% c(7777, 9999), NA, m)
    moderate <- per_week(d$PAD790Q, as.character(d$PAD790U)) * minutes(d$PAD800)
    vigorous <- per_week(d$PAD810Q, as.character(d$PAD810U)) * minutes(d$PAD820)
    moderate[d$PAD790Q %in% 0] <- 0
    vigorous[d$PAD810Q %in% 0] <- 0
  } else {
    days <- function(answer, n) ifelse(answer %in% 2, 0, ifelse(answer %in% 1 & n %in% 1:7, n, NA))
    minutes <- function(m) ifelse(m %in% c(7777, 9999), NA, m)
    moderate <- days(d$PAQ665, d$PAQ670) * ifelse(d$PAQ665 %in% 2, 0, minutes(d$PAD675))
    vigorous <- days(d$PAQ650, d$PAQ655) * ifelse(d$PAQ650 %in% 2, 0, minutes(d$PAD660))
  }
  data.frame(SEQN = d$SEQN, leisure_moderate = moderate, leisure_vigorous = vigorous,
             mvpa_equivalent = moderate + 2 * vigorous,
             sedentary_minutes = ifelse(d$PAD680 %in% c(7777, 9999), NA, d$PAD680))
}

# Liver elastography (2017 on): median CAP (dB/m), median stiffness (kPa), the stiffness
# IQR/median ratio, and exam status (1 complete, 2 partial, 3 ineligible, 4 not done).
elastography <- function(cycle) {
  d <- component("LUX", cycle)
  data.frame(SEQN = d$SEQN, cap = d$LUXCAPM, stiffness = d$LUXSMED, stiffness_iqr_ratio = d$LUXSIQRM,
             elastography_status = d$LUAXSTAT)
}

# Hepatitis B surface antigen (LBDHBG), hepatitis C antibody and RNA (LBXHCR): 1 positive.
hepatitis <- function(cycle) {
  b <- component("HEPBD", cycle)
  c <- component("HEPC", cycle)
  out <- data.frame(SEQN = union(b$SEQN, c$SEQN))
  out$hbsag <- b$LBDHBG[match(out$SEQN, b$SEQN)]
  # Hepatitis C antibody: LBDHCV through 2011-2012, LBDHCI from 2017 (2013-2016 released none).
  antibody <- intersect(c("LBDHCI", "LBDHCV"), names(c))
  out$hcv_antibody <- if (length(antibody) > 0) c[[antibody[1]]][match(out$SEQN, c$SEQN)] else NA
  out$hcv_rna <- if (has(c, "LBXHCR")) c$LBXHCR[match(out$SEQN, c$SEQN)] else NA
  out
}

crp <- function(cycle) {
  d <- component(if (cycle %in% c("2015-2016", "2017-2018", "2017-2020", REPLICATION_CYCLE)) "HSCRP" else "CRP", cycle)
  # hs-CRP in mg/L from 2015; CRP in mg/dL before 2011 (LBXCRP), converted to mg/L.
  value <- if (has(d, "LBXHSCRP")) d$LBXHSCRP else d$LBXCRP * 10
  data.frame(SEQN = d$SEQN, crp = value)
}

ferritin <- function(cycle) {
  d <- component("FERTIN", cycle)
  data.frame(SEQN = d$SEQN, ferritin = d$LBXFER)
}

# A categorical covariate with missing values kept as their own level ("Unclear"), as
# EmpowerStats papers in this literature report them.
with_unclear <- function(x, levels = sort(unique(x[!is.na(x)]))) {
  factor(ifelse(is.na(x), "unclear", as.character(x)), levels = c(as.character(levels), "unclear"))
}

# Merges data frames by SEQN, keeping every row of the first.
merge_all <- function(...) Reduce(function(a, b) merge(a, b, by = "SEQN", all.x = TRUE), list(...))

# Arthritis and its type. The type question's codes changed: 1999-2008 ask MCQ190 and 2009-2010
# MCQ191, with 1 = rheumatoid arthritis and 2 = osteoarthritis; from 2011 on (2021-2023 included)
# MCQ195 has 1 = osteoarthritis and 2 = rheumatoid arthritis. Returns told_arthritis (MCQ160A)
# and osteoarthritis and rheumatoid flags (1 yes, 0 no among those answering, NA otherwise).
arthritis <- function(cycle) {
  d <- component("MCQ", cycle)
  told <- yes(d$MCQ160A)
  if (has(d, "MCQ195")) {
    type <- d$MCQ195; oa_code <- 1; ra_code <- 2
  } else {
    type <- if (has(d, "MCQ191")) d$MCQ191 else d$MCQ190
    oa_code <- 2; ra_code <- 1
  }
  answered <- told %in% 0 | type %in% c(1, 2, 3, 4)
  data.frame(SEQN = d$SEQN, told_arthritis = told,
             osteoarthritis = ifelse(told %in% 1 & type %in% oa_code, 1, ifelse(answered, 0, NA)),
             rheumatoid = ifelse(told %in% 1 & type %in% ra_code, 1, ifelse(answered, 0, NA)))
}