correct use of yromit to remove years from plots and trend analysis, …

…propagates to all functions that use seasonal metrics

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: wqtrends
 Title: Assess Water Quality Trends with Generalized Additive Models
 Version: 1.4.2
-Date: 2024-01-26
+Date: 2024-03-03
 Authors@R: c(
     person(given = "Marcus",
       family = "Beck",

R/anlz_avgseason.R

@@ -5,6 +5,7 @@
 #' @param mod input model object as returned by \code{\link{anlz_gam}}
 #' @param doystr numeric indicating start Julian day for extracting averages
 #' @param doyend numeric indicating ending Julian day for extracting averages
+#' @param yromit optional numeric vector for years to omit from the output
 #'
 #' @return A data frame of period averages
 #' @export
@@ -22,7 +23,7 @@
 #'
 #' mod <- anlz_gam(tomod, trans = 'log10')
 #' anlz_avgseason(mod, doystr = 90, doyend = 180)
-anlz_avgseason <- function(mod, doystr = 1, doyend = 364) {
+anlz_avgseason <- function(mod, doystr = 1, doyend = 364, yromit = NULL) {
 
   # transformation
   trans <- mod$trans
@@ -63,6 +64,11 @@ anlz_avgseason <- function(mod, doystr = 1, doyend = 364) {
     tibble::tibble() %>% 
     dplyr::mutate(dispersion = dispersion) %>% 
     dplyr::select(yr, met, se, bt_lwr, bt_upr, bt_met, dispersion)
+
+  # fill NA for yromit (do not remove)
+  if(!is.null(yromit))
+    out <- out %>% 
+      dplyr::mutate_at(dplyr::vars(-yr), ~ifelse(yr %in% yromit, NA, .)) 
 
   return(out)
 

R/anlz_metseason.R

@@ -7,6 +7,7 @@
 #' @param doystr numeric indicating start Julian day for extracting averages
 #' @param doyend numeric indicating ending Julian day for extracting averages
 #' @param nsim numeric indicating number of random draws for simulating uncertainty
+#' @param yromit optional numeric vector for years to omit from the output
 #' @param ... additional arguments passed to \code{metfun}, e.g., \code{na.rm = TRUE)}
 #'
 #' @return A data frame of period metrics
@@ -27,7 +28,7 @@
 #'
 #' mod <- anlz_gam(tomod, trans = 'log10')
 #' anlz_metseason(mod, mean, doystr = 90, doyend = 180, nsim = 100)
-anlz_metseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, nsim = 1e4, ...) {
+anlz_metseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, nsim = 1e4, yromit = NULL, ...) {
 
   # transformation
   trans <- mod$trans
@@ -88,7 +89,12 @@ anlz_metseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, nsim =
   mets$dispersion <- dispersion
 
   out <- mets
-
+
+  # fill NA for yromit (do not remove)
+  if(!is.null(yromit))
+    out <- out %>% 
+      dplyr::mutate_at(dplyr::vars(-yr), ~ifelse(yr %in% yromit, NA, .)) 
+
   return(out)
 
 }

R/anlz_mixmeta.R

@@ -5,7 +5,6 @@
 #' @param metseason output from \code{\link{anlz_metseason}}
 #' @param yrstr numeric for starting year
 #' @param yrend numeric for ending year
-#' @param yromit optional numeric vector for years to omit from, inherited from \code{\link{show_metseason}}
 #' @details Parameters are not back-transformed if the original GAM used a transformation of the response variable
 #' 
 #' @concept analyze
@@ -25,17 +24,19 @@
 #' mod <- anlz_gam(tomod, trans = 'log10')
 #' metseason <- anlz_metseason(mod, doystr = 90, doyend = 180)
 #' anlz_mixmeta(metseason, yrstr = 2016, yrend = 2019)
-anlz_mixmeta <- function(metseason, yrstr = 2000, yrend = 2019, yromit = NULL){
+anlz_mixmeta <- function(metseason, yrstr = 2000, yrend = 2019){
 
   # input
   totrnd <- metseason %>% 
     dplyr::mutate(S = se^2) %>% 
     dplyr::filter(yr %in% seq(yrstr, yrend))
 
-  if(nrow(totrnd) != length(seq(yrstr, yrend)) & is.null(yromit))
+  if(nrow(totrnd) != length(seq(yrstr, yrend)))
     return(NA)
 
-  out <- mixmeta::mixmeta(met ~ yr, S = S, random = ~1|yr, data = totrnd, method = 'reml')
+  out <- try(mixmeta::mixmeta(met ~ yr, S = S, random = ~1|yr, data = totrnd, method = 'reml'), silent = TRUE)
+  if(inherits(out, 'try-error'))
+    return(NA)
 
   return(out)
 

R/anlz_sumtrndseason.R

@@ -7,10 +7,15 @@
 #' @param doyend numeric indicating ending Julian day for extracting averages
 #' @param justify chr string indicating the justification for the trend window
 #' @param win numeric vector indicating number of years to use for the trend window
+#' @param yromit optional numeric vector for years to omit from the plot, see details
 #'
 #' @return A data frame of slope estimates and p-values for each year
 #' @export
 #'
+#' @details
+#'  
+#' The optional \code{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.
+#' 
 #' @concept analyze
 #'
 #' @details This function is a wrapper to \code{\link{anlz_trndseason}} to loop across values in \code{win}, using \code{useave = TRUE} for quicker calculation of average seasonal metrics.  It does not work with any other seasonal metric calculations.
@@ -28,13 +33,14 @@
 #'
 #' mod <- anlz_gam(tomod, trans = 'log10')
 #' anlz_sumtrndseason(mod, doystr = 90, doyend = 180, justify = 'center', win = 2:3)
-anlz_sumtrndseason <- function(mod, doystr = 1, doyend = 364, justify = c('center', 'left', 'right'), win = 5:15){
+anlz_sumtrndseason <- function(mod, doystr = 1, doyend = 364, justify = c('center', 'left', 'right'), win = 5:15, 
+                               yromit = NULL){
 
   justify <- match.arg(justify)
 
   tmp <- NULL
   for(i in win){
-    res <- anlz_trndseason(mod, metfun = mean, doystr = doystr, doyend = doyend, justify = justify, win = i, useave = T)
+    res <- anlz_trndseason(mod, metfun = mean, doystr = doystr, doyend = doyend, justify = justify, win = i, useave = T, yromit = yromit)
     res$win <- i
     tmp <- rbind(tmp, res)
   }

R/anlz_trndseason.R

@@ -9,6 +9,7 @@
 #' @param justify chr string indicating the justification for the trend window
 #' @param win numeric indicating number of years to use for the trend window, see details
 #' @param nsim numeric indicating number of random draws for simulating uncertainty
+#' @param yromit optional numeric vector for years to omit from the output
 #' @param useave logical indicating if \code{anlz_avgseason} is used for the seasonal metric calculation
 #' @param ... additional arguments passed to \code{metfun}, e.g., \code{na.rm = TRUE)}
 #'
@@ -20,6 +21,8 @@
 #' @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 \code{\link{anlz_mixmeta}}) for the seasonal metrics (from \code{\link{anlz_metseason}}).
 #' 
 #' Note that for left and right windows, the exact number of years in \code{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., for 1990 it 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 \code{win = 8} and \code{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 \code{win = 7} and \code{justify = 'center'}, the estimate for 2000 will be centered on 1997 and 2003 (three years left, three years right, seven years total).
+#' 
+#' The optional \code{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.
 #'
 #' @family analyze
 #' 
@@ -34,7 +37,7 @@
 #'
 #' mod <- anlz_gam(tomod, trans = 'log10')
 #' anlz_trndseason(mod, doystr = 90, doyend = 180, justify = 'center', win = 4)
-anlz_trndseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, justify = c('center', 'left', 'right'), win = 5, nsim = 1e4, useave = FALSE, ...){
+anlz_trndseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, justify = c('center', 'left', 'right'), win = 5, nsim = 1e4, yromit = NULL, useave = FALSE, ...){
 
   justify <- match.arg(justify)
 
@@ -44,36 +47,39 @@ anlz_trndseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, justif
   # make sure user wants average and useave
   if(!chk & useave)
     stop('Specify metfun = mean if useave = T')
-  
+
   # estimate metrics
   if(useave)
-    metseason <- anlz_avgseason(mod, doystr = doystr, doyend = doyend)
+    metseason <- anlz_avgseason(mod, doystr = doystr, doyend = doyend, yromit = yromit)
   if(!useave)
-    metseason <- anlz_metseason(mod, metfun, doystr = doystr, doyend = doyend, nsim = nsim, ...)
+    metseason <- anlz_metseason(mod, metfun, doystr = doystr, doyend = doyend, nsim = nsim, yromit = yromit,...)
 
   tmp <- tibble::tibble(metseason)
   tmp$yrcoef <- NaN
   tmp$pval <- NaN
-  
+
   # iterate through years to get trend
   for(i in seq_along(tmp$yr)){
 
     yr <- tmp$yr[i]
 
-    if(justify == 'left')
-      mixmet <- anlz_mixmeta(tmp, yrstr = yr, yrend = yr + win - 1)
+    if(justify == 'left'){
+      yrstr <- yr
+      yrend <- yr + win - 1
+
+    }
+
+    if(justify == 'right'){
+      yrstr <- yr - win + 1
+      yrend <- yr
+    }
 
-    if(justify == 'right')
-      mixmet <- anlz_mixmeta(tmp, yrstr = yr - win + 1, yrend = yr)
-
     if(justify == 'center'){
-
       yrstr <- floor(yr - win / 2) + 1
       yrend <- floor(yr + win / 2)
-
-      mixmet <- anlz_mixmeta(tmp, yrstr = yrstr, yrend = yrend)
-
     }
+
+    mixmet <- anlz_mixmeta(tmp, yrstr = yrstr, yrend = yrend)
 
     if(inherits(mixmet, 'logical'))
       next

R/show_metseason.R

@@ -26,7 +26,7 @@
 #' 
 #' Setting \code{yrstr} or \code{yrend} to \code{NULL} will suppress plotting of the trend line for the meta-analysis regression model.
 #' 
-#' The optional \code{omityr} 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.
+#' The optional \code{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 \code{useave = T} to speed up calculations if \code{metfun = mean}.  This will use \code{\link{anlz_avgseason}} to estimate the seasonal summary metrics using a non-stochastic equation.
 #' 
@@ -52,8 +52,8 @@
 #'      ylab = 'Chlorophyll-a (ug/L)')
 #'      
 #' # omit years from the analysis
-#' show_metseason(mod, doystr = 90, doyend = 180, yrstr = 2015, yrend = 2019,
-#'      yromit = 2018, ylab = 'Chlorophyll-a (ug/L)')
+#' show_metseason(mod, doystr = 90, doyend = 180, yrstr = 2016, yrend = 2019,
+#'      yromit = 2017, ylab = 'Chlorophyll-a (ug/L)')
 #' }      
 show_metseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, yrstr = 2000, yrend = 2019, yromit = NULL, ylab, width = 0.9, size = 1.5, nsim = 1e4, useave = FALSE, base_size = 11, xlim = NULL, ylim = NULL, ...) {
 
@@ -63,17 +63,12 @@ show_metseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, yrstr =
   # make sure user wants average and useave
   if(!chk & useave)
     stop('Specify metfun = mean if useave = T')
-  
+
   # estimate metrics
   if(useave)
-    metseason <- anlz_avgseason(mod, doystr = doystr, doyend = doyend)
+    metseason <- anlz_avgseason(mod, doystr = doystr, doyend = doyend, yromit = yromit)
   if(!useave)
-    metseason <- anlz_metseason(mod, metfun, doystr = doystr, doyend = doyend, nsim = nsim, ...)
-
-  # omit years if yromit provided
-  if(!is.null(yromit))
-    metseason <- metseason %>% 
-      dplyr::filter(!yr %in% yromit)
+    metseason <- anlz_metseason(mod, metfun, doystr = doystr, doyend = doyend, nsim = nsim, yromit = yromit, ...)
 
   # transformation used
   trans <- mod$trans
@@ -93,8 +88,8 @@ show_metseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, yrstr =
 
   # plot output
   p <- ggplot2::ggplot(data = toplo1, ggplot2::aes(x = yr, y = bt_met)) + 
-    ggplot2::geom_point(colour = 'deepskyblue3', size = size) +
-    ggplot2::geom_errorbar(ggplot2::aes(ymin = bt_lwr, ymax = bt_upr), colour = 'deepskyblue3', width = width) +
+    ggplot2::geom_point(colour = 'deepskyblue3', size = size, na.rm = TRUE) +
+    ggplot2::geom_errorbar(ggplot2::aes(ymin = bt_lwr, ymax = bt_upr), colour = 'deepskyblue3', width = width, na.rm = TRUE) +
     ggplot2::theme_bw(base_size = base_size) + 
     ggplot2::theme(
       axis.title.x = ggplot2::element_blank()
@@ -104,8 +99,8 @@ show_metseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, yrstr =
   if(!is.null(yrstr) & !is.null(yrend)){
 
     # get mixmeta models
-    mixmet <- anlz_mixmeta(metseason, yrstr = yrstr, yrend = yrend, yromit = yromit)
-  
+    mixmet <- anlz_mixmeta(metseason, yrstr = yrstr, yrend = yrend)
+
     toplo2 <- data.frame(
       yr = seq(yrstr, yrend, length = 50)
       ) %>% 

R/show_mettrndseason.R

@@ -19,6 +19,8 @@
 #' 
 #' 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 \code{nms} and \code{cols} arguments, respectively.  The \code{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 \code{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.
+#' 
 #' @concept show
 #' 
 #' @return A \code{\link[ggplot2]{ggplot}} object
@@ -39,8 +41,8 @@
 show_mettrndseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, justify = c('center', 'left', 'right'), win = 5, nsim = 1e4, useave = FALSE, yromit = NULL, ylab, width = 0.9, size = 3, nms = NULL, cols = NULL, cmbn = F, base_size = 11, xlim = NULL, ylim = NULL, ...){
 
   # get seasonal metrics and trends
-  trndseason <- anlz_trndseason(mod = mod, metfun, doystr = doystr, doyend = doyend, justify = justify, win = win, useave = useave)
-  
+  trndseason <- anlz_trndseason(mod = mod, metfun, doystr = doystr, doyend = doyend, justify = justify, win = win, useave = useave, yromit = yromit)
+
   # handle nms and cols args if not combine (keep no trend and no estimate)
   if(!cmbn){
 
@@ -88,11 +90,6 @@ show_mettrndseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, jus
 
   names(cols) <- nms
 
-  # omit years if yromit provided
-  if(!is.null(yromit))
-    trndseason <- trndseason %>% 
-      dplyr::filter(!yr %in% yromit)
-
   # title
   dts <- as.Date(c(doystr, doyend), origin = as.Date("2000-12-31"))
   strt <- paste(lubridate::month(dts[1], label = T, abbr = T), lubridate::day(dts[1]))
@@ -104,11 +101,11 @@ show_mettrndseason <- function(mod, metfun = mean, doystr = 1, doyend = 364, jus
   subttl <- paste0('Points colored by trend for ', win, '-year, ', justify, '-justified window')
 
   toplo <- trndseason
-  
+
   # plot output
   p <- ggplot2::ggplot(data = toplo, ggplot2::aes(x = yr, y = bt_met)) + 
-    ggplot2::geom_errorbar(ggplot2::aes(ymin = bt_lwr, ymax = bt_upr), colour = 'black', width = width) +
-    ggplot2::geom_point(ggplot2::aes(fill = trnd), pch = 21, color = 'black', size = size) +
+    ggplot2::geom_errorbar(ggplot2::aes(ymin = bt_lwr, ymax = bt_upr), colour = 'black', width = width, na.rm = TRUE) +
+    ggplot2::geom_point(ggplot2::aes(fill = trnd), pch = 21, color = 'black', size = size, na.rm = TRUE) +
     ggplot2::theme_bw(base_size = base_size) +
     ggplot2::scale_fill_manual(values = cols, drop = F) +
     ggplot2::theme(

R/show_sumtrndseason.R

@@ -5,6 +5,7 @@
 #' @param mod input model object as returned by \code{\link{anlz_gam}}
 #' @param doystr numeric indicating start Julian day for extracting averages
 #' @param doyend numeric indicating ending Julian day for extracting averages
+#' @param yromit optional numeric vector for years to omit from the plot, see details
 #' @param justify chr string indicating the justification for the trend window
 #' @param win numeric vector indicating number of years to use for the trend window
 #' @param txtsz numeric for size of text labels inside the plot
@@ -18,6 +19,8 @@
 #'
 #' @details This function plots output from \code{\link{anlz_sumtrndseason}}.
 #' 
+#' The optional \code{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.
+#' 
 #' @family show
 #' 
 #' @examples
@@ -31,7 +34,8 @@
 #'
 #' mod <- anlz_gam(tomod, trans = 'log10')
 #' show_sumtrndseason(mod, doystr = 90, doyend = 180, justify = 'center', win = 2:3)
-show_sumtrndseason <- function(mod, doystr = 1, doyend = 364, justify = c('center', 'left', 'right'), 
+show_sumtrndseason <- function(mod, doystr = 1, doyend = 364, yromit = NULL, 
+                              justify = c('center', 'left', 'right'), 
                               win = 5:15, txtsz = 6, cols = c('lightblue', 'lightgreen'), 
                               base_size = 11){
 
@@ -41,7 +45,7 @@ show_sumtrndseason <- function(mod, doystr = 1, doyend = 364, justify = c('cente
   sig_vals <- c(-Inf, 0.005, 0.05, Inf)
 
   # get ests across all window widths
-  res <- anlz_sumtrndseason(mod, doystr = doystr, doyend = doyend, justify = justify, win = win)
+  res <- anlz_sumtrndseason(mod, doystr = doystr, doyend = doyend, justify = justify, win = win, yromit = yromit)
 
   # seasonal range for title
   dts <- as.Date(c(doystr, doyend), origin = as.Date("2000-12-31"))
@@ -66,8 +70,8 @@ show_sumtrndseason <- function(mod, doystr = 1, doyend = 364, justify = c('cente
     )
 
   p <- ggplot2::ggplot(toplo, ggplot2::aes(x = yr, y = win, fill = yrcoef)) +
-    ggplot2::geom_tile(color = 'black') + 
-    ggplot2::geom_text(ggplot2::aes(label = psig), size = txtsz) + 
+    ggplot2::geom_tile(color = 'black', na.rm = TRUE) + 
+    ggplot2::geom_text(ggplot2::aes(label = psig), size = txtsz, na.rm = TRUE) + 
     ggplot2::scale_x_continuous(expand = c(0, 0)) + 
     ggplot2::scale_y_continuous(expand = c(0, 0), breaks = win) +
     ggplot2::scale_fill_gradient2(low = cols[1], mid = 'white', high = cols[2], midpoint = 0) +