From b1f57140beed23fe8ee9b6db977f2e531f720ab0 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 22 May 2024 15:02:02 +0200 Subject: [PATCH 01/11] Bias Welch significance test --- R/Bias.R | 45 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 36 insertions(+), 9 deletions(-) diff --git a/R/Bias.R b/R/Bias.R index 5ec8d0f..ddf7b18 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -27,6 +27,7 @@ #' bias. The default value is FALSE. #'@param time_mean A logical value indicating whether to compute the temporal #' mean of the bias. The default value is TRUE. +#'@param alpha A numeric or NULL (default) to indicate the significance level using Weltch test. Only available when absolute is FALSE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -48,7 +49,7 @@ #'@importFrom ClimProjDiags Subset #'@export Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE, - absolute = FALSE, time_mean = TRUE, ncores = NULL) { + absolute = FALSE, time_mean = TRUE, alpha = NULL, ncores = NULL) { # Check inputs ## exp and obs (1) @@ -120,6 +121,12 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, if (!is.logical(time_mean) | length(time_mean) > 1) { stop("Parameter 'time_mean' must be one logical value.") } + ## alpha + if (!is.null(alpha)) { + if (!is.numeric(alpha) | length(alpha) > 1) { + stop("Parameter 'alpha' must be null or a numeric value.") + } + } ## ncores if (!is.null(ncores)) { if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | @@ -144,16 +151,19 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = na.rm, absolute = absolute, time_mean = time_mean, - ncores = ncores)$output1 - + alpha = alpha, + ncores = ncores) + + if (is.null(alpha)) { + bias <- bias$output1 + } return(bias) } .Bias <- function(exp, obs, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE, - absolute = FALSE, time_mean = TRUE) { + absolute = FALSE, time_mean = TRUE, alpha = NULL) { # exp and obs: [sdate, (dat)] - if (is.null(dat_dim)) { bias <- exp - obs @@ -164,15 +174,29 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, if (isTRUE(time_mean)) { bias <- mean(bias, na.rm = na.rm) } - + + if (!is.null(alpha)) { + if (!absolute) { + pval <- t.test(x = obs, y = exp, alternative = "two.sided")$p.value + sig <- pval <= alpha + } + } } else { nexp <- as.numeric(dim(exp)[dat_dim]) nobs <- as.numeric(dim(obs)[dat_dim]) bias <- array(dim = c(dim(exp)[time_dim], nexp = nexp, nobs = nobs)) - + pval <- array(dim = c(nexp = nexp, nobs = nobs)) + sig <- array(dim = c(nexp = nexp, nobs = nobs)) for (i in 1:nexp) { for (j in 1:nobs) { bias[, i, j] <- exp[, i] - obs[, j] + if (!is.null(alpha)) { + if (!absolute) { + pval[i,j] <- t.test(x = obs[,j], y = exp[,i], + alternative = "two.sided")$p.value + sig[i,j] <- pval <= alpha + } + } } } @@ -183,7 +207,10 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, if (isTRUE(time_mean)) { bias <- MeanDims(bias, time_dim, na.rm = na.rm) } + } + if (!is.null(alpha) && !absolute) { + return(list(bias = bias, sig = sig)) + } else { + return(bias) } - - return(bias) } -- GitLab From af589d3fbf8dd2e966d816da5c6a9232bbbaeddf Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 22 May 2024 15:28:32 +0200 Subject: [PATCH 02/11] describe return --- R/Bias.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Bias.R b/R/Bias.R index ddf7b18..b9292ca 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -35,7 +35,7 @@ #'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. +#'(i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. If alpha is specified, and absolute is FALSE, the result is a list with two elements, the bias as describe above and the significance as logical array with the same dimensions. #' #'@references #'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 -- GitLab From 6ac387e5afbf440acb3bdca952ef9294334307f6 Mon Sep 17 00:00:00 2001 From: vagudets Date: Thu, 13 Jun 2024 14:40:23 +0200 Subject: [PATCH 03/11] Rename sig to sign, minor corrections to the documentation --- R/Bias.R | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/R/Bias.R b/R/Bias.R index b9292ca..1763a7b 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -27,7 +27,8 @@ #' bias. The default value is FALSE. #'@param time_mean A logical value indicating whether to compute the temporal #' mean of the bias. The default value is TRUE. -#'@param alpha A numeric or NULL (default) to indicate the significance level using Weltch test. Only available when absolute is FALSE. +#'@param alpha A numeric or NULL (default) to indicate the significance level +#' using Welch's t-test. Only available when absolute is FALSE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -35,7 +36,10 @@ #'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. If alpha is specified, and absolute is FALSE, the result is a list with two elements, the bias as describe above and the significance as logical array with the same dimensions. +#'(i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. If +#'alpha is specified, and absolute is FALSE, the result is a list with two +#'elements: the bias as described above and the significance as a logical array +#'with the same dimensions. #' #'@references #'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7 @@ -48,8 +52,9 @@ #'@import multiApply #'@importFrom ClimProjDiags Subset #'@export -Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE, - absolute = FALSE, time_mean = TRUE, alpha = NULL, ncores = NULL) { +Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, + na.rm = FALSE, absolute = FALSE, time_mean = TRUE, + alpha = NULL, ncores = NULL) { # Check inputs ## exp and obs (1) @@ -178,7 +183,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, if (!is.null(alpha)) { if (!absolute) { pval <- t.test(x = obs, y = exp, alternative = "two.sided")$p.value - sig <- pval <= alpha + sign <- pval <= alpha } } } else { @@ -186,7 +191,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, nobs <- as.numeric(dim(obs)[dat_dim]) bias <- array(dim = c(dim(exp)[time_dim], nexp = nexp, nobs = nobs)) pval <- array(dim = c(nexp = nexp, nobs = nobs)) - sig <- array(dim = c(nexp = nexp, nobs = nobs)) + sign <- array(dim = c(nexp = nexp, nobs = nobs)) for (i in 1:nexp) { for (j in 1:nobs) { bias[, i, j] <- exp[, i] - obs[, j] @@ -194,7 +199,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, if (!absolute) { pval[i,j] <- t.test(x = obs[,j], y = exp[,i], alternative = "two.sided")$p.value - sig[i,j] <- pval <= alpha + sign[i,j] <- pval <= alpha } } } @@ -209,7 +214,7 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, } } if (!is.null(alpha) && !absolute) { - return(list(bias = bias, sig = sig)) + return(list(bias = bias, sign = sign)) } else { return(bias) } -- GitLab From 8b0d81e5fd5db460ea0d3669cd77fc1323e6f286 Mon Sep 17 00:00:00 2001 From: vagudets Date: Thu, 13 Jun 2024 16:37:53 +0200 Subject: [PATCH 04/11] =?UTF-8?q?Add=20example=C2=B7?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- R/Bias.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/Bias.R b/R/Bias.R index 1763a7b..c1a0336 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -48,6 +48,7 @@ #'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') +#'bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.05) #' #'@import multiApply #'@importFrom ClimProjDiags Subset -- GitLab From 8650eb334190a511836effc103d296abfb7e26ec Mon Sep 17 00:00:00 2001 From: vagudets Date: Thu, 13 Jun 2024 16:39:36 +0200 Subject: [PATCH 05/11] Modify .Rd file --- man/Bias.Rd | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/man/Bias.Rd b/man/Bias.Rd index 2a02f2d..4ebd761 100644 --- a/man/Bias.Rd +++ b/man/Bias.Rd @@ -13,6 +13,7 @@ Bias( na.rm = FALSE, absolute = FALSE, time_mean = TRUE, + alpha = NULL, ncores = NULL ) } @@ -44,6 +45,9 @@ 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{alpha}{A numeric or NULL (default) to indicate the significance level +using Welch's t-test. Only available when absolute is FALSE.} + \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } @@ -51,7 +55,10 @@ computation. The default value is NULL.} 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. +(i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted. If +alpha is specified, and absolute is FALSE, the result is a list with two +elements: the bias as described above and the significance as a logical array +with the same dimensions. } \description{ The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference @@ -66,6 +73,7 @@ pair of exp and obs data. 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') +bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.05) } \references{ -- GitLab From ef1cc11b5a385c3c5f52aa2417731c9c0ae5b532 Mon Sep 17 00:00:00 2001 From: vagudets Date: Tue, 18 Jun 2024 11:30:18 +0200 Subject: [PATCH 06/11] Set alpha = 0.05 as default, fix bug in significance computation, update unit tests --- R/Bias.R | 18 +++++--- man/Bias.Rd | 5 +- tests/testthat/test-Bias.R | 95 +++++++++++++++++++++++++------------- 3 files changed, 79 insertions(+), 39 deletions(-) diff --git a/R/Bias.R b/R/Bias.R index c1a0336..0f0a97e 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -48,14 +48,15 @@ #'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') -#'bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.05) +#'bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.01) +#'abs_bias <- Bias(exp = exp, obs = obs, absolute = TRUE, alpha = NULL) #' #'@import multiApply #'@importFrom ClimProjDiags Subset #'@export Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE, absolute = FALSE, time_mean = TRUE, - alpha = NULL, ncores = NULL) { + alpha = 0.05, ncores = NULL) { # Check inputs ## exp and obs (1) @@ -129,9 +130,14 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, } ## alpha if (!is.null(alpha)) { - if (!is.numeric(alpha) | length(alpha) > 1) { + if (any(!is.numeric(alpha) | alpha <= 0 | alpha >= 1 | length(alpha) > 1)) { stop("Parameter 'alpha' must be null or a numeric value.") } + if (absolute) { + alpha <- NULL + .warning("Parameter 'absolute' is TRUE, so 'alpha' has been set to", + "false and significance will not be returned.") + } } ## ncores if (!is.null(ncores)) { @@ -198,9 +204,9 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, bias[, i, j] <- exp[, i] - obs[, j] if (!is.null(alpha)) { if (!absolute) { - pval[i,j] <- t.test(x = obs[,j], y = exp[,i], - alternative = "two.sided")$p.value - sign[i,j] <- pval <= alpha + pval[i, j] <- t.test(x = obs[, j], y = exp[, i], + alternative = "two.sided")$p.value + sign[i, j] <- pval[i, j] <= alpha } } } diff --git a/man/Bias.Rd b/man/Bias.Rd index 4ebd761..a003d10 100644 --- a/man/Bias.Rd +++ b/man/Bias.Rd @@ -13,7 +13,7 @@ Bias( na.rm = FALSE, absolute = FALSE, time_mean = TRUE, - alpha = NULL, + alpha = 0.05, ncores = NULL ) } @@ -73,7 +73,8 @@ pair of exp and obs data. 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') -bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.05) +bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.01) +abs_bias <- Bias(exp = exp, obs = obs, absolute = TRUE, alpha = NULL) } \references{ diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R index 4c6cc99..490c275 100644 --- a/tests/testthat/test-Bias.R +++ b/tests/testthat/test-Bias.R @@ -67,7 +67,7 @@ test_that("1. Input checks", { ) # dat_dim expect_error( - Bias(exp1, obs1, dat_dim = TRUE, ), + Bias(exp1, obs1, dat_dim = TRUE), "Parameter 'dat_dim' must be a character string." ) expect_error( @@ -94,6 +94,11 @@ test_that("1. Input checks", { Bias(exp2, obs2, memb_dim = 'member', time_mean = 1.5), "Parameter 'time_mean' must be one logical value." ) + # alpha + expect_error( + Bias(exp1, obs1, alpha = TRUE), + "Parameter 'alpha' must be null or a numeric value." + ) # ncores expect_error( Bias(exp2, obs2, memb_dim = 'member', ncores = 1.5), @@ -106,25 +111,33 @@ test_that("1. Input checks", { test_that("2. Output checks: dat1", { expect_equal( - dim(Bias(exp1, obs1)), + dim(Bias(exp1, obs1)$bias), + c(lat = 2) + ) + expect_equal( + dim(Bias(exp1, obs1)$sign), c(lat = 2) ) expect_equal( - dim(Bias(exp1, obs1, time_mean = FALSE)), + dim(Bias(exp1, obs1, time_mean = FALSE)$bias), c(sdate = 10, lat = 2) ) expect_equal( - as.vector(Bias(exp1, obs1)), + as.vector(Bias(exp1, obs1)$bias), c(-0.07894886, 0.06907455), tolerance = 0.0001 ) expect_equal( - as.vector(Bias(exp1, obs1, absolute = TRUE)), + as.vector(Bias(exp1, obs1)$bias), + c(FALSE, FALSE), + ) + expect_equal( + as.vector(Bias(exp1, obs1, absolute = TRUE, alpha = NULL)), c(0.9557288, 0.8169118), tolerance = 0.0001 ) expect_equal( - as.vector(Bias(exp1, obs1, time_mean = FALSE, na.rm = TRUE))[1:5], + as.vector(Bias(exp1, obs1, time_mean = FALSE, na.rm = TRUE))$bias[1:5], c(0.27046074, -0.00120586, -2.42347394, 2.72565648, 0.40975953), tolerance = 0.0001 ) @@ -135,20 +148,24 @@ test_that("2. Output checks: dat1", { test_that("3. Output checks: dat2", { expect_equal( - dim(Bias(exp2, obs2, memb_dim = 'member')), + dim(Bias(exp2, obs2, memb_dim = 'member')$bias), c(lat = 2) ) expect_equal( - dim(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)), + dim(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)$bias), c(sdate = 10, lat = 2) ) expect_equal( - as.vector(Bias(exp2, obs2, memb_dim = 'member')), + dim(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)$sign), + c(lat = 2) + ) + expect_equal( + as.vector(Bias(exp2, obs2, memb_dim = 'member')$bias), 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]), + as.vector(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)$bias[1:2,1:2]), c(0.6755093, 0.1949769, 0.4329061, -1.9391461), tolerance = 0.0001 ) @@ -159,30 +176,39 @@ test_that("3. Output checks: dat2", { test_that("4. Output checks: dat3", { expect_equal( - dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')), + dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$bias), c(nexp = 2, nobs = 3, lat = 2) ) expect_equal( - dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)), + dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$sign), + c(nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)$bias), c(sdate = 10, nexp = 2, nobs = 3, lat = 2) ) expect_equal( - as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset'))[5:10], + dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)$sign), + c(nexp = 2, nobs = 3, lat = 2) + ) + expect_equal( + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$bias)[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], + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', + absolute = TRUE, alpha = NULL, 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,]) + as.vector(Bias(exp2, obs2, memb_dim = 'member')$bias), + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$bias[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,]) + as.vector(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)$bias), + as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)$bias[ ,1,1,]) ) }) @@ -190,49 +216,56 @@ test_that("4. Output checks: dat3", { ############################################## test_that("5. Output checks: dat4", { expect_equal( - dim(Bias(exp4, obs4)), + dim(Bias(exp4, obs4)$bias), NULL ) expect_equal( - dim(Bias(exp4, obs4, time_mean = F)), + dim(Bias(exp4, obs4, time_mean = F)$bias), c(sdate = 10) ) expect_equal( - as.vector(Bias(exp4, obs4, time_mean = F)), + dim(Bias(exp4, obs4, time_mean = F)$sign), + dim(Bias(exp4, obs4)$sign), + NULL + ) + expect_equal( + as.vector(Bias(exp4, obs4, time_mean = F)$bias), as.vector(exp4 - obs4) ) expect_equal( - as.vector(Bias(exp4, obs4, time_mean = F, absolute = T)), + as.vector(Bias(exp4, obs4, time_mean = F, absolute = T, alpha = NULL)), abs(as.vector(exp4 - obs4)) ) expect_equal( - as.vector(Bias(exp4, obs4, absolute = T)), + as.vector(Bias(exp4, obs4, absolute = T, alpha = NULL)), mean(abs(as.vector(exp4 - obs4))) ) - + expect_equal( - dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset')), + dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset')$bias), + dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset')$sign), c(nexp = 1, nobs = 1) ) expect_equal( - dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)), + dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)$bias), 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)) + as.vector(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)$bias), + as.vector(Bias(exp4, obs4, time_mean = F)$bias) ) # 4_2: NA + ## TODO: Significance expect_equal( - as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset')), + as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset')$bias), as.numeric(NA) ) expect_equal( - as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F))[c(1, 3)], + as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F)$bias)[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_1, obs4_2, dat_dim = 'dataset', time_mean = F)$bias)[c(2, 4:10)], as.vector(Bias(exp4, obs4, time_mean = F))[c(2, 4:10)] ) }) -- GitLab From ef5596054b3d7503f56ccd0557f5d0b50660d8ec Mon Sep 17 00:00:00 2001 From: vagudets Date: Tue, 18 Jun 2024 11:40:19 +0200 Subject: [PATCH 07/11] Add memb_dim to abs_bias example --- R/Bias.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Bias.R b/R/Bias.R index 0f0a97e..7b303de 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -49,7 +49,7 @@ #'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50)) #'bias <- Bias(exp = exp, obs = obs, memb_dim = 'member') #'bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.01) -#'abs_bias <- Bias(exp = exp, obs = obs, absolute = TRUE, alpha = NULL) +#'abs_bias <- Bias(exp = exp, obs = obs, memb_dim = 'member', absolute = TRUE, alpha = NULL) #' #'@import multiApply #'@importFrom ClimProjDiags Subset -- GitLab From 6e732beccd4d90d5726fbdd1487511f7c3e279ad Mon Sep 17 00:00:00 2001 From: vagudets Date: Tue, 18 Jun 2024 11:41:03 +0200 Subject: [PATCH 08/11] Update Rd file --- man/Bias.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/Bias.Rd b/man/Bias.Rd index a003d10..e94beb5 100644 --- a/man/Bias.Rd +++ b/man/Bias.Rd @@ -74,7 +74,7 @@ exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50)) bias <- Bias(exp = exp, obs = obs, memb_dim = 'member') bias2 <- Bias(exp = exp, obs = obs, memb_dim = 'member', alpha = 0.01) -abs_bias <- Bias(exp = exp, obs = obs, absolute = TRUE, alpha = NULL) +abs_bias <- Bias(exp = exp, obs = obs, memb_dim = 'member', absolute = TRUE, alpha = NULL) } \references{ -- GitLab From c6d4b72bdc9b60cf0a0113b9ce6c83a79e74ae88 Mon Sep 17 00:00:00 2001 From: vagudets Date: Tue, 18 Jun 2024 12:07:25 +0200 Subject: [PATCH 09/11] Fix unit tests --- tests/testthat/test-Bias.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R index 490c275..4024897 100644 --- a/tests/testthat/test-Bias.R +++ b/tests/testthat/test-Bias.R @@ -128,7 +128,7 @@ test_that("2. Output checks: dat1", { tolerance = 0.0001 ) expect_equal( - as.vector(Bias(exp1, obs1)$bias), + as.vector(Bias(exp1, obs1)$sign), c(FALSE, FALSE), ) expect_equal( @@ -266,6 +266,6 @@ test_that("5. Output checks: dat4", { ) expect_equal( as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F)$bias)[c(2, 4:10)], - as.vector(Bias(exp4, obs4, time_mean = F))[c(2, 4:10)] + as.vector(Bias(exp4, obs4, time_mean = F)$bias)[c(2, 4:10)] ) }) -- GitLab From b0abde50a3886a3e7f23f234657037707193c370 Mon Sep 17 00:00:00 2001 From: vagudets Date: Fri, 21 Jun 2024 11:44:41 +0200 Subject: [PATCH 10/11] Return NA for significance when result is NA --- R/Bias.R | 11 +++++++++-- tests/testthat/test-Bias.R | 9 ++++++++- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/R/Bias.R b/R/Bias.R index 7b303de..d488718 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -189,8 +189,12 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, if (!is.null(alpha)) { if (!absolute) { - pval <- t.test(x = obs, y = exp, alternative = "two.sided")$p.value - sign <- pval <= alpha + if (all(is.na(bias))) { + sign <- NA + } else { + pval <- t.test(x = obs, y = exp, alternative = "two.sided")$p.value + sign <- pval <= alpha + } } } } else { @@ -218,6 +222,9 @@ Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, if (isTRUE(time_mean)) { bias <- MeanDims(bias, time_dim, na.rm = na.rm) + if (!is.null(sign)) { + sign[which(is.na(bias))] <- NA + } } } if (!is.null(alpha) && !absolute) { diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R index 4024897..bab78ed 100644 --- a/tests/testthat/test-Bias.R +++ b/tests/testthat/test-Bias.R @@ -255,15 +255,22 @@ test_that("5. Output checks: dat4", { as.vector(Bias(exp4, obs4, time_mean = F)$bias) ) # 4_2: NA - ## TODO: Significance expect_equal( as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset')$bias), as.numeric(NA) ) + expect_equal( + as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset')$sign), + as.numeric(NA) + ) expect_equal( as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F)$bias)[c(1, 3)], as.numeric(c(NA, NA)) ) + expect_equal( + as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F)$sign), + FALSE + ) expect_equal( as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F)$bias)[c(2, 4:10)], as.vector(Bias(exp4, obs4, time_mean = F)$bias)[c(2, 4:10)] -- GitLab From 341f3a589f5b2dfc90543aad817826293fce4ab7 Mon Sep 17 00:00:00 2001 From: vagudets Date: Fri, 21 Jun 2024 12:01:32 +0200 Subject: [PATCH 11/11] Fix pipeline --- tests/testthat/test-Bias.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R index bab78ed..3dee73a 100644 --- a/tests/testthat/test-Bias.R +++ b/tests/testthat/test-Bias.R @@ -261,7 +261,7 @@ test_that("5. Output checks: dat4", { ) expect_equal( as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset')$sign), - as.numeric(NA) + as.logical(NA) ) expect_equal( as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F)$bias)[c(1, 3)], -- GitLab