This study utilizes three performance metrics: the root mean square error (RMSE), the coefficient of determination (R ²), and the a ₂₀ index. The equations for RMSE, R ², and a ₂₀ are expressed in Equations 15–17, respectively: R M S E = 1 N t ∑ i = 1 N t Y p i − Y i 2 (15) R 2 = 1 − ∑ i = 1 N t Y p i − Y i 2 ∑ i = 1 N t Y p i − Y ¯ p i 2 (16) a 20 − i n d e x = m 20 M (17) where Y _i and Y _pi represent the actual and predicted values for the ith observation, respectively; N _t denotes the number of testing models, and Y̅ _pi is the mean of the predicted values. The optimal ML model is characterized by the lowest RMSE and the highest R ² and a ₂₀ index values. The parameter M denotes the total number of samples in the dataset, while m20 refers to the number of samples for which the ratio of the observed value to the predicted value lies within the range of 0.80–1.20. This range indicates a ±20% tolerance, commonly accepted in engineering practices as a benchmark for acceptable prediction accuracy. The a ₂₀, defined as the ratio m20/M, therefore provides an intuitive and practical measure of a model’s predictive performance (Asteris et al., 2024a; Asteris et al., 2024b). In the case of a perfectly accurate model, the a ₂₀ would equal 1, indicating that all predictions fall within the acceptable range. The value of the proposed a ₂₀ lies in its interpretability and relevance to engineering applications, as it directly reflects the proportion of predictions that closely match the corresponding experimental results within a tolerable margin of error.

Once the dataset has been prepared and the ML technique selected, the next crucial step is defining the parameters of the model, which significantly influence its performance. In this study, hyperparameter optimization is carried out using a combination of grid search and K-fold CV to minimise the risk of overfitting. Initially, potential parameter ranges are defined as grids. The model is then trained iteratively, testing all possible parameter combinations, with performance evaluated using K-fold CV. This approach ensures a reliable assessment of predictive accuracy while reducing bias from the random division of training and testing data. The dataset is split into K equally sized subsets, and the model undergoes K iterations, where K−1 subsets are used for training and the remaining subset for validation. The final model performance is determined by averaging the results across all K iterations. In this study, a 10-fold CV approach is implemented, as it is a widely recognised standard that partitions the data into ten groups, effectively reducing the risk of overfitting.

Saltelli et al. (2008) proposed a variance-based SA method to evaluate how changes in model input values affect the corresponding output. This technique examines interactions between input variables and the output factor by keeping all input parameters constant except for one, which is systematically varied (Liu et al., 2020).

Friedman (2001) introduced the PDP to examine the marginal effect of a specific parameter on the output by depicting the average outcome values across different parameter settings. The PDP is useful for identifying the relationship between a feature and the target variable. The partial dependence, represented as f _S , for a subset of features x _S , is defined as expressed in Equation 20: f S x S = E X C f x S , x C = ∫ f x S , x C d P x C (20) where x _S represents the factors considered for the PDP, while x _C denotes the complementary factors. A PDP can be generated for a dataset {X _i , i = 1, … ,n}, as expressed in Equation 21: f ¯ S x S = 1 n ∑ i = 1 n f x S , x i C (21)

A critical assumption of PDP is the independence of input parameters. When features are highly correlated, the analysis may involve artificial data samples that are unrealistic, resulting in significant bias in the estimated effect of the feature.

ALE is an unbiased alternative to PDP that represents the average effect of a feature within a ML algorithm (Apley and Zhu, 2020). Unlike PDP, which averages f _S (x _S ) over all features and can produce biased results when features are highly correlated, ALE mitigates this issue by focusing on changes in predictions and constraining data to specific grids. ALE averages the changes in predictions by isolating the influence of correlated features, simplifying complex models by considering only one or two factors at a time. ALE can be expressed as Equation 22: f x s , A L E x S = ∫ z 0 , 1 x S E X C / X S f ^ S x S , x C x S = z S d z S − c = ∫ z 0 , 1 x S ∫ x C f ^ S z S , x C P x C z S d x C d z S − c (22)

where c is a fixed constant, f ^ S x S , x C = δ f ^ x S , x C δ x S denotes the local effects of x _S on f ^ . at (x ₁, x _S ); z _0,1 represents a value smaller than the minimum observation, and P (x _C |x _S ) indicates the density. ALE plots are centred around zero, making them particularly effective for correlated data and enhancing visualization. However, the choice of grid or interval can influence the plots and obscure data variability.

Lundberg and Lee (2017) introduced the SHAP method, a technique grounded in conditional expectations and game theory, for assessing model predictions. SHAP is used to analyse the influence of individual input features on each output. Broadly, SHAP helps rank features based on their contribution to interaction effects. SHAP employs an additive feature attribution approach to create an interpretable model. The resulting model, represented as a linear function, is the sum of the actual contributions associated with each parameter. The interpretable framework is formulated as expressed in Equation 23: f x = g x ′ = ϕ 0 + ∑ i = 1 M ϕ i x i ′ , (23) where x = (x ₁, x ₂, … , x _p ) represents the M input variables; and x’ _I denotes the simplified inputs. Using a mapping function x = hx (x’), the parameter x’ is transformed into x. Additionally, ϕ ₀ and ϕ _i represent a fixed value and the contribution of each parameter, respectively.

PDP, ALE, and SHAP each offer distinct strengths in interpreting complex ML models, as presented in Table 1. PDP is easy to understand and effective for identifying overall trends but may misrepresent effects when features are correlated. ALE addresses this by providing unbiased local effects even with correlated inputs, though its accuracy depends on appropriate interval selection. SHAP offers the most comprehensive insights, attributing both global and local feature importance while capturing interaction effects, albeit at a higher computational cost. Their combined use provides a robust interpretability framework that balances simplicity, accuracy, and depth of insight.

Effective ML model development requires rigorous data preprocessing and systematic hyperparameter tuning to ensure robust and generalizable performance. In this study, all input features were normalized with respect to the bolt diameter d ₀ or relevant physical properties to ensure comparability across different test cases and avoid scale-related bias in models sensitive to input magnitude. Specifically, geometric inputs such as end distance, edge distance, and pitch dimensions were divided by d ₀ to create dimensionless variables. This representation aligns with established bearing and shear-out design formulations and enables the model to generalise across connections with different bolt sizes. From a numerical standpoint, using normalized features mitigates scale-related bias in algorithms that are sensitive to feature magnitude, improving optimization stability and reducing the risk that variables with larger numerical ranges dominate the learning process. If raw, non-normalized geometric dimensions were used instead, the model would be forced to infer the relevant geometric ratios implicitly, potentially reducing generalisation quality across bolt diameters and degrading performance for scale-sensitive algorithms, without offering a clear advantage in predictive accuracy. This standardization improves convergence and stability across ML algorithms. The dataset compiled from 443 experimental results was reviewed for completeness. Since the source data was manually extracted from published studies, any entries with incomplete feature sets were excluded from analysis to avoid bias introduced by imputation. As a result, the dataset used in this study contained no missing values.

Tuning hyperparameters is a crucial process in refining ML models to achieve optimal performance. It entails modifying the settings of the model to obtain the best results for a given dataset. A grid search method is often used for this purpose, where a range of hyperparameter values is systematically explored to find the combination that provides the highest performance. 10-fold CV is incorporated to ensure the selected hyperparameters are reliable and to prevent overfitting. In this technique, the dataset is divided into 10 roughly equal-sized subsets. The model is trained and evaluated 10 times, with each iteration selecting a different subset as the validation set while the remaining subsets are used for training. This approach systematically explored combinations of hyperparameters and assessed their impact on model performance across multiple data splits. The robustness of this methodology lies not only in its breadth of parameter search but also in the multi-metric evaluation strategy used to determine the best-performing models. This approach guarantees that the performance of the model is assessed across various data splits, offering a more dependable measure of its ability to generalize. For each model and parameter combination, two primary metrics were computed: RMSE and R ². Both metrics were computed for each fold, and their average values across the 10 folds were used to evaluate each hyperparameter set. The following criteria were applied:

Primary selection metric: The optimal configuration minimized the average RMSE across the 10 folds, ensuring minimal prediction error.

Different models were tuned with specific parameter ranges, and selections reflected trade-offs between bias and variance. Table 3 presents the optimal hyperparameter values for each ML algorithm, which were selected based on their ability to maximize predictive accuracy while balancing bias and variance.

The sensitivity indices for the first-order (S _i ) and total effect (S _T ) of the input variables are determined using the regression models, as shown in Table 4 and Figure 5. These indices illustrate the contributions of the input variables to the variance of the model output. It can be seen that the variables e ₁/d ₀ and e ₂/d ₀ consistently show high contributions to the variance of the model, especially for S _i , across most models. The variable fu/fy generally has the lowest sensitivity indices, indicating a negligible influence on the variance of the model. The sensitivity of N _r varies across models, being significant in some cases (e.g., for SVR and RF) but less so in others. The LR and RR models exhibit similar sensitivity patterns, with e ₁/d ₀ and e ₂/d ₀ contributing the most to the output variance (S _i ≈ 0.35–0.39). For the SVR model, sensitivity indices are more evenly distributed among variables, with p ₂/d ₀ and N _r showing significant contributions. In KNN, S _T indicates high combined interactions (S _T values exceed 1.9). DT shows relatively even sensitivity, with moderate values for e ₁/d ₀ and e ₂/d ₀. RF highlights strong contributions from e ₁/d ₀, e ₂/d ₀, and interactions among other variables (S _T ≈ 1.1). For the AdaBoost model, variables e ₁/d ₀ and e ₂/d ₀ dominate both S _i and S _T . XGBoost, CatBoost, and LightGBM generally show high interaction effects, with S _T often exceeding 1. The e ₂/d ₀ variable consistently has the highest impact. The total variance contribution for first-order indices (∑S _i ) is less than or close to 1 in most cases, indicating that individual variables largely explain the model variance. For total effects (∑S _T ), values exceed 1 for certain models (e.g., KNN, LightGBM), suggesting significant interaction effects among the input variables.

In general, this SA is essential for identifying critical inputs and their interactions, guiding feature selection or model improvement strategies. The consistently high indices of e ₁/d ₀ and e ₂/d ₀ confirm them as primary design drivers, which can therefore be prioritised in the analyses. By contrast, the systematically low first-order and total-effect indices of p ₁/d ₀ indicate that it is a secondary descriptor; it can be fixed or omitted in simplified surrogate models with limited expected impact on predictive performance, although it is retained here for physical completeness. From a model-development standpoint, the pronounced gaps between S _i and S _T ​ for KNN, RF, XGBoost, CatBoost, and LightGBM confirm the presence of strong interaction effects, justifying the adoption of tree-based ensemble models that can capture higher-order interactions and guiding both the selection of CatBoost as the final surrogate and the design of the two-variable contour analyses.

Figure 6 presents a set of regression plots comparing predicted and measured values across different ML models. Each subfigure contains four plots:

Histogram of Residuals: Displays the distribution of residuals (measured-predicted), showing the error spread;

The diagonal line in the predicted vs. measured scatterplot serves as a reference to evaluate prediction accuracy. Points closer to the line indicate better performance. The spread in the residuals vs. predicted values plot reveals whether the model has systematic bias (e.g., over- or under-prediction for certain ranges). Differences in the residual and prediction histograms reflect variations in model accuracy and distribution alignment. The figure highlights how different models behave in terms of prediction accuracy, bias, and variance. Color coding (e.g., red and blue points) differentiates between subsets of data or categories, indicating performance consistency across different data types.

From visual inspection, XGBoost, CatBoost, and LightGBM appear to perform better. The residual histograms for these models are narrower and more centred around 0 compared to others. In the residual vs. predicted plots, no clear patterns or trends are visible, indicating reduced systematic error. The scatterplot of predicted vs. measured values shows that points are closer to the diagonal line, reflecting better prediction accuracy. LR and KNN seem to perform the worst, as they have wider residual distributions and systematic patterns in the residuals vs. predicted plots, suggesting poorer performance.

In general, XGBoost, CatBoost, and LightGBM leverage gradient boosting techniques, which combine the strengths of multiple weak learners (e.g., decision trees) to iteratively minimize error. When properly tuned, these models are robust to overfitting and effectively handle complex relationships and interactions within the data. In contrast, simpler models like LR or KNN may fail to capture non-linear relationships or interactions, leading to poorer performance.

In conclusion, the gradient boosting models (XGBoost, CatBoost, and LightGBM) perform better due to their superior ability to capture complex patterns, resulting in tighter residual distributions, minimal bias, and high alignment between predicted and measured values.

Figure 7 shows Taylor Diagrams used to evaluate and compare the performance of different ML models during (a) the training phase and (b) the testing phase. Taylor diagrams provide a concise summary of how well a model matches observations by plotting three statistics simultaneously:

Standard Deviation (horizontal axis): This indicates the variability of the predicted data relative to the observed data. Ideally, the standard deviation of the model (represented by the markers) should be close to that of the reference (red star).

Figure 7a shows the training phase. Most models cluster near the reference, indicating high correlation and standard deviation close to the reference during training. Certain models, such as XGBoost and LightGBM, performed better by showing higher correlation and smaller CRMSE compared to others. Figure 7b presents the testing phase. The models show more variation in their correlation and standard deviation compared to the training phase. The scatter suggests that some models generalize better to unseen data (higher correlation and smaller CRMSE), while others may overfit the training data. It can be seen that the CatBoost model outperforms the others, as it maintains relatively high correlation coefficients compared to the other models, indicating good generalization on unseen data. Its standard deviation remains closer to the reference compared to the other models, and its CRMSE value (distance to the red star) is smaller than those of the other models, suggesting lower prediction errors on the testing data.

Table 5 evaluates the performance of various ML-based regression models for predicting f(x). The metrics used to compare the models are Average R ², Average RMSE, and a ₂₀ index. Both LR and RR have low R ² values (∼0.286–0.320), indicating poor explanatory power. The RMSE is relatively high (0.554–0.556), suggesting less accurate predictions. Both SVR and KNN show significantly higher R ² values (∼0.685–0.686), indicating better predictive accuracy. KNN has the lowest RMSE (0.334), suggesting it provides the most precise predictions among the models. DT presents moderate performance with R ² = 0.320, RMSE = 0.526, and a ₂₀ = 0.805. RF shows better performance (R ² = 0.703, RMSE = 0.363, and a ₂₀ = 0.827). CatBoost has the highest R ² (0.711), showing the best explanatory power among all models. XGBoost and CatBoost have relatively low RMSE values (0.359), indicating good predictive precision. AdaBoost and LightGBM perform moderately, with lower R ² values and higher RMSE compared to CatBoost and XGBoost. As a result, CatBoost has the highest R ² (0.711), indicating it explains the target variable variance the best.

Among all models, the CatBoost model achieves the highest average R ² and one of the lowest RMSE values. Given that the normalised bearing capacity ranges from 0.28 to 5.44 with an average of 2.86, this RMSE corresponds to an error of approximately 12%–13% of the mean response and about 7% of the full range. Furthermore, the a ₂₀ index of 0.872 indicates that around 87% of predictions lie within ±20% of the experimental values. Although this does not eliminate variability, it represents a substantial improvement over simple linear models and is comparable to the scatter typically observed between design-code predictions and test results for bolted connections. On this basis, CatBoost is selected as the reference surrogate model for the interpretability and parametric analyses in the following sections, with the understanding that it is intended as a complementary, not stand-alone, design tool.

PDPs illustrate the relationship between a specific feature and the predicted outcome of a ML model, while keeping all other features fixed, as shown in Figure 9. Each subplot corresponds to a different input variable affecting the output of the model. The x-axis represents the values of the specific input variables, and the y-axis represents the Partial Dependence, or how the predictions of the model vary with changes in that specific variable. The black ticks indicate the distribution of data points (or unique values) for that variable in the training dataset. Denser areas suggest more frequent values in the data.

The observations for each variable are as follows:

e ₁/d ₀: Partial dependence increases almost monotonically, indicating a positive relationship between e ₁/d ₀ and the model output.

Variables like N _r and e ₁/d ₀ seem to have clearer trends, while fu/fy might require further investigation due to its variability.

ALE represents the average effect of the input variable on the prediction of the model, accounting for interactions with other variables, as shown in Figure 10. Positive ALE values suggest that the variable contributes positively to the prediction, while negative ALE values indicate a negative contribution. The x-axis represents the range of the respective input variable, normalized between 0 and 1.

Each plot corresponds to the effect of a specific input variable on the predictions of the model:

e ₁/d ₀: Shows a positive correlation; as e ₁/d ₀ increases, the ALE increases significantly, indicating a strong positive effect.

Tick marks on the x-axis represent the distribution of data points for the respective variable, indicating where the model is well-supported by data.

Figure 11a presents the SHAP summary plot, where each point corresponds to a Shapley value associated with a specific parameter. The plot arranges samples into rows, each containing an equal number. The Shapley values are plotted along the x-axis, while the input variables are listed on the y-axis in order of importance, with the most significant ones at the top. Samples with identical SHAP values for a factor are spread horizontally. The variable values are color-coded, with high values in red and low values in blue. Red indicates values that increase the SHAP score and the associated estimate.

The data suggests that increasing the normalized pitch perpendicular to loading and the number of bolt rows leads to a decrease in the SHAP value, and thus a reduction in the normalized bearing capacity. Conversely, increasing the normalized end distance, normalized edge distance, or the ultimate-to-yield stress ratio of the critical steel plate results in an increase in the normalized bearing capacity. Figure 11a ensures an even distribution of points across the rows. In Figure 11b, the global significance factor is represented by the average absolute SHAP value for each factor. SHAP analysis identifies normalized end distance and edge distance as the most significant factors, aligning with the variable importance results from the CatBoost model.

Figure 12 presents the SHAP dependence graph through scatter plots that show the SHAP value of one parameter against others. The colours in the graph signify how interactions with different variables affect the values on the horizontal axis, with many of these interactions being non-linear. In Figure 12a, the SHAP values increase as e ₁/d ₀ increases, indicating that higher values of e ₁/d ₀ have a more positive contribution to the predictions of the model. The colour gradient represents N _r , where higher N _r values tend to amplify the impact of e ₁/d ₀ on the predictions. In Figure 12b, the SHAP values generally increase as e ₂/d ₀ increases, although the relationship appears to level off or saturate at higher values. The colour gradient corresponds to e ₁/d ₀, showing that lower values of e ₁/d ₀ are associated with smaller SHAP contributions, while higher e ₁/d ₀ amplifies the SHAP values for e ₂/d ₀. In Figure 12c, for most of the range of p ₁/d ₀, the SHAP values are clustered near zero, but there is an increasing trend at higher p ₁/d ₀ values, indicating a positive contribution to predictions. The colour gradient for e ₁/d ₀ shows that higher e ₁/d ₀ values correspond to more positive SHAP values, suggesting an interaction between p ₁/d ₀ and e ₁/d ₀.

In Figure 12d, the SHAP values are distributed across a wide range, with noticeable clusters, suggesting non-linear relationships. There are regions where p ₂/d ₀ has little to no effect (SHAP values near zero) and others where it contributes positively or negatively. The feature e ₁/d ₀ interacts with p ₂/d ₀, with higher e ₁/d ₀ amplifying the SHAP values. In Figure 12e, the SHAP values exhibit a slight overall decreasing trend with N _r , indicating that larger N _r values generally contribute negatively to the predicted normalized bearing capacity. The colour gradient for e ₁/d ₀ further shows that this negative influence is more pronounced for lower e ₁/d ₀ and can be partially mitigated when e ₁/d ₀ is large. In Figure 12f, the SHAP values are scattered without a strong linear trend, indicating a complex or weak relationship between fu/fy and the predictions of the model. Certain clusters show positive or negative contributions. The colour bar for p ₂/d ₀ suggests that the interaction between fu/fy and p ₂/d ₀ is significant, with higher p ₂/d ₀ values leading to more pronounced positive or negative SHAP contributions.

In general, the behaviour of SHAP values indicates that e ₁/d ₀ and e ₂/d ₀ have the strongest and most consistent positive impact on predictions. Features like p ₂/d ₀ and fu/fy exhibit non-linear or scattered relationships, implying complex interactions within the model. The colour gradients (interaction effects) reveal secondary relationships between the displayed feature and auxiliary features, emphasizing the multi-feature dynamics of the model.

Figure 13 shows a SHAP force plot, which visualizes how individual feature values influence the prediction of the model. The plot represents how different input variables contribute to pushing the prediction away from its base value (the expected value of the model output across the dataset). The red section represents features that push the prediction higher. The blue section represents features that push the prediction lower. The grey vertical line in the middle indicates the base value of the prediction. 1.73 is the expected value of the model output before considering the influence of individual feature values. Features shown in red, including fu/fy, p ₂/d ₀, N _r , e ₂/d ₀, increase the output from the base value. Features shown in blue, including e1/d0 decrease the output from the base value. In general, the feature interpretations are as follows:

fu/fy = 0.8485: A feature that pushes the output upwards.

Figure 14 presents a SHAP decision plot. The order of features on the y-axis suggests their relative importance. Higher-ranked features, such as e ₁/d ₀, e ₂/d ₀, have a larger overall impact on the output. Features lower on the axis, such as p ₁/d ₀, contribute less overall, indicating they may be less influential. The blue-to-red gradient shows how feature values relate to their SHAP impact. For e ₁/d ₀, higher values (red) mostly result in positive SHAP values, meaning high e ₁/d ₀ pushes the output of the model to increase. For e ₂/d ₀, similarly, higher values tend to have a positive SHAP impact, but some variability suggests a non-linear relationship. Low feature values (blue) generally have a negative or neutral SHAP impact, pulling the output lower. Features like p ₂/d ₀ and N _r show tightly clustered SHAP values around 0. This could mean their contributions are smaller and more uniform across the dataset. Conversely, features like e ₁/d ₀ have a wider spread, indicating a more variable contribution to the output. In SHAP plots, overlapping colours and values suggest potential feature interactions. For example, the spread in N _r and fu/fy might indicate interactions with other features, where their SHAP values depend on the presence of other variables. The model output appears to centre near 2, suggesting that for most samples, predictions are not strongly influenced (in either positive or negative directions) by individual features. Features like e ₁/d ₀ and e ₂/d ₀ have notable impacts, pulling the output above or below this centre.

Figure 15 shows a SHAP waterfall plot, which is commonly used to explain the contribution of different input variables to a ML prediction for a single instance. The top-left part of the plot states f(x) = 1.729, which is the predicted output for the given instance. The bottom label shows E [f(X)] = 1.751, which represents the expected prediction across all instances in the dataset. The largest impact comes from e ₁/d ₀, which decreases the prediction by 0.53. The second major factor is e ₂/d ₀, which increases the prediction by 0.37. Smaller contributions are made by p ₂/d ₀ (+0.07), N _r (+0.03), fu/fy (+0.03), and p ₁/d ₀ (+0.01).

Figure 16 shows a series of contour plots that represent the output of a two-variable exponential decay function. Each subplot displays the relationship between two input variables (on the axis) and the corresponding output of the decay function (indicated by the contour levels and the colour gradient). The x- and y-axis in each subplot represent different parameter pairs or input variables. Each subplot investigates how these variable pairs interact and influence the output of the function. The contour lines connect points of equal output values, providing a visual representation of the decay behaviour. The colour map indicates the magnitude of the output, where typically:

Red shades signify higher output values.

The objective is likely to analyse and fit a two-variable exponential decay model to some data. The plots help visualize how different parameter combinations affect the output, which may be useful for identifying trends, dependencies, or sensitivity of the function to specific variables. Peaks (red areas) and troughs (blue areas) in the contour plots highlight regions of high and low output values, respectively. Smooth transitions or sharp gradients in colour and contour lines indicate the rate of change in the output with respect to changes in the input variables.

Figure 16a illustrates the relationship between the parameters e ₁/d ₀, e ₂/d ₀, and the normalized variable F _max/f _u ndt. As e ₁/d ₀ and e ₂/d ₀ increase, F _max/f _u ndt also increases, indicating a direct correlation between these variables. The gradient of the contour lines suggests a nonlinear interaction, with a steeper increase in F _max/f _u ndt for higher values of e ₂/d ₀. Figure 16b represents the variation of the parameter F _max/f _u ndt with respect to e ₁/d ₀ and p ₁/d ₀. The contours show that F _max/f _u ndt increases as both e ₁/d ₀ and p ₁/d ₀ increase, with a steeper gradient in regions of higher e ₁/d ₀. The behaviour suggests a nonlinear relationship, where changes in p ₁/d ₀ influence F _max/f _u ndt more strongly for larger e ₁/d ₀ values. Figure 16c illustrates the relationship between two normalized parameters, e ₁/d ₀ and p ₂/d ₀, and their influence on the response parameter F _max/f _u ndt. As e ₁/d ₀ increases, the response parameter F _max/f _u ndt grows significantly, particularly for higher values of p ₂/d ₀. The behaviour suggests a nonlinear dependency, with a steep gradient for lower e ₁/d ₀ transitioning to a smoother variation at higher values.

Figure 16d shows that the variable F _max/f _u ndt increases with both e ₁/d ₀ and N _r , indicating a direct relationship between these parameters and the maximum normalized force. The gradient is steeper along the e ₁/d ₀ axis compared to the N _r axis, suggesting that e ₁/d ₀ has a more significant influence on F _max/f _u ndt. Regions with higher F _max/f _u ndt values are concentrated towards the top-right corner, emphasizing the combined effect of high e ₁/d ₀ and N _r . Figure 16e displays the relationship between the parameters e ₁/d ₀, fu/fy, and the normalized value F _max/f _u ndt. As e ₁/d ₀ increases, F _max/f _u ndt generally rises, indicating a strong dependency on this ratio. Similarly, higher values of fu/fy correspond to increased F _max/f _u ndt, showcasing its influence on the performance of the system. Figure 16f shows the parametric relationship between the variables p ₁/d ₀, e ₂/d ₀, and F _max/f _u ndt, represented as contour levels. The behaviour indicates that as e ₂/d ₀ and p ₁/d ₀ increase, F _max/f _u ndt rises steadily, transitioning from lower values to higher values. The gradient suggests a nonlinear correlation, where F _max/f _u ndt is more sensitive to changes in e ₂/d ₀ at higher values of p ₁/d ₀.

Figure 16g shows that as e ₂/d ₀ and p ₂/d ₀ increase, F _max/f _u ndt generally increases, suggesting a direct influence of these parameters on the force ratio. The steep gradients at lower values of e ₂/d ₀ and p ₂/d ₀ indicate higher sensitivity of F _max/f _u ndt to changes in these regions. Figure 16h illustrates that as e ₂/d ₀ and N _r increase, the value of F _max/f _u ndt rises, transitioning from blue (low values) to red (high values), indicating a strong positive correlation. The plot reveals that F _max/f _u ndt grows more rapidly at higher values of both e ₂/d ₀ and N _r , suggesting nonlinear behaviour in this region. Figure 16i shows that as e ₂/d ₀ increases, the normalized F _max/f _u ndt also increases significantly, indicating a stronger structural response with higher ratios of eccentricity to diameter. Similarly, higher values of fu/fy correlate with increased F _max/f _u ndt, suggesting that higher ultimate-to-yield strength ratios lead to improved maximum force response.

Figure 16j indicates a sharp transition from higher values to lower values as p ₁/d ₀ increases beyond a threshold, particularly at lower p ₂/d ₀ values. The behaviour stabilizes at consistently lower values for larger p ₁/d ₀ and p ₂/d ₀, suggesting diminishing sensitivity in these regions. Figure 16k shows that as p ₁/d ₀ and N _r increase, the value of F _max/f _u ndt generally increases, transitioning from blue to red regions. The gradient patterns suggest a nonlinear interaction where both parameters significantly influence the outcome, with higher values concentrated in the upper-right corner. Figure 16l shows that the highest values (F _max/f _u ndt ≈ 1.8) occur at very low values of p ₁/d ₀ and slightly higher fu/fy values, suggesting sensitivity to these parameters. As p ₁/d ₀ increases, F _max/f _u ndt decreases significantly, indicating reduced structural performance in these regions.

Figure 16m indicates that the feature F _max/f _u ndt shows a significant gradient along the p ₂/d ₀ axis near zero, with values exceeding 2, indicating sensitivity to changes in p ₂/d ₀ at lower values. Beyond a certain threshold of p ₂/d ₀, F _max/f _u ndt stabilizes around 1 for most of the N _r range, suggesting a reduced influence of p ₂/d ₀ at higher values. Figure 16n shows that the strong gradient near p ₂/d ₀ = 0 suggests significant sensitivity in this region, with the ratio F _max/f _u ndt increasing rapidly. For larger values of p ₂/d ₀, the ratio stabilizes, and fu/fy appears to have a weaker influence overall, with the region dominated by a uniform lower value. Figure 16o shows that as N _r increases, F _max/f _u ndt decreases, indicating an inverse relationship. Similarly, as fu/fy increases, F _max/f _u ndt also increases, suggesting a direct relationship between these parameters.

From a practical design standpoint, the interpretable CatBoost model and the associated interpretability methods provide several qualitative recommendations for detailing double-shear, bearing-type bolted connections within the investigated domain. First, the contour plots confirm that increasing the normalized end and edge distances systematically enhances the normalized bearing capacity, with the highest values concentrated in regions where both ratios are simultaneously large. This observation highlights the benefit of providing end and edge distances that exceed minimum code requirements rather than merely satisfying them. Second, the transverse pitch ratio should not be selected too small: the response surfaces show a steep increase in normalized bearing capacity as the transverse pitch ratio increases from very low values, followed by a plateau at moderate spacings. This indicates that excessively tight transverse spacing should be avoided, while very large transverse pitch ratios offer limited additional benefit. Third, the longitudinal pitch ratio exhibits a non-monotonic influence. Relatively small longitudinal pitch ratios, when combined with suitable material properties, yield the highest normalized capacities, whereas excessively large longitudinal pitch ratios lead to a marked reduction in normalized bearing capacity. Finally, the effect of the number of bolt rows is secondary to, and interacts with, the geometric parameters: increasing the number of bolt rows can be beneficial when accompanied by sufficient end and edge distances and appropriate pitch values, but it does not guarantee higher normalized capacity on its own. These insights are intended as qualitative guidance to support detailing decisions and rapid design exploration using the proposed ML model. They remain constrained by the range of geometries and material properties represented in the experimental dataset and should be applied in conjunction with existing code provisions and checks for other governing failure modes.