diff --git a/DESCRIPTION b/DESCRIPTION index 0045265cfeeefd668ffd0fc634ae9e42cac29915..6780c83d1760c794b0289f671433edacbab07636 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -50,7 +50,8 @@ Imports: ncdf4, parallel, plyr, - SpecsVerification (>= 0.5.0) + SpecsVerification (>= 0.5.0), + multiApply (>= 2.0.0) Suggests: easyVerification, testthat diff --git a/NAMESPACE b/NAMESPACE index 4e8557898444268933dfffceb4f5d7d3c8460183..e8dbfa095eeff6ee5e92dce2d6fa6b88041f0fd7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,7 +6,6 @@ export(.RMS) export(.RMSSS) export(.RatioRMS) export(.RatioSDRMS) -export(.Trend) export(ACC) export(Alpha) export(AnimateMap) @@ -90,6 +89,7 @@ import(graphics) import(mapproj) import(maps) import(methods) +import(multiApply) import(ncdf4) import(parallel) import(plyr) @@ -116,6 +116,7 @@ importFrom(stats,kmeans) importFrom(stats,lm) importFrom(stats,mad) importFrom(stats,median) +importFrom(stats,na.fail) importFrom(stats,na.omit) importFrom(stats,na.pass) importFrom(stats,pf) diff --git a/R/Consist_Trend.R b/R/Consist_Trend.R index 2a82fd9c021ab1028379515fb94bbef325ffd3ff..31155bd6bcae60b927c771262132295fa93e59cb 100644 --- a/R/Consist_Trend.R +++ b/R/Consist_Trend.R @@ -109,8 +109,8 @@ Consist_Trend <- function(var_exp, var_obs, interval = 1) { detrendedobs <- array(dim = dimobs) tmp1 <- Trend(var_exp, 2, interval = interval) tmp2 <- Trend(var_obs, 2, interval = interval) - trendcat[1:dimexp[1], , , , , ] <- tmp1$trend - trendcat[(dimexp[1] + 1):dim(trendcat)[1], , , , , ] <- tmp2$trend + trendcat[1:dimexp[1], , , , , ] <- tmp1$trend[, 1:4, , , , ] + trendcat[(dimexp[1] + 1):dim(trendcat)[1], , , , , ] <- tmp2$trend[, 1:4, , , , ] detrendedmod[] <- tmp1$detrended detrendedobs[] <- tmp2$detrended # diff --git a/R/Eno.R b/R/Eno.R index 877710cb9a1769378ac81c480fcaff5b78459b46..21899e737c5c0e7cf0aaaac07d789fc34b0b01f8 100644 --- a/R/Eno.R +++ b/R/Eno.R @@ -6,15 +6,21 @@ #'statistical/inference tests.\cr #'Based on eno function from Caio Coelho from rclim.txt. #' -#'@param obs Matrix of any number of dimensions up to 10. -#'@param posdim Dimension along which to compute the effective sample size. +#'@param obs A numeric n-dimensional array with or without named dimensions. +#'@param posdim A character or an integer indicating the dimension along +#' which to compute the effective sample size. +#'@param na.action Function to be called to handle missing values, can be +#' na.pass or na.fail. Default: na.pass +#'@param ncores The number of cores to use for parallel computation. Default: NULL #' -#'@return Same dimensions as var but without the posdim dimension. +#'@return An array with the same dimensions as obs except the posdim dimension. #' #'@keywords datagen #'@author History:\cr #'0.1 - 2011-05 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Original code\cr -#'1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN +#'1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN +#'3.0 - 2019-08 (A. Ho, \email{an.ho@bsc.es}) - Adopt multiApply +#' #'@examples #'# See examples on Load() to understand the first lines in this example #' \dontrun{ @@ -47,70 +53,105 @@ #'eno <- Eno(sampleData$mod[1, 1, , 1, , ], 1) #'PlotEquiMap(eno, sampleData$lon, sampleData$lat) #' -#'@importFrom stats acf na.pass +#'@importFrom stats acf na.pass na.fail +#'@import multiApply #'@export -Eno <- function(obs, posdim) { - dimsvar <- dim(obs) - if (is.null(dimsvar)) { - dimsvar <- length(obs) +Eno <- function(obs, posdim, na.action = na.pass, ncores = NULL) { + + #Check params: + + if (is.null(obs)) { + stop("Parameter 'obs' cannot be NULL.") + } + if (!is.numeric(obs)) { + stop("Parameter 'obs' must be a numeric array.") + } + + if (is.null(posdim)) { + stop("Parameter 'posdim' cannot be NULL.") + } + + if (!is.numeric(posdim) && !is.character(posdim)) { + stop("Parameter 'posdim' must be an positive integer or character.") } - enlobs <- Enlarge(obs, 10) - outdim <- c(dimsvar, array(1, dim = (10 - length(dimsvar)))) - posaperm <- 1:10 - posaperm[posdim] <- 1 - posaperm[1] <- posdim - enlobs <- aperm(enlobs, posaperm) - dimsaperm <- outdim[posaperm] - # - # Loop on all dimensions to compute effective number of observations - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - enleno <- array(dim = c(1, dimsaperm[2:10])) - for (j2 in 1:dimsaperm[2]) { - for (j3 in 1:dimsaperm[3]) { - for (j4 in 1:dimsaperm[4]) { - for (j5 in 1:dimsaperm[5]) { - for (j6 in 1:dimsaperm[6]) { - for (j7 in 1:dimsaperm[7]) { - for (j8 in 1:dimsaperm[8]) { - for (j9 in 1:dimsaperm[9]) { - for (j10 in 1:dimsaperm[10]) { - tmp <- enlobs[, j2, j3, j4, j5, j6, j7, j8, j9, j10] - if (length(sort(tmp)) > 1 ) { - n <- length(sort(tmp)) - a <- acf(tmp, lag.max = n - 1, plot = FALSE, - na.action = na.pass)$acf[2:n, 1, 1] - s <- 0 - for (k in 1:(n - 1)) { - s <- s + (((n - k) / n) * a[k]) - } - enleno[1, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- min(n / (1 + (2 * s)), n) - } - } - } - } - } - } - } + + if (is.numeric(posdim)) { + if (posdim <= 0 | posdim %% 1 != 0) { + stop("Parameter 'posdim' must be a positive integer or character.") + } else if (posdim > length(dim(obs))) { + stop("Parameter 'posdim' excesses the dimension length of parameter 'obs'.") + } else { #posdim is okay for being integer. But change it to character. + if (!is.null(names(dim(obs)))) { + posdim <- names(dim(obs))[posdim] + } else { + names(dim(obs)) <- paste0("D", 1:length(dim(obs))) + posdim <- names(dim(obs))[posdim] } } } - # - # Back to the original dimensions - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - #dimsvar <- dimsvar[-posdim] - if (length(dimsvar) == 1) { - dimsvar <- 1 + + if (is.character(posdim)) { + if (is.null(names(dim(obs)))) { + stop(paste0("Parameter 'obs' doesn't have dimension names. ", + "Add names or specify parameter 'posdim' by number.")) + } else if (!any(posdim %in% names(dim(obs)))) { + stop(paste0("Parameter 'posdim' does not match any dimension name of ", + "parameter 'obs'.")) + } + } + + if (as.character(substitute(na.action)) != c("na.pass") && + as.character(substitute(na.action)) != c("na.fail")) { + stop("Parameter 'na.action' must be either 'na.pass' or 'na.fail'.") + } + + if(as.character(substitute(na.action))== c("na.fail") && any(is.na(obs))) { + stop(paste0("NA value in paratemter 'obs', which is not accepted when ", + "parameter 'na.action' = 'na.fail'.")) + } + + if (!is.null(ncores)) { + if (!is.numeric(ncores)) { + stop("Parameter 'ncores' must be a positive integer.") + } else if (ncores <= 0 | ncores %% 1 != 0) { + stop("Parameter 'ncores' must be a positive integer.") + } + } + + #Calculation: + + effnumobs <- Apply(data = list(obs), + target_dims = posdim, + fun = .Eno, + na.action = na.action, + ncores = ncores)$output1 + + return(effnumobs) +} + +.Eno <- function(x, na.action) { + # Checks: + if (!is.numeric(x)) { + stop("Parameter 'x' must be a numeric vector.") + } + + if (length(sort(x)) > 1) { + n <- length(sort(x)) +# if (n == 0) { +# n <- 1 +# } + a <- acf(x, lag.max = n - 1, plot = FALSE, + na.action = na.action)$acf[2:n, 1, 1] + s <- 0 + for (k in 1:(n - 1)) { + s <- s + (((n - k) / n) * a[k]) + } + + eno <- min(n / (1 + (2 * s)), n) + } else { - dimsvar <- dimsvar[-posdim] + eno <- NA } - effnumobs <- array(dim = dimsvar) - effnumobs[] <- aperm(enleno, posaperm) - # - # Outputs - # ~~~~~~~~~ - # - effnumobs + eno } + diff --git a/R/Trend.R b/R/Trend.R index ba35b4d7d7f4b2038b0fc9df783aadcd8a28f327..6bcfe5467f1e23c207b4d8072f9ecff62ecce076 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -1,167 +1,215 @@ -#'Computes the Trend of the Ensemble Mean +#'Compute the trend #' -#'Computes the trend along the forecast time of the ensemble mean by least -#'square fitting, and the associated error interval.\cr -#'Trend() also provides the time series of the detrended ensemble mean -#'forecasts.\cr -#'The confidence interval relies on a student-T distribution.\cr\cr -#'.Trend provides the same functionality but taking a matrix ensemble members -#'as input (exp). +#'Compute the linear trend or any degree of polynomial regression along the +#'forecast time. It returns the regression coefficients (including the intercept) +#'and the confidence intervals if needed. The detrended array is also provided.\cr +#'The confidence interval relies on the student-T distribution.\cr\cr #' -#'@param var An array of any number of dimensions up to 10. -#'@param posTR An integer indicating the position along which to compute the -#' trend. -#'@param interval A number of months/years between 2 points along posTR -#' dimension. Set 1 as default. -#'@param siglev A numeric value indicating the confidence level for the -#' computation of confidence interval. Set 0.95 as default. -#'@param conf A logical value indicating whether to compute the confidence -#' levels or not. Set TRUE as default. -#'@param exp An M by N matrix representing M forecasts from N ensemble members. +#'@param data An numeric array including the dimension along which the trend +#' is computed. +#'@param time_dim A character string indicating the dimension along which to +#' compute the trend. The default value is 'sdate'. +#'@param interval A positive numeric indicating the unit length between two +#' points along 'time_dim' dimension. The default value is 1. +#'@param polydeg A positive integer indicating the degree of polynomial +#' regression. The default value is 1. +#'@param conf A logical value indicating whether to retrieve the confidence +#' intervals or not. The default value is TRUE. +#'@param conf.lev A numeric indicating the confidence level for the +#' regression computation. The default value is 0.95. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. The default value is NULL. #' #'@return #'\item{$trend}{ -#' The intercept and slope coefficients for the least squares fitting of the -#' trend. -#' An array with same dimensions as parameter 'var' except along the posTR -#' dimension, which is replaced by a length 4 (or length 2 if conf = FALSE) -#' dimension, corresponding to the lower limit of the confidence interval -#' (only present if conf = TRUE), the slope, the upper limit of the confidence -#' interval (only present if conf = TRUE), and the intercept. +#' A numeric array with the first dimension 'stats', followed by the same +#' dimensions as parameter 'data' except the 'time_dim' dimension. The length +#' of the 'stats' dimension should be \code{polydeg + 1}, containing the +#' regression coefficients from the lowest order (i.e., intercept) to the +#' highest degree. +#'} +#'\item{$conf.lower}{ +#' A numeric array with the first dimension 'stats', followed by the same +#' dimensions as parameter 'data' except the 'time_dim' dimension. The length +#' of the 'stats' dimension should be \code{polydeg + 1}, containing the +#' lower limit of the \code{conf.lev}\% confidence interval for all the +#' regression coefficients with the same order as \code{$trend}. Only present +#' \code{conf = TRUE}. +#'} +#'\item{$conf.upper}{ +#' A numeric array with the first dimension 'stats', followed by the same +#' dimensions as parameter 'data' except the 'time_dim' dimension. The length +#' of the 'stats' dimension should be \code{polydeg + 1}, containing the +#' upper limit of the \code{conf.lev}\% confidence interval for all the +#' regression coefficients with the same order as \code{$trend}. Only present +#' \code{conf = TRUE}. #'} #'\item{$detrended}{ -#' Same dimensions as var with linearly detrended var along the posTR -#' dimension. +#' A numeric array with the same dimensions as paramter 'data', containing the +#' detrended values along the 'time_dim' dimension. #'} -#'Only in .Trend: -#'\item{$conf.int}{ -#' Corresponding to the limits of the \code{siglev}\% confidence interval -#' (only present if \code{conf = TRUE}) for the slope coefficient. -#'} #' #'@keywords datagen #'@author History:\cr #'0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@@ic3.cat}) - Original code\cr #'1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@@ic3.cat}) - Formatting to CRAN\cr #'2.0 - 2017-02 (A. Hunter, \email{alasdair.hunter@@bsc.es}) - Adapt to veriApply() +#'3.0 - 2019-12 (A. Ho, \email{an.ho@bsc.es}) - Adapt multiApply feature #'@examples #'# Load sample data as in Load() example: #'example(Load) #'months_between_startdates <- 60 -#'trend <- Trend(sampleData$obs, 3, months_between_startdates) -#' \donttest{ -#'PlotVsLTime(trend$trend, toptitle = "trend", ytitle = "K / (5 year)", -#' monini = 11, limits = c(-1,1), listexp = c('CMIP5 IC3'), -#' listobs = c('ERSST'), biglab = FALSE, hlines = 0, -#' fileout = 'tos_obs_trend.eps') -#'PlotAno(trend$detrended, NULL, startDates, -#' toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', -#' legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') -#' } +#'trend <- Trend(sampleData$obs, polydeg = 2) #' #'@rdname Trend +#'@import multiApply #'@export -Trend <- function(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) { - # - # Enlarge the size of var to 10 and move posTR to first position - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - dimsvar <- dim(var) - if (is.null(dimsvar)) { - dimsvar <- length(var) +Trend <- function(data, time_dim = 'sdate', interval = 1, polydeg = 1, + conf = TRUE, conf.lev = 0.95, ncores = NULL) { + + # Check inputs + ## data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") } - enlvar <- Enlarge(var, 10) - outdim <- c(dimsvar, array(1, dim = (10 - length(dimsvar)))) - if (conf) { - nvals <- 4 - poscoef2 <- 2 - poscoef1 <- 4 - } else { - nvals <- 2 - poscoef2 <- 1 - poscoef1 <- 2 + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") } - outdim[posTR] <- nvals - posaperm <- 1:10 - posaperm[posTR] <- 1 - posaperm[1] <- posTR - enlvar <- aperm(enlvar, posaperm) - dimsaperm <- outdim[posaperm] - # - # Loop on all dimensions to compute trends - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - enltrend <- array(dim = dimsaperm) - enldetrend <- array(dim = dim(enlvar)) - for (j2 in 1:dimsaperm[2]) { - for (j3 in 1:dimsaperm[3]) { - for (j4 in 1:dimsaperm[4]) { - for (j5 in 1:dimsaperm[5]) { - for (j6 in 1:dimsaperm[6]) { - for (j7 in 1:dimsaperm[7]) { - for (j8 in 1:dimsaperm[8]) { - for (j9 in 1:dimsaperm[9]) { - for (j10 in 1:dimsaperm[10]) { - tmp <- enlvar[, j2, j3, j4, j5, j6, j7, j8, j9, j10] - if (any(!is.na(tmp))) { - mon <- seq(tmp) * interval - lm.out <- lm(tmp ~ mon, na.action = na.omit) - enltrend[poscoef2, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- lm.out$coefficients[2] - enltrend[poscoef1, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- lm.out$coefficients[1] - if (conf) { - enltrend[c(1, 3), j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- confint(lm.out, level = siglev)[2, 1:2] - } - enldetrend[is.na(tmp) == FALSE, j2, j3, j4, j5, j6, j7, j8, - j9, j10] <- tmp[is.na(tmp) == FALSE] - lm.out$fitted.values - } - } - } - } - } - } - } - } + if (is.null(dim(data))) { #is vector + dim(data) <- c(length(data)) + names(dim(data)) <- time_dim + } + if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) { + stop("Parameter 'data' must have dimension names.") + } + ## 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(data))) { + stop("Parameter 'time_dim' is not found in 'data' dimension.") + } + ## interval + if (any(!is.numeric(interval) | interval <= 0 | length(interval) > 1)) { + stop("Parameter 'interval' must be a positive number.") + } + ## polydeg + if (!is.numeric(polydeg) | polydeg %% 1 != 0 | polydeg < 0 | + length(polydeg) > 1) { + stop("Parameter 'polydeg' must be a positive integer.") + } + ## conf + if (!is.logical(conf) | length(conf) > 1) { + stop("Parameter 'conf' must be one logical value.") + } + ## conf.lev + if (!is.numeric(conf.lev) | conf.lev < 0 | conf.lev > 1 | length(conf.lev) > 1) { + stop("Parameter 'conf.lev' must be a numeric number between 0 and 1.") + } + ## 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.") } } - # - # Back to the original dimensions - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - enldetrend <- aperm(enldetrend, posaperm) - dim(enldetrend) <- dimsvar - enltrend <- aperm(enltrend, posaperm) - dimsvar[posTR] <- nvals - dim(enltrend) <- dimsvar - # - # Outputs - # ~~~~~~~~~ - # - invisible(list(trend = enltrend, detrended = enldetrend)) + + ############################### + # Calculate Trend + dim_names <- names(dim(data)) + + if (conf) { + output_dims <- list(trend = 'stats', conf.lower = 'stats', + conf.upper = 'stats', detrended = time_dim) + } else if (!conf) { + output_dims <- list(trend = 'stats', detrended = time_dim) + } + + + output <- Apply(list(data), + target_dims = time_dim, + fun = .Trend, + output_dims = output_dims, + time_dim = time_dim, interval = interval, + polydeg = polydeg, conf = conf, + conf.lev = conf.lev, + ncores = ncores) + + #output <- lapply(output, .reorder, time_dim = time_dim, dim_names = dim_names) + + return(output) } -#'@rdname Trend -#'@importFrom stats lm na.omit confint -#'@export -.Trend <- function(exp, interval = 1, siglev = 0.95, conf = TRUE) { - - ensmean <- rowMeans(exp, na.rm = TRUE) - - if (any(!is.na(ensmean))) { - mon <- seq(ensmean) * interval - lm.out <- lm(ensmean ~ mon, na.action = na.omit) - trend <- c(lm.out$coefficients[2], lm.out$coefficients[1]) +.Trend <- function(x, time_dim = 'sdate', interval = 1, polydeg = 1, + conf = TRUE, conf.lev = 0.95) { + + mon <- seq(x) * interval + + # remove NAs for potential poly() + NApos <- 1:length(x) + NApos[which(is.na(x))] <- NA + x2 <- x[!is.na(NApos)] + mon2 <- mon[!is.na(NApos)] + + if (length(x2) > 0) { +# lm.out <- lm(x ~ mon, na.action = na.omit) + lm.out <- lm(x2 ~ poly(mon2, degree = polydeg, raw = TRUE), na.action = na.omit) + trend <- lm.out$coefficients #intercept, slope1, slope2,... + + if (conf) { + conf.lower <- confint(lm.out, level = conf.lev)[, 1] + conf.upper <- confint(lm.out, level = conf.lev)[, 2] + } + + detrended <- c() + detrended[is.na(x) == FALSE] <- x[is.na(x) == FALSE] - lm.out$fitted.values + } else { + trend <- rep(NA, polydeg + 1) + detrend <- NA if (conf) { - conf.int <- confint(lm.out, level = siglev)[2, 1:2] + conf.lower <- rep(NA, polydeg + 1) + conf.upper <- rep(NA, polydeg + 1) } - detrend <- ensmean[is.na(ensmean) == FALSE] - lm.out$fitted.values - } - - # - # Outputs - # ~~~~~~~~~ - # - invisible(list(trend = trend, conf.int = conf.int, detrended = detrend)) + } + + if (conf) { + return(list(trend = trend, conf.lower = conf.lower, conf.upper = conf.upper, + detrended = detrended)) + } else { + return(list(trend = trend, detrended = detrended)) + } + } + + .reorder <- function(output, time_dim, dim_names) { + + # Add dim name back + if (is.null(dim(output))) { + dim(output) <- c(stats = length(output)) + + } else { #is an array + if (length(dim(output)) == 1) { + if (!is.null(names(dim(output)))) { + dim(output) <- c(1, dim(output)) + names(dim(output))[1] <- time_dim + } else { + names(dim(output)) <- time_dim + } + } else { # more than one dim + if (names(dim(output))[1] != "") { + dim(output) <- c(1, dim(output)) + names(dim(output))[1] <- time_dim + } else { #regular case + names(dim(output))[1] <- time_dim + } + } + } + # reorder + pos <- match(dim_names, names(dim(output))) + output <- aperm(output, pos) + names(dim(output)) <- dim_names + names(dim(output))[names(dim(output)) == time_dim] <- 'stats' + + return(output) + } + diff --git a/man/Eno.Rd b/man/Eno.Rd index ba4f2088a3b8de234afe763c71bc401dfd70f4b1..b2a3123d81cd8727f783698b6c6b0f2d5e9dac66 100644 --- a/man/Eno.Rd +++ b/man/Eno.Rd @@ -4,15 +4,21 @@ \alias{Eno} \title{Computes Effective Sample Size With Classical Method} \usage{ -Eno(obs, posdim) +Eno(obs, posdim, na.action = na.pass, ncores = NULL) } \arguments{ -\item{obs}{Matrix of any number of dimensions up to 10.} +\item{obs}{A numeric n-dimensional array with or without named dimensions.} -\item{posdim}{Dimension along which to compute the effective sample size.} +\item{posdim}{A character or an integer indicating the dimension along +which to compute the effective sample size.} + +\item{na.action}{Function to be called to handle missing values, can be +na.pass or na.fail. Default: na.pass} + +\item{ncores}{The number of cores to use for parallel computation. Default: NULL} } \value{ -Same dimensions as var but without the posdim dimension. +An array with the same dimensions as obs except the posdim dimension. } \description{ Computes the effective number of independent values along the posdim @@ -58,6 +64,7 @@ PlotEquiMap(eno, sampleData$lon, sampleData$lat) History:\cr 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Original code\cr 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN +3.0 - 2019-08 (A. Ho, \email{an.ho@bsc.es}) - Adopt multiApply } \keyword{datagen} diff --git a/man/Trend.Rd b/man/Trend.Rd index 3b7f7bfdd541cac1246494570807500ac8718caa..f058009aa6b7fe56ce5591e7edc211a4d12e2477 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -1,74 +1,74 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/Trend.R \name{Trend} -\alias{.Trend} \alias{Trend} -\title{Computes the Trend of the Ensemble Mean} +\title{Compute the trend} \usage{ -Trend(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) - -.Trend(exp, interval = 1, siglev = 0.95, conf = TRUE) +Trend(data, time_dim = "sdate", interval = 1, polydeg = 1, conf = TRUE, + conf.lev = 0.95, ncores = NULL) } \arguments{ -\item{var}{An array of any number of dimensions up to 10.} +\item{data}{An numeric array including the dimension along which the trend +is computed.} + +\item{time_dim}{A character string indicating the dimension along which to +compute the trend. The default value is 'sdate'.} -\item{posTR}{An integer indicating the position along which to compute the -trend.} +\item{interval}{A positive numeric indicating the unit length between two +points along 'time_dim' dimension. The default value is 1.} -\item{interval}{A number of months/years between 2 points along posTR -dimension. Set 1 as default.} +\item{polydeg}{A positive integer indicating the degree of polynomial +regression. The default value is 1.} -\item{siglev}{A numeric value indicating the confidence level for the -computation of confidence interval. Set 0.95 as default.} +\item{conf}{A logical value indicating whether to retrieve the confidence +intervals or not. The default value is TRUE.} -\item{conf}{A logical value indicating whether to compute the confidence -levels or not. Set TRUE as default.} +\item{conf.lev}{A numeric indicating the confidence level for the +regression computation. The default value is 0.95.} -\item{exp}{An M by N matrix representing M forecasts from N ensemble members.} +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} } \value{ \item{$trend}{ - The intercept and slope coefficients for the least squares fitting of the - trend. - An array with same dimensions as parameter 'var' except along the posTR - dimension, which is replaced by a length 4 (or length 2 if conf = FALSE) - dimension, corresponding to the lower limit of the confidence interval - (only present if conf = TRUE), the slope, the upper limit of the confidence - interval (only present if conf = TRUE), and the intercept. + A numeric array with the first dimension 'stats', followed by the same + dimensions as parameter 'data' except the 'time_dim' dimension. The length + of the 'stats' dimension should be \code{polydeg + 1}, containing the + regression coefficients from the lowest order (i.e., intercept) to the + highest degree. } -\item{$detrended}{ - Same dimensions as var with linearly detrended var along the posTR - dimension. +\item{$conf.lower}{ + A numeric array with the first dimension 'stats', followed by the same + dimensions as parameter 'data' except the 'time_dim' dimension. The length + of the 'stats' dimension should be \code{polydeg + 1}, containing the + lower limit of the \code{conf.lev}\% confidence interval for all the + regression coefficients with the same order as \code{$trend}. Only present + \code{conf = TRUE}. } -Only in .Trend: -\item{$conf.int}{ - Corresponding to the limits of the \code{siglev}\% confidence interval - (only present if \code{conf = TRUE}) for the slope coefficient. +\item{$conf.upper}{ + A numeric array with the first dimension 'stats', followed by the same + dimensions as parameter 'data' except the 'time_dim' dimension. The length + of the 'stats' dimension should be \code{polydeg + 1}, containing the + upper limit of the \code{conf.lev}\% confidence interval for all the + regression coefficients with the same order as \code{$trend}. Only present + \code{conf = TRUE}. +} +\item{$detrended}{ + A numeric array with the same dimensions as paramter 'data', containing the + detrended values along the 'time_dim' dimension. } } \description{ -Computes the trend along the forecast time of the ensemble mean by least -square fitting, and the associated error interval.\cr -Trend() also provides the time series of the detrended ensemble mean -forecasts.\cr -The confidence interval relies on a student-T distribution.\cr\cr -.Trend provides the same functionality but taking a matrix ensemble members -as input (exp). +Compute the linear trend or any degree of polynomial regression along the +forecast time. It returns the regression coefficients (including the intercept) +and the confidence intervals if needed. The detrended array is also provided.\cr +The confidence interval relies on the student-T distribution.\cr\cr } \examples{ # Load sample data as in Load() example: example(Load) months_between_startdates <- 60 -trend <- Trend(sampleData$obs, 3, months_between_startdates) - \donttest{ -PlotVsLTime(trend$trend, toptitle = "trend", ytitle = "K / (5 year)", - monini = 11, limits = c(-1,1), listexp = c('CMIP5 IC3'), - listobs = c('ERSST'), biglab = FALSE, hlines = 0, - fileout = 'tos_obs_trend.eps') -PlotAno(trend$detrended, NULL, startDates, - toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', - legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') - } +trend <- Trend(sampleData$obs, polydeg = 2) } \author{ @@ -76,6 +76,7 @@ History:\cr 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@ic3.cat}) - Original code\cr 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@ic3.cat}) - Formatting to CRAN\cr 2.0 - 2017-02 (A. Hunter, \email{alasdair.hunter@bsc.es}) - Adapt to veriApply() +3.0 - 2019-12 (A. Ho, \email{an.ho@bsc.es}) - Adapt multiApply feature } \keyword{datagen} diff --git a/tests/testthat/test-Eno.R b/tests/testthat/test-Eno.R new file mode 100644 index 0000000000000000000000000000000000000000..e8c7e6cd30605b4c602e9f3ba5ed67ce09b6fe6f --- /dev/null +++ b/tests/testthat/test-Eno.R @@ -0,0 +1,102 @@ +context("Generic tests") +test_that("Sanity checks", { + +#error + expect_error( + Eno(NULL), + "Parameter 'obs' cannot be NULL." + ) + + expect_error( + Eno(c('a', 'b'), posdim = 1), + "Parameter 'obs' must be a numeric array." + ) + + expect_error( + Eno(1:10, c()), + "Parameter 'posdim' cannot be NULL." + ) + + expect_error( + Eno(1:10, posdim = FALSE), + "Parameter 'posdim' must be an positive integer or character." + ) + + expect_error( + Eno(1:10, posdim = -1), + "Parameter 'posdim' must be a positive integer or character." + ) + + expect_error( + Eno(1:10, posdim = 1.2), + "Parameter 'posdim' must be a positive integer or character." + ) + + expect_error( + Eno(array(1:10, dim = c(2, 5)), 3), + "Parameter 'posdim' excesses the dimension length of parameter 'obs'." + ) + + dat <- array(1:10, dim = c(2, 5)) + expect_error( + Eno(dat, 'a'), + paste0("Parameter 'obs' doesn't have dimension names. ", + "Add names or specify parameter 'posdim' by number.") + ) + + names(dim(dat)) <- c('a', 'b') + expect_error( + Eno(dat, 'c'), + paste0("Parameter 'posdim' does not match any dimension name of ", + "parameter 'obs'.") + ) + + expect_error( + Eno(dat, 'a', na.action = na.stop), + "Parameter 'na.action' must be either 'na.pass' or 'na.fail'." + ) + + dat[1] <- NA + expect_error( + Eno(dat, 'a', na.action = na.fail), + paste0("NA value in paratemter 'obs', which is not accepted when ", + "parameter 'na.action' = 'na.fail'.") + ) + + expect_error( + Eno(dat, 'a', ncores = FALSE), + "Parameter 'ncores' must be a positive integer." + ) + + expect_error( + Eno(dat, 'a', ncores = 1.5), + "Parameter 'ncores' must be a positive integer." + ) + +#equal + dat <- array(1:10, dim = c(2, 5)) + + expect_equal( + names(dim(Eno(dat, 1))), + c('D2') + ) + + names(dim(dat)) <- c('a', 'b') + expect_equal( + Eno(dat, 1), + array(rep(2, 5), dim = c(b = 5)) + ) + + dat[1] <- NA + expect_equal( + Eno(dat, 2)[1], + 4 + ) + + expect_equal( + Eno(dat, 1)[1], + as.numeric(NA) + ) + +}) + diff --git a/tests/testthat/test-Trend.R b/tests/testthat/test-Trend.R new file mode 100644 index 0000000000000000000000000000000000000000..8843b9ed75bcc651826880ced57d8a4d1010e193 --- /dev/null +++ b/tests/testthat/test-Trend.R @@ -0,0 +1,143 @@ +context("s2dverification::Trend tests") + +############################################## + # dat1 + dat1 <- array(c(-5, -7, -10:10, 12, 11, 7, 16), + dim = c(dat = 1, sdate = 13, ftime = 2)) + + # dat2 + set.seed(10) + dat2 <- c(1:10) + rnorm(10) + + # dat3 + set.seed(1) + dat3 <- array(c(1:60) + rnorm(60), + dim = c(sdate = 5, lat = 2, lon = 3, lev = 2)) + set.seed(1) + na <- floor(runif(5, min = 1, max = 60)) + dat3[na] <- NA + +############################################## + +test_that("1. Input checks", { + + expect_error( + Trend(c()), + "Parameter 'data' cannot be NULL." + ) + expect_error( + Trend(c(NA, NA)), + "Parameter 'data' must be a numeric array." + ) + expect_error( + Trend(list(a = array(rnorm(50), dim = c(dat = 5, sdate = 10)), b = c(1:4))), + "Parameter 'data' must be a numeric array." + ) + expect_error( + Trend(array(1:10, dim = c(2, 5))), + "Parameter 'data' must have dimension names." + ) + expect_error( + Trend(dat1, time_dim = 'a'), + "Parameter 'time_dim' is not found in 'data' dimension." + ) + expect_error( + Trend(array(c(1:25), dim = c(dat = 1, date = 5, ftime = 5))), + "Parameter 'time_dim' is not found in 'data' dimension." + ) + expect_error( + Trend(dat1, time_dim = 2), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + Trend(dat1, time_dim = c('a','sdate')), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + Trend(dat1, interval = 0), + "Parameter 'interval' must be a positive number." + ) + expect_error( + Trend(dat1, conf = 3), + "Parameter 'conf' must be one logical value." + ) + expect_error( + Trend(dat1, polydeg = 3.5), + "Parameter 'polydeg' must be a positive integer." + ) + expect_error( + Trend(dat1, ncore = 3.5), + "Parameter 'ncores' must be a positive integer." + ) +}) + +############################################## +test_that("2. Output checks: dat1", { + + expect_equal( + Trend(dat1)$trend, + array(c(-9.7692308, 0.6593407, 0.9615385, 0.7967033), + dim = c(stats = 2, dat = 1, ftime = 2)), + tolerance = 0.0001 + ) + expect_equal( + Trend(dat1)$conf.upper, + array(c(-7.4735367, 0.9485709, 3.0167860, 1.0556402), + dim = c(stats = 2, dat = 1, ftime = 2)), + tolerance = 0.0001 + ) + expect_equal( + summary(Trend(dat1)$detrended)[3], + c(Median = 0.1154), + tolerance = 0.001 + ) + +}) + +############################################## +test_that("3. Output checks: dat2", { + + expect_equal( + Trend(dat2), + list(trend = array(c(-0.182, 0.944), dim = c(stats = 2)), + conf.lower = array(c(-1.316, 0.761), dim = c(stats = 2)), + conf.upper = array(c(0.953, 1.127), dim = c(stats = 2)), + detrended = array(c(0.257, 0.110, -1.021, -0.193, 0.757, 0.909, + -0.633, 0.267, -0.939, 0.487), dim = c(sdate = 10))), + tolerance = 0.001 + ) + expect_equal( + Trend(dat2, interval = 2), + list(trend = array(c(-0.182, 0.472), dim = c(stats = 2)), + conf.lower = array(c(-1.316, 0.381), dim = c(stats = 2)), + conf.upper = array(c(0.953, 0.563), dim = c(stats = 2)), + detrended = array(c(0.257, 0.110, -1.021, -0.193, 0.757, 0.909, + -0.633, 0.267, -0.939, 0.487), dim = c(sdate = 10))), + tolerance = 0.001 + ) + expect_equal( + length(Trend(dat2, conf = F)), + 2 + ) + expect_equal( + names(Trend(dat2, conf = F)), + c('trend', 'detrended') + ) + +}) + +############################################## +test_that("4. Output checks: dat3", { + + expect_equal( + summary(Trend(dat3)$trend)[3], + c(Median = 1.3680), + tolerance = 0.0001 + ) + expect_equal( + dim(Trend(dat3, polydeg = 2)$trend), + c(stats = 3, lat = 2, lon = 3, lev = 2) + ) + +}) +