Biogeosciences Linear trends in seasonal vegetation time series and the modifiable temporal unit problem

Time series of vegetation indices (VI) derived from satellite imagery provide a consistent monitoring system for terrestrial plant productivity. They enable detection and quantification of gradual changes within the time frame covered, which are of crucial importance in global change studies, for example. However, VI time series typically contain a strong seasonal signal which complicates change detection. Commonly, trends are quantified using linear regression methods, while the effect of serial autocorrelation is remediated by temporal aggregation over bins having a fixed width. Aggregating the data in this way produces temporal units which are modifiable. Analogous to the wellknown Modifiable Area Unit Problem (MAUP), the way in which these temporal units are defined may influence the fitted model parameters and therefore the amount of change detected. This paper illustrates the effect of this Modifiable Temporal Unit Problem (MTUP) on a synthetic data set and a real VI data set. Large variation in detected changes was found for aggregation over bins that mismatched full lengths of vegetative cycles, which demonstrates that aperiodicity in the data may influence model results. Using 26 yr of VI data and aggregation over full-length periods, deviations in VI gains of less than 1 % were found for annual periods (with respect to seasonally adjusted data), while deviations increased up to 24 % for aggregation windows of 5 yr. This demonstrates that temporal aggregation needs to be carried out with care in order to avoid spurious model results.


Introduction
Vegetation systems provide a quick and measurable response to many environmental changes at a wide range of spatial and temporal scales.The availability of historical time series from satellite observations with daily global coverage makes operational monitoring of vegetation condition a matter of detecting and interpreting changes within these datasets.Change detection, however, is often complicated by a number of statistical preconditions that are intrinsic to time series of spectral vegetation indices with dense sampling intervals.Literature is replete with the most frequently used approach for detecting temporal trends, i.e. fitting linear regressions of a (temporally aggregated) vegetation index (VI) against time (Bai et al., 2008;Herrmann et al., 2005Herrmann et al., , 2007;;Olsson et al., 2005;Paruelo et al., 2004), but this needs to be done with care in order to avoid spurious trends.The detected slope (or gain) coefficient can be used to calculate the amount of change, but it is not always tested for significant deviation from zero, nor are standard statistical assumptions always respected (de Beurs and Henebry, 2005).Seasonal variation is an important cause for the data to violate assumptions like homogeneous variation and absence of serial correlation in the residuals.In few cases linear models were fitted directly to seasonal data (e.g.Pelkey et al., 2000), but seasonality is typically remediated using temporal aggregation, where the aggregation window (or bin size) corresponds to the length of a calendar year.The resulting bins can be regarded as temporal units, which, like spatial units, are modifiable (Taylor, 2010).In case of spatial units, it has been demonstrated that the size may influence the model results, which is known as the Modifiable Area Unit Problem (MAUP) (Openshaw and Taylor, 1979).This problem may affect a myriad of spatial studies in geography (Dark and Bram, 2007) and remote sensing (Marceau et al., 1994).Similarly, there is a Modifiable Temporal Unit Problem (MTUP) that is as troublesome as the MAUP.This is essentially a question of scale in the temporal dimension (C ¸öltekin et al., 2011).In analysis of time series of satellite vegetation indices this problem is easily disregarded, although it may result in misjudgements of temporal trends in the data.Aspects of the problem include the starting phase of a time series or segment, its extent and the level of temporal aggregation.Some studies which analyse temporal trends in vegetation activity have been debated in literature, in part because of such issues (Samanta et al., 2011;Wessels, 2009).The aim of this paper is to demonstrate possible MTUP effects in analysis of time series of satellite imagery using both real and simulated VI data and to provide, in this sense, a framework for linear time series regression.

Trend analysis and the modifiable temporal unit problem
Many trend analysis methods exist, including parametric and non-parametric approaches.The most common method to detect changes in cyclic time series is the use of a linear model (Eq. 1) obtained from ordinary least-squares (OLS) regression.The slope coefficient β, or gain, was used to calculate the change in Y as β times the number of bins.This number is determined by the aggregation level (or: number of observations per bin), which is equivalent to the sample interval.Given that the datasets consist of 24 observations per year, aggregation level 24 corresponds to yearly bins and so on.
Where ω t , the residual, is ideally independent and identically distributed (iid), i.e. white noise.The dependent variable Y can be any kind of VI or cyclical environmental parameter in general.The most common spectral vegetation indices are based on the rapid change in reflectance of chlorophyll between the red and near infrared (NIR) ranges.Here, we used the Normalized Difference Vegetation Index (NDVI), which is a commonly used proxy for terrestrial photosynthetic activity.A decrease over time is referred to as browning, whereas an increase indicates greening (de Jong et al., 2011).
Given this, we used a three-step approach to demonstrate the MTUP effect: 1.The influence of starting phase and data extent is illustrated using a perfectly harmonic model, without trend or noise components.Using a sample size of 24 observations per cycle this implies that the model residuals are far from independent, which invalidates linear regression by OLS.However, it provides a theoretical scenario to demonstrate our point that spurious slopes can be detected from cyclical data.
2. Next, the sample was aggregated into bins of fixed size.These bins represent different temporal units (or aggregation levels).The MTUP effect is demonstrated by calculating the change in NDVI from linear models (Eq. 1) fitted for different bin sizes.The minimum number of bins over the full length of the synthetic time series was set to 5.
3. Finally, step 2 was repeated using real data from Advanced Very High Resolution Radiometer (AVHRR) sensors.The detected changes were compared to the corresponding change obtained using seasonally adjusted data without temporal aggregation.The seasonal adjustment was carried out using a Fourier method (Roerink et al., 2000) with four components, following de Jong et al. (2011).The significance of the slopes in the seasonally adjusted data was assessed using generalized least-squares (GLS) in order to account for remaining short-lag serial correlation.A sample of 1000 pixels was used for calculation of deviations introduced by the MTUP and the state of Queensland in north-eastern Australia was used to illustrate possible spatial patterns introduced by different temporal aggregation schemes.
Provided that the level of serial autocorrelation can be disregarded after aggregation, the significance of the detected trends can be tested using OLS-based t-tests conform the hypotheses H0: β = 0 and HA: β = 0. Analogous to step 3 (described above), we used GLS-based tests to account for remaining serial autocorrelation.Slopes coefficients beyond the 0.05 confidence level were considered significant.All analyses were performed using standard R functionality (R Development Core Team, 2011).

Synthetic data
Synthetic time series were used to illustrate the effect of cyclic data on regression analysis.For this purpose, model parameters were chosen in such a way that they approximate the AVHRR time series described below for a temperate (non-forest) environment with a single growing season.
As such, the peak-to-peak amplitude (2A) was set to 0.6with a mean of 0.4 -and the number of observations per year to 24.NDVI (Y ) was simulated using a cosine model with no underlying positive or negative trend (Eq.2), in which a Biogeosciences, 9, 71-77, 2012 www.biogeosciences.net/9/71/2012/denotes the amplitude, time (t) the radians equivalent of the observation number (x) (Eq. 3) and the phase shift.

AVHRR data
The longest run of NDVI measurements is available from AVHRR sensors on board a series of National Oceanic and Atmospheric Administration (NOAA) satellites.Acquisition started in 1981 and is still on-going, but pre-processed datasets are for now available until the year 2006.Preprocessing includes correction for orbital decay, satellite changes and several atmospheric effects.The resulting NDVI values were aggregated into biweekly composites with ∼8 km spatial resolution.The full description of the dataset and the processing steps is provided by (Tucker et al., 2005).A sample of 1000 pixels was used for the MTUP analysis.
Other vegetation indices, including the Enhanced Vegetation Index (EVI) or the Soil Adjusted Vegetation Index (SAVI), have similar statistical characteristics and therefore the problem described here is not restricted to NDVI.

Synthetic data
In case of seasonal data with a dense sampling interval, both starting phase and extent influence the linear model.Figure 1a illustrates the ad absurdum case that a fit on a perfectly harmonic model without linear trend results in a range of slope coefficients -positive or negative, varying with the starting phase -but never zero.Zero slopes are obtained only if both sides of a minimum or a maximum are equally sampled.This might, however, result in over-or underestimation of the mean NDVI (Y ) from the model intercept (α).The slope coefficient is linearly related to the amplitude used for the seasonal model and is inversely related to the extent of the time series (or segments).The latter is illustrated in Fig. 1b, which shows 4 different extents (respectively 2, 4, 6 and 8 yr) and the associated linear models for a given starting phase.
The slope coefficient changed with the level of temporal aggregation, which influenced the detected change in NDVI.This is illustrated in Fig. 2 using a similar synthetic data set as above but now having a length of 26 full cycles, which is comparable to that of AVHRR data sets.The variation in detected NDVI change increased with the level of aggregation, which resulted in larger uncertainty in model predictions.The true change in NDVI (i.e.zero) is only obtained by aggregation over complete cyclic periods, i.e. one or more years.Given the perfect periodicity of this synthetic example, this is true by definition, but in reality calendar years may not fit the periodicity of VI time series because of shifts in vegetation phenology and variations in growing season length.For this reason, the same analysis was also performed on real AVHRR time series (see Sect. 3.2).
As mentioned, the starting phase and extent influence each segment of the linear analysis.The effect reduces with longer extents, but it appeared that the change for the longest available run of NDVI data (AVHRR) is of the same magnitude as changes found by trend studies which account for seasonality (e.g. de Jong et al., 2011;Wang et al., 2011;Zhou et al., 2001).Table 1 lists the calculated changes in NDVI for a perfectly periodic model without trend component (Fig. 1) for time spans of AVHRR (26 yr), MODIS (11 yr) and 1 yr.The seasonal peak-to-peak amplitude was set to 1 for ease of comparison to other amplitudes.Given that the order of maximum change reported in literature is ∼4.0 × 10 −2 for AVHRR data, trends induced by phase shift might introduce errors varying from 10 % to 90 % if seasonality is not accounted for.

AVHRR data
An example of the MTUP effect for an AVHRR pixel with a significant (p < 0.05) trend component is shown in Fig. 3.The change in NDVI was found to vary among aggregation levels: slight differences with respect to seasonally adjusted data occurred for 1 or 2 yr windows and substantial differences for all other cases.This illustrates that aperiodicity in the data or in the aggregation window might result in considerable deviations from the change in NDVI found using seasonal correction instead of temporal aggregation.If temporal aggregation is used to account for seasonality, the characteristics of the commonly used harmonic functions dictate the use of whole periods as aggregation windows, in order to force the deviation in slope coefficient to zero (Fig. 2).Aiming for bins holding complete seasons, the MTUP analysis was carried out on AVHRR data for aggregation levels of exactly 1 to 5 calendar years.It appeared that the significant changes in NDVI (i.e.significant for all aggregation levels) varied considerably with respect to trend analysis with seasonal adjustment instead of temporal aggregation.Table 2 lists the mean deviation and mean absolute deviation for the AVHRR sample and for the example in Fig. 3.This table shows that the mean deviation and the variation for the sample are low (around 1 %) for fine aggregation levels but considerably higher in case of coarser aggregation: a mean absolute deviation of 24 per cent was found for 5 yr aggregation.In few cases the detected change in NDVI switched between positive (greening) and negative (browning), although these trends could not be confirmed using significance tests.It was almost exclusively found that the amount of detected change increased with the level of aggregation, so greening and browning trends appeared stronger after aggregation.2 for the NDVI values corresponding to full-year aggregation levels.
panels show the change in NDVI with respect to the change detected without temporal aggregation.Spatial patterns are hardly perceivable at fine temporal aggregation levels, but they become apparent at coarser levels.Temporal aggregation will only lead to unbiased estimates of VI trends if the following requirements are met: 1.The regarded time series is perfectly periodic.Any type of aperiodicity, including phase shift or incomplete periods at the start or end, may result in incorrect model parameters.
2. The aggregation level corresponds exactly to the period of the seasonal signal (often one calendar year).Aggregating over multiple periods increases the risk of MTUP effects.
Other analysis methods may implicitly require similar conditions.For example, non-parametric trend tests such as the seasonal Mann-Kendall method (Hirsch and Slack, 1984) and related slope estimators like Theil-Sen selection (Birkes and Dodge, 1993) rely strongly on the absence of seasonal changes.In case of vegetation indices, this assumption often renders the method unsuitable for large spatial and temporal scales (de Jong et al., 2011).The results from this study indicate that similar effects may disturb linear regression.The effects, however, may not be as conspicuous because parametric models are less robust against this type of error than non-parametric models (McBride and Loftis, 1994).Seasonally adjusting the data using a decomposition method (e.g Cleveland et al., 1990) provides another approach for eliminating serial autocorrelation and the MTUP.
In that case, the seasonal model used should be appropriate for the growing regime and ideally should take possible seasonal changes into account.An example of such an approach is provided by Verbesselt et al. ( 2010).Yet other methods to describe seasonal time series include data generating processes without the use of deterministic functions.Changes in NDVI are likely to persist from one period to the next.As such, it may be reasonable to represent the process using an autoregressive (AR) model (Eq.4).
Where p represents the autoregressive order (or maximal lag), α i are stochastic model parameters and ω t is white noise (Cowpertwait and Metcalfe, 2009).Under these conditions, fitting a temporal trend to y t (e.g. using Eq. 1) will generate misleading results termed spurious regressions (Granger and Newbold, 1974).The spuriousness resides in the effect of autocorrelated residuals which bias the test towards rejection of the null hypothesis, even when the series are generated as statistically independent random walks (Phillips, 1986).This problem is analogous to the MTUP effects described in Sect.3.1 and will be subjected to the effects described in Sect.3.2, were y t temporally aggregated before fitting the temporal trend.
Using temporal aggregation, slopes may be found significant at a given aggregation level while not at another level.False positives (trend is found to be significant while it is not in reality) are likely to occur more frequently at low aggregation levels, but examples of the opposite case were found as well (not shown).If the significance of trends is not considered in the analysis, the MTUP may not only affect the amount but also the sign of detected changes in NDVI.These effects can be expected to be largest in regions with high seasonal amplitudes.MTUP as described here cannot affect analysis of time series without seasonal amplitude, e.g. in dense tropical forests or very sparsely vegetated areas.On the other hand, the use of vegetation indices is such regions www.biogeosciences.net/9/71/2012/Biogeosciences, 9, 71-77, 2012  is disputed; either because of signal saturation issues (Huete et al., 1997) or because of over-estimation of NDVI or related parameters (Fensholt et al., 2004).

Conclusions
Ordinary least squares (OLS) time series regression can be used to quantify trends in cyclic data but temporal aggregation needs to be carried out carefully in order to avoid spurious results.The risk of artefacts is minimal at an aggregation level corresponding to a full period, for instance a calendar year.Coarser aggregation levels tend to overestimate the amount of change and result in higher variation in model predictions, especially from 3 periods onwards.However, the use of a full-period window may be impractical because VI time series are hardly ever free of changes in seasonality.Aperiodicity within long-term time series of vegetation indices is intrinsic to certain land cover types and may arise from variations in start and length of growing seasons as a result of variations in temperature and/or precipitation.The starting phase and the choice of aggregation level -or temporal unit -affect the estimation of model parameters, including the slope coefficient or gain.In this study, the amount of absolute change that was attributed to the modifiable temporal unit problem (MTUP) varied between 1 % and 24 % for full-period aggregation levels.

Fig. 1 .
Fig. 1.The effect of phase shift (a) and extent (b) on linear regression of NDVI against time for seasonal data.The legend in (a) shows the phase shift in months.In (b), the colors each represent an additional extent of 2 yr.Accordingly, magenta, green, blue and red refer to extents of 2, 4, 6 and 8 yr respectively.

Fig. 2 .
Fig. 2.Change in NDVI detected using linear models depending on the temporal aggregation level.The plot illustrates the MTUP effect for a harmonic model of 26 yr (length of common AVHRR datasets) and an amplitude of 0.6 NDVI units.The points shown are the ones that resulted in a change in the number of bins, while increasing the aggregation level.The vertical dotted lines indicate aggregation levels corresponding to 1, 2, 3, 4 and 5 yr respectively and the grey line shows a LOESS curve indicating an average trend.

Fig. 3 .
Fig.3.Change in NDVI detected using linear models depending on the temporal aggregation level.The plot illustrates the MTUP effect for an AVHRR pixel with significant (p < 0.05) trend component (location: 46.21N/110.65E,Montana, USA).See the caption of Fig.2for additional information and Table2for the NDVI values corresponding to full-year aggregation levels.

Fig. 4 .
Fig. 4. Detected changes in NDVI for different levels of temporal aggregation -example for Queensland (Australia).The top-left map shows the location of Queensland and the 5 panels show the change in NDVI with respect to changes found using seasonal adjustedment without temporal aggregation.The aggregation levels correspond to 1, 2, 3, 4 and 5 yr (YR) respectively.

Table 1 .
Detected changes in NDVI using linear regression on synthetic seasonal data for different starting phases and extents.Change values were calculated as the slope coefficient β (Eq. 1) times the length of the time series and were multiplied by 100 for numerical convenience.The listed lengths of 26 yr and 11 yr correspond to the approximate length of AVHRR and MODIS time series respectively.The applied harmonic model had peak-to-peak amplitude one and therefore change values can be multiplied by the actual amplitude.Due to serial autocorrelation, t-values could not be calculated.

Table 2 .
The MTUP effect for a sample of 1000 AVHHR pixels and aggeregation levels (AL) of 1 to 5 yr.Detected changes in NDVI (dNDVIagg) were compared to those found using seasonally adjusted data (dNDVIref) and listed as dNDVIagg-dNDVIref (dNDVI), dND-VIagg/dNDVIref (pct) and square root of the sample variation of pct (sd).The last columns list the detected changes in NDVI for the pixel in Fig.3.In all cases, slopes were significant (p < 0.05).NDVI values were multiplied by 100 for numerical convenience.