From dd6e568014c3f96b7f498e6d240098e248a3666b Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 3 Jun 2020 13:55:06 +0200 Subject: [PATCH 1/2] Add p-value in Trend() --- R/Trend.R | 53 +++++++++++++++++++++++++++++-------- man/Trend.Rd | 16 ++++++++--- tests/testthat/test-Trend.R | 28 ++++++++++++++++++-- 3 files changed, 80 insertions(+), 17 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index c16b3c3..0f20efc 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 e9c890e..17076c6 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 ac9f14a..89052a0 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') ) -- GitLab From 17ab56068c8d37d81c93bea1bdc202fca0ca9a99 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 3 Jun 2020 13:57:46 +0200 Subject: [PATCH 2/2] Update NEWS.md --- NEWS.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NEWS.md b/NEWS.md index c3bf8ee..9f20054 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. -- GitLab