diff --git a/R/Bias.R b/R/Bias.R index 5ec8d0fd5d0a9c5c5891f31196c8135331a47f95..d4887189220daaad310cf7794870ffb108c6435d 100644 --- a/R/Bias.R +++ b/R/Bias.R @@ -27,6 +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 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. #' @@ -34,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. +#'(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 @@ -43,12 +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.01) +#'abs_bias <- Bias(exp = exp, obs = obs, memb_dim = 'member', 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, ncores = NULL) { +Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, + na.rm = FALSE, absolute = FALSE, time_mean = TRUE, + alpha = 0.05, ncores = NULL) { # Check inputs ## exp and obs (1) @@ -120,6 +128,17 @@ 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 (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)) { if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | @@ -144,16 +163,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 +186,33 @@ 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) { + if (all(is.na(bias))) { + sign <- NA + } else { + pval <- t.test(x = obs, y = exp, alternative = "two.sided")$p.value + sign <- 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)) + 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] + 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[i, j] <= alpha + } + } } } @@ -182,8 +222,14 @@ 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) { + return(list(bias = bias, sign = sign)) + } else { + return(bias) } - - return(bias) } diff --git a/man/Bias.Rd b/man/Bias.Rd index 2a02f2d52b9d39b8df0c7b8fa06ef02d702b1b11..e94beb54298a98aa07d3aadc878808a5ed2a9ca1 100644 --- a/man/Bias.Rd +++ b/man/Bias.Rd @@ -13,6 +13,7 @@ Bias( na.rm = FALSE, absolute = FALSE, time_mean = TRUE, + alpha = 0.05, 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,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.01) +abs_bias <- Bias(exp = exp, obs = obs, memb_dim = 'member', absolute = TRUE, alpha = NULL) } \references{ diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R index 4c6cc99f39cf5e981bad69ada0c44c599d97c5ee..3dee73a95b6d4dd8911b2e9e4ae9021bc3ed40e3 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)$sign), + 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')$sign), 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', 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,63 @@ 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 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')$sign), + as.logical(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))[c(2, 4:10)], - as.vector(Bias(exp4, obs4, time_mean = F))[c(2, 4:10)] + 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)] ) })