Version: | 1.1.6 |
Date: | 2023-10-10 |
Title: | Non-Parametric Trend Tests and Change-Point Detection |
Depends: | R (≥ 3.0) |
Description: | The analysis of environmental data often requires the detection of trends and change-points. This package includes tests for trend detection (Cox-Stuart Trend Test, Mann-Kendall Trend Test, (correlated) Hirsch-Slack Test, partial Mann-Kendall Trend Test, multivariate (multisite) Mann-Kendall Trend Test, (Seasonal) Sen's slope, partial Pearson and Spearman correlation trend test), change-point detection (Lanzante's test procedures, Pettitt's test, Buishand Range Test, Buishand U Test, Standard Normal Homogeinity Test), detection of non-randomness (Wallis-Moore Phase Frequency Test, Bartels rank von Neumann's ratio test, Wald-Wolfowitz Test) and the two sample Robust Rank-Order Distributional Test. |
Imports: | extraDistr (≥ 1.8.0) |
Suggests: | strucchange, Kendall, psych, datasets |
Classification/MSC-2010: | 62F03, 62G10, 62M10, 62P12 |
License: | GPL-3 |
Encoding: | UTF-8 |
ZipData: | yes |
NeedsCompilation: | yes |
RoxygenNote: | 7.2.2 |
Packaged: | 2023-10-10 07:57:22 UTC; thorsten |
Author: | Thorsten Pohlert |
Maintainer: | Thorsten Pohlert <thorsten.pohlert@gmx.de> |
Repository: | CRAN |
Date/Publication: | 2023-10-10 08:50:02 UTC |
Simulated data of Page (1955) as test-example for change-point detection
Description
Simulated data of Page (1955) as test-example for change-point detection taken from Table 1 of Pettitt (1979)
Usage
data(PagesData)
Format
a vector that contains 40 elements
Details
According to the publication of Pettitt (1979), the series comprise a
significant p = 0.014
change-point at i = 17
. The function
pettitt.test
computes the same U statistics as given by
Pettitt (1979) in Table1, row 4.
References
Page, E. S. (1954), A test for a change in a parameter occuring at an unknown point. Biometrika 41, 100–114.
Pettitt, A. N., (1979). A non-parametric approach to the change point problem. Journal of the Royal Statistical Society Series C, Applied Statistics 28, 126–135.
See Also
Examples
data(PagesData)
pettitt.test(PagesData)
Bartels Test for Randomness
Description
Performes a rank version of von Neumann's ratio test as proposed by Bartels. The null hypothesis of randomness is tested against the alternative hypothesis
Usage
bartels.test(x)
Arguments
x |
a vector of class "numeric" or a time series object of class "ts" |
Details
In this function, the test is implemented as given by Bartels (1982),
where the ranks r_1, \ldots, r_n
of
the X_i, \ldots, X_n
are used for the statistic:
T = \frac{\sum_{i=1}^n (r_i - r_{i+1})^2}{\sum_{i=1}^n (r_i - \bar{r})^2}
As proposed by Bartels (1982), the p
-value is calculated
for sample sizes in the range of (10 \le n < 100)
with
the non-standard beta distribution for the range 0 \le x \le 4
with parameters:
a = b = \frac{5 n \left( n + 1\right) \left(n - 1\right)^2}
{2 \left(n - 2\right) \left(5n^2 - 2n - 9\right)} - \frac{1}{2}
For sample sizes n \ge 100
a normal approximation with
N(2, 20/(5n + 7))
is used for p
-value calculation.
Value
A list with class "htest"
data.name |
character string that denotes the input data |
p.value |
the p-value |
statistic |
the test statistic |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
Note
The current function is for complete observations only.
References
R. Bartels (1982), The Rank Version of von Neumann's Ratio Test for Randomness, Journal of the American Statistical Association 77, 40–46.
See Also
Examples
# Example from Schoenwiese (1992, p. 113)
## Number of frost days in April at Munich from 1957 to 1968
##
frost <- ts(data=c(9,12,4,3,0,4,2,1,4,2,9,7), start=1957)
bartels.test(frost)
## Example from Sachs (1997, p. 486)
x <- c(5,6,2,3,5,6,4,3,7,8,9,7,5,3,4,7,3,5,6,7,8,9)
bartels.test(x)
## Example from Bartels (1982, p. 43)
x <- c(4, 7, 16, 14, 12, 3, 9, 13, 15, 10, 6, 5, 8, 2, 1, 11, 18, 17)
bartels.test(x)
Buishand Range Test for Change-Point Detection
Description
Performes the Buishand range test for change-point detection of a normal variate.
Usage
br.test(x, m = 20000)
Arguments
x |
a vector of class "numeric" or a time series object of class "ts" |
m |
numeric, number of Monte-Carlo replicates, defaults to 20000 |
Details
Let X
denote a normal random variate, then the following model
with a single shift (change-point) can be proposed:
x_i = \left\{
\begin{array}{lcl}
\mu + \epsilon_i, & \qquad & i = 1, \ldots, m \\
\mu + \Delta + \epsilon_i & \qquad & i = m + 1, \ldots, n \\
\end{array} \right.
with \epsilon \approx N(0,\sigma)
. The null hypothesis \Delta = 0
is tested against the alternative \Delta \ne 0
.
In the Buishand range test, the rescaled adjusted partial sums are calculated as
S_k = \sum_{i=1}^k \left(x_i - \hat{x}\right) \qquad (1 \le i \le n)
The test statistic is calculated as:
Rb = \left(\max S_k - \min S_k\right) / \sigma
.
The p.value
is estimated with a Monte Carlo simulation
using m
replicates.
Critical values based on m = 19999
Monte Carlo simulations
are tabulated for Rb / \sqrt{n}
by Buishand (1982).
Value
A list with class "htest" and "cptest"
data.name |
character string that denotes the input data |
p.value |
the p-value |
statistic |
the test statistic |
null.value |
the null hypothesis |
estimates |
the time of the probable change point |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
data |
numeric vector of Sk for plotting |
Note
The current function is for complete observations only.
References
T. A. Buishand (1982), Some Methods for Testing the Homogeneity of Rainfall Records, Journal of Hydrology 58, 11–27.
G. Verstraeten, J. Poesen, G. Demaree, C. Salles (2006), Long-term (105 years) variability in rain erosivity as derived from 10-min rainfall depth data for Ukkel (Brussels, Belgium): Implications for assessing soil erosion rates. Journal of Geophysical Research 111, D22109.
See Also
Examples
data(Nile)
(out <- br.test(Nile))
plot(out)
data(PagesData) ; br.test(PagesData)
Buishand U Test for Change-Point Detection
Description
Performes the Buishand U test for change-point detection of a normal variate.
Usage
bu.test(x, m = 20000)
Arguments
x |
a vector of class "numeric" or a time series object of class "ts" |
m |
numeric, number of Monte-Carlo replicates, defaults to 20000 |
Details
Let X
denote a normal random variate, then the following model
with a single shift (change-point) can be proposed:
x_i = \left\{
\begin{array}{lcl}
\mu + \epsilon_i, & \qquad & i = 1, \ldots, m \\
\mu + \Delta + \epsilon_i & \qquad & i = m + 1, \ldots, n \\
\end{array} \right.
with \epsilon \approx N(0,\sigma)
. The null hypothesis \Delta = 0
is tested against the alternative \Delta \ne 0
.
In the Buishand U test, the rescaled adjusted partial sums are calculated as
S_k = \sum_{i=1}^k \left(x_i - \bar{x}\right) \qquad (1 \le i \le n)
The sample standard deviation is
D_x = \sqrt{n^{-1} \sum_{i=1}^n \left(x_i - \bar{x}\right)}
The test statistic is calculated as:
U = \left[n \left(n + 1 \right) \right]^{-1} \sum_{k=1}^{n-1}
\left(S_k / D_x \right)^2
.
The p.value
is estimated with a Monte Carlo simulation
using m
replicates.
Critical values based on m = 19999
Monte Carlo simulations
are tabulated for U
by Buishand (1982, 1984).
Value
A list with class "htest" and "cptest"
data.name |
character string that denotes the input data |
p.value |
the p-value |
statistic |
the test statistic |
null.value |
the null hypothesis |
estimates |
the time of the probable change point |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
data |
numeric vector of Sk for plotting |
Note
The current function is for complete observations only.
References
T. A. Buishand (1982), Some Methods for Testing the Homogeneity of Rainfall Records, Journal of Hydrology 58, 11–27.
T. A. Buishand (1984), Tests for Detecting a Shift in the Mean of Hydrological Time Series, Journal of Hydrology 73, 51–69.
See Also
Examples
data(Nile)
(out <- bu.test(Nile))
plot(out)
data(PagesData)
bu.test(PagesData)
Cox and Stuart Trend Test
Description
Performes the non-parametric Cox and Stuart trend test
Usage
cs.test(x)
Arguments
x |
a vector or a time series object of class "ts" |
Details
First, the series is devided by three. It is compared, whether the
data of the first third of the series are larger or smaller than the
data of the last third of the series.
The test statistic of the Cox-Stuart trend test for n > 30
is calculated as:
z = \frac{|S - \frac{n}{6}|}{\sqrt{\frac{n}{12}}}
where S
denotes the maximum of the number of signs, i.e.
+
or -
, respectively. The z
-statistic is
normally distributed. For n \le 30
a continuity correction of
-0.5
is included in the denominator.
Value
An object of class "htest"
method |
a character string indicating the chosen test |
data.name |
a character string giving the name(s) of the data |
statistic |
the Cox-Stuart z-value |
alternative |
a character string describing the alternative hypothesis |
p.value |
the p-value for the test |
Note
NA values are omitted. Many ties in the series will lead to reject H0 in the present test.
References
L. Sachs (1997), Angewandte Statistik. Berlin: Springer.
C.-D. Schoenwiese (1992), Praktische Statistik. Berlin: Gebr. Borntraeger.
D. R. Cox and A. Stuart (1955), Quick sign tests for trend in location and dispersion. Biometrika 42, 80-95.
See Also
Examples
## Example from Schoenwiese (1992, p. 114)
## Number of frost days in April at Munich from 1957 to 1968
## z = -0.5, Accept H0
frost <- ts(data=c(9,12,4,3,0,4,2,1,4,2,9,7), start=1957)
cs.test(frost)
## Example from Sachs (1997, p. 486-487)
## z ~ 2.1, Reject H0 on a level of p = 0.0357
x <- c(5,6,2,3,5,6,4,3,7,8,9,7,5,3,4,7,3,5,6,7,8,9)
cs.test(x)
cs.test(Nile)
Correlated Seasonal Mann-Kendall Test
Description
Performs a Seasonal Mann-Kendall test under the presence of correlated seasons.
Usage
csmk.test(x, alternative = c("two.sided", "greater", "less"))
Arguments
x |
a time series object with class |
alternative |
the alternative hypothesis, defaults to |
Details
The Mann-Kendall scores are first computed for each season seperately.
The variance - covariance matrix is computed according to Libiseller and Grimvall (2002).
Finally the corrected Z-statistics for the entire series
is calculated as follows, whereas a continuity correction is employed
for n \le 10
:
z = \frac{\mathbf{1}^T \mathbf{S}}
{\sqrt{\mathbf{1}^T \mathbf{\Gamma}~\mathbf{1}}}
where
z
denotes the quantile of the normal distribution, 1 indicates a vector
with all elements equal to one, \mathbf{S}
is the vector of Mann-Kendall scores
for each season and \mathbf{\Gamma}
denotes the variance - covariance matrix.
Value
An object with class "htest"
data.name |
character string that denotes the input data |
p.value |
the p-value for the entire series |
statistic |
the z quantile of the standard normal distribution for the entire series |
null.value |
the null hypothesis |
estimates |
the estimates S and varS for the entire series |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
cov |
the variance - covariance matrix |
Note
Ties are not corrected. Current Version is for complete observations only.
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
cor
,
cor.test
,
mk.test
,
smk.test
Examples
csmk.test(nottem)
Monthly concentration of particle bound HCB, River Rhine
Description
Time series of monthly concentration of particle bound
Hexachlorobenzene (HCB) in \mu
g/kg at six different
monitoring sites at the River Rhine, 1995.1-2006.12
Usage
data(hcb)
Format
a time series object of class "mts"
we first column, series of station Weil (RKM 164.3)
ka second column, series of station Karlsruhe-Iffezheim (RKM 333.9)
mz third column, series of station Mainz (RKM 498.5)
ko fourth column, series of station Koblenz (RKM 590.3)
bh fith column, series of station Bad Honnef(RKM 645.8)
bi sixth column, series of station Bimmen (RKM 865.0)
Details
NO DATA values in the series were filled with estimated values
using linear interpolation (see approx
.
The Rhine Kilometer (RKM) is in increasing order from source to mouth of the River Rhine.
Source
International Commission for the Protection of the River Rhine
References
T. Pohlert, G. Hillebrand, V. Breitung (2011), Trends of persistent organic pollutants in the suspended matter of the River Rhine, Hydrological Processes 25, 3803–3817. doi:10.1002/hyp.8110
Examples
data(hcb)
plot(hcb)
mult.mk.test(hcb)
Lanzante's Test for Change Point Detection
Description
Performes a non-parametric test after Lanzante in order to test for a shift in the central tendency of a time series. The null hypothesis, no shift, is tested against the alternative, shift.
Usage
lanzante.test(x, method = c("wilcox.test", "rrod.test"))
Arguments
x |
a vector of class "numeric" or a time series object of class "ts" |
method |
the test method. Defaults to |
Details
Let X
denote a continuous random variable, then the following model
with a single shift (change-point) can be proposed:
x_i = \left\{
\begin{array}{lcl}
\theta + \epsilon_i, & \qquad & i = 1, \ldots, m \\
\theta + \Delta + \epsilon_i & \qquad & i = m + 1, \ldots, n \\
\end{array} \right.
with \theta(\epsilon) = 0
. The null hypothesis, H:\Delta = 0
is tested against the alternative A:\Delta \ne 0
.
First, the data are transformed into increasing ranks and for each time-step the adjusted rank sum is computed:
U_k = 2 \sum_{i=1}^k r_i - k \left(n + 1\right) \qquad k = 1, \ldots, n
The probable change point is located at the absolute maximum of the statistic:
m = k(\max |U_k|)
.
For method = "wilcox.test"
the Wilcoxon-Mann-Whitney two-sample
test is performed, using m
to split the series. Otherwise,
the robust rank-order distributional test (rrod.test
is
performed.
Value
A list with class "htest" and "cptest".
References
Lanzante, J. R. (1996), Resistant, robust and non-parametric techniques for the analysis of climate data: Theory and examples, including applications to historical radiosonde station data, Int. J. Clim., 16, 1197–1226.
See Also
Examples
data(maxau) ; plot(maxau[,"s"])
s.res <- lanzante.test(maxau[,"s"])
n <- s.res$nobs
i <- s.res$estimate
s.1 <- mean(maxau[1:i,"s"])
s.2 <- mean(maxau[(i+1):n,"s"])
s <- ts(c(rep(s.1,i), rep(s.2,(n-i))))
tsp(s) <- tsp(maxau[,"s"])
lines(s, lty=2)
print(s.res)
data(PagesData) ; lanzante.test(PagesData)
Annual suspended sediment concentration and flow data, River Rhine
Description
Annual time series of average suspended sediment concentration (s) in mg/l and average discharge (Q) in m^3 / s at the River Rhine, 1965.1-2009.1
Usage
data(maxau)
Format
a time series object of class "mts"
s. first column, suspended sediment concentration
Q. second column, average discharge
Source
Bundesanstalt für Gewässerkunde, Koblenz, Deutschland (Federal Institute of Hydrology, Koblenz, Germany)
Examples
data(maxau)
plot(maxau)
Mann-Kendall Trend Test
Description
Performs the Mann-Kendall Trend Test
Usage
mk.test(x, alternative = c("two.sided", "greater", "less"), continuity = TRUE)
Arguments
x |
a vector of class "numeric" or a time series object of class "ts" |
alternative |
the alternative hypothesis, defaults to |
continuity |
logical, indicates whether a continuity correction
should be applied, defaults to |
Details
The null hypothesis is that the data come from a population with independent realizations and are identically distributed. For the two sided test, the alternative hypothesis is that the data follow a monotonic trend. The Mann-Kendall test statistic is calculated according to:
S = \sum_{k = 1}^{n-1} \sum_{j = k + 1}^n
\mathrm{sgn}\left(x_j - x_k\right)
with \mathrm{sgn}
the signum function (see sign
).
The mean of S
is \mu = 0
. The variance including the
correction term for ties is
\sigma^2 = \left\{n \left(n-1\right)\left(2n+5\right) -
\sum_{j=1}^p t_j\left(t_j - 1\right)\left(2t_j+5\right) \right\} / 18
where p
is the number of the tied groups in the data set and
t_j
is the number of data points in the j
-th tied group.
The statistic S
is approximately normally distributed, with
z = S / \sigma
If continuity = TRUE
then a continuity correction will be employed:
z = \mathrm{sgn}(S) ~ \left(|S| - 1\right) / \sigma
The statistic S
is closely related to Kendall's \tau
:
\tau = S / D
where
D = \left[\frac{1}{2}n\left(n-1\right)-
\frac{1}{2}\sum_{j=1}^p t_j\left(t_j - 1\right)\right]^{1/2}
\left[\frac{1}{2}n\left(n-1\right) \right]^{1/2}
Value
A list with class "htest"
data.name |
character string that denotes the input data |
p.value |
the p-value |
statistic |
the z quantile of the standard normal distribution |
null.value |
the null hypothesis |
estimates |
the estimates S, varS and tau |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
Note
Current Version is for complete observations only.
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
cor.test
,
MannKendall
,
partial.mk.test
,
sens.slope
Examples
data(Nile)
mk.test(Nile, continuity = TRUE)
##
n <- length(Nile)
cor.test(x=(1:n),y=Nile, meth="kendall", continuity = TRUE)
Multivariate (Multisite) Mann-Kendall Test
Description
Performs a Multivariate (Multisite) Mann-Kendall test.
Usage
mult.mk.test(x, alternative = c("two.sided", "greater", "less"))
Arguments
x |
a time series object of class "ts" |
alternative |
the alternative hypothesis, defaults to |
Details
The Mann-Kendall scores are first computed for each variate (side) seperately.
S = \sum_{k = 1}^{n-1} \sum_{j = k + 1}^n
\mathrm{sgn}\left(x_j - x_k\right)
with \mathrm{sgn}
the signum function (see sign
).
The variance - covariance matrix is computed according to Libiseller and Grimvall (2002).
\Gamma_{xy} = \frac{1}{3}
\left[K + 4 \sum_{j=1}^n R_{jx} R_{jy} -
n \left(n + 1 \right) \left(n + 1 \right) \right]
with
K = \sum_{1 \le i < j \le n} \mathrm{sgn} \left\{ \left( x_j - x_i \right)
\left( y_j - y_i \right) \right\}
and
R_{jx} = \left\{ n + 1 + \sum_{i=1}^n
\mathrm{sgn} \left( x_j - x_i \right) \right\} / 2
Finally, the corrected z-statistics for the entire series
is calculated as follows, whereas a continuity correction is employed
for n \le 10
:
z = \frac{\sum_{i=1}^d S_i}{\sqrt{\sum_{j=1}^d \sum_{i=1}^d \Gamma_{ij}}}
where
z
denotes the quantile of the normal distribution
S
is the vector of Mann-Kendall scores
for each variate (site) 1 \le i \le d
and
\Gamma
denotes symmetric variance - covariance matrix.
Value
An object with class "htest"
data.name |
character string that denotes the input data |
p.value |
the p-value for the entire series |
statistic |
the z quantile of the standard normal distribution for the entire series |
null.value |
the null hypothesis |
estimates |
the estimates S and varS for the entire series |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
cov |
the variance - covariance matrix |
Note
Ties are not corrected. Current Version is for complete observations only.
References
Hipel, K.W. and McLeod, A.I. (1994), Time Series Modelling of Water Resources and Environmental Systems. New York: Elsevier Science.
Lettenmeier, D.P. (1988), Multivariate nonparametric tests for trend in water quality. Water Resources Bulletin 24, 505–512.
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
cor
,
cor.test
,
mk.test
,
smk.test
Examples
data(hcb)
mult.mk.test(hcb)
Partial Correlation Trend Test
Description
Performs a partial correlation trend test with either Pearson's or
Spearman's correlation coefficients (r(tx.z)
).
Usage
partial.cor.trend.test(x, z, method = c("pearson", "spearman"))
Arguments
x |
a "vector" or "ts" object that contains the variable, which is tested for trend (i.e. correlated with time) |
z |
a "vector" or "ts" object that contains the co-variate, which will be partialled out |
method |
a character string indicating which correlation coefficient is to be computed. One of "pearson" (default) or "spearman", can be abbreviated. |
Details
This function performs a partial correlation trend test using either the "pearson" correlation coefficient, or the "spearman" rank correlation coefficient (Hipel and McLoed (1994), p. 882). The partial correlation coefficient for the response variable "x" with time "t", when the effect of the explanatory variable "z" is partialled out, is defined as:
r_{tx.z} = \frac{r_{tx} - r_{tz}~r_{xz}}
{\sqrt{1 - r_{tz}^2} ~ \sqrt{1-r_{xz}^2}}
The H0: r_{tx.z} = 0
(i.e. no trend for "x", when
effect of "z" is partialled out) is tested against the
alternate Hypothesis, that there is a trend for "x", when the effect of
"z" is partialled out.
The partial correlation coefficient is tested for significance with
the student t distribution on df = n - 2
degree of freedom.
Value
An object of class "htest"
method |
a character string indicating the chosen test |
data.name |
a character string giving the name(s) of the data |
statistic |
the value of the test statistic |
estimate |
the partial correlation coefficient |
parameter |
the degrees of freedom of the test statistic in the case that it follows a t distribution |
alternative |
a character string describing the alternative hypothesis |
p.value |
the p-value of the test |
null.value |
The value of the null hypothesis |
Note
Current Version is for complete observations only.
References
Hipel, K.W. and McLeod, A.I. (1994), Time Series Modelling of Water Resources and Environmental Systems. New York: Elsevier Science.
Bahrenberg, G., Giese, E. and Nipper, J., (1992): Statistische Methoden in der Geographie, Band 2 Multivariate Statistik, Teubner, Stuttgart.
See Also
cor
,
cor.test
,
partial.r
,
partial.mk.test
,
Examples
data(maxau)
a <- tsp(maxau) ; tt <- a[1]:a[2]
s <- maxau[,"s"] ; Q <- maxau[,"Q"]
maxau.df <- data.frame(Year = tt, s =s, Q = Q)
plot(maxau.df)
partial.cor.trend.test(s,Q, method="pearson")
partial.cor.trend.test(s,Q, method="spearman")
Partial Mann-Kendall Trend Test
Description
Performs a partial Mann-Kendall Trend Test
Usage
partial.mk.test(x, y, alternative = c("two.sided", "greater", "less"))
Arguments
x |
a "vector" or "ts" object that contains the variable, which is tested for trend (i.e. correlated with time) |
y |
a "vector" or "ts" object that contains the variable, which effect on "x" is partialled out |
alternative |
character, the alternative method; defaults to "two.sided" |
Details
According to Libiseller and Grimvall (2002), the test statistic
for x
with its covariate y
is
z = \frac{S_x - r_{xy} S_y}{\left[ \left( 1 - r_{xy}^2 \right)
n \left(n - 1 \right) \left(2 n + 5 \right) / 18 \right]^{0.5}}
where the correlation r
is calculated as:
r_{xy} = \frac{\sigma_{xy}}
{n \left(n - 1\right) \left(2 n + 5 \right) / 18}
The conditional covariance between x
and y
is
\sigma_{xy} = \frac{1}{3}
\left[K + 4 \sum_{j=1}^n R_{jx} R_{jy} -
n \left(n + 1 \right) \left(n + 1 \right) \right]
with
K = \sum_{1 \le i < j \le n} \mathrm{sgn} \left\{ \left( x_j - x_i \right)
\left( y_j - y_i \right) \right\}
and
R_{jx} = \left\{ n + 1 + \sum_{i=1}^n
\mathrm{sgn} \left( x_j - x_i \right) \right\} / 2
Value
A list with class "htest"
method |
a character string indicating the chosen test |
data.name |
a character string giving the name(s) of the data |
statistic |
the value of the test statistic |
estimate |
the Mann-Kendall score S, the variance varS and the correlation between x and y |
alternative |
a character string describing the alternative hypothesis |
p.value |
the p-value of the test |
null.value |
the null hypothesis |
Note
Current Version is for complete observations only. The test statistic is not corrected for ties.
References
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
Examples
data(maxau)
s <- maxau[,"s"]; Q <- maxau[,"Q"]
partial.mk.test(s,Q)
Pettitt's Test for Change-Point Detection
Description
Performes a non-parametric test after Pettitt in order to test for a shift in the central tendency of a time series. The H0-hypothesis, no change, is tested against the HA-Hypothesis, change.
Usage
pettitt.test(x)
Arguments
x |
a vector of class "numeric" or a time series object of class "ts" |
Details
In this function, the test is implemented as given by Verstraeten et. al.
(2006), where the ranks r_1, \ldots, r_n
of
the X_i, \ldots, X_n
are used for the statistic:
U_k = 2 \sum_{i=1}^k r_i - k \left(n + 1\right) \qquad k = 1, \ldots, n
The test statistic is the maximum of the absolute value of the vector:
\hat{U} = \max |U_k|
.
The probable change-point K
is located where \hat{U}
has its maximum.
The approximate probability for a two-sided test is calculated
according to
p = 2 \exp^{-6K^2 / (T^3 + T^2)}
Value
A list with class "htest" and "cptest"
Note
The current function is for complete observations only.
The approximate probability is good for p \le 0.5
.
References
CHR (ed., 2010), Das Abflussregime des Rheins und seiner Nebenfluesse im 20. Jahrhundert, Report no I-22 of the CHR, p. 172.
Pettitt, A. N. (1979), A non-parametric approach to the change point problem. Journal of the Royal Statistical Society Series C, Applied Statistics 28, 126-135.
G. Verstraeten, J. Poesen, G. Demaree, C. Salles (2006), Long-term (105 years) variability in rain erosivity as derived from 10-min rainfall depth data for Ukkel (Brussels, Belgium): Implications for assessing soil erosion rates. Journal of Geophysical Research 111, D22109.
See Also
Examples
data(maxau) ; plot(maxau[,"s"])
s.res <- pettitt.test(maxau[,"s"])
n <- s.res$nobs
i <- s.res$estimate
s.1 <- mean(maxau[1:i,"s"])
s.2 <- mean(maxau[(i+1):n,"s"])
s <- ts(c(rep(s.1,i), rep(s.2,(n-i))))
tsp(s) <- tsp(maxau[,"s"])
lines(s, lty=2)
print(s.res)
data(PagesData) ; pettitt.test(PagesData)
Plotting cptest-objects
Description
Plotting method for objects inheriting from class "cptest"
Usage
## S3 method for class 'cptest'
plot(x, ...)
Arguments
x |
an object of class "cptest" |
... |
further arguments, currently ignored |
Examples
data(Nile)
(out <- br.test(Nile))
par(mfrow=c(2,1))
plot(Nile) ; plot(out)
Robust Rank-Order Distributional Test
Description
Performs Fligner-Pollicello robust rank-order distributional test for location.
Usage
rrod.test(x, ...)
## Default S3 method:
rrod.test(x, y, alternative = c("two.sided", "less", "greater"), ...)
## S3 method for class 'formula'
rrod.test(formula, data, subset, na.action, ...)
Arguments
x |
a vector of data values. |
... |
further arguments to be passed to or from methods. |
y |
an optional numeric vector of data values. |
alternative |
the alternative hypothesis. Defaults to |
formula |
a formula of the form |
data |
an optional matrix or data frame (or similar: see
|
subset |
an optional vector specifying a subset of observations to be used. |
na.action |
a function which indicates what should happen when
the data contain |
Details
The non-parametric RROD two-sample test can be used to test for differences in location, whereas it does not assume variance homogeneity.
Let X
and Y
denote two samples with sizes n_x
and n_y
of a continuous variable.First, the combined sample is transformed
into ranks in increasing order.
Let S_{xi}
and S_{yj}
denote the counts of Y
(X)
values having a lower rank than x_i
(y_j)
. The mean counts are:
\bar{S}_x = \sum_{i=1}^{n_x} S_{xi} / n_x
\bar{S}_y = \sum_{j=1}^{n_y} S_{yj} / n_y
The variances are:
s^2_{Sx} = \sum_{i=1}^{n_x} \left( S_{xi} - \bar{S}_x \right)^2
s^2_{Sy} = \sum_{j=1}^{n_y} \left( S_{yj} - \bar{S}_y \right)^2
The test statistic is:
z = \frac{1}{2}~
\frac{n_x \bar{S}_x - n_y \bar{S}_y}
{\left( \bar{S}_x \bar{S}_y + s^2_{Sx} + s^2_{Sy} \right)^{1/2}}
The two samples have significantly different location parameters,
if |z| > z_{1-\alpha/2}
.
The function calculates the p
-values of the null hypothesis
for the selected alternative than can be "two.sided"
, "greater"
or "less"
.
Value
A list with class "htest"
.
References
Fligner, M. A., Pollicello, G. E. III. (1981), Robust Rank Procedures for the Behrens-Fisher Problem, Journal of the American Statistical Association, 76, 162–168.
Lanzante, J. R. (1996), Resistant, robust and non-parametric techniques for the analysis of climate data: Theory and examples, including applications to historical radiosonde station data, Int. J. Clim., 16, 1197–1226.
Siegel, S. and Castellan, N. (1988), Nonparametric Statistics For The Behavioural Sciences, New York: McCraw-Hill.
See Also
Examples
## Two-sample test.
## Hollander & Wolfe (1973), 69f.
## Permeability constants of the human chorioamnion (a placental
## membrane) at term (x) and between 12 to 26 weeks gestational
## age (y). The alternative of interest is greater permeability
## of the human chorioamnion for the term pregnancy.
x <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)
y <- c(1.15, 0.88, 0.90, 0.74, 1.21)
rrod.test(x, y, alternative = "g")
## Formula interface.
boxplot(Ozone ~ Month, data = airquality)
rrod.test(Ozone ~ Month, data = airquality,
subset = Month %in% c(5, 8))
Seasonal Sen's Slope
Description
Computes seasonal Sen's slope for linear rate of change
Usage
sea.sens.slope(x)
Arguments
x |
a time series object of class "ts" |
Details
Acccording to Hirsch et al. (1982) the seasonal Sen's slope is calculated as follows:
d_{ijk} = \frac{x_{ij} - x_{ik}}{j - k}
for each (x_{ij}, x_{ik})
pair i = 1, \ldots, m
,
where (1 \leq k < j \leq n_i
and n_i
is
the number of known values in the i-th
season.
The seasonal slope estimator is the median of
the d_{ijk}
values.
Value
numeric, Seasonal Sen's slope.
Note
Current Version is for complete observations only.
References
Hipel, K.W. and McLeod, A.I. (1994), Time Series Modelling of Water Resources and Environmental Systems. New York: Elsevier Science.
Hirsch, R., J. Slack, R. Smith (1982), T echniques of Trend Analysis for Monthly Water Quality Data. Water Resources Research 18, 107-121.
Sen, P.K. (1968), Estimates of the regression coefficient based on Kendall's tau, Journal of the American Statistical Association 63, 1379–1389.
See Also
Examples
sea.sens.slope(nottem)
Sen's slope
Description
Computes Sen's slope for linear rate of change and corresponding confidence intervalls
Usage
sens.slope(x, conf.level = 0.95)
Arguments
x |
numeric vector or a time series object of class "ts" |
conf.level |
numeric, the level of significance |
Details
This test computes both the slope (i.e. linear rate of change) and confidence levels according to Sen's method. First, a set of linear slopes is calculated as follows:
d_{k} = \frac{x_j - x_i}{j - i}
for \left(1 \le i < j \le n \right)
, where d
is the slope, x denotes the variable, n is the number of data, and i,
j are indices.
Sen's slope is then calculated as the median from all slopes:
b_{Sen} = \textnormal{median}(d_k)
.
This function also computes the upper and lower confidence limits for sens slope.
Value
A list of class "htest".
estimates |
numeric, Sen's slope |
data.name |
character string that denotes the input data |
p.value |
the p-value |
statistic |
the z quantile of the standard normal distribution |
null.value |
the null hypothesis |
conf.int |
upper and lower confidence limit |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
Note
Current Version is for complete observations only.
References
Hipel, K.W. and McLeod, A.I. (1994), Time Series Modelling of Water Resources and Environmental Systems. New York: Elsevier Science.
Sen, P.K. (1968), Estimates of the regression coefficient based on Kendall's tau, Journal of the American Statistical Association 63, 1379–1389.
Examples
data(maxau)
sens.slope(maxau[,"s"])
mk.test(maxau[,"s"])
Seasonal Mann-Kendall Trend Test
Description
Performs a Seasonal Mann-Kendall Trend Test (Hirsch-Slack Test)
Usage
smk.test(x, alternative = c("two.sided", "greater", "less"), continuity = TRUE)
Arguments
x |
a time series object with class |
alternative |
the alternative hypothesis, defaults to |
continuity |
logical, indicates, whether a continuity correction
should be done; defaults to |
Details
The Mann-Kendall statistic for the $g$-th season is calculated as:
S_g = \sum_{i = 1}^{n-1} \sum_{j = i + 1}^n
\mathrm{sgn}\left(x_{jg} - x_{ig}\right), \qquad (1 \le g \le m)
with \mathrm{sgn}
the signum function (see sign
).
The mean of S_g
is \mu_g = 0
. The variance including the
correction term for ties is
\sigma_g^2 = \left\{n \left(n-1\right)\left(2n+5\right) -
\sum_{j=1}^p t_{jg}\left(t_{jg} - 1\right)\left(2t_{jg}+5\right) \right\} / 18
~~ (1 \le g \le m)
The seasonal Mann-Kendall statistic for the entire series is calculated according to
\begin{array}{ll}
\hat{S} = \sum_{g = 1}^m S_g &
\hat{\sigma}_g^2 = \sum_{g = 1}^m \sigma_g^2
\end{array}
The statistic S_g
is approximately normally distributed, with
z_g = S_g / \sigma_g
If continuity = TRUE
then a continuity correction will be employed:
z = \mathrm{sgn}(S_g) ~ \left(|S_g| - 1\right) / \sigma_g
Value
An object with class "htest" and "smktest"
data.name |
character string that denotes the input data |
p.value |
the p-value for the entire series |
statistic |
the z quantile of the standard normal distribution for the entire series |
null.value |
the null hypothesis |
estimates |
the estimates S and varS for the entire series |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
Sg |
numeric vector that contains S scores for each season |
varSg |
numeric vector that contains varS for each season |
pvalg |
numeric vector that contains p-values for each season |
Zg |
numeric vector that contains z-quantiles for each season |
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.
R. Hirsch, J. Slack, R. Smith (1982), Techniques of Trend Analysis for Monthly Water Quality Data, Water Resources Research 18, 107–121.
Examples
res <- smk.test(nottem)
## print method
res
## summary method
summary(res)
Standard Normal Homogeinity Test (SNHT) for Change-Point Detection
Description
Performes the Standard Normal Homogeinity Test (SNHT) for change-point detection of a normal variate.
Usage
snh.test(x, m = 20000)
Arguments
x |
a vector of class "numeric" or a time series object of class "ts" |
m |
numeric, number of Monte-Carlo replicates, defaults to 20000 |
Details
Let X
denote a normal random variate, then the following model
with a single shift (change-point) can be proposed:
x_i = \left\{
\begin{array}{lcl}
\mu + \epsilon_i, & \qquad & i = 1, \ldots, m \\
\mu + \Delta + \epsilon_i & \qquad & i = m + 1, \ldots, n \\
\end{array} \right.
with \epsilon \approx N(0,\sigma)
. The null hypothesis \Delta = 0
is tested against the alternative \Delta \ne 0
.
The test statistic for the SNHT test is calculated as follows:
T_k = k z_1^2 + \left(n - k\right) z_2^2 \qquad (1 \le k < n)
where
\begin{array}{l l}
z_1 = \frac{1}{k} \sum_{i=1}^k \frac{x_i - \bar{x}}{\sigma} &
z_2 = \frac{1}{n-k} \sum_{i=k+1}^n \frac{x_i - \bar{x}}{\sigma}. \\
\end{array}
The critical value is:
T = \max T_k.
The p.value
is estimated with a Monte Carlo simulation
using m
replicates.
Critical values based on m = 1,000,000
Monte Carlo simulations
are tabulated for T
by Khaliq and Ouarda (2007).
Value
A list with class "htest" and "cptest"
data.name |
character string that denotes the input data |
p.value |
the p-value |
statistic |
the test statistic |
null.value |
the null hypothesis |
estimates |
the time of the probable change point |
alternative |
the alternative hypothesis |
method |
character string that denotes the test |
data |
numeric vector of Tk for plotting |
Note
The current function is for complete observations only.
References
H. Alexandersson (1986), A homogeneity test applied to precipitation data, Journal of Climatology 6, 661–675.
M. N. Khaliq, T. B. M. J. Ouarda (2007), On the critical values of the standard normal homogeneity test (SNHT), International Journal of Climatology 27, 681–687.
G. Verstraeten, J. Poesen, G. Demaree, C. Salles (2006), Long-term (105 years) variability in rain erosivity as derived from 10-min rainfall depth data for Ukkel (Brussels, Belgium): Implications for assessing soil erosion rates. Journal of Geophysical Research 111, D22109.
See Also
Examples
data(Nile)
(out <- snh.test(Nile))
plot(out)
data(PagesData) ; snh.test(PagesData)
Object summaries
Description
Generic function "summary" for objects of class smktest
.
Usage
## S3 method for class 'smktest'
summary(object, ...)
Arguments
object |
an object of class |
... |
further arguments, currently ignored |
Wallis and Moore Phase-Frequency Test
Description
Performes the non-parametric Wallis and Moore phase-frequency test for testing the H0-hypothesis, whether the series comprises random data, against the HA-Hypothesis, that the series is significantly different from randomness (two-sided test).
Usage
wm.test(x)
Arguments
x |
a vector or a time series object of class "ts" |
Details
The test statistic of the phase-frequency test for n > 30
is calculated as:
z = \frac{| h - \frac{2 n - 7}{3}|}{\sqrt{\frac{16 n - 29}{90}}}
where h
denotes the number of phases, whereas the first and the
last phase is not accounted. The z
-statistic is
normally distributed. For n \le 30
a continuity correction of
-0.5
is included in the denominator.
Value
An object of class "htest"
method |
a character string indicating the chosen test |
data.name |
a character string giving the name(s) of the data |
statistic |
the Wallis and Moore z-value |
alternative |
a character string describing the alternative hypothesis |
p.value |
the p-value for the test |
Note
NA values are omitted. Many ties in the series will lead to reject H0 in the present test.
References
L. Sachs (1997), Angewandte Statistik. Berlin: Springer.
C.-D. Schoenwiese (1992), Praktische Statistik. Berlin: Gebr. Borntraeger.
W. A. Wallis and G. H. Moore (1941): A significance test for time series and other ordered observations. Tech. Rep. 1. National Bureau of Economic Research. New York.
See Also
Examples
## Example from Schoenwiese (1992, p. 113)
## Number of frost days in April at Munich from 1957 to 1968
## z = -0.124, Accept H0
frost <- ts(data=c(9,12,4,3,0,4,2,1,4,2,9,7), start=1957)
wm.test(frost)
## Example from Sachs (1997, p. 486)
## z = 2.56, Reject H0 on a level of p < 0.05
x <- c(5,6,2,3,5,6,4,3,7,8,9,7,5,3,4,7,3,5,6,7,8,9)
wm.test(x)
wm.test(nottem)
Wald-Wolfowitz Test for Independence and Stationarity
Description
Performes the non-parametric Wald-Wolfowitz test for independence and stationarity.
Usage
ww.test(x)
Arguments
x |
a vector or a time series object of class "ts" |
Details
Let x_1, x_2, ..., x_n
denote the sampled data, then the test
statistic of the Wald-Wolfowitz test is calculated as:
R = \sum_{i=1}^{n-1} x_i x_{i+1} + x_1 x_n
The expected value of R is:
E(R) = \frac{s_1^2 - s_2}{n - 1}
The expected variance is:
V(R) = \frac{s_2^2 - s_4}{n - 1} - E(R)^2 + \frac{s_1^4 -
4 s_1^2 s_2 + 4 s_1 s_3 + s_2^2 - 2 s_4}{(n - 1) (n - 2)}
with:
s_t = \sum_{i=1}^{n} x_i^t, ~~ t = 1, 2, 3, 4
For n > 10
the test statistic is normally distributed, with:
z = \frac{R - E(R)}{\sqrt{V(R)}}
ww.test calculates p-values from the standard normal distribution for the two-sided case.
Value
An object of class "htest"
method |
a character string indicating the chosen test |
data.name |
a character string giving the name(s) of the data |
statistic |
the Wald-Wolfowitz z-value |
alternative |
a character string describing the alternative hypothesis |
p.value |
the p-value for the test |
Note
NA values are omitted.
References
R. K. Rai, A. Upadhyay, C. S. P. Ojha and L. M. Lye (2013), Statistical analysis of hydro-climatic variables. In: R. Y. Surampalli, T. C. Zhang, C. S. P. Ojha, B. R. Gurjar, R. D. Tyagi and C. M. Kao (ed. 2013), Climate change modelling, mitigation, and adaptation. Reston, VA: ASCE. doi = 10.1061/9780784412718.
A. Wald and J. Wolfowitz (1943), An exact test for randomness in the non-parametric case based on serial correlation. Annual Mathematical Statistics 14, 378–388.
WMO (2009), Guide to Hydrological Practices. Volume II, Management of Water Resources and Application of Hydrological Practices, WMO-No. 168.
Examples
ww.test(nottem)
ww.test(Nile)
set.seed(200)
x <- rnorm(100)
ww.test(x)