Title:	Assess Water Quality Trends with Generalized Additive Models
Version:	1.5.1
Date:	2025-07-09
Description:	Assess Water Quality Trends for Long-Term Monitoring Data in Estuaries using Generalized Additive Models following Wood (2017) <doi:10.1201/9781315370279> and Error Propagation with Mixed-Effects Meta-Analysis following Sera et al. (2019) <doi:10.1002/sim.8362>. Methods are available for model fitting, assessment of fit, annual and seasonal trend tests, and visualization of results.
Depends:	R (≥ 3.5)
Imports:	dplyr, ggplot2, lubridate, mgcv, mixmeta, plotly, purrr, tibble, tidyr, viridisLite
License:	CC0
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.2
URL:	<https://github.com/tbep-tech/wqtrends/>, <https://tbep-tech.github.io/wqtrends/>
BugReports:	https://github.com/tbep-tech/wqtrends/issues
Suggests:	testthat (≥ 2.1.0), covr, english, knitr, rmarkdown
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2025-07-09 20:32:25 UTC; mbeck
Author:	Marcus Beck [aut, cre], Perry de Valpine [aut], Rebecca Murphy [aut], Ian Wren [aut], Ariella Chelsky [aut], Melissa Foley [aut], David Senn [aut]
Maintainer:	Marcus Beck <mbeck@tbep.org>
Repository:	CRAN
Date/Publication:	2025-07-09 20:50:02 UTC

Extract period (seasonal) averages from fitted GAM

Description

Extract period (seasonal) averages from fitted GAM

Usage

anlz_avgseason(mod, doystr = 1, doyend = 364, yromit = NULL)

Arguments

Value

A data frame of period averages

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
anlz_avgseason(mod, doystr = 90, doyend = 180)

Back-transform response variable

Description

Back-transform response variable after fitting GAM

Usage

anlz_backtrans(dat)

Arguments

Details

dat can be output from anlz_trans or anlz_prd

Value

dat with the value column back-transformed using info from the trans column

Examples

library(dplyr)

tomod <- rawdat %>% 
  filter(station %in% 34) %>% 
  filter(param %in% 'chl')
dat <- anlz_trans(tomod, trans = 'log10')
backtrans <- anlz_backtrans(dat)
head(backtrans)

mod <- anlz_gam(tomod, trans = 'log10')
dat <- anlz_prd(mod)
backtrans <- anlz_backtrans(dat)
head(backtrans)

Return summary statistics for GAM fits

Description

Return summary statistics for GAM fits

Usage

anlz_fit(mod)

Arguments

Details

Results show the overall summary of the model as Akaike Information Criterion (AIC), the generalized cross-validation score (GCV), and the R2 values. Lower values for AIC and GCV and higher values for R2 indicate improved model fit.

Value

A data.frame with summary statistics for GAM fits

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl')
  
mod <-  anlz_gam(tomod, trans = 'log10')
anlz_fit(mod)

Fit a generalized additive model to a water quality time series

Description

Fit a generalized additive model to a water quality time series

Usage

anlz_gam(moddat, kts = NULL, ...)

Arguments

Details

The model structure is as follows:

model S:

chl ~ s(cont_year, k = large)

The cont_year vector is measured as a continuous numeric variable for the annual effect (e.g., January 1st, 2000 is 2000.0, July 1st, 2000 is 2000.5, etc.) and doy is the day of year as a numeric value from 1 to 366. The function s models cont_year as a smoothed, non-linear variable. The optimal amount of smoothing on cont_year is determined by cross-validation as implemented in the mgcv package and an upper theoretical upper limit on the number of knots for k should be large enough to allow sufficient flexibility in the smoothing term. The upper limit of k was chosen as 12 times the number of years for the input data. If insufficient data are available to fit a model with the specified k, the number of knots is decreased until the data can be modelled, e.g., 11 times the number of years, 10 times the number of years, etc.

Value

a gam model object

Examples

library(dplyr)
tomod <- rawdat %>% 
  filter(station %in% 34) %>% 
  filter(param %in% 'chl')
anlz_gam(tomod, trans = 'log10')

Extract period (seasonal) metrics from fitted GAM

Description

Extract period (seasonal) metrics from fitted GAM

Usage

anlz_metseason(
  mod,
  metfun = mean,
  doystr = 1,
  doyend = 364,
  nsim = 10000,
  yromit = NULL,
  ...
)

Arguments

Details

This function estimates a metric of interest for a given seasonal period each year using results from a fitted GAM (i.e., from anlz_gam). The estimates are based on the predicted values for each seasonal period, with uncertainty of the metric based on repeated sampling of the predictions following uncertainty in the model coefficients.

Value

A data frame of period metrics

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
anlz_metseason(mod, mean, doystr = 90, doyend = 180, nsim = 100)

Fit a mixed meta-analysis regression model of trends

Description

Fit a mixed meta-analysis regression model of trends

Usage

anlz_mixmeta(metseason, yrstr = 2000, yrend = 2019)

Arguments

Details

Parameters are not back-transformed if the original GAM used a transformation of the response variable

Value

A list of mixmeta fitted model objects

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
metseason <- anlz_metseason(mod, doystr = 90, doyend = 180)
anlz_mixmeta(metseason, yrstr = 2016, yrend = 2019)

Estimate percent change trends from GAM results for selected time periods

Description

Estimate percent change trends from GAM results for selected time periods

Usage

anlz_perchg(mod, baseyr, testyr)

Arguments

Details

Working components of this function were taken from the gamDiff function in the baytrends package.

Value

A data frame of summary results for change between the years.

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl')
  
mod <- anlz_gam(tomod, trans = 'log10')
anlz_perchg(mod, baseyr = 1990, testyr = 2016)

Get predicted data from fitted GAMs across period of observation

Description

Get predicted data from fitted GAMs across period of observation

Usage

anlz_prd(mod, annual = FALSE)

Arguments

Value

a data.frame with predictions

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)
  
mod <- anlz_gam(tomod, trans = 'log10')
anlz_prd(mod)

Get predicted data from fitted GAMs across period of observation, every day

Description

Get predicted data from fitted GAMs across period of observation, every day

Usage

anlz_prdday(mod)

Arguments

Value

a data.frame with predictions

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
anlz_prdday(mod)

Get prediction matrix for a fitted GAM

Description

Get prediction matrix for a fitted GAM

Usage

anlz_prdmatrix(mod, doystr = 1, doyend = 364, avemat = FALSE)

Arguments

Details

Used internally by anlz_metseason, not to be used by itself

Value

a data.frame with predictors to use with the fitted GAM

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)
  
mod <- anlz_gam(tomod, trans = 'log10')
anlz_prdmatrix(mod, doystr = 90, doyend = 180)

Format p-values for show functions

Description

Format p-values for show functions

Usage

anlz_pvalformat(x)

Arguments

Value

p-value formatted as a text string, one of p < 0.001, 'p < 0.01', p < 0.05, or ns for not significant

Examples

anlz_pvalformat(0.05)

Return summary statistics for smoothers of GAMs

Description

Return summary statistics for smoothers of GAMs

Usage

anlz_smooth(mod)

Arguments

Details

Results show the individual effects of the modelled components of each model as the estimated degrees of freedom (edf), the reference degrees of freedom (Ref.df), the test statistic (F), and significance of the component (p-value). The significance of the component is in part based on the difference between edf and Ref.df.

Value

a data.frame with summary statistics for smoothers in each GAM

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl')
  
mod <- anlz_gam(tomod, trans = 'log10')
anlz_smooth(mod)

Retrieve summary statistics for seasonal metrics and trend results

Description

Retrieve summary statistics for seasonal metrics and trend results

Usage

anlz_sumstats(
  mod,
  metfun = mean,
  doystr = 1,
  doyend = 364,
  yrstr = 2000,
  yrend = 2019,
  yromit = NULL,
  nsim = 10000,
  confint = 0.95,
  useave = FALSE,
  ...
)

Arguments

Details

This function is primarily for convenience to return summary statistics of a fitted GAM from anlz_gam.

Note that confint only applies to the summary and coeffs list outputs. It does not apply to the metseason list element output that is default set to 95

Set useave = T to speed up calculations if metfun = mean. This will use anlz_avgseason to estimate the seasonal summary metrics using a non-stochastic equation.

Value

A list object with named elements:

• mixmet: mixmeta object of the fitted mixed-effects meta-analysis trend model

• metseason: tibble object of the fitted seasonal metrics as returned by anlz_metseason or anlz_avgseason

• summary: summary of the mixmet object

• coeffs: tibble object of the slope estimate coefficients from the mixmet model. An approximately linear slope estimate will be included as slope.approx if trans = 'log10' for the GAM used in mod.

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')

anlz_sumstats(mod, metfun = mean, doystr = 90, doyend = 180, yrstr = 2016, 
  yrend = 2019, nsim = 100)

Estimate seasonal rates of change based on average estimates for multiple window widths

Description

Estimate seasonal rates of change based on average estimates for multiple window widths

Usage

anlz_sumtrndseason(
  mod,
  doystr = 1,
  doyend = 364,
  justify = c("center", "left", "right"),
  win = 5:15,
  yromit = NULL
)

Arguments

Details

The optional yromit vector can be used to omit years from the plot and trend assessment. This may be preferred if seasonal estimates for a given year have very wide confidence intervals likely due to limited data, which can skew the trend assessments.

This function is a wrapper to anlz_trndseason to loop across values in win, using useave = TRUE for quicker calculation of average seasonal metrics. It does not work with any other seasonal metric calculations.

Value

A data frame of slope estimates and p-values for each year

See Also

Other analyze: anlz_trans(), anlz_trndseason()

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
anlz_sumtrndseason(mod, doystr = 90, doyend = 180, justify = 'center', win = 2:3)

Transform response variable

Description

Transform response variable prior to fitting GAM

Usage

anlz_trans(moddat, trans = c("log10", "ident"))

Arguments

Value

moddat with the value column transformed as indicated

See Also

Other analyze: anlz_sumtrndseason(), anlz_trndseason()

Examples

library(dplyr)
tomod <- rawdat %>% 
  filter(station %in% 34) %>% 
  filter(param %in% 'chl')
anlz_trans(tomod, trans = 'log10')

Estimate rates of change based on seasonal metrics

Description

Estimate rates of change based on seasonal metrics

Usage

anlz_trndseason(
  mod,
  metfun = mean,
  doystr = 1,
  doyend = 364,
  justify = c("center", "left", "right"),
  win = 5,
  nsim = 10000,
  yromit = NULL,
  useave = FALSE,
  ...
)

Arguments

Details

Trends are based on the slope of the fitted linear trend within the window, where the linear trend is estimated using a meta-analysis regression model (from anlz_mixmeta) for the seasonal metrics (from anlz_metseason). Set useave = T to speed up calculations if metfun = mean. This will use anlz_avgseason to estimate the seasonal summary metrics using a non-stochastic equation.

Note that for left and right windows, the exact number of years in win is used. For example, a left-centered window for 1990 of ten years will include exactly ten years from 1990, 1991, ... , 1999. The same applies to a right-centered window, e.g., 1990 would include 1981, 1982, ..., 1990 (if those years have data). However, for a centered window, picking an even number of years for the window width will create a slightly off-centered window because it is impossible to center on an even number of years. For example, if win = 8 and justify = 'center', the estimate for 2000 will be centered on 1997 to 2004 (three years left, four years right, eight years total). Centering for window widths with an odd number of years will always create a symmetrical window, i.e., if win = 7 and justify = 'center', the estimate for 2000 will be centered on 1997 and 2003 (three years left, three years right, seven years total).

The optional yromit vector can be used to omit years from the trend assessment. This may be preferred if seasonal estimates for a given year have very wide confidence intervals likely due to limited data, which can skew the trend assessments.

Value

A data frame of slope estimates and p-values for each year

See Also

Other analyze: anlz_sumtrndseason(), anlz_trans()

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
anlz_trndseason(mod, doystr = 90, doyend = 180, justify = 'center', win = 4)

Raw data from San Francisco Estuary (South Bay)

Description

Raw data from San Francisco Estuary (South Bay)

Usage

rawdat

Format

A data.frame object with 12411 rows and 8 columns

date

Date

station

int

param

chr

value

num

doy

num

cont_year

num

yr

num

mo

Ord.factor

Details

Data from datprc object in https://github.com/fawda123/SFbaytrends

Plot period (seasonal) averages from fitted GAM

Description

Plot period (seasonal) averages from fitted GAM

Usage

show_metseason(
  mod,
  metfun = mean,
  doystr = 1,
  doyend = 364,
  yrstr = 2000,
  yrend = 2019,
  yromit = NULL,
  ylab,
  width = 0.9,
  size = 1.5,
  seascol = "deepskyblue3",
  trndcol = "pink",
  nsim = 10000,
  useave = FALSE,
  base_size = 11,
  xlim = NULL,
  ylim = NULL,
  ...
)

Arguments

Details

Setting yrstr or yrend to NULL will suppress plotting of the trend line for the meta-analysis regression model.

The optional yromit vector can be used to omit years from the plot and trend assessment. This may be preferred if seasonal estimates for a given year have very wide confidence intervals likely due to limited data, which can skew the trend assessments.

Set useave = T to speed up calculations if metfun = mean. This will use anlz_avgseason to estimate the seasonal summary metrics using a non-stochastic equation.

Value

A ggplot object

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)
  
mod <- anlz_gam(tomod, trans = 'ident')

show_metseason(mod, doystr = 90, doyend = 180, yrstr = 2016, yrend = 2019, 
     ylab = 'Chlorophyll-a (ug/L)')


# show seasonal metrics without annual trend
show_metseason(mod, doystr = 90, doyend = 180, yrstr = NULL, yrend = NULL, 
     ylab = 'Chlorophyll-a (ug/L)')
     
# omit years from the analysis
show_metseason(mod, doystr = 90, doyend = 180, yrstr = 2016, yrend = 2019,
     yromit = 2017, ylab = 'Chlorophyll-a (ug/L)')
      

Plot seasonal metrics and rates of change

Description

Plot seasonal metrics and rates of change

Usage

show_mettrndseason(
  mod,
  metfun = mean,
  doystr = 1,
  doyend = 364,
  justify = c("center", "left", "right"),
  win = 5,
  nsim = 10000,
  useave = FALSE,
  yromit = NULL,
  ylab,
  width = 0.9,
  size = 3,
  nms = NULL,
  fils = NULL,
  cmbn = F,
  base_size = 11,
  xlim = NULL,
  ylim = NULL,
  ...
)

Arguments

Details

The plot is the same as that returned by show_metseason with the addition of points for the seasonal metrics colored by the trends estimated from anlz_trndseason for the specified window and justification.

Four colors are used to define increasing, decreasing, no trend, or no estimate (i.e., too few points for the window). The names and the colors can be changed using the nms and fils arguments, respectively. The cmbn argument can be used to combine the no trend and no estimate colors into one color and label. Although this may be desired for aesthetic reasons, the colors and labels may be misleading with the default names since no trend is shown for points where no estimates were made.

The optional yromit vector can be used to omit years from the plot and trend assessment. This may be preferred if seasonal estimates for a given year have very wide confidence intervals likely due to limited data, which can skew the trend assessments.

Value

A ggplot object

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
show_mettrndseason(mod, metfun = mean, doystr = 90, doyend = 180, justify = 'center', 
  win = 4, ylab = 'Chlorophyll-a (ug/L)')

Plot percent change trends from GAM results for selected time periods

Description

Plot percent change trends from GAM results for selected time periods

Usage

show_perchg(
  mod,
  baseyr,
  testyr,
  ylab,
  base_size = 11,
  xlim = NULL,
  ylim = NULL
)

Arguments

Value

A ggplot object

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl')
  
mod <- anlz_gam(tomod, trans = 'log10')

show_perchg(mod, baseyr = 1995, testyr = 2016, ylab = 'Chlorophyll-a (ug/L)')

Plot a 3-d surface of predictions

Description

Plot a 3-d surface of predictions

Usage

show_prd3d(mod, ylab)

Arguments

Value

a plotly surface

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)
  
mod <- anlz_gam(tomod, trans = 'log10')

show_prd3d(mod, ylab = 'Chlorophyll-a (ug/L)')

Plot predictions for GAMs against day of year

Description

Plot predictions for GAMs against day of year

Usage

show_prddoy(mod, ylab, yromit = NULL, size = 0.5, alpha = 1, base_size = 11)

Arguments

Details

The optional yromit vector can be used to omit years from the plot. This may be preferred if the predicted values from the model deviate substantially from other years likely due to missing data.

Value

A ggplot object

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl')

mod <- anlz_gam(tomod, trans = 'log10')

show_prddoy(mod, ylab = 'Chlorophyll-a (ug/L)')

Plot predictions for GAMs over time, by season

Description

Plot predictions for GAMs over time, by season

Usage

show_prdseason(
  mod,
  ylab,
  yromit = NULL,
  base_size = 11,
  xlim = NULL,
  ylim = NULL
)

Arguments

Value

A ggplot object

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)
  
mod <- anlz_gam(tomod, trans = 'log10')
show_prdseason(mod, ylab = 'Chlorophyll-a (ug/L)')

Plot predictions for GAMs over time series

Description

Plot predictions for GAMs over time series

Usage

show_prdseries(
  mod,
  ylab,
  yromit = NULL,
  alpha = 0.7,
  base_size = 11,
  xlim = NULL,
  ylim = NULL,
  col = "brown"
)

Arguments

Details

The optional yromit vector can be used to omit years from the plot. This may be preferred if the predicted values from the model deviate substantially from other years likely due to missing data.

Value

A ggplot object

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl')

mod <- anlz_gam(tomod, trans = 'log10')

show_prdseries(mod, ylab = 'Chlorophyll-a (ug/L)')

Plot seasonal rates of change based on average estimates for multiple window widths

Description

Plot seasonal rates of change based on average estimates for multiple window widths

Usage

show_sumtrndseason(
  mod,
  doystr = 1,
  doyend = 364,
  yromit = NULL,
  justify = c("center", "left", "right"),
  win = 5:15,
  txtsz = 6,
  cols = c("lightblue", "lightgreen"),
  base_size = 11
)

Arguments

Details

This function plots output from anlz_sumtrndseason.

The optional yromit vector can be used to omit years from the plot and trend assessment. This may be preferred if seasonal estimates for a given year have very wide confidence intervals likely due to limited data, which can skew the trend assessments.

Value

A ggplot2 plot

See Also

Other show: show_sumtrndseason2()

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
show_sumtrndseason(mod, doystr = 90, doyend = 180, justify = 'center', win = 2:3)

Plot seasonal rates of change in quarters based on average estimates for multiple window widths

Description

Plot seasonal rates of change in quarters based on average estimates for multiple window widths

Usage

show_sumtrndseason2(
  mod,
  yromit = NULL,
  justify = c("center", "left", "right"),
  win = 5:15,
  txtsz = 6,
  cols = c("lightblue", "lightgreen"),
  base_size = 11
)

Arguments

Details

This function is similar to show_sumtrndseason but results are grouped into seasonal quarters as four separate plots with a combined color scale.

The optional yromit vector can be used to omit years from the plot and trend assessment. This may be preferred if seasonal estimates for a given year have very wide confidence intervals likely due to limited data, which can skew the trend assessments.

Value

A ggplot2 plot

See Also

Other show: show_sumtrndseason()

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
show_sumtrndseason2(mod, justify = 'center', win = 2:3)

Plot rates of change based on seasonal metrics

Description

Plot rates of change based on seasonal metrics

Usage

show_trndseason(
  mod,
  metfun = mean,
  doystr = 1,
  doyend = 364,
  type = c("log10", "approx"),
  justify = c("left", "right", "center"),
  win = 5,
  ylab,
  nsim = 10000,
  yromit = NULL,
  useave = FALSE,
  base_size = 11,
  nms = NULL,
  fils = NULL,
  cols = NULL,
  xlim = NULL,
  ylim = NULL,
  ...
)

Arguments

Value

A ggplot object

Examples

library(dplyr)

# data to model
tomod <- rawdat %>%
  filter(station %in% 34) %>%
  filter(param %in% 'chl') %>% 
  filter(yr > 2015)

mod <- anlz_gam(tomod, trans = 'log10')
show_trndseason(mod, doystr = 90, doyend = 180, justify = 'left', win = 4,
     ylab = 'Slope Chlorophyll-a (ug/L/yr)')