Adding script for final analyses of total pesticide concentrations in…

… zooplankton and water for summer-fall period - script is unfinished

manuscript_contam/src/analysis_conc_water_zoop_summer_fall.R

@@ -0,0 +1,186 @@
+# NDFS Contaminants Manuscript
+# Purpose: Run statistical analyses for the total pesticide concentration data
+  # in water and zooplankton collected in summer-fall periods for the NDFS
+  # contaminants manuscript. Exploratory analyses and figures are in
+  # ND-FASTR/manuscript_contam/notebooks/explore_contam_conc_water_zoop.Rmd
+# Author: Dave Bosworth
+# Contacts: David.Bosworth@water.ca.gov
+
+# Load packages
+library(tidyverse)
+library(car)
+library(emmeans)
+library(broom)
+library(here)
+library(conflicted)
+
+# Declare package conflict preferences 
+conflicts_prefer(dplyr::filter())
+
+# Check if we are in the correct working directory
+i_am("manuscript_contam/src/analysis_conc_water_zoop_summer_fall.R")
+
+
+# 1. Import and Prepare Data ----------------------------------------------
+
+# Define root directory for pesticide data:
+fp_pest <- here("manuscript_contam/data/processed")
+
+# Total Pesticide Concentration data:
+df_conc_all <- read_rds(file.path(fp_pest, "contam_conc_water_zoop_2017-2019.rds"))
+
+# Pesticide application data:
+df_appl <- read_rds(file.path(fp_pest, "pesticide_use_daily_tot_2017-2020.rds"))
+
+# Daily average flow data:
+df_davg_flow <- read_rds(here("Water_Quality/Data_Processed/LIS_SR_DailyAvgFlow_2013-2022.rds"))
+
+# Prepare the pesticide application data to be joined to the concentration data
+  # by calculating monthly totals and standardizing to area of the watershed.
+df_appl_c <- df_appl %>% 
+  # Calculate monthly totals by region
+  summarize(TotalAppl = sum(TotalApplication), .by = c(Region, Year, Month)) %>% 
+  mutate(
+    # Standardize total application to area (in square miles) of the watershed 
+    TotalApplStd = if_else(
+      Region == "Sacramento River", 
+      TotalAppl/22526.79948,
+      TotalAppl/4217.503529
+    ),
+    # Define station that each region is assigned to
+    Station = case_match(
+      Region,
+      "Toe Drain" ~ "STTD", 
+      "Sacramento River" ~ "SHR"
+    )
+  )
+
+# Prepare the flow data to be joined to the concentration data using SR at
+  # Freeport flow data for SHR and LIS flow data for STTD. Use Daily Average Net
+  # flows for both.
+df_davg_flow_c <- df_davg_flow %>% 
+  select(-LIS_DailyAvgInstFlow) %>% 
+  pivot_longer(
+    cols = ends_with("NetFlow"),
+    names_to = "Station",
+    values_to = "DailyAvgNetFlow"
+  ) %>% 
+  mutate(
+    Station = case_match(
+      Station,
+      "SR_DailyAvgNetFlow" ~ "SHR",
+      "LIS_DailyAvgNetFlow" ~ "STTD"
+    )
+  )
+
+# Join pesticide concentration, application, and flow data together, and
+  # restrict the date ranges to summer-fall period (July 1 - Oct 31) for each year
+df_conc_all_c <- df_conc_all %>% 
+  mutate(
+    Year = year(Date),
+    Month = month(Date, label = TRUE)
+  ) %>% 
+  left_join(df_appl_c, by = join_by(Station, Year, Month)) %>% 
+  left_join(df_davg_flow_c, by = join_by(Date, Station)) %>% 
+  mutate(Month = as.numeric(Month)) %>%
+  filter(Month %in% 7:10) %>% 
+  select(Date, Year, Station, starts_with("TotalConc"), DailyAvgNetFlow, TotalApplStd)
+
+
+# 2. Zooplankton Analyses -------------------------------------------------
+
+# Prepare zooplankton concentration data for analysis
+df_zoop <- df_conc_all_c %>% 
+  select(-TotalConc_Water) %>% 
+  drop_na() %>% 
+  # Remove two samples with total concentrations equal to zero
+  filter(TotalConc_Zoop > 0) %>% 
+  mutate(
+    # Add log-transformed and sqrt-transformed zooplankton concentrations
+    TotalConc_Zoop_log = log(TotalConc_Zoop),
+    TotalConc_Zoop_sqrt = sqrt(TotalConc_Zoop),
+    # Convert year to character for categorical models
+    Year = as.character(Year)
+  )
+
+
+# 2.1 Linear models with flow and application -----------------------------
+
+# Run linear models separately for each station
+
+# SHR: 
+df_zoop_shr <- df_zoop %>% filter(Station == "SHR")
+
+# Run model for SHR using log transformed zooplankton concentrations since the
+  # model diagnostics looked best with these
+m_zoop_q_appl_shr_log <- df_zoop_shr %>% 
+  lm(TotalConc_Zoop_log ~ DailyAvgNetFlow + TotalApplStd, data = .)
+
+# Create ANOVA table for model and convert it to a tibble for export
+df_an_zoop_q_appl_shr <- tidy(anova(m_zoop_q_appl_shr_log)) %>% 
+  mutate(Response = "Zooplankton", Station = "SHR", .before = term)
+
+# STTD:
+df_zoop_sttd <- df_zoop %>% filter(Station == "STTD")
+
+# Run model for STTD using square root transformed zooplankton concentrations
+  # since the model diagnostics looked best with these
+m_zoop_q_appl_sttd_sqrt <- df_zoop_sttd %>% 
+  lm(TotalConc_Zoop_sqrt ~ DailyAvgNetFlow + TotalApplStd, data = .)
+
+# Create ANOVA table for model and convert it to a tibble for export
+df_an_zoop_q_appl_sttd <- tidy(anova(m_zoop_q_appl_sttd_sqrt)) %>% 
+  mutate(Response = "Zooplankton", Station = "STTD", .before = term)
+
+
+# 2.2 Linear model with station and year ----------------------------------
+
+# Run model using log transformed zooplankton concentrations since the model
+  # diagnostics looked best with these. Also run model without interaction between
+  # Station and Year since it wasn't significant.
+m_zoop_sta_yr_log <- df_zoop %>% lm(TotalConc_Zoop_log ~ Station + Year, data = .)
+
+# Create ANOVA table for model and convert it to a tibble for export
+# use type 2 sum of squares because of the unbalanced design
+df_an_zoop_sta_yr <- tidy(Anova(m_zoop_sta_yr_log, type = 2))
+
+# The main effect of station is not significant meaning that total pesticide
+  # concentrations in zooplankton don't differ between SHR and STTD. However, the
+  # main effect of Year is significant, so we'll export the pairwise contrasts to
+  # show how concentrations differ among the years. We'll also export the pairwise
+  # contrasts for station since it is close to significant.
+
+# Estimated marginal means for Year main effect
+em_zoop_yr <- emmeans(m_zoop_sta_yr_log, specs = "Year")
+
+# Estimated marginal means for Station main effect
+em_zoop_sta <- emmeans(m_zoop_sta_yr_log, specs = "Station")
+
+# Create table of contrasts and convert it to a tibble for export
+df_em_zoop_sta_yr <- bind_rows(tidy(pairs(em_zoop_yr)), tidy(pairs(em_zoop_sta)))
+
+
+# 3. Water Analyses -------------------------------------------------------
+
+# Prepare water concentration data for analysis
+df_water <- df_conc_all_c %>% 
+  select(-TotalConc_Zoop) %>% 
+  # only include 2019 since that was the only year water samples were collected
+    # at STTD
+  filter(Year == 2019) %>% 
+  drop_na() %>% 
+  # Add log-transformed water concentrations
+  mutate(TotalConc_Water_log = log(TotalConc_Water))
+
+# 3.1 Linear models with flow and application -----------------------------
+
+# Run linear models separately for each station
+
+# SHR: 
+df_water_shr <- df_water %>% filter(Station == "SHR")
+
+
+
+# STTD:
+df_water_sttd <- df_water %>% filter(Station == "STTD")
+