CST_MultiMetric.R

#'Multiple Metrics applied in Multiple Model Anomalies
#'
#'@author Mishra Niti, \email{niti.mishra@bsc.es}
#'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
#'@description This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations.
#'
#'@param data a CSTools object (list) giving as output of \code{Load} function from S2dverification package.
#'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.
#'@param multimodel a logical value indicating whether a Multi-Model Mean should be computed.
#'@param names a character vector indicating the names of the models or experiments to be compared.
#'
#'@return A CSTools object as the input parameter \code{data}, without the elements \code{mod} and \code{obs} and including:
#'\itemize{
#'\item\code{$ano_exp} {an array with the same dimensions as the \code{mod} label in the input \code{data}.}
#'\item\code{$ano_exp} {an array with the same dimensions as the \code{obs} label in the input \code{data}.}
#'\item\code{$metric} {An array with two datset dimensions from the \code{mod} and \code{obs} label in the input \code{data}. If \code{multimodel} is TRUE, the greatest first dimension correspons to the Multi-Model Mean. The third dimension contains the statistics selected. For metric \code{correlation}, the third dimension is of length four and they corresponds to the lower limit of the 95\% confidence interval, the statistics itselfs, the upper limit of the 95\% confidence interval and the 95\% significance level. For metric \code{rms}, the third dimension is length three and they corresponds to the lower limit of the 95\% confidence interval, the RMSE and the upper limit of the 95\% confidence interval. For metric \code{rmsss}, the third dimension is length two and they corresponds to the statistics itselfs and the p-value of the one-sided Fisher test with Ho: RMSSS = 0.}}
#'@details The output will include attributes in element \code{ano_exp} returning a vector with the names of the models or experiments.
#'@seealso \code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link[s2dverification]{Load}}
#'@references 
#' \url{http://link.springer.com/10.1007/s00382-018-4404-z}
#' 
#'@import s2dverification
#'@import stats
#'@examples
#'mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
#'dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
#'obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
#'dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
#'lon <- seq(0, 30, 5)
#'lat <- seq(0, 25, 5)
#'dat <- list(mod = mod, obs = obs, lat = lat, lon = lon)
#'a <- CST_MultiMetric(dat)
#'str(a)
#'@export
CST_MultiMetric <- function(data, metric = "correlation", multimodel = TRUE, names = NULL) {
  if (!is.list(data)) {
    stop("Parameter 'data' must be a list as output of Load function
from s2dverification package.")
  }
  if (!(all(c('mod', 'obs', 'lon', 'lat') %in% names(data)) && length(data) >= 4)) {
      stop("Parameter 'data' must contain at least 4 elements 'mod', 'obs', 
           'lon', 'lat'.")
  }
  if (!is.null(names(dim(data$mod))) & !is.null(names(dim(data$obs)))) {
    if (all(names(dim(data$mod)) %in% names(dim(data$obs)))) {
      dimnames <- names(dim(data$mod))
    } else {
      stop("Dimension names of elements 'mod' and 'obs' from parameter 'data' 
           should have the same name dimmension.")
    }
  } else {
    stop("Elements 'mod' and 'obs' from parameter 'data' should have 
         dimmension names.")
  }
  
  if (is.null(names)) {
    if (!is.null(names(data$Datasets$exp))) {
      names <-  names(data$Datasets$exp)
    } else {
      names <- paste0('mod', 1 : dim(data$mod)[1])
    }
  }
  if (!is.null(names)) {
    if (length(names) != dim(data$mod)[1]) {
      stop("Parameter 'names' must be the same length as the number of 
           models in the element 'mod' of parameter 'data'.")
    }
  }
  if (!is.logical(multimodel)) {
    stop("Parameter 'multimodel' must be a logical value.")
  }
  if (length(multimodel) > 1) {
    multimodel <- multimodel[1]
    warning("Parameter 'multimodel' has length > 1 and only the first element will be used.")
  }
  if (length(metric) > 1) {
    metric <- metric[1]
    warning("Parameter 'multimodel' has length > 1 and only the first element will be used.")
  }
  ano <- Ano_CrossValid(var_exp = data$mod, var_obs = data$obs)
  
  # seasonal average of anomalies per model
  AvgExp <- MeanListDim(ano$ano_exp, narm = T, c(2, 4))
  AvgObs <- MeanListDim(ano$ano_obs, narm = T, c(2, 4))
  
  # indv model correlation
  if (metric == 'correlation') {
    corr <- Corr(AvgExp, AvgObs, posloop = 1, poscor = 2)
  } else if (metric == 'rms') {
    corr <- RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2)
  } else if (metric == 'rmsss') {
    corr <- RMSSS(AvgExp, AvgObs, posloop = 1, posRMS = 2)
  } else {
    stop("Parameter 'metric' must be a character string indicating one of the 
         options: 'correlation', 'rms' or 'rmse'.")
  }
  if (multimodel == TRUE) {
    # seasonal avg of anomalies for multi-model
    AvgExp_MMM <- MeanListDim(AvgExp, narm = TRUE, 1)
    AvgObs_MMM <- MeanListDim(AvgObs, narm = TRUE, 1)
    # multi model correlation
    if (metric == 'correlation') {
      corr_MMM <- Corr(var_exp = InsertDim(AvgExp_MMM, 1, 1),
                       var_obs = InsertDim(AvgObs_MMM, 1, 1),
                       posloop = 1, poscor = 2)
    } else if (metric == 'rms') {
      corr_MMM <- RMS(var_exp = InsertDim(AvgExp_MMM, 1, 1),
                      var_obs = InsertDim(AvgObs_MMM, 1, 1),
                      posloop = 1, posRMS = 2)
    } else if (metric == 'rmsss') {
      corr_MMM <- RMSSS(var_exp = InsertDim(AvgExp_MMM, 1, 1),
                        var_obs = InsertDim(AvgObs_MMM, 1, 1), 
                        posloop = 1, posRMS = 2)
    } 
       corr <- abind::abind(corr, corr_MMM, along = 1)
    names <- c(names, 'MMM')
  }
  
  names(dim(corr)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6])
  
  data$ano_exp <- ano$ano_exp
  data$ano_obs <- ano$ano_obs
  data$metric <- corr
  data$obs <- NULL
  data$mod <- NULL
#  print(names)
 # attributes(data$exp)$names <- NULL
  data <- data[order(names(data))]
  attr(data[[1]], "experiments") <- names
  return(data)
}