From c770d4b2634805627bda95107d2373f9fc1ad9f1 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Tue, 30 Aug 2022 16:49:47 +0200 Subject: [PATCH 01/17] first version --- R/Bias.R | 32 +++++++++++++++++++++++++ R/BiasSS.R | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 100 insertions(+) create mode 100644 R/Bias.R create mode 100644 R/BiasSS.R diff --git a/R/Bias.R b/R/Bias.R new file mode 100644 index 0000000..e3d5242 --- /dev/null +++ b/R/Bias.R @@ -0,0 +1,32 @@ + +# exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +# obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) +# bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) + +Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { + + ## Checks + + + ## Ensemble mean + if (!is.null(memb_dim)) { + exp <- MeanDims(exp, memb_dim, na.rm = na.rm) + } + + ## Mean bias + bias <- multiApply::Apply(data = list(exp, obs), + target_dims = time_dim, + fun = .Bias, + ncores = ncores)$output1 + ## Return the mean bias + bias <- MeanDims(bias, time_dim, na.rm = na.rm) + + return(bias) +} + +.Bias <- function(exp, obs) { + + bias <- exp - obs + + return(bias) +} diff --git a/R/BiasSS.R b/R/BiasSS.R new file mode 100644 index 0000000..f71e5b5 --- /dev/null +++ b/R/BiasSS.R @@ -0,0 +1,68 @@ + + +# exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +# ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +# obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) +# biasSS1 <- BiasSS(exp = exp, obs = obs, ref = ref, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +# biasSS2 <- BiasSS(exp = exp, obs = obs, ref = NULL, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) + +BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { + + ## Checks + + + ## Ensemble mean + if (!is.null(memb_dim)) { + exp <- MeanDims(exp, memb_dim, na.rm = na.rm) + if (!is.null(ref)) { + ref <- MeanDims(ref, memb_dim, na.rm = na.rm) + } + } + + ## Mean bias skill score + if (is.null(ref)) { + ss <- multiApply::Apply(data = list(exp, obs), + target_dims = time_dim, + fun = .BiasSS, na.rm = na.rm, + ref = ref, ncores = ncores) + + } else { + ss <- multiApply::Apply(data = list(exp, obs, ref), + target_dims = time_dim, + fun = .BiasSS, na.rm = na.rm, + ncores = ncores) + } + return(ss) +} + +.BiasSS <- function(exp, obs, ref = NULL, na.rm = FALSE) { + + if (isTRUE(na.rm)) { + if (is.null(ref)) { + good_values <- !is.na(exp) & !is.na(obs) + exp <- exp[good_values] + obs <- obs[good_values] + } else { + good_values <- !is.na(exp) & !is.na(ref) & !is.na(obs) + exp <- exp[good_values] + ref <- ref[good_values] + obs <- obs[good_values] + } + } + + ## Bias of the exp + bias_exp <- .Bias(exp = exp, obs = obs) + + ## Bias of the ref + if (is.null(ref)) { ## Climatological forecast + ref <- mean(obs, na.rm = na.rm) + } + bias_ref <- .Bias(exp = ref, obs = obs) + + ## Skill score and significance + ss <- list() + ss$biasSS <- 1 - mean(abs(bias_exp), na.rm = na.rm) / mean(abs(bias_ref), na.rm = na.rm) + ss$sign <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif + + return(ss) +} -- GitLab From 1e3303ac34252b6461083e132b0073107971820c Mon Sep 17 00:00:00 2001 From: ALBA LLABRES Date: Wed, 31 Aug 2022 12:45:52 +0200 Subject: [PATCH 02/17] first version documentation added --- R/Bias.R | 39 ++++++++++++++++++++++++++++---- R/BiasSS.R | 65 ++++++++++++++++++++++++++++++++++++++++++++++++------ 2 files changed, 93 insertions(+), 11 deletions(-) diff --git a/R/Bias.R b/R/Bias.R index e3d5242..3f76dc7 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -1,7 +1,38 @@ - -# exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -# obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) -# bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +#'Compute the Mean Bias +#' +#'The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference +#'between the ensemble mean forecast and the observations. It is a deterministic +#'metric. Positive values indicate that the forecasts are on average too high +#'and negative values indicate that the forecasts are on average too low; however, +#'it gives no information about the typical magnitude of individual forecast errors. +#' +#'@param exp A named numerical array of the forecast with at least time +#' dimension. +#'@param obs A named numerical array of the observation with at least time +#' dimension. The dimensions must be the same as 'exp' except 'memb_dim' and +#' 'dat_dim'. +#'@param time_dim A character string indicating the name of the time dimension. +#' The default value is 'sdate'. +#'@param memb_dim A character string indicating the name of the member dimension +#' to compute the ensemble mean; it should be set to NULL if the parameter 'exp' +#' is already the ensemble mean. The default value is NULL. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. The default value is NULL. +#' +#'@return +#'A numerical array of Bias with dimensions the dimensions of +#''exp' except 'time_dim' and 'memb_dim' dimensions. +#' +#'@references +#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 +#' +#'@examples +#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) +#'bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +#' +#'@import multiApply +#'@export Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { diff --git a/R/BiasSS.R b/R/BiasSS.R index f71e5b5..ef94c85 100644 --- a/R/BiasSS.R +++ b/R/BiasSS.R @@ -1,10 +1,61 @@ - - -# exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -# ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -# obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) -# biasSS1 <- BiasSS(exp = exp, obs = obs, ref = ref, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) -# biasSS2 <- BiasSS(exp = exp, obs = obs, ref = NULL, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +#'Compute the Mean Bias Skill Score +#' +#'The Mean Bias Skill Score is based on the Mean Error (Wilks, 2011) +#'between the ensemble mean forecast and the observations. It measures +#'the accuracy of the forecast and it can be used to assess whether a +#'forecast presents an improvement or a worsening with respect to a reference +#'forecast. The Mean Bias Skill Score ranges between minus infinite and 1. +#'Positive values indicate that the forecast has higher skill than the +#'reference forecast, while negative values indicate that it has a lower skill. +#'Examples of reference forecasts are the climatological forecast (average of +#'the observations), a previous model version, or another model. It is computed +#'as BiasSS = 1 - Bias_exp / Bias_ref. The statistical significance is obtained +#'based on a Random Walk test at the 95% confidence level (DelSole and Tippett, +#'2016). +#' +#'@param exp A named numerical array of the forecast with at least time +#' dimension. +#'@param obs A named numerical array of the observation with at least time +#' dimension. The dimensions must be the same as 'exp' except 'memb_dim' and +#' 'dat_dim'. +#'@param ref A named numerical array of the reference forecast data with at +#' least time and member dimension. The dimensions must be the same as 'exp' +#' except 'memb_dim' and 'dat_dim'. If there is only one reference dataset, +#' it should not have dataset dimension. If there is corresponding reference +#' for each experiement, the dataset dimension must have the same length as in +#' 'exp'. If 'ref' is NULL, the climatological forecast is used as reference +#' forecast. The default value is NULL. +#'@param time_dim A character string indicating the name of the time dimension. +#' The default value is 'sdate'. +#'@param memb_dim A character string indicating the name of the member dimension +#' to compute the ensemble mean; it should be set to NULL if the parameter 'exp' +#' is already the ensemble mean. The default value is NULL. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. The default value is NULL. +#' +#'@return +#'\item{$biasSS}{ +#' A numerical array of BiasSS with dimensions the dimensions of +#' 'exp' except 'time_dim' and 'memb_dim' dimensions. +#'} +#'\item{$sign}{ +#' A logical array of the statistical significance of the BiasSS +#' with the same dimensions as $biasSS. +#'} +#' +#'@references +#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 +#'DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1 +#' +#'@examples +#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +#'ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) +#'biasSS1 <- BiasSS(exp = exp, obs = obs, ref = ref, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +#'biasSS2 <- BiasSS(exp = exp, obs = obs, ref = NULL, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +#' +#'@import multiApply +#'@export BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { -- GitLab From 95f5fa48c3fe3f0b08bbd2039df4e716010e117b Mon Sep 17 00:00:00 2001 From: ALBA LLABRES Date: Wed, 31 Aug 2022 13:23:28 +0200 Subject: [PATCH 03/17] documentation improvement --- R/BiasSS.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/BiasSS.R b/R/BiasSS.R index ef94c85..e59189b 100644 --- a/R/BiasSS.R +++ b/R/BiasSS.R @@ -2,9 +2,9 @@ #' #'The Mean Bias Skill Score is based on the Mean Error (Wilks, 2011) #'between the ensemble mean forecast and the observations. It measures -#'the accuracy of the forecast and it can be used to assess whether a -#'forecast presents an improvement or a worsening with respect to a reference -#'forecast. The Mean Bias Skill Score ranges between minus infinite and 1. +#'the accuracy of the forecast in comparison with a reference forecast to assess +#'whether the forecast presents an improvement or a worsening with respect to +#'that reference. The Mean Bias Skill Score ranges between minus infinite and 1. #'Positive values indicate that the forecast has higher skill than the #'reference forecast, while negative values indicate that it has a lower skill. #'Examples of reference forecasts are the climatological forecast (average of -- GitLab From 6b34657b30eac316a819bf9d52837194b36bb549 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Wed, 31 Aug 2022 15:09:33 +0200 Subject: [PATCH 04/17] added checks (dat_dim missing) --- R/Bias.R | 67 +++++++++++++++++++++++++++++++++--- R/BiasSS.R | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 159 insertions(+), 7 deletions(-) diff --git a/R/Bias.R b/R/Bias.R index 3f76dc7..15c3f4e 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -16,6 +16,8 @@ #'@param memb_dim A character string indicating the name of the member dimension #' to compute the ensemble mean; it should be set to NULL if the parameter 'exp' #' is already the ensemble mean. The default value is NULL. +#'@param na.rm A logical value indicating if NAs should be removed (TRUE) or +#' kept (FALSE) for computation. The default value is FALSE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -33,11 +35,68 @@ #' #'@import multiApply #'@export - Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { - ## Checks - + # Check inputs + ## exp and obs (1) + if (!is.array(exp) | !is.numeric(exp)) + stop('Parameter "exp" must be a numeric array.') + if (!is.array(obs) | !is.numeric(obs)) + stop('Parameter "obs" must be a numeric array.') + if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) | + any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) { + stop("Parameter 'exp' and 'obs' must have dimension names.") + } + ## time_dim + if (!is.character(time_dim) | length(time_dim) != 1) + stop('Parameter "time_dim" must be a character string.') + if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) { + stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.") + } + ## memb_dim + if (!is.character(memb_dim) | length(memb_dim) > 1) { + stop("Parameter 'memb_dim' must be a character string.") + } + if (!memb_dim %in% names(dim(exp))) { + stop("Parameter 'memb_dim' is not found in 'exp' dimension.") + } + ## dat_dim + # if (!is.null(dat_dim)) { + # if (!is.character(dat_dim) | length(dat_dim) > 1) { + # stop("Parameter 'dat_dim' must be a character string.") + # } + # if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) { + # stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.", + # " Set it as NULL if there is no dataset dimension.") + # } + # } + ## exp and obs (2) + name_exp <- sort(names(dim(exp))) + name_obs <- sort(names(dim(obs))) + name_exp <- name_exp[-which(name_exp == memb_dim)] + if (memb_dim %in% name_obs) { + name_obs <- name_obs[-which(name_obs == memb_dim)] + } + # if (!is.null(dat_dim)) { + # name_exp <- name_exp[-which(name_exp == dat_dim)] + # name_obs <- name_obs[-which(name_obs == dat_dim)] + # } + if (!identical(length(name_exp), length(name_obs)) | + !identical(dim(exp)[name_exp], dim(obs)[name_obs])) { + stop(paste0("Parameter 'exp' and 'obs' must have same length of ", + "all dimensions except 'memb_dim'")) # and 'dat_dim'.")) + } + ## na.rm + if (!is.logical(na.rm) | length(na.rm) > 1) { + stop("Parameter 'na.rm' must be one logical value.") + } + ## ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | + length(ncores) > 1) { + stop("Parameter 'ncores' must be either NULL or a positive integer.") + } + } ## Ensemble mean if (!is.null(memb_dim)) { @@ -56,8 +115,6 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, n } .Bias <- function(exp, obs) { - bias <- exp - obs - return(bias) } diff --git a/R/BiasSS.R b/R/BiasSS.R index e59189b..d456e8f 100644 --- a/R/BiasSS.R +++ b/R/BiasSS.R @@ -30,6 +30,8 @@ #'@param memb_dim A character string indicating the name of the member dimension #' to compute the ensemble mean; it should be set to NULL if the parameter 'exp' #' is already the ensemble mean. The default value is NULL. +#'@param na.rm A logical value indicating if NAs should be removed (TRUE) or +#' kept (FALSE) for computation. The default value is FALSE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -59,8 +61,101 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { - ## Checks - + # Check inputs + ## exp, obs, and ref (1) + if (!is.array(exp) | !is.numeric(exp)) { + stop("Parameter 'exp' must be a numeric array.") + } + if (!is.array(obs) | !is.numeric(obs)) { + stop("Parameter 'obs' must be a numeric array.") + } + if (any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) | + any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) { + stop("Parameter 'exp' and 'obs' must have dimension names.") + } + if (!is.null(ref)) { + if (!is.array(ref) | !is.numeric(ref)) + stop("Parameter 'ref' must be a numeric array.") + if (any(is.null(names(dim(ref))))| any(nchar(names(dim(ref))) == 0)) { + stop("Parameter 'ref' must have dimension names.") + } + } + ## time_dim + if (!is.character(time_dim) | length(time_dim) != 1) { + stop("Parameter 'time_dim' must be a character string.") + } + if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) { + stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.") + } + if (!is.null(ref) & !time_dim %in% names(dim(ref))) { + stop("Parameter 'time_dim' is not found in 'ref' dimension.") + } + ## memb_dim + if (!is.character(memb_dim) | length(memb_dim) > 1) { + stop("Parameter 'memb_dim' must be a character string.") + } + if (!memb_dim %in% names(dim(exp))) { + stop("Parameter 'memb_dim' is not found in 'exp' dimension.") + } + if (!is.null(ref) & !memb_dim %in% names(dim(ref))) { + stop("Parameter 'memb_dim' is not found in 'ref' dimension.") + } + ## dat_dim + # if (!is.null(dat_dim)) { + # if (!is.character(dat_dim) | length(dat_dim) > 1) { + # stop("Parameter 'dat_dim' must be a character string.") + # } + # if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) { + # stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.", + # " Set it as NULL if there is no dataset dimension.") + # } + # } + ## exp, obs, and ref (2) + name_exp <- sort(names(dim(exp))) + name_obs <- sort(names(dim(obs))) + name_exp <- name_exp[-which(name_exp == memb_dim)] + if (memb_dim %in% name_obs) { + name_obs <- name_obs[-which(name_obs == memb_dim)] + } + # if (!is.null(dat_dim)) { + # name_exp <- name_exp[-which(name_exp == dat_dim)] + # name_obs <- name_obs[-which(name_obs == dat_dim)] + # } + if (!identical(length(name_exp), length(name_obs)) | + !identical(dim(exp)[name_exp], dim(obs)[name_obs])) { + stop(paste0("Parameter 'exp' and 'obs' must have same length of ", + "all dimensions except 'memb_dim'")) # and 'dat_dim'.")) + } + if (!is.null(ref)) { + name_ref <- sort(names(dim(ref))) + name_ref <- name_ref[-which(name_ref == memb_dim)] + # if (!is.null(dat_dim)) { + # if (dat_dim %in% name_ref) { + # if (!identical(dim(exp)[dat_dim], dim(ref)[dat_dim])) { + # stop(paste0("If parameter 'ref' has dataset dimension, it must be", + # " equal to dataset dimension of 'exp'.")) + # } + # name_ref <- name_ref[-which(name_ref == dat_dim)] + # } + # } + if (!identical(length(name_exp), length(name_ref)) | + !identical(dim(exp)[name_exp], dim(ref)[name_ref])) { + stop(paste0("Parameter 'exp' and 'ref' must have the same length of ", + "all dimensions except 'memb_dim'")) #, + # " and 'dat_dim' if there is only one reference dataset.")) + } + } + ## na.rm + if (!is.logical(na.rm) | length(na.rm) > 1) { + stop("Parameter 'na.rm' must be one logical value.") + } + ## ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | + length(ncores) > 1) { + stop("Parameter 'ncores' must be either NULL or a positive integer.") + } + } ## Ensemble mean if (!is.null(memb_dim)) { -- GitLab From 73dfed032efb96e9b3aa8e2c0d767703255f3943 Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 5 Sep 2022 18:17:27 +0200 Subject: [PATCH 05/17] Refine doc and checks for memb_dim --- NAMESPACE | 2 ++ R/Bias.R | 51 ++++++++++++++++++++------------ R/BiasSS.R | 78 ++++++++++++++++++++++++++++--------------------- man/Bias.Rd | 57 ++++++++++++++++++++++++++++++++++++ man/BiasSS.Rd | 81 +++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 217 insertions(+), 52 deletions(-) create mode 100644 man/Bias.Rd create mode 100644 man/BiasSS.Rd diff --git a/NAMESPACE b/NAMESPACE index 032ebf9..6b21447 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,8 @@ export(AMV) export(AnimateMap) export(Ano) export(Ano_CrossValid) +export(Bias) +export(BiasSS) export(BrierScore) export(CDORemap) export(Clim) diff --git a/R/Bias.R b/R/Bias.R index 15c3f4e..f731484 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -3,8 +3,9 @@ #'The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference #'between the ensemble mean forecast and the observations. It is a deterministic #'metric. Positive values indicate that the forecasts are on average too high -#'and negative values indicate that the forecasts are on average too low; however, -#'it gives no information about the typical magnitude of individual forecast errors. +#'and negative values indicate that the forecasts are on average too low. +#'However, it gives no information about the typical magnitude of individual +#'forecast errors. #' #'@param exp A named numerical array of the forecast with at least time #' dimension. @@ -22,8 +23,8 @@ #' computation. The default value is NULL. #' #'@return -#'A numerical array of Bias with dimensions the dimensions of -#''exp' except 'time_dim' and 'memb_dim' dimensions. +#'A numerical array of bias with dimensions of 'exp' except 'time_dim' and +#''memb_dim' dimensions. #' #'@references #'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 @@ -34,15 +35,16 @@ #'bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) #' #'@import multiApply +#'@importFrom ClimProjDiags Subset #'@export Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { # Check inputs ## exp and obs (1) if (!is.array(exp) | !is.numeric(exp)) - stop('Parameter "exp" must be a numeric array.') + stop("Parameter 'exp' must be a numeric array.") if (!is.array(obs) | !is.numeric(obs)) - stop('Parameter "obs" must be a numeric array.') + stop("Parameter 'obs' must be a numeric array.") if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) | any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) { stop("Parameter 'exp' and 'obs' must have dimension names.") @@ -54,11 +56,21 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, n stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.") } ## memb_dim - if (!is.character(memb_dim) | length(memb_dim) > 1) { - stop("Parameter 'memb_dim' must be a character string.") - } - if (!memb_dim %in% names(dim(exp))) { - stop("Parameter 'memb_dim' is not found in 'exp' dimension.") + if (!is.null(memb_dim)) { + if (!is.character(memb_dim) | length(memb_dim) > 1) { + stop("Parameter 'memb_dim' must be a character string.") + } + if (!memb_dim %in% names(dim(exp))) { + stop("Parameter 'memb_dim' is not found in 'exp' dimension.") + } + if (memb_dim %in% names(dim(obs))) { + if (identical(as.numeric(dim(obs)[memb_dim]), 1)) { + obs <- ClimProjDiags::Subset(x = obs, along = memb_dim, indices = 1, drop = 'selected') + } else { + stop("Not implemented for observations with members ('obs' can have 'memb_dim', ", + "but it should be of length = 1).") + } + } } ## dat_dim # if (!is.null(dat_dim)) { @@ -73,9 +85,8 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, n ## exp and obs (2) name_exp <- sort(names(dim(exp))) name_obs <- sort(names(dim(obs))) - name_exp <- name_exp[-which(name_exp == memb_dim)] - if (memb_dim %in% name_obs) { - name_obs <- name_obs[-which(name_obs == memb_dim)] + if (!is.null(memb_dim)) { + name_exp <- name_exp[-which(name_exp == memb_dim)] } # if (!is.null(dat_dim)) { # name_exp <- name_exp[-which(name_exp == dat_dim)] @@ -97,17 +108,19 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, n stop("Parameter 'ncores' must be either NULL or a positive integer.") } } - + + ############################### + ## Ensemble mean if (!is.null(memb_dim)) { exp <- MeanDims(exp, memb_dim, na.rm = na.rm) } ## Mean bias - bias <- multiApply::Apply(data = list(exp, obs), - target_dims = time_dim, - fun = .Bias, - ncores = ncores)$output1 + bias <- Apply(data = list(exp, obs), + target_dims = time_dim, + fun = .Bias, + ncores = ncores)$output1 ## Return the mean bias bias <- MeanDims(bias, time_dim, na.rm = na.rm) diff --git a/R/BiasSS.R b/R/BiasSS.R index d456e8f..ee748ed 100644 --- a/R/BiasSS.R +++ b/R/BiasSS.R @@ -1,17 +1,17 @@ #'Compute the Mean Bias Skill Score #' -#'The Mean Bias Skill Score is based on the Mean Error (Wilks, 2011) -#'between the ensemble mean forecast and the observations. It measures -#'the accuracy of the forecast in comparison with a reference forecast to assess -#'whether the forecast presents an improvement or a worsening with respect to -#'that reference. The Mean Bias Skill Score ranges between minus infinite and 1. -#'Positive values indicate that the forecast has higher skill than the -#'reference forecast, while negative values indicate that it has a lower skill. -#'Examples of reference forecasts are the climatological forecast (average of -#'the observations), a previous model version, or another model. It is computed -#'as BiasSS = 1 - Bias_exp / Bias_ref. The statistical significance is obtained -#'based on a Random Walk test at the 95% confidence level (DelSole and Tippett, -#'2016). +#'The Mean Bias Skill Score is based on the Mean Error (Wilks, 2011) between the +#'ensemble mean forecast and the observations. It measures the accuracy of the +#'forecast in comparison with a reference forecast to assess whether the +#'forecast presents an improvement or a worsening with respect to that +#'reference. The Mean Bias Skill Score ranges between minus infinite and 1. +#'Positive values indicate that the forecast has higher skill than the reference +#'forecast, while negative values indicate that it has a lower skill. Examples +#'of reference forecasts are the climatological forecast (average of the +#'observations), a previous model version, or another model. It is computed as +#'\code{BiasSS = 1 - Bias_exp / Bias_ref}. The statistical significance is +#'obtained based on a Random Walk test at the 95% confidence level (DelSole and +#'Tippett, 2016). #' #'@param exp A named numerical array of the forecast with at least time #' dimension. @@ -19,17 +19,17 @@ #' dimension. The dimensions must be the same as 'exp' except 'memb_dim' and #' 'dat_dim'. #'@param ref A named numerical array of the reference forecast data with at -#' least time and member dimension. The dimensions must be the same as 'exp' -#' except 'memb_dim' and 'dat_dim'. If there is only one reference dataset, -#' it should not have dataset dimension. If there is corresponding reference -#' for each experiement, the dataset dimension must have the same length as in -#' 'exp'. If 'ref' is NULL, the climatological forecast is used as reference -#' forecast. The default value is NULL. +#' least time dimension. The dimensions must be the same as 'exp' except +#' 'memb_dim' and 'dat_dim'. If there is only one reference dataset, it should +#' not have dataset dimension. If there is corresponding reference for each +#' experiement, the dataset dimension must have the same length as in 'exp'. If +#' 'ref' is NULL, the climatological forecast is used as reference forecast. +#' The default value is NULL. #'@param time_dim A character string indicating the name of the time dimension. #' The default value is 'sdate'. #'@param memb_dim A character string indicating the name of the member dimension #' to compute the ensemble mean; it should be set to NULL if the parameter 'exp' -#' is already the ensemble mean. The default value is NULL. +#' and 'ref' are already the ensemble mean. The default value is NULL. #'@param na.rm A logical value indicating if NAs should be removed (TRUE) or #' kept (FALSE) for computation. The default value is FALSE. #'@param ncores An integer indicating the number of cores to use for parallel @@ -91,14 +91,21 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na stop("Parameter 'time_dim' is not found in 'ref' dimension.") } ## memb_dim - if (!is.character(memb_dim) | length(memb_dim) > 1) { - stop("Parameter 'memb_dim' must be a character string.") - } - if (!memb_dim %in% names(dim(exp))) { - stop("Parameter 'memb_dim' is not found in 'exp' dimension.") - } - if (!is.null(ref) & !memb_dim %in% names(dim(ref))) { - stop("Parameter 'memb_dim' is not found in 'ref' dimension.") + if (!is.null(memb_dim)) { + if (!is.character(memb_dim) | length(memb_dim) > 1) { + stop("Parameter 'memb_dim' must be a character string.") + } + if (!memb_dim %in% names(dim(exp))) { + stop("Parameter 'memb_dim' is not found in 'exp' dimension.") + } + if (memb_dim %in% names(dim(obs))) { + if (identical(as.numeric(dim(obs)[memb_dim]), 1)) { + obs <- ClimProjDiags::Subset(x = obs, along = memb_dim, indices = 1, drop = 'selected') + } else { + stop("Not implemented for observations with members ('obs' can have 'memb_dim', ", + "but it should be of length = 1).") + } + } } ## dat_dim # if (!is.null(dat_dim)) { @@ -113,9 +120,8 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na ## exp, obs, and ref (2) name_exp <- sort(names(dim(exp))) name_obs <- sort(names(dim(obs))) - name_exp <- name_exp[-which(name_exp == memb_dim)] - if (memb_dim %in% name_obs) { - name_obs <- name_obs[-which(name_obs == memb_dim)] + if (!is.null(memb_dim)) { + name_exp <- name_exp[-which(name_exp == memb_dim)] } # if (!is.null(dat_dim)) { # name_exp <- name_exp[-which(name_exp == dat_dim)] @@ -128,7 +134,9 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na } if (!is.null(ref)) { name_ref <- sort(names(dim(ref))) - name_ref <- name_ref[-which(name_ref == memb_dim)] + if (memb_dim %in% name_ref) { + name_ref <- name_ref[-which(name_ref == memb_dim)] + } # if (!is.null(dat_dim)) { # if (dat_dim %in% name_ref) { # if (!identical(dim(exp)[dat_dim], dim(ref)[dat_dim])) { @@ -156,7 +164,9 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na stop("Parameter 'ncores' must be either NULL or a positive integer.") } } - + + ############################ + ## Ensemble mean if (!is.null(memb_dim)) { exp <- MeanDims(exp, memb_dim, na.rm = na.rm) @@ -182,7 +192,9 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na } .BiasSS <- function(exp, obs, ref = NULL, na.rm = FALSE) { - + # exp and obs: [sdate] + # ref: [sdate] or NULL + if (isTRUE(na.rm)) { if (is.null(ref)) { good_values <- !is.na(exp) & !is.na(obs) diff --git a/man/Bias.Rd b/man/Bias.Rd new file mode 100644 index 0000000..a58d263 --- /dev/null +++ b/man/Bias.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Bias.R +\name{Bias} +\alias{Bias} +\title{Compute the Mean Bias} +\usage{ +Bias( + exp, + obs, + time_dim = "sdate", + memb_dim = NULL, + na.rm = FALSE, + ncores = NULL +) +} +\arguments{ +\item{exp}{A named numerical array of the forecast with at least time +dimension.} + +\item{obs}{A named numerical array of the observation with at least time +dimension. The dimensions must be the same as 'exp' except 'memb_dim' and +'dat_dim'.} + +\item{time_dim}{A character string indicating the name of the time dimension. +The default value is 'sdate'.} + +\item{memb_dim}{A character string indicating the name of the member dimension +to compute the ensemble mean; it should be set to NULL if the parameter 'exp' +is already the ensemble mean. The default value is NULL.} + +\item{na.rm}{A logical value indicating if NAs should be removed (TRUE) or +kept (FALSE) for computation. The default value is FALSE.} + +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} +} +\value{ +A numerical array of bias with dimensions of 'exp' except 'time_dim' and +'memb_dim' dimensions. +} +\description{ +The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference +between the ensemble mean forecast and the observations. It is a deterministic +metric. Positive values indicate that the forecasts are on average too high +and negative values indicate that the forecasts are on average too low. +However, it gives no information about the typical magnitude of individual +forecast errors. +} +\examples{ +exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) +bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) + +} +\references{ +Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 +} diff --git a/man/BiasSS.Rd b/man/BiasSS.Rd new file mode 100644 index 0000000..bce496e --- /dev/null +++ b/man/BiasSS.Rd @@ -0,0 +1,81 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BiasSS.R +\name{BiasSS} +\alias{BiasSS} +\title{Compute the Mean Bias Skill Score} +\usage{ +BiasSS( + exp, + obs, + ref = NULL, + time_dim = "sdate", + memb_dim = NULL, + na.rm = FALSE, + ncores = NULL +) +} +\arguments{ +\item{exp}{A named numerical array of the forecast with at least time +dimension.} + +\item{obs}{A named numerical array of the observation with at least time +dimension. The dimensions must be the same as 'exp' except 'memb_dim' and +'dat_dim'.} + +\item{ref}{A named numerical array of the reference forecast data with at +least time dimension. The dimensions must be the same as 'exp' except +'memb_dim' and 'dat_dim'. If there is only one reference dataset, it should +not have dataset dimension. If there is corresponding reference for each +experiement, the dataset dimension must have the same length as in 'exp'. If +'ref' is NULL, the climatological forecast is used as reference forecast. +The default value is NULL.} + +\item{time_dim}{A character string indicating the name of the time dimension. +The default value is 'sdate'.} + +\item{memb_dim}{A character string indicating the name of the member dimension +to compute the ensemble mean; it should be set to NULL if the parameter 'exp' +and 'ref' are already the ensemble mean. The default value is NULL.} + +\item{na.rm}{A logical value indicating if NAs should be removed (TRUE) or +kept (FALSE) for computation. The default value is FALSE.} + +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} +} +\value{ +\item{$biasSS}{ + A numerical array of BiasSS with dimensions the dimensions of + 'exp' except 'time_dim' and 'memb_dim' dimensions. +} +\item{$sign}{ + A logical array of the statistical significance of the BiasSS + with the same dimensions as $biasSS. +} +} +\description{ +The Mean Bias Skill Score is based on the Mean Error (Wilks, 2011) between the +ensemble mean forecast and the observations. It measures the accuracy of the +forecast in comparison with a reference forecast to assess whether the +forecast presents an improvement or a worsening with respect to that +reference. The Mean Bias Skill Score ranges between minus infinite and 1. +Positive values indicate that the forecast has higher skill than the reference +forecast, while negative values indicate that it has a lower skill. Examples +of reference forecasts are the climatological forecast (average of the +observations), a previous model version, or another model. It is computed as +\code{BiasSS = 1 - Bias_exp / Bias_ref}. The statistical significance is +obtained based on a Random Walk test at the 95% confidence level (DelSole and +Tippett, 2016). +} +\examples{ +exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) +obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) +biasSS1 <- BiasSS(exp = exp, obs = obs, ref = ref, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +biasSS2 <- BiasSS(exp = exp, obs = obs, ref = NULL, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) + +} +\references{ +Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 +DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1 +} -- GitLab From 0ffeef39a83272ab935e6f5ed86580cac5d2538d Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 5 Sep 2022 18:30:49 +0200 Subject: [PATCH 06/17] Correct if condition for ref ensemble mean --- R/BiasSS.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/BiasSS.R b/R/BiasSS.R index ee748ed..87edfdd 100644 --- a/R/BiasSS.R +++ b/R/BiasSS.R @@ -170,7 +170,7 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na ## Ensemble mean if (!is.null(memb_dim)) { exp <- MeanDims(exp, memb_dim, na.rm = na.rm) - if (!is.null(ref)) { + if (!is.null(ref) & memb_dim %in% names(dim(ref))) { ref <- MeanDims(ref, memb_dim, na.rm = na.rm) } } -- GitLab From 13cec946d79a8fcb660baae5be96f7ef74e4d759 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Tue, 6 Sep 2022 13:28:55 +0200 Subject: [PATCH 07/17] possibility of computing the time_mean and absolute_bias --- R/Bias.R | 48 +++++++++++++++++++++++++++++++++++------------- R/BiasSS.R | 30 +++++++++++++++--------------- 2 files changed, 50 insertions(+), 28 deletions(-) diff --git a/R/Bias.R b/R/Bias.R index f731484..316efc8 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -4,8 +4,7 @@ #'between the ensemble mean forecast and the observations. It is a deterministic #'metric. Positive values indicate that the forecasts are on average too high #'and negative values indicate that the forecasts are on average too low. -#'However, it gives no information about the typical magnitude of individual -#'forecast errors. +#'It also allows to compute the Absolute Mean Bias. #' #'@param exp A named numerical array of the forecast with at least time #' dimension. @@ -19,25 +18,30 @@ #' is already the ensemble mean. The default value is NULL. #'@param na.rm A logical value indicating if NAs should be removed (TRUE) or #' kept (FALSE) for computation. The default value is FALSE. +#'@param absolute A logical value indicating whether to compute the absolute bias. +#' The default value is FALSE. +#'@param time_dim A logical value indicating whether to compute the temporal +#' mean of the bias. The default value is TRUE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' #'@return -#'A numerical array of bias with dimensions of 'exp' except 'time_dim' and -#''memb_dim' dimensions. +#'A numerical array of bias with dimensions of 'exp' except 'time_dim' (if time_mean = T) +#'and 'memb_dim' dimensions. #' #'@references #'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 #' #'@examples -#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) -#'bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) +#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50)) +#'bias <- Bias(exp = exp, obs = obs, memb_dim = 'member') #' #'@import multiApply #'@importFrom ClimProjDiags Subset #'@export -Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { +Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, + absolute = FALSE, time_mean = TRUE, ncores = NULL) { # Check inputs ## exp and obs (1) @@ -101,6 +105,14 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, n if (!is.logical(na.rm) | length(na.rm) > 1) { stop("Parameter 'na.rm' must be one logical value.") } + ## absolute + if (!is.logical(absolute) | length(absolute) > 1) { + stop("Parameter 'absolute' must be one logical value.") + } + ## time_mean + if (!is.logical(time_mean) | length(time_mean) > 1) { + stop("Parameter 'time_mean' must be one logical value.") + } ## ncores if (!is.null(ncores)) { if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | @@ -116,18 +128,28 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, n exp <- MeanDims(exp, memb_dim, na.rm = na.rm) } - ## Mean bias + ## (Mean) Bias bias <- Apply(data = list(exp, obs), target_dims = time_dim, - fun = .Bias, + fun = .Bias, na.rm = na.rm, + absolute = absolute, + time_mean = time_mean, ncores = ncores)$output1 - ## Return the mean bias - bias <- MeanDims(bias, time_dim, na.rm = na.rm) return(bias) } -.Bias <- function(exp, obs) { +.Bias <- function(exp, obs, na.rm = FALSE, absolute = FALSE, time_mean = TRUE) { + bias <- exp - obs + + if (isTRUE(absolute)){ + bias <- abs(bias) + } + + if (isTRUE(time_mean)){ + bias <- mean(bias, na.rm = na.rm) + } + return(bias) } diff --git a/R/BiasSS.R b/R/BiasSS.R index 87edfdd..8e5da98 100644 --- a/R/BiasSS.R +++ b/R/BiasSS.R @@ -1,6 +1,6 @@ -#'Compute the Mean Bias Skill Score +#'Compute the Absolute Mean Bias Skill Score #' -#'The Mean Bias Skill Score is based on the Mean Error (Wilks, 2011) between the +#'The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks, 2011) between the #'ensemble mean forecast and the observations. It measures the accuracy of the #'forecast in comparison with a reference forecast to assess whether the #'forecast presents an improvement or a worsening with respect to that @@ -9,7 +9,7 @@ #'forecast, while negative values indicate that it has a lower skill. Examples #'of reference forecasts are the climatological forecast (average of the #'observations), a previous model version, or another model. It is computed as -#'\code{BiasSS = 1 - Bias_exp / Bias_ref}. The statistical significance is +#'\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance is #'obtained based on a Random Walk test at the 95% confidence level (DelSole and #'Tippett, 2016). #' @@ -50,16 +50,16 @@ #'DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1 #' #'@examples -#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -#'ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) -#'biasSS1 <- BiasSS(exp = exp, obs = obs, ref = ref, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) -#'biasSS2 <- BiasSS(exp = exp, obs = obs, ref = NULL, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) +#'ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) +#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50)) +#'biasSS1 <- AbsBiasSS(exp = exp, obs = obs, ref = ref, memb_dim = 'member') +#'biasSS2 <- AbsBiasSS(exp = exp, obs = obs, ref = NULL, memb_dim = 'member') #' #'@import multiApply #'@export -BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { +AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { # Check inputs ## exp, obs, and ref (1) @@ -179,19 +179,19 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na if (is.null(ref)) { ss <- multiApply::Apply(data = list(exp, obs), target_dims = time_dim, - fun = .BiasSS, na.rm = na.rm, + fun = .AbsBiasSS, na.rm = na.rm, ref = ref, ncores = ncores) } else { ss <- multiApply::Apply(data = list(exp, obs, ref), target_dims = time_dim, - fun = .BiasSS, na.rm = na.rm, + fun = .AbsBiasSS, na.rm = na.rm, ncores = ncores) } return(ss) } -.BiasSS <- function(exp, obs, ref = NULL, na.rm = FALSE) { +.AbsBiasSS <- function(exp, obs, ref = NULL, na.rm = FALSE) { # exp and obs: [sdate] # ref: [sdate] or NULL @@ -209,17 +209,17 @@ BiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na } ## Bias of the exp - bias_exp <- .Bias(exp = exp, obs = obs) + bias_exp <- .Bias(exp = exp, obs = obs, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) ## Bias of the ref if (is.null(ref)) { ## Climatological forecast ref <- mean(obs, na.rm = na.rm) } - bias_ref <- .Bias(exp = ref, obs = obs) + bias_ref <- .Bias(exp = ref, obs = obs, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) ## Skill score and significance ss <- list() - ss$biasSS <- 1 - mean(abs(bias_exp), na.rm = na.rm) / mean(abs(bias_ref), na.rm = na.rm) + ss$biasSS <- 1 - bias_exp / bias_ref ss$sign <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif return(ss) -- GitLab From eb376236ce9f8269a87f8bec364b0e5192f0f07d Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Tue, 6 Sep 2022 13:38:02 +0200 Subject: [PATCH 08/17] renamed to AbsBiasSS --- R/{BiasSS.R => AbsBiasSS.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{BiasSS.R => AbsBiasSS.R} (100%) diff --git a/R/BiasSS.R b/R/AbsBiasSS.R similarity index 100% rename from R/BiasSS.R rename to R/AbsBiasSS.R -- GitLab From 2d33532bf37b68c6701f6d74ce4565dcb4c8f87b Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Wed, 7 Sep 2022 17:14:52 +0200 Subject: [PATCH 09/17] Added dat_dim param, test-Bias and documentation --- NAMESPACE | 2 +- R/Bias.R | 103 ++++++++++++------ man/{BiasSS.Rd => AbsBiasSS.Rd} | 24 ++--- man/Bias.Rd | 33 ++++-- man/s2dv-package.Rd | 10 +- man/sampleDepthData.Rd | 6 +- man/sampleMap.Rd | 6 +- man/sampleTimeSeries.Rd | 6 +- tests/testthat/test-Bias.R | 182 ++++++++++++++++++++++++++++++++ 9 files changed, 300 insertions(+), 72 deletions(-) rename man/{BiasSS.Rd => AbsBiasSS.Rd} (81%) create mode 100644 tests/testthat/test-Bias.R diff --git a/NAMESPACE b/NAMESPACE index 6b21447..33fd850 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,11 +2,11 @@ export(ACC) export(AMV) +export(AbsBiasSS) export(AnimateMap) export(Ano) export(Ano_CrossValid) export(Bias) -export(BiasSS) export(BrierScore) export(CDORemap) export(Clim) diff --git a/R/Bias.R b/R/Bias.R index 316efc8..f7d6e31 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -13,21 +13,26 @@ #' 'dat_dim'. #'@param time_dim A character string indicating the name of the time dimension. #' The default value is 'sdate'. +#'@param dat_dim A character string indicating the name of dataset dimension. +#' The length of this dimension can be different between 'exp' and 'obs'. +#' The default value is NULL. #'@param memb_dim A character string indicating the name of the member dimension -#' to compute the ensemble mean; it should be set to NULL if the parameter 'exp' -#' is already the ensemble mean. The default value is NULL. +#' to compute the ensemble mean; it should be set to NULL if the parameter +#' 'exp' is already the ensemble mean. The default value is NULL. #'@param na.rm A logical value indicating if NAs should be removed (TRUE) or #' kept (FALSE) for computation. The default value is FALSE. -#'@param absolute A logical value indicating whether to compute the absolute bias. -#' The default value is FALSE. +#'@param absolute A logical value indicating whether to compute the absolute +#' bias. The default value is FALSE. #'@param time_dim A logical value indicating whether to compute the temporal #' mean of the bias. The default value is TRUE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' #'@return -#'A numerical array of bias with dimensions of 'exp' except 'time_dim' (if time_mean = T) -#'and 'memb_dim' dimensions. +#'A numerical array of bias with dimensions nexp, nobs and the rest dimensions +#'of 'exp' except 'time_dim' (if time_mean = T). nexp is the number of +#'experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation +#'(i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. #' #'@references #'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 @@ -40,7 +45,7 @@ #'@import multiApply #'@importFrom ClimProjDiags Subset #'@export -Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, +Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE, absolute = FALSE, time_mean = TRUE, ncores = NULL) { # Check inputs @@ -55,7 +60,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, } ## time_dim if (!is.character(time_dim) | length(time_dim) != 1) - stop('Parameter "time_dim" must be a character string.') + stop("Parameter 'time_dim' must be a character string.") if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) { stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.") } @@ -77,29 +82,29 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, } } ## dat_dim - # if (!is.null(dat_dim)) { - # if (!is.character(dat_dim) | length(dat_dim) > 1) { - # stop("Parameter 'dat_dim' must be a character string.") - # } - # if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) { - # stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.", - # " Set it as NULL if there is no dataset dimension.") - # } - # } + if (!is.null(dat_dim)) { + if (!is.character(dat_dim) | length(dat_dim) > 1) { + stop("Parameter 'dat_dim' must be a character string.") + } + if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) { + stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.", + " Set it as NULL if there is no dataset dimension.") + } + } ## exp and obs (2) name_exp <- sort(names(dim(exp))) name_obs <- sort(names(dim(obs))) if (!is.null(memb_dim)) { name_exp <- name_exp[-which(name_exp == memb_dim)] } - # if (!is.null(dat_dim)) { - # name_exp <- name_exp[-which(name_exp == dat_dim)] - # name_obs <- name_obs[-which(name_obs == dat_dim)] - # } + if (!is.null(dat_dim)) { + name_exp <- name_exp[-which(name_exp == dat_dim)] + name_obs <- name_obs[-which(name_obs == dat_dim)] + } if (!identical(length(name_exp), length(name_obs)) | !identical(dim(exp)[name_exp], dim(obs)[name_obs])) { stop(paste0("Parameter 'exp' and 'obs' must have same length of ", - "all dimensions except 'memb_dim'")) # and 'dat_dim'.")) + "all dimensions except 'memb_dim' and 'dat_dim'.")) } ## na.rm if (!is.logical(na.rm) | length(na.rm) > 1) { @@ -130,8 +135,10 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ## (Mean) Bias bias <- Apply(data = list(exp, obs), - target_dims = time_dim, - fun = .Bias, na.rm = na.rm, + target_dims = c(time_dim, dat_dim), + fun = .Bias, + dat_dim = dat_dim, + na.rm = na.rm, absolute = absolute, time_mean = time_mean, ncores = ncores)$output1 @@ -139,17 +146,43 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, return(bias) } -.Bias <- function(exp, obs, na.rm = FALSE, absolute = FALSE, time_mean = TRUE) { - - bias <- exp - obs - - if (isTRUE(absolute)){ - bias <- abs(bias) - } - - if (isTRUE(time_mean)){ - bias <- mean(bias, na.rm = na.rm) + +.Bias <- function(exp, obs, dat_dim = NULL, na.rm = FALSE, absolute = FALSE, time_mean = TRUE) { + + if (is.null(dat_dim)) { + bias <- exp - obs + + if (isTRUE(absolute)){ + bias <- abs(bias) + } + + if (isTRUE(time_mean)){ + bias <- mean(bias, na.rm = na.rm) + } + + } else { + nexp <- as.numeric(dim(exp)[dat_dim]) + nobs <- as.numeric(dim(obs)[dat_dim]) + bias <- array(dim = c(dim(exp)[1], nexp = nexp, nobs = nobs)) + + for (i in 1:nexp) { + for (j in 1:nobs) { + exp_data <- exp[ , i] + obs_data <- obs[ , j] + if (is.null(dim(exp_data))) dim(exp_data) <- c(dim(exp)[1]) + if (is.null(dim(obs_data))) dim(obs_data) <- c(dim(obs)[1]) + bias[ , i, j] <- exp_data - obs_data + } + } + + if (isTRUE(absolute)){ + bias <- abs(bias) + } + + if (isTRUE(time_mean)){ + bias <- apply(bias, c(2,3), mean, na.rm = na.rm) + } } - + return(bias) } diff --git a/man/BiasSS.Rd b/man/AbsBiasSS.Rd similarity index 81% rename from man/BiasSS.Rd rename to man/AbsBiasSS.Rd index bce496e..494982d 100644 --- a/man/BiasSS.Rd +++ b/man/AbsBiasSS.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/BiasSS.R -\name{BiasSS} -\alias{BiasSS} -\title{Compute the Mean Bias Skill Score} +% Please edit documentation in R/AbsBiasSS.R +\name{AbsBiasSS} +\alias{AbsBiasSS} +\title{Compute the Absolute Mean Bias Skill Score} \usage{ -BiasSS( +AbsBiasSS( exp, obs, ref = NULL, @@ -54,7 +54,7 @@ computation. The default value is NULL.} } } \description{ -The Mean Bias Skill Score is based on the Mean Error (Wilks, 2011) between the +The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks, 2011) between the ensemble mean forecast and the observations. It measures the accuracy of the forecast in comparison with a reference forecast to assess whether the forecast presents an improvement or a worsening with respect to that @@ -63,16 +63,16 @@ Positive values indicate that the forecast has higher skill than the reference forecast, while negative values indicate that it has a lower skill. Examples of reference forecasts are the climatological forecast (average of the observations), a previous model version, or another model. It is computed as -\code{BiasSS = 1 - Bias_exp / Bias_ref}. The statistical significance is +\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance is obtained based on a Random Walk test at the 95% confidence level (DelSole and Tippett, 2016). } \examples{ -exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) -biasSS1 <- BiasSS(exp = exp, obs = obs, ref = ref, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) -biasSS2 <- BiasSS(exp = exp, obs = obs, ref = NULL, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) +ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) +obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50)) +biasSS1 <- AbsBiasSS(exp = exp, obs = obs, ref = ref, memb_dim = 'member') +biasSS2 <- AbsBiasSS(exp = exp, obs = obs, ref = NULL, memb_dim = 'member') } \references{ diff --git a/man/Bias.Rd b/man/Bias.Rd index a58d263..fb0fafc 100644 --- a/man/Bias.Rd +++ b/man/Bias.Rd @@ -9,7 +9,10 @@ Bias( obs, time_dim = "sdate", memb_dim = NULL, + dat_dim = NULL, na.rm = FALSE, + absolute = FALSE, + time_mean = TRUE, ncores = NULL ) } @@ -21,35 +24,43 @@ dimension.} dimension. The dimensions must be the same as 'exp' except 'memb_dim' and 'dat_dim'.} -\item{time_dim}{A character string indicating the name of the time dimension. -The default value is 'sdate'.} +\item{time_dim}{A logical value indicating whether to compute the temporal +mean of the bias. The default value is TRUE.} \item{memb_dim}{A character string indicating the name of the member dimension -to compute the ensemble mean; it should be set to NULL if the parameter 'exp' -is already the ensemble mean. The default value is NULL.} +to compute the ensemble mean; it should be set to NULL if the parameter +'exp' is already the ensemble mean. The default value is NULL.} + +\item{dat_dim}{A character string indicating the name of dataset dimension. +The length of this dimension can be different between 'exp' and 'obs'. +The default value is NULL.} \item{na.rm}{A logical value indicating if NAs should be removed (TRUE) or kept (FALSE) for computation. The default value is FALSE.} +\item{absolute}{A logical value indicating whether to compute the absolute +bias. The default value is FALSE.} + \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ -A numerical array of bias with dimensions of 'exp' except 'time_dim' and -'memb_dim' dimensions. +A numerical array of bias with dimensions nexp, nobs and the rest dimensions +of 'exp' except 'time_dim' (if time_mean = T). nexp is the number of +experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation +(i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. } \description{ The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference between the ensemble mean forecast and the observations. It is a deterministic metric. Positive values indicate that the forecasts are on average too high and negative values indicate that the forecasts are on average too low. -However, it gives no information about the typical magnitude of individual -forecast errors. +It also allows to compute the Absolute Mean Bias. } \examples{ -exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50)) -obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50)) -bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1) +exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) +obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50)) +bias <- Bias(exp = exp, obs = obs, memb_dim = 'member') } \references{ diff --git a/man/s2dv-package.Rd b/man/s2dv-package.Rd index 3c98a95..ffb6783 100644 --- a/man/s2dv-package.Rd +++ b/man/s2dv-package.Rd @@ -6,7 +6,15 @@ \alias{s2dv-package} \title{s2dv: A Set of Common Tools for Seasonal to Decadal Verification} \description{ -The advanced version of package 's2dverification'. It is intended for 'seasonal to decadal' (s2d) climate forecast verification, but it can also be used in other kinds of forecasts or general climate analysis. This package is specially designed for the comparison between the experimental and observational datasets. The functionality of the included functions covers from data retrieval, data post-processing, skill scores against observation, to visualization. Compared to 's2dverification', 's2dv' is more compatible with the package 'startR', able to use multiple cores for computation and handle multi-dimensional arrays with a higher flexibility. +The advanced version of package 's2dverification'. It is + intended for 'seasonal to decadal' (s2d) climate forecast verification, but + it can also be used in other kinds of forecasts or general climate analysis. + This package is specially designed for the comparison between the experimental + and observational datasets. The functionality of the included functions covers + from data retrieval, data post-processing, skill scores against observation, + to visualization. Compared to 's2dverification', 's2dv' is more compatible + with the package 'startR', able to use multiple cores for computation and + handle multi-dimensional arrays with a higher flexibility. } \references{ \url{https://earth.bsc.es/gitlab/es/s2dv/} diff --git a/man/sampleDepthData.Rd b/man/sampleDepthData.Rd index 47e2f1b..77e4a7a 100644 --- a/man/sampleDepthData.Rd +++ b/man/sampleDepthData.Rd @@ -5,8 +5,7 @@ \alias{sampleDepthData} \title{Sample of Experimental Data for Forecast Verification In Function Of Latitudes And Depths} -\format{ -The data set provides with a variable named 'sampleDepthData'.\cr\cr +\format{The data set provides with a variable named 'sampleDepthData'.\cr\cr sampleDepthData$exp is an array that contains the experimental data and the dimension meanings and values are:\cr @@ -19,8 +18,7 @@ but in this sample is not defined (NULL).\cr\cr sampleDepthData$depths is an array with the 7 longitudes covered by the data.\cr\cr -sampleDepthData$lat is an array with the 21 latitudes covered by the data.\cr\cr -} +sampleDepthData$lat is an array with the 21 latitudes covered by the data.\cr\cr} \usage{ data(sampleDepthData) } diff --git a/man/sampleMap.Rd b/man/sampleMap.Rd index e4ec5a5..eaf8aa5 100644 --- a/man/sampleMap.Rd +++ b/man/sampleMap.Rd @@ -4,8 +4,7 @@ \name{sampleMap} \alias{sampleMap} \title{Sample Of Observational And Experimental Data For Forecast Verification In Function Of Longitudes And Latitudes} -\format{ -The data set provides with a variable named 'sampleMap'.\cr\cr +\format{The data set provides with a variable named 'sampleMap'.\cr\cr sampleMap$mod is an array that contains the experimental data and the dimension meanings and values are:\cr c(# of experimental datasets, # of members, # of starting dates, # of lead-times, # of latitudes, # of longitudes)\cr @@ -17,8 +16,7 @@ sampleMap$obs is an array that contains the observational data and the dimension sampleMap$lat is an array with the 2 latitudes covered by the data (see examples on Load() for details on why such low resolution).\cr\cr - sampleMap$lon is an array with the 3 longitudes covered by the data (see examples on Load() for details on why such low resolution). -} + sampleMap$lon is an array with the 3 longitudes covered by the data (see examples on Load() for details on why such low resolution).} \usage{ data(sampleMap) } diff --git a/man/sampleTimeSeries.Rd b/man/sampleTimeSeries.Rd index 7f058e2..05a8e79 100644 --- a/man/sampleTimeSeries.Rd +++ b/man/sampleTimeSeries.Rd @@ -4,8 +4,7 @@ \name{sampleTimeSeries} \alias{sampleTimeSeries} \title{Sample Of Observational And Experimental Data For Forecast Verification As Area Averages} -\format{ -The data set provides with a variable named 'sampleTimeSeries'.\cr\cr +\format{The data set provides with a variable named 'sampleTimeSeries'.\cr\cr sampleTimeSeries$mod is an array that contains the experimental data and the dimension meanings and values are:\cr c(# of experimental datasets, # of members, # of starting dates, # of lead-times)\cr @@ -17,8 +16,7 @@ sampleTimeSeries$obs is an array that contains the observational data and the di sampleTimeSeries$lat is an array with the 2 latitudes covered by the data that was area averaged to calculate the time series (see examples on Load() for details on why such low resolution).\cr\cr -sampleTimeSeries$lon is an array with the 3 longitudes covered by the data that was area averaged to calculate the time series (see examples on Load() for details on why such low resolution). -} +sampleTimeSeries$lon is an array with the 3 longitudes covered by the data that was area averaged to calculate the time series (see examples on Load() for details on why such low resolution).} \usage{ data(sampleTimeSeries) } diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R new file mode 100644 index 0000000..537625b --- /dev/null +++ b/tests/testthat/test-Bias.R @@ -0,0 +1,182 @@ +context("s2dv::Bias tests") + +############################################## + +# dat1 +set.seed(1) +exp1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2)) +set.seed(2) +obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2)) + +# dat2 +set.seed(1) +exp2 <- array(rnorm(40), dim = c(member = 2, sdate = 10, lat = 2)) +set.seed(2) +obs2 <- array(rnorm(20), dim = c(member = 1, sdate = 10, lat = 2)) + +# dat3 +set.seed(1) +exp3 <- array(rnorm(80), dim = c(member = 2, sdate = 10, lat = 2, dataset = 2)) +set.seed(2) +obs3 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3)) + +############################################## +test_that("1. Input checks", { + # exp and obs (1) + expect_error( + Bias(c()), + "Parameter 'exp' must be a numeric array." + ) + expect_error( + Bias(exp1, c()), + "Parameter 'obs' must be a numeric array." + ) + expect_error( + Bias(exp1, array(rnorm(20), dim = c(10, lat = 2))), + "Parameter 'exp' and 'obs' must have dimension names." + ) + # time_dim + expect_error( + Bias(exp1, obs1, time_dim = 1), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + Bias(exp1, obs1, time_dim = 'time'), + "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension." + ) + # memb_dim + expect_error( + Bias(exp1, obs1, memb_dim = TRUE), + "Parameter 'memb_dim' must be a character string." + ) + expect_error( + Bias(exp1, obs1, memb_dim = 'memb'), + "Parameter 'memb_dim' is not found in 'exp' dimension." + ) + expect_error( + Bias(exp2, array(rnorm(20), dim = c(member = 3, sdate = 10, lat = 2)), memb_dim = 'member'), + "Not implemented for observations with members ('obs' can have 'memb_dim', but it should be of length = 1).", fixed = TRUE + ) + # dat_dim + expect_error( + Bias(exp1, obs1, dat_dim = TRUE, ), + "Parameter 'dat_dim' must be a character string." + ) + expect_error( + Bias(exp1, obs1, dat_dim = 'dat_dim'), + "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension. Set it as NULL if there is no dataset dimension." + ) + # exp, ref, and obs (2) + expect_error( + Bias(exp1, array(1:9, dim = c(sdate = 9))), + "Parameter 'exp' and 'obs' must have same length of all dimensions except 'memb_dim' and 'dat_dim'." + ) + # na.rm + expect_error( + Bias(exp2, obs2, memb_dim = 'member', na.rm = 1.5), + "Parameter 'na.rm' must be one logical value." + ) + # absolute + expect_error( + Bias(exp2, obs2, memb_dim = 'member', ncores = 1.5), + "Parameter 'ncores' must be either NULL or a positive integer." + ) + # time_mean + expect_error( + Bias(exp2, obs2, memb_dim = 'member', time_mean = 1.5), + "Parameter 'time_mean' must be one logical value." + ) + # ncores + expect_error( + Bias(exp2, obs2, memb_dim = 'member', ncores = 1.5), + "Parameter 'ncores' must be either NULL or a positive integer." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + + expect_equal( + dim(Bias(exp1, obs1)), + c(lat = 2) + ) + expect_equal( + dim(Bias(exp1, obs1, time_mean = FALSE)), + c(sdate = 10, lat = 2) + ) + expect_equal( + as.vector(Bias(exp1, obs1)), + c(-0.07894886, 0.06907455), + tolerance = 0.0001 + ) + expect_equal( + as.vector(Bias(exp1, obs1, absolute = TRUE)), + c(0.9557288, 0.8169118), + tolerance = 0.0001 + ) + expect_equal( + as.vector(Bias(exp1, obs1, time_mean = FALSE, na.rm = TRUE))[1:5], + c(0.27046074, -0.00120586, -2.42347394, 2.72565648, 0.40975953), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("3. Output checks: dat2", { + + expect_equal( + dim(Bias(exp2, obs2, memb_dim = 'member')), + c(lat = 2) + ) + expect_equal( + dim(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)), + c(sdate = 10, lat = 2) + ) + expect_equal( + as.vector(Bias(exp2, obs2, memb_dim = 'member')), + c(-0.02062777, -0.18624194), + tolerance = 0.0001 + ) + expect_equal( + as.vector(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)[1:2,1:2]), + c(0.6755093, 0.1949769, 0.4329061, -1.9391461), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("4. Output checks: dat3", { + + expect_equal( + dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')), + c(nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)), + c(sdate = 10, nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset'))[5:10], + c(0.23519286, 0.18346575, -0.18624194, -0.07803352, 0.28918537, 0.39739379), + tolerance = 0.0001 + ) + expect_equal( + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', absolute = TRUE, time_mean = FALSE))[5:10], + c(0.2154482, 0.8183919, 2.1259250, 0.7796967, 1.5206510, 0.8463483), + tolerance = 0.0001 + ) + expect_equal( + as.vector(Bias(exp2, obs2, memb_dim = 'member')), + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')[1,1,]) + ) + expect_equal( + as.vector(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)), + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)[ ,1,1,]) + ) + +}) + + -- GitLab From 54d862b56d08893262023101f798666ccc7becca Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Fri, 9 Sep 2022 17:11:44 +0200 Subject: [PATCH 10/17] Added dat_dim parameter and test file --- R/AbsBiasSS.R | 238 +++++++++++++++++++++---------- R/Bias.R | 5 +- man/AbsBiasSS.Rd | 29 ++-- tests/testthat/test-AbsBiasSS.R | 240 ++++++++++++++++++++++++++++++++ 4 files changed, 424 insertions(+), 88 deletions(-) create mode 100644 tests/testthat/test-AbsBiasSS.R diff --git a/R/AbsBiasSS.R b/R/AbsBiasSS.R index 8e5da98..4585c8d 100644 --- a/R/AbsBiasSS.R +++ b/R/AbsBiasSS.R @@ -1,17 +1,17 @@ #'Compute the Absolute Mean Bias Skill Score #' -#'The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks, 2011) between the -#'ensemble mean forecast and the observations. It measures the accuracy of the -#'forecast in comparison with a reference forecast to assess whether the -#'forecast presents an improvement or a worsening with respect to that -#'reference. The Mean Bias Skill Score ranges between minus infinite and 1. +#'The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks, +#' 2011) between the ensemble mean forecast and the observations. It measures +#'the accuracy of the forecast in comparison with a reference forecast to assess +#'whether the forecast presents an improvement or a worsening with respect to +#'that reference. The Mean Bias Skill Score ranges between minus infinite and 1. #'Positive values indicate that the forecast has higher skill than the reference #'forecast, while negative values indicate that it has a lower skill. Examples #'of reference forecasts are the climatological forecast (average of the #'observations), a previous model version, or another model. It is computed as -#'\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance is -#'obtained based on a Random Walk test at the 95% confidence level (DelSole and -#'Tippett, 2016). +#'\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance +#'is obtained based on a Random Walk test at the 95% confidence level (DelSole +#'and Tippett, 2016). #' #'@param exp A named numerical array of the forecast with at least time #' dimension. @@ -30,6 +30,9 @@ #'@param memb_dim A character string indicating the name of the member dimension #' to compute the ensemble mean; it should be set to NULL if the parameter 'exp' #' and 'ref' are already the ensemble mean. The default value is NULL. +#'@param dat_dim A character string indicating the name of dataset dimension. +#' The length of this dimension can be different between 'exp' and 'obs'. +#' The default value is NULL. #'@param na.rm A logical value indicating if NAs should be removed (TRUE) or #' kept (FALSE) for computation. The default value is FALSE. #'@param ncores An integer indicating the number of cores to use for parallel @@ -37,12 +40,14 @@ #' #'@return #'\item{$biasSS}{ -#' A numerical array of BiasSS with dimensions the dimensions of -#' 'exp' except 'time_dim' and 'memb_dim' dimensions. +#' A numerical array of BiasSS with dimensions nexp, nobs and the rest +#' dimensions of 'exp' except 'time_dim' and 'memb_dim'. #'} #'\item{$sign}{ #' A logical array of the statistical significance of the BiasSS -#' with the same dimensions as $biasSS. +#' with the same dimensions as $biasSS. nexp is the number of +#' experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation +#' (i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. #'} #' #'@references @@ -58,8 +63,7 @@ #' #'@import multiApply #'@export - -AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) { +AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE, ncores = NULL) { # Check inputs ## exp, obs, and ref (1) @@ -108,49 +112,49 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, } } ## dat_dim - # if (!is.null(dat_dim)) { - # if (!is.character(dat_dim) | length(dat_dim) > 1) { - # stop("Parameter 'dat_dim' must be a character string.") - # } - # if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) { - # stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.", - # " Set it as NULL if there is no dataset dimension.") - # } - # } + if (!is.null(dat_dim)) { + if (!is.character(dat_dim) | length(dat_dim) > 1) { + stop("Parameter 'dat_dim' must be a character string.") + } + if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) { + stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.", + " Set it as NULL if there is no dataset dimension.") + } + } ## exp, obs, and ref (2) name_exp <- sort(names(dim(exp))) name_obs <- sort(names(dim(obs))) if (!is.null(memb_dim)) { name_exp <- name_exp[-which(name_exp == memb_dim)] } - # if (!is.null(dat_dim)) { - # name_exp <- name_exp[-which(name_exp == dat_dim)] - # name_obs <- name_obs[-which(name_obs == dat_dim)] - # } + if (!is.null(dat_dim)) { + name_exp <- name_exp[-which(name_exp == dat_dim)] + name_obs <- name_obs[-which(name_obs == dat_dim)] + } if (!identical(length(name_exp), length(name_obs)) | !identical(dim(exp)[name_exp], dim(obs)[name_obs])) { stop(paste0("Parameter 'exp' and 'obs' must have same length of ", - "all dimensions except 'memb_dim'")) # and 'dat_dim'.")) + "all dimensions except 'memb_dim' and 'dat_dim'.")) } if (!is.null(ref)) { name_ref <- sort(names(dim(ref))) - if (memb_dim %in% name_ref) { + if (!is.null(memb_dim) && memb_dim %in% name_ref) { name_ref <- name_ref[-which(name_ref == memb_dim)] } - # if (!is.null(dat_dim)) { - # if (dat_dim %in% name_ref) { - # if (!identical(dim(exp)[dat_dim], dim(ref)[dat_dim])) { - # stop(paste0("If parameter 'ref' has dataset dimension, it must be", - # " equal to dataset dimension of 'exp'.")) - # } - # name_ref <- name_ref[-which(name_ref == dat_dim)] - # } - # } + if (!is.null(dat_dim)) { + if (dat_dim %in% name_ref) { + if (!identical(dim(exp)[dat_dim], dim(ref)[dat_dim])) { + stop(paste0("If parameter 'ref' has dataset dimension, it must be", + " equal to dataset dimension of 'exp'.")) + } + name_ref <- name_ref[-which(name_ref == dat_dim)] + } + } if (!identical(length(name_exp), length(name_ref)) | !identical(dim(exp)[name_exp], dim(ref)[name_ref])) { stop(paste0("Parameter 'exp' and 'ref' must have the same length of ", - "all dimensions except 'memb_dim'")) #, - # " and 'dat_dim' if there is only one reference dataset.")) + "all dimensions except 'memb_dim' and 'dat_dim' if there is ", + "only one reference dataset.")) } } ## na.rm @@ -176,51 +180,135 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, } ## Mean bias skill score - if (is.null(ref)) { - ss <- multiApply::Apply(data = list(exp, obs), - target_dims = time_dim, - fun = .AbsBiasSS, na.rm = na.rm, - ref = ref, ncores = ncores) - + if (!is.null(ref)) { # use "ref" as reference forecast + if (!is.null(dat_dim) && dat_dim %in% names(dim(ref))) { + target_dims_ref <- c(time_dim, dat_dim) + } else { + target_dims_ref <- c(time_dim) + } + data <- list(exp = exp, obs = obs, ref = ref) + target_dims = list(exp = c(time_dim, dat_dim), + obs = c(time_dim, dat_dim), + ref = target_dims_ref) } else { - ss <- multiApply::Apply(data = list(exp, obs, ref), - target_dims = time_dim, - fun = .AbsBiasSS, na.rm = na.rm, - ncores = ncores) + data <- list(exp = exp, obs = obs) + target_dims = list(exp = c(time_dim, dat_dim), + obs = c(time_dim, dat_dim)) } - return(ss) + + output <- Apply(data, + target_dims = target_dims, + fun = .AbsBiasSS, + time_dim = time_dim, + dat_dim = dat_dim, + na.rm = na.rm, + ncores = ncores) + + return(output) } -.AbsBiasSS <- function(exp, obs, ref = NULL, na.rm = FALSE) { - # exp and obs: [sdate] - # ref: [sdate] or NULL +.AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE) { - if (isTRUE(na.rm)) { - if (is.null(ref)) { - good_values <- !is.na(exp) & !is.na(obs) - exp <- exp[good_values] - obs <- obs[good_values] - } else { - good_values <- !is.na(exp) & !is.na(ref) & !is.na(obs) - exp <- exp[good_values] - ref <- ref[good_values] - obs <- obs[good_values] + # exp and obs: [sdate (dat_dim)] + # ref: [sdate (dat_dim)] or NULL + + if (is.null(dat_dim)) { + if (isTRUE(na.rm)) { + if (is.null(ref)) { + good_values <- !is.na(exp) & !is.na(obs) + exp <- exp[good_values] + obs <- obs[good_values] + } else { + good_values <- !is.na(exp) & !is.na(ref) & !is.na(obs) + exp <- exp[good_values] + ref <- ref[good_values] + obs <- obs[good_values] + } + } + } else { + nexp <- as.numeric(dim(exp)[dat_dim]) + nobs <- as.numeric(dim(obs)[dat_dim]) + if (isTRUE(na.rm)) { + if (is.null(ref)) { + for (i in 1:nexp) { + for (j in 1:nobs) { + good_values <- !is.na(exp[ , i]) & !is.na(obs[ , j]) + exp[ , i] <- exp[ , i][good_values] + obs[ , j] <- obs[ , j][good_values] + } + } + } else if (length(dim(ref)) == 1) { + for (i in 1:nexp) { + for (j in 1:nobs) { + good_values <- !is.na(exp[ , i]) & !is.na(ref) & !is.na(obs[ , j]) + exp[ , i] <- exp[ , i][good_values] + obs[ , j] <- obs[ , j][good_values] + ref <- ref[good_values] + } + } + } else { + for (i in 1:nexp) { + for (j in 1:nobs) { + good_values <- !is.na(exp[ , i]) & !is.na(ref[ , i]) & !is.na(obs[ , j]) + exp[ , i] <- exp[ , i][good_values] + obs[ , j] <- obs[ , j][good_values] + ref[ , i] <- ref[ , i][good_values] + } + } + } } } - ## Bias of the exp - bias_exp <- .Bias(exp = exp, obs = obs, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + ## Bias of the exp + bias_exp <- .Bias(exp = exp, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + ## Bias of the ref - if (is.null(ref)) { ## Climatological forecast - ref <- mean(obs, na.rm = na.rm) + if (is.null(dat_dim)) { + if (is.null(ref)) { + ref <- mean(obs, na.rm = na.rm) + } + bias_ref <- .Bias(exp = ref, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + } else { + if (is.null(ref)) { + ref <- MeanDims(obs, time_dim, na.rm = na.rm) + bias_ref <- array(dim = c(nobs = nobs)) + for (j in 1:nobs) { + bias_ref[j] <- .Bias(exp = ref[j], obs = obs[ , j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + } + } else if (length(dim(ref)) == 1) { + bias_ref <- array(dim = c(nobs = nobs)) + for (j in 1:nobs) { + bias_ref[j] <- .Bias(exp = ref, obs = obs[ , j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + } + } else { + bias_ref <- .Bias(exp = ref, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + } } - bias_ref <- .Bias(exp = ref, obs = obs, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) - + ## Skill score and significance - ss <- list() - ss$biasSS <- 1 - bias_exp / bias_ref - ss$sign <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif - - return(ss) -} + if (is.null(dat_dim)) { + biasSS <- 1 - bias_exp / bias_ref + sign <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif + } else { + biasSS <- array(dim = c(nexp = nexp, nobs = nobs)) + sign <- array(dim = c(nexp = nexp, nobs = nobs)) + if (length(dim(bias_ref)) == 1) { + for (i in 1:nexp) { + for (j in 1:nobs) { + biasSS[i, j] <- 1 - bias_exp[i, j] / bias_ref[j] + sign[i, j] <- .RandomWalkTest(skill_A = bias_exp[i, j], skill_B = bias_ref[j])$signif + } + } + } else { + for (i in 1:nexp) { + for (j in 1:nobs) { + biasSS[i, j] <- 1 - bias_exp[i, j] / bias_ref[i, j] + sign[i, j] <- .RandomWalkTest(skill_A = bias_exp[i, j], skill_B = bias_ref[i, j])$signif + } + } + } + } + + return(list(biasSS = biasSS, sign = sign)) +} \ No newline at end of file diff --git a/R/Bias.R b/R/Bias.R index f7d6e31..e438115 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -137,6 +137,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, bias <- Apply(data = list(exp, obs), target_dims = c(time_dim, dat_dim), fun = .Bias, + time_dim = time_dim, dat_dim = dat_dim, na.rm = na.rm, absolute = absolute, @@ -147,7 +148,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, } -.Bias <- function(exp, obs, dat_dim = NULL, na.rm = FALSE, absolute = FALSE, time_mean = TRUE) { +.Bias <- function(exp, obs, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE, absolute = FALSE, time_mean = TRUE) { if (is.null(dat_dim)) { bias <- exp - obs @@ -180,7 +181,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, } if (isTRUE(time_mean)){ - bias <- apply(bias, c(2,3), mean, na.rm = na.rm) + bias <- MeanDims(bias, time_dim, na.rm = na.rm) } } diff --git a/man/AbsBiasSS.Rd b/man/AbsBiasSS.Rd index 494982d..89fab8a 100644 --- a/man/AbsBiasSS.Rd +++ b/man/AbsBiasSS.Rd @@ -10,6 +10,7 @@ AbsBiasSS( ref = NULL, time_dim = "sdate", memb_dim = NULL, + dat_dim = NULL, na.rm = FALSE, ncores = NULL ) @@ -37,6 +38,10 @@ The default value is 'sdate'.} to compute the ensemble mean; it should be set to NULL if the parameter 'exp' and 'ref' are already the ensemble mean. The default value is NULL.} +\item{dat_dim}{A character string indicating the name of dataset dimension. +The length of this dimension can be different between 'exp' and 'obs'. +The default value is NULL.} + \item{na.rm}{A logical value indicating if NAs should be removed (TRUE) or kept (FALSE) for computation. The default value is FALSE.} @@ -45,27 +50,29 @@ computation. The default value is NULL.} } \value{ \item{$biasSS}{ - A numerical array of BiasSS with dimensions the dimensions of - 'exp' except 'time_dim' and 'memb_dim' dimensions. + A numerical array of BiasSS with dimensions nexp, nobs and the rest + dimensions of 'exp' except 'time_dim' and 'memb_dim'. } \item{$sign}{ A logical array of the statistical significance of the BiasSS - with the same dimensions as $biasSS. + with the same dimensions as $biasSS. nexp is the number of + experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation + (i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. } } \description{ -The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks, 2011) between the -ensemble mean forecast and the observations. It measures the accuracy of the -forecast in comparison with a reference forecast to assess whether the -forecast presents an improvement or a worsening with respect to that -reference. The Mean Bias Skill Score ranges between minus infinite and 1. +The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks, +2011) between the ensemble mean forecast and the observations. It measures +the accuracy of the forecast in comparison with a reference forecast to assess +whether the forecast presents an improvement or a worsening with respect to +that reference. The Mean Bias Skill Score ranges between minus infinite and 1. Positive values indicate that the forecast has higher skill than the reference forecast, while negative values indicate that it has a lower skill. Examples of reference forecasts are the climatological forecast (average of the observations), a previous model version, or another model. It is computed as -\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance is -obtained based on a Random Walk test at the 95% confidence level (DelSole and -Tippett, 2016). +\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance +is obtained based on a Random Walk test at the 95% confidence level (DelSole +and Tippett, 2016). } \examples{ exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) diff --git a/tests/testthat/test-AbsBiasSS.R b/tests/testthat/test-AbsBiasSS.R new file mode 100644 index 0000000..2db29b5 --- /dev/null +++ b/tests/testthat/test-AbsBiasSS.R @@ -0,0 +1,240 @@ +context("s2dv::AbsBiasSS tests") + +############################################## + +# dat1 +set.seed(1) +exp1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2)) +set.seed(2) +obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2)) +set.seed(3) +ref1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2)) + +# dat2 +set.seed(1) +exp2 <- array(rnorm(60), dim = c(member = 2, sdate = 10, lat = 3)) +set.seed(2) +obs2 <- array(rnorm(30), dim = c(member = 1, sdate = 10, lat = 3)) +set.seed(3) +ref2 <- array(rnorm(60), dim = c(member = 2, sdate = 10, lat = 3)) + +# dat3 +set.seed(1) +exp3 <- array(rnorm(80), dim = c(member = 2, sdate = 10, lat = 2, dataset = 2)) +set.seed(2) +obs3 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3)) +set.seed(3) +ref3 <- array(rnorm(20), dim = c(sdate = 10, lat = 2)) + +# dat4 +set.seed(1) +exp4 <- array(rnorm(80), dim = c(member = 2, sdate = 10, lat = 2, dataset = 2)) +set.seed(2) +obs4 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3)) +set.seed(3) +ref4 <- array(rnorm(40), dim = c(sdate = 10, lat = 2, dataset = 2)) + +############################################## +test_that("1. Input checks", { + # exp and obs (1) + expect_error( + AbsBiasSS(c()), + "Parameter 'exp' must be a numeric array." + ) + expect_error( + AbsBiasSS(exp1, c()), + "Parameter 'obs' must be a numeric array." + ) + expect_error( + AbsBiasSS(exp1, array(rnorm(20), dim = c(10, lat = 2))), + "Parameter 'exp' and 'obs' must have dimension names." + ) + expect_error( + AbsBiasSS(exp1, array(rnorm(20), dim = c(10, lat = 2))), + "Parameter 'exp' and 'obs' must have dimension names." + ) + expect_error( + AbsBiasSS(exp1, obs1, ref = 'ref'), + "Parameter 'ref' must be a numeric array." + ) + expect_error( + AbsBiasSS(exp1, obs1, ref = array(rnorm(20), dim = c(10, lat = 2))), + "Parameter 'ref' must have dimension names." + ) + # time_dim + expect_error( + AbsBiasSS(exp1, obs1, time_dim = 1), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + AbsBiasSS(exp1, obs1, time_dim = 'time'), + "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension." + ) + expect_error( + AbsBiasSS(exp1, obs1, ref = array(rnorm(20), dim = c(lon = 10, lat = 2)), time_dim = 'time'), + "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension." + ) + # memb_dim + expect_error( + AbsBiasSS(exp1, obs1, memb_dim = TRUE), + "Parameter 'memb_dim' must be a character string." + ) + expect_error( + AbsBiasSS(exp1, obs1, memb_dim = 'memb'), + "Parameter 'memb_dim' is not found in 'exp' dimension." + ) + expect_error( + AbsBiasSS(exp2, array(rnorm(20), dim = c(member = 3, sdate = 10, lat = 2)), memb_dim = 'member'), + "Not implemented for observations with members ('obs' can have 'memb_dim', but it should be of length = 1).", fixed = TRUE + ) + # dat_dim + expect_error( + AbsBiasSS(exp1, obs1, dat_dim = TRUE, ), + "Parameter 'dat_dim' must be a character string." + ) + expect_error( + AbsBiasSS(exp1, obs1, dat_dim = 'dat_dim'), + "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension. Set it as NULL if there is no dataset dimension." + ) + # exp, ref, and obs (2) + expect_error( + AbsBiasSS(exp1, array(1:9, dim = c(sdate = 9))), + "Parameter 'exp' and 'obs' must have same length of all dimensions except 'memb_dim' and 'dat_dim'." + ) + expect_error( + AbsBiasSS(exp3, obs3, ref = obs3, dat_dim = 'dataset', memb_dim = 'member'), + "If parameter 'ref' has dataset dimension, it must be equal to dataset dimension of 'exp'." + ) + expect_error( + AbsBiasSS(exp1, obs1, ref = ref2), + "Parameter 'exp' and 'ref' must have the same length of all dimensions except 'memb_dim' and 'dat_dim' if there is only one reference dataset." + ) + # na.rm + expect_error( + AbsBiasSS(exp2, obs2, memb_dim = 'member', na.rm = 1.5), + "Parameter 'na.rm' must be one logical value." + ) + # ncores + expect_error( + AbsBiasSS(exp2, obs2, memb_dim = 'member', ncores = 1.5), + "Parameter 'ncores' must be either NULL or a positive integer." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + + expect_equal( + dim(AbsBiasSS(exp1, obs1)$biasSS), + c(lat = 2) + ) + expect_equal( + dim(AbsBiasSS(exp1, obs1, ref1)$biasSS), + c(lat = 2) + ) + expect_equal( + as.vector(AbsBiasSS(exp1, obs1)$biasSS), + c(-0.3103594, 0.0772921), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp1, obs1, ref1)$biasSS), + c(-0.07871642, 0.28868904), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp1, obs1)$sign), + c(FALSE, FALSE), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("3. Output checks: dat2", { + + expect_equal( + dim(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS), + c(lat = 3) + ) + expect_equal( + dim(AbsBiasSS(exp2, obs2, ref2, memb_dim = 'member')$biasSS), + c(lat = 3) + ) + expect_equal( + as.vector(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS), + c(-0.4743706, -0.2884077, -0.4064404), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp2, obs2, ref2, memb_dim = 'member')$biasSS), + c(-0.07319869, 0.06277502, -0.17321998), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("4. Output checks: dat3", { + + expect_equal( + dim(AbsBiasSS(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$biasSS), + c(nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + dim(AbsBiasSS(exp3, obs3, ref3, memb_dim = 'member', dat_dim = 'dataset')$biasSS), + c(nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + as.vector(AbsBiasSS(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$biasSS)[5:10], + c(-0.09794912, -0.11814710, -0.28840768, 0.02117376, -0.31282805, -0.08754112), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp3, obs3, ref3, memb_dim = 'member', dat_dim = 'dataset')$biasSS)[5:10], + c(0.264602089, 0.251073637, 0.006772886, 0.245427687, 0.168998894, 0.311602249), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS)[1:2], + as.vector(AbsBiasSS(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$biasSS[1,1,]) + ) + expect_equal( + as.vector(AbsBiasSS(exp3, obs3, ref3, memb_dim = 'member', dat_dim = 'dataset', na.rm = TRUE)$biasSS)[1:5], + c(-0.21373395, -0.34444526, 0.11039962, 0.05523797, 0.26460209) + ) + +}) +############################################## +test_that("5. Output checks: dat4", { + + expect_equal( + dim(AbsBiasSS(exp4, obs4, memb_dim = 'member', dat_dim = 'dataset')$biasSS), + c(nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + dim(AbsBiasSS(exp4, obs4, ref4, memb_dim = 'member', dat_dim = 'dataset')$biasSS), + c(nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + as.vector(AbsBiasSS(exp4, obs4, memb_dim = 'member', dat_dim = 'dataset', na.rm = TRUE)$biasSS)[5:10], + c(-0.09794912, -0.11814710, -0.28840768, 0.02117376, -0.31282805, -0.08754112), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp4, obs4, ref4, memb_dim = 'member', dat_dim = 'dataset')$biasSS)[5:10], + c(0.264602089, 0.133379718, 0.006772886, 0.242372951, 0.168998894, 0.272767238), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS)[1:2], + as.vector(AbsBiasSS(exp4, obs4, memb_dim = 'member', dat_dim = 'dataset')$biasSS[1,1,]) + ) + expect_equal( + as.vector(AbsBiasSS(exp4, obs4, ref4, memb_dim = 'member', dat_dim = 'dataset', na.rm = TRUE)$biasSS)[1:5], + c(-0.213733950, -0.214240924, 0.110399615, -0.009733463, 0.264602089) + ) + +}) + -- GitLab From e4cd6852d6c0ed54619621890759e9925873673d Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 14 Sep 2022 17:33:50 +0200 Subject: [PATCH 11/17] Improve document and code; add essential unit test --- R/AbsBiasSS.R | 49 +++++++++++++---------- R/Bias.R | 32 +++++++-------- man/AbsBiasSS.Rd | 3 +- man/Bias.Rd | 17 +++++--- tests/testthat/test-AbsBiasSS.R | 71 +++++++++++++++++++++++++++++++++ tests/testthat/test-Bias.R | 44 ++++++++++++++++++++ 6 files changed, 173 insertions(+), 43 deletions(-) diff --git a/R/AbsBiasSS.R b/R/AbsBiasSS.R index 4585c8d..d01b09e 100644 --- a/R/AbsBiasSS.R +++ b/R/AbsBiasSS.R @@ -11,7 +11,8 @@ #'observations), a previous model version, or another model. It is computed as #'\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance #'is obtained based on a Random Walk test at the 95% confidence level (DelSole -#'and Tippett, 2016). +#'and Tippett, 2016). If there is more than one dataset, the result will be +#'computed for each pair of exp and obs data. #' #'@param exp A named numerical array of the forecast with at least time #' dimension. @@ -63,7 +64,8 @@ #' #'@import multiApply #'@export -AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE, ncores = NULL) { +AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, + dat_dim = NULL, na.rm = FALSE, ncores = NULL) { # Check inputs ## exp, obs, and ref (1) @@ -209,8 +211,8 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, .AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE) { - # exp and obs: [sdate (dat_dim)] - # ref: [sdate (dat_dim)] or NULL + # exp and obs: [sdate, (dat_dim)] + # ref: [sdate, (dat_dim)] or NULL if (is.null(dat_dim)) { if (isTRUE(na.rm)) { @@ -232,27 +234,27 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, if (is.null(ref)) { for (i in 1:nexp) { for (j in 1:nobs) { - good_values <- !is.na(exp[ , i]) & !is.na(obs[ , j]) - exp[ , i] <- exp[ , i][good_values] - obs[ , j] <- obs[ , j][good_values] + good_values <- !is.na(exp[, i]) & !is.na(obs[, j]) + exp[, i] <- exp[, i][good_values] + obs[, j] <- obs[, j][good_values] } } - } else if (length(dim(ref)) == 1) { + } else if (length(dim(ref)) == 1) { #ref: [sdate] for (i in 1:nexp) { for (j in 1:nobs) { - good_values <- !is.na(exp[ , i]) & !is.na(ref) & !is.na(obs[ , j]) - exp[ , i] <- exp[ , i][good_values] - obs[ , j] <- obs[ , j][good_values] + good_values <- !is.na(exp[, i]) & !is.na(ref) & !is.na(obs[, j]) + exp[, i] <- exp[, i][good_values] + obs[, j] <- obs[, j][good_values] ref <- ref[good_values] } } } else { for (i in 1:nexp) { for (j in 1:nobs) { - good_values <- !is.na(exp[ , i]) & !is.na(ref[ , i]) & !is.na(obs[ , j]) - exp[ , i] <- exp[ , i][good_values] - obs[ , j] <- obs[ , j][good_values] - ref[ , i] <- ref[ , i][good_values] + good_values <- !is.na(exp[, i]) & !is.na(ref[, i]) & !is.na(obs[, j]) + exp[, i] <- exp[, i][good_values] + obs[, j] <- obs[, j][good_values] + ref[, i] <- ref[, i][good_values] } } } @@ -269,20 +271,26 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, ref <- mean(obs, na.rm = na.rm) } bias_ref <- .Bias(exp = ref, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + # bias_ref should be a number + } else { if (is.null(ref)) { ref <- MeanDims(obs, time_dim, na.rm = na.rm) + # ref should be [dat] bias_ref <- array(dim = c(nobs = nobs)) for (j in 1:nobs) { - bias_ref[j] <- .Bias(exp = ref[j], obs = obs[ , j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + bias_ref[j] <- .Bias(exp = ref[j], obs = obs[, j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) } - } else if (length(dim(ref)) == 1) { + # bias_ref should be [nobs] + } else if (length(dim(ref)) == 1) { # ref: [sdate] bias_ref <- array(dim = c(nobs = nobs)) for (j in 1:nobs) { - bias_ref[j] <- .Bias(exp = ref, obs = obs[ , j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + bias_ref[j] <- .Bias(exp = ref, obs = obs[, j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) } + # bias_ref should be [nobs] } else { bias_ref <- .Bias(exp = ref, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + # bias_ref should be [nexp, nobs] } } @@ -290,10 +298,11 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, if (is.null(dat_dim)) { biasSS <- 1 - bias_exp / bias_ref sign <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif + } else { biasSS <- array(dim = c(nexp = nexp, nobs = nobs)) sign <- array(dim = c(nexp = nexp, nobs = nobs)) - if (length(dim(bias_ref)) == 1) { + if (length(dim(bias_ref)) == 1) { # ref: [sdate] for (i in 1:nexp) { for (j in 1:nobs) { biasSS[i, j] <- 1 - bias_exp[i, j] / bias_ref[j] @@ -311,4 +320,4 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, } return(list(biasSS = biasSS, sign = sign)) -} \ No newline at end of file +} diff --git a/R/Bias.R b/R/Bias.R index e438115..0319a0f 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -4,7 +4,9 @@ #'between the ensemble mean forecast and the observations. It is a deterministic #'metric. Positive values indicate that the forecasts are on average too high #'and negative values indicate that the forecasts are on average too low. -#'It also allows to compute the Absolute Mean Bias. +#'It also allows to compute the Absolute Mean Bias or bias without temporal +#'mean. If there is more than one dataset, the result will be computed for each +#'pair of exp and obs data. #' #'@param exp A named numerical array of the forecast with at least time #' dimension. @@ -23,15 +25,15 @@ #' kept (FALSE) for computation. The default value is FALSE. #'@param absolute A logical value indicating whether to compute the absolute #' bias. The default value is FALSE. -#'@param time_dim A logical value indicating whether to compute the temporal +#'@param time_mean A logical value indicating whether to compute the temporal #' mean of the bias. The default value is TRUE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' #'@return -#'A numerical array of bias with dimensions nexp, nobs and the rest dimensions -#'of 'exp' except 'time_dim' (if time_mean = T). nexp is the number of -#'experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation +#'A numerical array of bias with dimensions c(nexp, nobs, the rest dimensions of +#''exp' except 'time_dim' (if time_mean = T) and 'memb_dim'). nexp is the number +#'of experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation #'(i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. #' #'@references @@ -148,39 +150,37 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, } -.Bias <- function(exp, obs, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE, absolute = FALSE, time_mean = TRUE) { +.Bias <- function(exp, obs, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE, + absolute = FALSE, time_mean = TRUE) { + # exp and obs: [sdate, (dat)] if (is.null(dat_dim)) { bias <- exp - obs - if (isTRUE(absolute)){ + if (isTRUE(absolute)) { bias <- abs(bias) } - if (isTRUE(time_mean)){ + if (isTRUE(time_mean)) { bias <- mean(bias, na.rm = na.rm) } } else { nexp <- as.numeric(dim(exp)[dat_dim]) nobs <- as.numeric(dim(obs)[dat_dim]) - bias <- array(dim = c(dim(exp)[1], nexp = nexp, nobs = nobs)) + bias <- array(dim = c(dim(exp)[time_dim], nexp = nexp, nobs = nobs)) for (i in 1:nexp) { for (j in 1:nobs) { - exp_data <- exp[ , i] - obs_data <- obs[ , j] - if (is.null(dim(exp_data))) dim(exp_data) <- c(dim(exp)[1]) - if (is.null(dim(obs_data))) dim(obs_data) <- c(dim(obs)[1]) - bias[ , i, j] <- exp_data - obs_data + bias[, i, j] <- exp[, i] - obs[, j] } } - if (isTRUE(absolute)){ + if (isTRUE(absolute)) { bias <- abs(bias) } - if (isTRUE(time_mean)){ + if (isTRUE(time_mean)) { bias <- MeanDims(bias, time_dim, na.rm = na.rm) } } diff --git a/man/AbsBiasSS.Rd b/man/AbsBiasSS.Rd index 89fab8a..029101d 100644 --- a/man/AbsBiasSS.Rd +++ b/man/AbsBiasSS.Rd @@ -72,7 +72,8 @@ of reference forecasts are the climatological forecast (average of the observations), a previous model version, or another model. It is computed as \code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance is obtained based on a Random Walk test at the 95% confidence level (DelSole -and Tippett, 2016). +and Tippett, 2016). If there is more than one dataset, the result will be +computed for each pair of exp and obs data. } \examples{ exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) diff --git a/man/Bias.Rd b/man/Bias.Rd index fb0fafc..2a02f2d 100644 --- a/man/Bias.Rd +++ b/man/Bias.Rd @@ -24,8 +24,8 @@ dimension.} dimension. The dimensions must be the same as 'exp' except 'memb_dim' and 'dat_dim'.} -\item{time_dim}{A logical value indicating whether to compute the temporal -mean of the bias. The default value is TRUE.} +\item{time_dim}{A character string indicating the name of the time dimension. +The default value is 'sdate'.} \item{memb_dim}{A character string indicating the name of the member dimension to compute the ensemble mean; it should be set to NULL if the parameter @@ -41,13 +41,16 @@ kept (FALSE) for computation. The default value is FALSE.} \item{absolute}{A logical value indicating whether to compute the absolute bias. The default value is FALSE.} +\item{time_mean}{A logical value indicating whether to compute the temporal +mean of the bias. The default value is TRUE.} + \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ -A numerical array of bias with dimensions nexp, nobs and the rest dimensions -of 'exp' except 'time_dim' (if time_mean = T). nexp is the number of -experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation +A numerical array of bias with dimensions c(nexp, nobs, the rest dimensions of +'exp' except 'time_dim' (if time_mean = T) and 'memb_dim'). nexp is the number +of experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation (i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. } \description{ @@ -55,7 +58,9 @@ The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference between the ensemble mean forecast and the observations. It is a deterministic metric. Positive values indicate that the forecasts are on average too high and negative values indicate that the forecasts are on average too low. -It also allows to compute the Absolute Mean Bias. +It also allows to compute the Absolute Mean Bias or bias without temporal +mean. If there is more than one dataset, the result will be computed for each +pair of exp and obs data. } \examples{ exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50)) diff --git a/tests/testthat/test-AbsBiasSS.R b/tests/testthat/test-AbsBiasSS.R index 2db29b5..6728726 100644 --- a/tests/testthat/test-AbsBiasSS.R +++ b/tests/testthat/test-AbsBiasSS.R @@ -34,6 +34,18 @@ obs4 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3)) set.seed(3) ref4 <- array(rnorm(40), dim = c(sdate = 10, lat = 2, dataset = 2)) +# dat5 +set.seed(1) +exp5 <- array(rnorm(10), dim = c(sdate = 10)) +exp5_1 <- array(exp5, dim = c(sdate = 10, dataset = 1)) +set.seed(2) +obs5 <- array(rnorm(10), dim = c(sdate = 10)) +obs5_1 <- array(obs5, dim = c(sdate = 10, dataset = 1)) +obs5_2 <- obs5_1 +obs5_2[c(1,3)] <- NA +set.seed(3) +ref5 <- array(rnorm(10), dim = c(sdate = 10)) + ############################################## test_that("1. Input checks", { # exp and obs (1) @@ -129,10 +141,18 @@ test_that("2. Output checks: dat1", { dim(AbsBiasSS(exp1, obs1)$biasSS), c(lat = 2) ) + expect_equal( + dim(AbsBiasSS(exp1, obs1)$sign), + c(lat = 2) + ) expect_equal( dim(AbsBiasSS(exp1, obs1, ref1)$biasSS), c(lat = 2) ) + expect_equal( + dim(AbsBiasSS(exp1, obs1, ref1)$sign), + c(lat = 2) + ) expect_equal( as.vector(AbsBiasSS(exp1, obs1)$biasSS), c(-0.3103594, 0.0772921), @@ -238,3 +258,54 @@ test_that("5. Output checks: dat4", { }) +############################################## +test_that("6. Output checks: dat5", { + expect_equal( + dim(AbsBiasSS(exp5, obs5)$biasSS), + NULL + ) + expect_equal( + dim(AbsBiasSS(exp5, obs5)$sign), + NULL + ) + expect_equal( + as.vector(AbsBiasSS(exp5, obs5)$biasSS), + -0.3103594, + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp5, obs5)$sign), + FALSE + ) + expect_equal( + as.vector(AbsBiasSS(exp5, obs5, ref5)$biasSS), + -0.07871642, + tolerance = 0.0001 + ) + #5_1 + expect_equal( + dim(AbsBiasSS(exp5_1, obs5_1, dat_dim = 'dataset')$biasSS), + c(nexp = 1, nobs = 1) + ) + expect_equal( + dim(AbsBiasSS(exp5_1, obs5_1, dat_dim = 'dataset')$sign), + c(nexp = 1, nobs = 1) + ) + expect_equal( + as.vector(AbsBiasSS(exp5_1, obs5_1, dat_dim = 'dataset')$biasSS), + as.vector(AbsBiasSS(exp5, obs5)$biasSS) + ) + expect_equal( + as.vector(AbsBiasSS(exp5_1, obs5_1, ref = ref5, dat_dim = 'dataset')$biasSS), + as.vector(AbsBiasSS(exp5, obs5, ref = ref5)$biasSS) + ) +# #5_2: NA test +# expect_equal( +# as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset')$biasSS), +# NA +# ) +# expect_equal( +# as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset', na.rm = T)$biasSS), +# +# ) +}) diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R index 537625b..1d81817 100644 --- a/tests/testthat/test-Bias.R +++ b/tests/testthat/test-Bias.R @@ -20,6 +20,14 @@ exp3 <- array(rnorm(80), dim = c(member = 2, sdate = 10, lat = 2, dataset = 2)) set.seed(2) obs3 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3)) +# dat4 +set.seed(1) +exp4 <- array(rnorm(10), dim = c(sdate = 10)) +exp4_1 <- array(exp4, dim = c(sdate = 10, dataset = 1)) +set.seed(2) +obs4 <- array(rnorm(10), dim = c(sdate = 10)) +obs4_1 <- array(obs4, dim = c(sdate = 10, dataset = 1)) + ############################################## test_that("1. Input checks", { # exp and obs (1) @@ -179,4 +187,40 @@ test_that("4. Output checks: dat3", { }) +############################################## +test_that("5. Output checks: dat4", { + expect_equal( + dim(Bias(exp4, obs4)), + NULL + ) + expect_equal( + dim(Bias(exp4, obs4, time_mean = F)), + c(sdate = 10) + ) + expect_equal( + as.vector(Bias(exp4, obs4, time_mean = F)), + as.vector(exp4 - obs4) + ) + expect_equal( + as.vector(Bias(exp4, obs4, time_mean = F, absolute = T)), + abs(as.vector(exp4 - obs4)) + ) + expect_equal( + as.vector(Bias(exp4, obs4, absolute = T)), + mean(abs(as.vector(exp4 - obs4))) + ) + + expect_equal( + dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset')), + c(nexp = 1, nobs = 1) + ) + expect_equal( + dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)), + c(sdate = 10, nexp = 1, nobs = 1) + ) + expect_equal( + as.vector(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)), + as.vector(Bias(exp4, obs4, time_mean = F)) + ) +}) -- GitLab From 785080990af0f97b8665411895b04b660fabc9e2 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 14 Sep 2022 17:34:01 +0200 Subject: [PATCH 12/17] Reformatting --- man/s2dv-package.Rd | 10 +--------- man/sampleDepthData.Rd | 6 ++++-- man/sampleMap.Rd | 6 ++++-- man/sampleTimeSeries.Rd | 6 ++++-- 4 files changed, 13 insertions(+), 15 deletions(-) diff --git a/man/s2dv-package.Rd b/man/s2dv-package.Rd index ffb6783..3c98a95 100644 --- a/man/s2dv-package.Rd +++ b/man/s2dv-package.Rd @@ -6,15 +6,7 @@ \alias{s2dv-package} \title{s2dv: A Set of Common Tools for Seasonal to Decadal Verification} \description{ -The advanced version of package 's2dverification'. It is - intended for 'seasonal to decadal' (s2d) climate forecast verification, but - it can also be used in other kinds of forecasts or general climate analysis. - This package is specially designed for the comparison between the experimental - and observational datasets. The functionality of the included functions covers - from data retrieval, data post-processing, skill scores against observation, - to visualization. Compared to 's2dverification', 's2dv' is more compatible - with the package 'startR', able to use multiple cores for computation and - handle multi-dimensional arrays with a higher flexibility. +The advanced version of package 's2dverification'. It is intended for 'seasonal to decadal' (s2d) climate forecast verification, but it can also be used in other kinds of forecasts or general climate analysis. This package is specially designed for the comparison between the experimental and observational datasets. The functionality of the included functions covers from data retrieval, data post-processing, skill scores against observation, to visualization. Compared to 's2dverification', 's2dv' is more compatible with the package 'startR', able to use multiple cores for computation and handle multi-dimensional arrays with a higher flexibility. } \references{ \url{https://earth.bsc.es/gitlab/es/s2dv/} diff --git a/man/sampleDepthData.Rd b/man/sampleDepthData.Rd index 77e4a7a..47e2f1b 100644 --- a/man/sampleDepthData.Rd +++ b/man/sampleDepthData.Rd @@ -5,7 +5,8 @@ \alias{sampleDepthData} \title{Sample of Experimental Data for Forecast Verification In Function Of Latitudes And Depths} -\format{The data set provides with a variable named 'sampleDepthData'.\cr\cr +\format{ +The data set provides with a variable named 'sampleDepthData'.\cr\cr sampleDepthData$exp is an array that contains the experimental data and the dimension meanings and values are:\cr @@ -18,7 +19,8 @@ but in this sample is not defined (NULL).\cr\cr sampleDepthData$depths is an array with the 7 longitudes covered by the data.\cr\cr -sampleDepthData$lat is an array with the 21 latitudes covered by the data.\cr\cr} +sampleDepthData$lat is an array with the 21 latitudes covered by the data.\cr\cr +} \usage{ data(sampleDepthData) } diff --git a/man/sampleMap.Rd b/man/sampleMap.Rd index eaf8aa5..e4ec5a5 100644 --- a/man/sampleMap.Rd +++ b/man/sampleMap.Rd @@ -4,7 +4,8 @@ \name{sampleMap} \alias{sampleMap} \title{Sample Of Observational And Experimental Data For Forecast Verification In Function Of Longitudes And Latitudes} -\format{The data set provides with a variable named 'sampleMap'.\cr\cr +\format{ +The data set provides with a variable named 'sampleMap'.\cr\cr sampleMap$mod is an array that contains the experimental data and the dimension meanings and values are:\cr c(# of experimental datasets, # of members, # of starting dates, # of lead-times, # of latitudes, # of longitudes)\cr @@ -16,7 +17,8 @@ sampleMap$obs is an array that contains the observational data and the dimension sampleMap$lat is an array with the 2 latitudes covered by the data (see examples on Load() for details on why such low resolution).\cr\cr - sampleMap$lon is an array with the 3 longitudes covered by the data (see examples on Load() for details on why such low resolution).} + sampleMap$lon is an array with the 3 longitudes covered by the data (see examples on Load() for details on why such low resolution). +} \usage{ data(sampleMap) } diff --git a/man/sampleTimeSeries.Rd b/man/sampleTimeSeries.Rd index 05a8e79..7f058e2 100644 --- a/man/sampleTimeSeries.Rd +++ b/man/sampleTimeSeries.Rd @@ -4,7 +4,8 @@ \name{sampleTimeSeries} \alias{sampleTimeSeries} \title{Sample Of Observational And Experimental Data For Forecast Verification As Area Averages} -\format{The data set provides with a variable named 'sampleTimeSeries'.\cr\cr +\format{ +The data set provides with a variable named 'sampleTimeSeries'.\cr\cr sampleTimeSeries$mod is an array that contains the experimental data and the dimension meanings and values are:\cr c(# of experimental datasets, # of members, # of starting dates, # of lead-times)\cr @@ -16,7 +17,8 @@ sampleTimeSeries$obs is an array that contains the observational data and the di sampleTimeSeries$lat is an array with the 2 latitudes covered by the data that was area averaged to calculate the time series (see examples on Load() for details on why such low resolution).\cr\cr -sampleTimeSeries$lon is an array with the 3 longitudes covered by the data that was area averaged to calculate the time series (see examples on Load() for details on why such low resolution).} +sampleTimeSeries$lon is an array with the 3 longitudes covered by the data that was area averaged to calculate the time series (see examples on Load() for details on why such low resolution). +} \usage{ data(sampleTimeSeries) } -- GitLab From 768e785696c1dde4ccffe6c44532b3add4976693 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 14 Sep 2022 17:42:37 +0200 Subject: [PATCH 13/17] Add unit test for NA --- tests/testthat/test-Bias.R | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R index 1d81817..842ecc2 100644 --- a/tests/testthat/test-Bias.R +++ b/tests/testthat/test-Bias.R @@ -27,6 +27,8 @@ exp4_1 <- array(exp4, dim = c(sdate = 10, dataset = 1)) set.seed(2) obs4 <- array(rnorm(10), dim = c(sdate = 10)) obs4_1 <- array(obs4, dim = c(sdate = 10, dataset = 1)) +obs4_2 <- obs4_1 +obs4_2[c(1, 3)] <- NA ############################################## test_that("1. Input checks", { @@ -222,5 +224,17 @@ test_that("5. Output checks: dat4", { as.vector(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)), as.vector(Bias(exp4, obs4, time_mean = F)) ) - + # 4_2: NA + expect_equal( + as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset')), + as.numeric(NA) + ) + expect_equal( + as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F))[c(1, 3)], + as.numeric(c(NA, NA)) + ) + expect_equal( + as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F))[c(2, 4:10)], + as.vector(Bias(exp4, obs4, time_mean = F))[c(2, 4:10)] + ) }) -- GitLab From 2de37faf3afbe9dd8178d9084b656548b449e94f Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Thu, 15 Sep 2022 10:54:02 +0200 Subject: [PATCH 14/17] Correct code for na.rm = T --- R/AbsBiasSS.R | 45 +-------------------------------- tests/testthat/test-AbsBiasSS.R | 18 ++++++------- 2 files changed, 10 insertions(+), 53 deletions(-) diff --git a/R/AbsBiasSS.R b/R/AbsBiasSS.R index d01b09e..effadb6 100644 --- a/R/AbsBiasSS.R +++ b/R/AbsBiasSS.R @@ -214,54 +214,11 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, # exp and obs: [sdate, (dat_dim)] # ref: [sdate, (dat_dim)] or NULL - if (is.null(dat_dim)) { - if (isTRUE(na.rm)) { - if (is.null(ref)) { - good_values <- !is.na(exp) & !is.na(obs) - exp <- exp[good_values] - obs <- obs[good_values] - } else { - good_values <- !is.na(exp) & !is.na(ref) & !is.na(obs) - exp <- exp[good_values] - ref <- ref[good_values] - obs <- obs[good_values] - } - } - } else { + if (!is.null(dat_dim)) { nexp <- as.numeric(dim(exp)[dat_dim]) nobs <- as.numeric(dim(obs)[dat_dim]) - if (isTRUE(na.rm)) { - if (is.null(ref)) { - for (i in 1:nexp) { - for (j in 1:nobs) { - good_values <- !is.na(exp[, i]) & !is.na(obs[, j]) - exp[, i] <- exp[, i][good_values] - obs[, j] <- obs[, j][good_values] - } - } - } else if (length(dim(ref)) == 1) { #ref: [sdate] - for (i in 1:nexp) { - for (j in 1:nobs) { - good_values <- !is.na(exp[, i]) & !is.na(ref) & !is.na(obs[, j]) - exp[, i] <- exp[, i][good_values] - obs[, j] <- obs[, j][good_values] - ref <- ref[good_values] - } - } - } else { - for (i in 1:nexp) { - for (j in 1:nobs) { - good_values <- !is.na(exp[, i]) & !is.na(ref[, i]) & !is.na(obs[, j]) - exp[, i] <- exp[, i][good_values] - obs[, j] <- obs[, j][good_values] - ref[, i] <- ref[, i][good_values] - } - } - } - } } - ## Bias of the exp bias_exp <- .Bias(exp = exp, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) diff --git a/tests/testthat/test-AbsBiasSS.R b/tests/testthat/test-AbsBiasSS.R index 6728726..4d67f4f 100644 --- a/tests/testthat/test-AbsBiasSS.R +++ b/tests/testthat/test-AbsBiasSS.R @@ -299,13 +299,13 @@ test_that("6. Output checks: dat5", { as.vector(AbsBiasSS(exp5_1, obs5_1, ref = ref5, dat_dim = 'dataset')$biasSS), as.vector(AbsBiasSS(exp5, obs5, ref = ref5)$biasSS) ) -# #5_2: NA test -# expect_equal( -# as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset')$biasSS), -# NA -# ) -# expect_equal( -# as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset', na.rm = T)$biasSS), -# -# ) + #5_2: NA test + expect_equal( + as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset')$biasSS), + as.numeric(NA) + ) + expect_equal( + as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset', na.rm = T)$biasSS), + c(-0.4636772) + ) }) -- GitLab From 0f01db04200bbacfcaebbfbd63aa2c9fff583cde Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Thu, 15 Sep 2022 15:56:04 +0200 Subject: [PATCH 15/17] Correct case na.rm = T for all possibilities --- R/AbsBiasSS.R | 102 +++++++++++++++----------------- tests/testthat/test-AbsBiasSS.R | 10 +++- 2 files changed, 58 insertions(+), 54 deletions(-) diff --git a/R/AbsBiasSS.R b/R/AbsBiasSS.R index effadb6..b5dbe77 100644 --- a/R/AbsBiasSS.R +++ b/R/AbsBiasSS.R @@ -201,7 +201,6 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, output <- Apply(data, target_dims = target_dims, fun = .AbsBiasSS, - time_dim = time_dim, dat_dim = dat_dim, na.rm = na.rm, ncores = ncores) @@ -209,72 +208,69 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, return(output) } -.AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE) { +.AbsBiasSS <- function(exp, obs, ref = NULL, dat_dim = NULL, na.rm = FALSE) { + # exp and obs: [sdate] + # ref: [sdate] or NULL - # exp and obs: [sdate, (dat_dim)] - # ref: [sdate, (dat_dim)] or NULL - - if (!is.null(dat_dim)) { - nexp <- as.numeric(dim(exp)[dat_dim]) - nobs <- as.numeric(dim(obs)[dat_dim]) - } - - ## Bias of the exp - bias_exp <- .Bias(exp = exp, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) - - ## Bias of the ref if (is.null(dat_dim)) { - if (is.null(ref)) { - ref <- mean(obs, na.rm = na.rm) + nexp <- 1 + nobs <- 1 + exp <- InsertDim(exp, posdim = 2, lendim = 1, name = 'dataset') + obs <- InsertDim(obs, posdim = 2, lendim = 1, name = 'dataset') + if (!is.null(ref)) { + ref <- InsertDim(ref, posdim = 2, lendim = 1, name = 'dataset') } - bias_ref <- .Bias(exp = ref, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) - # bias_ref should be a number - } else { - if (is.null(ref)) { - ref <- MeanDims(obs, time_dim, na.rm = na.rm) - # ref should be [dat] - bias_ref <- array(dim = c(nobs = nobs)) - for (j in 1:nobs) { - bias_ref[j] <- .Bias(exp = ref[j], obs = obs[, j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) - } - # bias_ref should be [nobs] - } else if (length(dim(ref)) == 1) { # ref: [sdate] - bias_ref <- array(dim = c(nobs = nobs)) - for (j in 1:nobs) { - bias_ref[j] <- .Bias(exp = ref, obs = obs[, j], dat_dim = NULL, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) - } - # bias_ref should be [nobs] + nexp <- as.numeric(dim(exp)[dat_dim]) + nobs <- as.numeric(dim(obs)[dat_dim]) + if (length(dim(ref)) == 1) { # ref: [sdate] + ref_dat_dim <- FALSE } else { - bias_ref <- .Bias(exp = ref, obs = obs, dat_dim = dat_dim, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) - # bias_ref should be [nexp, nobs] + ref_dat_dim <- TRUE } } - ## Skill score and significance - if (is.null(dat_dim)) { - biasSS <- 1 - bias_exp / bias_ref - sign <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif + biasSS <- array(dim = c(nexp = nexp, nobs = nobs)) + sign <- array(dim = c(nexp = nexp, nobs = nobs)) - } else { - biasSS <- array(dim = c(nexp = nexp, nobs = nobs)) - sign <- array(dim = c(nexp = nexp, nobs = nobs)) - if (length(dim(bias_ref)) == 1) { # ref: [sdate] - for (i in 1:nexp) { - for (j in 1:nobs) { - biasSS[i, j] <- 1 - bias_exp[i, j] / bias_ref[j] - sign[i, j] <- .RandomWalkTest(skill_A = bias_exp[i, j], skill_B = bias_ref[j])$signif - } - } + for (i in 1:nexp) { + exp_data <- exp[, i] + if (isTRUE(ref_dat_dim)) { + ref_data <- ref[, i] } else { - for (i in 1:nexp) { - for (j in 1:nobs) { - biasSS[i, j] <- 1 - bias_exp[i, j] / bias_ref[i, j] - sign[i, j] <- .RandomWalkTest(skill_A = bias_exp[i, j], skill_B = bias_ref[i, j])$signif + ref_data <- ref + } + for (j in 1:nobs) { + obs_data <- obs[, j] + if (isTRUE(na.rm)) { + if (is.null(ref)) { + good_values <- !is.na(exp_data) & !is.na(obs_data) + exp_data <- exp_data[good_values] + obs_data <- obs_data[good_values] + } else { + good_values <- !is.na(exp_data) & !is.na(ref_data) & !is.na(obs_data) + exp_data <- exp_data[good_values] + ref_data <- ref_data[good_values] + obs_data <- obs_data[good_values] } } + ## Bias of the exp + bias_exp <- .Bias(exp = exp_data, obs = obs_data, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + ## Bias of the ref + if (is.null(ref)) { ## Climatological forecast + ref_data <- mean(obs_data, na.rm = na.rm) + } + bias_ref <- .Bias(exp = ref_data, obs = obs_data, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + biasSS[i, j] <- 1 - bias_exp / bias_ref + sign[i, j] <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif } } + if (is.null(dat_dim)) { + dim(biasSS) <- NULL + dim(sign) <- NULL + } + + return(list(biasSS = biasSS, sign = sign)) } diff --git a/tests/testthat/test-AbsBiasSS.R b/tests/testthat/test-AbsBiasSS.R index 4d67f4f..a1fff9c 100644 --- a/tests/testthat/test-AbsBiasSS.R +++ b/tests/testthat/test-AbsBiasSS.R @@ -38,6 +38,8 @@ ref4 <- array(rnorm(40), dim = c(sdate = 10, lat = 2, dataset = 2)) set.seed(1) exp5 <- array(rnorm(10), dim = c(sdate = 10)) exp5_1 <- array(exp5, dim = c(sdate = 10, dataset = 1)) +exp5_2 <- exp5_1 +exp5_2[c(1,1)] <- NA set.seed(2) obs5 <- array(rnorm(10), dim = c(sdate = 10)) obs5_1 <- array(obs5, dim = c(sdate = 10, dataset = 1)) @@ -306,6 +308,12 @@ test_that("6. Output checks: dat5", { ) expect_equal( as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset', na.rm = T)$biasSS), - c(-0.4636772) + c(-0.4636772), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp5_2, obs5_2, ref = ref, dat_dim = 'dataset', na.rm = T)$biasSS), + c(-0.2542055), + tolerance = 0.0001 ) }) -- GitLab From 82c582393e763407199f0eacbf7c816e71987932 Mon Sep 17 00:00:00 2001 From: Eva Rifa Date: Thu, 15 Sep 2022 16:14:05 +0200 Subject: [PATCH 16/17] Fix pipeline --- R/AbsBiasSS.R | 2 ++ tests/testthat/test-AbsBiasSS.R | 30 +++++++++++++++--------------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/R/AbsBiasSS.R b/R/AbsBiasSS.R index b5dbe77..7c9e06d 100644 --- a/R/AbsBiasSS.R +++ b/R/AbsBiasSS.R @@ -220,6 +220,7 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, if (!is.null(ref)) { ref <- InsertDim(ref, posdim = 2, lendim = 1, name = 'dataset') } + ref_dat_dim <- FALSE } else { nexp <- as.numeric(dim(exp)[dat_dim]) nobs <- as.numeric(dim(obs)[dat_dim]) @@ -261,6 +262,7 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, ref_data <- mean(obs_data, na.rm = na.rm) } bias_ref <- .Bias(exp = ref_data, obs = obs_data, na.rm = na.rm, absolute = TRUE, time_mean = TRUE) + ## Skill score and significance biasSS[i, j] <- 1 - bias_exp / bias_ref sign[i, j] <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif } diff --git a/tests/testthat/test-AbsBiasSS.R b/tests/testthat/test-AbsBiasSS.R index a1fff9c..cfeeb09 100644 --- a/tests/testthat/test-AbsBiasSS.R +++ b/tests/testthat/test-AbsBiasSS.R @@ -301,19 +301,19 @@ test_that("6. Output checks: dat5", { as.vector(AbsBiasSS(exp5_1, obs5_1, ref = ref5, dat_dim = 'dataset')$biasSS), as.vector(AbsBiasSS(exp5, obs5, ref = ref5)$biasSS) ) - #5_2: NA test - expect_equal( - as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset')$biasSS), - as.numeric(NA) - ) - expect_equal( - as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset', na.rm = T)$biasSS), - c(-0.4636772), - tolerance = 0.0001 - ) - expect_equal( - as.vector(AbsBiasSS(exp5_2, obs5_2, ref = ref, dat_dim = 'dataset', na.rm = T)$biasSS), - c(-0.2542055), - tolerance = 0.0001 - ) + #5_2: NA test + expect_equal( + as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset')$biasSS), + as.numeric(NA) + ) + expect_equal( + as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset', na.rm = T)$biasSS), + c(-0.4636772), + tolerance = 0.0001 + ) + expect_equal( + as.vector(AbsBiasSS(exp5_2, obs5_2, ref = ref5, dat_dim = 'dataset', na.rm = T)$biasSS), + c(0.08069355), + tolerance = 0.0001 + ) }) -- GitLab From 49a2d17c004dbe9f5bab91affa84b9c20911aed2 Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 19 Sep 2022 09:40:50 +0200 Subject: [PATCH 17/17] Trivial format and comment adjustment --- R/AbsBiasSS.R | 6 ++++-- tests/testthat/test-AbsBiasSS.R | 4 ++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/R/AbsBiasSS.R b/R/AbsBiasSS.R index 7c9e06d..0838f25 100644 --- a/R/AbsBiasSS.R +++ b/R/AbsBiasSS.R @@ -209,9 +209,10 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, } .AbsBiasSS <- function(exp, obs, ref = NULL, dat_dim = NULL, na.rm = FALSE) { - # exp and obs: [sdate] - # ref: [sdate] or NULL + # exp and obs: [sdate, (dat_dim)] + # ref: [sdate, (dat_dim)] or NULL + # Adjust exp, obs, ref to have dat_dim temporarily if (is.null(dat_dim)) { nexp <- 1 nobs <- 1 @@ -243,6 +244,7 @@ AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, } for (j in 1:nobs) { obs_data <- obs[, j] + if (isTRUE(na.rm)) { if (is.null(ref)) { good_values <- !is.na(exp_data) & !is.na(obs_data) diff --git a/tests/testthat/test-AbsBiasSS.R b/tests/testthat/test-AbsBiasSS.R index cfeeb09..b2e70f3 100644 --- a/tests/testthat/test-AbsBiasSS.R +++ b/tests/testthat/test-AbsBiasSS.R @@ -39,12 +39,12 @@ set.seed(1) exp5 <- array(rnorm(10), dim = c(sdate = 10)) exp5_1 <- array(exp5, dim = c(sdate = 10, dataset = 1)) exp5_2 <- exp5_1 -exp5_2[c(1,1)] <- NA +exp5_2[1] <- NA set.seed(2) obs5 <- array(rnorm(10), dim = c(sdate = 10)) obs5_1 <- array(obs5, dim = c(sdate = 10, dataset = 1)) obs5_2 <- obs5_1 -obs5_2[c(1,3)] <- NA +obs5_2[c(1, 3)] <- NA set.seed(3) ref5 <- array(rnorm(10), dim = c(sdate = 10)) -- GitLab