From b1f57140beed23fe8ee9b6db977f2e531f720ab0 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
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 <nperez@bscearth372.int.bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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 <victoria.agudetse@bsc.es>
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