Research Regional Kendall Test for Trend DENNIS R. HELSEL* AND LONNA M. FRANS U.S. Geological Survey, P. O. Box 25046, MS 964, Lakewood, Colorado 80225
Trends in environmental variables are often investigated within a study region at more than one site. At each site, a trend analysis determines whether a trend has occurred. Yet often also of interest is whether a consistent trend is evident throughout the entire region. This paper adapts the Seasonal Kendall trend test to determine whether a consistent regional trend occurs in environmental variables.
Introduction Perhaps the most popular test for trend in environmental studies is the Seasonal Kendall test (1). First applied to concentrations of dissolved chemicals in surface waters (2), the test has been subsequently used to test for trends in biologic community structure (3), radionuclides in groundwaters (4), growth patterns of children (5), turbidity (6), estuarine salinity (7), lake water quality (8), atmospheric chemistry (9), flood and drought statistics (10), and precipitation chemistry (11). The test has been used by scientists in the United States (2), Japan (12), Sweden (13), Germany (14), The Netherlands (15), Norway (16), Australia (17), and Canada (18). The Seasonal Kendall test determines whether there is a monotonic (single-direction) trend over time. It is a nonparametric test, determining trend regardless of whether that trend is linear or whether data follow a normal distribution. Environmental data rarely are sufficiently symmetric to be modeled by a normal distribution. The flexibility to incorporate data of different shapes gives the Kendall procedure a distinct advantage over linear regression, where normality of residuals is required. The seasonal portion of the Kendall test has been important for testing trends in data with temporal variation such as concentrations in surface waters: trend patterns may differ from one month or quarter (“season”) to the next. Trends may occur in warm weather but not cold, or vice versa. To deal with differences between seasons, the Seasonal Kendall test separately tests for trend in each season, and then combines the results into one overall test. This overall test states whether there is a trend with time for that location, blocking out all seasonal differences in the pattern of change. The test is sometimes run after first subtracting off the influences of variables other than time, such as stream flow (19). Media such as soil and groundwater chemistry exhibit more spatial than seasonal variation. Where trends are of interest, they are usually of interest at multiple locations. For example, Stoline et al. (20) reported trends in water chemistry at 15 monitoring wells downgradient of a Superfund landfill site. However, trend results for the entire site were summarized only by inspection, rather than by a formal procedure. Kolpin et al. (21) analyzed trends in groundwater chemistry across Iowa using two-group tests, comparing early * Corresponding author phone: (303) 236-5340; fax: (303) 2361425; e-mail:
[email protected]. 4066 9 ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 13, 2006
versus later concentrations. They stated that “Few studies to date ... adequately address this issue [trend] on a regional scale”. They summarized regional trends by qualitatively describing results for wells grouped into three aquifer types. Again, no overall quantitative test was performed. In both cases, however, discussion pointed toward the need to summarize trend results across large portions of their study areas. Recently we (22) required a test for trend in groundwater nitrate concentrations across regions of the Columbia Basin Ground Water Management Area (GWMA) in Washington state. The management area was formed in 1998 with the primary goal of decreasing the high nitrate concentrations found in the area’s groundwater (23). The GWMA monitors progress toward meeting this goal by analyzing nitrate concentrations at over 400 wells every 2 years. A quantitative method was needed to summarize trends within the GWMA, and to compare trends among its three counties. To meet these objectives, we have adapted the Seasonal Kendall test to determine whether a consistent pattern of trend occurs across an entire area, at multiple locations. We call this test the Regional Kendall test for trend. Mann-Kendall and Seasonal Kendall Tests. Mann (24) first suggested using the nonparametric correlation coefficient Kendall’s tau (τ) to test for trend. The Y variable is evaluated for trend, while X is time. The Mann-Kendall test can be stated most generally as a test for whether Y values tend to increase or decrease with time (T). The null hypothesis (H0) and alternative hypothesis (H1) are stated as
H0: prob[Yj > Yi] ) 0.5,
where time Tj > Ti
H1: prob[Yj > Yi] * 0.5
(two-sided test)
No assumption of normality is required, but there must be no serial correlation in the Y data for the resulting p values to be correct. The test possesses the useful property of other nonparametric tests in that it is invariant to (monotonic) power transformations such as logarithms and square roots. Test results on the logarithms of Y will be identical to those for Y itself. Tau is most easily computed by first ordering all data pairs by increasing T. If a positive correlation exists, the Y values will increase more often than decrease as T increases. For a negative correlation, the Y values will decrease more often than increase. If no correlation exists, the Y values will increase and decrease about the same number of times. The test statistic (Kendall’s S) is calculated by subtracting the number of “discordant pairs” M, the number of (T,Y) pairs where Y decreases as T increases, from the number of “concordant pairs” P, the number of (T,Y) pairs where Y increases with increasing T:
S)P-M where
T j > Ti P ) “number of pluses”, the number of concordant pairs where Yj > Yi M ) “number of minuses”, the number of discordant pairs where Yj < Yi for all i ) 1, ..., (n - 1) and j ) (i + 1), ..., n 10.1021/es051650b Not subject to U.S. copyright. Publ. 2006 Am. Chem.Soc. Published on Web 06/06/2006
If there is no trend, so that the null hypothesis is true, about half of the comparisons would be expected to be concordant, and half discordant. Kendall’s S would then be expected to be close to 0. If S is significantly different from zero, the data indicate that a trend in Y occurs. Note that there are n(n - 1)/2 possible comparisons to be made among the n data pairs. If all Y values increased along with the T values, S ) n(n - 1)/2. In this situation the correlation coefficient τ should equal +1. When all Y values decrease with increasing T, S ) -n(n - 1)/2 and τ should equal -1. Therefore, dividing S by n(n - 1)/2 will give a value always falling between -1 and +1. This then is the definition of Kendall’s τ correlation coefficient, measuring the strength of the monotonic association between Y and time:
τ)
S n(n - 1)/2
To test for significance of τ, S is compared to what would be expected when the null hypothesis H0 is true. The p value summarizes the probability of getting the observed value of τ, or one more extreme, when the null hypothesis is true. When p is small, the likelihood that there is no trend is also small, and H0 is rejected. For n e 10, exact p values for the significance of τ are found in table B8 of ref 25 and elsewhere. For larger sample sizes S is converted to Z, a statistic that is approximated by a normal distribution. The formula for Z is given below. This approximation is used by most statistical software.
{
S-1 σS Ztau ) 0 S+1 σS and
if S > 0 if S ) 0 if S < 0
σS ) x(n/18)(n - 1)(2n + 5) The Seasonal Kendall test (1) accounts for seasonality by computing the Mann-Kendall test on each of m seasons separately, and then combining the results. So for monthly “seasons”, January data are compared only with January, February only with February, etc. No comparisons are made across season boundaries. The Seasonal Kendall statistic Sk is simply the sum of individual Kendall’s S statistics SL for each of the L ) 1 to m seasons. m
Sk )
∑S
L
L)1
When the product of number of seasons and number of years is more than about 25, the distribution of Sk can be approximated quite well by a normal distribution (1). Sk is standardized by subtracting off its expected value of 0 and dividing by its standard deviation σSk. A continuity correction of +1 or -1 is added to make the statistic more closely conform to a smooth curve. The resulting ZSk is evaluated against a table of the standard normal distribution.
{
Sk - 1 σS k
Z Sk ) 0 Sk + 1 σS k where
if Sk > 0 if Sk ) 0 if Sk < 0
µ Sk ) 0
x∑ m
σSk )
(nL/18)(nL - 1)(2nL + 5)
L)1
and
nL ) number of data in the Lth season When some of the Y and/or T values are tied, the formula for σSk must be modified (25). The overall Seasonal Kendall trend slope for Y over time T is computed as the median of all slopes between data pairs within the same season. No cross-season slopes contribute to the overall estimate of the Seasonal Kendall trend slope. This slope is the median rate of change of Y over time. Serial correlation among data within a season may cause the test statistic, which assumes independence among data, to be inaccurate. Hirsch and Slack (26) dealt with this by modifying the test statistic, accounting for the covariance among data points. They recommend that the adjustment for serial correlation be used when there are 10 or more observations within a season, resulting from simulations of monthly data (12 seasons). For fewer than 120 observations (10 values in sequence within a season), it is difficult to accurately estimate the covariance matrix needed for the adjustment, and they recommend using the original, unadjusted version of the test. The Seasonal Kendall test will find evidence of a trend over time when a trend of consistent direction occurs in many of the individual seasons. In contrast, assuming there are sufficient data to detect a trend, a finding of “no trend” may be caused by one of two patterns. First, there may be no trend in most or all of the seasons. The null hypothesis is, in fact, true. Second, trends within seasons may be present but occurring in opposite directions, canceling one another out. This heterogeneity of trend patterns results in a finding of “no consistent, overall trend” because values are increasing in some seasons while decreasing in others. This is the appropriate outcome of the Seasonal Kendall test, but trend patterns within individual seasons and their heterogeneity in direction may still be of interest. van Belle and Hughes (27) devised a test for trend heterogeneity. In their test, the Seasonal Kendall Z statistics for each season are squared and summed, forming a chi-square (χ2) statistic measuring the overall change in data regardless of direction. This statistic is partitioned into a chi-square due to heterogeneity and one due to trend. Each chi-square statistic can be separately tested to determine significance. The test of the trend chisquare is strongly related to the Seasonal Kendall results. If the test for heterogeneity is significant, van Belle and Hughes (27) state that the test for trend is “not appropriate”. Counterbalancing trends will likely result in an insignificant Seasonal Kendall test in this case, reflecting that there is no consistent trend in one direction throughout the time period. In other words, it would be “not appropriate” to state that there is “no trend” in the data. Trends exist, but in differing directions in different seasons of the year. However, it would be appropriate to state the results of the nonsignificant Seasonal Kendall test as “there is no consistent trend in one direction across all seasons”. Persons interested in determining whether trends exist in differing directions among the seasons are advised to perform the van Belle and Hughes test alongside the Seasonal Kendall test. Gilbert (4) provides a detailed example of the computation of the van Belle and Hughes test. Computation of the Regional Kendall Test. The Regional Kendall test substitutes location for season and computes the equivalent of the Seasonal Kendall test. Kendall’s S is VOL. 40, NO. 13, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
4067
computed for each of m locations, and the sum Sr of the individual statistics provides the overall test statistic for the regional test. m
Sr )
∑S
L
L)1
The same estimate of standard deviation as provided by (1) is used to compute the statistic and p value of the normal approximation for the overall Regional Kendall test.
{
Sr - 1 σS r
if Sr > 0 if Sr ) 0
Zr ) 0 Sr + 1 σS r where
if Sr < 0
µSr ) 0
σSr )
x
m
∑(n /18)(n L
L
- 1)(2nL + 5)
L)1
and
nL ) number of data at the Lth of m locations Like the Seasonal Kendall test, the Regional Kendall test is an “intrablock” test (27). Test statistics are computed on each block of data separately, and the overall test combines the individual test statistics so that no cross-block comparisons are made. For the Regional Kendall test, the blocking factor is location. If some locations exhibit an upward trend while others exhibit a downward trend, their Sr statistics will cancel out, and no consistent trend in the same direction across the locations will be found. The Regional Kendall test looks for consistency in the direction of trend at each location, and tests whether there is evidence for a general trend in a consistent direction throughout the region. Patterns at an individual location occurring in the same direction as the regional trend provide some evidence toward a significant regional trend, even if there is insufficient evidence of trend for that one location. Either the original Seasonal Kendall test (1) or the version modified (26) to correct for serial dependence (serial correlation) in the data can be used as the basis for the Regional Kendall test for trend. Groundwater concentrations may exhibit serial correlation between observations over time at a well, for example, especially if sampling frequencies are short in comparison to flow rates. However, following the guidelines given in ref 26, use of the modified test is not warranted unless there are at least 10 observations per location. One of the advantages of the Regional Kendall test is its “blocking” mechanism. Because individual Mann-Kendall tests are conducted at each location, all characteristics of that location are held constant, or “blocked out” of the test. Groundwater nitrate concentrations often decrease with well depth, for example. Methods that combine data from several wells must either take differing depths into account, or suffer from increased noise due to differences in depths among wells. Procedures for removing this variation include subtracting out the relation of nitrate with depth by computing the residuals from a LOWESS smooth (25). The subsequent trend test is then computed using the residuals as the Y variable. An advantage of the Regional Kendall test is that 4068
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 13, 2006
this extra step is not necessary, as depth does not change for data within a single well, and no cross-well comparisons are made. van Belle and Hughes (27) proposed chi-square statistics to test for homogeneity of trends among different locations, as well as among different seasons. Individual test statistics for combinations of seasons and locations are summed in a variety of ways in order to compute the tests, somewhat analogous to an analysis of variance procedure. Chapter 17 of ref 4 provides an example of using these tests to determine differences in the directions of trends among both seasons and locations, calling it “testing for global trends”. The basic requirement for performing this set of tests is sufficient numbers of data to compute a Kendall’s S statistic for each combination of season and location. For m locations and k seasons, a minimum sample size might be 5mk measured values. The Regional Kendall test has a somewhat simpler objective than the van Belle and Hughes set of tests. The Regional Kendall test will find evidence of a consistent spatial trend when a trend in the same direction occurs at many of the individual locations. Failing to reject the null hypothesis of “no regional trend” may be because there are no trends at most locations, or because observable trends have opposite directions at different locations, canceling one another out. The Regional Kendall test does not detect this heterogeneity of trend. Heterogeneity results in a finding of “no consistent, overall trend” by the test. Detecting heterogeneity of trend patterns is the objective of the van Belle and Hughes test. Persons interested in detecting whether trends exist in differing directions among the locations are advised to perform the van Belle and Hughes test alongside the Regional Kendall test. What Constitutes a “Region”? A region is a spatial area that can be defined by physical characteristics other than the trend results themselves. Although usually a contiguous area over space, it could be composed of a collection of disjoint areas all with the same physical characteristics, for example all areas with the same soil series. As with all sampling efforts, sampling for regional trend should be performed in a manner that ensures conditions over time and space within the region are fully covereds“representative sampling”. Sampling locations should not be clustered in only one portion of the region. To do so would incur the danger that differing trend patterns (including no trend) occur in the unsampled portions, and results from the cluster are being falsely projected throughout the entire region. Also, very closely spaced locations may be repeatedly sampling the same phenomena, i.e., spatial autocorrelation. Measurements are more like replicates in this case than independent measures of system behavior. If this occurs, repeated measures of trend (or no trend) will again falsely project into the entire region the results from a much smaller area. Gilbert (4) provides a good discussion of spatial sampling design. In the next section we provide an example of defining our region and avoiding spatial clustering. Application of the Regional Kendall Test. The Regional Kendall test was used to test for trends in groundwater nitrate concentrations in 474 wells distributed across three counties in the Columbia Basin GWMA (Figure 1). These wells were selected using a stratified random site selection (SRSS) technique (28) to ensure an equal-area random selection of wells across depth categories and areal extent of the GWMA. The resulting set of wells should be spatially unbiased and broadly representative of conditions throughout the GWMA. Use of SRSS avoids clustering of wells in small locations that would be spatially repetitive and therefore not independent. To additionally evaluate spatial independence, distances between well pairs were compared to distances that groundwaters would be able to travel over the timeframe of the
FIGURE 1. Study area for the GWMA. study. Median lateral hydraulic conductivities of 4.98 × 10-4 ft/s and 1.6 × 10-5 ft/s in the overburden and basalt aquifers, respectively (29), are small enough that among all possible pairs of the 474 wells only three well pairs are proximate enough to possibly have sampled the same water over the 4 year study. Spatial dependence is therefore unlikely to have affected our results. No significant trend was found considering all 474 wells together, as the Regional Kendall p value was 0.09. The median rate of change of -0.002 mg/L per year was not significantly different from zero. Wells within each of three counties, Franklin, Adams, and Grant, were then tested separately to detect whether there was a consistent concentration trend within each county. As seen in Figure 2, nitrate concentrations for the three counties generally decrease with depth, although concentrations for Franklin and Grant counties increase slightly with depth up to about 200 feet. The Regional Kendall test blocked out these differences due to depth by performing the Mann-Kendall test separately at each well. Still no significant trend in concentration was detected for any of the three counties: for more detail, see ref 22. There is no consistent regional trend in nitrate concentrations for all wells within the GWMA. Within the GWMA, 99 wells had a starting nitrate concentration in 1998 exceeding the maximum contaminant limit (MCL) for drinking water of 10 mg/L nitrate (see Figure 3). Logistic regression (not shown) indicated an increased
FIGURE 2. Nitrate concentrations for Adams, Franklin, and Grant counties. Nitrate concentrations generally decrease with depth. The lines drawn are LOWESS smooths (25), showing the general pattern of how nitrate changes with depth. probability of nitrate exceeding 10 mg/L in wells that had short well casings in areas of high fertilizer application on soils with high infiltration rates Therefore this group of wells is especially important as they are the most sensitive to changes in surface practices, and can be expected to most quickly show the effects of remediation. To test whether this group of 99 wells (group A) had a different trend pattern VOL. 40, NO. 13, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
4069
FIGURE 3. Wells within the GWMA exceeding 10 mg/L in 1998, shown as black dots. than the remaining 375 wells (group B), the test for a contrast between two subsets (27) was computed using the test statistic
(SA - SB)2 var(SA) + var(SB) where S is the sum of the Mann-Kendall statistics for all wells in that group, and var(S) is the corresponding variance normally placed in the denominator of the group’s Regional Kendall test. The test statistic is compared to a chi-square distribution with one degree of freedom. For these higher and lower nitrate groups, the test statistic is
(-89 - 18)2 ) 6.85, 361 + 1310
p ) 0.009
The null hypothesis of homogeneity is rejected, and therefore the wells with nitrate in 1998 above 10 mg/L have a different trend pattern than the wells with nitrate below 10 mg/L. The trend at these 99 high-nitrate wells was then investigated using the Regional Kendall test. The Regional Kendall test found a statistically significant downward trend in concentration from 1998 to 2002 (p value < 0.0001) for this group of 99 wells. The median rate of change was -0.4 mg/L per year. To illustrate the calculations for the Regional Kendall test, we performed the test using the 17 wells in Adams County with nitrate in 1998 exceeding 10 mg/L. A Mann-Kendall test was conducted at each well, each having three measured nitrate concentrations (in 1998, 2000, and 2002). For the n ) 3 observations, there were n(n - 1)/2 ) 3 possible comparisons at each well: 1998 to 2000, 4070
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 40, NO. 13, 2006
1998 to 2002, and 2000 to 2002. The maximum value that S could take was therefore +3, and the minimum was -3, at each well. If the null hypothesis were true, about half of the wells would be expected to have positive S values, and half would have negative S values. Rows in Table 1 present the results of the Mann-Kendall test for each of the 17 wells. With only three observations per well, none of the individual tests at a well were significant, even those where SL was +3 or -3, the most extreme values possible. However, evidence for trend accumulates as results from individual wells are combined. Out of the 17 tests, 14 have a negative slope. Only 8.5 would be expected to do so if the null hypothesis were true. For the strongest signals (where all comparisons are uptrends or all downtrends), there are several more -3 values than +3 values, so the strongest signals are more likely downtrends than uptrends. The Regional Kendall test incorporates this information when conducting the overall test for trend. Results for Adams County are shown in the last row in Table 1. Sr for the regional test equals -17, with a p value of 0.043. This indicates a significant downtrend, with an estimated slope of -0.35 mg/L of nitrate per year for the subset of wells within this county. Franklin County similarly showed a significant downward trend at high-nitrate wells of -0.46 mg/L per year (p value < 0.0001). Wells in Grant County did not show a statistically significant trend (p value ) 0.17). Differences for the three counties are pictured in Figure 4. Running the Regional Kendall test at all 99 wells together produced the significant test results cited earlier (p value < 0.0001, with a median rate of change of -0.4 mg/L per year) for this subset of wells within the GWMA.
FIGURE 4. Box plots of nitrate concentrations for wells exceeding 10 mg/L in 1998. Trend for the three sampling years is shown as a drop in the median (central line within the box) concentration across the three years. Trends for Adams and Franklin counties are significant, as measured by the Regional Kendall test. Trend for Grant County is not significant. Data from ref 22. Outliers as defined by the standard definition of a box plot (25) are plotted as individual asterisks.
TABLE 1. Mann-Kendall and Regional Kendall Test Results, Adams County, for Wells with Nitrate Concentrations Exceeding 10 mg/L in 1998 (22) Mann-Kendall Test Results ZL
well
N
SL
τ
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
1 -1 -1 -1 -3 -1 -3 -3 -1 3 -3 -1 1 3 -3 -3 -1
0.333 -0.333 -0.333 -0.333 -1.000 -0.333 -1.000 -1.000 -0.333 1.000 -1.000 -0.333 0.333 1.000 -1.000 -1.000 -0.333
N
Sr
τ
51
-17
-0.333
0.000 0.000 0.000 0.000 -1.044 0.000 -1.044 -1.044 0.000 1.044 -1.044 0.000 0.000 -1.044 -1.044 -1.044 0.000
Regional Kendall Test Results Zr
Discussion The power of a trend test is the ability to detect a trend when it is present. Several papers have evaluated the power of the Seasonal Kendall (SK) test for trend. Those results should apply directly to the Regional Kendall (RK) version of the test. Hirsch and Slack (26) presented their correction for serial correlation to the SK test by comparing its power to a previous test by Dietz and Killeen (30), demonstrating the greater power of the new test. van Belle and Hughes (27) compared the power of the SK test to an aligned-rank test by Farrell (31), finding the SK test to be generally about 96% as powerful. However, they noted that the SK test is much more amenable to missing data and nondetects, so that its greater utility outweighed the small difference in power. Taylor and Loftis (32) found the power of the Seasonal Kendall test to be greater than ANCOVA and regression alternatives for realistic condi-
-2.027
p value
intercept
slope
1.000 1.000 1.000 1.000 0.296 1.000 0.296 0.296 1.000 0.296 0.296 1.000 1.000 0.296 0.296 0.296 1.000
-76.1 3272.4 229.8 7336.7 15230.6 1023.7 416.6 1115.1 5567.4 -487.1 4896.2 8984.0 -236.9 -3632.2 2327.9 5636.5 4823.7
0.0500 -1.6250 -0.1100 -3.6650 -7.6000 -0.5000 -0.2000 -0.5525 -2.7800 0.2500 -2.4475 -4.4920 0.1250 1.8250 -1.1600 -2.8125 -2.4025
p value
intercept
slope
0.043
713.1
-0.35
tions of environmental data, particularly differences in variance between seasons, both with and without mild serial correlation. The SK test also had better power characteristics than Farrell’s test when the variances among seasons differed. Thas et al. (33) compared the power of four versions of trend tests, including the SK test with and without correction for serial correlation, using data that contained serial correlation. They demonstrated the value of the correction, showing that the uncorrected test rejected the null hypothesis too often (is “liberal”), mistaking serial correlation for trend. Their findings support why Hirsch and Slack (26) developed the corrected version. One clear advantage of the SK/RK test over regression methods is that it can be directly applied to data containing “nondetect” values (34). The test itself requires only a determination of an increase or decrease between observaVOL. 40, NO. 13, 2006 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
4071
tions, not the magnitude of that change, in order to perform the test. Therefore, a
