## ----------------------------------------------------------------------------
## From Yale's Public Health Data Science and Data Equity (DSDE) Team
## 
## Workshop: Code Smarter, Not Harder: Unlocking AI-Assisted Coding
##  Authors: Shelby Golden, M.S.
##     Date: 2026-03-09
## 
##    R version: 4.5.2
## renv version: 1.1.5
## 
## 
## Description: Thomas et al. 2023 analyzed NHANES cycles from 1999-2017 
##              and prepared datasets with specific metrics for their analysis. 
##              This script is a worked-through example guiding you through
##              the process of preparing the Mediterranean diet scores and
##              regenerating the plot from the paper. The goal is to focus
##              on how AI can assist you in the coding process, and not on
##              the specific details of the analysis. Therefore, guides are
##              provided to help you outline your analysis plan, and a solution
##              script is provided for you to refer to if you get stuck at 
##              any point. 
##
##              Refer to the Book of Workshops webpage for this example for 
##              summary details about the paper, the data, and questions to 
##              consider as you work through the example.


## ----------------------------------------------------------------------------
## SET UP THE ENVIRONMENT
## renv() will install all of the packages and their correct version used here.
#renv::init()
renv::restore()

# Install development version of package (not yet on CRAN). Only necessary
# during the initial environment set up.
#install.packages("devtools")
#devtools::install_github("dayoonkwon/BioAge")

## Load in the R packages used in this script from the project library.
suppressPackageStartupMessages({
  library("nhanesdata") # Access compiled NHANES data
  library("readr")      # For reading in the data
  library("tidyr")      # For tidying data 
  library("dplyr")      # For data manipulation
  library("stringr")    # For string manipulation
  library("ggplot2")    # For creating static visualizations
  library("gridExtra")  # Format multiple plots together
})


## This function will check if an element is not in a vector.
"%!in%" <- function(x,y)!('%in%'(x,y))

# This function will return NA if all values in a vector are NA, otherwise 
# it will return the sum of the vector with NA values removed.
sum_or_na <- function(x) {
  if (all(is.na(x))) NA_real_ else sum(x, na.rm = TRUE)
}




## ----------------------------------------------------------------------------
## LOAD IN THE DATA

# Load demographics data and reference the Codebook (data dictionary) URL
demo <- read_nhanes("demo")
demo_var <- get_url("demo")

# Load dietary individual food data for day 1 and 2, and reference the Codebook URLs.
# Simultaneously select the variables of interest and merge with demographics.
diet1 <- read_nhanes("dr1iff") |>
  select(seqn, dr1ifdcd, dr1ikcal, dr1iprot, dr1icarb, dr1isugr, dr1itfat, 
         dr1imfat, dr1isfat, dr1icalc, dr1ialco) |>
  left_join(demo |> select(seqn, riagendr), by = c("seqn"))

diet2 <- read_nhanes("dr2iff") |>
  select(seqn, dr2ifdcd, dr2ikcal, dr2iprot, dr2icarb, dr2isugr, dr2itfat, 
         dr2imfat, dr2isfat, dr2icalc, dr2ialco) |>
  left_join(demo |> select(seqn, riagendr), by = c("seqn"))

diet1_var <- get_url("dr1iff_l")
diet2_var <- get_url("dr2iff_l")


# Load the Mediterranean diet food group reference data, based on the USDA Food and 
# Nutrient Database for Dietary Studies (FNDDS) codes.
mediterranean_groups <- read_csv(file = "https://raw.githubusercontent.com/ysph-dsde/Book-of-Workshops/refs/heads/main/Workshops/Code-Smarter-Not-Harder/Data/Mediterranean%20Diet%20Food%20Groups_Trichopoulou%20et%20al%202003.csv", show_col_types = FALSE,
                                 col_types = cols(code = col_character()))

# Load the provided solution for the Mediterranean diet scoring.
MeDi_score_soln <- read_csv(file = "https://raw.githubusercontent.com/ysph-dsde/Book-of-Workshops/refs/heads/main/Workshops/Code-Smarter-Not-Harder/Data/MeDi%20Scores.csv", show_col_types = FALSE,
                            col_types = readr::cols(
                              seqn    = readr::col_double(),
                              avg_kcal= readr::col_double(),
                              sd_kcal = readr::col_double(),
                              MeDi    = readr::col_double()
                            ))

# Load the already compiled dataset.
df <- read_csv(file = "https://raw.githubusercontent.com/ysph-dsde/Book-of-Workshops/refs/heads/main/Workshops/Code-Smarter-Not-Harder/Data/Complete%20Dataset.csv", show_col_types = FALSE)



## ----------------------------------------------------------------------------
## MEDITERRRANEAN DIET SCORING

## Thomas et al. (2023) assessed participants' Mediterranean diet adherence 
## using a "MeDi score":
##    1) For beneficial foods (vegetables, potatoes, fruits, legumes and nuts, 
##       fish, cereals), participants receive +1 point if their caloric intake 
##       meets or exceeds the median.
## 
##    2) For detrimental foods (meat, poultry, dairy products), participants 
##       receive +1 point if their caloric intake is below the median.
## 
##    3) For alcohol, +1 point is awarded for mild-to-moderate consumption (0–1 
##       drink per day for females; 0–2 drinks per day for males). The  National
##       Institute on Alcohol Abuse and Alcoholism (NIAAA) defines one standard
##       drink as 14 grams of pure alcohol.
## 
##    4) For fats, particiants with a high monosaturated to saturated ratio
##       receive +1 point if their caloric intake meets or exceeds the median.
##
##    5) Sex-specific medians are calculated for each food category as reference 
##       values. Alcohol and fat consumption is scored using different criteria.
## 
## Let's first examine the Mediterranean diet food group reference data to understand
## how the FNDDS codes are categorized. The "code" column contains the first two
## digits of the FNDDS codes, which correspond to broad food categories.

# Summarize the number of FNDDS groups that fall into each Mediterranean diet category 
# (beneficial, detrimental, or N/A).
scoring_groups <- mediterranean_groups |>
  group_by(med_group, benefice) |>
  summarise(count = n(), .groups = "drop") |>
  pivot_wider(names_from = benefice, values_from = count, values_fill = list(count = 0))

scoring_groups |> filter(Yes > 0)
scoring_groups |> filter(No > 0)
scoring_groups |> filter(`N/A` > 0)

# Some participants only completed one 24-hour dietary recall (either Day 1 or 
# Day 2). We process each day separately, then merge them when averaging dietary 
# records across days.
all(diet1$seqn %in% diet2$seqn == TRUE)
all(diet2$seqn %in% diet1$seqn == TRUE)

# Annotate individual food records using the Mediterranean diet classification 
# scheme described above.
diet1_annotated <- diet1 |>
  mutate(
    dr1ifdcd = as.character(format(dr1ifdcd, scientific = FALSE)),
    first_two_digits = str_sub(dr1ifdcd, 1, 2),
  ) |>
  left_join(mediterranean_groups, by = c("first_two_digits" = "code"))

diet2_annotated <- diet2 |>
  mutate(
    dr2ifdcd = as.character(format(dr2ifdcd, scientific = FALSE)),
    first_two_digits = str_sub(dr2ifdcd, 1, 2),
  ) |>
  left_join(mediterranean_groups, by = c("first_two_digits" = "code"))

# Confirm that there were no missed matches after joining the tables
diet1_annotated |> filter(med_group %in% NA) |> pull(first_two_digits) |> unique()
diet2_annotated |> filter(med_group %in% NA) |> pull(first_two_digits) |> unique()

# As a sanity check, we examine the distribution of key nutrients (protein, 
# carbohydrates, sugar, total fat, calcium, and alcohol) to confirm that 
# Mediterranean food groups are correctly associated. Obvious deviations would 
# indicate potential issues with the joining process that require investigation.
diet1_annotated |>
  group_by(med_group) |>
  summarize(
    avg_prot = mean(dr1iprot, na.rm = TRUE),
    avg_carb = mean(dr1icarb, na.rm = TRUE),
    avg_sugar = mean(dr1isugr, na.rm = TRUE),
    avg_fat = mean(dr1itfat, na.rm = TRUE),
    avg_cal = mean(dr1icalc, na.rm = TRUE),
    avg_alc = mean(dr1ialco, na.rm = TRUE),
    .groups = 'drop'
  )


## Now we're ready to start preparing the Mediterranean diet scores using the criteria 
## outlined above. Use the guides provided below to help outline your analysis plan.
## If at any point you get stuck, refer to the provided solution script for guidance 
## on how to proceed, or patch solutions in when you need it to proceed.
## 
## This section is intended for you to try solving these steps using the AI assistant.
## The solution is already provided at the end for use in the next section.

# a. Explore the dataset to determine what processing considerations are required to 
#    correctly prepare the data.



# b. Calculate the total monosaturated fat, saturated fat, and grams of alcohol each
#    participant consumed on each day a 24-hour recall was recorded. HINT: Use the
#    sum_or_na() function to ensure that participants with all missing values are
#    assigned NA instead of 0.



# c. Combine the two days into one table by = c("seqn", "riagendr") and calculate
#    the average and standard deviation of monosaturated fat, saturated fat, and 
#    grams of alcohol the participant consumed. HINT: You will need to calculate
#    these values over each row.
 


# d. Calculate the monosturated to saturated fat ratio for each participant and 
#    convert grams of alcohol to nearest whole drinks consumed using the NIAAA 
#    definition of a standard drink. 
# 
#    HINT: If you want to calculate the mean and standard deviation of ratio
#    assume that consumption of monosaturated and unsaturated fats are independent,
#    i.e. Var(X / Y) = E[X^2] * E[1/Y^2] - ( E[X] * E[1/Y] )^2 
#    and E(X/Y) ≈ E[X] * E[1/Y].
# 
#    If you calculate the standard deviation of drinks simply use SD(X) / 14.



# e. Calculate the total caloric intake for each participant by beneficial and 
#    detrimental food groups, i.e. by = c("seqn", "riagendr", "med_group", "benefice").
#    HINT: Use the sum_or_na() function to ensure that participants with all missing 
#    values are assigned NA instead of 0.



# g. Combine the two days into one table by = c("seqn", "riagendr", "med_group", 
#    "benefice") and calculate the average and standard deviation of total calories
#    consumed for each group. HINT: Simply calculate these values over each row.



# h. Calculate the threshold for scoring: mean monosturated to saturated fat ratio
#    and mean beneficial and detrimental food group by sex. HINT: Simply use the
#    median() function to calculate the threshold for scoring.



# i. Finally, score each participant's Mediterranean diet adherence using the criteria 
#    outlined above. Note that the provided solution assumes that you kept the
#    beneficial and detrimental food groups separate. HINT: Use functions like case_when() 
#    and if_else() to apply the scoring criteria for specific conditions.



# j. How well does your final MeDi score match the provided solution? Join your calculation
#    with the provided solution, MeDi_score_soln, and calculate their difference.




## ----------------------------------------------------------------------------
## SUMMARY STATISTICS

## In this section we would like to assess how evenly distributed possible
## confounders are across the different activity and MeDi groups. We will then
## create a balanced sample of participants across the different activity and
## MeDi groups for plotting.

## a. Age groupings and MeDi categories need to be defined and added to the dataset.
##    Age: [20,40), [40,60), [60,+) and MeDi: Low (0-3), Moderate (4-5), High (6-9).



## b. Calculate the number of people who fall into each stratification: age group, 
##    BMI category, and sex.



## c. Use the smallest n to randomly select the same number of participants from
##    each stratification group. HINT: To make it easier to plot with and without
##    the balanced sample, create a new variable that indicates whether the
##    participant is in the balanced sample or not, instead of subsetting the data.




## ----------------------------------------------------------------------------
## REGENERATE PLOT

## In order to regenerate the plot, we need to create a new variable that combines 
## the activity level and MeDi category for each participant. We also need to ensure 
## that the levels of this new variable are ordered correctly for plotting.

# Ordering the combined levels, used on the x-axis
desired_ordering <- c(
 "Sedentary - Low",  "Sedentary - Moderate",  "Sedentary - High",
 "Low - Low",  "Low - Moderate",  "Low - High",
 "Moderate - Low",  "Moderate - Moderate",  "Moderate - High",
 "High - Low",  "High - Moderate",  "High - High"
)

# Apply the factor levels and create the combined variable for plotting
df <- df |>
  mutate(
    activity_level = factor(activity_level, levels = c("Sedentary", "Low", "Moderate", "High")),
    activity_MeDi = {
      combo <- str_c(activity_level, MeDi_category, sep = " - ")
      factor(combo,
             levels = c(intersect(desired_ordering, unique(combo)),
                        setdiff(unique(combo), desired_ordering)))
    },
    bmi_category = factor(bmi_category, levels = c("Normal weight", "Overweight", "Obesity"))
  ) |>
  relocate(activity_MeDi, .before = in_balanced)


## Each plot that will be generated has a few key components:
##    1) The x-axis will show the combined activity level and MeDi category variable.
##    2) The y-axis will show the mean advancement_modified_std for each group,
##       with error bars representing the 95% confidence interval.
##    3) The points and error bars will be colored by activity level.

## a. Calculate the mean advancement_modified_std and 95% confidence interval for 
##    each group (by = c("activity_MeDi", "activity_level")). HINT: Use the 
##    t-distribution to calculate the confidence intervals.



## b. Do the same plot above, but facet by age group (by = c("age_group",
##    "activity_MeDi", "activity_level")).



## c. Do the same plot above, but facet by BMI category (by = c("bmi_category", 
##    "activity_MeDi", "activity_level")).

