# Monitoring Changes in the Enhanced Vegetation Index to Inform the Management of Forests
## 2 Materials and Methods
### Statistical Methods  
#### 2.3.1 EVI Change Detection  
We used Breaks for Additive Seasonal and Trend (BFAST, [54; 55]) statistical methods to detect abrupt (i.e., breakpoints or breaks) and gradual changes (i.e., trends) in the forest EVI across time series pixels while accounting for seasonal variation. (Figure A1 in Appendix A.2 shows a schematic of the different types of breaks and trends mentioned here). BFAST's algorithms have been used in several studies analyzing vegetation dynamics [56; 57; 58; 59; 60; 61]. We used the BFAST R package (version 1.6.1) and three algorithms-standard bfast (hereafter bfast), bfast01 and bfastclassify.  
Thus, bfast decomposes a time series into seasonal, trend, and remainder components and integrates these with methods for detecting and characterizing time series changes [54; 55; 62]. Trend breakpoints are detected using piece-wise linear regression over segments of the detrended time series. The bfast function finds one or more trend breakpoints (hereafter, breaks) in a pixel time series. Among the key outputs, bfast provides the magnitude, direction, or type (positive or negative) and the approximate timing of each break (together with 95% confidence intervals, hereafter 95% CIs (Figure 2).  
Furthermore, we ran the bfast01 algorithm, which provides input to the bfastclassify algorithm to categorize trends (Figure 2). Contrary to the standard bfast algorithm, bfast01 runs a seasonally adjusted regression model to find the single biggest break (if any) across the entire time series. If a statistically significant break is found, then the algorithm fits two trend segments based on piece-wise linear regression, one on each side of the break. If no statistically significant break is found, the algorithm fits a single trend line to the time series. Each pixel time series processed by bfast01 is classified by bfastclassify into a trend type as proposed by [PERSON] et al. [57].  
One of the key parameters for BFAST's algorithms is \(h\), which determines the temporal window size between potentially detected breaks. \(h\) should be set according to the frequency and magnitude of the particular disturbance events and vegetation responses attempted to be captured [60; 63]. Often, these details are not known in advance. Here, we set \(h\) to 0.05, which is equivalent to one year of data given our time series (23 time steps per year when using the MODIS EVI 16-day temporal resolution). This choice of \(h\) means that once a break is detected, a new break can only be detected after one year of observations has passed. A sensitivity analysis showed that \(h=0.05\) provided the best compromise to leverage the disturbance data used here and produce annual summaries (Tables A2 and A3).  
After prepossessing and filtering out poor-quality pixels (e.g., pixels with high cloud cover, high aerosols, possible shadow, etc.), we ran BFAST algorithms on two decades of forest EVI time series (2003-2022). Analyses were run using a 0.01 significance level for bfast and bfast01 and 0.05 for bfastclassify. All BFAST analyses presented were conducted using \(h=0.05\). (For sensitivity analyses of \(h\), see Tables A2 and A3). To attribute forest EVI negative breaks to disturbances, we overlapped detected break and disturbance pixels in space and time. We assume that an overlap meant that a particular disturbance was linked to a particular EVI break.