## 2 Materials and Methods
### Methods
#### 2.3.2 XGBoost-SHAP Model  
To develop an urban expansion model, this study combines the eXtreme Gradient Boosting (XGBoost) algorithm and Shapley Additive ExPlanations (SHAP).  
The eXtreme Gradient Boosting (XGBoost) algorithm is an open-source library that supports an efficient implementation of gradient-boosted decision trees [34]. Since the package includes efficient linear model solver and tree learning algorithms, users can compute much faster than other existing gradient boosting tools. For the analysis, we used the version 1.5.2.1 of \"xgboost\" package in R, which was released in February 2022.  
The gradient boosting decision tree (GBDT) is an ensemble learning technique that combines a series of weak decision trees to build a strong learner [35]. In this algorithm, each decision tree is trained from the residuals of the previous one, and iteratively constructs a more accurate model until the loss function is minimized. Due to its high predictive precision and ability to deal with both categorical and continuous variables, GBDT has been widely used in various fields of research [36].  
For a given training dataset \(\{(x_{i},y_{i})\}\overset{N}{i=1}\), the initialized model with constant value is defined as below:  
\[f_{0}(x)=argmin_{\gamma}\sum_{i=1}^{N}L(y_{i}\,,\,\gamma) \tag{1}\]  
where \(\gamma\) denotes the constant, \(L(y,F(x))\) denotes a differentiable loss function of \(\gamma\), and \(argmin_{\gamma}\) indicates the value \(\gamma\) that minimizes the function.  
For \(m\) number of iterations, the negative gradient of the loss function is calculated as  
\[g_{m}(x_{i})=-\bigg{[}\frac{\partial L(y_{i},f(x_{i}))}{\partial t(x_{i})} \bigg{]}_{f=f_{m-1}} \tag{2}\]  
Here, \(g_{m}(x_{i})\) is calculated by taking a derivative of previous loss function \(f_{m-1}(x)\).  
Figure 3: Research procedure.  
Then, a base learner (or weak learner) solves the optimization problem, as follows:  
\[\theta_{m}=argmin_{\theta}\sum_{i=1}^{N}L(y_{i},F_{m-1}(x_{i})+\theta t(x;\;\mu_{ m})) \tag{3}\]  
where \(L(y_{i},F_{m-1}(x_{i})+\theta t(x;\;\mu_{m}))\) indicates the loss function on each node \(i\).  
Lastly, the model is updated as  
\[f_{m}(x)=f_{m-1}(x)+\theta_{m}t(x;\;\mu_{m}) \tag{4}\]  
Here, \(t(x;\;\mu_{m})\) denotes the selected node and \(\theta_{m}\) denotes the learning rate.  
In the XGBoost model, there are several parameters that require to be designated to maximize performance of the model, while preventing overfitting problems [37]. More specifically, the model needs to select the suitable number of iterations, maximum depth, the fraction of observations, and learning rate. In addition, the parameters including 'colsample_bytree', 'alpha', and 'lambda' determine the weights and fitness of the model.  
The Shapley Additive exPlanations (SHAP) was first proposed by [PERSON] and [PERSON] [38], and has been used to evaluate the relative importance of features in machine learning models. In SHAP, the importance of each independent variable on the model outcome is calculated based on its marginal contribution [39]. For an XGBoost model of group N with n features, the SHAP value \(\phi_{i}\) assigned to each feature \(i\) is represented as  
\[\phi_{i}=\sum_{S\in N}\frac{|S|!(n-|S|-1)!}{n!}[v(S\cup\{i\})-v(S)] \tag{5}\]  
where \(S\) represents feature subsets derived from \(n\), and \(v\) represents the input features within the set \(S\).