diff --git a/NAMESPACE b/NAMESPACE
index 9214a1a54c1b7c4bc7e1e81203b51ac5e674cc26..587aec211d4a7bb13e0ef6e9e7dd9a4a99f6eadd 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -77,6 +77,7 @@ export(Season)
 export(SignalNoiseRatio)
 export(Smoothing)
 export(Spectrum)
+export(SprErr)
 export(Spread)
 export(StatSeasAtlHurr)
 export(TPI)
diff --git a/R/SprErr.R b/R/SprErr.R
new file mode 100644
index 0000000000000000000000000000000000000000..f89d37a10fd41dd7f19ec03a848d5b8daa39b639
--- /dev/null
+++ b/R/SprErr.R
@@ -0,0 +1,220 @@
+#'Compute the ratio between the ensemble spread and RMSE
+#'
+#'Compute the ratio between the spread of the members around the 
+#'ensemble mean in experimental data and the RMSE between the ensemble mean of 
+#'experimental and observational data. The p-value and/or the statistical 
+#'significance is provided by a two-sided Fisher's test.
+#'
+#'@param exp A named numeric array of experimental data with at least two 
+#'  dimensions 'memb_dim' and 'time_dim'.
+#'@param obs A named numeric array of observational data with at least two 
+#'  dimensions 'memb_dim' and 'time_dim'. It should have the same dimensions as
+#'  parameter 'exp' except along 'dat_dim' and 'memb_dim'.
+#'@param dat_dim A character string indicating the name of dataset (nobs/nexp) 
+#'  dimension. The default value is NULL (no dataset).
+#'@param memb_dim A character string indicating the name of the member 
+#'  dimension. It must be one dimension in 'exp' and 'obs'. The default value 
+#'  is 'member'.
+#'@param time_dim A character string indicating the name of dimension along  
+#'  which the ratio is computed. The default value is 'sdate'.
+#'@param pval A logical value indicating whether to compute the p-value 
+#'  of the test Ho : SD/RMSE = 1 or not. The default value is TRUE.
+#'@param sign A logical value indicating whether to retrieve the statistical
+#'  significance of the test Ho: ACC = 0 based on 'alpha'. The default value is
+#'  FALSE.
+#'@param alpha A numeric indicating the significance level for the statistical
+#'  significance test. The default value is 0.05. 
+#'@param na.rm A logical value indicating whether to remove NA values. The 
+#'  default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return A list of two arrays with dimensions c(nexp, nobs, the rest of 
+#'  dimensions of 'exp' and 'obs' except memb_dim and time_dim), which nexp is
+#'  the length of dat_dim of 'exp' and nobs is the length of dat_dim of 'obs'. 
+#'  If dat_dim is NULL, nexp and nobs are omitted. \cr
+#'\item{$ratio}{
+#'  The ratio of the ensemble spread and RMSE.
+#'}
+#'\item{$p_val}{
+#'  The p-value of the two-sided Fisher's test with Ho: Spread/RMSE = 1. Only 
+#'  present if \code{pval = TRUE}.
+#'}
+#'
+#'@examples
+#'exp <- array(rnorm(30), dim = c(lat = 2, sdate = 3, member = 5))
+#'obs <- array(rnorm(30), dim = c(lat = 2, sdate = 3))
+#'sprerr1 <- SprErr(exp, obs)
+#'sprerr2 <- SprErr(exp, obs, pval = FALSE, sign = TRUE)
+#'sprerr3 <- SprErr(exp, obs, pval = TRUE, sign = TRUE)
+#'
+#'@import multiApply
+#'@export
+SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member', 
+                   time_dim = 'sdate', pval = TRUE, sign = FALSE, 
+                   alpha = 0.05, na.rm = FALSE, ncores = NULL) {
+  
+  # Check inputs 
+  ## exp and obs (1)
+  if (is.null(exp) | is.null(obs)) {
+    stop("Parameter 'exp' and 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(exp) | !is.numeric(obs)) {
+    stop("Parameter 'exp' and 'obs' must be a numeric array.")
+  }
+  if (is.null(dim(exp)) | is.null(dim(obs))) {
+    stop(paste0("Parameter 'exp' and 'obs' must be array with as least two ",
+                "dimensions memb_dim and time_dim."))
+  }
+  if (any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
+      any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  ## dat_dim
+  if (!is.null(dat_dim)) {
+    if (!is.character(dat_dim) | length(dat_dim) > 1) {
+      stop("Parameter 'dat_dim' must be a character string.")
+    }
+    if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
+      stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
+    }
+  }
+  ## memb_dim
+  if (!is.character(memb_dim) | length(memb_dim) > 1) {
+    stop("Parameter 'memb_dim' must be a character string.")
+  }
+  if (!memb_dim %in% names(dim(exp))) {
+    stop("Parameter 'memb_dim' is not found in 'exp' dimensions. ",
+         "'exp' must have the member dimension to compute the spread.")
+  }
+  # Add [member = 1] 
+  if (memb_dim %in% names(dim(exp)) & !memb_dim %in% names(dim(obs))) {
+    dim(obs) <- c(dim(obs), 1)
+    names(dim(obs))[length(dim(obs))] <- memb_dim
+  }
+  if (!memb_dim %in% names(dim(exp)) & memb_dim %in% names(dim(obs))) { ## check no longer needed?
+    dim(exp) <- c(dim(exp), 1)
+    names(dim(exp))[length(dim(exp))] <- memb_dim
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) > 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
+    stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  if (!is.null(dat_dim)) {
+    name_exp <- name_exp[-which(name_exp == dat_dim)]
+    name_obs <- name_obs[-which(name_obs == dat_dim)]
+  }
+  name_exp <- name_exp[-which(name_exp == memb_dim)]
+  name_obs <- name_obs[-which(name_obs == memb_dim)]
+  if (!identical(dim(exp)[name_exp], dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all the dimensions except 'dat_dim' and 'memb_dim'."))
+  }
+  ## pval
+  if (!is.logical(pval) | length(pval) > 1) {
+    stop("Parameter 'pval' must be one logical value.")
+  }
+  ## sign
+  if (!is.logical(sign) | length(sign) > 1) {
+    stop("Parameter 'sign' must be one logical value.")
+  }
+  # alpha
+  if (!is.numeric(alpha) | any(alpha < 0) | any(alpha > 1) | length(alpha) > 1) {
+    stop("Parameter 'alpha' must be a numeric number between 0 and 1.")
+  }
+  # na.rm
+  if (!na.rm %in% c(TRUE, FALSE)) {
+    stop("Parameter 'na.rm' must be TRUE or FALSE")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  } 
+  
+  
+  ###############################
+  # Calculate RatioSDRMS
+  
+  # If dat_dim = NULL, insert dat dim
+  remove_dat_dim <- FALSE
+  if (is.null(dat_dim)) {
+    dat_dim <- 'dataset'
+    exp <- InsertDim(exp, posdim = 1, lendim = 1, name = 'dataset')
+    obs <- InsertDim(obs, posdim = 1, lendim = 1, name = 'dataset')
+    remove_dat_dim <- TRUE
+  }
+  
+  res <- Apply(list(exp, obs), 
+               target_dims = list(c(dat_dim, memb_dim, time_dim), 
+                                  c(dat_dim, memb_dim, time_dim)),
+               pval = pval,  
+               sign = sign,
+               alpha = alpha,
+               na.rm = na.rm,
+               fun = .SprErr, 
+               ncores = ncores)
+  
+  if (remove_dat_dim) {
+    if (length(dim(res[[1]])) > 2) {
+      res <- lapply(res, Subset, c('nexp', 'nobs'), list(1, 1), drop = 'selected')
+    } else {
+      res <- lapply(res, as.vector)
+    }
+  }
+  
+  return(res)
+}
+
+.SprErr <- function(exp, obs, pval = TRUE, sign = FALSE, alpha = 0.05, na.rm = FALSE) {
+  
+  # exp: [dat_exp, member, sdate]
+  # obs: [dat_obs, member, sdate]
+  nexp <- dim(exp)[1]
+  nobs <- dim(obs)[1]
+  
+  # ensemble mean
+  ens_exp <- MeanDims(exp, 2, na.rm = na.rm) # [dat, sdate]
+  ens_obs <- MeanDims(obs, 2, na.rm = na.rm)
+  
+  # Create empty arrays
+  ratio <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
+  p.val <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
+  
+  for (jexp in 1:nexp) {
+    for (jobs in 1:nobs) {
+      
+      # spread and error
+      spread <- sqrt(mean(apply(exp[jexp,,], 2, var, na.rm = na.rm), na.rm = na.rm))
+      error <- sqrt(mean((ens_obs - ens_exp[jexp,])^2, na.rm = na.rm))
+      ratio[jexp, jobs] <- spread/error
+      
+      # effective sample size
+      enospr <- sum(Eno(apply(exp[jexp,,], 2, var, na.rm = na.rm), names(dim(exp))[3]))
+      enodif <- .Eno((ens_exp[jexp, ] - ens_obs[jobs, ])^2, na.action = na.pass)
+      if (pval | sign) {
+        f_statistic <- (enospr * spread^2 / (enospr - 1)) / (enodif * error^2 / (enodif - 1))
+        if (!is.na(f_statistic) & !is.na(enospr) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
+          p.val[jexp, jobs] <- pf(f_statistic, enospr - 1, enodif - 1)
+          p.val[jexp, jobs] <- 2 * min(p.val[jexp, jobs], 1 - p.val[jexp, jobs])
+        } else {
+          p.val[jexp, jobs] <- NA
+        }
+      }
+    }
+  }
+  
+  res <- list(ratio = ratio)
+  if (pval) {res$p.val <- p.val}
+  if (sign) {res$sign <- p.val <= alpha}
+  
+  return(res)
+}
diff --git a/man/SprErr.Rd b/man/SprErr.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..cdc647eb86373a69c6ff2fedd4a1b15dbc90b584
--- /dev/null
+++ b/man/SprErr.Rd
@@ -0,0 +1,80 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/SprErr.R
+\name{SprErr}
+\alias{SprErr}
+\title{Compute the ratio between the ensemble spread and RMSE}
+\usage{
+SprErr(
+  exp,
+  obs,
+  dat_dim = NULL,
+  memb_dim = "member",
+  time_dim = "sdate",
+  pval = TRUE,
+  sign = FALSE,
+  alpha = 0.05,
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named numeric array of experimental data with at least two 
+dimensions 'memb_dim' and 'time_dim'.}
+
+\item{obs}{A named numeric array of observational data with at least two 
+dimensions 'memb_dim' and 'time_dim'. It should have the same dimensions as
+parameter 'exp' except along 'dat_dim' and 'memb_dim'.}
+
+\item{dat_dim}{A character string indicating the name of dataset (nobs/nexp) 
+dimension. The default value is NULL (no dataset).}
+
+\item{memb_dim}{A character string indicating the name of the member 
+dimension. It must be one dimension in 'exp' and 'obs'. The default value 
+is 'member'.}
+
+\item{time_dim}{A character string indicating the name of dimension along
+which the ratio is computed. The default value is 'sdate'.}
+
+\item{pval}{A logical value indicating whether to compute the p-value 
+of the test Ho : SD/RMSE = 1 or not. The default value is TRUE.}
+
+\item{sign}{A logical value indicating whether to retrieve the statistical
+significance of the test Ho: ACC = 0 based on 'alpha'. The default value is
+FALSE.}
+
+\item{alpha}{A numeric indicating the significance level for the statistical
+significance test. The default value is 0.05.}
+
+\item{na.rm}{A logical value indicating whether to remove NA values. The 
+default value is TRUE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A list of two arrays with dimensions c(nexp, nobs, the rest of 
+ dimensions of 'exp' and 'obs' except memb_dim and time_dim), which nexp is
+ the length of dat_dim of 'exp' and nobs is the length of dat_dim of 'obs'. 
+ If dat_dim is NULL, nexp and nobs are omitted. \cr
+\item{$ratio}{
+ The ratio of the ensemble spread and RMSE.
+}
+\item{$p_val}{
+ The p-value of the two-sided Fisher's test with Ho: Spread/RMSE = 1. Only 
+ present if \code{pval = TRUE}.
+}
+}
+\description{
+Compute the ratio between the spread of the members around the 
+ensemble mean in experimental data and the RMSE between the ensemble mean of 
+experimental and observational data. The p-value and/or the statistical 
+significance is provided by a two-sided Fisher's test.
+}
+\examples{
+exp <- array(rnorm(30), dim = c(lat = 2, sdate = 3, member = 5))
+obs <- array(rnorm(30), dim = c(lat = 2, sdate = 3))
+sprerr1 <- SprErr(exp, obs)
+sprerr2 <- SprErr(exp, obs, pval = FALSE, sign = TRUE)
+sprerr3 <- SprErr(exp, obs, pval = TRUE, sign = TRUE)
+
+}
diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
new file mode 100644
index 0000000000000000000000000000000000000000..65ff728388ff671e2230fa0876396fc01da8de4c
--- /dev/null
+++ b/tests/testthat/test-SprErr.R
@@ -0,0 +1,597 @@
+library(s2dv)
+library(testthat)
+library(multiApply)
+library(ClimProjDiags)
+
+
+##############################################
+# data
+##############################################
+
+# dat1
+set.seed(1)
+exp1 <- array(rnorm(60), dim = c(member = 3, sdate = 10, lat = 2))
+set.seed(2)
+obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+
+# dat1_2
+exp1_2 <-  exp1 + 1
+obs1_2 <-  obs1 + 1
+
+# dat1_3: NAs
+exp1_3 <- exp1; exp1_3[1, 2, 1] <- NA
+obs1_3 <- obs1; obs1_3[2, 1] <- NA
+
+# dat2
+set.seed(1)
+exp2 <- array(rnorm(20), dim = c(member = 2, sdate = 10))
+set.seed(2)
+obs2 <- array(rnorm(10), dim = c(sdate = 10))
+set.seed(2)
+obs2_1 <- array(rnorm(10), dim = c(member = 1, sdate = 10))
+
+# dat3
+set.seed(1)
+exp3 <- array(rnorm(60), dim = c(member = 2, sdate = 10, dataset = 3))
+set.seed(2)
+obs3 <- array(rnorm(20), dim = c(member = 1, sdate = 10, dataset = 2))
+
+# dat3_2
+set.seed(1)
+exp3_2 <- array(rnorm(80), dim = c(member = 4, sdate = 5, dataset = 4))
+set.seed(2)
+obs3_2 <- array(rnorm(30), dim = c(member = 2, sdate = 5, dataset = 3))
+
+# dat4
+set.seed(1)
+exp4 <- array(rnorm(60), dim = c(member = 2, sdate = 10, dataset = 3, lat = 2))
+set.seed(2)
+obs4 <- array(rnorm(20), dim = c(member = 1, sdate = 10, dataset = 2, lat = 2))
+
+# dat4_2: NAs
+exp4_2 <- exp4; exp4_2[1, 2, 1, 1] <- NA
+obs4_2 <- obs4; obs4_2[1, 1:4, 1, 1] <- NA
+
+
+##############################################
+# tests
+##############################################
+
+test_that("1. Input checks", {
+
+  # exp and obs (1)
+  expect_error(
+    SprErr(c(), obs1), 
+    "Parameter 'exp' and 'obs' cannot be NULL." 
+  )
+  expect_error(
+    SprErr(obs1, c()), 
+    "Parameter 'exp' and 'obs' cannot be NULL." 
+  )
+  expect_error(
+    SprErr("", obs1),
+    "Parameter 'exp' and 'obs' must be a numeric array."
+  )
+  expect_error(
+    SprErr(exp1, ""), 
+    "Parameter 'exp' and 'obs' must be a numeric array."
+  )
+  expect_error(
+    SprErr(1, obs1), 
+    "Parameter 'exp' and 'obs' must be array with as least two dimensions memb_dim and time_dim."
+  )
+  expect_error(
+    SprErr(exp1, 1), 
+    "Parameter 'exp' and 'obs' must be array with as least two dimensions memb_dim and time_dim." 
+  )
+  expect_error(
+    SprErr(array(1), obs1), 
+    "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  expect_error(
+    SprErr(exp1, array(1)), 
+    "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  # dat_dim
+  expect_error(
+    SprErr(exp1, obs1, dat_dim = 1),
+    "Parameter 'dat_dim' must be a character string."
+  ) 
+  expect_error(
+    SprErr(exp1, obs1, dat_dim = 'dat'), 
+    "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
+  )  
+  # memb_dim
+  expect_error(
+    SprErr(exp1, obs1, memb_dim = 1),
+    "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+    SprErr(exp1, obs1, memb_dim = 'check'), 
+    "Parameter 'memb_dim' is not found in 'exp' dimensions. 'exp' must have the member dimension to compute the spread."
+  )
+  # time_dim
+  expect_error(
+    SprErr(exp1, obs1, time_dim = 1),
+    "Parameter 'time_dim' must be a character string."
+  )  
+  expect_error(
+    SprErr(exp1, obs1, time_dim = 'time'),
+    "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # exp and obs (2)
+  expect_error(
+    SprErr(exp1, array(1:9, dim = c(sdate = 9))),
+    "Parameter 'exp' and 'obs' must have same length of all the dimensions except 'dat_dim' and 'memb_dim'."
+  )
+  # pval
+  expect_error(
+    SprErr(exp1, obs1, pval = 1),
+    "Parameter 'pval' must be one logical value."
+  )  
+  # sign
+  expect_error(
+    SprErr(exp1, obs1, sign = 1),
+    "Parameter 'sign' must be one logical value."
+  )  
+  # alpha
+  expect_error(
+    SprErr(exp1, obs1, alpha = -0.05),
+    "Parameter 'alpha' must be a numeric number between 0 and 1."
+  )   
+  # na.rm
+  expect_error(
+    SprErr(exp1, obs1, na.rm = ""), 
+    "Parameter 'na.rm' must be TRUE or FALSE"
+  )  
+  # ncores
+  expect_error(
+    SprErr(exp2, obs2, ncores = 1.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+})
+
+##############################################
+
+
+test_that("2. Output checks: dat1", {
+  
+  # element names
+  expect_equal(
+    names(SprErr(exp1, obs1)),
+    c("ratio", "p.val")
+  )
+  expect_equal(
+    names(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)),
+    c("ratio", "sign")
+  )
+  expect_equal(
+    names(SprErr(exp1, obs1, pval = TRUE, sign = TRUE)),
+    c("ratio", "p.val", "sign")
+  )
+  # dimensions
+  expect_equal(
+    dim(SprErr(exp1, obs1)$ratio),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp1, obs1)$p.val),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$ratio), 
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$sign), 
+    c(lat = 2)
+  )
+  # values
+  expect_equal(
+    as.vector(SprErr(exp1, obs1)$ratio),
+    c(1.0646692, 0.6657522),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1)$p.val),
+    c(0.8549593, 0.2412730),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE)$sign),
+    c(FALSE, FALSE)
+  )
+  # pval = FALSE
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$ratio),
+    c(01.0646692, 0.6657522),  
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$sign),
+    c(FALSE, FALSE)  
+  )
+  # na.rm = TRUE
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, na.rm = TRUE)$ratio),
+    c(1.0646692, 0.6657522),  
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, na.rm = TRUE)$p.val),
+    c(0.8549593, 0.2412730),  
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE, na.rm = TRUE)$sign),
+    c(FALSE, FALSE)
+  )
+  # alpha
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE, alpha = 0.5)$sign),
+    c(FALSE, TRUE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE, alpha = 0.99)$sign),
+    c(TRUE, TRUE)
+  )
+  
+  # dat1_2
+  expect_equal(
+    SprErr(exp1, obs1),
+    SprErr(exp1_2, obs1_2)
+  )
+  expect_equal(
+    SprErr(exp1, obs1, sign = TRUE)$ratio,
+    SprErr(exp1_2, obs1_2)$ratio
+  )
+  expect_equal(
+    SprErr(exp1, obs1, sign = TRUE)$p.val,
+    SprErr(exp1_2, obs1_2)$p.val
+  )
+  
+  # dat1_3
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3)$ratio),
+    c(NA, 0.6657522),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3)$p.val),
+    c(NA, 0.241273),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, sign = TRUE)$sign),
+    c(NA, FALSE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, na.rm = TRUE)$ratio),
+    c(0.9656329, 0.6657522),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, na.rm = TRUE)$p.val),
+    c(0.896455, 0.241273),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, sign = TRUE, na.rm = TRUE)$sign),
+    c(FALSE, FALSE)
+  )
+  
+})
+
+##############################################
+
+
+test_that("3. Output checks: dat2", {
+  
+  expect_equal(
+    names(SprErr(exp2, obs2)),
+    c("ratio", "p.val")
+  )
+  expect_equal(
+    names(SprErr(exp2, obs2, sign = TRUE)),
+    c("ratio", "p.val","sign")
+  )
+  expect_equal(
+    dim(SprErr(exp2, obs2)$ratio),
+    NULL
+  )
+  expect_equal(
+    dim(SprErr(exp2, obs2)$p.val),
+    NULL
+  )
+  expect_equal(
+    dim(SprErr(exp2, obs2, sign = TRUE)$sign),
+    NULL
+  )
+  # values
+  expect_equal(
+    as.vector(SprErr(exp2, obs2)$ratio),
+    c(0.6866402),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp2, obs2)$p.val),
+    c(0.2779936),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp2, obs2, sign = TRUE)$sign), 
+    FALSE 
+  )
+  # sign = TRUE
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, sign = TRUE)$ratio),
+    c(0.6866402),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, sign = TRUE)$p.val),
+    c(0.2779936),
+    tolerance = 0.0001
+  )
+  # alpha
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, sign = TRUE, alpha = 0.99)$sign),
+    TRUE 
+  )
+  # na.rm = TRUE
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, na.rm = TRUE)$ratio),
+    c(0.6866402),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, na.rm = TRUE)$p.val),
+    c(0.2779936),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, na.rm = TRUE, sign = TRUE)$sign), 
+    FALSE 
+  )
+  # other
+  expect_equal(
+    SprErr(exp2, obs2),
+    SprErr(exp2, obs2_1)
+  )
+  
+})
+
+##############################################
+
+
+test_that("4. Output checks: dat3", {
+  
+  expect_equal(
+    dim(SprErr(exp3, obs3, dat_dim = 'dataset')$ratio),
+    c('nexp' = 3, 'nobs' = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp3, obs3, dat_dim = 'dataset')$p.val),
+    c('nexp' = 3, 'nobs' = 2)
+  )
+  # values
+  expect_equal(
+    mean(SprErr(exp3, obs3, dat_dim = 'dataset')$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset')$ratio),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset')$p.val)[1:3],
+    c(0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE)$sign)[1:3],
+    c(FALSE, FALSE, TRUE)
+  )
+  expect_equal(
+    mean(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE)$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE)$p.val)[1:3],
+    c(0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  # alpha
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE, alpha = 0.99)$sign)[1:3],
+    c(TRUE, TRUE, TRUE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE, alpha = 0.20)$sign)[1:3],
+    c(FALSE, TRUE, TRUE)
+  )
+  # na.rm = TRUE
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', na.rm = TRUE)$ratio),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', na.rm = TRUE)$p.val)[1:3],
+    c(0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', na.rm = TRUE, sign = TRUE)$sign)[1:3], 
+    c(FALSE, FALSE, TRUE) 
+  )
+  
+  # dat3_2
+  expect_equal(
+    dim(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$ratio),
+    c('nexp' = 4, 'nobs' = 3)
+  )
+  expect_equal(
+    dim(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$p.val),
+    c('nexp' = 4, 'nobs' = 3)
+  )
+  # values
+  expect_equal(
+    mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$ratio),
+    c(1.25586),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$p.val)[1:4],
+    c(0.6927309, 0.7390035, 0.8834023, 0.4421531),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE)$sign)[1:4],
+    c(FALSE, FALSE, FALSE, FALSE)
+  )
+  expect_equal(
+    mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE)$ratio),
+    c(1.25586),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE)$p.val)[1:4],
+    c(0.6927309, 0.7390035, 0.8834023, 0.4421531),
+    tolerance = 0.0001
+  )
+  # alpha
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE, alpha = 0.99)$sign)[1:3],
+    c(TRUE, TRUE, TRUE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE, alpha = 0.70)$sign)[1:4],
+    c(TRUE, FALSE, FALSE, TRUE)
+  )
+  # na.rm = TRUE
+  expect_equal(
+    mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', na.rm = TRUE)$ratio),
+    c(1.25586),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', na.rm = TRUE)$p.val)[1:4],
+    c(0.6927309, 0.7390035, 0.8834023, 0.4421531),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', na.rm = TRUE, 
+                     alpha = 0.70, sign = TRUE)$sign)[1:4], 
+    c(TRUE, FALSE, FALSE, TRUE) 
+  )
+  
+})
+
+############################################## 
+
+
+test_that("5. Output checks: dat4", {
+  
+  expect_equal(
+    dim(SprErr(exp4, obs4, dat_dim = 'dataset')$ratio),
+    c('nexp' = 3, 'nobs' = 2, 'lat' = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp4, obs4, dat_dim = 'dataset', sign = TRUE)$p.val),
+    c('nexp' = 3, 'nobs' = 2, 'lat' = 2)
+  )
+  # values
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset')$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset')$p.val),
+    c(0.1674805),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[, , 2]),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4, obs4, dat_dim = 'dataset')$p.val[, , 2]),
+    c(0.30405979, 0.18162950, 0.01675207, 0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[, , 1],
+    SprErr(exp3, obs3, dat_dim = 'dataset')$ratio
+  )
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset', sign = TRUE)$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset', sign = TRUE)$p.val),
+    c(0.1674805),
+    tolerance = 0.0001
+  )
+  
+  # dat4_2: NAs
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$ratio[, , 1]),
+    c(NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$p.val[, , 1]),
+    c(NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', sign = TRUE)$sign[, , 1]),
+    c(NA, NA, NA, NA, NA, NA),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$ratio[, , 2]),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$p.val[, , 2]),
+    c(0.30405979, 0.18162950, 0.01675207, 0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$ratio)),
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$p.val))
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$ratio[, , 1]),
+    c(0.4648097, 0.6571888, 0.3975804, 0.4648097, 0.6571888, 0.3975804),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$p.val[, , 1]),
+    c(0.04674814, 0.21983481, 0.01350139, 0.04351013, 0.22700305, 0.01124875),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm =TRUE, sign = TRUE)$sign[, , 1]),
+    c(TRUE, FALSE, TRUE, TRUE, FALSE, TRUE),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$ratio)),
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$p.val))
+  )
+  expect_equal(
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$ratio)),
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$p.val))
+  )
+  
+})