R/lnres_err.R

In fawda123/WRTDStidal: Weighted Regression for Water Quality Evaluation in Tidal Waters

#' Simulate random errors from a time series
#' 
#' Simulate random errors of a water quality time series by modelling the statistical properties of an observed dataset
#'
#' @param dat_in input \code{\link[base]{data.frame}} that must include discharge and decimal time columns, see example dataset \code{\link{daydat}}
#' @param yr numeric year value to use for the stationary model, defaults to the median year
#' @param comps logical indicating if the WRTDS model used to get response error measures is also returned, see value.
#' @param seed optional numeric value for random generation seed
#' 
#' @details Random errors for a stationary seasonal water quality time series on a daily time step are generated by modelling residuals from an observed dataset.  First, a stationary seasonal model is created by fitting a \code{\link{wrtds}} model and estimating an error distribution the residuals using the \code{\link[forecast]{auto.arima}} function.  Accumulated standard errors from the regression are also retained for each residual.  Random errors using the estimated auto-regressive structures are simulated using \code{\link[stats]{arima.sim}} for the entire year and multiplied by the corresponding standard error estimate from the regression.  The entire year is then repeated for every year in the observed time series.  The final simulated errors are rescaled to the range of the original residuals that were used to estimate the distribution.   
#' 
#' @return The original data frame with additional columns for the random errors (\code{errs}) and the standard error estimates for each residual (\code{scls}).  If \code{comps = TRUE}, a two-element list is returned that also includes the WRTDS model used as a basis for errors and scale values.
#' 
#' @export
#' 
#' @import dplyr forecast
#' 
#' @seealso \code{\link{daydat}}
#' 
#' @examples 
#' \dontrun{
#' ## example data
#' data(daydat)
#' 
#' ## get errors
#' lnres_err(daydat)
#'}
lnres_err <- function(dat_in, yr = NULL, comps = FALSE, seed = NULL) {

  # pick yr for stationary mod, defaults to median
  if(is.null(yr))
    yr <- median(unique(dat_in$year)) %>% 
      round(., 0)

  # set seed
  set.seed(seed)
      
  # data for the year to get stationary mod
  stat_dat <- filter(dat_in, year == yr)
  
  # use wrtds for the year
  tomod <- select(stat_dat, date, lnres, lnQ) %>% 
    rename(
      res = lnres,
      flo = lnQ
    ) %>% 
    mutate(lim = -1e6)
  cat('\nEstimating stationary, seasonal model with WRTDS...\n')
  resmod <- modfit(tomod, resp_type = 'mean')
  resmod <- resscls(resmod)

  # get model residuals, scale parameter, and decimal time minus year
  stat_toerr <- with(resmod,
    data.frame(
      dec_mo = stat_dat$dec_time - stat_dat$year,
      resids = res - fits,
      scls = scls
      )
    )

  # get arma model from resids
  errmod <- auto.arima(stat_toerr$resids, d = 0, seasonal = FALSE)
  ars <- coef(errmod)[grep('^ar', names(coef(errmod)))]
  mas <- coef(errmod)[grep('^ma', names(coef(errmod)))]
  
  # simulate rnorm errors using arma(p,q) process 
  errs <- arima.sim(list(ar = ars, ma = mas, order = c(length(ars), 0, length(mas))), n = nrow(dat_in), 
    rand.gen = function(x) rnorm(x, 0, 1))
  
  # add errors and scl values to dat_in
  dat_in$errs <- as.numeric(errs)
  dat_in$dec_mo <- with(dat_in, dec_time - year)
  dat_in <- left_join(dat_in, stat_toerr, by = 'dec_mo')
  
  # linear transform errrs to match those in the from the observed data
  rng <- range(dat_in$res, na.rm = T)
  rngerr <- range(errs, na.rm = TRUE)
  errs <- (errs - rngerr[1])/diff(rngerr) * diff(rng) + rng[1]
  
  # remove extra cols, sort on date
  tosel <- c('date', 'flo', 'lnres', 'Q', 'lnQ', 'jday', 'year', 'day', 'dec_time', 'scls', 'errs', 'lnQ_sim')
  dat_in <- dat_in[, names(dat_in) %in% tosel]
  dat_in <- arrange(dat_in, date)
  
  # also return wrtds model used to get scale values and errors for error simulations
  if(comps){
    
    out_ls <- list(
      resmod = resmod, 
      dat_in = dat_in
    )
    
    return(out_ls)
    
  }
      
  return(dat_in)
  
}