Uncertainty Quantification in Shear Wave Velocity Predictions: Integrating Explainable Machine Learning and Bayesian Inference

Abstract

The accurate prediction of shear wave velocity (Vs) is critical for earthquake engineering applications. However, the prediction is inevitably influenced by geotechnical variability and various sources of uncertainty. This paper investigates the effectiveness of integrating explainable machine learning (ML) model and Bayesian generalized linear model (GLM) to enhance both predictive accuracy and uncertainty quantification in Vs prediction. The study utilizes an Extreme Gradient Boosting (XGBoost) algorithm coupled with Shapley Additive Explanations (SHAPs) and partial dependency analysis to identify key geotechnical parameters influencing Vs predictions. Additionally, a Bayesian GLM is developed to explicitly account for uncertainties arising from geotechnical variability. The effectiveness and predictive performance of the proposed models were validated through comparison with real case scenarios. The results highlight the unique advantages of each model. The XGBoost model demonstrates good predictive performance, achieving high coefficient of determination ($R^2$), index of agreement (IA), Kling–Gupta efficiency (KGE) values, and low error values while effectively explaining the impact of input parameters on Vs. In contrast, the Bayesian GLM provides probabilistic predictions with 95% credible intervals, capturing the uncertainty associated with the predictions. The integration of these two approaches creates a comprehensive framework that combines the strengths of high-accuracy ML predictions with the uncertainty quantification of Bayesian inference. This hybrid methodology offers a powerful and interpretable tool for Vs prediction, providing engineers with the confidence to make informed decisions.

1. Introduction

Shear wave velocity (Vs) is a critical parameter in geotechnical and earthquake engineering. It has the greatest effect on the determination of dynamic behavior of soil and rock materials during seismic events [1,2,3]. Vs can be measured with either invasive (e.g., seismic cone penetration test, SCPT, cross-hole test, CHT, downhole test, DHT) or noninvasive (e.g., multichannel analysis of surface waves, MASW) techniques [4,5,6,7,8]. Integrating these tests with laboratory tests enables site characterization and the identification of potentially liquefiable lithological units, which is essential for hazard assessment [9,10]. While these methods offer the benefit of testing the soil in its natural environment, they tend to be time-consuming and economically infeasible for routine geotechnical exploration programs. Consequently, researchers have developed empirical correlations to estimate Vs using CPT parameters (cone tip resistance qc and sleeve friction fs). Many notable examples of such correlations can be found in the literature (e.g., see the works of Robertson and others [11,12,13]). However, these empirical models often exhibit limitations in generalization, as they are based on site-specific data and are unable to capture the complex, non-linear relationships inherent in the soil properties influencing Vs. Regardless of the techniques, there is unavoidable uncertainty in the final Vs profile, which contributes to uncertainty in site response predictions [14].

Accounting for Vs uncertainty is of paramount importance in engineering applications. While there has been a great deal of research on the significance of incorporating uncertainty of Vs into engineering practices, no universally accepted methods exist for realistically addressing this critical issue [14]. Matasovic and Hashash [15] highlighted the use of bounding-type profiles, which is developing upper and lower bounds (e.g., $\pm 20\%$) around a base Vs profile to represent potential variability. Additionally, statistical approaches, such as Monte Carlo simulations and the Toro [16] Vs randomization model, have been utilized to generate a range of Vs profiles that reflect the inherent uncertainty in subsurface conditions.

With the increasing availability of high-quality, freely accessible data (e.g., CPT data), the use of machine learning (ML) algorithms for soil characterization has become increasingly popular in recent years [17,18]. Research related to Vs prediction using various algorithms shows promising results. Notable examples include predictions of Vs from CPT parameters (e.g., see authors’ previous research and many others) [19,20,21,22]. ML algorithms excel at learning complex relationships within data without relying on the prior knowledge of the underlying mathematical models. This capability makes them highly effective for analyzing highly variable geotechnical data. However, the development of ML models must be approached with caution, as their performances heavily depend on several factors, including data preprocessing, feature selection, and hyperparameter tuning.

With the growing interest in predictive modeling for soil characterization, Bayesian generalized linear models (GLMs), have also emerged as a focus of study. Bayesian GLMs can be considered valuable tools alongside ML algorithms, particularly in geotechnical applications where data variability and uncertainty are prevalent. To the best of the authors’ knowledge, no previous study has reported on assessment of Vs predictions by integrating ML and Bayesian GLM approaches. The aim of this study is to develop Bayesian GLM and ML models to predict Vs using CPT parameters. For the ML model, the Extreme Gradient Boosting (XGBoost) algorithm coupled with Shapley Additive Explanations (SHAP) analysis was chosen due to its high performance in predictive models. Following the development of Bayesian GLM and XGBoost models, a comparative analysis was conducted against real SCPT data from different locations to validate the predictive capabilities of the proposed models.

2. Data-Driven Shear Wave Velocity (Vs) Prediction Models

In the following subsections, we will discuss Bayesian GLM and XGBoost algorithms for the development of Vs prediction models.

2.1. Bayesian Generalized Linear Model

The Bayesian framework has been utilized in geotechnical parameter estimations for some time now. This approach offers a powerful paradigm to estimate soil parameters in a probabilistic manner. However, studies on Vs prediction using the Bayesian framework are limited. Bayesian GLM provides a robust framework for integrating CPT data to enhance the accuracy and reliability of Vs predictions. Bayesian GLM are readily available through rstanarm [23] package and can be implemented through R programming language [23]. In this study, we formulate a Bayesian approach to model the relationship between Vs and CPT parameters using a GLM framework. The general form of Bayesian GLM can be expressed as follows [24]:

Let Y represent Vs and X represent CPT parameters (e.g., cone tip resistance, qc; sleeve friction, fs; friction ratio, Rf). Then, the relationship between Y and predictor variables X can be expressed as

$$Y=M\left(X\right)+E$$

(1)

where $Y$ is predicted Vs, $M\left(X\right)$ is the linear predictor function based on the predictor variables and $E$ is an error term.

The linear predictor can be defined as

$$M\left(X\right)={\beta }_0+{\beta }_{1}qc+{\beta }_{2}fs+\dots +{\beta }_p{X}_p$$

(2)

where ${\beta }_0$ is the intercept and ${\beta }_i(i=1,2,\dots ,p)$ is the coefficient of $X_i$.

The Bayesian GLM requires specifying prior distributions for the intercept and coefficients. The values of these parameters are equally likely to be positive or negative, but they are highly unlikely to be far from zero. These beliefs can be represented using normal distributions with a mean of zero and a small standard deviation. The default priors in rstanarm package perform well for most cases. Prior distributions for the coefficients can be defined as

$${\beta }_i~Normal\left(0,{\sigma }_{{\beta }_i}^2\right)$$

(3)

where ${\sigma }_{{\beta }_i}^2$ is variance of prior distribution of the coefficients. Similarly, the prior distribution for the intercept can be defined as

$${\beta }_0~Normal\left(0,{\sigma }_{{\beta }_0}^2\right)$$

(4)

where ${\sigma }_{{\beta }_0}^2$ is variance of prior distribution of the intercept.

The error term is assumed to follow normal distribution with a mean of 0 and variance ${\sigma }^2$).

$$E~Normal\left(0,{\sigma }^2\right)$$

(5)

Common choices for the prior distribution of the variance of error include the exponential, inverse gamma, and half-normal distributions. In this study, the exponential prior distribution (the default for GLM in the rstanarm package, with rate parameter λ = 1) was considered. It is worth noting that the priors are internally adjusted in rstanarm to account for data scaling and centering.

$$\sigma ~Exponential\left(\lambda \right)$$

(6)

Assuming normally distributed errors, the likelihood function can be expressed as

$$L\left(Y|X,\beta ,{\sigma }^2\right)={\prod }_{i=1}^n{\displaystyle \frac{1}{\sqrt{{2\pi \sigma }^2}}}exp\left(-{\displaystyle \frac{{\left(Y_i-M\left(X_i\right)\right)}^2}{{2\sigma }^2}}\right)$$

(7)

By using Bayes’ theorem, the posterior distribution of the parameters can be derived as

$$P\left(\beta |Y,X\right)\propto L\left(Y|X,\beta ,{\sigma }^2\right)P\left(\beta \right)$$

(8)

The actual values of posterior distributions of the parameters are rarely obtainable through analytical methods. Therefore, Markov Chain Monte Carlo (MCMC) simulation is commonly used to explore the posterior target distribution through random sampling [25,26]. In MCMC simulations, the Markov chains should converge to stationary distributions. To ensure that the chains mix well within parameter space, visual inspections of the chains and statistical tests, such as the Gelman and Rubin Rhat statistics [27], should be conducted. At convergence, Rhat is expected to be below 1.2, which indicates that the between-chain and within-chain variances for the model parameters are in good agreement [28,29,30]. This threshold serves as a guideline to confirm that the Markov chains have mixed well and that the posterior estimates are reliable and representative of the target distribution.

2.2. Extreme Gradient Boosting (XGBoost) Algorithm

XGBoost algorithm is a powerful machine learning algorithm known for its efficiency, scalability, and ability to handle large datasets with high-dimensional feature spaces and non-linear relationships [31]. In this study, XGBoost was selected to develop a Vs predictive model using seismic cone penetration tests (SCPTs) and SCPTs with a pore water pressure measurement (SCPTu) dataset. Each observation in the dataset represents a measurement from SCPT or SCPTu (e.g., cone resistance (qc), sleeve friction (fs), and depth, among others). These observations (features) are used as an input to the XGBoost model for the prediction of output (in this case the Vs). The XGBoost algorithm has the capability to discover patterns and relationships in the data without requiring pre-defined mathematical models. For instance, the relationships between Vs and geotechnical parameters (e.g., qc and fs) may involve non-linearities and are difficult to express mathematically. XGBoost overcomes this limitation by learning directly from the data, making it particularly suited for SCPT-based Vs predictions where the physical or mathematical relationships may be unknown or too complex to model explicitly. The model uses an ensemble of decision trees, built iteratively to minimize the prediction error by optimizing a loss function. Specifically, XGBoost combines weak learners (individual trees) in a sequential manner, where each tree corrects the errors of its predecessors.

In XGBoost, a regularization parameter was incorporated to reduce model complexity and prevent overfitting. The objective function of XGBoost is defined as follows [31]:

$$\mathcal{L}\left(\varphi \right)=\sum _{i=1}^{N}l\left(y_i,{\widehat{y}}_i\right)+\sum _{i=1}^T\Omega \left(f_t\right)$$

(9)

where $\mathcal{L}\left(\varphi \right)$ is objective function, $l\left(.\right)$ is loss function that measures difference between actual value $y_i$ and predicted value ${\widehat{y}}_i$, $N$ is number of observations, $T$ is number of estimators, and $\Omega \left(f_t\right)$ is regularization parameter defined as

$$\Omega \left(f_t\right)=\gamma T+{\displaystyle \frac{1}{2}}\lambda {\omega }_t$$

(10)

where $T$ is the number of nodes, and ${\omega }_t$ weight of the leaf nodes, and $\gamma$ and $\lambda$ are regularization coefficients that control the model’s complexity.

The predictive performance and generalization capability of XGBoost depends on the values of its hyperparameter. Therefore, the optimal hyperparameter values must be determined through hyperparameters tuning techniques. A next-generation optimization framework, Optuna, is widely used [32]. To prevent overfitting during the tuning process, k-fold cross validation techniques are commonly employed. This process helps the identification of the optimal parameter configuration that maximizes the model’s coefficient of determination ($R^2$) while avoiding overfitting.

To analyze the influence of input features on the optimized XGBoost model, partial dependence plots (PDPs) and SHAPs were utilized. SHAP, introduced by Lundberg [33], offers a unified framework for interpreting feature importance in the model. SHAP assigns a unique Shapley value (${\varphi }_{\mathrm{j}}$) to each feature based on its contribution to the model output. For a prediction $f\left(x\right)$, SHAP is expressed as [33]:

$$f\left(x\right)={\varphi }_0+\sum _{j=1}^N{\varphi }_{\mathrm{j}}$$

(11)

where ${\varphi }_0$ is base value (average model output for the dataset), ${\varphi }_{\mathrm{j}}$ is SHAP value for feature $j$and it is defined as

$${\varphi }_{\mathrm{j}}=\sum _{S\subseteq N\setminus \left\{j\right\}}{\displaystyle \frac{\left|S\right|!(\left|N\right|-\left|S\right|-1)!}{\left|N\right|!}}[f\left(S\cup \left\{j\right\}\right)-f\left(S\right)]$$

(12)

where $N$ is total number of input features, $S$ is a subset of features excluding $j$, $f\left(S\right)$is the model prediction for the subset $S$, $\left|S\right|$is the number of features in subset $S$, and $S\subseteq N\setminus \left\{j\right\}$ means that $S$ is a subset of the set of features $N$ excluding feature $j$. The feature with the highest absolute SHAP value is considered the most influential in the model’s decision-making process.

PDPs can be used to visualize the interactions between the input features and the target variable in machine learning models [34]. In this study, PDPs were employed to analyze the relationship between input features and the target variable, Vs. By marginalizing over all other features, PDPs enable a clear depiction of the isolated effect of a single feature on Vs. The PDPs can reveal whether the relationship between Vs and input features is linear, monotonic, or non-linear. From the PDP analysis, it is possible to identify key features that contribute significantly to the predictions of Vs.

3. Methodology

This study focuses on developing XGBoost and Bayesian GLMs for predicting Vs. Figure 1 presents the methodology used for the development and evaluation of these predictive models. The process begins with the preprocessing of SCPT and SCPTu datasets. The preprocessed dataset was used for the development of the XGBoost model, conducting partial dependence analysis to identify impactful features, and fitting the Bayesian GLM. Model evaluations include PDPs and SHAP analysis for the XGBoost model and a comparison of the predictive performance of both models. The following sub-sections provide a detailed description of the research data, performance metrics, and the assessment of the models’ predictive capabilities.

3.1. Training and Testing Dataset

The data used in this study were provided by Graz University of Technology [35] (Graz, Austria). The complete dataset consists of 1339 cone penetration tests (CPT, CPTu, SCPT, and SCPTu) conducted in Austria and Germany by Premstaller Geotechnik ZT GmbH [35] (Hallein, Austria). For this study, we utilized only the SCPT (46 tests in total) and SCPTu (50 tests in total) datasets. The SCPT and SCPTu can measure Vs at various intervals (e.g., at every 1 m) along the penetration depth, in addition to cone penetration data (i.e., cone tip resistance, qc and sleeve friction, fs). Before training the ML algorithm, the dataset underwent preprocessing. Raw data are rarely error-free or perfectly structured for ML training. Data often arrive in raw, messy, or incomplete forms, rendering them unsuitable for direct use in ML models. Initially, the dataset was filtered to include only SCPT and SCPTu tests. Unnecessary columns, and non-relevant measurements, were removed. Rows with zero or negative values for critical parameters were excluded. The differences in sampling intervals between CPT and Vs measurements were reconciled as per procedures described in [36]. This involved computing interval-based statistics (mean and coefficient of variation, cv) for relevant input features within the depth intervals where Vs data were recorded. Initial input features considered include mean cone tip resistance (qc_mean) and its coefficient of variation (qc_cv); mean sleeve friction (fs_mean) and its coefficient of variation (fs_cv); mean friction ratio (Rf_mean) and its coefficient of variation (Rf_cv); mean of soil behavior type index (Ic_mean); total overburden stress (σ,v); and depth of soil while Vs served as the target variable. Outliers in Vs were identified and removed using the interquartile range (IQR) method, with a threshold of 2.5 times the IQR above the third quartile. The cleaned and processed dataset was then saved for further modeling and analysis.

Table 1. Statistical summary of input features.

Metrics	qc_mean (MPa)	qc_cv (-)	fs_mean (kPa)	fs_cv (-)	Rf_mean (%)	Ic_mean (-)	Depth (m)	σ,v (kPa)	Vs (m/s)
Mean	4.24	0.24	50.53	0.26	2.29	2.74	13.38	254	237.6
STD	7.12	0.25	56.57	0.21	3.16	0.65	9.01	171	91.1
Minimum	0.02	0	0.35	0	0.09	0	0.5	9.5	22
Maximum	74.75	2.59	820.25	1.79	112.95	4.06	49.5	941	547
Count	3600	3600	3600	3600	3600	3600	3600	3600	3600

presents a statistical summary of the dataset after preprocessing, while Figure 2 displays a pair plot of the input features and the target variable, Vs. The diagonal of the matrix contains histograms that illustrate the distribution of each individual feature. The off-diagonal scatter plots depict the pairwise relationships between features. Additionally, the upper triangle includes linear correlation coefficients, highlighting the strength and direction of these relationships.

The XGBoost model was trained on 80% of the dataset, while the remaining 20% of the dataset was used to evaluate its performance. To ensure the model’s robustness and minimize the risk of overfitting, hyperparameter optimization was conducted using the Optuna framework. This optimization process was paired with 10-fold cross-validation, where the training dataset was divided into 10 equal subsets. In each trial, the model was trained on nine folds and validated on the remaining fold. The primary objective of the optimization was to maximize the coefficient of determination ($R^2$). The optimized model was used to conduct partial dependence analysis to examine the influence of each feature on Vs. Features with negligible importance were excluded, and the model was retrained using the updated feature set. For the Bayesian GLM, only features exhibiting a linear or nearly linear relationship with Vs were selected. The performance of the final XGBoost model was assessed using multiple statistical metrics. Additionally, the predictive performance of both the Bayesian GLM and XGBoost models was validated on a separate validation dataset.

3.2. Hungarian Seismic Cone Penetration Test (SCPT)

The cone penetration test (CPT) data have been in use by Hungarian geotechnical engineers for quite some time. A wealth of CPT data is readily available for research purposes. This study used SCPT data (see Figure 3) from Hungary (Pak site) to validate the predictive performance of Bayesian GLM and XGBoost models. The local soil deposit at the site is predominantly characterized by fluvial sediments from the Danube River. Details about the SCPT procedures and geological characteristics are available elsewhere and readers are encouraged to refer to [3,36,37,38,39].

3.3. Performance Measurements

The predictive performance of the Bayesian GLM and XGBoost models was evaluated by computing various performance and error metrics (

Table 2. Performance and error metrics used to evaluate Vs predictive capabilities of XGBoost and Bayesian GLM.

Metrics	Formula	Ideal Value	Equation. No
Correlation coefficient	$r={\displaystyle \frac{{\sum }_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=}^n{\left(x_i-\overline{x}\right)}^2}\sqrt{{\sum }_{i=}^n{\left(y_i-\overline{y}\right)}^2}}}$	1	(13)
Coefficient of determination	$R^2=1-{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{{\left(x_i-y_i\right)}^2}{{\left(x_i-\overline{x}\right)}^2}}$	1	(14)
Index of Agreement	$IA=1-{\displaystyle \frac{{\sum }_{i=1}^N{\left(y_i-{\widehat{y}}_i\right)}^2}{{\sum }_{i=}^N{\left(\left|y_i-\overline{y}\right|+|{\widehat{y}}_i-\overline{y}|\right)}^2}}$	1	(15)
Kling–Gupta efficiency	$KGE=1-\sqrt{{\left(r-1\right)}^2+{\left(\alpha -1\right)}^2+{\left(\beta -1\right)}^2}$	1	(16)
Mean squared error	$MSE={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}{\left(x_i-y_i\right)}^2$	0	(17)
Root mean squared relative error	$RMSRE=\sqrt{{\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}{\left({\displaystyle \frac{x_i-y_i}{x_i}}\right)}^2}$	0	(18)
Mean absolute error	$MAE={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}\left|x_i-y_i\right|$	0	(19)
Mean absolute relative error	$MARE={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}\left|{\displaystyle \frac{x_i-y_i}{x_i}}\right|$	0	(20)
Mean square relative error	$MSRE={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}{\left|{\displaystyle \frac{x_i-y_i}{x_i}}\right|}^2$	0	(21)
Mean bias error	$MBE={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}(x_i-y_i)$	0	(22)
Maximum absolute relative error	$MaxARE=Max\left(\left|{\displaystyle \frac{x_i-y_i}{x_i}}\right|\right)$	0	(23)

Let $X$ represent the measured Vs values and $Y$ represent the predicted Vs values. Then, $x_i$ represents the $i^{th}$ value of variable $X$, $\overline{x}$ is mean of $X$, $y_i$ is the $i^{th}$value of variable $Y$, $\overline{y}$ is the mean of $Y,$ $r$ represents the linear correlation coefficients, and $n$ is the number of data points.

). For Bayesian GLM, trace plot and R-hat, mean prediction, 95% credible intervals were used to assess the performance. Additionally, the linear correlation coefficients $\left(r\right)$ between the measured Vs values and the predicted values for the validation dataset were computed for both models.

The predictive performance of XGBoost model was further evaluated using several metrics, including the coefficient of determination ${(R}^2)$, index of agreement (IA), Kling–Gupta efficiency (KGE), mean squared error (MSE), root mean squared relative error (RMSRE), mean absolute error (MAE), mean absolute relative error (MARE), mean square relative error (MSRE), mean bias error (MBE), and maximum absolute relative error (MaxARE) [40,41,42]. These metrics provide a comprehensive quantitative evaluation of the models in predicting the Vs. IA, KGE, and $R^2$ indicate overall agreement and explanatory power, with ideal values close to unity, while error metrics quantify prediction errors (lower values indicate high performance).

4. Discussion

In the following subsections, the Vs prediction capabilities of the proposed models are evaluated. The Bayesian GLM was implemented in R while the XGBoost model was implemented in Python (version 3.11.7).

4.1. Evaluation of the Models

Table 3. Optimized hyperparameters values of XGBoost model.

Hyperparameters	Search Span	Optimized Values
Number of trees (n_estimators)	50–600	444
Learning rate	0.001–0.5	0.0094
Maximum depth of trees	1–10	10
Subsampling ratio	0.05–1	0.640
L1 Regularization term on weights (reg_alpha)	0.01–1	0.208
L2 Regularization term on weights (reg_lambda)	0.01–1	0.603
Column subsampling ratio per tree (colsample_bytree)	0.5–1	0.707
Minimum loss reduction required to make a further split (gamma)	0–10	8.79
Maximum number of bins for feature quantization (max_bin)	128–512	499
Balancing weight for positive and negative classes (scale_pos_weight)	0.1–10	6.09

presents the optimized critical hyperparameters of the XGBoost model. These hyperparameters were fine-tuned through Optuna with a 10-fold cross-validation technique to enhance the model’s predictive performance. The Optuna framework iteratively trained and evaluated models using the mean $R^2$ score across the cross-validation folds as the objective metric. By focusing on maximizing this metric, Optuna identified the optimal set of hyperparameters (see Table 3) that balanced accuracy and generalization. The results of the hyperparameter importance analysis are illustrated in Figure 4. As can be observed, the relative importance of each parameter is ranked based on its contribution to model performance. Among the hyperparameters, the learning rate (eta) emerged as the most influential factor. This underscores its pivotal role in determining the pace of model updates during training. Following the learning rate, the subsampling ratio and maximum tree depth (max_depth) were found to significantly affect model performance. The max_depth governs the complexity of the decision trees, allowing the model to capture complicated patterns in the dataset. Notably, the regularization term (lambda) exhibited relatively lower importance indicating that the other hyperparameters already provided sufficient control over model complexity.

The partial dependence analysis was conducted to explore the relationships between key input features and the predicted Vs. The results (see Figure 5) indicate that Vs exhibits positive correlations with σ,v, qc_mean, fs_mean, depth and qc_cv. The relationship between Ic_mean and Vs, however, is non-linear, indicating a complex influence on Vs predictions. This non-linear interaction can be effectively captured by the XGBoost model but could not be accommodated by the linear framework of the Bayesian GLM. Consequently, Ic_mean was excluded from the Bayesian GLM. In contrast, fs_cv, Rf_mean and Rf_cv show negligible influence on Vs, as observed from their flat PDP curves. These parameters were excluded from both the XGBoost and Bayesian GLMs to optimize performance. As a result, the final optimized XGBoost model was trained using σ,v, qc_mean, fs_mean, depth, qc_cv, and Ic_mean, while Bayesian GLM was fitted using σ,v, qc_mean, fs_mean, depth, and qc_cv.

The predictive performance of the optimized XGBoost model is presented in Figure 6. In Figure 6a, the Vs values predicted by XGBoost model for both the training and testing dataset are plotted against measured Vs values. The diagonal dashed green line (y = x) represents perfect prediction (i.e., the predicted Vs values match the measured Vs values). A good agreement is observed between predicted and measured Vs values. The residual (the difference between predicted and measured Vs) plots for the predictions are illustrated in Figure 6b. It can be observed that random scatters of residuals around zero (horizontal dashed green line) for both training and testing dataset, indicating unbiased predictions. On the right side of the scatter plots, histograms of the residuals are presented. These histograms offer a visual summary of the distribution of residuals. For the training dataset, the residuals appear to be symmetrically distributed around zero, with a relatively narrow spread, indicating strong predictive performance and minimal bias. In contrast, for the testing dataset, the residuals exhibit a slightly wider spread, reflecting the model’s reduced performance on unseen data. However, the residuals are still centered around zero, suggesting that the model generalizes reasonably well without significant overfitting. The results demonstrate the robustness of the XGBoost model and emphasize the importance of data preprocessing, feature selection, and hyperparameter optimization in improving model performance.

The predictive performance of the model across various performance metrics ($R^2$, IA, KGE) and error metrics are summarized in

Table 4. Performance of XGBoost model on test and train dataset.

Performance Metrics	Test Dataset	Train Dataset
$R^2$	0.54	0.91
IA	0.84	0.97
KGE	0.65	0.82
MSE	3792	781
RMSRE	0.39	0.22
MAE	41	19.6
MARE	0.21	0.11
MSRE	0.15	0.05
MBE	1.12	0.13
MaxARE	3.30	3.95

for both the training and testing dataset. For instance, the model achieved an $R^2$ of 0.54, an IA of 0.84 and KGE of 0.65 on the test dataset. For the training dataset, the model achieved an $R^2$ of 0.91, an IA of 0.97 and KGE of 0.82. The difference in performance between the training and testing datasets can be attributed to the complexity and variability inherent in the test data, which were not fully captured during model training. This difference is a common occurrence in ML models, particularly when the training data does not encompass all the potential variations present in the test set. Additionally, while the model was designed to generalize well and hyperparameter optimization techniques were employed to mitigate overfitting, the relatively lower performance on the test dataset reflects the challenges in accurately modeling highly complex or non-linear relationships between input features and the target variable Vs. However, these results are superior to those reported in [20]. The improvements can be attributable to the feature selection and data preprocessing step.

Figure 7 shows SHAP analysis results, where the input features are displayed on the y-axis in descending order based on their contribution to the Vs prediction. The most influential feature appears at the top. The x-axis represents SHAP values, which quantify the contribution of each feature to the model’s prediction. A positive SHAP value indicates that the feature increases the predicted Vs, whereas a negative SHAP value indicates a decrease in the predicted Vs. The SHAP values are represented by colored points that indicate the influence of the features on the model (Figure 7a). Each point reflects the impact of specific features on an observation from the entire database. The red color in the figure indicates high feature values while the blue color indicates low feature values. Positive SHAP values (red points) increase the Vs prediction while the negative SHAP value (blue points) decrease Vs prediction. For instance, for cone tip resistance (qc_mean), the red points predominantly on the positive side indicate higher qc_mean value generally increase Vs. On the other hand, the blue points on the negative side indicate lower qc_mean values generally decrease Vs. Interestingly, the instances of red points on the left or blue points on the right indicate a non-linear relationship between the feature and Vs. The global impact of features on Vs is shown by the mean absolute SHAP values of each feature (Figure 7b). Features with higher bar length indicate significant impact on the Vs. It can be observed that total overburden stress σ,v is the most influential feature, followed by cone tip resistance and sleeve frictions. On the other hand, the depth and soil behavior type index have lower impact on Vs. It is important to note that less impactful features were removed from the model, and only these features were incorporated into the final model.

The Bayesian GLM was fitted on the complete dataset using the final predictors identified through partial dependence analysis. It is important to note that the term “predictors” is more appropriate for the Bayesian GLM, while the term “input features” is suited for the ML models. In this study, both terms refer to the same geotechnical parameters. The model was developed using rstanarm package with four independent Markov chains, each consisting of 10,000 iterations, of which the first 1000 iterations were discarded as burn-in period. Weakly informative prior (normal prior (0, 2.5) for the intercept and coefficients, and an exponential prior (λ = 1) for the variance of the error term) were used. The performance of the model was evaluated using Rhat convergence diagnostic (evaluation of the agreement between between-chain and within-chain estimates for the model parameters), Monte Carlo Standard Error, MCSE (the ratio of standard deviation of the model parameters to the square root of the effective sample size), and trace plots of the parameters.

Figure 8 illustrates the convergence of Rhat values across the first 6000 iterations (burn-in period excluded) for the model fitted using four independent Markov chains. As can be observed, Rhat values remained consistently close to one (0.996–1.016) across iterations for each model parameter. This indicates the effective convergence of Markov chains. The trace plots (Figure 9) show the sampling behavior of Markov chains for each model parameter. It can be observed that the sampling demonstrates stable posterior distributions and reliable inference for all model parameters, as evidenced by well-mixed chains with no apparent trends.

Table 5. MCMC diagnostics for the Bayesian GLM.

Parameters	Mean	STD	MCSE (%)	Rhat	n_eff
Intercept	5.387	0.006	0.004	1.000	32,562
$\mathrm{q}\mathrm{c}\_\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}$	0.079	0.010	0.007	0.999	19,341
$\mathrm{q}\mathrm{c}\_\mathrm{c}\mathrm{v}$	0.015	0.007	0.004	1.000	28,913
$\mathrm{f}\mathrm{s}\_\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}$	0.092	0.010	0.007	1.000	19,305
$\mathsf{\sigma }\mathrm{v}$	0.092	0.786	0.65	0.999	14,770
$\mathrm{D}\mathrm{e}\mathrm{p}\mathrm{t}\mathrm{h}$	0.109	0.787	0.65	1.000	14,771

summarizes the MCMC diagnostics (mean, standard deviation, MCSE, Rhat, and effective sample size n_eff) for the model parameters. The final Rhat values (computed based on all iteration) are nearly one, further confirming the convergence of the MCMC chains. Additionally, high effective sample sizes and minimal MCSE values validate the robustness of the parameter estimates.

4.2. Validation of the Models

The Bayesian GLM was evaluated using validation SCPT data from Hungary to assess its predictive performance. Figure 10 presents the comparison between the measured and predicted Vs values across varying depths. Figure 10a displays the mean predictions (blue diamonds) with 95% credible intervals (2.5–97.5%), along with the measured Vs values (red stars). The close alignment of the predicted mean values with the measured values across the entire depths demonstrates the model’s accuracy and reliability. The lower plots (Figure 10b,c) provide a detailed analysis of predictions at a depth of 10 m (selected for illustration purpose only). The histogram (Figure 10c) shows the posterior distribution of Vs prediction at this depth. The blue curve represents the fitted probability density curve. It can be observed that the distribution follows lognormal distribution (attributable to data transformation during fitting the model). The vertical red arrow marks the measured value of Vs, while the blue arrows mark the 2.5% (left), mean prediction (middle), and 97.5% (right) credible interval. The 95% credible interval for this depth (142–627 m/s) encompasses the measured value (258 m/s). Additionally, Figure 10b summarizes the predictors used for this prediction. In general, the results demonstrate robustness of Bayesian GLM in capturing uncertainty and its capability to provide reliable predictions for Vs across varying depths. The use of credible intervals and posterior distributions further underscores its significance in quantifying prediction reliability.

In Figure 11, the predictive performance of the XGBoost model was compared against measured Vs values and Bayesian GLM predictions for the same validation dataset. The analysis involved comparing the outputs of the XGBoost model and the mean predictions of the Bayesian GLM across the Vs profile depth. For the Bayesian GLM, the 95% credible interval was also included to capture the uncertainty in the Vs predictions. It can be observed that both the XGBoost model’s predictions and the Bayesian GLM’s mean predictions align well with the measured Vs values across depths. Notably, the XGBoost model’s predictions fall within the 95% credible interval across the entire depth, showcasing its excellent predictive performance of the Vs. The uncertainty captured by the Bayesian GLM (the credible intervals) is crucial for site response analysis, where accounting for variability in Vs is essential. These intervals can be used to conduct sensitivity studies by varying the Vs profile within the predicted bounds, enabling the estimation of a range of possible seismic responses for the site. This approach minimizes the risk of underestimating hazards and provides a foundation for prioritizing additional site investigations in areas with significant uncertainty. Furthermore, to quantitatively assess the predictive performance of the two models, linear correlation coefficients $\left(r\right)$ were computed. Both models demonstrated good agreement with the measured Vs values, with the XGBoost model achieving a relatively higher $r$ value (0.39).

5. Conclusions

This study integrates Bayesian GLM and the XGBoost algorithm, coupled with SHAP and partial dependence analysis, to robustly predict Vs. Six key features were identified as the most influential factors for Vs prediction: cone tip resistance and its coefficient of variation, sleeve friction, soil behavior type index, depth, and overburden stress. These features formed the foundation for the final development of the XGBoost model. For the Bayesian GLM, cone tip resistance and its coefficient of variation, sleeve friction, depth, and overburden stress were identified as key predictor variables to fit the model on the same preprocessed datasets as the XGBoost model.

The XGBoost model demonstrated strong predictive performance, achieving high scores across various performance metrics and maintaining low prediction errors. SHAP analysis enhanced the model’s interpretability by quantifying the contribution of each feature to the Vs predictions. Overburden stress was identified as the most influential factor, followed by cone tip resistance and sleeve friction. The validation results on the SCPT dataset revealed a strong agreement between measured and predicted Vs values for both models, with the Bayesian GLM offering credible intervals representing prediction uncertainty.

The integration of Bayesian GLM and explainable ML algorithms signifies a robust framework for predicting Vs while leveraging the strengths of both approaches. Explainable ML model excel at identifying critical input geotechnical parameters and explaining their contributions to the prediction. This capability not only enhances the interpretation of the underlying factors influencing Vs but also delivers high predictive performance. On the other hand, Bayesian GLM offers complementary advantages by explicitly capturing the uncertainties associated with geotechnical variability, measurement errors, and other sources of uncertainty. This approach enhances confidence in Vs predictions and their applicability to seismic response assessments.

To further enhance the applicability of this work, developing a user-friendly computational platform for Vs prediction is recommended. Such a platform could integrate both Bayesian GLM and ML approaches, enabling geotechnical engineers to easily interpret predictions and uncertainties. Additionally, site-specific seismic response analyses using Vs profiles obtained from the two approaches should be conducted to further validate their performance in practical engineering scenarios.