From c6d7e525868e0fb00dfda2a52e8bfc80948a81f0 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Fri, 19 Jun 2020 15:48:22 +0200 Subject: [PATCH 1/8] added the RandomWalkTest function --- R/RandomWalkTest.R | 72 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 R/RandomWalkTest.R diff --git a/R/RandomWalkTest.R b/R/RandomWalkTest.R new file mode 100644 index 0000000..c39062b --- /dev/null +++ b/R/RandomWalkTest.R @@ -0,0 +1,72 @@ +#'Random walk test for skill differences +#' +#'@description Forecast comparison of 2 forecasts based on Random Walks (given skill of the 2 forecasts with respect to a common reference), +#'with significance estimate at the 5% confidence level, as in DelSole and Tippett (2015). +#' +#'@author Andrea Manrique \email{andrea.manrique@bsc.es} +#'@author Carlos Delgado-Torres \email{carlos.delgado@bsc.es} +#' +#'@param skill_A Time serie of the skill with the forecaster A's. +#'@param skill_B Time serie of the skill with the forecaster B's. +#'@param time_dim A character string indicating the name of the dimension along which the tests are computed. The default value is 'sdate'. +#'@param ncores An integer indicating the number of cores to use for parallel computation. The default value is NULL. +#' +#'@return A list of 2: +#'\item{$score}{ +#'The number of times that forecaster A has been better than forecaster B minus +#'the number of times that forecaster B has been better than forecaster A (for skill positively oriented). +#'If $score is positive forecaster A is better than forecaster B, and if $score is negative forecaster B is better than forecaster B. +#'} +#'\item{$signif}{ +#'Whether the difference is significant or not at the 5% significance level. +#'} +#' +#'@examples +#' dat1 <- array(c(1:40), dim = c(sdate = 10, lat = 2, lon = 2)) +#' dat2 <- array(c(11:50), dim = c(sdate = 10, lat = 2, lon = 2)) +#' print(RandomWalkTest(skill_A = dat1, skill_B = dat2, time_dim = 'sdate', ncores = 1)) +#' +#'@import multiApply +#'@export +RandomWalkTest <- function(skill_A, skill_B, time_dim = 'sdate', ncores = NULL){ + + ## Check inputs + if (is.null(skill_A) | is.null(skill_B)){ + stop("Parameters 'skill_A' and 'skill_B' cannot be NULL.") + } + if(!is.numeric(skill_A) | !is.numeric(skill_B)){ + stop("Parameter 'skill_A' and 'skill_B' must be a numerical array.") + } + if (!identical(dim(skill_A),dim(skill_B))){ + stop("Parameters 'skill_A' and 'skill_B' must have the same dimensions.") + } + if(!is.character(time_dim)){ + stop("Parameter 'time_dim' must be a character string.") + } + if(!time_dim %in% names(dim(skill_A)) | !time_dim %in% names(dim(skill_B))){ + stop("Parameter 'time_dim' is not found in 'skill_A' or 'skill_B' dimensions.") + } + if (!is.null(ncores)){ + if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores < 0 | length(ncores) > 1){ + stop("Parameter 'ncores' must be a positive integer.") + } + } + + ## Compute the Random Walk Test + res <- multiApply::Apply(data = list(skill_A, skill_B), + target_dims = time_dim, + fun = .RandomWalkTest, + ncores = ncores) + return(res) +} + +.RandomWalkTest <- function(skill_A, skill_B){ + + score <- cumsum(skill_A > skill_B) - cumsum(skill_A < skill_B) + + ## TRUE if significant (if last value is above or below 2*sqrt(N)) + signif<- ifelse(test = (score[length(skill_A)] < (-2)*sqrt(length(skill_A))) | (score[length(skill_A)] > 2*sqrt(length(skill_A))), + yes = TRUE, no = FALSE) + + return(list("score"=score[length(skill_A)],"signif"=signif)) +} \ No newline at end of file -- GitLab From c7170c32e8c9c057fbe169c23bde3d86a101c8f3 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Thu, 23 Jul 2020 16:08:20 +0200 Subject: [PATCH 2/8] improved documentation --- R/RandomWalkTest.R | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/R/RandomWalkTest.R b/R/RandomWalkTest.R index c39062b..724dbb3 100644 --- a/R/RandomWalkTest.R +++ b/R/RandomWalkTest.R @@ -1,6 +1,6 @@ #'Random walk test for skill differences #' -#'@description Forecast comparison of 2 forecasts based on Random Walks (given skill of the 2 forecasts with respect to a common reference), +#'@description Forecast comparison of the skill obtained with 2 forecasts (and using the same reference) based on Random Walks (given skill of the 2 forecasts with respect to a common reference), #'with significance estimate at the 5% confidence level, as in DelSole and Tippett (2015). #' #'@author Andrea Manrique \email{andrea.manrique@bsc.es} @@ -22,9 +22,12 @@ #'} #' #'@examples -#' dat1 <- array(c(1:40), dim = c(sdate = 10, lat = 2, lon = 2)) -#' dat2 <- array(c(11:50), dim = c(sdate = 10, lat = 2, lon = 2)) -#' print(RandomWalkTest(skill_A = dat1, skill_B = dat2, time_dim = 'sdate', ncores = 1)) +#' fcst_A <- array(c(11:50), dim = c(sdate = 10, lat = 2, lon = 2)) +#' fcst_B <- array(c(21:60), dim = c(sdate = 10, lat = 2, lon = 2)) +#' reference <- array(1:40, dim = c(sdate = 10, lat = 2, lon = 2)) +#' skill_A <- abs(fcst_A - reference) +#' skill_B <- abs(fcst_B - reference) +#' RandomWalkTest(skill_A = skill_A, skill_B = skill_B, time_dim = 'sdate', ncores = 1) #' #'@import multiApply #'@export -- GitLab From 3eb9f6a4d574664215441183cb971656e4ada5bf Mon Sep 17 00:00:00 2001 From: cdelgado Date: Thu, 23 Jul 2020 16:15:26 +0200 Subject: [PATCH 3/8] improved doc --- R/RandomWalkTest.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/RandomWalkTest.R b/R/RandomWalkTest.R index 724dbb3..230112f 100644 --- a/R/RandomWalkTest.R +++ b/R/RandomWalkTest.R @@ -1,6 +1,6 @@ #'Random walk test for skill differences #' -#'@description Forecast comparison of the skill obtained with 2 forecasts (and using the same reference) based on Random Walks (given skill of the 2 forecasts with respect to a common reference), +#'@description Forecast comparison of the skill obtained with 2 forecasts (using the same reference) based on Random Walks (given the skill of the 2 forecasts with respect to a common reference), #'with significance estimate at the 5% confidence level, as in DelSole and Tippett (2015). #' #'@author Andrea Manrique \email{andrea.manrique@bsc.es} -- GitLab From 44340fca565080a8fb35091547bbb722b2265f5b Mon Sep 17 00:00:00 2001 From: cdelgado Date: Thu, 23 Jul 2020 16:22:27 +0200 Subject: [PATCH 4/8] improved documentation --- R/RandomWalkTest.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/RandomWalkTest.R b/R/RandomWalkTest.R index 230112f..0ad8ea1 100644 --- a/R/RandomWalkTest.R +++ b/R/RandomWalkTest.R @@ -1,6 +1,6 @@ #'Random walk test for skill differences #' -#'@description Forecast comparison of the skill obtained with 2 forecasts (using the same reference) based on Random Walks (given the skill of the 2 forecasts with respect to a common reference), +#'@description Forecast comparison of the skill obtained with 2 forecasts (with recpect to a common reference) based on Random Walks, #'with significance estimate at the 5% confidence level, as in DelSole and Tippett (2015). #' #'@author Andrea Manrique \email{andrea.manrique@bsc.es} -- GitLab From 6d04dbf05195aa2f10500e67607fb3a7e127343b Mon Sep 17 00:00:00 2001 From: cdelgado Date: Thu, 23 Jul 2020 16:22:50 +0200 Subject: [PATCH 5/8] improved doc --- R/RandomWalkTest.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/RandomWalkTest.R b/R/RandomWalkTest.R index 0ad8ea1..31c2e61 100644 --- a/R/RandomWalkTest.R +++ b/R/RandomWalkTest.R @@ -1,6 +1,6 @@ #'Random walk test for skill differences #' -#'@description Forecast comparison of the skill obtained with 2 forecasts (with recpect to a common reference) based on Random Walks, +#'@description Forecast comparison of the skill obtained with 2 forecasts (with respect to a common reference) based on Random Walks, #'with significance estimate at the 5% confidence level, as in DelSole and Tippett (2015). #' #'@author Andrea Manrique \email{andrea.manrique@bsc.es} -- GitLab From efb3c2cbb94210d4eb43628d1f2db8eb92e89a79 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 28 Oct 2020 15:01:33 +0100 Subject: [PATCH 6/8] Format improvement of RandomWalkTest.R, and create unit test --- NAMESPACE | 1 + R/RandomWalkTest.R | 40 ++++++++------ man/RandomWalkTest.Rd | 57 ++++++++++++++++++++ tests/testthat/test-RandomWalkTest.R | 81 ++++++++++++++++++++++++++++ 4 files changed, 164 insertions(+), 15 deletions(-) create mode 100644 man/RandomWalkTest.Rd create mode 100644 tests/testthat/test-RandomWalkTest.R diff --git a/NAMESPACE b/NAMESPACE index 3e6f876..b196f73 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -32,6 +32,7 @@ export(PlotSection) export(PlotStereoMap) export(RMS) export(RMSSS) +export(RandomWalkTest) export(Regression) export(Reorder) export(Season) diff --git a/R/RandomWalkTest.R b/R/RandomWalkTest.R index 31c2e61..0771331 100644 --- a/R/RandomWalkTest.R +++ b/R/RandomWalkTest.R @@ -1,24 +1,31 @@ #'Random walk test for skill differences #' -#'@description Forecast comparison of the skill obtained with 2 forecasts (with respect to a common reference) based on Random Walks, -#'with significance estimate at the 5% confidence level, as in DelSole and Tippett (2015). +#'Forecast comparison of the skill obtained with 2 forecasts (with respect to a +#'common reference) based on Random Walks, with significance estimate at the 5% +#'confidence level, as in DelSole and Tippett (2015). #' -#'@author Andrea Manrique \email{andrea.manrique@bsc.es} -#'@author Carlos Delgado-Torres \email{carlos.delgado@bsc.es} -#' -#'@param skill_A Time serie of the skill with the forecaster A's. -#'@param skill_B Time serie of the skill with the forecaster B's. -#'@param time_dim A character string indicating the name of the dimension along which the tests are computed. The default value is 'sdate'. -#'@param ncores An integer indicating the number of cores to use for parallel computation. The default value is NULL. +#'@param skill_A A numerical array of the time series of the skill with the +#' forecaster A's. +#'@param skill_B A numerical array of the time series of the skill with the +#' forecaster B's. The dimensions should be identical as parameter 'skill_A'. +#'@param time_dim A character string indicating the name of the dimension along +#' which the tests are computed. The default value is 'sdate'. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. The default value is NULL. #' #'@return A list of 2: #'\item{$score}{ -#'The number of times that forecaster A has been better than forecaster B minus -#'the number of times that forecaster B has been better than forecaster A (for skill positively oriented). -#'If $score is positive forecaster A is better than forecaster B, and if $score is negative forecaster B is better than forecaster B. +#' A numerical array with the same dimensions as the input arrays except +#' 'time_dim'. The number of times that forecaster A has been better than +#' forecaster B minus the number of times that forecaster B has been better +#' than forecaster A (for skill positively oriented). If $score is positive +#' forecaster A is better than forecaster B, and if $score is negative +#' forecaster B is better than forecaster B. #'} #'\item{$signif}{ -#'Whether the difference is significant or not at the 5% significance level. +#' A logical array with the same dimensions as the input arrays except +#' 'time_dim'. Whether the difference is significant or not at the 5% +#' significance level. #'} #' #'@examples @@ -29,6 +36,9 @@ #' skill_B <- abs(fcst_B - reference) #' RandomWalkTest(skill_A = skill_A, skill_B = skill_B, time_dim = 'sdate', ncores = 1) #' +#'@author Andrea Manrique \email{andrea.manrique@bsc.es} +#'@author Carlos Delgado-Torres \email{carlos.delgado@bsc.es} +#' #'@import multiApply #'@export RandomWalkTest <- function(skill_A, skill_B, time_dim = 'sdate', ncores = NULL){ @@ -38,7 +48,7 @@ RandomWalkTest <- function(skill_A, skill_B, time_dim = 'sdate', ncores = NULL){ stop("Parameters 'skill_A' and 'skill_B' cannot be NULL.") } if(!is.numeric(skill_A) | !is.numeric(skill_B)){ - stop("Parameter 'skill_A' and 'skill_B' must be a numerical array.") + stop("Parameters 'skill_A' and 'skill_B' must be a numerical array.") } if (!identical(dim(skill_A),dim(skill_B))){ stop("Parameters 'skill_A' and 'skill_B' must have the same dimensions.") @@ -72,4 +82,4 @@ RandomWalkTest <- function(skill_A, skill_B, time_dim = 'sdate', ncores = NULL){ yes = TRUE, no = FALSE) return(list("score"=score[length(skill_A)],"signif"=signif)) -} \ No newline at end of file +} diff --git a/man/RandomWalkTest.Rd b/man/RandomWalkTest.Rd new file mode 100644 index 0000000..2946517 --- /dev/null +++ b/man/RandomWalkTest.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RandomWalkTest.R +\name{RandomWalkTest} +\alias{RandomWalkTest} +\title{Random walk test for skill differences} +\usage{ +RandomWalkTest(skill_A, skill_B, time_dim = "sdate", ncores = NULL) +} +\arguments{ +\item{skill_A}{A numerical array of the time series of the skill with the +forecaster A's.} + +\item{skill_B}{A numerical array of the time series of the skill with the +forecaster B's. The dimensions should be identical as parameter 'skill_A'.} + +\item{time_dim}{A character string indicating the name of the dimension along +which the tests are computed. The default value is 'sdate'.} + +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} +} +\value{ +A list of 2: +\item{$score}{ + A numerical array with the same dimensions as the input arrays except + 'time_dim'. The number of times that forecaster A has been better than + forecaster B minus the number of times that forecaster B has been better + than forecaster A (for skill positively oriented). If $score is positive + forecaster A is better than forecaster B, and if $score is negative + forecaster B is better than forecaster B. +} +\item{$signif}{ + A logical array with the same dimensions as the input arrays except + 'time_dim'. Whether the difference is significant or not at the 5% + significance level. +} +} +\description{ +Forecast comparison of the skill obtained with 2 forecasts (with respect to a +common reference) based on Random Walks, with significance estimate at the 5% +confidence level, as in DelSole and Tippett (2015). +} +\examples{ +fcst_A <- array(c(11:50), dim = c(sdate = 10, lat = 2, lon = 2)) +fcst_B <- array(c(21:60), dim = c(sdate = 10, lat = 2, lon = 2)) +reference <- array(1:40, dim = c(sdate = 10, lat = 2, lon = 2)) +skill_A <- abs(fcst_A - reference) +skill_B <- abs(fcst_B - reference) +RandomWalkTest(skill_A = skill_A, skill_B = skill_B, time_dim = 'sdate', ncores = 1) + +} +\author{ +Andrea Manrique \email{andrea.manrique@bsc.es} + +Carlos Delgado-Torres \email{carlos.delgado@bsc.es} +} + diff --git a/tests/testthat/test-RandomWalkTest.R b/tests/testthat/test-RandomWalkTest.R new file mode 100644 index 0000000..d50b9ce --- /dev/null +++ b/tests/testthat/test-RandomWalkTest.R @@ -0,0 +1,81 @@ +context("s2dv::RandomWalkTest tests") + +############################################## + set.seed(1) + dat1_A <- array(rnorm(64), dim = c(sdate = 4, ftime = 4, lat = 2, lon = 2)) + set.seed(2) + dat1_B <- array(rnorm(64), dim = c(sdate = 4, ftime = 4, lat = 2, lon = 2)) + + +############################################## +test_that("1. Input checks", { + + expect_error( + RandomWalkTest(c(), dat1_B), + "Parameters 'skill_A' and 'skill_B' cannot be NULL." + ) + expect_error( + RandomWalkTest(dat1_A, 'a'), + "Parameters 'skill_A' and 'skill_B' must be a numerical array." + ) + expect_error( + RandomWalkTest(array(1:10, dim = c(2, 5)), array(1:12, dim = c(2, 6))), + "Parameters 'skill_A' and 'skill_B' must have the same dimensions." + ) + expect_error( + RandomWalkTest(1:10, 2:11, time_dim = 12), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + RandomWalkTest(dat1_A, dat1_B, time_dim = 'date'), + "Parameter 'time_dim' is not found in 'skill_A' or 'skill_B' dimension." + ) + expect_error( + RandomWalkTest(dat1_A, dat1_B, ncores = T), + "Parameter 'ncores' must be a positive integer." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + res <- RandomWalkTest(dat1_A, dat1_B) + + expect_equal( + length(res), + 2 + ) + expect_equal( + names(res), + c("score", "signif") + ) + expect_equal( + dim(res$score), + c(ftime = 4, lat = 2, lon = 2) + ) + expect_equal( + dim(res$signif), + c(ftime = 4, lat = 2, lon = 2) + ) + expect_equal( + is.numeric(res$score), + TRUE + ) + expect_equal( + is.logical(res$signif), + TRUE + ) + expect_equal( + mean(res$score, na.rm = T), + 0 + ) + expect_equal( + res$score[, 1, 1], + c(0, 0, -2, -2) + ) + expect_equal( + res$score[, 1, 2], + c(0, 4, 2, 0) + ) + +}) -- GitLab From 630ba721cabb9352568defec262c1d67e5d099b7 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 28 Oct 2020 15:05:58 +0100 Subject: [PATCH 7/8] Add the RandomWalkTest authors to DESCRIPTION --- DESCRIPTION | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index f88f28d..9d44334 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,7 +4,9 @@ Version: 0.0.1 Authors@R: c( person("BSC-CNS", role = c("aut", "cph")), person("An-Chi", "Ho", , "an.ho@bsc.es", role = c("aut", "cre")), - person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = "aut")) + person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = "aut")), + person("Andrea", "Manrique", , "andrea.manrique@bsc.es", role = c("ctb")), + person("Carlos", "Delgado", , "carlos.delgado@bsc.es", role = c("ctb")) Description: The advanced version of package 's2dverification'. It is intended for 'seasonal to decadal' (s2d) climate forecast verification, but it can also be used in other kinds of forecasts or general climate analysis. -- GitLab From eda488931f88a1a21bbc77ed1d4297de5e6be6f8 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 28 Oct 2020 15:34:06 +0100 Subject: [PATCH 8/8] Typo fix --- DESCRIPTION | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9d44334..f61512a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,9 +4,9 @@ Version: 0.0.1 Authors@R: c( person("BSC-CNS", role = c("aut", "cph")), person("An-Chi", "Ho", , "an.ho@bsc.es", role = c("aut", "cre")), - person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = "aut")), - person("Andrea", "Manrique", , "andrea.manrique@bsc.es", role = c("ctb")), - person("Carlos", "Delgado", , "carlos.delgado@bsc.es", role = c("ctb")) + person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = "aut"), + person("Andrea", "Manrique", , "andrea.manrique@bsc.es", role = "ctb"), + person("Carlos", "Delgado", , "carlos.delgado@bsc.es", role = "ctb")) Description: The advanced version of package 's2dverification'. It is intended for 'seasonal to decadal' (s2d) climate forecast verification, but it can also be used in other kinds of forecasts or general climate analysis. -- GitLab