Probabilistic indirect models for undrained shear strength: addressing significant data missing and variability with advanced imputation and machine learning techniques

Probabilistic indirect models for undrained shear strength: addressing significant data missing and variability with advanced imputation and machine learning techniques Haibin Xiong1; Shaoheng Dai2; Peng Lan3; Xuzhen He4; Chenxi Tong5; Sheng Zhang6; and Daichao Sheng7

Haibin Xiong, Ph.D. candidate School of Civil and Environmental Engineering, University of Technology Sydney, Ultimo, NSW 2007, Australia. Email: haibin.xiong@student.uts.edu.au Shaoheng Dai, Postdoctoral School of Civil and Environmental Engineering, University of Technology Sydney, Ultimo, NSW 2007, Australia. Email: shaoheng.dai@uts.edu.au Peng Lan, Postdoctoral School of Infrastructure Engineering, Nanchang University, Xuefu Road, 999, Nanchang 330031, China. Xuzhen He, Senior Lecturer (Corresponding author) School of Civil and Environmental Engineering, University of Technology Sydney, Ultimo, NSW 2007, Australia. Email: xuzhen.he@uts.edu.au Chenxi Tong, Associate Professor School of Civil Engineering, Central South University, Changsha, 410075, China. Email: cxtong@csu.edu.cn

Sheng Zhang, Professor School of Civil Engineering, Central South University, Changsha, 410075, China.

Daichao Sheng, Distinguished Professor School of Civil and Environmental Engineering, University of Technology Sydney, Ultimo, NSW 2007, Australia.

ABSTRACT: Accurate prediction of undrained shear strength (su) is crucial for geotechnical design, but is often hampered by substantial uncertainty in traditional empirical methods. This study uses the CLAY/10/7490 global database to develop probabilistic indirect models to predict su based on Atterberg limits and piezocone cone penetration (CPTU) measurements. Firstly, the dataset has a high missing data rate and variability. We test three imputation methods – multivariate normal (MN), multiple imputation by chained equations (MICE), and miss forest (MF) – to fill the missing values. To validate their effectiveness, a Probabilistic Extreme Gradient Boosting (PXGB) model is developed, and the imputation methods are evaluated by comparing the PXGB’s performance when trained on the imputed datasets against that on the original incomplete data. Secondly, the indirect model is built by integrating a multihead attention (MHA) mechanism into an artificial neural network (ANN) to enhance information extraction from limited data, which leads to the MHA-based probabilistic neural networks (MHA-PNN) model. The models’ performance, alongside a conventional MN-based prediction model, was evaluated using root mean square error (RMSE), coefficient of determination (R²), mean absolute percentage error (MAPE), conditional interval width (wCI), and coverage rate (CR). Results demonstrate that the proposed MN-enhanced MHA-PNN model substantially outperforms other models in both prediction accuracy and uncertainty quantification. These findings highlight the potential of this integrated strategy for building robust probabilistic indirect models in geotechnical applications, particularly when confronted with sparse and incomplete datasets.

Keywords: Uncertainty quantification; Multivariate normal; Multiple imputation by chained equations; Miss forest; Multi-head attention

1. Introduction Undrained shear strength (su), defined as the maximum shear resistance of soil under conditions where pore water pressure does not dissipate, plays a crucial role in determining the immediate shear capacity of soil. Its accurate assessment remains a significant challenge in geotechnical design, particularly in transient loading scenarios such as seismic events (Belkhatir et al., 2012; Mun et al., 2016), tunnelling construction operations (Soe & Ukritchon, 2023; Ukritchon & Keawsawasvong, 2020), and slope stability assessments (Griffiths & Yu, 2015; Jiang & Huang, 2018; Shogaki & Kumagai, 2008). Current engineering practice primarily employs two approaches for determining su: (1) laboratory-based measurements such as unconfined compression (UC) and unconsolidated undrained (UU) triaxial testing, and (2) empirical correlations combined with in-situ measurements including field vane shear (FV) and piezocone cone penetration (CPTU) (Yoshimine et al., 1999). Additionally, back-analysis is widely used to estimate unknown soil properties based on laboratory and/or in-situ measurements (Hu et al., 2025b; Zhang et al., 2024).But it requires repeated execution of the forward models, which can be time-consuming and computationally intensive.

All these methods are subject to various sources of geotechnical uncertainty. These uncertainties can be broadly categorized into four types (Hu et al., 2025a; Phoon & Kulhawy, 1999b): (1) inherent variability, resulting from the natural heterogeneity of geological formations; (2) statistical uncertainty, due to the limited data available for estimating soil parameters; (3) measurement error, caused by disturbances during sampling, testing procedures, and other influencing factors (Ding & Loehr, 2019; Lunne et al., 2006; Tanaka, 2000); and (4) transformation uncertainty, introduced when laboratory or in-situ measurements are converted into design parameters using

empirical model (Phoon & Kulhawy, 1999a). These compounded uncertainties in su are reflected in the highly variable coefficients of variation (COVs) and biases. As reported by Ding and Loehr (2019), COVs range from 0.06-0.56, while biases (λ=su/su_ref, where su_ref represents the reference value at the same depth) vary from 0.40 to 1.22 across different testing methods. Such variability fundamentally challenges conventional deterministic approaches in building empirical correlations, underscoring the need for probabilistic frameworks that can effectively quantify uncertainty in su predictions.

Recent advances in uncertainty quantification have leveraged Bayesian frameworks with Markov Chain Monte Carlo (MCMC) sampling to update prior distributions using laboratory and field data for geotechnical parameters (Jong & Ong, 2024; Li et al., 2025). However, repeated forward model evaluations in MCMC are computationally expensive. To address this issue, data-driven indirect models have been widely recognized as an effective approach to lighten this computational burden (Lan et al., 2025; Lan et al., 2024; Nguyen et al., 2025; Sun et al., 2023; Wang et al., 2023). Despite their advantages, indirect models are often constructed using relatively small datasets, primarily due to the limited availability of site-specific geotechnical data (Bozorgzadeh et al., 2019; Zhang et al., 2023). To improve the predictive performance, researchers have increasingly turned to multi-source geotechnical datasets to expand the available training data (Beesley & Vardanega, 2020; Yang & Liu, 2024; Zhang et al., 2022; Zhang et al., 2021). For instance, Zhang et al. (2021) utilized the F-CLAY/7/216 and SCLAY/7/168 datasets (D’Ignazio et al., 2016) to develop extreme gradient boosting and random forest models as indirect models for predicting su in soft sensitive clays. Naturally, this has led to growing interest in integrating geotechnical data from around the world for more generalized modelling. One notable effort is by Ching and Phoon

(2014b), who compiled the CLAY/10/7490 database with 10 geotechnical features and 7490 rows of data, which are from 251 studies conducted across 30 countries or regions. However, the database bears many missing values and its completion rate is only 34.1%, posing a significant challenge for data-driven modelling. Therefore, effectively utilizing such highly incomplete and highly variable geotechnical big data remains a key and unresolved challenge.

While discarding incomplete records is straightforward, it results in substantial data loss and forfeits valuable information. To address this issue, a range of imputation techniques to assign a value to the missing data have been developed. Basic statistical approaches include mean, median, and mode imputation, whereas more advanced methods encompass linear, polynomial, and spline interpolation (Ye et al., 2024). Interpolation methods, in particular, are widely applied in geotechnical engineering to handle missing data (Kim & Ji, 2022; Xu et al., 2024; Ye et al., 2024). In addition, more sophisticated machine learning (ML) approaches have been explored, such as extreme gradient boosting (XGBoost) (Chian et al., 2023; Xie et al., 2022), k-nearest neighbours (KNN) (Chian et al., 2023; Xie et al., 2022), Bayesian ridge regression (Piao et al., 2025), and miss forest (MF) imputation (Chian et al., 2023), among others. While these methods have shown improved performance in imputing missing data and capturing nonlinear relationships, most of them primarily focus on point estimation. In contrast, statistical imputation methods, which often rely on global data distributions, have the potential to reduce uncertainty more effectively by preserving overall statistical characteristics.

In this paper, we chose: (1) multivariate normal (MN) as a simple imputation baseline method; (2) multiple imputation by chained equations (MICE), which is widely adopted in practice and represents a standard statistical imputation approach; and (3) MF as a representative of advanced ML-based imputation techniques. These methods were compared to provide a comprehensive evaluation across distinct methodological paradigms. To validate their effectiveness, a Probabilistic Extreme Gradient Boosting (PXGB) model is developed, and the imputation methods are evaluated by comparing the PXGB’s performance when trained on the imputed datasets against that on the original incomplete data. Additionally, since the imputed data are inferred from limited available information, their quality is highly dependent on the amount and relevance of the observed data. To enhance the model’s capability to extract meaningful patterns from such limited inputs, a multi-head attention (MHA) mechanism is incorporated into a conventional artificial neural network (ANN) architecture. Building upon this foundation, we propose a probabilistic modelling framework that integrates imputation strategies with ML techniques. Finally, the proposed framework is applied to the CLAY/10/7490 database and compared with a conventional MN-based prediction model to validate its superiority in terms of both prediction accuracy and uncertainty quantification.

The structure of this paper is as follows: Section 2 introduces the CLAY/10/7490 database, focusing on the prediction of su and the imputation of missing data using the MN-based model. Section 3 analyses the missing data patterns within the CLAY/10/7490 database and compares the statistical characteristics of the MN, MICE, and MF methods. In Section 4, the PXGB and MHA-PNN models are developed, and

the impact of different imputation techniques on prediction accuracy and uncertainty quantification is assessed. Finally, Section 5 presents the conclusions of this study. 2. Prediction and Imputation Based on Multivariate Normal (MN) Model 2.1 Missing data distribution in CLAY/10/7490 database The CLAY/10/7490 database is a comprehensive clay parameters database and the largest publicly available global database for clay, incorporating site data from 251 studies across 30 countries or regions worldwide (Ching & Phoon, 2014b). However, its completeness rate is around 34.1%. To improve the database's utility, it is essential to first analyse the distribution of the missing data. Figure 1 provides a visualization of data completeness across 8 parameters: (1) Atterberg limit parameters: liquid limit (LL), plasticity index (PI), liquidity index (LI); (2) CPTU parameters: normalized cone tip resistance ((qt–σv)/σ'v), effective cone tip resistance ((qt–u2)/σ'v), normalized excess pore pressure ((u2-u0)/σv), and pore pressure ratio Bq; and (3) normalized mobile undrained shear strength (su(mob)/σv), which represents the mobilized in situ shear strength under undrained failure, converted from su obtained from various tests using transformation models (Ching & Phoon, 2014a, 2014b). The figure illustrates the distribution of missing data, with white areas representing missing values and black areas indicating available data. The extent of missing data varies across different parameters, with CPTU parameters exhibiting particularly high levels of incompleteness. Other parameters in the database, such as preconsolidation stress (σp) and sensitivity (st), are often difficult to obtain in practice. Therefore, to develop a probabilistic indirect model using easily accessible variables for predicting su(mob)/σv, we deliberately focused on these seven input features (LL, PI, LI, and the four CPTU parameters). Notably, while this global database offers valuable insight into general

trends and uncertainty across a wide range of clays, it is intended as a basis for firstorder estimates rather than a substitute for site-specific predictions.

Fig. 1. Missing data distribution across 8 parameters for the CLAY/10/7490 database

Table 1 provides a more detailed statistical summary, revealing that su(mob)/σv and Atterberg limit parameters have relatively lower missing rates (approximately 40% to 50%), whereas CPTU parameters suffer from severe data missing, with some variables missing over 90% of their values. Table 1. The statistics of the CLAY/10/7490 database Type Atterberg Limit

CPTU su

Variable

Mean

COV

Min

Max

LL 68.649 0.811 18.067 550.000 PI 40.210 1.095 1.895 493.000 LI 1.005 0.795 -0.775 6.455 (qt-σv)/σv 9.450 1.233 -1.041 110.561 (qt-u2)/σv 5.655 1.460 0.357 108.201 (u2-u0)/σv 4.732 1.342 -0.065 117.327 Bq 0.569 0.358 0.000 1.172 su(mob)/σv 0.541 1.670 0.019 27.755

Points 4057 4503 3795 880 752 668 1017 3779

Missing (%) 45.8 39.9 49.3 88.3 90.1 91.1 86.4 49.5

Beyond data completeness, the table reveals significant variability across parameters, with CPTU parameters showing particularly high variation, as (qt–σv)/σ', (qt–u2)/σ'v, and (u2-u0)/σv exceeding a COV of 1.2. The highest variability is observed in su(mob)/σv (COV = 1.670), reflecting substantial variability in these parameters. Such variations pose challenges for analysis and modelling. To enhance the utility of the database, implementing strategic data imputation and adopting robust approaches to handle variable uncertainty are crucial for ensuring reliable geotechnical applications. 2.2 MN-based prediction model for su(mob)/σv Given the high missingness and variability observed in the database, an effective probabilistic framework is essential for estimating missing values while preserving the underlying statistical relationships among parameters. The MN-based model offers a structured approach to achieving this; however, as all parameters exhibit non-Gaussian distribution, as shown in Figure 2, a transformation to approximate normality is necessary before constructing the MN-based model. The Johnson system of distributions, comprising three families (log-normal (SL), unbounded (SU), and bounded (SB)) provides a flexible method for this transformation (Slifker & Shapiro, 1980).

Fig. 2. Distributions of all parameters and their Johnson fitting

The mathematical formulations are as follows:

(

(

−1 ax  bx + ax  sin ( x − by ) / a y Johnson − SU：f ( x) =  2 ay 1 + 1 + ( x − by ) / a y

(

Johnson − SB：f ( x ) =

(

(

))

(1)

)

)

(

) ) (2)

  ax  bx + ax  ln ( x − by ) / a y − ln 1 − ( x − by ) / a y   ay ( x − by ) / a y 1 − ( x − b y ) / a y

Johnson − SL：f ( x ) =

(

(

)(

(

)

ax  bx + ax ln(a y ) + ax ln ( x − by ) / a y  ay ( x − by ) / a y

))

(3)

where ax , bx , ay , and by are the parameters of Johnson distributions, x is the real value of each variable, and () is the value of the probability density function (PDF) for the standard normal distribution.

Following the approach outlined by Slifker and Shapiro (1980), we classify all parameters under the SU family and estimate the eight transformation model parameters, as summarized in Table 2. Figure 2 illustrates the fitted results alongside the histograms of all parameters, including the corresponding means (μ) and standard deviations (σ). Despite the varying degrees of skewness and kurtosis, the fitted curves align well with the true distributions, verifying that the Johnson-SU distribution is appropriate for transforming the parameters into normal space. Table 2. Fitted parameters of the Johnson distributions for CLAY/10/7490 database Variable ax bx ay by Family LL 0.75 -1.18 7.29 33.99 SU PI 0.70 -1.10 4.88 15.05 SU LI 1.50 -1.19 0.66 0.28 SU SU (qt-σv)/σv 0.64 -1.10 0.71 3.68 0.77 -1.14 0.71 1.96 SU (qt-u2)/σv SU (u2-u0)/σv 0.78 -0.87 0.82 2.15 Bq 2.52 0.26 0.47 0.63 SU SU su(mob)/σv 0.80 -1.44 0.06 0.14 Next, to develop the MN-based prediction model, we adopt the framework for constructing the distribution as outlined in previous studies (Ching et al., 2014; He et

al., 2025). This approach involves transforming the original parameters into a normal space, followed by estimating the statistical correlations among the transformed parameters to predict the probability distributions of the target variables.

The first step involves selecting appropriate transformation functions to map the original parameter into a normal space. Based on the distribution patterns observed in Figure 2 and the statistical results in Table 1, we make the following determinations: Since su(mob)/σv is always positive and predominantly concentrated at lower values, a logarithmic transformation is the most suitable choice. Similarly, as Bq falls within the range (0, 1.2), an inverse sigmoid function is appropriate. For all other parameters, the inverse hyperbolic function is applied. Accordingly, the transformation functions are defined as follows: sinh −1 ( X i ),  Ti = sigmoid −1 ( X 7 / 1.2), ln( X ), 8 

for i = 1, 2, , 6 for i = 7

(4)

for i = 8

where T1 to T8 represent the transformed parameters, with X1 = LL, X2 = PI, X3 = LI, X4 = (qt-σv)/σv, X5 = (qt-u2)/σv, X6 = (u2-u0)/σv, X7 = Bq, and X8 = su(mob)/σv. To further evaluate the effectiveness of the transformation, Figure 3 presents the distributions of the transformed parameters (T1 to T8) and their correlations. The diagonal figures demonstrate that T1 to T8 exhibit normality. Furthermore, normality tests, including the quantile-quantile (Q-Q) plots and Anderson-Darling tests, confirm that the transformed parameters follow a marginal normal distribution, validating the suitability of the transformations. The estimated mean (µi) and standard deviations (i) of these transformed variables are listed in Column 2-3 of Table 3. The Pearson correlation coefficients jk between the Tj and Tk is estimated using the bootstrapping method suggested by Ching et al. (2014) and are listed in Table 3. The last column of Figure 3

shows that the Atterberg limit parameters exhibit little to no correlation with su(mob)/σv, whereas the CPTU parameters have a strong linear relationship with it. Table 3. Parameters of the MN-based prediction model for su(mob)/σv prediction Ti

µi

i

T1 T2 T3 T4 T5 T6 T7 T8

4.71 4.09 0.78 2.46 2.00 1.90 -0.10 -1.17

0.35 0.48 0.31 0.41 0.43 0.38 0.45 0.52

ij T1 1.00 0.96 -0.23 0.00 0.02 0.10 0.16 0.07

T2 0.96 1.00 -0.33 -0.07 -0.04 0.05 0.18 0.03

T3 -0.23 -0.33 1.00 -0.11 -0.17 -0.14 0.01 -0.13

T4 0.00 -0.07 -0.11 1.00 0.86 0.77 -0.39 0.80

T5 0.02 -0.04 -0.17 0.86 1.00 0.49 -0.61 0.69

T6 0.10 0.05 -0.14 0.77 0.49 1.00 0.17 0.65

Fig. 3. Distributions and correlations of transformed parameters (T1-T8)

T7 0.16 0.18 0.01 -0.39 -0.61 0.17 1.00 -0.27

T8 0.07 0.03 -0.13 0.80 0.69 0.65 -0.27 1.00

All the transformed parameters, T1 to T8, follow a multivariate normal distribution, then the conditional distribution p(Tpre i |Tj, Tk ,…) is also a normal distribution and can serve as a predictive model. Here, Tpre i

represents the predicted distribution for the

transformed variable i. When i = 8, it corresponds to the predicted distribution of the transformed su(mob)/σv. Tj, Tk ,… are known parameters. For example, for a soil sample, three values, X2 = PI, X5 = (qt-u2)/σv, and X6 = (u2-u0)/σv are available. Then, T2, T5, and T6 can be integrated into the conditional distribution to have a predictive distribution for Tpre 8 , i.e., prediction regarding su(mob)/σv. The conditional distribution can be explicitly expressed as

(

Ti pre ~ N ipre , ( ipre )



pre i

= i + cij

cik

( ) = c − c pre i

2

ii

ij

 c jj c jk  c jk ckk  

cik

2

)    

 c jj c jk  c jk ckk  

 c11 c12   1 0   1 =0    C = c21 c22 2     21      

(5) −1

T j −  j  T −   k   k  

   

−1

(6)

 cij  c   ik   

(7)

12   1

0

1

2

 0   

   

(8)

where C denotes the covariance matrix, which can be calculated from the standard deviations (i) and the Pearson correlation coefficients jk as in Eq.(8). The estimated mean (µi) and standard deviations (i) are those in Table 3. µpre and  pre are just i i calculated intermediate values. This predictive model is then used for su(mob)/σv on the CLAY/10/7490 dataset, and the results are depicted in Figure 4. Because for different soil samples, we have different numbers of available (Tj, Tk,…) as shown in Figure 1. The results are plotted on different figures depending on how many inputs are

available. The vertical lines indicate the 95% confidence intervals, and scatter plots illustrate the comparison between predicted mean values and corresponding actual values on a logarithmic scale.

Fig. 4. Prediction for su(mob)/σv based on the MN-based prediction model (Labels 17 represent the number of available inputs, ranging from 1 to 7)

Notably, in Figure 4(a), when only a single input feature is available, the model tends to generate nearly identical predictions, which is inaccurate. This suggests that with an insufficient number of input features, the model lacks the necessary information to distinguish between varying cases, resulting in oversimplified and inaccurate predictions. Figures 4(b)-(c) illustrate the prediction results when the number of input features increases to two and three. Their predicted mean values mostly fall between 0.1 and 0.8, while the real values cover a wider range, which still has low accuracy. The main reason for this issue is that when the number of available input features does not exceed three, the database mainly consists of Atterberg limit parameters, with very few CPTU parameters. Consequently, when the model predominantly relies on the

Atterberg limit parameters as the main input feature, it fails to capture the true probabilistic relationship of su(mob)/σv, leading to inaccurate predictions.

Figures 4(d)-(g) illustrate the prediction results of the MN-based prediction model using four to seven input features. Compared to Figures 4(a)-(c), the predictions demonstrate a stronger alignment with the diagonal reference line. Additionally, the confidence intervals are generally narrower, indicating reduced uncertainty and greater reliability in the model’s predictions. This improvement is due to incorporating CPTU parameters, highlighting the importance of selecting informative features. Without them, the model struggles to accurately capture the true underlying distribution of the target variable. Note that the number of data points in Figures 4(a)-(g) varies due to differences in sample sizes associated with each label (available input features).

Overall, the MN-based prediction model provides a statistically rigorous approach for predicting su(mob)/σv, reducing dependence on less reliable empirical methods, especially in the presence of missing data. The results underscore the importance of incorporating relevant input features to enhance predictive accuracy. As more informative variables are included, the model better captures probabilistic relationships, leading to narrower confidence intervals and improved agreement with observations. Additionally, its ability to explicitly account for uncertainty makes it well-suited for geotechnical applications where data sparsity is a challenge. 2.3 MN-based imputation method As described in Section 2.2, the MN-based prediction model utilizes the linear statistical relationships among transformed variables to construct their multivariate normal distribution, which is then applied to predict the probabilistic distributions of

the target variables. While the results in Figures 4(d)-(g) demonstrate strong predictive performance, this section aims to assess whether imputing missing input parameters leads to improved results compared to using the original incomplete database. It also serves as a foundation for developing new models that require a complete database for training. Notably, only input data will be imputed, while target data will remain unchanged to avoid potential issues in subsequent prediction and validation.

In the previous section, we used the conditional distribution as a predictive model (MNbased prediction model) for su(mob)/σv only, i.e., Tpre with i = 8. If we choose i as the i indices of missing parameters, we can also have predictions for the missing parameters, referred to as the MN-based imputation method, specifically:

(

pre pre pre Tmiss ~ N miss , ( miss )

2

)

(9)

pre −1 miss = miss + C miss,non-miss C non-miss,non-miss (Tnon-miss − non-miss )

(10)

pre −1  miss =  miss + C miss,non-missC non-miss,non-miss C non-miss,miss

(11)

where the subscript miss denotes the indices of missing values, and the subscript nonmiss represents the indices of non-missing (observed) values. Consequently, for a given sample, if a non-missing value is present in the input, all missing values for that sample can be inferred. This prediction is a distribution; we choose the mean as the imputed value. After imputation, the imputed dataset is shown in Figure 5 and compared with the original dataset. Notably, the imputation method effectively preserves the statistical relationships among all variables. This aligns with the fundamental characteristics of the multivariate normal model, providing an initial validation of the effectiveness of the imputation process. The figure presents the pairwise relationships among T1-T7, where the upper triangle displays scatter plots of the original database, the lower triangle

shows scatter plots of the database obtained using imputation method based on multivariate normal model (hereafter referred to as the MN database), and the diagonal elements display histograms of the MN database. Comparatively, the variables in the MN database retain the original distribution characteristics while slightly expanding the range of variable distributions, all while maintaining a high degree of consistency with the original database. Furthermore, the MN database accurately reproduces the relationships among variables, thereby reinforcing our confidence in utilizing this database for further analyses.

Fig. 5. Distribution of transformed variables in the original (upper triangular) and imputed (lower triangular) database

After comparing their distribution characteristics, the next focus is on whether the MN database can be effectively used to build a similar multivariate model and enhance the predictive performance. Based on the results in Figure 4, which show prediction differences with varying available inputs, we categorized the MN database by the original number of available input features and compared its performance with the original database. Figure 6 illustrates this, with the original database shown by blue squares and the new model by red circles. As the number of input features increases, the predictive performance of both models gradually improves. This observation aligns with the trend observed in Figure 4, where an increase in the number of input features led to improved model performance.

Fig. 6. Comparison of prediction results using the original and MN database (Labels 17 represent the number of available inputs, ranging from 1 to 7)

Compared to the original results, the predictions obtained using the MN database in Figures 6(a)-(c) exhibit a more concentrated distribution. However, the results presented in Figures 6(a)-(g) remain largely consistent overall. This suggests that

imputing the original data with the MN-based imputation method does not lead to a significant improvement in prediction accuracy. The primary reason for this is that both the MN-based imputation method and the MN-based prediction model are based on multivariate normal distribution theory. Therefore, the MN-based imputation method does not provide additional independent information to improve the model's predictive capability. Consequently, the model produces similar outcomes, which is expected due to the overlap in the theoretical foundations of both approaches. Hence, the subsequent chapters will explore alternative methods and models aimed at further improving prediction accuracy. 3. Imputation Based on Multiple Imputation and Miss Forest 3.1 Missing data patterns for CLAY/10/7490 Database As described in Section 2.3, we have successfully imputed missing data in the Clay/10/7490 database using the MN-based imputation method, while preserving the statistical relationships between different variables. However, the predictive performance has not improved. Therefore, in this chapter, we will explore alternative data imputation methods, including various statistical and machine learning techniques. Given the high missing data rates characteristic of the Clay/10/7490 database, it is crucial to first analyse the underlying missing data mechanisms before investigating alternative approaches. Common missing data mechanisms include missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). The distinction between these mechanisms lies in the relationship between the missing and observed (or unobserved) data: MCAR is entirely random, MAR depends on observed data, and MNAR is influenced by unobserved data.

Figure 7 illustrates the correlation matrix of missing data patterns among the transformed variables (T1-T8), visualized as a heatmap. The colour scale ranges from deep red (indicating strong positive correlation) to deep blue (indicating strong negative correlation), representing the strength of relationships between the missing data patterns of different variables. It is observed that the transformed Atterberg limit parameters (T1, T2, and T3) exhibit a strong positive correlation (> 0.6), suggesting that their missing data patterns are likely influenced by similar factors, such as sampling limitations or measurement inadequacies. Since these factors are easily observable in field, their missing data mechanism is best characterized as MAR. Similarly, the transformed CPTU parameters (T4, T5, T6, and T7) display stronger correlations (> 0.8) and, as presented in Table 3, exhibit a high missing rate (> 86%). This suggests that their missing values are primarily influenced by common factors, such as equipment limitations, further supporting the assumption that their missing data mechanism follows the MAR.

Fig. 7. Correlation matrix of missing data patterns among transformed variables

In contrast, T8 shows weaker correlations with other variables, and in some cases, even negative correlations (e.g., -0.2 with T1, T2, and T3). This suggests that the data collection method or source of T8 is more independent, and its missing data mechanism may be closer to MNAR or have more complex characteristics. However, as T8 serves as the target variable, imputing its missing values could introduce bias. Therefore, consistent with the previous method, we choose not to impute T8. As a result, its more complex missing data mechanism does not influence the overall validity of our findings.

Given that the missing data mechanism for these variables (T1-T7) aligns with the MAR assumption, we employ both statistical and machine learning-based imputation techniques to improve data completeness. Specifically, we will apply multiple imputation by chained equations (MICE) and miss forest (MF), as both methods effectively utilize inter-variable correlations to enhance imputation accuracy. By these methods, we aim to mitigate the impact of missing data while preserving the underlying statistical structure of the database. 3.2 Multiple Imputation by Chained Equations (MICE) MICE is a robust statistical technique for addressing multivariate missing data, widely used due to its flexibility and effectiveness in imputing missing values. This section introduces the MICE method and provides a brief comparison with the MN-based imputation method outlined in Section 2.3. The MN-based imputation method relies on the assumption that all data follow a multivariate normal distribution, using this theoretical framework to estimate and impute missing data. In contrast, MICE offers greater flexibility by modelling each variable's missing values based on its own distribution, without assuming a multivariate normal distribution. Moreover, MICE do

not require variables to be transformed into an approximately normal distribution before imputation, making it more adaptable for handling diverse data types.

The MICE process is iterative, using regression models to impute missing data sequentially. In each iteration, other observed variables serve as predictors to estimate the missing values of the current variable. This approach models missing data based on the specific distribution of each variable, providing a more accurate imputation method. The main steps of the MICE procedure for imputing missing data are as follows: (1) Use a simple imputation method (e.g., mean imputation) to fill in the missing values as initial guesses. (2) Select a variable with missing values as the target variable and use other variables as predictors. Apply a regression model (such as linear regression, logistic regression, or other more complex models) to predict and impute the missing values of the target variable. (3) Repeat step (2) for all variables with missing values, using newly imputed values to refine subsequent imputations. (4) Iterate steps (2) and (3) multiple times until the imputed values converge, ensuring stability in the results. (5) Repeat this process several times (typically 5 to 10 iterations), generating multiple imputed databases. These databases are then analysed independently, and the results are combined to reduce uncertainty in the imputation process.

Initially, based on the results from the MN-based imputation method, each sample in the database was labelled according to the number of available input features, with labels ranging from 0 to 7. The original database, containing 4,973 samples, excluded data points with a label of 0. For comparative analysis, three subsets were created based on label thresholds: The first subset (Subset 1), with 4,316 samples, had labels greater

than 1; the second subset (Subset 2), with 3,747 samples, contained labels greater than 2; and the third subset (Subset 3), with 809 samples, consisted of labels greater than 3.

The MICE method was applied for 10 iterations to each of these subsets. The results are presented in Figure 8 as histograms. The top row of Figure 8 corresponds to the original dataset (excluding data points labelled as 0), while the second to fourth rows correspond to Subsets 1 to 3, respectively.

Fig. 8. Comparison of the MICE datasets (red dotted lines with markers) and the original datasets (black histograms). From top to bottom are the original dataset and Subsets 1 to 3.

As shown in the histograms of Figure 8, the overall shapes of the distributions remain consistent across the original dataset and Subsets 1 to 3. Minor differences can be observed in the distributions of variables T4, T5, and T6 for the original dataset and Subsets 1 and 2 (Figures 8(a4)-(a6), 8(b4)-(b6), and 8(c4)-(c6)), where slight shifts in the density and a slight increase in peak height are noticeable. For other variables, the MICE and original data distributions are nearly indistinguishable.

Notably, as the data subsets progress from the full dataset to Subset 1 and then to Subset 3, the difference between the MICE and original data distributions decreases. This trend suggests that the MICE method becomes increasingly robust in imputing missing values as the number of available features increases. Across all subsets, the MICE data consistently resemble the original data, with only minimal distributional differences observed. 3.3 Miss Forest (MF) MF is a variant of MICE and, like MICE, follows an iterative modelling approach, using known features to impute missing values. The key difference between MF and MICE lies in the imputation model: MICE uses traditional regression models, while MF replaces these with random forest models. As a non-parametric machine learning method, random forest is capable of capturing complex non-linear relationships between variables, making MF particularly effective for data with complicated interactions. Unlike MICE, which generates multiple imputed datasets, MF typically performs a single imputation. Both methods are widely used for missing data imputation. In this section, we will discuss the performance of the MF method when applied to impute the Clay/10/7490 database.

Like the analysis process in the previous section, we applied the MF method to impute missing values for each dataset, including the original database and Subsets 1-3. These datasets were used to investigate how different known data features impact the inference of missing data. Initially, missing values were imputed using the means of the variables to ensure the dataset was complete for model training. The algorithm then iteratively updates the missing values by training a random forest model on the available data. The process stops when the maximum change in imputed values between

iterations falls below a predefined convergence threshold, indicating convergence. Key hyperparameters include the maximum number of iterations (set as 10 in this study), the convergence threshold (10-4), the number of trees in the random forest (80), and the maximum depth of the trees (100).

Figure 9 presents a series of histograms comparing different variables (T1 to T7) across each subset. The black bars correspond to the original data, while the red bars represent the data imputed using the MF method. In Figures 9(a4), 9(b4), and 9(c4), the imputed data appear more concentrated than the original data distribution. For variables T4, T5, and T6, slight differences in the histogram density are observed between the original database and Subsets 1 and 2. Notably, in Subset 3, the distribution of the imputed data closely matches that of the original data, consistent with the results obtained using the MICE method. Overall, in most cases, the distribution of imputed values aligns closely with the original data, indicating that the imputation process is generally successful.

Fig. 9. Comparison of the MF datasets (red dotted lines with markers) and the original datasets (black histograms). From top to bottom are the original dataset and Subsets 1 to 3.

For better visualization and comparison, Figure 10 displays a series of violin plots showcasing the performance of the MN, MICE, and MF imputation methods on the Clay/10/7490 database. Each violin plot represents the probability distribution of the dataset after imputation, with the central black bar representing the quartiles of a corresponding box plot, while the white bar denotes the median.

In Figure 10, the distribution plots of variables T4-T7 clearly show that the database imputed by the MN database is more concentrated compared to the original database. Additionally, the database imputed by the MF method (hereafter referred to as the MF database) exhibits some extreme values and outliers, whereas the MICE method produces a more reliable database. This is primarily because, in the MICE method, the database is imputed 10 times, generating 10 different versions that are subsequently combined to form the final database.

Fig. 10. Comparison of the distributions of original and imputed values (MN, MICE, MF) for each variable in the CLAY/10/7490 database

To quantitatively assess the similarity between the original and imputed data distributions, we performed the Kolmogorov-Smirnov (K-S) test, which evaluates the

maximum difference between two cumulative distribution functions (CDFs). In the KS test, smaller D values and larger p-values (typically p > 0.05) indicate a closer match between distributions. As shown in Table 4, the MN method yields significantly higher D values for T4-T7 (0.222-0.290, p < 0.001), indicating overly concentrated distributions. The MF method shows moderate deviation (D = 0.166-0.211, p < 0.05). In contrast, MICE achieves the best performance with low D values (e.g., 0.044 for T7) and high p-values (p > 0.05), effectively preserving the original distribution characteristics. The poorer results for T4-T6 across all methods (p < 0.0001) are mainly due to their high missing rates (>85%), making accurate recovery difficult. Table 4. K-S test results comparing original and imputed distributions Ti T1 T2 T3 T4 T5 T6 T7

MN D=0.040, p=0.002 D=0.033, p=0.011 D=0.064, p=0.000 D=0.222, p=0.000 D=0.256, p=0.000 D=0.289, p=0.000 D=0.274, p=0.000

MICE D= 0.009, p= 0.986 D= 0.005, p= 1.000 D= 0.005, p= 1.000 D= 0.136, p= 0.000 D= 0.097, p= 0.000 D= 0.117, p= 0.000 D= 0.044, p= 0.070

MF D=0.029, p=0.048 D=0.023, p=0.157 D=0.030, p=0.036 D=0.200, p=0.000 D=0.185, p=0.000 D=0.211, p=0.000 D=0.166, p=0.000

Overall, as shown in Figure 10 and Table 4, the MICE database most closely aligns with the original database, followed by the MF database and then the MN database. However, these findings must be interpreted with caution, as statistical similarity represents only a necessary (not sufficient) criterion for evaluating imputation quality. Since the true values of missing data are unknown, it is impossible to determine which method is truly better. While preserving distributions and correlations offers useful insights, it does not guarantee that individual imputed values are accurate. This issue is especially critical for variables with high missingness (e.g., CPTU parameters with over 90% missing), where the imputed values reflect the predictions of the imputation method rather than the true data characteristics. Given these limitations, a more

meaningful evaluation is the effect of each imputation method on target‑variable prediction, which is examined in the following section. 4. Two Probabilistic Indirect Models 4.1 Probabilistic extreme gradient boosting (PXGB) Extreme gradient boosting (XGBoost) is a widely used and efficient ensemble learning algorithm built on a flexible gradient boosting framework, well-suited for structured data prediction tasks. A key advantage of XGBoost is its native ability to handle missing data during training. This eliminates the need for explicit imputation and simplifies the data preprocessing pipeline. Therefore, we can directly compare the performance of XGBoost applied to the original datasets versus the imputed datasets to evaluate the effectiveness of the imputation strategies (MN, MICE, and MF) discussed in Sections 2 and 3.

In this study, we extend XGBoost into a probabilistic framework, referred to as probabilistic XGBoost (PXGB), by assuming that the target variable follows a normal distribution. To enable uncertainty quantification in predictions, we adopt the negative log-likelihood (NLL) as the loss function for model optimization. This process can be formally expressed as:

(

T8pre ~ N  pre , ( pre )

(  , ) =f pre

pre

PXGB

2

)

(T1 , T2 , ..., T7 )

(12) (13)

N

L = − log ( p( y∣ i xi ; ) )

(14)

i =1

where Tpre represents the predicted probabilistic distribution for the transformed target 8 variable su(mob)/σv, L is the loss function (NLL), N is the number of data points,

p(yi∣xi;θ) is the probability of the target variable yi given the input xi and model parameters θ. 4.2 Multi-head attention-based probabilistic neural network (MHA-PNN) Attention mechanism is inspired by how humans focus on important details in a scene while ignoring irrelevant information. Multi-head attention (MHA) extends this concept by projecting queries (Q), keys (K), and values (V) into multiple lowdimensional subspaces. This allows for parallelized feature learning and enhances the extraction of information from multiple perspectives. Therefore, in this study, the MHA layer is incorporated into a traditional neural network model, with the probabilistic distributions of the target as the output. This model is referred to as the multi-head attention-based probabilistic neural network (MHA-PNN). The architecture of the MHA-PNN model comprises four key modules, including an input layer, a hidden layer, an MHA layer, and an output layer. The process begins with the input layer, which receives a series of transformed input features [T1, T2, ..., T7]. These features are subsequently mapped into the hidden layer, where they are represented as corresponding vectors [d1, d2, ..., dd]. Next, the hidden representations are processed through the MHA layer, yielding output vectors [o1, o2, ..., oj]. These outputs are then normalized to obtain the final representations [d1, d2, ..., dj], which serve as the predictors for estimating the probability distribution of the target variable. Additionally， to mitigate overfitting, dropout regularization is incorporated into both the hidden layer and the MHA layer. The overall architecture of the MHA-PNN model is illustrated in Figure 11.

Fig. 11. The architecture of the MHA-PNN model

As shown in Figure 11, the key aspect of the model lies in the implementation of the MHA layer and the subsequent probabilistic distribution prediction. For the MHA layer, linear transformations are applied to generate the Q, K, and V matrices. Specifically, each attention head independently projects the hidden representations using different trainable weight matrices, ensuring diverse feature extraction and capturing distinct aspects of the input information. The transformation for each attention head is given by the following equations: Qh = DWQ ,h + bQ ,h =  qh ,1 , qh ,2 , , qh,d  K h = DWK ,h + bK ,h =  kh ,1 , kh ,2 , , kh ,d 

(15)

Vh = DWV ,h + bV ,h = vh ,1 , vh ,2 , , vh ,d 

where WQ,h, WK,h, and WV,h represent the weight matrices for the Qh, Kh, and Vh, respectively, while bQ,h, bK,h, and bV,h are the corresponding bias terms. Here, h refers to the number of attention heads, d is the number of neurons in the hidden layer, and D represents the hidden representations. Each attention head functions independently, capturing distinct aspects of the input information. It leverages scaled dot-product

attention to compute the relevance between features, with the specific computation given as follows:

Q K T  H h = Attention(Qh , K h ,Vh ) = softmax  h h  Vh  d k  

(16)

where dk is the dimension of Kh, which is used to scale the dot product between Qh and Kh, the softmax function is the normalized exponential function, and Hh represents the output of the attention mechanism for the h-th attention head. The outputs from these heads are then concatenated (Concat) and linearly transformed to yield the final multihead attention output. To preserve the original features, we incorporate a residual connection in the multi-head attention mechanism. The final formula is as follows:

O = Multihead( H ,Wo , Bo , T ) = concat ( H1 , H 2 , , H h )Wo + Bo + T

(17)

where Wo is the weight matrix used to combine the information from different attention heads, and Bo is the corresponding bias term. The outputs of the MHA layer are normalized and then passed to the output layer, where a probabilistic modelling approach is used to obtain the probability distribution of the target variable. In the output layer, as shown in Eqs (12)-(14) and consistent with the PXGB model, the output is assumed to follow a normal distribution, and the NLL is adopted as the loss function. Notably, the softplus function is used to ensure that the predicted standard deviation remains positive, with a small constant (10-6) added for numerical stability. 4.3 Hyperparameter configuration Hyperparameter optimization is essential for improving the performance of machine learning models and ensuring accurate predictions of the target variable. To minimize the influence of hyperparameter settings on the performance of the PXGB and MHAPNN models, we employ Optuna, a hyperparameter optimization technique, to automate the identification of the optimal hyperparameters. Specifically, for the MHA-

PNN model, we focus on optimizing the following hyperparameters: the number of attention heads (Natt), the key dimension (dk), the number of units in the dense layer (Nunits), the dropout rate (rd) in both the dense and MHA layers, the learning rate (η), and the activation function (act). The initial search ranges for these hyperparameters are as follows: Natt and dk from 2 to 20, Nunits from 20 to 200, rd from 0.2 to 0.5, η from 10-4 to 10-1, and act options being ReLU, Sigmoid, and Tanh. For the PXGB model, the optimized hyperparameters include the learning rate (η), maximum tree depth (Dmax), minimum child weight (wmin), regularization term (γ), feature subsampling ratio (rfea), and instance subsample ratio (rins). Their corresponding search spaces are: η from 10-4 to 10-1, Dmax from 1 to 20, wmin from 1 to 20, γ from 0 to 0.2, rfea from 0.5 to 1, and rins from 0.8 to 1. For each dataset, a total of 10 trials were performed to identify the bestperforming configurations. The optimal hyperparameter configurations are summarized in Table 5. Table 5. Optimal hyperparameter configurations for the PXGB and MHA-PNN models Ti

Ori. Natt dk Nunits rd act η (10-2) 4.29 Dmax 17 wmin 10 γ 0.16 rfea 0.97 rins 0.86

PXGB MN MICE 2.86 6.87 19 11 10 8 0.14 0.16 0.78 0.79 0.86 0.88

MHA-PNN MF MN MICE MF 20 7 9 19 2 13 184 84 152 0.31 0.23 0.41 ReLU Tanh ReLU 8.15 0.012 0.21 0.10 20 14 0.18 0.87 0.86 -

The results show that the optimal hyperparameters vary across datasets, emphasizing the necessity of dataset-specific tuning. Leveraging these optimal models enhances the reliability and credibility of subsequent comparative analyses, ensuring that the findings are not affected by hyperparameter configurations.

4.4 Performance comparison and uncertainty quantification Based on the optimal hyperparameters determined in the previous section, we constructed the corresponding optimal model for each dataset. To ensure reliable evaluation, we adopted stratified 5-fold cross-validation, which maintains consistent label distributions across all folds. This approach prevents the absence of certain label types in any given fold, ensuring that model performance is assessed fairly and comprehensively. The random seed was fixed at 42 during data splitting and training to ensure reproducibility. In each fold, the model was trained on four folds and validated on the remaining one. Figures 12-13 present the validation results for each fold, comparing the PXGB and MHA-PNN models across different datasets.

Figure 12 presents the prediction results of the PXGB model under different labels, ranging from 1 to 7 for different rows. From left to right, each column displays the validation results for folds 1 through 5. In this figure, green, blue, red, and purple represent the predictions of models trained on the original dataset, the MN dataset, the MICE dataset, and the MF dataset, respectively. It is evident that the predictions across the five folds are visually consistent, demonstrating the model’s robustness and stability across different data splits. Moreover, the predictions based on the imputed datasets (MN, MICE, and MF) generally exhibit improved accuracy, as the results using the original dataset (in green) deviate more noticeably from the diagonal line. This visual deviation indicates larger prediction errors, thereby validating the effectiveness of the data imputation methods in enhancing model performance.

Fig. 12. PXGB model predictions across five folds using original and imputed datasets (left to right: folds 1-5; green = original dataset, blue = MN dataset, red = MICE dataset, purple = MF dataset).

Fig. 13. MHA-PNN model predictions across five folds using imputed datasets (left to right: folds 1-5; blue = MN dataset, red = MICE dataset, purple = MF dataset).

Figure 13 presents the prediction performance of the MHA-PNN model under five-fold cross-validation, using only the imputed datasets (blue for MN, red for MICE, and purple for MF). Similar to the PXGB model, the prediction results appear visually consistent across all five folds, demonstrating the stability of both proposed models. Notably, compared to PXGB, the MHA-PNN model shows fewer points deviating from the diagonal line, indicating improved prediction accuracy. However, as illustrated in Figures 12 and 13, as well as in the earlier MN-based prediction results shown in Figure 4, when the number of available input features is less than or equal to three, the model exhibits consistently poor predictive performance, regardless of whether the original or imputed dataset is used. To ensure a fair comparison between models, these lowinformation input cases (i.e., with ≤3 input features) are excluded from the subsequent evaluations, and only results with more than three input features are considered.

In this study, three commonly used evaluation metrics are employed to assess prediction accuracy: root mean square error (RMSE), coefficient of determination (R²), and mean absolute percentage error (MAPE). Additionally, to further evaluate the quality of uncertainty quantification, two additional metrics are employed: the width of the conditional interval (wCI) and the coverage rate (CR), which are defined as follows:

(

1 N pre  CR =  1 yireal   yipre ,L , yi ,U  N i =1 wCI =

1 N pre pre   yi,U − yi,L  N i =1 

)

(18)

(19)

where N is the total number of dataset samples, yreal denotes the real value of the i-th i pre pre sample, ypre i,M , yi,U and yi,L represent the predicted median, upper, and lower bounds of the

95% confidence interval, respectively, and 1(⋅) is the indicator function, which equals 1 if the real value falls within the predicted 95% confidence interval, and 0 otherwise.

Among these metrics, CR shows how often the true values fall inside the predicted range. A higher CR means better coverage. wCI indicates the average width of the predicted range. A smaller wCI means the model is more confident, as long as CR remains high enough. Table 6. Evaluation metrics of PXGB, MHA-PNN, and MN-based prediction (MNpred.) models on the original and imputed datasets Model (Dataset) MAPE RMSE R² wCI CR MN-pred. (Ori.) 0.48 0.59 -0.06 0.72 0.80 PXGB (Ori.) 1.72 1.17 -5.62 434.62 0.97 PXGB (MN) 1.22 0.78 -1.31 24.86 0.97 PXGB (MICE) 1.00 0.76 -1.29 28.81 0.93 PXGB (MF) 1.19 0.78 -1.40 44.95 0.95 MHA-PNN (MN) 0.34 0.37 0.55 1.03 0.94 MHA-PNN (MICE) 0.46 0.45 0.32 1.58 0.98 MHA-PNN (MF) 0.46 0.39 0.50 1.47 0.97 Table 6 provides a comprehensive comparison of the PXGB, MHA-PNN, and MNbased prediction models across both original and imputed datasets. Among the evaluated metrics, the PXGB model using the original dataset shows the poorest performance, with the least favourable metrics (MAPE = 1.72, RMSE = 1.17, R2 = -5.62, wCI = 434.62), reflecting its limited predictive capability when handling highly incomplete data. In contrast, notable improvements are observed for PXGB when applied to imputed datasets, particularly the MICE dataset, where the performance improves significantly (MAPE = 1.00, RMSE = 0.76, R2 = -1.29, wCI = 28.81), demonstrating the positive impact of effective data imputation. Specifically, compared to the PXGB model on the original dataset, the PXGB model on the MICE dataset reduced MAPE by approximately 42%, RMSE by 35%, increased R² from -5.62 to 1.29, and decreased wCI by about 93%, while CR remained stable.

Additionally, the best overall results are achieved by the MHA-PNN model trained on the MN dataset, with the lowest MAPE (0.34) and RMSE (0.37), as well as the highest R2 (0.55) and strong uncertainty metrics (wCI = 1.03, CR = 0.94). On the MN dataset, the MHA-PNN model substantially outperformed the PXGB model, reducing MAPE by approximately 72%, RMSE by 53%, increasing R2 from -1.31 to 0.55, and decreasing wCI by about 96%, while CR remained stable. The results collectively underscore that predictive performance is highly sensitive to both the choice of prediction model and the imputation method. Advanced architectures such as MHA-PNN, when paired with appropriate imputation strategies, can substantially enhance the accuracy and reliability of predictions, particularly in data-scarce geotechnical settings.

To maintain consistency, the three metrics (MAPE, RMSE, and wCI) were transformed by taking their reciprocals, ensuring that higher values uniformly indicate better performance. Subsequently, all five metrics were normalized using min-max scaling. This allows for a more intuitive comparison, as shown in Figure 14, which presents a radar chart of the normalized evaluation metrics across various models and datasets. In this figure, the blue colour scheme, ranging from dark to light, corresponds to the PXGB model results on the original (Ori.), MN, MICE, and MF datasets, respectively; the purple colour scheme, also from dark to light, represents the MHA-PNN model applied to the MN, MICE, and MF datasets; the green dotted line with markers denotes the MNbased prediction (MN-pred.) model on the original dataset. It is evident that the dark blue line with markers (PXGB on the original dataset) exhibits the poorest overall performance, indicated by the smallest radar area. Conversely, the purple dotted lines with markers representing MHA-PNN consistently achieve superior performance, occupying the outermost regions of the chart.

Fig. 14. Radar chart of normalized evaluation metrics for different models and datasets

Overall, these findings confirm that MHA-PNN significantly outperforms both PXGB and MN-pred. models in predictive accuracy and uncertainty quantification. Among the imputation methods applied to the MHA-PNN model, the MN approach offers the best trade-off between prediction accuracy and reliable uncertainty quantification. 5. Conclusions The prediction of su is essential for applications such as slope stability analysis and foundation bearing capacity evaluation. However, traditional estimation methods, including laboratory tests and in-situ techniques, are often subject to significant uncertainty. Moreover, the sparsity and variability of geotechnical data greatly limit the development of robust predictive models. To address these issues, this study proposes a probabilistic modelling framework that combines statistical and machine learning techniques to generate more reliable predictions of su along with uncertainty quantification. The framework is developed using the CLAY/10/7490 database, which presents substantial challenges due to its low completion rate (around 34.1%) and high

variability. To mitigate data incompleteness, we systematically explore the effects of different imputation strategies (MN, MICE, and MF) and develop an MN-based prediction model, a probabilistic XGBoost (PXGB) model, and a neural network variant for predicting probabilistic distributions of su.

Initially, the MN-based prediction model was applied to the CLAY/10/7490 database to predict the probabilistic distribution of su(mob)/σv using seven input features, including Atterberg limit parameters and CPTU parameters. The results revealed that the number of input features had a significant impact on prediction accuracy. Particularly, when fewer than four features were used, the model predictions became largely unreliable. Subsequently, the MN, MICE, and MF imputation methods were applied to address the data incompleteness. Statistical analyses revealed that although the imputed results remained sensitive to the number of input features, the distributions of the completed databases closely matched those of the original database. However, we emphasize that imputed values should not be interpreted as accurate reconstructions of the true, unobserved data; therefore, site-specific investigations remain essential for reliable decision-making.

The effectiveness of the MN, MICE, and MF imputation methods was validated by comparing the performance of the PXGB model trained on imputed datasets against its performance on the original incomplete dataset. To enhance the model’s capability in extracting information from limited inputs, an MHA mechanism was incorporated into the traditional ANN architecture. By using the predicted probability distribution of the target variable as the output, a novel MHA-PNN model was developed. Comprehensive evaluations using RMSE, R², MAPE, CR, and wCI demonstrate that the proposed MHA-

PNN consistently achieves lower prediction errors, higher coverage rates, and narrower uncertainty intervals. For example, on the MN dataset, the MHA-PNN model outperformed PXGB, reducing MAPE by approximately 72%, RMSE by 53%, increasing R2 from -1.31 to 0.55, and decreasing wCI by about 96%, while CR remained stable.

In summary, the MHA-PNN model exhibits clear advantages over the MN-based prediction and PXGB models in both predictive accuracy and uncertainty quantification across various datasets. This model effectively captures complex feature interactions, improves prediction accuracy under sparse input conditions, and provides wellcalibrated uncertainty estimates. However, we emphasize that this model is not a substitute for site-specific predictions. We explicitly recommend that engineers treat the model outputs as preliminary guidance and supplement them with local testing whenever feasible. In particular, when applying the model to geologically distinct or sensitive conditions (e.g., structured or soft clays), incorporating site-specific data is mandatory for reliable analysis. Acknowledgements The authors would like to thank the members of the TC304 Committee on Engineering Practice of Risk Assessment & Management of the International Society of Soil Mechanics and Geotechnical Engineering for developing the database 304dB used in this study and making it available for scientific inquiry. We also wish to thank Dr. Jianye Ching for contributing this CLAY/10/7490 database to the TC304 compendium of databases.

References Beesley, M. E., & Vardanega, P. J. (2020). Parameter variability of undrained shear strength and strain using a database of reconstituted soil tests. Canadian Geotechnical Journal, 57(8), 1247-1255. Belkhatir, M., Arab, A., Della, N., & Schanz, T. (2012). Experimental study of undrained shear strength of silty sand: effect of fines and gradation. Geotechnical and Geological Engineering, 30, 1103-1118. Bozorgzadeh, N., Harrison, J. P., & Escobar, M. D. (2019). Hierarchical Bayesian modelling of geotechnical data: application to rock strength. Géotechnique, 69(12), 1056-1070. Chian, S., Bi, J., Halai, H., Lee, N., & Leung, C. (2023). Filling missing soil data using machine learning. In Smart Geotechnics for Smart Societies (pp. 372-377). CRC Press. Ching, J., & Phoon, K.-K. (2014a). Correlations among some clay parameters — the multivariate distribution. Canadian Geotechnical Journal, 51(6), 686-704. https://doi.org/10.1139/cgj-2013-0353 Ching, J., & Phoon, K.-K. (2014b). Transformations and correlations among some clay parameters — the global database. Canadian Geotechnical Journal, 51(6), 663685. https://doi.org/10.1139/cgj-2013-0262 Ching, J., Phoon, K.-K., & Chen, C.-H. (2014). Modeling piezocone cone penetration (CPTU) parameters of clays as a multivariate normal distribution. Canadian Geotechnical Journal, 51(1), 77-91. D’Ignazio, M., Phoon, K.-K., Tan, S. A., & Länsivaara, T. T. (2016). Correlations for undrained shear strength of Finnish soft clays. Canadian Geotechnical Journal, 53(10), 1628-1645.

Ding, D., & Loehr, J. E. (2019). Variability and bias in undrained shear strength from different sampling and testing methods. Journal of Geotechnical and Geoenvironmental Engineering, 145(10), 04019082. Griffiths, D., & Yu, X. (2015). Another look at the stability of slopes with linearly increasing undrained strength. Géotechnique, 65(10), 824-830. He, X., Cai, G., & Sheng, D. (2025). Indirect models for SWCC parameters: reducing prediction uncertainty with machine learning. Computers and Geotechnics, 177, 106823. Hu, L., Kasama, K., Wang, G., & Takahashi, A. (2025a). Assessing the influence of geotechnical uncertainty on existing tunnel settlement caused by new tunneling underneath. Tunnelling and Underground Space Technology, 155, 106189. Hu, L., Wang, G., & Kasama, K. (2025b). Probabilistic back-analysis of earthquakeinduced 3D landslide model parameters and risk assessment for secondary slide. Engineering Geology, 352, 108062. Jiang, S.-H., & Huang, J. (2018). Modeling of non-stationary random field of undrained shear strength of soil for slope reliability analysis. Soils and Foundations, 58(1), 185-198. Jong, S., & Ong, D. (2024). A Bayesian inference framework for geomaterial characterization and evaluation of complex soil-structure interactions. Computers and Geotechnics, 172, 106452. Kim, H.-S., & Ji, Y. (2022). Three-dimensional geotechnical-layer mapping in Seoul using borehole database and deep neural network-based model. Engineering Geology,

297,

https://doi.org/https://doi.org/10.1016/j.enggeo.2021.106489

106489.

Lan, P., Huang, J., Su, J., Zhu, S., Chen, J., Zhang, S., & Jiang, S.-H. (2025). A rapid unsaturated infiltration prediction method for slope stability analysis considering uncertainties: Deep operator networks. Engineering Geology, 347, 107886. Lan, P., Su, J., Zhu, S., Huang, J., & Zhang, S. (2024). Reconstructing unsaturated infiltration behavior with sparse data via physics-informed deep learning. Computers and Geotechnics, 168, 106162. Li, Z., Kasama, K., Hu, L., Yu, J., He, Y., & Yahiro, Y. (2025). The strength evaluation of cement-mixed soil considering mixing quality using image analysis and Bayesian

approach.

Computers

and

Geotechnics,

183,

107225.

https://doi.org/https://doi.org/10.1016/j.compgeo.2025.107225 Lunne, T., Berre, T., Andersen, K. H., Strandvik, S., & Sjursen, M. (2006). Effects of sample disturbance and consolidation procedures on measured shear strength of soft marine Norwegian clays. Canadian Geotechnical Journal, 43(7), 726-750. https://doi.org/10.1139/t06-040 Mun, W., Teixeira, T., Balci, M., Svoboda, J., & McCartney, J. S. (2016). Rate effects on the undrained shear strength of compacted clay. Soils and Foundations, 56(4), 719-731. Nguyen, H. C., Xu, H., Nazem, M., Sousa, R., Kowalski, J., & Zhao, Q. (2025). Deeplearning surrogate models for the stability of a wide rectangular tunnel. Computers

and

Geotechnics,

179,

106946.

https://doi.org/https://doi.org/10.1016/j.compgeo.2024.106946 Phoon, K.-K., & Kulhawy, F. H. (1999a). Characterization of geotechnical variability. Canadian Geotechnical Journal, 36(4), 612-624.

Phoon, K.-K., & Kulhawy, F. H. (1999b). Evaluation of geotechnical property variability. Canadian Geotechnical Journal, 36(4), 625-639. Piao, S., Huang, S., Wei, Y., Tan, J., & Ma, B. (2025). A Novel Hybrid Bayesian-GroupBased Machine Learning (HB-GML) Method for Predicting Uniaxial Compressive Strength (UCS) of Rock. Rock Mechanics and Rock Engineering, 1-27. Shogaki, T., & Kumagai, N. (2008). A slope stability analysis considering undrained strength anisotropy of natural clay deposits. Soils and Foundations, 48(6), 805819. Slifker, J. F., & Shapiro, S. S. (1980). The Johnson system: selection and parameter estimation. Technometrics, 22(2), 239-246. Soe, T. E. E., & Ukritchon, B. (2023). Three-dimensional undrained face stability of circular tunnels in non-homogeneous and anisotropic clays. Computers and Geotechnics, 159, 105422. Sun, R., Pan, B., & Duan, Q. (2023). A surrogate modeling method for distributed land surface hydrological models based on deep learning. Journal of Hydrology, 624, 129944. Tanaka, H. (2000). Sample Quality of Cohesive Soils: Lessons from Three Sites, Ariake, Bothkennar

and

Drammen.

Soils

and

Foundations,

40(4),

57-74.

https://doi.org/https://doi.org/10.3208/sandf.40.4_57 Ukritchon, B., & Keawsawasvong, S. (2020). Undrained stability of unlined square tunnels in clays with linearly increasing anisotropic shear strength. Geotechnical and Geological Engineering, 38, 897-915.

Wang, T., Chen, W., Li, T., Connolly, D. P., Luo, Q., Liu, K., & Zhang, W. (2023). Surrogate-assisted uncertainty modeling of embankment settlement. Computers and Geotechnics, 159, 105498. Xie, J., Huang, J., Zeng, C., Huang, S., & Burton, G. J. (2022). A generic framework for geotechnical subsurface modeling with machine learning. Journal of Rock Mechanics

and

Geotechnical

Engineering,

14(5),

1366-1379.

https://doi.org/https://doi.org/10.1016/j.jrmge.2022.08.001 Xu, B., Rong, Z., Pang, R., Tan, W., & Wei, B. (2024). A novel method for settlement imputation and monitoring of earth-rockfill dams subjected to large-scale missing

data.

Advanced

Engineering

Informatics,

62,

102642.

https://doi.org/https://doi.org/10.1016/j.aei.2024.102642 Yang, H., & Liu, Z. (2024). A fused sampling method integrating geotechnical and geophysical data for assessing three-dimensional soil liquefaction-induced damage

capacity.

Computers

and

Geotechnics,

166,

106024.

https://doi.org/https://doi.org/10.1016/j.compgeo.2023.106024 Ye, J., Mao, Y., Cheng, K., Tan, X., Du, B., & Chen, W. (2024). A spatiotemporal recurrent neural network for missing data imputation in tunnel monitoring. Journal

of

Rock

Mechanics

and

Geotechnical

Engineering.

https://doi.org/https://doi.org/10.1016/j.jrmge.2024.11.019 Yoshimine, M., Robertson, P., & Wride, C. (1999). Undrained shear strength of clean sands to trigger flow liquefaction. Canadian Geotechnical Journal, 36(5), 891906. Zhang, P., Yin, Z.-Y., & Jin, Y.-F. (2022). Bayesian neural network-based uncertainty modelling: application to soil compressibility and undrained shear strength prediction. Canadian Geotechnical Journal, 59(4), 546-557.

Zhang, P., Yin, Z.-Y., & Sheil, B. (2023). Interpretable data-driven constitutive modelling of soils with sparse data. Computers and Geotechnics, 160, 105511. https://doi.org/https://doi.org/10.1016/j.compgeo.2023.105511 Zhang, S., Xiong, H.-B., Tong, C.-X., He, X., & Sheng, D. (2024). Modelling the evolving critical state of crushable soils and parameter optimization with an improved genetic algorithm (AIS–RCGA). Canadian Geotechnical Journal, 62, 1-19. https://doi.org/10.1139/cgj-2024-0131 Zhang, W., Wu, C., Zhong, H., Li, Y., & Wang, L. (2021). Prediction of undrained shear strength using extreme gradient boosting and random forest based on Bayesian optimization.

Geoscience

Frontiers,

https://doi.org/https://doi.org/10.1016/j.gsf.2020.03.007

12(1),

469-477.

Related documents

• Goal-Oriented Weighting of Reynolds-Stress Data for Learning Turbulence Models in Complex Flows#1105595
• Enhancing clinical reliability in pressure injury prediction: A conformal prediction approach with machine learning models.#196074
• Deep Learning Imputation of Missing Radius of Maximum Winds (Rmax) Values in Tropical Cyclone Best-Track Data#441807
• Explainable Machine Learning Using Sensor-Derived Biomechanical Features to Classify Elevated VALR-Related Loading Across Midsole Hardness Conditions in School-Aged Boys.#317770
• Embedding composite failure mechanics into neural networks: A physics-constrained framework for dual-criteria failure assessment of Type IV pressure vessels#246829