NDFS Synthesis Manuscript: Chlorophyll analysis
Models using categorical predictors - Revised Analysis
Dave Bosworth
February 29, 2024
Purpose
Explore and analyze the continuous chlorophyll data to be included in the NDFS synthesis manuscript. We will attempt to fit multiple models to the data set and use a model selection process to determine the best one. At a minimum the model will contain the three categorical variables (Year, Station, and Flow Pulse Period) and optionally a term to account for seasonality (either daily average water temperature or a GAM smooth for DOY). These models will only include representative stations for 4 habitat types - upstream (RD22), lower Yolo Bypass (STTD), Cache Slough complex (LIB), and downstream (RVB).
Global code and functions
# Load packages
library(tidyverse)
library(scales)
library(knitr)
library(mgcv)
library(car)
library(gratia)
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.2.3 (2023-03-15 ucrt)
## os Windows 10 x64 (build 19045)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate English_United States.utf8
## ctype English_United States.utf8
## tz America/Los_Angeles
## date 2024-02-29
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
## bslib 0.4.2 2022-12-16 [1] CRAN (R 4.2.2)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.2.3)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.2)
## car * 3.1-2 2023-03-30 [1] CRAN (R 4.2.3)
## carData * 3.0-5 2022-01-06 [1] CRAN (R 4.2.1)
## cli 3.6.2 2023-12-11 [1] CRAN (R 4.2.3)
## coda 0.19-4 2020-09-30 [1] CRAN (R 4.2.1)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.2.3)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.2)
## conflicted * 1.2.0 2023-02-01 [1] CRAN (R 4.2.2)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.1)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.2.1)
## digest 0.6.33 2023-07-07 [1] CRAN (R 4.2.3)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.2.3)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.1)
## emmeans * 1.8.5 2023-03-08 [1] CRAN (R 4.2.2)
## estimability 1.4.1 2022-08-05 [1] CRAN (R 4.2.1)
## evaluate 0.21 2023-05-05 [1] CRAN (R 4.2.3)
## fansi 1.0.6 2023-12-08 [1] CRAN (R 4.2.3)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.2)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.2.2)
## fs * 1.6.3 2023-07-20 [1] CRAN (R 4.2.3)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.1)
## ggplot2 * 3.4.3 2023-08-14 [1] CRAN (R 4.2.3)
## glue 1.7.0 2024-01-09 [1] CRAN (R 4.2.3)
## gratia * 0.8.1 2023-02-02 [1] CRAN (R 4.2.2)
## gtable 0.3.4 2023-08-21 [1] CRAN (R 4.2.3)
## here * 1.0.1 2020-12-13 [1] CRAN (R 4.2.1)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.2.3)
## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.2.3)
## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.2.3)
## httpuv 1.6.9 2023-02-14 [1] CRAN (R 4.2.2)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.1)
## jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.2.3)
## knitr * 1.42 2023-01-25 [1] CRAN (R 4.2.2)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.2.1)
## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.2.3)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.2.3)
## lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.2.3)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.1)
## MASS * 7.3-58.2 2023-01-23 [2] CRAN (R 4.2.3)
## Matrix 1.5-4 2023-04-04 [1] CRAN (R 4.2.3)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.1)
## mgcv * 1.8-42 2023-03-02 [1] CRAN (R 4.2.2)
## mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.1)
## multcomp * 1.4-23 2023-03-09 [1] CRAN (R 4.2.2)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.1)
## mvnfast 0.2.8 2023-02-23 [1] CRAN (R 4.2.2)
## mvtnorm * 1.1-3 2021-10-08 [1] CRAN (R 4.2.0)
## nlme * 3.1-162 2023-01-31 [2] CRAN (R 4.2.3)
## patchwork * 1.1.2 2022-08-19 [1] CRAN (R 4.2.1)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3)
## pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.2.3)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.1)
## pkgload 1.3.2.1 2023-07-08 [1] CRAN (R 4.2.3)
## prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.2.3)
## processx 3.8.2 2023-06-30 [1] CRAN (R 4.2.3)
## profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.1)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.1)
## ps 1.7.5 2023-04-18 [1] CRAN (R 4.2.3)
## purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.2.3)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.1)
## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.2.3)
## readr * 2.1.5 2024-01-10 [1] CRAN (R 4.2.3)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.1)
## rlang * 1.1.3 2024-01-10 [1] CRAN (R 4.2.3)
## rmarkdown 2.21 2023-03-26 [1] CRAN (R 4.2.3)
## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.2.1)
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.1)
## sandwich 3.0-2 2022-06-15 [1] CRAN (R 4.2.1)
## sass 0.4.6 2023-05-03 [1] CRAN (R 4.2.3)
## scales * 1.2.1 2022-08-20 [1] CRAN (R 4.2.1)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.1)
## shiny 1.7.4 2022-12-15 [1] CRAN (R 4.2.2)
## stringi 1.8.3 2023-12-11 [1] CRAN (R 4.2.3)
## stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.2.3)
## survival * 3.5-5 2023-03-12 [1] CRAN (R 4.2.3)
## TH.data * 1.1-2 2023-04-17 [1] CRAN (R 4.2.3)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.2.3)
## tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.2.3)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.1)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.2.2)
## timechange 0.3.0 2024-01-18 [1] CRAN (R 4.2.3)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.2.3)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.1)
## usethis 2.1.6 2022-05-25 [1] CRAN (R 4.2.1)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.2)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.2.3)
## withr 3.0.0 2024-01-16 [1] CRAN (R 4.2.3)
## xfun 0.39 2023-04-20 [1] CRAN (R 4.2.3)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.1)
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.2)
## zoo 1.8-12 2023-04-13 [1] CRAN (R 4.2.3)
##
## [1] C:/R/win-library/4.2
## [2] C:/Program Files/R/R-4.2.3/library
##
## ──────────────────────────────────────────────────────────────────────────────Create functions for document:
# Create summary figure for model results and observed data
plot_model_summary <- function(df_data, em_tuk_obj) {
# Calculate min and max values of observed data for each Station - Year group to
# determine vertical positioning of letters for figure
df_data_summ <- df_data %>%
summarize(
max_val = max(Chla),
min_val = min(Chla),
.by = c(StationCode, Year)
)
# Add significance grouping letters to the Tukey post-hoc results
df_tuk <- em_tuk_obj %>%
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),
Year = as.numeric(as.character(Year_fct))
) %>%
# Add min and max values of observed data to the Tukey post-hoc results and
# calculate vertical positioning of letters
left_join(df_data_summ, by = join_by(StationCode, Year)) %>%
mutate(max_val = if_else(upper.CL > max_val, upper.CL, max_val)) %>%
group_by(StationCode, Year) %>%
mutate(max_val = max(max_val)) %>%
ungroup() %>%
mutate(y_pos = max_val + (max_val - min_val) / 10) %>%
select(
StationCode,
Year,
FlowActionPeriod,
emmean,
lower.CL,
upper.CL,
group,
y_pos
)
# Create summary figure
df_tuk %>%
ggplot(
aes(
x = FlowActionPeriod,
y = emmean,
ymin = lower.CL,
ymax = upper.CL
)
) +
geom_boxplot(
data = df_data,
aes(x = FlowActionPeriod, 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) +
facet_wrap(
vars(StationCode, Year),
scales = "free_y",
nrow = 4,
labeller = labeller(.multi_line = FALSE)
) +
xlab("Flow Pulse Period") +
ylab(expression(Chlorophyll~Fluoresence~(mu*g~L^{-1}))) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
}Import Data
# Define file path for processed data
fp_data <- here("manuscript_synthesis/data/processed")
# Import daily average water quality data
df_wq <- readRDS(file.path(fp_data, "wq_daily_avg_2013-2019.rds"))Prepare Data
# Create a vector for the factor order of StationCode
sta_order <- c("RD22", "STTD", "LIB", "RVB")
# Prepare chlorophyll and water temperature data for exploration and analysis
df_chla_wt_c <- df_wq %>%
select(StationCode, Date, Chla, WaterTemp) %>%
drop_na(Chla) %>%
# Filter to only include representative stations for 3 habitat types - RD22, STTD, LIB
filter(StationCode %in% sta_order) %>%
mutate(
# Scale and log transform chlorophyll values
Chla_log = log(Chla * 1000),
# Add Year and DOY variables
Year = year(Date),
DOY = yday(Date)
) %>%
# Add Flow Action Periods
ndfa_action_periods() %>%
mutate(
# Apply factor orders to FlowActionPeriod and StationCode
FlowActionPeriod = factor(FlowActionPeriod, levels = c("Before", "During", "After")),
StationCode = factor(StationCode, levels = sta_order),
# Add a column for Year as a factor for the model
Year_fct = factor(Year)
) %>%
arrange(StationCode, Date)
# Create another dataframe that has up to 3 lag variables for chlorophyll to be
# used in the linear models
df_chla_wt_lag <- df_chla_wt_c %>%
# Fill in missing days for each StationCode-Year combination
group_by(StationCode, Year) %>%
complete(Date = seq.Date(min(Date), max(Date), by = "1 day")) %>%
# Create lag variables of scaled log transformed chlorophyll values
mutate(
lag1 = lag(Chla_log),
lag2 = lag(Chla_log, 2),
lag3 = lag(Chla_log, 3)
) %>%
ungroup()Explore sample counts by Station
df_chla_wt_c %>%
summarize(
min_date = min(Date),
max_date = max(Date),
num_samples = n(),
.by = c(StationCode, Year, FlowActionPeriod)
) %>%
complete(StationCode, Year, FlowActionPeriod) %>%
arrange(StationCode, Year, FlowActionPeriod) %>%
kable()| StationCode | Year | FlowActionPeriod | min_date | max_date | num_samples |
|---|---|---|---|---|---|
| RD22 | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| RD22 | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| RD22 | 2013 | After | 2013-10-03 | 2013-11-04 | 33 |
| RD22 | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| RD22 | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| RD22 | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| RD22 | 2015 | Before | 2015-07-24 | 2015-08-20 | 28 |
| RD22 | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RD22 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RD22 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| RD22 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RD22 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RD22 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RD22 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RD22 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RD22 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RD22 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RD22 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RD22 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RD22 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RD22 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| STTD | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| STTD | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| STTD | 2013 | After | 2013-10-03 | 2013-11-04 | 33 |
| STTD | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| STTD | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| STTD | 2014 | After | 2014-09-24 | 2014-11-08 | 29 |
| STTD | 2015 | Before | 2015-07-27 | 2015-08-20 | 25 |
| STTD | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| STTD | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| STTD | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| STTD | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| STTD | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| STTD | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| STTD | 2017 | During | 2017-08-29 | 2017-09-01 | 4 |
| STTD | 2017 | After | 2017-09-20 | 2017-09-25 | 6 |
| STTD | 2018 | Before | 2018-07-20 | 2018-08-27 | 39 |
| STTD | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| STTD | 2018 | After | 2018-09-27 | 2018-10-15 | 19 |
| STTD | 2019 | Before | 2019-07-26 | 2019-08-25 | 31 |
| STTD | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| STTD | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIB | 2013 | Before | 2013-07-16 | 2013-08-21 | 37 |
| LIB | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| LIB | 2013 | After | 2013-10-03 | 2013-11-17 | 46 |
| LIB | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| LIB | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| LIB | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| LIB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIB | 2015 | During | 2015-08-21 | 2015-10-01 | 38 |
| LIB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIB | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIB | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| LIB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIB | 2018 | Before | 2018-08-14 | 2018-08-27 | 14 |
| LIB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIB | 2019 | Before | 2019-07-11 | 2019-07-30 | 20 |
| LIB | 2019 | During | 2019-09-05 | 2019-09-21 | 17 |
| LIB | 2019 | After | 2019-09-22 | 2019-11-06 | 39 |
| RVB | 2013 | Before | 2013-07-07 | 2013-08-21 | 46 |
| RVB | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| RVB | 2013 | After | 2013-10-03 | 2013-11-17 | 46 |
| RVB | 2014 | Before | 2014-07-25 | 2014-09-07 | 45 |
| RVB | 2014 | During | 2014-09-10 | 2014-09-23 | 14 |
| RVB | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| RVB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| RVB | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RVB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| RVB | 2016 | Before | 2016-05-29 | 2016-07-13 | 39 |
| RVB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RVB | 2016 | After | 2016-08-02 | 2016-09-16 | 37 |
| RVB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RVB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RVB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RVB | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RVB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RVB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RVB | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RVB | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RVB | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
Except for a few Station-Year-FlowPeriod combinations with low sample counts, it appears we have adequate data.
Plots
df_chla_wt_c %>%
ggplot(aes(x = FlowActionPeriod, y = Chla, fill = FlowActionPeriod)) +
geom_boxplot() +
facet_grid(cols = vars(Year), rows = vars(StationCode)) +
scale_y_log10(labels = trans_format("log10", math_format(10^.x))) +
annotation_logticks(sides = "l") +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "top"
)
Models
Model 1: GAM 3-way interactions with s(DOY)
Initial Model
We’ll try running a GAM using a three-way interaction between Year, Flow Action Period, and Station, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_cat3 <- gam(
Chla_log ~ Year_fct * FlowActionPeriod * StationCode + s(DOY, k = 5),
data = df_chla_wt_c,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_cat3)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * FlowActionPeriod * StationCode + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 9.08117 0.11195 81.120
## Year_fct2014 0.19256 0.11869 1.622
## Year_fct2015 0.63992 0.12388 5.166
## Year_fct2016 0.30110 0.13277 2.268
## Year_fct2017 0.34181 0.11904 2.871
## Year_fct2018 0.25789 0.11910 2.165
## Year_fct2019 0.22433 0.11923 1.881
## FlowActionPeriodDuring 0.10196 0.12080 0.844
## FlowActionPeriodAfter 1.15670 0.12833 9.013
## StationCodeSTTD -0.25842 0.15631 -1.653
## StationCodeLIB -1.57515 0.12106 -13.011
## StationCodeRVB -1.69525 0.11955 -14.180
## Year_fct2014:FlowActionPeriodDuring -0.65640 0.14774 -4.443
## Year_fct2015:FlowActionPeriodDuring -0.89633 0.13934 -6.433
## Year_fct2016:FlowActionPeriodDuring -0.40149 0.15420 -2.604
## Year_fct2017:FlowActionPeriodDuring -0.26619 0.14252 -1.868
## Year_fct2018:FlowActionPeriodDuring -0.69788 0.13815 -5.052
## Year_fct2019:FlowActionPeriodDuring -0.98774 0.13946 -7.083
## Year_fct2014:FlowActionPeriodAfter -1.67549 0.13620 -12.301
## Year_fct2015:FlowActionPeriodAfter -1.08617 0.14172 -7.664
## Year_fct2016:FlowActionPeriodAfter -1.79581 0.15340 -11.707
## Year_fct2017:FlowActionPeriodAfter -1.54569 0.13661 -11.315
## Year_fct2018:FlowActionPeriodAfter -1.23303 0.13656 -9.029
## Year_fct2019:FlowActionPeriodAfter -0.72868 0.13671 -5.330
## Year_fct2014:StationCodeSTTD -0.05655 0.16778 -0.337
## Year_fct2015:StationCodeSTTD -1.23388 0.17581 -7.018
## Year_fct2016:StationCodeSTTD -0.39579 0.18049 -2.193
## Year_fct2017:StationCodeSTTD -0.50185 0.16778 -2.991
## Year_fct2018:StationCodeSTTD -1.39676 0.16880 -8.275
## Year_fct2019:StationCodeSTTD -1.62927 0.17055 -9.553
## Year_fct2014:StationCodeLIB 0.20981 0.13555 1.548
## Year_fct2015:StationCodeLIB -0.60544 0.13956 -4.338
## Year_fct2016:StationCodeLIB 0.89636 0.14555 6.159
## Year_fct2017:StationCodeLIB -0.04825 0.13555 -0.356
## Year_fct2018:StationCodeLIB -1.58967 0.15175 -10.476
## Year_fct2019:StationCodeLIB -0.63431 0.14377 -4.412
## Year_fct2014:StationCodeRVB 0.14898 0.13433 1.109
## Year_fct2015:StationCodeRVB -0.29352 0.13803 -2.126
## Year_fct2016:StationCodeRVB 0.02945 0.14538 0.203
## Year_fct2017:StationCodeRVB -0.19731 0.13420 -1.470
## Year_fct2018:StationCodeRVB -0.18456 0.13420 -1.375
## Year_fct2019:StationCodeRVB -0.85616 0.13420 -6.380
## FlowActionPeriodDuring:StationCodeSTTD 0.80425 0.16883 4.764
## FlowActionPeriodAfter:StationCodeSTTD -0.98260 0.17209 -5.710
## FlowActionPeriodDuring:StationCodeLIB 0.21389 0.13685 1.563
## FlowActionPeriodAfter:StationCodeLIB -0.79587 0.13840 -5.750
## FlowActionPeriodDuring:StationCodeRVB 0.26545 0.13551 1.959
## FlowActionPeriodAfter:StationCodeRVB -1.12033 0.13711 -8.171
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.18071 0.20886 0.865
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 1.01021 0.19762 5.112
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.11730 0.21366 0.549
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.39969 0.24026 1.664
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 0.97183 0.19561 4.968
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 1.01560 0.19873 5.111
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 1.69478 0.19535 8.676
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 1.21670 0.20024 6.076
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 1.91011 0.20366 9.379
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.59088 0.22273 2.653
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 1.22365 0.20032 6.109
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.74205 0.19490 3.807
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.47435 0.18398 2.578
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 1.00159 0.16683 6.004
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.01240 0.18509 0.067
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.09854 0.17490 -0.563
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 0.79121 0.18111 4.369
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.52898 0.18151 2.914
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 1.31915 0.16307 8.089
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 0.80891 0.16713 4.840
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 1.10508 0.17155 6.442
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 1.01109 0.16307 6.200
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.28694 0.17679 -1.623
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB 0.18793 0.17087 1.100
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB 0.08281 0.18420 0.450
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.40198 0.16492 2.438
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.26354 0.18495 -1.425
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.37809 0.17386 -2.175
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.04862 0.16668 0.292
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 0.97561 0.16857 5.787
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 1.60050 0.16208 9.875
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 0.84464 0.16587 5.092
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 1.70707 0.17278 9.880
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 1.20633 0.16197 7.448
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 0.97985 0.16197 6.050
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 0.78129 0.16197 4.824
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2014 0.104840
## Year_fct2015 2.56e-07 ***
## Year_fct2016 0.023410 *
## Year_fct2017 0.004116 **
## Year_fct2018 0.030441 *
## Year_fct2019 0.060007 .
## FlowActionPeriodDuring 0.398757
## FlowActionPeriodAfter < 2e-16 ***
## StationCodeSTTD 0.098385 .
## StationCodeLIB < 2e-16 ***
## StationCodeRVB < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring 9.22e-06 ***
## Year_fct2015:FlowActionPeriodDuring 1.47e-10 ***
## Year_fct2016:FlowActionPeriodDuring 0.009273 **
## Year_fct2017:FlowActionPeriodDuring 0.061896 .
## Year_fct2018:FlowActionPeriodDuring 4.66e-07 ***
## Year_fct2019:FlowActionPeriodDuring 1.77e-12 ***
## Year_fct2014:FlowActionPeriodAfter < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter 2.45e-14 ***
## Year_fct2016:FlowActionPeriodAfter < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter < 2e-16 ***
## Year_fct2018:FlowActionPeriodAfter < 2e-16 ***
## Year_fct2019:FlowActionPeriodAfter 1.06e-07 ***
## Year_fct2014:StationCodeSTTD 0.736090
## Year_fct2015:StationCodeSTTD 2.80e-12 ***
## Year_fct2016:StationCodeSTTD 0.028398 *
## Year_fct2017:StationCodeSTTD 0.002803 **
## Year_fct2018:StationCodeSTTD < 2e-16 ***
## Year_fct2019:StationCodeSTTD < 2e-16 ***
## Year_fct2014:StationCodeLIB 0.121786
## Year_fct2015:StationCodeLIB 1.49e-05 ***
## Year_fct2016:StationCodeLIB 8.38e-10 ***
## Year_fct2017:StationCodeLIB 0.721881
## Year_fct2018:StationCodeLIB < 2e-16 ***
## Year_fct2019:StationCodeLIB 1.06e-05 ***
## Year_fct2014:StationCodeRVB 0.267522
## Year_fct2015:StationCodeRVB 0.033551 *
## Year_fct2016:StationCodeRVB 0.839460
## Year_fct2017:StationCodeRVB 0.141600
## Year_fct2018:StationCodeRVB 0.169155
## Year_fct2019:StationCodeRVB 2.06e-10 ***
## FlowActionPeriodDuring:StationCodeSTTD 2.00e-06 ***
## FlowActionPeriodAfter:StationCodeSTTD 1.25e-08 ***
## FlowActionPeriodDuring:StationCodeLIB 0.118179
## FlowActionPeriodAfter:StationCodeLIB 9.85e-09 ***
## FlowActionPeriodDuring:StationCodeRVB 0.050231 .
## FlowActionPeriodAfter:StationCodeRVB 4.54e-16 ***
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.386993
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 3.40e-07 ***
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.583046
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.096318 .
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 7.16e-07 ***
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 3.42e-07 ***
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 1.39e-09 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.008026 **
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 1.14e-09 ***
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.000143 ***
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.009980 **
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 2.17e-09 ***
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.946589
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB 0.573217
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 1.29e-05 ***
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.003591 **
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 8.78e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 1.37e-06 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 1.38e-10 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 6.45e-10 ***
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB 0.104698
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB 0.271498
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB 0.653052
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.014850 *
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB 0.154298
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB 0.029734 *
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.770545
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 7.92e-09 ***
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 3.77e-07 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 1.25e-13 ***
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 1.64e-09 ***
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 1.48e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.674 3.947 21.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.917 Deviance explained = 91.9%
## -REML = 690.81 Scale est. = 0.085513 n = 2932appraise(m_gam_cat3)
k.check(m_gam_cat3)## k' edf k-index p-value
## s(DOY) 4 3.674445 0.969169 0.0425draw(m_gam_cat3, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_cat3, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_cat3))
Box.test(residuals(m_gam_cat3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_cat3)
## X-squared = 4457.9, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_cat3 <- as.formula("Chla_log ~ Year_fct * FlowActionPeriod * StationCode + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat3 <- gamm(
f_gam_cat3,
data = df_chla_wt_c,
method = "REML"
)
m_gamm_cat3_ar1 <- gamm(
f_gam_cat3,
data = df_chla_wt_c,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat3_ar2 <- gamm(
f_gam_cat3,
data = df_chla_wt_c,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat3_ar3 <- gamm(
f_gam_cat3,
data = df_chla_wt_c,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_cat3$lme,
m_gamm_cat3_ar1$lme,
m_gamm_cat3_ar2$lme,
m_gamm_cat3_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio
## m_gamm_cat3$lme 1 87 1555.622 2073.622 -690.8109
## m_gamm_cat3_ar1$lme 2 88 -1932.180 -1408.226 1054.0899 1 vs 2 3489.802
## m_gamm_cat3_ar2$lme 3 89 -1930.971 -1401.063 1054.4856 2 vs 3 0.791
## m_gamm_cat3_ar3$lme 4 90 -1930.016 -1394.154 1055.0080 3 vs 4 1.045
## p-value
## m_gamm_cat3$lme
## m_gamm_cat3_ar1$lme <.0001
## m_gamm_cat3_ar2$lme 0.3737
## m_gamm_cat3_ar3$lme 0.3067It looks like the AR(1) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(1) Model
# Pull out the residuals and the GAM model
resid_cat3_ar1 <- residuals(m_gamm_cat3_ar1$lme, type = "normalized")
m_gamm_cat3_ar1_gam <- m_gamm_cat3_ar1$gam
summary(m_gamm_cat3_ar1_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * FlowActionPeriod * StationCode + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error
## (Intercept) 9.463591 0.275846
## Year_fct2014 -0.305994 0.362737
## Year_fct2015 -0.039067 0.366802
## Year_fct2016 -0.583469 0.386983
## Year_fct2017 -0.389788 0.360791
## Year_fct2018 -0.083020 0.357441
## Year_fct2019 -0.273120 0.358662
## FlowActionPeriodDuring 0.011851 0.165058
## FlowActionPeriodAfter 0.116296 0.226462
## StationCodeSTTD -0.262267 0.389697
## StationCodeLIB -1.922016 0.359134
## StationCodeRVB -1.972146 0.353906
## Year_fct2014:FlowActionPeriodDuring -0.100874 0.229879
## Year_fct2015:FlowActionPeriodDuring -0.185263 0.230420
## Year_fct2016:FlowActionPeriodDuring 0.089598 0.231706
## Year_fct2017:FlowActionPeriodDuring 0.001276 0.229562
## Year_fct2018:FlowActionPeriodDuring -0.318753 0.229152
## Year_fct2019:FlowActionPeriodDuring 0.081118 0.229281
## Year_fct2014:FlowActionPeriodAfter -0.308063 0.310774
## Year_fct2015:FlowActionPeriodAfter 0.057162 0.314450
## Year_fct2016:FlowActionPeriodAfter -0.055906 0.315211
## Year_fct2017:FlowActionPeriodAfter 0.241987 0.310816
## Year_fct2018:FlowActionPeriodAfter -0.513149 0.310892
## Year_fct2019:FlowActionPeriodAfter 0.232759 0.310865
## Year_fct2014:StationCodeSTTD 0.121376 0.516288
## Year_fct2015:StationCodeSTTD -0.505324 0.519139
## Year_fct2016:StationCodeSTTD 0.473585 0.541212
## Year_fct2017:StationCodeSTTD 0.069676 0.529637
## Year_fct2018:StationCodeSTTD -1.169014 0.513617
## Year_fct2019:StationCodeSTTD -1.226391 0.514099
## Year_fct2014:StationCodeLIB 0.667012 0.489863
## Year_fct2015:StationCodeLIB 0.244754 0.486757
## Year_fct2016:StationCodeLIB 1.670439 0.503574
## Year_fct2017:StationCodeLIB 0.570171 0.486294
## Year_fct2018:StationCodeLIB -1.652782 0.502839
## Year_fct2019:StationCodeLIB -0.116978 0.505596
## Year_fct2014:StationCodeRVB 0.591941 0.486869
## Year_fct2015:StationCodeRVB 0.286957 0.481961
## Year_fct2016:StationCodeRVB 0.822577 0.505290
## Year_fct2017:StationCodeRVB 0.360495 0.482447
## Year_fct2018:StationCodeRVB -0.036471 0.477673
## Year_fct2019:StationCodeRVB -0.404221 0.479189
## FlowActionPeriodDuring:StationCodeSTTD 0.071263 0.233323
## FlowActionPeriodAfter:StationCodeSTTD -0.146511 0.319783
## FlowActionPeriodDuring:StationCodeLIB 0.097456 0.229446
## FlowActionPeriodAfter:StationCodeLIB 0.344124 0.312268
## FlowActionPeriodDuring:StationCodeRVB 0.025445 0.228706
## FlowActionPeriodAfter:StationCodeRVB -0.040395 0.311021
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.086674 0.325522
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 0.247102 0.325979
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD -0.157989 0.327673
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.182295 0.327414
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 0.383379 0.324930
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD -0.120457 0.325001
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 0.390595 0.441383
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 0.100276 0.444479
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 0.005327 0.445767
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD -2.160500 0.445992
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 0.493127 0.443459
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD -0.271730 0.441102
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB -0.041254 0.322247
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 0.111258 0.321784
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB -0.310961 0.323368
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.115647 0.321768
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 0.167424 0.323637
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB -0.514504 0.324512
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB -0.211145 0.433931
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB -0.563615 0.436438
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB -0.635066 0.437052
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB -0.656717 0.433931
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.123449 0.438709
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB -0.951986 0.438102
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB -0.104058 0.321761
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.087262 0.321125
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.183629 0.323383
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.096613 0.321241
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.239050 0.320609
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB -0.054098 0.320809
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 0.092093 0.433081
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB -0.218715 0.435504
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB -0.092014 0.437800
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB -0.429726 0.433034
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 0.359997 0.433034
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB -0.249515 0.433034
## t value Pr(>|t|)
## (Intercept) 34.307 < 2e-16 ***
## Year_fct2014 -0.844 0.398981
## Year_fct2015 -0.107 0.915188
## Year_fct2016 -1.508 0.131733
## Year_fct2017 -1.080 0.280068
## Year_fct2018 -0.232 0.816352
## Year_fct2019 -0.761 0.446423
## FlowActionPeriodDuring 0.072 0.942767
## FlowActionPeriodAfter 0.514 0.607617
## StationCodeSTTD -0.673 0.501000
## StationCodeLIB -5.352 9.40e-08 ***
## StationCodeRVB -5.573 2.75e-08 ***
## Year_fct2014:FlowActionPeriodDuring -0.439 0.660830
## Year_fct2015:FlowActionPeriodDuring -0.804 0.421451
## Year_fct2016:FlowActionPeriodDuring 0.387 0.699016
## Year_fct2017:FlowActionPeriodDuring 0.006 0.995566
## Year_fct2018:FlowActionPeriodDuring -1.391 0.164330
## Year_fct2019:FlowActionPeriodDuring 0.354 0.723520
## Year_fct2014:FlowActionPeriodAfter -0.991 0.321636
## Year_fct2015:FlowActionPeriodAfter 0.182 0.855764
## Year_fct2016:FlowActionPeriodAfter -0.177 0.859237
## Year_fct2017:FlowActionPeriodAfter 0.779 0.436307
## Year_fct2018:FlowActionPeriodAfter -1.651 0.098937 .
## Year_fct2019:FlowActionPeriodAfter 0.749 0.454072
## Year_fct2014:StationCodeSTTD 0.235 0.814153
## Year_fct2015:StationCodeSTTD -0.973 0.330443
## Year_fct2016:StationCodeSTTD 0.875 0.381623
## Year_fct2017:StationCodeSTTD 0.132 0.895346
## Year_fct2018:StationCodeSTTD -2.276 0.022917 *
## Year_fct2019:StationCodeSTTD -2.386 0.017120 *
## Year_fct2014:StationCodeLIB 1.362 0.173422
## Year_fct2015:StationCodeLIB 0.503 0.615125
## Year_fct2016:StationCodeLIB 3.317 0.000921 ***
## Year_fct2017:StationCodeLIB 1.172 0.241101
## Year_fct2018:StationCodeLIB -3.287 0.001025 **
## Year_fct2019:StationCodeLIB -0.231 0.817046
## Year_fct2014:StationCodeRVB 1.216 0.224158
## Year_fct2015:StationCodeRVB 0.595 0.551626
## Year_fct2016:StationCodeRVB 1.628 0.103651
## Year_fct2017:StationCodeRVB 0.747 0.454991
## Year_fct2018:StationCodeRVB -0.076 0.939144
## Year_fct2019:StationCodeRVB -0.844 0.398991
## FlowActionPeriodDuring:StationCodeSTTD 0.305 0.760064
## FlowActionPeriodAfter:StationCodeSTTD -0.458 0.646874
## FlowActionPeriodDuring:StationCodeLIB 0.425 0.671054
## FlowActionPeriodAfter:StationCodeLIB 1.102 0.270548
## FlowActionPeriodDuring:StationCodeRVB 0.111 0.911422
## FlowActionPeriodAfter:StationCodeRVB -0.130 0.896670
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.266 0.790058
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 0.758 0.448496
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD -0.482 0.629733
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.557 0.577727
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 1.180 0.238145
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD -0.371 0.710936
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 0.885 0.376267
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 0.226 0.821526
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 0.012 0.990467
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD -4.844 1.34e-06 ***
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 1.112 0.266232
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD -0.616 0.537928
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB -0.128 0.898142
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 0.346 0.729554
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB -0.962 0.336316
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.359 0.719315
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 0.517 0.604971
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB -1.585 0.112971
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB -0.487 0.626589
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB -1.291 0.196671
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB -1.453 0.146315
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB -1.513 0.130286
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.281 0.778430
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB -2.173 0.029864 *
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB -0.323 0.746415
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.272 0.785842
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.568 0.570189
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.301 0.763627
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.746 0.455963
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB -0.169 0.866099
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 0.213 0.831618
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB -0.502 0.615558
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB -0.210 0.833547
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB -0.992 0.321107
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 0.831 0.405853
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB -0.576 0.564525
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 1 1 6.113 0.0135 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.844
## Scale est. = 0.2084 n = 2932appraise(m_gamm_cat3_ar1_gam)
k.check(m_gamm_cat3_ar1_gam)## k' edf k-index p-value
## s(DOY) 4 1.000003 0.8948027 0draw(m_gamm_cat3_ar1_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat3_ar1_gam, pages = 1, all.terms = TRUE)
acf(resid_cat3_ar1)
Box.test(resid_cat3_ar1, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_cat3_ar1
## X-squared = 81.422, df = 20, p-value = 2.249e-09The AR(1) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat3_ar1_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * FlowActionPeriod * StationCode + s(DOY,
## k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 6 0.674 0.671
## FlowActionPeriod 2 0.214 0.808
## StationCode 3 17.439 3.22e-11
## Year_fct:FlowActionPeriod 12 1.275 0.226
## Year_fct:StationCode 18 3.643 3.10e-07
## FlowActionPeriod:StationCode 6 0.788 0.579
## Year_fct:FlowActionPeriod:StationCode 36 3.809 2.52e-13
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 1 1 6.113 0.0135The 3-way interaction term is significant in this model, so we’ll use this one for the model selection process.
rm(m_gam_cat3, m_gamm_cat3, m_gamm_cat3_ar2, m_gamm_cat3_ar3)Model 2: LM 3-way interactions with water temperature
Initial Model
We’ll try running a linear model using a three-way interaction between Year, Flow Action Period, and Station, and a covariate of water temperature to account for seasonality and its effect of primary production. First we’ll run the linear model without accounting for serial autocorrelation.
m_lm_cat3_wt <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode + WaterTemp, data = .)Lets look at the model summary and diagnostics:
summary(m_lm_cat3_wt)##
## Call:
## lm(formula = Chla_log ~ Year_fct * FlowActionPeriod * StationCode +
## WaterTemp, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.94130 -0.15240 -0.01843 0.12832 1.67114
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 9.091520 0.146902
## Year_fct2014 0.198374 0.120464
## Year_fct2015 0.697112 0.125473
## Year_fct2016 0.499141 0.129685
## Year_fct2017 0.399696 0.120555
## Year_fct2018 0.322851 0.120469
## Year_fct2019 0.296743 0.120521
## FlowActionPeriodDuring -0.012914 0.121685
## FlowActionPeriodAfter 1.023870 0.127559
## StationCodeSTTD -0.258571 0.158710
## StationCodeLIB -1.498069 0.122869
## StationCodeRVB -1.599870 0.120735
## WaterTemp 0.001499 0.003927
## Year_fct2014:FlowActionPeriodDuring -0.690117 0.150044
## Year_fct2015:FlowActionPeriodDuring -0.951150 0.141255
## Year_fct2016:FlowActionPeriodDuring -0.368134 0.153919
## Year_fct2017:FlowActionPeriodDuring -0.320170 0.144456
## Year_fct2018:FlowActionPeriodDuring -0.767423 0.139823
## Year_fct2019:FlowActionPeriodDuring -1.053011 0.141021
## Year_fct2014:FlowActionPeriodAfter -1.683995 0.138584
## Year_fct2015:FlowActionPeriodAfter -1.134384 0.143849
## Year_fct2016:FlowActionPeriodAfter -1.898621 0.147625
## Year_fct2017:FlowActionPeriodAfter -1.616375 0.138232
## Year_fct2018:FlowActionPeriodAfter -1.290452 0.138359
## Year_fct2019:FlowActionPeriodAfter -0.806678 0.138200
## Year_fct2014:StationCodeSTTD -0.055493 0.170381
## Year_fct2015:StationCodeSTTD -1.239396 0.178545
## Year_fct2016:StationCodeSTTD -0.393709 0.183343
## Year_fct2017:StationCodeSTTD -0.500097 0.170420
## Year_fct2018:StationCodeSTTD -1.411304 0.171393
## Year_fct2019:StationCodeSTTD -1.662689 0.173064
## Year_fct2014:StationCodeLIB 0.134381 0.137309
## Year_fct2015:StationCodeLIB -0.640971 0.141666
## Year_fct2016:StationCodeLIB 0.878693 0.145294
## Year_fct2017:StationCodeLIB -0.121064 0.137150
## Year_fct2018:StationCodeLIB -1.741056 0.152288
## Year_fct2019:StationCodeLIB -0.648100 0.145964
## Year_fct2014:StationCodeRVB 0.057426 0.135687
## Year_fct2015:StationCodeRVB -0.347710 0.139954
## Year_fct2016:StationCodeRVB -0.004647 0.144935
## Year_fct2017:StationCodeRVB -0.288873 0.135453
## Year_fct2018:StationCodeRVB -0.277759 0.135461
## Year_fct2019:StationCodeRVB -0.948070 0.135443
## FlowActionPeriodDuring:StationCodeSTTD 0.803834 0.171429
## FlowActionPeriodAfter:StationCodeSTTD -0.982891 0.174735
## FlowActionPeriodDuring:StationCodeLIB 0.137972 0.138700
## FlowActionPeriodAfter:StationCodeLIB -0.840357 0.140413
## FlowActionPeriodDuring:StationCodeRVB 0.170461 0.136965
## FlowActionPeriodAfter:StationCodeRVB -1.184447 0.138447
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.180445 0.212074
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 1.016561 0.200656
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.117289 0.216941
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.440122 0.243874
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 0.986785 0.198609
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 1.049215 0.201689
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 1.673604 0.198443
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 1.238548 0.203291
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 1.909035 0.206807
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.574813 0.225989
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 1.193435 0.203360
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.775968 0.197801
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.549126 0.186513
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 1.031172 0.169338
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.034492 0.186073
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.023296 0.177319
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 0.943190 0.182408
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.521109 0.184439
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 1.386480 0.165554
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 0.826730 0.169756
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 1.092401 0.171920
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 1.052577 0.165167
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.168268 0.177917
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB 0.173850 0.173553
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB 0.172082 0.186487
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.456077 0.167277
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.224505 0.185728
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.283141 0.176032
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.142863 0.168659
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 1.070085 0.170604
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 1.658774 0.163731
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 0.882146 0.168205
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 1.707304 0.172795
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 1.268057 0.163613
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 1.041350 0.163607
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 0.842850 0.163608
## t value Pr(>|t|)
## (Intercept) 61.888 < 2e-16 ***
## Year_fct2014 1.647 0.099722 .
## Year_fct2015 5.556 3.02e-08 ***
## Year_fct2016 3.849 0.000121 ***
## Year_fct2017 3.315 0.000926 ***
## Year_fct2018 2.680 0.007406 **
## Year_fct2019 2.462 0.013869 *
## FlowActionPeriodDuring -0.106 0.915490
## FlowActionPeriodAfter 8.027 1.45e-15 ***
## StationCodeSTTD -1.629 0.103381
## StationCodeLIB -12.192 < 2e-16 ***
## StationCodeRVB -13.251 < 2e-16 ***
## WaterTemp 0.382 0.702751
## Year_fct2014:FlowActionPeriodDuring -4.599 4.42e-06 ***
## Year_fct2015:FlowActionPeriodDuring -6.734 2.00e-11 ***
## Year_fct2016:FlowActionPeriodDuring -2.392 0.016834 *
## Year_fct2017:FlowActionPeriodDuring -2.216 0.026745 *
## Year_fct2018:FlowActionPeriodDuring -5.489 4.41e-08 ***
## Year_fct2019:FlowActionPeriodDuring -7.467 1.08e-13 ***
## Year_fct2014:FlowActionPeriodAfter -12.151 < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter -7.886 4.41e-15 ***
## Year_fct2016:FlowActionPeriodAfter -12.861 < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter -11.693 < 2e-16 ***
## Year_fct2018:FlowActionPeriodAfter -9.327 < 2e-16 ***
## Year_fct2019:FlowActionPeriodAfter -5.837 5.91e-09 ***
## Year_fct2014:StationCodeSTTD -0.326 0.744674
## Year_fct2015:StationCodeSTTD -6.942 4.78e-12 ***
## Year_fct2016:StationCodeSTTD -2.147 0.031847 *
## Year_fct2017:StationCodeSTTD -2.934 0.003368 **
## Year_fct2018:StationCodeSTTD -8.234 2.72e-16 ***
## Year_fct2019:StationCodeSTTD -9.607 < 2e-16 ***
## Year_fct2014:StationCodeLIB 0.979 0.327822
## Year_fct2015:StationCodeLIB -4.525 6.30e-06 ***
## Year_fct2016:StationCodeLIB 6.048 1.66e-09 ***
## Year_fct2017:StationCodeLIB -0.883 0.377466
## Year_fct2018:StationCodeLIB -11.433 < 2e-16 ***
## Year_fct2019:StationCodeLIB -4.440 9.33e-06 ***
## Year_fct2014:StationCodeRVB 0.423 0.672162
## Year_fct2015:StationCodeRVB -2.484 0.013032 *
## Year_fct2016:StationCodeRVB -0.032 0.974424
## Year_fct2017:StationCodeRVB -2.133 0.033040 *
## Year_fct2018:StationCodeRVB -2.050 0.040410 *
## Year_fct2019:StationCodeRVB -7.000 3.19e-12 ***
## FlowActionPeriodDuring:StationCodeSTTD 4.689 2.87e-06 ***
## FlowActionPeriodAfter:StationCodeSTTD -5.625 2.04e-08 ***
## FlowActionPeriodDuring:StationCodeLIB 0.995 0.319942
## FlowActionPeriodAfter:StationCodeLIB -5.985 2.44e-09 ***
## FlowActionPeriodDuring:StationCodeRVB 1.245 0.213396
## FlowActionPeriodAfter:StationCodeRVB -8.555 < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.851 0.394919
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 5.066 4.32e-07 ***
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.541 0.588790
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 1.805 0.071226 .
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 4.968 7.15e-07 ***
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 5.202 2.11e-07 ***
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 8.434 < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 6.092 1.26e-09 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 9.231 < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 2.544 0.011026 *
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 5.869 4.90e-09 ***
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 3.923 8.95e-05 ***
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 2.944 0.003265 **
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 6.089 1.29e-09 ***
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.185 0.852953
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.131 0.895485
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 5.171 2.49e-07 ***
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 2.825 0.004756 **
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 8.375 < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 4.870 1.18e-06 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 6.354 2.43e-10 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 6.373 2.16e-10 ***
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.946 0.344347
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB 1.002 0.316570
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB 0.923 0.356214
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 2.726 0.006441 **
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -1.209 0.226847
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -1.608 0.107844
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.847 0.397037
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 6.272 4.10e-10 ***
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 10.131 < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 5.244 1.68e-07 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 9.881 < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 7.750 1.27e-14 ***
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 6.365 2.27e-10 ***
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 5.152 2.76e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2969 on 2843 degrees of freedom
## Multiple R-squared: 0.9168, Adjusted R-squared: 0.9143
## F-statistic: 372.8 on 84 and 2843 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp) %>%
plot_lm_diag(Chla_log, m_lm_cat3_wt)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat3_wt$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat3_wt$residuals
## W = 0.94155, p-value < 2.2e-16acf(residuals(m_lm_cat3_wt))
Box.test(residuals(m_lm_cat3_wt), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3_wt)
## X-squared = 4569.2, df = 20, p-value < 2.2e-16Hmmm, the residuals deviate from a normal distribution according to visual inspection and the Shapiro-Wilk normality test. Also, 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, 2, and 3 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_lm_cat3_wt_lag1 <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode + WaterTemp + lag1, data = .)
acf(residuals(m_lm_cat3_wt_lag1))
Box.test(residuals(m_lm_cat3_wt_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3_wt_lag1)
## X-squared = 72.103, df = 20, p-value = 8.23e-08Lag 2
m_lm_cat3_wt_lag2 <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1, lag2) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode + WaterTemp + lag1 + lag2, data = .)
acf(residuals(m_lm_cat3_wt_lag2))
Box.test(residuals(m_lm_cat3_wt_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3_wt_lag2)
## X-squared = 86.229, df = 20, p-value = 3.358e-10Lag 3
m_lm_cat3_wt_lag3 <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1, lag2, lag3) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode + WaterTemp + lag1 + lag2 + lag3, data = .)
acf(residuals(m_lm_cat3_wt_lag3))
Box.test(residuals(m_lm_cat3_wt_lag3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3_wt_lag3)
## X-squared = 61.572, df = 20, p-value = 4.056e-06The model with 3 lag terms looks to have slightly better ACF and Box-Ljung test results than the model with 1 lag term. Let’s compare the 3 models using AIC and BIC to see which model those prefer.
Compare models
AIC(m_lm_cat3_wt_lag1, m_lm_cat3_wt_lag2, m_lm_cat3_wt_lag3)## df AIC
## m_lm_cat3_wt_lag1 87 -2236.312
## m_lm_cat3_wt_lag2 88 -2197.327
## m_lm_cat3_wt_lag3 89 -2156.275BIC(m_lm_cat3_wt_lag1, m_lm_cat3_wt_lag2, m_lm_cat3_wt_lag3)## df BIC
## m_lm_cat3_wt_lag1 87 -1716.947
## m_lm_cat3_wt_lag2 88 -1673.096
## m_lm_cat3_wt_lag3 89 -1627.216Both AIC and BIC prefer the model with 1 lag term. Let’s take a closer look at that one.
Lag 1 model summary
summary(m_lm_cat3_wt_lag1)##
## Call:
## lm(formula = Chla_log ~ Year_fct * FlowActionPeriod * StationCode +
## WaterTemp + lag1, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.53589 -0.05777 -0.00825 0.04708 1.32305
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 1.5433727 0.1254039
## Year_fct2014 -0.0182347 0.0704237
## Year_fct2015 0.1095640 0.0734276
## Year_fct2016 0.0559093 0.0756029
## Year_fct2017 0.0218486 0.0705498
## Year_fct2018 0.0187731 0.0704786
## Year_fct2019 0.0192824 0.0704744
## FlowActionPeriodDuring -0.0292730 0.0709500
## FlowActionPeriodAfter 0.1510594 0.0747653
## StationCodeSTTD -0.0524052 0.0935312
## StationCodeLIB -0.2784185 0.0732843
## StationCodeRVB -0.2988072 0.0723946
## WaterTemp 0.0006079 0.0021523
## lag1 0.8325996 0.0101960
## Year_fct2014:FlowActionPeriodDuring -0.0558688 0.0859745
## Year_fct2015:FlowActionPeriodDuring -0.1786278 0.0817444
## Year_fct2016:FlowActionPeriodDuring -0.0493003 0.0880415
## Year_fct2017:FlowActionPeriodDuring 0.0130437 0.0828428
## Year_fct2018:FlowActionPeriodDuring -0.0997253 0.0807290
## Year_fct2019:FlowActionPeriodDuring -0.1448880 0.0816976
## Year_fct2014:FlowActionPeriodAfter -0.2306039 0.0816099
## Year_fct2015:FlowActionPeriodAfter -0.1874236 0.0833563
## Year_fct2016:FlowActionPeriodAfter -0.3199999 0.0868441
## Year_fct2017:FlowActionPeriodAfter -0.2393589 0.0812347
## Year_fct2018:FlowActionPeriodAfter -0.1830236 0.0807085
## Year_fct2019:FlowActionPeriodAfter -0.0898510 0.0799562
## Year_fct2014:StationCodeSTTD 0.0213002 0.0995410
## Year_fct2015:StationCodeSTTD -0.2171526 0.1047109
## Year_fct2016:StationCodeSTTD -0.0668027 0.1067229
## Year_fct2017:StationCodeSTTD -0.0468161 0.0996974
## Year_fct2018:StationCodeSTTD -0.2103017 0.1011273
## Year_fct2019:StationCodeSTTD -0.2690573 0.1024073
## Year_fct2014:StationCodeLIB 0.0767255 0.0792581
## Year_fct2015:StationCodeLIB -0.1146177 0.0818759
## Year_fct2016:StationCodeLIB 0.1699111 0.0841153
## Year_fct2017:StationCodeLIB 0.0166567 0.0791619
## Year_fct2018:StationCodeLIB -0.2440582 0.0896066
## Year_fct2019:StationCodeLIB -0.0620441 0.0843249
## Year_fct2014:StationCodeRVB 0.0609742 0.0783740
## Year_fct2015:StationCodeRVB -0.0555450 0.0807707
## Year_fct2016:StationCodeRVB 0.0076851 0.0835560
## Year_fct2017:StationCodeRVB -0.0185683 0.0782872
## Year_fct2018:StationCodeRVB -0.0112261 0.0783011
## Year_fct2019:StationCodeRVB -0.1160530 0.0788572
## FlowActionPeriodDuring:StationCodeSTTD 0.1521880 0.1002826
## FlowActionPeriodAfter:StationCodeSTTD -0.1794923 0.1021000
## FlowActionPeriodDuring:StationCodeLIB 0.0466299 0.0798234
## FlowActionPeriodAfter:StationCodeLIB -0.1196394 0.0811585
## FlowActionPeriodDuring:StationCodeRVB 0.0533654 0.0788926
## FlowActionPeriodAfter:StationCodeRVB -0.1763340 0.0805508
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD -0.0137482 0.1210623
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 0.1967960 0.1157705
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.0205039 0.1239926
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.2718959 0.1377188
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 0.1457438 0.1145147
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 0.1824220 0.1161909
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 0.2637300 0.1153907
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 0.2500200 0.1173263
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 0.3573180 0.1202704
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.0286247 0.1317845
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 0.1651458 0.1171817
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.1043081 0.1139902
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.0219463 0.1051488
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 0.2156308 0.0966431
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.0238801 0.1049769
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.0680988 0.1001231
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 0.1074204 0.1039223
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.0431558 0.1047747
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 0.1679720 0.0950500
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 0.1468685 0.0966081
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 0.1581214 0.0981718
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 0.1552268 0.0943469
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.0911976 0.1010821
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB -0.0286357 0.0985368
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB -0.0275969 0.1056301
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.0905538 0.0950143
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.0038011 0.1049124
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.0951350 0.0994458
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB -0.0086012 0.0955471
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 0.1185233 0.0972339
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 0.2326492 0.0944915
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 0.1635210 0.0957945
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 0.2900616 0.0996603
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 0.1983165 0.0937732
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 0.1488385 0.0934833
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 0.0859221 0.0933149
## t value Pr(>|t|)
## (Intercept) 12.307 < 2e-16 ***
## Year_fct2014 -0.259 0.795710
## Year_fct2015 1.492 0.135776
## Year_fct2016 0.740 0.459658
## Year_fct2017 0.310 0.756819
## Year_fct2018 0.266 0.789977
## Year_fct2019 0.274 0.784405
## FlowActionPeriodDuring -0.413 0.679941
## FlowActionPeriodAfter 2.020 0.043432 *
## StationCodeSTTD -0.560 0.575322
## StationCodeLIB -3.799 0.000148 ***
## StationCodeRVB -4.127 3.77e-05 ***
## WaterTemp 0.282 0.777636
## lag1 81.660 < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring -0.650 0.515855
## Year_fct2015:FlowActionPeriodDuring -2.185 0.028957 *
## Year_fct2016:FlowActionPeriodDuring -0.560 0.575547
## Year_fct2017:FlowActionPeriodDuring 0.157 0.874901
## Year_fct2018:FlowActionPeriodDuring -1.235 0.216819
## Year_fct2019:FlowActionPeriodDuring -1.773 0.076260 .
## Year_fct2014:FlowActionPeriodAfter -2.826 0.004751 **
## Year_fct2015:FlowActionPeriodAfter -2.248 0.024624 *
## Year_fct2016:FlowActionPeriodAfter -3.685 0.000233 ***
## Year_fct2017:FlowActionPeriodAfter -2.947 0.003240 **
## Year_fct2018:FlowActionPeriodAfter -2.268 0.023423 *
## Year_fct2019:FlowActionPeriodAfter -1.124 0.261214
## Year_fct2014:StationCodeSTTD 0.214 0.830575
## Year_fct2015:StationCodeSTTD -2.074 0.038186 *
## Year_fct2016:StationCodeSTTD -0.626 0.531402
## Year_fct2017:StationCodeSTTD -0.470 0.638691
## Year_fct2018:StationCodeSTTD -2.080 0.037655 *
## Year_fct2019:StationCodeSTTD -2.627 0.008653 **
## Year_fct2014:StationCodeLIB 0.968 0.333104
## Year_fct2015:StationCodeLIB -1.400 0.161655
## Year_fct2016:StationCodeLIB 2.020 0.043480 *
## Year_fct2017:StationCodeLIB 0.210 0.833361
## Year_fct2018:StationCodeLIB -2.724 0.006496 **
## Year_fct2019:StationCodeLIB -0.736 0.461930
## Year_fct2014:StationCodeRVB 0.778 0.436640
## Year_fct2015:StationCodeRVB -0.688 0.491706
## Year_fct2016:StationCodeRVB 0.092 0.926725
## Year_fct2017:StationCodeRVB -0.237 0.812533
## Year_fct2018:StationCodeRVB -0.143 0.886007
## Year_fct2019:StationCodeRVB -1.472 0.141218
## FlowActionPeriodDuring:StationCodeSTTD 1.518 0.129230
## FlowActionPeriodAfter:StationCodeSTTD -1.758 0.078856 .
## FlowActionPeriodDuring:StationCodeLIB 0.584 0.559157
## FlowActionPeriodAfter:StationCodeLIB -1.474 0.140555
## FlowActionPeriodDuring:StationCodeRVB 0.676 0.498823
## FlowActionPeriodAfter:StationCodeRVB -2.189 0.028671 *
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD -0.114 0.909592
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 1.700 0.089264 .
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.165 0.868669
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 1.974 0.048448 *
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 1.273 0.203227
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 1.570 0.116523
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 2.286 0.022356 *
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 2.131 0.033177 *
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 2.971 0.002994 **
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.217 0.828062
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 1.409 0.158853
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.915 0.360238
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.209 0.834685
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 2.231 0.025746 *
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.227 0.820067
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.680 0.496465
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 1.034 0.301384
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.412 0.680450
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 1.767 0.077304 .
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 1.520 0.128561
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 1.611 0.107366
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 1.645 0.100025
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.902 0.367021
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB -0.291 0.771371
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB -0.261 0.793911
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.953 0.340645
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.036 0.971100
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.957 0.338825
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB -0.090 0.928278
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 1.219 0.222966
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 2.462 0.013872 *
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 1.707 0.087933 .
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 2.911 0.003637 **
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 2.115 0.034530 *
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 1.592 0.111466
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 0.921 0.357247
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1619 on 2806 degrees of freedom
## Multiple R-squared: 0.9753, Adjusted R-squared: 0.9746
## F-statistic: 1305 on 85 and 2806 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1) %>%
plot_lm_diag(Chla_log, m_lm_cat3_wt_lag1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat3_wt_lag1$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat3_wt_lag1$residuals
## W = 0.82377, p-value < 2.2e-16Hmmm, the residuals deviate from a normal distribution particularly at the both tails of the distribution. I’m not so sure about using this model. Let’s take a look at its ANOVA table using type 3 sum of squares.
Anova(m_lm_cat3_wt_lag1, 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) 3.972 1 151.4675 < 2.2e-16 ***
## Year_fct 0.298 6 1.8941 0.07815 .
## FlowActionPeriod 0.518 2 9.8844 5.276e-05 ***
## StationCode 0.678 3 8.6139 1.087e-05 ***
## WaterTemp 0.002 1 0.0798 0.77764
## lag1 174.861 1 6668.2887 < 2.2e-16 ***
## Year_fct:FlowActionPeriod 1.550 12 4.9252 4.070e-08 ***
## Year_fct:StationCode 2.086 18 4.4201 1.458e-09 ***
## FlowActionPeriod:StationCode 0.981 6 6.2339 1.622e-06 ***
## Year_fct:FlowActionPeriod:StationCode 2.516 36 2.6656 3.306e-07 ***
## Residuals 73.581 2806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All the interaction terms are significant, so we can use this model in our model selection process; however, I’m somewhat leery of using this model since the model residuals don’t look right and the water temperature term isn’t significant.
rm(m_lm_cat3_wt, m_lm_cat3_wt_lag2, m_lm_cat3_wt_lag3)Model 3: LM 3-way interactions without seasonal term
Initial Model
We’ll try running a linear model using a three-way interaction between Year, Flow Action Period, and Station without including a term to account for seasonality. First we’ll run the linear model without accounting for serial autocorrelation.
m_lm_cat3 <- df_chla_wt_lag %>%
drop_na(Chla_log) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode, data = .)Lets look at the model summary and diagnostics:
summary(m_lm_cat3)##
## Call:
## lm(formula = Chla_log ~ Year_fct * FlowActionPeriod * StationCode,
## data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.94020 -0.15272 -0.01851 0.12896 1.67404
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 9.127697 0.112233
## Year_fct2014 0.198038 0.120470
## Year_fct2015 0.697361 0.125480
## Year_fct2016 0.501078 0.129596
## Year_fct2017 0.401506 0.120470
## Year_fct2018 0.322344 0.120470
## Year_fct2019 0.298192 0.120470
## FlowActionPeriodDuring -0.016983 0.121226
## FlowActionPeriodAfter 1.011768 0.123565
## StationCodeSTTD -0.258422 0.158722
## StationCodeLIB -1.502245 0.122390
## StationCodeRVB -1.602973 0.120470
## Year_fct2014:FlowActionPeriodDuring -0.688171 0.149969
## Year_fct2015:FlowActionPeriodDuring -0.949838 0.141224
## Year_fct2016:FlowActionPeriodDuring -0.363312 0.153412
## Year_fct2017:FlowActionPeriodDuring -0.317224 0.144261
## Year_fct2018:FlowActionPeriodDuring -0.766901 0.139827
## Year_fct2019:FlowActionPeriodDuring -1.051709 0.140990
## Year_fct2014:FlowActionPeriodAfter -1.680057 0.138209
## Year_fct2015:FlowActionPeriodAfter -1.131101 0.143602
## Year_fct2016:FlowActionPeriodAfter -1.890869 0.146232
## Year_fct2017:FlowActionPeriodAfter -1.615221 0.138209
## Year_fct2018:FlowActionPeriodAfter -1.287915 0.138209
## Year_fct2019:FlowActionPeriodAfter -0.806499 0.138209
## Year_fct2014:StationCodeSTTD -0.056554 0.170371
## Year_fct2015:StationCodeSTTD -1.240850 0.178518
## Year_fct2016:StationCodeSTTD -0.395790 0.183276
## Year_fct2017:StationCodeSTTD -0.501852 0.170371
## Year_fct2018:StationCodeSTTD -1.412505 0.171377
## Year_fct2019:StationCodeSTTD -1.663228 0.173071
## Year_fct2014:StationCodeLIB 0.136902 0.137161
## Year_fct2015:StationCodeLIB -0.638991 0.141581
## Year_fct2016:StationCodeLIB 0.877051 0.145241
## Year_fct2017:StationCodeLIB -0.121157 0.137161
## Year_fct2018:StationCodeLIB -1.740702 0.152297
## Year_fct2019:StationCodeLIB -0.647366 0.145962
## Year_fct2014:StationCodeRVB 0.059310 0.135607
## Year_fct2015:StationCodeRVB -0.346447 0.139925
## Year_fct2016:StationCodeRVB -0.006954 0.144820
## Year_fct2017:StationCodeRVB -0.289593 0.135450
## Year_fct2018:StationCodeRVB -0.276845 0.135450
## Year_fct2019:StationCodeRVB -0.948437 0.135450
## FlowActionPeriodDuring:StationCodeSTTD 0.804249 0.171439
## FlowActionPeriodAfter:StationCodeSTTD -0.982601 0.174747
## FlowActionPeriodDuring:StationCodeLIB 0.140988 0.138485
## FlowActionPeriodAfter:StationCodeLIB -0.837443 0.139886
## FlowActionPeriodDuring:StationCodeRVB 0.173168 0.136791
## FlowActionPeriodAfter:StationCodeRVB -1.181285 0.138209
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.180713 0.212088
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 1.017187 0.200664
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.117299 0.216957
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.441666 0.243859
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 0.987575 0.198613
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 1.049563 0.201702
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 1.676516 0.198311
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 1.238481 0.203307
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 1.910108 0.206803
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.578552 0.225793
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 1.197049 0.203154
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.776012 0.197816
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.547249 0.186463
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 1.029646 0.169304
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.031705 0.185944
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.025634 0.177227
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 0.942241 0.182404
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.517986 0.184272
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 1.360721 0.165032
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 0.825930 0.169574
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 1.093050 0.171806
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 1.052668 0.165032
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.167234 0.177811
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB 0.173617 0.173411
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB 0.170454 0.186452
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.454910 0.167262
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.227130 0.185615
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.285808 0.175906
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.140900 0.168594
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 1.067887 0.170520
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 1.658838 0.163743
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 0.881042 0.168193
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 1.707706 0.172804
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 1.267278 0.163613
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 1.040801 0.163613
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 0.842243 0.163613
## t value Pr(>|t|)
## (Intercept) 81.328 < 2e-16 ***
## Year_fct2014 1.644 0.100313
## Year_fct2015 5.558 2.99e-08 ***
## Year_fct2016 3.866 0.000113 ***
## Year_fct2017 3.333 0.000871 ***
## Year_fct2018 2.676 0.007500 **
## Year_fct2019 2.475 0.013373 *
## FlowActionPeriodDuring -0.140 0.888593
## FlowActionPeriodAfter 8.188 3.96e-16 ***
## StationCodeSTTD -1.628 0.103605
## StationCodeLIB -12.274 < 2e-16 ***
## StationCodeRVB -13.306 < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring -4.589 4.65e-06 ***
## Year_fct2015:FlowActionPeriodDuring -6.726 2.10e-11 ***
## Year_fct2016:FlowActionPeriodDuring -2.368 0.017940 *
## Year_fct2017:FlowActionPeriodDuring -2.199 0.027961 *
## Year_fct2018:FlowActionPeriodDuring -5.485 4.51e-08 ***
## Year_fct2019:FlowActionPeriodDuring -7.459 1.15e-13 ***
## Year_fct2014:FlowActionPeriodAfter -12.156 < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter -7.877 4.74e-15 ***
## Year_fct2016:FlowActionPeriodAfter -12.931 < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter -11.687 < 2e-16 ***
## Year_fct2018:FlowActionPeriodAfter -9.319 < 2e-16 ***
## Year_fct2019:FlowActionPeriodAfter -5.835 5.97e-09 ***
## Year_fct2014:StationCodeSTTD -0.332 0.739955
## Year_fct2015:StationCodeSTTD -6.951 4.48e-12 ***
## Year_fct2016:StationCodeSTTD -2.160 0.030892 *
## Year_fct2017:StationCodeSTTD -2.946 0.003249 **
## Year_fct2018:StationCodeSTTD -8.242 2.55e-16 ***
## Year_fct2019:StationCodeSTTD -9.610 < 2e-16 ***
## Year_fct2014:StationCodeLIB 0.998 0.318308
## Year_fct2015:StationCodeLIB -4.513 6.64e-06 ***
## Year_fct2016:StationCodeLIB 6.039 1.76e-09 ***
## Year_fct2017:StationCodeLIB -0.883 0.377135
## Year_fct2018:StationCodeLIB -11.430 < 2e-16 ***
## Year_fct2019:StationCodeLIB -4.435 9.55e-06 ***
## Year_fct2014:StationCodeRVB 0.437 0.661879
## Year_fct2015:StationCodeRVB -2.476 0.013346 *
## Year_fct2016:StationCodeRVB -0.048 0.961707
## Year_fct2017:StationCodeRVB -2.138 0.032601 *
## Year_fct2018:StationCodeRVB -2.044 0.041057 *
## Year_fct2019:StationCodeRVB -7.002 3.13e-12 ***
## FlowActionPeriodDuring:StationCodeSTTD 4.691 2.84e-06 ***
## FlowActionPeriodAfter:StationCodeSTTD -5.623 2.06e-08 ***
## FlowActionPeriodDuring:StationCodeLIB 1.018 0.308730
## FlowActionPeriodAfter:StationCodeLIB -5.987 2.41e-09 ***
## FlowActionPeriodDuring:StationCodeRVB 1.266 0.205643
## FlowActionPeriodAfter:StationCodeRVB -8.547 < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD 0.852 0.394251
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 5.069 4.25e-07 ***
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.541 0.588787
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 1.811 0.070223 .
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 4.972 7.01e-07 ***
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 5.204 2.09e-07 ***
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 8.454 < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 6.092 1.27e-09 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 9.236 < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 2.562 0.010449 *
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 5.892 4.26e-09 ***
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 3.923 8.96e-05 ***
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 2.935 0.003363 **
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 6.082 1.35e-09 ***
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.171 0.864623
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.145 0.885006
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 5.166 2.56e-07 ***
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 2.811 0.004973 **
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 8.245 2.49e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 4.871 1.17e-06 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 6.362 2.31e-10 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 6.379 2.08e-10 ***
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.941 0.347033
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB 1.001 0.316819
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB 0.914 0.360692
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 2.720 0.006573 **
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -1.224 0.221181
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -1.625 0.104321
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB 0.836 0.403371
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 6.263 4.36e-10 ***
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 10.131 < 2e-16 ***
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 5.238 1.74e-07 ***
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 9.882 < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 7.746 1.31e-14 ***
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 6.361 2.32e-10 ***
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 5.148 2.82e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2969 on 2848 degrees of freedom
## Multiple R-squared: 0.9167, Adjusted R-squared: 0.9142
## F-statistic: 377.4 on 83 and 2848 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log) %>%
plot_lm_diag(Chla_log, m_lm_cat3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat3$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat3$residuals
## W = 0.94201, p-value < 2.2e-16acf(residuals(m_lm_cat3))
Box.test(residuals(m_lm_cat3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3)
## X-squared = 4580.6, df = 20, p-value < 2.2e-16Hmmm, the residuals deviate from a normal distribution according to visual inspection and the Shapiro-Wilk normality test. Also, 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, 2, and 3 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_lm_cat3_lag1 <- df_chla_wt_lag %>%
drop_na(Chla_log, lag1) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode + lag1, data = .)
acf(residuals(m_lm_cat3_lag1))
Box.test(residuals(m_lm_cat3_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3_lag1)
## X-squared = 72.215, df = 20, p-value = 7.889e-08Lag 2
m_lm_cat3_lag2 <- df_chla_wt_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode + lag1 + lag2, data = .)
acf(residuals(m_lm_cat3_lag2))
Box.test(residuals(m_lm_cat3_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3_lag2)
## X-squared = 86.357, df = 20, p-value = 3.191e-10Lag 3
m_lm_cat3_lag3 <- df_chla_wt_lag %>%
drop_na(Chla_log, lag1, lag2, lag3) %>%
lm(Chla_log ~ Year_fct * FlowActionPeriod * StationCode + lag1 + lag2 + lag3, data = .)
acf(residuals(m_lm_cat3_lag3))
Box.test(residuals(m_lm_cat3_lag3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat3_lag3)
## X-squared = 61.887, df = 20, p-value = 3.622e-06The model with 3 lag terms looks to have slightly better ACF and Box-Ljung test results than the model with 1 lag term. Let’s compare the 3 models using AIC and BIC to see which model those prefer.
Compare models
AIC(m_lm_cat3_lag1, m_lm_cat3_lag2, m_lm_cat3_lag3)## df AIC
## m_lm_cat3_lag1 86 -2244.328
## m_lm_cat3_lag2 87 -2205.247
## m_lm_cat3_lag3 88 -2163.924BIC(m_lm_cat3_lag1, m_lm_cat3_lag2, m_lm_cat3_lag3)## df BIC
## m_lm_cat3_lag1 86 -1730.815
## m_lm_cat3_lag2 87 -1686.851
## m_lm_cat3_lag3 88 -1640.684Both AIC and BIC prefer the model with 1 lag term. Let’s take a closer look at that one.
Lag 1 model summary
summary(m_lm_cat3_lag1)##
## Call:
## lm(formula = Chla_log ~ Year_fct * FlowActionPeriod * StationCode +
## lag1, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.53572 -0.05749 -0.00790 0.04703 1.32211
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 1.555707 0.113899
## Year_fct2014 -0.018563 0.070371
## Year_fct2015 0.109392 0.073378
## Year_fct2016 0.056486 0.075510
## Year_fct2017 0.022339 0.070470
## Year_fct2018 0.018340 0.070421
## Year_fct2019 0.019667 0.070408
## FlowActionPeriodDuring -0.031037 0.070627
## FlowActionPeriodAfter 0.145762 0.072563
## StationCodeSTTD -0.052358 0.093469
## StationCodeLIB -0.279822 0.072957
## StationCodeRVB -0.299776 0.072177
## lag1 0.832869 0.010180
## Year_fct2014:FlowActionPeriodDuring -0.054753 0.085855
## Year_fct2015:FlowActionPeriodDuring -0.177763 0.081659
## Year_fct2016:FlowActionPeriodDuring -0.047178 0.087690
## Year_fct2017:FlowActionPeriodDuring 0.014467 0.082655
## Year_fct2018:FlowActionPeriodDuring -0.099169 0.080665
## Year_fct2019:FlowActionPeriodDuring -0.143954 0.081610
## Year_fct2014:FlowActionPeriodAfter -0.228415 0.081320
## Year_fct2015:FlowActionPeriodAfter -0.185703 0.083147
## Year_fct2016:FlowActionPeriodAfter -0.316281 0.086029
## Year_fct2017:FlowActionPeriodAfter -0.238324 0.081149
## Year_fct2018:FlowActionPeriodAfter -0.181508 0.080546
## Year_fct2019:FlowActionPeriodAfter -0.089433 0.079899
## Year_fct2014:StationCodeSTTD 0.020969 0.099467
## Year_fct2015:StationCodeSTTD -0.217358 0.104624
## Year_fct2016:StationCodeSTTD -0.067518 0.106612
## Year_fct2017:StationCodeSTTD -0.047307 0.099606
## Year_fct2018:StationCodeSTTD -0.210328 0.101048
## Year_fct2019:StationCodeSTTD -0.268761 0.102336
## Year_fct2014:StationCodeLIB 0.077840 0.079102
## Year_fct2015:StationCodeLIB -0.113562 0.081760
## Year_fct2016:StationCodeLIB 0.169094 0.084033
## Year_fct2017:StationCodeLIB 0.016771 0.079109
## Year_fct2018:StationCodeLIB -0.243330 0.089539
## Year_fct2019:StationCodeLIB -0.061421 0.084254
## Year_fct2014:StationCodeRVB 0.061880 0.078255
## Year_fct2015:StationCodeRVB -0.054830 0.080687
## Year_fct2016:StationCodeRVB 0.006864 0.083450
## Year_fct2017:StationCodeRVB -0.018640 0.078233
## Year_fct2018:StationCodeRVB -0.010627 0.078227
## Year_fct2019:StationCodeRVB -0.115784 0.078804
## FlowActionPeriodDuring:StationCodeSTTD 0.152225 0.100212
## FlowActionPeriodAfter:StationCodeSTTD -0.179035 0.102029
## FlowActionPeriodDuring:StationCodeLIB 0.047928 0.079631
## FlowActionPeriodAfter:StationCodeLIB -0.115393 0.080800
## FlowActionPeriodDuring:StationCodeRVB 0.054557 0.078719
## FlowActionPeriodAfter:StationCodeRVB -0.174594 0.080337
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD -0.013777 0.120982
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 0.196731 0.115691
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.020454 0.123911
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 0.272393 0.137614
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 0.145720 0.114435
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 0.182219 0.116113
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 0.264473 0.115236
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 0.249619 0.117247
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 0.357228 0.120179
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.029780 0.131615
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 0.166206 0.117001
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.104044 0.113914
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.020903 0.105032
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 0.214645 0.096544
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.022669 0.104820
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.069168 0.099986
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 0.106665 0.103838
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.041545 0.104576
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 0.159782 0.094671
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 0.143514 0.096436
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 0.155277 0.098042
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 0.152139 0.094204
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.093581 0.100951
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB -0.031651 0.098381
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB -0.028437 0.105523
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.089854 0.094928
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.004904 0.104761
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.096288 0.099287
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB -0.009588 0.095425
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 0.117176 0.097097
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 0.232073 0.094425
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 0.162732 0.095709
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 0.289664 0.099591
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 0.197522 0.093695
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 0.148185 0.093411
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 0.085282 0.093241
## t value Pr(>|t|)
## (Intercept) 13.659 < 2e-16 ***
## Year_fct2014 -0.264 0.791966
## Year_fct2015 1.491 0.136127
## Year_fct2016 0.748 0.454486
## Year_fct2017 0.317 0.751268
## Year_fct2018 0.260 0.794549
## Year_fct2019 0.279 0.780013
## FlowActionPeriodDuring -0.439 0.660365
## FlowActionPeriodAfter 2.009 0.044659 *
## StationCodeSTTD -0.560 0.575414
## StationCodeLIB -3.835 0.000128 ***
## StationCodeRVB -4.153 3.37e-05 ***
## lag1 81.816 < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring -0.638 0.523693
## Year_fct2015:FlowActionPeriodDuring -2.177 0.029572 *
## Year_fct2016:FlowActionPeriodDuring -0.538 0.590617
## Year_fct2017:FlowActionPeriodDuring 0.175 0.861073
## Year_fct2018:FlowActionPeriodDuring -1.229 0.219031
## Year_fct2019:FlowActionPeriodDuring -1.764 0.077853 .
## Year_fct2014:FlowActionPeriodAfter -2.809 0.005006 **
## Year_fct2015:FlowActionPeriodAfter -2.233 0.025598 *
## Year_fct2016:FlowActionPeriodAfter -3.676 0.000241 ***
## Year_fct2017:FlowActionPeriodAfter -2.937 0.003343 **
## Year_fct2018:FlowActionPeriodAfter -2.253 0.024307 *
## Year_fct2019:FlowActionPeriodAfter -1.119 0.263097
## Year_fct2014:StationCodeSTTD 0.211 0.833047
## Year_fct2015:StationCodeSTTD -2.078 0.037844 *
## Year_fct2016:StationCodeSTTD -0.633 0.526587
## Year_fct2017:StationCodeSTTD -0.475 0.634869
## Year_fct2018:StationCodeSTTD -2.081 0.037482 *
## Year_fct2019:StationCodeSTTD -2.626 0.008679 **
## Year_fct2014:StationCodeLIB 0.984 0.325178
## Year_fct2015:StationCodeLIB -1.389 0.164955
## Year_fct2016:StationCodeLIB 2.012 0.044291 *
## Year_fct2017:StationCodeLIB 0.212 0.832124
## Year_fct2018:StationCodeLIB -2.718 0.006616 **
## Year_fct2019:StationCodeLIB -0.729 0.466060
## Year_fct2014:StationCodeRVB 0.791 0.429160
## Year_fct2015:StationCodeRVB -0.680 0.496854
## Year_fct2016:StationCodeRVB 0.082 0.934448
## Year_fct2017:StationCodeRVB -0.238 0.811696
## Year_fct2018:StationCodeRVB -0.136 0.891953
## Year_fct2019:StationCodeRVB -1.469 0.141872
## FlowActionPeriodDuring:StationCodeSTTD 1.519 0.128867
## FlowActionPeriodAfter:StationCodeSTTD -1.755 0.079412 .
## FlowActionPeriodDuring:StationCodeLIB 0.602 0.547303
## FlowActionPeriodAfter:StationCodeLIB -1.428 0.153367
## FlowActionPeriodDuring:StationCodeRVB 0.693 0.488333
## FlowActionPeriodAfter:StationCodeRVB -2.173 0.029842 *
## Year_fct2014:FlowActionPeriodDuring:StationCodeSTTD -0.114 0.909347
## Year_fct2015:FlowActionPeriodDuring:StationCodeSTTD 1.700 0.089150 .
## Year_fct2016:FlowActionPeriodDuring:StationCodeSTTD 0.165 0.868901
## Year_fct2017:FlowActionPeriodDuring:StationCodeSTTD 1.979 0.047869 *
## Year_fct2018:FlowActionPeriodDuring:StationCodeSTTD 1.273 0.202986
## Year_fct2019:FlowActionPeriodDuring:StationCodeSTTD 1.569 0.116685
## Year_fct2014:FlowActionPeriodAfter:StationCodeSTTD 2.295 0.021803 *
## Year_fct2015:FlowActionPeriodAfter:StationCodeSTTD 2.129 0.033341 *
## Year_fct2016:FlowActionPeriodAfter:StationCodeSTTD 2.972 0.002979 **
## Year_fct2017:FlowActionPeriodAfter:StationCodeSTTD 0.226 0.821010
## Year_fct2018:FlowActionPeriodAfter:StationCodeSTTD 1.421 0.155557
## Year_fct2019:FlowActionPeriodAfter:StationCodeSTTD 0.913 0.361132
## Year_fct2014:FlowActionPeriodDuring:StationCodeLIB 0.199 0.842263
## Year_fct2015:FlowActionPeriodDuring:StationCodeLIB 2.223 0.026275 *
## Year_fct2016:FlowActionPeriodDuring:StationCodeLIB 0.216 0.828796
## Year_fct2017:FlowActionPeriodDuring:StationCodeLIB -0.692 0.489134
## Year_fct2018:FlowActionPeriodDuring:StationCodeLIB 1.027 0.304402
## Year_fct2019:FlowActionPeriodDuring:StationCodeLIB 0.397 0.691196
## Year_fct2014:FlowActionPeriodAfter:StationCodeLIB 1.688 0.091570 .
## Year_fct2015:FlowActionPeriodAfter:StationCodeLIB 1.488 0.136815
## Year_fct2016:FlowActionPeriodAfter:StationCodeLIB 1.584 0.113355
## Year_fct2017:FlowActionPeriodAfter:StationCodeLIB 1.615 0.106423
## Year_fct2018:FlowActionPeriodAfter:StationCodeLIB -0.927 0.354012
## Year_fct2019:FlowActionPeriodAfter:StationCodeLIB -0.322 0.747686
## Year_fct2014:FlowActionPeriodDuring:StationCodeRVB -0.269 0.787573
## Year_fct2015:FlowActionPeriodDuring:StationCodeRVB 0.947 0.343953
## Year_fct2016:FlowActionPeriodDuring:StationCodeRVB -0.047 0.962665
## Year_fct2017:FlowActionPeriodDuring:StationCodeRVB -0.970 0.332229
## Year_fct2018:FlowActionPeriodDuring:StationCodeRVB -0.100 0.919970
## Year_fct2019:FlowActionPeriodDuring:StationCodeRVB 1.207 0.227614
## Year_fct2014:FlowActionPeriodAfter:StationCodeRVB 2.458 0.014041 *
## Year_fct2015:FlowActionPeriodAfter:StationCodeRVB 1.700 0.089190 .
## Year_fct2016:FlowActionPeriodAfter:StationCodeRVB 2.909 0.003660 **
## Year_fct2017:FlowActionPeriodAfter:StationCodeRVB 2.108 0.035108 *
## Year_fct2018:FlowActionPeriodAfter:StationCodeRVB 1.586 0.112766
## Year_fct2019:FlowActionPeriodAfter:StationCodeRVB 0.915 0.360457
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1618 on 2811 degrees of freedom
## Multiple R-squared: 0.9753, Adjusted R-squared: 0.9746
## F-statistic: 1323 on 84 and 2811 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log, lag1) %>%
plot_lm_diag(Chla_log, m_lm_cat3_lag1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat3_lag1$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat3_lag1$residuals
## W = 0.82379, p-value < 2.2e-16Hmmm, again the residuals deviate from a normal distribution particularly at the both tails of the distribution. I’m not so sure about using this model. Let’s take a look at its ANOVA table using type 3 sum of squares.
Anova(m_lm_cat3_lag1, 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) 4.886 1 186.5576 < 2.2e-16 ***
## Year_fct 0.299 6 1.9058 0.07623 .
## FlowActionPeriod 0.546 2 10.4251 3.084e-05 ***
## StationCode 0.690 3 8.7860 8.495e-06 ***
## lag1 175.299 1 6693.8549 < 2.2e-16 ***
## Year_fct:FlowActionPeriod 1.547 12 4.9233 4.106e-08 ***
## Year_fct:StationCode 2.081 18 4.4137 1.525e-09 ***
## FlowActionPeriod:StationCode 0.980 6 6.2376 1.606e-06 ***
## Year_fct:FlowActionPeriod:StationCode 2.519 36 2.6719 3.074e-07 ***
## Residuals 73.614 2811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All the interaction terms are significant, so we can use this model in our model selection process; however, I’m somewhat leery of using this model since the model residuals don’t look right.
rm(m_lm_cat3, m_lm_cat3_lag2, m_lm_cat3_lag3)Model 4: GAM 2-way interactions with s(DOY)
Initial Model
We’ll try running a GAM using all two-way interactions between Year, Flow Action Period, and Station, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_cat2 <- gam(
Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5),
data = df_chla_wt_c,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_cat2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.64975 0.06314 152.819 < 2e-16 ***
## Year_fct2014 -0.50817 0.06713 -7.570 5.00e-14 ***
## Year_fct2015 0.01230 0.06728 0.183 0.854968
## Year_fct2016 -0.39765 0.08211 -4.843 1.35e-06 ***
## Year_fct2017 -0.17198 0.06670 -2.579 0.009972 **
## Year_fct2018 -0.24509 0.06655 -3.683 0.000235 ***
## Year_fct2019 -0.29125 0.06660 -4.373 1.27e-05 ***
## FlowActionPeriodDuring -0.36091 0.06217 -5.806 7.11e-09 ***
## FlowActionPeriodAfter 0.32490 0.07297 4.452 8.81e-06 ***
## StationCodeSTTD -1.10053 0.06546 -16.811 < 2e-16 ***
## StationCodeLIB -2.14605 0.05927 -36.207 < 2e-16 ***
## StationCodeRVB -2.29306 0.05826 -39.358 < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring -0.41511 0.06845 -6.064 1.50e-09 ***
## Year_fct2015:FlowActionPeriodDuring -0.23567 0.05960 -3.954 7.86e-05 ***
## Year_fct2016:FlowActionPeriodDuring -0.36454 0.08077 -4.513 6.64e-06 ***
## Year_fct2017:FlowActionPeriodDuring -0.27042 0.06710 -4.030 5.72e-05 ***
## Year_fct2018:FlowActionPeriodDuring -0.12610 0.06249 -2.018 0.043693 *
## Year_fct2019:FlowActionPeriodDuring -0.24824 0.06405 -3.876 0.000109 ***
## Year_fct2014:FlowActionPeriodAfter -0.46589 0.05896 -7.902 3.86e-15 ***
## Year_fct2015:FlowActionPeriodAfter -0.30129 0.05891 -5.115 3.35e-07 ***
## Year_fct2016:FlowActionPeriodAfter -0.53779 0.08552 -6.288 3.69e-10 ***
## Year_fct2017:FlowActionPeriodAfter -0.63694 0.06061 -10.508 < 2e-16 ***
## Year_fct2018:FlowActionPeriodAfter -0.64681 0.06101 -10.602 < 2e-16 ***
## Year_fct2019:FlowActionPeriodAfter -0.19776 0.05999 -3.297 0.000991 ***
## Year_fct2014:StationCodeSTTD 0.90359 0.07400 12.211 < 2e-16 ***
## Year_fct2015:StationCodeSTTD -0.26654 0.06977 -3.821 0.000136 ***
## Year_fct2016:StationCodeSTTD 0.69399 0.07466 9.295 < 2e-16 ***
## Year_fct2017:StationCodeSTTD 0.22832 0.08222 2.777 0.005521 **
## Year_fct2018:StationCodeSTTD -0.50639 0.07249 -6.986 3.50e-12 ***
## Year_fct2019:StationCodeSTTD -0.90290 0.07074 -12.763 < 2e-16 ***
## Year_fct2014:StationCodeLIB 1.02330 0.06869 14.898 < 2e-16 ***
## Year_fct2015:StationCodeLIB 0.13438 0.06589 2.039 0.041498 *
## Year_fct2016:StationCodeLIB 1.55910 0.07022 22.202 < 2e-16 ***
## Year_fct2017:StationCodeLIB 0.51381 0.06752 7.610 3.70e-14 ***
## Year_fct2018:StationCodeLIB -1.44531 0.06880 -21.006 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.31303 0.07052 -4.439 9.39e-06 ***
## Year_fct2014:StationCodeRVB 0.94342 0.06862 13.749 < 2e-16 ***
## Year_fct2015:StationCodeRVB 0.27341 0.06541 4.180 3.01e-05 ***
## Year_fct2016:StationCodeRVB 0.81490 0.07118 11.448 < 2e-16 ***
## Year_fct2017:StationCodeRVB 0.32877 0.06731 4.885 1.09e-06 ***
## Year_fct2018:StationCodeRVB 0.31730 0.06584 4.819 1.52e-06 ***
## Year_fct2019:StationCodeRVB -0.21713 0.06630 -3.275 0.001069 **
## FlowActionPeriodDuring:StationCodeSTTD 1.49868 0.05041 29.731 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeSTTD 0.22609 0.04601 4.914 9.44e-07 ***
## FlowActionPeriodDuring:StationCodeLIB 0.81415 0.04847 16.798 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeLIB -0.08208 0.04215 -1.947 0.051590 .
## FlowActionPeriodDuring:StationCodeRVB 0.55212 0.04650 11.874 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeRVB -0.04128 0.04056 -1.018 0.308841
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.527 3.892 15.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.889 Deviance explained = 89.1%
## -REML = 1070.5 Scale est. = 0.11366 n = 2932appraise(m_gam_cat2)
k.check(m_gam_cat2)## k' edf k-index p-value
## s(DOY) 4 3.526722 1.026702 0.9175draw(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 = 6617.7, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_cat2 <- as.formula("Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat2 <- gamm(
f_gam_cat2,
data = df_chla_wt_c,
method = "REML"
)
m_gamm_cat2_ar1 <- gamm(
f_gam_cat2,
data = df_chla_wt_c,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat2_ar2 <- gamm(
f_gam_cat2,
data = df_chla_wt_c,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat2_ar3 <- gamm(
f_gam_cat2,
data = df_chla_wt_c,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_cat2$lme,
m_gamm_cat2_ar1$lme,
m_gamm_cat2_ar2$lme,
m_gamm_cat2_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio
## m_gamm_cat2$lme 1 51 2243.041 2547.337 -1070.521
## m_gamm_cat2_ar1$lme 2 52 -1907.890 -1597.627 1005.945 1 vs 2 4152.931
## m_gamm_cat2_ar2$lme 3 53 -1906.021 -1589.792 1006.011 2 vs 3 0.132
## m_gamm_cat2_ar3$lme 4 54 -1905.019 -1582.823 1006.509 3 vs 4 0.998
## p-value
## m_gamm_cat2$lme
## m_gamm_cat2_ar1$lme <.0001
## m_gamm_cat2_ar2$lme 0.7167
## m_gamm_cat2_ar3$lme 0.3179It looks like the AR(1) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(1) Model
# Pull out the residuals and the GAM model
resid_cat2_ar1 <- residuals(m_gamm_cat2_ar1$lme, type = "normalized")
m_gamm_cat2_ar1_gam <- m_gamm_cat2_ar1$gam
summary(m_gamm_cat2_ar1_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.400509 0.229508 40.959 < 2e-16
## Year_fct2014 -0.363899 0.299792 -1.214 0.224908
## Year_fct2015 -0.020134 0.299762 -0.067 0.946453
## Year_fct2016 -0.536829 0.319613 -1.680 0.093139
## Year_fct2017 -0.093840 0.297754 -0.315 0.752665
## Year_fct2018 -0.214040 0.294097 -0.728 0.466802
## Year_fct2019 -0.142536 0.295918 -0.482 0.630073
## FlowActionPeriodDuring 0.028091 0.098052 0.286 0.774521
## FlowActionPeriodAfter 0.319889 0.133116 2.403 0.016320
## StationCodeSTTD -0.228035 0.309700 -0.736 0.461602
## StationCodeLIB -1.754353 0.287274 -6.107 1.15e-09
## StationCodeRVB -1.972157 0.283865 -6.948 4.58e-12
## Year_fct2014:FlowActionPeriodDuring -0.123204 0.114744 -1.074 0.283035
## Year_fct2015:FlowActionPeriodDuring -0.078610 0.114363 -0.687 0.491900
## Year_fct2016:FlowActionPeriodDuring -0.088402 0.116556 -0.758 0.448243
## Year_fct2017:FlowActionPeriodDuring -0.001409 0.114895 -0.012 0.990214
## Year_fct2018:FlowActionPeriodDuring -0.126989 0.114737 -1.107 0.268483
## Year_fct2019:FlowActionPeriodDuring -0.070220 0.115683 -0.607 0.543896
## Year_fct2014:FlowActionPeriodAfter -0.248509 0.153022 -1.624 0.104482
## Year_fct2015:FlowActionPeriodAfter -0.119569 0.153648 -0.778 0.436515
## Year_fct2016:FlowActionPeriodAfter -0.251766 0.158549 -1.588 0.112412
## Year_fct2017:FlowActionPeriodAfter -0.519339 0.154236 -3.367 0.000769
## Year_fct2018:FlowActionPeriodAfter -0.345890 0.154314 -2.241 0.025071
## Year_fct2019:FlowActionPeriodAfter -0.106045 0.154616 -0.686 0.492857
## Year_fct2014:StationCodeSTTD 0.316809 0.411046 0.771 0.440925
## Year_fct2015:StationCodeSTTD -0.387762 0.400573 -0.968 0.333117
## Year_fct2016:StationCodeSTTD 0.480368 0.420543 1.142 0.253442
## Year_fct2017:StationCodeSTTD -0.549525 0.432411 -1.271 0.203888
## Year_fct2018:StationCodeSTTD -0.897047 0.407500 -2.201 0.027791
## Year_fct2019:StationCodeSTTD -1.320089 0.401163 -3.291 0.001012
## Year_fct2014:StationCodeLIB 0.620377 0.387973 1.599 0.109926
## Year_fct2015:StationCodeLIB 0.065587 0.379254 0.173 0.862712
## Year_fct2016:StationCodeLIB 1.360499 0.396236 3.434 0.000604
## Year_fct2017:StationCodeLIB 0.302380 0.383749 0.788 0.430783
## Year_fct2018:StationCodeLIB -1.624289 0.391022 -4.154 3.36e-05
## Year_fct2019:StationCodeLIB -0.591885 0.398336 -1.486 0.137416
## Year_fct2014:StationCodeRVB 0.661451 0.386592 1.711 0.087192
## Year_fct2015:StationCodeRVB 0.241599 0.375957 0.643 0.520520
## Year_fct2016:StationCodeRVB 0.753178 0.400012 1.883 0.059816
## Year_fct2017:StationCodeRVB 0.189157 0.381636 0.496 0.620180
## Year_fct2018:StationCodeRVB 0.182456 0.375881 0.485 0.627424
## Year_fct2019:StationCodeRVB -0.475199 0.377726 -1.258 0.208474
## FlowActionPeriodDuring:StationCodeSTTD 0.180882 0.087052 2.078 0.037810
## FlowActionPeriodAfter:StationCodeSTTD -0.319820 0.116710 -2.740 0.006176
## FlowActionPeriodDuring:StationCodeLIB 0.018425 0.086802 0.212 0.831920
## FlowActionPeriodAfter:StationCodeLIB -0.097061 0.115326 -0.842 0.400068
## FlowActionPeriodDuring:StationCodeRVB 0.019857 0.086058 0.231 0.817531
## FlowActionPeriodAfter:StationCodeRVB -0.113430 0.114426 -0.991 0.321626
##
## (Intercept) ***
## Year_fct2014
## Year_fct2015
## Year_fct2016 .
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring
## FlowActionPeriodAfter *
## StationCodeSTTD
## StationCodeLIB ***
## StationCodeRVB ***
## Year_fct2014:FlowActionPeriodDuring
## Year_fct2015:FlowActionPeriodDuring
## Year_fct2016:FlowActionPeriodDuring
## Year_fct2017:FlowActionPeriodDuring
## Year_fct2018:FlowActionPeriodDuring
## Year_fct2019:FlowActionPeriodDuring
## Year_fct2014:FlowActionPeriodAfter
## Year_fct2015:FlowActionPeriodAfter
## Year_fct2016:FlowActionPeriodAfter
## Year_fct2017:FlowActionPeriodAfter ***
## Year_fct2018:FlowActionPeriodAfter *
## Year_fct2019:FlowActionPeriodAfter
## Year_fct2014:StationCodeSTTD
## Year_fct2015:StationCodeSTTD
## Year_fct2016:StationCodeSTTD
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD *
## Year_fct2019:StationCodeSTTD **
## Year_fct2014:StationCodeLIB
## Year_fct2015:StationCodeLIB
## Year_fct2016:StationCodeLIB ***
## Year_fct2017:StationCodeLIB
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB
## Year_fct2014:StationCodeRVB .
## Year_fct2015:StationCodeRVB
## Year_fct2016:StationCodeRVB .
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## FlowActionPeriodDuring:StationCodeSTTD *
## FlowActionPeriodAfter:StationCodeSTTD **
## FlowActionPeriodDuring:StationCodeLIB
## FlowActionPeriodAfter:StationCodeLIB
## FlowActionPeriodDuring:StationCodeRVB
## FlowActionPeriodAfter:StationCodeRVB
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 2.493 2.493 6.19 0.00119 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.849
## Scale est. = 0.19107 n = 2932appraise(m_gamm_cat2_ar1_gam)
k.check(m_gamm_cat2_ar1_gam)## k' edf k-index p-value
## s(DOY) 4 2.492999 0.9032606 0draw(m_gamm_cat2_ar1_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat2_ar1_gam, pages = 1, all.terms = TRUE)
acf(resid_cat2_ar1)
Box.test(resid_cat2_ar1, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_cat2_ar1
## X-squared = 69.446, df = 20, p-value = 2.243e-07The AR(1) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat2_ar1_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 6 0.821 0.55369
## FlowActionPeriod 2 4.643 0.00970
## StationCode 3 25.306 3.72e-16
## Year_fct:FlowActionPeriod 12 2.563 0.00223
## Year_fct:StationCode 18 5.428 1.03e-12
## FlowActionPeriod:StationCode 6 6.812 3.45e-07
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 2.493 2.493 6.19 0.00119All the interaction terms are significant in this model, so we’ll use this one for the model selection process.
rm(m_gam_cat2, m_gamm_cat2, m_gamm_cat2_ar2, m_gamm_cat2_ar3)Model 5: LM 2-way interactions with water temperature
Initial Model
We’ll try running a linear model using all two-way interactions between Year, Flow Action Period, and Station, and a covariate of water temperature to account for seasonality and its effect of primary production. First we’ll run the linear model without accounting for serial autocorrelation.
m_lm_cat2_wt <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + WaterTemp, data = .)Lets look at the model summary and diagnostics:
summary(m_lm_cat2_wt)##
## Call:
## lm(formula = Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 +
## WaterTemp, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24466 -0.20454 -0.01048 0.18046 1.60342
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.593474 0.121134 79.197 < 2e-16
## Year_fct2014 -0.547554 0.067470 -8.116 7.09e-16
## Year_fct2015 0.031436 0.068067 0.462 0.644235
## Year_fct2016 -0.223022 0.071997 -3.098 0.001969
## Year_fct2017 -0.168358 0.067841 -2.482 0.013133
## Year_fct2018 -0.231084 0.067214 -3.438 0.000594
## Year_fct2019 -0.269036 0.067612 -3.979 7.09e-05
## FlowActionPeriodDuring -0.508660 0.057813 -8.798 < 2e-16
## FlowActionPeriodAfter 0.199869 0.064973 3.076 0.002117
## StationCodeSTTD -1.107208 0.066171 -16.733 < 2e-16
## StationCodeLIB -2.102456 0.060478 -34.764 < 2e-16
## StationCodeRVB -2.242817 0.058911 -38.072 < 2e-16
## WaterTemp 0.005785 0.004457 1.298 0.194396
## Year_fct2014:FlowActionPeriodDuring -0.391221 0.068867 -5.681 1.47e-08
## Year_fct2015:FlowActionPeriodDuring -0.255486 0.060197 -4.244 2.26e-05
## Year_fct2016:FlowActionPeriodDuring -0.324340 0.068169 -4.758 2.05e-06
## Year_fct2017:FlowActionPeriodDuring -0.265526 0.067783 -3.917 9.16e-05
## Year_fct2018:FlowActionPeriodDuring -0.123330 0.063191 -1.952 0.051071
## Year_fct2019:FlowActionPeriodDuring -0.262255 0.064611 -4.059 5.06e-05
## Year_fct2014:FlowActionPeriodAfter -0.441405 0.060053 -7.350 2.57e-13
## Year_fct2015:FlowActionPeriodAfter -0.323046 0.059840 -5.399 7.27e-08
## Year_fct2016:FlowActionPeriodAfter -0.668491 0.064841 -10.310 < 2e-16
## Year_fct2017:FlowActionPeriodAfter -0.662727 0.060473 -10.959 < 2e-16
## Year_fct2018:FlowActionPeriodAfter -0.661922 0.061961 -10.683 < 2e-16
## Year_fct2019:FlowActionPeriodAfter -0.230354 0.060291 -3.821 0.000136
## Year_fct2014:StationCodeSTTD 0.902569 0.074778 12.070 < 2e-16
## Year_fct2015:StationCodeSTTD -0.256846 0.070606 -3.638 0.000280
## Year_fct2016:StationCodeSTTD 0.702398 0.075613 9.289 < 2e-16
## Year_fct2017:StationCodeSTTD 0.243279 0.083021 2.930 0.003413
## Year_fct2018:StationCodeSTTD -0.519259 0.073214 -7.092 1.65e-12
## Year_fct2019:StationCodeSTTD -0.909889 0.071513 -12.723 < 2e-16
## Year_fct2014:StationCodeLIB 1.000167 0.069817 14.326 < 2e-16
## Year_fct2015:StationCodeLIB 0.126561 0.066662 1.899 0.057722
## Year_fct2016:StationCodeLIB 1.562525 0.071106 21.975 < 2e-16
## Year_fct2017:StationCodeLIB 0.495189 0.068360 7.244 5.57e-13
## Year_fct2018:StationCodeLIB -1.480353 0.069471 -21.309 < 2e-16
## Year_fct2019:StationCodeLIB -0.318747 0.071316 -4.470 8.14e-06
## Year_fct2014:StationCodeRVB 0.907615 0.069319 13.093 < 2e-16
## Year_fct2015:StationCodeRVB 0.260582 0.066115 3.941 8.30e-05
## Year_fct2016:StationCodeRVB 0.811548 0.072098 11.256 < 2e-16
## Year_fct2017:StationCodeRVB 0.305045 0.068161 4.475 7.92e-06
## Year_fct2018:StationCodeRVB 0.287415 0.066435 4.326 1.57e-05
## Year_fct2019:StationCodeRVB -0.242106 0.067040 -3.611 0.000310
## FlowActionPeriodDuring:StationCodeSTTD 1.507099 0.051010 29.545 < 2e-16
## FlowActionPeriodAfter:StationCodeSTTD 0.227376 0.046646 4.874 1.15e-06
## FlowActionPeriodDuring:StationCodeLIB 0.787165 0.048955 16.079 < 2e-16
## FlowActionPeriodAfter:StationCodeLIB -0.103418 0.043469 -2.379 0.017419
## FlowActionPeriodDuring:StationCodeRVB 0.524680 0.046880 11.192 < 2e-16
## FlowActionPeriodAfter:StationCodeRVB -0.066924 0.041547 -1.611 0.107337
##
## (Intercept) ***
## Year_fct2014 ***
## Year_fct2015
## Year_fct2016 **
## Year_fct2017 *
## Year_fct2018 ***
## Year_fct2019 ***
## FlowActionPeriodDuring ***
## FlowActionPeriodAfter **
## StationCodeSTTD ***
## StationCodeLIB ***
## StationCodeRVB ***
## WaterTemp
## Year_fct2014:FlowActionPeriodDuring ***
## Year_fct2015:FlowActionPeriodDuring ***
## Year_fct2016:FlowActionPeriodDuring ***
## Year_fct2017:FlowActionPeriodDuring ***
## Year_fct2018:FlowActionPeriodDuring .
## Year_fct2019:FlowActionPeriodDuring ***
## Year_fct2014:FlowActionPeriodAfter ***
## Year_fct2015:FlowActionPeriodAfter ***
## Year_fct2016:FlowActionPeriodAfter ***
## Year_fct2017:FlowActionPeriodAfter ***
## Year_fct2018:FlowActionPeriodAfter ***
## Year_fct2019:FlowActionPeriodAfter ***
## Year_fct2014:StationCodeSTTD ***
## Year_fct2015:StationCodeSTTD ***
## Year_fct2016:StationCodeSTTD ***
## Year_fct2017:StationCodeSTTD **
## Year_fct2018:StationCodeSTTD ***
## Year_fct2019:StationCodeSTTD ***
## Year_fct2014:StationCodeLIB ***
## Year_fct2015:StationCodeLIB .
## Year_fct2016:StationCodeLIB ***
## Year_fct2017:StationCodeLIB ***
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB ***
## Year_fct2014:StationCodeRVB ***
## Year_fct2015:StationCodeRVB ***
## Year_fct2016:StationCodeRVB ***
## Year_fct2017:StationCodeRVB ***
## Year_fct2018:StationCodeRVB ***
## Year_fct2019:StationCodeRVB ***
## FlowActionPeriodDuring:StationCodeSTTD ***
## FlowActionPeriodAfter:StationCodeSTTD ***
## FlowActionPeriodDuring:StationCodeLIB ***
## FlowActionPeriodAfter:StationCodeLIB *
## FlowActionPeriodDuring:StationCodeRVB ***
## FlowActionPeriodAfter:StationCodeRVB
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3409 on 2879 degrees of freedom
## Multiple R-squared: 0.8889, Adjusted R-squared: 0.8871
## F-statistic: 480 on 48 and 2879 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp) %>%
plot_lm_diag(Chla_log, m_lm_cat2_wt)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat2_wt$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat2_wt$residuals
## W = 0.97383, p-value < 2.2e-16acf(residuals(m_lm_cat2_wt))
Box.test(residuals(m_lm_cat2_wt), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2_wt)
## X-squared = 6771.9, df = 20, p-value < 2.2e-16Hmmm, the residuals deviate from a normal distribution according to visual inspection and the Shapiro-Wilk normality test. Also, 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, 2, and 3 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_lm_cat2_wt_lag1 <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + WaterTemp + lag1, data = .)
acf(residuals(m_lm_cat2_wt_lag1))
Box.test(residuals(m_lm_cat2_wt_lag1), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2_wt_lag1)
## X-squared = 58.222, df = 20, p-value = 1.337e-05Lag 2
m_lm_cat2_wt_lag2 <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1, lag2) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + WaterTemp + lag1 + lag2, data = .)
acf(residuals(m_lm_cat2_wt_lag2))
Box.test(residuals(m_lm_cat2_wt_lag2), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2_wt_lag2)
## X-squared = 70.972, df = 20, p-value = 1.263e-07Lag 3
m_lm_cat2_wt_lag3 <- df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1, lag2, lag3) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + WaterTemp + lag1 + lag2 + lag3, data = .)
acf(residuals(m_lm_cat2_wt_lag3))
Box.test(residuals(m_lm_cat2_wt_lag3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2_wt_lag3)
## X-squared = 51.51, df = 20, p-value = 0.0001342The model with 3 lag terms looks to have slightly better ACF and Box-Ljung test results than the model with 1 lag term. Let’s compare the 3 models using AIC and BIC to see which model those prefer.
Compare models
AIC(m_lm_cat2_wt_lag1, m_lm_cat2_wt_lag2, m_lm_cat2_wt_lag3)## df AIC
## m_lm_cat2_wt_lag1 51 -2211.062
## m_lm_cat2_wt_lag2 52 -2170.794
## m_lm_cat2_wt_lag3 53 -2130.150BIC(m_lm_cat2_wt_lag1, m_lm_cat2_wt_lag2, m_lm_cat2_wt_lag3)## df BIC
## m_lm_cat2_wt_lag1 51 -1906.607
## m_lm_cat2_wt_lag2 52 -1861.021
## m_lm_cat2_wt_lag3 53 -1815.092Both AIC and BIC prefer the model with 1 lag term. Let’s take a closer look at that one.
Lag 1 model summary
summary(m_lm_cat2_wt_lag1)##
## Call:
## lm(formula = Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 +
## WaterTemp + lag1, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.65329 -0.05723 -0.00618 0.04969 1.32146
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.191711 0.104387 11.416 < 2e-16
## Year_fct2014 -0.089576 0.033321 -2.688 0.007225
## Year_fct2015 -0.011064 0.033238 -0.333 0.739248
## Year_fct2016 -0.054635 0.035234 -1.551 0.121100
## Year_fct2017 -0.050797 0.033159 -1.532 0.125661
## Year_fct2018 -0.038796 0.032919 -1.179 0.238689
## Year_fct2019 -0.034912 0.033137 -1.054 0.292178
## FlowActionPeriodDuring -0.070377 0.028621 -2.459 0.013995
## FlowActionPeriodAfter 0.028026 0.031660 0.885 0.376103
## StationCodeSTTD -0.148629 0.033741 -4.405 1.10e-05
## StationCodeLIB -0.263814 0.035030 -7.531 6.72e-14
## StationCodeRVB -0.289045 0.035073 -8.241 2.57e-16
## WaterTemp 0.001135 0.002150 0.528 0.597647
## lag1 0.875941 0.008971 97.639 < 2e-16
## Year_fct2014:FlowActionPeriodDuring -0.039124 0.033709 -1.161 0.245884
## Year_fct2015:FlowActionPeriodDuring -0.036631 0.029396 -1.246 0.212812
## Year_fct2016:FlowActionPeriodDuring -0.020959 0.033251 -0.630 0.528531
## Year_fct2017:FlowActionPeriodDuring 0.009967 0.032981 0.302 0.762512
## Year_fct2018:FlowActionPeriodDuring -0.017099 0.030781 -0.555 0.578599
## Year_fct2019:FlowActionPeriodDuring -0.032189 0.031554 -1.020 0.307764
## Year_fct2014:FlowActionPeriodAfter -0.043399 0.029494 -1.471 0.141271
## Year_fct2015:FlowActionPeriodAfter -0.026004 0.029250 -0.889 0.374067
## Year_fct2016:FlowActionPeriodAfter -0.088595 0.032078 -2.762 0.005784
## Year_fct2017:FlowActionPeriodAfter -0.064878 0.030044 -2.159 0.030899
## Year_fct2018:FlowActionPeriodAfter -0.083066 0.030733 -2.703 0.006917
## Year_fct2019:FlowActionPeriodAfter -0.025718 0.029442 -0.874 0.382458
## Year_fct2014:StationCodeSTTD 0.130814 0.037118 3.524 0.000431
## Year_fct2015:StationCodeSTTD -0.012610 0.034173 -0.369 0.712157
## Year_fct2016:StationCodeSTTD 0.110179 0.037042 2.974 0.002960
## Year_fct2017:StationCodeSTTD 0.082878 0.040517 2.046 0.040894
## Year_fct2018:StationCodeSTTD -0.051509 0.035681 -1.444 0.148964
## Year_fct2019:StationCodeSTTD -0.110216 0.035476 -3.107 0.001910
## Year_fct2014:StationCodeLIB 0.135296 0.034940 3.872 0.000110
## Year_fct2015:StationCodeLIB 0.022980 0.032246 0.713 0.476120
## Year_fct2016:StationCodeLIB 0.208361 0.037105 5.615 2.15e-08
## Year_fct2017:StationCodeLIB 0.072614 0.033330 2.179 0.029438
## Year_fct2018:StationCodeLIB -0.182207 0.036080 -5.050 4.70e-07
## Year_fct2019:StationCodeLIB -0.031421 0.034726 -0.905 0.365637
## Year_fct2014:StationCodeRVB 0.130204 0.034546 3.769 0.000167
## Year_fct2015:StationCodeRVB 0.044107 0.032022 1.377 0.168504
## Year_fct2016:StationCodeRVB 0.116439 0.035649 3.266 0.001103
## Year_fct2017:StationCodeRVB 0.047096 0.033054 1.425 0.154322
## Year_fct2018:StationCodeRVB 0.048087 0.032206 1.493 0.135521
## Year_fct2019:StationCodeRVB -0.028393 0.032458 -0.875 0.381778
## FlowActionPeriodDuring:StationCodeSTTD 0.206326 0.028160 7.327 3.05e-13
## FlowActionPeriodAfter:StationCodeSTTD 0.012461 0.022847 0.545 0.585504
## FlowActionPeriodDuring:StationCodeLIB 0.097623 0.024806 3.935 8.50e-05
## FlowActionPeriodAfter:StationCodeLIB -0.028083 0.021099 -1.331 0.183283
## FlowActionPeriodDuring:StationCodeRVB 0.064487 0.023212 2.778 0.005502
## FlowActionPeriodAfter:StationCodeRVB -0.009380 0.020136 -0.466 0.641374
##
## (Intercept) ***
## Year_fct2014 **
## Year_fct2015
## Year_fct2016
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring *
## FlowActionPeriodAfter
## StationCodeSTTD ***
## StationCodeLIB ***
## StationCodeRVB ***
## WaterTemp
## lag1 ***
## Year_fct2014:FlowActionPeriodDuring
## Year_fct2015:FlowActionPeriodDuring
## Year_fct2016:FlowActionPeriodDuring
## Year_fct2017:FlowActionPeriodDuring
## Year_fct2018:FlowActionPeriodDuring
## Year_fct2019:FlowActionPeriodDuring
## Year_fct2014:FlowActionPeriodAfter
## Year_fct2015:FlowActionPeriodAfter
## Year_fct2016:FlowActionPeriodAfter **
## Year_fct2017:FlowActionPeriodAfter *
## Year_fct2018:FlowActionPeriodAfter **
## Year_fct2019:FlowActionPeriodAfter
## Year_fct2014:StationCodeSTTD ***
## Year_fct2015:StationCodeSTTD
## Year_fct2016:StationCodeSTTD **
## Year_fct2017:StationCodeSTTD *
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD **
## Year_fct2014:StationCodeLIB ***
## Year_fct2015:StationCodeLIB
## Year_fct2016:StationCodeLIB ***
## Year_fct2017:StationCodeLIB *
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB
## Year_fct2014:StationCodeRVB ***
## Year_fct2015:StationCodeRVB
## Year_fct2016:StationCodeRVB **
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## FlowActionPeriodDuring:StationCodeSTTD ***
## FlowActionPeriodAfter:StationCodeSTTD
## FlowActionPeriodDuring:StationCodeLIB ***
## FlowActionPeriodAfter:StationCodeLIB
## FlowActionPeriodDuring:StationCodeRVB **
## FlowActionPeriodAfter:StationCodeRVB
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1636 on 2842 degrees of freedom
## Multiple R-squared: 0.9745, Adjusted R-squared: 0.974
## F-statistic: 2215 on 49 and 2842 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log, WaterTemp, lag1) %>%
plot_lm_diag(Chla_log, m_lm_cat2_wt_lag1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat2_wt_lag1$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat2_wt_lag1$residuals
## W = 0.82111, p-value < 2.2e-16Hmmm, the residuals deviate from a normal distribution particularly at the both tails of the distribution. I’m not so sure about using this model. Let’s take a look at its ANOVA table using type 3 sum of squares.
Anova(m_lm_cat2_wt_lag1, 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) 3.490 1 130.3319 < 2.2e-16 ***
## Year_fct 0.305 6 1.8963 0.0777883 .
## FlowActionPeriod 0.404 2 7.5440 0.0005399 ***
## StationCode 1.947 3 24.2426 1.735e-15 ***
## WaterTemp 0.007 1 0.2786 0.5976472
## lag1 255.265 1 9533.3675 < 2.2e-16 ***
## Year_fct:FlowActionPeriod 0.650 12 2.0224 0.0190524 *
## Year_fct:StationCode 3.076 18 6.3821 8.929e-16 ***
## FlowActionPeriod:StationCode 1.964 6 12.2279 1.262e-13 ***
## Residuals 76.097 2842
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All the interaction terms are significant, so we can use this model in our model selection process; however, I’m somewhat leery of using this model since the model residuals don’t look right and the water temperature term isn’t significant.
rm(m_lm_cat2_wt, m_lm_cat2_wt_lag2, m_lm_cat2_wt_lag3)Model 6: LM 2-way interactions without seasonal term
Initial Model
We’ll try running a linear model using all two-way interactions between Year, Flow Action Period, and Station without including a term to account for seasonality. First we’ll run the linear model without accounting for serial autocorrelation.
m_lm_cat2 <- df_chla_wt_lag %>%
drop_na(Chla_log) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2, data = .)Lets look at the model summary and diagnostics:
summary(m_lm_cat2)##
## Call:
## lm(formula = Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2,
## data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23917 -0.20797 -0.01003 0.18141 1.59954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.72909 0.05993 162.332 < 2e-16 ***
## Year_fct2014 -0.54396 0.06742 -8.069 1.03e-15 ***
## Year_fct2015 0.03685 0.06794 0.542 0.587604
## Year_fct2016 -0.21205 0.07150 -2.966 0.003045 **
## Year_fct2017 -0.15756 0.06735 -2.339 0.019386 *
## Year_fct2018 -0.22989 0.06719 -3.422 0.000631 ***
## Year_fct2019 -0.25859 0.06717 -3.850 0.000121 ***
## FlowActionPeriodDuring -0.51916 0.05715 -9.084 < 2e-16 ***
## FlowActionPeriodAfter 0.15654 0.05490 2.851 0.004385 **
## StationCodeSTTD -1.11036 0.06611 -16.795 < 2e-16 ***
## StationCodeLIB -2.11223 0.05969 -35.389 < 2e-16 ***
## StationCodeRVB -2.25084 0.05854 -38.450 < 2e-16 ***
## Year_fct2014:FlowActionPeriodDuring -0.38921 0.06883 -5.655 1.71e-08 ***
## Year_fct2015:FlowActionPeriodDuring -0.25461 0.06017 -4.231 2.40e-05 ***
## Year_fct2016:FlowActionPeriodDuring -0.31300 0.06757 -4.632 3.78e-06 ***
## Year_fct2017:FlowActionPeriodDuring -0.26051 0.06764 -3.851 0.000120 ***
## Year_fct2018:FlowActionPeriodDuring -0.12634 0.06315 -2.001 0.045514 *
## Year_fct2019:FlowActionPeriodDuring -0.26402 0.06459 -4.088 4.47e-05 ***
## Year_fct2014:FlowActionPeriodAfter -0.43149 0.05866 -7.356 2.46e-13 ***
## Year_fct2015:FlowActionPeriodAfter -0.31598 0.05947 -5.313 1.16e-07 ***
## Year_fct2016:FlowActionPeriodAfter -0.64009 0.06079 -10.529 < 2e-16 ***
## Year_fct2017:FlowActionPeriodAfter -0.66061 0.06037 -10.943 < 2e-16 ***
## Year_fct2018:FlowActionPeriodAfter -0.65331 0.06145 -10.632 < 2e-16 ***
## Year_fct2019:FlowActionPeriodAfter -0.23425 0.06019 -3.892 0.000102 ***
## Year_fct2014:StationCodeSTTD 0.90351 0.07475 12.087 < 2e-16 ***
## Year_fct2015:StationCodeSTTD -0.26118 0.07051 -3.704 0.000216 ***
## Year_fct2016:StationCodeSTTD 0.69668 0.07547 9.232 < 2e-16 ***
## Year_fct2017:StationCodeSTTD 0.24158 0.08299 2.911 0.003629 **
## Year_fct2018:StationCodeSTTD -0.51797 0.07319 -7.077 1.84e-12 ***
## Year_fct2019:StationCodeSTTD -0.91093 0.07149 -12.742 < 2e-16 ***
## Year_fct2014:StationCodeLIB 0.99784 0.06933 14.392 < 2e-16 ***
## Year_fct2015:StationCodeLIB 0.12799 0.06656 1.923 0.054591 .
## Year_fct2016:StationCodeLIB 1.55122 0.07055 21.988 < 2e-16 ***
## Year_fct2017:StationCodeLIB 0.48912 0.06816 7.176 9.06e-13 ***
## Year_fct2018:StationCodeLIB -1.48134 0.06939 -21.349 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.32261 0.07119 -4.532 6.08e-06 ***
## Year_fct2014:StationCodeRVB 0.91193 0.06921 13.175 < 2e-16 ***
## Year_fct2015:StationCodeRVB 0.26138 0.06609 3.955 7.85e-05 ***
## Year_fct2016:StationCodeRVB 0.79955 0.07150 11.183 < 2e-16 ***
## Year_fct2017:StationCodeRVB 0.29718 0.06789 4.377 1.24e-05 ***
## Year_fct2018:StationCodeRVB 0.28667 0.06642 4.316 1.64e-05 ***
## Year_fct2019:StationCodeRVB -0.24807 0.06688 -3.709 0.000212 ***
## FlowActionPeriodDuring:StationCodeSTTD 1.51157 0.05088 29.706 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeSTTD 0.23411 0.04635 5.051 4.67e-07 ***
## FlowActionPeriodDuring:StationCodeLIB 0.79121 0.04878 16.220 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeLIB -0.09403 0.04234 -2.221 0.026453 *
## FlowActionPeriodDuring:StationCodeRVB 0.52855 0.04676 11.303 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeRVB -0.05663 0.04074 -1.390 0.164610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3408 on 2884 degrees of freedom
## Multiple R-squared: 0.8888, Adjusted R-squared: 0.887
## F-statistic: 490.7 on 47 and 2884 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log) %>%
plot_lm_diag(Chla_log, m_lm_cat2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat2$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat2$residuals
## W = 0.97416, p-value < 2.2e-16acf(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 = 6777.5, df = 20, p-value < 2.2e-16Hmmm, the residuals deviate from a normal distribution according to visual inspection and the Shapiro-Wilk normality test. Also, 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, 2, and 3 lag terms and compare how well they correct for autocorrelation.
Lag 1
m_lm_cat2_lag1 <- df_chla_wt_lag %>%
drop_na(Chla_log, lag1) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + 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 = 58.272, df = 20, p-value = 1.314e-05Lag 2
m_lm_cat2_lag2 <- df_chla_wt_lag %>%
drop_na(Chla_log, lag1, lag2) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + 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 = 70.963, df = 20, p-value = 1.267e-07Lag 3
m_lm_cat2_lag3 <- df_chla_wt_lag %>%
drop_na(Chla_log, lag1, lag2, lag3) %>%
lm(Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + lag1 + lag2 + lag3, data = .)
acf(residuals(m_lm_cat2_lag3))
Box.test(residuals(m_lm_cat2_lag3), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_lm_cat2_lag3)
## X-squared = 51.531, df = 20, p-value = 0.0001332The model with 3 lag terms looks to have slightly better ACF and Box-Ljung test results than the model with 1 lag term. Let’s compare the 3 models using AIC and BIC to see which model those prefer.
Compare models
AIC(m_lm_cat2_lag1, m_lm_cat2_lag2, m_lm_cat2_lag3)## df AIC
## m_lm_cat2_lag1 50 -2218.889
## m_lm_cat2_lag2 51 -2178.457
## m_lm_cat2_lag3 52 -2137.377BIC(m_lm_cat2_lag1, m_lm_cat2_lag2, m_lm_cat2_lag3)## df BIC
## m_lm_cat2_lag1 50 -1920.334
## m_lm_cat2_lag2 51 -1874.570
## m_lm_cat2_lag3 52 -1828.190Both AIC and BIC prefer the model with 1 lag term. Let’s take a closer look at that one.
Lag 1 model summary
summary(m_lm_cat2_lag1)##
## Call:
## lm(formula = Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 +
## lag1, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.65247 -0.05772 -0.00654 0.04937 1.31988
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.216102 0.092117 13.202 < 2e-16
## Year_fct2014 -0.088344 0.033278 -2.655 0.007981
## Year_fct2015 -0.009516 0.033158 -0.287 0.774134
## Year_fct2016 -0.051913 0.034969 -1.485 0.137771
## Year_fct2017 -0.048184 0.032910 -1.464 0.143268
## Year_fct2018 -0.038091 0.032891 -1.158 0.246923
## Year_fct2019 -0.032346 0.032904 -0.983 0.325669
## FlowActionPeriodDuring -0.072114 0.028315 -2.547 0.010921
## FlowActionPeriodAfter 0.020326 0.026837 0.757 0.448887
## StationCodeSTTD -0.149099 0.033701 -4.424 1.00e-05
## StationCodeLIB -0.264370 0.034736 -7.611 3.68e-14
## StationCodeRVB -0.290098 0.034937 -8.303 < 2e-16
## lag1 0.876122 0.008961 97.774 < 2e-16
## Year_fct2014:FlowActionPeriodDuring -0.038814 0.033673 -1.153 0.249140
## Year_fct2015:FlowActionPeriodDuring -0.036636 0.029373 -1.247 0.212412
## Year_fct2016:FlowActionPeriodDuring -0.018906 0.032956 -0.574 0.566237
## Year_fct2017:FlowActionPeriodDuring 0.010850 0.032896 0.330 0.741560
## Year_fct2018:FlowActionPeriodDuring -0.017809 0.030748 -0.579 0.562510
## Year_fct2019:FlowActionPeriodDuring -0.032701 0.031529 -1.037 0.299748
## Year_fct2014:FlowActionPeriodAfter -0.042201 0.028798 -1.465 0.142926
## Year_fct2015:FlowActionPeriodAfter -0.025477 0.029063 -0.877 0.380767
## Year_fct2016:FlowActionPeriodAfter -0.083856 0.030153 -2.781 0.005454
## Year_fct2017:FlowActionPeriodAfter -0.065266 0.029978 -2.177 0.029551
## Year_fct2018:FlowActionPeriodAfter -0.082123 0.030463 -2.696 0.007063
## Year_fct2019:FlowActionPeriodAfter -0.027320 0.029381 -0.930 0.352516
## Year_fct2014:StationCodeSTTD 0.130915 0.037089 3.530 0.000422
## Year_fct2015:StationCodeSTTD -0.013397 0.034117 -0.393 0.694580
## Year_fct2016:StationCodeSTTD 0.108933 0.036949 2.948 0.003222
## Year_fct2017:StationCodeSTTD 0.082475 0.040482 2.037 0.041710
## Year_fct2018:StationCodeSTTD -0.051148 0.035653 -1.435 0.151505
## Year_fct2019:StationCodeSTTD -0.110237 0.035452 -3.109 0.001893
## Year_fct2014:StationCodeLIB 0.133608 0.034692 3.851 0.000120
## Year_fct2015:StationCodeLIB 0.022159 0.032186 0.688 0.491212
## Year_fct2016:StationCodeLIB 0.204751 0.036804 5.563 2.89e-08
## Year_fct2017:StationCodeLIB 0.070248 0.033213 2.115 0.034510
## Year_fct2018:StationCodeLIB -0.183223 0.036022 -5.086 3.89e-07
## Year_fct2019:StationCodeLIB -0.033229 0.034651 -0.959 0.337653
## Year_fct2014:StationCodeRVB 0.130827 0.034487 3.794 0.000152
## Year_fct2015:StationCodeRVB 0.044126 0.032000 1.379 0.168021
## Year_fct2016:StationCodeRVB 0.113843 0.035334 3.222 0.001288
## Year_fct2017:StationCodeRVB 0.045416 0.032909 1.380 0.167671
## Year_fct2018:StationCodeRVB 0.047818 0.032185 1.486 0.137457
## Year_fct2019:StationCodeRVB -0.029589 0.032373 -0.914 0.360807
## FlowActionPeriodDuring:StationCodeSTTD 0.206980 0.028103 7.365 2.31e-13
## FlowActionPeriodAfter:StationCodeSTTD 0.013755 0.022690 0.606 0.544429
## FlowActionPeriodDuring:StationCodeLIB 0.098119 0.024725 3.968 7.41e-05
## FlowActionPeriodAfter:StationCodeLIB -0.025885 0.020547 -1.260 0.207847
## FlowActionPeriodDuring:StationCodeRVB 0.065104 0.023152 2.812 0.004956
## FlowActionPeriodAfter:StationCodeRVB -0.007414 0.019735 -0.376 0.707191
##
## (Intercept) ***
## Year_fct2014 **
## Year_fct2015
## Year_fct2016
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring *
## FlowActionPeriodAfter
## StationCodeSTTD ***
## StationCodeLIB ***
## StationCodeRVB ***
## lag1 ***
## Year_fct2014:FlowActionPeriodDuring
## Year_fct2015:FlowActionPeriodDuring
## Year_fct2016:FlowActionPeriodDuring
## Year_fct2017:FlowActionPeriodDuring
## Year_fct2018:FlowActionPeriodDuring
## Year_fct2019:FlowActionPeriodDuring
## Year_fct2014:FlowActionPeriodAfter
## Year_fct2015:FlowActionPeriodAfter
## Year_fct2016:FlowActionPeriodAfter **
## Year_fct2017:FlowActionPeriodAfter *
## Year_fct2018:FlowActionPeriodAfter **
## Year_fct2019:FlowActionPeriodAfter
## Year_fct2014:StationCodeSTTD ***
## Year_fct2015:StationCodeSTTD
## Year_fct2016:StationCodeSTTD **
## Year_fct2017:StationCodeSTTD *
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD **
## Year_fct2014:StationCodeLIB ***
## Year_fct2015:StationCodeLIB
## Year_fct2016:StationCodeLIB ***
## Year_fct2017:StationCodeLIB *
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB
## Year_fct2014:StationCodeRVB ***
## Year_fct2015:StationCodeRVB
## Year_fct2016:StationCodeRVB **
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## FlowActionPeriodDuring:StationCodeSTTD ***
## FlowActionPeriodAfter:StationCodeSTTD
## FlowActionPeriodDuring:StationCodeLIB ***
## FlowActionPeriodAfter:StationCodeLIB
## FlowActionPeriodDuring:StationCodeRVB **
## FlowActionPeriodAfter:StationCodeRVB
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1635 on 2847 degrees of freedom
## Multiple R-squared: 0.9745, Adjusted R-squared: 0.9741
## F-statistic: 2265 on 48 and 2847 DF, p-value: < 2.2e-16df_chla_wt_lag %>%
drop_na(Chla_log, lag1) %>%
plot_lm_diag(Chla_log, m_lm_cat2_lag1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
shapiro.test(m_lm_cat2_lag1$residuals)##
## Shapiro-Wilk normality test
##
## data: m_lm_cat2_lag1$residuals
## W = 0.82107, p-value < 2.2e-16Hmmm, again the residuals deviate from a normal distribution particularly at the both tails of the distribution. I’m not so sure about using this model. Let’s take a look at its ANOVA table using type 3 sum of squares.
Anova(m_lm_cat2_lag1, 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) 4.661 1 174.2870 < 2.2e-16 ***
## Year_fct 0.301 6 1.8764 0.0811165 .
## FlowActionPeriod 0.428 2 8.0014 0.0003426 ***
## StationCode 1.981 3 24.6957 9.017e-16 ***
## lag1 255.643 1 9559.7382 < 2.2e-16 ***
## Year_fct:FlowActionPeriod 0.654 12 2.0381 0.0179620 *
## Year_fct:StationCode 3.065 18 6.3671 9.975e-16 ***
## FlowActionPeriod:StationCode 1.958 6 12.2005 1.360e-13 ***
## Residuals 76.133 2847
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All the interaction terms are significant, so we can use this model in our model selection process; however, I’m somewhat leery of using this model since the model residuals don’t look right.
rm(m_lm_cat2, m_lm_cat2_lag2, m_lm_cat2_lag3)Model selection - GAM and LM models
As a summary, here are the 6 models we are comparing:
- Model 1 -
m_gamm_cat3_ar1- GAM 3-way interactions with s(DOY)
Formula: Chla_log ~ Year_fct * FlowActionPeriod * StationCode + s(DOY, k = 5) - Model 2 -
m_lm_cat3_wt_lag1- LM 3-way interactions with water temperature
Formula: Chla_log ~ Year_fct * FlowActionPeriod * StationCode + WaterTemp + lag1 - Model 3 -
m_lm_cat3_lag1- LM 3-way interactions without seasonal term
Formula: Chla_log ~ Year_fct * FlowActionPeriod * StationCode + lag1 - Model 4 -
m_gamm_cat2_ar1- GAM 2-way interactions with s(DOY)
Formula: Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5) - Model 5 -
m_lm_cat2_wt_lag1- LM 2-way interactions with water temperature
Formula: Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + WaterTemp + lag1 - Model 6 -
m_lm_cat2_lag1- LM 2-way interactions without seasonal term
Formula: Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + lag1
# AIC values
df_m_aic <-
AIC(
m_gamm_cat3_ar1$lme,
m_lm_cat3_wt_lag1,
m_lm_cat3_lag1,
m_gamm_cat2_ar1$lme,
m_lm_cat2_wt_lag1,
m_lm_cat2_lag1
) %>%
as_tibble(rownames = "Model")
# BIC values
df_m_bic <-
BIC(
m_gamm_cat3_ar1$lme,
m_lm_cat3_wt_lag1,
m_lm_cat3_lag1,
m_gamm_cat2_ar1$lme,
m_lm_cat2_wt_lag1,
m_lm_cat2_lag1
) %>%
as_tibble(rownames = "Model")
# Combine AIC and BIC
df_m_aic_bic <- left_join(df_m_aic, df_m_bic, by = join_by(Model, df))
# Sort by AIC
df_m_aic_bic %>% arrange(AIC)## # A tibble: 6 × 4
## Model df AIC BIC
## <chr> <dbl> <dbl> <dbl>
## 1 m_lm_cat3_lag1 86 -2244. -1731.
## 2 m_lm_cat3_wt_lag1 87 -2236. -1717.
## 3 m_lm_cat2_lag1 50 -2219. -1920.
## 4 m_lm_cat2_wt_lag1 51 -2211. -1907.
## 5 m_gamm_cat3_ar1$lme 88 -1932. -1408.
## 6 m_gamm_cat2_ar1$lme 52 -1908. -1598.# Sort by BIC
df_m_aic_bic %>% arrange(BIC)## # A tibble: 6 × 4
## Model df AIC BIC
## <chr> <dbl> <dbl> <dbl>
## 1 m_lm_cat2_lag1 50 -2219. -1920.
## 2 m_lm_cat2_wt_lag1 51 -2211. -1907.
## 3 m_lm_cat3_lag1 86 -2244. -1731.
## 4 m_lm_cat3_wt_lag1 87 -2236. -1717.
## 5 m_gamm_cat2_ar1$lme 52 -1908. -1598.
## 6 m_gamm_cat3_ar1$lme 88 -1932. -1408.According to AIC, model 3 (LM 3-way interactions without seasonal term) was the best model, while BIC preferred model 6 (LM 2-way interactions without seasonal term). We’ll select model 3 (LM 3-way interactions without seasonal term) because the 3-way interaction between categorical variables was significant in this model indicating that there are possible significant differences between Flow Pulse Periods within each Station-Year grouping.
Model Results
Model 3
Lets take a closer look at model 3: LM 3-way interactions without
seasonal term
Formula: Chla_log ~ Year_fct * FlowActionPeriod *
StationCode + lag1
Pairwise Contrasts
Tukey pairwise contrasts of Flow Pulse Period for each Station - Year combination:
em_lm_cat3_lag1 <- emmeans(m_lm_cat3_lag1, ~ FlowActionPeriod | StationCode * Year_fct)
pairs(em_lm_cat3_lag1)## StationCode = RD22, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During 0.03104 0.0706 2811 0.439 0.8990
## Before - After -0.14576 0.0726 2811 -2.009 0.1103
## During - After -0.17680 0.0391 2811 -4.527 <.0001
##
## StationCode = STTD, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During -0.12119 0.0711 2811 -1.705 0.2035
## Before - After 0.03327 0.0718 2811 0.463 0.8885
## During - After 0.15446 0.0384 2811 4.027 0.0002
##
## StationCode = LIB, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During -0.01689 0.0368 2811 -0.459 0.8903
## Before - After -0.03037 0.0361 2811 -0.842 0.6767
## During - After -0.01348 0.0345 2811 -0.390 0.9195
##
## StationCode = RVB, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During -0.02352 0.0348 2811 -0.677 0.7771
## Before - After 0.02883 0.0340 2811 0.849 0.6727
## During - After 0.05235 0.0347 2811 1.509 0.2869
##
## StationCode = RD22, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During 0.08579 0.0488 2811 1.757 0.1844
## Before - After 0.08265 0.0347 2811 2.385 0.0452
## During - After -0.00314 0.0481 2811 -0.065 0.9977
##
## StationCode = STTD, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During -0.05266 0.0483 2811 -1.090 0.5206
## Before - After -0.00278 0.0390 2811 -0.071 0.9972
## During - After 0.04987 0.0518 2811 0.962 0.6009
##
## StationCode = LIB, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During 0.01696 0.0482 2811 0.352 0.9342
## Before - After 0.03826 0.0340 2811 1.127 0.4975
## During - After 0.02131 0.0481 2811 0.443 0.8977
##
## StationCode = RVB, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During 0.05967 0.0512 2811 1.165 0.4743
## Before - After 0.02518 0.0342 2811 0.736 0.7418
## During - After -0.03450 0.0509 2811 -0.678 0.7762
##
## StationCode = RD22, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During 0.20880 0.0410 2811 5.092 <.0001
## Before - After 0.03994 0.0403 2811 0.991 0.5829
## During - After -0.16886 0.0367 2811 -4.601 <.0001
##
## StationCode = STTD, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During -0.14016 0.0423 2811 -3.311 0.0027
## Before - After -0.03064 0.0408 2811 -0.752 0.7327
## During - After 0.10951 0.0353 2811 3.100 0.0056
##
## StationCode = LIB, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During -0.05377 0.0360 2811 -1.495 0.2934
## Before - After 0.01182 0.0340 2811 0.348 0.9354
## During - After 0.06559 0.0359 2811 1.828 0.1607
##
## StationCode = RVB, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During 0.06439 0.0349 2811 1.846 0.1549
## Before - After 0.05180 0.0342 2811 1.514 0.2845
## During - After -0.01259 0.0346 2811 -0.364 0.9295
##
## StationCode = RD22, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During 0.07822 0.0520 2811 1.505 0.2888
## Before - After 0.17052 0.0442 2811 3.857 0.0003
## During - After 0.09230 0.0444 2811 2.078 0.0945
##
## StationCode = STTD, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During -0.09446 0.0522 2811 -1.811 0.1662
## Before - After -0.00767 0.0434 2811 -0.177 0.9829
## During - After 0.08679 0.0444 2811 1.954 0.1239
##
## StationCode = LIB, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During 0.00762 0.0443 2811 0.172 0.9839
## Before - After 0.13064 0.0345 2811 3.791 0.0005
## During - After 0.12302 0.0443 2811 2.778 0.0152
##
## StationCode = RVB, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During 0.02856 0.0459 2811 0.622 0.8080
## Before - After 0.05545 0.0381 2811 1.457 0.3119
## During - After 0.02689 0.0459 2811 0.586 0.8278
##
## StationCode = RD22, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During 0.01657 0.0429 2811 0.386 0.9212
## Before - After 0.09256 0.0345 2811 2.683 0.0201
## During - After 0.07599 0.0427 2811 1.780 0.1764
##
## StationCode = STTD, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During -0.40805 0.0847 2811 -4.820 <.0001
## Before - After 0.24182 0.0769 2811 3.143 0.0048
## During - After 0.64986 0.1097 2811 5.922 <.0001
##
## StationCode = LIB, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During 0.03781 0.0428 2811 0.883 0.6511
## Before - After 0.05582 0.0342 2811 1.634 0.2317
## During - After 0.01801 0.0427 2811 0.422 0.9065
##
## StationCode = RVB, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During 0.05830 0.0430 2811 1.355 0.3648
## Before - After 0.06964 0.0344 2811 2.027 0.1059
## During - After 0.01133 0.0426 2811 0.266 0.9618
##
## StationCode = RD22, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During 0.13021 0.0390 2811 3.340 0.0024
## Before - After 0.03575 0.0341 2811 1.049 0.5457
## During - After -0.09446 0.0383 2811 -2.466 0.0366
##
## StationCode = STTD, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During -0.16774 0.0408 2811 -4.110 0.0001
## Before - After 0.04857 0.0455 2811 1.068 0.5339
## During - After 0.21631 0.0486 2811 4.453 <.0001
##
## StationCode = LIB, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During -0.02439 0.0538 2811 -0.453 0.8931
## Before - After 0.24472 0.0524 2811 4.672 <.0001
## During - After 0.26911 0.0412 2811 6.529 <.0001
##
## StationCode = RVB, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During 0.08524 0.0384 2811 2.218 0.0684
## Before - After 0.06216 0.0342 2811 1.817 0.1642
## During - After -0.02308 0.0380 2811 -0.608 0.8159
##
## StationCode = RD22, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During 0.17499 0.0409 2811 4.278 0.0001
## Before - After -0.05633 0.0340 2811 -1.658 0.2217
## During - After -0.23132 0.0413 2811 -5.608 <.0001
##
## StationCode = STTD, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During -0.15945 0.0436 2811 -3.657 0.0008
## Before - After 0.01866 0.0380 2811 0.491 0.8754
## During - After 0.17811 0.0399 2811 4.461 <.0001
##
## StationCode = LIB, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During 0.08552 0.0551 2811 1.553 0.2665
## Before - After 0.09072 0.0457 2811 1.985 0.1161
## During - After 0.00520 0.0482 2811 0.108 0.9936
##
## StationCode = RVB, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During 0.00326 0.0395 2811 0.083 0.9962
## Before - After 0.03298 0.0340 2811 0.971 0.5950
## During - After 0.02972 0.0394 2811 0.755 0.7307
##
## P value adjustment: tukey method for comparing a family of 3 estimatesSummary Figure
df_chla_wt_lag %>%
drop_na(Chla, lag1) %>%
plot_model_summary(em_lm_cat3_lag1)
Observed daily average chlorophyll fluorescence values (boxplots) and model results (model means as red points ±95% confidence intervals as gray boxes) for the Flow Pulse Period comparisons by Station and Year. Different letters above boxplots identify statistically significant (p < 0.05) differences from a Tukey post-hoc test.
The model results don’t match the observed values well at all. Unfortunately, this linear model is not a good candidate to use for our analysis.
Model 1
Since the linear model using 3-way interactions doesn’t fit the
observed data that well, let’s take a closer look at the GAM version of
this model, model 1: GAM 3-way interactions with s(DOY)
Formula:
Chla_log ~ Year_fct * FlowActionPeriod * StationCode + s(DOY, k = 5)
Pairwise Contrasts
Tukey pairwise contrasts of Flow Pulse Period for each Station - Year combination:
# Add the model call back to the gam object so it works with emmeans
m_gamm_cat3_ar1_gam$call <- quote(
gamm(
f_gam_cat3,
data = df_chla_wt_c,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1),
method = "REML"
)
)
em_gamm_cat3_ar1 <- emmeans(m_gamm_cat3_ar1_gam, ~ FlowActionPeriod | StationCode * Year_fct)
pairs(em_gamm_cat3_ar1)## StationCode = RD22, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During -0.01185 0.165 2847 -0.072 0.9972
## Before - After -0.11630 0.226 2847 -0.514 0.8648
## During - After -0.10445 0.161 2847 -0.647 0.7942
##
## StationCode = STTD, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During -0.08311 0.165 2847 -0.504 0.8696
## Before - After 0.03022 0.226 2847 0.133 0.9902
## During - After 0.11333 0.161 2847 0.702 0.7625
##
## StationCode = LIB, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During -0.10931 0.160 2847 -0.684 0.7727
## Before - After -0.46042 0.216 2847 -2.128 0.0845
## During - After -0.35111 0.159 2847 -2.209 0.0699
##
## StationCode = RVB, Year_fct = 2013:
## contrast estimate SE df t.ratio p.value
## Before - During -0.03730 0.159 2847 -0.235 0.9700
## Before - After -0.07590 0.215 2847 -0.353 0.9335
## During - After -0.03861 0.159 2847 -0.243 0.9679
##
## StationCode = RD22, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During 0.08902 0.160 2847 0.555 0.8437
## Before - After 0.19177 0.214 2847 0.896 0.6428
## During - After 0.10274 0.160 2847 0.641 0.7975
##
## StationCode = STTD, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During -0.06891 0.161 2847 -0.427 0.9042
## Before - After -0.05232 0.218 2847 -0.240 0.9687
## During - After 0.01660 0.162 2847 0.102 0.9942
##
## StationCode = LIB, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During 0.03282 0.160 2847 0.205 0.9772
## Before - After 0.05879 0.214 2847 0.275 0.9593
## During - After 0.02597 0.160 2847 0.162 0.9856
##
## StationCode = RVB, Year_fct = 2014:
## contrast estimate SE df t.ratio p.value
## Before - During 0.16764 0.161 2847 1.044 0.5493
## Before - After 0.14007 0.214 2847 0.653 0.7905
## During - After -0.02757 0.160 2847 -0.172 0.9839
##
## StationCode = RD22, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During 0.17341 0.161 2847 1.077 0.5287
## Before - After -0.17346 0.219 2847 -0.791 0.7086
## During - After -0.34687 0.160 2847 -2.170 0.0766
##
## StationCode = STTD, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During -0.14495 0.161 2847 -0.898 0.6416
## Before - After -0.12722 0.219 2847 -0.580 0.8306
## During - After 0.01773 0.159 2847 0.111 0.9932
##
## StationCode = LIB, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During -0.03530 0.159 2847 -0.222 0.9732
## Before - After 0.04603 0.215 2847 0.214 0.9750
## During - After 0.08133 0.159 2847 0.512 0.8657
##
## StationCode = RVB, Year_fct = 2015:
## contrast estimate SE df t.ratio p.value
## Before - During 0.06070 0.159 2847 0.382 0.9226
## Before - After 0.08565 0.215 2847 0.399 0.9161
## During - After 0.02495 0.159 2847 0.157 0.9865
##
## StationCode = RD22, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During -0.10145 0.163 2847 -0.623 0.8075
## Before - After -0.06039 0.220 2847 -0.274 0.9593
## During - After 0.04106 0.161 2847 0.254 0.9649
##
## StationCode = STTD, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During -0.01472 0.163 2847 -0.090 0.9955
## Before - After 0.08079 0.220 2847 0.367 0.9284
## During - After 0.09552 0.161 2847 0.592 0.8245
##
## StationCode = LIB, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During 0.11206 0.160 2847 0.700 0.7633
## Before - After 0.23055 0.214 2847 1.077 0.5285
## During - After 0.11850 0.160 2847 0.741 0.7394
##
## StationCode = RVB, Year_fct = 2016:
## contrast estimate SE df t.ratio p.value
## Before - During 0.05674 0.161 2847 0.352 0.9339
## Before - After 0.07202 0.217 2847 0.331 0.9413
## During - After 0.01528 0.161 2847 0.095 0.9951
##
## StationCode = RD22, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During -0.01313 0.160 2847 -0.082 0.9963
## Before - After -0.35828 0.214 2847 -1.673 0.2157
## During - After -0.34516 0.160 2847 -2.159 0.0786
##
## StationCode = STTD, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During -0.26668 0.165 2847 -1.615 0.2396
## Before - After 1.94873 0.228 2847 8.548 <.0001
## During - After 2.21541 0.168 2847 13.223 <.0001
##
## StationCode = LIB, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During 0.00506 0.160 2847 0.032 0.9994
## Before - After -0.04569 0.214 2847 -0.213 0.9752
## During - After -0.05075 0.160 2847 -0.317 0.9460
##
## StationCode = RVB, Year_fct = 2017:
## contrast estimate SE df t.ratio p.value
## Before - During 0.05804 0.160 2847 0.363 0.9299
## Before - After 0.11184 0.214 2847 0.522 0.8605
## During - After 0.05380 0.160 2847 0.336 0.9395
##
## StationCode = RD22, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During 0.30690 0.159 2847 1.926 0.1315
## Before - After 0.39685 0.214 2847 1.851 0.1534
## During - After 0.08995 0.159 2847 0.565 0.8390
##
## StationCode = STTD, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During -0.14774 0.161 2847 -0.917 0.6295
## Before - After 0.05024 0.222 2847 0.227 0.9721
## During - After 0.19798 0.163 2847 1.215 0.4441
##
## StationCode = LIB, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During 0.04202 0.164 2847 0.257 0.9643
## Before - After 0.17618 0.222 2847 0.792 0.7080
## During - After 0.13416 0.161 2847 0.834 0.6818
##
## StationCode = RVB, Year_fct = 2018:
## contrast estimate SE df t.ratio p.value
## Before - During 0.04241 0.159 2847 0.266 0.9617
## Before - After 0.07725 0.214 2847 0.360 0.9309
## During - After 0.03484 0.159 2847 0.219 0.9740
##
## StationCode = RD22, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During -0.09297 0.159 2847 -0.583 0.8293
## Before - After -0.34906 0.214 2847 -1.629 0.2336
## During - After -0.25609 0.159 2847 -1.606 0.2434
##
## StationCode = STTD, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During -0.04377 0.161 2847 -0.272 0.9601
## Before - After 0.06919 0.217 2847 0.318 0.9457
## During - After 0.11296 0.160 2847 0.706 0.7602
##
## StationCode = LIB, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During 0.32408 0.170 2847 1.911 0.1358
## Before - After 0.25881 0.227 2847 1.138 0.4906
## During - After -0.06527 0.162 2847 -0.403 0.9146
##
## StationCode = RVB, Year_fct = 2019:
## contrast estimate SE df t.ratio p.value
## Before - During -0.06432 0.159 2847 -0.403 0.9143
## Before - After -0.05915 0.214 2847 -0.276 0.9589
## During - After 0.00517 0.159 2847 0.032 0.9994
##
## P value adjustment: tukey method for comparing a family of 3 estimatesSummary Figure
plot_model_summary(df_chla_wt_c, em_gamm_cat3_ar1)
Observed daily average chlorophyll fluorescence values (boxplots) and model results (model means as red points ±95% confidence intervals as gray boxes) for the Flow Pulse Period comparisons by Station and Year. Different letters above boxplots identify statistically significant (p < 0.05) differences from a Tukey post-hoc test.
The model means seem to match the observed values better than Model 3, but there is a large amount of uncertainty around the means resulting in the model not seeing any significant differences among the pairwise comparisons. Unfortunately, this linear model does not seem like a good candidate to use for our analysis as well.