Evaluation of Contributing Factors Affecting Number of Vehicles Involved in Crashes Using Machine Learning Techniques in Rural Roads of Cosenza, Italy

: The evaluation of road safety is a critical issue having to be conducted for successful safety management in road transport systems, whereas safety management is considered in road transportation systems as a challenging task according to the dynamic of this issue and the presence of a large number of effective parameters on road safety. Therefore, the evaluation and analysis of important contributing factors affecting the number of vehicles involved in crashes play a key role in increasing the efﬁciency of road safety. For this purpose, in this research work, two machine learning algorithms, including the group method of data handling (GMDH)-type neural network and a combination of support vector machine (SVM) and the grasshopper optimization algorithm (GOA), are employed. Hence, the number of vehicles involved in an accident is considered to be the output, and the seven factors affecting transport safety, including Daylight (DL), Weekday (W), Type of accident (TA), Location (L), Speed limit (SL), Average speed (AS), and Annual average daily trafﬁc (AADT) of rural roads in Cosenza, southern Italy, are selected as the inputs. In this study, 564 data sets from rural areas were investigated, and the relevant, effective parameters were measured. In the next stage, several models were developed to investigate the parameters affecting the safety management of road transportation in rural areas. The results obtained demonstrated that the “Type of accident” has the highest level and “Location” has the lowest importance in the investigated rural area. Finally, although the results of both algorithms were the same, the GOA-SVM model showed a better degree of accuracy and robustness than the GMDH model.


Introduction
Road transport is one of the oldest methods of transporting goods and individuals. With the increase in the population and the development of towns, the expansion of road transport has become an inevitable issue that plays a significant role in increasing and improving economic and social development. Numerous issues affect the quality and quantity of road transport, and safety management is one of the most significant subjects. Therefore, identifying, studying, and evaluating the contributing factors affecting road safety is inevitable to increase the level of safety management. Several parameters affect road safety, and valuable studies have been conducted to investigate road safety, such as Driver Behavior [1][2][3][4][5][6][7], Age of Driver and Vehicle [8][9][10][11][12][13], Weather Conditions [14][15][16][17][18], Road Geometry [19][20][21][22][23], and Lighting Conditions [24][25][26][27][28][29]. Siliquini et al. (2010) investigated the link between driving performance and the use of psychoactive drugs in teen drivers who were enrolled at common consuming locations. For this purpose, they designed and implemented the TEND by Night project in several European countries. Based on their results, the key findings were the change in response time between entering and exiting the recreation site, as well as its relationship with psychoactive drug usage [30]. Alonso et al. (2017) studied a relatively large statistical community of drivers to explore the behavioral and representational characteristics of drivers that modulate the smokingaccidents relationship. They discovered that, despite being aware of the consequences of smoking while driving, few drivers thought it was a risky behavior. Finally, they suggested various recommendations based on their findings, such as increasing awareness and control of this habit [31]. Zoe et al. (2018) carried out an efficient analysis of the expansion trend in road safety research between 2000 and 2018. Their results indicated that road safety research focused on five major areas, including accident frequency data investigation, driver behavior questionnaire, safety in numbers for walkers and bicyclists, injury and prevention of road traffic, and driving speed and road accidents [32].
Moreover, other comprehensive studies have been performed to examine the role of other factors on road safety. Gichaga (2017) reviewed the historical and cultural backgrounds involving road development and road safety features in Kenya. Based on his results, he made some recommendations for improvements in aspects of road safety [33]. Elvik et al. (2019) conducted a review of the relationship between speed and road safety. They supported two mathematical models. Their results showed that the speed of individual drivers has a similar relationship with safety as the mean speed of traffic [34]. In another study, the effect of air quality on road safety was evaluated by Sager (2019). He investigated the impact of increased air pollution on the amount of road traffic crashes. Then, he found out that there is a relationship between the number of accidents and the amount of PM2.5 [35].
On the other hand, several techniques exist to evaluate the parameters affecting road safety. Traditionally, classical before-after research works, statistical modeling, and personal judgment-based approaches are applied to chronological data [36]. Zheng et al. (2018) investigated the proposal of a bivariate extreme value model. Their proposed model could reduce the uncertainty in crash estimates [37]. Moreover, according to a bivariate extreme value theory (EVT) framework, Wang et al. (2019) introduced an accident forecast method. Their results prove that the proposed crash prediction method can provide very promising results compared to univariate models [38]. In another research work, Zheng et al. (2020) carried out a comprehensive review of research based on traffic conflicts in the road safety analysis. Then, they discussed conceptual and methodological matters related to traffic conflict modeling. Based on their results, it has been determined that although suitable research studies have been performed on this issue, the need for more studies is necessary [39]. An overview of road traffic accidents was conducted in the Eastern Province, KSA, from 2009 to 2016 by Jamal et al. (2020). They developed logistic regression models to forecast crash severity. Finally, they recommended some suggestions to prevent road accidents [40].
Due to the uncertainty and unforeseen conditions of parameters affecting road safety and the high ability of machine learning algorithms to solve complex problems, recently, the use of these methods in combination with classical methods or, in many studies, alone has been widely used [41][42][43][44][45][46][47][48][49]. Xu et al. (2018) evaluated the effect of road lighting on road safety using an artificial neural network. The results clearly showed that road lighting greatly contributed to road safety levels [50]. In another study, two types of artificial intelligence methods were applied to estimate the severity of crashes by Amiri et al. (2020). They used the crash information from California in 2012. Their results indicated that both artificial intelligence models were reliable for predicting the severity of crashes, and also, the light condition played an essential role in the severity level compared to the other contributing factors [51]. Guido et al. (2020) investigated some potential factors for road safety using two artificial intelligence methods. Their results indicated that the PSO algorithm had a superior function compared to the GA algorithm for evaluating factors in road accidents [52]. Shiran et al. (2021) carried out crash severity analysis by applying data mining approaches and multinomial logistic regression. They used an accident dataset of State Highways in California, USA. They found that the C5.0 model provides a higher performance capacity in evaluating crash severity analysis than other models [53]. A brief review of the literature is shown in Table 1. Table 1. A brief overview of some studies that employed machine learning techniques [41][42][43][44][45][46][47][48][49][50][51][52][53].

Researcher(s) Type of Techniques Description
Mussone et al. [41] ANN Modeling of urban vehicular accidents using ANN with an assessment of the main circumstances and causes of accidents.
Halim et al. [42] Review of AI techniques such as: GA, GP, CRF, ANN, PCA, Fuzzy Logic, TD Learning, SVM Assessment of studies based on AI approaches for accident prediction and identifying dangerous driving situations.
Castro et al. [43] Bayesian Network, Decision Trees, and ANN Evaluation of the impact of various factors on injury risk in order to improve the road safety level.
De Luca [44] MVA, ANN A comparison of road safety management prediction models on two-lane highways.
Shah et al. [45] DEA-ANN Identification and evaluation of the most important criteria in determining the level of road risk.
Liu et al. [46] ANFIS, Logistic Regression, Decision Tree, and SVM Examining real-time crash risk for urban freeways as a means of assessing road safety and traffic control decisions.
Guido et al. [47] GMDH Assessment of the effective parameters affecting accidents for the urban and rural areas.
Mokhtarimousavi et al. [48] SVM, CS-SVM, and Logit Model Temporal examination of accident severity determinants in worker-involved work zone crashes based on random parameters and machine learning methodologies.
Kitali et al. [49] SVM-FA Examination of the elements that influence the severity of injuries in crashes on express lanes facilities.
Xu et al. [50] ANN Study on the impact of road lighting on traffic safety.
Amiri et al. [51] ANN, GA-ANN Forecasting the severity of fixed object accidents among elderly drivers using two types of AI techniques.
Guido et al. [52] GA, PSO The use of clustering models to evaluate potential safety factors.
Studies of the past literature show that although valuable studies have been done, more research is needed to improve the quality of road transport based on the increasing expansion. Therefore, in this study, two machine learning algorithms, including the group method of data handling (GMDH)-type neural network and a combination of support vector machine (SVM) and the grasshopper optimization algorithm (GOA), are applied. Then, the obtained results of two algorithms are compared based on performance indicators to determine the performance of the models regarding the conditions and characteristics of a case study. Finally, a sensitivity analysis was conducted to prioritize the contributing factors affecting the number of vehicles involved in crashes in rural areas.

Site Description and Accident Monitoring
To evaluate the performance and usefulness of the proposed approaches for prioritizing the contributing factors affecting the number of vehicles involved in crashes in a rural area, a sample of 564 accident data was acquired from 2017 to 2018 on the rural road network of Cosenza province and has been analyzed ( Figure 1). It should be noted that the most comprehensive data sets were available for modeling in this period.
The accident data has been acquired from the Automobile Club Italia (ACI) database, which collaborates with the National Institute of Statistics (ISTAT) in collecting road accident data in Italy. This database provides a lot of information about accidents, such as date and place, road category, pavement conditions, weather conditions, crash dynamics, the type of vehicle involved, the causes of the accident, and the consequences for the people involved (injuries or deaths). The information does not include Property Damage Only (PDO) events under existing Italian legislation, defining road accidents as accidents only when they cause at least one injury [54]. According to ISTAT, accidents are classified as fatal or with injuries; therefore, it is not possible to distinguish the injured according to the level of severity. Based on the previous comment, in this paper, the measure of the number of vehicles involved in an accident, as referred to in Section 4, is expressed as a dichotomous variable (0,1): -When there is a vehicle involved in an accident, label 0 is considered. -When there is more than one vehicle involved in an accident, label 1 is considered.
Other information on roads and traffic flows was considered to ensure a more detailed analysis, such as speed limits, average speed, and average annual daily traffic (AADT). The speed limits are regulated by the "Nuovo Codice della Strada" [55,56], and their values have been acquired from the database of Azienda Nazionale Autonoma delle Strade (ANAS). The average speed was obtained through the available data of the historical traffic statistics of TomTom (TomTom Move) and Octo Telematics (Octo IoT Cloud), referring to the road sections with the observed accidents. The AADT was acquired from the PANAMA system of ANAS (ANAS Platform for Monitoring and Analyzing).
This study considers seven independent variables (i.e., the factors affecting the number and severity of accidents), which were selected based on all the available data for the case study. These variables include four qualitative variables (Daylight (DL), Weekday (W), Type of accident (TA), and Location (L)) and three quantitative variables (Speed limit (SL), Average speed (AS), and Annual Average Daily Traffic (AADT)). Table 2 shows the variables mentioned above, which are classified into various categories (codes) with their statements. The method of grouping variables into different categories depends on their characteristics and the number of observations involved in each study.

Methodology
The study of road safety is one of the inseparable issues of transportation engineering, which is usually defined by the absence of accidents and casualties. Lack of attention to road safety could impose irreparable financial and physical damages. Therefore, it is necessary to have a deep understanding of road safety, know all the effective components, and estimate the impact of each on this issue. Studying the literature reveals that most investigations in the field of road safety are based on logit models or regression models with artificial neural networks. On the one hand, the complexity and uncertainty of the factors affecting road safety, and on the other hand, the ability of machine learning algorithms to predict and navigate in the face of unexpected and uncertain issues, has resulted in the successful application of machine learning methods in road safety in recent years [57][58][59][60][61][62][63].
Accordingly, the main aim of this study is to investigate the factors affecting road safety in rural areas by using two machine training approaches, namely the GMDH model and the hybrid GOA-SVM model. For this purpose, the GMDH model was developed to achieve the best binary classification model by determining the best control parameters for GMDH. Moreover, the SVM was hybridized with the grasshopper optimization algorithm as a suitable evolutionary algorithm and was developed to optimize its three parameters. Finally, a sensitivity analysis was performed to evaluate and rank factors affecting the number of vehicles involved in an accident in rural areas. More discussions regarding machine learning models will be given in the next sections.

Group Method of Data Handling-Type Neural Network
In today's science, ANNs (as one of the branches of artificial intelligence) play a valuable role in the development of new technologies. Therefore, the use of ANNs to solve many complex problems in various fields of science is inevitable [64][65][66][67][68][69][70][71]. The group method of data handling (GMDH) is one of the artificial neural networks that was introduced by Ivakhnenko (1971). The GMDH has been successfully used for computer-based mathematical modeling in complex systems, as well as data mining, optimization, and pattern recognition problems [72]. The process of GMDH is similar to a type of self-organizing network [73,74]. The mapping between the input and output variables in a GMDH neural network is a nonlinear function. The Polynomial Neural Network (PNN) is one of the basic algorithms used to construct GMDH models. In GMDH modeling, input data are entered into the initial layer, and, after preparation, they are considered input for the second layer, and this process continues until the algorithm converges and stops. Finally, in the convergence process of the algorithm, if the results in the layer (n + 1) are better than the layer (n), then the algorithm converges. Equations (1) and (2) indicate the relationship between the approximate function ( ∧ f ) with the multi-input and single-output ( ∧ y) dataset and the least possible error between actual and predicted values [73,74].
where a i , a ij , a ijk , a ijkl , . . . are considered to be coefficients of the polynomial, and m is the amount of data. Furthermore, Equation (3) can be considered as a quadratic polynomial for 2 inputs according to Equation (4) [77].
The total error (E) is considered by minimizing the difference between the actual output (y) and predicted output ŷ = G (x i , x j ) for each pair of input variables x i and x j based on Equation (5). Moreover, the coefficients of each quadratic function are optimized [78].
Out of a total of n input variables, all alternatives for two independent variables are provided in the elementary form of the GMDH algorithm for providing the regression polynomial in the form of Equation (4) [79]. Therefore, n 2 = n (n−1) 2 neurons will be made in the primary layer of the feedforward neural network from the observations . The matrix form of Equation (4) can be considered to indicate the main form of the GMDH based on Equation (6) [80]. Y = Aa (6) Y = (y 1 , y 2 , y 3 , . . . , y m ) and a = (a 1 , a 2 , a 3 , a 4 , a 5 ) are considered the observed output vector and the vector coefficient of the quadratic polynomial, respectively. A is computed based on Equation (7) [81,82].
Finally, using the least-squares process from multiple regression analysis, a normal equation is achieved based on Equation (8), which calculates the vector of the best coefficients for Equation (4) [83].

Support Vector Machine
A support vector machine (SVM) is an effective machine learning method that was introduced by Cortes and Vapnik (1995) [84]. SVMs are a kind of supervised learning algorithms that are used in a wide range of modelings, such as regression and classification. The SVM presents a linear two-class classifier and aims to maximize the margin amongst two classes so that a classification hyperplane is formed in the center of the maximum margin. It provides many hyperplanes, while the goal of the support vector machine is to find the best hyperplane in n-dimensional space. Two labels are considered for this classification: label +1 is considered for cases that are above the hyperplane, and label −1 belongs to cases that are under the hyperplane. Equation (9) shows a group of the sample set that is used in classification learning data [85,86].
where y i is the target variable for the observed i-th sample (the sample category). It is also assumed that x i presents the i-th sample data. After the formation of hyperplanes, one of the hyperplanes has the highest margin, which is called the optimal hyperplane. This optimal hyperplane is determined by the existing support vectors and constraints.
where w and b are the weight vector and the bias vector, correspondingly. Then, considering an error coefficient, the constraints are rewritten and corrected according to Equations (12) and (13). This error coefficient is intended to ensure a more accurate classification [88]. where c is the penalty coefficient. Then, using the Lagrange method, SVM classification problems are considered as the following dual optimization problem based on Equation (14) [88].
where K is a mathematical function that is called the kernel function. There are different types of kernel functions, including linear (LIN), the radial basis function (RBF), and polynomial (POL), and their relationships are shown in Table 3. Gamma (γ) and d are necessary to define the kernel types. Gamma (γ) is used for RBF and POL and "d" represents the term of polynomial degree only for the POL kernel function [89,90]. The most important role of the kernel function is to take the dataset as an input and convert it into the required form. Knowledge of the use of various kernel functions in related situations can affect the quality of the category. Table 3. Equations of different kernel functions.

Grasshopper Optimization Algorithm
In recent years, meta-heuristic algorithms have played a very important role in dealing with complex and uncertain problems [91][92][93][94][95][96]. These algorithms have made significant progress in both academia and industry. The grasshopper optimization algorithm (GOA) is one of the newest meta-heuristic algorithms, which was presented by Saremi et al. (2017) [97]. This algorithm is based on swarm intelligence and is population-based, which was inspired by the group behavior of grasshoppers. Grasshoppers are usually seen individually or in groups. One of the essential features of the group of grasshoppers is their type of movement, which has a slow movement with small steps. In the GOA, two sections are defined for the search, including exploration and exploitation. In exploration, search agents are persuaded to move suddenly, whereas they want to move locally in exploitation. A mathematical model was introduced to simulate the swarming behavior of locusts according to Equation (15) [97,98].
where X i represents the position of the i-th grasshopper. S i , G i , and A i are the social interaction, the gravity force on the i-th grasshopper, and the wind advection, respectively. Then, to show the random behavior of grasshoppers in Equation (15), the random factors are used including r 1 , r 2 , and r 3 , which are random values within [0, 1]. Equation (15) is rewritten according to Equation (16) [97,99].
The gravity force on the i-th grasshopper (S i ) is presented based on Equation (17) [97]. where d ij andd ij present the distance between the i-th and j-th grasshopper and a unit vector from the i-th grasshopper to the j-th grasshopper, which is computed as d ij = x j − x i and , respectively. Based on Equation (18), S is a function to describe the strength of social forces. The motion of the grasshopper is affected by the repulsion and attraction factors between them, which is defined by s in the mathematical model of GOA based on Figure 2. According to Figure 1, there is a comfort zone (comfortable distance) in which neither the attraction nor the repulsion action takes place between two grasshoppers [97].
where f and l are the intensity of attraction and the attractive length scale, correspondingly. Changes in the amount of f and l indicate a change in the behavior of the grasshopper, so by changing them, the number of S changes, and the final results can change.
Equations (19) and (20) show the gravity force on the i-th grasshopper and the wind advection, respectively [97].
whereê g is a unity vector towards the center of earth, and g is the constant of gravity. Moreover, u andê w present a constant drift and a unity vector in the direction of the wind, correspondingly. It should be noted that nymph grasshoppers do not have wings, so their movements are strongly influenced by the wind. By inserting the values of each of these definitions, the mathematical model of the algorithm expands based on Equation (21) [97][98][99][100].  (19) and (20) show the gravity force on the i-th grasshopper and the wind advection, respectively [97]. It should be noted that if the locust population reaches the comfort zone quickly, the swarm will not converge to a specific point, and Equation (21) is not able to solve optimization problems directly. Therefore, the mathematical model presented in Equation (21) considers an upper bound and lower bound as well as two coefficients to balance the motion between the comfort, gravity, and repulsion zones, and Equation (21) is modified as Equation (22) [97,101].
where ub d and ub d represent the upper bound and lower bound in the D-th dimension, respectively.T d presents the value of the D-th dimension in the target, and c introduces a reducing coefficient to shrink the comfort zone, repulsion zone, and attraction zone that is shown in Equation (23) [97,102].
where cmax and cmin are the maximum value and the minimum value, correspondingly. Moreover, l demonstrates the current iteration, and L shows the maximum number of iterations. Parameter c needs to be reduced proportionally to the number of iterations needed to balance exploration and exploitation. For more information and explanations regarding the grasshopper optimization algorithm, refer to Saremi et al. (2017) [97].

Results and Discussion
As mentioned earlier, for the binary classification modeling, two machine learning algorithms, namely GMDH and the hybrid GOA-SVM, were used and developed, and the number of vehicles involved in the accidents was evaluated in this study. For this purpose, a valuable dataset was collected and, as described in Section 2, the seven factors affecting the number of vehicles involved in the crashes, including the Daylight (DL), Weekday (W), Type of accident (TA), Location (L), Speed limit (SL), Average speed (AS), and Annual average daily traffic (AADT), of rural areas of Cosenza in southern Italy, were considered. It should be noted that, initially, a comprehensive study was conducted on the literature, and a set of parameters affecting road safety was identified; about 18 parameters. Then, due to the limitations of data access, including incomplete data, lack of data, and incorrect data, these seven parameters affecting road safety were selected. In this study, all the data has been classified into two classes. To check the number of vehicles involved in an accident, the first class with the label "0" was considered for cases where at most one car was involved in an accident. The second class was labeled "1" for cases in which at least two or more vehicles were involved in a crash. This classification was based on the assumption that the main criterion for class separation is to consider the minimum number of vehicles involved in an accident. This was conducted by developing and constructing the best classification model to determine the correct classes with the highest possible accuracy by determining a mapping between the input and output data.
By developing models, the best models for each method are determined, and then the results obtained are compared. Finally, by performing a sensitivity analysis, the importance of the effect of each of the factors is determined. It is necessary to mention that in binary classification modeling, the use of accuracy and error in the confusion matrix are considered the most practical performance indicators. Therefore, to compare the performance evaluation of the models, the confusion matrix is used according to Figure 3 and Equations (24) and (25). Moreover, due to the range of changes and the scale of measurement of each of the studied parameters, the normalization of this data has a significant role in the data-driven system modeling approaches with appropriate accuracy because if they are not normalized, a factor on a larger scale may cause a computational deviation. Therefore, all data are normalized using the min-max normalization before modeling.

GMDH Modeling
In this work, GMDH was employed to construct a binary classification model f assessment of the number of vehicles involved in an accident using the MATLAB ronment. The number of vehicles involved in an accident was considered the depe variable, and DL, W, TA, L, SL, AS, and AADT were set as the independent var Hence, the optimal architecture of GMDH contributes greatly to the high performa GMDH models. Thus, the accurate determination of control parameters of the G model is an essential issue. Although there are no specific relationships to accurate termine control parameters, most of these parameters are obtained from past studi pert opinions, and trial and error. Therefore, according to experts' opinions and pre studies, a range is considered for some control parameters of the GMDH model, includes the set of the maximum number of layers (MNL) equal to 5, 10, 20, 40, a and the maximum number of neurons in a layer (MNNL) equal to 5, 10, 20, 40, a Moreover, the Selection Pressure (SP) is another control parameter of the GMDH m and it is a dimensionless number where the sensitivity of the modeling error is af by the SP, and it was considered equal to 0.5. Different rates for training and testin are considered in modeling, and several studies have been conducted in this [103,104]. Based on the suggestions and studies of Looney, 75% of the data (423 were used for training, and 25% of the data (141 cases) were applied as test data for eling [105]. In total, 25 models were constructed, and their results are demonstra Table 4.

GMDH Modeling
In this work, GMDH was employed to construct a binary classification model for the assessment of the number of vehicles involved in an accident using the MATLAB environment. The number of vehicles involved in an accident was considered the dependent variable, and DL, W, TA, L, SL, AS, and AADT were set as the independent variables. Hence, the optimal architecture of GMDH contributes greatly to the high performance of GMDH models. Thus, the accurate determination of control parameters of the GMDH model is an essential issue. Although there are no specific relationships to accurately determine control parameters, most of these parameters are obtained from past studies, expert opinions, and trial and error. Therefore, according to experts' opinions and previous studies, a range is considered for some control parameters of the GMDH model, which includes the set of the maximum number of layers (MNL) equal to 5, 10, 20, 40, and 50, and the maximum number of neurons in a layer (MNNL) equal to 5, 10, 20, 40, and 50. Moreover, the Selection Pressure (SP) is another control parameter of the GMDH model, and it is a dimensionless number where the sensitivity of the modeling error is affected by the SP, and it was considered equal to 0.5. Different rates for training and testing data are considered in modeling, and several studies have been conducted in this case [103,104]. Based on the suggestions and studies of Looney, 75% of the data (423 cases) were used for training, and 25% of the data (141 cases) were applied as test data for modeling [105]. In total, 25 models were constructed, and their results are demonstrated in Table 4.
After constructing different models and determining the accuracy performance of each model, all models were ranked based on a simple ranking method suggested by Zorlu et al. (2008), Table 5 indicates the ranking results [106].  According to Table 5, all models were ranked, and the seven developed models had the lowest rank, with a rank equal to 31 among the 25 developed models. While the 15th model had the best performance in comparison with other developed models. Its training accuracy and testing accuracy were 83.2% and 81.6%, correspondingly. The structure of the 15th developed model consisted of MNL, MNNL, and SP equal to 20, 50, and 0.5, respectively. The results of the confusion matrices of the training, the testing, and total data are indicated in Figure 4a-c, respectively. Table 5, all models were ranked, and the seven developed models had the lowest rank, with a rank equal to 31 among the 25 developed models. While the 15th model had the best performance in comparison with other developed models. Its training accuracy and testing accuracy were 83.2% and 81.6%, correspondingly. The structure of the 15th developed model consisted of MNL, MNNL, and SP equal to 20, 50, and 0.5, respectively. The results of the confusion matrices of the training, the testing, and total data are indicated in Figure 4a-c, respectively.  As mentioned earlier, 75% of the data set (423) was considered for the training dataset, and the rest (141) was assigned to the testing dataset. According to Figure 4a, the 15th binary classification model correctly identified 80 cases of the first class with the label "0" (at most one car was involved in an accident), while 18 cases of the second class (at least two or more vehicles were involved in an accident) were incorrectly estimated in the first class with the label "0". Note that it could classify the training dataset with an accuracy equal to 83.2%. Meanwhile, the results of the binary classification for the testing data according to Figure 4b show that 36 cases of the first class with label "0" and 79 cases of the second class with label "1" were correctly considered, while the 15th developed model wrongly predicted 5 and 21 cases of the first and second classes with labels "0" and "1", respectively. Consequently, the 15th model was able to obtain an acceptable accuracy of classifying test data of 81.6%. Finally, according to the binary classification results for the whole data set according to Figure 4c, it is clear that the 15th binary classification developed model could correctly classify 467 cases of two classes, and 97 cases of both classes were wrongly classified. Consequently, the total accuracy of the binary classification modeling was 82.8%. This analysis shows that GMDH is a reliable modeling method for predicting the number of vehicles involved in an accident.

GOA-SVM Modeling
GOA and SVM were combined to develop a predictive model in the MATLAB environment. The GOA algorithm was used to optimize some parameters of the SVM so that the SVM model shows the highest performance. In the GOA-SVM modeling of this study, the same datasets performed in the analysis of GMDH were used. For modeling, 75% of the dataset (423 data) is randomly defined as the training dataset, and the remaining 25% (141 data) is considered as the testing dataset [105]. Like modeling with the GMDH algorithm, the two classes were considered for all data with labels, including "0" and "1". To develop and optimize the parameters of the SVM model by the GOA, the control parameters of the GOA must be determined, which play an important role in the rapid and appropriate convergence of the model. Although there are no specific relationships to determine these parameters, based on previous studies and experts' opinions, a range was determined for each of them, such as the number of grasshopper populations (5, 10, 15, 20, 30, and 40) and the number of iterations (10, 20, 40, 50, and 100). Then, the most appropriate ones were selected by a trial-and-error approach [107].
Moreover, to further evaluate the model, k-fold cross-validation was used, in which the data were subdivided into K subsets. In this system, one of them was used at each time for validation, while the other K − 1 was applied for training. This procedure was carried out K times, with each data set being used exactly once for training and once for validation. Finally, the average result of this K validation was chosen as the final estimate. There is no specific method for determining the amount of k-fold, and it is determined according to the number of data and the opinion of experts. Hence, k-fold was considered to be equal to 3 in this study [107]. In addition, the three different types of kernel functions, including RBF, POL, and LIN, were used. According to the number of control parameters, 30 models were built for each kernel function of GOA-SVM. That means, in total, 90 models were made for GOA-SVM. Preliminary analyses were performed, and a comparison of the performance indicators for the best-developed models with three different kernels is shown in Figure 5.  Although the GOA algorithm was able to train SVM very well with different kernels, according to Figure 5, it is clear that the best-developed model with an RBF kernel had the highest degree of accuracy in training, validation, and testing in comparison with the POL and LIN kernel functions. Therefore, the best developed GOA-SVM model is considered with the RBF kernel function and the optimal control parameters of the model are shown in Table 6. Moreover, the value of error in each iteration was calculated based on Equation (25), and the result is indicated in Figure 6. According to Figure 6, modeling started with an error of about 0.195 and had different values until the 18th iteration, reaching 0.153 in the 19th iteration, which remained constant until the last iteration (40th). The obtained values of accuracy and error indicated the proper convergence of the model. Although the GOA algorithm was able to train SVM very well with different kernels, according to Figure 5, it is clear that the best-developed model with an RBF kernel had the highest degree of accuracy in training, validation, and testing in comparison with the POL and LIN kernel functions. Therefore, the best developed GOA-SVM model is considered with the RBF kernel function and the optimal control parameters of the model are shown in Table 6. Moreover, the value of error in each iteration was calculated based on Equation (25), and the result is indicated in Figure 6. According to Figure 6, modeling started with an error of about 0.195 and had different values until the 18th iteration, reaching 0.153 in the 19th iteration, which remained constant until the last iteration (40th). The obtained values of accuracy and error indicated the proper convergence of the model.   According to Figure 7, the best GOA-SVM model indicated higher performance than the best GMDH model for predicting the number of accidents by 84.6% and 83.4% for

Comparison of Models' Performance and Sensitivity Analysis
Each year, many people pass away due to road traffic accidents. Therefore, knowing the impact of various contributing factors on road accidents and taking the necessary measures to reduce accidents can significantly increase the level of road safety. This research employed two machine learning methods, namely GMDH-type neural network and GOA-SVM, to conduct the binary classification modeling. Based on the accuracy of the modeling performance, the best GMDH and GOA-SVM models were chosen after multiple modeling. A comparison was made between the best model of GMDH and the best model of GOA-SVM based on the accuracy of training and testing that is shown in Figure 7. According to the explanations in the previous section, it should be noted that the value of the validation accuracy model is considered instead of the value of training accuracy in the GOA-SVM model. Hence the value of validation accuracy of the GOA-SVM model should be compared to the value of training accuracy of the GMDH model.
According to Figure 7, the best GOA-SVM model indicated higher performance than the best GMDH model for predicting the number of accidents by 84.6% and 83.4% for training and testing accuracies, in comparison with 83.2% and 81.6% for training and testing accuracies, respectively. Although it is worth mentioning that both models had acceptable degrees of accuracy and robustness, it can be concluded that they are reliable systems of modeling for predicting the number of vehicles involved in an accident and can be used as useful tools for modeling road safety involved in transportation engineering.  According to Figure 7, the best GOA-SVM model indicated higher performance than the best GMDH model for predicting the number of accidents by 84.6% and 83.4% for training and testing accuracies, in comparison with 83.2% and 81.6% for training and testing accuracies, respectively. Although it is worth mentioning that both models had acceptable degrees of accuracy and robustness, it can be concluded that they are reliable systems of modeling for predicting the number of vehicles involved in an accident and can be used as useful tools for modeling road safety involved in transportation engineering.
Road accidents can cause considerable economic and human losses to society. Therefore, assessing the impact of parameters affecting the number of vehicles involved in an Road accidents can cause considerable economic and human losses to society. Therefore, assessing the impact of parameters affecting the number of vehicles involved in an accident can provide an in-depth insight for engineers involved in road safety management. A sensitivity analysis was performed to assess the impact of DL, W, TA, L, SL, AS, and AADT on the predicted number of vehicles involved in the accident. This sensitivity analysis is based on the cosine amplitude method according to Equation (26), in which r ij is the strength of the relationship, n shows the number of datasets, and x ik and y ij explain the input variables and the predicted output, correspondingly [108,109].
According to Figure 8, it is clear that the analyses of both models had similar results from the impact of the factors under consideration, which indicates the reliability of the results. Furthermore, the following remarks can be concluded: - The type of accident was the most significant factor among other contributing factors that affected the number of vehicles involved in the crashes. In general, certain types of accidents can be caused by a variety of issues, including a lack of traffic signs and poor road quality. The type of accident has an important effect on the number of vehicles involved in an accident. - The next factor is the average speed, which can increase the risk of accidents. Some researchers discovered that controlling other factors, such as traffic volume, road geometry, and the number of lanes, can reduce or eliminate the effects of average speed [110,111]. When the average speed is higher, the driver's response time is shorter, which can lead to an accident. Therefore, it is possible to control the impact of average speed by providing some types of measures, such as improvement to the location of road signs, speed limit enforcement methods, pavement markings, and vertical centerline treatments. - The third factor influencing the number of vehicles involved in an accident after the type of accident and the average speed is the annual average daily traffic that plays a key role in the development needs and priorities of road development for transportation planning. Moreover, some studies indicate that increasing the amount of AADT can lead to an increase in the frequency of accidents [112,113]. Therefore, to reduce the effects of AADT on this case study, it is recommended to consider other intercity transportation systems, such as trains, which are being considered in the review of coming urban development plans. - The subsequent contributing factor is the speed limit, which has a significant role in the behavior and decisions of drivers. Generally, it is considered that the speed limit is determined by the road conditions. If the speed limit is selected incorrectly on a part of the route and the driver is aware of this error due to the road conditions, they may lose confidence in the speed limits in other sections of the road and increase or decrease the speed based on their interpretation [114,115]. Hence, given the impact of this factor in this case study, it is suggested that a general review be considered in selecting the speed limit for rural roads in Cosenza. -Several appropriate studies have been conducted on the relationship between weekdays and accidents [116,117]. In line with national statistics, road accidents are more concentrated on holidays on the road network in Cosenza's province. In the present study, out of the 564 accidents, 65% (367) occurred on holidays. Due to the geographical location of Cosenza, the amount of traffic on holidays has experienced a relative increase. To reduce the exposure to the risk, the intensification of controls and monitoring of roads during the holidays would mitigate the effect of this factor by increasing police enforcement. -Extensive studies have been conducted on the effect of daylight on the number of vehicles involved in accidents, which shows this factor's high importance. This factor has been given priority in many studies, among other factors [50,118]. The amount of impact this factor has is heavily influenced by its geographic location and road lighting systems. This factor was determined as the sixth most effective factor out of seven factors, and this result was matched with the location and road lighting system of rural roads in Cosenza. - The last studied factor was the location that had the most negligible impact on the rate of vehicles involved in an accident, based on the results of both artificial intelligence models. Based on the type of structure of the rural roads in Cosenza, the location has not had much effect on the rate of accidents. Therefore, in future studies on this case study, the effect of this factor can be ignored. has been given priority in many studies, among other factors [50,118]. The amount of impact this factor has is heavily influenced by its geographic location and road lighting systems. This factor was determined as the sixth most effective factor out of seven factors, and this result was matched with the location and road lighting system of rural roads in Cosenza. - The last studied factor was the location that had the most negligible impact on the rate of vehicles involved in an accident, based on the results of both artificial intelligence models. Based on the type of structure of the rural roads in Cosenza, the location has not had much effect on the rate of accidents. Therefore, in future studies on this case study, the effect of this factor can be ignored. Human behavior is jointly responsible for road accidents, together with the vehicle and the environment (of which the infrastructure is part). It is strictly correlated to the speed and type of accident, and the author's analysis of the importance and the proposed methodology determined the impacts. Moreover, these two factors have the highest impact compared to other parameters on the rate of vehicles involved in accidents on the Human behavior is jointly responsible for road accidents, together with the vehicle and the environment (of which the infrastructure is part). It is strictly correlated to the speed and type of accident, and the author's analysis of the importance and the proposed methodology determined the impacts. Moreover, these two factors have the highest impact compared to other parameters on the rate of vehicles involved in accidents on the rural roads of Cosenza.
Therefore, based on the results obtained and after consulting experts, it can be concluded that several measures, such as improving roadway geometry, pavement markings, and vertical centerline treatments, are the most significant measures that should be considered in reviewing the plan for rural roads in Cosenza. Finally, it should be noted that the presented models with specific structures are location-sensitive and cannot be directly used on other rural roads.

Conclusions
The evaluation and analysis of contributing factors affecting the number of vehicles involved in the crashes can lead to a deep understanding of the existing situation and increase road safety by planning and taking a series of necessary measures. Therefore, this study made an attempt to predict the number of vehicles involved in the crashes using the GMDH-type neural network and GOA-SVM methods. This study was accomplished using 564 accident cases from rural roads in Cosenza in southern Italy. Several crash-related parameters, including DL, W, TA, L, SL, AS, and AADT, were set as input variables, and the number of vehicles involved in the crash was considered as the output data. According to each technique's number of control parameters, 25 models were made for GMDH and 30 models for each kernel function of GOA-SVM, so a total of 90 models were made for GOA-SVM. After modeling, the models for each technique were compared with each other in terms of the model's performance. Among the binary classification of GMDH models, the 15th model was selected with the highest score, and the GOA-SVM model with RBF kernel function was chosen among the binary classification of GOA-SVM models. In the comparison between the GMDH and GOA-SVM models, the GOA-SVM model had a higher capability for the prediction of the number of vehicles involved in the accident, although the performance of both models indicates that both can be used as useful tools for modeling the number of vehicles involved in an accident in transportation engineering. Consequently, a sensitivity analysis was performed based on the results obtained from both models. In both models, the type of accident and location had the highest and lowest impact compared to other parameters on the rate of vehicles involved in accidents on the rural roads of Cosenza, respectively. The results of this study showed the role of humans in causing accidents as the most contributing factor in road safety, which was fully consistent with previous studies in this field. As a result, it is recommended that organizations concerned with road safety, in addition to taking the required steps to enhance road condition, establish a long-term strategy for raising awareness and assessing driver performance.
For future work, it is recommended to examine other factors that can impact the number of vehicles involved in an accident, such as the age of drivers, the age of cars, gender, and the geometry of roads.