INTRODUCTION
In general asphalt binder properties are mostly described by investigating the rheology of the asphalt material. Rheology is a very powerful tool for establishing the properties of different materials that affect the pavement performance^{1}. Rheology is the study of the matter flow, which describes the elastic and viscous behavior of asphalt when loaded. Dynamic shear rheometer (DSR), which is one of the superpave test procedures, is a powerful tool used to characterize the viscoelastic behavior of asphalt binders^{2}.
As the asphalt binder becomes stiffer, the complex shear modulus goes higher and this allows binders to resist deformation and consequently hot mix asphalt (HMA) rutting. On the other hand, to resist fatigue cracking in asphalt pavement, the binder should not be too stiff because extremely stiff asphalts tend to crack rather than deform. The Superpave grading system uses the DSR outcomes; complex shear modulus (G*) and phase angle (δ) in the form of rutting (G*/sinδ) and fatigue cracking (G*×sinδ) parameters^{3}. Rutting parameter needs to be maximized by decreasing the phase angle and increasing the complex shear modulus in order to prevent rutting. As a result, the asphalt binder is able to recover deformation as it occurs and eventually becomes more rutting resistant. Conversely, fatigue cracking prevention requires the fatigue cracking parameter to be minimized by decreasing the phase angle and the complex shear modulus. This leads to a more elastic asphalt material that is not too stiff, which in turn makes the material more cracking resistant.
To fully characterize asphalt binders’ rheological properties, its performance needs to be tested over a wide range of temperatures that the asphalt is expected to be subjected to during its life. This is due to the fact that one of the main problems associated with asphalt pavement is asphalt binders’ sensitivity to temperature variation^{4}. Asphalt binders is also influenced by aging which can be split into two periods; shortterm aging that occurs during mixing, hauling and construction processes of the asphalt binder life cycle. This is simulated in the lab using the rolling thin film oven (RTFO) test. The second period is the longterm aging that occurs over the pavement service life and simulated using the pressure aging vessel (PAV) testing procedure^{5,6}.
Asphalt binders are generally graded by three grading systems; these are: Penetration grading system, viscosity grading system and superpave performance grading (PG) system. The PG grading system is focused on the expected climatic conditions where the asphalt binder is going to be used. A PG grade is defined by two numbers (both in °C): the average sevenday maximum pavement design temperature (PGhigh) and the minimum pavement design temperature (PGlow).
Based on the above, several factors that influence the rheological properties of asphalt binders were used as predictors; these factors are mostly related to testing conditions and asphalt binders’ properties. The factors are: test temperature (T), type of aging (Ag), low performance grade (PGlow), high performance grade (PGhigh), penetration (Pen), kinematic viscosity (KinV) and absolute or dynamic viscosity (AbsV). Literature lacks a comprehensive study that takes into consideration all of these factors in predicting the dynamic shear rheometer (DSR) complex shear modulus outcome. This study aims at filling this gap in the literature using the principle component analysis (PCA) approach.
The main objective of this study was to develop predictive models to predict the DSR complex shear modulus (G*) based on data collected from the longterm pavement performance (LTPP) website. This study discovers the relationship between the DSR complex shear modulus and a complete list of controlling factors that can be beneficial for researchers in the paving materials domain to save effort and time to acquire laboratory results and consequently cost. This study will also help the researcher to uncover the critical factors in governing the viscoelastic behavior of asphalt materials that many researchers were not able to fully explore. Thus a new robust predictive models of the DSR G* may be arrived at.
MATERIALS AND METHODS
Asphalt binders produced in the United States were used in this study. These asphalt binders were subjected to the Superpave DSR test method to obtain their complex shear modulus then results were stored in the LTPP website. Results of many years of laboratory testing were compiled from the LTPP website and analyzed using different analytical techniques such as the nonlinear regression (NLR), artificial neural network (ANN) and principle component analysis (PCA).
Data collection and processing: Data used in this study was obtained from the LTPP website, which is considered as part of the strategic highway research program (SHRP). Data compiled in this website represent the outcome of a large number of laboratory tests throughout the United States over the past years.
Table 1:  Variables definition and description for the complex shear modulus 


Dataset of complex shear modulus was linked and combined with the input variables using specific fields that are state code, sharp ID, construction number, sample number and layer number to end up having a complete dataset that could be used in the modeling of the complex shear modulus. After that the new dataset was filtered so that no missing data points and outliers and to obtain a proper dataset for the subsequent steps. The final dataset contained 627 records of complex shear modulus values. Table 1 lists the seven input factors that were selected based on previous work and the availability of data on LTPP website. The DSR test was conducted at several temperatures and a fixed frequency of 10 rad/sec. The type of aging was converted into a quantitative measure so that it can be used in the modeling of the complex shear modulus.
Research procedure: Four predictive models were developed in this study to investigate the Superpave DSR complex shear modulus. Modeling was performed using three analytical techniques (NLR, ANN and PCA). Neural Designer software (version 3) was used to perform the ANN analysis. The NLR was performed using the Statistical Package for the Social Sciences (SPSS; version 25) software and the PCA was performed using the JMP (version 14) data analysis software.
Nonlinear regression (NLR): In this process, the best mathematical function from a pool of existing functions is selected between each predictor and the response variable (in this case G*) and then the initial parameters are suggested. Afterwards, he best fit functions with their parameter estimates for each predictor are multiplied by each other using the SPSS nonlinear tool to predict G*. The MarquardtLevenberg algorithm was utilized here; this process is iterative and works to minimize the value of the squared sum of the difference between observed and predicted data^{7}. The target sum of squared adopted in this study is 1.00E8 or lower.

Fig. 1:  Details of the uses of all the instances in the dataset for complex modulus 
Artificial neural network (ANN): ANNs are braininspired systems that contain several interconnecting neurons which are programmed to mimic the human nervous system which can learn and interact with the environment and show responses. The process of the ANN with connection weights and biases is explained mathematically by Eq. 1, where f_{h} and f_{o} are the activation functions for hidden layer and output layer, respectively. B_{o} and B_{Hj} are the biases for output layer and hidden layer, respectively. W_{ij} is the weight connection between each one of inputs and neurons of hidden layer, W_{j} is the weight connection between the output layer and each neurons of hidden layer, X_{i} is the parameter input and Y_{p} is the predicted output:
The dataset used to build the ANN models was split into training, validation (selection) and testing subsets as shown in Fig. 1.
Typically, 65% of the dataset is used for training and 20% is used for validation and the remaining 15% for testing. The training subset is used to compute errors and gradients and to optimize the weights and biases at each interaction. The validation subset is used to test the performance of the model fit while training the model and evaluating the generalization of the running neural network. The testing subset is used after the training is finished to check the performance ability of the trained neural network.
The network architecture consists of input layer with eight neurons that are the type of aging that was split into two variables due to its qualitative nature: Original (ORG) and RTFO, test temperature (T), penetration (Pen), kinematic viscosity (KinV), absolute viscosity (AbsV), high performance grade (PGH) and low performance grade (PGL). Also, one neural hidden layer and one output layer were included in the network architecture as shown in Fig. 2.
Principle component analysis (PCA): PCA is one of the statistical analysis methods that can be used in data analysis and in building predictive models. Although it is very useful, this analytical technique is underutilized. PCA is a mathematical procedure that convert a number of correlated variables into a smaller number of uncorrelated variables called principle components (PCs)^{8}. PCA can serve in different ways in the science of data analysis, such as simplification, data reduction, modeling, outlier detection, variable selection, classification, prediction, etc.^{9}.

Fig. 3:  Architecture of PCR and PCNN models for the prediction of G* 
A study enumerated some usages of the obtained PCs including the establishment of a new set of orthogonal variables that are linear combinations of the original variables and at the same time contain exactly the same information as the original dataset, solving multicollinearity problem in the original dataset and the reduction of the ndimensions of the problem by finding out the most effective number of dimensions over which the data set exhibits variation^{10}.
Equation 2 is a mathematical expression of a principal component, a linear combination of the original input variables. Where is the original input variable, is the coefficient for each input, is a constant and is the number of the original input variables:
Using the selected PCs as input variables with the complex shear modulus as the response variable, the PCR analysis and the PCNN methods were used to develop the predictive models. The architecture of the PCAbased regression analysis and PCAbased neural network analysis are shown in Fig. 3. JMP software (version 14) was used for the PCA modeling.
RESULTS AND DISCUSSION
In this study, seven factors were used to investigate the Superpave DSR rheological properties of asphalt binders in terms of the complex shear modulus. Predictive models of the complex shear modulus were developed using three techniques; NLR, ANN and PCA.
Nonlinear regression (NLR): Descriptive statistics in terms of mean, standard deviation and number of data points for all seven factors are shown in Table 2. Aging factor is a qualitative variable that needs to be converted or coded so that it can be used in models’ development; 0 for Original (ORG) and 1 for RTFO.
Table 2:  Descriptive statistics of the complex shear modulus’s model 


Table 3:  Person’s correlation matrix of the complex shear modulus’s model 


Table 4:  Model evaluation of complex shear modulus model 


A matrix of Person’s correlation was utilized to estimate the correlation between each pair of independent variables and G*. Table 3 shows this correlation matrix; the testing temperature has the highest linear correlation among all the seven factors, followed by aging then penetration, PGhigh, absolute viscosity, kinematic viscosity and last is PGlow. In general, the correlation coefficients are classified from weak to moderate.
The curve estimation tool of the SPSS (version 25) was used to select the best fit between each independent variable and G* and then suggest the initial parameters. Table 4 represents each function used based on the highest R^{2} for each independent variable. Several nonlinear functions were used such as growth, power, cubic and compound. Also, the initial parameter estimate is shown for each independent variable and labeled in the table as “Model Evaluation #1”.
After many iterations utilizing the MarquardtLevenberg algorithm, the best fit functions with their parameter estimates for each independent variable are multiplied by each other using the SPSS nonlinear tool to predict the response variable (G*). The process in this study started with 3.211E+13 sum of squared residual and stopped after 92 iterations at sum of squared residual of 1.24E3 as shown in Table 4.
Equation 3 represents the developed NLR model for the prediction of obtained by the nonlinear regression analysis:
Figure 4 graphically displays the effectiveness of the developed NLR G* model by plotting the observed versus predicted values of G*.

Fig. 4:  Observed versus predicted complex shear modulus using nonlinear regression 

Fig. 5:  Normal probability plot of the residuals for complex shear modulus NLR model 

Fig. 6:  Network losses for each iteration 
The best fit (solid) line, the equality line (red dotted line) and R^{2} are all labelled on the figure. As shown, R^{2} was found to be 74.4%. The normal probability plot of the residuals for the NLR G* model is presented in Fig. 5, the figure indicates a normal distribution of the residuals.
Artificial neural network (ANN): After splitting the dataset into 65% training, 20% validation and 15% testing, an operation of minimummaximum scaling method was carried out resulting in a dataset with [μ = 0, σ = 1]. The dataset was normalized to be in the range of [1, 1] using the activation function called the hyperbolic tangent function. A supervised multilayer feedforward neural network model with hyperbolic tangent function in the hidden layers and linear function in the output layer was applied and the quasiNewton method was used here as the training algorithm.
The best performance of the neural network model was obtained by a repeated process. Thus, the best performance after much iteration was given by a neural network model that consists of one hidden layer containing 5 neurons. The architecture of this neural network can be written as 8:5:1, where the number of inputs in the scaling layer is 8, the number of outputs is 1 and the complexity is 5.
The loss index plays an important role in the use of a neural network. It provides a measure of the quality of the model. Figure 6 shows the losses in each iteration.

Fig. 7(ab): 
ANN model, (a) Observed versus predicted complex shear modulus and (b) Observed scaled versus predicted scaled complex shear modulus 
Table 5:  Error method used for complex shear modulus predictive model accuracy 


The initial value of the training loss is 4.41 and the final value after 208 iterations is 0.19. The initial value of the selection loss is 4.72 and the final value after 208 iterations is 0.20. The stopping criteria was the gradient norm goal with elapsed time of 3 sec. The final gradient norm parameter, which is the gradient of the loss function, has a value of 0.000758 and that indicates stability in the neural network model’s performance.
The normalized squared error (NSE) is used here and it has a value close to zero as shown in Table 5, which indicates that the model fits the dataset. Other types of error are also shown in Table 5, these are sum squared error (SSE), mean squared error (MSE) and root mean squared error (RMSE). Also, the model represents good generalization ability in which training, validation and testing dataset portions have close prediction accuracy.
This analysis leads to a high R^{2} of 80.74% between the observed and predicted G* values. Figure 7 shows the effectiveness obtained using the neural network model, once with the actual data points (Fig. 7a) and once with the scaled (normalized) data points (Fig. 7b). The normal probability plot of the residuals is presented in Fig. 8 which indicates a normal distribution of the residuals.

Fig. 8:  Normal probability plot of the residuals for complex shear modulus ANN model 
Table 6:  Mathematical expression of one hidden layer for the ANN G* model 


Table 7: 
Basic network validation results for nontraining dataset 

The mathematical expressions represented by the neural network are written below. It propagates the eight inputs in a feedforward model through the scaling layer into the perceptron layers to predict the G*. Table 6 shows the expressions for one hidden layer containing 5 neurons, all neurons were treated by a hyperbolic tangent activation function.
Equation 4 and 5 show the final scaled complex shear modulus (G*_{s}) and predicted complex shear modulus (G*_{p}) mathematical expressions; where y_{1}y_{5} represent the outputs of the five neurons in the hidden layer:
G*_{s} = (1.438560.100435×y_{1}+0.613896×y_{2}+0.246846×
y_{3}1.40126×y_{4}0.265329×y_{5})
 (4) 
G*_{p} = (0.5×(G*_{s}+1)×(15.50.3)+0.3 (5)  (5) 
Two nontraining data points were used after the model has been built to ensure the effectiveness of the developed neural network model. Both data points have low error value as shown in Table 7, indicating a reliable neural network model.
Principle component analysis (PCA) results: Starting with the G* dataset, the Person’s correlation matrix of the seven input variables is calculated and presented in Table 3. As shown, correlation coefficients with absolute value of 0.4 or greater are in boldface. Five of the numbers in the matrix were equal to 0.4 or greater indicating a high intercorrelation. So, to predict the G* accurately, PCA was used to replace the original input variables with a smaller set of principal components that are orthogonal (uncorrelated).
Using the JMP software, the eigenvalues and their corresponding eigenvectors are calculated and presented in Table 8. It is shown that 90.915% of the variation in the original data is explained by the first five principal components. The scree plot in Fig. 9 shows the location of the break (elbow) point which occurs after selecting the first five principal components. Based on the PCA results, the first five principal components were selected to create new uncorrelated input variables that are linear combination of the original variables.
Table 8:  Eigenvalues and eigenvectors of the complex shear modulus dataset Z matrix 


Table 9:  Model evaluation of the PCR complex shear modulus model 


The five selected principal components to create new input variables are shown in Eq. 6 through 10:
PC1 = 6.11360.05717×T+0.0632×KinV+0.03833×AbsV0.00813
×Pen +0.058094×PGL+0.0794×PGH+0.3086×Ag (6)
 (6) 
PC2 = 2.32050.01565×T+0.0868×KinV+0.05116×AbsV0.00333
×Pen+0.091313×PGL+0.04691×PGH+0.6116×Ag (7)
 (7) 
PC3 = 2.567230.02537×T+0.0472×KinV+0.02988×AbsV+
0.0043×Pen+0.01214×PGL+0.012225×PGH+1.7773×Ag (8)
 (8) 
PC4 = 2.562530.080023×T+0.0415×KinV+0.03077×AbsV+
0.01073×Pen+0.057858×PGL+0.04625×PGH+0.18335×Ag (9)
 (9) 
PC5 = 6.0144+0.07813×T+0.0667×KinV+0.034551×AbsV
+0.0115×Pen+0.0055142×PGL+0.01547×PGH+0.58802×Ag (10)
 (10) 

Fig. 9:  Scree plot of the complex shear modulus dataset 
These principal components will be used as new predictors while keeping the complex shear modulus as the response variable for modeling using two different techniques; the principle component regression (PCR) and the principal component neural network (PCNN).
Principle component regression (PCR) model: Results of the PCR model are presented in this section. The curve estimation tool of the SPSS program was used in which the most appropriate mathematical function was selected between each principal component and the response variable (G*). Table 9 represents each function used based on the highest R^{2} for each principal component. The initial parameter estimate is shown for each principal component and labeled in the table as “Model evaluation #1”.

Fig. 10:  Observed versus predicted complex shear modulus of the PCR model 

Fig. 11:  Normal probability plot of the residuals for complex shear modulus PCR model 
After many iterations utilizing the MarquardtLevenberg algorithm, the best fit functions with their parameter estimates for each principal component are multiplied by each other using the SPSS nonlinear tool to predict the response variable (G*). The process in this study started with 6.5E+7sum of squared residual and stopped after 150 iterations at sum of squared residual of 1.6E+3 as shown in Table 9.
Equation 11 represents the predicted complex shear modulus which was obtained by the nonlinear regression of the selected five principal components:
G* _{p} = [2.033+0.309×PC10.032×PC1 ^{2}0.021×PC1 ^{3}]×
[0.668+0.034×PC20.039×PC2^{2}]×[EXP
[0.0480.563×PC3)]×[EXP(0.156+0.371×PC4)]× (11)
[2.8922.693×PC5+0.206×PC5+0.229×PC5^{3}
 (11) 
The results showed that R^{2} between the predicted and observed complex shear moduli was 67.44%. Figure 10 graphically displays the ability of the developed PCR G*model to explain the variability in the data. The normal probability plot of the residuals for the PCR G* model is presented in Fig. 11, that indicates a normal distribution of the residuals.
Principle component neural network (PCNN) model: Results of the PCNN method are presented below. Least squares method that works to minimize the sum of squared errors is used for the artificial neural network. The network architecture of the developed PCNN model consists of a threelayer supervised feed forward network (5:5:1) that was developed using the JMP software as shown in Fig. 12. A hyperbolic tangent function in the hidden layers and linear function in the output layer were applied.
The connection weights and biases of the PCNN hidden layer of five neurons are presented by the following matrices.

Fig. 12:  Network architecture of the complex shear modulus PCNN model 

Fig. 13:  Observed versus predicted complex shear modulus of the PCNN model 
Where B_{o} and B_{Hj} are the biases for output layer and hidden layer, respectively. W_{ij} is the weight connection between the principal components and neurons of the hidden layer and W_{j} is the weight connection between the output layer and each one of the five neurons of the hidden layer:
Performance results for the training and validation datasets are presented by R^{2} and RMSE. The 418 data points that represent the training dataset has R^{2} of 81.49% with RMSE of 1.238, whereas the 209 data points that represent the validation dataset has R^{2} of 76.05% with RMSE of 1.257.

Fig. 14:  Plot of the complex shear modulus model residuals of the PCNN model 

Fig. 15:  Normal probability plot of the residuals for complex shear modulus PCNN model 
The coefficient of determination obtained using the PCNN model is 80.04% which indicates the amount of variability in the dataset that can be explained by this model. This is graphically shown in Fig. 13. Also, the residuals versus predicted plot is shown in Fig. 14 which represents the check of the homoscedasticity assumption; since there is no pattern in the residuals plot and it seems to be scattered, the assumption of the equality of variances is satisfied. Figure 15 shows the normal probability of the PCNN residuals for checking the normality assumption, the data points are not too far away from a straight line. Therefore, the normality assumption is satisfied. All modeling assumptions pertaining to the PCNN model are satisfied in this study.
The developed models using the previously mentioned techniques in this study revealed that as the testing temperature increases the complex shear modulus decreases. This finding is in agreement with several other studies^{2,3,4,1113}.
For aging variation, the complex shear modulus increases with aging. These results point out the increase in stiffness of binders which effect the asphalt flow properties and this leads to improving the rutting resistance (G*/sinδ) in contrast with fatigue resistance (G*×sinδ). This finding was also pointed out by several other studies^{1,5,6,1418}.
Kinematic and absolute viscosities have the same influence on the complex shear modulus; that is the increase of both viscosities leads to an increase in the complex shear modulus and this causes the binder to become stiff^{13,1921}. Complex shear modulus decreases with higher penetration values, which is in line results from other studies^{21,22}.
CONCLUSION
Results show that all techniques give good to very good predictive power, but the best to predict the G* is the ANN with an R^{2} of 80.74%. The prediction of the complex shear modulus based on the PCA is capable of producing more robust and reliable predictive models than other techniques when the original input parameters are highly correlated. Finally, the developed models may help in explaining the performance of asphalt pavements with less effort and time and consequently cost by limiting the use of fullscale lab measurements.
ACKNOWLEDGMENT
Authors would like to extend thanks and appreciation to the those who are responsible for upkeeping the LTPP website to help researchers investigate specific pavement related areas that are critical to pavement performance.