usgs_metab_long-term_analysis_v1.Rmd

---
output:
  pdf_document: default
  html_document: default
---
---
title: 'Data for: Are annual river productivity regimes changing over time?'
author: "Nick Marzolf"
date: "2023-05-08"
output: html_document




```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Introduction

RMarkdown file to accompany: *Are annual river productivity regimes changing over time?,* a manuscript submitted to *Limnology and Oceanography: Letters* Special Issue on Changing Phenology in Aquatic Ecosystems.

Authors and affiliations: Nicholas S. Marzolf^1,\*^, Michael J. Vlah^1^, Heili E. Lowman^2^, Weston M. Slaughter^1^, Emily S. Bernhardt^1^

^1^ Department of Biology, Duke University, Durham, NC, USA

^2^ Department of Natural Resources and Environmental Science, University of Nevada, Reno, Reno, Nevada, USA

* Corresponding author: nicholas.marzolf@duke.edu (ORCiD: 0000-0001-9146-1643)

To recreate this HTML, load the data files into a folder named 'Data' at the same path as the .rmd file, open the .rmd file, and click 'Knit'.

# Load packages

```{r load-packages, echo = TRUE, message = FALSE, warning = FALSE}
# data manipulation
library(tidyverse)
library(dplyr)
library(lubridate)
library(purrr)
library(tibble)

# plotting
library(ggplot2)
library(ggpubr)
library(ggExtra)
library(GGally)
source("C:/Users/Nick Marzolf/Desktop/Research/R code/theme_nick.R")
ggplot2::theme_set(theme_nick())

# spatial
library(sf)
library(spData)


# statistics
library(lme4)
library(nlme)
library(performance)
library(EnvStats)
library(EflowStats)
library(car)
library(merTools)

# reproducability
library(pander)

# define constants and axis labels
metab_units_area <- (expression(paste('GPP (g C ', m^-2, ' ',y^-1,')', sep = ' ')))

low_p <- '#103801'
high_p <- '#43fa00'

# quotient to convert GPP from O2 to C
# g C = g O2 * (1 mol O2/32 g O2) * (2 mol O/1 mol O2) * (1 mol C/1 mol O) * (12 g C/1 mol C)
rq <- (1/32)*(2/1)*(1/1)*(12/1) 
```

# Explore sites

```{r sites, warning = FALSE, echo = TRUE}
# histogram of watershed areas
table_s1 <- readr::read_csv('data/data_citation/1_site_info.csv')
table_s1

readr::write_csv(table_s1,
                 'data/data_citation/ms_usgs_annual_metab/table_s1_site_info.csv')

table_s1_plot <- table_s1 %>% 
  sf::st_as_sf(.,
               coords = c('Longitude', 'Latitude'),
               crs = 4326)

fig_s1_map <- ggplot()+
  geom_sf(data = world %>% 
            dplyr::filter(name_long == 'United States') %>% 
            sf::st_set_crs(4326),
          fill = 'white')+
  geom_sf(data = table_s1_plot)+ 
  coord_sf(xlim = c(-130, -65), ylim = c(23, 55), expand = FALSE)


fig_s1_hist <- ggplot(table_s1,
                      aes(x = `Drainage Area (km2)`*2.58999))+ # convert to km2
  geom_histogram()+
  labs(x = expression(paste('Watershed area (', km^2,')')),
       y = 'Count')+
  scale_x_log10()
fig_s1_hist

fig_s1 <- cowplot::plot_grid(fig_s1_map,
                             fig_s1_hist,
                             ncol = 2, 
                             rel_widths = c(2,1),
                             axis = 'b',
                             align = 'v')
fig_s1

ggsave(plot = fig_s1,
       'manuscript/long-term_GPP/figures/fig_s1.png',
       dpi = 1200,
       width = 8, height = 3)
```

# Load GPP data

```{r load-data, warning = FALSE, echo = TRUE}

river_metab <- readr::read_csv('data/data_citation/5_usgs_metabolism.csv') 
fluxnet_metab <- readr::read_csv('data/data_citation/6_fluxnet_metabolism.csv')

```

# Calculate Annual Productivity

```{r calc-annual, warning = FALSE, echo = TRUE}

metab_all_daily <- rbind(river_metab %>% 
                           dplyr::select(site, date, GPP = GPP) %>% 
                           dplyr::mutate(source = 'USGS',
                                         biome = NA),
                         fluxnet_metab %>% 
                           dplyr::select(site, date, GPP, biome) %>% 
                           dplyr::mutate(source = 'Fluxnet'))

metab_all_annual <- metab_all_daily %>% 
  dplyr::mutate(year = lubridate::year(date)) %>% 
  dplyr::group_by(source, site, year, biome) %>% 
  dplyr::summarise(GPP_ann_C = dplyr::if_else('USGS' %in% source,
                                              sum(GPP*rq, na.rm = TRUE),
                                              sum(GPP, na.rm = TRUE)),
                   GPP_daily_cv = EnvStats::cv(GPP, 
                                               na.rm = TRUE, method = 'l.moments')*100)


# plot of annual productivity time-series
fig_s2 <- ggplot2::ggplot(metab_all_annual,
                          aes(x = year,
                              y = GPP_ann_C,
                              group = site,
                              color = GPP_ann_C))+
  ggplot2::geom_line()+
  ggplot2::geom_point()+
  ggplot2::scale_color_gradient(name = metab_units_area,
                                low = low_p,
                                high = high_p)+
  ggplot2::facet_grid(source ~ .)+
  ggplot2::ylab(metab_units_area)+
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = 'none')
fig_s2

ggsave(plot = fig_s2,
       'manuscript/long-term_GPP/figures/fig_s2.png',
       dpi = 1200, width = 6, height = 4)
```

```{r viz-ann-prod, warning = FALSE, echo = TRUE}
fig_1_top <- ggplot2::ggplot(data = metab_all_annual,
                             ggplot2::aes(x = reorder(site, 
                                                      desc(GPP_ann_C)), 
                                          y = GPP_ann_C,
                                          color = GPP_ann_C))+
  #ggplot2::geom_point()+
  ggplot2::geom_line()+
  ggplot2::ylab(metab_units_area)+
  ggplot2::facet_wrap(. ~ source,
                      scales = 'free')+
  ggplot2::scale_color_gradient(name = metab_units_area,
                                low = low_p,
                                high = high_p)+
  ggplot2::labs(x = 'Site')+
  ggplot2::scale_y_log10(limits = c(10, 10000))+
  ggplot2::theme(axis.text.x = ggplot2::element_blank(),
                 panel.grid.major = ggplot2::element_blank(), 
                 panel.grid.minor = ggplot2::element_blank(),
                 legend.position = 'none')
fig_1_top

# calculate long term CV
lt_cv_df <- metab_all_annual %>% 
  dplyr::group_by(site, source, biome) %>% 
  dplyr::summarise(lt_mean_GPP = mean(GPP_ann_C, na.rm = TRUE),
                   n_obs = length(GPP_ann_C),
                   lt_sd_GPP = sd(GPP_ann_C, na.rm = TRUE),
                   lt_se_GPP = lt_sd_GPP/n_obs,
                   lt_cv = lt_mean_GPP/sd(GPP_ann_C, na.rm = TRUE),
                   yr_range = paste(min(year), max(year), sep = '-'))  

# ID the 5 most variable FLUXNET sites
lt_cv_df_sub_flux <- data.frame(lt_cv_df) %>% 
  dplyr::arrange(desc(lt_cv)) %>% 
  dplyr::filter(source == 'Fluxnet',
                between(row_number(), 1, 5))
lt_cv_df_sub_flux


# ID the 5 most variable river sites
river_names <- readr::read_csv('data/data_citation/1_site_info.csv')  
river_names <- river_names %>% 
  dplyr::mutate(site = ifelse(nchar(river_names$`Site Number`) == 8,
                              paste0('nwis_', `Site Number`),
                              paste0('nwis_0', `Site Number`))) %>% 
  dplyr::select(site, 
                station_nm = Name) %>%
  dplyr::mutate(site_nm = str_sub(station_nm, 0, nchar(station_nm) - 3),
                state = str_sub(station_nm, -2, -1),
                clean_name = paste(site_nm, state, sep = ' '))

lt_cv_df_sub_flux_usgs <- data.frame(lt_cv_df) %>% 
  dplyr::filter(source == 'USGS') %>% 
  dplyr::arrange(desc(lt_cv)) %>% 
  dplyr::filter(between(row_number(), 1, 5)) %>% 
  dplyr::left_join(river_names, by = 'site')
lt_cv_df_sub_flux_usgs

# and plot
fig_1_bottom <- ggplot2::ggplot()+
  ggplot2::geom_dotplot(data = lt_cv_df,
                        aes(x = source, y = lt_cv,
                            fill = lt_mean_GPP,
                            group = lt_mean_GPP),
                        dotsize = 0.8,
                        binaxis='y', 
                        stackdir='center',
                        stackgroups=TRUE, 
                        method="histodot",
                        binwidth = 3)+
  ggplot2::scale_fill_gradient(name = metab_units_area,
                               low = low_p, # #855E09
                               high = high_p)+ # #40E304
  ggplot2::ylab(expression(paste(CV[GPP],' (%)')))+
  ggrepel::geom_label_repel(data = lt_cv_df_sub_flux,
                            aes(x = 1, 
                                y = lt_cv,
                                label = biome),
                            nudge_x = c(0.3,-0.3,0.3,-0.3,-0.3),
                            nudge_y = c(1,1,-5,0,-1),
                            size = 2)+
  ggrepel::geom_label_repel(data = lt_cv_df_sub_flux_usgs,
                            aes(x = 2, 
                                y = lt_cv,
                                label = clean_name),
                            size = 2,
                            nudge_x = 0.5,
                            nudge_y = c(50,40,30,20,10))+
  ggplot2::theme(axis.title.x = element_blank(),
                 legend.position = 'none')
fig_1_bottom


fig_1 <- cowplot::plot_grid(fig_1_top,
                            fig_1_bottom,
                            nrow = 2, 
                            ncol = 1,
                            common.legend = TRUE,
                            align = 'hv',
                            axis = 'l',
                            labels = 'auto')
fig_1

ggsave(plot = fig_1,
       'manuscript/long-term_GPP/figures/fig_1.png',
       dpi = 1000,
       width = 6.5, height = 5)
```

```{r fig-2, warning = FALSE, echo = TRUE}
extreme_prod_sites <- lt_cv_df %>%
  data.frame() %>% 
  dplyr::filter(source == 'USGS') %>% 
  dplyr::arrange(lt_mean_GPP) %>% 
  dplyr::slice(c(which.min(lt_mean_GPP), which.max(lt_mean_GPP))) %>% 
  dplyr::pull(site)


extreme_site_years <- metab_all_annual %>% 
  dplyr::filter(site %in% extreme_prod_sites) %>% 
  dplyr::group_by(site) %>% 
  dplyr::arrange(GPP_ann_C) %>% 
  dplyr::slice(c(which.min(GPP_ann_C), which.max(GPP_ann_C))) %>% 
  dplyr::mutate(site_year = paste(site, year, sep = '-'))

extreme_years <- river_metab %>% 
  dplyr::mutate(site_year = paste(site, lubridate::year(date), sep = '-')) %>% 
  dplyr::filter(site_year %in% extreme_site_years$site_year) %>% 
  dplyr::mutate(year = factor(lubridate::year(date)),
                jday = lubridate::yday(date)) %>% 
  dplyr::left_join(river_names, 'site') %>% 
  dplyr::mutate(clean_name_case = stringr::str_to_upper(clean_name, 'en'))

fig_2 <- ggplot(extreme_years,
                aes(x = jday, 
                    y = GPP*rq,
                    color = year))+
  geom_line(alpha = 0.7)+
  geom_point()+
  ggplot2::scale_color_manual(values = c(high_p, low_p,
                                         low_p,high_p))+
  facet_grid(clean_name_case ~ .,
             scales = 'free')+
  labs(x = 'Day of Year',
       y = expression(paste('GPP (g C ', m^-2,' ',d^-1,')')))+
  theme(legend.position = 'none',
        strip.text = element_text(size = 6),
        strip.background = element_rect(fill = 'white'))
fig_2

ggsave(plot = fig_2,
       'manuscript/long-term_GPP/figures/fig_2.png',
       dpi = 1000, 
       width = 5, height = 6)
```

# Productivity trends

```{r prod-trends, warning = FALSE, echo = TRUE}
river_metab_ann <- metab_all_annual %>%
  dplyr::filter(source == 'USGS') 

annual_C_sum_trends <- river_metab_ann %>% 
  dplyr::select(site, year, GPP_ann_C) %>% 
  dplyr::group_by(site) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(sens = purrr::map(data, ~trend::sens.slope(x = .$GPP_ann_C)),
                mk_test = purrr::map(data, ~trend::mk.test(x = .$GPP_ann_C))) %>% 
  dplyr::mutate(sens_p = purrr::map(sens, ~.$p.value), 
                sens_s = purrr::map_dbl(sens, ~.$estimates),
                mk_p = purrr::map_dbl(mk_test, ~.$p.value),
                mk_s = purrr::map_dbl(mk_test, ~.$estimates['S']),
                sens_sig = ifelse(sens_p <= 0.05, 'significant', 'non-significant'),
                sens_slope = ifelse(sens_s > 0, 'increasing', 'decreasing'),
                mk_sig = ifelse(mk_p <= 0.05, 'significant', 'non-significant'),
                mk_slope = ifelse(mk_s > 0, 'increasing', 'decreasing')) 


annual_C_sum_sigs <- annual_C_sum_trends %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(sens_p = unlist(sens_p)) %>% 
  dplyr::select(site, sens_s, sens_p, sens_sig, sens_slope,
                mk_p, mk_s, mk_sig, mk_slope)


n_sig <- annual_C_sum_sigs %>% 
  dplyr::filter(sens_p <= 0.05) %>% 
  dplyr::group_by(sens_slope) %>% 
  summarise(n = n()) %>% 
  data.frame()


river_trends <- river_metab_ann %>% 
  dplyr::left_join(annual_C_sum_sigs) 

fig_3 <- ggplot2::ggplot(river_trends,
                         ggplot2::aes(x = year, 
                                      y = GPP_ann_C,
                                      color = sens_slope))+
  ggplot2::geom_line(aes(group = site),
                     linewidth = 1.25)+
  scale_y_log10(limits = c(10, 3000),
                breaks = c(10,30,100,300,1000,3000),
                labels = c(10,30,100,300,1000,3000))+
  scale_x_continuous(n.breaks = 10)+
  gghighlight::gghighlight(sens_sig == 'significant',
                           use_direct_label = FALSE,
                           keep_scales = TRUE,
                           unhighlighted_params = list(alpha("grey", 0.4),
                                                       linewidth = 0.1))+
  ggplot2::labs(y = metab_units_area)+
  scale_color_manual(name = element_blank(),
                     values = c(low_p,high_p),
                     labels = c(glue::glue('Decreasing (n = ', n_sig[1,2],')'),
                                glue::glue('Increasing (n = ', n_sig[2,2],')')))+
  theme(axis.title.x = element_blank(),
        legend.justification = c(0.95,0.025),
        legend.position = c(0.95, 0.025),
        legend.text = element_text(size = 10),
        legend.background = element_rect(color = 'black',
                                         linewidth = 0.25),
        legend.margin = margin(0, 2, 0.5, 0.1))
fig_3

ggsave(plot = fig_3,
       'manuscript/long-term_GPP/figures/fig_3.png',
       dpi = 1200,
       width = 6, height = 4)

table_s2 <- annual_C_sum_sigs %>% 
  dplyr::select(site, sens_s, sens_p, mk_s,mk_p)
colnames(table_s2) <- c('Site',
                        "Sen's Slope",
                        "Sen's p-value",
                        "Mann-Kendall S",
                        "Mann-Kendall p-value")
readr::write_csv(table_s2,
                 'data/data_citation/ms_usgs_annual_metab/table_s2_trend_stats.csv')

```

# Drivers of annual productivity

```{r prod-drivers, warning = FALSE, echo = TRUE}
river_light <- readr::read_csv('data/data_citation/7_river_light.csv')

library(BernhardtMetabolism)
source('code/functions/fill_and_normalize_metab.R')
river_metab_filled <- fill_and_normalize_metab(river_metab)

river_metab_wLight <- river_metab_filled %>% 
  dplyr::left_join(.,
                   river_light, 
                   by = c('site', 'date'))


river_metab_ann_drivers <- river_metab_wLight %>% 
  dplyr::select(date, site, GPP_filled, discharge, daily_PAR, temp.water) %>% 
  dplyr::group_by(site,
                  year = lubridate::year(date)) %>% 
  dplyr::summarise(GPP_ann_C = sum(GPP_filled*rq, na.rm = TRUE),         # g C m-2 y-1
                   Q_ann_mean = mean(discharge, na.rm = TRUE),
                   Q_ann_cv = EnvStats::cv(discharge, method = 'l.moments'),
                   light_ann_tot = sum(daily_PAR, na.rm = TRUE),
                   temp_ann_mean = mean(temp.water, na.rm = TRUE)) %>%      # mol m-2 y-1
  dplyr::filter(light_ann_tot != 0)  

# check the normality of inputs data
hist(river_metab_ann_drivers$GPP_ann_C) # log transform me!
hist(log10(river_metab_ann_drivers$GPP_ann_C))
hist(river_metab_ann_drivers$Q_ann_cv)
hist(river_metab_ann_drivers$temp_ann_mean)
hist(river_metab_ann_drivers$light_ann_tot)

# check for site-specific outliers
ggplot(river_metab_ann_drivers, aes(x = site, y = GPP_ann_C))+
  geom_boxplot()


# co-linearity of drivers
ggpairs(data = river_metab_ann_drivers %>% 
          dplyr::select(site, Q_ann_cv, light_ann_tot, temp_ann_mean),
        cardinality_threshold = 47)



# all 3 drivers
gpp_mult_reg <- lme4::lmer(data = river_metab_ann_drivers,
                           log10(GPP_ann_C) ~ Q_ann_cv + light_ann_tot + temp_ann_mean + (1|site))

# summarise and evaluate outputs
summary(gpp_mult_reg)
car::Anova(gpp_mult_reg)


fig_s3_a <- ggplot(data = data.frame(resid(gpp_mult_reg)))+
  geom_histogram(aes(x = resid.gpp_mult_reg.))+
  labs(x = 'Residuals',
       y = 'Count',
       title = 'Light + CVQ + Temp')
fig_s3_a

qqnorm(resid(gpp_mult_reg))

## extract random effects
mult_ranef <- nlme::ranef(gpp_mult_reg)[[1]] %>% 
  data.frame() %>% 
  tibble::rownames_to_column() %>% 
  dplyr::rename(site = rowname,
                ranef = X.Intercept.) %>% 
  dplyr::mutate(effect = 'interaction')

## fixed effects
nlme::fixef(gpp_mult_reg)


# light only
gpp_reg_light <- lme4::lmer(data = river_metab_ann_drivers,
                            log10(GPP_ann_C) ~ light_ann_tot + (1|site))
summary(gpp_reg_light)

fig_s3_b <- ggplot(data = data.frame(resid(gpp_reg_light)))+
  geom_histogram(aes(x = resid.gpp_reg_light.))+
  labs(x = 'Residuals',
       y = 'Count',
       title = 'Light')
fig_s3_b

## light only random effects
light_ranef <- nlme::ranef(gpp_reg_light)[[1]] %>% 
  data.frame() %>% 
  tibble::rownames_to_column() %>% 
  dplyr::rename(site = rowname,
                ranef = X.Intercept.) %>% 
  dplyr::mutate(effect = 'light')


# cv Q only
gpp_reg_cvQ <- lme4::lmer(data = river_metab_ann_drivers,
                          log10(GPP_ann_C) ~ Q_ann_cv + (1|site))
summary(gpp_reg_cvQ)

fig_s3_c <- ggplot(data = data.frame(resid(gpp_reg_cvQ)))+
  geom_histogram(aes(x = resid.gpp_reg_cvQ.))+
  labs(x = 'Residuals',
       y = 'Count',
       title = 'CVQ')
fig_s3_c

## CV_Q random effects
cvQ_ranef <- nlme::ranef(gpp_reg_cvQ)[[1]] %>% 
  data.frame() %>% 
  tibble::rownames_to_column() %>% 
  dplyr::rename(site = rowname,
                ranef = X.Intercept.) %>% 
  dplyr::mutate(effect = 'cv_Q')


# temperature only
gpp_reg_temp <- lme4::lmer(data = river_metab_ann_drivers,
                           log10(GPP_ann_C) ~ temp_ann_mean + (1|site))

summary(gpp_reg_temp)

fig_s3_d <- ggplot(data = data.frame(resid(gpp_reg_temp)))+
  geom_histogram(aes(x = resid.gpp_reg_temp.))+
  labs(x = 'Residuals',
       y = 'Count',
       title = 'Temp')
fig_s3_d



temp_ranef <- nlme::ranef(gpp_reg_temp)[[1]] %>% 
  data.frame() %>% 
  tibble::rownames_to_column() %>% 
  dplyr::rename(site = rowname,
                ranef = X.Intercept.) %>% 
  dplyr::mutate(effect = 'temperature')


# compile and compare model outputs

## print residual plots
fig_s3 <- cowplot::plot_grid(fig_s3_a,
                             fig_s3_b,
                             fig_s3_c,
                             fig_s3_d,
                             nrow = 2, ncol = 2)

ggsave(plot = fig_s3,
       'manuscript/long-term_GPP/figures/fig_s3.png',
       dpi = 1200, 
       width = 5, height = 5)

## Random effects (ie site-level effects)
table_s3 <- list(mult_ranef,
                 light_ranef,
                 temp_ranef,
                 cvQ_ranef) %>% 
  reduce(rbind) %>% 
  pivot_wider(names_from = effect,
              values_from = ranef)
table_s3[,-1] <- round(table_s3[,-1], 1)
colnames(table_s3) <- c('Site', 'Interaction', 'Light', 'Temperature', 'CVQ')

table_s3

readr::write_csv(table_s3,
                 'data/data_citation/ms_usgs_annual_metab/table_s3_random_effects.csv')

# Compare models by AICw
table_1 <- performance::compare_performance(gpp_reg_temp,
                                            gpp_reg_cvQ, 
                                            gpp_reg_light,
                                            gpp_mult_reg) 
table_1

# predict from best model
pred_river_metab_ann_drivers <- river_metab_ann_drivers %>% 
  dplyr::group_by(site, Q_ann_cv, light_ann_tot, temp_ann_mean) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(GPP_ann_C_pred = 10^(stats::predict(gpp_mult_reg, 
                                                    river_metab_ann_drivers,
                                                    re.form = ~(1|site))),
                intervals = merTools::predictInterval(gpp_mult_reg,
                                                      level = 0.95,
                                                      newdata = river_metab_ann_drivers,
                                                      n.sims = 1000,
                                                      stat = 'mean',
                                                      type = 'linear.prediction'),
                GPP_ann_C = river_metab_ann_drivers$GPP_ann_C)

cond_r2 <- round(performance::performance(gpp_mult_reg)$R2_conditional, 2)

line_1 <- expression(paste('Predicted GPP (g C ', m^-2, ' ',y^-1,')'))  
line_2 <- expression(paste(log[10],'(GPP) ~ MAT + ',CV[Q],' + PAR + (1|site)'))


fig_4 <- ggplot(pred_river_metab_ann_drivers,
                aes(x = 10^(intervals$fit)))+
  geom_ribbon(aes(ymax = 10^(intervals$upr),
                  ymin = 10^(intervals$lwr)),
              alpha = 0.2)+
  geom_point(aes(y = GPP_ann_C))+
  geom_abline(slope = 1,
              intercept = 0,
              linetype = 'dashed',
              color = 'black')+
  labs(x = "",
       y = expression(paste('Estimated GPP (g C ', m^-2, ' ',y^-1,')')))+
  scale_x_log10(limits = c(50, 3000),
                breaks = c(50, 100, 500, 1000, 3000),
                labels = c(50, 100, 500, 1000, 3000))+
  scale_y_log10(limits = c(30, 5000),
                breaks = c(30, 100, 1000, 5000),
                labels = c(30, 100, 1000, 5000))+
  # lims(x = c(0, 2500),
  #      y = c(0, 4000))+
  annotate('text',
           x = 2000, y = -400, 
           label = bquote(R[Conditional]^2 == .(cond_r2)))

fig_4 <- cowplot::ggdraw(fig_4)+
  cowplot::draw_label(line_1, x = 0.55, y = 0.06, size = 12)+
  cowplot::draw_label(line_2, x = 0.55, y = 0.02, size = 10)

fig_4

ggsave(plot = fig_4,
       'manuscript/long-term_GPP/figures/fig_4.png',
       dpi = 1200,
       width = 5, height = 4)

```

# Are the drivers changing?

```{r changers, warning=FALSE, echo = TRUE}
changing_sites <- river_trends %>% 
  dplyr::filter(sens_sig == 'significant') %>% 
  dplyr::group_by(site) %>% 
  dplyr::select(site, sens_slope) %>% 
  unique()


changers_drivers <- river_metab_wLight %>% 
  dplyr::filter(site %in% changing_sites$site) %>% 
  group_by(site, year = lubridate::year(date)) %>% 
  summarise(GPP_ann_C = sum(GPP*rq, na.rm = TRUE),
            mat = mean(temp.water, na.rm = TRUE),
            cv_q = EnvStats::cv(discharge, method = 'l.moments'),
            ann_light = sum(daily_PAR, na.rm = TRUE)) %>% 
  dplyr::left_join(changing_sites, 'site')

mat_trends <- changers_drivers %>% 
  dplyr::select(site, year, mat) %>% 
  dplyr::group_by(site) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(sens = purrr::map(data, ~trend::sens.slope(x = .$mat)),
                mk_test = purrr::map(data, ~trend::mk.test(x = .$mat))) %>% 
  dplyr::mutate(sens_p = purrr::map(sens, ~.$p.value), 
                sens_s = purrr::map_dbl(sens, ~.$estimates),
                mk_p = purrr::map_dbl(mk_test, ~.$p.value),
                mk_s = purrr::map_dbl(mk_test, ~.$estimates['S']),
                sens_sig = ifelse(sens_p <= 0.05, 'significant', 'non-significant'),
                sens_slope = ifelse(sens_s > 0, 'increasing', 'decreasing'),
                mk_sig = ifelse(mk_p <= 0.05, 'significant', 'non-significant'),
                mk_slope = ifelse(mk_s > 0, 'increasing', 'decreasing')) 


mat_sigs <- mat_trends %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(sens_p = unlist(sens_p)) %>% 
  dplyr::select(site, sens_sig, sens_slope, sens_s, sens_p)

n_sig_mat <- mat_sigs %>% 
  dplyr::filter(sens_p <= 0.05) %>% 
  dplyr::group_by(sens_slope) %>% 
  summarise(n = n()) %>% 
  data.frame()


cvq_trends <- changers_drivers %>% 
  dplyr::select(site, year, cv_q) %>% 
  dplyr::group_by(site) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(sens = purrr::map(data, ~trend::sens.slope(x = .$cv_q)),
                mk_test = purrr::map(data, ~trend::mk.test(x = .$cv_q))) %>% 
  dplyr::mutate(sens_p = purrr::map(sens, ~.$p.value), 
                sens_s = purrr::map_dbl(sens, ~.$estimates),
                mk_p = purrr::map_dbl(mk_test, ~.$p.value),
                mk_s = purrr::map_dbl(mk_test, ~.$estimates['S']),
                sens_sig = ifelse(sens_p <= 0.05, 'significant', 'non-significant'),
                sens_slope = ifelse(sens_s > 0, 'increasing', 'decreasing'),
                mk_sig = ifelse(mk_p <= 0.05, 'significant', 'non-significant'),
                mk_slope = ifelse(mk_s > 0, 'increasing', 'decreasing')) 


cvq_sigs <- cvq_trends %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(sens_p = unlist(sens_p)) %>% 
  dplyr::select(site, sens_sig, sens_slope, sens_s, sens_p)

n_sig_cvq <- cvq_trends %>% 
  dplyr::filter(sens_p <= 0.05) %>% 
  dplyr::group_by(sens_slope) %>% 
  summarise(n = n()) %>% 
  data.frame()


light_trends <- changers_drivers %>% 
  dplyr::select(site, year, ann_light) %>% 
  dplyr::group_by(site) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(sens = purrr::map(data, ~trend::sens.slope(x = .$ann_light)),
                mk_test = purrr::map(data, ~trend::mk.test(x = .$ann_light))) %>% 
  dplyr::mutate(sens_p = purrr::map(sens, ~.$p.value), 
                sens_s = purrr::map_dbl(sens, ~.$estimates),
                mk_p = purrr::map_dbl(mk_test, ~.$p.value),
                mk_s = purrr::map_dbl(mk_test, ~.$estimates['S']),
                sens_sig = ifelse(sens_p <= 0.05, 'significant', 'non-significant'),
                sens_slope = ifelse(sens_s > 0, 'increasing', 'decreasing'),
                mk_sig = ifelse(mk_p <= 0.05, 'significant', 'non-significant'),
                mk_slope = ifelse(mk_s > 0, 'increasing', 'decreasing')) 


table_s4 <- changing_sites %>% 
  left_join(mat_trends %>% 
              dplyr::select(site, 
                            mat_mk_p = mk_p), 
            'site') %>% 
  left_join(cvq_trends %>% 
              dplyr::select(site, 
                            cvq_mk_p = mk_p), 
            'site') %>% 
  left_join(light_trends %>% 
              dplyr::select(site, 
                            light_mk_p = mk_p), 
            'site') %>% 
  rename(ann_GPP_trend = sens_slope)
table_s4

readr::write_csv(table_s4,
                 'data/data_citation/ms_usgs_annual_metab/table_s4_driver_trends.csv')
```

# Session Info

```{r session-info}
pander(sessionInfo())
```