NDFS Synthesis Manuscript: Chlorophyll analysis
Combined models using weekly average flow as continuous predictor
Dave Bosworth
October 02, 2026
Purpose
Explore and analyze the continuous chlorophyll data to be included in the NDFS synthesis manuscript. We will attempt to fit multiple models to predict weekly average chlorophyll fluorescence values. All models will only include representative stations for 4 habitat types - upstream (RD22), lower Yolo Bypass (STTD), Cache Slough complex (LIB), and downstream (RVB). At a minimum, the models will contain the two categorical variables - Year and Station - as predictor variables. In some of the models, we will add weekly average flow as a continuous predictor which replaces the categorical predictor - flow action period - in the original analysis. Additionally, we’ll add a GAM smooth for Week number term to account for seasonality in some of the models. After fitting multiple models, we’ll use a model selection process to determine the best one.
Global code and functions
# Load packages
library(tidyverse)
library(scales)
library(knitr)
library(mgcv)
library(car)
library(gratia)
library(ggeffects)
library(emmeans)
library(multcomp)
library(here)
library(conflicted)
# Source functions
source(here("manuscript_synthesis/src/global_functions.R"))
# Declare package conflict preferences
conflicts_prefer(dplyr::filter(), dplyr::lag(), dplyr::select())Display current versions of R and packages used for this analysis:
devtools::session_info()## ─ Session info ────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.6.1 (2026-06-24 ucrt)
## os Windows 11 x64 (build 26200)
## system x86_64, mingw32
## ui Positron
## language (EN)
## collate English_United States.utf8
## ctype English_United States.utf8
## tz America/Los_Angeles
## date 2026-10-02
## pandoc 3.10 @ c:\\Users\\dboswort\\AppData\\Local\\Programs\\Positron\\resources\\app\\quarto\\bin\\tools/ (via rmarkdown)
## quarto 1.9.38 @ C:\\Users\\dboswort\\AppData\\Local\\Programs\\Quarto\\bin\\quarto.exe
##
## ─ Packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-8 2024-09-12 [1] CRAN (R 4.6.1)
## bit 4.6.0 2025-03-06 [1] CRAN (R 4.6.1)
## bit64 4.8.6 2026-09-01 [1] CRAN (R 4.6.1)
## bslib 0.12.0 2026-08-04 [1] CRAN (R 4.6.1)
## cachem 1.1.0 2024-05-16 [1] CRAN (R 4.6.1)
## car * 3.1-5 2026-02-03 [1] CRAN (R 4.6.1)
## carData * 3.0-6 2026-01-30 [1] CRAN (R 4.6.1)
## cli 3.6.6 2026-04-09 [1] CRAN (R 4.6.1)
## codetools 0.2-20 2024-03-31 [1] CRAN (R 4.6.1)
## conflicted * 1.2.0 2023-02-01 [1] CRAN (R 4.6.1)
## crayon 1.5.3 2024-06-20 [1] CRAN (R 4.6.1)
## datawizard 1.4.0 2026-09-10 [1] CRAN (R 4.6.1)
## devtools 2.5.2 2026-04-30 [1] CRAN (R 4.6.1)
## digest 0.6.39 2025-11-19 [1] CRAN (R 4.6.1)
## dplyr * 1.2.1 2026-04-03 [1] CRAN (R 4.6.1)
## ellipsis 0.3.3 2026-04-04 [1] CRAN (R 4.6.1)
## emmeans * 2.0.4 2026-07-15 [1] CRAN (R 4.6.1)
## estimability 2.0.0 2026-06-26 [1] CRAN (R 4.6.1)
## evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.6.1)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.6.1)
## fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.6.1)
## forcats * 1.0.1 2025-09-25 [1] CRAN (R 4.6.1)
## Formula 1.2-6 2026-08-03 [1] CRAN (R 4.6.1)
## fs * 2.1.0 2026-04-18 [1] CRAN (R 4.6.1)
## generics 0.1.4 2025-05-09 [1] CRAN (R 4.6.1)
## ggeffects * 2.3.2 2025-12-16 [1] CRAN (R 4.6.1)
## ggokabeito 0.1.0 2021-10-18 [1] CRAN (R 4.6.1)
## ggplot2 * 4.0.3 2026-04-22 [1] CRAN (R 4.6.1)
## glue 1.8.1 2026-04-17 [1] CRAN (R 4.6.1)
## gratia * 0.11.2 2026-02-07 [1] CRAN (R 4.6.1)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.6.1)
## haven 2.5.5 2025-05-30 [1] CRAN (R 4.6.1)
## here * 1.0.2 2025-09-15 [1] CRAN (R 4.6.1)
## hms 1.1.4 2025-10-17 [1] CRAN (R 4.6.1)
## htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.6.1)
## insight 1.5.4 2026-09-05 [1] CRAN (R 4.6.1)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.6.1)
## jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.6.1)
## knitr * 1.52 2026-09-06 [1] CRAN (R 4.6.1)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.6.1)
## lattice 0.23-1 2026-08-12 [1] CRAN (R 4.6.1)
## lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.6.1)
## lubridate * 1.9.5 2026-02-04 [1] CRAN (R 4.6.1)
## magrittr 2.0.5 2026-04-04 [1] CRAN (R 4.6.1)
## MASS * 7.3-66 2026-07-15 [1] CRAN (R 4.6.1)
## Matrix 1.7-6 2026-07-25 [1] CRAN (R 4.6.1)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.6.1)
## mgcv * 1.9-4 2025-11-07 [1] CRAN (R 4.6.1)
## mirai 2.7.3 2026-09-24 [1] CRAN (R 4.6.1)
## multcomp * 1.4-32 2026-08-21 [1] CRAN (R 4.6.1)
## multcompView 0.1-12 2026-07-26 [1] CRAN (R 4.6.1)
## mvnfast 0.2.8 2023-02-23 [1] CRAN (R 4.6.1)
## mvtnorm * 1.4-2 2026-07-12 [1] CRAN (R 4.6.1)
## nanonext 1.10.3 2026-09-23 [1] CRAN (R 4.6.1)
## nlme * 3.1-171 2026-09-01 [1] CRAN (R 4.6.1)
## otel 0.2.0 2025-08-29 [1] CRAN (R 4.6.1)
## patchwork * 1.3.2 2025-08-25 [1] CRAN (R 4.6.1)
## pillar 1.11.1 2025-09-17 [1] CRAN (R 4.6.1)
## pkgbuild 1.4.8 2025-05-26 [1] CRAN (R 4.6.1)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.6.1)
## pkgload 1.5.3 2026-06-15 [1] CRAN (R 4.6.1)
## purrr * 1.2.2 2026-04-10 [1] CRAN (R 4.6.1)
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.6.1)
## ragg 1.5.2 2026-03-23 [1] CRAN (R 4.6.1)
## RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.6.1)
## Rcpp 1.1.2 2026-07-05 [1] CRAN (R 4.6.1)
## readr * 2.2.0 2026-02-19 [1] CRAN (R 4.6.1)
## renv 1.2.4 2026-08-03 [1] CRAN (R 4.6.1)
## rlang * 1.3.0 2026-07-05 [1] CRAN (R 4.6.1)
## rmarkdown 2.32 2026-09-01 [1] CRAN (R 4.6.1)
## rprojroot 2.1.1 2025-08-26 [1] CRAN (R 4.6.1)
## S7 0.2.2 2026-04-22 [1] CRAN (R 4.6.1)
## sandwich 3.1-3 2026-08-03 [1] CRAN (R 4.6.1)
## sass 0.4.10 2025-04-11 [1] CRAN (R 4.6.1)
## scales * 1.4.0 2025-04-24 [1] CRAN (R 4.6.1)
## sessioninfo 1.2.4 2026-06-04 [1] CRAN (R 4.6.1)
## statmod 1.5.2 2026-05-17 [1] CRAN (R 4.6.1)
## stringi 1.8.9 2026-08-04 [1] CRAN (R 4.6.1)
## stringr * 1.6.0 2025-11-04 [1] CRAN (R 4.6.1)
## survival * 3.8-12 2026-09-09 [1] CRAN (R 4.6.1)
## systemfonts 1.3.2 2026-03-05 [1] CRAN (R 4.6.1)
## textshaping 1.0.5 2026-03-06 [1] CRAN (R 4.6.1)
## TH.data * 1.1-5 2025-11-17 [1] CRAN (R 4.6.1)
## tibble * 3.3.1 2026-01-11 [1] CRAN (R 4.6.1)
## tidyr * 1.3.2 2025-12-19 [1] CRAN (R 4.6.1)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.6.1)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.6.1)
## timechange 0.4.0 2026-01-29 [1] CRAN (R 4.6.1)
## tweedie 3.1.0 2026-05-19 [1] CRAN (R 4.6.1)
## tzdb 0.5.0 2025-03-15 [1] CRAN (R 4.6.1)
## usethis 3.2.2 2026-09-10 [1] CRAN (R 4.6.1)
## utf8 1.2.6 2025-06-08 [1] CRAN (R 4.6.1)
## vctrs 0.7.3 2026-04-11 [1] CRAN (R 4.6.1)
## viridisLite 0.4.3 2026-02-04 [1] CRAN (R 4.6.1)
## vroom 1.7.1 2026-03-31 [1] CRAN (R 4.6.1)
## withr 3.0.3 2026-06-19 [1] CRAN (R 4.6.1)
## xfun 0.61 2026-09-16 [1] CRAN (R 4.6.1)
## xtable 1.8-8 2026-02-22 [1] CRAN (R 4.6.1)
## yaml 2.3.12 2025-12-10 [1] CRAN (R 4.6.1)
## zoo 1.9-0 2026-07-31 [1] CRAN (R 4.6.1)
##
## [1] C:/R/win-library/4.6
## [2] C:/Program Files/R/R-4.6.1/library
## * ── Packages attached to the search path.
##
## ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Import Data
# Define file path for processed data
fp_data <- here("manuscript_synthesis/data/processed")
# Import weekly average water quality data
df_wq <- readRDS(file.path(fp_data, "wq_week_avg_2013-2019.rds"))
# Import weekly average flow data
df_flow <- readRDS(file.path(fp_data, "flow_week_avg_2013-2019.rds"))Prepare Data
# Create a vector for the factor order of StationCode
sta_order <- c("RD22", "STTD", "LIB", "RVB")
# We will use LIS flow data for all stations except for RD22
df_flow_c <- bind_rows(
expand_grid(
StationCode = c("STTD", "LIB", "RVB"),
df_flow |> filter(StationCode == "LIS") |> select(-StationCode)
),
df_flow |> filter(StationCode == "RD22")
)
# Prepare chlorophyll and flow data for exploration and analysis
df_chla_c1 <- df_wq %>%
select(StationCode, Year, Week, Chla) %>%
drop_na(Chla) %>%
# Filter to only include representative stations for 4 habitat types - RD22, STTD, LIB, RVB
filter(StationCode %in% sta_order) %>%
# Join flow data to chlorophyll data
left_join(df_flow_c, by = join_by(StationCode, Year, Week)) %>%
# Remove all NA flow values
drop_na(Flow) %>%
mutate(
# Scale and log transform chlorophyll values
Chla_log = log(Chla * 1000),
# Apply factor order to StationCode
StationCode = factor(StationCode, levels = sta_order),
# Add a column for Year as a factor for the model
Year_fct = factor(Year)
) %>%
arrange(StationCode, Year, Week)Explore sample counts by Station
df_chla_c1 %>%
summarize(
min_week = min(Week),
max_week = max(Week),
num_samples = n(),
.by = c(StationCode, Year)
) %>%
arrange(StationCode, Year) %>%
kable()| StationCode | Year | min_week | max_week | num_samples |
|---|---|---|---|---|
| RD22 | 2014 | 39 | 45 | 7 |
| RD22 | 2015 | 30 | 45 | 16 |
| RD22 | 2016 | 25 | 38 | 14 |
| RD22 | 2017 | 28 | 44 | 17 |
| RD22 | 2018 | 28 | 45 | 18 |
| RD22 | 2019 | 28 | 45 | 18 |
| STTD | 2013 | 33 | 44 | 12 |
| STTD | 2014 | 30 | 45 | 14 |
| STTD | 2015 | 30 | 46 | 17 |
| STTD | 2016 | 25 | 38 | 14 |
| STTD | 2017 | 28 | 39 | 10 |
| STTD | 2018 | 29 | 42 | 14 |
| STTD | 2019 | 30 | 45 | 16 |
| LIB | 2013 | 29 | 46 | 18 |
| LIB | 2014 | 30 | 45 | 16 |
| LIB | 2015 | 27 | 46 | 20 |
| LIB | 2016 | 22 | 38 | 17 |
| LIB | 2017 | 28 | 44 | 17 |
| LIB | 2018 | 33 | 45 | 13 |
| LIB | 2019 | 28 | 45 | 14 |
| RVB | 2013 | 27 | 46 | 20 |
| RVB | 2014 | 30 | 45 | 16 |
| RVB | 2015 | 27 | 46 | 20 |
| RVB | 2016 | 22 | 38 | 16 |
| RVB | 2017 | 28 | 44 | 17 |
| RVB | 2018 | 28 | 45 | 18 |
| RVB | 2019 | 28 | 45 | 18 |
Looking at the sample counts and date ranges, we’ll only include years 2015-2019 for the analysis.
df_chla_c2 <- df_chla_c1 %>%
filter(Year %in% 2015:2019) %>%
mutate(Year_fct = fct_drop(Year_fct))We’ll create another dataframe that has up to 2 lag variables for chlorophyll to be used in the models to help with serial autocorrelation.
df_chla_c2_lag <- df_chla_c2 %>%
# Fill in missing weeks for each StationCode-Year combination
group_by(StationCode, Year) %>%
# Create lag variables of scaled log transformed chlorophyll values
mutate(
lag1 = lag(Chla_log),
lag2 = lag(Chla_log, 2)
) %>%
ungroup()Plots
Let’s explore the data with some plots. First, lets plot the data in scatter plots of chlorophyll and flow faceted by Station and grouping all years together.
df_chla_c2 %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(vars(StationCode), scales = "free") +
theme_bw()
At first glance, I’m not sure how well flow is going to be able to predict chlorophyll concentrations. At the furthest upstream station - RD22 - chlorophyll appears to be highest at the lowest flows, but the variation is at its maximum at the lowest flows. There may be some dilution effect going on here at the higher flows. At STTD, there does seem to be a modest increase in chlorophyll concentrations at the mid-range flows. There appears to be no effect of flow on chlorophyll at RVB, but the range of chlorophyll concentrations is narrow at this station (between 0 and 5).
Let’s break these scatterplots apart by year to see how these patterns vary annually.
df_chla_c2 %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(
vars(StationCode, Year),
ncol = 5,
scales = "free",
labeller = labeller(.multi_line = FALSE)
) +
theme_bw()
The patterns appear to vary annually at each station, which may justify using a 3-way interaction.
Model 1: GAM Model with Flow and 3-way interactions
First, we will attempt to fit a generalized additive model (GAM) to the data set to help account for seasonality in the data. We’ll try running a GAM using a three-way interaction between Year, Weekly Average Flow, and Station, and a cyclic penalized cubic regression spline smooth term for week number to account for seasonality (restricting the k-value to 5 to reduce overfitting). Initially, we’ll run the GAM without accounting for serial autocorrelation.
Initial Model
m_gam_flow3 <- gam(
Chla_log ~ Year_fct * Flow * StationCode + s(Week, bs = "cc", k = 5),
data = df_chla_c2,
method = "REML",
knots = list(Week = c(0, 52))
)Lets look at the model summary and diagnostics:
summary(m_gam_flow3)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * Flow * StationCode + s(Week, bs = "cc",
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.806e+00 1.020e-01 96.165 < 2e-16 ***
## Year_fct2016 -9.021e-01 1.483e-01 -6.084 3.81e-09 ***
## Year_fct2017 -9.328e-01 2.269e-01 -4.112 5.15e-05 ***
## Year_fct2018 -4.196e-01 1.361e-01 -3.082 0.002260 **
## Year_fct2019 -1.494e-01 1.337e-01 -1.117 0.264848
## Flow -2.024e-03 4.522e-04 -4.476 1.11e-05 ***
## StationCodeSTTD -1.150e+00 1.282e-01 -8.965 < 2e-16 ***
## StationCodeLIB -2.113e+00 1.240e-01 -17.042 < 2e-16 ***
## StationCodeRVB -2.159e+00 1.240e-01 -17.411 < 2e-16 ***
## Year_fct2016:Flow 2.363e-03 7.198e-04 3.284 0.001155 **
## Year_fct2017:Flow 2.789e-02 1.167e-02 2.390 0.017499 *
## Year_fct2018:Flow 5.834e-04 6.248e-04 0.934 0.351234
## Year_fct2019:Flow -1.589e-05 5.498e-04 -0.029 0.976969
## Year_fct2016:StationCodeSTTD 1.224e+00 1.885e-01 6.492 3.83e-10 ***
## Year_fct2017:StationCodeSTTD 1.137e+00 3.090e-01 3.680 0.000279 ***
## Year_fct2018:StationCodeSTTD -2.514e-01 1.801e-01 -1.396 0.163694
## Year_fct2019:StationCodeSTTD -9.074e-01 1.770e-01 -5.126 5.52e-07 ***
## Year_fct2016:StationCodeLIB 1.700e+00 1.818e-01 9.350 < 2e-16 ***
## Year_fct2017:StationCodeLIB 1.023e+00 2.585e-01 3.957 9.61e-05 ***
## Year_fct2018:StationCodeLIB -1.667e+00 1.812e-01 -9.198 < 2e-16 ***
## Year_fct2019:StationCodeLIB -5.322e-01 1.764e-01 -3.017 0.002788 **
## Year_fct2016:StationCodeRVB 9.008e-01 1.834e-01 4.913 1.52e-06 ***
## Year_fct2017:StationCodeRVB 7.195e-01 2.585e-01 2.784 0.005738 **
## Year_fct2018:StationCodeRVB 1.184e-01 1.713e-01 0.691 0.489984
## Year_fct2019:StationCodeRVB -6.463e-01 1.710e-01 -3.779 0.000193 ***
## Flow:StationCodeSTTD 5.137e-03 6.881e-04 7.465 1.05e-12 ***
## Flow:StationCodeLIB 3.292e-03 6.715e-04 4.902 1.60e-06 ***
## Flow:StationCodeRVB 2.230e-03 6.715e-04 3.321 0.001017 **
## Year_fct2016:Flow:StationCodeSTTD -4.360e-03 1.051e-03 -4.148 4.45e-05 ***
## Year_fct2017:Flow:StationCodeSTTD -1.928e-02 1.304e-02 -1.478 0.140575
## Year_fct2018:Flow:StationCodeSTTD -8.881e-04 9.253e-04 -0.960 0.337983
## Year_fct2019:Flow:StationCodeSTTD -1.556e-03 8.188e-04 -1.900 0.058447 .
## Year_fct2016:Flow:StationCodeLIB -3.309e-03 1.035e-03 -3.198 0.001543 **
## Year_fct2017:Flow:StationCodeLIB -2.818e-02 1.184e-02 -2.381 0.017942 *
## Year_fct2018:Flow:StationCodeLIB 1.085e-03 9.277e-04 1.170 0.243134
## Year_fct2019:Flow:StationCodeLIB -1.224e-03 8.155e-04 -1.500 0.134624
## Year_fct2016:Flow:StationCodeRVB -3.139e-03 1.037e-03 -3.027 0.002697 **
## Year_fct2017:Flow:StationCodeRVB -2.708e-02 1.184e-02 -2.288 0.022899 *
## Year_fct2018:Flow:StationCodeRVB -1.100e-03 9.063e-04 -1.214 0.225781
## Year_fct2019:Flow:StationCodeRVB 3.505e-04 8.031e-04 0.436 0.662882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.73 3 27.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.91 Deviance explained = 92.1%
## -REML = 231.59 Scale est. = 0.10018 n = 324appraise(m_gam_flow3)
shapiro.test(residuals(m_gam_flow3))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_flow3)
## W = 0.98047, p-value = 0.0002157k.check(m_gam_flow3)## k' edf k-index p-value
## s(Week) 3 2.729607 0.9315452 0.1025draw(m_gam_flow3, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_flow3, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_flow3))
Box.test(residuals(m_gam_flow3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_flow3)
## X-squared = 77.528, df = 20, p-value = 1.027e-08Besides the Shapiro-Wilk normality test showing that the residuals aren’t normal, the diagnostic plots look pretty good. However, the residuals are autocorrelated.
Model with lag terms
Now, we’ll try to deal with the residual autocorrelation. We’ll run a series of models adding 1 and 2 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_gam_flow3_lag1 <- gam(
Chla_log ~ Year_fct * Flow * StationCode + s(Week, bs = "cc", k = 5) + lag1,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_flow3_lag1))
Box.test(residuals(m_gam_flow3_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_flow3_lag1)
## X-squared = 31.415, df = 20, p-value = 0.04994Lag 2
m_gam_flow3_lag2 <- gam(
Chla_log ~ Year_fct *
Flow *
StationCode +
s(Week, bs = "cc", k = 5) +
lag1 +
lag2,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_flow3_lag2))
Box.test(residuals(m_gam_flow3_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_flow3_lag2)
## X-squared = 18.253, df = 20, p-value = 0.5708Let’s use AIC to see how they compare.
Compare Models
AIC(m_gam_flow3, m_gam_flow3_lag1, m_gam_flow3_lag2)## df AIC
## m_gam_flow3 43.94127 216.0624
## m_gam_flow3_lag1 44.81882 142.2910
## m_gam_flow3_lag2 45.78607 125.2611It looks like the lag2 model has the best fit according to the AIC values. Let’s take a closer look at that one.
Lag 2 model summary
summary(m_gam_flow3_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * Flow * StationCode + s(Week, bs = "cc",
## k = 5) + lag1 + lag2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.792e+00 5.382e-01 12.619 < 2e-16 ***
## Year_fct2016 -6.729e-01 1.508e-01 -4.462 1.25e-05 ***
## Year_fct2017 -7.214e-01 2.257e-01 -3.197 0.001577 **
## Year_fct2018 -3.369e-01 1.333e-01 -2.528 0.012116 *
## Year_fct2019 -8.753e-02 1.299e-01 -0.674 0.501234
## Flow -1.475e-03 4.264e-04 -3.459 0.000642 ***
## StationCodeSTTD -7.815e-01 1.356e-01 -5.763 2.54e-08 ***
## StationCodeLIB -1.502e+00 1.593e-01 -9.430 < 2e-16 ***
## StationCodeRVB -1.526e+00 1.626e-01 -9.382 < 2e-16 ***
## lag1 4.259e-01 6.179e-02 6.892 4.84e-11 ***
## lag2 -1.185e-01 5.577e-02 -2.125 0.034583 *
## Year_fct2016:Flow 1.492e-03 6.661e-04 2.240 0.026033 *
## Year_fct2017:Flow 2.404e-02 1.267e-02 1.898 0.058963 .
## Year_fct2018:Flow 6.154e-04 5.707e-04 1.078 0.281984
## Year_fct2019:Flow 1.071e-04 5.050e-04 0.212 0.832153
## Year_fct2016:StationCodeSTTD 8.949e-01 1.914e-01 4.675 4.90e-06 ***
## Year_fct2017:StationCodeSTTD 8.825e-01 3.069e-01 2.875 0.004397 **
## Year_fct2018:StationCodeSTTD -9.336e-02 1.764e-01 -0.529 0.597170
## Year_fct2019:StationCodeSTTD -6.547e-01 1.792e-01 -3.653 0.000318 ***
## Year_fct2016:StationCodeLIB 1.211e+00 1.965e-01 6.160 3.05e-09 ***
## Year_fct2017:StationCodeLIB 7.927e-01 2.547e-01 3.113 0.002080 **
## Year_fct2018:StationCodeLIB -1.260e+00 2.041e-01 -6.172 2.85e-09 ***
## Year_fct2019:StationCodeLIB -3.854e-01 1.749e-01 -2.204 0.028478 *
## Year_fct2016:StationCodeRVB 6.156e-01 1.822e-01 3.378 0.000851 ***
## Year_fct2017:StationCodeRVB 5.531e-01 2.518e-01 2.196 0.029043 *
## Year_fct2018:StationCodeRVB 1.609e-01 1.649e-01 0.975 0.330323
## Year_fct2019:StationCodeRVB -4.362e-01 1.697e-01 -2.570 0.010763 *
## Flow:StationCodeSTTD 3.549e-03 6.689e-04 5.305 2.56e-07 ***
## Flow:StationCodeLIB 2.629e-03 6.255e-04 4.203 3.71e-05 ***
## Flow:StationCodeRVB 1.718e-03 6.201e-04 2.770 0.006044 **
## Year_fct2016:Flow:StationCodeSTTD -3.022e-03 9.752e-04 -3.098 0.002177 **
## Year_fct2017:Flow:StationCodeSTTD -1.777e-02 1.382e-02 -1.286 0.199744
## Year_fct2018:Flow:StationCodeSTTD -8.478e-04 8.539e-04 -0.993 0.321815
## Year_fct2019:Flow:StationCodeSTTD -1.063e-03 7.541e-04 -1.410 0.159740
## Year_fct2016:Flow:StationCodeLIB -2.382e-03 9.462e-04 -2.517 0.012488 *
## Year_fct2017:Flow:StationCodeLIB -2.445e-02 1.277e-02 -1.916 0.056606 .
## Year_fct2018:Flow:StationCodeLIB 9.517e-04 8.607e-04 1.106 0.269909
## Year_fct2019:Flow:StationCodeLIB -1.202e-03 7.458e-04 -1.611 0.108444
## Year_fct2016:Flow:StationCodeRVB -2.025e-03 9.502e-04 -2.131 0.034120 *
## Year_fct2017:Flow:StationCodeRVB -2.342e-02 1.277e-02 -1.834 0.067832 .
## Year_fct2018:Flow:StationCodeRVB -1.097e-03 8.253e-04 -1.329 0.185223
## Year_fct2019:Flow:StationCodeRVB -2.692e-05 7.335e-04 -0.037 0.970749
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.519 3 15.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.93 Deviance explained = 94.1%
## -REML = 190.54 Scale est. = 0.078179 n = 284appraise(m_gam_flow3_lag2)
shapiro.test(residuals(m_gam_flow3_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_flow3_lag2)
## W = 0.96943, p-value = 9.622e-06k.check(m_gam_flow3_lag2)## k' edf k-index p-value
## s(Week) 3 2.519218 0.9474706 0.1325draw(m_gam_flow3_lag2, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_flow3_lag2, pages = 1, all.terms = TRUE)
anova(m_gam_flow3_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * Flow * StationCode + s(Week, bs = "cc",
## k = 5) + lag1 + lag2
##
## Parametric Terms:
## df F p-value
## Year_fct 4 7.368 1.30e-05
## Flow 1 11.963 0.000642
## StationCode 3 33.991 < 2e-16
## lag1 1 47.497 4.84e-11
## lag2 1 4.517 0.034583
## Year_fct:Flow 4 2.449 0.046957
## Year_fct:StationCode 12 9.572 3.84e-15
## Flow:StationCode 3 10.730 1.21e-06
## Year_fct:Flow:StationCode 12 2.731 0.001745
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.519 3.000 15.38 <2e-16The model diagnostics look pretty good. We’ll use
m_gam_flow3_lag2 in the model selection process.
rm(m_gam_flow3, m_gam_flow3_lag1)Model 2: GAM Model with Flow and 2-way interactions
Now we’ll try running a GAM using all two-way interactions between Year, Flow, and Station.
Initial Model
m_gam_flow2 <- gam(
Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(Week, bs = "cc", k = 5),
data = df_chla_c2,
method = "REML",
knots = list(Week = c(0, 52))
)Lets look at the model summary and diagnostics:
summary(m_gam_flow2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(Week, bs = "cc",
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.7205112 0.0967498 100.471 < 2e-16 ***
## Year_fct2016 -0.5993849 0.1358856 -4.411 1.45e-05 ***
## Year_fct2017 -0.4210647 0.1288681 -3.267 0.001214 **
## Year_fct2018 -0.3812812 0.1255067 -3.038 0.002596 **
## Year_fct2019 -0.0783928 0.1246313 -0.629 0.529840
## Flow -0.0013679 0.0003299 -4.146 4.43e-05 ***
## StationCodeSTTD -1.0301452 0.1249632 -8.244 5.65e-15 ***
## StationCodeLIB -2.0236206 0.1216646 -16.633 < 2e-16 ***
## StationCodeRVB -2.0728187 0.1216470 -17.040 < 2e-16 ***
## Year_fct2016:Flow -0.0003023 0.0004073 -0.742 0.458519
## Year_fct2017:Flow 0.0011291 0.0016908 0.668 0.504805
## Year_fct2018:Flow 0.0003558 0.0003546 1.003 0.316528
## Year_fct2019:Flow -0.0005446 0.0003159 -1.724 0.085763 .
## Year_fct2016:StationCodeSTTD 0.8354066 0.1778523 4.697 4.06e-06 ***
## Year_fct2017:StationCodeSTTD 0.3708812 0.1977241 1.876 0.061683 .
## Year_fct2018:StationCodeSTTD -0.3116290 0.1722282 -1.809 0.071413 .
## Year_fct2019:StationCodeSTTD -1.0239195 0.1690948 -6.055 4.28e-09 ***
## Year_fct2016:StationCodeLIB 1.3567886 0.1726755 7.857 7.49e-14 ***
## Year_fct2017:StationCodeLIB 0.5463886 0.1826319 2.992 0.003009 **
## Year_fct2018:StationCodeLIB -1.5897173 0.1723422 -9.224 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.6547315 0.1704938 -3.840 0.000151 ***
## Year_fct2016:StationCodeRVB 0.5541429 0.1739615 3.185 0.001601 **
## Year_fct2017:StationCodeRVB 0.2230409 0.1825689 1.222 0.222810
## Year_fct2018:StationCodeRVB 0.0524243 0.1650727 0.318 0.751027
## Year_fct2019:StationCodeRVB -0.6484061 0.1646652 -3.938 0.000103 ***
## Flow:StationCodeSTTD 0.0036332 0.0003164 11.483 < 2e-16 ***
## Flow:StationCodeLIB 0.0025281 0.0003230 7.827 9.15e-14 ***
## Flow:StationCodeRVB 0.0017230 0.0003131 5.503 8.14e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.709 3 29.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.897 Deviance explained = 90.6%
## -REML = 188.07 Scale est. = 0.11471 n = 324appraise(m_gam_flow2)
shapiro.test(residuals(m_gam_flow2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_flow2)
## W = 0.98509, p-value = 0.001963k.check(m_gam_flow2)## k' edf k-index p-value
## s(Week) 3 2.709144 0.9617908 0.2525draw(m_gam_flow2, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_flow2, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_flow2))
Box.test(residuals(m_gam_flow2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_flow2)
## X-squared = 126.85, df = 20, p-value < 2.2e-16Besides the Shapiro-Wilk normality test showing that the residuals aren’t normal, the diagnostic plots look really good. However, the residuals are autocorrelated.
Model with lag terms
Now, we’ll try to deal with the residual autocorrelation. We’ll run a series of models adding 1 and 2 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_gam_flow2_lag1 <- gam(
Chla_log ~ (Year_fct + Flow + StationCode)^2 +
s(Week, bs = "cc", k = 5) +
lag1,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_flow2_lag1))
Box.test(residuals(m_gam_flow2_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_flow2_lag1)
## X-squared = 28.41, df = 20, p-value = 0.1Lag 2
m_gam_flow2_lag2 <- gam(
Chla_log ~ (Year_fct + Flow + StationCode)^2 +
s(Week, bs = "cc", k = 5) +
lag1 +
lag2,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_flow2_lag2))
Box.test(residuals(m_gam_flow2_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_flow2_lag2)
## X-squared = 19.353, df = 20, p-value = 0.499The model with 1 lag term already seems to address the serial autocorrelation, but the lag2 model is even better. Let’s use AIC to see how they compare.
Compare Models
AIC(m_gam_flow2, m_gam_flow2_lag1, m_gam_flow2_lag2)## df AIC
## m_gam_flow2 31.93221 249.4926
## m_gam_flow2_lag1 32.77541 149.3128
## m_gam_flow2_lag2 33.74833 138.0220Again, it looks like the lag2 model has the best fit according to the AIC values. Let’s take a closer look at that one.
Lag 2 model summary
summary(m_gam_flow2_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(Week, bs = "cc",
## k = 5) + lag1 + lag2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.6927620 0.5088819 11.187 < 2e-16 ***
## Year_fct2016 -0.3997008 0.1305545 -3.062 0.002441 **
## Year_fct2017 -0.2694646 0.1205579 -2.235 0.026286 *
## Year_fct2018 -0.2518784 0.1180629 -2.133 0.033858 *
## Year_fct2019 -0.0039353 0.1160004 -0.034 0.972964
## Flow -0.0008715 0.0003024 -2.882 0.004291 **
## StationCodeSTTD -0.5722347 0.1242222 -4.607 6.51e-06 ***
## StationCodeLIB -1.2125300 0.1496448 -8.103 2.34e-14 ***
## StationCodeRVB -1.2321433 0.1526991 -8.069 2.92e-14 ***
## lag1 0.5309148 0.0609174 8.715 3.98e-16 ***
## lag2 -0.1173388 0.0563545 -2.082 0.038339 *
## Year_fct2016:Flow -0.0002401 0.0003722 -0.645 0.519412
## Year_fct2017:Flow 0.0005429 0.0014803 0.367 0.714116
## Year_fct2018:Flow 0.0003463 0.0003167 1.093 0.275306
## Year_fct2019:Flow -0.0003186 0.0002831 -1.125 0.261462
## Year_fct2016:StationCodeSTTD 0.5152892 0.1679207 3.069 0.002386 **
## Year_fct2017:StationCodeSTTD 0.2218232 0.1835982 1.208 0.228105
## Year_fct2018:StationCodeSTTD -0.1434021 0.1587806 -0.903 0.367312
## Year_fct2019:StationCodeSTTD -0.6524006 0.1625663 -4.013 7.91e-05 ***
## Year_fct2016:StationCodeLIB 0.8102839 0.1728554 4.688 4.53e-06 ***
## Year_fct2017:StationCodeLIB 0.3481908 0.1672413 2.082 0.038357 *
## Year_fct2018:StationCodeLIB -0.9897797 0.1795194 -5.513 8.71e-08 ***
## Year_fct2019:StationCodeLIB -0.4536142 0.1601279 -2.833 0.004988 **
## Year_fct2016:StationCodeRVB 0.3165400 0.1623957 1.949 0.052384 .
## Year_fct2017:StationCodeRVB 0.1199746 0.1659188 0.723 0.470295
## Year_fct2018:StationCodeRVB 0.0725750 0.1508970 0.481 0.630965
## Year_fct2019:StationCodeRVB -0.3942603 0.1538694 -2.562 0.010981 *
## Flow:StationCodeSTTD 0.0021542 0.0003131 6.880 4.74e-11 ***
## Flow:StationCodeLIB 0.0017400 0.0003015 5.771 2.31e-08 ***
## Flow:StationCodeRVB 0.0010239 0.0002855 3.586 0.000403 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.454 3 13.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.924 Deviance explained = 93.3%
## -REML = 136.22 Scale est. = 0.084737 n = 284appraise(m_gam_flow2_lag2)
shapiro.test(residuals(m_gam_flow2_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_flow2_lag2)
## W = 0.96549, p-value = 2.591e-06k.check(m_gam_flow2_lag2)## k' edf k-index p-value
## s(Week) 3 2.454221 0.978034 0.3175draw(m_gam_flow2_lag2, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_flow2_lag2, pages = 1, all.terms = TRUE)
anova(m_gam_flow2_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(Week, bs = "cc",
## k = 5) + lag1 + lag2
##
## Parametric Terms:
## df F p-value
## Year_fct 4 3.998 0.00366
## Flow 1 8.307 0.00429
## StationCode 3 25.413 2.15e-14
## lag1 1 75.957 3.98e-16
## lag2 1 4.335 0.03834
## Year_fct:Flow 4 1.764 0.13668
## Year_fct:StationCode 12 7.377 1.34e-11
## Flow:StationCode 3 18.037 1.24e-10
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.454 3.000 13.2 <2e-16The model diagnostics look pretty good. Note that the 2-way
interaction between Year and Flow isn’t significant. We’ll use
m_gam_flow2_lag2 in the model selection process.
rm(m_gam_flow2, m_gam_flow2_lag1)Model 3: GAM Model with 2-way interaction between Station and Year but without Flow
Next we’ll try running a GAM using a two-way interaction between Year and Station but not including flow as a predictor.
Initial Model
m_gam_cat2 <- gam(
Chla_log ~ Year_fct * StationCode + s(Week, bs = "cc", k = 5),
data = df_chla_c2,
method = "REML",
knots = list(Week = c(0, 52))
)Lets look at the model summary and diagnostics:
summary(m_gam_cat2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + s(Week, bs = "cc", k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.52235 0.10443 91.186 < 2e-16 ***
## Year_fct2016 -0.59040 0.15473 -3.816 0.000165 ***
## Year_fct2017 -0.23453 0.14502 -1.617 0.106871
## Year_fct2018 -0.30480 0.14296 -2.132 0.033810 *
## Year_fct2019 -0.15134 0.14296 -1.059 0.290635
## StationCodeSTTD -0.75861 0.14490 -5.236 3.09e-07 ***
## StationCodeLIB -1.80328 0.13971 -12.907 < 2e-16 ***
## StationCodeRVB -1.86323 0.13971 -13.336 < 2e-16 ***
## Year_fct2016:StationCodeSTTD 0.84525 0.21375 3.954 9.57e-05 ***
## Year_fct2017:StationCodeSTTD 0.01425 0.22080 0.065 0.948601
## Year_fct2018:StationCodeSTTD -0.34732 0.20766 -1.673 0.095448 .
## Year_fct2019:StationCodeSTTD -0.84806 0.20352 -4.167 4.04e-05 ***
## Year_fct2016:StationCodeLIB 1.37570 0.20491 6.714 9.42e-11 ***
## Year_fct2017:StationCodeLIB 0.26636 0.19964 1.334 0.183133
## Year_fct2018:StationCodeLIB -1.58046 0.20674 -7.645 2.82e-13 ***
## Year_fct2019:StationCodeLIB -0.54733 0.20374 -2.686 0.007623 **
## Year_fct2016:StationCodeRVB 0.55807 0.20654 2.702 0.007281 **
## Year_fct2017:StationCodeRVB -0.02388 0.19964 -0.120 0.904881
## Year_fct2018:StationCodeRVB -0.01596 0.19679 -0.081 0.935423
## Year_fct2019:StationCodeRVB -0.60593 0.19679 -3.079 0.002267 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.412 3 17.56 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.844 Deviance explained = 85.4%
## -REML = 196.13 Scale est. = 0.17285 n = 324appraise(m_gam_cat2)
shapiro.test(residuals(m_gam_cat2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_cat2)
## W = 0.98188, p-value = 0.0004132k.check(m_gam_cat2)## k' edf k-index p-value
## s(Week) 3 2.412098 0.9343303 0.0975draw(m_gam_cat2, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_cat2, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_cat2))
Box.test(residuals(m_gam_cat2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_cat2)
## X-squared = 285.11, df = 20, p-value < 2.2e-16Besides the Shapiro-Wilk normality test showing that the residuals aren’t normal, the diagnostic plots look really good. However, the residuals are autocorrelated.
Model with lag terms
Now, we’ll try to deal with the residual autocorrelation. We’ll run a series of models adding 1 and 2 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_gam_cat2_lag1 <- gam(
Chla_log ~ Year_fct * StationCode + s(Week, bs = "cc", k = 5) + lag1,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_cat2_lag1))
Box.test(residuals(m_gam_cat2_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_cat2_lag1)
## X-squared = 29.793, df = 20, p-value = 0.07328Lag 2
m_gam_cat2_lag2 <- gam(
Chla_log ~ Year_fct * StationCode + s(Week, bs = "cc", k = 5) + lag1 + lag2,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_cat2_lag2))
Box.test(residuals(m_gam_cat2_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_cat2_lag2)
## X-squared = 15.574, df = 20, p-value = 0.7427The model with 1 lag term already seems to address the serial autocorrelation, but the lag2 model is even better. Let’s use AIC to see how they compare.
Compare Models
AIC(m_gam_cat2, m_gam_cat2_lag1, m_gam_cat2_lag2)## df AIC
## m_gam_cat2 23.73452 374.9887
## m_gam_cat2_lag1 24.37544 200.7548
## m_gam_cat2_lag2 25.49185 184.6154Again, it looks like the lag2 model has the best fit according to the AIC values. Let’s take a closer look at that one.
Lag 2 model summary
summary(m_gam_cat2_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + s(Week, bs = "cc", k = 5) +
## lag1 + lag2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.065399 0.505647 8.040 3.18e-14 ***
## Year_fct2016 -0.304260 0.131386 -2.316 0.021350 *
## Year_fct2017 -0.105866 0.119535 -0.886 0.376628
## Year_fct2018 -0.127785 0.118521 -1.078 0.281958
## Year_fct2019 -0.025808 0.117799 -0.219 0.826753
## StationCodeSTTD -0.258138 0.124366 -2.076 0.038912 *
## StationCodeLIB -0.763750 0.146773 -5.204 3.97e-07 ***
## StationCodeRVB -0.788523 0.150028 -5.256 3.08e-07 ***
## lag1 0.744321 0.059661 12.476 < 2e-16 ***
## lag2 -0.174329 0.060677 -2.873 0.004401 **
## Year_fct2016:StationCodeSTTD 0.385743 0.180505 2.137 0.033530 *
## Year_fct2017:StationCodeSTTD -0.039684 0.184574 -0.215 0.829932
## Year_fct2018:StationCodeSTTD -0.128853 0.171712 -0.750 0.453691
## Year_fct2019:StationCodeSTTD -0.419929 0.172216 -2.438 0.015424 *
## Year_fct2016:StationCodeLIB 0.604215 0.182603 3.309 0.001069 **
## Year_fct2017:StationCodeLIB 0.109151 0.163821 0.666 0.505820
## Year_fct2018:StationCodeLIB -0.717763 0.188519 -3.807 0.000175 ***
## Year_fct2019:StationCodeLIB -0.290952 0.169497 -1.717 0.087251 .
## Year_fct2016:StationCodeRVB 0.234250 0.171967 1.362 0.174320
## Year_fct2017:StationCodeRVB -0.034849 0.163505 -0.213 0.831385
## Year_fct2018:StationCodeRVB 0.007175 0.160844 0.045 0.964456
## Year_fct2019:StationCodeRVB -0.287943 0.163206 -1.764 0.078858 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.134 3 6.191 7.45e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.909 Deviance explained = 91.6%
## -REML = 106.2 Scale est. = 0.10243 n = 284appraise(m_gam_cat2_lag2)
shapiro.test(residuals(m_gam_cat2_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_cat2_lag2)
## W = 0.9567, p-value = 1.809e-07k.check(m_gam_cat2_lag2)## k' edf k-index p-value
## s(Week) 3 2.133943 0.9481868 0.19draw(m_gam_cat2_lag2, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_cat2_lag2, pages = 1, all.terms = TRUE)
anova(m_gam_cat2_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + s(Week, bs = "cc", k = 5) +
## lag1 + lag2
##
## Parametric Terms:
## df F p-value
## Year_fct 4 1.657 0.160485
## StationCode 3 11.575 3.81e-07
## lag1 1 155.647 < 2e-16
## lag2 1 8.255 0.004401
## Year_fct:StationCode 12 3.433 0.000107
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.134 3.000 6.191 7.45e-05The model diagnostics look pretty good but not quite as good as with
the initial model. We’ll use m_gam_cat2_lag2 in the model
selection process.
rm(m_gam_cat2, m_gam_cat2_lag1)Model 4: Linear Model with Flow and 3-way interactions
Let’s try the weekly average model as a linear model with a three-way interaction between Year, Weekly Average Flow, and Station but without the smooth term for week number. Initially, we’ll run the model without accounting for serial autocorrelation.
Initial Model
m_lm_flow3 <- lm(Chla_log ~ Year_fct * Flow * StationCode, data = df_chla_c2)Lets look at the model summary and diagnostics:
summary(m_lm_flow3)##
## Call:
## lm(formula = Chla_log ~ Year_fct * Flow * StationCode, data = df_chla_c2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.97374 -0.19702 -0.02017 0.16626 1.33105
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.7615689 0.1149026 84.955 < 2e-16 ***
## Year_fct2016 -0.7803006 0.1662325 -4.694 4.17e-06 ***
## Year_fct2017 -1.2750955 0.2533518 -5.033 8.59e-07 ***
## Year_fct2018 -0.3716521 0.1542539 -2.409 0.016617 *
## Year_fct2019 -0.1064280 0.1515132 -0.702 0.482984
## Flow -0.0021678 0.0005067 -4.278 2.58e-05 ***
## StationCodeSTTD -1.1482782 0.1453232 -7.902 6.10e-14 ***
## StationCodeLIB -2.0659490 0.1403787 -14.717 < 2e-16 ***
## StationCodeRVB -2.1117161 0.1403787 -15.043 < 2e-16 ***
## Year_fct2016:Flow 0.0028066 0.0008103 3.464 0.000615 ***
## Year_fct2017:Flow 0.0498038 0.0128457 3.877 0.000131 ***
## Year_fct2018:Flow 0.0004445 0.0007085 0.627 0.530893
## Year_fct2019:Flow -0.0000748 0.0006226 -0.120 0.904452
## Year_fct2016:StationCodeSTTD 1.2499596 0.2136261 5.851 1.34e-08 ***
## Year_fct2017:StationCodeSTTD 1.7155377 0.3407645 5.034 8.53e-07 ***
## Year_fct2018:StationCodeSTTD -0.2869259 0.2037980 -1.408 0.160256
## Year_fct2019:StationCodeSTTD -0.9572118 0.2007104 -4.769 2.96e-06 ***
## Year_fct2016:StationCodeLIB 1.7614390 0.2054580 8.573 6.59e-16 ***
## Year_fct2017:StationCodeLIB 1.3626781 0.2886090 4.722 3.68e-06 ***
## Year_fct2018:StationCodeLIB -1.8236205 0.2044106 -8.921 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.6043253 0.1995777 -3.028 0.002688 **
## Year_fct2016:StationCodeRVB 0.9752324 0.2069080 4.713 3.82e-06 ***
## Year_fct2017:StationCodeRVB 1.0593695 0.2886090 3.671 0.000289 ***
## Year_fct2018:StationCodeRVB 0.0493574 0.1940940 0.254 0.799451
## Year_fct2019:StationCodeRVB -0.7029752 0.1938728 -3.626 0.000341 ***
## Flow:StationCodeSTTD 0.0049297 0.0007794 6.325 9.84e-10 ***
## Flow:StationCodeLIB 0.0028763 0.0007599 3.785 0.000187 ***
## Flow:StationCodeRVB 0.0018144 0.0007599 2.388 0.017609 *
## Year_fct2016:Flow:StationCodeSTTD -0.0043111 0.0011917 -3.618 0.000352 ***
## Year_fct2017:Flow:StationCodeSTTD -0.0352894 0.0145777 -2.421 0.016114 *
## Year_fct2018:Flow:StationCodeSTTD -0.0007185 0.0010472 -0.686 0.493211
## Year_fct2019:Flow:StationCodeSTTD -0.0012877 0.0009279 -1.388 0.166310
## Year_fct2016:Flow:StationCodeLIB -0.0032468 0.0011734 -2.767 0.006027 **
## Year_fct2017:Flow:StationCodeLIB -0.0483336 0.0131090 -3.687 0.000272 ***
## Year_fct2018:Flow:StationCodeLIB 0.0017105 0.0010495 1.630 0.104231
## Year_fct2019:Flow:StationCodeLIB -0.0008137 0.0009237 -0.881 0.379107
## Year_fct2016:Flow:StationCodeRVB -0.0031157 0.0011756 -2.650 0.008493 **
## Year_fct2017:Flow:StationCodeRVB -0.0472314 0.0131090 -3.603 0.000371 ***
## Year_fct2018:Flow:StationCodeRVB -0.0007736 0.0010268 -0.753 0.451830
## Year_fct2019:Flow:StationCodeRVB 0.0007609 0.0009095 0.837 0.403500
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3591 on 284 degrees of freedom
## Multiple R-squared: 0.8977, Adjusted R-squared: 0.8837
## F-statistic: 63.92 on 39 and 284 DF, p-value: < 2.2e-16df_chla_c2 %>% plot_lm_diag(Chla_log, m_lm_flow3)
shapiro.test(residuals(m_lm_flow3))##
## Shapiro-Wilk normality test
##
## data: residuals(m_lm_flow3)
## W = 0.98191, p-value = 0.0004191acf(residuals(m_lm_flow3))
Box.test(residuals(m_lm_flow3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_flow3)
## X-squared = 151.84, df = 20, p-value < 2.2e-16The residuals deviate from a normal distribution according to visual inspection and the Shapiro-Wilk normality test. Also, model definitely has residual autocorrelation as indicated by the ACF plot and the Box-Ljung test.
Model with lag terms
Now, we’ll try to deal with the residual autocorrelation and the non-normal residuals. We’ll run a series of linear models adding 1 and 2 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_lm_flow3_lag1 <- df_chla_c2_lag %>%
drop_na(Chla_log, lag1) %>%
lm(Chla_log ~ Year_fct * Flow * StationCode + lag1, data = .)
acf(residuals(m_lm_flow3_lag1))
Box.test(residuals(m_lm_flow3_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_flow3_lag1)
## X-squared = 25.245, df = 20, p-value = 0.1922Lag 2
m_lm_flow3_lag2 <- df_chla_c2_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
lm(Chla_log ~ Year_fct * Flow * StationCode + lag1 + lag2, data = .)
acf(residuals(m_lm_flow3_lag2))
Box.test(residuals(m_lm_flow3_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_flow3_lag2)
## X-squared = 13.923, df = 20, p-value = 0.8344The model with 1 lag term already has better ACF and Box-Ljung test results than the initial model. Let’s use AIC to see how they compare.
Compare Models
AIC(m_lm_flow3, m_lm_flow3_lag1, m_lm_flow3_lag2)## df AIC
## m_lm_flow3 41 295.0642
## m_lm_flow3_lag1 42 179.7469
## m_lm_flow3_lag2 43 172.2352Again, it looks like the lag2 model has the best fit according to the AIC values. Let’s take a closer look at that one.
Lag 2 model summary
summary(m_lm_flow3_lag2)##
## Call:
## lm(formula = Chla_log ~ Year_fct * Flow * StationCode + lag1 +
## lag2, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.03686 -0.14697 -0.00613 0.11959 1.00420
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4975717 0.5381025 10.217 < 2e-16 ***
## Year_fct2016 -0.5272351 0.1585453 -3.325 0.001020 **
## Year_fct2017 -0.7714772 0.2455699 -3.142 0.001890 **
## Year_fct2018 -0.2640250 0.1445395 -1.827 0.068982 .
## Year_fct2019 -0.0528143 0.1412730 -0.374 0.708846
## Flow -0.0014995 0.0004555 -3.292 0.001143 **
## StationCodeSTTD -0.6512369 0.1458385 -4.465 1.23e-05 ***
## StationCodeLIB -1.2126820 0.1650472 -7.347 3.09e-12 ***
## StationCodeRVB -1.2247878 0.1681936 -7.282 4.60e-12 ***
## lag1 0.5301472 0.0652367 8.127 2.28e-14 ***
## lag2 -0.0907567 0.0601615 -1.509 0.132718
## Year_fct2016:Flow 0.0017664 0.0007138 2.475 0.014027 *
## Year_fct2017:Flow 0.0299815 0.0135893 2.206 0.028306 *
## Year_fct2018:Flow 0.0005670 0.0006222 0.911 0.363076
## Year_fct2019:Flow 0.0001914 0.0005492 0.348 0.727804
## Year_fct2016:StationCodeSTTD 0.7826555 0.2073799 3.774 0.000202 ***
## Year_fct2017:StationCodeSTTD 1.0521849 0.3299013 3.189 0.001614 **
## Year_fct2018:StationCodeSTTD -0.0991403 0.1918570 -0.517 0.605809
## Year_fct2019:StationCodeSTTD -0.5693466 0.1949134 -2.921 0.003819 **
## Year_fct2016:StationCodeLIB 1.0737971 0.2106794 5.097 6.96e-07 ***
## Year_fct2017:StationCodeLIB 0.8020328 0.2761937 2.904 0.004026 **
## Year_fct2018:StationCodeLIB -1.1176227 0.2212095 -5.052 8.60e-07 ***
## Year_fct2019:StationCodeLIB -0.3518422 0.1901208 -1.851 0.065442 .
## Year_fct2016:StationCodeRVB 0.6035354 0.1967708 3.067 0.002406 **
## Year_fct2017:StationCodeRVB 0.5994637 0.2727511 2.198 0.028908 *
## Year_fct2018:StationCodeRVB 0.1036231 0.1796164 0.577 0.564534
## Year_fct2019:StationCodeRVB -0.3803675 0.1848640 -2.058 0.040705 *
## Flow:StationCodeSTTD 0.0029446 0.0007232 4.072 6.32e-05 ***
## Flow:StationCodeLIB 0.0021381 0.0006779 3.154 0.001814 **
## Flow:StationCodeRVB 0.0012279 0.0006719 1.827 0.068869 .
## Year_fct2016:Flow:StationCodeSTTD -0.0025056 0.0010605 -2.363 0.018942 *
## Year_fct2017:Flow:StationCodeSTTD -0.0193301 0.0149430 -1.294 0.197039
## Year_fct2018:Flow:StationCodeSTTD -0.0006183 0.0009282 -0.666 0.505995
## Year_fct2019:Flow:StationCodeSTTD -0.0007664 0.0008211 -0.933 0.351515
## Year_fct2016:Flow:StationCodeLIB -0.0020124 0.0010295 -1.955 0.051756 .
## Year_fct2017:Flow:StationCodeLIB -0.0296984 0.0137441 -2.161 0.031692 *
## Year_fct2018:Flow:StationCodeLIB 0.0010813 0.0009386 1.152 0.250423
## Year_fct2019:Flow:StationCodeLIB -0.0009727 0.0008129 -1.197 0.232628
## Year_fct2016:Flow:StationCodeRVB -0.0016700 0.0010336 -1.616 0.107443
## Year_fct2017:Flow:StationCodeRVB -0.0286477 0.0137460 -2.084 0.038203 *
## Year_fct2018:Flow:StationCodeRVB -0.0007750 0.0008989 -0.862 0.389463
## Year_fct2019:Flow:StationCodeRVB 0.0001112 0.0008000 0.139 0.889543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3051 on 242 degrees of freedom
## Multiple R-squared: 0.929, Adjusted R-squared: 0.917
## F-statistic: 77.26 on 41 and 242 DF, p-value: < 2.2e-16df_chla_c2_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
plot_lm_diag(Chla_log, m_lm_flow3_lag2)
shapiro.test(residuals(m_lm_flow3_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_lm_flow3_lag2)
## W = 0.9644, p-value = 1.828e-06Anova(
m_lm_flow3_lag2,
type = 3,
contrasts = list(topic = contr.sum, sys = contr.sum)
)## Anova Table (Type III tests)
##
## Response: Chla_log
## Sum Sq Df F value Pr(>F)
## (Intercept) 9.7165 1 104.3787 < 2.2e-16 ***
## Year_fct 1.8414 4 4.9454 0.0007538 ***
## Flow 1.0089 1 10.8376 0.0011431 **
## StationCode 5.8077 3 20.7963 5.065e-12 ***
## lag1 6.1477 1 66.0404 2.279e-14 ***
## lag2 0.2118 1 2.2757 0.1327181
## Year_fct:Flow 1.0592 4 2.8446 0.0247760 *
## Year_fct:StationCode 7.2387 12 6.4801 5.404e-10 ***
## Flow:StationCode 1.7700 3 6.3378 0.0003741 ***
## Year_fct:Flow:StationCode 2.3512 12 2.1048 0.0172994 *
## Residuals 22.5276 242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model diagnostics look okay, but not as good as with the GAM
models. We’ll use m_lm_flow3_lag2 in the model selection
process.
rm(m_lm_flow3, m_lm_flow3_lag1)Model 5: Linear Model with Flow and 2-way interactions
Let’s try a linear model using all two-way interactions between Year, Weekly Average Flow, and Station. Initially, we’ll run the model without accounting for serial autocorrelation.
Initial Model
m_lm_flow2 <- lm(
Chla_log ~ (Year_fct + Flow + StationCode)^2,
data = df_chla_c2
)Lets look at the model summary and diagnostics:
summary(m_lm_flow2)##
## Call:
## lm(formula = Chla_log ~ (Year_fct + Flow + StationCode)^2, data = df_chla_c2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99534 -0.22307 -0.03153 0.20115 1.36296
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.6943716 0.1094244 88.594 < 2e-16 ***
## Year_fct2016 -0.4689162 0.1527744 -3.069 0.002344 **
## Year_fct2017 -0.4662132 0.1463633 -3.185 0.001600 **
## Year_fct2018 -0.3651423 0.1429793 -2.554 0.011156 *
## Year_fct2019 -0.0701638 0.1420134 -0.494 0.621627
## Flow -0.0016930 0.0003628 -4.666 4.65e-06 ***
## StationCodeSTTD -1.0540539 0.1423626 -7.404 1.38e-12 ***
## StationCodeLIB -1.9989312 0.1385827 -14.424 < 2e-16 ***
## StationCodeRVB -2.0471430 0.1385528 -14.775 < 2e-16 ***
## Year_fct2016:Flow 0.0001952 0.0004545 0.430 0.667779
## Year_fct2017:Flow 0.0036316 0.0018887 1.923 0.055456 .
## Year_fct2018:Flow 0.0004829 0.0004020 1.201 0.230592
## Year_fct2019:Flow -0.0003297 0.0003562 -0.926 0.355401
## Year_fct2016:StationCodeSTTD 0.8736753 0.2027976 4.308 2.24e-05 ***
## Year_fct2017:StationCodeSTTD 0.5774444 0.2236000 2.582 0.010289 *
## Year_fct2018:StationCodeSTTD -0.3143766 0.1963384 -1.601 0.110400
## Year_fct2019:StationCodeSTTD -1.0367475 0.1929058 -5.374 1.56e-07 ***
## Year_fct2016:StationCodeLIB 1.4274612 0.1961874 7.276 3.10e-12 ***
## Year_fct2017:StationCodeLIB 0.6143627 0.2080527 2.953 0.003400 **
## Year_fct2018:StationCodeLIB -1.7055558 0.1958498 -8.708 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.6916575 0.1940127 -3.565 0.000424 ***
## Year_fct2016:StationCodeRVB 0.6391321 0.1974165 3.237 0.001343 **
## Year_fct2017:StationCodeRVB 0.2879781 0.2080017 1.384 0.167249
## Year_fct2018:StationCodeRVB 0.0149704 0.1882767 0.080 0.936678
## Year_fct2019:StationCodeRVB -0.6628422 0.1878057 -3.529 0.000483 ***
## Flow:StationCodeSTTD 0.0035768 0.0003609 9.910 < 2e-16 ***
## Flow:StationCodeLIB 0.0024150 0.0003682 6.560 2.40e-10 ***
## Flow:StationCodeRVB 0.0015361 0.0003566 4.308 2.24e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3865 on 296 degrees of freedom
## Multiple R-squared: 0.8765, Adjusted R-squared: 0.8652
## F-statistic: 77.81 on 27 and 296 DF, p-value: < 2.2e-16df_chla_c2 %>% plot_lm_diag(Chla_log, m_lm_flow2)
shapiro.test(residuals(m_lm_flow2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_lm_flow2)
## W = 0.9802, p-value = 0.0001911acf(residuals(m_lm_flow2))
Box.test(residuals(m_lm_flow2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_flow2)
## X-squared = 239.47, df = 20, p-value < 2.2e-16The residuals deviate from a normal distribution according to visual inspection and the Shapiro-Wilk normality test. Also, model definitely has residual autocorrelation as indicated by the ACF plot and the Box-Ljung test.
Model with lag terms
Now, we’ll try to deal with the residual autocorrelation and the non-normal residuals. We’ll run a series of linear models adding 1 and 2 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_lm_flow2_lag1 <- df_chla_c2_lag %>%
drop_na(Chla_log, lag1) %>%
lm(Chla_log ~ (Year_fct + Flow + StationCode)^2 + lag1, data = .)
acf(residuals(m_lm_flow2_lag1))
Box.test(residuals(m_lm_flow2_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_flow2_lag1)
## X-squared = 29.04, df = 20, p-value = 0.08698Lag 2
m_lm_flow2_lag2 <- df_chla_c2_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
lm(Chla_log ~ (Year_fct + Flow + StationCode)^2 + lag1 + lag2, data = .)
acf(residuals(m_lm_flow2_lag2))
Box.test(residuals(m_lm_flow2_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_flow2_lag2)
## X-squared = 16.346, df = 20, p-value = 0.6949The model with 2 lag terms seems to be okay in terms of serial autocorrelation. Let’s use AIC to see how they compare.
Compare Models
AIC(m_lm_flow2, m_lm_flow2_lag1, m_lm_flow2_lag2)## df AIC
## m_lm_flow2 29 332.1787
## m_lm_flow2_lag1 30 181.2203
## m_lm_flow2_lag2 31 176.4292Again, it looks like the lag2 model has the best fit according to the AIC values. Let’s take a closer look at that one.
Lag 2 model summary
summary(m_lm_flow2_lag2)##
## Call:
## lm(formula = Chla_log ~ (Year_fct + Flow + StationCode)^2 + lag1 +
## lag2, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.20397 -0.15989 -0.01927 0.11946 0.96043
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4998924 0.4956888 9.078 < 2e-16 ***
## Year_fct2016 -0.2943646 0.1351127 -2.179 0.030277 *
## Year_fct2017 -0.2628023 0.1284627 -2.046 0.041811 *
## Year_fct2018 -0.2050746 0.1259745 -1.628 0.104785
## Year_fct2019 0.0077718 0.1241524 0.063 0.950135
## Flow -0.0010099 0.0003095 -3.263 0.001252 **
## StationCodeSTTD -0.4756410 0.1317070 -3.611 0.000367 ***
## StationCodeLIB -0.9610523 0.1524165 -6.305 1.27e-09 ***
## StationCodeRVB -0.9695374 0.1550963 -6.251 1.71e-09 ***
## lag1 0.6313481 0.0631115 10.004 < 2e-16 ***
## lag2 -0.0943825 0.0599941 -1.573 0.116918
## Year_fct2016:Flow 0.0002841 0.0003772 0.753 0.451999
## Year_fct2017:Flow 0.0014882 0.0015620 0.953 0.341608
## Year_fct2018:Flow 0.0004459 0.0003377 1.320 0.187957
## Year_fct2019:Flow -0.0001049 0.0002986 -0.351 0.725789
## Year_fct2016:StationCodeSTTD 0.4590973 0.1797473 2.554 0.011230 *
## Year_fct2017:StationCodeSTTD 0.2611445 0.1965932 1.328 0.185255
## Year_fct2018:StationCodeSTTD -0.1222489 0.1703977 -0.717 0.473766
## Year_fct2019:StationCodeSTTD -0.5428917 0.1736903 -3.126 0.001980 **
## Year_fct2016:StationCodeLIB 0.7212841 0.1828336 3.945 0.000103 ***
## Year_fct2017:StationCodeLIB 0.3150986 0.1796582 1.754 0.080657 .
## Year_fct2018:StationCodeLIB -0.8442328 0.1911986 -4.415 1.49e-05 ***
## Year_fct2019:StationCodeLIB -0.3984735 0.1712576 -2.327 0.020766 *
## Year_fct2016:StationCodeRVB 0.3436105 0.1728820 1.988 0.047936 *
## Year_fct2017:StationCodeRVB 0.1211150 0.1783372 0.679 0.497671
## Year_fct2018:StationCodeRVB 0.0466405 0.1619971 0.288 0.773651
## Year_fct2019:StationCodeRVB -0.3254659 0.1649252 -1.973 0.049532 *
## Flow:StationCodeSTTD 0.0018283 0.0003316 5.514 8.62e-08 ***
## Flow:StationCodeLIB 0.0014442 0.0003202 4.511 9.87e-06 ***
## Flow:StationCodeRVB 0.0007276 0.0003025 2.405 0.016874 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.313 on 254 degrees of freedom
## Multiple R-squared: 0.9216, Adjusted R-squared: 0.9127
## F-statistic: 103 on 29 and 254 DF, p-value: < 2.2e-16df_chla_c2_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
plot_lm_diag(Chla_log, m_lm_flow2_lag2)
shapiro.test(residuals(m_lm_flow2_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_lm_flow2_lag2)
## W = 0.96256, p-value = 1.027e-06Anova(
m_lm_flow2_lag2,
type = 3,
contrasts = list(topic = contr.sum, sys = contr.sum)
)## Anova Table (Type III tests)
##
## Response: Chla_log
## Sum Sq Df F value Pr(>F)
## (Intercept) 8.0720 1 82.4112 < 2.2e-16 ***
## Year_fct 0.9908 4 2.5290 0.041132 *
## Flow 1.0431 1 10.6493 0.001252 **
## StationCode 4.5471 3 15.4744 2.824e-09 ***
## lag1 9.8020 1 100.0738 < 2.2e-16 ***
## lag2 0.2424 1 2.4749 0.116918
## Year_fct:Flow 0.5183 4 1.3228 0.261880
## Year_fct:StationCode 5.7798 12 4.9174 2.643e-07 ***
## Flow:StationCode 3.4705 3 11.8107 2.875e-07 ***
## Residuals 24.8788 254
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model diagnostics look somewhat worse than those for the 3-way
interaction model. Note that the 2-way interaction between Year and Flow
isn’t significant. We’ll use m_lm_flow2_lag2 in the model
selection process.
rm(m_lm_flow2, m_lm_flow2_lag1)Model 6: Linear Model with 2-way interaction between Station and Year but without Flow
We’ll try running a linear model using a two-way interaction between Year and Station but not including flow as a predictor. Initially, we’ll run the model without accounting for serial autocorrelation.
Initial Model
m_lm_cat2 <- lm(Chla_log ~ Year_fct * StationCode, data = df_chla_c2)Lets look at the model summary and diagnostics:
summary(m_lm_cat2)##
## Call:
## lm(formula = Chla_log ~ Year_fct * StationCode, data = df_chla_c2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.35952 -0.25178 -0.05221 0.27546 1.31786
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.45473 0.11258 83.981 < 2e-16 ***
## Year_fct2016 -0.40045 0.16480 -2.430 0.015684 *
## Year_fct2017 -0.19510 0.15686 -1.244 0.214519
## Year_fct2018 -0.26862 0.15473 -1.736 0.083569 .
## Year_fct2019 -0.11515 0.15473 -0.744 0.457316
## StationCodeSTTD -0.75644 0.15686 -4.822 2.25e-06 ***
## StationCodeLIB -1.74964 0.15104 -11.584 < 2e-16 ***
## StationCodeRVB -1.80959 0.15104 -11.980 < 2e-16 ***
## Year_fct2016:StationCodeSTTD 0.84307 0.23146 3.642 0.000317 ***
## Year_fct2017:StationCodeSTTD 0.10665 0.23835 0.447 0.654872
## Year_fct2018:StationCodeSTTD -0.35014 0.22440 -1.560 0.119719
## Year_fct2019:StationCodeSTTD -0.88642 0.22033 -4.023 7.25e-05 ***
## Year_fct2016:StationCodeLIB 1.38154 0.22188 6.227 1.58e-09 ***
## Year_fct2017:StationCodeLIB 0.21272 0.21604 0.985 0.325578
## Year_fct2018:StationCodeLIB -1.72136 0.22289 -7.723 1.66e-13 ***
## Year_fct2019:StationCodeLIB -0.62118 0.22038 -2.819 0.005138 **
## Year_fct2016:StationCodeRVB 0.57425 0.22355 2.569 0.010683 *
## Year_fct2017:StationCodeRVB -0.07752 0.21604 -0.359 0.719984
## Year_fct2018:StationCodeRVB -0.06960 0.21295 -0.327 0.744018
## Year_fct2019:StationCodeRVB -0.65957 0.21295 -3.097 0.002135 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4503 on 304 degrees of freedom
## Multiple R-squared: 0.8278, Adjusted R-squared: 0.817
## F-statistic: 76.92 on 19 and 304 DF, p-value: < 2.2e-16df_chla_c2 %>% plot_lm_diag(Chla_log, m_lm_cat2)
shapiro.test(residuals(m_lm_cat2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_lm_cat2)
## W = 0.98776, p-value = 0.007767acf(residuals(m_lm_cat2))
Box.test(residuals(m_lm_cat2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2)
## X-squared = 374.61, df = 20, p-value < 2.2e-16Besides the Shapiro-Wilk normality test showing that the residuals aren’t normal, the diagnostic plots look pretty good. However, the residuals are autocorrelated.
Model with lag terms
Now, we’ll try to deal with the residual autocorrelation and the non-normal residuals. We’ll run a series of linear models adding 1 and 2 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_lm_cat2_lag1 <- df_chla_c2_lag %>%
drop_na(Chla_log, lag1) %>%
lm(Chla_log ~ Year_fct * StationCode + lag1, data = .)
acf(residuals(m_lm_cat2_lag1))
Box.test(residuals(m_lm_cat2_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2_lag1)
## X-squared = 33.232, df = 20, p-value = 0.03182Lag 2
m_lm_cat2_lag2 <- df_chla_c2_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
lm(Chla_log ~ Year_fct * StationCode + lag1 + lag2, data = .)
acf(residuals(m_lm_cat2_lag2))
Box.test(residuals(m_lm_cat2_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2_lag2)
## X-squared = 16.358, df = 20, p-value = 0.6942The model with 2 lag terms seems to be okay in terms of serial autocorrelation. Let’s use AIC to see how they compare.
Compare Models
AIC(m_lm_cat2, m_lm_cat2_lag1, m_lm_cat2_lag2)## df AIC
## m_lm_cat2 21 423.8688
## m_lm_cat2_lag1 22 212.6171
## m_lm_cat2_lag2 23 201.4048Again, it looks like the lag2 model has the best fit according to the AIC values. Let’s take a closer look at that one.
Lag 2 model summary
summary(m_lm_cat2_lag2)##
## Call:
## lm(formula = Chla_log ~ Year_fct * StationCode + lag1 + lag2,
## data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.32499 -0.17269 -0.01809 0.14998 0.85380
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.372162 0.481718 7.000 2.14e-11 ***
## Year_fct2016 -0.181269 0.131159 -1.382 0.168132
## Year_fct2017 -0.077050 0.123280 -0.625 0.532516
## Year_fct2018 -0.087196 0.122108 -0.714 0.475808
## Year_fct2019 0.004655 0.121586 0.038 0.969491
## StationCodeSTTD -0.204739 0.127949 -1.600 0.110768
## StationCodeLIB -0.610005 0.145883 -4.181 3.95e-05 ***
## StationCodeRVB -0.628436 0.148825 -4.223 3.33e-05 ***
## lag1 0.797550 0.060317 13.223 < 2e-16 ***
## lag2 -0.158015 0.061978 -2.550 0.011357 *
## Year_fct2016:StationCodeSTTD 0.330530 0.186247 1.775 0.077111 .
## Year_fct2017:StationCodeSTTD -0.018475 0.190484 -0.097 0.922808
## Year_fct2018:StationCodeSTTD -0.126243 0.177047 -0.713 0.476452
## Year_fct2019:StationCodeSTTD -0.386387 0.177700 -2.174 0.030572 *
## Year_fct2016:StationCodeLIB 0.523514 0.186858 2.802 0.005463 **
## Year_fct2017:StationCodeLIB 0.063625 0.169164 0.376 0.707135
## Year_fct2018:StationCodeLIB -0.674882 0.194480 -3.470 0.000608 ***
## Year_fct2019:StationCodeLIB -0.297316 0.175374 -1.695 0.091202 .
## Year_fct2016:StationCodeRVB 0.220236 0.177556 1.240 0.215946
## Year_fct2017:StationCodeRVB -0.061200 0.169078 -0.362 0.717671
## Year_fct2018:StationCodeRVB -0.025173 0.166265 -0.151 0.879774
## Year_fct2019:StationCodeRVB -0.281198 0.168751 -1.666 0.096840 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3312 on 262 degrees of freedom
## Multiple R-squared: 0.9095, Adjusted R-squared: 0.9022
## F-statistic: 125.3 on 21 and 262 DF, p-value: < 2.2e-16df_chla_c2_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
plot_lm_diag(Chla_log, m_lm_cat2_lag2)
shapiro.test(residuals(m_lm_cat2_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_lm_cat2_lag2)
## W = 0.95081, p-value = 3.589e-08Anova(
m_lm_cat2_lag2,
type = 3,
contrasts = list(topic = contr.sum, sys = contr.sum)
)## Anova Table (Type III tests)
##
## Response: Chla_log
## Sum Sq Df F value Pr(>F)
## (Intercept) 5.3755 1 49.0040 2.140e-11 ***
## Year_fct 0.3048 4 0.6947 0.596226
## StationCode 2.5266 3 7.6776 6.165e-05 ***
## lag1 19.1789 1 174.8376 < 2.2e-16 ***
## lag2 0.7130 1 6.5002 0.011357 *
## Year_fct:StationCode 3.5215 12 2.6752 0.002074 **
## Residuals 28.7402 262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model diagnostics don’t look that great. However, we’ll use
m_lm_cat2_lag2 in the model selection process.
rm(m_lm_cat2, m_lm_cat2_lag1)Model 7: GAM model using smooths for Flow
Finally, we’ll try running a GAM model using smooths for weekly average flow by Station and Year and a smooth term for week number to account for seasonality. Initially, we’ll run the model without accounting for serial autocorrelation.
Initial Model
m_gam_sflow <- gam(
Chla_log ~ s(Flow, by = StationCode) +
s(Flow, by = Year_fct) +
Year_fct * StationCode +
s(Week, bs = "cc", k = 5),
data = df_chla_c2,
method = "REML",
knots = list(Week = c(0, 52))
)Lets look at the model summary and diagnostics:
summary(m_gam_sflow)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ s(Flow, by = StationCode) + s(Flow, by = Year_fct) +
## Year_fct * StationCode + s(Week, bs = "cc", k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.67880 0.08731 110.859 < 2e-16 ***
## Year_fct2016 -0.62257 0.13173 -4.726 3.58e-06 ***
## Year_fct2017 -0.41450 0.13578 -3.053 0.002477 **
## Year_fct2018 -0.38170 0.11622 -3.284 0.001148 **
## Year_fct2019 -0.16065 0.11673 -1.376 0.169815
## StationCodeSTTD -0.85927 0.12346 -6.960 2.29e-11 ***
## StationCodeLIB -1.94019 0.11878 -16.334 < 2e-16 ***
## StationCodeRVB -2.02615 0.11644 -17.400 < 2e-16 ***
## Year_fct2016:StationCodeSTTD 0.80733 0.17935 4.502 9.79e-06 ***
## Year_fct2017:StationCodeSTTD 0.35252 0.19820 1.779 0.076356 .
## Year_fct2018:StationCodeSTTD -0.21879 0.16956 -1.290 0.197966
## Year_fct2019:StationCodeSTTD -0.89203 0.17236 -5.175 4.26e-07 ***
## Year_fct2016:StationCodeLIB 1.36901 0.17332 7.899 5.92e-14 ***
## Year_fct2017:StationCodeLIB 0.57926 0.18074 3.205 0.001502 **
## Year_fct2018:StationCodeLIB -1.56010 0.16797 -9.288 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.52393 0.17176 -3.050 0.002497 **
## Year_fct2016:StationCodeRVB 0.57276 0.17177 3.335 0.000966 ***
## Year_fct2017:StationCodeRVB 0.24077 0.17725 1.358 0.175416
## Year_fct2018:StationCodeRVB 0.08339 0.16018 0.521 0.603050
## Year_fct2019:StationCodeRVB -0.59061 0.16021 -3.686 0.000271 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Flow):StationCodeRD22 1.000e+00 1.0000188 0.293 0.58892
## s(Flow):StationCodeSTTD 2.990e+00 3.7013885 4.651 0.00172 **
## s(Flow):StationCodeLIB 2.416e+00 3.0027006 2.616 0.05135 .
## s(Flow):StationCodeRVB 1.000e+00 1.0001210 0.255 0.61427
## s(Flow):Year_fct2015 1.000e+00 1.0000747 0.228 0.63362
## s(Flow):Year_fct2016 1.697e+00 2.1101506 0.386 0.62924
## s(Flow):Year_fct2017 4.749e-05 0.0000943 0.065 0.99803
## s(Flow):Year_fct2018 1.000e+00 1.0000992 0.027 0.86978
## s(Flow):Year_fct2019 1.001e+00 1.0029186 0.200 0.65737
## s(Week) 2.763e+00 3.0000000 33.756 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 103/104
## R-sq.(adj) = 0.903 Deviance explained = 91.3%
## -REML = 141.49 Scale est. = 0.10764 n = 324appraise(m_gam_sflow)
shapiro.test(residuals(m_gam_sflow))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_sflow)
## W = 0.98081, p-value = 0.0002518k.check(m_gam_sflow)## k' edf k-index p-value
## s(Flow):StationCodeRD22 9 1.000009e+00 1.0488673 0.7825
## s(Flow):StationCodeSTTD 9 2.989703e+00 1.0488673 0.7875
## s(Flow):StationCodeLIB 9 2.415561e+00 1.0488673 0.8225
## s(Flow):StationCodeRVB 9 1.000061e+00 1.0488673 0.8075
## s(Flow):Year_fct2015 9 1.000038e+00 1.0488673 0.7875
## s(Flow):Year_fct2016 9 1.696835e+00 1.0488673 0.7900
## s(Flow):Year_fct2017 9 4.748618e-05 1.0488673 0.8025
## s(Flow):Year_fct2018 9 1.000050e+00 1.0488673 0.8300
## s(Flow):Year_fct2019 9 1.001492e+00 1.0488673 0.7875
## s(Week) 3 2.763257e+00 0.9642118 0.2825concurvity(m_gam_sflow, full = FALSE)$worst## para s(Flow):StationCodeRD22 s(Flow):StationCodeSTTD s(Flow):StationCodeLIB
## para 1.000000e+00 2.561680e-01 7.487812e-02 7.679159e-02
## s(Flow):StationCodeRD22 2.561680e-01 1.000000e+00 1.074363e-26 9.371473e-27
## s(Flow):StationCodeSTTD 7.487812e-02 1.101331e-26 1.000000e+00 4.131051e-28
## s(Flow):StationCodeLIB 7.679159e-02 9.271762e-27 5.679579e-28 1.000000e+00
## s(Flow):StationCodeRVB 8.759218e-02 8.234788e-27 5.735916e-28 4.420274e-28
## s(Flow):Year_fct2015 2.252320e-01 4.236437e-01 2.465837e-01 2.573537e-01
## s(Flow):Year_fct2016 1.881476e-01 6.035297e-01 1.543428e-01 2.446523e-01
## s(Flow):Year_fct2017 1.882716e-01 3.384493e-01 6.714084e-02 1.873644e-01
## s(Flow):Year_fct2018 1.937259e-01 3.067056e-01 1.811522e-01 1.938793e-01
## s(Flow):Year_fct2019 2.011684e-01 8.259018e-01 2.906021e-01 2.818269e-01
## s(Week) 9.912595e-31 1.512455e-01 4.081059e-02 6.396698e-02
## s(Flow):StationCodeRVB s(Flow):Year_fct2015 s(Flow):Year_fct2016 s(Flow):Year_fct2017
## para 8.759218e-02 2.252320e-01 1.881476e-01 1.882716e-01
## s(Flow):StationCodeRD22 8.837323e-27 4.236437e-01 6.035297e-01 3.384493e-01
## s(Flow):StationCodeSTTD 3.737793e-28 2.465837e-01 1.543428e-01 6.714084e-02
## s(Flow):StationCodeLIB 2.333734e-28 2.573537e-01 2.446523e-01 1.873644e-01
## s(Flow):StationCodeRVB 1.000000e+00 2.449126e-01 1.563972e-01 1.839605e-01
## s(Flow):Year_fct2015 2.449126e-01 1.000000e+00 2.238503e-26 4.918031e-14
## s(Flow):Year_fct2016 1.563972e-01 2.043733e-26 1.000000e+00 1.534799e-14
## s(Flow):Year_fct2017 1.839605e-01 1.259481e-13 1.595872e-14 1.000000e+00
## s(Flow):Year_fct2018 1.810226e-01 1.362544e-26 4.730280e-27 1.875936e-14
## s(Flow):Year_fct2019 2.894419e-01 1.200993e-25 3.535738e-27 1.950144e-14
## s(Week) 6.093752e-02 1.188708e-01 1.887421e-01 1.683704e-01
## s(Flow):Year_fct2018 s(Flow):Year_fct2019 s(Week)
## para 1.937259e-01 2.011684e-01 4.547140e-31
## s(Flow):StationCodeRD22 3.067056e-01 8.259018e-01 1.512455e-01
## s(Flow):StationCodeSTTD 1.811522e-01 2.906021e-01 4.081059e-02
## s(Flow):StationCodeLIB 1.938793e-01 2.818269e-01 6.396698e-02
## s(Flow):StationCodeRVB 1.810226e-01 2.894419e-01 6.093752e-02
## s(Flow):Year_fct2015 1.358798e-26 1.248024e-25 1.188708e-01
## s(Flow):Year_fct2016 6.224435e-27 2.970225e-27 1.887421e-01
## s(Flow):Year_fct2017 5.521515e-14 6.662451e-14 1.683704e-01
## s(Flow):Year_fct2018 1.000000e+00 3.023201e-27 1.080063e-01
## s(Flow):Year_fct2019 3.103856e-27 1.000000e+00 9.053392e-02
## s(Week) 1.080063e-01 9.053392e-02 1.000000e+00draw(m_gam_sflow, select = 10, residuals = TRUE, rug = FALSE)
plot(m_gam_sflow, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_sflow))
Box.test(residuals(m_gam_sflow), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_sflow)
## X-squared = 118.48, df = 20, p-value = 5.551e-16The diagnostic plots look really good. However, the residuals are autocorrelated.
Model with lag terms
Now, we’ll try to deal with the residual autocorrelation. We’ll run a series of models adding 1 and 2 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_gam_sflow_lag1 <- gam(
Chla_log ~ s(Flow, by = StationCode) +
s(Flow, by = Year_fct) +
Year_fct * StationCode +
s(Week, bs = "cc", k = 5) +
lag1,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_sflow_lag1))
Box.test(residuals(m_gam_sflow_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_sflow_lag1)
## X-squared = 27.978, df = 20, p-value = 0.1099Lag 2
m_gam_sflow_lag2 <- gam(
Chla_log ~ s(Flow, by = StationCode) +
s(Flow, by = Year_fct) +
Year_fct * StationCode +
s(Week, bs = "cc", k = 5) +
lag1 +
lag2,
data = df_chla_c2_lag,
method = "REML",
knots = list(Week = c(0, 52))
)
acf(residuals(m_gam_sflow_lag2))
Box.test(residuals(m_gam_sflow_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_sflow_lag2)
## X-squared = 25.085, df = 20, p-value = 0.1982The model with 2 lag terms seems to be okay in terms of serial autocorrelation. Let’s use AIC to see how they compare.
Compare Models
AIC(m_gam_sflow, m_gam_sflow_lag1, m_gam_sflow_lag2)## df AIC
## m_gam_sflow 37.77154 235.9396
## m_gam_sflow_lag1 36.39619 141.1780
## m_gam_sflow_lag2 43.29041 118.4925It looks like the lag2 model has the best fit according to the AIC values. Let’s take a closer look at that one.
Lag 2 model summary
summary(m_gam_sflow_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ s(Flow, by = StationCode) + s(Flow, by = Year_fct) +
## Year_fct * StationCode + s(Week, bs = "cc", k = 5) + lag1 +
## lag2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.93613 0.49719 11.939 < 2e-16 ***
## Year_fct2016 -0.46750 0.15129 -3.090 0.002233 **
## Year_fct2017 -0.29715 0.13290 -2.236 0.026268 *
## Year_fct2018 -0.25798 0.10614 -2.431 0.015794 *
## Year_fct2019 -0.06591 0.10535 -0.626 0.532112
## StationCodeSTTD -0.50788 0.11588 -4.383 1.74e-05 ***
## StationCodeLIB -1.19970 0.14229 -8.431 3.00e-15 ***
## StationCodeRVB -1.25900 0.14378 -8.756 3.43e-16 ***
## lag1 0.48127 0.05898 8.160 1.79e-14 ***
## lag2 -0.09558 0.05472 -1.747 0.081933 .
## Year_fct2016:StationCodeSTTD 0.57426 0.19182 2.994 0.003039 **
## Year_fct2017:StationCodeSTTD 0.28317 0.17741 1.596 0.111744
## Year_fct2018:StationCodeSTTD -0.09540 0.15250 -0.626 0.532168
## Year_fct2019:StationCodeSTTD -0.56643 0.16020 -3.536 0.000486 ***
## Year_fct2016:StationCodeLIB 0.84535 0.19645 4.303 2.44e-05 ***
## Year_fct2017:StationCodeLIB 0.41287 0.16229 2.544 0.011579 *
## Year_fct2018:StationCodeLIB -1.01914 0.17242 -5.911 1.14e-08 ***
## Year_fct2019:StationCodeLIB -0.35713 0.15839 -2.255 0.025040 *
## Year_fct2016:StationCodeRVB 0.33325 0.18543 1.797 0.073555 .
## Year_fct2017:StationCodeRVB 0.15440 0.15792 0.978 0.329208
## Year_fct2018:StationCodeRVB 0.10157 0.14341 0.708 0.479493
## Year_fct2019:StationCodeRVB -0.37451 0.14662 -2.554 0.011247 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Flow):StationCodeRD22 1.000e+00 1.000e+00 0.348 0.5559
## s(Flow):StationCodeSTTD 2.357e+00 2.929e+00 2.533 0.0587 .
## s(Flow):StationCodeLIB 2.456e+00 3.049e+00 2.597 0.0525 .
## s(Flow):StationCodeRVB 1.000e+00 1.000e+00 0.065 0.7989
## s(Flow):Year_fct2015 1.000e+00 1.000e+00 0.052 0.8196
## s(Flow):Year_fct2016 5.612e+00 6.468e+00 2.625 0.1007
## s(Flow):Year_fct2017 3.691e-06 7.327e-06 0.010 0.9998
## s(Flow):Year_fct2018 1.000e+00 1.000e+00 0.008 0.9296
## s(Flow):Year_fct2019 1.000e+00 1.000e+00 0.036 0.8489
## s(Week) 2.545e+00 3.000e+00 16.088 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 105/106
## R-sq.(adj) = 0.932 Deviance explained = 94.1%
## -REML = 91.294 Scale est. = 0.076243 n = 284appraise(m_gam_sflow_lag2)
shapiro.test(residuals(m_gam_sflow_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_sflow_lag2)
## W = 0.9633, p-value = 1.295e-06k.check(m_gam_sflow_lag2)## k' edf k-index p-value
## s(Flow):StationCodeRD22 9 1.000010e+00 1.123814 0.9700
## s(Flow):StationCodeSTTD 9 2.356924e+00 1.123814 0.9925
## s(Flow):StationCodeLIB 9 2.456281e+00 1.123814 0.9675
## s(Flow):StationCodeRVB 9 1.000030e+00 1.123814 0.9900
## s(Flow):Year_fct2015 9 1.000031e+00 1.123814 0.9750
## s(Flow):Year_fct2016 9 5.612204e+00 1.123814 0.9825
## s(Flow):Year_fct2017 9 3.691465e-06 1.123814 0.9900
## s(Flow):Year_fct2018 9 1.000020e+00 1.123814 0.9775
## s(Flow):Year_fct2019 9 1.000049e+00 1.123814 0.9750
## s(Week) 3 2.544953e+00 1.009239 0.5375concurvity(m_gam_sflow_lag2, full = FALSE)$worst## para s(Flow):StationCodeRD22 s(Flow):StationCodeSTTD s(Flow):StationCodeLIB
## para 1.000000e+00 2.570357e-01 6.467291e-02 6.873111e-02
## s(Flow):StationCodeRD22 2.570357e-01 1.000000e+00 6.832305e-27 5.903932e-27
## s(Flow):StationCodeSTTD 6.467291e-02 6.800983e-27 1.000000e+00 3.215191e-28
## s(Flow):StationCodeLIB 6.873111e-02 6.780406e-27 3.701507e-28 1.000000e+00
## s(Flow):StationCodeRVB 8.179021e-02 5.058301e-27 5.964115e-28 7.386603e-28
## s(Flow):Year_fct2015 2.288070e-01 4.243169e-01 2.480268e-01 2.881836e-01
## s(Flow):Year_fct2016 1.865535e-01 6.038788e-01 1.530147e-01 2.425054e-01
## s(Flow):Year_fct2017 1.866197e-01 3.603111e-01 6.561655e-02 1.890768e-01
## s(Flow):Year_fct2018 1.930861e-01 3.083022e-01 1.807200e-01 1.938443e-01
## s(Flow):Year_fct2019 2.019358e-01 8.269739e-01 2.900501e-01 2.806141e-01
## s(Week) 1.672167e-29 1.645987e-01 7.498632e-02 9.984918e-02
## s(Flow):StationCodeRVB s(Flow):Year_fct2015 s(Flow):Year_fct2016 s(Flow):Year_fct2017
## para 8.179021e-02 2.288070e-01 1.865535e-01 1.866197e-01
## s(Flow):StationCodeRD22 5.535822e-27 4.243169e-01 6.038788e-01 3.603111e-01
## s(Flow):StationCodeSTTD 3.328124e-28 2.480268e-01 1.530147e-01 6.561655e-02
## s(Flow):StationCodeLIB 3.558085e-28 2.881836e-01 2.425054e-01 1.890768e-01
## s(Flow):StationCodeRVB 1.000000e+00 2.467403e-01 1.556821e-01 1.858040e-01
## s(Flow):Year_fct2015 2.467403e-01 1.000000e+00 4.157153e-26 1.978233e-12
## s(Flow):Year_fct2016 1.556821e-01 2.858411e-26 1.000000e+00 7.668165e-14
## s(Flow):Year_fct2017 1.858040e-01 1.650521e-12 3.566287e-15 1.000000e+00
## s(Flow):Year_fct2018 1.920115e-01 1.077666e-25 8.938628e-28 3.281103e-13
## s(Flow):Year_fct2019 2.889471e-01 5.409089e-27 1.689988e-27 4.718115e-13
## s(Week) 7.686187e-02 1.550599e-01 3.097224e-01 1.944685e-01
## s(Flow):Year_fct2018 s(Flow):Year_fct2019 s(Week)
## para 1.930861e-01 2.019358e-01 9.486679e-30
## s(Flow):StationCodeRD22 3.083022e-01 8.269739e-01 1.645987e-01
## s(Flow):StationCodeSTTD 1.807200e-01 2.900501e-01 7.498632e-02
## s(Flow):StationCodeLIB 1.938443e-01 2.806141e-01 9.984918e-02
## s(Flow):StationCodeRVB 1.920115e-01 2.889471e-01 7.686187e-02
## s(Flow):Year_fct2015 1.118029e-25 4.642002e-27 1.550599e-01
## s(Flow):Year_fct2016 7.576510e-28 1.173662e-27 3.097224e-01
## s(Flow):Year_fct2017 2.499652e-13 2.531812e-13 1.944685e-01
## s(Flow):Year_fct2018 1.000000e+00 2.888272e-27 1.118373e-01
## s(Flow):Year_fct2019 2.569263e-27 1.000000e+00 1.047812e-01
## s(Week) 1.118373e-01 1.047812e-01 1.000000e+00draw(m_gam_sflow_lag2, select = 1:4, residuals = TRUE, rug = FALSE)
draw(m_gam_sflow_lag2, select = 5:9, residuals = TRUE, rug = FALSE)
draw(m_gam_sflow_lag2, select = 10, residuals = TRUE, rug = FALSE)
anova(m_gam_sflow_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ s(Flow, by = StationCode) + s(Flow, by = Year_fct) +
## Year_fct * StationCode + s(Week, bs = "cc", k = 5) + lag1 +
## lag2
##
## Parametric Terms:
## df F p-value
## Year_fct 4 3.673 0.00633
## StationCode 3 29.298 3.26e-16
## lag1 1 66.583 1.79e-14
## lag2 1 3.051 0.08193
## Year_fct:StationCode 12 7.868 2.23e-12
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Flow):StationCodeRD22 1.000e+00 1.000e+00 0.348 0.5559
## s(Flow):StationCodeSTTD 2.357e+00 2.929e+00 2.533 0.0587
## s(Flow):StationCodeLIB 2.456e+00 3.049e+00 2.597 0.0525
## s(Flow):StationCodeRVB 1.000e+00 1.000e+00 0.065 0.7989
## s(Flow):Year_fct2015 1.000e+00 1.000e+00 0.052 0.8196
## s(Flow):Year_fct2016 5.612e+00 6.468e+00 2.625 0.1007
## s(Flow):Year_fct2017 3.691e-06 7.327e-06 0.010 0.9998
## s(Flow):Year_fct2018 1.000e+00 1.000e+00 0.008 0.9296
## s(Flow):Year_fct2019 1.000e+00 1.000e+00 0.036 0.8489
## s(Week) 2.545e+00 3.000e+00 16.088 <2e-16The model diagnostics look a little worse than with the initial model
but they still look pretty good. We’ll use m_gam_sflow_lag2
in the model selection process.
rm(m_gam_sflow, m_gam_sflow_lag1)Model selection
Now we’ll compare the seven candidate models with AIC to select the one with the best fit. As a summary, here are the 7 models we are comparing:
- Model 1 -
m_gam_flow3_lag2- GAM 3-way interactions with s(Week)
Formula: Chla_log ~ Year_fct * Flow * StationCode + s(Week, bs = “cc”, k = 5) + lag1 + lag2 - Model 2 -
m_gam_flow2_lag2- GAM 2-way interactions with s(Week)
Formula: Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(Week, bs = “cc”, k = 5) + lag1 + lag2 - Model 3 -
m_gam_cat2_lag2- GAM 2-way interaction between Station and Year with s(Week) but without Flow
Formula: Chla_log ~ Year_fct * StationCode + s(Week, bs = “cc”, k = 5) + lag1 + lag2 - Model 4 -
m_lm_flow3_lag2- LM 3-way interactions without seasonal term
Formula: Chla_log ~ Year_fct * Flow * StationCode + lag1 + lag2 - Model 5 -
m_lm_flow2_lag2- LM 2-way interactions without seasonal term
Formula: Chla_log ~ (Year_fct + Flow + StationCode)^2 + lag1 + lag2 - Model 6 -
m_lm_cat2_lag2- LM 2-way interaction between Station and Year but without Flow and seasonal term
Formula: Chla_log ~ Year_fct * StationCode + lag1 + lag2 - Model 7 -
m_gam_sflow_lag2- GAM using smooths for Flow with s(Week)
Formula: Chla_log ~ s(Flow, by = StationCode) + s(Flow, by = Year_fct) + Year_fct * StationCode + s(Week, bs = “cc”, k = 5) + lag1 + lag2
# AIC values
df_m_aic <-
AIC(
m_gam_flow3_lag2,
m_gam_flow2_lag2,
m_gam_cat2_lag2,
m_lm_flow3_lag2,
m_lm_flow2_lag2,
m_lm_cat2_lag2,
m_gam_sflow_lag2
) %>%
as_tibble(rownames = "Model") %>%
mutate(Model_Number = 1:7, .before = Model)
# BIC values
df_m_bic <-
BIC(
m_gam_flow3_lag2,
m_gam_flow2_lag2,
m_gam_cat2_lag2,
m_lm_flow3_lag2,
m_lm_flow2_lag2,
m_lm_cat2_lag2,
m_gam_sflow_lag2
) %>%
as_tibble(rownames = "Model")
# Combine AIC and BIC and calculate differences from lowest value
df_m_aic_bic <-
left_join(df_m_aic, df_m_bic, by = join_by(Model, df)) %>%
mutate(across(c(AIC, BIC), ~ .x - min(.x), .names = "{.col}_delta")) %>%
select(starts_with("Model"), df, starts_with("AIC"), starts_with("BIC"))
# Sort by AIC
df_m_aic_bic %>% arrange(AIC)## # A tibble: 7 × 7
## Model_Number Model df AIC AIC_delta BIC BIC_delta
## <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 7 m_gam_sflow_lag2 43.3 118. 0 276. 15.3
## 2 1 m_gam_flow3_lag2 45.8 125. 6.77 292. 31.2
## 3 2 m_gam_flow2_lag2 33.7 138. 19.5 261. 0
## 4 4 m_lm_flow3_lag2 43 172. 53.7 329. 68.0
## 5 5 m_lm_flow2_lag2 31 176. 57.9 290. 28.4
## 6 3 m_gam_cat2_lag2 25.5 185. 66.1 278. 16.5
## 7 6 m_lm_cat2_lag2 23 201. 82.9 285. 24.2Export the AIC/BIC table.
df_m_aic_bic %>%
arrange(AIC) %>%
mutate(across(
where(is.numeric),
~ paste0(formatC(.x, digits = 1, format = "f"), "##")
)) %>%
write_csv(here(
"manuscript_synthesis/results/tables/chl_aic_weekly_models.csv"
))According to AIC, Model 7 (GAM using smooths for Flow with s(Week)) was the model with the best fit. BIC preferred Model 2 (GAM 2-way interactions with s(Week)). BIC prefers simpler models, so I’m inclined to go with Model 2, especially since its difficult to interpret the smooths for Flow in Model 7. Before we proceed with Model 2, let’s revisit the model diagnostics and take a closer look at how the back-transformed fitted values from the model match the observed values.
Model 2 diagnostics
appraise(m_gam_flow2_lag2)
shapiro.test(residuals(m_gam_flow2_lag2))##
## Shapiro-Wilk normality test
##
## data: residuals(m_gam_flow2_lag2)
## W = 0.96549, p-value = 2.591e-06k.check(m_gam_flow2_lag2)## k' edf k-index p-value
## s(Week) 3 2.454221 0.978034 0.36draw(m_gam_flow2_lag2, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_flow2_lag2, pages = 1, all.terms = TRUE)
anova(m_gam_flow2_lag2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(Week, bs = "cc",
## k = 5) + lag1 + lag2
##
## Parametric Terms:
## df F p-value
## Year_fct 4 3.998 0.00366
## Flow 1 8.307 0.00429
## StationCode 3 25.413 2.15e-14
## lag1 1 75.957 3.98e-16
## lag2 1 4.335 0.03834
## Year_fct:Flow 4 1.764 0.13668
## Year_fct:StationCode 12 7.377 1.34e-11
## Flow:StationCode 3 18.037 1.24e-10
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 2.454 3.000 13.2 <2e-16Model 2 observed vs fitted values
df_chla_c2_lag2 <- df_chla_c2_lag %>% drop_na(lag1, lag2)
df_m_gam_flow2_lag2_fit <- df_chla_c2_lag2 %>%
fitted_values(m_gam_flow2_lag2, data = .) %>%
mutate(fitted_bt = exp(.fitted) / 1000)
plt_m_gam_flow2_lag2_fit <- df_m_gam_flow2_lag2_fit %>%
ggplot(aes(x = fitted_bt, y = Chla)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") +
theme_bw() +
labs(
x = "Back-transformed Fitted Values",
y = "Observed Values"
)
plt_m_gam_flow2_lag2_fit
Let’s group by station.
plt_m_gam_flow2_lag2_fit + facet_wrap(vars(StationCode), scales = "free")
Now, group by year.
plt_m_gam_flow2_lag2_fit + facet_wrap(vars(Year_fct), scales = "free")
Everything looks pretty decent. Not perfect, but pretty good given the number of data points. Note that variability does increase as the chlorophyll values increase. Before proceeding with Model 2, let’s look more closely at its results.
Model 2 Results
Effect of Flow by Station
# Calculate min and max flows for each station to narrow down x-axis in the plot
df_chla_flow_sta_summ <- df_chla_c2_lag2 %>%
summarize(
Flow_min = min(Flow),
Flow_max = max(Flow),
.by = c(StationCode)
) %>%
mutate(
Flow_buffer = (Flow_max - Flow_min) * 0.05,
Flow_min = Flow_min - Flow_buffer,
Flow_max = Flow_max + Flow_buffer
)
# Calculate effects of flow on chlorophyll for each station holding the
# non-focal variables constant - marginal effects/adjusted predictions
df_gam_flow2_sta_eff <-
as.data.frame(
predict_response(
m_gam_flow2_lag2,
terms = c("Flow", "StationCode"),
margin = "marginalmeans"
),
terms_to_colnames = TRUE
) %>%
as_tibble() %>%
# Narrow down range of flow values for each station
left_join(df_chla_flow_sta_summ, by = join_by(StationCode)) %>%
filter(Flow >= Flow_min & Flow <= Flow_max) %>%
transmute(
StationCode,
Flow,
# Back calculate model fits and confidence levels
across(c(predicted, conf.low, conf.high), ~ exp(.x) / 1000)
)
# Create effects plot
plt_gam_flow2_sta_eff <- df_gam_flow2_sta_eff %>%
ggplot(aes(x = Flow, y = predicted)) +
geom_point(
data = df_chla_c2_lag2,
aes(y = Chla, color = Year_fct),
alpha = 0.6
) +
geom_line(linewidth = 1) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.25) +
facet_wrap(vars(StationCode), scales = "free") +
theme_bw() +
labs(
x = "Flow (cfs)",
y = expression(
Chlorophyll ~ Fluoresence ~ (mu * g ~ L^{
-1
})
)
) +
scale_x_continuous(breaks = breaks_extended(6)) +
scale_color_viridis_d(name = "Year", option = "C")
plt_gam_flow2_sta_eff
Effect of Flow by Year
# Calculate min and max flows for each station to narrow down x-axis in the plot
df_chla_flow_yr_summ <- df_chla_c2_lag2 %>%
summarize(
Flow_min = min(Flow),
Flow_max = max(Flow),
.by = c(Year_fct)
) %>%
mutate(
Flow_buffer = (Flow_max - Flow_min) * 0.05,
Flow_min = Flow_min - Flow_buffer,
Flow_max = Flow_max + Flow_buffer
)
# Calculate effects of flow on chlorophyll for each year holding the
# non-focal variables constant - marginal effects/adjusted predictions
df_gam_flow2_yr_eff <-
as.data.frame(
predict_response(
m_gam_flow2_lag2,
terms = c("Flow", "Year_fct"),
margin = "marginalmeans"
),
terms_to_colnames = TRUE
) %>%
as_tibble() %>%
# Narrow down range of flow values for each station
left_join(df_chla_flow_yr_summ, by = join_by(Year_fct)) %>%
filter(Flow >= Flow_min & Flow <= Flow_max) %>%
transmute(
Year_fct,
Flow,
# Back calculate model fits and confidence levels
across(c(predicted, conf.low, conf.high), ~ exp(.x) / 1000)
)
# Create effects plot
plt_gam_flow2_yr_eff <- df_gam_flow2_yr_eff %>%
ggplot(aes(x = Flow, y = predicted)) +
geom_point(
data = df_chla_c2_lag2,
aes(y = Chla, color = StationCode),
alpha = 0.6
) +
geom_line(linewidth = 1) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.25) +
facet_wrap(vars(Year_fct), scales = "free") +
theme_bw() +
labs(
x = "Flow (cfs)",
y = expression(
Chlorophyll ~ Fluoresence ~ (mu * g ~ L^{
-1
})
)
) +
scale_x_continuous(breaks = breaks_extended(6)) +
scale_color_viridis_d(name = "Station", option = "C") +
theme(
legend.margin = margin(0, 0, 0, 0),
legend.position = "inside",
legend.position.inside = c(0.8, 0.3)
)
plt_gam_flow2_yr_eff
These results look reasonable, but not interesting.
Station by Year Contrasts
# Estimated marginal means for station by year
em_gam2_sta_yr <- emmeans(m_gam_flow2_lag2, ~ StationCode | Year_fct)
# Tukey post-hoc contrasts
pairs(em_gam2_sta_yr)## Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## RD22 - STTD 0.43227 0.1190 252 3.641 0.0019
## RD22 - LIB 1.09948 0.1440 252 7.617 <0.0001
## RD22 - RVB 1.16562 0.1480 252 7.887 <0.0001
## STTD - LIB 0.66720 0.1170 252 5.717 <0.0001
## STTD - RVB 0.73334 0.1190 252 6.150 <0.0001
## LIB - RVB 0.06614 0.0979 252 0.676 0.9062
##
## Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## RD22 - STTD -0.08302 0.1210 252 -0.689 0.9014
## RD22 - LIB 0.28919 0.1190 252 2.427 0.0746
## RD22 - RVB 0.84908 0.1410 252 6.035 <0.0001
## STTD - LIB 0.37221 0.1160 252 3.204 0.0083
## STTD - RVB 0.93209 0.1370 252 6.791 <0.0001
## LIB - RVB 0.55988 0.1190 252 4.691 <0.0001
##
## Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## RD22 - STTD 0.21045 0.1460 252 1.442 0.4742
## RD22 - LIB 0.75128 0.1460 252 5.161 <0.0001
## RD22 - RVB 1.04564 0.1560 252 6.693 <0.0001
## STTD - LIB 0.54083 0.1410 252 3.832 0.0009
## STTD - RVB 0.83519 0.1500 252 5.581 <0.0001
## LIB - RVB 0.29436 0.1110 252 2.644 0.0429
##
## Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## RD22 - STTD 0.57568 0.1250 252 4.606 <0.0001
## RD22 - LIB 2.08925 0.2110 252 9.919 <0.0001
## RD22 - RVB 1.09304 0.1430 252 7.653 <0.0001
## STTD - LIB 1.51358 0.1740 252 8.707 <0.0001
## STTD - RVB 0.51736 0.1200 252 4.320 0.0001
## LIB - RVB -0.99621 0.1400 252 -7.109 <0.0001
##
## Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## RD22 - STTD 1.08467 0.1390 252 7.792 <0.0001
## RD22 - LIB 1.55309 0.1680 252 9.259 <0.0001
## RD22 - RVB 1.55988 0.1670 252 9.332 <0.0001
## STTD - LIB 0.46842 0.1220 252 3.851 0.0009
## STTD - RVB 0.47520 0.1160 252 4.089 0.0003
## LIB - RVB 0.00679 0.1120 252 0.060 0.9999
##
## P value adjustment: tukey method for comparing a family of 4 estimates# Create table of contrasts and convert it to a tibble for plot
df_gam2_sta_yr <- em_gam2_sta_yr %>%
cld(sort = FALSE, Letters = letters) %>%
as_tibble() %>%
mutate(
group = str_remove_all(.group, fixed(" ")),
# back transform log-transformed results
across(c(emmean, lower.CL, upper.CL), ~ exp(.x) / 1000)
) %>%
# Add min and max values of observed data to the Tukey post-hoc results and
# calculate vertical positioning of letters
left_join(
df_chla_c2_lag2 %>%
summarize(
max_val = max(Chla),
min_val = min(Chla),
.by = Year_fct
),
by = join_by(Year_fct)
) %>%
mutate(max_val = if_else(upper.CL > max_val, upper.CL, max_val)) %>%
group_by(Year_fct) %>%
mutate(max_val = max(max_val)) %>%
ungroup() %>%
mutate(y_pos = max_val + (max_val - min_val) / 10) %>%
select(
StationCode,
Year_fct,
emmean,
lower.CL,
upper.CL,
group,
y_pos
)
# Create boxplot showing Tukey post-hoc results
plt_gam2_sta_yr <- df_gam2_sta_yr %>%
ggplot(
aes(
x = StationCode,
y = emmean,
ymin = lower.CL,
ymax = upper.CL
)
) +
geom_boxplot(
data = df_chla_c2_lag2,
aes(x = StationCode, y = Chla),
inherit.aes = FALSE
) +
geom_crossbar(color = "grey82", fill = "grey", alpha = 0.7, linewidth = 0.1) +
geom_point(color = "red") +
geom_text(aes(y = y_pos, label = group), size = 3.5) +
facet_wrap(vars(Year_fct), scales = "free_y") +
xlab("Station") +
ylab(expression(
Chlorophyll ~ Fluoresence ~ (mu * g ~ L^{
-1
})
)) +
theme_bw()
plt_gam2_sta_yr
Many of the contrasts are significant. The model under predicts RD22, but the uncertainty looks reasonable for all values.
Year by Station Contrasts
# Estimated marginal means for year by station
em_gam2_yr_sta <- emmeans(m_gam_flow2_lag2, ~ Year_fct | StationCode)
# Tukey post-hoc contrasts
pairs(em_gam2_yr_sta)## StationCode = RD22:
## contrast estimate SE df t.ratio p.value
## Year_fct2015 - Year_fct2016 0.41530 0.124 252 3.358 0.0080
## Year_fct2015 - Year_fct2017 0.23419 0.138 252 1.697 0.4377
## Year_fct2015 - Year_fct2018 0.22938 0.112 252 2.057 0.2425
## Year_fct2015 - Year_fct2019 0.02464 0.110 252 0.223 0.9994
## Year_fct2016 - Year_fct2017 -0.18111 0.137 252 -1.322 0.6781
## Year_fct2016 - Year_fct2018 -0.18593 0.115 252 -1.614 0.4897
## Year_fct2016 - Year_fct2019 -0.39067 0.117 252 -3.340 0.0085
## Year_fct2017 - Year_fct2018 -0.00481 0.131 252 -0.037 1.0000
## Year_fct2017 - Year_fct2019 -0.20956 0.133 252 -1.580 0.5115
## Year_fct2018 - Year_fct2019 -0.20474 0.105 252 -1.941 0.2982
##
## StationCode = STTD:
## contrast estimate SE df t.ratio p.value
## Year_fct2015 - Year_fct2016 -0.09999 0.116 252 -0.859 0.9115
## Year_fct2015 - Year_fct2017 0.01237 0.190 252 0.065 1.0000
## Year_fct2015 - Year_fct2018 0.37278 0.119 252 3.140 0.0161
## Year_fct2015 - Year_fct2019 0.67704 0.124 252 5.481 <0.0001
## Year_fct2016 - Year_fct2017 0.11235 0.190 252 0.593 0.9761
## Year_fct2016 - Year_fct2018 0.47277 0.128 252 3.680 0.0026
## Year_fct2016 - Year_fct2019 0.77702 0.136 252 5.714 <0.0001
## Year_fct2017 - Year_fct2018 0.36041 0.195 252 1.844 0.3508
## Year_fct2017 - Year_fct2019 0.66467 0.202 252 3.287 0.0101
## Year_fct2018 - Year_fct2019 0.30426 0.118 252 2.576 0.0779
##
## StationCode = LIB:
## contrast estimate SE df t.ratio p.value
## Year_fct2015 - Year_fct2016 -0.39498 0.114 252 -3.461 0.0057
## Year_fct2015 - Year_fct2017 -0.11400 0.171 252 -0.666 0.9634
## Year_fct2015 - Year_fct2018 1.21916 0.152 252 8.004 <0.0001
## Year_fct2015 - Year_fct2019 0.47825 0.117 252 4.073 0.0006
## Year_fct2016 - Year_fct2017 0.28098 0.178 252 1.583 0.5098
## Year_fct2016 - Year_fct2018 1.61414 0.188 252 8.581 <0.0001
## Year_fct2016 - Year_fct2019 0.87323 0.143 252 6.106 <0.0001
## Year_fct2017 - Year_fct2018 1.33316 0.213 252 6.265 <0.0001
## Year_fct2017 - Year_fct2019 0.59225 0.187 252 3.171 0.0146
## Year_fct2018 - Year_fct2019 -0.74091 0.137 252 -5.397 <0.0001
##
## StationCode = RVB:
## contrast estimate SE df t.ratio p.value
## Year_fct2015 - Year_fct2016 0.09876 0.109 252 0.908 0.8936
## Year_fct2015 - Year_fct2017 0.11422 0.171 252 0.668 0.9630
## Year_fct2015 - Year_fct2018 0.15680 0.102 252 1.537 0.5396
## Year_fct2015 - Year_fct2019 0.41890 0.109 252 3.838 0.0015
## Year_fct2016 - Year_fct2017 0.01545 0.175 252 0.088 1.0000
## Year_fct2016 - Year_fct2018 0.05804 0.112 252 0.517 0.9856
## Year_fct2016 - Year_fct2019 0.32013 0.117 252 2.731 0.0523
## Year_fct2017 - Year_fct2018 0.04259 0.173 252 0.246 0.9992
## Year_fct2017 - Year_fct2019 0.30468 0.177 252 1.726 0.4201
## Year_fct2018 - Year_fct2019 0.26209 0.106 252 2.467 0.1016
##
## P value adjustment: tukey method for comparing a family of 5 estimates# Create table of contrasts and convert it to a tibble for plot
df_gam2_yr_sta <- em_gam2_yr_sta %>%
cld(sort = FALSE, Letters = letters) %>%
as_tibble() %>%
mutate(
group = str_remove_all(.group, fixed(" ")),
# back transform log-transformed results
across(c(emmean, lower.CL, upper.CL), ~ exp(.x) / 1000)
) %>%
# Add min and max values of observed data to the Tukey post-hoc results and
# calculate vertical positioning of letters
left_join(
df_chla_c2_lag2 %>%
summarize(
max_val = max(Chla),
min_val = min(Chla),
.by = StationCode
),
by = join_by(StationCode)
) %>%
mutate(max_val = if_else(upper.CL > max_val, upper.CL, max_val)) %>%
group_by(StationCode) %>%
mutate(max_val = max(max_val)) %>%
ungroup() %>%
mutate(y_pos = max_val + (max_val - min_val) / 10) %>%
select(
StationCode,
Year_fct,
emmean,
lower.CL,
upper.CL,
group,
y_pos
)
# Create boxplot showing Tukey post-hoc results
plt_gam2_yr_sta <- df_gam2_yr_sta %>%
ggplot(
aes(
x = Year_fct,
y = emmean,
ymin = lower.CL,
ymax = upper.CL
)
) +
geom_boxplot(
data = df_chla_c2_lag2,
aes(x = Year_fct, y = Chla),
inherit.aes = FALSE
) +
geom_crossbar(color = "grey82", fill = "grey", alpha = 0.7, linewidth = 0.1) +
geom_point(color = "red") +
geom_text(aes(y = y_pos, label = group), size = 3.5) +
facet_wrap(vars(StationCode), scales = "free_y") +
xlab("Year") +
ylab(expression(
Chlorophyll ~ Fluoresence ~ (mu * g ~ L^{
-1
})
)) +
theme_bw()
plt_gam2_yr_sta
Again, the model under predicts RD22, but the uncertainty looks reasonable for all values.
Conclusions
After looking at the results more closely from Model 2, I think Model 2 (GAM 2-way interactions with s(Week)) looks pretty good.