Setup and background

This document explores how variation in sampling frequency and timing may influence conclusions for trend analyses. We use seasonal metrics of chlorophyll concentration as estimated from mixed-effects meta-analysis analysis of GAM results as described in Beck et al. 2022. Observational chlorophyll data are used for station 30 in South Bay, as shown below in the chldat data object. All source materials are here.

The primary question addressed is:

  • How does variation in sample effort affect conclusions on trends, given that effort varied by month and season for the duration of the time series?

First, some setup material for the analysis:

library(tidyverse)
library(wqtrends)
library(patchwork)
library(lubridate)
library(colorspace)
library(here)
library(knitr)

opts_chunk$set(warning = FALSE, message = FALSE, echo = T)

load(file = here('data/datprc.RData'))

# chl data station 30
chldat <- datprc %>% 
  filter(station == 30) %>% 
  filter(param == 'chl')

# year groups
delim <- 3
delims <- unique(year(chldat$date)) %>% 
  as.numeric %>% 
  range %>% 
  {seq(.[1], .[2], length.out = delim + 1)} %>% 
  round(0)

# add temporal groupings to chldat 
chldat <- chldat %>% 
  mutate(
    doy = yday(date - 31),
    yr = ifelse(mo == 'Jan', yr - 1, yr),
    yrgrp = cut(yr, breaks = delims, include.lowest = T, right = F, dig.lab = 4), 
    qrt = quarter(date, fiscal_start = 2),
    qrt = factor(qrt, levels = c('1', '2', '3', '4'), labels = c('FMA', 'MJJ', 'ASO', 'NDJ'))
  ) %>% 
  filter(yr > 1977)

This is the dataset that is evaluated. Note the addition of different columns for temporal and seasonal groupings. The yr and doy variables are offset by one month, where year starts on February to more closely follow the water year.

head(chldat)
##         date station param value doy cont_year   yr  mo       yrgrp qrt
## 1 1978-02-16      30   chl  2.20  16  1978.126 1978 Feb [1978,1992) FMA
## 2 1978-03-14      30   chl  4.65  42  1978.197 1978 Mar [1978,1992) FMA
## 3 1978-03-20      30   chl  9.70  48  1978.214 1978 Mar [1978,1992) FMA
## 4 1978-04-05      30   chl  2.90  64  1978.258 1978 Apr [1978,1992) FMA
## 5 1978-05-24      30   chl  3.90 113  1978.392 1978 May [1978,1992) MJJ
## 6 1978-05-25      30   chl  2.10 114  1978.395 1978 May [1978,1992) MJJ

Down-sampling methods

Here, it is shown how effort varied by month and seasonal groupings at station 30 for the approximate forty year period of record.

smpeff_plo <- function(datin){

  toplo <- datin 
  
  p1 <- ggplot(toplo, aes(x = mo, group = yrgrp, fill = yrgrp)) + 
    geom_bar(position = position_dodge(), color = 'black', alpha = 0.7) + 
    theme_minimal() + 
    labs(
      fill = 'Year group', 
      x = NULL ,
      y = 'Observations'
    )
  
  p2 <- ggplot(toplo, aes(x = qrt, group = yrgrp, fill = yrgrp)) + 
    geom_bar(position = position_dodge(), color = 'black', alpha = 0.7) + 
    theme_minimal() + 
    labs(
      fill = 'Year group', 
      x = NULL ,
      y = 'Observations'
    )
  
  p <- p1 + p2 + plot_layout(ncol = 1, guides = 'collect')
  
  return(p)
  
}

smpeff_plo(chldat)

It is clearly shown that sampling in the spring occurred more frequently earlier in the period of record, whereas more sampling in the fall occurred later in the period of record. Regardless of time period, more sampling occurs in the spring.

To evaluate an effect on trend conclusions, the observed time series can be modeled with a GAM and seasonal conclusions can be evaluated for different time periods. This can then be compared with different “down-sampled” versions of the original time series to test different scenarios. The two tested scenarios include:

  1. Equal sampling between months
  2. Equal sampling between year groups

The down-sampling for both scenarios is demonstrated below.

Scenario 1: equal sampling between months

The minimum number of samples in each month, year group combination is identified across all groups. The down-sampled dataset takes a random sample for each month, year group combination that is equal in size to the minimum criteria. A function is made to use later. Note the use of a random seed, which becomes important later.

set.seed(123)
downsamp_fun1 <- function(datin){
  
  tosmp <- datin %>% 
    group_by(mo, yrgrp) %>% 
    summarise(
      cnt = n(), 
      .groups = 'drop'
    ) %>% 
    group_by(mo) %>% 
    filter(cnt == min(cnt)) %>% 
    ungroup() %>% 
    select(-yrgrp)
  
  minmo <- min(tosmp$cnt)
  
  out <- datin %>%
    group_by(yrgrp, mo) %>% 
    sample_n(size = minmo) %>% 
    ungroup()
  
  return(out)  
  
}

chlsmp1 <- downsamp_fun1(chldat)

Now, evaluating the counts in each temporal category as before shows that the down-sampled time series has a more equal distribution of samples for each month.

smpeff_plo(chlsmp1)

Scenario 2: equal sampling across year groups

This down-sampling scenario is similar as above, except the approximate sampling bias towards spring is more balanced between the year groups. The minimum number of samples for each month across year groups is used as the down-sampling criteria.

downsamp_fun2 <- function(datin){

  tosmp <- datin %>% 
    group_by(mo, yrgrp) %>% 
    summarise(
      cnt = n(), 
      .groups = 'drop'
    ) %>% 
    group_by(mo) %>% 
    filter(cnt == min(cnt)) %>% 
    ungroup() %>% 
    select(-yrgrp)
  
  out <- datin %>%
    left_join(tosmp, by = 'mo') %>% 
    group_by(yrgrp, mo) %>% 
    sample_n(unique(cnt))  
  
  return(out)
  
}

chlsmp2 <- downsamp_fun2(chldat)
smpeff_plo(chlsmp2)

Evaluate effect on maximum peak estimates

The above was repeated for the February to April period using five random down-sampled estimates and extracting the estimate of the maximum chlorophyll peak. This tests the hypothesis that less frequent sampling earlier in the period of record may not capture high chlorophyll events and that should be reflected in lower GAM estimates in the down-sampled data. The down-sampling assumes equal month.

# random sample label
rnds <- paste('Equal month', 1:5)

# get max seasonal estimates from random down-sampled data
ests <- tibble(chlsmp = rnds) %>% 
  group_by(chlsmp) %>% 
  nest() %>% 
  mutate(
    data = list(chldat), 
    data = purrr::map(data, downsamp_fun1), 
    mods = purrr::map(data, anlz_gam, trans = 'log10'),
    vals = purrr::map(mods, anlz_metseason, metfun = max, doystr = 1, doyend = 89, nsim = 1e4, na.rm = T)
  ) %>% 
  select(chlsmp, vals) %>% 
  unnest('vals')

# get max seasonal estimates from full dataset
origests <- anlz_metseason(modorig, metfun = max, doystr = 1, doyend = 89, nsim = 1e4, na.rm = T) %>% 
  mutate(chlsmp = 'Full dataset')

# combine all results
allests <- bind_rows(ests, origests) %>% 
  mutate(
    col = case_when(
      chlsmp == 'Full dataset' ~ 'Full dataset', 
      T ~ 'Down-sampled'
    ), 
    col = factor(col, levels = c('Full dataset', 'Down-sampled')),
    chlsmp = factor(chlsmp, levels = c('Full dataset', rnds))
  )

head(allests)

Now make the plot for comparison.

wd <- 0.5
ggplot(data = allests, aes(x = yr, y = bt_met, color = col, group = chlsmp)) + 
  geom_point(position = position_dodge(width = wd)) +
  geom_errorbar(aes(ymin = bt_lwr, ymax = bt_upr), position = position_dodge(width = wd), width = 0) +
  scale_y_log10() +
  scale_color_manual(values = c('tomato1', 'black')) +
  theme_bw() + 
  theme(
    axis.title.x = element_blank(), 
    legend.position = 'top'
  ) + 
  labs(
    y = 'Chlorophyll-a (ug/L)', 
    color = NULL
  )

The above was repeated using 1000 random samples, then taking the average of the seasonal metrics for each random sample in each year. The average of the random samples was estimated mixed-effects meta-analysis to account for the uncertainty of each seasonal metric for each random sample. The mixed-effects meta-analysis model was used to calculate an intercept only model.

The random samples and their seasonal metrics take about 3 hours to run, but the code is shown here for reproducibility.

rnds <- paste('Equal month', 1:1000)

# get max seasonal estimates from random down-sampled data
str <- Sys.time()
smpests <- tibble(chlsmp = rnds) %>% 
  group_by(chlsmp) %>% 
  nest() %>% 
  mutate(
    data = list(chldat), 
    data = purrr::map(data, downsamp_fun1), 
    mods = purrr::map(data, anlz_gam, trans = 'log10'),
    vals = purrr::map(mods, anlz_metseason, metfun = max, doystr = 1, doyend = 89, nsim = 1e4, na.rm = T)
  ) %>% 
  select(chlsmp, vals) %>% 
  unnest('vals')
print(Sys.time() - str)

save(smpests, file = here('data/smpests.RData'))
load(file = here('data/smpests.RData'))

cmb <- smpests %>%
  group_by(yr) %>%
  nest() %>% 
  mutate(
    mixmet = purrr::map(data, function(x) mixmeta::mixmeta(met ~ 1, S = se, random = ~1|chlsmp, data = x, method = 'reml')), 
    met = purrr::map(mixmet, function(x) summary(x)$coefficients[1]),
    se = purrr::map(mixmet, function(x) summary(x)$coefficients[2])
  ) %>% 
  select(yr, met, se) %>% 
  unnest(cols = c('met', 'se')) %>% 
  mutate(
    bt_lwr = 10^(met - 1.96 * se),
    bt_upr = 10^(met + 1.96 * se),
    bt_met = 10^met,
    chlsmp = '1000 sim avg'
  )
  
# combine all results
allests <- bind_rows(cmb, origests) %>% 
  mutate(chlsmp = factor(chlsmp, levels = c('Full dataset', '1000 sim avg')))

wd <- 0.5
ggplot(data = allests, aes(x = yr, y = bt_met, color = chlsmp, group = chlsmp)) + 
  geom_point(position = position_dodge(width = wd)) +
  geom_errorbar(aes(ymin = bt_lwr, ymax = bt_upr), position = position_dodge(width = wd), width = 0) +
  scale_y_log10() +
  scale_color_manual(values = c('tomato1', 'black')) +
  theme_bw() + 
  theme(
    axis.title.x = element_blank(), 
    legend.position = 'top'
  ) + 
  labs(
    y = 'Chlorophyll-a (ug/L)', 
    color = NULL
  )