ahernmultiobjectiveextensivehypothesis2024.md 204 KB


category: literaturenote citekey: ahernmultiobjectiveextensivehypothesis2024 title: Multi-objective extensive hypothesis testing for the estimation of advanced crash frequency models authors: "Ahern, Zeke; Corry, Paul; Rabbani, Wahi; Paz, Alexander" year: 2024 date: 2024-00-00 2024 doi: 10.1016/j.aap.2024.107690 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2024.107690" zotero_key: TXZ5BBFW zotero_storage: GMSQBKXX collections: transpordi_planeerimine folder: 001_artiklid/Road Safety Indicators firstAuthor: "Ahern, Zeke"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

journal homepage: www.elsevier.com/locate/aap

Multi-objective extensive hypothesis testing for the estimation of advanced crash frequency models

Zeke Ahern a,∗ , Paul Corry b , Wahi Rabbani c , Alexander Paz a

  • a School of Civil & Environment Engineering, Queensland University of Technology, 2 George Street, Brisbane, 4000 QLD, Australia
  • b School of Mathematical Sciences, Queensland University of Technology, 2 George Street, Brisbane, 4000 QLD, Australia
  • c Department of Transport and Main Roads, Brisbane, 4000 QLD, Australia

A R T I C L E I N F O

Keywords: Metaheuristic Regression Random parameters Hypothesis testing Crash data Prediction

A B S T R A C T

Analyzing crash data is a complex and labor-intensive process that requires careful consideration of multiple interdependent modeling aspects, such as functional forms, transformations, likely contributing factors, correlations, and unobserved heterogeneity. Limited time, knowledge, and experience may lead to oversimplified, over-fitted, or misspecified models overlooking important insights. This paper proposes an extensive hypothesis testing framework including a multi-objective mathematical programming formulation and solution algorithms to estimate crash frequency models considering simultaneously likely contributing factors, transformations, non-linearities, and correlated random parameters. The mathematical programming formulation minimizes both in-sample fit and out-of-sample prediction. To address the complexity and non-convexity of the mathematical program, the proposed solution framework utilizes a variety of metaheuristic solution algorithms. Specifically, Harmony Search demonstrated minimal sensitivity to hyperparameters, enabling an efficient search for solutions without being influenced by the choice of hyperparameters. The effectiveness of the framework was evaluated using two real-world datasets and one synthetic dataset. Comparative analyses were performed using the two real-world datasets and the corresponding models published in literature by independent teams. The proposed framework showed its capability to pinpoint efficient model specifications, produce accurate estimates, and provide valuable insights for both researchers and practitioners. The proposed approach allows for the discovery of numerous insights while minimizing the time spent on model development. By considering a broader set of contributing factors, models with varied qualities can be generated. For instance, when applied to crash data from Queensland, the proposed approach revealed that the inclusion of medians on sharp curved roads can effectively reduce the occurrence of crashes, when applied to crash data from Washington, the simultaneous consideration of traffic volume and road curvature resulted in a notable reduction in crash variances but an increase in crash means.

1. Introduction

Crash count models are essential tools for examining likely crash contributing factors and underlying unobserved heterogeneity (Geedipally and Lord, 2010; Mannering and Bhat, 2014). Analysis using these models offer valuable insights for decision-making regarding risk management, policy development, and resource allocation (AASHTO, 2011). This, in turn, contributes to the development of effective countermeasures that play a pivotal role in preventing crashes, saving lives, and mitigating crash severity (Lozano et al., 2012; Mohammed et al., 2019; Wijnen et al., 2019; Cerwick et al., 2014). The field of crash count modeling offers a diverse range of analytical and data-driven models, each with its own strengths and limitations, influenced by the characteristics of the data (Mannering et al., 2016). Recent literature in this field encompasses various methodologies, including machine learning (Mannering et al., 2020), advanced heterogeneous statistical models (Behara et al., 2021), causality models (Tanishita and van Wee, 2017) and simplified statistical methods (Srinivasan and Bauer, 2013; Das et al., 2021b). All of these methodologies strive to prioritize objectives by considering factors such as statistical fit and causality. While statistical fit and causality are important in model development, there is no universally applicable approach that satisfies both criteria, as there is often a trade-off between them. Therefore, model estimation requires careful consideration of this trade-off (Washington et al., 2020).

Data-driven methods, such as machine learning (ML), including support vector machines (SVMs) (Corazza et al., 2010; Xu et al., 2012; Mokhtarimousavi et al., 2019) and neural networks (NNs) (Delen et al.,

E-mail address: zeke.ahern@hdr.qut.edu.au (Z. Ahern).

https://doi.org/10.1016/j.aap.2024.107690

∗ Corresponding author.

2006; Cansiz, 2010; Silva et al., 2021), offer an alternative to theorydriven approaches by capturing relationships between potential contributing factors and observed outcomes without relying on specific statistical modeling specifications. These methods are renowned for their prediction accuracy (García et al., 2010; Nassiri et al., 2014). However, their primary focus is on predictive accuracy rather than interpretability; they encounter challenges revealing causality and associated inferences. In the context of developing strategies to reduce crashes and improve road safety, the interpretability of models becomes crucial for understanding the underlying factors and relationships that contribute to crash occurrences (Mannering et al., 2020). SVM models have illustrated superior performance in effectively and accurately predicting crashes compared to traditional statistical models (Li et al., 2008). They can even offer comparable or better performance than neural networks while being faster to implement. However, a drawback of SVMs is their lack of interpretability, making it challenging to understand the underlying reasons behind their predictions. This is why SVMs and other machine learning models are often referred to as ''black boxes''. While SVM and NN models are primarily recommended for studies solely focused on predicting motor vehicle crashes when interpretability is not considered a priority, there are alternative ML approaches with greater potential for interpretability.

One approach that provides interpretability and feature selection is random forrest (rf), which utilizes importance scores to identify relevant features (Breiman, 2001; Li et al., 2019; Al-Moqri et al., 2020; Das et al., 2018). However, in count regression problems, categorical variables complicate feature extraction and understanding the joint effects of multiple variables. Another alternative, offering interpretability, is Bayesian networks, which explicitly model probabilities, dependencies and correlations between variables (Zou and Yue, 2017). However, constructing Bayesian networks for crash count prediction requires substantial data and prior knowledge, and their predictive accuracy may be limited as crashes involve various factors that do not strictly conform to an 'if-then' structure (Lin et al., 2015; Ma et al., 2018). Moreover, Bayesian networks may not be well-suited for road analysis and can face stability issues with small sample size and low mean values in crash data (Song et al., 2022; Xie et al., 2007).

In contrast, statistical methods are often preferred in practice due to their ability to handle limited data, address challenges such as under-reporting and false records, and establish causality (Lord and Mannering, 2010). Furthermore, ML methods, which are fundamentally learning models, have their own susceptibilities. These include overfitting, sensitivity to modeling parameters, and vulnerability to outliers. R1cParticularly, these methods may overfit to the commonality of zero crashes, leading to incorrect and inconsistent estimates (Li et al., 2016; Silva et al., 2021). A major limitation of machine learning techniques in crash data analysis is their inability to reveal causal inferences, which are crucial for understanding the underlying factors contributing to crashes (Ding et al., 2018). Therefore, statistical models remain steadfast, considering their ability to provide meaningful and intuitive insights.

Within the realm of statistical methods, Poisson regression is a frequently employed approach for analyzing crash data. This popularity is primarily due to the tendency of crashes to follow a Poisson process under homogeneous conditions (Lord et al., 2005). However, complications arise when the data exhibits an excessive number of zero observations or when the observations are not evenly distributed — the counts are over-dispersed. In such scenarios, employing a Negative Binomial (NB) or Generalized Poisson (GP) model is typically considered sufficient for handling over-dispersion (Yu et al., 2019; Ash et al., 2019; Ackaah and Salifu, 2011) . However, not all road conditions necessarily conform to the assumptions of the NB or GP model. Instances of under-dispersion in road crashes, where the observed variation is less than what would be expected under a Poisson distribution, can occur. In such cases, the Conway-Maxwell Poisson model could be a suitable alternative (Simmachan et al., 2022). However, implementing the Conway-Maxwell Poisson model often necessitates Bayesian methods due to the lack of a closed-form solution, leading to less frequent use given the rarity of under-dispersion in crash data and long computation times. Despite their complexity, Bayesian models are reserved for instances without closed-form options, as outlined in Ma et al. (2008) and Lord and Mannering (2010). The NB model is favored for over-dispersed data because of its closed-form and ability to handle extra variation through gamma, exponential, and Poisson-Lognormal mixtures (Lord et al., 2005; Ash et al., 2020; Ma et al., 2008). These mixtures, however, often lack closed-form distributions. Bayesian methods facilitate the use of non-closed form functions, but minimal differences are noted between maximum simulated loglikelihood (MSL) estimation and Bayesian approaches for available closed-form models (Farrell and Ludwig, 2008). This highlights the importance of prioritizing computational efficiency by using closedform functions, as the MSL approach offers significant computational advantages, enabling the testing of a greater number of alternative hypotheses (Greene, 2007). This computational efficiency leads to superior results, as it enables the consideration of many alternative and meaningful hypotheses, facilitating formal, spatial, and mechanistic modeling decisions crucial for prediction and interpretation (Fay and Punt, 2006).

Alternative statistical methods, such as advanced heterogeneity models, excel in dealing with unobserved heterogeneity in crash data. One example is latent class Tobit regression models, which censors the data by assuming a limit on the number of crashes (Anderson et al., 2017). A prevailing and more conventional approach to address unobserved heterogeneity involves extending probabilistic models by incorporating random parameters (Moomen et al., 2020). This extension allows the distributions to be estimated as functions of observed contributing factors, thereby offering enhanced flexibility and sensitivity in capturing crash conditions (Behara et al., 2021; Behnood and Mannering, 2017; Seraneeprakarn et al., 2017; Venkataraman et al., 2014). This approach provides a more flexible and accurate representation of the data compared to models with fixed or average terms (Alnawmasi and Mannering, 2022; Anastasopoulos and Mannering, 2009; Mannering et al., 2016; Washington et al., 2020; Behara et al., 2021). Several studies have proposed advanced variations of random parameter models to effectively account for unobserved heterogeneity and capture correlation between different sources of heterogeneity (Venkataraman et al., 2011b; Huo et al., 2020). Given the multifaceted interplay of factors, acknowledging potential correlation structures in unobserved heterogeneity from random parameters is key. Ignoring these could undervalue parameter variances, especially when significant correlation exists (Saeed et al., 2019). To overcome this challenge, reviercorcorrelated random-parameter count models have been developed (Hou et al., 2018; Huo et al., 2020; Saeed et al., 2019; Venkataraman et al., 2011b; Ahmed et al., 2021, 2022). These models offer a framework to capture complex interactions and provide accurate estimates of model parameters. For instance, Yu et al. (2015) illustrated the superiority of the correlated random parameters approach in capturing the complex interactions between weather and geometric road conditions, which were shown to be interconnected and have related effects on crashes. Utilizing such approaches helps ensure accurate estimation of coefficients and better understanding of the amplified effects of these factors on crash outcomes, and thus ensures that the coefficients of the models are identifiable (Hou et al., 2018). Additionally, incorporating spatial correlation through random parameter models has been suggested (Quddus, 2008; Li et al., 2018) because it provides insight into how different road factors can influence each other. Although the use of correlated random parameters is generally considered to provide extra flexibility in capturing the different sources of unobserved heterogeneity, the approach is largely data-dependent (Fountas et al., 2018; Tang et al., 2020). For example, Saeed et al. (2019) found that both uncorrelated and correlated random parameters models produced comparable results. Consequently, it may be unnecessary to use correlated random parameters because they might introduce insignificant or mis-specified aspects into the model. For instance, if the road segments of interest share very similar explanatory characteristics, employing complex modeling structures involving random parameters and correlation could potentially lead to overfitting. This is due to the inclusion of additional model parameters for estimation, making the models overly complex. Therefore, complex models may not always be the appropriate solution (Hausman, 2001; Hutcheon et al., 2010).

While random parameter models have been extensively utilized in road safety studies, they are limited by their inherent parametric assumptions (Mannering and Bhat, 2014). The typical practice is to employ a Normal distribution due to its simplicity (Haque and Washington, 2014). However, there exist more unique distributions that are better suited for crash data analysis (Ghitany et al., 2008; Geedipally et al., 2012). These parametric assumptions may restrict the flexibility of random parameter models and their variants (Mannering and Bhat, 2014; Beeramoole et al., 2023; Ahern et al., 2024). Therefore, exploring alternative priors beyond the commonly used Normal distribution becomes necessary to capture the unobserved heterogeneity present in the data (Heydari, 2018). To adequately represent this heterogeneity, testing various distributional assumptions may be required , including specialized distributions such as Lindley, Gamma, and Triangular, among others (Vellaisamy and Upadhye, 2009).

Even when the estimated model is well-fitted, extensive predictive checks are imperative to assess generalization (Fouskakis and Draper, 2008; Cunto and Saccomanno, 2008). However, due to time constraints or limited understanding of data partitioning techniques, these predictive checks can often be neglected. Quantitative predictive checks play a vital role in verifying the transferability and broad applicability of the estimated model (Persaud, 2001; Washington et al., 2005). In selecting a predictive assessment metric, root-mean squared error (RMSE) and mean absolute error (MAE) are notably effective for count data, outperforming metrics such as mean absolute prediction error, particularly when the data has a high density of zero values (Hadji Hosseinlou et al., 2018). Both RMSE and MAE offer advantages by distilling model performance into a single quantifiable metric (Abdel-Aty and Haleem, 2011). Furthermore, they are relatively easy to implement and interpret, and address overfitting (Lu et al., 2014; Gomes et al., 2019). Therefore, the adoption of a robust predictive assessment metric is paramount to ensure model credibility, as neglecting adequate predictive checks can yield biased results (Rashidi et al., 2022).

2. Research gap and proposed contribution

Despite extensive research and development to incorporate the above model development efforts, there are inherent limitations stemming from the associated demanding task (Hauer, 2015). Analysts must navigate a systematic approach, involving multiple stages such as data collection and processing, model selection and fitting, evaluation, and implementation. Each stage necessitates numerous decisions, including hypotheses about likely contributing factors, accounting for unobserved heterogeneity, determining distributional assumptions, data transformations, statistical model assumptions, and implementing techniques for handling excessive zeros, among other considerations. Moreover, ensuring model accuracy is frequently neglected, leaving insufficient evidence about performance beyond estimation (Chai and Draxler, 2014). As a result, there is potential risk regarding interpretability and generalizability (Sifringer et al., 2020). Model development faces limitations imposed by available resources, human bias, and insufficient exploration, thus potentially missing out on valuable extensive hypothesis testing. These limitations are not unique to count data and road safety data; similar difficulties have been investigated in the context of discrete outcome models and cluster analysis (Beeramoole et al., 2023; Arteaga et al., 2022; Rodrigues et al., 2020; Paz et al., 2019). testR1In such cases, optimization methods have enabled the generation and simultaneous testing of numerous hypotheses, resulting in the development of parsimonious models that effectively alleviate the aforementioned limitations. For example, Ahern et al. (2024) proposed a framework that considers simultaneously multiple modeling decisions but evaluates the specifications using only in-sample goodness-of-fit. There is no literature yet that provides a detailed investigation regarding trade-offs between model in-sample fitness and out-of-sample prediction accuracy, including how these performance measures affect the specification search while considering simultaneously multiple interdependent modeling decisions. Section 5 of this paper provides some critical insights into these trade-offs. Furthermore, Ahern et al. (2024) did not consider correlated random parameters, which could help address unobserved heterogeneity and relieve potential biased estimates.

Road safety analysis will also benefit by a framework that enables the simultaneous consideration of multiple important modeling aspects and associated hypotheses testing (Shin and Washington, 2012; Paz et al., 2019). Such a framework would assist analysts in overcoming challenges associated with ad hoc estimation, heterogeneity, distributional assumptions, and model prediction, leading to estimates with less chances of modeling errors relative to the current state-of-the-art.

The testing of diverse hypotheses, considering complex modeling characteristics, and providing an impartial exploration should be enabled by this framework (Wegman, 2017). testR1arHowever, such framework should not be intended to replace the analyst. Instead, it should aim to support the estimation process by alleviating subjective selection of hypotheses and enabling extensive strategic and unbiased search of the solution space. The experience and knowledge of the analyst remain crucial to ensure the correctness of road-safety interpretation and associated outcomes derived from the estimated models. The proposed framework is designed to serve as a decision-assistance tool, enhancing modeling and rigorous hypothesis testing to identify potential crash contributing factors and their interdependencies. Analysts are encouraged to remain actively involved in the estimation process, utilizing the framework to save time and potentially develop more accurate models than those achievable with current state-of-the-art methods (Beeramoole et al., 2023, 2024).

Contributions to literature are summarized as follows:

  • A multi-objective framework for the estimation of crash frequency models that generates and tests a broad range of meaningful hypotheses considering simultaneously potential likely contributing factors, non-linearities, correlated and non-correlated random parameters involving unique combinations of distributions, and transformations, all contained within a selected closed-form probabilistic model.
  • A mathematical programming formulation including a formal description of many crash count data modeling aspects. The formulation provides details involved in model estimation including selection of likely contributing factors, modeling them with uncorrelated or correlated random parameters, considering their distributional assumptions, and transforming them as needed. All these aspects are seamlessly considered for statistical selection of the most suitable probabilistic close-form, including its relevant assumptions, such as a dispersion parameter if applicable.
  • The formulation involves minimizing two objectives, an in sample goodness of fit metric and an out sample prediction metric, to seek balance between model simplicity and prediction accuracy.
  • A multi-objective search algorithm to optimize model performance while ensuring transparent record-keeping of Pareto-Efficient solutions. This allows analysts to select the most suitable solution using their knowledge and experience, offering flexibility and informed decision-making among multiple options.
  • The proposed optimization framework employs a diverse range of metaheuristic algorithms and alternative search strategies to handle unique data sets with complex and non-convex problem characteristics.

  • The framework enables capture of valuable insights that conventional methods often overlook. Through the simultaneous exploration of numerous hypotheses, it facilitates the discovery of significant insights that analysts might not have been able to identify due to limited time, knowledge, experience, or bias.

  • The proposed framework enables the investigation and estimation of specifications for a wide range of modeling applications. This includes scenarios involving many contributing factors, over and equally dispersed data sources, transformations, inter-factor correlations, and a diverse range of distributions to capture unobserved heterogeneity.

3. Methodology

This section presents a methodology that encompasses a mathematical programming formulation and associated metaheuristic solution algorithm.

3.1. Mathematical programming formulation

Consider the following simple fixed effects model for the expected number of crashes for observation (road segment) , as illustrated in Eq. (1). The term exp( ) is the model error term, where conforms to the assumptions of an analyst's selected model. For example, if the model conforms to Negative Binomial (NB), then the error follows a gamma distribution with shape parameter 1 and scale parameter . Therefore, in the formulation presented later, the error term is conditioned by the assumptions from a selected closed-form model from a set , which may include Poisson and NB models commonly used in the literature (Lord et al., 2005; Geedipally and Lord, 2010). This error term is appropriately conditioned in Section 3.4.

$$\lambda_i = \exp \left[ \beta^F X_i^F + \varsigma_i \right] \quad \forall i \in \mathcal{I}$$ (1)

In this model, a vector of likely contributing factors for observation is selected to form , with the coefficients for the fixed effects represented by . By utilizing the model in Eq. (1), access to the necessary information for analyzing and predicting the frequency of crashes in a specific road segment is obtained by the analyst. However, it should be noted that this is subject to appropriate selection of likely contributing factors to form . The widely adopted simple fixed effects model finds extensive usage in the field to extract valuable insights regarding patterns and trends concerning crash occurrences. Its popularity arises from its simplicity, making it easy to implement and interpret. However, it is essential to recognize that this simplicity may potentially oversimplify the underlying crash processes, limiting its ability to capture inherent heterogeneity as a function of the likely contributing factors. For instance, while the error term accounts for heterogeneity across segments, it does not consider any variability among likely contributing factors themselves. As a result, the model explanatory power might be constrained when attempting to capture nuanced interactions and complexities that could exist between different contributing factors. To address this limitation, many studies have utilized a traditional non-correlated random parameters (RP) or a more complex correlated random parameters (CRP) approach, which allows for correlations to exist among the random parameters (Huo et al., 2020; Hou et al., 2018). Despite the increased complexity of the CRP model, the same set of likely contributing factors are used in both RP and CRP approaches. Therefore, the simple fixed effects model in Eq. (1) can be extended by incorporating both RP and CRP components, as shown by generalized equation in Eq. (2).

$$\lambda_{i|\omega_i} = \exp\left[\beta^F X_i^F + \beta_i^R X_i^R + \zeta_i\right] \quad \forall i \in \mathcal{I}$$ (2)

where, Eq. (2) is conditioned by a vector of random parameters , which is × 1 in size, is the number of random parameters in the final model, enables observations to be subject to variations based on a predetermined distribution; it is commonly assumed that all random parameters follow a Normal distribution, as is often practiced by analysts (Greene, 2007). However, in Section 3.3, the utilization of alternative distributions is illustrated by introducing a set of random parameter distributions , which better capture and reflect the variations of . Furthermore, in Eq. (2), there exists , which represents the likely contributing factors utilized to incorporate random parameters. The coefficients of the means and random parameters, along with potential correlations between them, are represented by . However, to adequately account for these correlations, it is essential to construct appropriately (Saeed et al., 2019). The construction of with correlated random parameters is presented in Eq. (3).

$$\beta_i^R = \beta^R + C\omega_i \quad \forall i \in I$$ (3)

where, is the variance–covariance matrix (Cholesky matrix) which can capture correlations among random parameters. This Cholesky matrix can capture the covariance matrix ̂ such that . = ̂. The Cholesky matrix is assumed to be a lower-triangular matrix with the following structure when all random parameters have correlations to each other:

$$C = \begin{bmatrix} (\sigma{1,1})^2 \ \sigma{2,1}^2 & (\sigma{2,2})^2 \ \vdots & \vdots & \ddots \ \sigma{i-1,1}^2 & \sigma{i-1,2}^2 & \cdots & (\sigma{i-1,p-1})^2 \ \sigma{i,1}^2 & \sigma{i,2}^2 & \cdots & \sigma{i,p-1}^2 & (\sigma{i,p})^2 \end{bmatrix}$$

$$(4)$$

Here, is the Cholesky matrix of the random parameter vector It is a square matrix of size × . The diagonal elements of represent the variance of each individual random parameter, while the off-diagonal elements represent the covariance between each pair of random parameters. Non-zero off-diagonal elements in indicate that the random parameters are correlated with each other, which means that the values of one parameter can provide information about the values of another parameter. The correlation between random parameters can be positive or negative, according to the sign of corresponding off-diagonal element. In other words, captures the interdependence among the random parameters and is used to account for the potential correlation among them.

However, merely presenting a random parameters and a correlated random parameters model in isolation fails to acknowledge the potential benefits of their combined application. Specifically, assuming that all parameters are estimable with correlations, as illustrated in Eq. (4), can result in computationally intensive processes that may not correctly capture the underlying crash process. Hence, it becomes crucial to construct the Cholesky matrix appropriately, reflecting the most suitable model for accurate representation.

While may contain all random correlated parameters, some of those random parameters in need to be determined as strictly noncorrelated, which is defined as the set , and others as correlated, which is defined as the set , such that | ∪ | = and ∪ = . This construction is necessary, because a mixture of correlated and non-correlated random parameters can constitute the most suitable model structure, yielding enhanced outcomes particularly in coefficient estimation. The highest hierarchical structure that the proposed framework can consider is a fully correlated random parameter model, where = . This structure is what is being captured in Eq. (4) .

By adjusting the structure of with the assistance of decision variables, it is possible to create models with different correlation patterns among the random parameters. The specific values and arrangements in parameter matrix determine the characteristics and extent of the correlations. For example, if is an empty matrix (which would be imposed by decisions later in Section 3.4), the model simplifies to a fixed-effects model. This approach can help alleviate subjective hypotheses, considering the primary focus on achieving a balance between fit and complexity, effectively addressing the bias–variance trade-off (Kuha, 2004). The dimensionality of changes in relation to , and similarly, the dimensionality of changes in relation to . Therefore, the determination of and in Eq. (2) is required. Selecting and from the available data ,, along with any necessary transformations, is also a required modeling decision. Additionally, it may be desirable to consider reduced specifications that capture strictly random, correlated random, fixed, or combinations thereof. In this case, can be treated as a decision variable with a structure encompassing potential correlation and non-correlation elements. By adopting this approach, the most suitable hierarchical model for the given data can be identified.

The formulation of a mathematical programming problem aims to determine the most suitable values for , , , , and , enabling the effective utilization of the CRP model and reduced hierarchies while minimizing the objective functions. The mathematical approach for estimating crash-count data models is presented in Section 3.3, which involves evaluating different modeling decisions, such as the type of probabilistic model, relevant contributing factors and their transformations, and the potential for correlated or non-correlated random parameters and their associated parametric distributions. , , , and are decision variables, and thus calculated by solving the mathematical program. The notation used in this formulation is explained in Table 1.

3.2. Notation and definition

In Table 1, the mathematical notation is described, divided into multiple sections including: subscripts, sets, data, and binary and discrete decision variables (Structural). Further, continuous decision variables estimated through the lower-level maximum log-likelihood function, along with the subsequent logic derived from the implied coefficients are sectioned off in its own subheading (Decision Variables - Solver Dependant).

The estimation of crash count models is formulated as a multiobjective optimization mathematical program. The first objective is to minimize an in sample goodness of fit metric from , where = {AIC, BIC, HQIC}. The analyst can choose among these metrics considering the study objectives, data characteristics, and experience, as illustrated in Eqs. (5)–(8). As Sawalha and Sayed (2006) suggests, model generality should adhere to the principle of parsimony. This principle advocates for explaining as much data variability as possible with the least number of explanatory variables. The goal is to capture the actual underlying interdependence among dependent and independent variables, an objective that becomes increasingly attainable as the sample size grows. Therefore evaluations should focus on selection of models that best represent the data (Burnham and Anderson, 1998). There are divergent opinions about model comparison with some favoring complex models (Lin et al., 2014), while others prefer a balance between uncertainty, stemming from overly complex models, and bias, arising from overly simplistic models. According to Aho et al. (2014), a model should be parsimonious, as simple as possible while still effectively capturing the fundamental characteristics of the data under study. The model should be based on a set of parameters that the researcher identifies as meaningful and relevant to the system being studied. These parameters should include essential covariates, interactions, and higher-order terms, while maintaining minimal complexity.

In the realm of model selection, the Bayesian Information Criterion (BIC) has earned its reputation for consistently selecting the model that aligns best with the underlying data generation process (Vrieze, 2012), while the Aikike's Information Criterion (AIC) excels as a minimax-rate optimal tool for estimating a regression function (Yang, 2005). The BIC, in particular, is often the preferred criterion because of its roots in Bayesian approximations (Kuha, 2004; Beeramoole et al., 2023). It combines the likelihood of observing the data (expressed as in the formulation) with a penalty term linked to the model size or the number of parameters. As the number of parameters increases, the penalty term rises, actively discouraging overfitting and making the BIC adaptable to larger and more complex datasets. On the other hand, although the AIC has been a cornerstone in model selection for years, it may encounter difficulties in capturing the intricacies of complex models where data is extensive (Kuha, 2004). Different goodness of fit metrics, each presenting unique perspectives, are primarily designed to encapsulate the complexities of model selection. A balance between these criteria might provide an effective model selection approach. In this context, the Hanan-Quinn Criterion (HQIC), proposed by Hannan and Quinn (1979), serves as a potential blend that combines the advantages of both the AIC and BIC.

$$AIC = -2 \times ln(L) + 2\zeta \tag{5}$$

$$BIC = -2 \times ln(L) + \zeta \times ln(|\mathcal{I}^{in}|)$$ (6)

$$HQIC = -2 \times ln(L) + 2\zeta \times ln(ln(|\mathcal{I}^{in}|))$$ (7)

$$CAIC = -2 \times ln(L) + \zeta ln(|\mathcal{I}^{in}| + 1)$$ (8)

where, is the likelihood of the model estimation for all in-sample observations.

The second objective is a choice from , where = {MAE, RMSE}, which is calculated from the out-of sample observations to support model estimation, prevent overfitting and make accurate predictions (Xu et al., 2021; Hauer, 2015). These objectives are expressed in Eqs. (9) and (10), which was selected based on its frequent utilization in prior research, its suitability for quantifying a single measurable output, its ease of implementation, and its ability to mitigate risk of overfitting (Abdel-Aty and Haleem, 2011; Lu et al., 2014; Gomes et al., 2019). Furthermore, the incorporation of multiple objectives offers the potential to generate multiple final models, inheriting a trade-off between the two objectives. As a result, analysts need to prioritize between objectives, recognizing that one objective may hold greater importance than the other. This enables assessment of model performance and facilitates well-informed decision-making.

$$RMSE = \sqrt{\frac{1}{|I^{out}|}} \sum_{i \in I^{out}} (Y_i - y_i)^2$$ (9)

$$MAE = \frac{1}{|I^{out}|} \sum_{i \in I^{out}} |Y_i - y_i|$$

$$\tag{10}$$

3.3. Mathematical program - Objective function

minimize 1 (℧ ) subject to 2 (℧ ) ≤ 2 , = 1, 2, …,

minimize 2 (℧ ) subject to 1 (℧ ) ≤ 1 , = 1, 2, …,

where:

1 (℧ ) is the 1 objective function, which is defined based on a selected equation from Eqs. (5)–(8).

2 (℧ ) is the 2 objective function, which is defined based on a selected equation from Eqs. (9) and (10).

℧ encompasses all of the decision variables, as defined in Table 1, where index refers to a Pareto-efficient solution in .

1 and 2 are the upper bounds for and objectives, respectively for all Pareto-efficient solutions in .

This objective function uses the epsilon constraint method because it is simple and easy to implement, and it is able to generate a set of Pareto-efficient solutions, which cannot be improved in one objective without worsening at least one of the other objectives (Deb et al., 2002; Dutta and Das, 2019).

Table 1 Notation: Multi-objective extensive hypothesis testing for the estimation of crash count data models.

Sets, data and subscripts Subscript to denote road segment observation. Set of all observations, = 1, 2, … *, ,*…. , Set of in sample and out-of sample observations, = 1, 2, …*, ,…. Set of all potential contributing factors available from the data = 1,*… *, ,*… , . Note || depicts the cardinality of set Set of likely contributing factors, which are strictly deemed as fixed effects. Set of likely contributing factors which are non-correlated random parameters. Set of likely contributing factors which are correlated random parameters. Set of likely contributing factors which are random parameters such that = ∪ ★ Set of all likely contributing factors, including, fixed and random parameters such that ★ = ∪ 1 Set of in-sample objective measures that measure goodness of fit {AIC, BIC, HQIC, …} 2 Set of out of sample objective measure that measure model prediction {MAE, RMSE, … } Set of model types, where = {Poisson, NB, … } Set of model assumptions capturing the error terms imposed by the model, where = {0, (1, ), … } Set of distribution types, where = {Uniform, Normal, Triangular, LogNormal… }. Set of mathematical transformations, where = {(), (), (), (), (), (), ℎ()*,*… }. Set of Pareto-efficient solutions. , Road data for observation ∈ , including all potential contributing factors ∈ . Set of real numbers within the allowable range determined by the analyst. Where = { ∈ ∶ ≤ ≤ }, where a and b are the lower and upper bounds. Correlation threshold to exclude highly correlated fixed effects. Small coefficient, to ensure model estimation is not under-determined Decision variables - Structural , Data that has been transformed using transformations from for each data column ∈ ★ over all observations ∈ . are the selected likely contributing factors for the fixed effects, where , is the value at observation for element ∈ are the selected likely contributing factors for the random parameters, where , is the value at observation for element ∈ *,*′ Binary variables to denote if element *,*′ of is to be estimated, allowing for *,*′ , Indicator variable, selecting distribution in for random parameter ; ∀ ∈ , ∈ . Is a vector of the selected random parameters, to be used for the purpose of random sample draws. Where , ∼ ∈ subject to , = 1, ∀ ∈ , ∈ . Indicator variable taking value 1 if likely contributing factor has an associated fixed parameter; 0 otherwise. Indicator variable taking value 1 if likely contributing factor has an associated random parameter; 0 otherwise. Indicator variable taking value 1 if likely contributing factor has an associated correlated random parameter; 0 otherwise. , Indicator variable taking value 1 if likely contributing factor is transformed according to subscript in ; 0 otherwise Vector of indicator variables, taking value 1 if the crash model is from selected set ; = {1 , … , , … , }, ∀ ∈ . Number of parameters in the estimated crash model. Decision variables - Solver dependant Vector of coefficients for likely contributing factors with fixed effects; = {1 *,*… , | }. Note, this is dependant on Vector of coefficients for random effects within the means associated with likely contributing factors; = {1 *,*… , | }. Note, this is dependant on Cholesky matrix of estimable parameters where ∈ R| |×| ̂ Covariance matrix, where ̂ = *. ,*′ Estimable parameter in matrix , ∀ ∈ Global parameters for all coefficients in the final model that encompass , .

3.4. Mathematical program - Constraints

Expected number of counts by the model

These two objectives are dependent on a lower-level estimation to determine the values of MSL and the coefficients for the parameters in our model, denoted by , where denotes a global set of parameters that accommodate the varying levels of fixed, random, and correlated random parameters (ie. , ). The simulated log-likelihood function described in Eq. (11) is used as the estimation method, including an error term determined by which is selected from the set to account for different probabilities of closed-form count models, such as Poisson, NB or GP. For example if = 1 then = 0.

$$ln(L) = \sum_{i \in \mathcal{I}} (ln \left[ P(Y = y_i | \lambda_i, \kappa_m, \omega_i, \varsigma_i, \Theta) \right])$$ (11)

In this paper, the , representing the observed crashes, is predicted by Eq. (2) which includes multiple conditional requirements, including , , , , and . Forming involves two components: the structural decisions and the solver-dependent decisions. The structural decisions are formalized as decisions that are constrained by the mathematical program, while the solver-dependent decisions refer to the coefficients that depend on the MSL. Therefore, the structural decisions can be perceived as hypotheses that an analyst could develop . For instance, an analyst may explore the possibility of introducing additional variables into the equation as detailed in Hauer (2015) before subjecting it to the influence of MSL.

After a hypothesis has been defined by the structural decisions, the model can be deemed as strictly fixed (ie. ), in which case the dependent vector of random parameter distributional assumptions , along with , would be empty ( = ∅). However, when the is non-empty, the determination of , with respect to the distribution of each relevant element from and the corresponding distribution from the set becomes necessary.

Alternatively, when is non-empty, the inclusion of occurs, which corresponds to random parameter coefficients . However, the inclusion of correlated random parameters implies the presence of additional beyond the number of likely contributing factors in To handle this, the sets (set of non-correlated random parameters) and (set of strictly correlated random parameters) are utilized for assistance. The likely contributing factors that exhibit correlation among random parameters is denoted by which is represented by Eq. (3). can vary in structure, allowing for flexibility, and as a result, the associated can be constrained to explicitly represent different types of random parameter structures. These structures can include strictly non-correlated, strictly correlated, or a mixture of the two in terms of random parameter relationships. Therefore, Eq. (11) is heavily influenced by numerous decisions, particularly those involved in the construction of Eq. (2).

In Eq. (12), the selection of likely contributing factors , , and is made to ensure that each variable exhibits relevant characteristics. Consider this to be equivalent to the act of incorporating an additional variable into the model equation during the process of model development (Hauer, 2015). The value of indicates whether likely contributing factor is estimated as a fixed effect, while represents a random parameter effect that is uncorrelated, and represents a random parameter effect that is correlated with other correlated random parameters. These decision variables ( , , and ) determine which likely contributing factors are included in the final model, and identify the variables that are strictly fixed, strictly random, or strictly correlated random parameters. Additionally, the model allows for various data transformations, which can further enhance model estimation. Constraints on the data transformations, related to likely contributing factors, are formulated in Eq. (13). Similarly, in Eq. (15), relevant transformations are applied using indicator variables ,, which determine the appropriate function to be applied to data column ∈ , resulting in a transformed vector of observations for each . The set includes various transformations that can be used for likely contributing factor within ,. The final set of potential contributing factors from is considered for the new transformed dataset only if Eqs. (13) and (14) are satisfied. Note in Eq. (14), a transformation should be disregarded if it leads to values that fall outside the set of real numbers within an analysts desired range (

$$1 \ge \alpha_k + r_k + \iota_k \qquad \forall k \in \mathcal{K} \tag{12}$$

$$\sum{l \in \mathcal{T}} \tau{k,l} = \alpha_k + r_k + \iota_k \qquad \forall k \in \mathcal{K}$$ (13)

$$\tau{k,t} = 0 : t(X{i,k}) \notin \mathcal{R}_c \qquad \forall t \in \mathcal{T}, k \in \mathcal{K}$$ (14)

$$X{i,k}^{\dagger} = t \cdot (X{i,k}) \qquad \forall k \in \mathcal{K}, t \in \mathcal{T} : \tau_{k,t} = 1, i \in I$$ (15)

To ensure that the random parameters ( of length | ) are sampled from the correct distribution, certain constraints must be imposed. Specifically, should be distribution-specific for each likely contributing factor in . Managing the dimensionality of requires updating the random parameter set to include the dependencies of and . As a result, subsets of are defined for fixed parameters as , random parameters as , and correlated random parameters as . These updates are depicted in Eqs. (16)–(20), respectively.

Additionally, to handle the issue of numerous zero elements in the transformed data, variables and are utilized to extract only relevant transformed characteristics, as described in Eq. (2). This separation is achieved by reducing , into and as illustrated in Eqs. (21) and (22). Specifically, is obtained using Eq. (21), while is derived from Eq. (22).

$$\mathcal{K}_F = {k \in \mathcal{K} | \alpha_k = 1} \tag{16}$$

$$\mathcal{K}_R = {k \in \mathcal{K} | r_k = 1} \tag{17}$$

$$\mathcal{K}_C = {k \in \mathcal{K} | i_k = 1} \tag{18}$$

$$\mathcal{K}{P} = \mathcal{K}{R} \cup \mathcal{K}_{C} \tag{19}$$

$$\mathcal{K}_{\star} = \mathcal{K}_P \cup \mathcal{K}_F$$ (20)

Now that the sets are defined, the construction of ,, and , is achieved by using sets in Eqs. (21) and (22).

$$X{i,k}^F = X{i,k}^{\dagger}$$ $\forall i \in I, k \in \mathcal{K}_F$ (21)

$$X{i,k}^R = X{i,k}^{\dagger} \qquad \forall i \in I, k \in \mathcal{K}_P$$ (22)

With , appropriately defined, can be determined; distribution specific. In each observation, a draw is made from a unique distribution according to the random parameter ; this is enforced in Eq. (23). Similarly, the most appropriate error terms are enforced by Eq. (24), the application of the most appropriate distributions for each random parameter is achieved in Eq. (26). Random parameters are based on the parametric assumption of a distributional form contained within set . Thus, enforcement that at most one distribution can be selected is achieved by Eq. (25).

$$\omega{i,k} \sim d$$ : $\xi{k,d} = 1$ $\forall d \in D, k \in K_P, \forall i \in I$ (23)

$$\varsigma_i \sim m : \kappa_m = 1 \qquad \forall m \in \mathcal{M}^E, i \in \mathcal{I}$$ (24)

Furthermore, to accommodate random parameters, an indicator variable is introduced in Eq. (25) to select the appropriate distribution.

$$\sum{d \in D} \xi{k,d} = 1 \qquad \forall k \in \mathcal{K}_{\mathcal{P}}$$ (25)

$$\sum_{m \in \mathcal{M}} \kappa_m = 1 \tag{26}$$

Hence, through Eq. (12), the likely contributing factors can only be defined as fixed, random, or correlated random. For the fixed effect terms, highly correlated terms are removed, which is enforced in Eq. (27).

$$\alphak \times |X{i,k}^{\dagger}, X_{i,k'}^{\dagger}| \le c \forall k, k' \in \mathcal{K}$$ (27)

To ensure consistency in both correlated and non-correlated parameter settings, all diagonal elements are set to 1, as mandated by Eq. (28). Subsequently, Eq. (29) guarantees the estimability of lower triangular correlated random parameters. Finally, to exclude undesired estimation avoiding correlated random parameters, Eqs. (30) and (31) effectively disable and prevent estimation.

$$\gamma_{k,k} = 1 \quad \forall k \in \mathcal{K}_P$$ (28)

$$\gamma_{k,k'} = 1 \quad \forall k, k' \in \mathcal{K}_C : k' < k$$ (29)

$$\gamma_{k,k'} = 0 \quad \forall k, k' \in \mathcal{K}_R, : k' \neq k$$ (30)

$$\gamma_{k,k'} = 0 \quad \forall k, k' \in \mathcal{K}_C, : k' > k$$ (31)

The objective function is influenced by its total number of parameters, distinguished as . This is achieved in Eq. (32), which includes the parameters encountered by fixed, random, correlated and dispersion.

$$\zeta = \sum_{k \in \mathcal{K}_F} (\alphak) + \sum{k \in \mathcal{K}_R} (rk) + \sum{k \in \mathcal{K}_C} (\iotak) + \sum{k \in \mathcal{K}P} \sum{j \in \mathcal{K}P} \gamma{k,j} + (1 - \kappa_m) \forall m \in \mathcal{M}$$ (32)

4. Multi-objective extensive hypothesis testing solution algorithm

4.1. Motivation

The above multi-objective mathematical program exhibits nonlinearity, intricate hierarchical structures, and a large number of interdependent decision variables, necessitating an efficient solution algorithm. Due to the intricacy of the proposed estimation problem, multiple iterations may be required to find solutions that meet the study requirements, resulting in considerable computational time. For example, the analysis might involve considering total or fatal crash frequency. To enhance results, the problem can be solved multiple times with gained knowledge incorporated through constraints to test alternative hypotheses and objective functions. Enforcing contributing factors or imposing restrictions based on previous findings can lead to more reliable solutions. Consequently, leveraging past results and constraints can guide the design of future experiments and research endeavors.

Given that safety strategies are influenced by data characteristics, essential contributing factors and relevant exposure effects must be carefully considered in the final model. Finding an efficient solution algorithm is imperative to effectively address the complexity of the proposed multi-objective mathematical program and to achieve meaningful results that align with the intended goals.

One of the most promising approaches that satisfies the criteria of efficient model selection while maintaining interpretability was introduced by Veeramisti et al. (2020). They employed a Simulated Annealing (SA) algorithm to estimate clusterwise safety performance functions (SPFs), enabling simultaneous estimation of the optimal number of clusters and associated safety performance functions. However, their framework was limited to estimating SPFs using only exposure, neglecting other likely contributing factors, non-linearities, correlations, and

1 Equivalently ∑ ∈ ( + + ).

random parameters including correlations and their distributional assumptions. Despite these limitations, the use of a metaheuristic-search approach to determine an optimal number of clusters was effective.

In the context of combinatorial optimization approaches to support modeling decisions, Paz et al. (2019) proposed a mathematical programming formulation and a SA solution algorithm to identify mixed-Logit specifications and associated parameter distributions. Addressing the decision environment context, Beeramoole et al. (2023) proposed a Harmony Search (HS) algorithm to test hypotheses that encompass likely contributing factors, non-linearities, correlations, and random parameters. However, these problems are based on mixedlogit models, leaving a gap for a similar approach to address count data problems, even though they share similar characteristics. Research has consistently illustrated the effectiveness of metaheuristics in solving a wide spectrum of optimization problems. These methods typically required reasonable resources and time through strategic search techniques, especially when alternative approaches are ineffective or excessively complex (Chew et al., 2013; Koza, 1994; Chakroborty, 2003).

Among the available metaheuristic algorithms that can be employed to solve the proposed mathematical program, there are algorithms that have demonstrated suitability due to their unique and efficient search strategies, as well as their limited and easily-tuned hyperparameters (Oliveira et al., 2017). Notably, HS (Alia and Mandava, 2011; Beeramoole et al., 2023), SA (Kirkpatrick et al., 1983), and Differential Evolution (DE) (Dia and Panwai, 2009) have exhibited exceptional performance in solving discrete optimization problems, evading local optima and converging to efficient solutions that surpass local search and exact methods for non-linear and non-convex problems (Durillo et al., 2010; Ommen et al., 2014).

While metaheuristics cannot guarantee a global optimum, they can escape local optima and intensify the search to achieve convergence towards near-optimal solutions with reasonable resource requirements (Oliveira et al., 2017). Hence, three metaheuristic-search based approaches were devised and implemented to solve the proposed mathematical programming formulation. Each metaheuristic leverages a unique search strategy, inexorably seeking convergence towards a set of solutions that conform to the Pareto-frontier.

4.2. Solution algorithms

The three metaheuristic algorithms developed in this study and in Ahern et al. (2024) follow a discrete framework. Discrete DE and HS have been inspired by the works of Wu et al. (2018) and Kattan et al. (2010), respectively, and SA is naturally discrete. All algorithms have been adjusted to include a time-based stopping criterion in the initial testing phase to ensure fair comparison between algorithms. For the final proposed models in Sections 5.4–5.6, an alternative stopping criterion is utilized, which is based on the number of iterations without a Pareto-efficient improvement. This criterion, denoted as , is employed as referenced in Ahmadian et al. (2021) and Badillo et al. (2013). See Table 6 for reference to the algorithms hyperparameters.

HS is a metaheuristic algorithm inspired by the improvisation process of jazz musicians, wherein initial harmonies are adjusted over time through pitch changes to achieve an optimal harmony (Geem et al., 2001). Three main hyperparameters define HS which include: Harmony Memory Size (), Pitch Adjustment Rate (), and Harmony Memory Considering Rate (). In the discrete implementation of HS, the Pitch Adjustment Index () governs the permissible displacement of an ordered sequence of decisions when activated by the pitch adjustment operator. The Harmony Memory set (), is initially filled up to a size equal to . Each Harmony Memory embodies a prospective hypothesis specification for the mathematical program. In this context, a specification refers to a crash count model that can be analyzed, and includes goodness of fit metrics and predictive accuracy information. A single element of , represented as Harmony , signifies a solution to the mathematical program. The quality of the Harmony is evaluated based on the multiple objectives Eqs. (9) and (10).

In our implementation of Multi-objective HS, a dynamic is utilized which enables the to be equal to Pareto-efficient set and avoid sensitivity issues while determining an optimal size. Dynamic, because we allow to be equal to if all solutions within are mutually non-dominated and the size of the Paretoefficient set is greater than the initial population size. Nonetheless, remains a hyperparameter, as we enforce that the size of cannot fall below its initial population size. Thus, the size of || will invariably be at least equal to but is free to exceed this size, if all solutions are Pareto efficient. This dynamic modification mechanism ensures a diverse range of solutions, even in scenarios where || = 1.

To sort the set, is arranges non-dominated solutions first, followed by those that tend towards balancing the non-dominated, thereby maintaining a controlled population size. This arrangement strategy draws inspiration from the Strength Based Pareto Evolutionary Algorithm 2 (SPEA2) (Ordoñez et al., 2018).

Throughout the algorithm, at least two sets are continually monitored including the Pareto-efficient , which will be analyzed at the end, and the , to ensure exploration. If the Pareto-efficient set is larger than the input harmony memory size, these sets are equivalent.

A non-dominated strength Pareto weighted trailing approach is utilized to derive the generated harmony set, with new solutions improvised based on this set. This developed HS algorithm is illustrated in Algorithm 1.

SA incorporates hyperparameters such as an initial temperature (0 ), a cooling factor ( ), temperature steps ( ) The initial temperature 0 is calculated with a heuristic to start the algorithm with a desired initial acceptation rate (), thus reducing hyperparameter sensitivity, (see Johnson et al. (1989)). 0 is determined with = 0*.*5 so that approximately half of the candidate solutions are accepted at the begging of the search. This approach is based on the variance of the objective value for a random population of initial solutions, denoted as ̄. It has been found to be less sensitive to the results compared to defining 0 in other studies (Simon, 2013; Ahern et al., 2022; Baxendale et al., 2021; Ahern et al., 2024).

The initial temperature is calculated according to Eq. (33). In the multi-objective version of SA, each objective is assigned an energy. Therefore, for analysis, two energies, specifically 1 and 2 , are considered. If the proposed solution dominates the incumbent, new solutions are accepted. Otherwise, the energies of the proposed solution are compared based on their sum.

$$T_0 = \frac{\Delta_1 \bar{E}_1 + \Delta_2 \bar{E}_2}{\log(\chi)}$$ (33)

The SA algorithm is illustrated in Algorithm 2. The multi-objective SA exhibits a similar capability to track all Pareto-efficient solutions within set , when compared to the multi-objective HS approach. As a result, at termination, all solutions can still be effectively observed, demonstrating the trade-off between the two objectives.

Finally, discrete DE bears similarity to HS, except that it generates solutions based on three subsets of the population rather than the entire solution space. Generating new solutions from this mechanism and possessing one less hyperparameter than HS. This reduction in hyperparameters can be advantageous for future problems, provided the control mechanism is distinctive enough to effectively explore the entire search space, thereby minimizing unnecessary complexity. The main parameters that require tuning in discrete DE are the crossover rate () and the population size ().

A mutation function, represented by Eq. (34), depends on three solutions in the current population. It governs the mutations of elements in our decision variables, where represents a selected decision from ℧. Variable , represents the chosen decision for solution in

Algorithm 1 Multi-Objective Harmony Search

1: function MO-HS(, ,  ,  )
2: HM ← Generate_Random_Specifications(HMS) a
3: HM ← SPEA2_Sort(HM) b
4:  ← NonDominatedSort(HM) c ⊳ Track Pareto Set
5: while current_runtime() ≤   or termination(WIC) do d
6: for  = 1, 2, ..,  do
7: if (0, 1) e≤  then
8: r ← rand_int(1, HMS)f ⊳ Select random specification
9: H[j] ← HM[r, j]
10: else
11: H[j] ← rand_decision(j)g ⊳ Generate random
  specification
12: end if
13: if (0, 1) ≥   then
14: [] ←Conduct_pitch_adjustment() h
15: end if
16: end for
17: 1
         (℧), 2
             (℧) of H ⊳ Evaluate Objectives Eqs. (9) and (10)
18: if 1
          (℧), 2
               (℧) ≤ 1
                     (℧), 2
                         (℧) of ([]) then ⊳
  Compare the solution to worst in HM.
19: [] ←  ⊳ Replace the worst solution.
20:  ← SPEA2_Sort()
21:  ← NonDominatedSort( ∪ )
22: end if
23: end while
24: end function
  • a Generates a population of random specifications to the size of the input.
  • b Sort by the SPEA2 algorithm.
  • c Breaks down to the Pareto Efficient Set, such that each solution is non-dominating to each other.
  • d Find the number of iterations without a Pareto-Efficient improvement, terminate on threshold parameter .
  • e Generates a random continuous value uniformly distributed between 0 and 1.
  • f Generates a random integer, from 1 to HMS.
  • g For element j of the chromosome, generate any possible decision randomly available.
  • h Discrete Pitch Adjustment.

the population for element and ,, correspond to current indexed element for a solutions at indexed element , such that the amount of decisions are in a in order. An equation for , utilizes the modulo argument to ensure that new solutions can be generated within the available decision space, without exceeding its bounds. The equation is presented in Eq. (34).

$$v{p,j} \leftarrow x{p,j}[(\ell{p,j} + (\ell{r1,j} - \ell{r2,j})) \mod (|v{p,j}|)]$$ (34)

where, is the index of the current solution observed from the population, and where 1 and 2 are distinct, randomly chosen solutions from the population.

The DE, illustrated in Algorithm 3, and HS algorithms utilize the same population control with Pareto-efficient tracking. In our implementation of DE, the population set is defined as because it corresponds to the same population based-control method which was implemented in the HS algorithm.

4.3. Repair method

In each iteration of all tested metaheuristic algorithms, the solution representation may require modification to adhere to the constraints of the mathematical program. This process entails generating and

Algorithm 2 Multi-Objective Simulated Annealing Algorithm

1: function MOSA(0
                ,  , , , 
2: 0 = {1
          , …, 
                } ← Generate_Random_Specifications(
3:  ←  (, 0
              a ⊳ Eq. (33)
4:  ← NonDominatedSort(0
5:  ← SPEA2_Sort()[1] ⊳ Gets the best solution from the Pareto
6: 1
      , 2 ← 1
             (℧) and 2
                    (℧) of  ⊳ Calculate the two objectives
7: while current_runtime() < MaxTime or termination(WIC) do
8: for q = 1,2, … , step do
9: 
          ′ ← randomly chosen neighboring state of 
10: 
           ′
           , ′
             2 ← 1
                  (℧) and 2
                         (℧) of 
                              ′
11: if 
            ′
             < 1 or 
                    ′
                     < 2
                         then
12:  ← 
13:  ← NonDominatedSort( ∪ )
14: 1
             , 2 ← 
                   ′
                   1
                    , ′
15: else
16:  ← 
               −(
                  −1+
                       −2
                         )∕
17: if (0, 1) <  then
18:  ← 
19: 1
                , 2 ← 
                      ′
                      , ′
20: end if
21: end if
22: end for
23:  ←  × 
24: end while
25: end function

Algorithm 3 Multi-Objective Differential Evolution

1: function MO-DE(,  ,  )
2: HM ← RandomiseModelSpecifications( )
3: HM ← SPEA2_Sort(HM)
4:  ← NonDominatedSort(HM)
5: while current_runtime ≤   do
6: for  = 1, 2, …,  do
7:  ← []
8: 1 ← rand_int(1, ) ∶ 1 ≠ 
9: 2 ← rand_int(1, ) ∶ 2 ≠ {, 1
10: for each dimension  ∈  do
11: , ← ,[(, +  (1
                           , − 2
                                ,)) mod (|,|)] ⊳
  Offspring created as per Eq. (34)
12: if (0, 1) <  then
13:  ←  ⊳ Crossover Occurs on element j
14: else
15:  ←  ⊳ Element j remains the same for
  solution 
16: end if
17: end for
18: end for
19: for  = 1, 2, … ,  do
20: if 1
            (℧), 2
                (℧) of  ≤ 1
                         (℧), 2
                             (℧) of ([]) then ⊳
  Compare solution to worst in HM
21: [] ← 
22:  ←SPEA2_Sort() ⊳ update population
23:  ← NonDominatedSort( ∪ 
24: end if
25: end for
26: end while
27: end function

a Function to calculate initial temperature, takes population of solutions, and the desired acceptance ratio .

Fig. 1. Hypothesis testing repair procedure.

conducting the appropriate significance level testing for potential contributing factors. Likely contributing factors are those associated with a rejected null hypothesis. Similarly, specific factors within its proposed hierarchical structure associated with a null hypothesis that was not rejected are removed from the specification. This iterative process persists until the contributing factors are determined to be significant. Non-significant factors are retained only if their exclusion would result in the removal of significant ones, particularly when significance is present in the random parameters. A summary of this procedure is provided in Fig. 1.

4.4. Analyst input

To assist an analyst, the proposed framework takes into account all decisions and constraints in the formulated mathematical program. The analyst can select objective functions and constraints using knowledge and experience as well as insights gained during model development.

The flowchart in Fig. 2 visually depicts the integration of the metaheuristic solution algorithms within the framework.

5. Numerical experiments

The proposed framework and solution algorithms were evaluated by performing numerical tests to assess performance. Three different datasets were utilized to showcase various aspects of the solution algorithm, including its ability to:

  • 1. Estimate crash frequency models considering multiple objectives including minimization of the BIC and RMSE;
  • 2. Simultaneously test multiple hypotheses involving random parameters (correlated and non-correlated), associated distributions, data transformations, probabilistic model type;
  • 3. Provide insights regarding likely crash contributing factors; and
  • 4. Show how correlated random parameters may yield more reliable estimates, ensuring better statistical fit and support the prediction of out of sample data.

5.1. Data

In order to conduct benchmarking analyses, two datasets previously analyzed in literature, as well as one synthetic dataset, were used. The first real dataset, as analyzed in Behara et al. (2021), was obtained from the Queensland Department of Transport and Main Roads (TMR). This dataset encompasses 1848 unique roadway segments in Queensland, focusing on factors that could potentially contribute to head-on collisions involving Heavy Vehicles. It comprises 67 potential contributing factors, including, but not limited to, annual average daily travel (AADT), multi-combination vehicular traffic (MCV), number of lanes (Nlanes), and Roadtype is motorway (M). For a comprehensive description of the data, refer to Behara et al. (2021). This experiment was primarily aimed at examining the multi objective benefits of the proposed framework by exhaustively considering a Pareto-frontier of solutions.

The second dataset, derived from Washington et al. (2010), was utilized to evaluate the significance of correlated random parameters. Spanning from 1990 to 1995, this dataset incorporates information from 275 roadway segments in Washington State. It includes the predictor variable, representing the number of crashes, along with 30 potential contributing factors.2

The synthetic dataset was specifically crafted to evaluate the ability of the proposed framework to capture a known data generation process with relative complex characteristics including correlated random parameters, which followed a Poisson process. Complete knowledge of the data generation process is only possible using synthetic data. The synthetic dataset was generated comprising of 2000 observations, encapsulating both probable and potential contributing factors. The latter were generated but not linked to the count generative process. All three probable contributing factors were composed based on a correlated random parameters model, utilizing a normal distribution in every case. The creation of this synthetic dataset involved conducting a Monte Carlo simulation, adhering to a methodology akin to the one proposed by Hilbe (2010) and in our previous study (Ahern et al., 2024).

The model generation proceeded as follows:

  • 1. Incorporate three variables (X1, X2, and X3) as independent variables in the model, along with five other non-significant variables (X4-X8) intended to introduce random noise. This ensures that the three key contributing factors are integrated into the model, which will then influence the count response, conforming to the correlated random parameters model. Factors X1 and X2 were modeled as normal random parameters, exhibiting mutual correlations, forming complete matrix of .
  • 2. Simulate crash occurrences using the coefficients and explanatory variables as outlined in Eq. (2). This method is designed to highlight the shared correlated random parameters between factors X1 and X2, following the normal distribution.
  • 3. Repeat steps 1 and 2, 100 times, where the coefficients from Eq. (2) remain the same.3 100 was selected because it has been considered sufficient for reliable estimates, to achieve credible count predictions when averaged (Lord, 2006; Roque and Cardoso, 2014)
  • 4. Compute the average counts from these repeated trials to form the 'true' counts, alongside an average of the explanatory variables.

5.2. Experimental details

For simulating the maximum likelihood estimate, 200 Halton draws ( = 200) were utilized as per (Halton, 1960). The choice of Halton draws was influenced by their proven superiority over pure

3 The parameters to generate this model are $$\beta{constant} = -2$$ , $\beta{x1} = 1$ , $\beta{x2} = -1$ , $\beta{x3} = -3$ C = $\begin{bmatrix} \sqrt{2} & \ -\frac{1}{2} & \sqrt{4} \ -\frac{1}{2} & -\frac{1}{2} \end{bmatrix}$

$$\beta_{x3} = -3, C = \begin{bmatrix} \frac{1}{\sqrt{7}} & \frac{1}{\sqrt{7}} \ \frac{3}{\sqrt{14}} & \frac{14}{\sqrt{14}} & \frac{14}{\sqrt{1}} \end{bmatrix}$$

2 This problem has also been analyzed by Zhan et al. (2015) and Milton et al. (2008).

Fig. 2. Multi-objective extensive hypothesis testing framework.

random draws, and the number 200 has been deemed sufficient in various studies (Bhat, 2003; Train, 2009, 2016). Consequently, with a simulated maximum likelihood estimate, the estimation of Eq. (11) will be simulated to assume the form represented in Eq. (35).

$$ln(L) = \sum{i \in I} ln \left[ \frac{1}{H{draws}} \sum{h=1}^{H{draws}} P(Y_i = y_i | \lambda_i, \omega_i, \zeta_i, \Theta) \right]$$ (35)

To demonstrate the ability of the proposed framework, Python software implementing the extensive hypothesis testing was developed. The software facilitates the efficient discovery of likely contributing factors and model structure through the proposed extensive hypothesis testing. The software, implemented in Python, is available as an installable package and has been used in prior research (Ahern et al., 2023, 2024). All computational experiments were conducted using this software, with each experiment assigned to a single thread on an SGI Altix XE Cluster. The analyst has the opportunity to define the set of decisions to consider, provide initial solutions, select variables to include in or exclude from the model (if any), and define the objective functions. The following sets of decisions were included in the software to perform the experiments in this paper:

  • = {Poisson, NB}
  • = {0, (1, )} (1, )
  • • = {Uniform, Normal, Triangular, Tn_Normal, Ln_Normal}
  • = {(), (), (), (), ()}
  • = , where 70% of the data was used for in-sample goodness of fit
  • = , where 30% of the data was used for out-of sample prediction

Set is defined with commonly used count data models, Poisson or NB,4 Subsequently, for each model, the error term captured in the data varies depending on the distribution used. Therefore, it is essential to account for these errors appropriately using . This term helps capture the model error for each mean equation . Assuming a

4 For a complete formulation of Poisson and NB see Mannering et al. (2016).

$$\begin{aligned} \text{Given: } \kappa_{\textit{Poisson}} &= 1, & \text{then,} & Pr(Y_i = y_i | \lambda_i, \text{Poisson}) &= \frac{\exp\left(-\lambda_i\right) \lambda_i^{y_i}}{yi!} \ \text{Given: } \kappa{\textit{NB}} &= 1, & \text{then,} & Pr(Y_i = y_i | \lambda_i, \text{NB}) &= \left(\frac{1/\nu}{1/\nu + \lambda_i}\right)^{1/\nu} \left(\frac{\Gamma(1/\nu + y_i)}{\Gamma(1/\nu) + y_i!}\right) \left(\frac{\lambda_i}{1/\nu + \lambda_i}\right)^{y_i} \end{aligned}$$

Box I.

Table 2 Hyperparameters for each algorithm.

| Algorithm | Hyperparameters covering array | | | | |------------------------|------------------------------------------|-----------------------------------------------------|---------------------------------|--| | Harmony search | 𝐻𝑀𝑆 | 𝐻𝑀𝐶𝑅 | 𝑃 𝐴𝑅 | | | | {20, 25, 30, 35, 40, 45, 50} | {0.5, 0.65, 0.8} | {0.3, .5, 0.7} | | | Simulated annealing | 𝑇𝑠 | 𝑇𝛼 | 𝐶𝑅 | | | | {2, 5, 10} | {0.95, 0.975, 0.99} | {0.05, 12.5, 20.27.5, 35, 42.5} | | | Differential evolution | 𝑃 𝑜𝑝𝑆𝑖𝑧𝑒
{20, 25, 30, 35, 40, 45, 50} | 𝐶𝑅
{0.1, 0.1625, 0.2875, 0.35, .4125, .5275, .6} | | |

dispersion parameter which captures the assumptions of ∈ , the models can be generalized as the following from the equations in Box I.

5.3. Covering arrays

Achieving efficient performance requires the essential task of tuning the hyperparameters in either HS, SA, or DE, as the effectiveness of metaheuristics often depends on the data used. Covering arrays were employed to assess the sensitivity of the metaheuristic hyperparameters and ensure their proper tuning for each algorithm (Izquierdo-Marquez et al., 2018). The selection of hyperparameters for the following experiments was based on the guidelines presented in Simon (2013). The specific hyperparameters considered for analysis are outlined in Table 2. Furthermore, during the initial testing phase, in order to meet time-restricted modeling requirements and achieve a sufficient range of hypotheses, the value of was configured to 10,000 s of CPU time, Ahmadian et al. (2021). Later, in Sections 5.4 to 5.6 when the final models are proposed, the time limit was not restricted and was replaced with the termination criterion of , which facilitates the intended exploration of the selected algorithm.

Hyperparameters were selected and evaluated against all datasets using a weighted approach. A weighted approach was used to facilitate a single measurable output from each algorithm, facilitating a direct comparative analysis of the algorithms, on all datasets. The outcomes of this analysis are detailed in Fig. 3.

Initial results using the three datasets indicated that the chosen metaheuristics are capable of solving the proposed extensive hypothesis testing problem. These results were obtained using a set of covering arrays depicted in Fig. 3. Performance analysis indicated that HS provides superior convergence. This results aligns with existing literature, where HS dominated 16 alternative multi-objective algorithms (Cobos et al., 2020; Beeramoole et al., 2023).

HS, while showing less sensitivity, still benefited from careful hyperparameter tuning. Notably, an increased Harmony memory and pitch adjustment rate bolstered performance across all datasets. Harmony memory size, specifically, showed minimal drawbacks, with a size of 30 appearing to be a viable choice. Both performance evaluation and hyperparameter analyses underscored Harmony Search as the best option across all three datasets. A heatmap analysis showcases the identification of superior regions in Fig. 3.

SA likely underperformed, primarily due to its inability to track a population of solutions, which in turn constrained its capacity to effectively explore both objectives. Frequently chosen for single-objective problems, Simulated Annealing may encounter performance difficulties due to hyperparameter sensitivity, struggling to identify the entire Pareto-Frontier set of solutions. Its propensity to minimize over a specific weighted metric among the two tested objectives might limit its performance. In contrast, DE possibly falls short due to an imbalance between exploration and exploitation stemming from its hyperparameters, leading to solutions becoming trapped in local optima.

In multi-objective optimization, the Pareto-efficient set necessitates an algorithm specifically designed for multi-objective optimization, such as HS. Future research should consider exploring alternative optimization algorithms and further validating the effectiveness of Harmony Search across other datasets.

5.4. Experiment using data from Queensland

After conducting an extensive hypothesis testing using data from Queensland, a Pareto-frontier was obtained, including mainly uncorrelated and correlated random parameter solutions. Only solutions that are in close proximity to the Pareto-frontier are highlighted and can be observed in Fig. 4. From Fig. 4, a solution was chosen for analysis, and the results are displayed in Table 3. 5 Note, certain variables may not be statistically significant on their own. However, they have been included in the final model due to the significance of their counterparts with random or correlated random parameters. These variables cannot be dropped from the model without affecting their significance.

Choosing a solution from the array of possibilities presented by the frontier necessitates an analyst who leverages knowledge and experience for decision-making. The framework proposed thus simplifies the modeling process, requiring minimal initial input and aiming to gather as much resourceful information as possible. Acting on behalf of an analyst, the chosen model was selected because of its balance between in-sample fit and out-of-sample prediction (Kuha, 2004). This was supported by considering both RMSE and BIC among all Paretoefficient solutions, thereby providing excellent balance between fit and prediction.

This model aligns closely with the prevailing literature and is based on reliable parameter estimates, including a commonly used offset expression, with AADT and road length, which together vehicle density (Markon and Krueger, 2004; Ackaah and Salifu, 2011; Polus and Cohen, 2012; Agbelie, 2016). By integrating vehicle density in to the model, a consistent offset term is ensured for all proposed models. This approach ensures that each hypothesized solution, with the offset

5 This is only one solutions reported from the extensive hypothesis search, for alternative proposed models see https://github.com/zahern/data/tree/ main/Results.

Fig. 3. Initial hyperparameter testing: objectives are weighted in this experiment. Each combination represents 5 seeded trials with 10,000 seconds of CPU time.

Fig. 4. Minimizing BIC and RMSE via extensive hypothesis testing: Pareto-frontier.

term in mind, delivers additional valuable insights from the other contributing factors.

All solutions existing on the Pareto-frontier surpass the benchmark model in terms of both objectives. Some solutions, not previously

Table 3 Models compared: Benchmark model from Behara et al. (2021) compared to the extensive hypothesis testing framework.

NB model found by Behara et al. (2021) NB model found using the HS algorithm
Effect 𝜏 Coeff Prob z > Z Effect 𝜏 Coeff Prob z > Z
Constanta no −12.41 0.00*** Constant no −14.82 0.00***
AADT log 2.99 0.00*** AADT arcsinh 5.95 0.00***
MCV log 5.70 0.00*** MCV arcsinh 5.85 0.00***
Urban single carriageways (US) no 2.29 0.02* Shoulder width (SW) no −1.81 0.07.
Rural single and medium speed
(RSMS)
no 2.69 0.01* Speed limit in km/h (SP) sqrt 4.23 0.00***
Significant curvature (moderate to
very sharp which is longer than
50% of the segment) (Curve50)
no −1.15 0.06. Length of road in km (Length) arcsinh 1.85 0.00***
Rural single and high speed
(RSHS)
no 4.11 0.00*** Rural dual carriageway (RD) no −27.01 0.00***
Very sharp curvature (VS
Curve)
no −14.61 0.00***
Median is present (PM) no −3.04 0.00***
Scale parameters for the distributions of random parameters
AADT (Std. Dev.) normal −0.69 0.49 RD (Std. Dev.) tn_normal −4.83 0.00***
Curve50 (Std. Dev.) normal −0.02 0.98 Cholb
: VS Curve (Std. Dev.
ln_normal)
−11.70 0.00***
RSHS (Std. Dev.) normal 0.04 0.97 Chol: PM (Std. Dev.
tn_normal) . VS Curve (Std.
Dev. ln_normal)
−1.30 0.19
Chol: PM (Std. Dev.
tn_normal)
−0.16 0.87
Model type
Dispersion value
Log-Likelihood
BIC
2.04
−198.57
481.80
Negative binomial Negative binomial
2.98
−171.44
462.51

a Bolded text to indicate shared across both models.

identified as Pareto-efficient by the extensive hypothesis testing framework, have remarkably outperformed the benchmark model as well. This result underscores the potential for more efficient models to be proposed by analysts if this approach was adopted (Ahern et al., 2024). Furthermore, even when faced with datasets burdened by excessive zeros, adopting a complex approach that balances numerous contributing factors with correlated random parameters proves to be a viable strategy to address the issue of excessive zeros, but more research is required to address this problem more effectively (Ali et al., 2022).

The hypothesis put forth by the proposed framework posits that the utilization of MCV heightens the probability of crashes, and notably affirmed by existing literature (Ackaah and Salifu, 2011; Mannering et al., 2016). Factors such as vehicle size and weight, limited visibility, extended stopping distance, likelihood of driver fatigue, and mechanical issues inherently render these MCVs riskier than other vehicles (AASHTO, 2011; Akgüngör and Doan, 2009). The framework has also drawn a correlation between MCV usage and crash risk, utilizing the arcsinh transformation. This outcome implies that a substantial increase in MCV occupancy has a waning influence on crash risk, substantiating the sustainability of current standards governing MCV use. Behara et al. (2021) reported a similar finding, where a logarithmic transformation exhibited traits akin to arcsinh. However, within the MCV context, the arcsinh function proves more fitting than the logarithmic function. The arcsinh function, unlike the commonly used logarithm, is particularly adept at handling skewed data. Additionally, it maintains the sign of values through a more gradual growth pattern. These attributes make it a better match for the specific traits of MCV data.

The proposed framework reveled that increasing speed on a highway has the potential to increase the frequency of crashes. This finding is consistent with previous studies that have identified various mechanisms through which speed contributes to crashes, including increased risky driver behavior, reduced driver control, longer stopping distance, and slower reaction time (Garber and Ehrhart, 2000; Islam et al., 2014; Yasmin et al., 2022; Cai et al., 2021; Ponnaluri, 2018). Although addressing speed can be challenging due to its impact on traffic flow and driver satisfaction, it is important to learn about a potential speeding problem (Das and Sun, 2016). The results suggests that rural double highways are significantly less prone to crashes. Therefore, future road design, should prioritize mitigation strategies for other road classes before focusing on rural double highways.

Additionally, very sharp curvature in the road has been found to reduce crashes, which may initially seem counter-intuitive because curves are typically associated with increased crash risk (Gooch et al., 2016; Yahaya et al., 2016; Behara et al., 2021). However, highways with sharp curvature could benefit from risk-adjusting behaviors, such as drivers reducing their speed (Oviedo-Trespalacios et al., 2017). Moreover, the implementation of effective mitigation strategies on Queensland highways, such as enhanced signage, delineation, highfriction pavements, and wider clear zones, is likely contributing to the reduction in observed crashes (Antonucci et al., 2004; Gooch et al., 2016). This is evident despite the historical occurrence of crashes on sharply curved highways compared to other regions (Onate-Vega et al., 2020).

Lastly, the estimated model suggests a correlation between roadway curvature and the presence of a median. This indicates that significant curvature, when complemented by a median, furnishes an additional safety layer, resulting in a decrease in crashes. The negative skew associated with the proposed lognormal distribution reinforces the notion that the safety risks of curves can be effectively mitigated by the inclusion of a median. The median, denoted as PM, serves several safety-enhancing functions. It deters crossover-roadway crashes, provides a protected waiting area for turning vehicles, restricts direct access from external vehicles, and aids in preventing head-on collisions (Donnell and Mason, 2006; Zou and Tarko, 2016).

b Chol: Cholesky element which has captured the correlations between the two contributing factors.

Fig. 5. Minimizing BIC and RMSE via extensive hypothesis testing algorithm: Pareto-frontier.

While our analysis only examined one model in the Pareto-frontier, it is important to note that the analyst should ultimately determine the final specification. Other solutions within the Pareto-frontier may provide further insights, nuances, and features that could be useful. For example, the chosen model highlights the importance of homogeneous sections along rural double highways, incorporating medians, and probes the feasibility of including medians in their current state on sharply curved highways.

5.5. Experiment using data from washington

After conducting an extensive hypothesis testing using data from Washington, a Pareto-frontier was obtained. Results are provided in Fig. 5.

Numerous studies, including those by Miaou (1994) and Behara et al. (2021), have indicated AADT and segment length as significant predictors of road crashes. As such, the product of ( ) and ℎ6 was included in this analysis. Both AADT and segment length rank among the most relevant variables for predicting crash frequency on road segments and capture vehicular density on roads (Mehta and Lou, 2013; AASHTO, 2011; Haghani et al., 2019; Kim et al., 2013). By specifying an offset term, the final model ensures the inclusion of vehicular traffic density, which guides other modeling specifications and helps alleviate bias in the estimates (Mannering et al., 2016). Moreover, the offset term employed in the original study was ( ) × ℎ. Our extensive hypothesis framework, therefore, provides an opportunity to uncover insights akin to those disclosed in Mannering et al. (2016). The framework also provides alternative insights that can be compared against the initial study by Washington et al. (2020), leading to a better understanding of the relationship between exposure factors and crash frequencies.

Beyond exposure, AADT per lane was found to correlate positively with crash frequency. While this effect could, to some extent, be attributed to the offset, it nonetheless underscores the significance of high-density traffic per lane. In essence, as lanes are added, safety potentially decreases due to reduced traffic load per lane. Our findings, which align with previous research (Mecheri et al., 2017), suggest that adding lanes could decrease the risk of crashes, due to the decrease in traffic density per lane.

The model derived using the proposed framework indicates that at peak-hour for traffic, the probability of crashes is more likely. Milton and Mannering (1998) posited that highways with higher traffic volumes tend to have a larger share of activity throughout the day. As a result, these roads are generally more prone to crashes. Conversely, roads with lower volumes tend to have a high proportion of their traffic during peak hour traffic. Thus, Milton and Mannering (1998) suggested that despite high percentage of peak hour traffic, roads with lower volumes roads in AADT are likely to witness a reduced overall crash frequency. This is attributed to their overall low traffic volume, which presents fewer opportunities for crashes. However, our framework observed an inverse relationship, and is bolstered by other research (Das et al., 2021a). For instance, Chan et al. (2010) proposed that during peak hour traffic, drivers on low volume roads may be likely to exhibit aggressive behaviors that heighten the chances of conflicts with neighboring vehicles. This could include decelerating for rutting or changing lanes to bypass highway cracking. This is likely to occur when high proportion of peak hour traffic causes congestion.

Several factors associated with grade breaks, have been observed to heighten crash frequency. Intriguingly, this finding contrasts the assertions in Milton et al. (2008), which postulate that grade breaks contribute to the prevention of serious crashes by prompting riskadjusting behavior. Nevertheless, it is crucial to weigh the potential influence of risk compensation, which arises from the perception of worsening conditions (Chen et al., 2017). Factors such as diminished visibility, vehicle control difficulties, increased wear and tear, drainage issues causing hydroplaning, and sudden, unexpected alterations in road grade resulting from grade breaks can all amplify the likelihood of crashes (Rusli et al., 2018; Das et al., 2018). Grade breaks, being natural obstacles, challenge drivers with an additional hurdle that necessitates a change in driving behavior (Rezapour et al., 2019).

Similarly, the number of interchanges has been observed to marginally increase the probability of crashes, contradicting a counterintuitive finding in Venkataraman et al. (2011a) where the same data was employed. Supporting our findings, a study by Feknssa et al. (2023) suggested that interchanges increase the number of conflict points. Therefore, it is imperative to thoroughly consider these factors and their implications during road design. Minimizing and clearly marking grade breaks, ensuring gradual changes in road grade, and reducing interchanges can afford drivers ample time to adjust their speed and driving style, thereby mitigating the risk of crashes.

HISNOW is included with random parameters in both the models from Table 4. Inclement weather conditions generally increase driving hazards, prompting most drivers to adopt more risk-averse behaviors. This could, in some cases, reduce the frequency of crashes as people may choose not to drive or drive slower than during normal conditions. However, the reduced visibility during inclement weather conditions also has been found to increase the risk of crashes (Das et al., 2021a). This complex interplay is well-captured by the use of random parameters. These parameters account for the variability in conditions and behaviors; even if drivers are more risk-averse, some roads may still present excessive danger due to poor traction. The same finding is supported by Jackson and Sharif (2014); drivers tend to travel at lower speeds when it rains, thereby reducing the likelihood of crashes.

6 Also known as the exposure index (AASHTO, 2011).

Table 4 Specification found by Washington et al. (2020) compared to the extensive hypothesis testing framework.

Specification found By Washington et al. (2020) Poisson model found through the HS algorithm
Effect 𝜏 Coeff Pr 𝑧 > 𝑍 Effect 𝜏 Coef Prob 𝑧 > 𝑍
Constant no 2.82 0.00*** Constant no −3.84 0.00***
Curvature per mile (CPM) no −0.21 0.02* CPM
no
0.72 0.00***
Segments where average
snowfall per month is above 1
inch (HISNOW)
no −0.10 0.68 HISNOW no 0.03 0.67
Friction value (0 to 100 with 100
being high) (FRICTION)
no −0.41 0.58 Minimum grade in the segment
no
(MIGRADE)
0.78 0.00***
Average precipitation in a month
is less than 1.5 inches (LOWPE)
no −0.36 0.78 Maximum grade difference in
the segment (MXGRDIFF)
arcsinh 0.86 0.00***
Number of grade breaks per mile
(GBRPM)
no 0.07 0.82 Number of interchanges in the
segment (INTECHAG)
no 0.65 0.00***
Interchanges per mile (INTPM) no −0.34 0.22 The percentage of daily transit
in peak hour (PEAKHR)
arcsinh 2.72 0.00***
Number of grade breaks in the
segment (GRADEBR)
arcsinh 2.00 0.00***
GBRPM arcsinh 0.84 0.00***
Average daily travel per lane
(ADTLANE)
arcsinh 0.97 0.00***
Median width (MEDWIDTH) no −0.14 0.15
Scale parameters for distributions of random parameters
INTPM (Std. Dev. normal) 0.93 0.00*** PEAKHR (Std. Dev. tn_normal) 0.23 0.00***
CPM (Std. Dev. normal) 0.89 0.00*** GRADEBR (Std. Dev. normal) 0.17 0.03*
HISNOW (Std. Dev. normal) 0.63 0.00*** GBRPM (Std. Dev. ln_normal) −0.35 0.00***
HISNOW (Std. Dev. tn_normal) 0.29 0.00***
Chol: MEDWIDTH (Std. Dev.
normal)
0.17 0.02*
Chol: ADTLANE (Std. Dev.
ln_normal) . MEDWIDTH (Std.
Dev. normal)
1.22 0.00***
Chol: ADTLANE (Std. Dev.
ln_normal)
0.52 0.00***
Chol: CPM (Std. Dev.
ln_normal) . MEDWIDTH (Std.
Dev. normal)
1.67 0.00***
Chol: CPM (Std. Dev.
ln_normal) . ADTLANE (Std.
Dev. ln_normal)
−1.40 0.00***
Chol: CPM (Std. Dev. ln_normal) −1.17 0.00***
Model type
Dispersion value
Log-Likelihood
BIC
1.42
−701.86
1467.55
Negative binomial Poisson
–
−617.62
1367.42

Median width, was also included in the model. With this it was found that the transition from larger to smaller median widths did not significantly influence the number of crashes, with only a slight uptick in crash frequency as median widths diminished. This outcome aligns with the typical use of median barriers, primarily installed to avert cross-median crashes. The choice of these barriers is governed by factors such as median width and traffic volumes, a notion supported by numerous studies (Donnell and Mason, 2006; Abdel-Aty and Radwan, 2000; Park and Kim, 2010). However, the dynamics change when considering correlated random parameters for effects curvature and lane density together, under which the crash frequency can fluctuate considerably.

Large median widths may instill a sense of safety in drivers. Paradoxically, this could inadvertently prompt reckless driving; an effect that becomes especially noticeable when considering road curvature and high lane traffic density. The implication is that the perceived safety afforded by large median widths, particularly considering challenging scenarios such as driving along curves and high traffic densities, may inadvertently lead to an increase in crashes. Some research suggests that larger medians might inadvertently foster less safe conditions due to their perceived safety connotation (Telford et al., 1953). According to Anastasopoulos et al. (2012), median width impacts safety in subtle ways. Narrow medians typically feature barriers, but as the median width expands, the slopes become gentler, potentially mitigating the risk of severe overturn crashes. This concurs with our finding that increasing median width can help alleviate challenging conditions such as driving along curves, which could escalate the risk of overturn crashes. If the median width does not meet standards, it may fall short in preventing severe median crossover crashes (Shankar et al., 1996).

The original benchmark model by Washington et al. (2020) registered a decrease in crash likelihood as curves per mile increased. In contrast, our model, alongside (Behara et al., 2021; Gooch et al., 2016; Kvasnes et al., 2021), suggests a different, and arguably more intuitive, effect. However, by using lognormal random parameters, negative skew in the standard deviations across road segments was discerned. This skew implies a trend towards a decreased mean crash probability as curves per mile rises — an insight somewhat aligned

Fig. 6. Minimizing BIC and RMSE via Extensive hypothesis testing algorithm: Pareto-frontier.

with the benchmark model (Washington et al., 2020). Interestingly, this negative skew is further intensified by the contributing factor AADT per lane, which are correlated with each other, indicating that under specific conditions, lanes with high density and significant curvature might actually reduce crash likelihood more than the mean effect of increasing the number curves per mile alone. Although this may initially seem counterintuitive, it underscores the complex interplay of factors influencing crash frequencies and reiterates the nuanced role that road characteristics play in driving safety. As learned using data from Queensland, significant curvature can result in fewer crashes, likely due to driver behavior and effective mitigation strategies. The enhanced safety associated with lane AADT amidst curvature is likely due to compelling drivers to exercise caution, so as to avoid encroaching on neighboring lanes, which are densely occupied.

The selected model integrates correlated random parameters, a trend also present across the Pareto frontier solutions. This aligns with Saeed et al. (2019), where both correlated and non-correlated random parameters yielded similar model estimations. However, the focus on uncovering new insights surpasses the importance of estimation alone. For example, the identification of the increased crash rate on roads with increasing median widths, significant curvature and high traffic density, would not have been identified without estimating and analyzing the correlated random parameter model. This highlights the significance of an extensive hypothesis testing that can consider correlated random parameters for obtaining deeper insights that cannot be captured without considering correlations. It is essential to recognize that the presence of correlations between the random parameters relies on the characteristics of the data and may not be universally present in other cases. Therefore, testing for the existence of such correlations becomes necessary. The significance of the proposed framework lies in its extensive search for the most suitable model, whether it is a uncorrelated fixed, random or correlated random parameter model, along with other characteristics tested within the mathematical program.

5.6. Experiment using synthetic data

The extensive hypothesis testing using synthetic data yielded a Pareto-frontier comprising fourteen distinct solutions. One of these solutions within the Pareto-frontier accurately identified the synthetic model. Only solutions that lie relatively close to the Pareto-frontier have been highlighted and are presented with a level of visual transparency to indicate that they are not Pareto-efficient. synthreviewerThe results are presented in Fig. 6, while the output of the compared models are presented in Table 5.

The model search was driven by the RMSE and the BIC. However, the search results would likely differ if alternative objectives were employed, such as HQIC. In such a scenario, even if the search stumbled upon the synthetic specification, it may not have been deemed Paretoefficient given the use of HQIC instead of BIC. This underscores the role of an analyst to apply their expertise in selecting objectives and using constraints to guide model selection.

Considering the RMSE and BIC objectives, the synthetic model selected for analysis showed to be Pareto-efficient, despite our inability to capture the true model. This outcome highlights the importance of utilizing both objectives during model search, as it aids in identifying models that offer a balance between fit and precision. The BIC, known for favoring simpler models that effectively explain the data, addresses the bias–variance trade-off and reinforces its utility in model selection (Burnham and Anderson, 2004). Our chosen model, which mirrors the structure of the final model but lacks one variable, exemplifies this. Although the omitted variable was used in model generalization, this shows that a estimated model can outperform the synthetic true model in one objective, indicating a trade-off, subject to analyst preference.

The use of the epsilon-constraint method was highlighted as crucial, as demonstrated by the Pareto-frontier. Without the Pareto-frontier approach, a singular minimization of the BIC objective function may have misleadingly indicated that the proposed model in Table 5 is indeed the superior statistical fit. This is because the proposed model, despite representing the smallest BIC on the Pareto-frontier, inevitably yields a worst prediction accuracy as opposed to the true model. If a single-objective method were adopted optimizing BIC, this critical detail could be overlooked, and even if our framework tested for the true model, it would have been rejected. BIC, considers the number of parameters and penalizes excessive complexity, typically favoring parsimonious models over those with greater generalization. However, while this supports parsimonious models, generalization can still be a concern.

Therefore, to address this issue, it is crucial to recognize that generalizability does not always guarantee prediction accuracy. Hence, it was necessary to incorporate the second objective to mitigate this potential problem (Beeramoole et al., 2023). This is further supported by the Pareto-frontier. Which reveals that an excessive emphasis on goodness of fit can compromise the second objective; RMSE, which is a crucial metric for out-of-sample prediction. In this instance, a trade-off is discerned in the synthetic model between these two objectives. If an additional objective function were not incorporated, the accurate identification of the data-generating model would have been unattainable under strict minimization of BIC or RMSE alone.

Lastly, the insights drawn from this experiment underscore the challenges inherent in recreating the synthetic model, particularly when additional factors escalate complexity and amplify the likelihood of encountering local optima. Nevertheless, the proposed framework has showed its capability in capturing most of the significant effects in the model, and further identified the correlation structure between

Table 5 Comparison of results between true and estimated model.

Effect True model specification Model found through the proposed framework
𝜏 Coeff Prob 𝜏 Coeff Prob
Constant no −1.53 0.00*** no −2.54 0.00***
X1 no 1.15 0.00*** no 1.27 0.00***
X2 no −1.55 0.00***
X3 no 2.52 0.00*** no 2.99 0.00***
Scale parameters for distributions of random parameters
X1 (Std. Dev. normal) 0.78 0.00*** 1.21 0.00***
Chol: X2 (Std. Dev. normal) 0.61 0.00***
Chol: X3 (Std. Dev. normal) 0.51 0.00*** 1.35 0.00***
Chol: X3 (Std. Dev. normal). X1 (Std. Dev. normal) 0.58 0.00*** 1.42 0.00***
Chol: X2 (Std. Dev. normal). X1 (Std. Dev. normal) 1.02 0.00***
Chol: X3 (Std. Dev. normal). X2 (Std. Dev. normal) 0.95 0.00***
Model type Poisson Poisson
log-likelihood −2429.106 −2430.166
BIC 4930.66 4903.80
RMSE 3.11 3.97

factors ''X1'' and ''X3''. Moreover, the choice of a Poisson distribution over other alternatives exemplifies success pinpointing the most fitting closed-form model.

6. Conclusion

This paper proposed an unbiased multi-objective optimization framework that significantly enhances the development of crash count data models. Traditional estimation methods, as highlighted in previous studies (Mannering, 2018; Lord and Geedipally, 2011; Saeed et al., 2019), may be prone to bias and overlook unique specifications and insights present in the data. In contrast, the proposed framework generates and tests a broad range of hypothesis considering simultaneously likely contributing factors, unobserved heterogeneity, transformations, and correlations. The framework includes a mathematical programming formulation that minimizes two objectives: the goodness of fit for in-sample data and the assessment of out-of-sample predictions. This approach seeks a balance between model generalizability and accuracy. The framework encompasses various metaheuristic solution algorithms, enabling efficient search for parsimonious specifications and facilitating its applicability to various datasets. The effectiveness of the framework was evaluated using two real-world datasets and one synthetic dataset. Results were compared with those published using the same data. This evaluation and comparative analysis illustrated the ability of the proposed framework to estimate parsimonious crash count data models efficiently, capable of proposing complex models whilst reducing barriers in model development such as time, bias', and ad hoc processes (Beeramoole et al., 2023; Arteaga et al., 2022; Paz et al., 2019).

The experiment utilizing synthetic data demonstrated the potential of multi-objective optimization and chosen functions to capture complex specifications involving correlation. Results confirmed that the proposed approach can accurately capture the synthetic model used to create the data while avoiding local optima regions.

Results using data from Queensland yielded a model consistent with an early study (Behara et al., 2021). The model adequately captured the dispersion in the data while identifying similar contributing factors. Notably, the framework also discovered additional factors that were not previously identified, demonstrating its ability to uncover more comprehensive insights. The results underscored the importance of mitigating grade breaks, and recognizing the significance of the median. Further exploration of other solutions within the Pareto-frontier may reveal additional likely contributing factors that influence the underlying crash generation process.

An analysis of data from Washington uncovered a range of potential contributing factors to traffic incidents, including the percentage of daily traffic during peak hours, heavy snowfall, grade breaks, and various road characteristics such as the number of interchanges, curvature per mile, and median width. The results underscore the intricacies of these elements, demonstrating that, contrary to expectations, higher peak-hour traffic, extreme weather conditions, and wider medians could actually lead to an increase in crashes, potentially due to driver behavior and varying road conditions. Additionally, the study highlighted the significance of incorporating correlated random parameters to gain deeper insights because of the interplay between modeling effects, which non-correlated random parameters model overlook (Amiri et al., 2020b; Afghari et al., 2018). As such, consideration of correlated random parameters enhances the capability of the proposed framework to identify critical factors affecting road safety, and also helps to alleviate biased estimates.

Further research could explore ways to improve the multi-objective framework as well as testing with alternative models and data sets. There are several extensions that could further enhance the framework such as considering heterogeneity in the means (Behara et al., 2021; Ali et al., 2022). Another extension worth considering is the inclusion of grouped random parameters (Amiri et al., 2020a; Heydari et al., 2016), which can account for group-specific effects by considering random effects specific to different groups. Considering simultaneously heterogeneity in the means and grouped random parameters would enable a more extensive hypothesis testing and search for effects associated with different subgroups or factors.

CRediT authorship contribution statement

Zeke Ahern: Writing – review & editing, Writing – original draft, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Paul Corry: Writing – review & editing, Supervision, Methodology. Wahi Rabbani: Writing – review & editing, Supervision. Alexander Paz: Writing – review & editing, Validation, Supervision, Methodology, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Appendix

See Table 6.

Table 6 Notation: Representation of the solution and the properties of metaheuristics.

Experimental default parameters
𝑀𝑎𝑥𝑇 𝑖𝑚𝑒 Max allowable time for each algorithm.
𝑊 𝐼𝐶 Number of iterations, without an improvement in the Pareto-Efficient Set 
Solution representation
𝑣𝑝 Solution vector to the mathematical program, where 𝑝 depicts an indexed solution.
𝑣𝑝,𝑘 Specification of decision element 𝑝 , ∀𝑘 ∈ ℧.
𝓁𝑝,𝑘 Index value of the solution 𝑝 for decision 𝑘.
Harmony search
𝑃 𝐴𝐼 Pitch adjustment index.
𝐻𝑀 Set of distinct solutions (𝑣𝑝
) in the harmony memory.
𝐻𝑀𝑆 Harmony memory size (controls the set size of 𝐻𝑀).
𝐻𝑀𝐶𝑅 The probability of retaining the previous solution
Differential evolution
𝐶𝑅 The probability of a crossover occurring on 𝑣𝑝,𝑗
𝑃 𝑜𝑝𝑆𝑖𝑧𝑒 The population size.
𝐴𝐼 Adjustment index.
Simulated annealing
𝑇𝑠𝑡𝑒𝑝 Number of iterations before the temperature is reduced.
𝑇𝛼 Temperature reduction factor.

A.1. Acronyms

AADT Annual average daily travel

ADTLANE Average daily travel per lane

AIC Aikike's Information Criterion

BIC Bayesian Information Criterion

CPM Curvature per mile

CRP Correlated random parameters

Curve50 Significant curvature (moderate to very sharp which is longer than 50% of the segment)

DE Differential Evolution

FRICTION Friction value (0 to 100 with 100 being high)

GBRPM Number of grade breaks per mile

GP Generalized Poisson

GRADEBR Number of grade breaks in the segment

HISNOW Segments where average snowfall per month is above 1 inch

HQIC Hannan–Quinn Criterion

HS Harmony Search

INTECHAG Number of interchanges in the segment

INTPM Interchanges per mile

LENGTH Segment length in miles

LOWPE Average precipitation in a month is less than 1.5 inches

M Roadtype is motorway

MAE Mean absolute error

MCV Multi-combination vehicular traffic

MEDWIDTH Median width

MIGRADE Minimum grade in the segment

ML Machine learning

MSL Maximum simulated log-likelihood

MXGRDIFF Maximum grade difference in the segment

NB Negative Binomial

NLANES Number of lanes

NN Neural network

PEAKHR The percentage of daily transit in peak hour

PM Median is present

RD Rural dual carriageway

RF Random forest

RMSE Root-mean squared error

RP Random parameters

RSHS Rural single and high speed

RSMS Rural single and medium speed

SA Simulated Annealing

SP Speed limit in Km/Hr

SPEA2 Strength Based Pareto Evolutionary Algorithm 2

SPF Safety performance function

SVM Support vector machine

SW Shoulder width

US Urban single carriageways

VS Curve Very sharp curvature

References