# Indices that papers in this literature use as exposures, by their published formulas.

# Systemic immune-inflammation index: platelets x neutrophils / lymphocytes (counts in 10^3/uL).
sii <- function(platelets, neutrophils, lymphocytes) platelets * neutrophils / lymphocytes

# Weight-adjusted-waist index: waist (cm) / sqrt(weight (kg)).
wwi <- function(waist, weight) waist / sqrt(weight)

# Triglyceride-glucose index: ln(triglycerides (mg/dL) x fasting glucose (mg/dL) / 2).
tyg <- function(tg, glucose) log(tg * glucose / 2)

# Visceral adiposity index (Amato et al. 2010), triglycerides and HDL in mmol/L.
vai <- function(sex, waist, bmi, tg_mmol, hdl_mmol) {
  ifelse(sex == 1, waist / (39.68 + 1.88 * bmi) * (tg_mmol / 1.03) * (1.31 / hdl_mmol),
         waist / (36.58 + 1.89 * bmi) * (tg_mmol / 0.81) * (1.52 / hdl_mmol))
}

# Unit conversions to mmol/L.
tg_mmol <- function(tg_mg) tg_mg * 0.01129
cholesterol_mmol <- function(c_mg) c_mg * 0.02586

# Dietary Inflammatory Index (Shivappa et al. 2014, Public Health Nutr 17:1689): for each food
# parameter, z = (intake - global mean) / global SD, turned into a centered percentile
# 2 * pnorm(z) - 1, times the parameter's overall inflammatory effect score, summed. These are
# the 28 of its 45 parameters that NHANES's total-nutrient files supply, with the published
# constants (as the dietaryindex R package, MIT license, tabulates them).
DII_PARAMETERS <- data.frame(
  parameter = c("ALCOHOL", "VITB12", "VITB6", "BCAROTENE", "CAFFEINE", "CARB", "CHOLES", "KCAL", "TOTALFAT", "FIBER",
                "FOLICACID", "IRON", "MG", "MUFA", "NIACIN", "N3FAT", "N6FAT", "PROTEIN", "PUFA", "RIBOFLAVIN",
                "SATFAT", "SE", "THIAMIN", "VITA", "VITC", "VITD", "VITE", "ZN"),
  effect = c(-0.278, 0.106, -0.365, -0.584, -0.11, 0.097, 0.11, 0.18, 0.298, -0.663,
             -0.19, 0.032, -0.484, -0.009, -0.246, -0.436, -0.159, 0.021, -0.337, -0.068,
             0.373, -0.191, -0.098, -0.401, -0.424, -0.446, -0.419, -0.313),
  mean = c(13.98, 5.15, 1.47, 3718, 8.05, 272.2, 279.4, 2056, 71.4, 18.8,
           273, 13.35, 310.1, 27, 25.9, 1.06, 10.8, 79.4, 13.88, 1.7,
           28.6, 67, 1.7, 983.9, 118.2, 6.26, 8.73, 9.84),
  sd = c(3.72, 2.7, 0.74, 1720, 6.67, 40, 51.2, 338, 19.4, 4.9,
         70.7, 3.71, 139.4, 6.1, 11.77, 1.06, 7.5, 13.9, 3.76, 0.79,
         8, 25.1, 0.66, 518.6, 43.46, 2.21, 1.49, 2.19),
  stringsAsFactors = FALSE
)

# The NHANES nutrient (the suffix after DR1T/DR2T) behind each DII parameter, with its unit
# conversion: caffeine is in grams in the DII and in milligrams in NHANES; n-3 fat is ALA (18:3),
# stearidonic acid (18:4), EPA (20:5), DPA (22:5), and DHA (22:6); n-6 fat is linoleic (18:2) and
# arachidonic (20:4) acid; folic acid is DR?TFA; vitamin A is retinol activity equivalents.
DII_NUTRIENTS <- c("ALCO", "VB12", "VB6", "BCAR", "CAFF", "CARB", "CHOL", "KCAL", "TFAT", "FIBE", "FA", "IRON", "MAGN",
                   "MFAT", "NIAC", "P183", "P184", "P205", "P225", "P226", "P182", "P204", "PROT", "PFAT", "VB2", "SFAT",
                   "SELE", "VB1", "VARA", "VC", "VD", "ATOC", "ZINC")

dii_parameters <- function(diet) {
  get <- function(v) if (v %in% names(diet)) diet[[v]] else NA_real_
  data.frame(ALCOHOL = get("ALCO"), VITB12 = get("VB12"), VITB6 = get("VB6"), BCAROTENE = get("BCAR"),
             CAFFEINE = get("CAFF") / 1000, CARB = get("CARB"), CHOLES = get("CHOL"), KCAL = get("KCAL"),
             TOTALFAT = get("TFAT"), FIBER = get("FIBE"), FOLICACID = get("FA"), IRON = get("IRON"), MG = get("MAGN"),
             MUFA = get("MFAT"), NIACIN = get("NIAC"),
             N3FAT = get("P183") + get("P184") + get("P205") + get("P225") + get("P226"),
             N6FAT = get("P182") + get("P204"), PROTEIN = get("PROT"), PUFA = get("PFAT"), RIBOFLAVIN = get("VB2"),
             SATFAT = get("SFAT"), SE = get("SELE"), THIAMIN = get("VB1"), VITA = get("VARA"), VITC = get("VC"),
             VITD = get("VD"), VITE = get("ATOC"), ZN = get("ZINC"))
}

# The DII from a data frame of nutrient intakes named as DII_NUTRIENTS (for example the output of
# dietary_totals(cycle, DII_NUTRIENTS)), over `parameters` (by default every parameter the data
# have). A parameter missing for a person is left out of that person's sum.
dii <- function(diet, parameters = DII_PARAMETERS$parameter) {
  values <- dii_parameters(diet)
  score <- rep(0, nrow(values))
  for (p in parameters) {
    k <- DII_PARAMETERS[DII_PARAMETERS$parameter == p, ]
    contribution <- (2 * stats::pnorm((values[[p]] - k$mean) / k$sd) - 1) * k$effect
    score <- score + ifelse(is.na(contribution), 0, contribution)
  }
  score
}

# Composite Dietary Antioxidant Index (Wright et al. 2004, Am J Epidemiol 160:68): the sum over
# vitamins A, C, E, zinc, selenium, and carotenoids of (intake - mean) / SD, with the mean and SD
# from the study's own sample. `stats` is a list with mean and sd for each component, computed on
# the paper's cycles (constants) and reused unchanged in 2021-2023.
CDAI_COMPONENTS <- c("vitamin_a", "vitamin_c", "vitamin_e", "zinc", "selenium", "carotenoids")

cdai_components <- function(diet) {
  get <- function(v) if (v %in% names(diet)) diet[[v]] else NA_real_
  data.frame(vitamin_a = get("VARA"), vitamin_c = get("VC"), vitamin_e = get("ATOC"), zinc = get("ZINC"),
             selenium = get("SELE"), carotenoids = get("ACAR") + get("BCAR") + get("CRYP") + get("LYCO") + get("LZ"))
}

cdai <- function(components, stats) {
  Reduce(`+`, lapply(CDAI_COMPONENTS, function(k) (components[[k]] - stats[[k]]$mean) / stats[[k]]$sd))
}