09-start_to_finish.Rmd

# From start to finish {#start-finish}


> Let's pull all these steps together into a few functions.


The script below shows an example workflow that starts with 24-hour air toxics data from MPCA's Air Toxics database, walks through each of the analysis steps outlined in the previous sections, and produces a suite of data validation tools, POC analysis tools, completeness checks, annual summaries, site comparisons, and pollution roses.


## Using functions  

The first step is to put our analysis methods into a few reusable functions, which can then be re-used for all of our monitoring conditions. 


_Click below to view the example functions._

<div class="toggle">
<button class="btn_code">Show __R__ code</button>

<br>

```{r, eval=F}
# Functions

install.MPCAair.packages = function() {

install.packages(c(
  "easypackages",
  "tidyverse",
  "data.table",
  "nortest",
  "car",
  "DT",
  "RcppRoll",
  "shiny",
  "rsconnect",
  "EnvStats",
  "openair",
  "ggbeeswarm",
  "reshape",
  "corrplot",
  "RODBC"
))

}

MPCA_air_libraries = function() {
  
library(easypackages)
libraries(c(
  "tidyverse",
  "lubridate",
  "stats",
  "nortest",
  "car",
  "DT",
  "RcppRoll",
  "shiny",
  "rsconnect",
  "EnvStats",
  "openair",
  "ggbeeswarm",
  "reshape",
  "corrplot",
  "RODBC"
))

}

flag <- function(data) {
  
  library(tidyverse)
  library(lubridate)
  
  data <- mutate(data,
                AQSID_flag = is.na(AQSID),
                POC_flag = is.na(as.numeric(POC) ) | as.numeric(POC) < 0 | as.numeric(POC) > 4,
                Parameter_flag = is.na(as.numeric(Parameter)) | nchar(as.character(Parameter)) != 5,
                Date_flag = is.na(ymd(Date)),
                Result_flag = !is.na(Result) & (abs(as.numeric(Result) ) >= 900 | is.na(as.numeric(Result) ) ),
                Null_flag = (is.na(Result) & is.na(Null_Data_Code) ) | (!is.na(Result) & !is.na(Null_Data_Code) ),
                MDL_flag = is.character(MDL) | MDL < 0,
                Pollutant_flag = is.na(Pollutant),
                any_flag = (AQSID_flag + POC_flag + Parameter_flag + Date_flag + Result_flag + Null_flag +
                              MDL_flag + Pollutant_flag) > 0
                )
}

remove_flagged <- function(data) {
  
  library(tidyverse)
  
  return(filter(data, !any_flag) %>% 
           select(-contains("flag") ) %>%
           mutate(Result = as.numeric(Result) ) )
}

flag_duplicates = function(data) {
  library(tidyverse)
  return(data %>% group_by(AQSID, POC, Parameter, Date) %>% mutate(Count = n(), duplicate_flag = Count > 1))
}

average_duplicates = function(data) {
  library(tidyverse)
  
  dupe_averaging = function(Result, Censored) {
    
    Result = Result[!is.na(Result)]
    Censored = Censored[!is.na(Censored)]
    
    if(all(Censored, na.rm = T)) {
      return (mean(Result, na.rm = T))
    }
    else {
      return (mean(Result[!Censored], na.rm = T ) )
    }
    
  }
  
  data = data %>% group_by(AQSID, POC, Parameter, Pollutant, Date, MDL) %>% mutate(Result = dupe_averaging(Result, Censored), Censored = all(Censored, na.rm = T) ) %>% slice(1) %>% ungroup() %>%
    mutate(Result = ifelse(is.na(Result), NA, Result), Censored = ifelse(is.na(Result), NA, Censored))
  return (data)
}



time_series_plots = function(data) {
  
  library(shiny)
  library(tidyverse)
  library(RcppRoll)
  library(DT)
  library(rsconnect)


data <- mutate(data, sitePOC = paste0(AQSID,"-", POC) )
pollutant <- unique(data$Pollutant)
site <- unique(data$sitePOC)


shinyApp(
  ui = fluidPage(
                 fluidRow(
                   column(3,
                          style = "padding-bottom: 20px;",
                          inputPanel(
                            selectInput("pollutant", label="Choose a pollutant", choices = pollutant, 
                                        selected="Benzene"),
                            selectInput("site", label="Choose a site", choices = site, selected=270535501),
                            dateRangeInput("date", label = "Select date range", start = "2009-01-01", end = 
                                             "2013-12-31", min = "2009-01-01", max = "2013-12-31") ) ),
                   column(9,
                          plotOutput('detlim', height = "400px")))),
  
  
  server = function(input, output) {
    
    
    output$detlim <- renderPlot({
      data_sub = filter(data, Pollutant==input$pollutant, sitePOC == input$site, Date >= input$date[1], Date <= input$date[2])
      ggplot(data=data_sub, aes(x= Date, y=Result)) +
        geom_point(aes(color=Censored), size =3, alpha=0.55) +
        geom_line() +
        scale_x_date(date_labels = "%D") +
        xlab(NULL) +
        ylab("Result (ug/m3)") +
        expand_limits(y=c(0, max(data_sub$Result))) +
        scale_colour_manual(values= c("#197519"[FALSE %in% unique(data_sub$Censored)], "#0000FF"[TRUE %in% unique(data_sub$Censored)]), breaks=c(FALSE, TRUE)) +
        theme(text = element_text(size=15), axis.text.x = element_text(angle = -90, vjust = 0.3,  size=14)) +
        ggtitle(paste0("Time series for ", input$pollutant, " at site ", input$site))
    })
    
  })

}




POC_compare = function(data) {

library(shiny)
library(tidyverse)
library(rsconnect)  
  
data <- distinct(data, AQSID, POC, Parameter, Date, Pollutant, Year, .keep_all = T) #replace with better cleaning function

data <- spread(data, POC, Result) %>% 
        mutate(Status = ifelse(`1` < MDL  & `2` < MDL, "POCs 1 and 2 below MDL", ifelse(`1` < MDL, "POC 1 below MDL", ifelse(`2` < MDL, "POC 2 below MDL", "POCs 1 and 2 above MDL") ) ) ) %>% 
       drop_na(Status)


Pollutant <- unique(data$Pollutant)

Site <- unique(data$AQSID)

shinyApp(
  ui = fluidPage(
                 fluidRow(
                   column(3,
                          style = "padding-bottom: 20px;",
                          inputPanel(
                            selectInput("Pollutant", label="Choose a pollutant", choices = Pollutant),
                            selectInput("Site", label="Choose a site", choices = Site),
                            dateRangeInput("date", label = "Select date range", start = "2009-01-01", end = 
                                             "2013-12-31", min = "2009-01-01", max = "2013-12-31"))),
                   column(9,
                          plotOutput('normviz', height = "500px")))),
  
  
  server = function(input, output) {
    
    
    
    
    output$normviz <- renderPlot({
      data_sub = filter(data, Pollutant==input$Pollutant, AQSID == input$Site, Date >= input$date[1], Date <= input$date[2])
      ggplot(data_sub, aes(x = `1`, y = `2`, color = Status)) +
        geom_point(size = 3) +
        geom_segment(x=-1000, xend=1000, y=-1000, yend=1000, color="red", size=1) +
        labs(title = "POC comparison chart", x = "POC 1", y = "POC 2", subtitle = paste("Correlation =", round(cor(data_sub$`1`,data_sub$`2`, use = "complete"), 2) ) )
    })
    
  })

}

POC_average = function(data) {
  library(tidyverse)
  
  POC_averaging = function(Result, Censored) {
    
    Result = Result[!is.na(Result)]
    Censored = Censored[!is.na(Censored)]
    
    if(all(Censored, na.rm = T)) {
      return (mean(Result, na.rm = T))
    }
    else {
      return (mean(Result[!Censored], na.rm = T ) )
    }
    
  }
  
  data = data %>% group_by(AQSID, Parameter, Pollutant, Date, MDL) %>% mutate(Result = POC_averaging(Result, Censored), Censored = all(Censored, na.rm = T) ) %>% slice(1) %>% select(-POC) %>% ungroup() %>%
    mutate(Result = ifelse(is.na(Result), NA, Result), Censored = ifelse(is.na(Result), NA, Censored))
  return (data)
}


completeness_check = function(data) {

# Create a sampling calendar based on EPA's air toxics monitoring schedule
sample_calendar <- function(start         = "2012-01-01", 
                            end           = "2016-12-31", 
                            day_interval  = 6,
                            type          = "air_toxics") {
  
  library(lubridate)
  
  # Convert 'start' and 'end' to class date
  start <- ymd(start)
  end   <- ymd(end)
  
  # Set official start date to selected EPA calendar
  if(type == "air_toxics") {
      epa_start <- ymd("1989-12-24")
  } else {
      epa_start <- start
  }
  
  # Create full table of sampling dates
  calendar <- seq(from = epa_start, 
                  to   = end, 
                  by   = paste(day_interval, "days"))
  
  
  # Subset to user's date range
  calendar <- calendar[calendar >= start & calendar <= end]
  
  return(calendar)
  
}


# Find the year range of your data
date_range <- range(data$Date)

# Create expected sample calendar
epa_schedule <- tibble(Date = sample_calendar(start = format(date_range[1], "%Y-01-01"), #Extend range to first day of the year
                                                  end   = format(date_range[2], "%Y-12-31"), #Extend range to last day of the year
                                                  day_interval = 6))

# Add year and calendar quarter columns
epa_schedule <- epa_schedule %>% mutate(Year        = year(Date),
                                        cal_quarter = quarter(Date))


# Count the expected number of samples per quarter and year.
epa_schedule <- epa_schedule %>% 
                group_by(Year, cal_quarter) %>%
                summarize(expected_quarter_samples = length(unique(Date))) %>%
                group_by(Year) %>%
                mutate(expected_annual_samples = sum(expected_quarter_samples))


# Assign each date to a calendar quarter
data <- data %>% mutate(cal_quarter = quarter(Date))

# Count the number of sampling dates for each quarter and year.
data <- data %>% 
        group_by(AQSID, Parameter, Pollutant, Year, cal_quarter) %>%
        mutate(valid_quarter_samples = length(unique(Date[!is.na(Result)]))) %>%
        group_by(AQSID, Parameter, Pollutant, Year) %>%
        mutate(valid_annual_samples = length(unique(Date[!is.na(Result)])))


# Join expected sample table to data by quarter and year columns
data <- left_join(data, epa_schedule, by = c("Year", "cal_quarter"))


# Divide valid samples by expected samples
data <- data %>% 
        group_by(AQSID, Parameter, Pollutant, Year, cal_quarter) %>%
        summarise(pct_quarter_samples = round(valid_quarter_samples[1] / expected_quarter_samples[1], 2)) %>%
        mutate(Complete = pct_quarter_samples >= 0.75) %>%
        group_by(AQSID, Parameter, Pollutant, Year) %>% 
        summarise(Complete = sum(Complete, na.rm = T) == 4,
                  lowest_quarter = ifelse(n() == 4, min(pct_quarter_samples, na.rm = T), 0) )

return(data %>% select(AQSID, Parameter, Pollutant, Year, Complete, lowest_quarter) %>% ungroup() )

}
        



UCL_95 = function(data, Boot_Repeats = 1000) {
  library(EnvStats)
  set.seed(2017)
  
  annual_AT_means = function(air_toxics) {
  
  air_toxics = mutate(air_toxics, Year = year(ymd(Date)), Quarter = quarter(ymd(Date)) )
  sample_complete = air_toxics %>% completeness_check()
  
  enough_detects = air_toxics %>% group_by(AQSID, Parameter, Pollutant, Year) %>% summarise(Detected = mean(Censored, na.rm = T) <= 0.8 )
  
  site_means = air_toxics %>% group_by(AQSID, Parameter, Pollutant, Year) %>% summarise(Mean = ifelse(length(unique(Result[!is.na(Result) & !Censored] ) ) < 2, NA,
     ifelse (any(Censored, na.rm = T), elnormAltCensored(Result, Censored, method = "impute.w.mle", ci = F)$parameters[[1]], mean(Result, na.rm = T) ) ) )
  
  site_means = left_join(site_means, sample_complete, by = c("AQSID", "Parameter", "Pollutant", "Year") ) %>%
    left_join(enough_detects, by = c("AQSID", "Parameter", "Pollutant", "Year") ) %>% mutate(Mean = ifelse(Complete & Detected, Mean, NA), ID = paste(AQSID, Parameter, Pollutant, Year) )

  return(site_means)
}
  
  MLE_est <- function(data){
    
    results = data$Result
    censored = data$Censored
    n = sum(!is.na(results))
    
    if (length(unique(results[!is.na(results) & !censored] ) ) < 2 ) {
      MLE_means = NA
    }
    
    else {
      random.rows = NULL  
      random.rows = sample(which(!is.na(censored) & (!censored) & !duplicated(results) ), 2, replace = FALSE)
      random.rows = c(random.rows, sample(which(!is.na(censored)), n-2, replace = TRUE))
      
      MLE_means = ifelse(sum(censored[random.rows], na.rm = T) == 0, mean(results[random.rows]), elnormAltCensored(results[random.rows], censored[random.rows], method = "impute.w.mle", ci = F)$parameters[[1]] )
    }
    
    return(MLE_means)
  }
  
  data = mutate(data, ID = paste(AQSID, Parameter, Pollutant, Year), Result = ifelse(Censored, MDL, Result) )
  
  Bootstrap_means = replicate(Boot_Repeats, (by(data, data$ID, MLE_est) ) )
  
  CL = apply(Bootstrap_means, 1, function(x) sort(x)[ceiling(0.95 * Boot_Repeats)] )
  CL = data.frame(ID = names(CL), UCL95 = unname(CL))
  
  annual_summary = left_join(annual_AT_means(data), CL, by = "ID") %>% mutate(UCL95 = ifelse(Complete & Detected, UCL95, NA) ) %>% select(-ID)
  
  return(annual_summary)
  
}

site_compare = function(data, site_number, Boot_Repeats = 1000) {
  library(EnvStats)
  set.seed(2017)
  
  annual_AT_means = function(air_toxics) {
  
  air_toxics = mutate(air_toxics, Year = year(ymd(Date)), Quarter = quarter(ymd(Date)) )
  sample_complete = air_toxics %>% completeness_check()
  
  enough_detects = air_toxics %>% group_by(AQSID, Parameter, Pollutant, Year) %>% summarise(Detected = mean(Censored, na.rm = T) <= 0.8 )
  
  site_means = air_toxics %>% group_by(AQSID, Parameter, Pollutant, Year) %>% summarise(Mean = ifelse(length(unique(Result[!is.na(Result) & !Censored] ) ) < 2, NA,
     ifelse (any(Censored, na.rm = T), elnormAltCensored(Result, Censored, method = "impute.w.mle", ci = F)$parameters[[1]], mean(Result, na.rm = T) ) ) )
  
  site_means = left_join(site_means, sample_complete, by = c("AQSID", "Parameter", "Pollutant", "Year") ) %>%
    left_join(enough_detects, by = c("AQSID", "Parameter", "Pollutant", "Year") ) %>% mutate(Mean = ifelse(Complete & Detected, Mean, NA), ID = paste(AQSID, Parameter, Pollutant, Year) )

  return(site_means)
}
  
  MLE_est <- function(data){
    
    results = data$Result
    censored = data$Censored
    n = sum(!is.na(results))
    
    if (length(unique(results[!is.na(results) & !censored] ) ) < 2 ) {
      MLE_means = NA
    }
    
    else {
      random.rows = NULL  
      random.rows = sample(which(!is.na(censored) & (!censored) & !duplicated(results) ), 2, replace = FALSE)
      random.rows = c(random.rows, sample(which(!is.na(censored)), n-2, replace = TRUE))
      
      MLE_means = ifelse(sum(censored[random.rows], na.rm = T) == 0, mean(results[random.rows]), elnormAltCensored(results[random.rows], censored[random.rows], method = "impute.w.mle", ci = F)$parameters[[1]] )
    }
    
    return(MLE_means)
  }
  
  data = mutate(data, Result = ifelse(Censored, MDL, Result), ID = paste(AQSID, Parameter, Pollutant, Year))
  
  Bootstrap_means = replicate(Boot_Repeats, (by(data, data$ID, MLE_est) ) )
  
  Bootstrap_means = rownames_to_column(as.data.frame(Bootstrap_means), "ID" )
  
  Bootstrap_means = right_join(annual_AT_means(data), Bootstrap_means, by = "ID")
  
  Bootstrap_means = Bootstrap_means %>% group_by(Parameter, Pollutant, Year) %>% arrange(desc(AQSID == site_number), .by_group = T ) %>%
    group_by(Parameter, Pollutant, Year) %>% mutate_at(vars(num_range("V", 1:Boot_Repeats)), funs(c(first(.), (. - first(.))[-1])) ) %>% ungroup()
  
  LB = select(Bootstrap_means, num_range("V", 1:Boot_Repeats) ) %>% apply(1, function(x) sort(-x)[floor(0.025 * Boot_Repeats)] )
  UB = select(Bootstrap_means, num_range("V", 1:Boot_Repeats) ) %>% apply(1, function(x) sort(-x)[ceiling(0.975 * Boot_Repeats)] )
  
  CI = data.frame(Lower = LB, Upper = UB)
  CI = bind_cols(CI, Bootstrap_means) %>% select(Lower:ID) %>% group_by(Pollutant, Year) %>%
    mutate(Lower = ifelse(any(AQSID == site_number & Complete & Detected) & AQSID != site_number & Complete & 
        Detected, Lower, NA), Upper = ifelse(any(AQSID == site_number & Complete & Detected) & AQSID != 
        site_number & Complete & Detected, Upper, NA), Comparison = ifelse(Lower > 0, "Higher", ifelse(Upper < 0, 
        "Lower", "Same") ) )
  
  return(CI %>% ungroup() )
  
}

correlation_plots = function(data, site) {
  library(tidyverse)
  library(corrplot)
  
  data_site <- filter(data, AQSID %in% site) %>% select(Date, Pollutant, Result)
  analytes <- spread(data_site, Pollutant, Result, drop=T)
  analytes$Date <- NULL

  coranalytes <- cor(analytes, method="kendall", use="pairwise.complete.obs") %>% as.data.frame()
  coranalytes <- select_if(coranalytes, function(x) !all(is.na(x))) %>% filter_all(any_vars(!is.na(.))) %>% 
    as.matrix()
  rownames(coranalytes) <- colnames(coranalytes)

  return(corrplot(coranalytes, method = "circle", type="lower", tl.cex=0.6) )#plot matrix
}




correlation_plots = function(data) {
  library(tidyverse)
  library(corrplot)
  
  site <- unique(data$AQSID)

  shinyApp(
  ui = fluidPage(
                 fluidRow(
                   column(3,
                          style = "padding-bottom: 20px;",
                          inputPanel(
                            selectInput("site", label="Choose a site", choices = site),
                            dateRangeInput("date", label = "Select date range", start = "2009-01-01", end = 
                                             "2013-12-31", min = "2009-01-01", max = "2013-12-31"))),
                            column(9,
                            plotOutput('normviz', height = "500px")))),
  
  
  server = function(input, output) {
    
    
    
    
    output$normviz <- renderPlot({
      data_sub = filter(data, AQSID==input$site, Date >= input$date[1], Date <= 
                          input$date[2]) %>% select(Date, Pollutant, Result)
      analytes <- spread(data_sub, Pollutant, Result, drop=T)
      analytes$Date <- NULL

      coranalytes <- cor(analytes, method="kendall", use="pairwise.complete.obs") %>% as.data.frame()
      coranalytes <- select_if(coranalytes, function(x) !all(is.na(x))) %>% filter_all(any_vars(!is.na(.))) %>% 
      as.matrix()
      rownames(coranalytes) <- colnames(coranalytes)
      dummy_obj = corrplot(coranalytes, method = "circle", type="lower", tl.cex=0.6)
      
    })
    
  })
  
}



beeswarm_plot = function(data) {
  library(dplyr)
  library(ggbeeswarm)
  library(ggplot2)
  
  data = mutate(data, AQSID = as.character(AQSID) )
    
  Pollutant <- unique(data$Pollutant)
  Year <- unique(data$Year)

  shinyApp(
  ui = fluidPage(
                 fluidRow(
                   column(3,
                          style = "padding-bottom: 20px;",
                          inputPanel(
                            selectInput("Pollutant", label="Choose a pollutant", choices = Pollutant),
                             dateRangeInput("date", label = "Select date range", start = "2009-01-01", end = 
                                             "2013-12-31", min = "2009-01-01", max = "2013-12-31") ) ),
                            column(9,
                            plotOutput('normviz', height = "500px")))),
  
  
  server = function(input, output) {
    
    
    
    
    output$normviz <- renderPlot({
      data_sub = filter(data, Pollutant==input$Pollutant, Date >= input$date[1], Date <= 
                          input$date[2], !is.na(Result) )
      ggplot(data_sub, aes(y = AQSID, x = Result, color = Censored) ) + 
      geom_quasirandom(groupOnX=F) +   
      labs(title = paste(data_sub$Pollutant[1]), x = "Result (ug/m^3)" )
      
    })
    
  })
  
}




read_met_data_tableau <- function(years = 2009:2017, stations = "MSP") {
  
  library(tidyverse)
  
  years <- paste0(years, collapse = ",")
  
  stations <- paste0(stations, collapse = ",")
  
  url <- paste0("http://tableau.pca.state.mn.us/views/WeatherObservations2009-2017/HourTable.csv?Station=",
               stations, "&Year=", years, "&Month=1,2,3,4,5,6,7,8,9,10,11,12")
  
  met_data <- read_csv(url, col_types = "??????????-") %>% 
              select(Station, Year, Month, Day, Hour, everything()) %>% 
              arrange(Station, Year, Month, Day, Hour)
  
  return(met_data)
}

pollution_roses = function(data, met_data, num_breaks = 5) { 
  # Met data must be in Tableau format
  
  library(tidyverse)
  library(openair)
  library(reshape)
  library(shiny)
  library(rsconnect)
  
  data$Date <- ymd(data$Date)
  
  met_data <- dplyr::rename(met_data, wd = `Wind Dir`, ws = `Wind Spd MPH`) %>%
              mutate(date = paste0(Year,"/",Month,"/",Day," ",Hour,":00"), date = ymd_hm(date)) %>%
              select(-Day, -Month, -Hour, -Year) %>%
              timeAverage(avg.time = "day") %>%
              mutate(date = ymd(date))
  
  data <- left_join(data, met_data, by = c("Date" = "date"))
  
  Pollutant <- unique(data$Pollutant)
  
  Site <- unique(data$AQSID)

  shinyApp(
  ui = fluidPage(
                 fluidRow(
                   column(3,
                          style = "padding-bottom: 20px;",
                          inputPanel(
                            selectInput("Pollutant", label="Choose a pollutant", choices = Pollutant),
                            selectInput("Site", label="Choose a site", choices = Site),
                            dateRangeInput("date", label = "Select date range", start = "2009-01-01", end = 
                                             "2013-12-31", min = "2009-01-01", max = "2013-12-31"))),
                   column(9,
                          plotOutput('normviz', height = "500px")))),
  
  
  server = function(input, output) {
    
    
    
    
    output$normviz <- renderPlot({
      data_sub = filter(data, Pollutant==input$Pollutant, AQSID == input$Site, Date >= input$date[1], Date <= 
                          input$date[2], !is.na(Result))
      data_sub = data_sub %>% mutate(MDL = max(MDL), minimum = min(Result), maximum = max(Result), Result = 
                                       ifelse(Censored, 1e-16, Result))
      breaks_site = NULL
      if(!all(data_sub$Censored)){
        
        breaks_site = c(breaks_site, 0,
                      round_any( c(data_sub$MDL[1], data_sub$MDL[1] + (data_sub$maximum[1] - data_sub$MDL[1]) * 
                                     (1:(num_breaks-1) / (num_breaks-1) ) ), 0.0001, ceiling ) )
        pollutionRose(data_sub, statistic = "abs.count", pollutant = "Result", breaks = breaks_site, 
                  key.footer="ug/m3", main=paste("Daily Average Pollution Rose for",
                  data_sub$Pollutant[1],"\n", data_sub$AQSID[1]) )
      }
      else {
        breaks_site = c(breaks_site, c(0, round_any( c(data_sub$MDL[1], 2*data_sub$MDL[1] ), 0.0001, ceiling ) ) )
        
        pollutionRose(data_sub, statistic = "abs.count", pollutant = "Result", breaks = breaks_site, 
                  key.footer="ug/m3", main=paste("Daily Average Pollution Rose for",
                  data_sub$Pollutant[1],"\n", data_sub$AQSID[1]) )
        
      }
      
    })
    
  })

}

```

</div>

<br>


## A simpler analysis 

Now we use our functions above to write a simpler and easier to read analysis script.


_Click below to view the example._

<div class="toggle">
<button class="btn_code">Show __R__ code</button>

```{r, eval=F, message=F, warning=F}

library(tidyverse)
MPCA_air_libraries()

#Import data
data <- read_csv('https://raw.githubusercontent.com/MPCA-air/air-methods/master/airtoxics_data_2009_2013.csv')

names(data)[1:10] <- c("AQSID", "POC", "Parameter", "Date","Result",
                       "Null_Data_Code", "MDL", "Pollutant", "Year", "CAS")

# Data cleaning

flagged <- data %>% flag() #Check with QA about flagged values if necessary.

duplicates_flagged <- data %>% flag_duplicates()

data <- data %>% mutate(Censored = Result < MDL) %>% average_duplicates() #We decide to average the duplicates here

# Data Validation
time_series <- data %>% time_series_plots()


# Collocated monitors
poc_comparisons <- data %>% POC_compare()

data <- data %>% POC_average()


# Completeness checks
complete <- data %>% completeness_check()


# Summary statistics
annual_summary <- data %>% filter(AQSID == 270370020) %>% UCL_95(100)


# Site Comparisons
site_number <- 270370020 #Our favorite site

comparisons_to_FH <- data %>% 
                    filter(AQSID %in% c(270370020, 270370470, 271230871) ) %>%
                    site_compare(site_number, 50)

pollutant_correlations <- data %>% correlation_plots()

beeswarms <- data %>% beeswarm_plot()


# Pollution Roses

poll_rose_data <- data

met_data <- read_met_data_tableau(2009:2013, "MSP")

num_breaks <- 5

pol_roses <- poll_rose_data %>% pollution_roses(met_data, num_breaks)
```

</div>

<br> 

<br>[Back to top](#start-finish)