---
title: "Obesity Rates in USA Analysis of 10,612 Adults from 2011-2018 National Health and Nutrition Examination Survey (NHANES)"
author: "Made with R Studio by Sobhan Mardan"
output: 
  flexdashboard::flex_dashboard:
    source_code: embed
    
      
---

```{r setup, include=FALSE}
library(flexdashboard)
library(magrittr) 
library(dplyr)
library(ggplot2)
library(tidyr)
library(ggprism)
library(nhanesA)
library(gtsummary)
library(survey)
library(trelliscopejs)
```


```{r, echo=FALSE, error=FALSE, include=FALSE, cache=FALSE, dev='jpeg'}
# #Demographic data cycles 2011-2018
# demo_2011_18 <- full_join(nhanes("Demo_G"),nhanes("Demo_H")) %>% 
#   full_join(full_join(nhanes("Demo_I"),nhanes("Demo_J"))) %>% 
#   select(SEQN, RIAGENDR, RIDAGEYR, RIDRETH1) 
# 
# 
# #Body measures data cycles 2011-2018
# body_2011_18 <- full_join(nhanes("BMX_G"), nhanes("BMX_H")) %>% 
#   full_join(full_join(nhanes("BMX_I"),nhanes("BMX_J"))) %>% 
#   select(SEQN, BMXWT, BMXHT, BMXBMI)
# 
# #Body Composition data cycles 2011-2018
# dxa_2011_18 <- full_join(nhanes("DXX_G"),nhanes("DXX_H")) %>% 
#   full_join(full_join(nhanes("DXX_I"),nhanes("DXX_J"))) %>% 
#    rowwise() %>% 
#   mutate(ALM = sum(c(DXDRLLE, DXDLLLE, DXDRALE, DXDLALE))/1000) %>% 
#   select(SEQN, DXDTOLE, DXDTOPF, ALM, DXXLSBMD) 
# 
# # Lab works data cycles 2011-2018
# lab_2011_18 <- full_join(nhanes("GHB_G"), nhanes("GHB_H")) %>% 
#   full_join(full_join(nhanes("GHB_I"),nhanes("GHB_J"))) %>%
#   mutate(diabetes_status = case_when(
#     LBXGH < 5.7 ~ "Normal",
#     LBXGH >= 5.7 & LBXGH <6.4 ~ "Prediabetic", 
#     LBXGH >= 6.5 ~ "Diabetic"
#   ), 
#   diabetes_status = factor(diabetes_status, 
#                            levels = c("Normal", "Prediabetic",
#                                       "Diabetic"))) %>% 
#   rename(a1c = "LBXGH")
# 
# # Individual cycles 
# 
# # 2011-2012
# cycle_2011_12<- 
#   nhanes("Demo_G") %>%
#   dplyr::pull(var = SEQN)
# 
# # 2013-14
# cycle_2013_14 <-
#   nhanes("Demo_H") %>% 
#   dplyr::pull(var = SEQN)
# 
# # 2015-16
# cycle_2015_16 <- 
#   nhanes("Demo_I") %>% 
#   dplyr::pull(var = SEQN)
# 
# # 2017-18
# cycle_2017_18 <- 
#   nhanes("Demo_J") %>% 
#   dplyr::pull(var = SEQN)
# 
# #Joining of data and Mutating data cycles 2011-2018
# df_2011_18<- full_join(demo_2011_18, body_2011_18, by = "SEQN") %>% 
#   full_join(dxa_2011_18, by = "SEQN") %>% 
#   full_join(lab_2011_18, by = "SEQN") %>% 
#   mutate(BMXHT = BMXHT/100, 
#          SMI = ALM/BMXHT^2, 
#          Gender = case_when(
#            RIAGENDR == 1 ~ "Male" ,
#            RIAGENDR ==2 ~ "Female"
#          ), 
#          Age_Group = case_when(
#            RIDAGEYR >= 8 & RIDAGEYR  <= 19 ~ "8-19", 
#            RIDAGEYR >= 20 & RIDAGEYR <= 29 ~ "20-29",
#            RIDAGEYR >= 30 & RIDAGEYR <= 39 ~ "30-39", 
#            RIDAGEYR >= 40 & RIDAGEYR <= 49 ~ "40-49", 
#            RIDAGEYR >= 50 & RIDAGEYR <= 59 ~ "50-59" 
#            ), 
#          BMI_class = case_when(
#            BMXBMI <18.5 ~ "Underweight", 
#            BMXBMI >=18.5 & BMXBMI <= 24.9 ~ "Normal", 
#            BMXBMI >= 25.0 & BMXBMI <= 29.9 ~ "Overweight", 
#            BMXBMI >= 30.0 & BMXBMI <= 34.9 ~ "Class I obese",
#            BMXBMI >=35 ~ "Class II-III obese" 
#          ), 
#          BMI_class = factor(BMI_class,levels=c("Underweight","Normal",
#                                                  "Overweight", 
#                                                "Class I obese",
#                                                "Class II-III obese")), 
#         Gender = factor(Gender, 
#                             levels = c("Male", "Female")),
#         Age_Group = factor(Age_Group,
#                                levels = c("8-19", "20-29", "30-39",
#                                           "40-49", "50-59")),
#         cycle_year = case_when(
#           SEQN %in% c(cycle_2011_12) ~ "2011-2012 wave", 
#           SEQN %in% c(cycle_2013_14) ~ "2013-2014 wave", 
#           SEQN %in% c(cycle_2015_16) ~ "2015-2016 wave", 
#           SEQN %in% c(cycle_2017_18) ~ "2017-2018 wave")
#         )%>% 
#   drop_na() 


df_2011_18<-
read.csv(file = "Data/df_2011_2018.csv") %>% 
   filter(RIDAGEYR>= 18) %>% 
  mutate(BMI_class = factor(BMI_class,levels=c("Underweight","Normal",
                                                 "Overweight",
                                               "Class I obese",
                                               "Class II-III obese")),
        Gender = factor(Gender,
                            levels = c("Male", "Female")),
        Age_Group = case_when(
    RIDAGEYR >= 18 & RIDAGEYR<= 19 ~ "18-19",
    TRUE ~ as.character(Age_Group)),
        Age_Group = factor(Age_Group,
                               levels = c("18-19", 
                                          "20-29", "30-39",
                                          "40-49", "50-59"))
  )


```



Column {data-width=650, .tabset} 
-----------------------------------------------------------------------
### **Trends in Obesity**

```{r, dev='jpeg', fig.width=20, fig.height=10,  fig.cap= "NHANES Data 2011-2018", dpi = 500}
 
attach(df_2011_18)
a<-
df_2011_18 %>% 
  group_by(BMI_class, cycle_year) %>% 
  summarise(n=n()) %>%
  arrange(BMI_class, cycle_year) %>% 
  bind_cols(
    round(prop.table(table(BMI_class, cycle_year), c(2))*100, 1) %>% 
  data.frame() %>%
    arrange(BMI_class, cycle_year) %>%
    select(Freq) 
  ) %>%
  ggplot(aes(x = cycle_year, y = n, fill = BMI_class))+ 
  geom_col(position = position_dodge2(), color = "black", size = 1) + 
  geom_text(aes(label = paste0(Freq, "%")), position = position_dodge(0.9), vjust = -1, size = 5)+ 
  scale_fill_brewer(name = "BMI Class")+ 
  scale_y_continuous(breaks = seq(0,1500,500))+
  theme_minimal(base_size = 30) + 
  theme(axis.text = element_text(colour = "black"),
        legend.title = element_blank()
        )+  
  labs(y = "Number of People", 
       x = "NHANES Data Collection Waves")

a



  
```


### **BMI vs Body Fat**

```{r, dev='jpeg', fig.width=20, fig.height=10,  fig.cap= "NHANES Data 2011-2018", dpi = 500}

attach(df_2011_18)
b<-
df_2011_18 %>% 
  ggplot(aes(x = BMXBMI, y = DXDTOPF))+ 
  geom_point(aes(fill = Gender), pch = 21, size= 7, alpha = .75)+ 
  scale_fill_manual(values = c("#3F6E9A", "#ACD4EC"))+
  geom_smooth(lty = 2, se = F, size = 2, color = "black" , method = "lm") + 
  labs(x = "Body Mass Index (kg/m2)", 
       y = "Body Fat (%)") +
  theme_minimal(base_size = 30) +
  theme(axis.text = element_text(color = "black"), 
        legend.text = element_text(face = "bold"), 
        legend.key.size = unit(1, "cm"),
        legend.title = element_blank()
        )
b


```


### **Body Fat and Age Group**

``````{r, dev='jpeg', fig.width=20, fig.height=10,  fig.cap= "NHANES Data 2011-2018", dpi = 500}
attach(df_2011_18)
c<-
df_2011_18 %>% 
  ggplot(aes(x = Age_Group, y = DXDTOPF, fill = Gender))+
  geom_boxplot(size = 2, outlier.alpha = 0) +
  scale_fill_manual(values = c("#3F6E9A", "#ACD4EC"))+
  labs(x = "Age Group (years)", 
       y = "Body Fat (%)") +
  theme_prism(base_size = 30) +
  scale_x_discrete(guide = "prism_bracket") + 
  scale_y_continuous(guide = "prism_offset", 
                     breaks = seq(0,60,10), 
                     expand = c(.1,2)
                     )+
  theme(axis.text = element_text(color = "black"), 
        legend.text = element_text(face = "bold"), 
        legend.key.size = unit(1, "cm"),
        legend.key.height = unit(3,"cm"),
        panel.grid = element_blank(),
        legend.title = element_blank()
        )
c
  
```


### **Summary Table**
```{r, echo=FALSE, error=FALSE, cache=FALSE, dev='jpeg'}

attach(df_2011_18)
# df_2011_18 %>%
#   select(Gender, BMXBMI, DXDTOPF, Age_Group) %>% 
#   tbl_summary(by = Age_Group) 
# 
tbl1<- 
df_2011_18 %>% 
  select(DXDTOPF,Gender, BMI_class, BMXBMI) %>% 
  tbl_strata(strata = Gender, ~.x %>% 
  tbl_summary(by = BMI_class, 
              type = all_continuous() ~ "continuous2" ,
              statistic = 
                list(all_continuous() ~ c("{mean} ± {sd}",
                                          "{median} ({p25}, {p75})",
                                          "{min} - {max}")),
               label = list(DXDTOPF ~ "Body Fat (%)",
                            BMXBMI ~ "BMI (kg/m2)"
                            ) 
  ) %>%
  modify_header(all_stat_cols() ~ "**{level}**<br>N = {n} <br> ({style_percent(p)}%)" ) %>%
  
  italicize_levels()
  ) 

tbl1  

             
#  des1<-
#  svydesign(ids = ~0 , data = df_2011_18)
#  summary(
# svyglm(formula = DXDTOPF ~ BMXBMI, design = des1)
#  ) 
#mod1<-
#cor.test(BMXBMI, DXDTOPF , method = "pearson")
# 
# mod1<- svyglm(formula = DXDTOPF ~ BMXBMI, design = des1)
#  
# tbl_regression(mod1, exponentiate = T) 
# mod1

```



Column {data-width=350, .tabset}
-----------------------------------------------------------------------

### **Purpose**

Obesity is a chronic medical condition that affects millions of people in the United States. Obesity is associated with many metabolic diseases including type 2 diabetes, high blood pressure, fatty liver etc. Using National Health and Nutrition Examination Survey from 2011 to 2018, I wanted to examine the following questions: 

* What is the prevalence of obesity by gender, age and the year in which data was collected? 

* Is Body Mass Index (BMI) a good indicator of body fat percentage as measured by the gold standard method dual energy x-ray absorptiometry (DEXA)? 

* What are the differences between body fat composition between those in different weight classes?
 
By using this dashboard my hope is to create an output for visualization and analysis of the publicly available data. 

The findings of this dashboard is for the purposes of personal interests, projects and teaching only. 




### **Gender and Obesity class** 

```{r, dev='jpeg', fig.cap= "NHANES Data 2011-2018", fig.width= 8, fig.height= 8, dpi=500}
attach(df_2011_18)
d<-
df_2011_18 %>% 
   group_by(BMI_class, Gender) %>% 
  summarise(n=n()) %>%
  arrange(BMI_class, Gender) %>% 
  bind_cols(
    round(prop.table(table(BMI_class, Gender), c(1))*100, 1) %>% 
  data.frame() %>%
    arrange(BMI_class, Gender) %>% 
    select(Freq)
  ) %>% 
  ggplot(aes(x = BMI_class, y = n, fill = Gender))+ 
  geom_col(color = "black", position = position_dodge(), size= 1)+ 
  scale_fill_manual(values = c("#3F6E9A", "#ACD4EC")) + 
  labs(x = "Body Mass Index Category",
       y = "Number of People") + 
  geom_text(aes(label = paste0(Freq, "%")), vjust = 1.5, position = position_dodge(.9), color = "white", fontface = "bold") +
  theme_prism(base_size = 15)+ 
  scale_x_discrete(guide = "prism_bracket", 
                   ) + 
  scale_y_continuous(guide = "prism_offset", 
                     breaks = seq(0,5000,500)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
       # plot.background = element_rect(fill = "#EEDBBD")
        legend.text = element_text(face="bold"),
        legend.title = element_blank())
d

```


### **Gender Proportions** 

```{r, dev='jpeg', fig.cap= "NHANES Data 2011-2018", fig.width= 8,dpi=500}
attach(df_2011_18)
e<-
df_2011_18 %>% 
  group_by(Age_Group, Gender) %>% 
  summarise(n=n()) %>% 
  arrange(Age_Group, Gender) %>% 
  bind_cols(
    round(prop.table(table(Age_Group, Gender), c(1))*100, 1) %>% 
  data.frame() %>% 
    arrange(Age_Group, Gender) %>% 
    select(Freq)
  ) %>% 
  #count(Age_Group, Gender, name = "n")#%>%  
  ggplot(aes(x = Age_Group,
             y = ifelse(test = Gender == "Male", 
                        yes = -n, 
                        no = n) , 
             fill = Gender)) + 
  geom_col(color = "black", width = .75, size = 1) + 
  coord_flip() + 
  scale_fill_manual(values = c("#3F6E9A", "#ACD4EC")) +  
  labs(x = "Age Group (years)", 
       y = "Number of People") + 
  geom_text(data = ~subset(., Gender=="Male"), aes(label = paste0(Freq, "%")), hjust = -.3, color = "white", fontface = "bold") +
  geom_text(data = ~subset(., Gender=="Female"), aes(label = paste0(Freq, "%")), hjust = 1.3, color = "white", fontface = "bold") +
  theme_prism(base_size = 15)+ 
  scale_x_discrete(guide = guide_prism_bracket(outside = FALSE)) + 
  scale_y_continuous(guide = "prism_offset" ,
                     breaks = seq(-5000,5000,500), 
                     labels = paste0(as.character(c(seq(5000, 0, -500), 
                                                     seq(500,5000,500)))), 
                     expand = c(0 , 500)) +  
  theme(strip.text = element_text(size = 15), 
        #plot.background = element_rect(fill = "#EEDBBD")
        legend.text = element_text(face="bold")) 
e

```

### **Body Fat and BMI Class and Gender Anova Table**
```{r, echo=FALSE, error=FALSE, cache=FALSE, dev='jpeg'}

attach(df_2011_18)
bmi_gender<-
  aov(DXDTOPF~BMI_class*Gender)
#  summary(bmi_gender)

poshoc<-
TukeyHSD(bmi_gender)
poshocbmigender<-
poshoc$`BMI_class:Gender` %>% 
  data.frame() %>% 
  mutate(diff = round(diff, 1), 
         lwr = round(lwr, 1), 
         upr = round(upr, 1), 
         p.adj = round(p.adj , 3), 
         ci = paste0("(", lwr, "-",upr, ")"), 
         p = if_else(p.adj <0.0001, true = as.character("<0.001"), false = as.character(p.adj))) %>% 

  select(diff, ci, p) %>% 
  rename("Difference Body Fat (%)" = diff, 
         "Confidence Interval" = ci) %>% 
  DT::datatable(rownames = T, options = list(columnDefs = list(list(className = 'dt-center', 
                                     targets = "_all"))))

poshocbmigender


```


### **Summary**

Based on the data collected from 2011 to 2018 I have answered the following questions: 

* What is the prevalence and distribution of obesity rates by gender, age and the year in which data was collected? 

From 2011 to 2018 about 2 in 3 American were considered overweight or living with obesity.

Based on BMI, there were more men who were considered overweight (58.4% vs 41.6%) and women tended to experience higher obesity class II-III rates versus men (61.1% vs 38.9%). Additionally, on average females tended to have higher body fat percentages across all age groups (*p* <0.001). 

* Is Body Mass Index a good indicator of body fat as measured by the gold standard method dual energy x-ray absorptiometry (DEXA)? 

I correlated BMI against body fat and found that there was a positive relationship between the two with *r* = 0.59, *p* <0.001. This suggests that BMI as an estimate of fat percentage is adequate but not perfect.  

* What are the differences between body fat composition between those in different weight classes?

It was found that body fat percentages are significantly different across all weight classes and genders. An anova table has been provided. The only interaction that was not significant, it was found that on average underweight females had roughly 1% higher body fat compared to overweight males.