---
title: "Run Summary Statistics"
author: "Vy Nguyen"
date: "1/12/2023"
output:
  pdf_document: default
  html_document: default
  word_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Goal 
We will learn how to calculate summary statistics using the chemical dataset. We will specifically learn to 
1. calculate the distribution statistics for one variable
2. stratify the distribution of one variable by a categorical variable (e.g. sex)
3. calculate the distribution statistics for one variable while accounting for NHANES sampling design
4. stratify the distribution of one variable by a categorical variable (e.g. sex) while accounting for NHANES sampling design
5. calculate the distribution statistics for multiple variables
6. calculate the distribution statistics for multiple variables while accounting for NHANES sampling design

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Install packages to run analysis
You only need to run this ONCE!
```{r, echo = FALSE}
install.packages("survey") # Needed to run a regression model that accounts for the NHANES sampling design
install.packages("tidyverse") # Needed to pipe functions together
install.packages("broom") # Needed to easily extract regression statistics
```

## Load packages
Run this so you have access to the functions in these packages.
```{r, echo = FALSE}
library("survey") 
library("tidyverse") 
library("broom") 
```

## Loading NHANES datasets
Make sure you set your working directory to be the folder containing file for the NHANES data. 
```{r, echo = FALSE}
load("./w - nhanes_1988_2018.RData")
```


## Merge NHANES datasets
Please use any dataset that ends with "_clean". 
Let's merge the demographics, chemicals, and weights datasets together using tidyverse syntax.
```{r, echo = FALSE}
nhanes_merged <- full_join(demographics_clean, 
                           chemicals_clean, 
                           by = c("SEQN",
                                  "SEQN_new", 
                                  "SDDSRVYR")) %>%
  full_join(., 
            weights_clean, 
            by = c("SEQN",
                   "SEQN_new", 
                   "SDDSRVYR")) 
```

## Data dictionary
Take a look at the dictionary to know the codename and description of the variables
```{r, echo = FALSE}
View(df_dictionary)
```

## Ensure that sex and race/ethnicity variables are categorical
factor() is a base R function that encode a variable as categorical. 
relevel() is a base R function that set the reference group. 

We ensured that sex (RIAGENDR) is a categorical variable.
We ensured that race/ethnicity (RIDRETH1) is a categorical variable with the reference group as Non-Hispanic Whites (3). 
```{r, echo = FALSE}
nhanes_merged <- nhanes_merged %>%
  mutate(RIAGENDR = factor(RIAGENDR)) %>%
  mutate(RIDRETH1 = factor(RIDRETH1) %>%
           relevel(., ref = 3))
```

## Choose only variables that you will use

For the portion on calculating summary statistics for multiple variables, we'll iterate over the metal biomarkers to get distribution statistics for each metal biomarkers. 
So let's include the metal biomarkers and their corresponding weight codenames in this NHANES subset to prepare for subsequential analysis.
```{r, echo = FALSE}
metals_codename <- df_dictionary %>%
  filter(chemical_family == "Metals") %>%
  filter(grepl("replicate", variable_description_use) == FALSE) %>%
  pull(variable_codename_use) 
  
# Corresponding weight codenames for the chemical biomarkers are the concatenation of "WT_" and the codename of the chemical biomarker
weights_codename_for_metals <- paste("WT_"
                                     , metals_codename
                                     , sep = "")

nhanes_subset <- nhanes_merged %>%
  select("SEQN",
         "SEQN_new", 
         "SDDSRVYR",
         "SDMVPSU",
         "SDMVSTRA",
         all_of(metals_codename),
         all_of(weights_codename_for_metals),
         "RIDAGEYR",
         "RIAGENDR",
         "RIDRETH1")

```

## Distribution statistics for one variable

Let's calculate the distribution statistics for blood lead
```{r, echo = FALSE}

dataset_lead <- nhanes_subset %>%
  select("SEQN",
         "SEQN_new", 
         "SDDSRVYR",
         "SDMVPSU",
         "SDMVSTRA",
         "WT_LBXBPB",
         "LBXBPB",
         "RIDAGEYR",
         "RIAGENDR",
         "RIDRETH1") %>%
  na.omit(.)

df_stats_lead <- dataset_lead %>%
  summarise(min = min(LBXBPB),
            perc_10 = quantile(LBXBPB, probs = 0.1),
            perc_25 = quantile(LBXBPB, probs = 0.25),
            median = median(LBXBPB),
            mean = mean(LBXBPB),
            perc_75 = quantile(LBXBPB, probs = 0.75),
            perc_90 = quantile(LBXBPB, probs = 0.9),
            max = max(LBXBPB),
            num_with_measurements = n())
```

Let's stratify the distribution statistics for blood lead by a categorical variable (e.g. sex)
```{r, echo = FALSE}

df_stats_lead_sex <- dataset_lead %>%
  mutate(sex = ifelse(RIAGENDR == 1,
                      "males", 
                      "females")) %>%
  group_by(sex) %>%
  summarise(min = min(LBXBPB),
            perc_10 = quantile(LBXBPB, probs = 0.1),
            perc_25 = quantile(LBXBPB, probs = 0.25),
            median = median(LBXBPB),
            mean = mean(LBXBPB),
            perc_75 = quantile(LBXBPB, probs = 0.75),
            perc_90 = quantile(LBXBPB, probs = 0.9),
            max = max(LBXBPB),
            num_with_measurements = n()) %>%
  ungroup(.)

```


## Distribution statistics for one variable accounting for NHANES sampling design

The survey weights for a given two-year cycle in the set of cycles 3-10 (2003-2018) is representative of a national sample (i.e. the survey weights in a given cycle should sum to the population size of the US). 
The survey weights for cycles 1 and 2 (1999-2002) TOGETHER is representative of a national sample. 

When we combine survey cycles, let's suppose for two cycles (cycle 3 and 4), then the survey weights are going to be representative of two national samples, so any estimates are going to be representative of two US population. This doesn't make sense! So we multiply the survey weights by 1/2 (i.e. 1 over the number of studied cycles), and then the adjusted survey weights are representative of ONLY ONE national sample. 

More examples can be found under the tab "When and How to Construct Weights When Combining Survey Cycles" at https://wwwn.cdc.gov/nchs/nhanes/tutorials/weighting.aspx 

Sanity check:  sum the adjusted survey weights and it should be close to the US population size. 

```{r, echo = FALSE}

# Determine the cycles with measurements for blood lead
unique_cycles_lead <- dataset_lead %>%
  pull(SDDSRVYR) %>%
  unique(.)

# Determine the total number of cycles with measurements for lead
num_cycles_lead <- length(unique_cycles_lead)

# Calculate the adjusted survey weights
dataset_lead <- dataset_lead %>%
  mutate(adjusted_weights = ifelse(SDDSRVYR %in% c(1,2), 
                                   (2/num_cycles_lead)*WT_LBXBPB,
                                   (1/num_cycles_lead)*WT_LBXBPB))

# Checking our calculations
dataset_lead %>%
  select(SDDSRVYR,
         WT_LBXBPB,
         adjusted_weights) %>%
  unique(.) %>%
  mutate(multiplier = adjusted_weights/WT_LBXBPB) %>%
  View(.)

# Account for the sampling design
dsn_lead <- svydesign(ids = ~SDMVPSU, 
                      strata = ~SDMVSTRA,
                      weights = ~adjusted_weights, 
                      nest = TRUE, 
                      data = dataset_lead)

# Calculate distribution statistics for lead and these statistics are generalizable to the entire US population
df_stats_lead_svy <- dataset_lead %>%
  summarise(min = svyquantile(~LBXBPB, dsn_lead, quantiles = 0)[1],
            perc_10 = svyquantile(~LBXBPB, dsn_lead, quantiles = 0.1)[1],
            perc_25 = svyquantile(~LBXBPB, dsn_lead, quantiles = 0.25)[1],
            median = svyquantile(~LBXBPB, dsn_lead, quantiles = 0.5)[1],
            mean = svymean(~LBXBPB, dsn_lead)[1],
            perc_75 = svyquantile(~LBXBPB, dsn_lead, quantiles = 0.75)[1],
            perc_90 = svyquantile(~LBXBPB, dsn_lead, quantiles = 0.90)[1],
            max = svyquantile(~LBXBPB, dsn_lead, quantiles = 1)[1]) %>%
  ungroup(.)



```

Let's stratify the distribution statistics for blood lead by a categorical variable (e.g. sex)
```{r, echo = FALSE}

df_stats_leads_svy_sex <- dataset_lead %>%
  summarise(min = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svyquantile, quantiles = 0, keep.var=FALSE)[,2],
            perc_10 = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svyquantile, quantiles = 0.1, keep.var=FALSE)[,2],
            perc_25 = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svyquantile, quantiles = 0.25, keep.var=FALSE)[,2],
            median = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svyquantile, quantiles = 0.5, keep.var=FALSE)[,2],
            mean = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svymean)[,2],
            perc_75 = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svyquantile, quantiles = 0.75, keep.var=FALSE)[,2],
            perc_90 = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svyquantile, quantiles = 0.9, keep.var=FALSE)[,2],
            max = svyby(~LBXBPB, ~RIAGENDR, dsn_lead, svyquantile, quantiles = 1, keep.var=FALSE)[,2]) %>%
  mutate(RIAGENDR = unique(dataset_lead$RIAGENDR)) %>%
  relocate(RIAGENDR)

```

## Distribution statistics for multiple variables

We will iterate over each metal biomarker, calculate the statistics for each metal biomarker, and stitch the dataset of statistics together. 

To iterate over each metal biomarkers, we will use map() on the codenames of the metal biomarkers.

To calculate the statistics for each metal biomarker, we will customize a function (e.g. calculate_distribution_stats) to calculate the percentiles.

To stitch the dataset of statistics together, we will use bind_rows().
```{r, echo = FALSE}

calculate_distribution_stats <- function(x, 
                                         df_nhanes,
                                         dataset_dictionary)
{
  # Print out the codename of a given chemical biomarker
  # This helps with debugging to detect which chemical biomarker is causing problems.
  print(x)
  
  # Use the data dictionary to extract the name of the chemical
  chemical_name <- dataset_dictionary %>%
    filter(variable_codename_use == x) %>%
    pull(variable_description_use) %>%
    unique(.)
  
  # Define a subset to include any participants with measurements for the chemical biomarker
  subset_x <- df_nhanes %>%
    select("SEQN",
           "SEQN_new", 
           "SDDSRVYR",
           all_of(x)) %>%
    na.omit(.)
  
  # Let's change the column name of the chemical biomarker to "chem" to make it easier to calculate the statistics. 
  index_chem <- which(colnames(subset_x) == x)
  colnames(subset_x)[index_chem] <- "chem"
  
  df_stats_x <- subset_x %>%
    summarise(min = min(chem),
              perc_10 = quantile(chem, probs = 0.1),
              perc_25 = quantile(chem, probs = 0.25),
              median = median(chem),
              mean = mean(chem),
              perc_75 = quantile(chem, probs = 0.75),
              perc_90 = quantile(chem, probs = 0.9),
              max = max(chem),
              num_with_measurements = n()) %>%
    mutate(variable_codename_use = x) %>% # Include the codename to know which statistics belongs to which chemical when we stitch the datasets together
    mutate(variable_description_use = chemical_name) %>% # Include the chemical name
    relocate(variable_codename_use) %>% # Make the column of chemical codenames be the 1st column
    relocate(variable_description_use, 
             .after = variable_codename_use) # Make the column of chemical names be the column after the column of codenames (or the 2nd column)
  
  # For map(), we must return a data frame.
  return(df_stats_x)
}

df_stats_metals <- metals_codename %>%
    map(., 
        calculate_distribution_stats, 
        nhanes_subset,
        df_dictionary) %>%
    bind_rows(.) 

```

## Distribution statistics for multiple variables accounting for NHANES sampling design

Again, we will iterate over each metal biomarker, calculate the statistics for each metal biomarker, and stitch the dataset of statistics together. 

Now the difference is that when we customize the function for each metal biomarker, we need to adjust the survey weights, so that our estimates will be representative of ONLY ONE national sample.

```{r, echo = FALSE}

calculate_distribution_stats_svy <- function(x, 
                                         df_nhanes,
                                         dataset_dictionary)
{
  # Print out the codename of a given chemical biomarker
  print(x)
  
  # Define the weight codename for a given chemical biomarker
  weight_codename <- paste("WT_",
                           x,
                           sep = "")
  
  # Use the data dictionary to extract the name of the chemical
  chemical_name <- dataset_dictionary %>%
    filter(variable_codename_use == x) %>%
    pull(variable_description_use) %>%
    unique(.)
  
  # Define a subset to include any participants with measurements for the chemical biomarker
  subset_x <- df_nhanes %>%
    select("SEQN",
           "SEQN_new", 
           "SDDSRVYR",
           "SDMVPSU",
           "SDMVSTRA",
           all_of(x),
           all_of(weight_codename)) %>%
    na.omit(.)
  
  # Extract the unadjusted survey weights from the subset
   unadjusted_weights <- subset_x %>%
    pull(all_of(weight_codename))
  
  # Determine the cycles that the metal biomarker was measured
  unique_cycles_x <- subset_x %>%
    pull(SDDSRVYR) %>%
    unique(.)
  # print(unique_cycles_x)

  # Determine the total number of cycles
  num_cycles <- length(unique_cycles_x)
  # print(num_cycles)
  
  # Determine if the chemical biomarker had measurements in both cycle 1 and cycle 2
  indicator_cycles <- ifelse(1 %in% unique_cycles_x & 2 %in% unique_cycles_x,
                            "yes", 
                            "no")
  # print(indicator_cycles)
  
  # Calculate the adjusted weights
  if(indicator_cycles == "yes")
  {
    adjusted_weights <- ifelse(subset_x$SDDSRVYR %in% c(1,2), 
                               (2/num_cycles)*unadjusted_weights,
                               (1/num_cycles)*unadjusted_weights)
  } else {
    
    adjusted_weights <- (1/num_cycles)*unadjusted_weights
    
  }
  
  # Include the adjusted weights in the subset and exclude any weights that are 0 or else the code to calculate the statistics may not run
  subset_x <- subset_x %>%
    mutate(adjusted_weights = adjusted_weights) %>%
    filter(adjusted_weights != 0)
  
    # Let's change the column name of the chemical biomarker to "chem" to make it easier to calculate the statistics.
  index_chem <- which(colnames(subset_x) == x)
  colnames(subset_x)[index_chem] <- "chem"

  # Define the survey object to account for the sampling design
  dsn_x <- svydesign(ids = ~SDMVPSU, 
                      strata = ~SDMVSTRA,
                      weights = ~adjusted_weights, 
                      nest = TRUE, 
                      data = subset_x)
  
  # Calculate the distribution statistics
  # The [1] is necessary or else the column names will be changed to some less readable
  df_stats_x <- data.frame("min" = svyquantile(~chem, dsn_x, quantiles = 0)[1],
                           "perc_10" = svyquantile(~chem, dsn_x, quantiles = 0.1)[1],
                           "perc_25" = svyquantile(~chem, dsn_x, quantiles = 0.25)[1],
                           "median" = svyquantile(~chem, dsn_x, quantiles = 0.5)[1],
                           "mean" = svymean(~chem, dsn_x)[1] %>% unlist(.),
                           "perc_75" = svyquantile(~chem, dsn_x, quantiles = 0.75)[1],
                           "perc_90" = svyquantile(~chem, dsn_x, quantiles = 0.90)[1],
                           "max" = svyquantile(~chem, dsn_x, quantiles = 1)[1]) %>%
    mutate(variable_codename_use = x) %>%
    mutate(variable_description_use = chemical_name) %>%
    relocate(variable_codename_use) %>%
    relocate(variable_description_use,
             .after = variable_codename_use)
  
  return(df_stats_x)
}

df_stats_metals_svy <- metals_codename %>%
    map(., 
        calculate_distribution_stats_svy, 
        nhanes_subset,
        df_dictionary) %>%
    bind_rows(.) 

```

