Temperature, Salinity, and Secchi analyses
Sam Bashevkin
August 01, 2023
Load packages and functions
require(conflicted)
require(dplyr)
require(ggplot2)
require(patchwork)
require(knitr)
require(car)
require(multcomp)
require(emmeans)
require(stringr)
require(readr)
require(tidyr)
library(kableExtra)
library(here)
source(here("src/analysis/global_analysis_functions.R"))
conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr")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)
## 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)
## conflicted * 1.2.0 2023-02-01 [1] CRAN (R 4.2.2)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.1)
## devtools 2.4.5 2022-10-11 [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)
## 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)
## 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)
## 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)
## 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)
## 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)
## 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)
## 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)
## 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)
## 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 water quality measurements
wq_data <- read_rds(here("data/processed/wq/lt_avg_wq_meas.rds"))
# Prepare data for analyses
wq_data_c <- wq_data %>%
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),
Drought = factor(Drought, levels = c("D", "N", "W")),
YearAdj_s = (YearAdj - mean(YearAdj) / sd(YearAdj)),
across(c(Secchi, Salinity), list(l = log))
)Plots
Note that each facet has a different y-axis scale
By year
Temperature
ggplot(wq_data_c, aes(x=YearAdj, y=Temperature, fill=Drought))+
geom_point(shape=21, color="black")+
facet_wrap(Season~Region, scales="free_y")+
color_pal_drought()+
ylab("Temperature (°C)")+
xlab("Adjusted Water Year") +
theme_bw()+
theme(axis.text.x=element_text(angle=45, hjust=1))## Warning: Removed 4 rows containing missing values (`geom_point()`).
Salinity
ggplot(wq_data_c, aes(x=YearAdj, y=Salinity, fill=Drought))+
geom_point(shape=21, color="black")+
facet_wrap(Season~Region, scales="free_y")+
color_pal_drought()+
ylab("Salinity")+
xlab("Adjusted Water Year") +
theme_bw()+
theme(axis.text.x=element_text(angle=45, hjust=1))## Warning: Removed 4 rows containing missing values (`geom_point()`).
Secchi
Secchi is also more regionally than seasonally variable
plt_se_point <- ggplot(wq_data_c, aes(x=YearAdj, y=Secchi, fill=Drought))+
geom_point(shape=21, color="black")+
facet_wrap(Season~Region, scales="free_y")+
color_pal_drought()+
ylab("Secchi depth (cm)")+
xlab("Adjusted Water Year") +
theme_bw()+
theme(axis.text.x=element_text(angle=45, hjust=1))
plt_se_point## Warning: Removed 14 rows containing missing values (`geom_point()`).
Save Secchi depth figure for Supplemental Info.
ggsave(
here("results/figures/secchi_point_region_season.jpg"),
plot = plt_se_point,
width = 11,
height = 9,
units = "in",
dpi = 300
)## Warning: Removed 14 rows containing missing values (`geom_point()`).By Drought index
Temperature
ggplot(wq_data_c, aes(x=Drought, y=Temperature, fill=Drought))+
geom_boxplot()+
facet_wrap(Season~Region, scales="free_y")+
color_pal_drought()+
ylab("Temperature (°C)")+
theme_bw()+
theme(axis.text.x=element_text(angle=45, hjust=1))## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
Salinity
ggplot(wq_data_c, aes(x=Drought, y=Salinity, fill=Drought))+
geom_boxplot()+
facet_wrap(Season~Region, scales="free_y")+
color_pal_drought()+
ylab("Salinity")+
theme_bw()+
theme(axis.text.x=element_text(angle=45, hjust=1))## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
Secchi
ggplot(wq_data_c, aes(x=Drought, y=Secchi, fill=Drought))+
geom_boxplot()+
facet_wrap(Season~Region, scales="free_y")+
color_pal_drought()+
ylab("Secchi depth (cm)")+
theme_bw()+
theme(axis.text.x=element_text(angle=45, hjust=1))## Warning: Removed 14 rows containing non-finite values (`stat_boxplot()`).
Analyses
Temperature
m_t<-aov(Temperature ~ (Drought + Season + Region)^2, data=wq_data_c)Check assumptions:
model_plotter(m_t, wq_data_c, "Temperature")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Check results
m_t_Anova<-Anova(m_t, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
m_t_Anova## Anova Table (Type III tests)
##
## Response: Temperature
## Sum Sq Df F value Pr(>F)
## (Intercept) 3033.95 1 4799.0185 < 2.2e-16 ***
## Drought 10.61 2 8.3885 0.0002457 ***
## Season 2582.95 3 1361.8787 < 2.2e-16 ***
## Region 5.43 4 2.1482 0.0730497 .
## Drought:Season 60.77 6 16.0197 < 2.2e-16 ***
## Drought:Region 31.80 8 6.2876 6.062e-08 ***
## Season:Region 178.04 12 23.4679 < 2.2e-16 ***
## Residuals 568.98 900
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Partial R2
How much variability is explained by the Drought index?
partial.r2(m_t_Anova)## [1] 0.1534957Post-hoc test
m_t_contrasts_r <- pub_figure_plotter(df_data=wq_data_c, param=Temperature, y_label="Temperature (°C)",
fct_grp=Region, model=m_t, plt_title = "Water temperature")
m_t_contrasts_r$diffs## # A tibble: 5 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Marsh 16.4 16.2 0.260
## 2 Suisun Bay 16.3 15.8 0.448
## 3 Confluence 16.5 16.0 0.540
## 4 South-Central Delta 17.6 16.6 0.902
## 5 North Delta 16.4 15.0 1.42m_t_contrasts_s <- pub_figure_plotter(df_data=wq_data_c, param=Temperature, y_label="Temperature (°C)",
fct_grp=Season, model=m_t, plt_title = "Water temperature")
m_t_contrasts_s$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 10.1 10.2 -0.106
## 2 Spring 16.4 14.9 1.47
## 3 Summer 21.8 21.3 0.536
## 4 Fall 18.3 17.3 0.952Salinity
Salinity data had to be log transformed to fit the assumption of normality
m_sa<-aov(Salinity_l ~ (Drought + Season + Region)^2, data=wq_data_c)Check assumptions:
model_plotter(m_sa, wq_data_c, "Salinity_l")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Check results
m_sa_Anova<-Anova(m_sa, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
m_sa_Anova## Anova Table (Type III tests)
##
## Response: Salinity_l
## Sum Sq Df F value Pr(>F)
## (Intercept) 61.12 1 238.5221 < 2.2e-16 ***
## Drought 30.41 2 59.3430 < 2.2e-16 ***
## Season 24.91 3 32.3992 < 2.2e-16 ***
## Region 536.06 4 523.0062 < 2.2e-16 ***
## Drought:Season 8.41 6 5.4731 1.414e-05 ***
## Drought:Region 46.45 8 22.6613 < 2.2e-16 ***
## Season:Region 45.49 12 14.7934 < 2.2e-16 ***
## Residuals 230.62 900
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Partial R2
How much variability is explained by the Drought index?
partial.r2(m_sa_Anova)## [1] 0.2699633Post-hoc test
m_sa_contrasts_r <- pub_figure_plotter(df_data=wq_data_c, param=Salinity, y_label="Salinity (PSU)",
fct_grp=Region, model=m_sa, log_trans = TRUE, plt_title = "Salinity")
m_sa_contrasts_r$diffs## # A tibble: 5 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Marsh 4.68 1.22 3.46
## 2 Suisun Bay 8.51 2.10 6.41
## 3 Confluence 1.19 0.234 0.956
## 4 South-Central Delta 0.264 0.142 0.123
## 5 North Delta 0.0782 0.0636 0.0147m_sa_contrasts_s <- pub_figure_plotter(df_data=wq_data_c, param=Salinity, y_label="Salinity (PSU)",
fct_grp=Season, model=m_sa, log_trans = TRUE, plt_title = "Salinity")
m_sa_contrasts_s$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 0.999 0.375 0.623
## 2 Spring 0.706 0.199 0.507
## 3 Summer 1.06 0.362 0.694
## 4 Fall 1.32 0.569 0.753Secchi Depth
Secchi depth had to be log-transformed to comply with the assumption of normality
m_se<-aov(Secchi_l ~ (Drought + Season + Region)^2, data=wq_data_c)Check assumptions:
model_plotter(m_se, wq_data_c, "Secchi_l")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Check results
m_se_Anova<-Anova(m_se, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
m_se_Anova## Anova Table (Type III tests)
##
## Response: Secchi_l
## Sum Sq Df F value Pr(>F)
## (Intercept) 405.41 1 5467.3158 < 2.2e-16 ***
## Drought 3.90 2 26.3122 7.904e-12 ***
## Season 5.64 3 25.3312 9.897e-16 ***
## Region 18.56 4 62.5914 < 2.2e-16 ***
## Drought:Season 0.58 6 1.2997 0.25454
## Drought:Region 1.41 8 2.3745 0.01563 *
## Season:Region 5.32 12 5.9819 4.112e-10 ***
## Residuals 65.99 890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Partial R2
How much variability is explained by the Drought index?
partial.r2(m_se_Anova)## [1] 0.08192418Post-hoc test
m_se_contrasts_r <- pub_figure_plotter(df_data=wq_data_c, param=Secchi, y_label="Secchi depth (cm)",
fct_grp=Region, model=m_se, log_trans = TRUE, plt_title = "Secchi depth")
m_se_contrasts_r$diffs## # A tibble: 5 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Marsh 34.4 24.4 10.0
## 2 Suisun Bay 45.7 34.3 11.4
## 3 Confluence 58.8 47.8 11.1
## 4 South-Central Delta 84.9 57.8 27.1
## 5 North Delta 99.5 61.1 38.4m_se_contrasts_s <- pub_figure_plotter(df_data=wq_data_c, param=Secchi, y_label="Secchi depth (cm)",
fct_grp=Season, model=m_se, log_trans = TRUE, plt_title = "Secchi depth")
m_se_contrasts_s$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 59.7 38.3 21.4
## 2 Spring 49.4 36.4 13.0
## 3 Summer 56.6 42.3 14.3
## 4 Fall 78.1 56.0 22.1Secchi Depth with annual trend
Secchi depth had to be log-transformed to comply with the assumption of normality
m_se2<-aov(Secchi_l ~ (Drought + Season + Region)^2 + YearAdj*Season*Region, data=wq_data_c)Check assumptions:
model_plotter(m_se2, wq_data_c, "Secchi_l")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Check results
m_se2_Anova<-Anova(m_se2, type=3, contrasts=list(topic=contr.sum, sys=contr.sum))
m_se2_Anova## Anova Table (Type III tests)
##
## Response: Secchi_l
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.001 1 0.0286 0.8658303
## Drought 3.551 2 38.4174 < 2.2e-16 ***
## Season 0.171 3 1.2308 0.2973327
## Region 0.843 4 4.5615 0.0011916 **
## YearAdj 0.015 1 0.3338 0.5635566
## Drought:Season 1.001 6 3.6109 0.0015231 **
## Drought:Region 1.272 8 3.4401 0.0006562 ***
## Season:Region 1.207 12 2.1758 0.0111877 *
## Season:YearAdj 0.163 3 1.1742 0.3184950
## Region:YearAdj 0.895 4 4.8438 0.0007229 ***
## Season:Region:YearAdj 1.213 12 2.1870 0.0107196 *
## Residuals 40.203 870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Partial R2
How much variability is explained by the Drought index?
partial.r2(m_se2_Anova)## [1] 0.126524Post-hoc test
m_se2_contrasts_r <- pub_figure_plotter(df_data=wq_data_c, param=Secchi, y_label="Secchi depth (cm)",
fct_grp=Region, model=m_se2, log_trans = TRUE, plt_title = "Secchi depth")
m_se2_contrasts_r$diffs## # A tibble: 5 × 4
## Region D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Suisun Marsh 33.6 24.8 8.72
## 2 Suisun Bay 45.1 34.9 10.2
## 3 Confluence 56.0 50.5 5.55
## 4 South-Central Delta 78.8 62.9 15.9
## 5 North Delta 95.3 65.2 30.1m_se2_contrasts_s <- pub_figure_plotter(df_data=wq_data_c, param=Secchi, y_label="Secchi depth (cm)",
fct_grp=Season, model=m_se2, log_trans = TRUE, plt_title = "Secchi depth")
m_se2_contrasts_s$diffs## # A tibble: 4 × 4
## Season D W difference
## <fct> <dbl> <dbl> <dbl>
## 1 Winter 58.3 39.4 18.9
## 2 Spring 47.4 38.3 9.12
## 3 Summer 54.3 44.5 9.86
## 4 Fall 73.7 59.8 13.9Save 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(m_t_Anova, rownames = "Parameter"), Model="Temperature"),
mutate(as_tibble(m_sa_Anova, rownames = "Parameter"), Model="Salinity"),
mutate(as_tibble(m_se2_Anova, rownames = "Parameter"), Model="Secchi depth")
) %>%
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) |
|---|---|---|---|---|---|
| Temperature | (Intercept) |
|
1 |
|
< 0.001 |
| Temperature | Drought | 10.61 | 2 | 8.389 | < 0.001 |
| Temperature | Season |
|
3 |
|
< 0.001 |
| Temperature | Region | 5.432 | 4 | 2.148 | 0.07305 |
| Temperature | Drought:Season | 60.77 | 6 | 16.02 | < 0.001 |
| Temperature | Drought:Region | 31.80 | 8 | 6.288 | < 0.001 |
| Temperature | Season:Region | 178.0 | 12 | 23.47 | < 0.001 |
| Temperature | Residuals | 569.0 | 900 | NA | NA |
| Salinity | (Intercept) | 61.12 | 1 | 238.5 | < 0.001 |
| Salinity | Drought | 30.41 | 2 | 59.34 | < 0.001 |
| Salinity | Season | 24.91 | 3 | 32.40 | < 0.001 |
| Salinity | Region | 536.1 | 4 | 523.0 | < 0.001 |
| Salinity | Drought:Season | 8.414 | 6 | 5.473 | < 0.001 |
| Salinity | Drought:Region | 46.45 | 8 | 22.66 | < 0.001 |
| Salinity | Season:Region | 45.49 | 12 | 14.79 | < 0.001 |
| Salinity | Residuals | 230.6 | 900 | NA | NA |
| Secchi depth | (Intercept) | 0.001320 | 1 | 0.02856 | 0.8658 |
| Secchi depth | Drought | 3.551 | 2 | 38.42 | < 0.001 |
| Secchi depth | Season | 0.1706 | 3 | 1.231 | 0.2973 |
| Secchi depth | Region | 0.8432 | 4 | 4.561 | 0.001192 |
| Secchi depth | YearAdj | 0.01543 | 1 | 0.3338 | 0.5636 |
| Secchi depth | Drought:Season | 1.001 | 6 | 3.611 | 0.001523 |
| Secchi depth | Drought:Region | 1.272 | 8 | 3.440 | < 0.001 |
| Secchi depth | Season:Region | 1.207 | 12 | 2.176 | 0.01119 |
| Secchi depth | Season:YearAdj | 0.1628 | 3 | 1.174 | 0.3185 |
| Secchi depth | Region:YearAdj | 0.8953 | 4 | 4.844 | < 0.001 |
| Secchi depth | Season:Region:YearAdj | 1.213 | 12 | 2.187 | 0.01072 |
| Secchi depth | Residuals | 40.20 | 870 | 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, "TempSalSecchi_anovas.csv"))Post-hoc tests
tukeys <-
bind_rows(
mutate(as_tibble(m_t_contrasts_r$contrasts), Model="Temperature", `Contrast type`="region"),
mutate(as_tibble(m_t_contrasts_s$contrasts), Model="Temperature", `Contrast type`="season"),
mutate(as_tibble(m_sa_contrasts_r$contrasts), Model="Salinity", `Contrast type`="region"),
mutate(as_tibble(m_sa_contrasts_s$contrasts), Model="Salinity", `Contrast type`="season"),
mutate(as_tibble(m_se2_contrasts_r$contrasts), Model="Secchi depth", `Contrast type`="region"),
mutate(as_tibble(m_se2_contrasts_s$contrasts), Model="Secchi depth", `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 |
|---|---|---|---|---|---|---|---|---|
| Temperature | region | Suisun Marsh | D - N | -0.09981 | 0.1379 | 900 | -0.7240 | 0.7493 |
| Temperature | region | Suisun Marsh | D - W | 0.2605 | 0.1465 | 900 | 1.778 | 0.1776 |
| Temperature | region | Suisun Marsh | N - W | 0.3603 | 0.1583 | 900 | 2.277 | 0.05962 |
| Temperature | region | Suisun Bay | D - N | 0.2205 | 0.1372 | 900 | 1.608 | 0.2429 |
| Temperature | region | Suisun Bay | D - W | 0.4481 | 0.1439 | 900 | 3.114 | 0.005395 |
| Temperature | region | Suisun Bay | N - W | 0.2275 | 0.1564 | 900 | 1.455 | 0.3134 |
| Temperature | region | Confluence | D - N | 0.2060 | 0.1372 | 900 | 1.501 | 0.2907 |
| Temperature | region | Confluence | D - W | 0.5403 | 0.1439 | 900 | 3.756 | < 0.001 |
| Temperature | region | Confluence | N - W | 0.3343 | 0.1564 | 900 | 2.138 | 0.08292 |
| Temperature | region | South-Central Delta | D - N | 0.6159 | 0.1372 | 900 | 4.490 | < 0.001 |
| Temperature | region | South-Central Delta | D - W | 0.9019 | 0.1439 | 900 | 6.269 | < 0.001 |
| Temperature | region | South-Central Delta | N - W | 0.2860 | 0.1564 | 900 | 1.829 | 0.1608 |
| Temperature | region | North Delta | D - N | 0.7882 | 0.1372 | 900 | 5.746 | < 0.001 |
| Temperature | region | North Delta | D - W | 1.421 | 0.1439 | 900 | 9.876 | < 0.001 |
| Temperature | region | North Delta | N - W | 0.6326 | 0.1564 | 900 | 4.045 | < 0.001 |
| Temperature | season | Winter | D - N | -0.1974 | 0.1227 | 900 | -1.609 | 0.2424 |
| Temperature | season | Winter | D - W | -0.1055 | 0.1294 | 900 | -0.8155 | 0.6935 |
| Temperature | season | Winter | N - W | 0.09187 | 0.1405 | 900 | 0.6537 | 0.7903 |
| Temperature | season | Spring | D - N | 0.5911 | 0.1232 | 900 | 4.798 | < 0.001 |
| Temperature | season | Spring | D - W | 1.475 | 0.1298 | 900 | 11.36 | < 0.001 |
| Temperature | season | Spring | N - W | 0.8836 | 0.1405 | 900 | 6.287 | < 0.001 |
| Temperature | season | Summer | D - N | 0.1863 | 0.1227 | 900 | 1.518 | 0.2829 |
| Temperature | season | Summer | D - W | 0.5365 | 0.1287 | 900 | 4.169 | < 0.001 |
| Temperature | season | Summer | N - W | 0.3502 | 0.1399 | 900 | 2.504 | 0.03334 |
| Temperature | season | Fall | D - N | 0.8047 | 0.1227 | 900 | 6.559 | < 0.001 |
| Temperature | season | Fall | D - W | 0.9516 | 0.1287 | 900 | 7.395 | < 0.001 |
| Temperature | season | Fall | N - W | 0.1469 | 0.1399 | 900 | 1.050 | 0.5453 |
| Salinity | region | Suisun Marsh | D - N | 0.5693 | 0.08777 | 900 | 6.487 | < 0.001 |
| Salinity | region | Suisun Marsh | D - W | 1.342 | 0.09327 | 900 | 14.39 | < 0.001 |
| Salinity | region | Suisun Marsh | N - W | 0.7730 | 0.1007 | 900 | 7.672 | < 0.001 |
| Salinity | region | Suisun Bay | D - N | 0.7397 | 0.08733 | 900 | 8.470 | < 0.001 |
| Salinity | region | Suisun Bay | D - W | 1.401 | 0.09159 | 900 | 15.30 | < 0.001 |
| Salinity | region | Suisun Bay | N - W | 0.6617 | 0.09957 | 900 | 6.646 | < 0.001 |
| Salinity | region | Confluence | D - N | 0.9571 | 0.08733 | 900 | 10.96 | < 0.001 |
| Salinity | region | Confluence | D - W | 1.626 | 0.09159 | 900 | 17.75 | < 0.001 |
| Salinity | region | Confluence | N - W | 0.6691 | 0.09957 | 900 | 6.720 | < 0.001 |
| Salinity | region | South-Central Delta | D - N | 0.4048 | 0.08733 | 900 | 4.636 | < 0.001 |
| Salinity | region | South-Central Delta | D - W | 0.6239 | 0.09159 | 900 | 6.811 | < 0.001 |
| Salinity | region | South-Central Delta | N - W | 0.2190 | 0.09957 | 900 | 2.200 | 0.07175 |
| Salinity | region | North Delta | D - N | 0.09576 | 0.08733 | 900 | 1.097 | 0.5164 |
| Salinity | region | North Delta | D - W | 0.2076 | 0.09159 | 900 | 2.267 | 0.06104 |
| Salinity | region | North Delta | N - W | 0.1119 | 0.09957 | 900 | 1.124 | 0.4997 |
| Salinity | season | Winter | D - N | 0.3919 | 0.07811 | 900 | 5.017 | < 0.001 |
| Salinity | season | Winter | D - W | 0.9788 | 0.08238 | 900 | 11.88 | < 0.001 |
| Salinity | season | Winter | N - W | 0.5869 | 0.08948 | 900 | 6.560 | < 0.001 |
| Salinity | season | Spring | D - N | 0.9004 | 0.07843 | 900 | 11.48 | < 0.001 |
| Salinity | season | Spring | D - W | 1.268 | 0.08266 | 900 | 15.34 | < 0.001 |
| Salinity | season | Spring | N - W | 0.3680 | 0.08948 | 900 | 4.113 | < 0.001 |
| Salinity | season | Summer | D - N | 0.5279 | 0.07811 | 900 | 6.759 | < 0.001 |
| Salinity | season | Summer | D - W | 1.070 | 0.08192 | 900 | 13.07 | < 0.001 |
| Salinity | season | Summer | N - W | 0.5425 | 0.08906 | 900 | 6.091 | < 0.001 |
| Salinity | season | Fall | D - N | 0.3932 | 0.07811 | 900 | 5.034 | < 0.001 |
| Salinity | season | Fall | D - W | 0.8436 | 0.08192 | 900 | 10.30 | < 0.001 |
| Salinity | season | Fall | N - W | 0.4504 | 0.08906 | 900 | 5.057 | < 0.001 |
| Secchi depth | region | Suisun Marsh | D - N | 0.1779 | 0.03810 | 870 | 4.669 | < 0.001 |
| Secchi depth | region | Suisun Marsh | D - W | 0.3007 | 0.04096 | 870 | 7.341 | < 0.001 |
| Secchi depth | region | Suisun Marsh | N - W | 0.1228 | 0.04288 | 870 | 2.864 | 0.01191 |
| Secchi depth | region | Suisun Bay | D - N | 0.1444 | 0.03767 | 870 | 3.834 | < 0.001 |
| Secchi depth | region | Suisun Bay | D - W | 0.2571 | 0.04020 | 870 | 6.396 | < 0.001 |
| Secchi depth | region | Suisun Bay | N - W | 0.1127 | 0.04243 | 870 | 2.656 | 0.02190 |
| Secchi depth | region | Confluence | D - N | 0.06421 | 0.03767 | 870 | 1.704 | 0.2041 |
| Secchi depth | region | Confluence | D - W | 0.1042 | 0.04020 | 870 | 2.592 | 0.02621 |
| Secchi depth | region | Confluence | N - W | 0.04000 | 0.04243 | 870 | 0.9428 | 0.6133 |
| Secchi depth | region | South-Central Delta | D - N | 0.1476 | 0.03767 | 870 | 3.918 | < 0.001 |
| Secchi depth | region | South-Central Delta | D - W | 0.2250 | 0.04020 | 870 | 5.597 | < 0.001 |
| Secchi depth | region | South-Central Delta | N - W | 0.07741 | 0.04243 | 870 | 1.824 | 0.1622 |
| Secchi depth | region | North Delta | D - N | 0.1507 | 0.04013 | 870 | 3.757 | < 0.001 |
| Secchi depth | region | North Delta | D - W | 0.3798 | 0.04045 | 870 | 9.389 | < 0.001 |
| Secchi depth | region | North Delta | N - W | 0.2290 | 0.04363 | 870 | 5.250 | < 0.001 |
| Secchi depth | season | Winter | D - N | 0.2162 | 0.03370 | 870 | 6.415 | < 0.001 |
| Secchi depth | season | Winter | D - W | 0.3911 | 0.03607 | 870 | 10.84 | < 0.001 |
| Secchi depth | season | Winter | N - W | 0.1749 | 0.03811 | 870 | 4.590 | < 0.001 |
| Secchi depth | season | Spring | D - N | 0.08475 | 0.03486 | 870 | 2.431 | 0.04038 |
| Secchi depth | season | Spring | D - W | 0.2138 | 0.03647 | 870 | 5.863 | < 0.001 |
| Secchi depth | season | Spring | N - W | 0.1291 | 0.03852 | 870 | 3.350 | 0.002423 |
| Secchi depth | season | Summer | D - N | 0.1237 | 0.03455 | 870 | 3.580 | 0.001055 |
| Secchi depth | season | Summer | D - W | 0.2003 | 0.03604 | 870 | 5.556 | < 0.001 |
| Secchi depth | season | Summer | N - W | 0.07657 | 0.03836 | 870 | 1.996 | 0.1138 |
| Secchi depth | season | Fall | D - N | 0.1233 | 0.03369 | 870 | 3.658 | < 0.001 |
| Secchi depth | season | Fall | D - W | 0.2083 | 0.03596 | 870 | 5.793 | < 0.001 |
| Secchi depth | season | Fall | N - W | 0.08502 | 0.03795 | 870 | 2.240 | 0.06515 |
# 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, "TempSalSecchi_tukeys.csv"))Manuscript Figures
Combine all figures of model results for the water quality parameters into larger figures for the manuscript, one for the regional analyses and one for the seasonal analyses.
# Region plots
plt_wq_region <- comb_pub_figs(list(m_t_contrasts_r$plt, m_sa_contrasts_r$plt, m_se2_contrasts_r$plt))
ggsave(
file.path(fp_plots, "water_quality_region.jpg"),
plot = plt_wq_region,
width = 8,
height = 9,
units = "in",
dpi = 300
)
# Season plots
plt_wq_season <- comb_pub_figs(list(m_t_contrasts_s$plt, m_sa_contrasts_s$plt, m_se2_contrasts_s$plt))
ggsave(
file.path(fp_plots, "water_quality_season.jpg"),
plot = plt_wq_season,
width = 8,
height = 9,
units = "in",
dpi = 300
)