diff --git a/NEWS.md b/NEWS.md index c3bf8eeefa6bb8ac85cfae641d0f7d23d6403fb3..9f20054b9199276760eacdf58106a2ce0c520e12 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,8 @@ # s2dv 0.0.2 (Release date: 2020-) - Change the default of Season() parameter 'time_dim' from 'sdate' to 'ftime'. +- Add p-value by ANOVA in Trend(). +- Change MeanDims() na.rm default to FALSE to be in line with mean() +- Remove unecessary parameter checks in Clim() # s2dv 0.0.1 (Release date: 2020-02-07) - The package is the advanced version of package 's2dverification', adopting the regime of package 'multiApply' for all the analytic functions. Most of the other functions for plotting and data retrieval in 's2dverification' are also preserved in this package. diff --git a/R/Trend.R b/R/Trend.R index c16b3c31c37cbae925563439ef33aefc4db60cf4..0f20efc0cb07ffd5dec036b77835b19c2c81a7d3 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -2,8 +2,10 @@ #' #'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 +#'and the detrended array. The confidence intervals and p-value are also +#'provided if needed.\cr +#'The confidence interval relies on the student-T distribution, and the p-value +#'is calculated by ANOVA. #' #'@param data An numeric array including the dimension along which the trend #' is computed. @@ -17,6 +19,8 @@ #' 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 pval A logical value indicating whether to compute the p-value or not. +#' The default value is TRUE. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' @@ -45,6 +49,9 @@ #' regression coefficients with the same order as \code{$trend}. Only present #' \code{conf = TRUE}. #'} +#'\item{$p.val}{ +#' The p-value calculated by anova(). Only present if \code{pval = TRUE}. +#'} #'\item{$detrended}{ #' A numeric array with the same dimensions as paramter 'data', containing the #' detrended values along the 'time_dim' dimension. @@ -54,13 +61,13 @@ #'# Load sample data as in Load() example: #'example(Load) #'months_between_startdates <- 60 -#'trend <- Trend(sampleData$obs, polydeg = 2) +#'trend <- Trend(sampleData$obs, polydeg = 2, interval = months_between_startdates) #' #'@rdname Trend #'@import multiApply #'@export Trend <- function(data, time_dim = 'sdate', interval = 1, polydeg = 1, - conf = TRUE, conf.lev = 0.95, ncores = NULL) { + conf = TRUE, conf.lev = 0.95, pval = TRUE, ncores = NULL) { # Check inputs ## data @@ -101,6 +108,10 @@ Trend <- function(data, time_dim = 'sdate', interval = 1, polydeg = 1, 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.") } + ## pval + if (!is.logical(pval) | length(pval) > 1) { + stop("Parameter 'pval' must be one logical value.") + } ## ncores if (!is.null(ncores)) { if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores < 0 | @@ -113,10 +124,15 @@ Trend <- function(data, time_dim = 'sdate', interval = 1, polydeg = 1, # Calculate Trend dim_names <- names(dim(data)) - if (conf) { + if (conf & pval) { + output_dims <- list(trend = 'stats', conf.lower = 'stats', + conf.upper = 'stats', p.val = 'stats', detrended = time_dim) + } else if (conf & !pval) { output_dims <- list(trend = 'stats', conf.lower = 'stats', conf.upper = 'stats', detrended = time_dim) - } else if (!conf) { + } else if (!conf & pval) { + output_dims <- list(trend = 'stats', p.val = 'stats', detrended = time_dim) + } else { output_dims <- list(trend = 'stats', detrended = time_dim) } @@ -127,16 +143,14 @@ Trend <- function(data, time_dim = 'sdate', interval = 1, polydeg = 1, output_dims = output_dims, time_dim = time_dim, interval = interval, polydeg = polydeg, conf = conf, - conf.lev = conf.lev, + conf.lev = conf.lev, pval = pval, ncores = ncores) - #output <- lapply(output, .reorder, time_dim = time_dim, dim_names = dim_names) - return(output) } .Trend <- function(x, time_dim = 'sdate', interval = 1, polydeg = 1, - conf = TRUE, conf.lev = 0.95) { + conf = TRUE, conf.lev = 0.95, pval = TRUE) { mon <- seq(x) * interval @@ -156,20 +170,37 @@ Trend <- function(data, time_dim = 'sdate', interval = 1, polydeg = 1, conf.upper <- confint(lm.out, level = conf.lev)[, 2] } + if (pval) { + p.val <- as.array(anova(lm.out)$'Pr(>F)'[1]) + } + 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.lower <- rep(NA, polydeg + 1) conf.upper <- rep(NA, polydeg + 1) } + + if (pval) { + p.val <- NA + } + } - if (conf) { + if (conf & pval) { return(list(trend = trend, conf.lower = conf.lower, conf.upper = conf.upper, + p.val = p.val, detrended = detrended)) + } else if (conf & !pval) { + return(list(trend = trend, conf.lower = conf.lower, conf.upper = conf.upper, detrended = detrended)) + } else if (!conf & pval) { + return(list(trend = trend, p.val = p.val, detrended = detrended)) } else { return(list(trend = trend, detrended = detrended)) } diff --git a/man/Trend.Rd b/man/Trend.Rd index e9c890e4d6e7baa2a17a529b053baa952a39547b..17076c6e690e8063bb85b0bbf17550b24c9342f4 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -5,7 +5,7 @@ \title{Compute the trend} \usage{ Trend(data, time_dim = "sdate", interval = 1, polydeg = 1, conf = TRUE, - conf.lev = 0.95, ncores = NULL) + conf.lev = 0.95, pval = TRUE, ncores = NULL) } \arguments{ \item{data}{An numeric array including the dimension along which the trend @@ -26,6 +26,9 @@ intervals or not. The default value is TRUE.} \item{conf.lev}{A numeric indicating the confidence level for the regression computation. The default value is 0.95.} +\item{pval}{A logical value indicating whether to compute the p-value or not. +The default value is TRUE.} + \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } @@ -54,6 +57,9 @@ A list containing: regression coefficients with the same order as \code{$trend}. Only present \code{conf = TRUE}. } +\item{$p.val}{ + The p-value calculated by anova(). Only present if \code{pval = TRUE}. +} \item{$detrended}{ A numeric array with the same dimensions as paramter 'data', containing the detrended values along the 'time_dim' dimension. @@ -62,14 +68,16 @@ A list containing: \description{ 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 +and the detrended array. The confidence intervals and p-value are also +provided if needed.\cr +The confidence interval relies on the student-T distribution, and the p-value +is calculated by ANOVA. } \examples{ # Load sample data as in Load() example: example(Load) months_between_startdates <- 60 -trend <- Trend(sampleData$obs, polydeg = 2) +trend <- Trend(sampleData$obs, polydeg = 2, interval = months_between_startdates) } diff --git a/tests/testthat/test-Trend.R b/tests/testthat/test-Trend.R index ac9f14a1811ec19feab5913b5d16077610f958f6..89052a0b119d09c02682f8a0b822caa3f3e722f8 100644 --- a/tests/testthat/test-Trend.R +++ b/tests/testthat/test-Trend.R @@ -61,6 +61,18 @@ test_that("1. Input checks", { Trend(dat1, conf = 3), "Parameter 'conf' must be one logical value." ) + expect_error( + Trend(dat1, conf.lev = 3), + "Parameter 'conf.lev' must be a numeric number between 0 and 1." + ) + expect_error( + Trend(dat1, conf.lev = TRUE), + "Parameter 'conf.lev' must be a numeric number between 0 and 1." + ) + expect_error( + Trend(dat1, pval = 0.95), + "Parameter 'pval' must be one logical value." + ) expect_error( Trend(dat1, polydeg = 3.5), "Parameter 'polydeg' must be a positive integer." @@ -86,6 +98,12 @@ test_that("2. Output checks: dat1", { dim = c(stats = 2, dat = 1, ftime = 2)), tolerance = 0.0001 ) + expect_equal( + Trend(dat1)$p.val, + array(c(0.0003916111, 3.065828e-05), + dim = c(stats = 1, dat = 1, ftime = 2)), + tolerance = 0.0001 + ) expect_equal( median(Trend(dat1)$detrended, na.rm = TRUE), 0.1153846, @@ -102,6 +120,7 @@ test_that("3. Output checks: 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)), + p.val = array(2.277774e-06, dim = c(stats = 1)), 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 @@ -111,16 +130,21 @@ test_that("3. Output checks: dat2", { 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)), + p.val = array(2.277774e-06, dim = c(stats = 1)), 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 + 3 + ) + expect_equal( + length(Trend(dat2, pval = F)), + 4 ) expect_equal( - names(Trend(dat2, conf = F)), + names(Trend(dat2, conf = F, pval = F)), c('trend', 'detrended') )