How BFAST Trend and Seasonal Model Components Affect Disturbance Detection in Tropical Dry Forest and Temperate Forest

: Time series analysis has gained popularity in forest disturbance monitoring thanks to the availability of satellite and airborne remote sensing images and the development of different time series methods for change detection. Previous research has focused on time series data noise reduction, the magnitude of breakpoints, and accuracy assessment; however, few have looked in detail at how the trend and seasonal model components contribute to disturbance detection in different forest types. Here, we use Landsat time series images spanning 1994–2018 to map forest disturbance in a western Paciﬁc area of Mexico, where both temperate and tropical dry forests have been subject to severe deforestation and forest degradation processes. Since these two forest types have distinct seasonal characteristics, we investigate how trend and seasonal model components, such as the goodness-of-ﬁt ( R 2 ), magnitude of change, amplitude, and model length in a stable historical period, affect forest disturbance detection. We applied the Breaks For Additive Season and Trend Monitor (BFAST) algorithm and after accuracy assessment by stratiﬁed random sample points, and we obtained 68% and 86% of user accuracy and 75.6% and 86% of producer’s accuracy in disturbance detection, in tropical dry forests and temperate forests, respectively. We extracted the noncorrelated trend and seasonal model components R 2 , magnitude, amplitude, length of the stable historical period, and percentage of pixels with NA and tested their effects on disturbance detection employing forest-type speciﬁc logistic regression. Our results showed that, for all forests combined, the amplitude and stable historical period length contributed to disturbance detection. While for tropical dry forest alone, amplitude was the main predictor, and for the temperate forest alone, the stable historical period length contributed most to the prediction, although it was not statistically signiﬁcant. These ﬁndings provide insights for improving the results of forest disturbance detection in different forest types.


Introduction
Forests play a vital role in ecosystem goods and services to humanity, by providing energy, shelter, and livelihoods [1]. Human-induced land use and management practices, such as deforestation for agriculture, logging, plantation, or transitional subsistence farming, such as shifting cultivation, have led to forest cover loss [2]. Reliable information on forest cover and its changes is crucial for policymakers to design effective plans in forest conservation.
Forest disturbances have been detected and quantified using multitemporal spaceborne optical remote sensing, such as Landsat, which has been providing images for (1) Do TDF and TF differ in accuracy in forest disturbance detection by BFAST trend and seasonal model? (2) Is the difference in accuracy related to the BFAST components, such as the magnitude, goodness-of-fit, amplitude, and length of the stable historical period? (3) Within the same type of forest, do the BFAST trend and seasonal model components affect the correct detection of forest disturbance? (4) How does the percentage of pixels with no data values in a TS, caused by clouds and shadow removals in those pixels, affect forest disturbance detection?
We aim to answer these questions by exploring the potential of Landsat NDVI TS from 1994 to 2018 with the BFAST model for forest disturbance detection in both TF and TDF. These two types of forest in the Ayuquila River Basin have been subject to severe disturbances, with agriculture conversion and shifting cultivation in the TDF and selective logging in the TF. We compare the disturbance detection in these two types of forest by examining components in their trend and seasonal model, such as goodness-of-fit, magnitude, amplitude, length of stable period, and percentage of no data value in the TS, and associate these factors to the difference in disturbance detection. ate forests (TF), encompassing four research questions:

Materials and Methods
(1) Do TDF and TF differ in accuracy in forest disturbance detection by BFAST trend and seasonal model? (2) Is the difference in accuracy related to the BFAST components, such as the magnitude, goodness-of-fit, amplitude, and length of the stable historical period? (3) Within the same type of forest, do the BFAST trend and seasonal model components affect the correct detection of forest disturbance? (4) How does the percentage of pixels with no data values in a TS, caused by clouds and shadow removals in those pixels, affect forest disturbance detection?
We aim to answer these questions by exploring the potential of Landsat NDVI TS from 1994 to 2018 with the BFAST model for forest disturbance detection in both TF and TDF. These two types of forest in the Ayuquila River Basin have been subject to severe disturbances, with agriculture conversion and shifting cultivation in the TDF and selective logging in the TF. We compare the disturbance detection in these two types of forest by examining components in their trend and seasonal model, such as goodness-of-fit, magnitude, amplitude, length of stable period, and percentage of no data value in the TS, and associate these factors to the difference in disturbance detection. Figure 1 presents the overall data and data processing flow.

Materials and Methods
. Figure 1. Flowchart of the data and data processing methods applied in this study.

Forest in Ayuquila River Basin and Main Drivers of Forest Disturbances
The study area, the watershed of the Ayuquila River, is located in the western Pacific area of Mexico ( Figure 2). It is characterized by large topographic variations, with the elevation ranging from 260 to 2500 m above mean sea level; the average annual precipitation is 800-1200 mm and occurs mostly between June and October; the average monthly temperature ranges from 18 to 22 °C [58]. There are both temperate forests (TF) and tropical dry forests (TDF) in this area. TF is distributed in the higher elevation of the watershed, occupying about 12% of the watershed, including pine (Pinus spp.), fir (Abies spp.), and oak forests (Quercus spp.). TDF is distributed in the lower topographic areas, occupying about 24% of the watershed and composed of deciduous and semideciduous woodlands with much lower biomass density, canopy cover, and height than the TF. While TDF is used mainly for shifting cultivation, fuelwood extraction, cattle grazing, poles extraction for constructing fences, and harvesting mushrooms and medicinal plants [59], TF is exploited mainly for timber, although recently, avocado (Persea americana) plantations have

Forest in Ayuquila River Basin and Main Drivers of Forest Disturbances
The study area, the watershed of the Ayuquila River, is located in the western Pacific area of Mexico ( Figure 2). It is characterized by large topographic variations, with the elevation ranging from 260 to 2500 m above mean sea level; the average annual precipitation is 800-1200 mm and occurs mostly between June and October; the average monthly temperature ranges from 18 to 22 • C [58]. There are both temperate forests (TF) and tropical dry forests (TDF) in this area. TF is distributed in the higher elevation of the watershed, occupying about 12% of the watershed, including pine (Pinus spp.), fir (Abies spp.), and oak forests (Quercus spp.). TDF is distributed in the lower topographic areas, occupying about 24% of the watershed and composed of deciduous and semideciduous woodlands with much lower biomass density, canopy cover, and height than the TF. While TDF is used mainly for shifting cultivation, fuelwood extraction, cattle grazing, poles extraction for constructing fences, and harvesting mushrooms and medicinal plants [59], TF is exploited

Data Collection and Preprocessing
We employed Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI level-1 terraincorrected (L1T) surface reflectance scenes at path/row (29/65) from 1994 to 2018 available from Google Earth Engine (GEE, [60]). We applied cloud-and-shadow masks, with the FMASK algorithm [61] built into the GEE to mask out clouds and cloud shadows. Intercalibration of the different Landsat sensors has been recommended to reduce the sensor effects on the radiometry of the scenes [62]. In our case, we considered using uncalibrated scenes that would have minor effects on the change detection, since the selected historical period  includes at least two complete years where Landsat-7 ETM and Landsat-8 OLI were operational. Therefore, we assume that the BFAST algorithm can fit a stable model, considering the NDVI variations between the two different sensors.
We computed the NDVI from the near infrared (NIR) and red (RED) bands [6]: where the RED and NIR are band 3 (0.626-0.693 μm) and band 4 (0.776-0.904 μm), respectively, of Landsat TM and ETM+, and the band 4 (0.630-0.680 μm) and band 5 (0.845-0.885 μm), respectively, of Landsat 8 OLI. The compiled NDVI TS images have data values ranging from −1 to 1. Figure 3 shows an NDVI TS for TDF and TF from January 2013 to June 2018. TDF demonstrates a more pronounced seasonal pattern than TF.

Data Collection and Preprocessing
We employed Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI level-1 terraincorrected (L1T) surface reflectance scenes at path/row (29/65) from 1994 to 2018 available from Google Earth Engine (GEE, [60]). We applied cloud-and-shadow masks, with the FMASK algorithm [61] built into the GEE to mask out clouds and cloud shadows. Intercalibration of the different Landsat sensors has been recommended to reduce the sensor effects on the radiometry of the scenes [62]. In our case, we considered using uncalibrated scenes that would have minor effects on the change detection, since the selected historical period  includes at least two complete years where Landsat-7 ETM and Landsat-8 OLI were operational. Therefore, we assume that the BFAST algorithm can fit a stable model, considering the NDVI variations between the two different sensors.
We computed the NDVI from the near infrared (NIR) and red (RED) bands [6]: where the RED and NIR are band 3 (0.626-0.693 µm) and band 4 (0.776-0.904 µm), respectively, of Landsat TM and ETM+, and the band 4 (0.630-0.680 µm) and band 5 (0.845-0.885 µm), respectively, of Landsat 8 OLI. The compiled NDVI TS images have data values ranging from −1 to 1. Figure 3 shows an NDVI TS for TDF and TF from January 2013 to June 2018. TDF demonstrates a more pronounced seasonal pattern than TF. In the Landsat TS, the pixels with clouds and cloud shadows removed have values of "no data" (NA). Pixels with NA also include those in the data gap in Landsat-7 ETM+ and those with missing values, due to the difference in acquisition swath extent in different dates. The monthly distribution of the NA observations for the forest area is shown in Table 1. The percentage of NA observations is lower in the dry season from January to May and from November to December, and higher in the rainy season from June to October.

Forest Disturbance Detection
To monitor forest disturbances, we applied the BFAST Monitor algorithm [63]. In this study, we define "forest disturbance" as negative changes that occur in the forest canopy induced by human activities, such as logging and shifting cultivation, or natural events such as fires and hurricanes. A detailed description of the algorithm can be found in [34,47]. Briefly, this algorithm works by first defining a stable historical period using a reverse-ordered CUSUM test (ROC) [46]. Afterward, it derives a first-order harmonic seasonal model from data in the stable historical period (Equation (2)) and projects it into the monitoring period: where the intercept , slope , amplitude , and phases are the unknown parameters to describe the trend and seasonal patterns, which can be estimated and tested using ordinary least squares (OLS) techniques. The frequency that describes the number of annual observations is and represents the harmonic terms, and we adopt = 1. In the Landsat TS, the pixels with clouds and cloud shadows removed have values of "no data" (NA). Pixels with NA also include those in the data gap in Landsat-7 ETM+ and those with missing values, due to the difference in acquisition swath extent in different dates. The monthly distribution of the NA observations for the forest area is shown in Table 1. The percentage of NA observations is lower in the dry season from January to May and from November to December, and higher in the rainy season from June to October.

Forest Disturbance Detection
To monitor forest disturbances, we applied the BFAST Monitor algorithm [63]. In this study, we define "forest disturbance" as negative changes that occur in the forest canopy induced by human activities, such as logging and shifting cultivation, or natural events such as fires and hurricanes. A detailed description of the algorithm can be found in [34,47]. Briefly, this algorithm works by first defining a stable historical period using a reverse-ordered CUSUM test (ROC) [46]. Afterward, it derives a first-order harmonic seasonal model from data in the stable historical period (Equation (2)) and projects it into the monitoring period: where the intercept α 1 , slope α 2 , amplitude γ j , and phases δ j are the unknown parameters to describe the trend and seasonal patterns, which can be estimated and tested using ordinary least squares (OLS) techniques. The frequency that describes the number of annual observations is f and k represents the harmonic terms, and we adopt k = 1.
We define 1994-2015 as the historical period and January 2016-June 2018 as the monitoring period, considering that a disturbance would probably still be recognizable in the field verification (October 2018), within two and a half years after the disturbance occurred. The data in the monitoring period are checked for consistency with the model built in the stable historical period by way of moving sums (MOSUM) of the residuals [64]. The definition for each time point (t) in a monitoring period is expressed in Equation (3): wherev is the predicted variance from the trend in the historical period, n is the number of observations in the sample, and y t andŷ t are the observed and predicted values at time point (t), respectively; y t −ŷ t are the residuals in the MOSUM window, and h is the MOSUM bandwidth, defined as a fraction of the number of observations in the sample, often assigned as 0.25n for the disturbance detection [34,47,63]. When the null hypothesis of the stability of the seasonality pattern is rejected, a breakpoint is registered. The decision to reject this null hypothesis is based on a boundary condition; more detail is provided in [34].
If the TS pattern in the monitoring period is stable compared to the data in the stable historical period, all of the residuals and their MOSUM values should be near zero. Otherwise, a breakpoint is registered in the monitoring period. In addition, the residual (y t −ŷ t ) is calculated as a magnitude of change. Disturbances are often indicated by negative magnitudes; however, not all the breakpoints indicate disturbances because of interannual variability [47]. Therefore, a magnitude threshold can be applied to reduce false positive detections [22,47,48].
We chose the magnitude threshold by field verification carried out in 15-24 October 2018. We verified 12 breakpoints with a magnitude of change ranging from (|0.05|) to (|0.2|) and identified false detection corresponding to magnitude < (|0.2|). Note that, for forest disturbance detection, the magnitude value was negative, but we used its absolute value to facilitate interpretation of the results. We eliminated all breakpoints with magnitudes < |0.2| and took the earliest remaining breakpoint for each pixel as the "final" forest disturbance. During the field campaign, anthropogenic causes of forest disturbances including deforestation for Agave spp. plantations and fire-induced forest damage were identified. Since only forest pixels were analyzed, the changes detected are very likely the result of forest clearing for cattle ranching, shifting cultivation, or permanent agriculture such as avocado plantations.
The BFAST Monitor approach requires a historical period from which a stable TS can be modeled to detect the forest disturbance in the monitoring period [47]. We applied the BFAST algorithm to all the pixels in the watershed and applied a forest mask to keep only the forest disturbances detected in the forest area. For this, we visually updated a land cover map from 2011 to January 2016 to make certain the existing forests before the beginning of the monitoring period was being analyzed. We merged the non-forest categories, including agriculture, grassland, human settlements, riparian vegetation, and water bodies, into one non-forest category, and merged all forest types, including evergreen, deciduous, and shrubland, into a single category of forest. Landsat NDVI TS images were compiled in GEE [60]. The BFAST Monitor algorithm was applied through the bfastspatial package [37] (https://github.com/loicdtx/bfastSpatial, accessed on 20 May 2021) using R [64]. It is worth mentioning that BFAST implementation in GEE is now available, which might be of interest for future users [65].

Accuracy Assessment
We evaluated the accuracy of forest disturbances detection by adopting the protocol described in [66,67]. The protocol includes a sampling and response design that defines a representative set of samples and uses it to initially form a traditional confusion matrix [67], which was then corrected by a mapped area to estimate the unbiased area for each class. This procedure also had the purpose of generating variables of the seasonal model components from the correctly and falsely identified disturbances for the analysis of logistic regression; more details can be found in Section 2.5. We calculated the total number of samples by applying the equation for stratified random sampling (Equation (4)): where S Ô is the standard error of the estimated overall accuracy, and we adopt 0.01 as suggested by [66]. W i is the mapped proportion of area of class i and S i is the standard deviation of class i, and S i = U i (1 − U i ). The variables are listed in Table 2. We distributed the total sampling points in the four classes as follows: 50 points each to the two disturbance classes, TDF-disturbance and TF-disturbance, with the rest of the points distributed proportionally to the forest classes with no disturbance detected, TDF-no change and TF-no change. For the response design, we visually interpreted the points employing TS Sentinel-2 images available in GEE and high spatial resolution Google Earth images, considering both spatial and temporal accuracy. During the interpretation process, we deleted 25 points from the 327 points of TDF-no change, and 12 points from the 235 points of TF-no change because of the errors in the mask, which correspond to non-forest areas at the beginning of the monitoring period (2016).
We summarized the results of accuracy assessment in a confusion matrix using both the number of points and the proportion of the categories [66][67][68]. Using the area-weighted confusion matrix, we calculated the estimated accuracy indices, including overall accuracy and user's and producer's accuracies. We also estimated the areas of disturbance and permanent forests with 95% confidence intervals. The estimation assumes that the reference validation data used for the computation of the confusion matrix can represent a proportion of the whole area [66].

Logistic Regression with BFAST Model Components
We assume that certain BFAST model components, such as the magnitude of change, goodness-of-fit measured by coefficient of determination (R 2 ), the amplitude and the length of the stable historical period, and the percentage of NA in a TS play a role in the performance of change detection in TDF and TF. Therefore, we used the verified random points for correct and false detection to extract the NDVI TS and calculated the former mentioned BFAST model components.
We incorporated these components as independent variables in a logistic regression to test their importance in labeling the disturbance as correct or false (Equation (5)).
where ∝ 0 is the intercept, ∝ 1 . . . ..∝ 4 are the expected change in the logarithm odds of the dependent variable according to changes in each independent variable, x 1 . . . x 4 are the independent variables (predictors), and y is the dependent variable (correct or false positive detected disturbance). The forest type variable was encoded as binary, defining TF as 0 and TDF as 1, while all the other variables were numerical variables. The dependent variable was also encoded as binary, defining 0 as incorrect detection and 1 as correct detection.
As we are interested in the role of forest type in disturbance detection, we applied three logistic regressions, each with a different focus of forest type: all forests, only TDF, and only TF. From a statistical point of view, it is important to carry out the analysis on pixels with the same length of data availability. This is to exclude the possibility that the analysis is affected by the uneven data distribution within the sampled pixels. In satellite data, especially from the tropics, the contamination from clouds and cloud shadows cause invalid observations in a TS (Table 1) and possibly uneven distributions of the observations. However, in our case, since the pixels are from the same study area with a relatively small size, and the TS data include a long historical period (1994-2015), we assume that the difference in the distribution of observations in the TS data will be relatively insignificant. Nevertheless, the percentage of NA was included among the analyzed variables, as an indicator of the ratio of invalid/valid observations. Prior to fitting the logistic regression model, we conducted a variance inflation factor (VIF) statistical analysis to assess the degree of collinearity among the variables. According to [69], the VIF can be derived from the multiple correlation coefficient (Equation (6)).
where R i 2 is the multiple correlation coefficient of independent variable X i regressed on the remaining independent variables. For example, a variable X i with a coefficient R i 2 of 0.6 would result in a VIF of 2.5, and a higher coefficient corresponds to higher VIF. Following the recommendations from [70], we excluded highly redundant variables that have VIF > 2.5.

Change Detection with BFAST
In total, 496.7 ha were detected as forest disturbances in a period of two and a half years, January 2016-June 2018, with about 329.22 ha in TDF and 167.5 ha in TF. The definition of disturbance categories and areas of detected disturbances and forest covers with no change are summarized in Table 3. Figure 4 presents the forest disturbance detection with BFAST in NDVI TS of both TF and TDF. The yellow vertical line indicates the breakpoint that was registered as disturbance with a magnitude of change equal or lower than |0.2|. In Figure 4A, a disturbance was detected on 02 December 2016 in a TF, and in Figure 4B, a disturbance was detected on 10 December 2016 in a TDF.  Figure 4 presents the forest disturbance detection with BFAST in NDVI TS of both TF and TDF. The yellow vertical line indicates the breakpoint that was registered as disturbance with a magnitude of change equal or lower than |0.2|. In Figure 4A, a disturbance was detected on 02 December 2016 in a TF, and in Figure 4B, a disturbance was detected on 10 December 2016 in a TDF.

Accuracy Assessment of Detected Changes in TDF and in TF
Based on the independent reference sample, the labeled points were used to quantify the accuracy of the map of forest disturbances using confusion metrics weighted by area [68]. Table 4 presents the confusion matrices in both the number of samples and the area weighted probabilities. The change detection obtained a 96.52% overall accuracy estimated by the area weighted confusion matrix, and 68% and 86% of user's accuracy for disturbances in TDF and in TF, respectively. Because each of the two disturbance classes occupied a small proportion of the map area, 0.13% and 0.07% for disturbances in tropical dry forest and in temperate forest, respectively, omissions in these classes had a large negative impact on the precision of the two forest disturbance classes. In the case of TDF-disturbance, the confusion matrix with number of samples shows that 11 points that were selected in the TDF-no change stratum but were identified as disturbances (TDF-disturbance). This omission error for TDF-disturbance has a large impact on the producer's accuracy, since it occurs in a large TDF-no change stratum with 58% of the total map area. Because of this impact, the producer's accuracy for TDF-disturbance decreased from 68% in a sample-based error matrix to barely 4% in an area-weighted matrix. Similarly, the producer's accuracy of TF disturbances decreased from 86% to 4.1%. The results of accuracy assessment show that TF-disturbance had higher accuracy in both user's and producer's accuracies than TDF-disturbance. Table 4. Confusion matrix for accuracy assessment. The upper half of the confusion matrix was built with counts of samples and the lower half was based on the sample counts, but weighted by area proportions of the change detection map with categories of disturbances in tropical dry forest (TDF) and temperate forest (TF) and categories of no change in both forests.

Statistical Description of the BFAST Model Components
We extracted the BFAST trend and seasonal model components from the 100 verified pixels of disturbances (Table 4), including the amplitude, stable historical period length, magnitude, NA percentage, and R 2 . Since we are interested in the role of forest type in disturbance detection, we applied three models of logistic regressions, including all forests, only TDF and only TF.
We balanced the number of observations. In the forest model, we have 46 observations, 23 in each category (i.e., correct, or false). In the TDF-only model, we have 32 observations, 16 in each category, and in the TF-only model, we have 14 observations, 7 in each category. The statistics of the variables are presented in Table 5. Table 5. Statistics of BFAST model components in correct and false disturbance detections as explanatory variables in the logistic regression model. The statistics include the mean and standard deviations of each variable. In the TDF-only dataset, the correct detections correspond to higher R 2 , slightly higher NA percentage, same magnitude, higher stable historical period length, and higher amplitude value. This seems to suggest that within the same TDF, the correct detections correspond to the BFAST model, with high values in R 2 , stable historical period length and amplitude, and perhaps slightly higher NA percentage. In the TF-only dataset, the correct detections correspond to higher R 2 , lower NA percentage, higher magnitude, higher stable historical period length, and similar but slightly lower amplitude. This seems to suggest that within the same TF, the correct detections correspond to BFAST model, with high values in R 2 and stable historical period length. For both TDF and TF, when measured separately, higher values in R 2 , and stable historical period length are important for correct detection, while NA percentage does not seem to affect the correct or false detections.

Detection of Multi-Collinearity among the Variables
The VIF was calculated for each dataset (Table 6). For all forests, three variables showed a VIF > 2.5: forest type, amplitude and R 2 . We removed forest type and R 2 from the dataset for the logistic regression model (LRM). Although amplitude also has VIF higher than 2.5, after removing forest type, the value of VIF for amplitude was lowered, and therefore, we included amplitude in the regression. In the case of TDF, the only two redundant variables were amplitude and NA percentage, so we eliminated NA percentage, as all forest models included amplitude as a significant variable (Section 3.3.3). Finally, for TF, the redundant variables were amplitude and R 2 and we opted to include amplitude for the same reason as the case of TDF.

Logistic Regression Models
The three fitted logistic regression models are presented in Table 7. The first model with all forests has the amplitude and stable historical period length as independent variables to predict forest disturbance. Both are statistically significant (p < 0.05) for the correct disturbance detection. The TDF-only model has amplitude as a significant predictor, and the TF-only model has stable historical period length as a predictor, however, it was not statistically significant. The odds ratio is calculated as the Euler's constant (e) raised to the power of the estimate (β), and it represents the change in the odds of correctly identifying a disturbance depending on the independent variables. The odds ratio ranges from 0 to infinite. When its value is between 0 and 1, the independent variable decreases the odds of a correct detection; when it is larger than 1, the independent variable increases the odds of correct detection. All the independent variables have an odds ratio larger than 1, therefore all the independent variables in the regression model increases the odds of correct detection.

Forest Disturbances in TDF and TF
The disturbance detection in TF obtained a higher accuracy than in TDF. This is in line with the finding from [38] showing that the trend and seasonal model such as BFAST tends to yield higher accuracy in forests with less seasonality, such as those dominated by conifers, because this facilitates the discrimination between phenological changes and disturbances. Indeed, in our study area, the TDF has a very pronounced seasonality controlled mainly by humidity, with a rapid response to the onset of the rainy season, reflected in the drastic changes in NDVI. The variation in the precipitation from year to year might cause BFAST flags breakpoints that are triggered by interannual variations in the TDF phenology rather than disturbances [38].
Additionally, in ARB, disturbances in TF are caused mainly by selective logging, forest fires, and recently by permanent conversion to avocado plantations. These activities happen on a larger scale and cause more severe damage to the forest than the subtle changes, such as those caused by shifting cultivation, which is the main type of forest disturbance in TDF. In our experiment, most detected disturbances had led to complete forest clearing instead of vegetation thinning.

Logistic Regression with All Forests
The logistic regression for all forests shows that amplitude and stable historical period length are significant variables with positive coefficient. Although the odds ratio of amplitude (342.09) is much higher than that of stable historical period length (1.17), since amplitude has a small value range (i.e., 0-1), the effect of its high odds ratio is balanced out. In fact, stable historical period length has higher significance value (p = 0.003) in the regression than amplitude (p = 0.029).
Since amplitude often acts as an indicator of forest type: higher amplitude usually corresponds to more seasonal forest such as TDF (see Figure 3). At a first glance, the positive relation between amplitude and correct disturbance detection seems contradictory knowing that TDF has lower accuracy than TF. Nevertheless, a closer look at our data shows two key aspects: (1) the false disturbances in TDF demonstrate a similar value and variation as TF false and correct; and (2) correctly identified disturbances had higher amplitude values, especially in TDF (see Table 5). Therefore, the significant effect of amplitude corresponds to the effect of TDF. This topic will be further discussed in the next section.

Forest-Type Specific Logistic Regression
The forest-type specific logistic regression shows that for TDF, the single most important predictor of correct or false detection is amplitude. As typical areas covered by TDF correspond to very pronounced seasonal NDVI patterns (i.e., high amplitude values, see Figure 3); low amplitude values could be used as an indicator of models that were fitted rather inaccurately, and therefore, are more prone to false detections. We consider that a main reason for this inaccurate fit is related to the presence of artifacts in the TS, which is not reflected in the NA percentage.
In the case of TF, stable historical period length was the most important variable, however, it was not statistically significant with p = 0.057. Since we have a low number of observations in the false disturbance category (n = 7), the effect of stable historical period length on correct/false disturbance identification might be obscured. Although it would be desirable to follow a sampling procedure to increase the number of false disturbances, it is difficult to identify these areas prior to the verification procedure.

General Considerations
The other two BFAST components, NA percentage and magnitude, did not show as significant variables in the disturbance detection. NA percentage was used as an indicator of the data quality in the entire TS; however, it does not necessarily represent the data quality in the modelled period (covered by the stable historical period length). We believe that a more relevant variable would be an indicator for the data quality in the stable historical period length, that will ultimately affect the BFAST model fit and disturbance detection. Additionally, other data quality indicators, such as those related to the presence of artifacts, might be of greater importance to the fit of the BFAST model and therefore affects the disturbance detection.
Magnitude is the variable that is often used in previous studies for eliminating false disturbances. Our study adopted the magnitude threshold of |0.2|, which resulted from field verification. The lack of explanatory power of magnitude in the logistic regression is related to the fact that a threshold was already applied, and therefore, magnitude values higher than |0.2| did not gave more information for discriminating between false and correct disturbances.
Previous studies [38] showed that integrating components such as magnitude, slope, and amplitude for forest-type specific models improved disturbance detection. Based on our results, we found that, in addition, stable historic period length is an important factor as well in refining the results of forest disturbance detection with BFAST model.

Conclusions
This study evaluated the contribution of the components of BFAST model and the percentage of NAs in a time series to forest disturbance detection. Time series NDVI from Landsat spanning 1994-2018 was applied, employing BFAST model in both temperate forest and tropical dry forest and its accuracy was evaluated by 624 random sample points. Afterwards, BFAST model components, including goodness-of-fit (R 2 ), magnitude, amplitude, length of the stable historical period, and the percentage of NA, were extracted from the time series data of the verified disturbances. The main findings are the following: (1) The disturbance detection yielded higher accuracy in the temperate forest than in the tropical dry forest.
(2) The difference in accuracy is associated with the BFAST model components, especially the amplitude and stable historical period length. (3) Within the same forest type, the components affect disturbance detection. For TDF, the amplitude was the statistically significant factor, and for TF, the stable historical period length contributed most to the detection, however, with no statistical significance. (4) The percentage of NAs in a time series was not an important factor in disturbance detection.
This study confirmed that BFAST model components, though varying with forest type, play an important role in forest disturbance detection. To explore the full potential of these methods to map and quantify forest disturbances, future studies should evaluate different BFAST model components in combination with techniques that help reduce the signalto-noise ratio. For example, variables that represent data quality in the stable historical period might be more relevant as a data quality indicator in the BFAST model fitting, and therefore, help distinguish between correct and false detected disturbances. Funding: This research was funded by the Consejo Nacional de Ciencia y Tecnologia (CONACyT) grant number "Ciencia Basica" SEP-285349.

Data Availability Statement:
The time series Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI data compilation and preprocessing were carried out in GEE. The scripts in R for disturbance detection by BFAST monitor are available through the link https://github.com/alequech/Change-forest-cover (accessed on 20 May 2021). The extraction of the trend and seasonal model components is available through the link https://github.com/JonathanVSV/BFAST_params_extr (accessed on 20 May 2021).