memb_dim and ref arguments in RMSSS

We are opening this issue to propose some developments for s2dv::RMSSS, but they are not urgent.

include memb_dim parameter to allow the computation of the ensemble mean inside the function. The RMSSS is a metric to evaluate deterministic forecasts. Then, it has to be computed with the ensemble mean rather than with the individual members. The function now fails if the dimensions of exp and obs are not identical. The idea would be to add the memb_dim parameter. If it is NULL, the function will do nothing. If it is a character (e.g., 'member'), the ensemble mean will be computed, and that dimension will disappear.
include ref parameter to allow the computation of the RMSSS of the forecast with respect to a reference forecast (other than the climatological forecast). This would be similar to RPSS and CRPSS. If ref = NULL, the climatological forecast is used (defined as no anomaly). If ref is an array, the RMS of both exp and ref are computed individually, and the RMSSS is computed as RMSSS = 1 - RMS_exp / RMS_ref.

I already have a version of the function with these parameters (please see below). The problem that we see is that the current function returns the p-value based on a one-sided Fisher test. My function returns the significance based on a two-sided Random Walk test (sorry for coming back to significance issues).

library(multiApply)
library(s2dv)

RMSSS_edited <- function(exp, obs, ref = NULL, ref_value = 0, time_dim = 'year', member_dim = 'member', na.rm = FALSE, ncores = 1){
  
  ## ref_value is only used if ref = NULL
  ## for typical climatological forecast, ref_value should be 0
  ## for normalised climatological forecast, ref_value should be 1
  
  ## Ensemble means
  if (!is.null(member_dim)){
    exp <- multiApply::Apply(data = exp, target_dims = member_dim, fun = mean, na.rm = na.rm, ncores = ncores)$output1
  }
  if (!is.null(ref)){
    if (!is.null(member_dim)){
      ref <- multiApply::Apply(data = ref, target_dims = member_dim, fun = mean, na.rm = na.rm, ncores = ncores)$output1
    }
  } else { ## Climatological forecast
    ref <- array(data = ref_value, dim = dim(obs))
  }
  
  output <- multiApply::Apply(data = list(exp = exp, obs = obs, ref = ref),
                              target_dims = time_dim, fun = .RMSSS_edited,
                              time_dim = time_dim, ncores = ncores)
  return(output)
}

.RMSSS_edited <- function(exp, obs, ref, time_dim){
  
  ## RMS
  exp <- s2dv::InsertDim(data = exp, posdim = 1, lendim = 1, name = 'dat')
  ref <- s2dv::InsertDim(data = ref, posdim = 1, lendim = 1, name = 'dat')
  obs <- s2dv::InsertDim(data = obs, posdim = 1, lendim = 1, name = 'dat')
  rms_exp <- drop(s2dv::RMS(exp = exp, obs = obs, time_dim = time_dim, dat_dim = 'dat', comp_dim = NULL, conf = FALSE, ncores = NULL)$rms)
  rms_ref <- drop(s2dv::RMS(exp = ref, obs = obs, time_dim = time_dim, dat_dim = 'dat', comp_dim = NULL, conf = FALSE, ncores = NULL)$rms)
  exp <- drop(exp)
  ref <- drop(ref)
  obs <- drop(obs)
  
  ## Error
  error_exp <- array(data = abs(exp - obs), dim = c(time = length(obs)))
  error_ref <- array(data = abs(ref - obs), dim = c(time = length(obs)))
  
  ## RMSSS and significance
  output <- NULL
  output$rmsss <- 1 - rms_exp / rms_ref
  output$sign <- s2dv::RandomWalkTest(skill_A = error_exp, skill_B = error_ref, time_dim = 'time', ncores = 1)$signif
  
  return(output)
}

Do we have a meeting to discuss this all together?

Thank you!
@nmilders @vagudets Carlos