Driver Injury Prediction and Factor Analysis in Passenger Vehicle-to-Passenger Vehicle Collision Accidents Using Explainable Machine Learning

1. Introduction

Passenger vehicles (PVs), as motor vehicles designed to transport people and their belongings, have significantly enhanced the convenience of human mobility [1,2,3]. PVs are primarily categorized into two major groups: passenger cars and light trucks [4]. Among various vehicle types, PVs not only constitute the majority of the vehicle fleet but also account for the highest proportion of traffic accidents [5]. According to the 2024 analysis conducted by the National Highway Traffic Safety Administration on the 2022 Crash Report Sampling System data [6], more than 64% of police-reported crashes involved two-vehicle collisions, with PV-PV crashes accounting for the highest proportion at 73.4%. Drivers are the most severely affected individuals in traffic accidents, representing 67.7% of all road traffic casualties in 2021. Therefore, understanding the factors influencing collision occurrences and injury severity is crucial for enhancing traffic safety.

PV-PV collisions constitute a complex dynamic process, as illustrated in Figure 1, encompassing dual-vehicle and driver characteristics, environmental conditions, causal factors, and pre-collision and collision dynamics [7]. accident severity prediction and feature analysis involve four key steps [8]: data acquisition, preprocessing, feature selection, and model selection. First, obtaining a reliable and comprehensive accident dataset is foundational [9]. Second, preprocessing steps like outlier removal, missing value imputation, encoding, and balancing ensure data suitability for analysis [10]. Third, feature selection addresses redundancy and enhances prediction performance and interpretability [11]. Finally, selecting an appropriate model tailored to the dataset and research objectives is critical for effective analysis.

Recent studies have focused on identifying and quantifying factors influencing accident injury severity, primarily through statistical methods. For instance, Dongdong Song et al. (2023) utilized a random parameters bivariate probability model to analyze UK truck-car collision data (2017–2019), identifying key factors affecting driver injury severity and proposing safety improvements [12]. Similarly, Hongren Gong et al. (2022) applied a Bayesian multivariate random parameters logit model to U.S. GES and CRSS two-vehicle collision data (2010–2018), uncovering significant associations between factors and injury risk, along with heterogeneity and temporal instability in effects [13]. In another study, Joon-Ki Kim et al. (2013) analyzed California single-vehicle crash data (2003–2004) using a mixed logit model, highlighting the influence of age and gender on injury severity and the presence of heterogeneity [14]. Donald Mathew Cerwick et al. (2014) compared mixed logit and latent class models using multi-vehicle heavy truck crash data from Iowa (2007–2012), finding that while the latent class model fit slightly better, the mixed logit model yielded probability predictions closer to observed values, with similar marginal effects for most variables [15].

Statistical methods have been widely used due to their strong interpretability, aiding in understanding the mechanisms behind accident heterogeneity and severity correlations. However, their dependence on linear assumptions and independence constraints severely limits their ability to model the complex, nonlinear interactions among accident-related factors. [16]. Traffic accidents, as complex systems, involve multiple collinear and higher-order interactions among influencing factors [17,18], such as vehicle dynamics and environmental conditions. The results obtained from traditional methods heavily depend on the validity of their assumptions. If these assumptions deviate from reality—such as neglecting nonlinear interactions—the reliability and generalizability of the analysis are significantly compromised [19,20,21]. In contrast, machine learning (ML) methods excel in nonlinear modeling and can effectively capture intricate relationships among accident features without relying on strict statistical assumptions. They can extract latent patterns from high-dimensional data, substantially improving the accuracy of accident severity prediction [22,23].

In recent years, ML techniques have emerged as a promising alternative to traditional data analysis methods for injury severity analysis, offering the ability to extract valuable insights from large-scale, complex, and heterogeneous datasets. Ensemble machine learning (EML) and deep learning (DL) models are widely used in traffic accident severity prediction, but their applications differ. DL models are mainly used for macroscopic accident analysis, utilizing spatial-temporal data from images, videos, or sequences to predict accident occurrence and severity, aiding traffic management in intervention and emergency response. Yu et al. (2021) developed a Deep Spatio-Temporal Graph Convolutional Network using heterogeneous traffic data in Beijing to improve accident risk prediction, the model enhances prediction accuracy, enabling early warnings of potential hazards and supporting safer route selection [24].Fares Alhaek et al. (2024) proposed a CNN-BiLSTM-based DL method for traffic accident severity prediction. Experiments on UK city data showed it outperformed baselines [25]. However, due to the structured nature of accident severity datasets, Tabular is the most basic form of traffic accident datasets. EML models has demonstrated its superiority in handling datasets with tabular data format. EML’s strength lies in its ability to combine multiple base models, enabling it to capture complex patterns and relationships within the data more effectively than single models. In the field of accident severity prediction, this advantage is particularly evident. For example, various ensemble machine learning models have outperformed other machine learning models. Models like Random Forest [26], XGBoost [27,28], LightGBM [29], CatBoost [30] and AdaBoost [31] are more suitable for microscopic accident severity prediction. These models offer better interpretability, which means it is easier to understand how different factors contribute to accident severity.

EML algorithms improve predictive performance by integrating multiple decision trees. Although a single decision tree offers inherent interpretability, the depth and number of nodes in ensemble models significantly increase in high-dimensional accident feature spaces, leading to excessive model complexity and a notable reduction in interpretability [32]. Despite their superior predictive performance, the "black-box" nature of ensemble models hinders the interpretation of accident causation mechanisms. Explainable ML techniques, SHAP [33], provide a viable solution to this issue. SHAP, based on Shapley values from game theory, quantifies feature contributions and reveals interaction effects. Its additive feature attribution property aligns well with tree-based ensemble models, supporting both global and local interpretability [34]. Widely integrated with ensemble machine learning, SHAP has become a preferred tool for analyzing the relationship between accident severity and contributing factors, offering a new paradigm for passenger vehicle collision analysis. It not only ensures the accuracy of the prediction model but also enables in - depth interpretation of the model and feature analysis.

Current research on accident severity prediction faces several challenges. First, natural accident datasets suffer from severe class imbalance [6,12,13], where fatal accident cases are significantly underrepresented compared to non-injury cases. This imbalance poses challenges for traditional resampling methods in effectively balancing severity levels [35]. Moreover, accident datasets often contain a mix of categorical and numerical features, which limits the applicability of Euclidean distance-based oversampling techniques such as Synthetic Minority Over-sampling Technique (SMOTE) [36]. Thus, addressing data imbalance at the data collection stage is crucial. Second, existing studies primarily rely on empirical heuristics for feature selection [12], lacking rigorous scientific justification. Consequently, the representativeness and rationality of selected features remain inadequately explained.

To address the issue of class imbalance in accident severity datasets, this study first integrates the CRSS and FARS datasets from the NHTSA. This integration reduces the severity imbalance ratio from 175.1 in CRSS and 10.2 in FARS to 5.1 in the combined dataset. Furthermore, a combination of oversampling and undersampling techniques is employed to achieve a balanced distribution of severity categories. For feature selection, traditional methods often rely on Pearson correlation coefficients to eliminate highly correlated and irrelevant variables. However, this approach is not suitable for accident datasets dominated by categorical features, as categorical values represent discrete labels rather than numerical relationships. o enhance the rationality of feature extraction, this study adopts the mRMR algorithm based on mutual information [37]. This algorithm has been widely applied in Disease Diagnosis [38,39]. Compared with other feature extraction methods, it demonstrates the best performance while requiring the fewest number of features [40]. Subsequently, this study applies ensemble machine learning algorithms—including CatBoost, AdaBoost, XGBoost, Random Tree, LightGBM, and decision trees—to accident severity prediction. By comparing the predictive performance of different models, the best-performing model is selected for SHAP-based interpretability analysis, facilitating insights into feature importance and improving the explainability of accident severity predictions.

The contributions of this paper can be summarized as follows:

(1) Integrated two distinct accident datasets to address imbalance issues and applied the mRMR algorithm for key feature extraction, reducing redundancy and improving model accuracy across all severity levels.

(2) Proposed a novel micro-level analytical approach for PV-PV accidents, incorporating dual-vehicle factors to simultaneously enhance prediction accuracy and model interpretability, advancing beyond traditional macro-level methods.

(3) Performed a systematic evaluation of state-of-the-art EML models for PV-PV accident severity prediction and integrated the optimal ensemble model with SHAP interpretation, enabling comprehensive feature analysis and addressing a critical research gap in traffic safety analytics.

The remainder of this paper is organized as follows. Section 2 describes the methodologies, encompassing data preprocessing techniques, feature extraction methods, predictive modeling approaches, and explainable AI algorithms. Section 3 presents the dataset’s statistical characteristics and feature analysis. Section 4 details the experimental results, including model performance evaluation and critical factor analysis. Section 5 discusses the research implications and study limitations. Finally, Section 6 concludes the paper and proposes future research directions.

2. Methods

The technical route of this study is shown in Figure 2, which is mainly divided into three parts: the first part is data collection and preprocessing, which aims to provide a structured data set that meets research needs to improve the accuracy and reliability of the model; the second part is model training and tuning, which aims to improve the predictive performance and generalization ability of the model by optimizing model parameters; the third part is to interpret and analyze the characteristics of the optimal model to reveal the key factors behind the model decision and enhance the interpretability and transparency of the model.

2.1. Accident Data Collection and Preprocessing

2.1.1. Accident Data Collection

In accident severity prediction research, class imbalance is a common issue, particularly due to the relatively low number of fatal accident samples, which limits models’ ability to predict severe accidents. Additionally, these datasets typically contain both categorical and numerical features, making traditional Euclidean distance-based SMOTE algorithms less effective for data balancing. To address dataset imbalance, this study integrates two complementary NHTSA datasets: CRSS and FARS. The FARS dataset, containing comprehensive fatal accident records, effectively compensates for the scarcity of severe cases in CRSS. This integration enhances the representation of high-severity accidents at the data source level, which improves the model’s predictive accuracy for severe accident classification.

2.1.2. Dataset Preprocessing

1) Dataset Compilation:

The accident-related data are extracted from multiple relational tables within both CRSS and FARS datasets, including accident records, person information, event sequences, and vehicle damage reports. These heterogeneous data sources are consolidated into a unified structured dataset through a systematic integration process. A composite primary key, comprising Accident ID and Driver ID, is implemented to ensure precise data linkage and maintain referential integrity across all integrated tables.

2) Case Selection:

This study focuses on predicting driver injury severity in PV-PV collisions. Accordingly, vehicle collision cases involving exactly two passenger vehicles with corresponding driver characteristics were extracted from the integrated dataset, ensuring alignment with the research focus.

3) Data Cleaning:

Features with excessive missing values were removed. Features with minimal missing values were imputed using a context-based approach. Cases lacking critical information, such as driver injury severity levels, were excluded to ensure data quality.

4) Encoding and Normalization:

Categorical feature encoding in the original dataset was not directly related to injury severity. Although decision tree-based ensemble models are unaffected by encoding order, it influences the interpretability of SHAP-based feature analysis. Therefore, target encoding was applied, ranking categorical features in ascending order based on their mean injury severity levels. Numerical features were standardized to eliminate scale differences, ensuring model stability and convergence. Minor injury categories (e.g., Possible injury and Suspected minor injury) lacked sufficient distinction. Thus, injury severity levels were reclassified as follows:

No apparent injury → Level0

Possible injury & Suspected minor injury → Level1

Suspected serious injury → Level2

Fatal injury → Level3

5) Data Balancing:

The integration of CRSS and FARS inherently mitigated class imbalance. To further balance the dataset, a combination of random oversampling and undersampling was applied [41]. The average sample size across injury severity levels was used as a reference: oversampling was conducted for underrepresented classes, while undersampling was applied to overrepresented classes. To prevent overfitting and enhance model generalization, categorical features remained unchanged during oversampling, whereas numerical features were perturbed by adding random noise in the range of -5% to +5%. The total dataset size remained constant before and after balancing.

2.2. Feature Selection

Feature selection is a crucial step in accident severity prediction, as redundant and irrelevant features can negatively impact model performance. In this study, the mRMR [37] method was employed for feature selection in accident severity prediction. Compared to traditional correlation-based feature selection methods (e.g., Pearson correlation), mRMR is more effective for datasets with categorical features, as it does not assume linear relationships. The mRMR algorithm is a widely used method for selecting an optimal subset of features by balancing two key objectives:

1. Maximum Relevance: The selected features should have the highest relevance to the target variable.

2. Minimum Redundancy: The selected features should be minimally correlated with each other to reduce redundancy.

2.2.1. Mutual Information (MI) Definition

The mutual information $I(X;Y)$ quantifies the dependency between two variables X and Y, and is defined as:

$$I(X;Y)=\sum _{x\in X}\sum _{y\in Y}p(x,y)log{\displaystyle \frac{p(x,y)}{p\left(x\right)p\left(y\right)}}$$

(1)

where $p(x,y)$ represents the joint probability distribution, and $p\left(x\right)$ and $p\left(y\right)$ are the marginal distributions of X and Y, respectively.

2.2.2. Maximum Relevance Criterion

The relevance of a feature set S with respect to the target variable c is computed as the average mutual information between each feature $X_i$ and c:

$$\mathcal{R}(S,c)={\displaystyle \frac{1}{\left|S\right|}}\sum _{X_i\in S}I(X_i;c)$$

(2)

where $\left|S\right|$ is the number of selected features. The goal is to maximize this relevance measure.

2.2.3. Minimum Redundancy Criterion

To ensure that the selected features are not highly correlated with each other, the redundancy term is defined as the average mutual information among selected features:

$$\mathcal{D}\left(S\right)={\displaystyle \frac{1}{{\left|S\right|}^2}}\sum _{X_i,X_j\in S}I(X_i;X_j)$$

(3)

The goal is to minimize this redundancy measure.

2.2.4. mRMR Optimization Objective

The final objective of mRMR is to simultaneously maximize relevance and minimize redundancy, which can be formulated as:

$$max\left[\mathcal{R}(S,c)-\mathcal{D}\left(S\right)\right]$$

(4)

This ensures that the selected feature subset retains the most informative variables while reducing redundancy.

2.3. Prediction Model

Given the superior EML models in accident severity prediction, this study adopts them as the primary predictive approach. Specifically, five top-performing EML models—Random Forest, LightGBM, XGBoost, CatBoost, and AdaBoost—are selected, with the Decision Tree model serving as the baseline. A brief introduction to the principles and core formulas of each model is provided below.

2.3.1. Decision Tree

A Decision Tree is a tree-structured model for classification and regression [42], which recursively partitions the feature space to construct decision rules. Its core idea is to select the optimal splitting feature based on Information $Gain$ or $Gini$ Impurity.

Information $Gain$ measures the improvement in dataset purity before and after splitting, calculated as:

$$\mathrm{Information}\mathrm{Gain}=H\left(D\right)-\sum _{v\in \mathrm{Values}\left(A\right)}{\displaystyle \frac{|D_v|}{\left|D\right|}}H\left(D_v\right),$$

(5)

Where $H\left(D\right)$ is the entropy of dataset D, and $D_v$ is the subset where feature A takes the value v.

$GiniImpurity$ measures the impurity of a dataset, calculated as:

$$\mathrm{Gini}\left(D\right)=1-\sum _{i=1}^k{p}_i^2,$$

(6)

Where $p_i$ is the proportion of samples of the $i-th$ class in dataset D.

2.3.2. Random Forest

Random Forest is an ensemble learning method based on decision trees, designed to enhance generalization by constructing multiple decision trees and aggregating their predictions [43]. Its core principle involves building diverse decision trees using Bootstrap Sampling and Random Feature Selection. For classification tasks, Random Forest determines the final prediction through majority voting:

$$\widehat{y}=\mathrm{mode}\left(\{{\widehat{y}}_1,{\widehat{y}}_2,\dots ,{\widehat{y}}_T\}\right),$$

(7)

Where $\widehat{y}$ is the prediction of the $i-th$ decision tree, and T is the number of trees.

2.3.3. LightGBM

LightGBM is an efficient implementation of Gradient Boosting Trees that leverages histogram-based decision tree algorithms and Gradient-based One-Side Sampling (GOSS) to accelerate training [44]. The objective function of LightGBM consists of a loss function and a regularization term:

$$\mathcal{L}=\sum _{i=1}^{n}L(y_i,{\widehat{y}}_i)+\sum _{k=1}^K\Omega \left(f_k\right),$$

(8)

Where $\mathcal{L}$ is the loss function, and $\Omega \left(f_k\right)$ is the complexity regularization term for the k-th tree.

2.3.4. XGBoost

XGBoost is a high-performance algorithm based on Gradient Boosting Trees [45], which improves model accuracy by using a second-order Taylor expansion to approximate the objective function and introduces regularization terms to prevent overfitting. The objective function of XGBoost is:

$$\mathcal{L}=\sum _{i=1}^{n}L(y_i,{\widehat{y}}_i)+\sum _{k=1}^K\left(\gamma T_k+{\displaystyle \frac{1}{2}}\lambda {\parallel w_k\parallel }^2\right),$$

(9)

Where $T_k$ is the number of leaf nodes in the k-th tree, $w_k$ represents the vector of leaf node weights in the k-th tree, and $\gamma$ and $\lambda$ are regularization coefficients. The term $\parallel w_k{\parallel }^2$ is the L2 norm (squared) of the leaf node weights, serving as a regularization to penalize large weights and reduce model complexity, thus preventing overfitting.

2.3.5. CatBoost

CatBoost is a gradient boosting algorithm optimized for categorical features [46], which processes them using Target Statistics for encoding and Ordered Boosting to prevent target leakage. For a categorical feature x, CatBoost encodes it using Target Statistics as follows:

$${\widehat{x}}_i={\displaystyle \frac{{\sum }_{j<i}{\mathbb{I}}_{x_j=x_i}·y_j+a·p}{{\sum }_{j<i}{\mathbb{I}}_{x_j=x_i}+a}},$$

(10)

where a is a smoothing parameter to control the impact of prior knowledge, and p is the global mean target value. The term ${\mathbb{I}}_{x_j=x_i}$ is an indicator function that equals 1 if $x_j=x_i$, and 0 otherwise.

2.3.6. AdaBoost

AdaBoost is an ensemble learning method based on weighted voting [47], which iteratively adjusts the sample weights to train multiple weak classifiers and aggregates their predictions. In the t-th iteration, the sample weights are updated as follows:

$$w_i^{(t+1)}=w_i^{\left(t\right)}·exp\left(-{\alpha }_t{y}_i{h}_t\left(x_i\right)\right),$$

(11)

where $w_i^{(t+1)}$ is the weight of sample i in the $t+1$-th iteration, $w_i^{\left(t\right)}$ is the weight of sample i in the t-th iteration, ${\alpha }_t$ is the weight of the t-th weak classifier, $y_i$ is the true label of sample i (with values $\pm 1$), and $h_t\left(x_i\right)$ is the prediction of the t-th weak classifier for sample $x_i$. The weight ${\alpha }_t$ of the weak classifier is computed as

$${\alpha }_t={\displaystyle \frac{1}{2}}ln\left({\displaystyle \frac{1-{ϵ}_t}{{ϵ}_t}}\right),$$

(12)

where ${ϵ}_t$ represents the error rate of the weak classifier in the t-th iteration.

2.4. Model Interpretation and Feature Analysis

SHAP is a unified framework for interpreting machine learning models by assigning each feature a contribution value based on its impact on the model’s prediction [48]. It is rooted in cooperative game theory and Shapley values, ensuring fair and consistent feature attribution. SHAP offers several key advantages: Additivity: SHAP values decompose the model’s prediction into the sum of contributions from each feature, providing a clear and interpretable explanation. Consistency: If a feature’s impact on the model’s prediction increases, its SHAP value will also increase, ensuring reliable interpretations. Efficiency: For tree-based models, Tree SHAP computes SHAP values in polynomial time, making it significantly faster than traditional methods. Interpretability: SHAP values provide both local (individual predictions) and global (overall feature importance) interpretability.

2.4.1. Shapley Values

Shapley values distribute the contribution of each feature fairly by considering all possible feature subsets. For a model f and feature i, the Shapley value ${\varphi }_i$ is computed as:

$${\varphi }_i=\sum _{S\subseteq F\setminus \left\{i\right\}}{\displaystyle \frac{\left|S\right|!\left(\right|F|-|S|-1)!}{\left|F\right|!}}\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right),$$

(13)

where F is the set of all features, $\left|F\right|$ is the total number of features, S is a subset of features excluding feature i, $\left|S\right|$ is the number of features in subset S, $v\left(S\right)$ and $v(S\cup \{i\left\}\right)$ are the model outputs for feature subsets S and $S\cup \left\{i\right\}$, respectively, and the weight ${\displaystyle \frac{\left|S\right|!\left(\right|F|-|S|-1)!}{\left|F\right|!}}$ represents the probability of the feature subset S in all possible permutations.

2.4.2. Additivity Property

SHAP values satisfy the additivity property, meaning the model’s prediction can be decomposed into the sum of the SHAP values of all features plus a baseline value:

$$f\left(x\right)=E\left(f\left(x\right)\right)+\sum _{i=1}^M{\varphi }_i,$$

(14)

where $E\left(f\right(x\left)\right)$ is the baseline value (typically the mean prediction over the training data), and M is the total number of features. This property ensures that SHAP values provide a clear and interpretable explanation of the model’s output.

2.4.3. Consistency

SHAP values are consistent: if a feature’s impact on the model’s prediction increases, its SHAP value will also increase. This ensures reliable interpretations across different models and datasets.

2.4.4. Efficiency for Tree Models (Tree SHAP)

For tree-based models (e.g., decision trees, random forests, gradient boosting), Tree SHAP computes SHAP values efficiently in polynomial time by leveraging the tree structure. The SHAP value for feature i in a single tree is:

$${\varphi }_i=\sum _{v\in \mathrm{Paths}\left(i\right)}\left(\mathrm{PathProb}\left(v\right)·\left(\mathrm{Value}(v\cup \{i\left\}\right)-\mathrm{Value}\left(v\right)\right)\right),$$

(15)

where Paths(i) is the set of paths where feature i is used in a split, PathProb(v) is the probability of path v, and Value(v) is the prediction value at the end of path v.

Figure 3. Vehicle features in the accident dataset. (a) Diagram illustrating vehicle collision locations based on NHTSA classification [49]. (b) Distribution of vehicle brands (VEH_MKE) in the dataset. The number inside the circle represents the percentage of the corresponding brand in total accidents, while the number next to the brand name indicates its assigned label.

Figure 3. Vehicle features in the accident dataset. (a) Diagram illustrating vehicle collision locations based on NHTSA classification [49]. (b) Distribution of vehicle brands (VEH_MKE) in the dataset. The number inside the circle represents the percentage of the corresponding brand in total accidents, while the number next to the brand name indicates its assigned label.

3. Data Description

This study uses crash data from the CRSS and FARS databases maintained by the NHTSA from 2016 to 2022. The CRSS provides a nationally representative sample of police-reported crashes, including those resulting in property damage, injuries, or fatalities, while the FARS records fatal crashes where at least one fatality occurs within 30 days of the accident. The two databases collectively include 594,773 crash cases, covering single-vehicle, two-vehicle, and multi-vehicle collisions, with PV-PV collisions being the most frequent, accounting for 225,461 cases. A detailed breakdown of the raw data is presented in Table 1. More information on these databases can be found at NHTSA CRSS and NHTSA FARS.

3.1. Descriptive Statistics of Accident Data

To ensure data completeness and reliability, accident cases missing critical feature information were removed. Since this study focuses on predicting the injury severity of the subject vehicle’s driver, certain attributes of the other party in the collision were considered less relevant. For instance, the seat belt usage of the other driver is unlikely to significantly impact the subject driver’s injury severity. Consequently, some cases were retained despite incomplete information on the other driver, resulting in instances where only one driver in a collision was analyzed.

After preprocessing, the final dataset comprises 116,998 accident cases and 182,296 drivers, ensuring that each retained case contains sufficient information for injury severity prediction. A detailed summary of the dataset after preprocessing is presented in Table 2, Table 3 and Table 4.

After applying the mRMR algorithm, the feature set was reduced from 65 to 33, achieving a balance between relevance and redundancy. At this selection point, feature redundancy is minimized, and further increasing the number of features does not yield noticeable improvements in predictive accuracy. The retained features encompass key accident-related factors, ensuring both robust model performance and meaningful interpretability in accident analysis. The correlation among the selected features, assessed using Cramér’s V [50], is visualized in 6. This refinement enhances model performance by ensuring the inclusion of only the most relevant and non-redundant features, which are categorized into three groups:

1) Environmental Features (7 features): These features are related to weather conditions and accident location. Table 2 presents all 7 environmental features.

2) Collision and Pre-Collision Features (8 features): These features describe the driver’s evasive actions before impact and the collision dynamics. Table 3 lists six of these features, while the initial impact areas of both vehicles (IMP_ARE1 and IMP_ARE2) are shown in Figure 4. For better visualization, the definition of the initial impact area of a vehicle is illustrated in Figure 3(a).

3) Participant Features (18 features): These features include driver and vehicle-related characteristics for both parties involved in the collision. Table 4 presents 14 of these features, the age distributions of the two drivers (DRV_AGE1 and DRV_AGE2) are illustrated in Figure 5, and the vehicle brand statistics of the involved vehicles are shown in Figure 3 (b).

This selection process ensures that the retained features contribute meaningfully to accident severity prediction while minimizing redundancy and noise in the dataset.

3.2. Feature Statistical Description

Building on the statistical analysis of environmental factors, accidents most frequently occur during daytime, under clear weather, on straight and level roads, and on two-way undivided roadways outside intersections—scenarios where vehicles spend most of their travel time. While crash frequency is high due to greater exposure, severe injuries are less common, as good visibility, well-maintained pavement, and predictable traffic flow enhance safety. However, two-way undivided roads stand out as an exception, exhibiting a relatively higher fatality risk despite their prevalence, likely due to the increased likelihood of high-speed head-on collisions.

Shifting focus to collision dynamics, pre-collision features play a crucial role in accident severity. Braking upon detecting danger is associated with a lower risk of fatal injuries (8%), whereas vehicle instability before impact—such as skidding laterally—correlates with a substantially higher severity risk (70.9%), likely due to factors such as slippery roads, excessive speed, and driver errors. Accident type (ACC_TYPE) also influences injury outcomes: rear-end collisions, though the most frequent, are linked to a relatively lower fatal injury risk (4.4%), while head-on collisions present the highest risk due to greater collision energy. Impact location further affects injury severity, with side impacts at the 9, 10, and 3 o’clock positions exhibiting significantly higher fatal injury risks (43.5%, 38.1%, and 35.3%, respectively), likely due to limited energy absorption and intrusion into the driver’s space.

It is important to note that since this study integrates two datasets with different class distributions—one covering all accident types and the other focusing on high-severity cases—the reported percentages indicate relative risk levels rather than absolute probabilities. Being struck by the front end (12 o’clock position) of the opposing vehicle poses the highest risk (28%). Similarly, higher pre-crash speeds (TRV_SPD1, TRV_SPD2) are strongly correlated with increased injury severity, while rollovers and loss-of-control crashes—often involving multiple impacts with objects such as the ground, trees, or guardrails—further elevate the risk of severe injuries.

An analysis of driver and vehicle characteristics identifies additional factors influencing injury severity. Male drivers exhibit a higher risk of fatal injury than female drivers, with Cramér’s V [50] analysis revealing significant associations between gender and vehicle attributes, including VEH_AGE1, VEH_BDY1, VEH_MKE1, and TRV_SPD. Furthermore, impairment due to alcohol consumption or drug influence is a critical factor associated with increased injury severity.

Finally, protective measures and vehicle attributes play a key role in reducing injury severity. Seat belt use and airbag deployment significantly lower injury severity, with non-deployed airbags typically seen in low-energy collisions that result in minor injuries. Heavier vehicles generally offer better protection due to lower collision accelerations, with vehicle weight also linked to structural attributes like height, width, and length that affect crash safety. Additionally, collisions involving minivans or cargo vans tend to result in less severe injuries, possibly due to factors such as professional drivers and lower urban speeds.

4. Experimental Results

To predict accident severity, five EML models and a Decision Tree model were employed, treating the prediction as a multi-class classification task. Model performance was evaluated using multiple metrics, as shown in Figure 6. Since accuracy alone does not provide a complete assessment, particularly in imbalanced classification problems, precision, recall, and F1-score were also used to ensure a more comprehensive evaluation.

Despite merging the CRSS and FARS databases, class imbalance persisted, with severity levels distributed as follows: Level 0 (45.3%), Level 1 (26.9%), Level 2 (8.8%), and Level 3 (19%). To mitigate this, random oversampling and undersampling techniques were applied, using the average class size as a benchmark. Classes exceeding the benchmark were downsampled, while those below it were upsampled. For synthetic data generation, categorical features remained unchanged, while numerical features were modified by adding random noise within ±10% of the original values to enhance data variability while preserving statistical integrity.

Bayesian optimization was employed to determine the optimal hyperparameters for each model, offering advantages such as reduced computational cost and efficient exploration of the hyperparameter space. Model performance was evaluated both before and after data balancing, with hyperparameters tuned separately for each scenario. The dataset was split into 60% training, 20% validation, and 20% testing sets.

Following data balancing, all models demonstrated significant performance improvements. As shown in Table 5, prior to balancing, XGBoost achieved the highest accuracy (79.52%) and precision (71.8%), while CatBoost excelled in recall (68.6%) and F1-score (68.69%). After balancing, XGBoost outperformed all models across all four metrics (Accuracy = 84.9%, Precision = 84.85%, Recall = 84.9%, F1-score = 84.87%), followed by Random Forest (Accuracy = 83.74%, Precision = 83.63%, Recall = 83.77%, F1-score = 83.66%). However, Decision Tree exhibited a slight accuracy decrease of 0.73% after balancing, despite notable improvements in precision, recall, and F1-score. This suggests that the model initially favored predicting the majority class in the imbalanced dataset. The redistribution of class instances reduced this bias, enhancing its ability to classify minority classes correctly, while causing a slight accuracy drop.

4.1. Accident Feature Analysis

Based on the model performance evaluation, XGBoost demonstrated superior performance across all metrics and was therefore selected as the target model for interpretability analysis in this study. The SHAP method, based on Shapley values (as detailed in Section 2), was employed to systematically interpret its decision-making mechanisms.

Global SHAP analysis quantitatively assessed the contributions of individual features to predictions across the four accident severity classes. Figure 7(a-d) illustrates the SHAP value distributions for different severity categories. The horizontal axis represents the range of SHAP values, indicating feature influence intensity (positive values support the prediction of the corresponding severity class, while negative values indicate opposing effects). The vertical axis lists feature variables, with scatter point colors representing feature magnitudes (for numerical variables) or encoded categorical labels (for nominal variables)

Analytical results indicate that TRV_SPD1 and TRV_SPD2 make the greatest contribution to Level 0 predictions, suggesting that the travel speeds of both collision parties play a crucial role in determining accident severity by influencing kinetic energy transfer. Lower velocities (cool-toned colors) increase the likelihood of Level 0 predictions, whereas higher velocities (warm-toned colors) lower the probability of being classified as Level 0.

Further investigation reveals that most features exhibit nonlinear relationships between their values and SHAP contributions. Most features exhibit a bimodal distribution, where lower-value ranges correspond to predominantly negative SHAP values, transitioning to positive values at higher ranges. Additionally, certain features display polarity reversal across different severity classes. This phenomenon suggests that the model captures context-dependent variations in feature contributions. Notably, vehicle speed parameters are strongly correlated with accident severity, supporting the validity of the model’s decision-making process and its consistency with fundamental collision dynamics principles.

4.1.1. Feature Importance Analysis

Feature importance analysis is essential for understanding model decision-making and guiding accident analysis. While SHAP value distributions show how individual feature values impact predictions, they do not quantify the overall significance of features. To address this, the mean absolute SHAP value (Mean(|SHAP|)) is used to measure the average contribution of each feature across all predictions. The feature importance ranking results are shown in Figure 8, with the feature importance for another eight types of accidents presented in Figure A1. TRV_SPD1, TRV_SPD2, ACC_TYPE, and IMP_ARE2 are identified as the most influential features in the model. However, globally important features do not always align with those most critical for specific severity levels. For instance, the primary predictors for Level 3 severity are TRV_SPD2, DRV_AGE1, and SAF_ARBG.

Moreover, to assess the impact of feature importance on model performance, features were progressively removed in ascending order of importance, followed by model retraining and evaluation. The relationship between feature count and model accuracy, illustrated in Figure 9, demonstrates that retaining high-importance features results in significant performance improvements, while low-importance features contribute less substantially.

Although global rankings identify key predictors, they may overlook features with strong localized effects. For example, DRV_ALC1 ranks 22nd globally but exhibits widely dispersed SHAP values, indicating that despite its low frequency, it has a substantial impact on accident severity. This underscores the necessity of considering both overall and scenario-specific importance for a more comprehensive assessment of accident severity determinants.

4.1.2. Feature Dependency Analysis

In this study, feature dependency analysis is performed using the SHAP dependence plot, which captures the marginal effect of individual features on the target variable while illustrating the distribution and variability of SHAP values. Each point on the plot represents the predicted target value for a specific observation in the dataset. The analysis focuses on Level3 severity, representing fatal accidents—the most critical crash severity category. Thus, feature dependency analysis primarily examines feature behavior at this level to provide deeper insights into their contributions to fatal crash outcomes. Specifically, four key features—TRAV_SPD1, IMP_ARE1, VEH_MKE1, and VEH_BDY1—are analyzed. Additionally, the system autonomously identifies the feature with the strongest interaction within the dataset, as shown in Figure 10.

Figure 10(a) presents the interaction between TRV_SPD1 and DRV_SEX1. As TRV_SPD1 increases, the SHAP value shifts from negative to positive along the y-axis, with a transition threshold between 40 mph and 55 mph, indicating a critical speed range for predicting fatal accidents. The scatter point colors represent driver gender, revealing nearly identical SHAP distributions for male and female drivers. Although DRV_SEX1 exhibits a globally bimodal impact, its narrow SHAP distribution suggests a negligible overall influence compared to speed, whose broader distribution signifies a more substantial effect.

Figure 10(b) illustrates the interaction between IMP_ARE1 and IMP_ARE2, underscoring the strong association between impact location and fatal crash outcomes. Collisions at the 9 o’clock and 10 o’clock positions contribute most significantly to fatal injury predictions, aligning with real-world crash dynamics. These zones, located near the driver, absorb less impact energy, leading to severe structural deformation and reduced survival space. The effect is exacerbated when these areas are struck by another vehicle’s 11 o’clock or 12 o’clock positions, where impact forces are typically highest.

Figure 10(c) examines the influence of VEH_MKE1 and VEH_WGT1, demonstrating brand-specific variations in vehicle safety performance. The x-axis represents encoded vehicle brand labels (refer to Figure 4), allowing for a direct comparison of brand influence on accident severity. Scatter point colors denote vehicle curb weight, highlighting how certain manufacturers specialize in distinct vehicle categories. For example, Hummer (15) focuses on off-road vehicles, while RAM (10) primarily produces light trucks—both featuring high curb weights and reinforced body structures. Conversely, Volkswagen (13) and Tesla primarily manufacture passenger cars with lower curb weights but high safety ratings, aligning with their real-world safety reputations.

Figure 10(d) explores the relationship between VEH_BDY1 and VEH_WGT1, emphasizing the significant variation in crash survivability across vehicle body types. Light trucks exhibit lower susceptibility to fatal crashes compared to passenger vehicles (PVs), which are more prone to severe injuries. Scatter point colors indicate that light trucks generally have higher curb weights than PVs. Among body types, pickup trucks demonstrate the highest safety performance, whereas convertibles rank lowest, underscoring the critical role of vehicle structure and weight in crash outcomes.

4.1.3. Bivariate SHAP Contribution Analysis

Given that SHAP values exhibit additivity, the sum of individual feature SHAP values equals their cumulative contribution to the model’s prediction. This property allows for the evaluation of the combined impact of two or more features on the target variable by aggregating their SHAP values. In this study, Bivariate SHAP Contribution Analysis investigates the joint influence of feature pairs on model predictions by analyzing the sum of their SHAP values. This approach uncovers feature interactions and their collective impact on accident severity classification. By visualizing the aggregated SHAP contributions through a scatter heatmap, patterns emerge that show how specific feature combinations lead to higher or lower predicted severity levels, providing deeper insights into interdependent factors driving model decisions.

Figure 11 (a-d) displays the joint contribution of TRV_SPD1 & TRV_SPD2 to injury severity levels ("Level 0" to "Level 3"). The color of each scatter point reflects the magnitude of the SHAP value, illustrating how increasing vehicle speeds lead to higher predicted severity. Similarly, Figure 11 (e-h) depicts the combined effect of DRV_AGE1 & VEH_AGE1, showing that drivers over 80 years old are more likely to experience Level 3 injuries, with the age threshold decreasing as vehicle age increases, reflecting improvements in vehicle safety technologies. Figure 11 (i-l) demonstrates the joint influence of VEH_WGT1 & VEH_WGT2, revealing that heavier vehicles typically sustain lower injury severity compared to lighter vehicles in collisions.

This analysis effectively reveals the joint influence of feature pairs on accident severity outcomes. The clustering of colors within the heatmaps highlights distinct support, rejection, and transition regions. Red regions indicate the support domain, where predictions are more likely to align with the corresponding category, while deep blue regions represent the rejection domain, where predictions are less likely to belong to the specified category. Light green and light yellow regions correspond to transition zones, where predictions exhibit intermediate tendencies and neither strongly support nor reject the predicted outcomes. The clarity of these clusters underscores the contribution of these feature interactions in enhancing model prediction accuracy.

4.1.4. Instance-Level SHAP Interpretation

Instance-Level SHAP Interpretation decomposes the SHAP values for a single prediction instance, revealing the specific contribution of each feature to the final prediction and explaining why the model made that particular decision. Figure 12 presents the SHAP value waterfall plot for a randomly selected accident case from the dataset, where the actual injury severity is classified as "Level 3." The x-axis represents the SHAP values, while the y-axis lists the features in descending order of importance. The feature labeled ’Others’ represents the cumulative SHAP values of lower-contributing features, ensuring a more concise visualization.

The SHAP value waterfall plot provides an intuitive and transparent representation of the model’s decision-making process, illustrating how individual features contribute to the final prediction. The baseline value, $E\left(f\right(x)$, representing the expected model output across the dataset, is computed as 0.360, 0.568, 0.555, and 0.465 for injury severity levels 0 to 3, respectively. By aggregating the SHAP values of each feature relative to this baseline, the cumulative SHAP value for each severity level is obtained. The final predicted probabilities for each severity category are then derived using the Softmax function.

This instance-level interpretation enhances the explainability of the model by highlighting the key factors that influence the classification outcome. Such insights are particularly valuable for understanding the underlying risk factors in individual accident cases and can aid in data-driven decision-making for traffic safety analysis and policy development.

5. Implication and Limitations

5.1. Implications

This study focuses on driver injury severity prediction and factor analysis in PV-PV collisions, addressing a critical research gap in this domain. By leveraging interpretable machine learning techniques, we not only enhance prediction accuracy but also provide deeper insights into the key factors influencing injury severity, thereby informing traffic safety improvements and accident prevention strategies.

Moreover, by integrating the CRSS and FARS datasets, we have significantly mitigated the issue of class imbalance, reducing it by two orders of magnitude. This enhancement strengthens the model’s ability to predict high-severity injury cases, which are often underrepresented in imbalanced datasets. Given that accurately identifying severe injury cases is essential for emergency response planning and policy formulation, this improvement substantially increases the practical applicability of our findings.

Furthermore, this study employs SHAP to interpret the machine learning model, enabling a precise quantification of each feature’s impact, magnitude, and direction on injury severity. This interpretability not only enhances the model’s credibility but also provides actionable insights for stakeholders, such as policymakers and vehicle manufacturers, facilitating the development of targeted safety interventions.

5.2. Limitations

Despite the promising results, this study has certain limitations. Traffic accident analysis and injury severity prediction are inherently complex, and several critical factors cannot be accurately captured in structured datasets. For instance, driver psychological state and fatigue levels play a significant role in accident outcomes but are difficult to quantify. Additionally, the collision process itself is highly dynamic—secondary collisions triggered by the initial impact can significantly affect injury severity, yet these cascading effects are not well represented in traditional datasets. As a result, while our model provides valuable insights, its predictive capability remains limited by the absence of real-world complexity factors.

5.3. Future Work

Future research could enhance injury severity prediction by incorporating vision-based methods to automatically extract both macro-level (e.g., vehicle trajectories, collision impact forces) and micro-level (e.g., driver facial expressions, pre-crash behavior) features from accident videos. This approach would enable a more comprehensive understanding of the accident process, overcoming the limitations of relying solely on structured datasets that capture only a single moment in time. By integrating video-based accident reconstruction with machine learning, future studies could achieve a more dynamic and holistic analysis of crash events, which would lead to improved predictive performance and more effective safety interventions.

6. Conclusions

This study employed interpretable machine learning algorithms to predict driver injury severity and analyze accident characteristics in PV-PV collisions. The proposed prediction model achieved high accuracy, with superior precision, recall, and F1-score across all severity levels compared to existing studies. Notably, the model demonstrated significantly improved accuracy in predicting severe injury and fatality levels, which are the most critical categories in accident prediction and analysis. This improvement stems from the integration of the CRSS and FARS datasets, which mitigated class imbalance in high-severity cases, and the application of the mRMR algorithm, which reduced redundancy and refined feature selection.

By leveraging SHAP-based interpretable machine learning, this study provides a high-performance predictive model for injury analysis and accident causation studies in PV-PV collisions. Focusing on micro-level accident cases, the study identified the most influential features associated with accident occurrence and progression, ensuring a more precise and reliable prediction framework.

Global and local SHAP analyses revealed the contribution and directional impact of various accident-related features on injury severity prediction. Through multiple visualization techniques, the study examined both individual feature effects and combined feature interactions from different perspectives. These insights enhance the understanding of key accident factors, supporting more effective injury prevention strategies, risk assessment, and safety interventions. The proposed approach offers a robust framework for future research on accident severity prediction and interpretable machine learning applications in traffic safety analysis.

Appendix A.1

Figure A1. This is a figure. Schemes follow the same formatting.

Figure A1. This is a figure. Schemes follow the same formatting.