rtrend: Trend Estimating Tools
Description
The traditional linear regression trend, Modified Mann-Kendall (MK) non-parameter trend and bootstrap trend are included in this package. Linear regression trend is rewritten by '.lm.fit'. MK trend is rewritten by 'Rcpp'. Finally, those functions are about 10 times faster than previous version in R. Reference: Hamed, K. H., & Rao, A. R. (1998). A modified Mann-Kendall trend test for autocorrelated data. Journal of hydrology, 204(1-4), 182-196. doi:10.1016/S0022-1694(97)00125-X.
Author(s)
Maintainer: Dongdong Kong kongdd.sysu@gmail.com (ORCID)
Authors:
Heyang Song 7800556@gmail.com (ORCID)
See Also
Useful links:
faster autocorrelation based on ffw
Description
This function is 4-times faster than stats::acf()
Usage
acf.fft(x, lag.max = NULL)
Arguments
lag.max
maximum lag at which to calculate the acf.
Default is 10\log_{10}(N/m) where N is the
number of observations and m the number of series. Will
be automatically limited to one less than the number of observations
in the series.
Value
An array with the same dimensions as x containing the estimated autocorrelation.
References
https://github.com/santiagobarreda/phonTools/blob/main/R/fastacf.R
https://gist.github.com/FHedin/05d4d6d74e67922dfad88038b04f621c
https://gist.github.com/ajkluber/f293eefba2f946f47bfa
http://www.tibonihoo.net/literate_musing/autocorrelations.html#wikispecd
https://lingpipe-blog.com/2012/06/08/autocorrelation-fft-kiss-eigen
Examples
set.seed(1)
x <- rnorm(100)
r_fast <- acf.fft(x)
r <- acf(x, plot=FALSE, lag.max=100)$acf[,,1]
apply function for 3d array
Description
NA values will be removed automatically
Usage
apply_3d(
array,
dim = 3,
FUN = rowMeans2,
by = NULL,
scale = 1,
na.rm = TRUE,
...
)
Arguments
array
A 3d array
dim
giving the subscripts to split up data by.
FUN
function, should only be row applied function, e.g. matrixStats::rowMeans2,
matrixStats::rowMins, matrixStats::rowRanges.
Because 3d array will be convert to matrix first, with the aggregated dim in
the last dimension.
by
If not provided (
NULL), the aggregated dim will be disappear. For example, daily precipitation[nrow, ncol, 31-days]aggregate into monthly[nrow, ncol].If provided,
byshould be equal to the aggregateddim. For example, daily precipitation[nrow, ncol, 365-days]aggregate into monthly[nrow, ncol, 12-months]. In that situation,byshould be equal to 365, and beformat(date, '%Y%m').
scale
in the same length of by, or a const value,
value_returned = FUN(x)*scale. This parameter is designed for converting
monthly to yearly, meanwhile multiply days in month.
Currently, same group should have the same scale factor. Otherwise, only the
first is used.
See Also
apply_row matrixStats::rowRanges
Examples
set.seed(1)
size <- c(10, 8, 31)
arr <- array(rnorm(10 * 8 * 31), dim = size)
by <- c(rep(1, 10), rep(2, 21))
r2 <- apply_3d(arr, 3, by = by, FUN = rowMeans)
## Not run:
arr_yearly <- apply_3d(arr, by = year(dates), scale = days_in_month(dates))
## End(Not run)
apply_col
Description
-
apply_col: aggregate by col, return a[ngrp, ncol]matrix -
apply_row: aggregate by row, return a[nrow, ngrp]matrix
Usage
apply_col(mat, by, FUN = colMeans2, scale = 1, ...)
apply_row(mat, by, FUN = rowMeans2, scale = 1, ...)
Arguments
mat
matrix, [nrow, ncol]
by
integer vector, with the dim of [ntime]
scale
in the same length of by, or a const value,
value_returned = FUN(x)*scale. This parameter is designed for converting
monthly to yearly, meanwhile multiply days in month.
Currently, same group should have the same scale factor. Otherwise, only the
first is used.
Details
For example, setting the dimension of mat is [ngrid, ntime],
if you want to aggregate by time, apply_row should be used here;
if you want to aggregate by region (grids), apply_col should be used.
Note
This function also suits for big.matrix object.
Examples
mat <- matrix(rnorm(4 * 6), 4, 6)
mat_bycol <- apply_col(mat, c(1, 1, 2, 2), colMeans)
mat_byrow <- apply_row(mat, c(1, 1, 2, 2, 3, 3), rowMeans)
array_3dTo2d
Description
array_3dTo2d
Usage
array_3dTo2d(array, I_grid = NULL)
array_2dTo3d(array, I_grid = NULL, dim)
Arguments
array
array with the dimension of [nlon, nlat, ntime]
I_grid
subindex of [nrow, ncol]
dim
[nrow, ncol]
Value
[nlat*nlon, ntime]
chunk
Description
chunk
Usage
chunk(x, nchunk = 6)
Arguments
x
a vector or list
nchunk
the number of chunks to be splitted
References
https://stackoverflow.com/questions/3318333/split-a-vector-into-chunks-in-r
parallel apply and llply
Description
parallel apply and llply
Usage
llply_par(X, FUN, ..., byrow = TRUE, .combine = c)
parLapply2(X, FUN, ..., byrow = TRUE, .combine = c)
apply_par(X, .margins = 1, FUN, ..., .progress = "text")
Modified Mann Kendall
Description
If valid observations <= 5, NA will be returned.
Usage
mkTrend_r(y, ci = 0.95, IsPlot = FALSE)
mkTrend(y, x = seq_along(y), ci = 0.95, IsPlot = FALSE)
Arguments
y
numeric vector
ci
critical value of autocorrelation
IsPlot
boolean
x
(optional) numeric vector
Details
mkTrend is 4-fold faster with .lm.fit.
Value
-
Z0: The original (non corrected) Mann-Kendall test Z statistic. -
pval0: The original (non corrected) Mann-Kendall test p-value -
Z: The new Z statistic after applying the correction -
pval: Corrected p-value after accounting for serial autocorrelationN/n*sValue of the correction factor, representing the quotient of the number of samples N divided by the effective sample sizen*s -
slp: Sen slope, The slope of the (linear) trend according to Sen test
Note
slp is significant, if pval < alpha.
Author(s)
Dongdong Kong
References
Hipel, K.W. and McLeod, A.I. (1994), Time Series Modelling of Water Resources and Environmental Systems. New York: Elsevier Science.
Libiseller, C. and Grimvall, A., (2002), Performance of partial Mann-Kendall tests for trend detection in the presence of covariates. Environmetrics 13, 71–84, doi:10.1002/env.507.
See Also
fume::mktrend and trend::mk.test
Examples
x <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69)
r <- mkTrend(x)
r_cpp <- mkTrend(x, IsPlot = TRUE)
movmean
Description
NA and Inf values in the y will be ignored automatically.
Usage
movmean(y, halfwin = 1L, SG_style = FALSE, w = NULL)
movmean2(y, win_left = 1L, win_right = 0L, w = NULL)
movmean_2d(mat, win_left = 3L, win_right = 0L)
Arguments
y
A numeric vector.
halfwin
Integer, half of moving window size
SG_style
If true, head and tail values will be in the style of SG (more weights on the center point), else traditional moving mean style.
w
Corresponding weights of y, with the same length.
win_left, win_right
windows size in the left and right
mat
numeric matrix
Examples
x <- 1:100
x[50] <- NA; x[80] <- Inf
s1 <- movmean(x, 2, SG_style = TRUE)
s2 <- movmean(x, 2, SG_style = FALSE)
movmean2(c(4, 8, 6, -1, -2, -3, -1), 2, 0)
movmean2(c(4, 8, NA, -1, -2, Inf, -1), 2, 0)
Set dimensions of an Object
Description
Set dimensions of an Object
Usage
set_dim(x, dim)
set_dimnames(x, value)
Arguments
x
an R object, for example a matrix, array or data frame.
dim
integer vector, see also base::dim()
value
For the default method, either NULL or
a numeric vector, which is coerced to integer (by truncation).
See Also
Examples
x <- 1:12
set_dim(x, c(3, 4))
slope_arr
Description
slope_arr
Usage
slope_arr(
arr,
fun = rtrend::slope_mk,
return.list = FALSE,
.progress = "text",
...
)
Value
t, A 3d array, with the dim of [nx, ny, 2].
-
t[,,1]: slope -
t[,,2]: pvalue
calculate slope of rast object
Description
calculate slope of rast object
Usage
slope_rast(
r,
period = c(2001, 2020),
outfile = NULL,
fun = rtrend::slope_mk,
...,
overwrite = FALSE,
.progress = "text"
)
rast_filter_time(r, period = c(2001, 2020))
Arguments
r
A yearly rast object, which should have time attribute
period
c(year_begin, year_end)
outfile
The path of outputed tiff file. If specified, slope and
pvalue will be written into outfile.
fun
the function used to calculate slope, see slope() for details.
...
other parameters ignored
overwrite
logical. If TRUE, outfile is overwritten.
.progress
name of the progress bar to use, see
create_progress_bar
Value
A terra rast object, with bands of slope and pvalue.
See Also
Examples
library(rtrend)
library(terra)
f <- system.file("rast/MOD15A2_LAI_China_G050_2001-2020.tif", package = "rtrend")
r <- rast(f)
r
time(r)
slp <- slope_rast(r,
period = c(2001, 2020),
outfile = "LAI_trend.tif", overwrite = TRUE,
fun = rtrend::slope_mk, .progress = "none"
)
# if you want to show progress, set `.progress = "text"`
slp
plot(slp)
file.remove("LAI_trend.tif")
slope
Description
-
slope: linear regression slope -
slope_p: linear regression slope and p-value -
slope_mk: mann kendall Sen's slope and p-value -
slope_sen: same asslope_mk, but with no p-value -
slope_boot: bootstrap slope and p-value
Usage
slope_sen(y, x = NULL)
slope(y, x, ...)
slope_p(y, x, fast = TRUE)
slope_sen_r(y, x = seq_along(y), ...)
slope_mk(y, x = NULL, ...)
slope_boot(y, x = NULL, slope_FUN = slope, times = 100, alpha = 0.1, seed, ...)
Arguments
y
vector of observations of length n, or a matrix with n rows.
x
vector of predictor of length n, or a matrix with n rows.
...
ignored.
fast
Boolean. If true, stats::.lm.fit() will be used, which is 10x
faster than stats::lm() .
slope_FUN
one of slope() , slope_p() , slope_mk()
times
The number of bootstrap replicates.
alpha
significant level, defalt 0.1
seed
a single value, interpreted as an integer, or NULL
(see ‘Details’).
Value
-
slope: linear regression coefficient -
pvalue:p-value <= 0.05`` means that corresponding slope' is significant. -
sd:Std. Error
For slope_boot, slope is estimated in many times. The lower, mean, upper
and standard deviation (sd) are returned.
Examples
y <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69)
r <- slope(y)
r_p <- slope_p(y)
r_mk <- slope_mk(y)
r_boot <- slope_boot(y)
Weighted Savitzky-Golay
Description
NA and Inf values in the y has been ignored automatically.
Usage
smooth_wSG(y, halfwin = 1L, d = 1L, w = NULL)
smooth_SG(y, halfwin = 1L, d = 1L)
Arguments
y
colvec
halfwin
halfwin of Savitzky-Golay
d
polynomial of degree. When d = 1, it becomes moving average.
w
colvec of weight
Examples
y <- c(1, 3, 2, 5, 6, 8, 10, 1)
w <- seq_along(y)/length(y)
halfwin = 2
d = 2
s1 <- smooth_wSG(y, halfwin, d, w)
s2 <- smooth_SG(y, halfwin, d)
split_data
Description
split_data
Usage
split_data(x, nchunk = 6, byrow = TRUE)
Arguments
x
a vector, list, or matrix (not support 3d array)
nchunk
the number of chunks to be splitted
byrow
If TRUE, split by row, otherwise by column.
summary_lm
Description
summary method for class ".lm.fit".. It's 200 times faster than traditional lm.
Usage
summary_lm(obj, ...)
Arguments
obj
Object returned by .lm.fit .
...
ignored
Value
a p x 4 matrix with columns for the estimated coefficient, its standard error, t-statistic and corresponding (two-sided) p-value. Aliased coefficients are omitted.
Examples
set.seed(129)
n <- 100
p <- 2
X <- matrix(rnorm(n * p), n, p) # no intercept!
y <- rnorm(n)
obj <- .lm.fit (x = cbind(1, X), y = y)
info <- summary_lm(obj)