explore_phys_activity_questions_solutions.Rmd

---
title: "NHANES BioData Exploratory Data Analysis - Physical Activity"
author: "JM"
date: '`r Sys.Date()`'
output:   
  github_document:
    toc: yes
---


```{r setup, include=FALSE}
library(ggplot2)
library(dplyr)
library(DT)
library(data.table)
library(naniar)
library(visdat)
library(skimr)
library(readr)
library(NHANES)
library(janitor)
library(tidyverse)
#### Global chunk options -----------------------------

knitr::opts_chunk$set(
  eval       = TRUE,    # whether to run code in code chunk
  include    = TRUE,    # whether to include the chunk output
  echo       = TRUE,   # Whether to show code chunk in final output
  error      = TRUE,    # whether to display error messages
  message    = FALSE,   # whether to preserve messages
  warning    = FALSE   # whether to preserve warnings
)
```

```{r data, include=FALSE}

source("R/helper.R")

#data_dictionary <- read.csv("data/shhs-data-dictionary-0.13.1-variables.csv")
#myDataFrame <- data.table(read_rds("data/common_data_small.rds"))
#load("data/bmi_diabetes.rda")
#myDataFrame <- data.table(bmi_diabetes)

data("NHANES")

##specify outcome variable here
outcome_var <- c("PhysActive")
## specify covariates here (including outcome variable)
covariates <- c("Gender", "Age", "SurveyYr", "Race1", "Race3" ,"MaritalStatus", 
                "BMI", "HHIncome", "Education",
                "BMI_WHO", "BPSysAve", "TotChol", "Depressed", "LittleInterest", 
                "Pulse", "Diabetes", "DiabetesAge",
                "PhysActive","PhysActiveDays","PhysActiveDaysAtLeast3",
                "SleepHrsNight", "SleepTrouble", "SleepHrsNightCat","TVHrsDay", "AlcoholDay", "SmokeNow","Marijuana")

NHANES <- NHANES %>% mutate(
  PhysActiveDaysAtLeast3=factor(1*(PhysActiveDays>=3),levels=c(0,1),labels=c("No","Yes")),
  SleepHrsNightCat=case_when(SleepHrsNight<6 ~ "<6",
                             dplyr::between(SleepHrsNight,6,9) ~ "6-9",
                             SleepHrsNight>9 ~ ">9",
                             TRUE ~ as.character(NA)))

myDataFrame <- data.table(NHANES)[,covariates,with=FALSE]

remove_categories <- outcome_var
categoricalVars <- sort(names(get_category_variables(myDataFrame)))
cat_no_outcome <- setdiff(categoricalVars, remove_categories)

remove_numeric <- c("nsrrid")
numericVars <- sort(get_numeric_variables(myDataFrame))
numericVars <- setdiff(numericVars, remove_numeric)

theme_set(theme_classic(base_size = 15))
data_dictionary <- readr::read_csv("data/data_dictionary.csv")

data_dictionary <- data_dictionary %>%
  add_row(VariableName = "PhysActiveDaysAtLeast3",
          Definition = "PhysActiveDays>=3 ~ Yes, PhysActiveDays < 3 ~ No") %>%
  add_row(VariableName = "SleepHrsNightCat",
          Definition = "SleepHrsNight categorized into <6hrs, [6-9]hrs, >9hrs") %>%
  arrange(VariableName)

```


This is a report showing possible answers posed in the lesson plan on [github](https://github.com/laderast/nhanes_explore/blob/master/lesson_plan.md).

# Data

Variables: 

```{r}
names(myDataFrame)
```

# What is the percent missing for PhysActive (yes/no physically active)?

```{r}
myDataFrame %>% tabyl(PhysActive) %>% adorn_pct_formatting() %>% kable
```

# What are the characteristics of people who have missing responses for PhysActive?


- tend to be missing in other lifestyle questionnaires (smoking, education, etc)
- seem to be mostly children

```{r}
myDataFrame %>% filter(is.na(PhysActive)) %>% skim()

myDataFrame %>%
      data.frame() %>%
      gg_miss_fct(fct = PhysActive) + 
      theme(axis.text = element_text(size = 15))

ggplot(myDataFrame, aes(x=PhysActive, y=Age, fill=PhysActive)) + 
      geom_boxplot() + theme(text=element_text(size=20), axis.text.x = element_text(angle=90))

```

# How does missingness and responses in PhysActive relate to responses and missingness in PhysActiveDays and PhysActiveDaysAtLeast3? 

How does the proportion of missingness in PhysActiveDaysAtLeast3 (yes = PhysActiveDays>= 3, no < 3) relate to PhysActive response? Is this what you would expect?

We would expect if `PhysActive`="No" then `PhysActiveDaysAtLeast3` would be either "No" always, or missing. However, there is 77% missingness and 16.7% "Yes" Similarly, when `PhysActive` is NA, there are still ~25% responses to `PhysActiveDaysAtLeast3`.

I would use `PhysActiveDays` and `PhysActiveDaysAtLeast3` very cautiously in analyses!

```{r}
data_dictionary %>% filter(str_detect(VariableName,"PhysAct")) %>% kable
```

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(PhysActive) %>%
      count(PhysActiveDaysAtLeast3) %>% mutate(ratio=scales::percent(n/sum(n)))
    
myDataFrame %>% 
      ggplot(aes(x=PhysActive, fill=PhysActiveDaysAtLeast3)) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))
```

Also, if we look at median `PhysActiveDays` between people who are physically active and who are not, it is roughly the same. Huh??

```{r}
ggplot(myDataFrame, aes(x=PhysActive, y=PhysActiveDays, fill=PhysActive)) + 
      geom_boxplot() + theme(text=element_text(size=20), axis.text.x = element_text(angle=90))

```

A more granular look at how many physically active days people reported (note how we can see the differences here, but not with using PhysActiveDays as a continuous variable and looking at the median):

```{r}

percent_table <- myDataFrame %>% data.frame() %>% group_by(PhysActive) %>%
      count(PhysActiveDays) %>% mutate(ratio=scales::percent(n/sum(n)))

myDataFrame %>% 
      ggplot(aes(x=PhysActive, fill=factor(PhysActiveDays))) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))
```


# Is there something strange about Age?

Age is truncated at 80 for privacy reasons.

```{r}
ggplot(myDataFrame, aes(x=Age, y=BMI)) + geom_point()
```

# What is the median and interquartile range of PhysActiveDays?

- This can be found in the skim output: Overview -> Tabular Summary of Data
- Median 3, IQR [2,5]

```{r}
myDataFrame %>% select(contains("PhysAct")) %>% skim
```

# Is number of sleep hours associated with physical activity? Consider using the categorized variable `SleepHrsNightCat`.

If we compare median number of hours slept between physical activity categories, we see no difference. However, if we look at sleep as a categorical variable, we see that people with "normal" sleep (6-9 hours sa night) have a slightly higher proportion of physical activity. However, this is a risky analysis since the result really depends on how we categorize our sleep variable.

```{r}
ggplot(myDataFrame, aes(x=PhysActive, y=SleepHrsNight, fill=PhysActive)) + 
      geom_boxplot() + theme(text=element_text(size=20), axis.text.x = element_text(angle=90))

ggplot(myDataFrame, aes(x=PhysActiveDaysAtLeast3, y=SleepHrsNight, fill=PhysActive)) + 
      geom_boxplot() + theme(text=element_text(size=20), axis.text.x = element_text(angle=90))
```


```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(SleepHrsNightCat) %>%
      count(PhysActive) %>% mutate(ratio=scales::percent(n/sum(n)))
    
myDataFrame %>% 
      ggplot(aes(x=SleepHrsNightCat, fill=PhysActive)) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))

```

If we use `PhysActiveDaysAtLeast3` we don't see as strong of an association:

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(SleepHrsNightCat) %>%
      count(PhysActiveDaysAtLeast3) %>% mutate(ratio=scales::percent(n/sum(n)))

myDataFrame %>%
      ggplot(aes(x=SleepHrsNightCat, fill=PhysActiveDaysAtLeast3)) +
      geom_bar(position="fill", color="black") +
  theme(text=element_text(size=20),
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio),
                position=position_fill(vjust=0.5))
```

# What about levels of education?


Higher levels of education are associated with increase proportions of physical activity (yes/no) as well as higher proportions of physical activity days >=3.

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(Education) %>%
      count(PhysActive) %>% mutate(ratio=scales::percent(n/sum(n)))

myDataFrame %>%
      ggplot(aes(x=Education, fill=PhysActive)) +
      geom_bar(position="fill", color="black") +
  theme(text=element_text(size=20),
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio),
                position=position_fill(vjust=0.5))
```

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(Education) %>%
      count(PhysActiveDaysAtLeast3) %>% mutate(ratio=scales::percent(n/sum(n)))

myDataFrame %>%
      ggplot(aes(x=Education, fill=PhysActiveDaysAtLeast3)) +
      geom_bar(position="fill", color="black") +
  theme(text=element_text(size=20),
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio),
                position=position_fill(vjust=0.5))
```

# What about BMI? Is the association of BMI with physical activity dependent on which measure of activity you use -- PhysActive (yes/no) or PhysActiveDaysAtLeast3?

There seems to be a small difference in BMI comparing people in physical activity yes/no groups, but when comparing physically active >=3 days yes or no, the difference is nonexistent. Either way, the differences are actually quite small so it is unlikely to be statistically significant.

```{r}
ggplot(myDataFrame, aes(x=PhysActive, y=BMI, fill=PhysActive)) + 
      geom_boxplot() + theme(text=element_text(size=20), axis.text.x = element_text(angle=90))

ggplot(myDataFrame, aes(x=PhysActiveDaysAtLeast3, y=BMI, fill=PhysActiveDaysAtLeast3)) + 
      geom_boxplot() + theme(text=element_text(size=20), axis.text.x = element_text(angle=90))

ggplot(myDataFrame, aes(x=factor(PhysActiveDays), y=BMI, fill=factor(PhysActiveDays))) + 
      geom_boxplot() + theme(text=element_text(size=20), axis.text.x = element_text(angle=90))
```

# Are there other covariates you expect to be associated with physical activity? How do they relate to each other?

We've explored hours of sleep, BMI, and education. There are other  lifestyle variables such as hours of TV watched, depression, smoking, Marijuana use, that might be related to physical activity. These are likely associated with each other through a myriad of causal effect pathways.

First, we can note how education and number of hours of sleep per night are also associated, so education may be a confounder in the association of hours per sleep and physical activity.

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(Education) %>%
      count(SleepHrsNightCat) %>% mutate(ratio=scales::percent(n/sum(n)))
    
myDataFrame %>% 
      ggplot(aes(x=Education, fill=SleepHrsNightCat)) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))
```

Physical activity is also associated with several other health outcomes including depression and diabetes:

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(PhysActive) %>%
      count(Depressed) %>% mutate(ratio=scales::percent(n/sum(n)))
    
myDataFrame %>% 
      ggplot(aes(fill=Depressed, x=PhysActive)) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))
```

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(PhysActive) %>%
      count(Diabetes) %>% mutate(ratio=scales::percent(n/sum(n)))
    
myDataFrame %>% 
      ggplot(aes(x=PhysActive, fill=Diabetes)) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))
```

Risk factors such as smoking and TV hours per day are also related. People who smoke are less likely to be physically active.

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(SmokeNow) %>%
      count(PhysActive) %>% mutate(ratio=scales::percent(n/sum(n)))
    
myDataFrame %>% 
      ggplot(aes(x=SmokeNow, fill=PhysActive)) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))
```

The relationship between number of TV hours and physical activity is somewhat strange, since those who said 0 hours of TV have a large proportion of not physically active, but 0-1 hours has a much lower proportion of non-physically active and that proportion increases with the number of hours of TV watched. Perhaps people who say no TV are somehow different than people who admit to watching a small amount of TV.

```{r}
percent_table <- myDataFrame %>% data.frame() %>% group_by(TVHrsDay) %>%
      count(PhysActive) %>% mutate(ratio=scales::percent(n/sum(n)))
    
myDataFrame %>% 
      ggplot(aes(x=TVHrsDay, fill=PhysActive)) + 
      geom_bar(position="fill", color="black") + 
  theme(text=element_text(size=20), 
        axis.text.x = element_text(angle = 90)) +
      geom_text(data = percent_table, mapping = aes(y=n, label=ratio), 
                position=position_fill(vjust=0.5))
```


# Data Dictionary

```{r}
data_dictionary %>% filter(VariableName %in% covariates) %>% kable
```