datafest_data_cleaning_and_EDA.Rmd

---
title: "Data aggregation, cleaning, and EDA"
output:
  html_document:
    df_print: paged
  html_notebook: default
  pdf_document: default
date: "April 14th 2020"
---

```{r setup, include=FALSE}
library(fs) #file utils
library(lubridate) # date utils 
library(DT) # pretty interactive data tables
library(tidyverse)
```

## Datasets used: 

(Obviously can be adjusted)

Arizona: 
- PM 2.5 data for 2020 
- Carbon Monoxide data for 2020
- Covid cases by county

California:
- PM 2.5 data for 2020 
- Carbon Monoxide data for 2020
- Covid cases by county

### Initial import & cleaning (Note that many locations lack PM 2.5 *and* CO data): 

```{r datasets}
data_directory <- "/Users/elizabethhaderer/Documents/github/duke_datafest_covid19/datasets" # change appropriately

# Covid data from NYT: https://github.com/nytimes/covid-19-data 
covid_by_county_url <- "https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-counties.csv"
covid_by_county_dfs <- read_csv(url(covid_by_county_url)) %>% 
  filter(state %in% c("California", "Arizona")) %>% 
  mutate_if(is.character, toupper) %>% 
  rename_all(toupper) %>% 
  group_by(STATE) %>% # split by state to be able to combine w/EPA data more easily 
  group_split()

# Arizona EPA data
az_files <- fs::dir_ls(data_directory, regexp = "ad_viz_plotval_data_AZ*")
az_df <- az_files %>% 
  map(read_csv) %>% 
  map(~ .x %>% 
        mutate_if(is.character, toupper)) %>% 
  reduce(full_join, 
         by=c("Date","STATE_CODE", "STATE", "COUNTY_CODE", "COUNTY", "Site ID", "Site Name", "SITE_LATITUDE", "SITE_LONGITUDE", "CBSA_CODE", "CBSA_NAME"), 
         suffix = c("_co", "_pm25")) %>% 
  rename_all(toupper) %>% 
  mutate_at("DATE", function(x) as.Date(x, format = "%m/%d/%Y")) %>% 
  left_join(covid_by_county_dfs[[1]], by=c("DATE", "STATE", "COUNTY"))

# California EPA data 
ca_files <- fs::dir_ls(data_directory, regexp = "ad_viz_plotval_data_CA*")
ca_df <- ca_files %>% 
  map(read_csv) %>% 
  map(~ .x %>% 
        mutate_if(is.character, toupper)) %>% 
  reduce(full_join, 
         by=c("Date","STATE_CODE", "STATE", "COUNTY_CODE", "COUNTY", "Site ID", "Site Name", "SITE_LATITUDE", "SITE_LONGITUDE", "CBSA_CODE", "CBSA_NAME"), 
         suffix = c("_co", "_pm25")) %>% 
  rename_all(toupper) %>% 
  mutate_at("DATE", function(x) as.Date(x, format = "%m/%d/%Y")) %>% 
  left_join(covid_by_county_dfs[[2]], by=c("DATE", "STATE", "COUNTY"))
```

### Look @ data sets:

```{r}
# A simplified list of columns to look at. I don't know what a lot of them mean so this is just a guess. 
relevant_initial_cols <- c("DATE", "STATE", "COUNTY", "SITE NAME", "DAILY MAX 8-HOUR CO CONCENTRATION", "DAILY MEAN PM2.5 CONCENTRATION", "CASES", "DEATHS")

#Note. According to [this site](covid19.healthdata.org), California's stay at home order was March 19th, 2020. 
ca_stay_at_home_start <- as.Date("2020-03-19")
```

```{r}
az_df %>% 
  select(relevant_initial_cols) %>% 
  filter(DATE > ca_stay_at_home_start) %>% 
  datatable(filter="top", 
            caption="AZ data on CO, PM 2.5, & COVID Cases/Deaths")
```

```{r}
ca_df %>% 
  select(relevant_initial_cols) %>% 
  filter(DATE > ca_stay_at_home_start) %>% 
  datatable(filter="top", 
            caption="CA data on CO, PM 2.5, & COVID Cases/Deaths")
```



##Comparing Pre & Post Shut Down

California Shut Down Date: 03/19/2020
Pre Shut Down Date Subset:03/05/2020-3/19/2020
Post Shut Down Date Subset : 3/20/2020-4/4/2020

Assumptions: PM 2.5 measurements are independent, and each day's mean is a representative sample for that region so we will treat as a random sample.

Arizona

```{r az}
#Subset Pre & Post and calculate mean PM 2.5 for that date range
az_pre <- subset(az_df, DATE > "2020-02-18" & DATE < "2020-03-20")

az_pre_mean <- mean(az_pre$`DAILY MEAN PM2.5 CONCENTRATION`)

az_post <- subset(az_df, DATE> "2020-03-19" & DATE < "2020-04-21")

az_post_mean <- mean(az_post$`DAILY MEAN PM2.5 CONCENTRATION`)
#Conduct a t.test to see if there is a significant difference

var.test(az_pre$`DAILY MEAN PM2.5 CONCENTRATION`, az_post$`DAILY MEAN PM2.5 CONCENTRATION`)

t.test(az_pre$`DAILY MEAN PM2.5 CONCENTRATION`, az_post$`DAILY MEAN PM2.5 CONCENTRATION`, var.equal = FALSE)


```


California
```{r ca}
#Subset Pre & Post PM 2.5 for that date range
ca_pre <- subset(ca_df, DATE > "2020-02-18" & DATE < "2020-03-20")

ca_pre_mean <- mean(ca_pre$`DAILY MEAN PM2.5 CONCENTRATION`)

ca_post <- subset(ca_df, DATE> "2020-03-19" & DATE < "2020-04-21")

ca_post_mean <- mean(ca_post$`DAILY MEAN PM2.5 CONCENTRATION`)

#Conduct a t.test to see if there is a significant difference

var.test(ca_pre$`DAILY MEAN PM2.5 CONCENTRATION`, ca_post$`DAILY MEAN PM2.5 CONCENTRATION`)

t.test(ca_pre$`DAILY MEAN PM2.5 CONCENTRATION`, ca_post$`DAILY MEAN PM2.5 CONCENTRATION`, var.equal = FALSE)



```


### Animation of PM 2.5 values in CA and AZ 

```{r}
### Setup for animation ###

library(maps)
library(ggmap)
library(mapdata)
library(ggthemes)
library(gganimate)
library(gifski)
library(viridis)

# Add more steps to color scale 
custom_colors <- viridis(20)

combined_ca_az <- az_df %>% 
  full_join(ca_df) %>% 
  select(DATE,SITE_LONGITUDE, SITE_LATITUDE, `DAILY MEAN PM2.5 CONCENTRATION`, STATE) %>% 
  filter(`DAILY MEAN PM2.5 CONCENTRATION` < 50) # remove outlier points just from visualization because they skew the color scale too much 

# Define plot 
combined_map_plot <- ggplot() + 
  geom_polygon(data=az_and_ca, aes(x=long, y = lat, group = group), 
               fill="lightgrey", color="white") +
  coord_fixed(1.3) +
  geom_point(data= combined_ca_az, 
             aes(x=SITE_LONGITUDE, y=SITE_LATITUDE, 
                 size=40, 
                 fill=`DAILY MEAN PM2.5 CONCENTRATION`),
             colour="black",
             alpha=0.8,
             pch=21) + 
  scale_fill_gradientn(colours=custom_colors,
        name = "Daily mean PM 2.5 Concentration",
    guide = guide_colorbar(
      direction = "horizontal",
      barheight = unit(5, units = "mm"),
      barwidth = unit(50, units = "mm"),
      draw.ulim = F,
      title.position = 'top',
      # some shifting around
      title.hjust = 0.5,
      label.hjust = 0.5
  )) + 
  labs(title="PM 2.5 levels: {closest_state}") + 
  theme_map() + 
  guides(size=FALSE) + 
  theme(legend.position="bottom", 
        plot.title = element_text(size=27),
        legend.title = element_text(size=16),
        legend.text = element_text(size=16)) + 
  transition_states(DATE, transition_length=1, state_length =1) 

# animate plot by date 
animate(combined_map_plot, 
        nframes=300,
        renderer=gifski_renderer())
```

```{r}
# Save plot to gif 
anim_save("ca_and_az_animation.gif", 
          animation = last_animation())
```