From a94ea1678aabf7c17d67d0fea42b73b648bae4e2 Mon Sep 17 00:00:00 2001 From: "Neven S. Fuckar" Date: Thu, 19 Mar 2015 19:00:52 +0100 Subject: [PATCH 01/18] Upgrade of Trend.R to include polynomial fit of an arbitrary degree/order. Fixes #36 --- R/Trend.R | 75 +++++++++++++++++++++++++++++++--------------------- man/Trend.Rd | 44 ++++++++++++++++++------------ 2 files changed, 72 insertions(+), 47 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index ade5e76b..a3d2d795 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -1,43 +1,55 @@ -Trend <- function(var, posTR = 2, interval = 1) { - # Compute trends along the (posTR)th dimension of matrix var by least square - # fitting, and the associated an error interval. - # Provide also the detrended data. +Trend <- function(var, posTR = 2, interval = 1, pldegree = 1) { + # Computes trend or higher order polynomial along the posTR dimension of tensor var + # with the method of least squares, and estimates the associated error interval. + # Provides also the detrended or residual (after removing fitted polynomial) data. # The confidence interval relies on a student-T distribution. # # Args: - # var: Matrix of any number of dimensions up to 10. - # posTR: Position along which to compute trends. + # var: Tensor of any number of dimensions up to 10. + # posTR: Position along which to compute trends or polynomial fits. # interval: Number of months between 2 points along posTR dimension. - # Default = 1. + # Default = 1 (assumes monthly data). + # pldegree: The degree/order of the polynomial to be fitted and removed + # (keep in mind that the degree must be smaller than the number + # of points along posTR dimension). Default = 1 (linear fit). # # Returns: - # $trend: Same dimensions as var except along posTR dimension which is - # replaced by a length 3 dimension, corresponding to the lower - # limit of the 95% confidence interval, the computed trends and - # the upper limit of the 95% confidence interval for each point of - # the matrix along all the other dimensions. - # $detrended: Same dimensions as var with linearly detrended var along the - # posTR dimension. - # + # $trend: Same dimensions as var except along posTR dimension which is replaced + # by a dimension with length 3*(pldegree + 1) showing for + # i) pldegree=1: a^(lo)_1, a_1 (slope or trend), a^(hi)_1, a_0 (intercept), + # a^(lo)_0, and a^(hi)_0 due to backward compatiliblity requirement, + # ii) pldegree>1: a^(lo)_pldegree, a_pldegree, a^(hi)_pldegree, + # a^(lo)_(pldegree-1), a_(pldegree-1), a^(hi)_(pldegree-1), ... + # a^(lo)_0, a_0, and a^(hi)_0, where (lo) is the lower limit of the + # 95% confidence interval and (hi) is the upper limit of the 95% confidence + # interval for each point of the tensor var along all other dimensions. + # $detrended: Same dimensions as var with residual values (after removing polynomial + # fit of pldegree order) along the posTR dimension. + # # History: # 1.0 # 2011-05 (V. Guemas, vguemas@ic3.cat) # Original code + # 1.1 modified to include fit of a polynomial with an arbitrary (but constrained) + # degree/order and a more complete $trend output by N.S. Fuckar, 2015-01 # # Enlarge the size of var to 10 and move posTR to first position # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # dimsvar <- dim(var) - if (is.null(dimsvar)) { - dimsvar <- length(var) + if (is.null(dimsvar)) { + dimsvar <- length(var) } enlvar <- Enlarge(var, 10) outdim <- c(dimsvar, array(1, dim = (10 - length(dimsvar)))) - outdim[posTR] <- 4 + outdim[posTR] <- 3*(pldegree + 1) posaperm <- 1:10 posaperm[posTR] <- 1 posaperm[1] <- posTR enlvar <- aperm(enlvar, posaperm) dimsaperm <- outdim[posaperm] + if ( dim(enlvar)[1]<(pldegree + 1) ) { + stop("selected polynomial degree (pldegree) is too high") + } # # Loop on all dimensions to compute trends # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -56,17 +68,20 @@ Trend <- function(var, posTR = 2, interval = 1) { tmp <- enlvar[, j2, j3, j4, j5, j6, j7, j8, j9, j10] if (length(sort(tmp)) > 0) { mon <- seq(tmp) * interval - lm.out <- lm(tmp ~ mon, na.action = na.omit) - enltrend[1, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- confint(lm.out)[2, 1] - enltrend[2, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- lm.out$coefficients[2] - enltrend[3, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- confint(lm.out)[2, 2] - enltrend[4, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- lm.out$coefficients[1] - enldetrend[is.na(tmp) == FALSE, j2, j3, j4, j5, j6, j7, j8, - j9, j10] <- tmp[is.na(tmp) == FALSE] - lm.out$fitted.values + lm.out <- lm(tmp ~ poly(mon, degree=pldegree, raw=TRUE), na.action = na.omit) + for (jj in (pldegree + 1):1) { + jj1 <- 1 + (pldegree + 1) - jj + if ( (jj==1) & (pldegree==1) ) { + enltrend[3*jj1-2, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- lm.out$coefficients[jj] + enltrend[3*jj1-1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 1] + enltrend[3*jj1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 2] + } else { + enltrend[3*jj1-2, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 1] + enltrend[3*jj1-1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- lm.out$coefficients[jj] + enltrend[3*jj1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 2] + } + } + enldetrend[is.na(tmp) == FALSE, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- tmp[is.na(tmp) == FALSE] - lm.out$fitted.values } } } @@ -82,7 +97,7 @@ Trend <- function(var, posTR = 2, interval = 1) { # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # detrended <- array(dim = dimsvar) - dimsvar[posTR] <- 4 + dimsvar[posTR] <- 3*(pldegree + 1) trend <- array(dim = dimsvar) trend[] <- aperm(enltrend, posaperm) detrended[] <- aperm(enldetrend, posaperm) diff --git a/man/Trend.Rd b/man/Trend.Rd index b4bd0080..b6a88d7b 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -4,37 +4,46 @@ Computes Trends } \description{ - Computes the trend along the posTR dimension of the matrix var by least - square fitting, and the associated an error interval. - Provide also the detrended data. + Computes trend or higher order polynomial along the posTR dimension of tensor var + with the method of least squares, and estimates the associated error interval. + Provides also the detrended or residual (after removing fitted polynomial) data. The confidence interval relies on a student-T distribution. } \usage{ - Trend(var, posTR = 2, interval = 1) + Trend(var, posTR = 2, interval = 1, pldegree=1) } \arguments{ \item{var}{ - Matrix of any number of dimensions up to 10. + Tensor of any number of dimensions up to 10. } \item{posTR}{ - Position along which to compute the trend. + Position along which to compute trends or polynomial fits. } \item{interval}{ - Number of months/years between 2 points along posTR dimension. - Default = 1. - The trend would be provided in number of units per month or year. + Number of months between 2 points along posTR dimension. + Default = 1 (assumes monthly data). + } + \item(pldgree){ + The degree/order of the polynomial to be fitted and removed + (keep in mind that the degree must be smaller than the number + of points along posTR dimension). Default = 1 (linear fit) } } \value{ - \item{$trend}{ - Same dimensions as var except along the posTR dimension which is replaced - by a length 3 dimension, corresponding to the lower limit of the 95\% - confidence interval, trends and the upper limit of the 95\% confidence - interval for each point of the matrix along all the other dimensions. + \item{$trend}{ + Same dimensions as var except along posTR dimension which is replaced + by a dimension with length 3*(pldegree + 1) showing for + i) pldegree=1: a^(lo)_1, a_1 (slope or trend), a^(hi)_1, a_0 (intercept), + a^(lo)_0, and a^(hi)_0 due to backward compatiliblity requirement, + ii) pldegree>1: a^(lo)_pldegree, a_pldegree, a^(hi)_pldegree, + a^(lo)_(pldegree-1), a_(pldegree-1), a^(hi)_(pldegree-1), ... + a^(lo)_0, a_0, and a^(hi)_0, where (lo) is the lower limit of the + 95% confidence interval and (hi) is the upper limit of the 95% confidence + interval for each point of the tensor var along all other dimensions. } \item{$detrended}{ - Same dimensions as var with linearly detrended var along the posTR - dimension. + Same dimensions as var with residual values (after removing polynomial + fit of pldegree order) along the posTR dimension. } } \examples{ @@ -54,7 +63,8 @@ } \author{ History: - 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@ic3.cat}) - Original code + 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@ic3.cat}) - Original code 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@ic3.cat}) - Formatting to CRAN + 1.1 - 2015-01 (N.S. Fuckar, \email{neven.fuckar@ic3.cat}) - Upgrade to polynomial fit } \keyword{datagen} -- GitLab From e08f5d8588fde1d21ec2d83ef2e2f3b533db2999 Mon Sep 17 00:00:00 2001 From: "Neven S. Fuckar" Date: Mon, 23 Mar 2015 14:58:37 +0100 Subject: [PATCH 02/18] more polished help --- R/Trend.R | 16 +++++++++------- man/Trend.Rd | 17 +++++++++++------ 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index a3d2d795..4277125f 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -1,12 +1,12 @@ Trend <- function(var, posTR = 2, interval = 1, pldegree = 1) { - # Computes trend or higher order polynomial along the posTR dimension of tensor var - # with the method of least squares, and estimates the associated error interval. + # Computes linear trend or higher order polynomial along the posTR dimension of tensor + # var with the method of least squares, and estimates the associated error interval. # Provides also the detrended or residual (after removing fitted polynomial) data. # The confidence interval relies on a student-T distribution. # # Args: - # var: Tensor of any number of dimensions up to 10. - # posTR: Position along which to compute trends or polynomial fits. + # var: Matrix of any number of dimensions up to 10. + # posTR: Position along which to compute trends. # interval: Number of months between 2 points along posTR dimension. # Default = 1 (assumes monthly data). # pldegree: The degree/order of the polynomial to be fitted and removed @@ -27,9 +27,11 @@ Trend <- function(var, posTR = 2, interval = 1, pldegree = 1) { # fit of pldegree order) along the posTR dimension. # # History: - # 1.0 # 2011-05 (V. Guemas, vguemas@ic3.cat) # Original code - # 1.1 modified to include fit of a polynomial with an arbitrary (but constrained) - # degree/order and a more complete $trend output by N.S. Fuckar, 2015-01 + # 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@ic3.cat}) - Original code + # 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@ic3.cat}) - Formatting to CRAN + # 1.1 - 2015-01 (N.S. Fuckar, \email{neven.fuckar@ic3.cat}) - Upgrade to polynomial fit + # i.e., modified to include fit of a raw polynomial with an arbitrary + # (but constrained) degree/order and a more complete $trend output. # # Enlarge the size of var to 10 and move posTR to first position diff --git a/man/Trend.Rd b/man/Trend.Rd index b6a88d7b..89d323dd 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -1,11 +1,11 @@ \name{Trend} \alias{Trend} \title{ - Computes Trends + Computes Trends and more } \description{ - Computes trend or higher order polynomial along the posTR dimension of tensor var - with the method of least squares, and estimates the associated error interval. + Computes linear trend or higher order polynomial along the posTR dimension of tensor + var with the method of least squares, and estimates the associated error interval. Provides also the detrended or residual (after removing fitted polynomial) data. The confidence interval relies on a student-T distribution. } @@ -60,11 +60,16 @@ PlotAno(trend$detrended, NULL, startDates, toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') + +... (add one example of nonlinear fit) + } \author{ History: - 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@ic3.cat}) - Original code - 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@ic3.cat}) - Formatting to CRAN - 1.1 - 2015-01 (N.S. Fuckar, \email{neven.fuckar@ic3.cat}) - Upgrade to polynomial fit + 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@ic3.cat}) - Original code + 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@ic3.cat}) - Formatting to CRAN + 1.1 - 2015-01 (N.S. Fuckar, \email{neven.fuckar@ic3.cat}) - Upgrade to polynomial fit + i.e., modified to include fit of a raw polynomial with an arbitrary + (but constrained) degree/order and a more complete $trend output. } \keyword{datagen} -- GitLab From 5b3d710ed5130a73f3cf44edbea865d83c9f8725 Mon Sep 17 00:00:00 2001 From: "Neven S. Fuckar" Date: Mon, 23 Mar 2015 15:19:10 +0100 Subject: [PATCH 03/18] bug fix in Trend.Rd --- man/Trend.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/Trend.Rd b/man/Trend.Rd index 89d323dd..71d3b55e 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -23,7 +23,7 @@ Number of months between 2 points along posTR dimension. Default = 1 (assumes monthly data). } - \item(pldgree){ + \item{pldgree}{ The degree/order of the polynomial to be fitted and removed (keep in mind that the degree must be smaller than the number of points along posTR dimension). Default = 1 (linear fit) -- GitLab From 5a57b613df39cd681955148cf94edccfda9becfa Mon Sep 17 00:00:00 2001 From: "Neven S. Fuckar" Date: Mon, 23 Mar 2015 15:23:28 +0100 Subject: [PATCH 04/18] more help fixes --- man/Trend.Rd | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/man/Trend.Rd b/man/Trend.Rd index 71d3b55e..b4c9046b 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -23,7 +23,7 @@ Number of months between 2 points along posTR dimension. Default = 1 (assumes monthly data). } - \item{pldgree}{ + \item{pldegree}{ The degree/order of the polynomial to be fitted and removed (keep in mind that the degree must be smaller than the number of points along posTR dimension). Default = 1 (linear fit) @@ -61,8 +61,7 @@ toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') -... (add one example of nonlinear fit) - +#... (add one example of nonlinear fit) } \author{ History: -- GitLab From 83fc08bba3d385e33e8a13003ce9505d9781b7e6 Mon Sep 17 00:00:00 2001 From: "Neven S. Fuckar" Date: Mon, 23 Mar 2015 17:59:39 +0100 Subject: [PATCH 05/18] added example of quadratic fit --- man/Trend.Rd | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/man/Trend.Rd b/man/Trend.Rd index b4c9046b..d41e45a0 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -47,21 +47,27 @@ } } \examples{ + # A test of the default linear fit: startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101') sampleData <- Load('tos', c('i00k'), c('ERSST'), startDates, nleadtime = 124, leadtimemin = 1, leadtimemax = 60, output = 'areave', latmin = 30, latmax = 45, lonmin = 0, lonmax = 40) months_between_startdates <- 60 trend <- Trend(sampleData$obs, 3, months_between_startdates) - 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') +##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') -> there is an error in this step? PlotAno(trend$detrended, NULL, startDates, toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', - legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') - -#... (add one example of nonlinear fit) + legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') # these are warnings in this step! + # A simple test of a nonlinear fit: + y0 <- c(10.1, 8.8, 7.3, 4.12, -0.314, -5, -10.6, -18.2) + y1 <- Trend(var=y0, posTR=1) # linear fit + y2 <- Trend(var=y0, posTR=1, pldegree=2) # quadratic fit + plot(y0) + lines(y0 - y1$detrended, col='blue') # linear fit + lines(y0 - y2$detrended, col='red') # quadratic fit } \author{ History: -- GitLab From 4f3038153b06866a43f23fab7a1f8bab87f89a51 Mon Sep 17 00:00:00 2001 From: "Neven S. Fuckar" Date: Mon, 23 Mar 2015 21:09:58 +0100 Subject: [PATCH 06/18] more work on help --- man/Trend.Rd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/man/Trend.Rd b/man/Trend.Rd index d41e45a0..ebab9a6f 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -10,7 +10,7 @@ The confidence interval relies on a student-T distribution. } \usage{ - Trend(var, posTR = 2, interval = 1, pldegree=1) + Trend(var, posTR = 2, interval = 1, pldegree = 1) } \arguments{ \item{var}{ @@ -57,10 +57,10 @@ ##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') -> there is an error in this step? +## fileout = 'tos_obs_trend.eps') -> there is an error in this step, but eps file is produced? PlotAno(trend$detrended, NULL, startDates, toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', - legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') # these are warnings in this step! + legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') # there are warnings in this step! # A simple test of a nonlinear fit: y0 <- c(10.1, 8.8, 7.3, 4.12, -0.314, -5, -10.6, -18.2) y1 <- Trend(var=y0, posTR=1) # linear fit @@ -74,7 +74,7 @@ 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas@ic3.cat}) - Original code 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@ic3.cat}) - Formatting to CRAN 1.1 - 2015-01 (N.S. Fuckar, \email{neven.fuckar@ic3.cat}) - Upgrade to polynomial fit - i.e., modified to include fit of a raw polynomial with an arbitrary + -> modified to include fit of a raw polynomial with an arbitrary (but constrained) degree/order and a more complete $trend output. } \keyword{datagen} -- GitLab From 81223b862836de2f79432e8ed6fa4dd0e3d7167e Mon Sep 17 00:00:00 2001 From: "Neven S. Fuckar" Date: Mon, 20 Jul 2015 16:36:38 +0200 Subject: [PATCH 07/18] minor changes --- R/Trend.R | 2 +- man/Trend.Rd | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index 4277125f..5f02cc8b 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -2,7 +2,7 @@ Trend <- function(var, posTR = 2, interval = 1, pldegree = 1) { # Computes linear trend or higher order polynomial along the posTR dimension of tensor # var with the method of least squares, and estimates the associated error interval. # Provides also the detrended or residual (after removing fitted polynomial) data. - # The confidence interval relies on a student-T distribution. + # The confidence interval relies on a student-t distribution. # # Args: # var: Matrix of any number of dimensions up to 10. diff --git a/man/Trend.Rd b/man/Trend.Rd index ebab9a6f..fb173def 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -7,10 +7,10 @@ Computes linear trend or higher order polynomial along the posTR dimension of tensor var with the method of least squares, and estimates the associated error interval. Provides also the detrended or residual (after removing fitted polynomial) data. - The confidence interval relies on a student-T distribution. + The confidence interval relies on a student-t distribution. } \usage{ - Trend(var, posTR = 2, interval = 1, pldegree = 1) + Trend(var, posTR=2, interval=1, pldegree=1) } \arguments{ \item{var}{ -- GitLab From 23ec3d40acdc1d4e5f85385975ee44fc79acb45c Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 12 Aug 2019 13:15:58 +0200 Subject: [PATCH 08/18] Eno() + Apply --- NAMESPACE | 1 + R/Eno.R | 166 +++++++++++++++++++++++--------------- man/Eno.Rd | 15 +++- tests/testthat/test-Eno.R | 102 +++++++++++++++++++++++ 4 files changed, 217 insertions(+), 67 deletions(-) create mode 100644 tests/testthat/test-Eno.R diff --git a/NAMESPACE b/NAMESPACE index 3e35710e..21307590 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -115,6 +115,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/Eno.R b/R/Eno.R index 877710cb..207e0007 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,104 @@ #'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 #'@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/man/Eno.Rd b/man/Eno.Rd index ba4f2088..b2a3123d 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/tests/testthat/test-Eno.R b/tests/testthat/test-Eno.R new file mode 100644 index 00000000..e8c7e6cd --- /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) + ) + +}) + -- GitLab From a3b6635fc3fd58d7de8aa3cdf177af9396239bff Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 12 Aug 2019 14:24:43 +0200 Subject: [PATCH 09/18] Add multiApply to import list --- DESCRIPTION | 3 ++- NAMESPACE | 1 + R/Eno.R | 1 + 3 files changed, 4 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 78935b85..87339be5 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 21307590..6c9b99a1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -89,6 +89,7 @@ import(graphics) import(mapproj) import(maps) import(methods) +import(multiApply) import(ncdf4) import(parallel) import(plyr) diff --git a/R/Eno.R b/R/Eno.R index 207e0007..21899e73 100644 --- a/R/Eno.R +++ b/R/Eno.R @@ -54,6 +54,7 @@ #'PlotEquiMap(eno, sampleData$lon, sampleData$lat) #' #'@importFrom stats acf na.pass na.fail +#'@import multiApply #'@export Eno <- function(obs, posdim, na.action = na.pass, ncores = NULL) { -- GitLab From 4707ede91f63bd975cde1fcd0b055bfbc273d093 Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 13 Aug 2019 14:10:56 +0200 Subject: [PATCH 10/18] Update to master (v2.8.5). Add polynomial feature in Trend() and correct its documentation. Trend() Example bugfix. --- R/Consist_Trend.R | 4 ++-- R/Trend.R | 35 +++++++++++++++++------------------ man/Trend.Rd | 35 +++++++++++++++++------------------ 3 files changed, 36 insertions(+), 38 deletions(-) diff --git a/R/Consist_Trend.R b/R/Consist_Trend.R index 827440fc..4932c955 100644 --- a/R/Consist_Trend.R +++ b/R/Consist_Trend.R @@ -107,8 +107,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/Trend.R b/R/Trend.R index 2dd16180..cfbc7ad6 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -65,27 +65,26 @@ #'2.8.6 - 2019-08 (A. Ho, \email{an.ho@bsc.es}) - Adapt polynomial regression #' #'@examples -#' # A test of the default linear fit: -#' startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101') -#' sampleData <- Load('tos', c('i00k'), c('ERSST'), startDates, nleadtime = 124, -#' leadtimemin = 1, leadtimemax = 60, output = 'areave', -#' latmin = 30, latmax = 45, lonmin = 0, lonmax = 40) -#' months_between_startdates <- 60 -#' trend <- Trend(sampleData$obs, 3, months_between_startdates) -#'##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') -> there is an error in this step, but eps file is produced? -#' PlotAno(trend$detrended, NULL, startDates, -#' toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', -#' legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') # there are warnings in this step! +#'example(Load) +#'months_between_startdates <- 60 +#'trend <- Trend(sampleData$obs, 3, months_between_startdates) +#'tmp <- Subset(trend$trend, c('dataset', 'member', 'sdate', 'ftime'), +#' list(1, 1, 1:4, 1), drop = FALSE) +#'PlotVsLTime(tmp, 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') +#' #' # A simple test of a nonlinear fit: #' y0 <- c(10.1, 8.8, 7.3, 4.12, -0.314, -5, -10.6, -18.2) -#' y1 <- Trend(var=y0, posTR=1) # linear fit -#' y2 <- Trend(var=y0, posTR=1, pldegree=2) # quadratic fit +#' y1 <- Trend(var = y0, posTR = 1) # linear fit +#' y2 <- Trend(var = y0, posTR = 1, pldegree = 2) # quadratic fit #' plot(y0) -#' lines(y0 - y1$detrended, col='blue') # linear fit -#' lines(y0 - y2$detrended, col='red') # quadratic fit +#' lines(y0 - y1$detrended, col = 'blue') # linear fit +#' lines(y0 - y2$detrended, col = 'red') # quadratic fit #' #'@rdname Trend #'@export diff --git a/man/Trend.Rd b/man/Trend.Rd index 2529ae10..17784565 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -76,27 +76,26 @@ The confidence interval relies on a student-T distribution.\cr\cr as input (exp). } \examples{ - # A test of the default linear fit: - startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101') - sampleData <- Load('tos', c('i00k'), c('ERSST'), startDates, nleadtime = 124, - leadtimemin = 1, leadtimemax = 60, output = 'areave', - latmin = 30, latmax = 45, lonmin = 0, lonmax = 40) - months_between_startdates <- 60 - trend <- Trend(sampleData$obs, 3, months_between_startdates) -##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') -> there is an error in this step, but eps file is produced? - PlotAno(trend$detrended, NULL, startDates, - toptitle = 'detrended anomalies (along the startdates)', ytitle = 'K', - legends = 'ERSST', biglab = FALSE, fileout = 'tos_detrended_obs.eps') # there are warnings in this step! +example(Load) +months_between_startdates <- 60 +trend <- Trend(sampleData$obs, 3, months_between_startdates) +tmp <- Subset(trend$trend, c('dataset', 'member', 'sdate', 'ftime'), + list(1, 1, 1:4, 1), drop = FALSE) +PlotVsLTime(tmp, 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') + # A simple test of a nonlinear fit: y0 <- c(10.1, 8.8, 7.3, 4.12, -0.314, -5, -10.6, -18.2) - y1 <- Trend(var=y0, posTR=1) # linear fit - y2 <- Trend(var=y0, posTR=1, pldegree=2) # quadratic fit + y1 <- Trend(var = y0, posTR = 1) # linear fit + y2 <- Trend(var = y0, posTR = 1, pldegree = 2) # quadratic fit plot(y0) - lines(y0 - y1$detrended, col='blue') # linear fit - lines(y0 - y2$detrended, col='red') # quadratic fit + lines(y0 - y1$detrended, col = 'blue') # linear fit + lines(y0 - y2$detrended, col = 'red') # quadratic fit } \author{ -- GitLab From a13bc53e551aed92941bcc04e71162466b5c5ffa Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 4 Dec 2019 17:51:17 +0100 Subject: [PATCH 11/18] Use multiApply --- R/Trend.R | 271 +++++++++++++++++++++++++++--------------------------- 1 file changed, 136 insertions(+), 135 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index badb1f98..81c22de1 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -9,7 +9,7 @@ #'as input (exp). #' #'@param var An array of any number of dimensions up to 10. -#'@param posTR An integer indicating the position along which to compute the +#'@param posdim 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. @@ -21,14 +21,16 @@ #' #'@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 +#' An array with same dimensions as parameter 'var' except along the posdim #' 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. #'} +#'\item{$conf.int}{ +#' Corresponding to the limits of the \code{siglev}\% confidence interval +#' (only present if \code{conf = TRUE}) for the slope coefficient. +#'} #'\item{$detrended}{ #' Same dimensions as var with linearly detrended var along the posTR #' dimension. @@ -44,6 +46,7 @@ #'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 to multiApply() #'@examples #'# Load sample data as in Load() example: #'example(Load) @@ -61,146 +64,144 @@ #' #'@rdname Trend #'@export -Trend <- function(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE, pldegree = 1) { - # - # Enlarge the size of var to 10 and move posTR to first position - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - dimsvar <- dim(var) - if (is.null(dimsvar)) { - dimsvar <- length(var) - } - enlvar <- Enlarge(var, 10) - outdim <- c(dimsvar, array(1, dim = (10 - length(dimsvar)))) - - - if (conf) { - nvals <- 3 * (pldegree + 1) #6, 9, 12, ... - } else { - nvals <- 3 + pldegree #4, 5, 6, ... - } - - outdim[posTR] <- nvals - - posaperm <- 1:10 - posaperm[posTR] <- 1 - posaperm[1] <- posTR - enlvar <- aperm(enlvar, posaperm) - dimsaperm <- outdim[posaperm] - - if (dim(enlvar)[1] < (pldegree + 1)) { - stop("Parameter 'pldegree' is too high.") - } - - # - # 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] +Trend <- function(var, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = TRUE, pldegree = 1) { - if (any(!is.na(tmp))) { - mon <- seq(tmp) * interval - lm.out <- lm(tmp ~ poly(mon, degree = pldegree, raw=TRUE), - na.action = na.omit) - - if (conf) { - - for (jj in (pldegree + 1):1) { - - jj1 <- 1 + (pldegree + 1) - jj - - if (jj == 1) { #last, deg = 0 - - enltrend[3 * jj1 - 2, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- lm.out$coefficients[jj] - enltrend[3 * jj1 - 1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 1] - enltrend[3 * jj1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 2] - - } else { - - enltrend[3 * jj1 - 2, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 1] - enltrend[3 * jj1 - 1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- lm.out$coefficients[jj] - enltrend[3 * jj1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 2] - - } - } - } else { #conf = FALSE - - for (jj in (pldegree + 1):1) { + # Check inputs + ## var + if (is.null(var)) { + stop("Parameter 'var' cannot be NULL.") + } + if (!is.numeric(var)) { + stop("Parameter 'var' must be a numeric array.") + } + if (is.null(dim(var))) { #is vector + dim(var) <- c(length(var)) + names(dim(var)) <- sdate_dim + } + if(any(is.null(names(dim(var))))) { + stop("Parameter 'var' must have dimension names.") + } + ## sdate_dim + if (!is.character(sdate_dim)) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% names(dim(var))) { + stop("Parameter 'sdate_dim' is not found in 'var' dimension.") + } + if (length(sdate_dim) > 1) { + sdate_dim <- sdate_dim[1] + warning("Parameter 'sdate_dim' has length > 1 and only the first ", + "element will be used.") + } + ## interval + if (!is.numeric(interval) | interval %% 1 != 0 | interval < 0) { + stop("Parameter 'interval' must be a positive integer.") + } + if (length(interval) > 1) { + interval <- interval[1] + warning("Parameter 'interval' has length > 1 and only the first ", + "element will be used.") + } + ## siglev + if (!is.numeric(siglev) | siglev < 0 | siglev > 1) { + stop("Parameter 'siglev' must be a numeric number between 0 and 1.") + } + if (length(siglev) > 1) { + siglev <- siglev[1] + warning("Parameter 'siglev' has length > 1 and only the first ", + "element will be used.") + } + ## conf + if (!is.logical(conf)) { + stop("Parameter 'conf' must be logical.") + } + if (length(conf) > 1) { + conf <- conf[1] + warning("Parameter 'conf' has length > 1 and only the first ", + "element will be used.") + } + ## pldegree + if (!is.numeric(pldegree) | pldegree %% 1 != 0 | pldegree < 0) { + stop("Parameter 'pldegree' must be a positive integer.") + } + if (length(pldegree) > 1) { + pldegree <- pldegree[1] + warning("Parameter 'pldegree' has length > 1 and only the first ", + "element will be used.") + } - jj1 <- 1 + (pldegree + 1) - jj + ############################### + # Calculate Trend + dim_names <- names(dim(var)) + + output <- Apply(list(var), + target_dims = sdate_dim, + fun = .Trend, + output_dims = , + interval = interval, siglev = siglev, conf = conf, + pldegree = pldegree) + + reorder <- function(output) { + + # Add dim name back + if (is.null(dim(output))) { + dim(output) <- c(sdate = 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] <- sdate_dim + } else { + names(dim(output)) <- sdate_dim + } + } else { # more than one dim + if (names(dim(output))[1] != "") { + dim(output) <- c(1, dim(output)) + names(dim(output))[1] <- sdate_dim + } else { #regular case + names(dim(output))[1] <- sdate_dim + } + } + } + # reorder + pos <- match(dim_names, names(dim(output))) + output <- aperm(output, pos) + names(dim(output)) <- dim_names - if (jj == 1) { + return(output) + } + output <- lapply(output, reorder) - enltrend[jj1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- lm.out$coefficients[jj] - enltrend[jj1 + 1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 1] - enltrend[jj1 + 2, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- confint(lm.out)[jj, 2] + return(output) +} - } else { - enltrend[jj1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- lm.out$coefficients[jj] - } - } - } +.Trend <- function(x, interval = interval, siglev = siglev, conf = conf, pldegree = pldegree) { + if (any(!is.na(x))) { + mon <- seq(x) * interval + lm.out <- lm(x ~ poly(mon, degree = pldegree, raw = TRUE), na.action = na.omit) + trend <- lm.out$coefficients #intercept, slope1, slope2,... - enldetrend[is.na(tmp) == FALSE, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- tmp[is.na(tmp) == FALSE] - lm.out$fitted.values + if (conf) { + conf.lower <- confint(lm.out, level = siglev)[-1, 1] #limits for slope + conf.upper <- confint(lm.out, level = siglev)[-1, 2] #limits for slope + } - } - } - } - } - } - } - } - } + detrended <- c() + detrended[is.na(x) == FALSE] <- x[is.na(x) == FALSE] - lm.out$fitted.values + } else { + trend <- NA + detrend <- NA + if (conf) { + conf.lower <- NA + conf.upper <- NA } } - # - # 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)) -} + if (conf) { + return(list(trend = trend, conf.lower = conf.lower, conf.upper = conf.upper, detrended = detrended)) + } else { + return(list(trend = trend, detrended = detrended)) + } -#'@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]) - if (conf) { - conf.int <- confint(lm.out, level = siglev)[2, 1:2] - } - detrend <- ensmean[is.na(ensmean) == FALSE] - lm.out$fitted.values - } - - # - # Outputs - # ~~~~~~~~~ - # - invisible(list(trend = trend, conf.int = conf.int, detrended = detrend)) } -- GitLab From 9b4f0700266be7e2a9172ecb7109d798049e680e Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 5 Dec 2019 14:09:06 +0100 Subject: [PATCH 12/18] Finish Trend() multiApply-development. update Trend documentation except examples --- NAMESPACE | 1 - R/Trend.R | 218 +++++++++++++++++++++++---------------------------- man/Trend.Rd | 76 +++++++++--------- 3 files changed, 138 insertions(+), 157 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 4e855789..d4ff8d01 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) diff --git a/R/Trend.R b/R/Trend.R index 81c22de1..804d6a18 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -1,45 +1,46 @@ -#'Computes the Trend of the Ensemble Mean +#'Compute the Trend of the Ensemble Mean #' -#'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). -#' -#'@param var An array of any number of dimensions up to 10. -#'@param posdim 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. +#'Compute the polynomial regression along the forecast time with any degree of +#'polynomial. 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 data An numeric array including the dimension along which the trend is computed. +#'@param sdate_dim A character string indicating the dimension along which to +#' compute the trend. +#'@param interval A positive numeric indicating the unit length between two +#' points along 'sdate_dim' dimension. Set 1 as default. +#'@param siglev A numeric indicating the confidence level for the +#' regression computation. Set 0.95 as default. +#'@param conf A logical value indicating whether to retrieve the confidence +#' intervals or not. Set TRUE as default. +#'@param polydeg A positive integer indicating the degree of polynomial +#' regression. Set 1 as default. #'@return #'\item{$trend}{ -#' An array with same dimensions as parameter 'var' except along the posdim -#' 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 same dimensions as parameter 'data' except the +#' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing +#' the regression coefficients. The length of the 'stats' dimension should be +#' \code{polydeg + 1}. +#'} +#'\item{$conf.lower}{ +#' A numeric array with same dimensions as parameter 'data' except the +#' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing +#' the lower limit of the \code{siglev}\% confidence interval for all +#' the regression coefficients expect the intercept term. The length of the +#' 'stats' dimension should be \code{polydeg}. Only present if \code{conf = TRUE}. +#'} +#'\item{$conf.upper}{ +#' A numeric array with same dimensions as parameter 'data' except the +#' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing +#' the upper limit of the \code{siglev}\% confidence interval for all +#' the regression coefficients expect the intercept term. The length of the +#' 'stats' dimension should be \code{polydeg}. Only present if \code{conf = TRUE}. #'} -#'\item{$conf.int}{ -#' Corresponding to the limits of the \code{siglev}\% confidence interval -#' (only present if \code{conf = TRUE}) for the slope coefficient. -#'} #'\item{$detrended}{ -#' Same dimensions as var with linearly detrended var along the posTR -#' dimension. +#' A numeric array with the same dimensions as paramter 'data', the detrended +#' values along the 'sdate_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 @@ -64,127 +65,103 @@ #' #'@rdname Trend #'@export -Trend <- function(var, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = TRUE, pldegree = 1) { +Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = TRUE, polydeg = 1) { # Check inputs - ## var - if (is.null(var)) { - stop("Parameter 'var' cannot be NULL.") - } - if (!is.numeric(var)) { - stop("Parameter 'var' must be a numeric array.") - } - if (is.null(dim(var))) { #is vector - dim(var) <- c(length(var)) - names(dim(var)) <- sdate_dim - } - if(any(is.null(names(dim(var))))) { - stop("Parameter 'var' must have dimension names.") - } + ## data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") + } + if (is.null(dim(data))) { #is vector + dim(data) <- c(length(data)) + names(dim(data)) <- sdate_dim + } + if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) { + stop("Parameter 'data' must have dimension names.") + } ## sdate_dim - if (!is.character(sdate_dim)) { - stop("Parameter 'sdate_dim' must be a character string.") - } - if (!sdate_dim %in% names(dim(var))) { - stop("Parameter 'sdate_dim' is not found in 'var' dimension.") - } - if (length(sdate_dim) > 1) { - sdate_dim <- sdate_dim[1] - warning("Parameter 'sdate_dim' has length > 1 and only the first ", - "element will be used.") - } + if (!is.character(sdate_dim) | length(sdate_dim) > 1) { + stop("Parameter 'sdate_dim' must be a character string.") + } + if (!sdate_dim %in% names(dim(data))) { + stop("Parameter 'sdate_dim' is not found in 'data' dimension.") + } ## interval - if (!is.numeric(interval) | interval %% 1 != 0 | interval < 0) { - stop("Parameter 'interval' must be a positive integer.") - } - if (length(interval) > 1) { - interval <- interval[1] - warning("Parameter 'interval' has length > 1 and only the first ", - "element will be used.") - } + if (any(!is.numeric(interval) | interval <= 0 | length(interval) > 1)) { + stop("Parameter 'interval' must be a positive number.") + } ## siglev - if (!is.numeric(siglev) | siglev < 0 | siglev > 1) { - stop("Parameter 'siglev' must be a numeric number between 0 and 1.") - } - if (length(siglev) > 1) { - siglev <- siglev[1] - warning("Parameter 'siglev' has length > 1 and only the first ", - "element will be used.") - } + if (!is.numeric(siglev) | siglev < 0 | siglev > 1 | length(siglev) > 1) { + stop("Parameter 'siglev' must be a numeric number between 0 and 1.") + } ## conf - if (!is.logical(conf)) { - stop("Parameter 'conf' must be logical.") - } - if (length(conf) > 1) { - conf <- conf[1] - warning("Parameter 'conf' has length > 1 and only the first ", - "element will be used.") - } - ## pldegree - if (!is.numeric(pldegree) | pldegree %% 1 != 0 | pldegree < 0) { - stop("Parameter 'pldegree' must be a positive integer.") - } - if (length(pldegree) > 1) { - pldegree <- pldegree[1] - warning("Parameter 'pldegree' has length > 1 and only the first ", - "element will be used.") - } + if (!is.logical(conf) | length(conf) > 1) { + stop("Parameter 'conf' must be one logical value.") + } + ## polydeg + if (!is.numeric(polydeg) | polydeg %% 1 != 0 | polydeg < 0 | + length(polydeg) > 1) { + stop("Parameter 'polydeg' must be a positive integer.") + } - ############################### - # Calculate Trend - dim_names <- names(dim(var)) + ############################### + # Calculate Trend + dim_names <- names(dim(data)) - output <- Apply(list(var), - target_dims = sdate_dim, - fun = .Trend, - output_dims = , - interval = interval, siglev = siglev, conf = conf, - pldegree = pldegree) + output <- Apply(list(data), + target_dims = sdate_dim, + fun = .Trend, + output_dims = , + interval = interval, siglev = siglev, conf = conf, + polydeg = polydeg) reorder <- function(output) { - # Add dim name back + # Add dim name back if (is.null(dim(output))) { - dim(output) <- c(sdate = length(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] <- sdate_dim + names(dim(output))[1] <- 'stats' } else { - names(dim(output)) <- sdate_dim + names(dim(output)) <- 'stats' } } else { # more than one dim if (names(dim(output))[1] != "") { dim(output) <- c(1, dim(output)) - names(dim(output))[1] <- sdate_dim + names(dim(output))[1] <- 'stats' } else { #regular case - names(dim(output))[1] <- sdate_dim + names(dim(output))[1] <- 'stats' } } } - # reorder - pos <- match(dim_names, names(dim(output))) - output <- aperm(output, pos) - names(dim(output)) <- dim_names + # reorder + pos <- match(dim_names, names(dim(output))) + output <- aperm(output, pos) + names(dim(output)) <- dim_names - return(output) - } - output <- lapply(output, reorder) + return(output) + } + output <- lapply(output, reorder) return(output) } -.Trend <- function(x, interval = interval, siglev = siglev, conf = conf, pldegree = pldegree) { +.Trend <- function(x, interval = interval, siglev = siglev, conf = conf, polydeg = polydeg) { if (any(!is.na(x))) { mon <- seq(x) * interval - lm.out <- lm(x ~ poly(mon, degree = pldegree, raw = TRUE), na.action = na.omit) + lm.out <- lm(x ~ poly(mon, degree = polydeg, raw = TRUE), na.action = na.omit) trend <- lm.out$coefficients #intercept, slope1, slope2,... if (conf) { - conf.lower <- confint(lm.out, level = siglev)[-1, 1] #limits for slope - conf.upper <- confint(lm.out, level = siglev)[-1, 2] #limits for slope + conf.lower <- confint(lm.out, level = siglev)[-1, 1] #limits for slope + conf.upper <- confint(lm.out, level = siglev)[-1, 2] #limits for slope } detrended <- c() @@ -205,3 +182,4 @@ Trend <- function(var, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = } } + diff --git a/man/Trend.Rd b/man/Trend.Rd index 46be5c89..8d3e9a20 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -1,59 +1,61 @@ % 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 of the Ensemble Mean} \usage{ -Trend(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) - -.Trend(exp, interval = 1, siglev = 0.95, conf = TRUE) +Trend(data, sdate_dim = "sdate", interval = 1, siglev = 0.95, + conf = TRUE, polydeg = 1) } \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{posTR}{An integer indicating the position along which to compute the -trend.} +\item{sdate_dim}{A character string indicating the dimension along which to +compute the trend.} -\item{interval}{A number of months/years between 2 points along posTR -dimension. Set 1 as default.} +\item{interval}{A positive numeric indicating the unit length between two +points along 'sdate_dim' dimension. Set 1 as default.} -\item{siglev}{A numeric value indicating the confidence level for the -computation of confidence interval. Set 0.95 as default.} +\item{siglev}{A numeric indicating the confidence level for the +regression computation. Set 0.95 as default.} -\item{conf}{A logical value indicating whether to compute the confidence -levels or not. Set TRUE as default.} +\item{conf}{A logical value indicating whether to retrieve the confidence +intervals or not. Set TRUE as default.} -\item{exp}{An M by N matrix representing M forecasts from N ensemble members.} +\item{polydeg}{A positive integer indicating the degree of polynomial +regression. Set 1 as default.} } \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 same dimensions as parameter 'data' except the + 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing + the regression coefficients. The length of the 'stats' dimension should be + \code{polydeg + 1}. } -\item{$detrended}{ - Same dimensions as var with linearly detrended var along the posTR - dimension. +\item{$conf.lower}{ + A numeric array with same dimensions as parameter 'data' except the + 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing + the lower limit of the \code{siglev}\% confidence interval for all + the regression coefficients expect the intercept term. The length of the + 'stats' dimension should be \code{polydeg}. Only present if \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 same dimensions as parameter 'data' except the + 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing + the upper limit of the \code{siglev}\% confidence interval for all + the regression coefficients expect the intercept term. The length of the + 'stats' dimension should be \code{polydeg}. Only present if \code{conf = TRUE}. +} +\item{$detrended}{ + A numeric array with the same dimensions as paramter 'data', the detrended + values along the 'sdate_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 polynomial regression along the forecast time with any degree of +polynomial. 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: @@ -76,5 +78,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 to multiApply() } \keyword{datagen} + -- GitLab From 7d173fe9df6fa6254cd24e5bfd208165763da128 Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 5 Dec 2019 17:27:57 +0100 Subject: [PATCH 13/18] Revise Trend() and create unit test --- NAMESPACE | 1 + R/Trend.R | 39 +++++------ man/Trend.Rd | 22 ++----- tests/testthat/test-Trend.R | 127 ++++++++++++++++++++++++++++++++++++ 4 files changed, 152 insertions(+), 37 deletions(-) create mode 100644 tests/testthat/test-Trend.R diff --git a/NAMESPACE b/NAMESPACE index d4ff8d01..a44102f8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -89,6 +89,7 @@ import(graphics) import(mapproj) import(maps) import(methods) +import(multiApply) import(ncdf4) import(parallel) import(plyr) diff --git a/R/Trend.R b/R/Trend.R index 804d6a18..fdbbe5b3 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -20,22 +20,23 @@ #'\item{$trend}{ #' A numeric array with same dimensions as parameter 'data' except the #' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing -#' the regression coefficients. The length of the 'stats' dimension should be +#' the regression coefficients from the lowest order (i.e., intercept) to +#' the highest degree. The length of the 'stats' dimension should be #' \code{polydeg + 1}. #'} #'\item{$conf.lower}{ #' A numeric array with same dimensions as parameter 'data' except the #' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing #' the lower limit of the \code{siglev}\% confidence interval for all -#' the regression coefficients expect the intercept term. The length of the -#' 'stats' dimension should be \code{polydeg}. Only present if \code{conf = TRUE}. +#' the regression coefficients with the same order as $trend. The length of the +#' 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. #'} #'\item{$conf.upper}{ #' A numeric array with same dimensions as parameter 'data' except the #' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing #' the upper limit of the \code{siglev}\% confidence interval for all -#' the regression coefficients expect the intercept term. The length of the -#' 'stats' dimension should be \code{polydeg}. Only present if \code{conf = TRUE}. +#' the regression coefficients with the same order as $trend. The length of the +#' 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. #'} #'\item{$detrended}{ #' A numeric array with the same dimensions as paramter 'data', the detrended @@ -52,18 +53,10 @@ #'# 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(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = TRUE, polydeg = 1) { @@ -128,16 +121,16 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = if (length(dim(output)) == 1) { if (!is.null(names(dim(output)))) { dim(output) <- c(1, dim(output)) - names(dim(output))[1] <- 'stats' + names(dim(output))[1] <- sdate_dim } else { - names(dim(output)) <- 'stats' + names(dim(output)) <- sdate_dim } } else { # more than one dim if (names(dim(output))[1] != "") { dim(output) <- c(1, dim(output)) - names(dim(output))[1] <- 'stats' + names(dim(output))[1] <- sdate_dim } else { #regular case - names(dim(output))[1] <- 'stats' + names(dim(output))[1] <- sdate_dim } } } @@ -145,6 +138,7 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = pos <- match(dim_names, names(dim(output))) output <- aperm(output, pos) names(dim(output)) <- dim_names + names(dim(output))[names(dim(output)) == sdate_dim] <- 'stats' return(output) } @@ -160,8 +154,8 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = trend <- lm.out$coefficients #intercept, slope1, slope2,... if (conf) { - conf.lower <- confint(lm.out, level = siglev)[-1, 1] #limits for slope - conf.upper <- confint(lm.out, level = siglev)[-1, 2] #limits for slope + conf.lower <- confint(lm.out, level = siglev)[, 1] + conf.upper <- confint(lm.out, level = siglev)[, 2] } detrended <- c() @@ -176,7 +170,8 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = } if (conf) { - return(list(trend = trend, conf.lower = conf.lower, conf.upper = conf.upper, detrended = detrended)) + return(list(trend = trend, conf.lower = conf.lower, conf.upper = conf.upper, + detrended = detrended)) } else { return(list(trend = trend, detrended = detrended)) } diff --git a/man/Trend.Rd b/man/Trend.Rd index 8d3e9a20..50ecdca4 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -29,22 +29,23 @@ regression. Set 1 as default.} \item{$trend}{ A numeric array with same dimensions as parameter 'data' except the 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing - the regression coefficients. The length of the 'stats' dimension should be + the regression coefficients from the lowest order (i.e., intercept) to + the highest degree. The length of the 'stats' dimension should be \code{polydeg + 1}. } \item{$conf.lower}{ A numeric array with same dimensions as parameter 'data' except the 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing the lower limit of the \code{siglev}\% confidence interval for all - the regression coefficients expect the intercept term. The length of the - 'stats' dimension should be \code{polydeg}. Only present if \code{conf = TRUE}. + the regression coefficients with the same order as $trend. The length of the + 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. } \item{$conf.upper}{ A numeric array with same dimensions as parameter 'data' except the 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing the upper limit of the \code{siglev}\% confidence interval for all - the regression coefficients expect the intercept term. The length of the - 'stats' dimension should be \code{polydeg}. Only present if \code{conf = TRUE}. + the regression coefficients with the same order as $trend. The length of the + 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. } \item{$detrended}{ A numeric array with the same dimensions as paramter 'data', the detrended @@ -61,16 +62,7 @@ The confidence interval relies on the student-T distribution.\cr\cr # 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{ diff --git a/tests/testthat/test-Trend.R b/tests/testthat/test-Trend.R new file mode 100644 index 00000000..a9f4c631 --- /dev/null +++ b/tests/testthat/test-Trend.R @@ -0,0 +1,127 @@ +context("s2dverification::Trend tests") + +############################################## + dat1 <- array(c(-5, -7, -10:10, 12, 11, 7, 16), + dim = c(dat = 1, sdate = 13, ftime = 2)) + + set.seed(10) + dat2 <- c(1:10) + rnorm(10) +############################################## + +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, sdate_dim = 'a'), + "Parameter 'sdate_dim' is not found in 'data' dimension." + ) + + expect_error( + Trend(array(c(1:25), dim = c(dat = 1, date = 5, ftime = 5))), + "Parameter 'sdate_dim' is not found in 'data' dimension." + ) + + expect_error( + Trend(dat1, sdate_dim = 2), + "Parameter 'sdate_dim' must be a character string." + ) + + expect_error( + Trend(dat1, sdate_dim = c('a','sdate')), + "Parameter 'sdate_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(data, polydeg = 3.5), + "Parameter 'polydeg' must be a positive integer." + ) + +}) + +############################################## +test_that("2. Output checks: dat1", { + + expect_equal( + Trend(dat1)$trend, + array(c(-9.769, 0.659, 0.962, 0.797), + dim = c(dat = 1, stats = 2, ftime = 2)), + tolerance = 0.001 + ) + + expect_equal( + dim(Trend(dat1)$conf.upper), + c(dat = 1, stats = 2, ftime = 2) + ) + + expect_equal( + as.vector(summary(Trend(dat1)$detrended)), + c(-4.319, -0.707, 0.115, 0.000, 0.809, 4.110), + 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(stats = 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(stats = 10))), + tolerance = 0.001 + ) + + expect_equal( + length(Trend(dat2, conf = F)), + 2 + ) + + expect_equal( + names(Trend(dat2, conf = F)), + c('trend', 'detrended') + ) + +}) + -- GitLab From a0ea78827f42761f7d4a68bcceeb8ef6d8bea1af Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 9 Dec 2019 14:39:19 +0100 Subject: [PATCH 14/18] Add multiApply to DESCRIPTION --- DESCRIPTION | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0045265c..2fd9876e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -50,7 +50,8 @@ Imports: ncdf4, parallel, plyr, - SpecsVerification (>= 0.5.0) + SpecsVerification (>= 0.5.0), + multiApply Suggests: easyVerification, testthat -- GitLab From 9da101a272fcb9552b95ec15cb10efc644790b4b Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 18 Dec 2019 10:08:26 +0100 Subject: [PATCH 15/18] Trend small fix --- R/Trend.R | 70 +++++++++++++++++++++++++++------------------------- man/Trend.Rd | 3 ++- 2 files changed, 38 insertions(+), 35 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index fdbbe5b3..388afd69 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -5,7 +5,8 @@ #'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 data An numeric array including the dimension along which the trend is computed. +#'@param data An numeric array including the dimension along which the trend +#' is computed. #'@param sdate_dim A character string indicating the dimension along which to #' compute the trend. #'@param interval A positive numeric indicating the unit length between two @@ -107,42 +108,11 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = output <- Apply(list(data), target_dims = sdate_dim, fun = .Trend, - output_dims = , + #output_dims = , interval = interval, siglev = siglev, conf = conf, polydeg = polydeg) - reorder <- function(output) { - - # 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] <- sdate_dim - } else { - names(dim(output)) <- sdate_dim - } - } else { # more than one dim - if (names(dim(output))[1] != "") { - dim(output) <- c(1, dim(output)) - names(dim(output))[1] <- sdate_dim - } else { #regular case - names(dim(output))[1] <- sdate_dim - } - } - } - # reorder - pos <- match(dim_names, names(dim(output))) - output <- aperm(output, pos) - names(dim(output)) <- dim_names - names(dim(output))[names(dim(output)) == sdate_dim] <- 'stats' - - return(output) - } - output <- lapply(output, reorder) + output <- lapply(output, .reorder) return(output) } @@ -178,3 +148,35 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = } + .reorder <- function(output) { + + # 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] <- sdate_dim + } else { + names(dim(output)) <- sdate_dim + } + } else { # more than one dim + if (names(dim(output))[1] != "") { + dim(output) <- c(1, dim(output)) + names(dim(output))[1] <- sdate_dim + } else { #regular case + names(dim(output))[1] <- sdate_dim + } + } + } + # reorder + pos <- match(dim_names, names(dim(output))) + output <- aperm(output, pos) + names(dim(output)) <- dim_names + names(dim(output))[names(dim(output)) == sdate_dim] <- 'stats' + + return(output) + } + diff --git a/man/Trend.Rd b/man/Trend.Rd index 50ecdca4..bc2d9526 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -8,7 +8,8 @@ Trend(data, sdate_dim = "sdate", interval = 1, siglev = 0.95, conf = TRUE, polydeg = 1) } \arguments{ -\item{data}{An numeric array including the dimension along which the trend is computed.} +\item{data}{An numeric array including the dimension along which the trend +is computed.} \item{sdate_dim}{A character string indicating the dimension along which to compute the trend.} -- GitLab From 68ae7fd57db4b9112782266bb811b87c6ad10800 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 18 Dec 2019 10:34:02 +0100 Subject: [PATCH 16/18] Trend() fix for .reorder() --- R/Trend.R | 4 ++-- tests/testthat/test-Trend.R | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index 388afd69..dce13528 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -112,7 +112,7 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = interval = interval, siglev = siglev, conf = conf, polydeg = polydeg) - output <- lapply(output, .reorder) + output <- lapply(output, .reorder, sdate_dim = sdate_dim, dim_names = dim_names) return(output) } @@ -148,7 +148,7 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = } - .reorder <- function(output) { + .reorder <- function(output, sdate_dim, dim_names) { # Add dim name back if (is.null(dim(output))) { diff --git a/tests/testthat/test-Trend.R b/tests/testthat/test-Trend.R index a9f4c631..73e888f2 100644 --- a/tests/testthat/test-Trend.R +++ b/tests/testthat/test-Trend.R @@ -61,7 +61,7 @@ test_that("1. Input checks", { ) expect_error( - Trend(data, polydeg = 3.5), + Trend(dat1, polydeg = 3.5), "Parameter 'polydeg' must be a positive integer." ) -- GitLab From e580d006e5405421097b1702d13ed5cef78e541b Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 18 Dec 2019 16:42:24 +0100 Subject: [PATCH 17/18] Add param 'ncores' --- R/Trend.R | 15 +++++++++++++-- tests/testthat/test-Trend.R | 4 ++++ 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index dce13528..6148f84c 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -17,6 +17,9 @@ #' intervals or not. Set TRUE as default. #'@param polydeg A positive integer indicating the degree of polynomial #' regression. Set 1 as default. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. Default value is NULL. +#' #'@return #'\item{$trend}{ #' A numeric array with same dimensions as parameter 'data' except the @@ -59,7 +62,7 @@ #'@rdname Trend #'@import multiApply #'@export -Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = TRUE, polydeg = 1) { +Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = TRUE, polydeg = 1, ncores = NULL) { # Check inputs ## data @@ -100,6 +103,13 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = length(polydeg) > 1) { stop("Parameter 'polydeg' must be a positive integer.") } + ## 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 Trend @@ -110,7 +120,8 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = fun = .Trend, #output_dims = , interval = interval, siglev = siglev, conf = conf, - polydeg = polydeg) + polydeg = polydeg, + ncores = ncores) output <- lapply(output, .reorder, sdate_dim = sdate_dim, dim_names = dim_names) diff --git a/tests/testthat/test-Trend.R b/tests/testthat/test-Trend.R index 73e888f2..98f5d22a 100644 --- a/tests/testthat/test-Trend.R +++ b/tests/testthat/test-Trend.R @@ -65,6 +65,10 @@ test_that("1. Input checks", { "Parameter 'polydeg' must be a positive integer." ) + expect_error( + Trend(dat1, ncore = 3.5), + "Parameter 'ncores' must be a positive integer." + ) }) ############################################## -- GitLab From 394be4cc32dadf1bcf1f3a17f5d8f219b8d557c5 Mon Sep 17 00:00:00 2001 From: aho Date: Wed, 29 Jan 2020 18:24:10 +0100 Subject: [PATCH 18/18] Fix NA bug. Add one more unit test --- R/Trend.R | 150 +++++++++++++++++++++--------------- man/Trend.Rd | 67 ++++++++-------- tests/testthat/test-Trend.R | 76 ++++++++++-------- 3 files changed, 166 insertions(+), 127 deletions(-) diff --git a/R/Trend.R b/R/Trend.R index 6148f84c..6bcfe546 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -1,50 +1,52 @@ -#'Compute the Trend of the Ensemble Mean +#'Compute the trend #' -#'Compute the polynomial regression along the forecast time with any degree of -#'polynomial. It returns the regression coefficients (including the intercept) +#'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 data An numeric array including the dimension along which the trend #' is computed. -#'@param sdate_dim A character string indicating the dimension along which to -#' compute the trend. +#'@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 'sdate_dim' dimension. Set 1 as default. -#'@param siglev A numeric indicating the confidence level for the -#' regression computation. Set 0.95 as default. -#'@param conf A logical value indicating whether to retrieve the confidence -#' intervals or not. Set TRUE as default. +#' points along 'time_dim' dimension. The default value is 1. #'@param polydeg A positive integer indicating the degree of polynomial -#' regression. Set 1 as default. +#' 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. Default value is NULL. +#' computation. The default value is NULL. #' #'@return #'\item{$trend}{ -#' A numeric array with same dimensions as parameter 'data' except the -#' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing -#' the regression coefficients from the lowest order (i.e., intercept) to -#' the highest degree. The length of the 'stats' dimension should be -#' \code{polydeg + 1}. +#' 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 same dimensions as parameter 'data' except the -#' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing -#' the lower limit of the \code{siglev}\% confidence interval for all -#' the regression coefficients with the same order as $trend. The length of the -#' 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. +#' 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 same dimensions as parameter 'data' except the -#' 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing -#' the upper limit of the \code{siglev}\% confidence interval for all -#' the regression coefficients with the same order as $trend. The length of the -#' 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. +#' 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', the detrended -#' values along the 'sdate_dim' dimension. +#' A numeric array with the same dimensions as paramter 'data', containing the +#' detrended values along the 'time_dim' dimension. #'} #' #'@keywords datagen @@ -52,7 +54,7 @@ #'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 to multiApply() +#'3.0 - 2019-12 (A. Ho, \email{an.ho@bsc.es}) - Adapt multiApply feature #'@examples #'# Load sample data as in Load() example: #'example(Load) @@ -62,7 +64,8 @@ #'@rdname Trend #'@import multiApply #'@export -Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = TRUE, polydeg = 1, ncores = NULL) { +Trend <- function(data, time_dim = 'sdate', interval = 1, polydeg = 1, + conf = TRUE, conf.lev = 0.95, ncores = NULL) { # Check inputs ## data @@ -74,34 +77,34 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = } if (is.null(dim(data))) { #is vector dim(data) <- c(length(data)) - names(dim(data)) <- sdate_dim + 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.") } - ## sdate_dim - if (!is.character(sdate_dim) | length(sdate_dim) > 1) { - stop("Parameter 'sdate_dim' must be a character string.") + ## time_dim + if (!is.character(time_dim) | length(time_dim) > 1) { + stop("Parameter 'time_dim' must be a character string.") } - if (!sdate_dim %in% names(dim(data))) { - stop("Parameter 'sdate_dim' is not found in 'data' dimension.") + 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.") } - ## siglev - if (!is.numeric(siglev) | siglev < 0 | siglev > 1 | length(siglev) > 1) { - stop("Parameter 'siglev' must be a numeric number between 0 and 1.") + ## 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.") } - ## polydeg - if (!is.numeric(polydeg) | polydeg %% 1 != 0 | polydeg < 0 | - length(polydeg) > 1) { - stop("Parameter 'polydeg' must be a positive integer.") + ## 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)) { @@ -115,38 +118,57 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = # 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 = sdate_dim, + target_dims = time_dim, fun = .Trend, - #output_dims = , - interval = interval, siglev = siglev, conf = conf, - polydeg = polydeg, + output_dims = output_dims, + time_dim = time_dim, interval = interval, + polydeg = polydeg, conf = conf, + conf.lev = conf.lev, ncores = ncores) - output <- lapply(output, .reorder, sdate_dim = sdate_dim, dim_names = dim_names) + #output <- lapply(output, .reorder, time_dim = time_dim, dim_names = dim_names) return(output) } -.Trend <- function(x, interval = interval, siglev = siglev, conf = conf, polydeg = polydeg) { - if (any(!is.na(x))) { - mon <- seq(x) * interval - lm.out <- lm(x ~ poly(mon, degree = polydeg, raw = TRUE), na.action = na.omit) +.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 = siglev)[, 1] - conf.upper <- confint(lm.out, level = siglev)[, 2] + 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 <- NA + trend <- rep(NA, polydeg + 1) detrend <- NA if (conf) { - conf.lower <- NA - conf.upper <- NA + conf.lower <- rep(NA, polydeg + 1) + conf.upper <- rep(NA, polydeg + 1) } } @@ -159,7 +181,7 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = } - .reorder <- function(output, sdate_dim, dim_names) { + .reorder <- function(output, time_dim, dim_names) { # Add dim name back if (is.null(dim(output))) { @@ -169,16 +191,16 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = if (length(dim(output)) == 1) { if (!is.null(names(dim(output)))) { dim(output) <- c(1, dim(output)) - names(dim(output))[1] <- sdate_dim + names(dim(output))[1] <- time_dim } else { - names(dim(output)) <- sdate_dim + 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] <- sdate_dim + names(dim(output))[1] <- time_dim } else { #regular case - names(dim(output))[1] <- sdate_dim + names(dim(output))[1] <- time_dim } } } @@ -186,7 +208,7 @@ Trend <- function(data, sdate_dim = 'sdate', interval = 1, siglev = 0.95, conf = pos <- match(dim_names, names(dim(output))) output <- aperm(output, pos) names(dim(output)) <- dim_names - names(dim(output))[names(dim(output)) == sdate_dim] <- 'stats' + names(dim(output))[names(dim(output)) == time_dim] <- 'stats' return(output) } diff --git a/man/Trend.Rd b/man/Trend.Rd index bc2d9526..f058009a 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -2,60 +2,65 @@ % Please edit documentation in R/Trend.R \name{Trend} \alias{Trend} -\title{Compute the Trend of the Ensemble Mean} +\title{Compute the trend} \usage{ -Trend(data, sdate_dim = "sdate", interval = 1, siglev = 0.95, - conf = TRUE, polydeg = 1) +Trend(data, time_dim = "sdate", interval = 1, polydeg = 1, conf = TRUE, + conf.lev = 0.95, ncores = NULL) } \arguments{ \item{data}{An numeric array including the dimension along which the trend is computed.} -\item{sdate_dim}{A character string indicating the dimension along which to -compute the trend.} +\item{time_dim}{A character string indicating the dimension along which to +compute the trend. The default value is 'sdate'.} \item{interval}{A positive numeric indicating the unit length between two -points along 'sdate_dim' dimension. Set 1 as default.} +points along 'time_dim' dimension. The default value is 1.} -\item{siglev}{A numeric indicating the confidence level for the -regression computation. Set 0.95 as default.} +\item{polydeg}{A positive integer indicating the degree of polynomial +regression. The default value is 1.} \item{conf}{A logical value indicating whether to retrieve the confidence -intervals or not. Set TRUE as default.} +intervals or not. The default value is TRUE.} -\item{polydeg}{A positive integer indicating the degree of polynomial -regression. Set 1 as default.} +\item{conf.lev}{A numeric indicating the confidence level for the +regression computation. The default value is 0.95.} + +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} } \value{ \item{$trend}{ - A numeric array with same dimensions as parameter 'data' except the - 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing - the regression coefficients from the lowest order (i.e., intercept) to - the highest degree. The length of the 'stats' dimension should be - \code{polydeg + 1}. + 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 same dimensions as parameter 'data' except the - 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing - the lower limit of the \code{siglev}\% confidence interval for all - the regression coefficients with the same order as $trend. The length of the - 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. + 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 same dimensions as parameter 'data' except the - 'sdate_dim' dimension, which is replaced by a 'stats' dimension containing - the upper limit of the \code{siglev}\% confidence interval for all - the regression coefficients with the same order as $trend. The length of the - 'stats' dimension should be \code{polydeg + 1}. Only present if \code{conf = TRUE}. + 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', the detrended - values along the 'sdate_dim' dimension. + A numeric array with the same dimensions as paramter 'data', containing the + detrended values along the 'time_dim' dimension. } } \description{ -Compute the polynomial regression along the forecast time with any degree of -polynomial. It returns the regression coefficients (including the intercept) +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 } @@ -71,7 +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 to multiApply() +3.0 - 2019-12 (A. Ho, \email{an.ho@bsc.es}) - Adapt multiApply feature } \keyword{datagen} diff --git a/tests/testthat/test-Trend.R b/tests/testthat/test-Trend.R index 98f5d22a..8843b9ed 100644 --- a/tests/testthat/test-Trend.R +++ b/tests/testthat/test-Trend.R @@ -1,11 +1,22 @@ 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", { @@ -14,57 +25,46 @@ test_that("1. Input checks", { 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, sdate_dim = 'a'), - "Parameter 'sdate_dim' is not found in 'data' dimension." + 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 'sdate_dim' is not found in 'data' dimension." + "Parameter 'time_dim' is not found in 'data' dimension." ) - expect_error( - Trend(dat1, sdate_dim = 2), - "Parameter 'sdate_dim' must be a character string." + Trend(dat1, time_dim = 2), + "Parameter 'time_dim' must be a character string." ) - expect_error( - Trend(dat1, sdate_dim = c('a','sdate')), - "Parameter 'sdate_dim' must be a character string." + 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." @@ -76,19 +76,19 @@ test_that("2. Output checks: dat1", { expect_equal( Trend(dat1)$trend, - array(c(-9.769, 0.659, 0.962, 0.797), - dim = c(dat = 1, stats = 2, ftime = 2)), - tolerance = 0.001 + array(c(-9.7692308, 0.6593407, 0.9615385, 0.7967033), + dim = c(stats = 2, dat = 1, ftime = 2)), + tolerance = 0.0001 ) - expect_equal( - dim(Trend(dat1)$conf.upper), - c(dat = 1, stats = 2, ftime = 2) + 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( - as.vector(summary(Trend(dat1)$detrended)), - c(-4.319, -0.707, 0.115, 0.000, 0.809, 4.110), + summary(Trend(dat1)$detrended)[3], + c(Median = 0.1154), tolerance = 0.001 ) @@ -103,25 +103,22 @@ test_that("3. Output checks: dat2", { 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(stats = 10))), + -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(stats = 10))), + -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') @@ -129,3 +126,18 @@ test_that("3. Output checks: dat2", { }) +############################################## +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) + ) + +}) + -- GitLab