Hi @vagudets, @nperez, @nmilders,
I have implemented the option to use the effective number of degrees of freedom (N.eff) to account for time series autocorrelation when assessing the statistical significance of the Skill Scores (MSSS, RMSSS, RPSS and CRPSS) with RandomWalktest.
RandomWalkTest: implemented in the atomic function, and included the parameter's documentation and checks. N.eff can be FALSE (default, the length of the time series will be used for the test, so the autocorrelation will not be accounted for), a number (that will be used in all the cases), or an array with the same dimensions as skill_A except time_dim. (I used FALSE instead of NA to avoid confusion with other functions like DiffCorr, where NA means that N.eff are computed automatically with Eno). If test.type == two.sided.approx, N.eff cannot be used, and a warning will arise. In cases with extremely high autocorrelation, N.eff might be lower than the difference of times that one forecast is better than the other, so binom.test fails (it does not make sense that, when flipping a coin 20 times, it shows one side more than 20 times). In such cases, such a difference in times is replaced with N.eff (Núria's idea, to be confirmed with Paco).
MSSS, RMSSS, RPSS and CRPSS: implemented in the big and atomic functions, and included the parameter's documentation and checks. N.eff can be NA (default, N.eff is automatically computed with Eno), FALSE (the length of the time series will be used for the test, so the autocorrelation will not be accounted for), a number (that will be used in all the cases), or an array with the same dimensions as skill_A except time_dim.
Please let me know if you have any questions or feedback. I have made some basic checks (see below), but Nadia has kindly volunteered to test these developments with real data.
Thank you!
Carlos
library(multiApply)
library(s2dv)
.Eno <- s2dv:::.Eno
.RPS <- s2dv:::.RPS
.CRPS <- s2dv:::.CRPS
.GetProbs <- s2dv:::.GetProbs
gitlab <- '/esarchive/scratch/cdelgado/gitlab/s2dv/R/'
source(paste0(gitlab,'RandomWalkTest.R'))
source(paste0(gitlab,'MSSS.R'))
source(paste0(gitlab,'RMSSS.R'))
source(paste0(gitlab,'RPSS.R'))
source(paste0(gitlab,'CRPSS.R'))
scores_A <- array(data = rnorm(1000), c(year = 50, lat = 2, lon = 3))
scores_B <- array(data = rnorm(1000), c(year = 50, lat = 2, lon = 3))
exp <- array(rnorm(1000), dim = c(year = 50, lat = 2, lon = 3, member = 10))
obs <- array(rnorm(1000), dim = c(year = 50, lat = 2, lon = 3))
N.eff <- array(10, dim = c(lat = 2, lon = 3))
## RandomWalkTest
rw_s2dv <- s2dv::RandomWalkTest(skill_A = scores_A, skill_B = scores_B, time_dim = 'year', pval = T, sign = T, test.type = 'two.sided')
rw_F <- RandomWalkTest(skill_A = scores_A, skill_B = scores_B, time_dim = 'year', pval = T, sign = T, test.type = 'two.sided', N.eff = FALSE)
rw_10 <- RandomWalkTest(skill_A = scores_A, skill_B = scores_B, time_dim = 'year', pval = T, sign = T, test.type = 'two.sided', N.eff = 10)
rw_array <- RandomWalkTest(skill_A = scores_A, skill_B = scores_B, time_dim = 'year', pval = T, sign = T, test.type = 'two.sided', N.eff = N.eff)
identical(rw_s2dv,rw_F) # TRUE
identical(rw_s2dv,rw_10) # FALSE
identical(rw_array,rw_10) # TRUE
## MSSS
msss_s2dv <- s2dv::MSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided')
msss_F <- MSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = FALSE)
msss_NA <- MSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = NA)
msss_10 <- MSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = 10)
msss_array <- MSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = N.eff)
identical(msss_s2dv,msss_F) # TRUE
identical(msss_10,msss_F) # FALSE
identical(msss_10,msss_array) # TRUE
## RMSSS
rmsss_s2dv <- s2dv::RMSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided')
rmsss_F <- RMSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = FALSE)
rmsss_NA <- RMSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = NA)
rmsss_10 <- RMSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = 10)
rmsss_array <- RMSSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', pval = T, sig_method = 'Random Walk', sig_method.type = 'two.sided', N.eff = N.eff)
identical(rmsss_s2dv,rmsss_F) # TRUE
identical(rmsss_10,rmsss_F) # FALSE
identical(rmsss_10,rmsss_array) # TRUE
## RPSS
rpss_s2dv <- s2dv::RPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05)
rpss_F <- RPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = FALSE)
rpss_NA <- RPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = NA)
rpss_10 <- RPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = 10)
rpss_array <- RPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = N.eff)
identical(rpss_s2dv,rpss_F) # TRUE
identical(rpss_10,rpss_F) # FALSE
identical(rpss_10,rpss_array) # TRUE
## CRPSS
crpss_s2dv <- s2dv::CRPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05)
crpss_F <- CRPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = FALSE)
crpss_NA <- CRPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = NA)
crpss_10 <- CRPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = 10)
crpss_array <- CRPSS(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', sig_method.type = 'two.sided', alpha = 0.05, N.eff = N.eff)
identical(crpss_s2dv,crpss_F) # TRUE
identical(crpss_10,crpss_F) # FALSE
identical(crpss_10,crpss_array) # TRUE