# Applying multiple testing procedures to detect change in East African vegetation

###### Abstract

The study of vegetation fluctuations gives valuable information toward effective land use and development. We consider this problem for the East African region based on the Normalized Difference Vegetation Index (NDVI) series from satellite remote sensing data collected between 1982 and 2006 over 8-kilometer grid points. We detect areas with significant increasing or decreasing monotonic vegetation changes using a multiple testing procedure controlling the mixed directional false discovery rate (mdFDR). Specifically, we use a three-stage directional Benjamini–Hochberg (BH) procedure with proven mdFDR control under independence and a suitable adaptive version of it. The performance of these procedures is studied through simulations before applying them to the vegetation data. Our analysis shows increasing vegetation in the Northern hemisphere as well as coastal Tanzania and generally decreasing Southern hemisphere vegetation trends, which are consistent with historical evidence.

10.1214/13-AOAS686 \volume8 \issue1 2014 \firstpage286 \lastpage308 \newproclaimproProcedure

Detecting change in East African vegetation

A]\fnmsNicolle \snmClements\corref\thanksreft1,m1label=e1]nclement@sju.edu, B]\fnmsSanat K. \snmSarkar\thanksreft2,m2label=e2]sanat@temple.edu, B]\fnmsZhigen \snmZhao\thanksreft3,m2label=e3]zhaozhg@temple.edu and C]\fnmsDong-Yun \snmKim\thanksrefm3,m4label=e4]kimd10@nhlbi.nih.gov \thankstextt1Supported by NSF Grant DMS-10-06344. \thankstextt2Supported by NSF Grants DMS-10-06344 and DMS-12-08735. \thankstextt3Supported by NSF Grant DMS-12-08735.

False discovery rate \kwddirectional false discovery rate \kwdNDVI \kwdEast Africa vegetation.

## 1 Introduction

The need to understand the Earth’s ecology and land cover is becoming increasingly important as the impacts of climate change start to affect animal and plant life, which ultimately affect human life. Knowledge of current vegetation trends and the ability to make accurate predictions is essential to minimize times of food scarcity in underdeveloped countries. Vegetation trends are also closely related to sustainability issues, such as management of conservation areas and wildlife habitats, precipitation and drought monitoring, improving land usage for livestock, and finding optimum agriculture seeding and harvest dates for crops.

The United Nations has given attention recently to precipitation and vegetation monitoring in East Africa, where a severe drought hit the entire region in mid-2011. The drought has caused a food crisis across Somalia, Ethiopia and Kenya, threatening the livelihood of over 10 million people. In many areas, the precipitation rate during the “long” rainy season from April to June 2011 was less than 30% of the average of 1995–2010. The lack of rain led to vegetation decline, crop failure and widespread loss of livestock, as high as 40%–60% in some areas [OCHA (2011)].

Droughts are commonly thought to occur from prolonged periods with less than average precipitation, which is then followed by a decline in vegetation growth. However, droughts can also arise independently from precipitation changes when soil conditions and erosion are triggered by poorly planned agricultural endeavors. Overfarming, excessive irrigation and deforestation can all adversely impact the ability of the land to capture and hold water. Thus, current efforts are geared toward monitoring agriculture and vegetation changes in hopes of minimizing the effects of low precipitation rates in future years.

Assessment of changes in a region’s vegetation structure is challenging, especially in topographically diverse areas, like East Africa. Forecasting future vegetation and agricultural planning become particularly difficult when unknown trends are occurring. However, the regions with vegetation changes are often the areas of most interest in land use management. For example, if a previously underdeveloped region is experiencing increasing trends in vegetation growth, meaning the land is able to sustain plant growth, local farmers could utilize this area to grow crops or raise livestock in future years. On the other hand, if a region is experiencing decreases in vegetation growth, this could be an indicator of overfarming, putting crops and livestock at risk of drought.

Data collection on vegetation and land cover are typically done through satellite remote sensing. The remote sensing imagery is used to convert the observed elements (i.e., the image color, texture, tone and pattern) into numeric quantities at each pixel in the image. The image pixels correspond to a square grid of land, the size of which depends on the satellite’s resolution. One such numeric value is the normalized difference vegetation index (NDVI). The NDVI has been shown to be highly correlated with vegetation parameters such as green-leaf biomass and green-leaf area, and hence is of considerable value for vegetation monitoring [Curran (1980), Jackson, Slater and Pinter (1983)]. The NDVI standard scale ranges from 1 to 1, indicating how much live green vegetation is contained in the targeted pixel. An NDVI value close to 1 indicates more abundant vegetation. Low values of NDVI (say, 0.1 and below) correspond to scarce vegetation consisting mostly of rock, sand and dirt, for example. A range of moderate values (0.2 to 0.3) indicates short vegetations such as shrub or grassland; larger NDVI values can be found in rainforests (0.6 to 0.8). Often, negative NDVI values are consolidated to be zero since negative values indicate nonvegetation and are of little use for vegetation monitoring.

Statistical and computational methods are needed to analyze remotely sensed data, like NDVI values, to determine trends in land condition and to predict areas at risk from degradation. Methodologies that detect land cover changes need to be sensitive as well as accurate, since it can be costly and risky to relocate human populations, agriculture or livestock to new regions of detected change. In such spatio-temporal data, existing change detection methodologies include geographically weighted regression [Foody (2003)], principal component analysis [Hayes and Sader (2001)] and smoothing polynomial regression [Chen, Jonsson and Tamura (2004)]. However, these methods are unable to provide an upper bound on false detections. Since there is large risk associated with falsely declaring an area to have significant vegetation changes, land use managers seek new methods that have a meaningful control over such errors.

In this article, we revisit the problem of detecting vegetation changes in East Africa based on the NDVI data and propose applying multiple testing methodologies. Such methodologies are very useful for detecting changes of statistical significance with a control over an overall measure of false detections and are being currently used in many other modern scientific investigations. We should point out that Vrieling, de Beurs and Brown (2008) did investigate this problem in a hypothesis testing framework, but, unlike ours, did not attempt to address the inherent multiplicity issue by controlling an overall false detection rate while making their final conclusions. Our proposed multiple testing methods have been developed by fine-tuning some existing ones in order to adequately capture the specific data structure and answer questions in the present context. In particular, there is a local dependency among nearby hypotheses (e.g., neighboring NDVI pixels) that should be taken into account and one should be able to identify the increasing/decreasing direction of a vegetation trend for an NDVI pixel once a significant change is detected. Our methods aim to incorporate such local dependency and control an error rate, the mixed directional false discovery rate (mdFDR), which is an overall measure of nondirectional as well as directional false detections.

We organize the paper as follows. In Section 2 we describe the East African NDVI data, its source and the associated multiple testing problem. In Section 3 we present our proposed mdFDR controlling procedures after providing some background information and notation related to multiple testing, and assessing spatial correlation. We consider dividing the East African region into subregions to adequately capture local dependencies among the NDVI values. The semivariogram plot [Cressie and Wikle (2011)], which is used to investigate the presence of spatial autocorrelation, helps determine the size of each subregion. We propose two procedures, Procedures 3.3 and 3.3, to control the mdFDR under such group or block dependence structure. Procedure 3.3 is referred to as a three-stage directional Benjamini–Hochberg (BH) procedure whose mdFDR control is theoretically shown (in an Appendix) assuming independence between but not within blocks (or subregions) and numerically examined under various dependence scenarios through simulations. Procedure 3.3 is an adaptive version of Procedure 3.3 designed to improve Procedure 3.3 through estimating the proportion of true null hypotheses within each subregion, and its mdFDR control is studied only through simulations. The findings of these simulations are reported in Section 4. In Section 5 we illustrate the applications of these proposed methods to the NDVI data collected in East Africa from 1982–2006. Discussions and concluding remarks are in Section 6.

## 2 Data description and the statistical problem

East Africa spans a wide variety of climate types and precipitation regimes which are reflected in its vegetation cover. To capture this, satellite imagery was collected over a sub-Saharan region of East Africa that includes five countries in their entirety (Kenya, Uganda, Tanzania, Burundi and Rwanda) and portions of seven countries (Somalia, Ethiopia, South Sudan, Democratic Republic of Congo, Malawi, Mozambique and Zimbabwe). This roughly “rectangular” region extends from 27.8E to 42.0E longitude and 15S to 6.2N latitude. Also included in the region are several East African Great Lakes such as Lake Victoria, Lake Malawi and Lake Tanganyika.

The remotely sensed images were recorded twice a month from 1982–2006 and then converted to NDVI values. Hence, the spatio-temporal data set consists of approximately 50,000 sites (pixels), each with 600 time series observations (24 observations per year over 25 years). The satellite’s resolution corresponds to each pixel spanning an 8 km km grid of land. This Global Inventory Modeling and Mapping Studies (GIMMS) data set is derived from imagery obtained from the Advanced Very High Resolution Radiometer (AVHRR) instrument onboard the National Oceanic and Atmospheric Administration (NOAA) satellite series 7, 9, 11, 14, 16 and 17. The NDVI values have been corrected for calibration, view geometry, volcanic aerosols, cloud coverage and other effects not related to vegetation change [Tucker et al. (2005)].

The cyclic/seasonal behavior of the NDVI time series at each pixel in the region is color-coded in Figure 1(a). For example, areas with dark green have a six-month periodicity while areas with light green have a twelve-month periodicity. Periodogram analysis indicates very strong peaks at six and twelve months, respectively, in these regions. The periodicities are reflected in bimodal and unimodal shapes in annual NDVI series, which in turn correspond to two and one rainy season each year. Figure 1(b) displays the average NDVI values for each grid point (site) over the region. Blue areas indicate regions containing only water (Indian Ocean, Lake Victoria, etc.), and thus no vegetation index was recorded. The light and dark green areas have more green vegetation on average compared to the drier areas, represented with yellow, orange and red. In this figure, one can see how this East African region spans the NDVI scale. Desert regions (with low NDVI) are within a few hundred kilometers of wetlands and rain forests (with extremely high NDVI), illustrating the large variability of climate types and precipitation regimes in this region. Figure 2 shows the time series plots for two selected pixels. The top series is a pixel selected from Southern Kenya and has a unimodal periodic pattern. The bottom series of Figure 2 is a pixel selected from the Democratic Republic of the Congo and has a bimodal periodic pattern.

We consolidated all the negative NDVI values to zero, as commonly done in vegetation monitoring, and then rescaled the remaining values by 1000. This is because negative values indicate nonvegetation areas, so they are of little use for our purpose. Prior to the analysis, we examined the data for quality assurance and eliminated a small number of pixels that were found to have several consecutive years with identical data values, which may be due to data entry errors or machine malfunction.

This data set was first examined in Vrieling, de Beurs and Brown (2008) where the interest was in studying several phenology indicators, including start of the season, length of season, time of maximum NDVI, maximum NDVI and cumulative NDVI over the season. After extracting these indicators for every year, trend tests were conducted to detect regions of significant changes in phenology indicators. The percentage of pixels with the trend test -value less than for each phenology indicator was reported separately for positive and negative slopes. The reported results indicate that much of the region has “significant” vegetation change. For example, the cumulative NDVI indicator detected 44.2% of sites with -values less than 0.10. However, this study fails to address the important statistical issue of multiplicity when making these claims about significant vegetation changes and their directions simultaneously for all the regions based on hypothesis testing.

When testing a single null hypothesis against a two-sided alternative, two types of error can occur when a directional decision is made following rejection of the null hypothesis. These are Type I error and Type III (or directional) errors. The Type I error occurs when the null hypothesis is falsely rejected, while the Type III error occurs when the null hypothesis is correctly rejected but a wrong directional decision is made about the alternative. For instance, when declaring a particular 8 km8 km grid of land as “significantly” changing in terms of vegetation, a Type I error is made if the area is not truly changing, and a Type III error is made if the area is truly changing but in the opposite direction of what is determined from the data. When such decisions are made simultaneously based on testing multiple hypotheses, as in Vrieling, de Beurs and Brown (2008), one should adjust for multiplicity and control an overall measure of Types I and III errors. Without such multiplicity adjustment, more Types I and III errors can occur than the desired level. It is particularly important to avoid these errors as much as possible in the present application. Land use managers, government and local farmers are looking to relocate East African populations of people, livestock and crops to areas of promising vegetation changes and avoid regions with decreasing changes. Since these migrations can be risky and costly, a careful consideration of the multiplicity issue seems essential when making declarations of significant vegetation changes.

In this paper, we revisit the work in Vrieling, de Beurs and Brown (2008) to adequately address the multiplicity issue. To test each 8 km8 km grid of land for vegetation change, we use the cumulative NDVI phenology indicator for each season or, equivalently, use the average NDVI per season. Since the East African region straddles the Equator, “seasons” are classified by precipitation changes rather than temperature, and it is quite probable that a particular site can have vegetation changes in one season and not another. This region receives rain in two distinct seasons, locally referred to as the “long rains” (April–June) and the “short rains” (November–December). The long rains provide more rainfall than the short rains, but generally the arrival of the short rains is more predictable. In hopes of capturing any seasonal changes, trend tests were conducted on the seasonal NDVI averages at each site for the first dry season (January–March), long rain season (April–June), second dry season (July–October) and short rain season (November–December).

To test for significant trend in each of the four seasons, we apply the monotonic trend test proposed by Brillinger (1989) for a time series consisting of a signal and stationary autocorrelated errors. We use the seasonal averages for each year as the observed time series. This test examines the null hypothesis that the series has a signal, that is, constant in time against the alternative hypothesis that the signal is monotonically increasing or decreasing in time. The test statistic is a standardized version of a linear combination of the time series values, with coefficients given in Abelson and Tukey (1963). This statistic is approximately normal with mean zero if and only if the null hypothesis is true. More specifically, given the 25-year data , , on the first dry (FD) season NDVI average in a particular site, Brillinger’s trend test can be applied for that season assuming the model.

(1) |

for , where is a deterministic signal, and is a zero mean stationary noise series, for . The test statistic is the ratio of the linear combination , with

, and the estimate of the standard error of this linear combination. The hypotheses of interest are the null and the two-sided alternative , where . This test can be similarly applied for testing the vegetation trend for the remaining three precipitation seasons. These tests were implemented by adapting R code written by Dr. Vito Muggeo, found at https://stat.ethz.ch/ pipermail/r-help/2002-December/027669.html.

Thus, for each site (8 km8 km grid of land), we have four -values, each providing an evidence of vegetation change occurring over the years in that particular season—the smaller the -value, the higher is the evidence of a significant vegetation change. Our goal is to do the following for each site: (i) combine the four seasonal -values to form a yearly -value, (ii) decide based on this yearly -value if a significant vegetation change has occurred over the years at that site, and (iii) if vegetation change is found significant, detect the season(s) that contributes to this change as well as the direction in which this change has taken place. We wish to accomplish this goal simultaneously for all sites (50,000) in the East African region in a multiple testing framework designed to ensure a control over a meaningful combined measure of statistical Types I and III errors.

## 3 The proposed multiple testing procedures

We propose two multiple testing procedures that would be suitable for applications to the vegetation data. Before that, we need to provide some background in multiple testing and a brief outline of our idea to determine the size of each subregion capturing local spatial dependencies using variograms.

### 3.1 Background in multiple testing

When simultaneously testing several null hypotheses, procedures have traditionally been developed to control the familywise error rate (FWER), which is the probability of at least one Type I error (i.e., rejecting true null hypothesis), at a desired level, say, . However, this notion of error rate is often too conservative when the number of hypotheses being tested becomes large, as in the present application. Therefore, there has been a recent surge in statistical research to define alternative, less stringent error rates and to develop multiple testing methods that control them. The false discovery rate (FDR), which is the expected proportion of Type I errors among all rejected null hypotheses, introduced by Benjamini and Hochberg (1995), is one of these alternative error rates that has received much attention.

In Benjamini and Hochberg (1995) a method was proposed, known as the BH method for short, for controlling the FDR. For testing null hypotheses , , using their respective -values , , it operates as follows: consider the ordered versions of the -values, , find , and reject the null hypotheses whose -values are less than or equal to , provided the maximum exists; otherwise, accept all null hypotheses. This procedure controls the FDR at level , under the assumption of independence or positive dependence (in a certain sense) of the -values. More specifically, the FDR of the BH method equals when the -values are independent, and is less than when the -values are positively dependent [Benjamini and Yekutieli (2001), Sarkar (2002)], where is the (true) proportion of null hypotheses. The difference between and the FDR gets larger with increasing (positive) dependence among the -values.

Often, it becomes essential for researchers, as in the present application, to determine the direction of significance, rather than significance alone, when testing multiple null hypotheses against two-sided alternatives. In other words, for each test, researchers have to decide whether or not the null hypothesis should be rejected and, if rejected, determine the direction of the alternative. Typically, this direction is determined based on the test statistic falling in the right- or left-hand side of the rejection region. Such decisions can potentially lead to one of two types of errors for each test, resulting in rejection of the null hypothesis—the Type I error if the null hypothesis is true or the directional error, also known as the Type III error, if the null hypothesis is not true but the direction of the alternative is falsely declared.

To deal with both Types I and III errors in an FDR framework, the notion of mixed directional FDR (mdFDR) has been introduced in Benjamini and Yekutieli (2005). It is defined as the expected proportion of Types I and III errors among all rejected null hypotheses. A method in Benjamini and Yekutieli (2005) was given for independent tests that controls the mdFDR when testing multiple simple hypotheses against two-sided alternatives. They proved that the original BH method controlling the FDR at can be augmented to make directional decision upon rejecting a null hypothesis according to the corresponding test statistic falling in the right- or left-hand side of the rejection region without causing the mdFDR to exceed . This so-called augmented method is referred to as the directional BH procedure.

This directional BH procedure will be extended in this paper to a situation, conforming more to the present application, where the -values are not all independent but can be grouped in such a way that they are mostly dependent within but not between the groups. Unlike in the case of multiple testing applications to genomics where gene pathways provide a natural way of grouping the -values, there are no clearly defined so-called “vegetation pathways” in the present application that we can consider for grouping the -values. Nevertheless, a statistically meaningful approach can be devised for grouping -values using variograms as outlined in the following section.

### 3.2 Variogram and its use in forming subregions capturing local spatial dependencies

The variogram is an important characteristic describing the degree of spatial dependence of a spatial random field or stochastic process . The variogram is defined as [Cressie and Wikle (2011)], and defined as if the spatial field has a constant mean. The function itself is called the semivariogram. The variogram (or the semivariogram) becomes a function of if the process is stationary, and of if it is also isotropic. For a second-order stationary process, that is, isotropic, the (theoretical) semivariogram rises from the origin to the upper asymptote Var which is called the sill of the semivariogram. The distance at which a certain fraction of the asymptote is reached is called the range of the semivariogram.

Given data , , on , the empirical semivariogram is given by , where is the set of pairs of observation at locations and such that , and is the cardinality of this set [Cressie and Wikle (2011)]. The range can be estimated by plotting the empirical semivariogram against . When extrapolated to zero distance, the empirical semivariogram reaches a nonzero value, called a nugget, caused from sampling error resulting in dissimilar values for samples at locations close to each other. See Figure 3 for illustration in the NDVI application.

In our present application, we implicitly assume that the dependency among the NDVI values is localized. In fact, as the first law of geography states in Tobler (1970), “everything is related to everything else, but near things are more related than distant things.” The range estimated from the empirical semivariogram based on the NDVI values can provide an idea of the size of neighborhoods or subregions capturing this local dependency. For instance, let be the estimated range. Then, one may consider creating subregions of a size grid box of land, where is the number of pixels greater than or equal to . The local dependency will be mostly concentrated within these subregions. Although spatial autocorrelations would still be present to some extent among the NDVI values for pixels in both sides of the boundaries, we will ignore them for the time being in order to theoretically develop our proposed procedures in the next subsection.

### 3.3 The proposed procedures

Suppose that the East African region is divided into subregions of a size grid box of land, where is the number of pixels which is determined by the range of the semivariogram plot, as described above. Each pixel, an 8 km8 km grid box of land, will be referred to as a “location.” Note that some locations in a subregion may be missing in the sense of containing only water and hence producing no NDVI observations. However, we will only consider the subregions with at least one nonmissing location. Let be the number of such subregions and be the number of locations in the th subregion.

We use a two-sided monotonic trend test for each of the four seasons in a location using the Brillinger test as described in Section 2. With being the monotonic trend parameter as defined in the Brillinger test for the th subregion, th location and th season, where , , and , let and be, respectively, the test statistic and the corresponding -value for testing the null hypothesis against its two-sided alternative. Using a Bonferroni correction over the seasons, we define the so-called combined -value for the th location in the th subregion as . The combined -value for the th subregion is defined as , which is a Bonferroni correction over the seasons and locations. With representing the null hypothesis corresponding to , we consider as the null hypothesis corresponding to th subregion, and as the null hypothesis corresponding to the th location in the th subregion.

We propose a procedure that tests the ’s against their respective two-sided alternatives and detects the directions of the alternatives for the rejected ’s. It operates in three stages, by testing , , at the first stage; , , for each such that is rejected, at the second stage; and , , for each such that is rejected, at the third stage. More specifically, our first procedure is defined as follows:

[(Three-stage directional BH)] {longlist}[Stage 3.]

Apply the BH method to test , , based on their respective -values as follows: consider the (increasingly) ordered versions of the ’s, , find , and reject the ’s for which the -values are less than or equal to , provided this maximum exists, otherwise, accept all .

For every such that is rejected at stage 1, consider testing , , based on their respective -values , , as follows: reject if .

For every such that is rejected at stage 2, first consider testing , , based on their respective -values , , as follows: reject if ; then, for each rejected , decide the direction of the monotonic trend to be the same as that of sign().

The first two stages in Procedure 3.3 identify the locations with significant vegetation changes, while the third stage allows one to make a more detailed analysis for each significant location by specifying the seasons that contribute to those changes as well as the directions in which these changes have occurred.

###### Theorem 3.1

The three-stage directional BH procedure controls the mdFDR at level if the subregions are independent.

A proof of this theorem is given in an Appendix.

We should point out that our assumption of dependence within, but not between, the subregions is made only to provide a theoretical framework for the development of our procedure in Theorem 3.1, even though, as said above, there is some dependence among the subregions. It is important to verify that this procedure can continue to control the mFDR under a certain type of positive dependence condition among the subregions, like the one that would be similar to the positive regression dependence on the subset (PRDS) condition [Benjamini and Yekutieli (2001), Sarkar (2002)] in the present context. We will do that numerically, since it seems difficult to do theoretically.

It is seen when proving the above theorem that the mdFDR of Procedure 3.3 is , where is the proportion of true ’s (out of the null hypotheses) in the th subregion. If were known, one would have used , instead of , in Procedure 3.3, to get a tighter control over the mdFDR at . In reality, when is unknown, one can consider estimating it from the data. With that in mind, we propose our next procedure as an adaptive version of Procedure 3.3 by estimating using a Storey, Taylor and Siegmund (2004) type estimate.

[(Adaptive three-stage directional BH)] Consider Procedure 3.3 with replaced by , where

(2) |

for any . Typically is chosen to be 0.5.

It is important to note how each is being adjusted in this adaptive test based on the information shared by the other -values in each subregion. If gets larger (or smaller), indicating more (or less) nonsignificant (or significant) -values in the th subregion, then moves further away from (or gets shrunk toward) zero, making it more likely to be nonsignificant (or significant) also.

## 4 Simulation studies

We ran a number of simulation studies to examine the mdFDR control property and the power of our proposed procedures. Keep in mind that the proposed procedures were developed assuming arbitrary dependence among locations within each subregion () and among the four seasons at each location (). To account for the correlation between pixels, we assume spatial stationarity across the region. We also assume isotropic spatial autocorrelation, which means that the process causing the spatial autocorrelation acts in the same way in all directions. Isotropic correlations depend only on the distance between locations and , but not on the direction.

Some frequently used isotropic covariance functions are the exponential model (), the Gaussian model () and the spherical model (), where is a scaling factor which is related to the range of the variogram. For bounded variograms that reach the sill asymptotically (e.g., exponential model), in practice, is taken to be the distance where the model reaches of the sill (also called the practical range).

In this simulation study, we selected the exponential correlation model to simulate spatial dependence between pixels within a subregion. More specifically, we did the simulation studies under the following dependence scenario.

The -values arise from test statistics , , , , where are random signals that are i.i.d. Bernoulli (). The covariance matrix, , includes correlations on three levels: between season correlation, between pixel correlation, and between subregion correlation. More specifically, let , where is the distance between location and , , and is estimated by the size of the subregion, , where , and , where , and then define . In other words, we assume , for , with for , and for .

We simulated both mdFDR and (average) power, the expected proportion of correctly rejected among all the false null hypotheses, for both Procedures 3.3 and 3.3 (with , as often considered) by choosing or ; , 100 or 400; or ; , or ; and , 0.2 or 0.5. The reason for selecting a negative value for is that vegetation trends have been noted [Vrieling, de Beurs and Brown (2008)] to change in different directions in two successive seasons.

Num. | mdFDRpower (P1) | mdFDRpower (P2) | BY | |||
---|---|---|---|---|---|---|

1 | ||||||

2 | ||||||

3 | ||||||

4 | ||||||

5 | ||||||

6 | ||||||

7 | ||||||

8 | ||||||

9 | ||||||

10 | ||||||

11 | ||||||

12 | ||||||

13 | ||||||

14 | ||||||

15 | ||||||

16 | ||||||

17 | ||||||

18 | ||||||

19 | ||||||

20 | ||||||

21 | ||||||

22 | ||||||

23 | ||||||

24 | ||||||

25 | ||||||

26 | ||||||

27 | ||||||

28 | ||||||

29 | ||||||

30 | ||||||

31 | ||||||

32 | ||||||

33 | ||||||

34 | ||||||

35 | ||||||

36 |

The simulated values were obtained based on 1000 simulation runs using . Table 1 compares Procedures 3.3 and 3.3, and Benjamini and Yekutieli (2001) in terms of these simulated mdFDR and power at several combinations of the aforementioned chosen values. It is to be noted that the Benjamini–Yekutieli procedure (BY for short) is an FDR controlling procedure for arbitrarily correlated -values. Although Procedures 3.3 and 3.3 are designed to control the mdFDR, it is worth comparing them to existing procedures that have similar dependence assumptions, namely, the BY. It would not be fair to compare to methods such as that of Benjamini and Hochberg (1995), since it requires independence or positive dependence between all -values, while Procedures 3.3 and 3.3 have more relaxed assumptions.

As seen from Table 1, the simulated mdFDR of Procedure 3.3 remains stably controlled across all correlation combinations (, 0, 0.4, 0.8; , 0.2, 0.5). Procedure 3.3 can outperform Procedure 3.3, in the sense of having higher power while still maintaining control of the mdFDR at the desired level , with weakly to moderately correlated data within subregions. However, if the data are moderately to largely correlated with small group sizes (), Procedure 3.3 can lose control of the mdFDR. Although unfortunate, this is not surprising, knowing that this type of adaptive procedure considered in the contexts of FDR or FWER control also becomes unstable with large correlations among the underlying test statistics. We can also see that as the subregion size increases ( from 9 to 400), the power decreases. This can, however, be attributed to the fact that a Bonferroni type combination of -values had to be considered, because of arbitrary dependence, to define subregion and location specific -values. In comparison, the BY procedure also seems to maintain control of the mdFDR, even though it is an FDR procedure. However, the BY procedure has smaller power in every scenario.

## 5 East African NDVI results

We use the Brillinger test as described in Section 2 to test for significant vegetation trend over the years separately for the four seasons at each location. Specifically, with representing the NDVI average for the th subregion, th location and th season in the th year, where , , and , we consider, for each fixed , the following model:

(3) |

and test is a constant signal vs. is monotonic in time, using the Brillinger test statistic and the corresponding approximate -value. A negative significant test statistic provides evidence of a monotonic decreasing trend, while a positive significant test statistic suggests an increasing monotonic trend.

We applied Procedures 3.3 and 3.3 (with and ) based on the -values for the above tests to the region to screen for significant seasonal vegetation changes over the years and the directions in which these changes are taking place.

As mentioned before, we used the semivariogram plot to determine the grid size () for each subregion, where each site’s averaged NDVI value over all years is . The empirical semivariogram plot for the NDVI data is shown in Figure 3. As seen from this plot, the range of the semivariogram is approximately 15 pixels, meaning NDVI values for locations with a Euclidean distance greater than 15 pixels apart are uncorrelated. Thus, it would be appropriate to choose the group size .

In particular, each subregion’s minimum -value () is used to represent the corresponding group in the stage 1 BH method in our procedure. Thus, if the locations of and are at least 15 pixels apart, we can consider according to this semivariogram plot that the corresponding subregions are independent.

After carefully considering various group sizes of , we choose , yielding groups, each with locations. Using to group the locations, only 1.37% of the group minima were closer than 15 pixels to another group minimum. Thus, we are satisfied that each group, represented by the group minimum, is nearly independent from the others.

This region has several bodies of water, including the African Great Lakes and part of the Indian Ocean, where there is no vegetation. Thus, the remote sensing pixels covering entirely water will correspond to a missing location in a subregion’s grid, that is, the subregions that straddle land and water will have .

The results of Procedures 3.3 and 3.3 for each of the four seasons are shown in Figures 4 and 5, respectively. Sites with a significant seasonal increasing change in vegetation are plotted in green. Sites with significant seasonal negative vegetation change are plotted in red. The nonsignificant sites are represented by tan. Using Procedure 3.3, we detected , , and pixels with significant increasing or decreasing changes in their respective seasonal NDVI averages (first dry season, long rainy season, second dry season and short rainy season). The second dry season has concentrated locations in coastal and central Tanzania with increasing average NDVI and the short rainy season has concentrated decreasing vegetation changes directly South of Lake Victoria, both of which are potentially important findings for land use managers. Using Procedure 3.3, the number of pixels with significant -values increases to 44, 42, 99 and 417 in their respective seasonal NDVI averages in generally the same regions as found from Procedure 3.3.

Overall, the results show increasing vegetation trends in the Northern hemisphere as well as coastal Eastern Tanzania. Decreasing vegetation trends are mostly concentrated directly South of Lake Victoria. Another noticeable finding is that the second dry season and short rainy season (which make up the last 6 months of the calendar year) are the seasons that contain the majority of the significance. These findings are consistent with historical evidence and other climate change investigations done in this region, as described below. However, this is the first study to reach these findings while maintaining control of a meaningful level of Types I and III errors.

## 6 Discussion and concluding remarks

The motivation of this paper lies in the fact that the currently available statistical approach to detecting vegetation changes over different locations in a region, like in East Africa, developed so far in the framework of testing multiple hypotheses [e.g., in Vrieling, de Beurs and Brown (2008)], may be questionable. Current methodologies have not taken into account the multiplicity by guarding against an overall measure of false discoveries, directional and nondirectional, and hence can potentially produce too many falsely discovered vegetation locations, more than what is statistically acceptable. It is important to avoid falsely discovered locations in the present application since the data findings could be used to relocate East African populations of people, livestock and crops, which is risky and costly.

Our findings, in terms of the proportion of regions with discovered vegetation change, are notably a smaller subset of the conclusions from other studies in this region, including Cole et al. (2000), Vrieling, de Beurs and Brown (2008), Usongo and Nagahuedi (2008) and Duveiller et al. (2007). First to analyze this data, Vrieling, de Beurs and Brown (2008) used trend tests on several phenology indicators at every location in the region and found a substantial proportion of “significant” changes in all indicators—as high as 44%. However, their conclusions failed to address any control on the error rate while testing thousands of hypotheses simultaneously.

Addressing this multiplicity issue is important, since making false claims about significant vegetation change is costly when it involves risking the livelihood of entire populations of people and livestock. Our proposed methods not only address the multiplicity by controlling an overall measure of combined directional and nondirectional false discoveries, the mdFDR, but also are developed with the idea of making them as powerful as possible by adequately capturing spatial dependency present in the data. Since sites tend to be dependent more locally than globally, we consider grouping the hypotheses into suitable clusters before developing these methods to be a way of capturing such local dependency. The idea of using grouped hypotheses has been successfully used in Clements, Sarkar and Guo (2011) in a two-stage format and to control the FDR. By augmenting this procedure to include directional errors, we are able to detect important directional vegetation changes in all four precipitation seasons in the East African region, while maintaining control of the mdFDR. It is important to point out, however, that Procedure 3.3 is one that offers an mdFDR controlling procedure in the present setting, that is, robust against spatial dependency but does not explicitly use such dependency (quantitatively speaking). Its adaptive version, Procedure 3.3, attempts to explicitly use such spatial dependency.

To reiterate, controlling an FDR related error rate incorporating directional errors, like the mdDFR, is the rationale behind proposing our procedures, since detecting areas with significant increasing or decreasing vegetation change is one of the primary objectives in the present application. Our procedure is an augmented version of an FDR controlling procedure, similar to Benjamini and Yekutieli (2005). There are other FDR controlling procedures proposed in similar cluster settings [Benjamini and Heller (2007), Pacifico et al. (2004)]. However, it remains to be determined if these procedures can be augmented as in Procedure 3.3 without losing control over the mdFDR. Then, one can consider these procedures as relevant competitors of ours and evaluate the performance of our procedure relative to them.

Utilizing the information gleaned after applying Procedures 3.3 and 3.3, more sophisticated modeling techniques can be applied to the locations with seasonal changes. By first detecting locations of change using multiple testing, we are protected from investigating too many falsely discovered locations. Specifically, spatio-temporal modeling and forecasting may be of interest to land use management to optimize utilization of the land based on projected NDVI.

There are a few areas for improvement to this study. Although the proposed procedures do not require any dependence assumption for the -values within each subregion, the theoretical proof demands that subregions be independent to maintain control of the mdFDR. It would be interesting to theoretically investigate the performance of the proposed three-stage directional procedures under more complex subregion dependence structures. Second, selection of the optimal subregion size is debatable, as with any tuning parameter. This parameter needs to be large enough such that one can reasonably assume the subregions are independent, yet keeping in mind that, in light of the simulation studies, larger group sizes yield less powerful procedures due to the Bonferroni adjustments. The idea of estimating the range of a variogram is one such way to select the subregion size .

## Appendix: Proof of Theorem 3.1

{pf}Let be the total number of ’s that have been rejected, and and , respectively, be the numbers of Types I and III errors that occurred out of these rejections. Then

where is the FDR, and is the (pure) directional FDR.

Let us consider using also as an indicator variable with (or ), indicating that the null hypothesis is true (or false). Then,

where is the number of significant subregions in the first stage of the procedure. Hence,

since [borrowing the idea from Guo and Sarkar (2012)]. Let be the number of significant subregions that would have been obtained if we had completely ignored the th subregion and applied the first-stage BH method to the rest of the subregion -values using the critical values , . Then, it can be shown that

Since the subregions are assumed independent, taking expectation in (Appendix: Proof of Theorem 3.1) and using that in (Appendix: Proof of Theorem 3.1), we see that

where is the proportion of true null hypotheses among the total null hypotheses in the th subregion.

We now work with the dFDR. With representing the true sign of the Brillinger’s monotonic trend parameter , can be expressed as follows:

from which we first have

Making arguments similar to those used for the FDR, we then have

Notice that , where is the cumulative distribution function of the standard normal. Therefore, assuming without any loss of generality that when , we have, for such ,

(8) | |||

The last inequality follows from the fact that, when , the distribution of is stochastically increasing in . Using (Appendix: Proof of Theorem 3.1) in (Appendix: Proof of Theorem 3.1), we see that

(9) |

where is the proportion of false null hypotheses among the total null hypotheses in the th subregion.

Thus, we finally have

(10) |

proving the desired result.

## Acknowledgments

The NDVI data set was collected as part of a Michigan State University research project, namely, the “Dynamic Interactions among People, Livestock, and Savanna Ecosystems under Climate Change” project (funded by the National Science Foundation Biocomplexity of Coupled Human and Natural Systems Program, Award No. BCS/CNH 0709671). We thank the anonymous referees for their constructive comments which have helped to improve the quality of the paper.

## References

- Abelson and Tukey (1963) {barticle}[mr] \bauthor\bsnmAbelson, \bfnmRobert P.\binitsR. P. \AND\bauthor\bsnmTukey, \bfnmJohn W.\binitsJ. W. (\byear1963). \btitleEfficient utilization of non-numerical information in quantitative analysis: General theory and the case of simple order. \bjournalAnn. Math. Statist. \bvolume34 \bpages1347–1369. \bidissn=0003-4851, mr=0156411 \bptokimsref \endbibitem
- Benjamini and Heller (2007) {barticle}[mr] \bauthor\bsnmBenjamini, \bfnmYoav\binitsY. \AND\bauthor\bsnmHeller, \bfnmRuth\binitsR. (\byear2007). \btitleFalse discovery rates for spatial signals. \bjournalJ. Amer. Statist. Assoc. \bvolume102 \bpages1272–1281. \biddoi=10.1198/016214507000000941, issn=0162-1459, mr=2412549 \bptokimsref \endbibitem
- Benjamini and Hochberg (1995) {barticle}[mr] \bauthor\bsnmBenjamini, \bfnmYoav\binitsY. \AND\bauthor\bsnmHochberg, \bfnmYosef\binitsY. (\byear1995). \btitleControlling the false discovery rate: A practical and powerful approach to multiple testing. \bjournalJ. Roy. Statist. Soc. Ser. B \bvolume57 \bpages289–300. \bidissn=0035-9246, mr=1325392 \bptokimsref \endbibitem
- Benjamini and Yekutieli (2001) {barticle}[mr] \bauthor\bsnmBenjamini, \bfnmYoav\binitsY. \AND\bauthor\bsnmYekutieli, \bfnmDaniel\binitsD. (\byear2001). \btitleThe control of the false discovery rate in multiple testing under dependency. \bjournalAnn. Statist. \bvolume29 \bpages1165–1188. \biddoi=10.1214/aos/1013699998, issn=0090-5364, mr=1869245 \bptokimsref \endbibitem
- Benjamini and Yekutieli (2005) {barticle}[mr] \bauthor\bsnmBenjamini, \bfnmYoav\binitsY. \AND\bauthor\bsnmYekutieli, \bfnmDaniel\binitsD. (\byear2005). \btitleFalse discovery rate-adjusted multiple confidence intervals for selected parameters. \bjournalJ. Amer. Statist. Assoc. \bvolume100 \bpages71–93. \biddoi=10.1198/016214504000001907, issn=0162-1459, mr=2156820 \bptokimsref \endbibitem
- Brillinger (1989) {barticle}[mr] \bauthor\bsnmBrillinger, \bfnmDavid R.\binitsD. R. (\byear1989). \btitleConsistent detection of a monotonic trend superposed on a stationary time series. \bjournalBiometrika \bvolume76 \bpages23–30. \biddoi=10.1093/biomet/76.1.23, issn=0006-3444, mr=0991419 \bptokimsref \endbibitem
- Chen, Jonsson and Tamura (2004) {barticle}[author] \bauthor\bsnmChen, \bfnmJ.\binitsJ., \bauthor\bsnmJonsson, \bfnmP.\binitsP. \AND\bauthor\bsnmTamura, \bfnmM.\binitsM. (\byear2004). \btitleA simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. \bjournalRemote Sensing of Environment \bvolume91 \bpages332–344. \bptokimsref \endbibitem
- Clements, Sarkar and Guo (2011) {binproceedings}[author] \bauthor\bsnmClements, \bfnmN.\binitsN., \bauthor\bsnmSarkar, \bfnmS. K.\binitsS. K. \AND\bauthor\bsnmGuo, \bfnmW.\binitsW. (\byear2011). \btitleAstronomical transient detection controlling the false discovery rate. In \bbooktitleStatistical Challenges in Modern Astronomy V (\beditor\binitsE. D. \bsnmFeigelson \AND\beditor\binitsG. J. \bsnmBabu, eds.) \bpages383–396. \bpublisherSpringer, \blocationNew York. \bptokimsref \endbibitem
- Cole et al. (2000) {barticle}[author] \bauthor\bsnmCole, \bfnmJ. E.\binitsJ. E., \bauthor\bsnmDunbar, \bfnmR. B.\binitsR. B., \bauthor\bsnmMcClanahan, \bfnmT. R.\binitsT. R. \AND\bauthor\bsnmMuthiga, \bfnmN. A.\binitsN. A. (\byear2000). \btitleTropical pacific forcing of decadal SST variability in the western Indian Ocean over the past two centuries. \bjournalScience \bvolume287 \bpages617–619. \bptokimsref \endbibitem
- Cressie and Wikle (2011) {bbook}[mr] \bauthor\bsnmCressie, \bfnmNoel\binitsN. \AND\bauthor\bsnmWikle, \bfnmChristopher K.\binitsC. K. (\byear2011). \btitleStatistics for Spatio-Temporal Data. \bpublisherWiley, \blocationHoboken, NJ. \bidmr=2848400 \bptokimsref \endbibitem
- Curran (1980) {barticle}[author] \bauthor\bsnmCurran, \bfnmP. J.\binitsP. J. (\byear1980). \btitleMultispectral remote sensing of vegetation amount. \bjournalProgress in Physical Geography \bvolume4 \bpages315–341. \bptokimsref \endbibitem
- Duveiller et al. (2007) {barticle}[author] \bauthor\bsnmDuveiller, \bfnmG.\binitsG., \bauthor\bsnmDefourny, \bfnmP.\binitsP., \bauthor\bsnmDesclee, \bfnmB.\binitsB. \AND\bauthor\bsnmMayaux, \bfnmP.\binitsP. (\byear2007). \btitleDeforestation in Central Africa: Estimates at regional, national and landscape levels by advanced processing of systematically-disturbed Landsat extracts. \bjournalRemote Sensing of Environment \bvolume112 \bpages1969–1981. \bptokimsref \endbibitem
- Foody (2003) {barticle}[author] \bauthor\bsnmFoody, \bfnmG. M.\binitsG. M. (\byear2003). \btitleGeographical weighting as a further refinement to regression modeling: An example focused on the NDVI–rainfall relationship. \bjournalRemote Sensing of Environment \bvolume88 \bpages283–293. \bptokimsref \endbibitem
- Guo and Sarkar (2012) {bmisc}[author] \bauthor\bsnmGuo, \bfnmW.\binitsW. \AND\bauthor\bsnmSarkar, \bfnmS.\binitsS. (\byear2012). \bhowpublishedAdaptive controls of the FWER and FDR under block dependence. Unpublished manuscript. Available at http:// web.njit.edu/~wguo/research.html. \bptokimsref \endbibitem
- Hayes and Sader (2001) {barticle}[author] \bauthor\bsnmHayes, \bfnmD. J.\binitsD. J. \AND\bauthor\bsnmSader, \bfnmS. A.\binitsS. A. (\byear2001). \btitleComparison of change-detection techniques for monitoring tropical forest clearing and vegetation regrowth in a time series. \bjournalPhotogrammetric Engineering and Remote Sensing \bvolume67 \bpages1067–1075. \bptokimsref \endbibitem
- Jackson, Slater and Pinter (1983) {barticle}[author] \bauthor\bsnmJackson, \bfnmR. D.\binitsR. D., \bauthor\bsnmSlater, \bfnmP. N.\binitsP. N. \AND\bauthor\bsnmPinter, \bfnmP. J.\binitsP. J. (\byear1983). \btitleDiscrimination of growth and water stress in wheat by various vegetation indices through clear and turbid atmospheres. \bjournalRemote Sensing of Environment \bvolume13 \bpages187–208. \bptokimsref \endbibitem
- OCHA (2011) {bmisc}[author] \bauthor\bsnmOCHA (\byear2011). \bhowpublishedEastern Africa drought humanitarian report No. 3. OCHA, UN Office for the Coordination of Humanitarian Affairs reliefweb.int. \bptokimsref \endbibitem
- Pacifico et al. (2004) {barticle}[mr] \bauthor\bsnmPacifico, \bfnmM. P.\binitsM. P., \bauthor\bsnmGenovese, \bfnmC.\binitsC., \bauthor\bsnmVerdinelli, \bfnmI.\binitsI. \AND\bauthor\bsnmWasserman, \bfnmL.\binitsL. (\byear2004). \btitleFalse discovery control for random fields. \bjournalJ. Amer. Statist. Assoc. \bvolume99 \bpages1002–1014. \biddoi=10.1198/0162145000001655, issn=0162-1459, mr=2109490 \bptokimsref \endbibitem
- Sarkar (2002) {barticle}[mr] \bauthor\bsnmSarkar, \bfnmSanat K.\binitsS. K. (\byear2002). \btitleSome results on false discovery rate in stepwise multiple testing procedures. \bjournalAnn. Statist. \bvolume30 \bpages239–257. \biddoi=10.1214/aos/1015362192, issn=0090-5364, mr=1892663 \bptokimsref \endbibitem
- Storey, Taylor and Siegmund (2004) {barticle}[mr] \bauthor\bsnmStorey, \bfnmJohn D.\binitsJ. D., \bauthor\bsnmTaylor, \bfnmJonathan E.\binitsJ. E. \AND\bauthor\bsnmSiegmund, \bfnmDavid\binitsD. (\byear2004). \btitleStrong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: A unified approach. \bjournalJ. R. Stat. Soc. Ser. B Stat. Methodol. \bvolume66 \bpages187–205. \biddoi=10.1111/j.1467-9868.2004.00439.x, issn=1369-7412, mr=2035766 \bptokimsref \endbibitem
- Tobler (1970) {barticle}[author] \bauthor\bsnmTobler, \bfnmW.\binitsW. (\byear1970). \btitleA computer movie simulating urban growth in the Detroit region. \bjournalEconomic Geography \bvolume46 \bpages234–240. \bptokimsref \endbibitem
- Tucker et al. (2005) {barticle}[author] \bauthor\bsnmTucker, \bfnmC.\binitsC., \bauthor\bsnmPinzon, \bfnmJ.\binitsJ., \bauthor\bsnmBrown, \bfnmM.\binitsM., \bauthor\bsnmSlayback, \bfnmD.\binitsD., \bauthor\bsnmPak, \bfnmE.\binitsE., \bauthor\bsnmMahoney, \bfnmR.\binitsR., \bauthor\bsnmVermote, \bfnmE.\binitsE. \AND\bauthor\bsnmSaleous, \bfnmN.\binitsN. (\byear2005). \btitleAn extended AVHRR 8-km NDVI data set compatible with MODIS and SPOT vegetation NDVI data. \bjournalInternational Journal of Remote Sensing \bvolume26 \bpages4485–4498. \bptokimsref \endbibitem
- Usongo and Nagahuedi (2008) {barticle}[author] \bauthor\bsnmUsongo, \bfnmL.\binitsL. \AND\bauthor\bsnmNagahuedi, \bfnmJ.\binitsJ. (\byear2008). \btitleParticipatory land-use planning for priority landscapes of the Congo Basin. \bjournalUnasylva \bvolume230 \bpages17–24. \bptokimsref \endbibitem
- Vrieling, de Beurs and Brown (2008) {bmisc}[author] \bauthor\bsnmVrieling, \bfnmA.\binitsA., \bauthor\bparticlede \bsnmBeurs, \bfnmK. M.\binitsK. M. \AND\bauthor\bsnmBrown, \bfnmM. E.\binitsM. E. (\byear2008). \bhowpublishedRecent trends in agricultural production of Africa based on AVHRR NDVI time series. Proceedings of the SPIE Conference: Remote Sensing for Agriculture, Ecosystems and Hydrology X. \bptokimsref \endbibitem