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 ofexp
andobs
are not identical. The idea would be to add thememb_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. Ifref = NULL
, the climatological forecast is used (defined as no anomaly). Ifref
is an array, the RMS of bothexp
andref
are computed individually, and the RMSSS is computed asRMSSS = 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?