R/rfa.R

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

#' Perform regional frequency analysis using L-moments on pooling groups
#'
#' @param obj Output from \link{rfa}. Containt the index-Flood and the Lmoments
#' of all sites
#' @param dis If a vector of distances, the target must be included with distance zero.
#'  If a distance matrix or a distance object, the target must be specified.
#' @param indx A vector or string to force new index-flood value
#' @param rtn A vector of return period
#' @param nk The number of donor
#' @param distr Regioal distribution to use. If omited, a best distribution is selected
#' @param nsim Number of simulation used to compute Homegeneity and GOF statistics
#' @param nk.min  Minimal size of pooling groups
#' @param rm.target Should the target be removed from the
#'
#' @export
#'
#' @examples
#'
#'
#' library(lmomco)
#' coord <- replicate(2, runif(100))
#' dis <- dist(coord)
#'
#' para <- vec2par(c(100, 30, .1),'gev')
#' xdata <- lapply(1:100, function(ii) rlmomco(rpois(1,20)+20, para))
#' names(xdata) <- paste('site',1:100)
#'
#' obj <- rfa(xdata, 'median')
#'
#' rfaRegFit(obj)
#'
#' # Perform the analysis of the site 13th site
#' regGev <- poolGrp(obj, dis = dis, target = 13, distr = 'gev')
#' regGev
#'
#' # Could choose the distribution automatically
#' poolGrp(obj, dis = dis, target = 13)
#'
#' # To get information on the regional growth curve
#' poolGrp(obj, dis = dis, target = 13, indx = 1)
#'
poolGrp <- function(obj, dis, rtn = c(2,5,10,20,50,100),
                    nk = 25, sset = NULL, indx = NULL,
                    distr = NULL, rm.target = FALSE, target = NULL){


  ## select an initial set of donors
  nk <- min(nk,nrow(obj$lmom))

  ## Organise the distance depending if dis is vector or not
  if(is.null(target))
    target <- which(dis==0)
  else
    dis <- as.matrix(dis)[target,]

  ## Should the analysis be restricted to a giving subset
  if(!is.null(sset)){

    if(is.logical(sset)) sset <- which(sset)

    dis[-sset] <- Inf
  }

  ## extract the list of donor sites
  if(rm.target)
    donor <- sort(order(dis)[seq(nk)+1])
  else
    donor <- sort(order(dis)[seq(nk)])

  if(!is.null(distr))
    fit <- rfaRegFit(obj, sset = donor, rtn = rtn, distr = distr)
  else
    fit <- rfaRegFit(obj, sset = donor, rtn = rtn)

  ## Predict flood quantiles at target
  if(is.null(indx) & !is.null(target))
    indx <- obj$indx[target]

  if(is.null(indx))
    indx <- 1

  qhat <- indx %o% fit$rgc
  rownames(qhat) <- names(indx)
  colnames(qhat) <- as.character(rtn)

  rmom <- as.numeric(fit$rmom$ratios)
  rmom[1] <- fit$rmom$lambdas[1]

  ## Organize and return result
  ans <- list(donor = donor, rmom = rmom, para = fit$para, qhat = qhat,
              gof = fit$gof)

  class(ans) <- c('poolGrp','list')

  return(ans)
}



#' @export
print.poolGrp <-
function(obj, long = FALSE){
  cat("\n\n## RFA using pooling group and L-moments ##\n\n")

  ## Print only if we want a long report
  if(long & !is.null(obj$tst)){
    cat("\nRegional L-moments\n")
    print(obj$tst$rmom)

    cat('\n')
    print(obj$tst)
  }

  cat('\nRegional Distribution : ',obj$para$type,'\n')
  print(obj$para$para)

  cat("\nFlood quantiles\n")
  print(obj$qhat)
}

########################################################################
#' Regional Frequency analysis
#'
#' Prepare data for regional frequency analysis. The index-flood is extracted
#' and L-moments of the standardized observation are computed
#'
#' @param Form,lx,x A list of sites or the combinaison of a formula and a data frame.
#' The formula should be 'value ~ siteID'.
#'
#' @param indx Type of index-Flood. Either 'mean' or 'median'.
#'
#'
#' @export
rfa <- function(x,...) UseMethod('rfa',x)


#' @export
#' @rdname rfa
rfa.list <- function(lx, indx = 'median'){

  ## extract the index-flood
  if(indx == 'median')
    scaleFactor <- sapply(lx, median, na.rm = TRUE)
  else if(indx == 'mean')
    scaleFactor <- sapply(lx, mean, na.rm = TRUE)
  else
    scaleFactor <- indx

  ## Standardize the data
  lx <- mapply('/', lx, scaleFactor)
  n <- sapply(lx, length)

  ## Compute the L-moments ratios
  lmm <- t(sapply(lx, Lmoments::Lmoments, na.rm = TRUE, rmax = 5))
  lmm[,3] <- lmm[,3]/lmm[,2]
  lmm[,4] <- lmm[,4]/lmm[,2]
  lmm[,5] <- lmm[,5]/lmm[,2]
  lmm[,2] <- lmm[,2]/lmm[,1]

  colnames(lmm) <- c('l_1', 't', 't_3', 't_4', 't_5')

  ans <- list(indx = scaleFactor, n = n, lmom = lmm)

  class(ans) <- c('rfa','list')

  return(ans)

}

#' @export
rfa.formula <-
function(form,x, indx = 'median'){
  x <- model.frame(form, as.data.frame(x))
  lx <- split(x[,1], as.character(x[,2]))

  rfa.list(lx, indx)
}


#' @export
#' @rdname rfa
rfa.default <- function(lmom, indx = 'median'){
  ans <- list(indx = indx, lmom = lmm)
  class(ans) <- c('rfa','list')
}

########################################################################
#' Fit a regional growth curve
#'
#' Fit the regional growth curve of a subset of site. If the distribution
#' is not specified, the best between \code{gev}, \code{glo},\code{ln3} and
#' \code{pe3} is chosen according to the goodness of fit of the Kurtosis (L-ratios)
#'
#' @param obj Output from \code{\link{rfa}}
#' @param distr Familly of distributions. See \link{lmom2par}.
#' @param sset Subset of values to fit. Default all observations.
#' @param rtn Vector of return periods
#' @param w Vector of weights. Default record lengths
#'
#' @export
rfaRegFit <-
function(obj, distr = NULL, sset = NULL, gof = TRUE,
         rtn = c(2,5,10,20,50,100), w = NULL){

  if(is.null(sset))
    sset <- seq(nrow(obj$lmom))


  ## Extract subset of data
  l <- obj$lmom[sset,]

  ## Default use weights proportional to record lengths
  if(is.null(w)) w <- obj$n[sset]

  ww <- w / sum(w)

  ## compute regional L-moments
  rmom <- apply(l,2, function(z,w) sum(w*z), ww)
  rmom <- lmomco::vec2lmom(rmom)

  ## Compute parameters and goodness-of-fit for common familly of distribution


  ## If the distribution not specified, use the best fit according to kurtosis
  if(is.null(distr)){
    ldistr <- c('gev','glo','ln3','pe3')
    lpara <- lapply(ldistr, function(d) lmomco::lmom2par(rmom,d))
    lkur <- sapply(lpara, function(z) lmomco::par2lmom(z)$ratios[4])
    para <- lpara[[which.min(abs(lkur-rmom$ratios[4]))]]

  } else
    para <- lmomco::lmom2par(rmom,distr)

  ## Compute regional growth curve
  rgc <- lmomco::qlmomco(1-1/rtn,para)

  ans <- list(rmom = rmom, para = para, rgc = rgc)

  class(ans) <- c('rfaRegFit','list')

  return(ans)
}

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

 cat('\nRegional L-moments\n\n')
 rmom <- as.numeric(obj$rmom$ratios)
 rmom[1] <- obj$rmom$lambdas[1]
 names(rmom) <- c('l1','t','t_3','t_4','t_5')
 print(rmom)

 cat('\nDistribution: ',obj$para$type,'\n')
 cat('\nRegional parameters\n')
 print(obj$para$para)

 cat('\nRegional growth curve\n')
 print(obj$rgc)

}

#' @export
rfaTst <- function(obj, sset = NULL, nsim = 1000){

  if(is.null(sset))
    sset <- seq_along(obj$indx)

  xd <- as.regdata(data.frame(name = seq_along(obj$indx)[sset],
                              n = obj$n[sset],
                              obj$lmom[sset,]))

  return(lmomRFA::regtst(xd, nsim))
}