Nutrient and chlorophyll analysis
Rosemary Hartman
August 01, 2023
Load packages and functions
library(multcomp)
library(car)
library(broom)
library(tidyverse)
library(lubridate)
library(emmeans)
library(lme4)
library(lmerTest)
library(DHARMa)
library(effects)
library(visreg)
library(knitr)
library(kableExtra)
library(here)
source(here("src/analysis/global_analysis_functions.R"))Display R and package versions:
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.2.3 (2023-03-15 ucrt)
## os Windows 10 x64 (build 19044)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate English_United States.utf8
## ctype English_United States.utf8
## tz America/Los_Angeles
## date 2023-08-01
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
## backports 1.4.1 2021-12-13 [1] CRAN (R 4.2.0)
## boot 1.3-28.1 2022-11-22 [2] CRAN (R 4.2.3)
## broom * 1.0.4 2023-03-11 [1] CRAN (R 4.2.2)
## bslib 0.4.2 2022-12-16 [1] CRAN (R 4.2.2)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.2.3)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.2)
## car * 3.1-2 2023-03-30 [1] CRAN (R 4.2.3)
## carData * 3.0-5 2022-01-06 [1] CRAN (R 4.2.1)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.3)
## coda 0.19-4 2020-09-30 [1] CRAN (R 4.2.1)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.2.3)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.2)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.1)
## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.2.1)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.2.1)
## DHARMa * 0.4.6 2022-09-08 [1] CRAN (R 4.2.1)
## digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.2)
## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.2.3)
## effects * 4.2-2 2022-07-13 [1] CRAN (R 4.2.2)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.1)
## emmeans * 1.8.5 2023-03-08 [1] CRAN (R 4.2.2)
## estimability 1.4.1 2022-08-05 [1] CRAN (R 4.2.1)
## evaluate 0.21 2023-05-05 [1] CRAN (R 4.2.3)
## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.2.2)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.2)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.2.2)
## fs 1.6.2 2023-04-25 [1] CRAN (R 4.2.3)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.1)
## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.2.3)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.1)
## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.2.3)
## here * 1.0.1 2020-12-13 [1] CRAN (R 4.2.1)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.2.3)
## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.2.3)
## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.2.3)
## httpuv 1.6.9 2023-02-14 [1] CRAN (R 4.2.2)
## httr 1.4.6 2023-05-08 [1] CRAN (R 4.2.3)
## insight 0.19.1 2023-03-18 [1] CRAN (R 4.2.3)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.1)
## jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.2.2)
## kableExtra * 1.3.4 2021-02-20 [1] CRAN (R 4.2.1)
## knitr * 1.42 2023-01-25 [1] CRAN (R 4.2.2)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.2.1)
## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.2.3)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.1)
## lme4 * 1.1-32 2023-03-14 [1] CRAN (R 4.2.3)
## lmerTest * 3.1-3 2020-10-23 [1] CRAN (R 4.2.1)
## lubridate * 1.9.2 2023-02-10 [1] CRAN (R 4.2.2)
## magrittr * 2.0.3 2022-03-30 [1] CRAN (R 4.2.1)
## MASS * 7.3-58.2 2023-01-23 [2] CRAN (R 4.2.3)
## Matrix * 1.5-4 2023-04-04 [1] CRAN (R 4.2.3)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.1)
## mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.1)
## minqa 1.2.5 2022-10-19 [1] CRAN (R 4.2.2)
## mitools 2.4 2019-04-26 [1] CRAN (R 4.2.2)
## multcomp * 1.4-23 2023-03-09 [1] CRAN (R 4.2.2)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.1)
## mvtnorm * 1.1-3 2021-10-08 [1] CRAN (R 4.2.0)
## nlme 3.1-162 2023-01-31 [2] CRAN (R 4.2.3)
## nloptr 2.0.3 2022-05-26 [1] CRAN (R 4.2.1)
## nnet 7.3-18 2022-09-28 [2] CRAN (R 4.2.3)
## numDeriv 2016.8-1.1 2019-06-06 [1] CRAN (R 4.2.0)
## patchwork * 1.1.2 2022-08-19 [1] CRAN (R 4.2.1)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3)
## pkgbuild 1.4.0 2022-11-27 [1] CRAN (R 4.2.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.1)
## pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.2.2)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.1)
## processx 3.8.1 2023-04-18 [1] CRAN (R 4.2.3)
## profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.1)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.1)
## ps 1.7.5 2023-04-18 [1] CRAN (R 4.2.3)
## purrr * 1.0.1 2023-01-10 [1] CRAN (R 4.2.2)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.1)
## Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.2)
## readr * 2.1.4 2023-02-10 [1] CRAN (R 4.2.2)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.1)
## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.3)
## rmarkdown 2.21 2023-03-26 [1] CRAN (R 4.2.3)
## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.2.1)
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.1)
## rvest 1.0.3 2022-08-19 [1] CRAN (R 4.2.1)
## sandwich 3.0-2 2022-06-15 [1] CRAN (R 4.2.1)
## sass 0.4.6 2023-05-03 [1] CRAN (R 4.2.3)
## scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.1)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.1)
## shiny 1.7.4 2022-12-15 [1] CRAN (R 4.2.2)
## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.2.2)
## stringr * 1.5.0 2022-12-02 [1] CRAN (R 4.2.2)
## survey 4.1-1 2021-07-19 [1] CRAN (R 4.2.2)
## survival * 3.5-5 2023-03-12 [1] CRAN (R 4.2.3)
## svglite 2.1.1 2023-01-10 [1] CRAN (R 4.2.2)
## systemfonts 1.0.4 2022-02-11 [1] CRAN (R 4.2.1)
## TH.data * 1.1-2 2023-04-17 [1] CRAN (R 4.2.3)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.2.3)
## tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.2.2)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.1)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.2.2)
## timechange 0.2.0 2023-01-11 [1] CRAN (R 4.2.2)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.2.3)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.1)
## usethis 2.1.6 2022-05-25 [1] CRAN (R 4.2.1)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.2)
## vctrs 0.6.2 2023-04-19 [1] CRAN (R 4.2.3)
## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.2.3)
## visreg * 2.7.0 2020-06-04 [1] CRAN (R 4.2.1)
## webshot 0.5.4 2022-09-26 [1] CRAN (R 4.2.1)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.1)
## xfun 0.39 2023-04-20 [1] CRAN (R 4.2.3)
## xml2 1.3.4 2023-04-27 [1] CRAN (R 4.2.3)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.1)
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.2)
## zoo 1.8-12 2023-04-13 [1] CRAN (R 4.2.3)
##
## [1] C:/R/win-library/4.2
## [2] C:/Program Files/R/R-4.2.3/library
##
## ──────────────────────────────────────────────────────────────────────────────Load data
# Import long-term average data for the nutrients
nutr_data <- read_rds(here("data/processed/wq/lt_avg_nutr.rds"))
# Import long-term average data for chlorophyll
chla_data <- read_rds(here("data/processed/wq/lt_avg_chla.rds"))
# Combine nutrient and chlorophyll data and prepare for analyses
nutr_chla_data_c <-
full_join(nutr_data, chla_data, by = join_by(YearAdj, Season, Region, SVIndex, YearType, Drought)) %>%
mutate(
# Change region names to "South-Central Delta" and "North Delta"
Region = case_when(
Region == "SouthCentral" ~ "South-Central Delta",
Region == "North" ~ "North Delta",
TRUE ~ Region
),
Region = factor(
Region,
levels = c("Suisun Marsh", "Suisun Bay", "Confluence", "South-Central Delta", "North Delta")
),
Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Fall")),
Year_fac = factor(YearAdj),
YearType = factor(YearType, levels = c("Critical", "Dry", "Below Normal", "Above Normal", "Wet")),
Drought = factor(Drought, levels = c("D", "N", "W")),
YearAdj_s = (YearAdj - mean(YearAdj) / sd(YearAdj)),
across(c(starts_with("Diss"), Chlorophyll), list(l = log))
)Analyses
Ammonia
I’m going to start with Ammonia. My initial look at Ammonia made it seem pretty boring. I’d lean towards leaving it out of the paper, but I”ll see if including some sort of dummy variable for WWTP upgrades helps at all.
Cloern 2019 discussed a lot of Ammonia patterns: Ammonium While the silicate source to San Francisco Bay is river inflow, the ammonium source is primarily municipal wastewater (Jassby 2008). Comparison of silicate and ammonium variability provides an opportunity to explore the differing patterns of nutrients delivered by runoff and point sources. The trends of increasing ammonium were significant for most months (Fig. 6B), so I explored models to explain variability of annual mean ammonium concentration. The best-fitting simple model explained 59% of ammonium variability with two processes— Outflow and wastewater inputs of TKN (Fig. 8D). The size effects of Outflow (−1.0) and TKN loading (+0.8) were comparable (Table 4). The relationship with Outflow is nonlinear and negative, with highest ammonium concentrations during dry years and lowest during wet years (Fig. 8D). This is opposite the Outflow-silicate relationship because river inflow is a source of silicate but it dilutes wastewater-derived ammonium. The shape of nutrient-Outflow relationships can therefore give information about the relative importance of riverine and point-source inputs. Ammonium concentration in the estuary was significantly and linearly related to TKN loading (Fig. 8D), so the long-term trend of increasing ammonium concentration (Fig. 4; Table 2) can be attributed to increasing wastewater loading over time (Fig. 7E), but not to changes in Outflow and dilution.
Ammonium at the monthly scale Analyses of water-quality time series data can yield different insights into processes depending on the time scale considered. As an example, I explored regression models to identify the key drivers of ammonium variability at the monthly scale to compare against analyses at the annual scale. This was motivated by the strong monthly pattern of ammonium trends—largest increases in winter and smallest in summer (Fig. 6B), suggesting a possible influence of seasonal temperature variability. The data series included 346 monthly measurements beginning January 1987. The best fitting model described monthly ammonium variability as a function of three variables—Outflow and TKN loading (as above), and temperature. Monthly ammonium was a nonlinear and negative function of Outflow and a linear positive function of TKN loading (Fig. 9)—the same functional responses observed at the annual scale (Fig. 8D). The Outflow effect was −2.1 and the TKN-loading effect was +1.0. However, at the monthly time scale, temperature had the largest effect (−5.8) through a negative linear relationship. Thus, both annual and monthly variability of ammonium, including trends of increase over time, are associated with variability of TKN loading as a source and Outflow through its dilution effect. However, a stronger temperature effect was revealed by analyzing monthly data. This effect is not evident in analyses of the annual data because annual temperature fluctuations are small relative to the high-amplitude monthly variability (Fig. 4A). The seasonal temperature effect on ammonium concentration in the estuary is presumably a response to the near-linear relationship between nitrification rate and temperature (e.g., Sudarno et al. 2011). Highest nitrification rates during the warm season would explain the negative relationship between ammonium and temperature (Fig. 9A), the seasonal pattern of lowest ammonium concentrations during summer (Fig. 4B), and the seasonal pattern of ammonium increase over time—largest in winter and near-zero in summer (Fig. 6B). Thus, the long-term trends of ammonium increase in the estuary can be attributed to increased wastewater inputs. But the expression of increased loading is not evident in summer (Fig. 6B), implying that wastewater ammonium is nitrified in the time of transit from the wastewater source to the sampling site in the estuary during the warm, low-flow season. Direct measurements of nitrate production and travel time in the lower Sacramento River support this conclusion (Kraus et al. 2017).
Plots
#plot of Ammonia by water year type
ggplot(nutr_chla_data_c, aes(x = YearType, y = DissAmmonia_l, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Ammonia mg/L log transformed")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 43 rows containing non-finite values (`stat_boxplot()`).
#plot of Ammonia by water year type
ggplot(nutr_chla_data_c, aes(x = YearType, y = DissAmmonia, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Ammonia mg/L")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 43 rows containing non-finite values (`stat_boxplot()`).
#plot of Ammonia by Sac valley index
ggplot(nutr_chla_data_c, aes(x = SVIndex, y = DissAmmonia_l))+
geom_point(alpha = 0.7, aes(color =Drought))+
geom_smooth(method = "lm")+
facet_grid(Season~Region)+
color_pal_drought(aes_type = "color")+
ylab("Ammonia mg/L (log-transformed)")+
xlab("Sac Valley Index")+
theme_bw()## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 43 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 43 rows containing missing values (`geom_point()`).
#plot of ammonia by drought
ggplot(nutr_chla_data_c, aes(x = Drought, y = DissAmmonia_l, fill = Drought))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_drought()+
ylab("Ammonia mg/L log transformed")+
xlab("Year Type")+
theme_bw()## Warning: Removed 43 rows containing non-finite values (`stat_boxplot()`).
#Ammonia by year
ggplot(nutr_chla_data_c, aes(x = YearAdj, y = DissAmmonia, fill = YearType))+
geom_col(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Ammonia mg/L")+
theme_bw()## Warning: Removed 43 rows containing missing values (`position_stack()`).
#plot of ammonia by drought
ggplot(nutr_chla_data_c, aes(x = Drought, y = DissAmmonia, fill = Drought))+
geom_boxplot(alpha = 0.7)+
facet_wrap(~Season, scales = "free_y")+
color_pal_drought()+
ylab("Ammonia mg/L")+
xlab("Year Type")+
theme_bw()## Warning: Removed 43 rows containing non-finite values (`stat_boxplot()`).
Run Model
Now we’ll do the statistics. I’m going to take out 2021, because of the Sac WWTP upgrade. Stockton WWTP was upgraded in 2006
Ammonia <- nutr_chla_data_c %>%
select(!starts_with(c("DissNitrateNitrite", "DissOrthophos", "Chlorophyll"))) %>%
drop_na(DissAmmonia_l) %>%
filter(YearAdj < 2021) %>%
mutate(
StocktonUpgrade = case_when(
YearAdj < 2006 ~ "Before",
YearAdj >= 2006 ~ "After"
)
)Now add the treatment plant term
Am2 = lm(DissAmmonia_l ~ Drought*Region + Drought*Season + Region*Season + Region*StocktonUpgrade, data = Ammonia)Check assumptions
model_plotter(Am2, Ammonia, "DissAmmonia_l")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Check results
summary(Am2)##
## Call:
## lm(formula = DissAmmonia_l ~ Drought * Region + Drought * Season +
## Region * Season + Region * StocktonUpgrade, data = Ammonia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7791 -0.2205 0.0233 0.2491 1.0467
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -1.801157 0.080324 -22.424
## DroughtN -0.122208 0.092376 -1.323
## DroughtW -0.343996 0.099713 -3.450
## RegionConfluence 0.036890 0.105230 0.351
## RegionSouth-Central Delta -0.273997 0.105687 -2.593
## RegionNorth Delta -0.156065 0.105687 -1.477
## SeasonSpring -0.512492 0.094812 -5.405
## SeasonSummer -1.150789 0.095056 -12.106
## SeasonFall -0.766192 0.094997 -8.065
## StocktonUpgradeBefore -0.157276 0.065768 -2.391
## DroughtN:RegionConfluence 0.063716 0.098599 0.646
## DroughtW:RegionConfluence 0.002578 0.108454 0.024
## DroughtN:RegionSouth-Central Delta 0.033224 0.100301 0.331
## DroughtW:RegionSouth-Central Delta 0.193461 0.111051 1.742
## DroughtN:RegionNorth Delta -0.204596 0.100301 -2.040
## DroughtW:RegionNorth Delta -0.645709 0.111051 -5.815
## DroughtN:SeasonSpring -0.273769 0.098045 -2.792
## DroughtW:SeasonSpring -0.169713 0.102563 -1.655
## DroughtN:SeasonSummer 0.026420 0.098045 0.269
## DroughtW:SeasonSummer 0.408359 0.103293 3.953
## DroughtN:SeasonFall 0.098777 0.098045 1.007
## DroughtW:SeasonFall 0.337498 0.102563 3.291
## RegionConfluence:SeasonSpring 0.116499 0.115511 1.009
## RegionSouth-Central Delta:SeasonSpring -0.248050 0.117196 -2.117
## RegionNorth Delta:SeasonSpring 0.693275 0.117196 5.915
## RegionConfluence:SeasonSummer 0.118074 0.116166 1.016
## RegionSouth-Central Delta:SeasonSummer -0.535880 0.117517 -4.560
## RegionNorth Delta:SeasonSummer 1.246926 0.117517 10.611
## RegionConfluence:SeasonFall 0.096819 0.115511 0.838
## RegionSouth-Central Delta:SeasonFall -0.309729 0.117514 -2.636
## RegionNorth Delta:SeasonFall 1.302829 0.117514 11.087
## RegionConfluence:StocktonUpgradeBefore -0.102772 0.093058 -1.104
## RegionSouth-Central Delta:StocktonUpgradeBefore 0.943546 0.094415 9.994
## RegionNorth Delta:StocktonUpgradeBefore 0.450659 0.094415 4.773
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## DroughtN 0.186313
## DroughtW 0.000597 ***
## RegionConfluence 0.726024
## RegionSouth-Central Delta 0.009738 **
## RegionNorth Delta 0.140239
## SeasonSpring 9.04e-08 ***
## SeasonSummer < 2e-16 ***
## SeasonFall 3.44e-15 ***
## StocktonUpgradeBefore 0.017065 *
## DroughtN:RegionConfluence 0.518367
## DroughtW:RegionConfluence 0.981045
## DroughtN:RegionSouth-Central Delta 0.740566
## DroughtW:RegionSouth-Central Delta 0.081957 .
## DroughtN:RegionNorth Delta 0.041766 *
## DroughtW:RegionNorth Delta 9.47e-09 ***
## DroughtN:SeasonSpring 0.005385 **
## DroughtW:SeasonSpring 0.098457 .
## DroughtN:SeasonSummer 0.787657
## DroughtW:SeasonSummer 8.54e-05 ***
## DroughtN:SeasonFall 0.314081
## DroughtW:SeasonFall 0.001053 **
## RegionConfluence:SeasonSpring 0.313559
## RegionSouth-Central Delta:SeasonSpring 0.034673 *
## RegionNorth Delta:SeasonSpring 5.31e-09 ***
## RegionConfluence:SeasonSummer 0.309799
## RegionSouth-Central Delta:SeasonSummer 6.10e-06 ***
## RegionNorth Delta:SeasonSummer < 2e-16 ***
## RegionConfluence:SeasonFall 0.402233
## RegionSouth-Central Delta:SeasonFall 0.008594 **
## RegionNorth Delta:SeasonFall < 2e-16 ***
## RegionConfluence:StocktonUpgradeBefore 0.269827
## RegionSouth-Central Delta:StocktonUpgradeBefore < 2e-16 ***
## RegionNorth Delta:StocktonUpgradeBefore 2.23e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3861 on 661 degrees of freedom
## Multiple R-squared: 0.6898, Adjusted R-squared: 0.6743
## F-statistic: 44.54 on 33 and 661 DF, p-value: < 2.2e-16Aman = Anova(Am2, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
Aman## Anova Table (Type III tests)
##
## Response: DissAmmonia_l
## Sum Sq Df F value Pr(>F)
## (Intercept) 74.968 1 502.8181 < 2.2e-16 ***
## Drought 1.779 2 5.9676 0.00270 **
## Region 1.656 3 3.7025 0.01160 *
## Season 23.028 3 51.4846 < 2.2e-16 ***
## StocktonUpgrade 0.853 1 5.7186 0.01706 *
## Drought:Region 9.670 6 10.8096 1.719e-11 ***
## Drought:Season 7.447 6 8.3249 1.031e-08 ***
## Region:Season 46.160 9 34.4004 < 2.2e-16 ***
## Region:StocktonUpgrade 22.852 3 51.0915 < 2.2e-16 ***
## Residuals 98.552 661
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# visreg(Am2, xvar = "Drought", by = "Region")
# visreg(Am2, xvar = "Drought", by = "Season")
# visreg(Am2, xvar = "Region", by = "Season")
# visreg(Am2, xvar = "StocktonUpgrade", by = "Region")
# plot pairwise comparisons
# amp2 = emmeans(Am2, pairwise ~ Drought*Region)
# plot(amp2, comparisons = T)
#
# amp2x = emmeans(Am2, pairwise ~ Drought*Season)
# plot(amp2x, comparisons = T)Partial R2
How much variability is explained by the Drought index?
partial.r2(Aman)## [1] 0.1608936Post-hoc tests and figures
AmRegions = pub_figure_plotter(Ammonia, DissAmmonia, expression(Ammonium~(`mg-N`~L^{-1})), Region, Am2, log_trans = T,
plt_title = "Ammonium")
AmRegions$diffs## # A tibble: 4 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Bay 0.0831 0.0681 0.0151
## 2 Confluence 0.0890 0.0731 0.0159
## 3 South-Central Delta 0.0771 0.0766 0.000499
## 4 North Delta 0.200 0.0861 0.114AmSeasons = pub_figure_plotter(Ammonia, DissAmmonia, expression(Ammonium~(`mg-N`~L^{-1})), Season, Am2, log_trans = T,
plt_title = "Ammonium")
AmSeasons$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 0.163 0.103 0.0596
## 2 Spring 0.112 0.0599 0.0521
## 3 Summer 0.0633 0.0603 0.00297
## 4 Fall 0.0992 0.0881 0.0111# Run Post-hoc test for the treatment plant upgrade by region
amRegime = emmeans(Am2, pairwise ~ StocktonUpgrade|Region)
# Calculate effect sizes and match up p-values
amRegime_diff = as_tibble(amRegime$emmeans) %>%
mutate(trans = exp(emmean)) %>%
select(StocktonUpgrade, Region, trans) %>%
pivot_wider(names_from = StocktonUpgrade, values_from = trans) %>%
mutate(difference = After - Before) %>%
left_join(as_tibble(amRegime$contrasts) %>% select(Region, p.value))
amRegime_diff## # A tibble: 4 × 5
## Region After Before difference p.value
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Suisun Bay 0.0798 0.0682 0.0116 1.71e- 2
## 2 Confluence 0.0920 0.0709 0.0211 8.69e- 5
## 3 South-Central Delta 0.0498 0.109 -0.0595 1.77e-28
## 4 North Delta 0.116 0.155 -0.0394 1.72e- 5# ggplot(Ammonia, aes(x = StocktonUpgrade, y = DissAmmonia)) +
# facet_wrap(~Region) + geom_boxplot()
visreg(Am2, xvar = "StocktonUpgrade", by = "Region")
Well, the treatment plant term improves the model a bit, but it doesn’t actually change the drought story. However including the region*season interaction and taking out 2021 does make this a little more interesting than before. Droughts in the north are definitely higher in Ammonia.
Nitrate + Nitrite
Plots
#plot of Nitrate by water year type
ggplot(nutr_chla_data_c, aes(x = YearType, y = DissNitrateNitrite_l, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Nitrate Nitrite (log-transformed)")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
ggplot(nutr_chla_data_c, aes(x = YearType, y = DissNitrateNitrite, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Nitrate mg/L")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
#plot of Nitrate by Sac valley index
ggplot(nutr_chla_data_c, aes(x = SVIndex, y = DissNitrateNitrite_l))+
geom_point(alpha = 0.7, aes(color =Drought))+
geom_smooth(method = "lm")+
facet_grid(Season~Region)+
color_pal_drought(aes_type = "color")+
ylab("Nitrate nitrite mg/L (log-transformed)")+
xlab("Sac Valley Index")+
theme_bw()## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 4 rows containing missing values (`geom_point()`).
#plot ofNitrate by drought
ggplot(nutr_chla_data_c, aes(x = Drought, y = DissNitrateNitrite, fill = Drought))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_drought()+
ylab("Nitrate mg/L")+
xlab("Year Type")+
theme_bw()## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
#Nitrate by year
ggplot(nutr_chla_data_c, aes(x = YearAdj, y = DissNitrateNitrite, fill = YearType))+
geom_col(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Nitrate mg/L")+
theme_bw()## Warning: Removed 4 rows containing missing values (`position_stack()`).
Run Model
OK, now for the statistics. The treatment plant upgrades didn’t have as big an impact on nitrate, so I’ll leave it out for now, but I can put it back in if folks think it’s important.
Nitrate <- nutr_chla_data_c %>%
select(!starts_with(c("DissAmmonia", "DissOrthophos", "Chlorophyll"))) %>%
drop_na(DissNitrateNitrite_l)Nat2 = lm(DissNitrateNitrite_l ~ Drought*Region + Drought*Season + Region*Season, data = Nitrate)Check assumptions
model_plotter(Nat2, Nitrate, "DissNitrateNitrite_l")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Check results
summary(Nat2)##
## Call:
## lm(formula = DissNitrateNitrite_l ~ Drought * Region + Drought *
## Season + Region * Season, data = Nitrate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.36705 -0.18837 0.00195 0.18596 0.80554
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.75621 0.05398 -14.008 < 2e-16 ***
## DroughtN -0.14685 0.06718 -2.186 0.029143 *
## DroughtW -0.24473 0.07075 -3.459 0.000574 ***
## RegionConfluence -0.02516 0.07012 -0.359 0.719852
## RegionSouth-Central Delta 0.68025 0.06980 9.746 < 2e-16 ***
## RegionNorth Delta -0.83261 0.06980 -11.929 < 2e-16 ***
## SeasonSpring -0.04045 0.06970 -0.580 0.561836
## SeasonSummer -0.23860 0.06970 -3.423 0.000654 ***
## SeasonFall -0.16999 0.06970 -2.439 0.014974 *
## DroughtN:RegionConfluence 0.07389 0.07163 1.032 0.302612
## DroughtW:RegionConfluence 0.01110 0.07511 0.148 0.882545
## DroughtN:RegionSouth-Central Delta 0.17729 0.07154 2.478 0.013436 *
## DroughtW:RegionSouth-Central Delta 0.28100 0.07528 3.733 0.000205 ***
## DroughtN:RegionNorth Delta 0.31068 0.07154 4.343 1.61e-05 ***
## DroughtW:RegionNorth Delta 0.31236 0.07528 4.149 3.74e-05 ***
## DroughtN:SeasonSpring -0.27050 0.07163 -3.776 0.000172 ***
## DroughtW:SeasonSpring -0.29209 0.07564 -3.862 0.000123 ***
## DroughtN:SeasonSummer -0.25496 0.07163 -3.559 0.000396 ***
## DroughtW:SeasonSummer -0.13317 0.07564 -1.761 0.078731 .
## DroughtN:SeasonFall -0.24882 0.07163 -3.474 0.000544 ***
## DroughtW:SeasonFall -0.14134 0.07564 -1.869 0.062092 .
## RegionConfluence:SeasonSpring -0.11730 0.08590 -1.365 0.172525
## RegionSouth-Central Delta:SeasonSpring -0.23353 0.08591 -2.718 0.006723 **
## RegionNorth Delta:SeasonSpring -0.25038 0.08591 -2.914 0.003676 **
## RegionConfluence:SeasonSummer -0.17448 0.08590 -2.031 0.042612 *
## RegionSouth-Central Delta:SeasonSummer -0.57945 0.08591 -6.745 3.16e-11 ***
## RegionNorth Delta:SeasonSummer -0.62434 0.08591 -7.267 9.60e-13 ***
## RegionConfluence:SeasonFall -0.06774 0.08590 -0.789 0.430608
## RegionSouth-Central Delta:SeasonFall -0.29555 0.08591 -3.440 0.000615 ***
## RegionNorth Delta:SeasonFall -0.23123 0.08591 -2.691 0.007280 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2929 on 718 degrees of freedom
## Multiple R-squared: 0.8159, Adjusted R-squared: 0.8085
## F-statistic: 109.8 on 29 and 718 DF, p-value: < 2.2e-16NatAn =Anova(Nat2, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
NatAn## Anova Table (Type III tests)
##
## Response: DissNitrateNitrite_l
## Sum Sq Df F value Pr(>F)
## (Intercept) 16.830 1 196.2212 < 2.2e-16 ***
## Drought 1.099 2 6.4094 0.0017418 **
## Region 40.840 3 158.7157 < 2.2e-16 ***
## Season 1.315 3 5.1120 0.0016651 **
## Drought:Region 3.374 6 6.5566 9.357e-07 ***
## Drought:Season 2.380 6 4.6255 0.0001269 ***
## Region:Season 6.875 9 8.9057 8.543e-13 ***
## Residuals 61.584 718
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# visreg(Nat2, xvar = "Drought", by = "Region")
# visreg(Nat2, xvar = "Drought", by = "Season")
# visreg(Nat2, xvar = "Region", by = "Season")
# Pairwise comparisons
# pairsNN = emmeans(Nat2, pairwise ~ Drought*Region)
# plot(pairsNN, comparisons = T)
# pairsNN2 = emmeans(Nat2, pairwise ~ Drought*Season)
# plot(pairsNN2, comparisons = T)Partial R2
How much variability is explained by the Drought index?
partial.r2(NatAn)## [1] 0.1001508Post-hoc tests and figures
NNRegion = pub_figure_plotter(Nitrate, DissNitrateNitrite, expression(Nitrate+Nitrite~(`mg-N`~L^{-1})), Region, Nat2, log_trans = T, plt_title = "Nitrate + Nitrite")
NNRegion$diffs## # A tibble: 4 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Bay 0.420 0.285 0.134
## 2 Confluence 0.374 0.257 0.117
## 3 South-Central Delta 0.628 0.565 0.0628
## 4 North Delta 0.138 0.129 0.00987NNSeason = pub_figure_plotter(Nitrate, DissNitrateNitrite, expression(Nitrate+Nitrite~(`mg-N`~L^{-1})), Season, Nat2, log_trans = T, plt_title = "Nitrate + Nitrite")
NNSeason$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 0.449 0.409 0.0401
## 2 Spring 0.371 0.252 0.119
## 3 Summer 0.251 0.200 0.0509
## 4 Fall 0.327 0.258 0.0684OK, there is a big main effect of drought here, with droughts having higher nitrate, especially in the spring, summer and fall (no effect in the winter). If I’m reading this right, droughts mean increased nitrate in Suisun Bay and the Confluence, but not the North, whereas Ammonia increased in the North. Interesting.
Ortho-phosphate
According to Cloern, total phosphorus loading from the WWTP decreased a lot in the 80s, but there was no real phosphorus trend over time.
Phosphate Although phosphate concentration did not change over the long term (Table 2), there was substantial variability of phosphate concentration as an early period of increase during the 1980s, peak concentrations in the early 1990s (Fig. 3) followed by a period of significant decrease (Fig. 5B). Records of phosphorus loading from the Sacramento Regional Wastewater Treatment Facility are not available during the early era of phosphate increase. However, measurements since 1987 showed a steep decrease in wastewater TP loading in the early 1990s that paralleled the downstream phosphate decrease at Sta. D8. This suggests that wastewater effluent is a significant source of phosphorus, as well as ammonium nitrogen, to the upper estuary. The best-fitting regression model (Table 4) explained 73% of annual-mean phosphate concentration as a negative function of Outflow (effect size = −0.74) and a positive function of wastewater TP loading (effect size = +0.25), similar to model results for ammonium. Relative to the Outflow effect, the effect of wastewater TP loading was smaller than the effect of ammonium loading (Table 4). This suggests that the wastewater enrichment of riverine phosphate is smaller than its enrichment of riverine ammonium. The parallel steep declines of phosphorus loading and concentrations in the estuary during the 1990s (Fig. 3) occurred after it became evident that phosphorus enrichment was degrading water quality of U.S. lakes and rivers. State bans on phosphorus-containing detergents grew in the 1970s and 1980s, and by 1994 the manufacture of phosphorus-based detergents ended (Litke 1999). Thus, the oscillating patterns of phosphate concentrations in San Francisco Bay (and other estuaries and lakes) include eras of increase associated with population growth (Jassby 2008) and decrease after policies were implemented to remove phosphorus from household detergents.
Plots
#plot of Phosphorus by water year type
ggplot(nutr_chla_data_c, aes(x = YearType, y = DissOrthophos_l, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Orthophosphate mg/L (log-transformed)")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
ggplot(nutr_chla_data_c, aes(x = YearType, y = DissOrthophos, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Orthophosphate mg/L")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
#plot of phosphorus by Sac valley index
ggplot(nutr_chla_data_c, aes(x = SVIndex, y = DissOrthophos_l))+
geom_point(alpha = 0.7, aes(color =Drought))+
geom_smooth(method = "lm")+
facet_grid(Season~Region)+
color_pal_drought(aes_type = "color")+
ylab("Orthophosphate mg/L (log-transformed)")+
xlab("Sac Valley Index")+
theme_bw()## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 2 rows containing missing values (`geom_point()`).
ggplot(nutr_chla_data_c, aes(x = YearAdj, y = DissOrthophos)) + geom_point(aes(color = Region))+
geom_smooth(aes(color = Region))+
theme_bw()+ ylab("Dissolved Orthophosphate")## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite values (`stat_smooth()`).
## Removed 2 rows containing missing values (`geom_point()`).
#Wow. Phosphorus actually looks pretty good there.
ggplot(nutr_chla_data_c, aes(x = Drought, y = DissOrthophos, fill = Drought))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_drought()+
ylab("Phospohrus mg/L")+
xlab("Year Type")+
theme_bw()## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
#Phosphorus by year
ggplot(nutr_chla_data_c, aes(x = YearAdj, y = DissOrthophos, fill = YearType))+
geom_col(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Phosphorus mg/L")+
theme_bw()## Warning: Removed 2 rows containing missing values (`position_stack()`).
Huh. Cloern didn’t see any overall trends in Phosphorus, but I think that’s a pattern. Anyway.
Run Model
Log-transforming phosphorus didn’t look quite as nice as nitrogen and ammonium. Square-root transformation looked better, but the residuals actually looked pretty good for both, so I just stuck with the log-transformation to keep it consistent.
Phos <- nutr_chla_data_c %>%
select(!starts_with(c("DissAmmonia", "DissNitrateNitrite", "Chlorophyll"))) %>%
drop_na(DissOrthophos_l)Phos2 = lm(DissOrthophos_l ~ Drought*Region + Drought*Season + Region*Season, data = Phos)Check assumptions
# Check assumptions
model_plotter(Phos2, Phos, "DissOrthophos_l")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Check results
# Check results
summary(Phos2)##
## Call:
## lm(formula = DissOrthophos_l ~ Drought * Region + Drought * Season +
## Region * Season, data = Phos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02192 -0.20210 -0.00486 0.20653 0.92993
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.57751 0.05623 -45.840 < 2e-16 ***
## DroughtN -0.03398 0.07001 -0.485 0.627579
## DroughtW -0.26424 0.07339 -3.600 0.000340 ***
## RegionConfluence -0.06539 0.07307 -0.895 0.371149
## RegionSouth-Central Delta 0.40598 0.07256 5.595 3.14e-08 ***
## RegionNorth Delta -0.37363 0.07256 -5.149 3.38e-07 ***
## SeasonSpring 0.01983 0.07259 0.273 0.784845
## SeasonSummer 0.20571 0.07259 2.834 0.004729 **
## SeasonFall 0.30510 0.07259 4.203 2.97e-05 ***
## DroughtN:RegionConfluence 0.02206 0.07464 0.296 0.767648
## DroughtW:RegionConfluence -0.03546 0.07827 -0.453 0.650684
## DroughtN:RegionSouth-Central Delta 0.14423 0.07455 1.935 0.053431 .
## DroughtW:RegionSouth-Central Delta 0.19197 0.07818 2.455 0.014311 *
## DroughtN:RegionNorth Delta 0.10988 0.07455 1.474 0.140956
## DroughtW:RegionNorth Delta -0.08961 0.07818 -1.146 0.252122
## DroughtN:SeasonSpring -0.26884 0.07465 -3.602 0.000338 ***
## DroughtW:SeasonSpring -0.23109 0.07827 -2.952 0.003256 **
## DroughtN:SeasonSummer -0.18959 0.07465 -2.540 0.011299 *
## DroughtW:SeasonSummer -0.14533 0.07827 -1.857 0.063759 .
## DroughtN:SeasonFall -0.16136 0.07465 -2.162 0.030978 *
## DroughtW:SeasonFall -0.14173 0.07827 -1.811 0.070606 .
## RegionConfluence:SeasonSpring -0.05640 0.08952 -0.630 0.528910
## RegionSouth-Central Delta:SeasonSpring -0.13917 0.08928 -1.559 0.119502
## RegionNorth Delta:SeasonSpring -0.09148 0.08928 -1.025 0.305885
## RegionConfluence:SeasonSummer -0.13567 0.08952 -1.515 0.130088
## RegionSouth-Central Delta:SeasonSummer -0.42110 0.08928 -4.717 2.88e-06 ***
## RegionNorth Delta:SeasonSummer -0.15276 0.08928 -1.711 0.087513 .
## RegionConfluence:SeasonFall -0.11522 0.08952 -1.287 0.198510
## RegionSouth-Central Delta:SeasonFall -0.57882 0.08928 -6.483 1.67e-10 ***
## RegionNorth Delta:SeasonFall 0.04688 0.08928 0.525 0.599716
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3052 on 720 degrees of freedom
## Multiple R-squared: 0.5239, Adjusted R-squared: 0.5047
## F-statistic: 27.32 on 29 and 720 DF, p-value: < 2.2e-16Panova = Anova(Phos2, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
Panova## Anova Table (Type III tests)
##
## Response: DissOrthophos_l
## Sum Sq Df F value Pr(>F)
## (Intercept) 195.740 1 2101.3318 < 2.2e-16 ***
## Drought 1.298 2 6.9676 0.001007 **
## Region 11.066 3 39.5996 < 2.2e-16 ***
## Season 2.329 3 8.3326 1.880e-05 ***
## Drought:Region 1.814 6 3.2459 0.003731 **
## Drought:Season 1.550 6 2.7733 0.011324 *
## Region:Season 7.531 9 8.9827 6.395e-13 ***
## Residuals 67.068 720
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# visreg(Phos2, xvar = "Drought", by = "Region")
# visreg(Phos2, xvar = "Drought", by = "Season")
# visreg(Phos2, xvar = "Region", by = "Season")
# pairsP = emmeans(Phos2, pairwise ~ Drought|Region)
# plot(pairsP, comparisons = T)
# plot(emmeans(Phos2, pairwise ~ Drought), comparisons = T)Partial R2
How much variability is explained by the Drought index?
# Calculate partial R2 for drought effect
partial.r2(Panova)## [1] 0.06499649Post-hoc tests and figures
OPRegion = pub_figure_plotter(Phos, DissOrthophos, expression(`Ortho-phosphate`~(`mg-P`~L^{-1})), Region, Phos2, log_trans = T, plt_title = "Ortho-phosphate")
OPRegion$diffs## # A tibble: 4 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Bay 0.0867 0.0585 0.0282
## 2 Confluence 0.0752 0.0490 0.0263
## 3 South-Central Delta 0.0979 0.0800 0.0179
## 4 North Delta 0.0568 0.0350 0.0218OPSeason = pub_figure_plotter(Phos, DissOrthophos, expression(`Ortho-phosphate`~(`mg-P`~L^{-1})), Season, Phos2, log_trans = T, plt_title = "Ortho-phosphate")
OPSeason$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 0.0753 0.0588 0.0165
## 2 Spring 0.0715 0.0443 0.0272
## 3 Summer 0.0775 0.0523 0.0252
## 4 Fall 0.0869 0.0589 0.0280For phosphorus, the story is still really the main effect of drought. Lower Phosphorus in wet years.
Run Model with SV Index
We didn’t use this for the manuscript.
#now again as a linear model versus sac valley index
Phos2SV = lm(DissOrthophos_l ~ SVIndex*Region + SVIndex*Season + Region*Season, data = Phos)
# plot(Phos2SV)
visreg(Phos2SV, xvar = "SVIndex", by = "Region")
visreg(Phos2SV, xvar = "SVIndex", by = "Season")
visreg(Phos2SV, xvar = "Region", by = "Season")
summary(Phos2SV)##
## Call:
## lm(formula = DissOrthophos_l ~ SVIndex * Region + SVIndex * Season +
## Region * Season, data = Phos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19277 -0.18359 -0.00559 0.18146 0.79352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.335747 0.082874 -28.184 < 2e-16
## SVIndex -0.040839 0.009127 -4.474 8.89e-06
## RegionConfluence -0.042069 0.096245 -0.437 0.662168
## RegionSouth-Central Delta 0.241509 0.095586 2.527 0.011728
## RegionNorth Delta -0.249407 0.095586 -2.609 0.009261
## SeasonSpring 0.209556 0.095733 2.189 0.028916
## SeasonSummer 0.370184 0.095733 3.867 0.000120
## SeasonFall 0.416756 0.095733 4.353 1.53e-05
## SVIndex:RegionConfluence -0.003289 0.009707 -0.339 0.734808
## SVIndex:RegionSouth-Central Delta 0.032748 0.009684 3.382 0.000759
## SVIndex:RegionNorth Delta -0.014948 0.009684 -1.544 0.123122
## SVIndex:SeasonSpring -0.042530 0.009707 -4.381 1.35e-05
## SVIndex:SeasonSummer -0.033428 0.009707 -3.444 0.000607
## SVIndex:SeasonFall -0.025434 0.009707 -2.620 0.008973
## RegionConfluence:SeasonSpring -0.056628 0.082516 -0.686 0.492763
## RegionSouth-Central Delta:SeasonSpring -0.137324 0.082293 -1.669 0.095604
## RegionNorth Delta:SeasonSpring -0.089638 0.082293 -1.089 0.276401
## RegionConfluence:SeasonSummer -0.135902 0.082516 -1.647 0.099994
## RegionSouth-Central Delta:SeasonSummer -0.419258 0.082293 -5.095 4.46e-07
## RegionNorth Delta:SeasonSummer -0.150919 0.082293 -1.834 0.067075
## RegionConfluence:SeasonFall -0.115446 0.082516 -1.399 0.162214
## RegionSouth-Central Delta:SeasonFall -0.576979 0.082293 -7.011 5.40e-12
## RegionNorth Delta:SeasonFall 0.048719 0.082293 0.592 0.554023
##
## (Intercept) ***
## SVIndex ***
## RegionConfluence
## RegionSouth-Central Delta *
## RegionNorth Delta **
## SeasonSpring *
## SeasonSummer ***
## SeasonFall ***
## SVIndex:RegionConfluence
## SVIndex:RegionSouth-Central Delta ***
## SVIndex:RegionNorth Delta
## SVIndex:SeasonSpring ***
## SVIndex:SeasonSummer ***
## SVIndex:SeasonFall **
## RegionConfluence:SeasonSpring
## RegionSouth-Central Delta:SeasonSpring .
## RegionNorth Delta:SeasonSpring
## RegionConfluence:SeasonSummer .
## RegionSouth-Central Delta:SeasonSummer ***
## RegionNorth Delta:SeasonSummer .
## RegionConfluence:SeasonFall
## RegionSouth-Central Delta:SeasonFall ***
## RegionNorth Delta:SeasonFall
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2813 on 727 degrees of freedom
## Multiple R-squared: 0.5916, Adjusted R-squared: 0.5792
## F-statistic: 47.86 on 22 and 727 DF, p-value: < 2.2e-16Chlorophyll
This is what Cloern has to say:
Phytoplankton biomass in estuaries is regulated by river driven transport processes and the balance between growth and grazing rates. Chl a at Sta. D8 decreased abruptly in 1987 (Table 3), and the largest declines occurred in summer—May through September (Fig. 6C). The best-fitting model of summer Chl a had an interactive effect between Outflow and clam abundance. This interaction was expected because the relationship between Chl a and Outflow changed after the Potamocorbula introduction (Jassby 2008). Most (80%) of that variability was associated with Outflow and clam abundance treated as a categorical variable (Fig. 8B). No other variables, such as temperature, turbidity (Secchi depth), nutrient concentrations, or loadings had comparable effects on summer Chl a.
Summer phytoplankton blooms occurred regularly in this region of the estuary during the 1970s and 1980s. These were attributed to the accumulation of phytoplankton cells within an estuarine turbidity maximum that was positioned in Suisun Bay (Sta. D8) when Outflow was in the range 100–350 m3 s−1 (Cloern et al. 1983). Thus, flow-phytoplankton relationships were established early in the study of San Francisco Bay. Observations made over the subsequent 3+ decades confirm that summer Chl a decreases when Outflow exceeds about200–300 m3 s−1 (Fig. 8B), primarily through its effect on residence time (Jassby 2008). A summer bloom did not develop in 1977, a year of extreme drought when salinity increased high enough that the upper estuary was colonized by marine benthic filter feeders including M. arenaria (Nichols 1985). The 1977 anomaly was the first evidence that bivalve grazing can suppress blooms in this estuary. Summer blooms disappeared completely after P. amurensis became established as a permanent resident in 1987. The regression model (Table 4) shows that the effect of clam presence was large (−13.1) compared to the outflow effect (−2.5). The nearly fourfold loss of phytoplankton biomass (Table 2) and production restructured biological communities and altered nutrient dynamics, establishing what has become an iconic example of abrupt disturbance by the introduction of a non-native species (Cloern and Jassby 2012). Regression models allow us to compare the effects of that human disturbance with the effect of varying freshwater inflow.
Plots
#plot of chlorophyll by water year type
ggplot(nutr_chla_data_c, aes(x = YearType, y = Chlorophyll_l, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Chlorophyl ug/L (log-transformed)")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
#plot of chlorophyll by Sac valley index
ggplot(nutr_chla_data_c, aes(x = SVIndex, y = Chlorophyll_l))+
geom_point(alpha = 0.7, aes(color =Drought))+
geom_smooth(method = "lm")+
facet_grid(Season~Region)+
color_pal_drought(aes_type = "color")+
ylab("Chlorophyl ug/L (log-transformed)")+
xlab("Sac Valley Index")+
theme_bw()## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
#plot of chlorophyll by drought
ggplot(nutr_chla_data_c, aes(x = Drought, y = Chlorophyll_l, fill = Drought))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_drought()+
ylab("Chlorophyl ug/L (log-transformed)")+
xlab("Year Type")+
theme_bw()## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
ggplot(nutr_chla_data_c, aes(x = Drought, y = Chlorophyll, fill = Drought))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_drought()+
ylab("Chlorophyl ug/L")+
xlab("Year Type")+
theme_bw()## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
#plot of chlorophyll by water year type - just the more recent years
ggplot(filter(nutr_chla_data_c, YearAdj > 1989), aes(x = YearType, y = Chlorophyll_l, fill = YearType))+
geom_boxplot(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Chlorophyl ug/L (log-transformed)")+
xlab("Year Type")+
theme_bw()+
scale_x_discrete(labels = c("C", "D", "B", "A", "W"))## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
#Bleh, nothing going on there
#plot of chlorophyll by year
ggplot(nutr_chla_data_c, aes(x = YearAdj, y = Chlorophyll_l, fill = YearType))+
geom_col(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_yrtype()+
ylab("Chlorophyl ug/L (log-transformed)")+
theme_bw()## Warning: Removed 1 rows containing missing values (`position_stack()`).
ggplot(nutr_chla_data_c, aes(x = YearAdj, y = Chlorophyll_l, fill = Drought))+
geom_col(alpha = 0.7)+
facet_grid(Season~Region)+
color_pal_drought()+
ylab("Chlorophyl ug/L (log-transformed)")+
theme_bw()## Warning: Removed 1 rows containing missing values (`position_stack()`).
ggplot(nutr_chla_data_c, aes(x = YearAdj, y = Chlorophyll_l, color = Region))+
geom_point(alpha = 0.7)+
geom_smooth()+
# color_pal_drought()+
ylab("Chlorophyl ug/L (log-transformed)")+
theme_bw()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Removed 1 rows containing missing values (`geom_point()`).
Run Model
And now for some statistics
ChlaA <- nutr_chla_data_c %>%
select(!starts_with(c("DissAmmonia", "DissNitrateNitrite", "DissOrthophos"))) %>%
drop_na(Chlorophyll_l) %>%
mutate(
Regime = case_when(
YearAdj > 1987 ~ "post-clam",
YearAdj < 1988 ~ "pre-clam"
)
)# model including pre/post clam regime
Chl1 = lm(Chlorophyll_l ~ Drought*Region + Drought*Season + Region*Season + Region*Regime, data = ChlaA)Check assumptions
model_plotter(Chl1, ChlaA, "Chlorophyll_l")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
qqPlot(Chl1)
## [1] 38 727Check results
summary(Chl1)##
## Call:
## lm(formula = Chlorophyll_l ~ Drought * Region + Drought * Season +
## Region * Season + Region * Regime, data = ChlaA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.95802 -0.31743 -0.02147 0.28326 1.93480
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.26722 0.09340 2.861 0.004347
## DroughtN -0.10583 0.11769 -0.899 0.368807
## DroughtW -0.25845 0.12182 -2.122 0.034212
## RegionConfluence 0.13337 0.12076 1.104 0.269783
## RegionSouth-Central Delta 1.27937 0.12076 10.594 < 2e-16
## RegionNorth Delta 0.46710 0.12076 3.868 0.000120
## SeasonSpring 0.75421 0.11948 6.312 4.82e-10
## SeasonSummer 0.64672 0.11948 5.413 8.48e-08
## SeasonFall 0.23890 0.11948 1.999 0.045938
## Regimepre-clam 0.93698 0.08538 10.974 < 2e-16
## DroughtN:RegionConfluence -0.31701 0.12749 -2.487 0.013124
## DroughtW:RegionConfluence -0.26037 0.13076 -1.991 0.046831
## DroughtN:RegionSouth-Central Delta -0.61212 0.12749 -4.801 1.92e-06
## DroughtW:RegionSouth-Central Delta -0.62607 0.13076 -4.788 2.05e-06
## DroughtN:RegionNorth Delta -0.15155 0.12749 -1.189 0.234935
## DroughtW:RegionNorth Delta -0.03289 0.13076 -0.252 0.801467
## DroughtN:SeasonSpring 0.22364 0.12276 1.822 0.068913
## DroughtW:SeasonSpring 0.22809 0.12874 1.772 0.076869
## DroughtN:SeasonSummer 0.46433 0.12276 3.782 0.000168
## DroughtW:SeasonSummer 0.36279 0.12874 2.818 0.004965
## DroughtN:SeasonFall 0.40148 0.12276 3.270 0.001125
## DroughtW:SeasonFall 0.53279 0.12874 4.139 3.91e-05
## RegionConfluence:SeasonSpring 0.22717 0.14701 1.545 0.122729
## RegionSouth-Central Delta:SeasonSpring 0.11638 0.14701 0.792 0.428847
## RegionNorth Delta:SeasonSpring -0.50248 0.14701 -3.418 0.000667
## RegionConfluence:SeasonSummer 0.12858 0.14701 0.875 0.382091
## RegionSouth-Central Delta:SeasonSummer 0.51493 0.14701 3.503 0.000489
## RegionNorth Delta:SeasonSummer -0.83737 0.14701 -5.696 1.79e-08
## RegionConfluence:SeasonFall 0.18224 0.14701 1.240 0.215514
## RegionSouth-Central Delta:SeasonFall 0.22759 0.14701 1.548 0.122043
## RegionNorth Delta:SeasonFall -0.88572 0.14701 -6.025 2.71e-09
## RegionConfluence:Regimepre-clam -0.23619 0.12073 -1.956 0.050809
## RegionSouth-Central Delta:Regimepre-clam -0.56547 0.12073 -4.684 3.37e-06
## RegionNorth Delta:Regimepre-clam -0.70011 0.12073 -5.799 1.00e-08
##
## (Intercept) **
## DroughtN
## DroughtW *
## RegionConfluence
## RegionSouth-Central Delta ***
## RegionNorth Delta ***
## SeasonSpring ***
## SeasonSummer ***
## SeasonFall *
## Regimepre-clam ***
## DroughtN:RegionConfluence *
## DroughtW:RegionConfluence *
## DroughtN:RegionSouth-Central Delta ***
## DroughtW:RegionSouth-Central Delta ***
## DroughtN:RegionNorth Delta
## DroughtW:RegionNorth Delta
## DroughtN:SeasonSpring .
## DroughtW:SeasonSpring .
## DroughtN:SeasonSummer ***
## DroughtW:SeasonSummer **
## DroughtN:SeasonFall **
## DroughtW:SeasonFall ***
## RegionConfluence:SeasonSpring
## RegionSouth-Central Delta:SeasonSpring
## RegionNorth Delta:SeasonSpring ***
## RegionConfluence:SeasonSummer
## RegionSouth-Central Delta:SeasonSummer ***
## RegionNorth Delta:SeasonSummer ***
## RegionConfluence:SeasonFall
## RegionSouth-Central Delta:SeasonFall
## RegionNorth Delta:SeasonFall ***
## RegionConfluence:Regimepre-clam .
## RegionSouth-Central Delta:Regimepre-clam ***
## RegionNorth Delta:Regimepre-clam ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5026 on 717 degrees of freedom
## Multiple R-squared: 0.6796, Adjusted R-squared: 0.6648
## F-statistic: 46.08 on 33 and 717 DF, p-value: < 2.2e-16Ch1_Anova<-Anova(Chl1, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
Ch1_Anova## Anova Table (Type III tests)
##
## Response: Chlorophyll_l
## Sum Sq Df F value Pr(>F)
## (Intercept) 2.067 1 8.1848 0.0043472 **
## Drought 1.138 2 2.2519 0.1059389
## Region 34.549 3 45.5978 < 2.2e-16 ***
## Season 13.177 3 17.3915 6.571e-11 ***
## Regime 30.416 1 120.4313 < 2.2e-16 ***
## Drought:Region 9.933 6 6.5551 9.399e-07 ***
## Drought:Season 6.877 6 4.5384 0.0001578 ***
## Region:Season 30.337 9 13.3465 < 2.2e-16 ***
## Region:Regime 10.462 3 13.8075 9.212e-09 ***
## Residuals 181.086 717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# visreg(Chl1, xvar = "Drought", by = "Region")
# visreg(Chl1, xvar = "Drought", by = "Season")
# visreg(Chl1, xvar = "Region", by = "Season")
# pairsC = emmeans(Chl1, pairwise ~ Drought*Region)
# plot(pairsC, comparisons = T)
# pairsC2 = emmeans(Chl1, pairwise ~ Drought*Season)
# plot(pairsC2, comparisons = T)
#
# pairsC3 = emmeans(Chl1, pairwise ~Regime*Region)
# plot(pairsC3, comparisons = T)Partial R2
How much variability is explained by the Drought index?
partial.r2(Ch1_Anova)## [1] 0.09017665Post-hoc tests and figures
ChRegion = pub_figure_plotter(ChlaA, Chlorophyll, expression(Chlorophyll~(mu*g~L^{-1})), Region, Chl1, log_trans = T, plt_title = "Chlorophyll")
ChRegion$diffs## # A tibble: 4 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Bay 3.14 3.22 -0.0715
## 2 Confluence 3.65 2.88 0.773
## 3 South-Central Delta 10.6 5.77 4.79
## 4 North Delta 2.03 2.01 0.0210CHSeason = pub_figure_plotter(ChlaA, Chlorophyll, expression(Chlorophyll~(mu*g~L^{-1})), Season, Chl1, log_trans = T, plt_title = "Chlorophyll")
CHSeason$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 2.77 1.70 1.07
## 2 Spring 5.65 4.36 1.30
## 3 Summer 5.03 4.44 0.594
## 4 Fall 3.12 3.26 -0.142# Run Post-hoc test for Regime by region
emRegime = emmeans(Chl1, pairwise ~ Regime|Region)
# Calculate effect sizes and match up p-values
emRegime_diff = as_tibble(emRegime$emmeans) %>%
mutate(trans = exp(emmean)) %>%
select(Regime, Region, trans) %>%
pivot_wider(names_from = Regime, values_from = trans) %>%
mutate(difference = `post-clam` - `pre-clam`) %>%
left_join(as_tibble(emRegime$contrasts) %>% select(Region, p.value))
emRegime_diff## # A tibble: 4 × 5
## Region `post-clam` `pre-clam` difference p.value
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Suisun Bay 2.10 5.35 -3.25 5.22e-26
## 2 Confluence 2.26 4.56 -2.30 1.03e-15
## 3 South-Central Delta 6.18 8.96 -2.78 1.54e- 5
## 4 North Delta 1.80 2.28 -0.482 5.66e- 3visreg(Chl1, xvar = "Regime", by = "Region")
Well, I added the regime shift at the clam invasion and now my diagnostic plots look better, droughts definitely have higher chlorophyll in the south-central region, and definitely have higher winter and spring chlorophyll YAY!!
Run Alternate Models
Try the linear model with Sac valley index. We didn’t use this for the manuscript.
Chl2 = lm(Chlorophyll_l ~ SVIndex*Region + SVIndex*Season + Region*Season + Region*Regime, data = ChlaA)
summary(Chl2)##
## Call:
## lm(formula = Chlorophyll_l ~ SVIndex * Region + SVIndex * Season +
## Region * Season + Region * Regime, data = ChlaA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.83126 -0.31906 -0.04021 0.28741 1.89132
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29957 0.14977 2.000 0.045861
## SVIndex -0.01783 0.01649 -1.081 0.279921
## RegionConfluence 0.21763 0.17316 1.257 0.209228
## RegionSouth-Central Delta 1.72978 0.17316 9.989 < 2e-16
## RegionNorth Delta 0.51365 0.17316 2.966 0.003114
## SeasonSpring 0.76352 0.17236 4.430 1.09e-05
## SeasonSummer 0.49602 0.17236 2.878 0.004122
## SeasonFall 0.03521 0.17236 0.204 0.838177
## Regimepre-clam 0.96605 0.08336 11.589 < 2e-16
## SVIndex:RegionConfluence -0.02904 0.01758 -1.651 0.099083
## SVIndex:RegionSouth-Central Delta -0.09812 0.01758 -5.580 3.41e-08
## SVIndex:RegionNorth Delta -0.01180 0.01758 -0.671 0.502354
## SVIndex:SeasonSpring 0.01496 0.01746 0.857 0.391933
## SVIndex:SeasonSummer 0.04918 0.01746 2.817 0.004984
## SVIndex:SeasonFall 0.05917 0.01746 3.389 0.000740
## RegionConfluence:SeasonSpring 0.22642 0.14839 1.526 0.127487
## RegionSouth-Central Delta:SeasonSpring 0.11562 0.14839 0.779 0.436118
## RegionNorth Delta:SeasonSpring -0.50324 0.14839 -3.391 0.000733
## RegionConfluence:SeasonSummer 0.12782 0.14839 0.861 0.389303
## RegionSouth-Central Delta:SeasonSummer 0.51417 0.14839 3.465 0.000561
## RegionNorth Delta:SeasonSummer -0.83812 0.14839 -5.648 2.33e-08
## RegionConfluence:SeasonFall 0.18149 0.14839 1.223 0.221701
## RegionSouth-Central Delta:SeasonFall 0.22683 0.14839 1.529 0.126787
## RegionNorth Delta:SeasonFall -0.88648 0.14839 -5.974 3.63e-09
## RegionConfluence:Regimepre-clam -0.30359 0.11785 -2.576 0.010192
## RegionSouth-Central Delta:Regimepre-clam -0.67036 0.11785 -5.688 1.86e-08
## RegionNorth Delta:Regimepre-clam -0.72751 0.11785 -6.173 1.11e-09
##
## (Intercept) *
## SVIndex
## RegionConfluence
## RegionSouth-Central Delta ***
## RegionNorth Delta **
## SeasonSpring ***
## SeasonSummer **
## SeasonFall
## Regimepre-clam ***
## SVIndex:RegionConfluence .
## SVIndex:RegionSouth-Central Delta ***
## SVIndex:RegionNorth Delta
## SVIndex:SeasonSpring
## SVIndex:SeasonSummer **
## SVIndex:SeasonFall ***
## RegionConfluence:SeasonSpring
## RegionSouth-Central Delta:SeasonSpring
## RegionNorth Delta:SeasonSpring ***
## RegionConfluence:SeasonSummer
## RegionSouth-Central Delta:SeasonSummer ***
## RegionNorth Delta:SeasonSummer ***
## RegionConfluence:SeasonFall
## RegionSouth-Central Delta:SeasonFall
## RegionNorth Delta:SeasonFall ***
## RegionConfluence:Regimepre-clam *
## RegionSouth-Central Delta:Regimepre-clam ***
## RegionNorth Delta:Regimepre-clam ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5072 on 724 degrees of freedom
## Multiple R-squared: 0.6704, Adjusted R-squared: 0.6585
## F-statistic: 56.63 on 26 and 724 DF, p-value: < 2.2e-16# plot(Chl2)
visreg(Chl2, xvar ="SVIndex", by = "Region")
visreg(Chl2, xvar = "SVIndex", by = "Season")
visreg(Chl2, xvar = "Region", by = "Season")
The rest of this was an attempt to use Bayesian stats but I decided not to go there after all.
library(brms)
chlb = brm(Chlorophyll ~ Region*Drought+ Season*Drought + Region*Season, data = ChlaA,
family = lognormal, backend = "cmdstanr", silent = 2)
summary(chlb)
pp_check(chlb)
cex = conditional_effects(chlb, ask = FALSE)
cexSave all model outputs
# Define file paths for output and figures for the manuscript
fp_output <- here("results/tables")
fp_plots <- here("results/figures")Anovas
anovas <-
bind_rows(
mutate(as_tibble(Aman, rownames = "Parameter"), Model = "Ammonium"),
mutate(as_tibble(NatAn, rownames = "Parameter"), Model = "Nitrate + Nitrite"),
mutate(as_tibble(Panova, rownames = "Parameter"), Model = "Ortho-phosphate"),
mutate(as_tibble(Ch1_Anova, rownames = "Parameter"), Model = "Chlorophyll")
) %>%
mutate(
`Pr(>F)` = if_else(
`Pr(>F)` < 0.001,
"< 0.001",
as.character(formatC(signif(`Pr(>F)`, 4), digits = 4, format = "fg", flag = "#"))
),
across(c(`Sum Sq`, `F value`), ~formatC(signif(.x, 4), digits = 4, format = "fg", flag = "#"))
) %>%
relocate(Model)
kbl(anovas) %>% kable_styling(fixed_thead = TRUE)| Model | Parameter | Sum Sq | Df | F value | Pr(>F) |
|---|---|---|---|---|---|
| Ammonium | (Intercept) | 74.97 | 1 | 502.8 | < 0.001 |
| Ammonium | Drought | 1.779 | 2 | 5.968 | 0.002700 |
| Ammonium | Region | 1.656 | 3 | 3.702 | 0.01160 |
| Ammonium | Season | 23.03 | 3 | 51.48 | < 0.001 |
| Ammonium | StocktonUpgrade | 0.8526 | 1 | 5.719 | 0.01706 |
| Ammonium | Drought:Region | 9.670 | 6 | 10.81 | < 0.001 |
| Ammonium | Drought:Season | 7.447 | 6 | 8.325 | < 0.001 |
| Ammonium | Region:Season | 46.16 | 9 | 34.40 | < 0.001 |
| Ammonium | Region:StocktonUpgrade | 22.85 | 3 | 51.09 | < 0.001 |
| Ammonium | Residuals | 98.55 | 661 | NA | NA |
| Nitrate + Nitrite | (Intercept) | 16.83 | 1 | 196.2 | < 0.001 |
| Nitrate + Nitrite | Drought | 1.099 | 2 | 6.409 | 0.001742 |
| Nitrate + Nitrite | Region | 40.84 | 3 | 158.7 | < 0.001 |
| Nitrate + Nitrite | Season | 1.315 | 3 | 5.112 | 0.001665 |
| Nitrate + Nitrite | Drought:Region | 3.374 | 6 | 6.557 | < 0.001 |
| Nitrate + Nitrite | Drought:Season | 2.380 | 6 | 4.626 | < 0.001 |
| Nitrate + Nitrite | Region:Season | 6.875 | 9 | 8.906 | < 0.001 |
| Nitrate + Nitrite | Residuals | 61.58 | 718 | NA | NA |
| Ortho-phosphate | (Intercept) | 195.7 | 1 |
|
< 0.001 |
| Ortho-phosphate | Drought | 1.298 | 2 | 6.968 | 0.001007 |
| Ortho-phosphate | Region | 11.07 | 3 | 39.60 | < 0.001 |
| Ortho-phosphate | Season | 2.329 | 3 | 8.333 | < 0.001 |
| Ortho-phosphate | Drought:Region | 1.814 | 6 | 3.246 | 0.003731 |
| Ortho-phosphate | Drought:Season | 1.550 | 6 | 2.773 | 0.01132 |
| Ortho-phosphate | Region:Season | 7.531 | 9 | 8.983 | < 0.001 |
| Ortho-phosphate | Residuals | 67.07 | 720 | NA | NA |
| Chlorophyll | (Intercept) | 2.067 | 1 | 8.185 | 0.004347 |
| Chlorophyll | Drought | 1.138 | 2 | 2.252 | 0.1059 |
| Chlorophyll | Region | 34.55 | 3 | 45.60 | < 0.001 |
| Chlorophyll | Season | 13.18 | 3 | 17.39 | < 0.001 |
| Chlorophyll | Regime | 30.42 | 1 | 120.4 | < 0.001 |
| Chlorophyll | Drought:Region | 9.933 | 6 | 6.555 | < 0.001 |
| Chlorophyll | Drought:Season | 6.877 | 6 | 4.538 | < 0.001 |
| Chlorophyll | Region:Season | 30.34 | 9 | 13.35 | < 0.001 |
| Chlorophyll | Region:Regime | 10.46 | 3 | 13.81 | < 0.001 |
| Chlorophyll | Residuals | 181.1 | 717 | NA | NA |
# Excel is the worst, I had to add letters to these numbers so it wouldn't re-round them
anovas %>%
mutate(
across(
c(`Sum Sq`, `F value`, `Pr(>F)`),
~paste0(ifelse(str_detect(.x, fixed("NA")) | is.na(.x), "", .x), "ZZZ")
)
) %>%
write_csv(file.path(fp_output, "Nutr_Chla_anovas.csv"))Post-hoc tests
tukeys <-
bind_rows(
mutate(AmRegions$contrasts, Model = "Ammonium", `Contrast type` = "region"),
mutate(AmSeasons$contrasts, Model = "Ammonium", `Contrast type` = "season"),
mutate(NNRegion$contrasts, Model = "Nitrate + Nitrite", `Contrast type` = "region"),
mutate(NNSeason$contrasts, Model = "Nitrate + Nitrite", `Contrast type` = "season"),
mutate(OPRegion$contrasts, Model = "Ortho-phosphate", `Contrast type` = "region"),
mutate(OPSeason$contrasts, Model = "Ortho-phosphate", `Contrast type` = "season"),
mutate(ChRegion$contrasts, Model = "Chlorophyll", `Contrast type` = "region"),
mutate(CHSeason$contrasts, Model = "Chlorophyll", `Contrast type` = "season")
) %>%
rename(
`p-value` = p.value,
`t-ratio` = t.ratio
) %>%
mutate(
`p-value` = if_else(
`p-value` < 0.001,
"< 0.001",
as.character(formatC(signif(`p-value`, 4), digits = 4, format = "fg", flag = "#"))
),
Group = if_else(is.na(Region), as.character(Season), as.character(Region)),
across(c(estimate, SE, `t-ratio`), ~formatC(signif(.x, 4), digits = 4, format = "fg", flag = "#"))
) %>%
select(
Model,
`Contrast type`,
Group,
Contrast = contrast,
Estimate = estimate,
SE,
DF = df,
`t-ratio`,
`p-value`
)
kbl(tukeys) %>% kable_styling(fixed_thead = TRUE)| Model | Contrast type | Group | Contrast | Estimate | SE | DF | t-ratio | p-value |
|---|---|---|---|---|---|---|---|---|
| Ammonium | region | Suisun Bay | D - N | 0.1594 | 0.06957 | 661 | 2.290 | 0.05783 |
| Ammonium | region | Suisun Bay | D - W | 0.2000 | 0.07670 | 661 | 2.607 | 0.02527 |
| Ammonium | region | Suisun Bay | N - W | 0.04061 | 0.07791 | 661 | 0.5212 | 0.8610 |
| Ammonium | region | Confluence | D - N | 0.09564 | 0.06987 | 661 | 1.369 | 0.3580 |
| Ammonium | region | Confluence | D - W | 0.1974 | 0.07671 | 661 | 2.573 | 0.02774 |
| Ammonium | region | Confluence | N - W | 0.1017 | 0.07829 | 661 | 1.300 | 0.3959 |
| Ammonium | region | South-Central Delta | D - N | 0.1261 | 0.07225 | 661 | 1.746 | 0.1891 |
| Ammonium | region | South-Central Delta | D - W | 0.006499 | 0.08031 | 661 | 0.08093 | 0.9964 |
| Ammonium | region | South-Central Delta | N - W | -0.1196 | 0.08044 | 661 | -1.487 | 0.2978 |
| Ammonium | region | North Delta | D - N | 0.3639 | 0.07225 | 661 | 5.037 | < 0.001 |
| Ammonium | region | North Delta | D - W | 0.8457 | 0.08031 | 661 | 10.53 | < 0.001 |
| Ammonium | region | North Delta | N - W | 0.4817 | 0.08044 | 661 | 5.988 | < 0.001 |
| Ammonium | season | Winter | D - N | 0.1491 | 0.07006 | 661 | 2.128 | 0.08492 |
| Ammonium | season | Winter | D - W | 0.4564 | 0.07456 | 661 | 6.121 | < 0.001 |
| Ammonium | season | Winter | N - W | 0.3073 | 0.07781 | 661 | 3.949 | < 0.001 |
| Ammonium | season | Spring | D - N | 0.4229 | 0.06954 | 661 | 6.081 | < 0.001 |
| Ammonium | season | Spring | D - W | 0.6261 | 0.07364 | 661 | 8.502 | < 0.001 |
| Ammonium | season | Spring | N - W | 0.2032 | 0.07816 | 661 | 2.600 | 0.02576 |
| Ammonium | season | Summer | D - N | 0.1227 | 0.06954 | 661 | 1.764 | 0.1824 |
| Ammonium | season | Summer | D - W | 0.04805 | 0.07462 | 661 | 0.6440 | 0.7958 |
| Ammonium | season | Summer | N - W | -0.07465 | 0.07908 | 661 | -0.9440 | 0.6126 |
| Ammonium | season | Fall | D - N | 0.05035 | 0.06954 | 661 | 0.7240 | 0.7494 |
| Ammonium | season | Fall | D - W | 0.1189 | 0.07364 | 661 | 1.615 | 0.2401 |
| Ammonium | season | Fall | N - W | 0.06857 | 0.07816 | 661 | 0.8773 | 0.6548 |
| Nitrate + Nitrite | region | Suisun Bay | D - N | 0.3404 | 0.05065 | 718 | 6.721 | < 0.001 |
| Nitrate + Nitrite | region | Suisun Bay | D - W | 0.3864 | 0.05311 | 718 | 7.275 | < 0.001 |
| Nitrate + Nitrite | region | Suisun Bay | N - W | 0.04596 | 0.05761 | 718 | 0.7978 | 0.7045 |
| Nitrate + Nitrite | region | Confluence | D - N | 0.2665 | 0.05065 | 718 | 5.262 | < 0.001 |
| Nitrate + Nitrite | region | Confluence | D - W | 0.3753 | 0.05311 | 718 | 7.066 | < 0.001 |
| Nitrate + Nitrite | region | Confluence | N - W | 0.1087 | 0.05761 | 718 | 1.888 | 0.1429 |
| Nitrate + Nitrite | region | South-Central Delta | D - N | 0.1631 | 0.05052 | 718 | 3.229 | 0.003715 |
| Nitrate + Nitrite | region | South-Central Delta | D - W | 0.1054 | 0.05336 | 718 | 1.975 | 0.1192 |
| Nitrate + Nitrite | region | South-Central Delta | N - W | -0.05775 | 0.05795 | 718 | -0.9966 | 0.5792 |
| Nitrate + Nitrite | region | North Delta | D - N | 0.02974 | 0.05052 | 718 | 0.5887 | 0.8263 |
| Nitrate + Nitrite | region | North Delta | D - W | 0.07402 | 0.05336 | 718 | 1.387 | 0.3481 |
| Nitrate + Nitrite | region | North Delta | N - W | 0.04428 | 0.05795 | 718 | 0.7641 | 0.7251 |
| Nitrate + Nitrite | season | Winter | D - N | 0.006385 | 0.05077 | 718 | 0.1258 | 0.9913 |
| Nitrate + Nitrite | season | Winter | D - W | 0.09362 | 0.05397 | 718 | 1.734 | 0.1931 |
| Nitrate + Nitrite | season | Winter | N - W | 0.08723 | 0.05829 | 718 | 1.496 | 0.2932 |
| Nitrate + Nitrite | season | Spring | D - N | 0.2769 | 0.05052 | 718 | 5.480 | < 0.001 |
| Nitrate + Nitrite | season | Spring | D - W | 0.3857 | 0.05299 | 718 | 7.279 | < 0.001 |
| Nitrate + Nitrite | season | Spring | N - W | 0.1088 | 0.05761 | 718 | 1.889 | 0.1426 |
| Nitrate + Nitrite | season | Summer | D - N | 0.2613 | 0.05052 | 718 | 5.173 | < 0.001 |
| Nitrate + Nitrite | season | Summer | D - W | 0.2268 | 0.05299 | 718 | 4.280 | < 0.001 |
| Nitrate + Nitrite | season | Summer | N - W | -0.03456 | 0.05761 | 718 | -0.5999 | 0.8202 |
| Nitrate + Nitrite | season | Fall | D - N | 0.2552 | 0.05052 | 718 | 5.051 | < 0.001 |
| Nitrate + Nitrite | season | Fall | D - W | 0.2350 | 0.05299 | 718 | 4.434 | < 0.001 |
| Nitrate + Nitrite | season | Fall | N - W | -0.02026 | 0.05761 | 718 | -0.3516 | 0.9341 |
| Ortho-phosphate | region | Suisun Bay | D - N | 0.1889 | 0.05278 | 720 | 3.579 | 0.001070 |
| Ortho-phosphate | region | Suisun Bay | D - W | 0.3938 | 0.05535 | 720 | 7.115 | < 0.001 |
| Ortho-phosphate | region | Suisun Bay | N - W | 0.2048 | 0.06003 | 720 | 3.412 | 0.001965 |
| Ortho-phosphate | region | Confluence | D - N | 0.1669 | 0.05278 | 720 | 3.161 | 0.004656 |
| Ortho-phosphate | region | Confluence | D - W | 0.4292 | 0.05535 | 720 | 7.755 | < 0.001 |
| Ortho-phosphate | region | Confluence | N - W | 0.2624 | 0.06003 | 720 | 4.370 | < 0.001 |
| Ortho-phosphate | region | South-Central Delta | D - N | 0.04469 | 0.05265 | 720 | 0.8489 | 0.6727 |
| Ortho-phosphate | region | South-Central Delta | D - W | 0.2018 | 0.05522 | 720 | 3.654 | < 0.001 |
| Ortho-phosphate | region | South-Central Delta | N - W | 0.1571 | 0.06003 | 720 | 2.617 | 0.02454 |
| Ortho-phosphate | region | North Delta | D - N | 0.07904 | 0.05265 | 720 | 1.501 | 0.2909 |
| Ortho-phosphate | region | North Delta | D - W | 0.4834 | 0.05522 | 720 | 8.753 | < 0.001 |
| Ortho-phosphate | region | North Delta | N - W | 0.4043 | 0.06003 | 720 | 6.735 | < 0.001 |
| Ortho-phosphate | season | Winter | D - N | -0.03507 | 0.05291 | 720 | -0.6627 | 0.7852 |
| Ortho-phosphate | season | Winter | D - W | 0.2475 | 0.05547 | 720 | 4.462 | < 0.001 |
| Ortho-phosphate | season | Winter | N - W | 0.2826 | 0.06003 | 720 | 4.707 | < 0.001 |
| Ortho-phosphate | season | Spring | D - N | 0.2338 | 0.05265 | 720 | 4.440 | < 0.001 |
| Ortho-phosphate | season | Spring | D - W | 0.4786 | 0.05522 | 720 | 8.667 | < 0.001 |
| Ortho-phosphate | season | Spring | N - W | 0.2448 | 0.06003 | 720 | 4.078 | < 0.001 |
| Ortho-phosphate | season | Summer | D - N | 0.1545 | 0.05265 | 720 | 2.935 | 0.009634 |
| Ortho-phosphate | season | Summer | D - W | 0.3928 | 0.05522 | 720 | 7.114 | < 0.001 |
| Ortho-phosphate | season | Summer | N - W | 0.2383 | 0.06003 | 720 | 3.970 | < 0.001 |
| Ortho-phosphate | season | Fall | D - N | 0.1263 | 0.05265 | 720 | 2.399 | 0.04404 |
| Ortho-phosphate | season | Fall | D - W | 0.3892 | 0.05522 | 720 | 7.048 | < 0.001 |
| Ortho-phosphate | season | Fall | N - W | 0.2629 | 0.06003 | 720 | 4.380 | < 0.001 |
| Chlorophyll | region | Suisun Bay | D - N | -0.1665 | 0.09023 | 717 | -1.846 | 0.1556 |
| Chlorophyll | region | Suisun Bay | D - W | -0.02247 | 0.09255 | 717 | -0.2428 | 0.9680 |
| Chlorophyll | region | Suisun Bay | N - W | 0.1441 | 0.09919 | 717 | 1.452 | 0.3147 |
| Chlorophyll | region | Confluence | D - N | 0.1505 | 0.09006 | 717 | 1.671 | 0.2172 |
| Chlorophyll | region | Confluence | D - W | 0.2379 | 0.09237 | 717 | 2.575 | 0.02753 |
| Chlorophyll | region | Confluence | N - W | 0.08743 | 0.09919 | 717 | 0.8814 | 0.6522 |
| Chlorophyll | region | South-Central Delta | D - N | 0.4456 | 0.09006 | 717 | 4.947 | < 0.001 |
| Chlorophyll | region | South-Central Delta | D - W | 0.6036 | 0.09237 | 717 | 6.534 | < 0.001 |
| Chlorophyll | region | South-Central Delta | N - W | 0.1580 | 0.09919 | 717 | 1.593 | 0.2493 |
| Chlorophyll | region | North Delta | D - N | -0.01498 | 0.09006 | 717 | -0.1663 | 0.9849 |
| Chlorophyll | region | North Delta | D - W | 0.01042 | 0.09237 | 717 | 0.1128 | 0.9930 |
| Chlorophyll | region | North Delta | N - W | 0.02540 | 0.09919 | 717 | 0.2561 | 0.9645 |
| Chlorophyll | season | Winter | D - N | 0.3760 | 0.08775 | 717 | 4.285 | < 0.001 |
| Chlorophyll | season | Winter | D - W | 0.4883 | 0.09149 | 717 | 5.337 | < 0.001 |
| Chlorophyll | season | Winter | N - W | 0.1123 | 0.09894 | 717 | 1.135 | 0.4929 |
| Chlorophyll | season | Spring | D - N | 0.1524 | 0.08755 | 717 | 1.740 | 0.1910 |
| Chlorophyll | season | Spring | D - W | 0.2602 | 0.09129 | 717 | 2.850 | 0.01248 |
| Chlorophyll | season | Spring | N - W | 0.1078 | 0.09894 | 717 | 1.090 | 0.5206 |
| Chlorophyll | season | Summer | D - N | -0.08833 | 0.08755 | 717 | -1.009 | 0.5715 |
| Chlorophyll | season | Summer | D - W | 0.1255 | 0.09129 | 717 | 1.375 | 0.3548 |
| Chlorophyll | season | Summer | N - W | 0.2138 | 0.09894 | 717 | 2.161 | 0.07870 |
| Chlorophyll | season | Fall | D - N | -0.02548 | 0.08755 | 717 | -0.2911 | 0.9544 |
| Chlorophyll | season | Fall | D - W | -0.04451 | 0.09129 | 717 | -0.4876 | 0.8772 |
| Chlorophyll | season | Fall | N - W | -0.01903 | 0.09894 | 717 | -0.1923 | 0.9798 |
# Excel is the worst, I had to add letters to these numbers so it wouldn't re-round them
tukeys %>%
mutate(
across(
c(Estimate, SE, `t-ratio`, `p-value`),
~paste0(ifelse(str_detect(.x, fixed("NA")) | is.na(.x), "", .x), "ZZZ")
)
) %>%
write_csv(file.path(fp_output, "Nutr_Chla_tukeys.csv"))Manuscript Figures
Combine all figures of model results for the nutrient and chlorophyll parameters into larger figures for the manuscript, one for the regional analyses and one for the seasonal analyses.
# Region plots
plt_nutr_chla_region <- comb_pub_figs(list(AmRegions$plt, NNRegion$plt, OPRegion$plt, ChRegion$plt))
ggsave(
file.path(fp_plots, "nutr_chla_region.jpg"),
plot = plt_nutr_chla_region,
width = 8,
height = 11,
units = "in",
dpi = 300
)
# Season plots
plt_nutr_chla_season <- comb_pub_figs(list(AmSeasons$plt, NNSeason$plt, OPSeason$plt, CHSeason$plt))
ggsave(
file.path(fp_plots, "nutr_chla_season.jpg"),
plot = plt_nutr_chla_season,
width = 8,
height = 11,
units = "in",
dpi = 300
)