
kendallSeasonalTrendTest(y, ...)
## S3 method for class 'formula':
kendallSeasonalTrendTest(y, data = NULL, subset,
na.action = na.pass, ...)
## S3 method for class 'default':
kendallSeasonalTrendTest(y, season, year,
alternative = "two.sided", correct = TRUE, ci.slope = TRUE, conf.level = 0.95,
independent.obs = TRUE, data.name = NULL, season.name = NULL, year.name = NULL,
parent.of.data = NULL, subset.expression = NULL, ...)
## S3 method for class 'data.frame':
kendallSeasonalTrendTest(y, ...)
## S3 method for class 'matrix':
kendallSeasonalTrendTest(y, ...)
y
must be numeric vector of observations.
When y
is a data frame, all columns must be numeric.
When y
is a matrix, it mas.data.frame
to a data frame) containing the variables in the model.
If not found in data
, the variables are taken from environment(for
NA
s.
The default is na.pass
.y
were taken. The length of season
must equal the length of y
.y
were taken.
The length of year
must equal the length of y
."two.sided"
(tau not equal to 0; the default),
"less"
(tau less than 0), and "greater"
(tau greater than 0).TRUE
.TRUE
.0.95
.y
are seially independent.
The default value is TRUE
.deparse(substitute(y))
.deparse(substitute(season))
.deparse(substitute(year))
."htest"
containing the results of the hypothesis
test. See the help file for htest.object
for details.
In addition, the following components are part of the list returned by
kendallSeasonalTrendTest
:independent.obs=TRUE
.independent.obs=FALSE
.kendallTrendTest
) that allows for seasonality in observations collected over time.
They call this test the seasonal Kendall test. Their test is appropriate for testing for
trend in each season when the trend is always in the same direction across all seasons.
van Belle and Hughes (1984) introduced a heterogeneity test for trend which is appropriate for testing
for trend in any direction in any season. Hirsch and Slack (1984) proposed an extension to the seasonal
Kendall test that allows for serial dependence in the observations. The function
kendallSeasonalTrendTest
includes all of these tests, as well as an extension of the
van Belle-Hughes heterogeneity test to the case of serial dependence.
Testing for Trend Assuming Serial Independence
The Model
Assume observations are taken over two or more years, and assume a single year
can be divided into two or more seasons. Let $p$ denote the number of seasons.
Let $X$ and $Y$ denote two continuous random variables with some joint
(bivariate) distribution (which may differ from season to season). Let $N_j$
denote the number of bivariate observations taken in the $j$'th season (over two
or more years) ($j = 1, 2, \ldots, p$), so that
kendallTrendTest
).
Also, assume all of the $Y$ observations are independent.
The function kendallSeasonalTrendTest
assumes that the $X$ values always
denote the year in which the observation was taken. Note that within any season,
the $X$ values need not be unique. That is, there may be more than one
observation within the same year within the same season. In this case, the
argument y
must be a numeric vector, and you must supply the additional
arguments season
and year
.
If there is only one observation per season per year (missing values allowed), it is
usually easiest to supply the argument y
as an $n \times p$ matrix or
data frame, where $n$ denotes the number of years. In this case
sign
function:
kendallTrendTest
).
For each season, if the predictor variables (the $X$'s) are strictly increasing
(e.g., Equation (3) above), then Equation (10) simplifies to
kendallTrendTest
).
In the case of tied observations,
kendallTrendTest
).
Estimating $\tau$, Slope, and Intercept
The function kendall.SeasonalTrendTest
returns estimated values of
Kendall's $\tau$, the slope, and the intercept for each season, as well as a
single estimate for each of these three quantities combined over all seasons.
The overall estimate of $\tau$ is the weighted average of the p seasonal
$\tau$'s:
kendallTrendTest
).
We can compute the estimated slope for season $j$ as:
kendallTrendTest
).
The overall estimate of intercept is the median of the $p$ seasonal estimates of
intercept:
kendallTrendTest
).
Confidence Interval for the Slope
Gilbert (1987, p.227-228) extends his method of computing a confidence interval for
the slope to the case of seasonal observations. Let $N'$ denote the number of
defined two-point estimated slopes that are used in Equation (22) above and let
kendallSeasonalTrendTest
does.
The Van Belle-Hughes Heterogeneity Test for Trend
The seasonal Kendall test described above is appropriate for testing the null
hypothesis (6) against the alternative hypothesis of a trend in at least one season.
All of the trends in each season should be in the same direction.
The seasonal Kendall test is not appropriate for testing for trend when there are
trends in a positive direction in one or more seasons and also negative trends in
one or more seasons. For example, for the following set of observations, the
seasonal Kendall statistic $S'$ is 0 with an associated two-sided p-value of 1,
even though there is clearly a positive trend in season 1 and a negative trend in
season 2.
kendallTrendTest
, htest.object
, cor.test
.# Reproduce Example 14-10 on page 14-38 of USEPA (2009). This example
# tests for trend in analyte concentrations (ppm) collected monthly
# between 1983 and 1985.
head(EPA.09.Ex.14.8.df)
# Month Year Unadj.Conc Adj.Conc
#1 January 1983 1.99 2.11
#2 February 1983 2.10 2.14
#3 March 1983 2.12 2.10
#4 April 1983 2.12 2.13
#5 May 1983 2.11 2.12
#6 June 1983 2.15 2.12
tail(EPA.09.Ex.14.8.df)
# Month Year Unadj.Conc Adj.Conc
#31 July 1985 2.31 2.23
#32 August 1985 2.32 2.24
#33 September 1985 2.28 2.23
#34 October 1985 2.22 2.24
#35 November 1985 2.19 2.25
#36 December 1985 2.22 2.23
# Plot the data
#--------------
Unadj.Conc <- EPA.09.Ex.14.8.df$Unadj.Conc
Adj.Conc <- EPA.09.Ex.14.8.df$Adj.Conc
Month <- EPA.09.Ex.14.8.df$Month
Year <- EPA.09.Ex.14.8.df$Year
Time <- paste(substring(Month, 1, 3), Year - 1900, sep = "-")
n <- length(Unadj.Conc)
Three.Yr.Mean <- mean(Unadj.Conc)
dev.new()
par(mar = c(7, 4, 3, 1) + 0.1, cex.lab = 1.25)
plot(1:n, Unadj.Conc, type = "n", xaxt = "n",
xlab = "Time (Month)",
ylab = "ANALYTE CONCENTRATION (mg/L)",
main = "Figure 14-15. Seasonal Time Series Over a Three Year Period",
cex.main = 1.1)
axis(1, at = 1:n, labels = rep("", n))
at <- rep(c(1, 5, 9), 3) + rep(c(0, 12, 24), each = 3)
axis(1, at = at, labels = Time[at])
points(1:n, Unadj.Conc, pch = 0, type = "o", lwd = 2)
points(1:n, Adj.Conc, pch = 3, type = "o", col = 8, lwd = 2)
abline(h = Three.Yr.Mean, lwd = 2)
legend("topleft", c("Unadjusted", "Adjusted", "3-Year Mean"), bty = "n",
pch = c(0, 3, -1), lty = c(1, 1, 1), lwd = 2, col = c(1, 8, 1),
inset = c(0.05, 0.01))
# Perform the seasonal Kendall trend test
#----------------------------------------
kendallSeasonalTrendTest(Unadj.Conc ~ Month + Year,
data = EPA.09.Ex.14.8.df)
#Results of Hypothesis Test
#--------------------------
#
#Null Hypothesis: All 12 values of tau = 0
#
#Alternative Hypothesis: The seasonal taus are not all equal
# (Chi-Square Heterogeneity Test)
# At least one seasonal tau != 0
# and all non-zero tau's have the
# same sign (z Trend Test)
#
#Test Name: Seasonal Kendall Test for Trend
# (with continuity correction)
#
#Estimated Parameter(s): tau = 0.9722222
# slope = 0.0600000
# intercept = -131.7350000
#
#Estimation Method: tau: Weighted Average of
# Seasonal Estimates
# slope: Hirsch et al.'s
# Modification of
# Thiel/Sen Estimator
# intercept: Median of
# Seasonal Estimates
#
#Data: y = Unadj.Conc
# season = Month
# year = Year
#
#Data Source: EPA.09.Ex.14.8.df
#
#Sample Sizes: January = 3
# February = 3
# March = 3
# April = 3
# May = 3
# June = 3
# July = 3
# August = 3
# September = 3
# October = 3
# November = 3
# December = 3
# Total = 36
#
#Test Statistics: Chi-Square (Het) = 0.1071882
# z (Trend) = 5.1849514
#
#Test Statistic Parameter: df = 11
#
#P-values: Chi-Square (Het) = 1.000000e+00
# z (Trend) = 2.160712e-07
#
#Confidence Interval for: slope
#
#Confidence Interval Method: Gilbert's Modification of
# Theil/Sen Method
#
#Confidence Interval Type: two-sided
#
#Confidence Level: 95%
#
#Confidence Interval: LCL = 0.05786914
# UCL = 0.07213086
#==========
# Clean up
#---------
rm(Unadj.Conc, Adj.Conc, Month, Year, Time, n, Three.Yr.Mean, at)
graphics.off()
Run the code above in your browser using DataLab