R/trend.R

In martindurocher/floodRFA: A R-package for regional frequency of floods

###################################################
#' Mann-Kendall test for trend using block bootstrap
#'
#' Return the results of the Mann-Kendall test for trends. In case
#' of autocorrelation, the p-value of the test is obtained by
#' block bootstrap. This function is a wrapper for the function
#' \code{\link{MannKendall}}.
#' If not specified by \code{block}, an automatic selection of
#' the block sizes is provided according to the largest significant
#' autocorrelation lag.
#'
#' @param x Sample.
#'
#' @param nsim Number of bootstrap sample.
#'
#' @param block Size of the bootstrap blocks.
#'   If \code{block} is \code{NULL} it is automatically selected
#'
#' @param block.max Maximal size for for \code{block} if selected automatically.
#'
#' @param acf.alpha Determine the p-value of a significant autocorrelation.
#'
#' @param acf.tol Value added to the largest significant autocorrelation to
#'   determine the block size.
#'
#' @param ... Other parameters pass to \code{\link{tsboot}}.
#'
#' @section References:
#'
#'
#' @export
#'
#' @examples
#'
#'
#' set.seed(1)
#'
#' ## Case of independant data.
#' x <- rlnorm(100, 5, .1)
#' MKendall(x, block = 0)
#'
#' ## Case of autocorrelated data without trend.
#' x<- rep(0,100)
#' for(ii in seq_along(x)[-1])
#'   x[ii] <- 0.5 * x[ii-1] + rnorm(1,0,sd=15)
#'
#' MKendall(x, block = 2)
#' MKendall(x)
#'
#' ## Case with autocorrelation with trend
#' MKendall(x+seq_along(x)/3)
#'
#'  ## Case of trend in peaks magnitude
#'  pid <- which.floodPeaks(flow~date, canadaFlood$daily, u= 1000, r = 14)
#'  MKendall(canadaFlood$daily$flow[pid])


MKendall <- function(x, block = NULL, block.max = 5, nboot = 1000,
                     acf.ci = 0.95, acf.tol = 1, season = FALSE,
                     ...){

  ## If a seasonal test is performed, the data must be a time series
  if(season & class(x) != 'ts')
    stop()

  ## Test under serial independence
  stat <- Kendall::MannKendall(x)

  ## Automatic selection of the bootstrap block size
  if(is.null(block)){

    ## compute the ranks
    z <- rank(x, na.last = 'keep')/(length(x)+1)

    ## compute the autocorrelation
    a  <- acf(z, na.action = na.pass, plot = FALSE)

    ## found a bound for the autocorrelation
    b <- qnorm((1 + acf.ci)/2)/sqrt(a$n.used)

    ## Use the highest (+acf.tol) lag as block size
    a0 <- a$acf[seq(2,block.max+1)]
    if(all(a0 <= b)){
      block <- 0

    } else {
      block <- max(which(a0>b))
      block <- min(block+acf.tol, block.max)
    }

  }

  if(block > 0){
    ## Perform block bootstrap if necessary
    boot <- boot::tsboot(x, R = nboot, l = block, sim = "fixed",
              statistic = function(z) Kendall::MannKendall(z)$tau, ...)$t

    ## Compute the p-value
    pv <- mean(abs(boot) > abs(stat$tau))

  }else{
    pv <- stat$sl
  }

  ans <- list(stat = stat$tau, pvalue = pv, block = block)

  class(ans) <- c('mkendall')
  return(ans)
}

#' @export
print.mkendall <- function(obj){

  if(obj$block>0){
    cat('\nMann-Kendall Test for trend using block bootstrap\n')
    cat('\nBlock size:', obj$block)
  }else{
    cat('\n\nMann-Kendall Test for trend')
  }

  cat('\nTau: ', format(obj$stat, digit = 4))
  cat('\np-value: ', format(obj$pvalue, digit = 4))

}