R/trend.R

In martindurocher/floodStat: FloodNet Regional Frequency Analysis

###################################################
#' Mann-Kendall test for trend using block bootstrap
#'
#' Return 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.
#'
#'
#' @section References:
#'
#' Önöz, B., Bayazit, M. (2012). Block bootstrap for Mann–Kendall trend
#'   test of serially dependent data. Hydrol. Process. 26, 3552–3560.
#'   https://doi.org/10.1002/hyp.8438
#'
#' @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)
#'
#' acf(x)
#' MKendall(x, block = 2)
#' MKendall(x)
#'
#' ## Case with autocorrelation with trend
#' MKendall(x+seq_along(x)/3)


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

  if(!all(is.finite(x)))
    stop('X Contains non finite values')


  ## Test under serial independence
  stat <- MkStat(x)

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

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

    ## 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
    xboot <- replicate(nboot, SampleBlock(x, block))
    sboot <- apply(xboot,2,MkStat, id = combn(length(x),2), pval = FALSE)

    ## Compute the p-value
    pv <- mean(abs(sboot) > abs(stat[1]))

  }else{
    pv <- stat[3]
  }

  ans <- list(stat = stat[1], 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('\nS: ', obj$stat)
  cat('\np-value: ', format(obj$pvalue, digit = 4))

}

## Function that calculate the Statistis
MkStat <- function(x, id = NULL, pval = TRUE){

  n <- length(x)

  if(is.null(id))
   id <- combn(n,2)

  ## Compute the test statistic
  s <- sum(sign(x[id[2,]] -x[id[1,]]))

  ## If necessary compute the pvalue
  if(pval){

    ## Correction for ties
    tj <- table(x)
    tj <- tj[tj>1]
    tj <- sum(tj*(tj-1)*(2*tj+5))

    ## Variance of the statistic
    v <- (n*(n-1)*(2*n+5) - tj )/18

    ## Normal approximation
    if(s>0){
      z <- abs(s-1)/sqrt(v)
    } else if(s < 0){
      z <- abs(s+1)/sqrt(v)
    }else{
      z <- 0
    }

    ## Return p-value based on normal approximation
    ans <- c(s,v,2-2*pnorm(z))

  } else{
    ans <- s
  }

  return(ans)
}

SampleBlock <- function(x, k){
  n <- length(x)
  m <- ceiling(n/k)
  id <- sample(seq(n-k), m,replace = TRUE)
  id <- sapply(id, function(ii) seq(ii,ii+k-1))
  id <- as.numeric(id)[seq(n)]

  return(x[id])
}