Comparison of Modeling Grassland Degradation with and without Considering Localized Spatial Associations in Vegetation Changing Patterns

Grassland ecosystems worldwide are confronted with degradation. It is of great importance to understand long-term trajectory patterns of grassland vegetation by advanced analytical models. This study proposes a new approach called a binary logistic regression model with neighborhood interactions, or BLR-NIs, which is based on binary logistic regression (BLR), but fully considers the spatio-temporally localized spatial associations or characterization of neighborhood interactions (NIs) in the patterns of grassland vegetation. The BLR-NIs model was applied to a modeled vegetation degradation of grasslands in the Xilin river basin, Inner Mongolia, China. Residual trend analysis on the normalized difference vegetation index (RESTREND-NDVI), which excluded the climatic impact on vegetation dynamics, was adopted as a preprocessing step to derive three human-induced trajectory patterns (vegetation degradation, vegetation recovery, and no significant change in vegetation) during two consecutive periods, T1 (2000–2008) and T2 (2007–2015). Human activities, including livestock grazing intensity and transportation accessibility measured by road network density, were included as explanatory variables for vegetation degradation, which was defined for locations if vegetation recovery or no significant change in vegetation in T1 and vegetation degradation in T2 were observed. Our work compared the results of BLR-NIs and the traditional BLR model that did not consider NIs. The study showed that: (1) both grazing intensity and road density had a positive correlation to vegetation degradation based on the traditional BLR model; (2) only road density was found to positively correlate to vegetation degradation by the BLR-NIs model; NIs appeared to be critical factors to predict vegetation degradation; and (3) including NIs in the BLR model improved the model performance substantially. The study provided evidence for the importance of including localized spatial associations between the trajectory patterns for mapping vegetation degradation, which has practical implications for designing management policies to counterpart grassland degradation in arid and semi-arid areas.


Introduction
Grasslands are the most widely distributed natural ecosystems, and play a vital role in global carbon recycling, as well as serve as an important livestock husbandry base for agricultural activities [1,2].Grassland degradation has attracted wide attention worldwide.Grassland covers more than 40% of China's total land area, a large percentage of which is found in the Inner Mongolia Autonomous Region (IMAR) [3].Due to arid and semi-arid climate conditions, as well as the impact of human activities during the past several decades, the grassland ecosystems showed a limited capacity to withstand human interference; as a result, grassland degradation has been widely observed in IMAR [4][5][6].Grassland degradation is a complex process in which grassland productivity, vegetation coverage, grass height, soil biodiversity, or community composition decreases [7][8][9][10].Numerous studies focusing on the characteristics, processes, and mechanisms of grassland degradation and on the restoration of degraded vegetation have been conducted [1,7].Development in remote sensing technology in the past few decades has significantly advanced the work in vegetation cover monitoring [11][12][13][14][15]. Huge datasets were accumulated by remote sensing sensors, which greatly facilitated studies in grassland ecosystems.Based on the integration of remotely sensing data and other socio-economic data, a few important issues related to grassland degradation were analyzed.For example, steppe degradation in arid and semi-arid areas could be the result of both climate conditions, particularly due to limited precipitation, and intensified human activities [16,17].Furthermore, studies on the driving factors of grassland degradation, both natural and human-induced, presented some critical insights for understanding the pathways of vegetation deterioration and their mechanisms [18,19].Typical examples of natural driving factors included climatic conditions, local environments (e.g., topography such as elevation, slope, or aspect), and soil properties [20][21][22].Some human-induced factors, including overpopulation, overgrazing from livestock, and transportation, were also identified as being threatening factors for the grassland ecosystem [23][24][25][26].Intensified human disturbance due to increased accessibility to remote regions, as a result of the construction of more road networks, have led to vegetation degradation in many parts of the world [17,[27][28][29].
While the constraints from natural factors are difficult to remove, threats from human-induced factors are relatively easier to alleviate by means of designing more sound land use policies, e.g., rotational grazing, fenced grazing, or through vegetation recovery projects [30].Therefore, quantifying vegetation degradation due to human-induced factors attracted much attention from grassland policymakers and experts in academia [31,32].Evans and Geerken proposed a method called residual trend analysis, or simply RESTREND, to discriminate the roles on dryland degradation played from climatic and human-induced factors based on the residuals of regression between the normalized difference vegetation index (NDVI) and precipitation [33].RESTREND has since been applied widely to explore the associations between grassland degradation and human-induced vegetation changes [31, [34][35][36][37].
Remote sensing and spatial analysis can be combined to map the trajectory patterns (TPs) of grassland vegetation based on historical time-series data.It is of great value to understand the factors that explain TPs, and preferably, to predict future vegetation degradation.However, land use change is a complex process that is not easy to model [38][39][40].Pattern based land-use models have been widely applied to map land-use changes [41].Relying on observations from a list of potential driving variables, pattern based land-use models were implemented to quantify the relationships of land-use patterns and those variables of interest [42].In addition, localized spatial associations, or neighborhood interactions contained in the land use patterns, have proved to be helpful in studying land use changes [43].Cellular automata, usually referred to as a CA model, is a typical example that uses the links or characteristics of neighborhood interactions in land use patterns to model the dynamics of land use expansion [44].Neighborhood interactions (NIs) reveal the degree of dependence between TPs in a local space, which are usually defined by a neighborhood matrix, and NIs could be important factors for determining further directions of TPs [45].By characterizing neighborhood interactions, global as well as spatially varied, or localized associations between different TPs were measured [46].There have been a few methods proposed for quantifying such neighborhood interactions.For example, Verburg et al. [46] proposed an index called the enrichment factor (F); for every location defined in a grid dataset, the characteristics of interactions between land cover types could be measured.In the case of modeling grassland degradation, NIs, implying possible spatial interactions between TPs, might also be helpful.
The primary aim of this paper is to validate the effectiveness of including NIs along with other human-induced driving factors in modeling grassland vegetation degradation.Binary logistic regression (BLR) has been widely applied to studying cause-effect relationships.BLR is capable of predicting the odds of being a case or not (responses from a binary variable) based on the values of the independent variables (predictors), and thus is fit for studying vegetation degradation [18].We compared two models based on binary logistic regression: the traditional binary logistic regression model (simply referred to as BLR hereafter), and a binary logistic regression model with neighborhood interactions (referred to as BLR-NIs).With the support of historical data, it is possible to predict the probability of future vegetation degradation by the BLR models.

Study Area
The study area was in the Xilin river basin located between 43  1).It is one of the most representative steppe zones in China.Out of the total area about 10,000 km 2 , approximately 90% of the region is covered by grassland [44].The area belongs to the continental semi-arid grasslands of the central Asian steppe ecosystem [47].The annual mean temperature of the area ranges from 1 • C to 4 • C [48], while annual mean precipitation is about 350 mm, most of which concentrates between June and August [49].Needle grass (scientifically known as Stipa grandis) and false wheatgrass (Leymus chinensis) are two dominant species in the area [50], and dominant soil types include chestnut, chernozem, and Aeolian [51].The growing season (from May to September) of local plants usually starts in late April and ends in late September.Agriculture is the primary income for the local people.Breeding livestock on the grassland is one of the most important activities for the farmers.Fenced grazing and rotational grazing were adopted since the early 2000s in some parts of the region to avoid overgrazing from the grassland animals.
Sustainability 2018, 10, x FOR PEER REVIEW 3 of 15 neighborhood interactions.For example, Verburg et al. [46] proposed an index called the enrichment factor (F); for every location defined in a grid dataset, the characteristics of interactions between land cover types could be measured.In the case of modeling grassland degradation, NIs, implying possible spatial interactions between TPs, might also be helpful.The primary aim of this paper is to validate the effectiveness of including NIs along with other human-induced driving factors in modeling grassland vegetation degradation.Binary logistic regression (BLR) has been widely applied to studying cause-effect relationships.BLR is capable of predicting the odds of being a case or not (responses from a binary variable) based on the values of the independent variables (predictors), and thus is fit for studying vegetation degradation [18].We compared two models based on binary logistic regression: the traditional binary logistic regression model (simply referred to as BLR hereafter), and a binary logistic regression model with neighborhood interactions (referred to as BLR-NIs).With the support of historical data, it is possible to predict the probability of future vegetation degradation by the BLR models.

Study Area
The study area was in the Xilin river basin located between 43°26′ and 44°29′ North and 115°32′ and 117°12′ East in the Inner Mongolia Autonomous Region (IMAR), China (Figure 1).It is one of the most representative steppe zones in China.Out of the total area about 10,000 km 2 , approximately 90% of the region is covered by grassland [44].The area belongs to the continental semi-arid grasslands of the central Asian steppe ecosystem [47].The annual mean temperature of the area ranges from 1 °C to 4 °C [48], while annual mean precipitation is about 350 mm, most of which concentrates between June and August [49].Needle grass (scientifically known as Stipa grandis) and false wheatgrass (Leymus chinensis) are two dominant species in the area [50], and dominant soil types include chestnut, chernozem, and Aeolian [51].The growing season (from May to September) of local plants usually starts in late April and ends in late September.Agriculture is the primary income for the local people.Breeding livestock on the grassland is one of the most important activities for the farmers.Fenced grazing and rotational grazing were adopted since the early 2000s in some parts of the region to avoid overgrazing from the grassland animals.

Data Sources
In order to study the vegetation degradation in relationship to human activities (including grazing intensity and road network density), we acquired NDVI data, climatic data, grazing intensity data, and road network data.The NDVI data were obtained from the NASA MODIS MOD13Q1 vegetation index dataset for the years 2000-2015.MOD13Q1 were provided every 16 days at 250-meter spatial resolution.Maximum annual NDVI (NDVI max ) was computed and represented the annual climax productivity of the grassland vegetation.The climatic data were obtained from China's Meteorological Data Sharing Service System (http://cdc.cma.gov.cn), which included monthly precipitation data collected from 680 weather stations across the entire country.Ordinary Kriging spatial interpolation was adopted to generate monthly precipitation thematic maps, considering its capability of giving unbiased predictions with minimum variance for the target variable.Grazing intensity (GI) was one of the most representative human-induced factors that caused grassland degradation.GI data were collected from the local statistical yearbook during the study period.There were five recorded animal species, that is, sheep, horse, buffalo, camel, and donkey.Since different livestock have varied impacts on grass consumption and vegetation damage, each of the livestock was converted to a standard animal unit, called a sheep unit, through a conversion coefficient to compute the grazing intensity [18].By converting all of the animal species to sheep unit, it is more convenient to model the overall grazing effect.Local farmers helped determine the coefficients, which represented the magnitude of grass consumption and vegetation damage from a particular animal type relative to that from one sheep, through a questionnaire survey.Specifically, one horse is equivalent to six sheep, and thus was multiplied by six to transform into the sheep unit.Similarly, cattle, camel, and donkey were multiplied by five, seven, and three, respectively.The road network was collected from an open source project called Open Street Map (OSM), as well as traffic road maps from the local transportation bureau.The line density analysis from the ArcGIS Spatial Analyst toolbox was used to compute road density (RD) for both OSM and the traffic road map, separately.Then, the two RD maps were overlaid to get the maximum value as RD.In the end, all of the dataset were resampled to 250 m × 250 m, and clipped using the boundary of the study region.

Methodology
The NDVI-based residual trend (RESTREND-NDVI) method proposed by Evans and Geerken [33] was used to extract human-induced grassland degradation for two consecutive periods, T 1 (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and T 2 (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), respectively.Here, the purpose of designing two periods was to test the power of predicting patterns of vegetation change for the next stage (namely T 2 ) based on the previous stage (T 1 ), because vegetation change not only depends on spatial context, but also presents time-lag, which indicates the linkage of patterns of vegetation change in a temporal context.NDVI has a significant positive correlation with rainfall in arid and semi-arid drylands [34].By performing regression analysis between precipitation data and annual NDVI max at the pixel level, the associated residuals in each period could be mapped in a linear regression [33].A positive trend in residuals indicates an increasing NDVI signal relative to precipitation, i.e., vegetation recovery (RR), whereas a negative trend in residuals indicates a declining NDVI signal relative to precipitation, i.e., vegetation degradation (DD); if the residuals show no significant trend, then no significant change (NS) in vegetation is detected [33,34,52].
Based on the above patterns (RR, DD, or NS) of vegetation change identified for each period (T 1 and T 2 ), vegetation degradation (VD) was redefined based on the conversions matrix of TPs from T 1 to T 2 : if RR or NS was observed in T 1 and DD was observed in T 2 , suggesting significant decreased vegetation productivity in T 2 compared with that in T 1 .This definition emphasized the changes of TPs in two consecutive periods, and this conversion was attributed to human activities.Similarly, vegetation restoration (VR) was redefined for the cases when DD or NS was present in T 1 and RR observed in T 2 , indicating an improvement of vegetation productivity in T 2 .
The time-lag effect may influence patterns of vegetation change, and the characterization of neighborhood interactions was computed by Verburg's F index [46], which measures the localized impact of each pattern (RR, DD, or NS) on vegetation degradation in a local region relative to that in a whole study area.As there were three involved patterns, namely RR, DD, and NS, which were derived from the vegetation dynamics based on the historical NDVI and precipitation data, three F i,k,d variables corresponding to pattern k (k = RR, DD, or NS) at each pixel (location) i within a local window defined by a predefined size d were computed for period T 1 .The average neighborhood characteristic for a particular pattern k can be calculated by taking the average F for all pixels.We were interested in testing whether the neighborhood interactions at an earlier period (T 1 ) given by the three F indices helped to predict vegetation degradation.The parameter d was relatively difficult to determine.Verburg et al. [46] found that the average neighborhood characteristic for a land use type changed from other land use types, and varied with different window sizes, but usually leveled off after a threshold.We found that the average neighborhood characteristic for DD developed from RR and NS showed a negligible change to window size when it was beyond seven units.Considering this result, as well as the previous study on vegetation density interpolation [18], d was set to seven (equivalent to a window with 1.75 km in width).
Binary logistic regression (BLR) is a bivariate statistical analysis model that predicts the probability of a binary response, for example, the presence or absence of an event.The BLR model provides insights into the relationship between a dependent variable, and one or a set of independent explanatory variables.Formally, the probability of occurring vegetation degradation (y) is a function for a vector X(x 1 , x 2 , . . ., x n ): where y is a binary variable that takes 1 as meaning an occurrence of vegetation degradation, or 0 as meaning the opposite; P( y = 1|X) denotes the probability that the event occurred at a certain location; element x i (i = 1, 2, . . ., n) in X corresponds variable i explaining y, P( y = 1|X)/1 − P( y = 1|X) is the odds ratio, and ln(P( y = 1|X)/(1 − P( y = 1|X))) is a logit transformation.β 0 and β i (i = 1, 2, . . ., n) are the regression coefficients of the function, which are best fit by a maximum likelihood estimation algorithm.Climate impact was excluded from being a contributing factor, since the vegetation degradation was derived by the RESTREND-NDVI method.Therefore, explanatory variables from human activities, including RD and GI, were taken in the BLR analysis, while the post-processed vegetation degradation (VD) was taken as the dependent variable.The RD and GI in T 2 were temporally averaged and input to the BLR models.In addition, neighborhood interactions (NIs) are also possible factors.It is of great value to know whether including NIs in the patterns of vegetation change improves the performance of vegetation degradation modeling.To examine the contribution of NIs, two BLR models-BLR-NIs that considered NIs and the other that did not (simply BLR)-were built and compared.All of the explanatory variables were standardized into zero-mean, following the formula: where x * i denotes the standardized value of variable x for sample/location i, x i is the original value of the variable, x is the mean value, and σ is the standard deviation of the original variable.After substituting x i with x * i in Equation ( 1), the result of the regression can be expressed in terms of conditional probability at any spatial location (pixel i) to be predicted: The probability of vegetation degradation (VD) at a location i, indicated by p i , in the study region could be calculated according to Equation (3).It is a straightforward computation to make a probability map of vegetation degradation over a given region when the required data for the models (BLR and BLR-NIs) are available.
A total of 6.5% of pixels labeled with VD and non-VD were randomly selected from the processed dataset and split into two parts with equal size (5000 samples) to fit and test the BLR and BLR-NIs, respectively.A confusion matrix and relative operating characteristic (ROC) were adopted to evaluate the performance of the models [53].The confusion matrix method uses a contingency table (also known as a cross-tabulation) to display a multivariate frequency distribution of the variables, and evaluate the classification accuracy of the models.The ROC method is based on a curve relating the true-positive rate and the false-positive rate over a range of cutoff values in classifying the probability.The area under the curve (AUC) is used to evaluate the accuracy of the model.The AUC varies between 0.5 (completely random) and 1 (perfect discrimination), and values above 0.7 are generally assumed to indicate a good model fit [54].All of the statistics were performed using the Statistical Package for the Social Sciences (SPSS, v22).

Spatio-Temporal Patterns in Vegetation Change
While a great part of area showed no obvious change (NS) in vegetation cover for both T 1 and T 2 , a significant difference was observed in the spatial distribution of RR and DD between period T 1 and T 2 (Figure 2).The DD during T 1 was mostly observed in the northwestern part (covering 16.75% of the total area), and RR mainly dominated the southeastern part, covering 12.80% of the total study area.Conversely, during T 2 , DD mainly appeared in southeastern region, which only took 8.64% of the total area, but RR shifted to the northwestern part, which showed a much larger area (36.29%) when compared with that in its previous period.Furthermore, a great part of the area in the southwest was also categorized as RR during T 2 .
Sustainability 2018, 10, x FOR PEER REVIEW 6 of 15 A total of 6.5% of pixels labeled with VD and non-VD were randomly selected from the processed dataset and split into two parts with equal size (5000 samples) to fit and test the BLR and BLR-NIs, respectively.A confusion matrix and relative operating characteristic (ROC) were adopted to evaluate the performance of the models [53].The confusion matrix method uses a contingency table (also known as a cross-tabulation) to display a multivariate frequency distribution of the variables, and evaluate the classification accuracy of the models.The ROC method is based on a curve relating the true-positive rate and the false-positive rate over a range of cutoff values in classifying the probability.The area under the curve (AUC) is used to evaluate the accuracy of the model.The AUC varies between 0.5 (completely random) and 1 (perfect discrimination), and values above 0.7 are generally assumed to indicate a good model fit [54].All of the statistics were performed using the Statistical Package for the Social Sciences (SPSS, v22).

Spatio-Temporal Patterns in Vegetation Change
While a great part of area showed no obvious change (NS) in vegetation cover for both T1 and T2, a significant difference was observed in the spatial distribution of RR and DD between period T1 and T2 (Figure 2).The DD during T1 was mostly observed in the northwestern part (covering 16.75% of the total area), and RR mainly dominated the southeastern part, covering 12.80% of the total study area.Conversely, during T2, DD mainly appeared in southeastern region, which only took 8.64% of the total area, but RR shifted to the northwestern part, which showed a much larger area (36.29%) when compared with that in its previous period.Furthermore, a great part of the area in the southwest was also categorized as RR during T2.The patterns of vegetation change present in each of the two periods (T1 and T2) exhibited much spatial variation.After the patterns in T1 and T2 were overlaid, several combinations were obtained.There was a clear spatial shift in the distribution of the patterns between the two study periods.For example, locations that showed stable (NS) or even recovered (RR) in vegetation productivity at an earlier stage (T1), but significantly declined (denoted as DD) at a later stage (T2), were mostly located in the southeastern part of the study area (Figure 2).Conversely, most of the degraded (DD) or no significant (NS) change in vegetation cover at T1 appeared to recover (RR) at a later stage (Figure 2).This spatial shift in the distribution of the patterns in vegetation change indicated a revertible trend in vegetation changes.The result from the overlaying the layers of the patterns of vegetation change The patterns of vegetation change present in each of the two periods (T 1 and T 2 ) exhibited much spatial variation.After the patterns in T 1 and T 2 were overlaid, several combinations were obtained.There was a clear spatial shift in the distribution of the patterns between the two study periods.For example, locations that showed stable (NS) or even recovered (RR) in vegetation productivity at an earlier stage (T 1 ), but significantly declined (denoted as DD) at a later stage (T 2 ), were mostly located in the southeastern part of the study area (Figure 2).Conversely, most of the degraded (DD) or no significant (NS) change in vegetation cover at T 1 appeared to recover (RR) at a later stage (Figure 2).This spatial shift in the distribution of the patterns in vegetation change indicated a revertible trend in vegetation changes.The result from the overlaying the layers of the patterns of vegetation change indicated a relatively small percentage of vegetation degradation compared to that of vegetation recovery (7.71% for VD vs. 34.78%for VR, in Figure 3).indicated a relatively small percentage of vegetation degradation compared to that of vegetation recovery (7.71% for VD vs. 34.78%for VR, in Figure 3).

Comparison of Vegetation Degradation by the BLR Models
The modeling results for vegetation degradation by the BLR-with and without NIs-are listed in Table 1.The BLR model demonstrated that grazing intensity and road density were significant and positive contributors to the human-induced vegetation degradation.For grazing intensity, the odds ratio Exp(β) could be interpreted as an increase of grazing intensity by one standardized unit almost doubling the possibility of human-induced degradation (184.1%).Road density exhibited an even more significant impact on the predicted probability of grassland degradation.When road density increased by one standardized unit, the probability for vegetation degradation increased more than six times (i.e., 6.339).
When the localized association variables, or neighborhood interactions (NIs) in TPs, were added to the BLR model (BLR-NIs), some changes were observed.First, grazing intensity was excluded from the model, demonstrating that the impact of grazing intensity was diminished (from significant in the BLR to insignificant) with the inclusion of NIs.Road density in BLR-NIs was still significantly and positively correlated with grassland degradation, with the coefficient (β) reduced from 1.847 to 1.671, suggesting that a one unit increase in standardized road density would lead to a high probability (increased more than five times) of vegetation degradation.Though the neighborhood interactions (DD, RR, and NS) in TPs demonstrated their impact on vegetation degradation, they presented different roles.For example, a higher F value in RR (FRR) would lead to a

Comparison of Vegetation Degradation by the BLR Models
The modeling results for vegetation degradation by the BLR-with and without NIs-are listed in Table 1.The BLR model demonstrated that grazing intensity and road density were significant and positive contributors to the human-induced vegetation degradation.For grazing intensity, the odds ratio Exp(β) could be interpreted as an increase of grazing intensity by one standardized unit almost doubling the possibility of human-induced degradation (184.1%).Road density exhibited an even more significant impact on the predicted probability of grassland degradation.When road density increased by one standardized unit, the probability for vegetation degradation increased more than six times (i.e., 6.339).
When the localized association variables, or neighborhood interactions (NIs) in TPs, were added to the BLR model (BLR-NIs), some changes were observed.First, grazing intensity was excluded from the model, demonstrating that the impact of grazing intensity was diminished (from significant in the BLR to insignificant) with the inclusion of NIs.Road density in BLR-NIs was still significantly and positively correlated with grassland degradation, with the coefficient (β) reduced from 1.847 to 1.671, suggesting that a one unit increase in standardized road density would lead to a high probability (increased more than five times) of vegetation degradation.Though the neighborhood interactions (DD, RR, and NS) in TPs demonstrated their impact on vegetation degradation, they presented different roles.For example, a higher F value in RR (F RR ) would lead to a greater likelihood for vegetation degradation, meaning that the locations with RR at T 1 were much more likely to experience vegetation degradation in the next period.The VD probability with a unit increase in standardized F RR would increase more than eight times (Exp(β) = 8.812).Conversely, a higher F value in DD (F DD ) would cut the probability in vegetation degradation in half (Exp(β) = 0.473) with a standardized unit increase.Lastly, the F value for NS (F NS ) did not show a significant impact on the vegetation degradation.

Accuracy Evaluation
Table 2 presented the accuracy of the two models, comparing the results of the observed and predicted VD.The percentage correct for the positive case (VD) was 68.0% for the BLR model and 77.6% for the BLR-NIs model, while that for negative case (non-VD) was 91.8% and 93.5% for the BLR model and the BLR-NIs model, respectively, suggesting that both models had better prediction power for the negative cases.The overall accuracy for the BLR model and the BLR-NIs model was 79.9% and 85.5%, respectively.From the performance of both the percentage correct and overall accuracy, the BLR-NIs model showed improved power in predicting VD over the BLR model, indicating the importance of the inclusion of localized spatial associations in TPs.The ROC curves for the two BLR models were compared.Figure 4 and Table 3 show that the AUC for the BLR and BLR-NIs models was 0.850 and 0.919, confirming that the model BLR-NIs had a better accuracy than the BLR model in predicting the probability of vegetation degradation.Finally, pseudo Cox & Snell and Nagelkerke R 2 , both computed as the ratio of the variance explained by independent variables to the total of the dependent, but from different aspects, for the two models are summarized in Table 3.For both R 2 , the values for the BLR-NIs model were much higher than those for the BLR model, suggesting that the variables taken in the BLR-NIs model have more explanatory power.

Discussion
This paper tried to analyze human-induced grassland degradation through residual trend analysis on the normalized difference vegetation index (RESTREND-NDVI), which is often adopted for arid or semi-arid regions, and to quantify the degradation with two typical factors from human activities.Transportation (road density) and grazing (intensity) were selected in the study because they had been widely recognized as driving factors of vegetation productivity.Our emphasis on human-induced vegetation degradation was based on the consideration that they were relatively easier to be reduced through designing sound grassland management plans.Nevertheless, very few models accounted for the neighborhood interactions (NIs) in the patterns of vegetation change when studying grassland vegetation degradation.Here, we compared two binary logistic regression (BLR) models to explain the vegetation degradation.While both of them included the two human-induced factors, our initial hypothesis was that NIs could be important predictors, and might be able improve model performance.
Overgrazing has long been regarded as the principal factor driving vegetation degradation [8].Previous studies also indicated that overgrazing has been present in the Xilin river basin for a long time, triggering a series of ecological problems, including serious grassland degradation [23,48,55].This was also confirmed in the BLR model, which presented a significant and positive coefficient in the fit BLR function, suggesting that higher grazing intensity led to a greater likeliness of vegetation degradation.Under the influence of overgrazing, the structure, productivity, or biological characteristics of grassland vegetation could be affected [9,10].Furthermore, overgrazing on vegetation often resulted in the degraded succession of vegetation communities, causing degradation in vegetation productivity in the end [56].However, it is interesting that the grazing intensity was not listed as a significant contributor to vegetation degradation once the characteristics

Discussion
This paper tried to analyze human-induced grassland degradation through residual trend analysis on the normalized difference vegetation index (RESTREND-NDVI), which is often adopted for arid or semi-arid regions, and to quantify the degradation with two typical factors from human activities.Transportation (road density) and grazing (intensity) were selected in the study because they had been widely recognized as driving factors of vegetation productivity.Our emphasis on human-induced vegetation degradation was based on the consideration that they were relatively easier to be reduced through designing sound grassland management plans.Nevertheless, very few models accounted for the neighborhood interactions (NIs) in the patterns of vegetation change when studying grassland vegetation degradation.Here, we compared two binary logistic regression (BLR) models to explain the vegetation degradation.While both of them included the two human-induced factors, our initial hypothesis was that NIs could be important predictors, and might be able improve model performance.
Overgrazing has long been regarded as the principal factor driving vegetation degradation [8].Previous studies also indicated that overgrazing has been present in the Xilin river basin for a long time, triggering a series of ecological problems, including serious grassland degradation [23,48,55].This was also confirmed in the BLR model, which presented a significant and positive coefficient in the fit BLR function, suggesting that higher grazing intensity led to a greater likeliness of vegetation degradation.Under the influence of overgrazing, the structure, productivity, or biological characteristics of grassland vegetation could be affected [9,10].Furthermore, overgrazing on vegetation often resulted in the degraded succession of vegetation communities, causing degradation in vegetation productivity in the end [56].However, it is interesting that the grazing intensity was not listed as a significant contributor to vegetation degradation once the characteristics of neighborhood interactions were accounted for in the modified BLR model, i.e., the BLR-NIs model.This result indicated that the grazing effect could have been included and reflected in the characterization of the neighborhood interactions, which was captured by the F indices in the BLR-NIs model.
Compared to other factors, road accessibility showed the strongest impact on vegetation degradation.The significance was confirmed by both of the BLR models.While constructing more road networks expanded the domain of human activities, it also might lead to more damage to the grassland vegetation due to the easier access to vegetated areas.The side impact on vegetation due to road networks on grasslands has been documented by a few studies investigated in Inner Mongolia and Mongolia [18,57].The paper [57] also analyzed a few reasons, including the hardness of the top layer and soil erosion, for the damage imposed by road tracks on grasslands.The BLR models in the current study suggest that a one-unit increase in the standardized road network density led to a five to six-fold increase in the probability of vegetation degradation.Traditionally, local people relied more on fixed locations for livestock grazing; in the last few decades, more and more road networks were built with the economic development of the region, which greatly facilitated moving grazing and the formation of more village centers in the countryside [18].The current study indicated that the grassland vegetation adjacent to the road network was more likely to suffer degradation.
Many factors, ranging from human activities to climate change, have proved to be important for grassland vegetation dynamics.However, previous studies rarely consider how neighborhood interactions impact vegetation degradation.In fact, vegetation degradation is not only affected by exterior factors (e.g., precipitation or grazing), it also varies under different spatial and temporal contexts of the patterns of vegetation change.For example, in our previous study, from the perspective of vegetation succession patterns, we found that grassland degradation was closely related to the background vegetation types from which the successions developed [18].This study revealed that vegetation degradation was more likely to occur at places where vegetation showed a recovery pattern (denoted as RR) at an earlier time, but was less likely to happen if vegetation had already showed degradation (labeled as a DD pattern) previously in the BLR-NIs model.This finding suggests that the patterns of vegetation change in the study region tended to reverse/convert over time, making it necessary to continue to enforce vegetation protection policies in the already recovered vegetated area, while at the same time strengthening efforts to improve other degraded areas.Out of those policies, an adaptive grazing strategy has been introduced in the past decade to allow grassland pasture to be fed by livestock in a rotational manner, which prohibits continuous grazing at the same location, and instead encourages feeding livestock at different locations based on the availability of pasture.Adaptive grazing tends to balance grazing intensity over an entire region by interrupting grazing always in degraded (DD) areas, and looking for vegetation recovered (RR) areas or areas with no significant (NS) change in vegetation.Other enforced strategies such as fenced grazing also prevented grassland vegetation from severe degradation [4].Similar outcomes resulted from the enforcement of those policies (e.g., adaptive grazing) may suggest certain clues for the insignificance of the grazing effect once the neighborhood interactions in TPs were included in the model.
Regression analysis is a typical method for estimating the relationships among variables.Binary logistic regression analysis is based on a regression model particularly fit for investigating the dependence of a variable with two responses, one positive and one negative, on a set of explanatory variables.Understanding the probability of the occurrence (or non-occurrence) of vegetation degradation is very important for the sustainable development of the grassland ecosystem.It is reasonable to apply BLR-based models to investigate factors related to grassland degradation.With support from available statistical packages such as SPSS, BLR models could be quickly fit and evaluated.The two BLR models (BLR versus BLR-NIs) compared the results between including and excluding neighborhood interactions in the prediction of the probability of vegetation degradation.Compared with the BLR model, the BLR-NIs model had much improvement in the prediction power.The outcome from the BLR models, particularly the odd ratio (i.e., Exp(β)) provided a valuable reference for policymakers in regard to designing informed strategies for improving grassland management in the future.Along the lines of the regression model, non-stationary models such as geographically weighted regression may provide better explanation in exploring the relationships among the variables.However, the main purpose of the current work was to test the effectiveness of including NIs; as a result, the more general forms of regression models (i.e., BLR and BLR-NIs) were compared.Further studies that combine geographically weighted regression with binary logistic regression may be preferred in finding improved models for analyzing grassland degradation.
The current study confirmed the effectiveness of accounting for the neighborhood characteristics in predicting vegetation degradation in BLR-NIs.Nevertheless, a few limitations should be acknowledged.First, RESTREND-NDVI was recognized as being able to separate the impacts on vegetation productivity from climatic and human-induced factors, but their interaction effect was not addressed.Therefore, the human-induced vegetation degradation patterns derived from RESTREND-NDVI may contain some interaction components of the principal climatic and human-induced factors.Moreover, certain uncertainties regarding the results from RESTREND-NDVI have also been pointed out [58].Those scenarios will inevitably lead to accumulated uncertainties in the current study.Second, more detail work should be input to measure the localized spatial associations or the neighborhood interactions in different TPs.For example, a spatially varied rather than fixed size in the definition of neighborhood might be necessary in modeling NIs, because measurements of neighborhood interactions are possibly constrained by the definition of a neighborhood matrix [59].Lastly, the current study only took into account two factors, namely GI and RD, from human activities.The performance of the models might be improved by the inclusion of other factors from human activities, such as the distribution of village centers.However, compiling the qualified dataset required more effort than expected, and future studies need to examine more potential variables, which may account for more variance in the vegetation degradation.

Conclusions
Historical remote sensing archives and climate data are valuable resources for mapping long-term grassland vegetation dynamics.Modeling the trajectory patterns of vegetation is important for both policymakers of grassland management and academic researchers.Vegetation dynamics are attributed to many factors, including climate change and human activities.While climatic constraints have proved to restrict vegetation productivity in arid or semi-arid grassland, human responses to lower those constraints are often limited, and sometimes impossible.Thus, many practical measures focusing on improving the grassland ecosystem turn to the optimized management of human-induced factors by balancing the requirement of grazing from agriculture production and the protection of grassland pastures at the same time, possibly using fenced grazing or adaptive grazing.Therefore, modeling the vegetation degradation from human-induced factors seems more promising.
Excluding the impact on vegetation dynamics from climate factors has been tested, particularly for arid and semi-arid grassland regions.In this study, we subtracted the impact of the climate factor by RESTREND-NDVI, and then focused on the human-induced factors that exhibited a significant impact on grassland degradation.However, modeling grassland degradation should also take into account the spatial and temporal context of the patterns of vegetation change themselves.Vegetation changes are not only related to the vegetation properties at the exact location, they are also related to the vegetation properties within a spatial neighborhood.Possible links between the patterns of vegetation change also exist from a temporal context, since vegetation cover at a later stage is developed on the basis of a previous one.To analyze such neighborhood interactions, this study considered three F indices (measurements of characteristics of neighborhood interactions), along with factors from human activities (i.e., grazing intensity and road density) in the BLR-NIs model.The BLR model indicated a positive and significant contribution to the grassland degradation from grazing intensity and the density of the road network.However, the variable of grazing intensity was excluded from the list of significant factors in the BLR-NIs model; instead, neighborhood interactions from DD (indicated by F DD ) and RR (F RR ) on the vegetation degradation were found to be significant.With the inclusion of NIs, the study showed that vegetation degradation was found to occur in a recovered vegetated area (high F for RR) during an early period, but less likely to happen in a previously degraded area (high F for DD), suggesting an unstable status in the vegetation dynamics from a long-term and temporal perspective.Thus, continuous effort must be taken to improve the grassland ecosystem.Our study confirmed that neighborhood interactions should be taken into consideration when modeling the patterns of vegetation change because, compared to the BLR, the BLR-NIs model had a significantly improved model prediction performance.The proposed method, i.e., integrating BLR and neighborhood interactions, can be a useful tool for modeling grassland vegetation dynamics.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.Location of the study area.

Figure 1 .
Figure 1.Location of the study area.

Figure 2 .
Figure 2. Trajectory patterns (DD, RR, and NS) of vegetation during T1 (a) and T2 (b), in the Xilin river basin.Non-vegetated (built-up, in white) areas were excluded from further analysis.DD: degradation, RR: recovery, NS: no significant change.

Figure 2 .
Figure 2. Trajectory patterns (DD, RR, and NS) of vegetation during T 1 (a) and T 2 (b), in the Xilin river basin.Non-vegetated (built-up, in white) areas were excluded from further analysis.DD: degradation, RR: recovery, NS: no significant change.

Figure 3 .
Figure 3. Grassland trajectory patterns (vegetation degradation denoted as VD and vegetation recovery denoted as VR), from a long-term perspective, computed from the patterns of vegetation change during two consecutive periods (T1 and T2) (see Figure 2).

Figure 3 .
Figure 3. Grassland trajectory patterns (vegetation degradation denoted as VD and vegetation recovery denoted as VR), from a long-term perspective, computed from the patterns of vegetation change during two consecutive periods (T 1 and T 2 ) (see Figure 2).

Figure 4 .
Figure 4. Performance comparison by relative operating characteristic (ROC) curves for the two models.

Figure 4 .
Figure 4. Performance comparison by relative operating characteristic (ROC) curves for the two models.

Acknowledgments:
Funding support for this study included the National Natural Science Foundation of China (Nos.41371371 and 41501441), Open Fund of Key Laboratory of Geographic Information Science (Ministry of Education), East China Normal University (No. KLGIS2017A05), Hubei Provincial Natural Science Foundation of China (No. ZRMS2017000737), and Large Scale Environment Remote Sensing Platform (Nos.16000009, 16000011 and 16000012), Wuhan University.The authors thank the anonymous reviewers for their constructive comments that have helped us to improve the manuscript.Author Contributions: Yuwei Wang, Zongyao Sha and Zhenyu Wang conceived and designed the experiments; Yuwei Wang performed the experiments; Ruren Li and Xiaoliang Meng analyzed the data; Xingjun Ju and Yuguo Zhao contributed materials and analysis tools; Yuwei Wang and Zongyao Sha wrote the paper together.

Table 1 .
Results of two binary logistic regression (BLR) models.

Table 2 .
Classification table for two BLR models.VD: vegetation degradation.

Table 3 .
Statistics of the pseudo R 2 and area under the curve (AUC) for BLR models with and without neighborhood interactions (NIs).

Table 3 .
Statistics of the pseudo R 2 and area under the curve (AUC) for BLR models with and without neighborhood interactions (NIs).
a particular period (e.g., T 1 in this study) DD vegetation degradation during a particular period (e.g., T 1 in this study) NS no significant change in vegetation during a particular period (e.g., T 1 in this study) F a measurement for neighborhood interactions proposed in Verburg's paper F DD , F RR , F NS F value computed for DD, RR and NS, respectively based on the conversions matrix of TPs from T 1 to T 2 when RR or NS was observed in T 1 and DD was observed in T 2 .Note that it is different from DD VR redefined vegetation restoration based on the conversions matrix of TPs from T 1 to T 2 when DD or NS was present in T 1 and RR observed in T 2 .Note that it is different from RR AUC area under the curve ROC relative operating characteristic