From 8291ee5c1b9ee40cbd6e8dea88ba0ce60c4cb8b7 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Tue, 30 Aug 2022 12:37:01 +0200 Subject: [PATCH 1/6] included two-sided test --- R/DiffCorr.R | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index c250814..8fdcbf0 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -38,6 +38,8 @@ #' steps with no missing values in all "exp", "ref", and "obs" will be used. If #' "na.fail", an error will arise if any of "exp", "ref", or "obs" contains any #' NA. The default value is "return.na". +#'@param two.sided A logical indicating whether to perform a two-sided test. If +#' FALSE, a one.sided test is performed. The default value is FALSE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -72,7 +74,7 @@ #'@export DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', memb_dim = NULL, method = 'pearson', alpha = NULL, - handle.na = 'return.na', ncores = NULL) { + handle.na = 'return.na', two.sided = FALSE, ncores = NULL) { # Check inputs ## exp, ref, and obs (1) @@ -143,6 +145,11 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', if (!handle.na %in% c('return.na', 'only.complete.triplets', 'na.fail')) { stop('Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".') } + ## two.sided + if (!is.logical(two.sided)){ + stop('Parameter "two.sided" must be TRUE or FALSE') + } + ## ncores if (!is.null(ncores)) { if (any(!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | length(ncores) > 1)) { @@ -184,22 +191,24 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', ref = time_dim, N.eff = NULL), output_dims = output_dims, fun = .DiffCorr, method = method, - alpha = alpha, handle.na = handle.na, ncores = ncores) + alpha = alpha, handle.na = handle.na, + two.sided = two.sided, ncores = ncores) } else { output <- Apply(data = list(exp = exp, obs = obs, ref = ref), target_dims = list(exp = time_dim, obs = time_dim, ref = time_dim), output_dims = output_dims, N.eff = N.eff, fun = .DiffCorr, method = method, - alpha = alpha, handle.na = handle.na, ncores = ncores) + alpha = alpha, handle.na = handle.na, + two.sided = two.sided, ncores = ncores) } return(output) } -.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na') { +.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na', two.sided = FALSE) { - .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL) { + .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, two.sided = FALSE) { # Correlation difference cor.exp <- cor(x = exp, y = obs, method = method) @@ -220,7 +229,11 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', if (is.null(alpha)) { output$p.val <- p.value } else { - output$sign <- ifelse(!is.na(p.value) & p.value <= alpha, TRUE, FALSE) + if (isFALSE(two.sided)){ + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha, TRUE, FALSE) + } else if (isTRUE(two.sided)){ + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha/2, TRUE, FALSE) + } } return(output) } -- GitLab From dfaad430a399cd9a858c688f1d1451079b62f2b6 Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 30 Aug 2022 17:31:35 +0200 Subject: [PATCH 2/6] Trivial document and format fix; unit test for two-sided added --- R/DiffCorr.R | 16 ++++++++-------- man/DiffCorr.Rd | 12 ++++++++---- tests/testthat/test-DiffCorr.R | 13 +++++++++++++ 3 files changed, 29 insertions(+), 12 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index 8fdcbf0..9fb6067 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -4,10 +4,10 @@ #'Positive values of the correlation difference indicate that the forecast is #'more skillful than the reference forecast, while negative values mean that #'the reference forecast is more skillful. The statistical significance of the -#'correlation differences is computed with a one-sided test for equality of -#'dependent correlation coefficients (Steiger, 1980; Siegert et al., 2017) using -#'effective degrees of freedom to account for the autocorrelation of the time -#'series (von Storch and Zwiers, 1999). +#'correlation differences is computed with a one-sided or two-sided test for +#'equality of dependent correlation coefficients (Steiger, 1980; Siegert et al., +#'2017) using effective degrees of freedom to account for the autocorrelation of +#'the time series (von Storch and Zwiers, 1999). #' #'@param exp A named numerical array of the forecast data with at least time #' dimension. @@ -146,8 +146,8 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', stop('Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".') } ## two.sided - if (!is.logical(two.sided)){ - stop('Parameter "two.sided" must be TRUE or FALSE') + if (!is.logical(two.sided)) { + stop("Parameter 'two.sided' must be TRUE or FALSE.") } ## ncores if (!is.null(ncores)) { @@ -229,9 +229,9 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', if (is.null(alpha)) { output$p.val <- p.value } else { - if (isFALSE(two.sided)){ + if (isFALSE(two.sided)) { output$sign <- ifelse(!is.na(p.value) & p.value <= alpha, TRUE, FALSE) - } else if (isTRUE(two.sided)){ + } else if (isTRUE(two.sided)) { output$sign <- ifelse(!is.na(p.value) & p.value <= alpha/2, TRUE, FALSE) } } diff --git a/man/DiffCorr.Rd b/man/DiffCorr.Rd index d8ff65c..ca16cd4 100644 --- a/man/DiffCorr.Rd +++ b/man/DiffCorr.Rd @@ -14,6 +14,7 @@ DiffCorr( method = "pearson", alpha = NULL, handle.na = "return.na", + two.sided = FALSE, ncores = NULL ) } @@ -56,6 +57,9 @@ steps with no missing values in all "exp", "ref", and "obs" will be used. If "na.fail", an error will arise if any of "exp", "ref", or "obs" contains any NA. The default value is "return.na".} +\item{two.sided}{A logical indicating whether to perform a two-sided test. If +FALSE, a one.sided test is performed. The default value is FALSE.} + \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } @@ -81,10 +85,10 @@ Compute the correlation difference between two deterministic forecasts. Positive values of the correlation difference indicate that the forecast is more skillful than the reference forecast, while negative values mean that the reference forecast is more skillful. The statistical significance of the -correlation differences is computed with a one-sided test for equality of -dependent correlation coefficients (Steiger, 1980; Siegert et al., 2017) using -effective degrees of freedom to account for the autocorrelation of the time -series (von Storch and Zwiers, 1999). +correlation differences is computed with a one-sided or two-sided test for +equality of dependent correlation coefficients (Steiger, 1980; Siegert et al., +2017) using effective degrees of freedom to account for the autocorrelation of +the time series (von Storch and Zwiers, 1999). } \examples{ exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50)) diff --git a/tests/testthat/test-DiffCorr.R b/tests/testthat/test-DiffCorr.R index e0834e1..9a087fc 100644 --- a/tests/testthat/test-DiffCorr.R +++ b/tests/testthat/test-DiffCorr.R @@ -89,6 +89,11 @@ test_that("1. Input checks", { DiffCorr(exp2, obs2, ref2, handle.na = TRUE), 'Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".' ) + # two.sided + expect_error( + DiffCorr(exp2, obs2, ref2, two.sided = 2), + "Parameter 'two.sided' must be TRUE or FALSE." + ) # ncores expect_error( DiffCorr(exp2, obs2, ref2, ncores = 1.5), @@ -131,10 +136,18 @@ c(0.27347087, 0.50556882, 0.08855968, 0.24199701, 0.22935182, 0.88336336), tolerance = 0.0001 ) expect_equal( +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$diff.corr), +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, two.sided = T)$diff.corr) +) +expect_equal( as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$sign), rep(FALSE, 6) ) expect_equal( +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, two.sided = T)$sign), +rep(FALSE, 6) +) +expect_equal( suppressWarnings(as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "spearman")$diff.corr)), c(0.07272727, 0.54545455, 0.10909091, -0.01212121, -0.03636364, 1.01818182), tolerance = 0.0001 -- GitLab From e59119cff930f6569594f129fc5481a7b087c702 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Wed, 31 Aug 2022 12:23:28 +0200 Subject: [PATCH 3/6] posible tests: two-sided, one-sided-higher and one-sided-lower --- R/DiffCorr.R | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index 9fb6067..0775106 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -38,8 +38,11 @@ #' steps with no missing values in all "exp", "ref", and "obs" will be used. If #' "na.fail", an error will arise if any of "exp", "ref", or "obs" contains any #' NA. The default value is "return.na". -#'@param two.sided A logical indicating whether to perform a two-sided test. If -#' FALSE, a one.sided test is performed. The default value is FALSE. +#'@param type A character string indicating the type of test. It can be "two-sided" +#' (to assess whether the skill of "exp" and "ref" are significantly different), +#' "one-sided-higher" (to assess whether the skill of "exp" is significantly higher +#' than that of "ref"), or "one-sided-lower" (to assess whether the skill of "exp" +#' is significantly lower than that of "ref"). The default value is "two-sided". #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -74,7 +77,7 @@ #'@export DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', memb_dim = NULL, method = 'pearson', alpha = NULL, - handle.na = 'return.na', two.sided = FALSE, ncores = NULL) { + handle.na = 'return.na', type = "two-sided", ncores = NULL) { # Check inputs ## exp, ref, and obs (1) @@ -145,9 +148,9 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', if (!handle.na %in% c('return.na', 'only.complete.triplets', 'na.fail')) { stop('Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".') } - ## two.sided - if (!is.logical(two.sided)) { - stop("Parameter 'two.sided' must be TRUE or FALSE.") + ## type + if (!type %in% c('two-sided','one-sided-higher','one-sided-lower')) { + stop("Parameter 'type' must be 'two-sided', 'one-sided-higher' or 'one-sided-lower'.") } ## ncores if (!is.null(ncores)) { @@ -192,7 +195,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', output_dims = output_dims, fun = .DiffCorr, method = method, alpha = alpha, handle.na = handle.na, - two.sided = two.sided, ncores = ncores) + type = type, ncores = ncores) } else { output <- Apply(data = list(exp = exp, obs = obs, ref = ref), target_dims = list(exp = time_dim, obs = time_dim, @@ -200,15 +203,15 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', output_dims = output_dims, N.eff = N.eff, fun = .DiffCorr, method = method, alpha = alpha, handle.na = handle.na, - two.sided = two.sided, ncores = ncores) + type = type, ncores = ncores) } return(output) } -.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na', two.sided = FALSE) { +.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na', type = 'two.sided') { - .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, two.sided = FALSE) { + .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, type = 'two.sided') { # Correlation difference cor.exp <- cor(x = exp, y = obs, method = method) @@ -229,10 +232,12 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', if (is.null(alpha)) { output$p.val <- p.value } else { - if (isFALSE(two.sided)) { - output$sign <- ifelse(!is.na(p.value) & p.value <= alpha, TRUE, FALSE) - } else if (isTRUE(two.sided)) { + if (type == 'two-sided') { output$sign <- ifelse(!is.na(p.value) & p.value <= alpha/2, TRUE, FALSE) + } else if (type == 'one-sided-higher') { + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr > 0, TRUE, FALSE) + } else if (type == 'one-sided-lower') { + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr < 0, TRUE, FALSE) } } return(output) @@ -250,7 +255,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', ref <- ref[!nna] output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha) + N.eff = N.eff, alpha = alpha, type = type) } else if (handle.na == 'return.na') { # Data contain NA, return NAs directly without passing to .diff.corr @@ -263,7 +268,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', } else { ## There is no NA output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha) + N.eff = N.eff, alpha = alpha, type = type) } return(output) } -- GitLab From 831ef2d316c7f5245a273f453bcef23a28a6abf1 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 31 Aug 2022 17:17:03 +0200 Subject: [PATCH 4/6] Change param 'type' to 'test.type'; amend unit tests --- R/DiffCorr.R | 41 +++++++++++++++++++--------------- man/DiffCorr.Rd | 10 ++++++--- tests/testthat/test-DiffCorr.R | 26 ++++++++++++++------- 3 files changed, 48 insertions(+), 29 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index 0775106..f4ff2a9 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -38,11 +38,12 @@ #' steps with no missing values in all "exp", "ref", and "obs" will be used. If #' "na.fail", an error will arise if any of "exp", "ref", or "obs" contains any #' NA. The default value is "return.na". -#'@param type A character string indicating the type of test. It can be "two-sided" -#' (to assess whether the skill of "exp" and "ref" are significantly different), -#' "one-sided-higher" (to assess whether the skill of "exp" is significantly higher -#' than that of "ref"), or "one-sided-lower" (to assess whether the skill of "exp" -#' is significantly lower than that of "ref"). The default value is "two-sided". +#'@param test.type A character string indicating the type of significance test. +#' It can be "two-sided" (to assess whether the skill of "exp" and "ref" are +#' significantly different), "one-sided-higher" (to assess whether the skill of +#' "exp" is significantly higher than that of "ref"), or "one-sided-lower" (to +#' assess whether the skill of "exp" is significantly lower than that of +#' "ref"). The default value is "two-sided". #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -77,7 +78,7 @@ #'@export DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', memb_dim = NULL, method = 'pearson', alpha = NULL, - handle.na = 'return.na', type = "two-sided", ncores = NULL) { + handle.na = 'return.na', test.type = "two-sided", ncores = NULL) { # Check inputs ## exp, ref, and obs (1) @@ -148,10 +149,13 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', if (!handle.na %in% c('return.na', 'only.complete.triplets', 'na.fail')) { stop('Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".') } - ## type - if (!type %in% c('two-sided','one-sided-higher','one-sided-lower')) { - stop("Parameter 'type' must be 'two-sided', 'one-sided-higher' or 'one-sided-lower'.") + ## test.type + if (!test.type %in% c('two-sided', 'one-sided-higher', 'one-sided-lower')) { + stop("Parameter 'test.type' must be 'two-sided', 'one-sided-higher' or 'one-sided-lower'.") } + #NOTE: warning can be removed in the next release + warning("The significance test has changed after s2dv_1.2.0. The default method is 'two-sided' and the one-sided test is distinguished into the higher or lower side.") + ## ncores if (!is.null(ncores)) { if (any(!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | @@ -195,7 +199,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', output_dims = output_dims, fun = .DiffCorr, method = method, alpha = alpha, handle.na = handle.na, - type = type, ncores = ncores) + test.type = test.type, ncores = ncores) } else { output <- Apply(data = list(exp = exp, obs = obs, ref = ref), target_dims = list(exp = time_dim, obs = time_dim, @@ -203,15 +207,16 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', output_dims = output_dims, N.eff = N.eff, fun = .DiffCorr, method = method, alpha = alpha, handle.na = handle.na, - type = type, ncores = ncores) + test.type = test.type, ncores = ncores) } return(output) } -.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na', type = 'two.sided') { +.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, + handle.na = 'return.na', test.type = 'two.sided') { - .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, type = 'two.sided') { + .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, test.type = 'two.sided') { # Correlation difference cor.exp <- cor(x = exp, y = obs, method = method) @@ -232,11 +237,11 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', if (is.null(alpha)) { output$p.val <- p.value } else { - if (type == 'two-sided') { + if (test.type == 'two-sided') { output$sign <- ifelse(!is.na(p.value) & p.value <= alpha/2, TRUE, FALSE) - } else if (type == 'one-sided-higher') { + } else if (test.type == 'one-sided-higher') { output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr > 0, TRUE, FALSE) - } else if (type == 'one-sided-lower') { + } else if (test.type == 'one-sided-lower') { output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr < 0, TRUE, FALSE) } } @@ -255,7 +260,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', ref <- ref[!nna] output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha, type = type) + N.eff = N.eff, alpha = alpha, test.type = test.type) } else if (handle.na == 'return.na') { # Data contain NA, return NAs directly without passing to .diff.corr @@ -268,7 +273,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', } else { ## There is no NA output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha, type = type) + N.eff = N.eff, alpha = alpha, test.type = test.type) } return(output) } diff --git a/man/DiffCorr.Rd b/man/DiffCorr.Rd index ca16cd4..36b2fb9 100644 --- a/man/DiffCorr.Rd +++ b/man/DiffCorr.Rd @@ -14,7 +14,7 @@ DiffCorr( method = "pearson", alpha = NULL, handle.na = "return.na", - two.sided = FALSE, + test.type = "two-sided", ncores = NULL ) } @@ -57,8 +57,12 @@ steps with no missing values in all "exp", "ref", and "obs" will be used. If "na.fail", an error will arise if any of "exp", "ref", or "obs" contains any NA. The default value is "return.na".} -\item{two.sided}{A logical indicating whether to perform a two-sided test. If -FALSE, a one.sided test is performed. The default value is FALSE.} +\item{test.type}{A character string indicating the type of significance test. +It can be "two-sided" (to assess whether the skill of "exp" and "ref" are +significantly different), "one-sided-higher" (to assess whether the skill of +"exp" is significantly higher than that of "ref"), or "one-sided-lower" (to +assess whether the skill of "exp" is significantly lower than that of +"ref"). The default value is "two-sided".} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} diff --git a/tests/testthat/test-DiffCorr.R b/tests/testthat/test-DiffCorr.R index 9a087fc..60ac7d4 100644 --- a/tests/testthat/test-DiffCorr.R +++ b/tests/testthat/test-DiffCorr.R @@ -89,14 +89,14 @@ test_that("1. Input checks", { DiffCorr(exp2, obs2, ref2, handle.na = TRUE), 'Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".' ) - # two.sided + # test.type expect_error( - DiffCorr(exp2, obs2, ref2, two.sided = 2), - "Parameter 'two.sided' must be TRUE or FALSE." + DiffCorr(exp2, obs2, ref2, test.type = "two.sided"), + "Parameter 'test.type' must be 'two-sided', 'one-sided-higher' or 'one-sided-lower'." ) # ncores expect_error( - DiffCorr(exp2, obs2, ref2, ncores = 1.5), + suppressWarnings(DiffCorr(exp2, obs2, ref2, ncores = 1.5)), 'Parameter "ncores" must be either NULL or a positive integer.' ) @@ -104,6 +104,8 @@ test_that("1. Input checks", { ############################################## test_that("2. Output checks: dat1", { +suppressWarnings({ + expect_equal( names(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb')), c("diff.corr", "p.val") @@ -136,15 +138,19 @@ c(0.27347087, 0.50556882, 0.08855968, 0.24199701, 0.22935182, 0.88336336), tolerance = 0.0001 ) expect_equal( -as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$diff.corr), -as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, two.sided = T)$diff.corr) +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided-higher")$diff.corr), +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "two-sided")$diff.corr) ) expect_equal( as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$sign), rep(FALSE, 6) ) expect_equal( -as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, two.sided = T)$sign), +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided-higher")$sign), +rep(FALSE, 6) +) +expect_equal( +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided-lower")$sign), rep(FALSE, 6) ) expect_equal( @@ -194,12 +200,15 @@ DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', handle.na = 'only.complete.triplet "There is no complete set of forecasts and observations." ) - +}) # suppressWarnings }) ############################################## test_that("3. Output checks: dat2", { + +suppressWarnings({ + expect_equal( names(DiffCorr(exp2, obs2, ref2)), c("diff.corr", "p.val") @@ -223,4 +232,5 @@ DiffCorr(exp2, obs2, ref2)$diff, tolerance = 0.0001 ) +}) # suppressWarnings }) -- GitLab From 144d56ae10c62b70ad13a89db2f831037b1cac20 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Mon, 5 Sep 2022 13:32:58 +0200 Subject: [PATCH 5/6] Z-test for two-sided --- R/DiffCorr.R | 72 +++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 51 insertions(+), 21 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index f4ff2a9..277e4d2 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -40,10 +40,10 @@ #' NA. The default value is "return.na". #'@param test.type A character string indicating the type of significance test. #' It can be "two-sided" (to assess whether the skill of "exp" and "ref" are -#' significantly different), "one-sided-higher" (to assess whether the skill of -#' "exp" is significantly higher than that of "ref"), or "one-sided-lower" (to -#' assess whether the skill of "exp" is significantly lower than that of -#' "ref"). The default value is "two-sided". +#' significantly different with a z-test on Fisher z-transformed correlation +#' coefficients) or "one-sided" (to assess whether the skill of "exp" is +#' significantly higher than that of "ref" with the test described in +#' Steiger, 1980). The default value is "two-sided". #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -150,11 +150,11 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', stop('Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".') } ## test.type - if (!test.type %in% c('two-sided', 'one-sided-higher', 'one-sided-lower')) { - stop("Parameter 'test.type' must be 'two-sided', 'one-sided-higher' or 'one-sided-lower'.") + if (!test.type %in% c('two-sided', 'one-sided')) { + stop("Parameter 'test.type' must be 'two-sided' or 'one-sided'.") } #NOTE: warning can be removed in the next release - warning("The significance test has changed after s2dv_1.2.0. The default method is 'two-sided' and the one-sided test is distinguished into the higher or lower side.") + warning("The default significance test has changed after s2dv_1.2.0. The default method is 'two-sided'.") ## ncores if (!is.null(ncores)) { @@ -216,7 +216,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', .DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na', test.type = 'two.sided') { - .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, test.type = 'two.sided') { + .diff.corr_one.sided <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL) { # Correlation difference cor.exp <- cor(x = exp, y = obs, method = method) @@ -233,21 +233,40 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', } R <- (1 - r12 * r12 - r13 * r13 - r23 * r23) + 2 * r12 * r13 * r23 t <- (r12 - r13) * sqrt((N.eff - 1) * (1 + r23) / (2 * ((N.eff - 1) / (N.eff - 3)) * R + 0.25 * (r12 + r13)^2 * (1 - r23)^3)) - p.value <- 1 - pt(t, df = N.eff - 3) + p.value <- pt(t, df = N.eff - 3, lower.tail = FALSE) if (is.null(alpha)) { output$p.val <- p.value } else { - if (test.type == 'two-sided') { - output$sign <- ifelse(!is.na(p.value) & p.value <= alpha/2, TRUE, FALSE) - } else if (test.type == 'one-sided-higher') { - output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr > 0, TRUE, FALSE) - } else if (test.type == 'one-sided-lower') { - output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr < 0, TRUE, FALSE) - } + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr > 0, TRUE, FALSE) } return(output) } - + + .diff.corr_two.sided <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL) { + + # Correlation difference + cor.exp <- cor(x = exp, y = obs, method = method) + cor.ref <- cor(x = ref, y = obs, method = method) + output <- NULL + output$diff.corr <- cor.exp - cor.ref + + # Significance with two-sided z-test on Fisher z-transformed correlation coefficients (Hinkel et al, 1988) + if (is.na(N.eff)) { + N.eff <- .Eno(x = obs, na.action = na.pass) ## effective degrees of freedom + } + Z <- abs(0.5*log((1+cor.exp)/(1-cor.exp)) - 0.5*log((1+cor.ref)/(1-cor.ref))) / sqrt(1/(N.eff-3) + 1/(N.eff-3)) + # sign <- Z >= qnorm(p = alpha/2, lower.tail = FALSE) + p.value <- pnorm(q = Z, lower.tail = FALSE) + + if (is.null(alpha)) { + output$p.val <- p.value + } else { + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha/2, TRUE, FALSE) + } + + return(output) + } + #================================================== if (anyNA(exp) | anyNA(obs) | anyNA(ref)) { ## There are NAs @@ -259,8 +278,13 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', obs <- obs[!nna] ref <- ref[!nna] - output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha, test.type = test.type) + if (identical(test.type,'two-sided')){ + output <- .diff.corr_two.sided(exp = exp, obs = obs, ref = ref, method = method, + N.eff = N.eff, alpha = alpha) + } else if (identical(test.type,'one-sided')){ + output <- .diff.corr_one.sided(exp = exp, obs = obs, ref = ref, method = method, + N.eff = N.eff, alpha = alpha) + } } else if (handle.na == 'return.na') { # Data contain NA, return NAs directly without passing to .diff.corr @@ -272,8 +296,14 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', } } else { ## There is no NA - output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha, test.type = test.type) + + if (identical(test.type,'two-sided')){ + output <- .diff.corr_two.sided(exp = exp, obs = obs, ref = ref, method = method, + N.eff = N.eff, alpha = alpha) + } else if (identical(test.type,'one-sided')){ + output <- .diff.corr_one.sided(exp = exp, obs = obs, ref = ref, method = method, + N.eff = N.eff, alpha = alpha) + } } return(output) } -- GitLab From 83a0ea11a746388ab68bb813738b60913bb8d603 Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 22 Sep 2022 12:55:37 +0200 Subject: [PATCH 6/6] Combine two .diff.corr(); correct unit tests --- R/DiffCorr.R | 97 ++++++++++++++-------------------- man/DiffCorr.Rd | 11 ++-- tests/testthat/test-DiffCorr.R | 30 +++++++---- 3 files changed, 65 insertions(+), 73 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index 277e4d2..95996b6 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -41,9 +41,9 @@ #'@param test.type A character string indicating the type of significance test. #' It can be "two-sided" (to assess whether the skill of "exp" and "ref" are #' significantly different with a z-test on Fisher z-transformed correlation -#' coefficients) or "one-sided" (to assess whether the skill of "exp" is -#' significantly higher than that of "ref" with the test described in -#' Steiger, 1980). The default value is "two-sided". +#' coefficients, Hinkel et al, 1988) or "one-sided" (to assess whether the +#' skill of "exp" is significantly higher than that of "ref" with the test +#' described in Steiger, 1980). The default value is "two-sided". #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -72,7 +72,8 @@ #' exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50)) #' obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50)) #' ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50)) -#' res <- DiffCorr(exp, obs, ref, memb_dim = 'member') +#' res_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', alpha = 0.05) +#' res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided') #' #'@import multiApply #'@export @@ -216,54 +217,46 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', .DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na', test.type = 'two.sided') { - .diff.corr_one.sided <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL) { - + .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, test.type = 'two.sided') { # Correlation difference cor.exp <- cor(x = exp, y = obs, method = method) cor.ref <- cor(x = ref, y = obs, method = method) output <- NULL output$diff.corr <- cor.exp - cor.ref - - # Significance with one-sided test for equality of dependent correlation coefficients (Steiger, 1980) - r12 <- cor.exp - r13 <- cor.ref - r23 <- cor(exp, ref) - if (is.na(N.eff)) { - N.eff <- .Eno(x = obs, na.action = na.pass) ## effective degrees of freedom - } - R <- (1 - r12 * r12 - r13 * r13 - r23 * r23) + 2 * r12 * r13 * r23 - t <- (r12 - r13) * sqrt((N.eff - 1) * (1 + r23) / (2 * ((N.eff - 1) / (N.eff - 3)) * R + 0.25 * (r12 + r13)^2 * (1 - r23)^3)) - p.value <- pt(t, df = N.eff - 3, lower.tail = FALSE) - if (is.null(alpha)) { - output$p.val <- p.value - } else { - output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr > 0, TRUE, FALSE) - } - return(output) - } - - .diff.corr_two.sided <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL) { - - # Correlation difference - cor.exp <- cor(x = exp, y = obs, method = method) - cor.ref <- cor(x = ref, y = obs, method = method) - output <- NULL - output$diff.corr <- cor.exp - cor.ref - - # Significance with two-sided z-test on Fisher z-transformed correlation coefficients (Hinkel et al, 1988) if (is.na(N.eff)) { N.eff <- .Eno(x = obs, na.action = na.pass) ## effective degrees of freedom } - Z <- abs(0.5*log((1+cor.exp)/(1-cor.exp)) - 0.5*log((1+cor.ref)/(1-cor.ref))) / sqrt(1/(N.eff-3) + 1/(N.eff-3)) - # sign <- Z >= qnorm(p = alpha/2, lower.tail = FALSE) - p.value <- pnorm(q = Z, lower.tail = FALSE) - - if (is.null(alpha)) { - output$p.val <- p.value + + if (test.type == 'one-sided') { + # Significance with one-sided test for equality of dependent correlation coefficients (Steiger, 1980) + r12 <- cor.exp + r13 <- cor.ref + r23 <- cor(exp, ref) + R <- (1 - r12 * r12 - r13 * r13 - r23 * r23) + 2 * r12 * r13 * r23 + t <- (r12 - r13) * sqrt((N.eff - 1) * (1 + r23) / (2 * ((N.eff - 1) / (N.eff - 3)) * R + 0.25 * (r12 + r13)^2 * (1 - r23)^3)) + p.value <- pt(t, df = N.eff - 3, lower.tail = FALSE) + if (is.null(alpha)) { + output$p.val <- p.value + } else { + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr > 0, TRUE, FALSE) + } + + } else if (test.type == 'two-sided') { + # Significance with two-sided z-test on Fisher z-transformed correlation coefficients (Hinkel et al, 1988) + Z <- abs(0.5 * log((1 + cor.exp) / (1 - cor.exp)) - 0.5 * log((1 + cor.ref) / (1 - cor.ref))) / sqrt(1/ (N.eff - 3) + 1/(N.eff - 3)) + # sign <- Z >= qnorm(p = alpha/2, lower.tail = FALSE) + p.value <- pnorm(q = Z, lower.tail = FALSE) + + if (is.null(alpha)) { + output$p.val <- p.value + } else { + output$sign <- ifelse(!is.na(p.value) & p.value <= alpha / 2, TRUE, FALSE) + } + } else { - output$sign <- ifelse(!is.na(p.value) & p.value <= alpha/2, TRUE, FALSE) + stop("Parameter 'test.type' is not supported.") } - + return(output) } @@ -278,13 +271,8 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', obs <- obs[!nna] ref <- ref[!nna] - if (identical(test.type,'two-sided')){ - output <- .diff.corr_two.sided(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha) - } else if (identical(test.type,'one-sided')){ - output <- .diff.corr_one.sided(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha) - } + output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, + N.eff = N.eff, alpha = alpha, test.type = test.type) } else if (handle.na == 'return.na') { # Data contain NA, return NAs directly without passing to .diff.corr @@ -296,14 +284,9 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', } } else { ## There is no NA - - if (identical(test.type,'two-sided')){ - output <- .diff.corr_two.sided(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha) - } else if (identical(test.type,'one-sided')){ - output <- .diff.corr_one.sided(exp = exp, obs = obs, ref = ref, method = method, - N.eff = N.eff, alpha = alpha) - } + output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, + N.eff = N.eff, alpha = alpha, test.type = test.type) } + return(output) } diff --git a/man/DiffCorr.Rd b/man/DiffCorr.Rd index 36b2fb9..e0e5ec3 100644 --- a/man/DiffCorr.Rd +++ b/man/DiffCorr.Rd @@ -59,10 +59,10 @@ NA. The default value is "return.na".} \item{test.type}{A character string indicating the type of significance test. It can be "two-sided" (to assess whether the skill of "exp" and "ref" are -significantly different), "one-sided-higher" (to assess whether the skill of -"exp" is significantly higher than that of "ref"), or "one-sided-lower" (to -assess whether the skill of "exp" is significantly lower than that of -"ref"). The default value is "two-sided".} +significantly different with a z-test on Fisher z-transformed correlation +coefficients, Hinkel et al, 1988) or "one-sided" (to assess whether the +skill of "exp" is significantly higher than that of "ref" with the test +described in Steiger, 1980). The default value is "two-sided".} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} @@ -98,7 +98,8 @@ the time series (von Storch and Zwiers, 1999). exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50)) obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50)) ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50)) -res <- DiffCorr(exp, obs, ref, memb_dim = 'member') +res_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', alpha = 0.05) +res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided') } \references{ diff --git a/tests/testthat/test-DiffCorr.R b/tests/testthat/test-DiffCorr.R index 60ac7d4..b8c8029 100644 --- a/tests/testthat/test-DiffCorr.R +++ b/tests/testthat/test-DiffCorr.R @@ -20,7 +20,6 @@ set.seed(3) obs2 <- array(rnorm(10), dim = c(sdate = 10)) - ############################################## test_that("1. Input checks", { # exp and obs (1) @@ -92,7 +91,7 @@ test_that("1. Input checks", { # test.type expect_error( DiffCorr(exp2, obs2, ref2, test.type = "two.sided"), - "Parameter 'test.type' must be 'two-sided', 'one-sided-higher' or 'one-sided-lower'." + "Parameter 'test.type' must be 'two-sided' or 'one-sided'." ) # ncores expect_error( @@ -125,7 +124,7 @@ tolerance = 0.0001 ) expect_equal( as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb')$p), -c(0.26166060, 0.15899774, 0.39264452, 0.27959883, 0.34736305, 0.07479832), +c(0.28984782, 0.14395946, 0.41374677, 0.32253784, 0.32677872, 0.03659053), tolerance = 0.0001 ) expect_equal( @@ -138,7 +137,7 @@ c(0.27347087, 0.50556882, 0.08855968, 0.24199701, 0.22935182, 0.88336336), tolerance = 0.0001 ) expect_equal( -as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided-higher")$diff.corr), +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided")$diff.corr), as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "two-sided")$diff.corr) ) expect_equal( @@ -146,11 +145,7 @@ as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$sign), rep(FALSE, 6) ) expect_equal( -as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided-higher")$sign), -rep(FALSE, 6) -) -expect_equal( -as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided-lower")$sign), +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided")$sign), rep(FALSE, 6) ) expect_equal( @@ -160,7 +155,7 @@ tolerance = 0.0001 ) expect_equal( suppressWarnings(as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "spearman")$p)), -c(0.4358970, 0.1341575, 0.3448977, 0.5114262, 0.5264872, 0.0437861), +c(0.44479470, 0.11854232, 0.36667620, 0.49095264, 0.46962814, 0.01714977), tolerance = 0.0001 ) expect_equal( @@ -170,7 +165,7 @@ tolerance = 0.0001 ) expect_equal( as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', N.eff = Neff1)$p), -c(0.27841537, 0.15899774, 0.40096749, 0.27959883, 0.35889690, 0.07479832), +c(0.30406438, 0.14395946, 0.42005522, 0.32253784, 0.33887632, 0.03659053), tolerance = 0.0001 ) @@ -223,10 +218,23 @@ NULL ) expect_equal( DiffCorr(exp2, obs2, ref2)$p, +0.3201317, +tolerance = 0.0001 +) +expect_equal( +DiffCorr(exp2, obs2, ref2, test.type = 'one-sided')$p, 0.6577392, tolerance = 0.0001 ) expect_equal( +DiffCorr(exp2, obs2, ref2, test.type = 'one-sided', alpha = 0.5)$sign, +FALSE +) +expect_equal( +DiffCorr(exp2, obs2, ref2, test.type = 'two-sided', alpha = 0.65)$sign, +TRUE +) +expect_equal( DiffCorr(exp2, obs2, ref2)$diff, -0.2434725, tolerance = 0.0001 -- GitLab