From 2369358026689133be52eed48e679d9ba4ac4f0d Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Tue, 11 Oct 2022 16:24:22 +0200 Subject: [PATCH 1/3] two-sided test with Steiger (1980) --- R/DiffCorr.R | 38 ++++++++++++++++++++++---------------- 1 file changed, 22 insertions(+), 16 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index 95996b6..0fec699 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -7,7 +7,7 @@ #'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). +#'the time series (Zwiers and von Storch, 1995). #' #'@param exp A named numerical array of the forecast data with at least time #' dimension. @@ -40,10 +40,9 @@ #' 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 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". +#' significantly different) or "one-sided" (to assess whether the +#' skill of "exp" is significantly higher than that of "ref") following 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. #' @@ -66,7 +65,7 @@ #'@references #'Steiger, 1980; https://content.apa.org/doi/10.1037/0033-2909.87.2.245 #'Siegert et al., 2017; https://doi.org/10.1175/MWR-D-16-0037.1 -#'von Storch and Zwiers, 1999; https://doi.org/10.1017/CBO9780511612336 +#'Zwiers and von Storch, 1995; https://doi.org/10.1175/1520-0442(1995)008%3C0336:TSCIAI%3E2.0.CO;2 #' #'@examples #' exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50)) @@ -227,14 +226,20 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', N.eff <- .Eno(x = obs, na.action = na.pass) ## effective degrees of freedom } + # Significance with one-sided or two-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)) + 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)) + + ## H0: the skill of exp is not higher than that of ref + ## H1: the skill of exp is higher than that of ref + p.value <- pt(t, df = N.eff - 3, lower.tail = FALSE) + if (is.null(alpha)) { output$p.val <- p.value } else { @@ -242,10 +247,11 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', } } 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) + + ## H0: the skill difference of exp and ref is zero + ## H1: the skill difference of exp and ref is different from zero + + p.value <- pt(abs(t), df = N.eff - 3, lower.tail = FALSE) if (is.null(alpha)) { output$p.val <- p.value -- GitLab From 2830084f975e6b900380d789b017589e35ce195f Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Tue, 11 Oct 2022 16:30:10 +0200 Subject: [PATCH 2/3] modified tests --- R/DiffCorr.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index 0fec699..270853f 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -71,8 +71,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_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') +#' res_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'two-sided', alpha = 0.05) +#' res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided', alpha = 0.05) #' #'@import multiApply #'@export -- GitLab From f124d48d3d862abd7b6f687ec64a8dba4a23fb8d Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 13 Oct 2022 13:19:22 +0200 Subject: [PATCH 3/3] Update DiffCorr.Rd and correct unit test --- R/DiffCorr.R | 10 +++++----- man/DiffCorr.Rd | 15 +++++++-------- tests/testthat/test-DiffCorr.R | 34 +++++++++++++++++++++++++--------- 3 files changed, 37 insertions(+), 22 deletions(-) diff --git a/R/DiffCorr.R b/R/DiffCorr.R index 270853f..a1fbae3 100644 --- a/R/DiffCorr.R +++ b/R/DiffCorr.R @@ -40,9 +40,9 @@ #' 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) or "one-sided" (to assess whether the -#' skill of "exp" is significantly higher than that of "ref") following Steiger -#' (1980). The default value is "two-sided". +#' significantly different) or "one-sided" (to assess whether the skill of +#' "exp" is significantly higher than that of "ref") following 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. #' @@ -65,14 +65,14 @@ #'@references #'Steiger, 1980; https://content.apa.org/doi/10.1037/0033-2909.87.2.245 #'Siegert et al., 2017; https://doi.org/10.1175/MWR-D-16-0037.1 -#'Zwiers and von Storch, 1995; https://doi.org/10.1175/1520-0442(1995)008%3C0336:TSCIAI%3E2.0.CO;2 +#'Zwiers and von Storch, 1995; https://doi.org/10.1175/1520-0442(1995)008<0336:TSCIAI>2.0.CO;2 #' #'@examples #' 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_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'two-sided', alpha = 0.05) -#' res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided', alpha = 0.05) +#' res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided', alpha = NULL) #' #'@import multiApply #'@export diff --git a/man/DiffCorr.Rd b/man/DiffCorr.Rd index e0e5ec3..e907d3b 100644 --- a/man/DiffCorr.Rd +++ b/man/DiffCorr.Rd @@ -59,10 +59,9 @@ 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 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".} +significantly different) or "one-sided" (to assess whether the skill of +"exp" is significantly higher than that of "ref") following 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.} @@ -92,18 +91,18 @@ the reference forecast is more skillful. The statistical significance of the 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). +the time series (Zwiers and von Storch, 1995). } \examples{ 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_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') +res_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'two-sided', alpha = 0.05) +res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided', alpha = NULL) } \references{ Steiger, 1980; https://content.apa.org/doi/10.1037/0033-2909.87.2.245 Siegert et al., 2017; https://doi.org/10.1175/MWR-D-16-0037.1 -von Storch and Zwiers, 1999; https://doi.org/10.1017/CBO9780511612336 +Zwiers and von Storch, 1995; https://doi.org/10.1175/1520-0442(1995)008<0336:TSCIAI>2.0.CO;2 } diff --git a/tests/testthat/test-DiffCorr.R b/tests/testthat/test-DiffCorr.R index b8c8029..3c35818 100644 --- a/tests/testthat/test-DiffCorr.R +++ b/tests/testthat/test-DiffCorr.R @@ -18,6 +18,9 @@ set.seed(2) ref2 <- array(rnorm(10), dim = c(sdate = 10)) set.seed(3) obs2 <- array(rnorm(10), dim = c(sdate = 10)) +## generate time auto-correlation (Eno will change) +set.seed(3) +obs2_2 <- array(1:10 + rnorm(10), dim = c(sdate = 10)) ############################################## @@ -123,8 +126,8 @@ 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')$p), -c(0.28984782, 0.14395946, 0.41374677, 0.32253784, 0.32677872, 0.03659053), +as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', test.type = "two-sided")$p), +c(0.26166060, 0.15899774, 0.39264452, 0.27959883, 0.34736305, 0.07479832), tolerance = 0.0001 ) expect_equal( @@ -155,7 +158,7 @@ tolerance = 0.0001 ) expect_equal( suppressWarnings(as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "spearman")$p)), -c(0.44479470, 0.11854232, 0.36667620, 0.49095264, 0.46962814, 0.01714977), +c(0.4358970, 0.1341575, 0.3448977, 0.4885738, 0.4735128, 0.0437861), tolerance = 0.0001 ) expect_equal( @@ -165,11 +168,10 @@ tolerance = 0.0001 ) expect_equal( as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', N.eff = Neff1)$p), -c(0.30406438, 0.14395946, 0.42005522, 0.32253784, 0.33887632, 0.03659053), +c(0.27841537, 0.15899774, 0.40096749, 0.27959883, 0.35889690, 0.07479832), tolerance = 0.0001 ) - #--------------------------- exp1[1] <- NA expect_equal( @@ -217,8 +219,8 @@ dim(DiffCorr(exp2, obs2, ref2)$p), NULL ) expect_equal( -DiffCorr(exp2, obs2, ref2)$p, -0.3201317, +DiffCorr(exp2, obs2, ref2, test.type = 'two-sided')$p, +0.3422608, tolerance = 0.0001 ) expect_equal( @@ -227,11 +229,11 @@ DiffCorr(exp2, obs2, ref2, test.type = 'one-sided')$p, tolerance = 0.0001 ) expect_equal( -DiffCorr(exp2, obs2, ref2, test.type = 'one-sided', alpha = 0.5)$sign, +DiffCorr(exp2, obs2, ref2, test.type = 'one-sided', alpha = 0.7)$sign, FALSE ) expect_equal( -DiffCorr(exp2, obs2, ref2, test.type = 'two-sided', alpha = 0.65)$sign, +DiffCorr(exp2, obs2, ref2, test.type = 'two-sided', alpha = 0.7)$sign, TRUE ) expect_equal( @@ -240,5 +242,19 @@ DiffCorr(exp2, obs2, ref2)$diff, tolerance = 0.0001 ) +# obs2_2 +expect_equal( +DiffCorr(exp2, obs2_2, ref2, test.type = 'two-sided')$p, +0.4524316, +tolerance = 0.0001 +) +expect_equal( +DiffCorr(exp2, obs2_2, ref2, test.type = 'one-sided')$p, +0.5475684, +tolerance = 0.0001 +) + + + }) # suppressWarnings }) -- GitLab