# Prediction of California Bearing Ratio of Granular Soil by Multivariate Regression and Gene Expression Programming

Authors: Mudthir Bakri; Ibrahim Aldhari; Mohammed Saleh Alfawzan
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

This research demonstrates the results of an investigation into the California bearing ratio (CBR) of granular soils from Qassim region, Saudi Arabia, using multilinear regression (MLR), pure quadratic (PQ) models, and gene expression programming (GEP) methods utilized to develop mathematical models for estimating the CBR based on basic soil index properties. In this study, samples were collected from different borrowing pits in the Qassim area. Forty-three samples of soil were taken and transferred to a laboratory for examination. Seven multilinear regressions and seven PQ models were investigated, while four GEP models were made. The selection of each model variable depends on soil indices, grouping into grain size distribution, Atterberg limits, and compaction parameters. The results of this analysis showed that the PQ model had a higher accuracy [coefficient of determination (R2) = 0.89, root mean square error (RMSE) = 16.006, uncertainty (U95) = 16.17, and reliability = 57%] than the multilinear regression model, which has a lower accuracy model [R2 = 0.811, RMSE = 20.791, U95 = 15.569, and reliability = 51%]. The best GEP model yields [R2 = 0.776, RMSE = 22.552, U95 = 15.787, and reliability = 53%]. Furthermore, sensitivity analysis was conducted to distinguish the influences of different input variables on CBR; it was found that fines percentage (F200), maximum dry density (MDD), and optimum moisture content (OMC) are the most influential variables.


---

## Body

## 1. Introduction

Saudi Arabian road network consists of 73,171 km. Saudi Arabia has a vast and well-developed transportation network. Roads carry approximately 65 percent of freight and 80 percent of passenger transportation [11]. The total number of automobiles in Saudi Arabia exceeds 12 million, and the growth rate of the number of vehicles is approximately 5% per annum over the last five years [20]. Geotechnical engineering is one of the greatest crucial subjects in the initial stages of infrastructure planning and design, owing to the concept that poor geotechnics might consequence in unnecessarily high costs if not addressed properly. The CBR value is a critical parameter in the structural design of pavements. The CBR value may be determined directly in the laboratory using AASHTO T 193 and ASTM D 1883.Two broad correlations for predicting CBR were also presented in the Guide for Mechanistic and Empirical design (2001). The models are constructed using empirical parameter data for the following soils; D60 (passing 60% diameter), P200 (percent passing sieve no. 200; US sieve), and PI (plastic index). PI and % passing no. 200 US sieve (or 0.075 mm size sieve) are two parameters that are included in the suggested model for plastic fine-grained soils [15]. The recommended equation is(1)CBR=751+0.728F200PI,where F200 is passing no. 200 US sieve (%).PI is the plasticity index. For the nonplastic coarse-grained soil, the proposed equation is(2)CBR=28D600.358.Determine the CBR value of each soil sample that demands at least four days in the laboratory. Historically, civil engineers have struggled to determine a representative CBR value for pavement design [36]. The CBR value varies according to the kind of soil and its various qualities. Due to the variability in soil engineering features, conducting this experiment with soil samples gathered from a few places cannot be used for the entire road path. To overcome this, a substantial number of specimens must be collected for testing, which causes the method to be costly and time consuming. Black made the first attempt to predict CBR in 1962 [13]; his work depends on remolded soil and its friction angle. It is also shown that the suction of such a remolded soil and its true angle of friction can be inferred from data on the liquid limit, plastic limits, and water content of the soil. Agarwal and Ghanekar predicted the CBR deepening on the grain size distribution as well as the consistency parameters, specifically the PI [2]. These attempts were followed by many subsequent attempts using more advanced methods of prediction, each of these attempts expresses the classification, considerable types of index properties, and unique study area samples [3, 5, 7]. CBR was determined by Yildirim and Gunaydin [35] utilizing soft computing; they used soil index properties and compaction properties of subgrade soil in Turkey as input parameters. Their work revealed a high correlation (R2 = 0.80–0.95) employing soil indices. They recommend using the provided formulas in the early stages of the design process [35]. Kumar (2014) proposed a linear regression model for ML and MI soils on the basis of the variables PI, liquid limit (LL), plastic limit (PL), MDD, and OMC, for ML and CBR and PI have an inverse linear regression relation [32]. Several linear regression models for estimating soaked CBR in clay soils are compared by Ramasubbarao and Siva [26, 26]. Tenpe et al. constructed several mathematical models to predict CBR utilizing a variety of soil parameters as independent variables, with the assistance of the GEP and ANN models. They also conducted comparative research to demonstrate that there is little difference between the findings acquired from GEP and those from ANN [6]. Numerous scholars have presented such models [6, 29, 34], but there is still a knowledge gap concerning the assessment and comparison of their predictive capability. Aleksandra et al. investigated the effects of compaction parameters, namely MDD and OMC on CBR [10].Interestingly, no correlation formula has been created utilizing the soil of the Qassim area. The region’s geological formations are divided into two large landmasses: the Arabian Shield in the west and the Arabian Shelf in the east. The Arabian Shield is made up of igneous and metamorphic rocks, while the Arabian Shelf is made up of sedimentary rocks. The weathering products of these rocks have generated the soils. Predicting the value of the CBR test would shorten the time from 96 hours to minutes, achieving such a goal would reduce the overall cost of projects in various ways. Al-Qassim region is one of the most important logistical areas. Since at present, there is no published work correlating soil properties and CBR in the Qassim area and its environs. So that the purpose of the study is to determine if GEP and multivariate regression approaches are applicable to estimate CBR of flexible pavement subbases and subbase layers in the Qassim area.


### 1.1. Samples Collection

Qassim is a region of predominantly desert nature with a multigeological nature, and the so-called Arab Shield is located in the southwestern part of it. It is located in the middle of the Kingdom of Saudi Arabia [4]. Qassim is considered a vital region linking the north with the south and the west with the east, thus heavy carriers pass through its roads. In Saudi Arabia, roads are built in accordance with the Ministry of Transportation [8] (MOT, 1998) regulations. Tables 1 and 2 present the gradation of the subbase layer and the limit of indices. During the building of new roads in Qassim area, samples were taken from various places and borrowed pits. These samples were sieved to find the percentage of aggregate (gravel) and the percentage of aggregate passing through sieve no. 4 (4.75 mm) and were retained on sieve no. 200 (0.075 mm) (sand) and the percentage passing through sieve no. 200 (0.075 mm) (fines). The consistency tests were conducted on fine soil to determine the liquid and plastic limits. The relationship between moisture content and dry density was then determined by the proctor modified test. Then, each sample was made at OMC and at different densities, with the one with the highest dry unit weight being tested to see what its CBR was for each unit weight.Table 1
Ministry of transportation specification limits or gradation of the subbase layer.
Sieve size designationMOT specification limits (% passing)Upper limitLower limit2″1001001 1/2″ 1″100901″85553/4″80503/8″7040No. 46030No. 105020No. 403010No. 200150Table 2
MOT specifications of soil indices for materials of subbase layer.
Sand equivalentPlasticity indexAbrasion lossCBR25 min.6 max.50 max.50 min.


### 1.2. Experimental Work

In this research, to study the factors which explain CBR characteristics and its prediction, soil samples were gathered from 43 borrowed pits in the Qassim area as shown in Figure1, with a different number of specimens taken from each location. The test method used to determine the CBR is given by MOT (1998) speciﬁcations, in which it complies with MRDTM: 213 ASTM D 1883-AASHTO 193 [31]. The compaction parameters are found in accordance to ASTM D 698–1557, AASHTO T99-T180. The particle size distribution shall follow MRDTM: 204. MRDTM: 208/209 was used to conduct the test of LL, PL, and the PI. These limits were determined on the soil samples having a fine percentage greater than 5%. Table 3shows that in general the average numbers of fines (silt and clay), sand, and gravel were 15.32%, 28.3%, and 56.421%, respectively. The averages of LL, PL, and PI were 7.8%, 6.18%, and 1.62%, respectively, while the average water content OMC was 7.73%, and the average MDD was 2.15 g/cm3.Figure 1
Location of samples.Table 3
Descriptive statistics analysis of the tested samples.
Gravel %Fines %Sand %Liquid limitPlastic limitPlasticity indexMaximum dry densityOptimum moisture contentCBR %Mean56.42115.32728.37.86.1811.6232.157.7389.45Std. error of mean3.0281.0542.351.8411.480.3770.01240.2887.37Std. deviation19.8566.9115.4312.079.72.4730.081341.8948.32Variance394.2747.81238.17145.894.066.120.073.5912334.77Median63.41419.10002.157760Minimum10.705.312.630.000.000.001.915.6021Maximum80.3728.375.835.2028.706.802.2416.50172AASHTO classification system [(ASTM designation D-3282; AASHTO method M145)] was used to classify the soil samples, and it was found that 34 samples are classified as A-1a and A-1-b, and ten samples were classified as A-2.Grain size analyses were performed firstly to classify the soil in the studied area and to quantify the number of fines [F200], sand percent [S], and gravel percentage [G]; grain size distribution curve is shown in Figure 2. Then after analysis of all samples, the tests of LLs and plastic limits were conducted; the plasticity chart is shown in Figure 3. Thirty percent of soil samples had a liquid limit larger than zero, which was classed as inorganic clays of low plasticity to medium plasticity and inorganic silt of low compressibility. The majority of the samples depicted fall above line A and below line U. The existence of soil samples parallel to the line A on the plasticity chart shows that they have the same geological origin [19, 30]. Table 3 provides a full statistical summary of the soil. To see the distribution of the parameters more clearly, frequency histograms were created for each parameter in Figure 4.Figure 2
Representative grain size distribution of soil samples.Figure 3
Plasticity diagram showing sample results.Figure 4
GEP algorithm.


### 1.3. Multivariate Regression

Single and multivariate regression analysis are well-known predictive techniques that are frequently used to construct models for predicting desired output parameters based on specified input factors such as the physical and mechanical characteristics of rocks [17]. Single regression studies examine the connections between input and output variables using a linear, logarithmic power or exponential function. Multivariate regression studies are advantageous when several input parameters are present in complicated connections, and the common form of a multivariate regression equation is(3)CBR=a+∑i=1nbGpi+ε,where CBR is California bearing ratio (dependent variable). Gpi is soil index or mechanical (independent variable), n is the number of variables, bi is the coefficient of regression, a is the constant, and ε is the error estimate.CBR will be used as the dependent variable, while the soil parameters in Table4 will be used as independent variables individually, or they will be used as groups of parameters depending on the relationship between these data and the correlation matrix as shown in Table 5. In addition, the selection of input variables of multilinear equations depends on grouping the soil indices into three categories; the first one is the grain size distribution, in which the parameters are gravel percent, sand percent, and fines percent. The second category is Atterberg limits, which are LL, PL, and PI. The third category is “compaction variables,” which include “MDD” and optimum water content. This will help to find the highest correlation parameters with CBR.Table 4
Summary of sample results.
No.Gravel %Fines %Sand %LLPLPIMDD g/cm3OMC %CBR %1727.620.40002.1856.412324820.631.423.5194.52.128.4473685.326.70002.158.61144737.919.10002.196.8113549.921.228.920.815.15.72.1578.2123672.19.218.70002.2366.9172763.420.116.520.6164.62.1058.480843.213.443.40002.146.754975.397.517.110002.2356.61501015.526.957.60002.0926.5271140.323.236.50002.1357.343.81226.919.753.40002.156745.51380.37712.630002.2376.81551462.719.817.534.528.75.82.1447.3601510.713.575.80002.0427.7351676.119.414.490002.2366.91701775.397.517.110002.2356.61501862.82116.20002.1286.8511975.657.217.150002.236.61502059.921.518.60002.1478.2452165.717.31724.518.65.92.1668.3512243.213.443.40002.1846.6542340.715.443.90002.2056.7602473.119.617.290002.2257.11202513.624.861.617.915.32.62.1497.2212672.698.418.9110.98.72.22.2316.91542749.921.228.920.815.15.72.1578.2512873.38.817.90002.2375.61432976.649.116.260002.2366.51533061.619.518.92217.24.82.1158.8563171.11414.932.127.54.62.1567.9483263.89.127.10002.2366.91703345.128.326.634.227.86.42.1259.3493463.420.116.520.6164.62.1058.4803535.225.938.9235.228.46.81.9910.5373625.924.249.90001.91416.5483745.319.435.30002.0510.5543876.416.916.690002.246.21303922.924.952.20001.9311.8584037.521.5411812.45.62.059.2334170.9510.818.250002.2361264270.58.5210002.216.61234376.38.515.20002.2066119Table 5
Correlations matrix.
CBRGrF200SLLPLPIMDDOMCCBR10.773−0.82−0.623−0.42−0.42−0.4030.716−0.437Gr0.7731−0.738−0.953−0.121−0.118−0.1270.727−0.46F200−0.82−0.73810.4990.5240.5160.536−0.7410.565S−0.623−0.9530.4991−0.081−0.081−0.078−0.60.337LL−0.42−0.1210.524−0.08110.9980.968−0.2820.26PL−0.42−0.1180.516−0.0810.99810.951−0.2750.252PI−0.403−0.1270.536−0.0780.9680.9511−0.2960.28MDD0.7160.727−0.741−0.6−0.282−0.275−0.2961−0.844OMC−0.437−0.460.5650.3370.260.2520.28−0.8441


### 1.4. Multinonlinear Regression

Nonlinear regression estimates the dynamic connection between the dependent variable and the function. Nonlinear multivariable regression estimates the interaction between independent and dependent variables. There are various modeling techniques and regression approaches for specific aspects, especially when data are limited [17].It is possible that the greatest association between the answer and independent variables may be explained by step-by-step regression. This algorithmic technique selects the proper model subsets to filter forward or backward. First, choose a permanent model and use the model words before optimizing fitness. That is, step-by-step regression is a forward selection procedure that examines the significance of all prior variables. For variables whose partial square sums do not meet the minimal requirement, the retroactive exclusion procedure is improved and variables are excluded one at a time until the minimum condition is met. Step-by-step regression requires more computations than forward or backward detailed calculations, but it yields better results [33]. Scheffe proposed a second-order polynomial for MNR analysis of mechanical characteristics [21]. MNLR in PQ forms has been studied.


### 1.5. Gene Expression Programming (GEP)

GEP, which is an algorithm based on genetic algorithm (GA) and genetic programming (GP), was ﬁrst proposed by Ferreira [27]. The GEP technique is an integration of two evolutionary strategies that are broadly acknowledged in the scientific community. The first approach is called the GA methodology, and it simplifies complex relationships by depicting them with linear structures called chromosomes that are of a fixed length. In addition, the second method is known as GP, and it makes use of the expression tree (ET) configurations that may take on a wide range of forms and dimensions [22]. Developing a program code recorded in linear chromosomes of fixed length, it exhibits some characteristics akin to biological evolution. GEP’s primary goal is to create a mathematical function that fits a collection of data supplied to the GEP model. The GEP method performs symbolic regression on the mathematical equation using the majority of GA’s genetic operators. GEP and GA have some distinctions. Individuals are shown in GA as symbolic strings of fixed length (chromosomes), but they are depicted in GP as nonlinear entities of varying sizes and forms (parse trees), whereas in GEP, each individual is represented by a fixed-length string that is used to make ETs of different sizes and shapes [9]. The primary GEP process is presented in Figure 5. Five variables are employed in GEP analysis: function set, terminal set, fitness function, control parameters, and stop condition. After encoding the issue into the candidate solution and specifying the fitness function, a randomly generated population of viable people (chromosomes) is constructed and then translated into an ET matching to a mathematical equation. Following that, the projected values are compared to the actual ones, and the fitness score for each chromosome is calculated. The algorithm is ended when a suitable fitness score is reached. If this is not the case, some chromosomes are randomly selected and then mutated to generate new generations. This method is repeated until a suitable fitness score is produced, at which point noncoding chromosomal analysis is undertaken to determine the optimal solution to the problem [14].Figure 5
Frequency histograms of the index, MDD, and OMC.


### 1.6. GEP Elements

The ET and chromosomes are the two most important components of the GEP. The chromosome is made up of one or more genes and it is represented by a mathematical equation on the genome. To code any mathematical equation, two bilingual and conclusive languages named Karva languages (the language of genes and the language of ET) are utilized, both of which are bilingual and conclusive. This property of the GEP is important for inferring the genotype with more precision than it would otherwise be possible. The GEP genes are divided into two sections, which are referred to as the head and tail. It is necessary to use the head of a gene, which contains certain mathematical functions, variables, and constants, for the encoding of a particular function. However, the tail of a gene that contains just variables and constants is used as supplemental terminal symbols in addition to the main terminal symbols. These symbols are used in the event that the terminal symbols in the head are insufficient to express a particular function. Many other functions, such as fundamental arithmetic operators, trigonometric functions, and any other mathematical or user-defined functions, may be utilized in the head of a gene.


### 1.7. GEP Operators

The major operators in GEP are selection, mutation, transposition, and crossover (recombination), which are all based on natural selection. By using these operators, the chromosomes may be adjusted to improve the fitness score of the following generation for survival. The operator rates, which are specified at the outset of the model construction, indicate a particular chance of a chromosome appearing in the model. In general, a mutation rate of 0.001 to 0.1 is suggested, with the lower limit being 0.001. The transposition operator and the crossover operator, on the other hand, are recommended to be 0.1 and 0.4, respectively [24].


### 1.8. Model Performance

To examine the performance of the prediction models, the RMSE, the coefficient of determination (R2), as well as mean absolute percentage error (MAPE), the scatter index (SI) [23], Bias, and index of agreement (IOA) are being used as shown in (4)–(9).(4)RMSE=∑i=1NCBRCBRm−CBRp2N,R2=1−∑i=1NCBRm−CBRp2∑i=1NCBRp2,MAPE=1N∑i=1NCBRm−CBRpCBRp×100,SI=1/N∑i=1NCBRpi−CBRp¯+CBRmi−CBRm¯21/N∑i=1NCBRmi,Bias=∑i=1NCBRp−CBRmN,IOA=1−∑i=1NCBRpi−CBRmi2∑i=1NCBRpi−CBRm¯+CBRmi−CBRm¯2.n, CBRm, and CBRp indicate the total number of data measured experimentally, and the developed predicted value, respectively.In addition, the Durbin–Watson statistic was used to determine the existence of multicollinearity [1]. The Durbin–Watson statistic has a range of values between 0 and 4, and the midpoint (that is, 2) indicates that there is no association between the input variables. Therefore, a value between 1.5 and 2.5 is sufficient for obtaining models that are unaffected by multicollinearity. This section covers mathematical expressions and models. Each model’s R2 and P-value (p) are described in the Tables 6 and 7.(5)DW=∑t=2Tet−et−1∑t=1Tet2,whereas DW is the Durbin–Watson test and et and et−1 least-squares regression residuals.


### 1.9. Uncertainty and Reliability

The main purpose of uncertainty analysis is to narrow the predicted range in which the real value of an experiment’s result resides. This estimated range is expressed as an interval and is referred to as the uncertainty interval. It may be approximated using the calculated errors for the measuring method of the experiment in question [28]. U95 is one of the uncertainty analysis algorithms used to calculate the uncertainty interval. The value of U95 relating to a specific experiment result may be read as follows: if the provided experiment is performed again, the real value of the experiment’s outcome will fall within the specified uncertainty range around 95 times out of every 100 trials. Furthermore, the value of U95 is provided by(6)U95=1.96N∑i=1NCBRmi−CBRm¯2+∑i=1NCBRmi−CBRpi2.


### 1.10. Relative Absolute Error and Reliability

The relative absolute error is a relative statistic that compares a predictive model’s performance to that of a simple model. The predictive model’s performance is defined as the total absolute difference between the actual and anticipated values (i.e., the error). The basic model’s performance is defined as the total absolute difference between the realized value and the average of all realized values [25]. In other words, the relative absolute error determines if a model outperforms just forecasting the average (i.e., the simple model). The relative absolute error formula is(7)RAE=∑i=1NCBRmi−CBRpi∑i=1NCBRmi−CBRm¯.The meaning of the relative absolute error is straightforward: if the RAE is less than one, the model outperforms the basic model. The relative absolute error for a perfect model is 0. RAE should ideally be as near to zero as feasible.A statistical technique known as reliability analysis may be used to evaluate a model’s level of general consistency. To be more specific, it determines whether or not a proposed model reaches a level of performance that is within acceptable limits. A measurement that may be used for the examination of reliability is defined as(8)Reliability=100%N∑i=1Nki.Then after we check to see if RAEi is greater thanδ, if it is, we set ki equal to 1; otherwise, we set it equal to 0. We define ki as the number of times the value of RAE is either less than or equal to that of δ. According to Chinese standards, the ideal value of δ is 0.2 [28].


## 1.1. Samples Collection

Qassim is a region of predominantly desert nature with a multigeological nature, and the so-called Arab Shield is located in the southwestern part of it. It is located in the middle of the Kingdom of Saudi Arabia [4]. Qassim is considered a vital region linking the north with the south and the west with the east, thus heavy carriers pass through its roads. In Saudi Arabia, roads are built in accordance with the Ministry of Transportation [8] (MOT, 1998) regulations. Tables 1 and 2 present the gradation of the subbase layer and the limit of indices. During the building of new roads in Qassim area, samples were taken from various places and borrowed pits. These samples were sieved to find the percentage of aggregate (gravel) and the percentage of aggregate passing through sieve no. 4 (4.75 mm) and were retained on sieve no. 200 (0.075 mm) (sand) and the percentage passing through sieve no. 200 (0.075 mm) (fines). The consistency tests were conducted on fine soil to determine the liquid and plastic limits. The relationship between moisture content and dry density was then determined by the proctor modified test. Then, each sample was made at OMC and at different densities, with the one with the highest dry unit weight being tested to see what its CBR was for each unit weight.Table 1
Ministry of transportation specification limits or gradation of the subbase layer.
Sieve size designationMOT specification limits (% passing)Upper limitLower limit2″1001001 1/2″ 1″100901″85553/4″80503/8″7040No. 46030No. 105020No. 403010No. 200150Table 2
MOT specifications of soil indices for materials of subbase layer.
Sand equivalentPlasticity indexAbrasion lossCBR25 min.6 max.50 max.50 min.


## 1.2. Experimental Work

In this research, to study the factors which explain CBR characteristics and its prediction, soil samples were gathered from 43 borrowed pits in the Qassim area as shown in Figure1, with a different number of specimens taken from each location. The test method used to determine the CBR is given by MOT (1998) speciﬁcations, in which it complies with MRDTM: 213 ASTM D 1883-AASHTO 193 [31]. The compaction parameters are found in accordance to ASTM D 698–1557, AASHTO T99-T180. The particle size distribution shall follow MRDTM: 204. MRDTM: 208/209 was used to conduct the test of LL, PL, and the PI. These limits were determined on the soil samples having a fine percentage greater than 5%. Table 3shows that in general the average numbers of fines (silt and clay), sand, and gravel were 15.32%, 28.3%, and 56.421%, respectively. The averages of LL, PL, and PI were 7.8%, 6.18%, and 1.62%, respectively, while the average water content OMC was 7.73%, and the average MDD was 2.15 g/cm3.Figure 1
Location of samples.Table 3
Descriptive statistics analysis of the tested samples.
Gravel %Fines %Sand %Liquid limitPlastic limitPlasticity indexMaximum dry densityOptimum moisture contentCBR %Mean56.42115.32728.37.86.1811.6232.157.7389.45Std. error of mean3.0281.0542.351.8411.480.3770.01240.2887.37Std. deviation19.8566.9115.4312.079.72.4730.081341.8948.32Variance394.2747.81238.17145.894.066.120.073.5912334.77Median63.41419.10002.157760Minimum10.705.312.630.000.000.001.915.6021Maximum80.3728.375.835.2028.706.802.2416.50172AASHTO classification system [(ASTM designation D-3282; AASHTO method M145)] was used to classify the soil samples, and it was found that 34 samples are classified as A-1a and A-1-b, and ten samples were classified as A-2.Grain size analyses were performed firstly to classify the soil in the studied area and to quantify the number of fines [F200], sand percent [S], and gravel percentage [G]; grain size distribution curve is shown in Figure 2. Then after analysis of all samples, the tests of LLs and plastic limits were conducted; the plasticity chart is shown in Figure 3. Thirty percent of soil samples had a liquid limit larger than zero, which was classed as inorganic clays of low plasticity to medium plasticity and inorganic silt of low compressibility. The majority of the samples depicted fall above line A and below line U. The existence of soil samples parallel to the line A on the plasticity chart shows that they have the same geological origin [19, 30]. Table 3 provides a full statistical summary of the soil. To see the distribution of the parameters more clearly, frequency histograms were created for each parameter in Figure 4.Figure 2
Representative grain size distribution of soil samples.Figure 3
Plasticity diagram showing sample results.Figure 4
GEP algorithm.


## 1.3. Multivariate Regression

Single and multivariate regression analysis are well-known predictive techniques that are frequently used to construct models for predicting desired output parameters based on specified input factors such as the physical and mechanical characteristics of rocks [17]. Single regression studies examine the connections between input and output variables using a linear, logarithmic power or exponential function. Multivariate regression studies are advantageous when several input parameters are present in complicated connections, and the common form of a multivariate regression equation is(3)CBR=a+∑i=1nbGpi+ε,where CBR is California bearing ratio (dependent variable). Gpi is soil index or mechanical (independent variable), n is the number of variables, bi is the coefficient of regression, a is the constant, and ε is the error estimate.CBR will be used as the dependent variable, while the soil parameters in Table4 will be used as independent variables individually, or they will be used as groups of parameters depending on the relationship between these data and the correlation matrix as shown in Table 5. In addition, the selection of input variables of multilinear equations depends on grouping the soil indices into three categories; the first one is the grain size distribution, in which the parameters are gravel percent, sand percent, and fines percent. The second category is Atterberg limits, which are LL, PL, and PI. The third category is “compaction variables,” which include “MDD” and optimum water content. This will help to find the highest correlation parameters with CBR.Table 4
Summary of sample results.
No.Gravel %Fines %Sand %LLPLPIMDD g/cm3OMC %CBR %1727.620.40002.1856.412324820.631.423.5194.52.128.4473685.326.70002.158.61144737.919.10002.196.8113549.921.228.920.815.15.72.1578.2123672.19.218.70002.2366.9172763.420.116.520.6164.62.1058.480843.213.443.40002.146.754975.397.517.110002.2356.61501015.526.957.60002.0926.5271140.323.236.50002.1357.343.81226.919.753.40002.156745.51380.37712.630002.2376.81551462.719.817.534.528.75.82.1447.3601510.713.575.80002.0427.7351676.119.414.490002.2366.91701775.397.517.110002.2356.61501862.82116.20002.1286.8511975.657.217.150002.236.61502059.921.518.60002.1478.2452165.717.31724.518.65.92.1668.3512243.213.443.40002.1846.6542340.715.443.90002.2056.7602473.119.617.290002.2257.11202513.624.861.617.915.32.62.1497.2212672.698.418.9110.98.72.22.2316.91542749.921.228.920.815.15.72.1578.2512873.38.817.90002.2375.61432976.649.116.260002.2366.51533061.619.518.92217.24.82.1158.8563171.11414.932.127.54.62.1567.9483263.89.127.10002.2366.91703345.128.326.634.227.86.42.1259.3493463.420.116.520.6164.62.1058.4803535.225.938.9235.228.46.81.9910.5373625.924.249.90001.91416.5483745.319.435.30002.0510.5543876.416.916.690002.246.21303922.924.952.20001.9311.8584037.521.5411812.45.62.059.2334170.9510.818.250002.2361264270.58.5210002.216.61234376.38.515.20002.2066119Table 5
Correlations matrix.
CBRGrF200SLLPLPIMDDOMCCBR10.773−0.82−0.623−0.42−0.42−0.4030.716−0.437Gr0.7731−0.738−0.953−0.121−0.118−0.1270.727−0.46F200−0.82−0.73810.4990.5240.5160.536−0.7410.565S−0.623−0.9530.4991−0.081−0.081−0.078−0.60.337LL−0.42−0.1210.524−0.08110.9980.968−0.2820.26PL−0.42−0.1180.516−0.0810.99810.951−0.2750.252PI−0.403−0.1270.536−0.0780.9680.9511−0.2960.28MDD0.7160.727−0.741−0.6−0.282−0.275−0.2961−0.844OMC−0.437−0.460.5650.3370.260.2520.28−0.8441


## 1.4. Multinonlinear Regression

Nonlinear regression estimates the dynamic connection between the dependent variable and the function. Nonlinear multivariable regression estimates the interaction between independent and dependent variables. There are various modeling techniques and regression approaches for specific aspects, especially when data are limited [17].It is possible that the greatest association between the answer and independent variables may be explained by step-by-step regression. This algorithmic technique selects the proper model subsets to filter forward or backward. First, choose a permanent model and use the model words before optimizing fitness. That is, step-by-step regression is a forward selection procedure that examines the significance of all prior variables. For variables whose partial square sums do not meet the minimal requirement, the retroactive exclusion procedure is improved and variables are excluded one at a time until the minimum condition is met. Step-by-step regression requires more computations than forward or backward detailed calculations, but it yields better results [33]. Scheffe proposed a second-order polynomial for MNR analysis of mechanical characteristics [21]. MNLR in PQ forms has been studied.


## 1.5. Gene Expression Programming (GEP)

GEP, which is an algorithm based on genetic algorithm (GA) and genetic programming (GP), was ﬁrst proposed by Ferreira [27]. The GEP technique is an integration of two evolutionary strategies that are broadly acknowledged in the scientific community. The first approach is called the GA methodology, and it simplifies complex relationships by depicting them with linear structures called chromosomes that are of a fixed length. In addition, the second method is known as GP, and it makes use of the expression tree (ET) configurations that may take on a wide range of forms and dimensions [22]. Developing a program code recorded in linear chromosomes of fixed length, it exhibits some characteristics akin to biological evolution. GEP’s primary goal is to create a mathematical function that fits a collection of data supplied to the GEP model. The GEP method performs symbolic regression on the mathematical equation using the majority of GA’s genetic operators. GEP and GA have some distinctions. Individuals are shown in GA as symbolic strings of fixed length (chromosomes), but they are depicted in GP as nonlinear entities of varying sizes and forms (parse trees), whereas in GEP, each individual is represented by a fixed-length string that is used to make ETs of different sizes and shapes [9]. The primary GEP process is presented in Figure 5. Five variables are employed in GEP analysis: function set, terminal set, fitness function, control parameters, and stop condition. After encoding the issue into the candidate solution and specifying the fitness function, a randomly generated population of viable people (chromosomes) is constructed and then translated into an ET matching to a mathematical equation. Following that, the projected values are compared to the actual ones, and the fitness score for each chromosome is calculated. The algorithm is ended when a suitable fitness score is reached. If this is not the case, some chromosomes are randomly selected and then mutated to generate new generations. This method is repeated until a suitable fitness score is produced, at which point noncoding chromosomal analysis is undertaken to determine the optimal solution to the problem [14].Figure 5
Frequency histograms of the index, MDD, and OMC.


## 1.6. GEP Elements

The ET and chromosomes are the two most important components of the GEP. The chromosome is made up of one or more genes and it is represented by a mathematical equation on the genome. To code any mathematical equation, two bilingual and conclusive languages named Karva languages (the language of genes and the language of ET) are utilized, both of which are bilingual and conclusive. This property of the GEP is important for inferring the genotype with more precision than it would otherwise be possible. The GEP genes are divided into two sections, which are referred to as the head and tail. It is necessary to use the head of a gene, which contains certain mathematical functions, variables, and constants, for the encoding of a particular function. However, the tail of a gene that contains just variables and constants is used as supplemental terminal symbols in addition to the main terminal symbols. These symbols are used in the event that the terminal symbols in the head are insufficient to express a particular function. Many other functions, such as fundamental arithmetic operators, trigonometric functions, and any other mathematical or user-defined functions, may be utilized in the head of a gene.


## 1.7. GEP Operators

The major operators in GEP are selection, mutation, transposition, and crossover (recombination), which are all based on natural selection. By using these operators, the chromosomes may be adjusted to improve the fitness score of the following generation for survival. The operator rates, which are specified at the outset of the model construction, indicate a particular chance of a chromosome appearing in the model. In general, a mutation rate of 0.001 to 0.1 is suggested, with the lower limit being 0.001. The transposition operator and the crossover operator, on the other hand, are recommended to be 0.1 and 0.4, respectively [24].


## 1.8. Model Performance

To examine the performance of the prediction models, the RMSE, the coefficient of determination (R2), as well as mean absolute percentage error (MAPE), the scatter index (SI) [23], Bias, and index of agreement (IOA) are being used as shown in (4)–(9).(4)RMSE=∑i=1NCBRCBRm−CBRp2N,R2=1−∑i=1NCBRm−CBRp2∑i=1NCBRp2,MAPE=1N∑i=1NCBRm−CBRpCBRp×100,SI=1/N∑i=1NCBRpi−CBRp¯+CBRmi−CBRm¯21/N∑i=1NCBRmi,Bias=∑i=1NCBRp−CBRmN,IOA=1−∑i=1NCBRpi−CBRmi2∑i=1NCBRpi−CBRm¯+CBRmi−CBRm¯2.n, CBRm, and CBRp indicate the total number of data measured experimentally, and the developed predicted value, respectively.In addition, the Durbin–Watson statistic was used to determine the existence of multicollinearity [1]. The Durbin–Watson statistic has a range of values between 0 and 4, and the midpoint (that is, 2) indicates that there is no association between the input variables. Therefore, a value between 1.5 and 2.5 is sufficient for obtaining models that are unaffected by multicollinearity. This section covers mathematical expressions and models. Each model’s R2 and P-value (p) are described in the Tables 6 and 7.(5)DW=∑t=2Tet−et−1∑t=1Tet2,whereas DW is the Durbin–Watson test and et and et−1 least-squares regression residuals.


## 1.9. Uncertainty and Reliability

The main purpose of uncertainty analysis is to narrow the predicted range in which the real value of an experiment’s result resides. This estimated range is expressed as an interval and is referred to as the uncertainty interval. It may be approximated using the calculated errors for the measuring method of the experiment in question [28]. U95 is one of the uncertainty analysis algorithms used to calculate the uncertainty interval. The value of U95 relating to a specific experiment result may be read as follows: if the provided experiment is performed again, the real value of the experiment’s outcome will fall within the specified uncertainty range around 95 times out of every 100 trials. Furthermore, the value of U95 is provided by(6)U95=1.96N∑i=1NCBRmi−CBRm¯2+∑i=1NCBRmi−CBRpi2.


## 1.10. Relative Absolute Error and Reliability

The relative absolute error is a relative statistic that compares a predictive model’s performance to that of a simple model. The predictive model’s performance is defined as the total absolute difference between the actual and anticipated values (i.e., the error). The basic model’s performance is defined as the total absolute difference between the realized value and the average of all realized values [25]. In other words, the relative absolute error determines if a model outperforms just forecasting the average (i.e., the simple model). The relative absolute error formula is(7)RAE=∑i=1NCBRmi−CBRpi∑i=1NCBRmi−CBRm¯.The meaning of the relative absolute error is straightforward: if the RAE is less than one, the model outperforms the basic model. The relative absolute error for a perfect model is 0. RAE should ideally be as near to zero as feasible.A statistical technique known as reliability analysis may be used to evaluate a model’s level of general consistency. To be more specific, it determines whether or not a proposed model reaches a level of performance that is within acceptable limits. A measurement that may be used for the examination of reliability is defined as(8)Reliability=100%N∑i=1Nki.Then after we check to see if RAEi is greater thanδ, if it is, we set ki equal to 1; otherwise, we set it equal to 0. We define ki as the number of times the value of RAE is either less than or equal to that of δ. According to Chinese standards, the ideal value of δ is 0.2 [28].


## 2. Results and Discussion

### 2.1. Multivariate Models Development and Prediction of CBR

The data were fully analyzed for the relationship between CBR and the parameters, each one individually, followed by a multiple linear regression relationship between CBR and the parameters: grain size, the relationship between CBR and Atterberg limit, and the relationship between the compaction parameter and CBR. With the help of the correlation matrix, a certain parameter will be selected to get the best model representing this data. Each model will be evaluated based on the correlation coefficient of determination (R2) and its level of significance (P≤0.05).The correlation matrix in Table5 clearly illustrates the nature of the relationship between the CBR and the soil parameters individually. As was expected, the studied area is sandy in nature. Therefore, all Atterberg coefficients appear to have a poor correlation with CBR. On the other hand, the relationship between the granular size distribution parameters, especially F200, appears from the compaction parameters (MDD) such that the highest correlation is related to CBR.


### 2.2. Multilinear Regression

Several attempts were conducted to find the best combination of parameters that represent the CBR correlation with the soil index and compaction parameters based on the findings from the correlation matrix (Table5). Out of several models, a possible correlation was made among those who exceeded R2 of 0.5 and, at the same time, was statistically significant. Multilinear regression models (MLR) are presented herein. Table 6 presents the outcomes of the multilinear regressions. MLR 1 and 3 in Table 6 are controlled by the grain size distribution; those models revealed R2 of 0.753 and 0.74. Consequently, MLR 4 presents the compaction parameters MDD and OMC as input parameters to predict the corresponding CBR; it revealed R2 of 0.61. MLR 2, 5, and 7 incorporate the effect of all indices as predictors for CBR; a better performance is revealed from these models where their R2 is 0.811, 0.79, and 0.792 consequently. The combination of grain size parameters and compaction parameters in MLR.6 gives R2 of 0.775. Consistency parameters could be considered with no contribution in MLR.6; this may be attributed to the fact that the soil in this study is nonplastic due to low plasticity soil. The best model of multivariate linear regression model is shown in (8), Table 7.(9)CBR=−835+6.76LL−8.97PL+321MDD−15.6lnGr−32.6lnS−52.9lnF200,where R2 is 0.81, and P<0.05.Table 6
Multilinear regression analysis significant models∗.
Model idParametersModelR2DWMLR 1Ln (Gr), Ln (F200) and Ln (S)CBR=−60.3+10.8LnGr−66.60LnF200−17.8LnS0.7531.84738MLR 2Ln (Gr), Ln (F200), Ln (S), PL, LL, MDD and OMCCBR=−835+6.76LL−8.97PL+321MDD+9.76OMC−15.6LnGr−52.9LnF200−32.6LnS0.8111.80383MLR 3Gr, F200, and SCBR=−1039 + 12.2Gr+ 7.6F200+ 11.4S0.7401.82060MLR 4MDD and OMCCBR=−1571 + 717MDD + 14.84OMC0.6101.84245MLR 5Gr, F200, S, PL, LL, MDD and OMCCBR =−1483 + 9.5Gr+ 6.2F200+ 8.9S+ 4.26LL− 5.97PL+ 289MDD+ 8.92OMC0.7921.62133MLR 6Gr, F200, S, MDD and OMCCBR=−1588 + 10.5Gr + 6.5F200+ 10.0S + 294MDD + 9.25OMC0.7751.67606MLR 7Ln (Gr), Ln (F200), Ln (S), MDD, and OMCCBR=−813−7.6LnGr−58.4LnF200−21.3LnS+315MDD+9.59OMC0.7901.78415∗Significance whereas P<0.01.Table 7
PQ [multinonlinear] regression analysis significant models∗.
Model idParametersModelR2DWPQ 1Ln (Gr), Ln (F200) and Ln (S)CBR=−2213−972LnGr−673LnF200−1041LnS−222LnGr2−106.0LnF2002−246LnS20.7871.98443PQ 2Ln (Gr), Ln (F200), Ln (S) PL, LL, MDD and OMCCBR=11235+6.49LL−8.36PL−10817MDD+15.1OMC+150LnGr+50LnF200−80LnS+2616MDD2−0.766OMC2+53LnGr2+14.5LnF2002−35.8LnS20.851.32803PQ 3Gr, F200, and SCBR=−638 + 5.2Gr+ 0.2F200+ 10.8S+ 0.0366Gr2+ 0.136F2002− 0.0420S20.7941.85436PQ 4MDD and OMCCBR=13434 − 13348MDD+ 28.3OMC+ 3286MDD2− 1.265OMC20.7941.28769PQ 5Gr, F200, S, PL, LL, MDD, and OMCCBR=9184 + 5.1Gr− 4.8F200+ 2.6S+ 7.86LL− 10.01PL− 9284MDD+ 16.6OMC− 0.0074Gr2+ 0.211F2002+ 0.0145S2+ 2257MDD2− 0.790OMC20.8881.50942PQ 6Gr, F200, S, MDD and OMCCBR=9623+5.9Gr−3.7F200+3.5S−9764MDD+18.0OMC−0.0100Gr2+0.187F2002+0.0103S2+2370MDD2−0.867OMC20.8711.64803PQ 7Ln (Gr), Ln (F200), Ln (S), MDD, and OMCCBR=12059+236LnGr+89LnF200−8LnS−11425MDD+14.3OMC+72LnGr2+20.8LnF2002−22.6LnS2+2757MDD2−0.759OMC20.881.52370∗Significance whereas P<0.01.


### 2.3. PQ Models

The PQ which is none linear generally contains linear variables and constants; it is effective when such a pattern does not appear to be linear, and the relationships of the variables tend to be curvilinear [16]. Table 8 illustrates the PQ models for the prediction of CBR of granular soil. The model’s significance (P) value is zero, signifying robust models.Table 8
GEP parameters obtained for the best three models.
Chromosomes30Error type fitness functionRMSE andR2Count of the genes5Head size10Linking function+Function set+, -, ×, /, √Mutation rate0.044One-point recombination rate0.3The recombination rate of two points0.3Inversion rate and transposition rate0.1The best model determination coefficient (R2) is 0.888, which is better than MLR CBR model in (8). The best PQ regression model is shown in (9).(10)CBR=9184+5.1Gr−4.8F200+2.6S+7.86LL−10.01PL−9284MDD+16.6OMC000000−0.0074Gr2+0.211F2002+.0145S2+2257MDD2−0.79OMC2,where R2 is 0.888, and P<0.05.


### 2.4. GEP Model Development

The primary objective of developing the GEP model was to provide mathematical functions that might be used to predict the CBR value. Numerous GEP models with differing quantities of input variables were built in this study in order to estimate the CBR. Each model’s number of input variables was distinct. The GEP models use the same input parameters as the multivariate models, which are also shown in Table3. Similar to multivariate models, all GEP data were arranged into a format with seven input parameters, including percentages of gravel (Gr), sand (S), and fines (F200), as well as LL, PI, MDD, and optimal water content (OMC). To create the best GEP model, the number of chromosomes was altered between 20 and 35, the number of genes between 3 and 7, and the head size between 8 and 13. During the building of GEP models, it was observed that the same GEP parameters were gathered for the “four best models”; Table 8 lists these parameters.Results of the best models which have better statistical performance that are presented in (5)–(8), in which the input parameters and mathematical function are included. Model GEP1 has three inputs that include grain size distribution. Models incorporate the OMC, while model GEP3 includes the parameters compaction and grain size analysis with PI. Model 4 has all indices that can affect the value of CBR.(11)GEP.1.CBR=18.0824+2.073Gr−F200+GrF200−Gr×F200S,GEP.2.CBR=21.056+4GrOMC+2OMC−F200+7.135GrF200,GEP.3.CBR=3.13+Gr−S−PI+0.334S−8.964−6.1861−GrF200+OMC2×F200Gr,GEP.4.CBR=−148+98.8MDD+3.9PI−0.8603F200+Gr+OMC−0.89PL×MDD−MDD×F200+PIPL/OMC−3.882−23.61/F200.The IOA, the SI, and Bias are the three indices that are produced in order to carry out a full statistical evaluation of the models that have been offered for the GEP. Table9 displays the findings of the statistical analysis that was carried out. According to Table 10, GEP 1 and GEP 4 have the highest coefficient of determination (R2 of 0.776), although GEP4 has the least RMSE with a value of 22.552. According to Table 9, GEP 1 has the lowest SI of 0.252, which indicates that it has better GEP model performance than the other GEP models. Good indexes of agreement can be seen across all GEP models, with values ranging from 0.928 for GEP 4 to 0.984 for GEP 3. In addition, the predictions of CBR that were made available by GEP 1 had a higher reliability (53%) and uncertainty (15.78) in comparison to the predictions that were made available by GEP 4 (45% reliability and 16.766 uncertainty).Table 9
Statistical performances of the proposed GEP models.
Model idSIIOABiasGEP10.2522310.958626−33.5488GEP20.2563080.933823−33.7025GEP30.5226930.984099−46.5304GEP40.8198130.928315−16.7601Table 10
Performances of multivariate and GEP models.
ModelMAEMSERMSEMAPER2U95ReliabilityMLR118.177563.80323.74524.280.7515.94049MLR216.546432.26120.79124.510.8115.56951MLR319.208593.44324.36128.120.7416.03349MLR424.952890.20829.83641.580.6116.83343MLR517.783474.15921.77527.340.7915.69949MLR617.618513.07822.65126.830.7815.81649MLR716.698478.32621.87124.690.7915.70253PQ115.933485.93922.04418.940.7915.72355PQ211.560244.78315.64616.750.8915.02262PQ315.822470.64121.69421.350.7915.67857PQ416.843470.53021.69226.930.7915.67851PQ511.490256.20616.00616.170.8915.08457PQ611.744293.37217.12817.160.8715.17457PQ712.087273.34016.53318.000.8815.12762GEP118.251509.4422.55224.790.77615.78753GEP216.401523.86122.88820.310.77115.82853GEP316.878532.63623.078921.1950.76616.80345GEP418.471510.71522.59928.3930.77616.76645


### 2.5. Models Comparison

In this section, an examination, discussion, and comparison of the results were acquired from multilinear, PQ, and GEP models. The coefficient of determination (R2), RMSE, MAPE, and Durbin–Watson (DW) were employed as statistical verification tools to assess the accuracy of the created models' outputs. Table 10 summarizes the statistical performances of the multivariate and GEP models. According to Table 10, the majority of models exhibit acceptable agreement between the developed model and experimental values. Table 10 further shows that the performance of the majority of models derived by multilinear regression, PQ, and GEP techniques is acceptable in terms of R2, RMSE, MAPE, uncertainty, and reliability. This event can be regarded as the models' inputs being appropriately chosen. Table 10 presents the statistical performances of multivariate and GEP models.However, in statistics, the total error performance of a link between two groups may be assessed using coefficient of determination (R) values. According to Benjamin [12], if a given model provides an R value greater than 0.80, there is a strong correlation between the measured and predicted values for the whole available database of data. In addition, the statistical performance of any model is assessed using error criteria such as RMSE; MAPE, an important criterion; and R value, because a model with a high R2 value may have a high RMSE or MAPE value. Taking into consideration these criteria, the best multivariate models are MLR2, PQ2, and PQ5. However, the four GEP models have a higher R2 and a lower RMSE of 22.56. The error magnitudes of these models are sufficient for forecasting the CBR of soils. In addition, Model PQ5 is clearly the best model for both GEP and multivariate models based on all statistical measures. In Figures 6 and 7, the estimated CBR values generated by PQ5, of multivariate and GEP1 are visually contrasted with the observed CBR values. Figures 6 and 7 show that both the multivariate model and the GEP model make good predictions based on the experimental data.Figure 6
Comparison of the estimated PQ model and measured values of CBR.Figure 7
Comparison of the estimated GEP1 model and experimental values of CBR.An importance analysis of GEP is also performed to realize to what extent the input variables can affect the output CBR (see Figures8 and 9). Gravel% and Fines% are found to be the most contributing variables to CBR, followed by Sand%, as shown in model GEP 3 figure. While in model GEP 4, where the contribution of the consistency limit is incorporated into the model, it is found that the consistency parameters are affecting the CBR significantly, which is in agreement with the study performed by El-Ashwah et al. [18].Figure 8
Importance analysis for CBR GEP 3.Figure 9
Importance analysis for CBR GEP 4.


### 2.6. Sensitivity Analysis

To determine the quantitatively relative implication of each input variable on the CBR, the MLR, PQ, and GEP models were employed to carry out a sensitivity analysis. The analysis is executed such that one variable of 15(6) and (8) is taken out each time to assess the effect of that variable’s inclusion on output. For (5), the results of the analysis depict that F200 (R2 = 0.723, RMSE = 24.94, DW = 1.696) is the highest influence variable on the CBR, while GR in which (R2 = 0.808, RMSE = 20.95, DW = 1.718) has the least influence on the CBR equation based on multilinear regression. For (11), the MDD (R2 = 0.8162, RMSE = 20.474, DW = 1.805) is the highest influence variable on the CBR, as well as the Gr (R2 = 0.8868, RMSE = 16.06837, DW = 1.485). In (10), the most significant influencing parameter is OMC, where R2 = 0.721, RMSE = 25.522, DW = 1.767, while PI is the least influencing parameter, where R2 = 0.788, RMSE = 22.38, DW = 1.774. Statistical parameters resulting from sensitivity analysis of the best models are shown in Table 11.Table 11
Statistical results of the sensitivity analysis.
Best modelFactor missingRMSER2DWMLRGR20.949670.80751.71827F20024.939930.72731.69628S21.594080.79551.65002PL21.418790.79881.69307LL21.334640.80041.6908MDD22.761740.77281.81755OMC22.500850.7781.80682PQGR16.068370.88681.48538F20016.751090.8771.24835S16.100030.88631.4504PL16.951780.8741.55601LL16.876310.87511.53486MDD20.474080.81621.80475OMC16.455620.88131.68452GEPGR23.6490.7621.609F20023.55120.7621.609PI22.38140.7881.774MDD25.25130.73021.771OMC25.52210.7211.767


## 2.1. Multivariate Models Development and Prediction of CBR

The data were fully analyzed for the relationship between CBR and the parameters, each one individually, followed by a multiple linear regression relationship between CBR and the parameters: grain size, the relationship between CBR and Atterberg limit, and the relationship between the compaction parameter and CBR. With the help of the correlation matrix, a certain parameter will be selected to get the best model representing this data. Each model will be evaluated based on the correlation coefficient of determination (R2) and its level of significance (P≤0.05).The correlation matrix in Table5 clearly illustrates the nature of the relationship between the CBR and the soil parameters individually. As was expected, the studied area is sandy in nature. Therefore, all Atterberg coefficients appear to have a poor correlation with CBR. On the other hand, the relationship between the granular size distribution parameters, especially F200, appears from the compaction parameters (MDD) such that the highest correlation is related to CBR.


## 2.2. Multilinear Regression

Several attempts were conducted to find the best combination of parameters that represent the CBR correlation with the soil index and compaction parameters based on the findings from the correlation matrix (Table5). Out of several models, a possible correlation was made among those who exceeded R2 of 0.5 and, at the same time, was statistically significant. Multilinear regression models (MLR) are presented herein. Table 6 presents the outcomes of the multilinear regressions. MLR 1 and 3 in Table 6 are controlled by the grain size distribution; those models revealed R2 of 0.753 and 0.74. Consequently, MLR 4 presents the compaction parameters MDD and OMC as input parameters to predict the corresponding CBR; it revealed R2 of 0.61. MLR 2, 5, and 7 incorporate the effect of all indices as predictors for CBR; a better performance is revealed from these models where their R2 is 0.811, 0.79, and 0.792 consequently. The combination of grain size parameters and compaction parameters in MLR.6 gives R2 of 0.775. Consistency parameters could be considered with no contribution in MLR.6; this may be attributed to the fact that the soil in this study is nonplastic due to low plasticity soil. The best model of multivariate linear regression model is shown in (8), Table 7.(9)CBR=−835+6.76LL−8.97PL+321MDD−15.6lnGr−32.6lnS−52.9lnF200,where R2 is 0.81, and P<0.05.Table 6
Multilinear regression analysis significant models∗.
Model idParametersModelR2DWMLR 1Ln (Gr), Ln (F200) and Ln (S)CBR=−60.3+10.8LnGr−66.60LnF200−17.8LnS0.7531.84738MLR 2Ln (Gr), Ln (F200), Ln (S), PL, LL, MDD and OMCCBR=−835+6.76LL−8.97PL+321MDD+9.76OMC−15.6LnGr−52.9LnF200−32.6LnS0.8111.80383MLR 3Gr, F200, and SCBR=−1039 + 12.2Gr+ 7.6F200+ 11.4S0.7401.82060MLR 4MDD and OMCCBR=−1571 + 717MDD + 14.84OMC0.6101.84245MLR 5Gr, F200, S, PL, LL, MDD and OMCCBR =−1483 + 9.5Gr+ 6.2F200+ 8.9S+ 4.26LL− 5.97PL+ 289MDD+ 8.92OMC0.7921.62133MLR 6Gr, F200, S, MDD and OMCCBR=−1588 + 10.5Gr + 6.5F200+ 10.0S + 294MDD + 9.25OMC0.7751.67606MLR 7Ln (Gr), Ln (F200), Ln (S), MDD, and OMCCBR=−813−7.6LnGr−58.4LnF200−21.3LnS+315MDD+9.59OMC0.7901.78415∗Significance whereas P<0.01.Table 7
PQ [multinonlinear] regression analysis significant models∗.
Model idParametersModelR2DWPQ 1Ln (Gr), Ln (F200) and Ln (S)CBR=−2213−972LnGr−673LnF200−1041LnS−222LnGr2−106.0LnF2002−246LnS20.7871.98443PQ 2Ln (Gr), Ln (F200), Ln (S) PL, LL, MDD and OMCCBR=11235+6.49LL−8.36PL−10817MDD+15.1OMC+150LnGr+50LnF200−80LnS+2616MDD2−0.766OMC2+53LnGr2+14.5LnF2002−35.8LnS20.851.32803PQ 3Gr, F200, and SCBR=−638 + 5.2Gr+ 0.2F200+ 10.8S+ 0.0366Gr2+ 0.136F2002− 0.0420S20.7941.85436PQ 4MDD and OMCCBR=13434 − 13348MDD+ 28.3OMC+ 3286MDD2− 1.265OMC20.7941.28769PQ 5Gr, F200, S, PL, LL, MDD, and OMCCBR=9184 + 5.1Gr− 4.8F200+ 2.6S+ 7.86LL− 10.01PL− 9284MDD+ 16.6OMC− 0.0074Gr2+ 0.211F2002+ 0.0145S2+ 2257MDD2− 0.790OMC20.8881.50942PQ 6Gr, F200, S, MDD and OMCCBR=9623+5.9Gr−3.7F200+3.5S−9764MDD+18.0OMC−0.0100Gr2+0.187F2002+0.0103S2+2370MDD2−0.867OMC20.8711.64803PQ 7Ln (Gr), Ln (F200), Ln (S), MDD, and OMCCBR=12059+236LnGr+89LnF200−8LnS−11425MDD+14.3OMC+72LnGr2+20.8LnF2002−22.6LnS2+2757MDD2−0.759OMC20.881.52370∗Significance whereas P<0.01.


## 2.3. PQ Models

The PQ which is none linear generally contains linear variables and constants; it is effective when such a pattern does not appear to be linear, and the relationships of the variables tend to be curvilinear [16]. Table 8 illustrates the PQ models for the prediction of CBR of granular soil. The model’s significance (P) value is zero, signifying robust models.Table 8
GEP parameters obtained for the best three models.
Chromosomes30Error type fitness functionRMSE andR2Count of the genes5Head size10Linking function+Function set+, -, ×, /, √Mutation rate0.044One-point recombination rate0.3The recombination rate of two points0.3Inversion rate and transposition rate0.1The best model determination coefficient (R2) is 0.888, which is better than MLR CBR model in (8). The best PQ regression model is shown in (9).(10)CBR=9184+5.1Gr−4.8F200+2.6S+7.86LL−10.01PL−9284MDD+16.6OMC000000−0.0074Gr2+0.211F2002+.0145S2+2257MDD2−0.79OMC2,where R2 is 0.888, and P<0.05.


## 2.4. GEP Model Development

The primary objective of developing the GEP model was to provide mathematical functions that might be used to predict the CBR value. Numerous GEP models with differing quantities of input variables were built in this study in order to estimate the CBR. Each model’s number of input variables was distinct. The GEP models use the same input parameters as the multivariate models, which are also shown in Table3. Similar to multivariate models, all GEP data were arranged into a format with seven input parameters, including percentages of gravel (Gr), sand (S), and fines (F200), as well as LL, PI, MDD, and optimal water content (OMC). To create the best GEP model, the number of chromosomes was altered between 20 and 35, the number of genes between 3 and 7, and the head size between 8 and 13. During the building of GEP models, it was observed that the same GEP parameters were gathered for the “four best models”; Table 8 lists these parameters.Results of the best models which have better statistical performance that are presented in (5)–(8), in which the input parameters and mathematical function are included. Model GEP1 has three inputs that include grain size distribution. Models incorporate the OMC, while model GEP3 includes the parameters compaction and grain size analysis with PI. Model 4 has all indices that can affect the value of CBR.(11)GEP.1.CBR=18.0824+2.073Gr−F200+GrF200−Gr×F200S,GEP.2.CBR=21.056+4GrOMC+2OMC−F200+7.135GrF200,GEP.3.CBR=3.13+Gr−S−PI+0.334S−8.964−6.1861−GrF200+OMC2×F200Gr,GEP.4.CBR=−148+98.8MDD+3.9PI−0.8603F200+Gr+OMC−0.89PL×MDD−MDD×F200+PIPL/OMC−3.882−23.61/F200.The IOA, the SI, and Bias are the three indices that are produced in order to carry out a full statistical evaluation of the models that have been offered for the GEP. Table9 displays the findings of the statistical analysis that was carried out. According to Table 10, GEP 1 and GEP 4 have the highest coefficient of determination (R2 of 0.776), although GEP4 has the least RMSE with a value of 22.552. According to Table 9, GEP 1 has the lowest SI of 0.252, which indicates that it has better GEP model performance than the other GEP models. Good indexes of agreement can be seen across all GEP models, with values ranging from 0.928 for GEP 4 to 0.984 for GEP 3. In addition, the predictions of CBR that were made available by GEP 1 had a higher reliability (53%) and uncertainty (15.78) in comparison to the predictions that were made available by GEP 4 (45% reliability and 16.766 uncertainty).Table 9
Statistical performances of the proposed GEP models.
Model idSIIOABiasGEP10.2522310.958626−33.5488GEP20.2563080.933823−33.7025GEP30.5226930.984099−46.5304GEP40.8198130.928315−16.7601Table 10
Performances of multivariate and GEP models.
ModelMAEMSERMSEMAPER2U95ReliabilityMLR118.177563.80323.74524.280.7515.94049MLR216.546432.26120.79124.510.8115.56951MLR319.208593.44324.36128.120.7416.03349MLR424.952890.20829.83641.580.6116.83343MLR517.783474.15921.77527.340.7915.69949MLR617.618513.07822.65126.830.7815.81649MLR716.698478.32621.87124.690.7915.70253PQ115.933485.93922.04418.940.7915.72355PQ211.560244.78315.64616.750.8915.02262PQ315.822470.64121.69421.350.7915.67857PQ416.843470.53021.69226.930.7915.67851PQ511.490256.20616.00616.170.8915.08457PQ611.744293.37217.12817.160.8715.17457PQ712.087273.34016.53318.000.8815.12762GEP118.251509.4422.55224.790.77615.78753GEP216.401523.86122.88820.310.77115.82853GEP316.878532.63623.078921.1950.76616.80345GEP418.471510.71522.59928.3930.77616.76645


## 2.5. Models Comparison

In this section, an examination, discussion, and comparison of the results were acquired from multilinear, PQ, and GEP models. The coefficient of determination (R2), RMSE, MAPE, and Durbin–Watson (DW) were employed as statistical verification tools to assess the accuracy of the created models' outputs. Table 10 summarizes the statistical performances of the multivariate and GEP models. According to Table 10, the majority of models exhibit acceptable agreement between the developed model and experimental values. Table 10 further shows that the performance of the majority of models derived by multilinear regression, PQ, and GEP techniques is acceptable in terms of R2, RMSE, MAPE, uncertainty, and reliability. This event can be regarded as the models' inputs being appropriately chosen. Table 10 presents the statistical performances of multivariate and GEP models.However, in statistics, the total error performance of a link between two groups may be assessed using coefficient of determination (R) values. According to Benjamin [12], if a given model provides an R value greater than 0.80, there is a strong correlation between the measured and predicted values for the whole available database of data. In addition, the statistical performance of any model is assessed using error criteria such as RMSE; MAPE, an important criterion; and R value, because a model with a high R2 value may have a high RMSE or MAPE value. Taking into consideration these criteria, the best multivariate models are MLR2, PQ2, and PQ5. However, the four GEP models have a higher R2 and a lower RMSE of 22.56. The error magnitudes of these models are sufficient for forecasting the CBR of soils. In addition, Model PQ5 is clearly the best model for both GEP and multivariate models based on all statistical measures. In Figures 6 and 7, the estimated CBR values generated by PQ5, of multivariate and GEP1 are visually contrasted with the observed CBR values. Figures 6 and 7 show that both the multivariate model and the GEP model make good predictions based on the experimental data.Figure 6
Comparison of the estimated PQ model and measured values of CBR.Figure 7
Comparison of the estimated GEP1 model and experimental values of CBR.An importance analysis of GEP is also performed to realize to what extent the input variables can affect the output CBR (see Figures8 and 9). Gravel% and Fines% are found to be the most contributing variables to CBR, followed by Sand%, as shown in model GEP 3 figure. While in model GEP 4, where the contribution of the consistency limit is incorporated into the model, it is found that the consistency parameters are affecting the CBR significantly, which is in agreement with the study performed by El-Ashwah et al. [18].Figure 8
Importance analysis for CBR GEP 3.Figure 9
Importance analysis for CBR GEP 4.


## 2.6. Sensitivity Analysis

To determine the quantitatively relative implication of each input variable on the CBR, the MLR, PQ, and GEP models were employed to carry out a sensitivity analysis. The analysis is executed such that one variable of 15(6) and (8) is taken out each time to assess the effect of that variable’s inclusion on output. For (5), the results of the analysis depict that F200 (R2 = 0.723, RMSE = 24.94, DW = 1.696) is the highest influence variable on the CBR, while GR in which (R2 = 0.808, RMSE = 20.95, DW = 1.718) has the least influence on the CBR equation based on multilinear regression. For (11), the MDD (R2 = 0.8162, RMSE = 20.474, DW = 1.805) is the highest influence variable on the CBR, as well as the Gr (R2 = 0.8868, RMSE = 16.06837, DW = 1.485). In (10), the most significant influencing parameter is OMC, where R2 = 0.721, RMSE = 25.522, DW = 1.767, while PI is the least influencing parameter, where R2 = 0.788, RMSE = 22.38, DW = 1.774. Statistical parameters resulting from sensitivity analysis of the best models are shown in Table 11.Table 11
Statistical results of the sensitivity analysis.
Best modelFactor missingRMSER2DWMLRGR20.949670.80751.71827F20024.939930.72731.69628S21.594080.79551.65002PL21.418790.79881.69307LL21.334640.80041.6908MDD22.761740.77281.81755OMC22.500850.7781.80682PQGR16.068370.88681.48538F20016.751090.8771.24835S16.100030.88631.4504PL16.951780.8741.55601LL16.876310.87511.53486MDD20.474080.81621.80475OMC16.455620.88131.68452GEPGR23.6490.7621.609F20023.55120.7621.609PI22.38140.7881.774MDD25.25130.73021.771OMC25.52210.7211.767


## 3. Conclusion

This study’s primary objective is to investigate the applicability of multivariate regression analysis and gene expression programming (GEP) for predicting CBR. To do this, CBR test data for coarse-grained soils belonging to the A-1 and A-2 soil groups were supplied from various locations in the Qassim area. Seven MLR and PQ nonlinear regressions and four GEP models with different input variables were examined to determine the optimal connection between fundamental soil indices and the parameter CBR. The models’ performance was assessed using statistical verification criteria. The models MLR, PQ, and GEP with eight input parameters yielded the greatest results. It can be concluded that the MLR, PQ, and GEP are capable of acquiring the relationship between CBR and basic soil parameters and may be used to predict CBR values of soils. The results of this analysis showed that the PQ model had a higher accuracy [coefficient of determination (R2) = 0.89, RMSE = 16.006, uncertainty (U95) = 16.17, and reliability = 57%] than the multilinear regression model, which has a lower accuracy [R2 = 0.811, RMSE = 20.791, U95 = 15.569, and reliability = 51%]. The best GEP model yields [R2 = 0.776, RMSE = 22.552, U95 = 15.787, and reliability = 53%]. Furthermore, an important analysis revealed that gravel percent and fine percent have the greatest influence on CBR, followed by consistency parameters (LL, PL, PI) and OMC. Given how hard, time-consuming, and expensive CBR soil tests are, it can be emphasized that using MLR, PQ, and GEP models to estimate CBR for granular soils in terms of soil parameters could be a useful method for early stages of material development or as a way to judge the validity of CBR values. Compaction parameters (OMC and MDD) give a weaker correlation with the CBR of granular soils, but when used with gravel and fines percentage, they play a significant role in the determination of the CBR of granular soil in the Qassim area. A sensitivity analysis was conducted to yield that F200 is the most variable factor affecting the multilinear regression model, while MDD is the most variable factor affecting the PQ model, and the optimum moisture content (OMC) significant variable affects the best GEP model. The models used in this work to estimate CBR were based on granular-grained soils with low to nonplastic properties. As a result, it cannot accurately figure out CBR in fine soils or soils of high plasticity.



---

*Source: 7426962-2022-08-31.xml*