Reliability prediction of the axial ultimate bearing capacity of piles: A hierarchical Bayesian method

Abstract

The sea reclamation is one of the efficient ways to alleviate the shortage of land resources due to population growth, and the corresponding axial ultimate bearing capacity of piles has become one of the critical factors for evaluating the performance of the soil layer reclamation work. Many models are used to analyze the testing data. However, these models cannot describe the mean population bearing capacity and unit-to-unit variation simultaneously, and they cannot give the reliability of predicting the axial ultimate bearing capacity of piles. Thus, they are rarely used in practice. In this article, we propose a mixed-effects model, which could overcome the drawback of the models in the literature. A hierarchical Bayesian framework is developed to estimate the model parameters using Gibbs sampling. The proposed model is applied to a real pile dataset collected in silt-rock layer area, and we predict the mean axial bearing capacities under different reliability levels.

Introduction

With the rapid development of economy in coastal area and the increasing pressure of population growth, the contradiction between land resources and space resources is becoming more and more serious. In order to solve the problem of land deficit in the coastal area, Japan, Holland, Singapore, and China mainly use the way of sea reclamation, which can be found in Zhuang et al., ¹ Swinbanks, ² Fang et al., ³ and Lee. ⁴ The typical soil layer of reclamation work is shown in Figure 1. The geological conditions of such a case are very complicated, and general piles encounter construction difficulties in practice. Due to the drilled shaft technology, it has been widely used in many kinds of soil layers to meet the requirements of different bearing capacity, such as complicated formation of coastal areas ⁵ and silt-rock fill layer. ⁶ However, the compression of the silt-rock fill layer at the layer interface is extremely different. Due to the disturbance of the silt layer in the construction, the mechanical properties change greatly. In addition, because of the large water content, and the influence by the sea ebb and tide, the soil properties are obviously different from the conventional type. Therefore, a lot of geotechnical engineering problems could come out consequently.

Figure 1.

Slit-rock fill layer.

The ultimate bearing capacity is considered to be one of the significant factors that governs the design of pile foundations, denoted as $Q_u$ in this article. What engineering designers have always concerned is how to determine $Q_u$ reasonably, and how to make full use of the technical and economic benefits of the piles. In the past decades, numerous theoretical and experimental procedures have been proposed to predict axial bearing capacity of the piles, which can be found in O’Neill et al., ⁷ Peng et al.,^8,9 Randolph, ¹⁰ Tajimi, ¹¹ Xu and collegues,^12,13 and Zhu et al.^14,15 However, accurate evaluation of pile bearing capacity and certain interpretation of pile load transfer mechanism are still far from being accomplished due to the complexity of the problem. Still, the most direct and reliable approach to determine the pile axial bearing capacity is statically loading the pile until up to its failure based on some construction industry standard (CIS). ¹⁶ The CIS can refer to JGJ94-2008 of China ¹⁷ and FHWA NHI-10-016 of America. ¹⁸

The main drawbacks of following CIS are high cost and time-consuming. Therefore, most of the $Q-s$ curves obtained from the experiment are not complete, where Q represents the axial loading and s represents the settlement of the pile cap. Thus, it is difficult for engineers to determine the $Q_u$ value directly and accurately. In fact, the incomplete $Q-s$ curve can also partly reflect the influence of the pile load transfer mechanism, the size of the pile, the depth of soil, and soil mechanical properties on the pile foundation $Q_u$ , and it has a certain trend. Thus, many researchers used the incomplete $Q-s$ curve rather than the theory of pure mechanics to establish scientific and effective models to predict the value of $Q_u$ . In previous research, the $Q-s$ curve is usually assumed to be a certain known function, such as the exponential function^19–23 and the hyperbolic function model.^22,24–26 Then, based on the experimental data obtained, the least-squares method, variable scale optimization, and other mathematical methods are utilized to estimate soil parameters. Finally, the estimated function is used to predict $Q_u$ , according to some axial ultimate bearing capacity criterion. Obviously, the accuracy of predicting the ultimate loading is related to the mathematical model, and the estimation of the ultimate loading of this kind of model can be reached only if the displacement tends to infinity. Actually, these methods are very limited in practice. The reasons are as follows: firstly, the number of measurements is finite and subjected to the construction budget; secondly, the displacement should not exceed the ultimate threshold too much according to CIS, otherwise, the effect will be catastrophic for the whole engineering experiment; finally, it is meaningless for engineer that the displacement tends to infinity. Luo et al. ²⁷ viewed the settlement of the pile cap as the general time, and constructed G(1, 1) prediction model based on the first-order dynamic differential equations via the gray system theory proposed by Deng. ²⁸ Similar works can be found in Peng et al. ²⁹ When the method is used to solve the problem, the accuracy and reliability of the results mainly depend on the experts’ empirical knowledge. For other methods, Shahin ³⁰ and Momeni et al. ³¹ applied neural network model to predict the ultimate bearing capacity based on the dynamic test. Ju and Lee, ³² Lee and Ju, ³³ and Armaghani et al. ³⁴ studied the pile ultimate bearing capacity under the static test via neural network model. However, the results based on these methods still need to be verified in practice. Deng and Lu ³⁵ concluded that the accuracy of the prediction based on various models is not only related to the mathematical model but also closely related to the properties of the soil layer.

In addition, all the above models are only effective to analyze the bearing capacity of a single pile and cannot obtain the mean bearing capacity of the piles in the area. In fact, the selected piles from the project site may be very different due to unobservable endogenous factors, such as initial variations in raw materials, and due to exogenous factors, such as different structures of silt layers. The unit-to-unit variation among piles make the previous models fail to analyze the data well, and thus obtaining the bearing capacities of the sampled plies does not make any sense for practice. Random-effects models have been proved useful in dealing with these unobserved heterogeneities and could describe the mean population bearing capacity and the unit-to-unit variation of the piles simultaneously. Thus, for better analyzing the data, a mixed-effects model is proposed for the pile data in this article.

The data

General introduction of the project

The project site is located in a coastal city in China, with its border connecting to the city on the east and the south, and the sea on the west and the north. Geological exploration is divided into two stages: preliminary and detailed exploration. The methods of the exploration are drilling, in-situ testing, and soil lab test. The initial ground water level of the site is from 0.40 to 2.50 m, and the depth of the stable water level is from 0.30 to 2.30 m. The water is mainly the shallow phreatic water and the pore water in the middle and lower gravel soil. The original surface of the site is the weak silt soil layer with an average of 15 m thickness. This kind of soil layer is gray, fluid plastic, has high water content, large void ratio, and high compressibility, and it is easy to be disturbed and deformed. Moreover, it contains a small amount of shell debris and silt sand. The middle layer of the site mainly consists of clay with medium compressibility, the lower part is heavily weathered, and medium weathering quartz syenite porphyry, which has good mechanical properties. It is located between 60 and 80 m below the surface.

In order to conveniently facilitate the construction of the project, all the areas in the site need to be backfilled as the project was started 13 months earlier. The filling material is mainly composed of block stone, gravel, grail, and so on. The stone diameter is generally from 0.20 to 0.70 m, while the biggest one could be more than 1.50 m. The thickness of refilling treatment is from 6.50 to 15.80 m. The schematic diagram of the silt-rock fill layer is shown in Figure 1, and the diagram of refilling treatment in the site is shown in Figure 2.

Figure 2.

Backfill site.

The drilled shaft pile is used in this area. The diameter of the pile is 800 mm, and the effective length of the pile is 72 m. Pile concrete strength grade is C40. The whole length of the main reinforcement cage is 12ϕ25, which represents 12 reinforcement metals with 25 mm diameter. The thickness of longitudinal reinforcement is 55 mm. The full section of pile which enters the bearing layer is not less than 2400 mm, the sediment thickness is less than 50 mm, and the bearing layer of the pile is the medium weathering quartz syenite porphyry.

Experiments

All the piles were subjected to single pile vertical compressive static loading experiment according to the JGJ106-2003 technical code for testing of building foundation piles published by the Ministry of Construction of the People’s Republic of China. ¹⁷ Under axial compressive load with slow rate loading method, we add 600 kN load to each step until the pile’s axial deformation is greater than the threshold of 40 mm. The experiments have 15 min frequency and each loading has 1 h. The next stage of the loading is applied after the settlement is stable. The curves of the eight test piles are shown in Figure 3, and the eight test piles come from different sites. In order to conveniently compare the accuracy of the prediction model based on the measured curve, the eight test piles are conventionally loaded to 6200 kN. After that, they would not be unloaded, but continued to be loaded to s greater than 40 mm. The final axial deformation was found to be 43.8, 44.5, 42.6, 42.3, 47.9, 42.3, 41.2, and 42.5 mm for the eight piles, and the responding axial loading was 8600, 8600, 9200, 8600, 9200, 8600, 8000, and 8000 kN, respectively. The $Q-s$ curve of vertical static load test (SLT) of single pile with field measurement at the later stage is shown in Figure 3. For the pile data, three problems are concerned for the practicing engineers:

1. How to model the pile axial deformation path more accurately?

2. How to estimate the model parameters?

3. How to predict the ultimate bearing capacity according to different reliability levels, and thus can be used as the standard in the project site?

Figure 3.

Pile axial deformation trace in silt-rock layer area.

Model

Suppose that there are n different piles in the engineering experiment, where ${\mathit{y}}_i$ and ${\mathit{X}}_i$ ( $i=1,2,\dots ,n$ ) represent the axial deformation and axial loading, respectively. A mixed-effects model is proposed to describe the pile axial deformation ${\mathit{y}}_i$ over the axial loading ${\mathit{X}}_i$

$${\mathit{y}}_i={\mathit{X}}_i(\beta +{\mathit{b}}_i)+{\mathit{ϵ}}_i$$

(1)

here ${\mathit{y}}_i=(y_{i1},y_{i2},\dots ,y_{i{m}_i})\prime$ , where $m_i$ is the total number of measurements of the ith pile. ${\mathit{X}}_i=(x_{i1},x_{i2},\dots ,x_{i{m}_i})\prime$ is the design matrix, where ${\mathit{x}}_{\mathrm{ij}}=(x_{\mathrm{ij}},x_{\mathrm{ij}}^2,\dots ,x_{\mathrm{ij}}^p)\prime$ , $j=1,2,\dots \mathrm{,}{m}_i$ , and p is the degree of polynomial. ${\mathit{b}}_i=(b_{i1},b_{i2},\dots ,b_{\mathrm{ip}})\prime$ are the random effects that represent the intrinsic variation of single pile, where $b_{\mathrm{il}}$ follows a normal distribution with zero mean, ${\sigma }_l^2$ variance, and $\mathrm{cov}(b_{\mathrm{il}},b_{\mathrm{ih}})=0$ , where $l\ne h=1,2,\dots ,p$ . Thus, ${\mathit{b}}_i$ follows a p-dimensional normal distribution ${\mathbf{N}}_p(0,\mathbf{\sum })$ , where $\Sigma =\mathrm{diag}({\sigma }_1^2,{\sigma }_2^2,\dots {\sigma }_p^2)$ . $\mathit{\beta }=({\beta }_1,{\beta }_2,\dots ,{\beta }_p)\prime$ is the coefficient of the fixed effects, which usually represents the average deformation value of all tested piles. ${\mathit{\epsilon }}_i=({\mathit{\epsilon }}_{i1},{\mathit{\epsilon }}_{i2},\dots ,{\mathit{\epsilon }}_{i{m}_i})\prime$ are the measurement errors, and we assume that ${\mathit{\epsilon }}_{\mathrm{ij}}$ follows a normal distribution with zero mean, ${\sigma }_{\epsilon }^2$ variance, and $\mathrm{cov}({\epsilon }_{\mathrm{ij}},{\epsilon }_{\mathrm{ik}})=0$ , where $j\ne k=1,2,\dots ,m_i$ . Thus, $\mathit{\epsilon }~\mathit{N}(0,{\sigma }_{\epsilon }^2{\mathit{I}}_{m_i})$ , where ${\mathit{I}}_{m_i}$ is a $m_i\times m_i$ unit matrix.

With some simple calculations, the response variable ${\mathit{y}}_i$ conditioned on ${\mathit{b}}_i$ yields a $m_i$ -dimensional normal distribution ${\mathit{N}}_{m_i}({\mathit{X}}_i(\beta +{\mathit{b}}_{\mathit{i}}),{\sigma }_{\epsilon }^2{\mathit{I}}_{m_i})$ . In the model, the interesting parameters include the coefficient $\mathit{\beta }$ , the variance parameter ${\sigma }_{\epsilon }^2$ of the measurement errors, and the covariance matrix $\Sigma$ of the random effects. Denote $\mathbf{\Lambda }=(\mathit{\beta },\mathit{\Sigma },{\mathit{\sigma }}_{\epsilon }^2)$ . After integrating all the unobservable random effects $\mathit{B}=\{{\mathit{b}}_i,i=1,2,\dots ,n\}$ , the likelihood function of $\mathit{Y}=({\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_n)\prime$ is

$$L(\mathbf{\Lambda }|\mathit{y})=\underset{i=1}{\overset{n}{\Pi }}\{\underset{-\infty }{\overset{+\infty }{\int }}f(y_i|\beta ,{\sigma }_{\epsilon }^2,{\mathit{b}}_i)f\left({\mathit{b}}_i\right|\Sigma \left)\right\}d{\mathit{b}}_i$$

(2)

Bayesian approach

It is intractable to obtain the maximum likelihood estimations (MLEs) of the parameters $\mathbf{\Lambda }=(\mathit{\beta },\mathit{\Sigma },{\mathit{\sigma }}_{\epsilon }^2)$ through maximizing the above likelihood function, because the integration in terms of ${\mathit{b}}_{\mathit{i}}$ may be very complicated in high-dimensional space. Therefore, maximum likelihood method is not a good way for the mixed-effects model. An alternative way is Bayesian approach, which can avoid the complex integration using Markov chain Monte Carlo (MCMC) algorithm. In the reminder of this section, we explore the mixed-effects model using a hierarchical Bayesian framework. Prior distribution plays usually a key role in Bayesian method. The prior distribution should not only reflect the available prior information on the model parameters but also offer algebraic convenience.^36–39 When the prior distribution responding posterior distribution follows the same distribution family or similar kernel function (kernel function is a class of functions with same pattern or structure), the prior is called conjugate prior. Therefore, the conjugate prior can give a closed-form posterior and computational convenience. Before conducting Bayesian analysis for our proposed model (equation (1)), the priors of the polynomial coefficient $\beta$ , the normal variance parameter ${\sigma }_{\epsilon }^2$ , and the covariance matrix $\Sigma$ should be specified, respectively. Checking formula (1), the normal variance parameter ${\sigma }_{\epsilon }^2$ of the measurement errors has an inverse Gamma kernel function. Thus, we set the prior of ${\sigma }_{\epsilon }^2$ as inverse Gamma distribution. Second, considering the covariance matrix $\Sigma$ in the p-dimensional normal distribution, whose kernel function has same pattern with inverse Wishart distribution, the prior of $\Sigma$ is set as inverse Wishart distribution. Finally, considering the parameter $\beta$ , no distributions have the same structure. For this case, a noninformative prior $\pi (\mathit{\beta })\propto 1$ is widely used. ⁴⁰ Based on this analysis, set the prior of parameters Λ as follows

$$\mathit{\pi }({\sigma }_{\epsilon }^2)=\mathrm{IGa}({\tau }_{\epsilon },{\upsilon }_{\epsilon }),\mathit{\pi }(\mathit{\Sigma })=\mathrm{Iw}(\mathbf{\Psi },\upsilon ),\mathit{\pi }(\mathit{\beta })\propto 1$$

where ${\tau }_{\epsilon },{\upsilon }_{\epsilon },\Psi ,\mathrm{and}\upsilon$ are the predetermined constants based on the experience of engineers and experts. If this information is not available for the model parameter, small values of these hyper-parameters are suggested. ⁴¹

However, setting a noninformative prior for $\beta$ could also lead to the loss of available prior information on the model parameters. In our study, a special prior distribution for the parameter $\mathit{\beta }$ is presented as follows

$$\mathit{\pi }(\beta |{\mathit{\sigma }}_{\epsilon }^2)={\mathit{N}}_p(\mathit{\mu },g^2\mathit{V}{\sigma }_{\epsilon }^2)$$

where $\mathit{V}=\frac{1}{n}\sum _{i=1}^n{({\mathit{X}\prime }_i{\mathit{X}}_i)}^{-1}$ and $g^2$ is a hyper-parameter controlling the overall magnitude of the prior variance (e.g. $g^2\mathit{V}{\sigma }_{\epsilon }^2$ ). For the purpose of simplifying the posterior calculation, a noninformative prior $\pi (g)\propto 1$ is widely used.^42,43 However, this setting will lose the available data knowledge. A more attractive prior is inverse Gamma distribution, which is a conjugate prior for $g^2$ . Therefore, an inverse Gamma distribution is selected as the prior of the hyper-parameter $g^2$ in this article, that is

$$\mathit{\pi }\left(g^2\right)=\mathrm{IGa}({\tau }_g,{\upsilon }_g)$$

A more clear hierarchical Bayesian structure can be found in the blue doted square part in Figure 4, where the each level of hierarchical Bayesian structure is denoted through the left text in subfigure.

Figure 4.

Framework for hierarchical Bayesian model in this study.

According to the Bayesian theorem, given the prior and the likelihood function, we can conduct posterior inference through Bayesian formula. For the mixed-effects model, the posterior distribution is complicated. It is hard to directly perform statistical inference for its conditional posteriors. ⁴¹ However, MCMC algorithm is an alternative way to solve this problem efficiently. Gibbs sampling as a particular MCMC algorithm has been used frequently to conduct statistical inference in multidimensional problems. In this article, Gibbs sampling is also employed to estimate model parameters in the mixed-effects model. In order to calculate the axial ultimate bearing capacity of the pile and to obtain parameter estimation simultaneously, an improved Gibbs sampling algorithm is developed. This is based on the fact that the axial ultimate bearing capacity of the pile is a function of the model parameter. The detailed sampling procedure has been summarized in Appendix 1. The whole hierarchical Bayesian procedure is represented in the green dotted part in Figure 4.

The proposed Bayesian approach is based on the given polynomial degree p. Thus, p has an important impact on the mixed-effects model. Improper p will cause underfitting or overfitting. In fact, determination of p is a model-selection problem. The red dotted square part of Figure 4 describes this dynamic procedure. The common model-selection criteria include deviance information criterion (DIC), ⁴⁴ Bayesian information criterion (BIC), ⁴⁵ and Akaike’s information criterion (AIC), ⁴⁶ defined by $\mathrm{DIC}=-2\log (p(y|\theta ))+C$ , $\mathrm{BIC}=2\log (n)k+2\log (p(y|\hat{\theta }))$ , and $\mathrm{AIC}=2\log (p(y|\hat{\theta }))-2k$ , respectively. AIC and BIC are not appropriate for our model for two reasons. First, they require calculating the maximum likelihood estimation of the model parameters, which is not readily available from Gibbs sampling. Second, as Spiegelhalter et al. ⁴⁴ pointed out, BIC is not a suitable model-selection instrument for the hierarchical Bayesian model. Besides, it is important to note that DIC is appropriate for the motivation that we are interested in forecasting the future value of the ith tested product based on the early observations of this sample, which is consistent with our goals. According to this, DIC will be adopted for model-selection in this article. A direct usage of the above DIC is impropriate since C is an unknown constant, and an improved DIC is given by Spliegelhalter et al. ⁴⁴ as follows

$$\mathrm{DIC}=D(\hat{\theta })+2P_D$$

where $D(\hat{\theta })$ is two times the log-likelihood function; $\hat{\theta }$ is the estimated value of the model parameter; $p_D$ is the effective number of the model parameters and is defined as $p_D=\overline{D}(\theta )-D(\hat{\theta })$ , where $\overline{D}(\theta )={\mathit{E}}_{\theta }[D(\theta )]$ , which measures the extent of model fitting data. The larger it is, the worse the fit. Solving the exception $E_{\theta }[D(\theta )]$ directly is not easy, and a common usage approximate way is to calculate the mean of $D(\theta )$ over the posterior samples of $\theta$ , defined as ${\mathit{E}}_{\theta }[D(\theta )]\approx \frac{1}{N_w}\sum _{w=1}^{N_w}D({\theta }^{(w)})$ , where $N_w$ is the number of the posterior samples.

Besides, we also use the traditional $R^2$ statistics for model selection. To compute $R^2$ , we use the formula $R^2=1-\mathrm{RSS}/\mathrm{TSS}$ , where $\mathrm{RSS}=\sum _{i=1}^n\sum _{j=1}^{m_i}(y_{\mathrm{ij}}-{\hat{y}}_{\mathrm{ij}}{)}^2$ is the residual sum of square and $\mathrm{TSS}=\sum _{i=1}^n\sum _{j=1}^{m_i}(y_{\mathrm{ij}}-{\overline{y}}_{\mathrm{ij}}{)}^2$ is the total sum of square.

Application

In this section, we apply mixed-effects model to the piles data described in section “The data” via Bayesian approach. As we discussed in the previous section, p has an important effect on model inference since our Bayesian approach depends on p. Thus, p should be determined first based on DIC proposed in the previous section. In this article, we consider six candidate models under different polynomial degrees $p\in \{1,2,\dots ,6\}$ and omit the situation ( $p⩾7$ ), denoted as $M_p$ , because large p will increase model complexity and lead to overfitting. When we conduct Bayesian inference for the candidate model, the initial values of the hyper-parameters must be given first. With few prior knowledge in terms of the model parameters, we should set small values for the hyper-parameters in the prior. For model $M_p$ , we assign ${\tau }_{\epsilon }={\upsilon }_{\epsilon }={\tau }_{g^2}={\upsilon }_{g^2}=1$ , $\Psi ={\mathit{I}}_3$ , and $\mu =(-0.5,0.6,0.7)\prime$ . Gibbs sampling runs 150,000 iterations and the first 100,000 iterations are discarded.

The results of model-selection under different p can be found in Table 1. When the polynomial degree $p=3$ , the corresponding DIC and $\overline{D}$ are much smaller than the other models, which suggests that model $M_3$ is the best among the candidate models. Therefore, $M_3$ is the optimal candidate model for the piles data. $R^2$ statistics under different p are all close to 1 except for $M_1$ with $p=1$ , which suggests that the mixed-effects model with $2⩽p⩽6$ fit the pile data well. Figure 5 demonstrates the traces of the transformed DIC and $R^2$ as p increases from 1 to 6. The curve of DIC is low convex and the lowest point appears in $p=3$ . When p is greater than 2, the trace of $R^2$ becomes flat, which suggests that $M_2$ or $M_3$ is suitable for the piles data. There are two reasons for this selection. One reason is that greater p could lead to over-fit phenomena, which usually makes the model fit the experiment data effectively but cannot reflect the nature of the test set. Another reason is that greater p would make the model complex and increase the computational cost.

Table 1.

Model-selection index under different polynomial degree p.

Model	p	$\hat{D}$	n	DIC	$R^2$
$M_1$	1	718.33	726.60	735.00	0.8995
$M_2$	2	443.11	459.41	475.70	0.9794
$M_3$	3	266.42	290.22	313.92	0.9913
$M_4$	4	667.31	761.60	855.91	0.9936
$M_5$	5	1075.0	1168.0	1260.0	0.9948
$M_6$	6	1611.0	1651.0	1690.0	0.9953

Figure 5.

Trace of DIC and $R^2$ after taking the scaling and centering transformation.

According to the previous analysis, $M_3$ is the optimal model for this batch of piles data. The parameter estimation of the model and the responding standard deviation (SD), mean, median, and credible interval are presented in Table 2. Figure 6 displays the eight test piles’ axial deformation trace over the axial loading, where the real and the dashed lines represent the realistic and the fitting trace, respectively. Obviously, all the fitting lines are extremely close to the real lines.

Table 2.

Parameter estimation of the model based on the complete test data.

Parameter	Mean	SD	Median	95% CI
${\beta }_1$	1.456	0.427	1.411	(0.734, 2.077)
${\beta }_2$	−0.140	0.415	−0.205	(−0.731, 0.475)
${\beta }_3$	−0.418	0.521	−0.252	(−1.137, 0.124)
${\sigma }_1^2$	1.225	1.053	0.930	(0.083, 4.008)
${\sigma }_2^2$	1.292	1.110	0.989	(0.084, 4.176)
${\sigma }_3^2$	1.327	1.095	1.047	(0.092, 4.132)
${\sigma }_{\epsilon }^2$	1.289	0.188	1.280	(0.946, 1.686)

Figure 6.

Realistic and predict path under the complete test data.

The prediction of the axial bearing capacity makes more sense than fitting the pile trace since the complete engineer testing is usually expensive. In order to predict the axial bearing, the piles data are divided into training set and test set at axial loading 6200 kN. Solving the mixed-effects model under the training set, the result has been presented in Table 3. After this, a common practice is conducting point prediction at special test point via the training model. However, the practical sense of the point prediction does not make sense. In this article, the quantile approach is used to acquire the coverage rate. First step is extracting a sequence of the predicted values of the response variable ${\hat{y}}_{\mathrm{ij}}^{(w)}$ from the parameter sequence $L^{(w)}$ through Gibbs sampling, where $i\in \{1,2,3,\dots ,8\}$ and $j\in \{6.8,7.4,8.0,8.6,9.2\}$ . And then solving the quantiles using the quantile algorithm 7 provided by Hyndman and Fan, ⁴⁷ which can be written as

$$Q(p)=(1-\gamma )X_{(k)}+\gamma X_{(k+1)}$$

(3)

where $X_{(k)}$ is the kth-order statistic, $p\in [(k-1)/(n-1),k/(n-1)]$ , and n is the sample size. When n is an even number, $\gamma =1-p-\lfloor 1-p\rfloor$ ; when n is an odd number, $\gamma =3/2-p-\lfloor 3/2-p\rfloor$ , where $\lfloor x\rfloor$ denotes the largest integer not greater than x. Then, the 95% and 90% one-sided credible intervals are demonstrated in Figure 7. Most of the test points fall into the 95% one-sided credible interval except four points, which are slightly better than 90% one-sided confidence interval. Although some points are outside the 95% area, they are all measured at the relative lower loadings and below the minimum lower bound; the corresponding pile axial deformation is further smaller than the axial deformation threshold of 40 mm, which means that these points are safety measuring sites. Therefore, the quantile approach is appropriate for predicting the axial ultimate bearing capacity of piles.

Table 3.

Model parameter estimation based on the complete test data truncated at 6200 kN.

Parameter	Mean	SD	Median	95% CI
${\beta }_1$	−0.336	0.311	1.411	(−0.682, 0.369)
${\beta }_2$	0.507	0.363	0.446	(0.014, 1.055)
${\beta }_3$	−0.010	0.271	−0.121	(−0.291, 0.380)
${\sigma }_1^2$	0.236	0.432	0.294	(0.043, 0.540)
${\sigma }_2^2$	1.679	0.395	1.762	(0.811, 2.200)
${\sigma }_3^2$	0.832	0.470	0.887	(0.112, 1.603)
${\sigma }_{\epsilon }^2$	8.723	1.648	8.610	(5.896, 12.260)

Figure 7.

The one-sided credible interval for the test data.

We further present the one-sided credible intervals of the axial ultimate bearing capacity of piles under different combinations of guarantee rate and pile axial deformation values, and the results are listed in Table 4. According to JGJ94-2008 of China, ¹⁷ the axial ultimate deformation is 40 mm, and the design axial loading value is 6200 kN in this pile experiment. The safety of the building is guaranteed, while this would require much more construction budget. Actually, this design value is too conservative according to Table 4. From Table 4, we can see that when the the axial deformation reaches 40 mm, the axial ultimate loading under different guarantee rates from 80% to 100% with 2.5% increment is 9308, 9261, 9200, 9102, 8986, 8902, 8827, 8746, and 8362 kN, respectively. If we choose the axial ultimate loading as 8362 kN, then there is a 100% guarantee rate that the vertical displacement of the pile foundation will not exceed 40 mm. Such a value is 34.8% higher than the original design value of 6200 kN. Thus, a single pile can not only save cost up to US$15,000 but also shorten the construction time of the pile foundation.

Table 4.

Axial ultimate bearing capacity of piles under different guarantee rates and displacement.

Displacement (mm)	100%	97.5%	95%	92.5%	90%	87.5%	85%	82.5%	80%
40	8362	8746	8827	8902	8986	9102	9200	9261	9308
35	7820	8149	8219	8284	8358	8492	8603	8658	8699
30	7237	7508	7565	7619	7681	7838	7964	8011	8045
25	6591	6809	6855	6896	6964	7119	7260	7298	7324
20	5871	6035	6069	6099	6134	6273	6428	6468	6493
15	5021	5151	5176	5197	5222	5286	5364	5397	5426
10	3866	4062	4084	4099	4113	4125	4137	4151	4166
5	2222	2449	2496	2536	2582	2638	2692	2731	2766

Besides, we also compute the average values of ultimate bearing capacity based on empirical formula method ( $M_e$ ), ⁴⁸ polynomial prediction model ( $M_p$ ), ⁴⁹ our proposed mixed-effects model ( $M_3$ ), and the static load test (SLT). ¹⁷ The results are listed in Table 5. We can see from Table 5 that

• The axial ultimate bearing capacity of $M_3$ could reach 97.5% of the theoretical ultimate bearing capacity of piles which is 8570 kN (SLT), denoted as ${\tilde{Q}}_u$ , which is 1.31 times higher than the original design of 6200 kN.

• The performance of model $M_e$ is highly dependent on the sample size. In this dataset, the sample size is only eight, which could make the model $M_e$ fail, and the axial ultimate bearing capacity only reaches 74% of ${\tilde{Q}}_u$ .

• Model $M_p$ ignores the random effects and uses fixed errors instead of random errors. The prediction of the axial ultimate bearing capacity only reaches 89% of ${\tilde{Q}}_u$ .

• Thus, our proposed model performs the best due to considering the fixed effect, random effect, and random error simultaneously.

Table 5.

Ultimate bearing capacity values obtained by different models.

Methods	$M_e$	$M_p$	$M_3$	SLT
$Q_u(\mathrm{kN})$	6383	7629	8362	8570

Conclusion

In this article, as we mentioned in section “The data,” there are three problems that need to be solved based on the experimental data:

• For the first problem, we propose a mixed-effects model, which can describe the mean ultimate bearing capacity of the project site and the unit-to-unit variation among the piles. Besides, the model is very flexible due to the unknown degree of polynomial p, and p can be determined by the DIC. As we can see in the data analysis, the new model fits the data well when p = 3.

• For the second one, a hierarchical Bayesian method is presented to estimate the model parameters. The g prior and inverse-Gamma prior are assigned for the model parameters, which not only facilitates the posterior calculation but also reflects more available prior information.

• For the last one, we predict the ultimate bearing capacity based on the posterior samples, and use the one-sided interval estimates, which can be achieved by the quantile approach. Under different guarantee rates and the axial deformation, we compute the lower quantiles of the ultimate bearing capacity. As we can see in Table 4, the values are much larger than the design value of 6200 kN in practice, which means a much larger design value can be chosen, thus saving the cost.

Compared to the traditional vertical bearing capacity formula of single pile, the mixed-effects model takes many factors into consideration, which can predict the vertical ultimate bearing capacity of single pile much better. Due to utilizing hierarchical Bayesian method, the point estimation and interval evaluation of model parameters can be conducted simultaneously, and the improved Gibbs sampling algorithm can make the estimation more efficient.