Applying spatial autocorrelation and logistic regression to analyze land cover change trajectory in a forested watershed

Taiwan is more susceptible to earthquakes and typhoons because of the unique geographical position in the Northwestern Pacific Rim. It is necessary to understand the effects of such frequent natural disturbances on land cover change for watershed management and disaster mitigation. This study applied both global and local Moran’s I statistics to analyze the spatial autocorrelation of landslide patches in a natural disturbed watershed in eastern Taiwan. The land cover maps extracted from FORMOSAT-2 satellite images acquired in 2005, 2008, and 2011. A logistic regression model validated to predict occurrence probability of change trajectory. Results showed that spatial pattern of homogeneous landslide patches presented on small scales; while heterogeneous landslide pattern was on larger scales. Landslide patches had higher positive spatial autocorrelation and indicated that as regional hotspots in study area. After the trajectory calculation and classification, two unchanged and five changed trajectories dominated the study area. The most significant transformation of land cover was from forest to landslide and channel. In addition, a fittest logistic regression model predicted the occurrence probability of change trajectory in the study area. Among these environmental variables for logistic regression, lithology was the most important spatial determinant for the change trajectories. Curvature and aspect variables were also significant. This spatial statistical model was helpful for predicting the occurrence probabilities of the change trajectories.


INTRODUCTION
Land cover is one of the most fundamental variables for watershed management and relates to many parts of the human and physical environments (Southworth et al. 2004). Change in land cover can transform biophysical surface features and lead to direct impacts on hydrological cycle, atmosphere, and the ecosystem services (Zomlot et al. 2017;Homer et al. 2020;Talukdar et al. 2020;Yin et al. 2020). Detection of land cover change can contribute vital information for the sustainable development and environmental protection of the watershed. Hence, it is necessary to assess land cover change at various scales for sustainable land management.
Land cove change, a complex, dynamic process, may take place due to natural and/or anthropogenic reasons and be affected by the local conditions, regional contexts, govern-ment policies, and external influences (Mertens and Lambin 2000;Verburg et al. 2010). However, change in land cover is more complicated and unique in Taiwan than other countries because of its two geographic locations. Firstly, the Taiwan orogeny is forming along a complex plate boundary between the Eurasian and Philippine Sea plates where Longitudinal Valley of eastern Taiwan, serving as a suture zone on land, experiences intensive collision, resulting in frequent earthquakes (Hao et al. 2019). Strong motion of earthquake is a principal trigger of slope failures, landslides, and rock falls (Wu et al. 2006a, b;Chang et al. 2011;Mozziconacci et al. 2013). Secondly, located in the northwest Pacific, Taiwan, one of the most intense tropical cyclones (typhoons) prone areas in the world, is affected by an average number of 3 -4 typhoons striking the island per year (Janapati et al. 2019). Eastern Taiwan is more likely to be affected by typhoon since most of the typhoons hitting Taiwan make landfalls along the east coast of Taiwan (Su et al. 2012). Typhoons' extreme rainfall cause hazardous landslides, debris flows, Terr. Atmos. Ocean. Sci., Vol. 32, No. 1, 35-52, February 2021 and floods (Chen 2007;Janapati et al. 2019). The particular positions in eastern Taiwan represent one of the few places in the world where earthquakes and typhoons frequently happen and, therefore, provide a rare opportunity to evaluate the role of natural disturbances in dynamics of land cover change for environmental protection and disaster mitigation. Taimali watershed, a typical example of forested watershed in eastern Taiwan, was chosen for this research due to fierce typhoons causing losses or damage of property, infrastructure, and lives over the past decade.
Land properties vary greatly along the different environmental gradients, so land cover types often exhibit spatial dependency or autocorrelation (Li et al. 2009). Cliff and Ord (1981) pointed out that spatial autocorrelation exists where there is a systematic spatial variation in values across a given area. In general, measures of spatial autocorrelation include global and local indexes. Global measures of spatial autocorrelation assume spatial stationarity and produce a single value that summarizes the entire study area (Wulder et al. 2007). On the contrary, local measures of spatial autocorrelation put emphasis on local variations within patterns of spatial dependence and are useful for identifying spatial relationships (Anselin 1995). Several techniques have been developed to measure both global and local spatial autocorrelations. Moran's I statistic, one of the most popular measures, was widely used to estimate general patterns of spatial dependency and frequently applied in the study of land cover change. The global Moran's I is a global parameter for measuring spatial dependency in whole large area (Cliff and Ord 1981); whereas, the local Moran's I examines the individual locations and recognizes the hotspots based on a comparison with the neighboring samples (Zhang et al. 2008). Since landslide occurrence is the most substantial change of land cover in the study area, Moran's I index was employed to identify the hotspots.
There are several ways to monitor land cover change, including the proportion, annual rates and spatial distribution of change. Furthermore, the change trajectory is also applicable to describe spatially explicit representations of change at the pixel level based on pixel histories (Rindfuss et al. 2004;Mena 2008). Land cover change trajectories, a succession of land cover types for a pixel in the raster images of time series over more than two observation years, not only allow the illustration of change patterns in the temporal dimension, but also consider their continuity and direction (Mertens and Lambin 2000;Zhou et al. 2008a, b;Ruiz and Domon 2009). Thus, trajectory analysis is applied to identify a sequence of land cover types for a given sampling unit under the influences of natural disturbances over time in the present study.
Trajectories of land cover change can be partitioned into stable and dynamic types. Stable type is featured by pixel histories that are unchanging over time; whereas dynamic type is characterized by pixel histories with varied classes across the time series (Mena 2008). Owing to the features of the dependent variable, binary logistic regression is a favorable method for analyzing the probability of land cover change trajectories and is utilized in this research. Logistic regression is a method for estimating the probability of occurrence of an event from dichotomous data by use of a set of variables (Bui et al. 2011). Many studies applied the logistic regression to estimate the probability of change for each pixel, especially for detecting instability factors and mapping landslide susceptibility (Bai et al. 2010;Lin et al. 2010;Nandi and Shakoor 2010;Bui et al. 2011). However, few studies have assessed the probability of change on each pixel for land cover change trajectories (Braimoh and Vlek 2005;Mena 2008).
This research combines the applications of spatial autocorrelation, trajectory analysis, and logistic regression to provide more spatially explicit and detailed information on land cover change dynamics under frequent natural disturbances, by utilizing FORMOSAT-2 satellite images acquired in 2005, 2008, and 2011 within the Taimali watershed. The main objectives of this study are to (1) assess spatial autocorrelation of landslide patches using global and local Moran's I statistics; (2) detect the dynamics of land cover change trajectories under frequent natural disturbances; and (3) establish a spatial statistical model based on logistic regression for projecting the occurrence probabilities of land cover change trajectories.

Study Area and Natural Disturbances
The Taimali watershed principally lies in Jialan and Liqui, Jinfeng Township, Taitung County of eastern Taiwan. The research area, a total area of 211.5 km 2 , is a typical mountainous watershed dominated by forest, and its elevation ranges from 3 to 3090 m at 57.6% average slope. Buildup and farmland account for less than 1.8% of the total study area. They are mostly located in the lower Taimali watershed where Jialan is the only one major settlement with a total population of about 1400 from 2012 to 2015 based on the population statistics of Taimali Household Registration Office (2015). On the other hand, forested land, covering 83.7% of land area, is the predominant land cover type during the study period. Forest chiefly occupies the upper and middle Taimali watershed; where the most notable land cover conversion occurs in relation to forest converting into landslides and channels over the study period. Under the influence of powerful earthquakes and typhoons, the Taimali watershed underwent remarkable land change. Recently, there were three large earthquakes whose moment magnitudes (M w ) were larger than 6 in southeastern Taiwan, including the Chengkung Earthquake (M w 6.8) in 2003, the Taitung Earthquake (M w 6.1) in 2006, and the Jiashian Earthquake (M w 6.3) in 2010. These earthquakes resulted in rock falls and landslides (Wu et al. 2006b;Hu et al. 2007), substantial coseismic ground displacements (Wu et al. 2006a), and many people injured (Rau et al. 2012). In addition, based on the typhoon database established by the Central Weather Bureau (CWB) of Taiwan, there were 23 typhoons hitting Taiwan during the period 2005 -2011. Among them, it is reported that three typhoons striking the Taimali watershed produced heavy rainfalls and triggered severe debris flows, landslides, and floods, including the Typhoon Haitang in 2005, the Typhoon Morakot in 2009, and the Typhoon Fanapi in 2010 (Fig. 1).

Image Classification and Accuracy Assessment
The land cover in the Taimali watershed relied on FORMOSAT-2 images obtained from the Spatial Information Research Center (SIRC) of National Taiwan University. The image classification performance was achieved by applying the ERDAS IMAGINE version 9.2 software (Leica Geosystems Inc.). The supervised classification algorithm with the maximum likelihood decision rule was applied to classify all the images. Three remote sensing images with 8-m resolution, acquired in 2005, 2008, and 2011, were classified into four land cover classes: landslide patches, forest matrix, human-made patches, and channel corridors. More comprehensive information about the land cover classification can refer to the research by Yeh and Liaw (2015). In order to evaluate the accuracy of classification of remotely sensed data, the current study utilized stratified random sampling to produce 256 reference points, whose values conducted a consistency check with the class values of the classified image (Jensen 2005).

Spatial Autocorrelation Assessment
For quantitative variables, Moran's I index is the most commonly used index to assess the global level of spatial autocorrelation (Cliff and Ord 1981;Li et al. 2009). Moran's I index indicates the degree of similarity or dissimilarity between the values of the variable considered and ranges approximately from +1 to -1 (Uuemaa et al. 2008). Moran's I value near +1 indicates clustering, while a value near -1 indicates dispersion, and 0 or near to 0 represents no spatial autocorrelation, that means a random pattern (Fernandes et al. 2011). Global Moran's I index is defined as following (Uuemaa et al. 2008;Li et al. 2009):  where X i and X j are the values of the observed variable at sites i and j; n the total number of sites; i = 1, 2, . . ., n and j = 1, 2, . . ., n (i ≠ j), μ the mean of all X i and X j ; W ij the weights representing proximity relationships between sites i and j; they form together a spatially 0 -1 contiguity matrix.
A correlogram, a graph where spatial correlation values (in our case Moran's I) are plotted on the ordinate against distance interval (lag) among sites on the abscissa, can vary only from +1 to -1, depending upon whether the correlation between locations is positive or negative (Rossi et al. 1992;Mander et al. 2010;Legendre and Legendre 2012). The correlogram, one of the most commonly used spatial structure functions, allows one to quantify spatial dependence and partition it amongst distance classes (Legendre and Legendre 2012). In this study, we applied the PASSaGE software (Rosenberg and Anderson 2011) for calculating Moran's I values with 20 different lag distances from 0.5 to 10 km. Each interval of lag distance was 500 m. Results were used to construct the correlogram graphs of Moran's I values versus lag distances.
The Moran's I values can be expected to be normally distributed, which allowed testing their significance through the comparison of the standardized deviates against those of the standard normal distribution (Cliff and Ord 1981). Usually, the Moran's I could be standardized to Z for testing the significance, and supposes approximate normality of Z (Salvador 2000). For instance, Z score larger than 1.96 is at the significant level of 0.05. In this study, we standardized the Moran's I values to Z values, and calculated statistical significance (p < 0.05).
We also calculated the local Moran's I, which is a measure of contagion that includes the effect of the spatial neighborhood (Fernandes et al. 2011). Local Moran's I has a spatial autocorrelation value for each sampling pixel, rather than the single value of the global Moran's I for an entire area. Local Moran's I index is described as following (Anselin 1995;Zhang et al. 2008).
In general, the spatial clusters of local Moran's I values could be categorized into four groups (Ping et al. 2004;Zhang et al. 2008): (1) high-high cluster (high values in a high value neighborhood); (2) low-low cluster (low values in a low value neighborhood); (3) high-low cluster (a high value in a low value neighborhood); and (4) low-high cluster (a low value in a high value neighborhood). The highhigh spatial clusters can be regarded as "regional hotspots", while the low-low clusters are "cool spots", and both highlow cluster and low-high cluster are spatial outliers (Zhang et al. 2008).

Method of Trajectory Calculation
This study develops a categorical map that demonstrates the land cover change trajectories at the pixel level. Each land cover class is described as a code in the land cover raster layer as 1, 2, 3, and 4 to represent the landslide patches, forest matrix, human-made patches, and channel corridors, respectively. The trajectory layer shows a sequence of codes for each pixel, which can be called trajectory codes. When there are fewer than 10 classes of land cover, the trajectory codes can be calculated by the following equation (Mena 2008;Wang et al. 2012Wang et al. , 2013: where TC ij is the trajectory code of the pixel at row i and column j in the trajectory layer; n is the quantity of the time points; (C t1 ) ij , (C t2 ) ij , and (C tn ) ij are the codes in land cover layers of different time points at the given pixel, and they must be ranked by time points in ascending order. Trajectory codes are established through formula calculation by applying the raster calculator in ArcGIS 9.3.1 in this study. For example, the trajectory of land covers that converts from forest to landslide, and then to landslide on a pixel over three time points is marked as 211 (forest → landslide → landslide), and so on.

Environmental Variables
One major goal of this study is to build a spatial model of probability of land cover change trajectories. Because the most remarkable change trajectory was from forested land into landslides in this study area, eight environment variables were chosen chiefly according to data availability and projection models for landslide occurrence in reference to previous studies (Lin et al. 2010;Bui et al. 2011;Zhou et al. 2018). Eight environmental variables were selected for following analysis, including lithology, distance to faults, rainfall, distance to rivers, elevation, slope, aspect, and curvature (Fig. 2). To obtain a better understanding of land cover change prone areas in each environmental variable, the continuous environmental variables (distance, elevation, slope, and so on) were classified into several small intervals by breakpoints in reference to previous research by Bui et al. (2011) and Zhou et al. (2018). These variables are described as follows.
Lithology is considered as an underlying driving factor for landslide occurrence (Bai et al. 2010;Lin et al. 2010;Bui et al. 2011). Geologic formations in the study area are composed of three strata including the Pilushan Formation (code 1), the Lushan Formation (code 2), and the Tananao Schist (code 3).
There are two faults in the study area, and their spatial distribution is similar to the letter Y shape. One fault passes across the central part of the Taimali watershed in a northsouth direction, and the other fault is elongated in a northeast to southwest direction at the northern part of the Taimali watershed. The distance to faults is calculated by measuring the straight-line distance to faults with ArcGIS Spatial Analyst Tools and is categorized into five classes: 0 -2000 , 2000 -4000 m, 4000 -6000 m, 6000 -8000 m, and > 8000 m.
Because there is no weather station in the study area, Kriging interpolation is used to estimate the average annual rainfall of the study area. An average annual rainfall map is produced by applying the ordinary Kriging method with the spherical semivariogram model based on the rainfall dataset from the neighboring 12 weather stations. Average annual rainfall is divided into four classes: 2000 -3000 mm, 3000 -4000 mm, 4000 -5000 mm, and 5000 -6000 mm. Under the effects of the natural disturbances, sediment from landslides deposited in channels caused the channel expansion. Therefore, distance to rivers is calculated by measuring the straight-line distance to the Taimali stream based on the remote sensing images in 2005 with ArcGIS Spatial Analyst Tools and is categorized into five classes: 0 -300 m, 300 -600 m, 600 -900 m, 900 -1200 m, and > 1200 m. The slope gradient is extracted from the DTM, and then divided into six slope gradient categories: 0° -5°, 5° -15°, 15° -25°, 25° -35°, 35° -45°, and > 45°.
The aspect is also extracted from the DTM, and the aspect image is classified and labeled with nine identifiers based on Miliaresis (2008 Curvature is mathematically defined as the change in slope angle along a very small arc of the curve, and can be expressed as the inverse of the radius of a circle that is tangent over a small arc of the curve (Ohlmacher 2007). The curvature is positive for convex landforms (e.g., hills and ridges), the curvature is negative for concave ones (e.g., depressions and valleys), and the curvature is zero for flat ones (Florinsky 2012).

Logistic Regression Analysis
In this study, each pixel of trajectory layer was reset as either value 1 if "with change" or value 0 if "without change". Because this is a binary variable, either 1 or 0, the logistic regression model can be applied to describe the probability of land cover change. The logistic regression model, based on the logistic function, can be expressed in the equation (Bai et al. 2010;Kleinbaum and Klein 2010): where P is the probability of occurrence of land cover change for each pixel, which ranges between 0 and 1 on a S-shaped curve. The value z can be substituted for the linear sum expression of 0 b plus 1 b times X 1 plus 2 b times X 2 , and so on to k b times X k : where 0 b , a constant term, is the intercept of the model, i b (i = 1, 2, …, k) are the slope coefficients of the model, and X i (i = 1, 2, …, k) are independent variables of interest. To find the best fitting set of parameters, maximum likelihood (ML) estimation is used to estimate the parameters in the model. Logistic regression was performed by the SPSS software package using the forward stepwise method.
Model fitting via logistic regression is sensitive to collinearities among the independent variables (Hosmer and Lemeshow 2000). The tolerance (TOL) and variance inflation factor (VIF) can be applied to discern multicollinearity in this situation. The TOL is defined as 1 -R 2 , where R 2 is the coefficient of determination calculated by regressing each variable on all the other independent variables (Allison 1999). The lower the TOL value, the higher the multicollinearity. A TOL value less than 0.2 demonstrates the presence of multicollinearity between independent variables, and a value less than 0.1 reveals serious muliticollinearity among them (Menard 2002). The VIF is obtained from the reciprocal of tolerance and shows an inflated value of variance of the coefficient in comparison with its value in the absence of collinearity (Allison 1999). TOL and VIF have no formal cutoff values to determine whether multicollinearity is present between variables. Allison (1999) suggests that a TOL value less than 0.4, that is, the VIF value greater than 2.5, may reveal the presence of multicollinearity. In this study, variables with TOL < 0.4 and VIF > 2.5 are rejected from the logistic regression analysis.
The R 2 statistic, also called the coefficient of determination, which represents the proportion of variation in the data accounted for by a linear regression model, is not appropriate for use with the logistic regression model (Hilbe 2009). The -2 Log likelihood (-2LL) is generally a measure of lack of fit, conceptually the opposite of R 2 in multiple regression, which is an index of variance (Osborne 2015). Cox and Snell R 2 and Nagelkerke R 2 , a kind of the pseudo-R 2 statistics, were appropriate for logistic models (Bai et al. 2010;Lin et al. 2010;Bui et al. 2011). Cox and Snell R 2 and Nagelkerke R 2 could indicate how much the variability of the dependent variable may be explained by all included predictor variables (George and Mallery 2010). The model fit is better when the values of Cox and Snell R 2 and Nagelkerke R 2 are higher. Both two pseudo-R 2 statistics were applied in this study to assess the model fit.
Implementing a validation for the prediction results is regarded as one of the important and essential procedures in prediction modeling (Chung and Fabbri 2003). The accuracy of logistic regression analysis in modeling the probability of land cover change can be evaluated by the area under the curve (AUC) values and the relative change intensity (RCI) values. The area under the receiver operating characteristic (ROC) curve provides a measure of the logistic regression model's ability to classify cases into one of two groups (Hosmer and Lemeshow 2000). The AUC value ranges from 0 to 1, and higher values indicate better fit. As a general rule, an AUC value greater than 0.7 is considered acceptable (Hosmer and Lemeshow 2000). Besides, the RCI is used to understand change intensity in a region comparing with other regions in a watershed. The RCI has a value of 1 when the change ratio in the region equals that in the entire watershed, and increases when the change ratio in the region is larger than that in the entire watershed. The RCI is calculated by the following equations (Wang et al. 2012): Rec Aec Ae = where RCI can be expressed by the proportion of the change ratio in a region (Rzc) to that in the entire watershed (Rec). The (Rzc) equals the proportion of changed area in the region (Azc) to the entire area of the (Az). The (Rec) equals the proportion of changed area in the entire watershed (Aec) to the entire area of the watershed (Ae).

Spatial Autocorrelation
Three land cover maps of the study area (2005, 2008, and 2011) were produced from FORMOSAT-2 images classification (Fig. 3). For these three images, the overall classification accuracies range from 83.98 to 90.23%, with Kappa coefficients ranging from 0.80 to 0.88. The resulting four classes were used in the trajectory analysis to detect dynamic processes of land cover change. The area of the forestland, the most dominant type in this watershed, was 198.0 km 2 in 2005, but reduced to 187.8 km 2 in 2008, with a further reduction to 177.1 km 2 in 2011. On the contrary, the area of the landslide patches rapidly enlarged from 2.1 km 2 in 2005 to 10.7 km 2 in 2008, and further extended to 17.6 km 2 in 2011. The landslide area dramatically increased due to the natural disturbances of earthquakes and typhoons.
Since the amount of total landslide area increased significantly over the study period, this research focused on the spatial autocorrelation analysis of landslides using both global and local Moran's I statistics. The global Moran's I on three different dates versus the 20 different distance lags was shown in Fig. 4. These correlograms showed that higher positive autocorrelations were present for landslide occurrence at shorter lag distances for all three dates. Although the Moran's I values at first lag distance of > 0 to ≤ 500 m were lower than that at second lag distance of > 500 to ≤ 1000 m in 2005 and 2011, Moran's I obviously decreased faster before the critical scale (about 2.5 km) for all three dates. Hence, 2.5 km was discerned as a critical scale in this study. This result revealed that spatial pattern of homogeneous landslide patches occurred at small scales. Although negative spatial autocorrelations occasionally occurred over some distance lags, Moran's I values became lower and more stable after the critical scale. This result indicated that spatial pattern of heterogeneous landslide patches occurred at large scales. Moreover, compared with the correlograms of all three dates, highest Moran's I value (0.50) was detected in 2011 at a lag distance of > 500 to ≤ 1000 m. Besides, all Moran's I values at each lag distance in 2011 were higher than that in 2005 and 2008. The finding suggested that stronger spatial autocorrelation of landslide patches existed in 2011. The area percentage of landslide patches increased from 1.0% in 2005, 5.1% in 2008, to 8.3% in 2011. The total areas of landslide patches almost expanded 8 times between 2005 and 2011 (Yeh and Liaw 2015). Based on the correlogram analysis, the higher Moran's I value at shorter lag distances in 2011 may be resulted from the strong clustering effect of landslide occurrences on specific areas.
In this study, local Moran's I values were standardized to Z scores and mapped by applying the ArcGIS 9.3.1 (Fig. 5). Local Moran's I can identify the landslide patches with high values at its neighbors as "hotspots". Through the spatial analysis and calculation using the PASSaGE software (Rosenberg and Anderson 2011), landslide patches with high-high spatial cluster can be regarded as a regional hotspot. This spatial distribution of regional hotspots provided site-specific information for management of degraded areas. Totally there were 27 grids whose Z score was larger than 1.96 at the significant level of 0.05 in 2005, whereas there were 72 grids with Z score larger than 1.96 in 2011 (Fig. 5). This result indicated that the number and area of landslide hotspots were expanding dramatically from 2005 to 2011. A research identifying the landslide hotspots in the entire island of Taiwan also revealed that southern Taiwan and catchments in the southern portion of eastern Taiwan became landslide hotspots after 2009, mainly attributed to Typhoon Morakot (Lin et al. 2017).

Land Cover Change Trajectories
The trajectory layer of land cover change was produced through trajectory calculation based on three land cover maps. There 44 different trajectories were identified (Fig. 6), including four unchanged trajectories (84.47%) and 40 changed trajectories (15.53%). The selection of the main trajectories in the study area followed the criteria suggested by Mena (2008). The trajectories covering less than 1% of the total watershed were excluded from the following analysis because they might be derived from the artifacts of classification error. As a result, there were top seven trajectories and all of them were larger than 1% of the total study area. These trajectories included two unchanged and five changed trajectories, and occupied 95.08% of the watershed area. Among them, the most dominant trajectory was the unchanged forest cover (222) comprising 80.32% of the watershed area. The other unchanged trajectory was channel cover (444) with 3.01% of the watershed area. Moreover,    the other five changed trajectories occupied 11.75% of the watershed area, including the trajectories 221 (5.14%), 211 (2.13%), 224 (1.77%), 212 (1.60%), and 244 (1.11%). The two first ones, 221 and 211, denoted the conversions from forests to landslides; whereas the trajectories 224 and 244 represented the conversions from forests to channels. Besides, the trajectory 212 signified natural revegetation on landslides or channels. The Taimali watershed is a representative of forested mountainous watershed and located in typical humid tropical monsoon climate. Previous studies revealed that earthquakes and typhoons often induced the conversions from forests to landslides (Chang et al. 2007;Lin et al. 2010) and sediment from landslides was frequently deposited in channels (Tsai et al. 2013). Furthermore, landslide scars and disturbed, bare land are usually revegetated quickly in the humid tropical climate (Douglas et al. 1999). Due to the physical environmental features of eastern Taiwan, several kinds of land cover change trajectories are prone to happen in this study area.

Change Analysis in Connection with Environmental Variables
After the trajectory layer was acquired, spatial distribution of change trajectories revealed the relationship between change trajectories and environmental variables. In this study, the values of selected main change trajectories in the trajectory layer were reset as 1, and then the values of other trajectories were reset as 0. Consequently, a binary map with value 1 or 0 of land cover changes were produced (Fig. 7). This binary map overlaid with the maps of eight environmental variables to measure the percentage of the land cover changes and RCI value (Table 1).
Pilushan and Lushan Formations comprised almost the entire watershed and each one covered nearly half of the watershed. However, RCI values larger than 1 found in the formations of Pilushan and Tananao Schist indicated that land cover changes were highly concentrated in them. In fact, more than 80% of the total changed area was located in Pilushan Formation. It seemed that Pilushan Formation was prone to land cover change because of relatively fragile geological conditions. A case study conducted in the Luye catchment, north to our study site, also showed that landslide ratio in Pilushan Formation is higher than Tananao Schist due to its weaker rock strength and more joint numbers (Shih and Chen 2010).
Each of the five classes of distance to faults covered about 20% of the study area. The RCI values descended with the increment of distance to faults. This result suggested that the high intensity of land cover change was situated within the region of 0 -6000 m distance to faults. The region adjacent to faults underwent land cover changes, probably due to the fragile geologic structure in the neighborhood of the faults. A study conducted by Lee et al. (2018) also revealed that the area of high landslide susceptibility was located near the boundaries of fault zones.
The rainfall regions with 4000 -5000 and 3000 -4000 mm took up the two greatest proportions of the study area. However, the two biggest RCI values were found in the rainfall regions of 4000 -5000 and 5000 -6000 mm with 1.85 and 1.37, respectively. Moreover, nearly 88% of the total changed area was located in both rainfall regions. Land cover conversion from forested land into landslides was the most remarkable change trajectory and mainly located in the upper and middle Taimali watershed over the study period. The upper and middle Taimali watershed with higher rainfall experience concentrated and extensive land cover changes. A study carried out in northern Vietnam by Bui et al. (2011) also indicated that highest concentration of landslides occurred in the two highest rainfall classes and a large amount of rainfall considerably impacts the landslide activity.
The changed ratios presented a decreasing trend with the increasing distance to rivers. The RCI values also displayed the same trend and the greatest RCI value 1.74 was found within the region of 0 -300 m distance to rivers. This result revealed that the high intensity and large extent of Fig. 7. Binary map of land cover trajectories. land cover change were closely adjacent and tied to rivers. A previous study performed in central Taiwan by Chang et al. (2007) demonstrated that landslides triggered by typhoons were apt to be approaching to stream and landslides triggered by earthquakes had an inclination to be close to ridge lines. Due to typhoons exerting a stronger influence on landslide occurrence than earthquakes during the study period, land cover change occurred mainly near the stream channel in the current research. The elevation region of 500 -1000 m took up the greatest proportion of study area at about 45%, but the RCI values larger than 1 were found at elevations between 1000 -2500 m. When explaining the relationship between elevation and land cover change, the correlation between elevation and rainfall and slope should be considered. In our study area, the regions with higher elevations are featured by greater annual rainfall and steeper slopes, both of which are important factors triggering landslides (Bui et al. 2011).

Variable Class Area (km 2 ) Area (%) Changed area (km 2 ) Changed ratio (%) RCI
The proportions of study area covered by slope gradients between 15 -45° were all larger than 20%. Overall, the RCI values increased with the ascending slope grades. The gradient regions of 5 -15°, 35 -45°, > 45° had the higher RCI values of 1.13, 1.23, and 1.46, respectively. The gradients can be portioned into two types. One is less steep slope (5 -15°) and the other is steep (35 -45° and > 45°). As mentioned above, the most remarkable land cover change related to conversions from forest cover to landslides and channels over the study period. The former occurred mainly on the steep slopes of upper and middle Taimali watershed, and the latter took place chiefly along the less steep Taimali streamside.
The three greatest proportions of area based on aspect classes were northeast, east, and southeast in sum about 48%. RCI values greater than 1 were located in the slope of NE, E, SE, and S. The outlet of Taimali watershed is eastward and typhoons with heavy rainfall often strike Taiwan from the east. As a result, the eastward slope was the dominant region where land cover change easily occurred.
The curvature classes of the study area were primarily comprised of concave and convex landforms that occupied 48 and 49% of the entire watershed, respectively. However, the RCI value of concave landform was higher than convex and flat landforms. A case study conducted in northeastern Taiwan by Lee et al. (2018) developed an integrated landslide hazard assessment approach and also indicated that concave-shaped topographic features had an adverse effect and was the dominant factor triggering landslides.

Multicollinearity Diagnostics
Multicollinearity among the variables may influence the accuracy of logistic regression model. TOL and VIF were applied to test the multicollinearity among the eight environmental variables. The results showed that some variables were correlated to other variables based on the criteria (rejected by TOL < 0.4 and VIF > 2.5) ( Table 2), especially for the variable of average annual rainfall (AAR) with lowest TOL value (0.125) and highest VIF value (7.997). Therefore, AAR variable was removed, and there was no multicollinearity again among other predictor variables ( Table 2).
The AAR variable is expected to have a relation with the dependent variable because the study site is prone to be attacked by frequent typhoon disturbances. Typhoons may bring a great amount of precipitation and drive changes in land cover. However, the AAR variable excluded based on the multicollinearity diagnostics indicated that AAR variable strongly correlated with the other variables. Because AAR variable was estimated from the surrounding 12 weather stations using the Kriging interpolation technique, there was some uncertainty about the rainfall estimates. Moreover, we attempted to apply other variables related to rainfall conditions. Annual maximum 1-day and consecutive 3-day precipitation data in 2005, 2008, and 2011, respectively, were also estimated from the neighboring 12 weather stations using the Kriging interpolation method to represent extreme precipitation events, which were tested again in multicollinearity diagnostics. Nevertheless, all annual maximum 1-day and 3-day precipitation variables were rejected by the criteria. Because the AAR was correlated to elevation and slope variables, it was omitted in the logistic regression analysis.

The Logistic Regression Model
Logistic regression can demonstrate the relationship between a target variable and multiple environmental variables, and project the occurring probability of land cover change trajectories. The measure of -2 Log likelihood (-2LL) indicates how well a model fits the data, and smaller -2LL values reveal that the model fits the data better (George and Mallery 2010). The -2LL values in the training data set decreased from 11622.574 in step 1 to 10056.432 in step 7; meanwhile, higher Cox and Snell R 2 and Nagelkerke R 2 values, 0.260 and 0.347, respectively, demonstrate a better model fit (Table 3). Besides, the predicted classification accuracy of the training data set was 82.4% for changed pixels, 62.0% for unchanged pixels, and 72.3% for the overall predicted accuracy. These results for model fit analysis indicated that this logistic regression model was reasonable.
The coefficients of independent variables from the fittest logistic regression model were showed in Table 4. Larger absolute values of the coefficient revealed that the independent variables were the main determinant for occurrence probability of change trajectories. Among these variables, lithology with highest absolute value of the coefficient was the most important spatial determinant for the change trajectories. In this study area, the most notable change trajectory was from forest cover into landslides. Lithology, regarded as an underlying driving factor for landslide occurrence, exerted the greatest influence on trajectories of land cover change. Its negative logistic coefficient indicated that the occurrence probability of change trajectories decreased with an increase in code value of lithological classes. Thus, trajectories of change were dominated mainly by Pilushan Formation in the study area. In addition, curvature was the second critical variable with negative value of logistic coefficient for the model. The concave landform had higher occurrence probability of change trajectory; while the convex landform had lowest occurrence probability. Moreover, aspect was the third significant determinant for the change trajectory. A negative coefficient of aspect showed that eastward slopes raised the likelihood of occurrence of change trajectory, including NE, E, and SE. The eastward slopes were on the windward slope and were prone to be hit by typhoons due to facing the typhoon tracks.
The other environmental variables were significant in explaining the change trajectory, but their coefficient values were lower. Both variables of slope and elevation had positive coefficient values, and implied that the likelihood of change occurrence risen with ascending gradients and altitudes. In contrast, both variables of distance to faults and rivers had negative coefficients, indicating that the probability of change occurrence decreased with increasing distance to faults and rivers.

Model Validation
This fittest logistic regression model was validated in the testing data set. The correct classification percentage was calculated to evaluate the performance of the model. The predicted classification accuracy of the testing data set was 81.6% for changed pixels, 64.5% for unchanged pixels, and 72.9% for the overall predicted accuracy. The model performance of testing data set was better than that of training data set. Moreover, the AUC value was high up to 0.785, which was considered as an acceptable discrimination (Fig. 8). It also indicated that classification result based on the logistic regression model was satisfactory.
The probability map of change occurrence was produced by applying the logistic regression model. Some studies presented that the occurrence probability could be further categorized into several classes (Bai et al. 2010;Nandi and Shakoor 2010;Bui et al. 2011). In this study, the probability map was divided into five equal interval classes of change likelihood, including very low (0.0 -0.2), low (0.2 -0.4), medium (0.4 -0.6), high (0.6 -0.8), and very high (0.8 -1.0) (Fig. 9). Among the five classes of probability, the class of very high occurrence probability covered the lowest percentage as 5.43% of the study area. In contrast, the class of very low occurrence probability had the largest percentage as 32.79% of the study area. This result was in accordance to the study by Can et al. (2005), which pointed out that high probability class should cover a small area. In addition, RCI values, measuring the concentration of land cover change in a region of the study area, were increasing with ascending classes. The greatest RCI value 2.68 was in class 5, but the total area was smallest. These results revealed that change trajectory considerably coincided in the zones which have higher probability of change and occupy a small area (Table 5).

CONCLUSION
Analyzing the change trajectories of land cover under the effects of natural disturbances, such as earthquakes and typhoons, can provide essential information for watershed management. In recent years, the effect of global warming on typhoon frequency and intensity has also received growing attention for disaster protection and mitigation. Under these circumstances, this study selected a mountainous watershed experiencing frequent earthquakes and typhoons in eastern Taiwan for understanding the influences of natural disturbances on land cover change trajectory. The land cover maps were extracted from FORMOSAT-2 satellite images acquired in 2005, 2008, and 2011. The spatial autocorrelation of landslide patches was analyzed by using global and local Moran's I statistics. A spatial statistical model based on logistic regression was developed for predicting the occurrence probabilities of land cover change trajectories.
The result showed that higher global Moran's I values were present for landslide occurrence at shorter lag distances for all three dates. It indicated that spatial pattern of homogeneous landslide patches occurred at small scales. On the contrary, spatial pattern of landslide patches was more heterogeneous at large scales. Besides, all global Moran's I values at each lag distance in 2011 were higher than that in 2005 and 2008, which suggested that stronger spatial autocorrelation of landslide patches existed in 2011. Based on local Moran's I values, landslide patches had higher positive spatial autocorrelation and presented as regional hotspots in study area. This spatial distribution of regional hotspots provided site-specific information for management of degraded areas.
Through formula calculation, selection criteria, and trajectory classification, two unchanged and five changed trajectories dominated 95.08% of the study area. Individually, the most dominant trajectory was the unchanged forest cover (222) comprising 80.32% of the watershed. However, the most significant transformation of land cover was from forest to landslide and channel in study area. Considering environmental variables, these change trajectories were related to fragile geological condition, the vicinities of faults and rivers, and geomorphological characteristics.
Based on the exploration of the relationship between change trajectory and environmental variables, land cover change occurring was subject to the geologic condition of Pilushan Formation and the vicinities of faults and rivers. In addition, land cover experienced highly concentrated changes at elevations of 1000 -2500 m, gradients greater than 35°, and on the eastward slope. Besides, the change intensity with concave landform was higher.
For multicollinearity diagnostics, the variable of average annual rainfall was rejected by the criteria. A fittest logistic regression model, which accuracy was high up to 72.3% at all, predicted the occurrence probability of change trajectory in the study area. Among these environmental variables for logistic regression, lithology was the most important spatial determinant for the change trajectories. Pilushan Formation dominated the probability of change occurrence in the study area. In addition, curvature was the second critical variable, and concave landform had higher probability of change trajectory. Aspect was also significant to influence change trajectory, especially for eastward. This spatial statistical model was helpful for predicting the occurrence probabilities of the change trajectories.   (MOST 103-2410-H-003-104). We express our sincere appreciation to the anonymous reviewers for their constructive comments and suggestions that helpfully improve the quality of this article.