The area of interest in this study was Batu City, East Java Province, Indonesia, located at 7°46′–7°55′S and 112°30′–112°39′E, with an area of ±199.09 km(Ali et al., 2024). There were 276 sample points in the area, consisting of 138 landslide (red dots) and 138 nonlandslide (green dots) points (Figure 2). This area was selected because of its high susceptibility to landslides. This susceptibility is mainly due to its physiographic conditions, as the area forms a topography with steep slopes due to the surrounding Mount Panderman, Mount Arjuno, and Mount Welirang. This topographic factor is also amplified by the elevation variation, which ranges from 680 to 1,200 m above sea level.

This study utilized 13 conditioning factors, all of which were accessed and processed within the GEE cloud computing platform to ensure reproducibility and computational efficiency. All predictor variables used, along with their sources and spatial resolutions, are summarized in Table 1. Some data are available in raster format at varying resolutions (90–500 m), whereas other data are vector-based (e.g., distance from rivers, geological formations, and distance from lineaments). All vector data were converted to raster format, and all layers were resampled to a uniform spatial resolution (10 m) using the GEE platform to ensure that all variables share the same coverage and projection system. Once all conditioning factor maps were ready in the raster format, they were then combined into a single multiband image. This combination is necessary because machine learning algorithms require all input variables (predictors) for each pixel in the form of a single image with multiple channels (bands).

Supervised machine learning models require labeled datasets for training. This study assigned two class labels to the samples: landslide points (1) and nonlandslide points (0). The distribution of landslide and nonlandslide sample points was sourced from the National Disaster Management Agency (BNPB) landslide potential inventory map, accessed through the InaRISK platform (Website link). The dataset used consisted of 276 sample points, comprising 138 landslide points and 138 nonlandslide points. Both classes had the same number of points to produce an initial balance in the training data. To avoid class bias and ensure sample representativeness, stratified random sampling was applied to the dataset (Fu et al.,2025; Zhou et al., 2024). This approach ensured balanced representation of both classes. Finally, the entire dataset was divided into two parts at a 70:30 ratio, with 70% used as training data and 30% as testing data.

Landslide susceptibility mapping in this study used 13 conditioning factors as predictor variables. In machine learning, they must be mutually independent and not strongly correlated (Qi et al., 2021). Therefore, a Pearson correlation matrix of landslide predictor variables was used before model formulation. Multicollinearity occurs when two variables in a multiple regression model are strongly correlated (Liu et al., 2022). This process is essential in modeling to ensure the stability and reliability of the predictive model being formulated. The correlation matrix was calculated using all predictor variables. Correlation coefficient values above 0.7 indicated high multicollinearity (Ullah et al., 2024).

The machine learning algorithms used in this study consisted of five models: RF, CART, GTB, SVM, and KNN. RF builds several decision trees on a random subset of the data and features and combines the voting results from these trees to produce a final decision ((Le et al., 2024); (Khan et al., 2021)). Meanwhile, CART splits data using the Gini index, forming an easy-to-understand tree structure. The resulting tree model consists of internal nodes representing separation criteria and leaf nodes that contain class predictions (Zhou et al., 2024). GTB operates through an iterative process; it builds a series of decision trees sequentially, with each new tree correcting the errors of the previous one (Tao et al., 2022). All trees are combined using weighted ensemble prediction to gradually improve accuracy ((Kavzoglu & Teke, 2022); (Li et al., 2025)). SVM identifies the optimal hyperplane that maximizes the margin between the landslide and nonlandslide classes using support vectors ((James et al., 2013); (Meena et al., 2025)). Finally, KNN predicts the value of a new observation based on the class or average of the nearest neighbors in feature space, usually calculated using Euclidean distance (Le et al., 2024).

The optimal configurations for all machine learning models resulting from hyperparameter optimization are shown in Table 2. The RF model was configured with 100 decision trees and a seed value of 0. It produced an out-of-bag error of 0.114, indicating a good level of internal accuracy. The CART model adopted an ensemble approach comprising 193 trees, with a maximum structural restriction of 15 nodes and a maximum depth of 8 to control model complexity. This configuration balanced model complexity and prediction accuracy. The GTB model achieved optimal performance with 200 trees, a learning rate (shrinkage) of 0.05, and a sampling rate of 0.7. This configuration helped improve model stability. The SVM model used a linear kernel with a gamma of 0.005 and a cost (regularization parameter) of 1. It created an optimal decision boundary for class separation. Finally, the KNN model achieved the best performance with the number of neighbors (k) set to 10 and using the Manhattan metric. It was proven effective at handling high-dimensional data with 193 features. These optimized configurations ensured that each model could optimally extract patterns from the landslide distribution dataset.

The evaluation matrix used to evaluate model performance includes accuracy, precision, recall (sensitivity), F1-score, kappa coefficient, and AUC. Figure 3 comprehensively compares the model predictions and actual values. The main components are true positive (TP), true negative (TN), false positive (FP), and false negative (FN). TP represents the correctly predicted positive cases (landslide points), whereas TN represents correctly predicted negative cases (nonlandslide points). Meanwhile, FP represents negative cases (nonlandslide points) that are predicted as positive (landslide points), whereas FN represents positive cases (landslide points) predicted as negative (nonlandslide points).

The performances of the models were validated using various evaluation metrics: accuracy, precision, recall (sensitivity), and F1-score, as formulated in Equations (1)–(4) ((Khadka et al., 2025); (Qiu et al., 2024)).

(1)

(2)

(3)

(4)

In addition, the performances of the models obtained were validated using the kappa coefficient, as shown in Equation (5) (Khadka et al., 2025). This metric can provide a more reliable description of the model’s predictive performance, especially for unbalanced data. In some cases, it is considered the standard metric for model comparison.

(5)

Where

(6)

And

(7)

This metric is considered more reliable than accuracy because it accounts for agreement occurring by chance. Model performance is generally described on a scale from 0 to 1, as shown in Table 3(Landis & Koch, 1977).

Finally, model performance was also validated using ROC curves. The model with the highest AUC value and optimum balance between precision and recall was considered the most representative model. Spatial validation was applied using spatial cross-validation techniques to reduce the effects of spatial autocorrelation (spatial overfitting) and improve prediction reliability ((Le et al., 2024); (Tjahjono et al., 2024)).

Spatially, the maps of the 13 predictor variables were generated in QGIS 3.34 using the UTM 49S projection, which was appropriate for the study area. The type of spatial data and the grid resolution produced varied according to the data available on the GEE platform. They showed consistent geographic coverage, confirming the spatial alignment among predictor layers—a key precondition for reliable pixel-based modeling. Each map represented a predictor, with the distribution of values visualized through color gradients. This visual inventory facilitated a qualitative assessment of spatial heterogeneity across the study area. In addition, it enabled the identification of potential multicollinearity among predictor variables before performing quantitative hyperparameter optimization. The 13 predictor maps produced in the study are shown in Figure 4.

The significance of the correlations between the landslide predictor variables, as quantified using the Pearson correlation matrix, is presented in Figure 5. The elevation variable exhibited a strong positive correlation with annual rainfall (r = 0.88). Meanwhile, the strong correlation between elevation and annual rainfall exceeded the specified threshold ((Ullah et al., 2024); (Qi et al., 2021)). Therefore, to avoid multicollinearity problems due to strong correlation, one of them was eliminated. This study retained the elevation variable and eliminated the rainfall variable from the model. The elevation variable also showed a strong positive correlation with distance from rivers (r = 0.64) and distance from roads (r = 0.60). These suggest that areas with high elevation tend to experience greater precipitation intensity and are often spatially proximal to drainage systems and road infrastructure.

In contrast, geological formation with distance from rivers (r = −0.00) and slope (r = 0.03) showed very weak correlations. In addition, soil type and landforms exhibited a very weak correlation with most other predictors, as reflected in the correlation values of soil type with landforms (r = 0.01) and geological formation (r = −0.01), and landforms with distance from roads (r = 0.01). A similar pattern was also observed for NDVI and aspect, which showed no significant relationship with many other variables. Furthermore, significant inverse relationships were observed. The negative correlation between LULC and distance from roads (r = −0.59) suggests that certain land cover types tend to occur farther from roads. Similarly, the negative correlation between TWI and slope (r = −0.58) indicates that steeper slopes consistently correspond to lower TWI values. By eliminating the annual rainfall variable, the correlations obtained between the predictor variables were less than the threshold of 0.7, confirming the multivariable interdependence in the slope system. This was done to emphasize the importance of selecting predictor variables to avoid multicollinearity in landslide susceptibility mapping.

Landslide susceptibility mapping was conducted using the results of a Pearson correlation matrix analysis, where annual rainfall data were excluded to avoid multicollinearity. This resulted in 12 predictors used for mapping: elevation, slope, aspect, landforms, distance from rivers, TWI, soil types, LULC, NDVI, distance from roads, geological formations, and distance from lineaments. The maps generated by each algorithm are illustrated in Figure 6, showing varying distributions of susceptibility classes. The maps produced by RF and GTB exhibit a more concentrated distribution of areas with high and very high susceptibility. This indicates better consistency in susceptibility zone identification compared with the other models. In contrast, the maps produced by CART and KNN tend to display more scattered and fragmented patterns, with high susceptibility spread across several locations. In addition, SVM shows a relatively moderate distribution pattern compared with the other models. These pattern differences indicate variations in the algorithms’ ability to process nonlinear relationships between variables. The spatial distribution patterns on the landslide susceptibility maps indicate that areas with high and very high susceptibility classes in Batu City are concentrated in zones with steep slopes. These high-risk areas dominate the complex topography formed by the presence of young Quaternary volcanoes—namely, Mount Panderman, Mount Arjuno, and Mount Welirang.

The percentage areas of the susceptibility classes generated by the five models are presented in Figure 7. Based on the total area of Batu City (199.09 km²), the RF model showed the most balanced distribution, with significant portions in the very high class (60.70 km² or 30.49%) and very low class (70.16 km² or 35.24%), indicating its good ability to differentiate between each class. In contrast, CART and GTB showed strong polarization trends, with a dominance of the very high class of 89.94 km² (45.18%) and 83.74 km² (42.06%), respectively, and a lack of representation of the medium class, indicating their tendency to overpredict the very low or very high classes. SVM produced a relatively even distribution across all susceptibility classes, with peaks in the very low (80.38 km² or 40.37%) and very high (49.56 km² or 24.89%) classes, reflecting moderate characteristics in susceptibility classification. Furthermore, KNN showed the highest concentration in the very low (59.61 km² or 29.94%) and low (45.99 km² or 23.10%) classes, indicating that it was more conservative than the other models in predicting high-risk zones. These results show significant variations in classification across models. According to the distribution of landslide susceptibility classes, the RF model presents the best and most realistic distribution. Therefore, the selection of this model must consider the distribution of areas within each class to align with the decision-making context for risk-based mitigation.

The results for 83 testing datasets using the five algorithm models are shown in a confusion matrix (Table 4). These results indicate that RF achieved the best performance, with TP and TN of 35 and 43, respectively. This model also had the lowest total classification errors (FP and FN) of 5, followed by GTB with a total error of 7. In addition, CART and SVM tended to have higher FN—8 and 9 points, respectively. These indicate the model’s inability to detect positive instances. KNN had a high TN (43 points) but the highest FP (10 points). These indicate the error rate in predicting negative instances. Thus, RF and GTB are more reliable and balanced approaches for the tested datasets.

Based on the comprehensive analysis using 83 testing data points (representing 30% of the total data), RF demonstrated the best performance compared with the other four models. This model produced the highest values for accuracy (0.940), precision (0.946), recall (0.921), F1-score (0.933), and kappa coefficient (0.878), as shown in Table 5. Therefore, it is the most reliable and robust at identifying patterns in the testing data among all studied models (Melati et al., 2024). Its high accuracy and precision indicate that it can make accurate predictions overall and minimize false-positive classification errors. In addition, it produces high recall, indicating its ability to detect almost all truly positive instances (minimizing FNs). Moreover, the high kappa coefficient confirms that the level of agreement between the model’s predictions and the actual conditions falls into the almost perfect agreement category, exceeding what would be expected by chance (Landis & Koch, 1977).

The second-best model after RF was GTB. It performed strongly and was a balanced and effective predictor. In contrast, KNN, despite having a high recall of 0.900, had a relatively low precision of 0.730. This high FP rate makes its positive predictions less reliable. These results indicate its weakness: it tends to be overly aggressive in assigning positive labels, resulting in many FPs. Meanwhile, CART and SVM showed moderate capabilities but were less consistent than the other models. Therefore, apart from RF, the advantage of GTB over CART, SVM, and KNN lies in its ability to maintain the right balance between accuracy and reliability.

The ROC curves for landslide susceptibility prediction for the five models are shown in Figure 8. The results show that RF achieved the highest AUC of 0.976. It is regarded as the best classifier for mapping landslide susceptibility, particularly in the study area. The AUC value quantifies the probability that a model can randomly rank positive instances with higher scores than negative instances. Thus, a value of 0.976 indicates an excellent separation ability of 97.6%. The level of discrimination ensures that the resulting landslide susceptibility map can effectively differentiate between high- and low-susceptibility zones, providing a credible basis for field mitigation planning. The discriminatory performance of RF also exceeded that of the other models, such as GTB (0.962), CART (0.930), SVM (0.927), and KNN (0.900). This excellent AUC performance reinforces the conclusions drawn from the previous evaluation metrics, providing comprehensive confirmation that RF is the most optimized and consistent model in the overall performance assessment. Thus, the landslide susceptibility map and variable importance results of RF are selected for subsequent analyses.

The results of the variable importance analysis of RF in landslide susceptibility mapping are shown in Figure 9. Based on the variable importance, the total accumulated mean decrease in impurity value for all variables was 285.02. The slope variable exhibited the most significant contribution, with a mean decrease in impurity value of 43.77, equivalent to 15.36% of the total importance. This result indicates that slope is the main controlling factor of landslides, as it can affect slope stability. This is because it can affect shear stress and soil mass stability (Zeng et al., 2024). The mean decrease in impurity measures the total reduction in impurity achieved by each variable across the entire ensemble tree. A higher value generally indicates stronger predictive power (Breiman et al., 1984).

In addition, the other variables that contributed to landslide susceptibility after the slope variable were distance from lineaments (11.12%), NDVI (10.96%), distance from roads (10.42%), elevation (9.95%), TWI (9.79%), aspect (9.78%), and distance from rivers (9.43%). They showed the complexity of factors influencing landslide susceptibility in the study area. The influence of the lineament variable, which ranked second-highest, indicated the role of geological structures in underground water flow control and the formation of weak zones in rocks. In contrast, the NDVI variable reflected vegetation’s role in stabilizing slopes through root systems and rainfall interception ((Pappaka et al., 2025); (Hussain et al., 2022)). The influence of the distance from roads variable indicated that human activity, infrastructure, and vehicles can destabilize slopes by increasing loads in landslide-prone areas ((Ersayin & Uzun, 2025); (Qi et al., 2024)). TWI and distance from rivers both indicated the influence of hydrological conditions on pore water pressure, which is the main trigger of soil mass movement ((Nurwatik et al., 2022); (Xu et al., 2022)).

Finally, the remaining variables exhibited lower importance values, including LULC (5.69%), geological formation (3.59%), landforms (3.29%), and soil type (0.62%). These results also indicate that morphometric variables, including slope, TWI, aspect, and proximity features (distance from lineaments, roads, and rivers), have stronger predictive capabilities than inherent material characteristics (geological formation and soil type) and landform classification. This is probably due to the homogeneity of geological and soil characteristics in the study area and the low data resolution, which is insufficient to capture variations in material properties at this mapping scale. Furthermore, variables such as geological formation and soil type may be more static and less spatial than morphometric parameters that directly reflect slope stability conditions. These results indicate that the landslide susceptibility modeling approach in the study area is more effective when it focuses on topographic characteristics and proximity factors that are more responsive to spatial changes in slope conditions.

Based on model performance evaluation metrics, RF achieved the best overall performance for landslide susceptibility in Batu City. The advantage of RF’s ensemble nature is that it can build multiple decision trees from randomized subsets of data and features, thereby reducing the risk of overfitting and improving prediction stability (Liu et al., 2022). Furthermore, its ability to handle complex (nonlinear) relationships between landslide variables and its strong discriminatory power result in a balanced and proportional distribution of each susceptibility class. This aligns with previous studies confirming that RF outperforms other algorithms for landslide susceptibility mapping in tropical volcanic regions, as evidenced by its high AUC ((Melati et al., 2024); (Wu et al., 2022)). However, some previous studies reported that GTB was slightly superior to RF (Kavzoglu & Teke, 2022). This contradiction is likely due to differences in the study area’s characteristics, which featured steeper slopes and more heterogeneous land use. The RF bagging mechanism provides better stability against extreme gradient values, whereas the boosting approach in the GTB tends to amplify errors from outliers (Liu et al., 2022). Thus, tree-based ensemble algorithms, particularly RF, achieve excellent results with hyperparameter optimization, further strengthening the justification for the model selection in this study.

The variable importance analysis identified the most influential variables for landslide potential. In particular, the slope angle directly affects shear stress and soil stability (Zeng et al., 2024). Consequently, steep slopes amplify gravitational forces, which can trigger soil movement. The volcanic material in this area is highly weathered and poorly consolidated, making the slope composition vulnerable. If these conditions are present on a steep slope, the shear stress will naturally approach or exceed the soil’s shear strength (Muhardi & Perdhana, 2025). This situation becomes even riskier when triggered by external factors, such as rain or anthropogenic activities (Hasan et al., 2024). This explains why the slope variable had the highest importance compared with the other variables (Ersayin & Uzun, 2025). Therefore, mitigation efforts can be implemented by focusing on slope stabilization and land use in topographic zones with extreme slopes.

Although slope is often identified as a dominant factor in various landslide susceptibility studies, the dominance of a factor is highly determined by local geological, climatic, and anthropogenic characteristics ((Pasang & Kubíček, 2020); (Arabameri et al., 2020)). In areas with homogeneous volcanic material, structural parameters (e.g., distance to lineaments and NDVI) are more influential than soil type or geological formation. In contrast, in areas with high human activity, distance from roads and land cover become the primary factors, even surpassing slope (Pasang & Kubíček, 2020). In tectonically active regions, distance from roads and NDVI dominate, whereas elevation and distance from rivers play a greater role in reservoir areas because of periodic water level fluctuations ((Qi et al., 2024); (Zhou et al., 2024)). Furthermore, optimizing environmental factors reportedly improves modeling accuracy—that is, elevation and lithology can exhibit influence comparable to that of slope (Li et al., 2025). Thus, although slope is a dominant factor, it is not absolute and depends on environmental conditions, triggering mechanisms, and the quality of the data used.

The landslide susceptibility distribution pattern in the study area shows that high and very high classes are generally found on steep mountain slopes. This is consistent with the variable importance analysis, which identified slope as the dominant factor. This is because areas with high susceptibility are generally formed by the presence of a series of young Quaternary volcanoes, such as Mount Panderman, Mount Arjuno, and Mount Welirang. They can produce volcanic materials, such as ash, lapilli, breccia, and tuff, which are not yet properly consolidated (Noviyanto et al., 2020). These young volcanic materials have high porosity, facilitating rainwater absorption, which then increases pore water pressure and soil mass weight. These conditions produce loose, easily eroded materials susceptible to weathering, reducing slope stability in steep areas (Susilo et al., 2023). The combination of loose material properties and steep slope accelerates destabilization in the study area.

Classification on the GEE platform is inherently limited to pixel-level predictions and is unable to capture more detailed spatial information (e.g., slope connectivity or landslide paths). Consequently, pixel-based models, including those used in this study, tend to overlook the connectivity between upper and lower slopes. This limitation significantly affects the results, as failures in the upper slope often affect lower-slope areas (Wu et al., 2022). To address this, future studies should use slope units as the mapping unit instead of grid cells and include spatial validation. From a practical perspective, the susceptibility map generated in this study identifies an area of 70.16 km² (35.24%) as a very-high-susceptibility zone. These results can serve as a basis for disaster mitigation strategies, including slope stabilization, land-use zoning, and the development of early warning systems, in line with recommendations to integrate landslide susceptibility mapping into regional spatial planning (Tjahjono et al., 2024). However, the application of this model to other  volcanic regions remains dependent on recalibrating the relevance of predictors and the quality of local data, as machine learning performance is highly sensitive to spatial resolution and the completeness of the data inventory ((Hussain et al., 2025); (Zhou et al., 2024)). Furthermore, the landslide inventory in this study still relied on historical point data from the BNPB’s InaRISK platform. Therefore, direct field surveys are essential to ensure that landslide sample points are in areas that have experienced landslide events and to confirm and assess the potential for future landslides.

To improve data accuracy in landslide area inventories, future studies should integrate high-resolution datasets (e.g., drone-based digital elevation model mapping, direct field soil sampling, and laboratory soil testing) with dynamic factors, including multitemporal rainfall, land cover changes, and InSAR techniques. Additionally, this study identified slope, lineament density, and NDVI as dominant factors influencing landslides, enabling the implementation of phased, spatially tailored interventions. For example, areas with high lineament density require detailed field geological mapping to identify weak zones, whereas areas with low NDVI should be a priority for large-scale reforestation. Furthermore, susceptibility maps can be integrated into rainfall-based landslide early warning systems by setting rainfall thresholds in the study area based on landslide history. Finally, as a disaster mitigation strategy, these results can guide infrastructure zoning by avoiding the construction of new roads or the expansion of residential settlements in high-risk areas, thereby reducing future risks.