Driver injury severity analysis of work zone crashes: A Bayesian hierarchical generalized ordered probit approach

Abstract

Highway work zones are locations where severe traffic crashes tend to occur. Most of the extant research on work zone crash severity neglects the discrepancy in the injuries sustained by different drivers involved in the same crash. Admittedly, it is essential to analyse crash-level factors to their highest injury severity; but it is equally important to understand driver-level contributing factors to their injury severity to establish effective safety countermeasures to minimize drivers’ injury severity. Thus, this research aims to identify the factors with significant impacts on the driver injury severity of work zone crashes and estimate their effects on each severity level. Data on 3 880 drivers involved in 2 134 work zone crashes are obtained from the Crash Report Sampling System (CRSS) database of the United States and employed for the empirical investigation. A Bayesian hierarchical generalized ordered probit model is advocated for analysing the driver injury severity. Model performance indices suggest that the advocated hierarchical model is superior to the generalized ordered probit model, and considerable within-crash correlation is found across the observed driver injury severity. The estimated parameters show that driver age and sex, alcohol use, vehicle age and type, speeding and speed limit, weather conditions, lighting conditions and crash type have significant effects on the driver injury severity in work zone crashes. Marginal effects of the significant factors on each divided injury severity level are also estimated. Countermeasures are proposed from the results to reduce severe injuries sustained by drivers involved in work zone crashes.

Highlights

• The driver injury severity of work zone crashes is analysed.

• A Bayesian hierarchical generalized ordered probit model is proposed for the analysis.

• The effects of factors related to driver, vehicle and environment on driver injury severity are quantified.

• The proposed hierarchical model outperforms the generalized ordered probit model.

• Countermeasures for preventing driver injury in work zone crashes are provided.

1. Introduction

Highways need to be regularly maintained for good performance in operation and safety. For the maintenance, it is unavoidable to have to set up work zones occasionally on highways. Due to narrowed lanes, poor road condition, workers and complex arrangements of traffic control devices and signs, the risk of traffic crashes and severe injuries in work zones is usually higher than in other locations [1].

To develop effective countermeasures to mitigate the injury sustained by vehicle occupants, a number of studies have analysed factors contributing to the injury severity in work zone crashes [2]. Some of the studies covered all crash types in/near work zones [1, 3–7], while others focused on the injury severity of work zone crashes in particular collision types, such as rear-end [8–10], single-vehicle [11, 12] and truck-involved crashes [13–15]. Meanwhile, some studies concentrated on specific environments, including lighting conditions [16, 17], weather conditions [18] and work zone configurations [19, 20]. Nonetheless, most of the empirical investigations were at the crash level. That is, the injury severity of an observed crash, which is usually represented by the most severe injury sustained by occupants in the crash, is used as the dependent variable in a model. Although this approach is useful to find crash-level contributing factors, there are two main drawbacks. First, it neglects the differences in the injury severity of all driver-vehicle units/occupants within the same crash. Second, effects of the factors related to human and vehicles may not be precisely quantified, especially for multiple vehicle crashes. Therefore, exploring the driver injury severity at the driver/vehicle level provides a deeper insight into the factors contributing to crash casualties [21].

Analytical methods play a pivotal role in establishing a suitable relationship between crash injury severity and contributing factors [22, 23]. The methods for analysing crash injury severity can be divided into two groups: statistical models and machine-learning models. Although machine-learning methods, such as the support vector machine [4, 5], random forest [4, 6], classification and regression tree [16] and self-paced ensemble model [3], have been proposed for predicting the injury severity of work zone crashes, the black-box feature of these methods limits their ability to clearly interpret the effects of contributing factors on the crash severity. Therefore, statistical models (e.g. logit/probit model) are more prevalent, given their great interpretability. Specifically, Ghasemzadeh and Ahmed [18] used a binary probit regression for analysing the injury severity of work zone crashes, which is divided into two levels (i.e. property damage only and injury). For the work zone crashes with more than two severity outcomes, the multinomial logit model [1, 20], mixed logit model [13, 15, 17] and random parameters logit model with heterogeneity in means and variances [9, 11, 12] have been adopted for the analysis. The mixed/random parameters logit models can handle the unobserved heterogeneity which is an important issue in crash severity data [24]. The heterogeneity issue can also be addressed by finite mixture/latent class methods [25] or Markov switching approaches [26]. These models do not take the ordinal nature of the crash injury severity into consideration. To accommodate the ordinal nature, ordered logit/probit [7, 8] and generalized ordered logit/probit models [14, 19] have been employed. Compared with ordered logit/probit models, generalized ordered logit/probit models can remove the restriction imposed by the fixed thresholds on the effects of observed factors on crash severity [27], and thus yield a better model fit. Nonetheless, applying generalized ordered logit/probit models in analysing the driver injury severity may neglect the within-crash correlation, which denotes the potential correlation in the severity levels of injury sustained by the drivers who are involved in the same traffic crash. Omitting the within-crash correlation across driver severities may result in biased parameter estimates [23]. Huang and Aty [28] pointed out that the Bayesian hierarchical modelling framework can capture the within-crash correlation while improving model fit performance.

In light of the previous research, this study aims to analyse the driver injury severity of work zone crashes using a Bayesian hierarchical generalized ordered probit model. To validate its performance, the advocated hierarchical model is compared with a generalized ordered probit model in the context of Bayesian inference. The main contributions of this article are that: 1) the effects of factors related to driver, vehicle, roadway, environment and crash configuration on driver injury severity are explicitly revealed; 2) a Bayesian hierarchical generalized ordered probit model is proposed for examining crash injury severity, which can account for the ordinal nature and within-crash correlation simultaneously.

The remainder of the article is organized as follows. Section 2 introduces the work zone crash data for the empirical analysis. Section 3 presents the formulations of generalized ordered probit model and hierarchical generalized ordered probit model. The results of model comparison, estimated parameters and marginal effects are summarized and discussed in Section 4. Section 5 draws conclusions and provides suggestions for future research.

2. Data preparation

We collected three-year (from 2016 to 2018) crash data from the Crash Report Sampling System (CRSS), which is managed by the National Highway Traffic Safety Administration (NHTSA) of the USA. In the CRSS, all crash records are derived from law enforcement agencies and sampled nationwide. Excluding crash records with incomplete information, the data on 2 134 work zone crashes, where 3 880 drivers are involved, are extracted and used for the empirical analysis. In the original crash records, driver injury severity levels are divided into four: no apparent injury, minor injury, serious injury and fatality. Among the 3 880 drivers involved in work zone crashes, 2 939 had no apparent injury, 765 had minor injuries, 143 had serious injuries and 33 had fatal injuries. Due to the rarity of fatal injuries in the crash dataset, we combined fatal injury with serious injury to constitute the most severe injury level in the analysis.

For each observed driver, information on driver characteristics (e.g. age, sex, alcohol use), vehicle characteristics (e.g. age, type, speeding), roadway characteristics (e.g. horizontal and vertical alignments, speed limit), environmental characteristics (e.g. weather conditions, lighting conditions, rural/urban), presence of traffic control and collision type is extracted from the CRSS. These characteristics are categorized based on their original definitions, transportation engineering experience and previous studies [9, 29]. They are used as covariates in the models of driver injury severity. Table 1 displays the description and descriptive statistics of the response variable and covariates.

3. Modelling framework

We advocate a Bayesian hierarchical generalized ordered probit model for examining the driver injury severity of work zone crashes. To demonstrate its advantages, the advocated model is compared with its counterpart generalized ordered probit model. In this section, the structures of the two models are clearly described first (Section 3.1). The processes of Bayesian estimation and the criterion for model comparison are then introduced (Section 3.2). Finally, the calculation method for the maginal effects of covariates on driver injury severity is presented (Section 3.3).

3.1. Model formulation

3.1.1. Generalized ordered probit model

The generalized ordered probit model is a method to account for the ordinal nature of divided injury severity levels while relaxing the restriction imposed by the fixed thresholds in standard ordered logit/probit models [27]. For the |$j$|th driver involved in work zone crash i, the covariates X_i,j are assumed to be linearly linked to a latent variable U_i,j, which is usually called injury propensity:

where β is a vector of parameters (including a constant term) to be estimated, and ε_i,j is a residual term which follows a standard normal distribution.

The latent injury propensity U_i,j is matched to the recorded injury severity degree Y_i,j of the jth driver in work zone crash i, based on the driver-specific thresholds |$\mu _{i,j,k}$| (k=1, 2, …, K−1):

where 1, 2, …, k, …, K respectively denote the ordered injury severity levels which are numbered ascendingly (i.e. 1 represents ‘no apparent injury’, 2 represents ‘minor injury’ and 3 represents ‘serious injury or fatality’ in the current analysis). To improve the flexibility in quantifying the effects of the covariates while maintaining their order relationships (⁠|${\mu _{i,j,1}} < {\mu _{i,j,2}} < \cdots < {\mu _{i,j,k}} < \cdots < {\mu _{i,j,K - 1}}$|⁠), the variable thresholds are formulated as suggested by Eluru et al. [27]:

where Z_i,j,l is a vector of covariates related to the (l+1)th threshold; and α_l is a vector of estimable parameters (including a constant term). For the uniqueness of identification, the first threshold |$\mu _{i,j,}$|₁ is set as zero for each observation.

Based on the above assumptions, the cumulative probability of the jth driver in work zone crash i sustaining injury in the kth (k = 1, 2, …, K) severity degree, P_i,j,k, is computed as:

where Φ(·) is the cumulative distribution function of the standard normal distribution.

As a consequence, the probability of the jth driver in crash i sustaining the kth level injury, p_i,j,K, is calculated as:

3.1.2. Hierarchical generalized ordered probit model

A work zone crash often involves multiple driver-vehicle units. There may be some unobserved factors with common effects on the injury severities of drivers involved in the same crash, which result in the within-crash correlation. Neglecting the within-crash correlation in driver injury severity may lead to inaccurate, biased parameter estimations, or both [30]. To jointly account for the ordinal nature and within-crash correlation, a hierarchical generalized ordered probit model is developed for the analysis. Specifically, on the basis of the generalized ordered probit model, a crash-specific random term δ_i,l is added into the formulation of the threshold |$\mu _{i,j,}$|_l+1. That is, the Eq. (3) is modified as:

where δ_i,l is assumed to be normally distributed:

where σ_l (>0) is the standard deviation.

3.2. Bayesian estimation and comparison criterion

3.2.1. Bayesian estimation

As the dominate estimation method for hierarchical models, Bayesian inference is employed to estimate the parameters/hyper-parameters in the above models. It is required to specify a prior distribution for each parameter/hyper-parameter, which reflects the prior information on it [31, 32]. In the absence of useful prior knowledge, a diffused normal distribution, N(0, 10⁴), is used as the prior of each element in β and α_l; and a diffused gamma distribution, G(0.001,0.001), is employed as the prior of |$1/{\rm{\sigma }}_l^2$| (i.e. the precision of δ_i,l).

We conduct the Bayesian estimation of the models in the freeware WinBUGS [33]. For each model, we set a chain of 100 000 Markov chain Monte Carlo (MCMC) simulation iterations, where the first 50 000 iterations are removed as burn-in, while the remaining 50 000 iterations are employed to infer the posterior distributions of the parameters of interest. The convergence of the MCMC simulations is empirically determined by monitoring if the Monte Carlo simulation error of each parameter is less than 5% of its posterior standard deviation and visually inspecting the trace plots for the parameters.

3.2.2. Model comparison criterion

The deviance information criterion (DIC), which provides a hybrid measure of model fitting and complexity, is used for comparing the comprehensive performance of the estimated generalized ordered probit and hierarchical generalized ordered probit models. DIC is the most widely used criterion for Bayesian model comparison and can be directly obtained after the model estimation in WinBUGS. According to the definition in Ref. [34], it is computed as:

where |$\bar{D}$| represents the posterior mean deviance used to measure the model goodness of fit, and pD represents the effective number of parameters used to measure the model complexity. In general, we prefer a model with a lower DIC value. Following the suggestion of Spiegelhalter et al. [35], we can rule out the model which yields a higher DIC value with a difference over 10.

3.3. Marginal effect

From the estimation results of the coefficients (i.e. the elements in β and α_l) in the above models, we can directly learn whether a certain covariate is significantly associated with the driver injury severity of work zone crashes; however, the direction or magnitude of its effect on the probability of each severity level is unclear. Thus, the marginal effect of each covariate with significant effect is calculated. Specifically, for a continuous covariate, its marginal effect on an injury severity level is calculated according to the first-order derivative with respect to it. For more details on the calculation equations, see Ref. [36]. For a dummy covariate x_m, its marginal effect is calculated as the changed probability when the covariate varies from 0 to 1:

The calculations of the marginal effects are specific to the individual driver. For the whole dataset, each significant covariate's average marginal effects for all drivers involved in the observed work zone crashes are estimated.

4. Result analysis

4.1. Model comparison

Table 2 displays the results of Bayesian estimation and assessment for the two models, where only the covariates with significant effects at the 90% credibility level are included. Based on the results, we found that: 1) the hierarchical generalized ordered probit model significantly outperforms the generalized ordered probit model in term of goodness of fit, since the hierarchical model yields a lower |$\bar{D}$| value (4 630 for the hierarchical generalized ordered probit model versus 4 919 for the generalized ordered probit model); 2) the generalized ordered probit model is more parsimonious since there are fewer effective parameters in it, as indicated by the pD value (249 for the hierarchical generalized ordered probit model versus 52 for the generalized ordered probit model); 3) the overall performance of the hierarchical generalized ordered probit model is much better than that of the generalized ordered probit model, since the hierarchical model yields a lower DIC value and its difference to that of the generalized ordered probit model is 92 (4 879 for the hierarchical generalized ordered probit model versus 4 971 for the generalized ordered probit model). The above findings are generally consistent with those in the literature [28, 30, 37]: capturing within-crash correlation between driver injury severities using Bayesian hierarchical models can substantially improve model fit. In addition, considerable within-crash correlation is demonstrated by the Bayesian estimate of σ₁, which is statistically significant at the 95% credibility level.

According to the Bayesian parameter estimates, we revealed that there are certain differences in the significant covariates contributing to the latent injury propensity and the threshold between minor injury and serious injury/fatality in the two models. Specifically, ‘cloudy’ is significantly associated with the latent injury propensity in the hierarchical model, but its effect on the latent injury propensity or the threshold is not significant in the generalized ordered probit model. The hierarchical model suggests that ‘driver age’ >60 has significant effects on both of the latent injury propensity and the threshold, while the generalized ordered probit model suggests that it has a significant effect on the latent injury propensity only. The differences on the parameter estimates bring about the discrepancies on the marginal effects of these covariates in the two models, which are shown in Table 3. The extant literature [28, 30] concluded that accommodating the within-crash correlation across driver injury severity in the hierarchical models is helpful to reduce bias in the parameter estimation.

4.2. Interpretation of parameter estimates and marginal effects

Our empirical analysis of driver injury severity of work zone crashes is based on the estimation results in the hierarchical generalized ordered probit model, since it significantly outperforms its counterpart—the generalized ordered probit model.

According to the estimation results in Table 2, female drivers are associated with a higher injury propensity and a larger threshold between minor injury and serious injury/fatality as compared to male drivers (the reference category). Specifically, female drivers are 8.0% and 0.6% more likely to sustain minor injury and serious injury/fatality in work zone crashes, respectively. The results are generally consistent with the findings in Refs. [9, 20]. The authors argued that the effect of sex on crash injury severity may be attributed to the combination of physiological and behavioural factors which are considerably different between female and male drivers.

Driver age has a significant effect on the injuries sustained in work zone crashes. According to the results in Table 3, relative to young (<25) drivers (the reference category), middle-aged (25−60) and old (>60) drivers are 4.2% and 3.4% more likely to experience minor injury crashes, respectively; they are 1.4% and 4.7% more likely to experience serious injury/fatality crashes, respectively. The results are reasonable, because the vision of aged drivers is usually poorer, their reaction time to emergency is longer and their muscle strength is weaker [30, 38].

The significantly positive effect of ‘alcohol use’ on the latent injury propensity suggests that drunk drivers are prone to sustaining severe injury in work zone crashes. Specifically, the likelihood of resulting in minor injury and serious injury/fatality will increase by 13.3% and 5.7%, respectively, for drivers with alcohol use. The results are in line with engineering intuitions and the previous findings [9, 39]: under the influence of alcohol, drivers have more difficult in manipulating vehicles, perceiving potential dangers and taking proper actions to avoid oncoming collisions.

The negative posterior mean of the parameter for ‘truck’ in the link function of latent injury propensity suggests that truck drivers are less likely to be injured in work zone crashes than automobile drivers (the reference vehicle type). Specifically, the likelihoods of minor injury and serious injury/fatality sustained by truck drivers are 7.5% and 2.4% lower than those of automobile drivers, respectively. The results are attributable to the greater mass and structural rigidity of trucks, which provides stronger capacity to protect their occupants, i.e. higher crashworthiness [37, 40].

There is a significant association between driver injury severity and vehicle age. Relative to new (age <3) vehicles, middle-aged (3−13) and old (age>13) vehicles increase the likelihood of resulting in serious injury/fatality by 1.7% and 3.1%, respectively. As implied by a number of previous studies [18, 37, 41], the worse physical performance and the absence of advanced safety equipment of older vehicles may be responsible for the increased injury/fatality likelihood of their drivers.

The positive and negative effects of ‘speeding’ on the latent injury propensity and the threshold between minor injury and serious injury/fatality, respectively, suggest that vehicle speeding increases the severity level of driver injury. According to the results in Table 3, if the speed of a vehicle exceeds the posted limit when it collides with other vehicle(s) or fixed object(s), the likelihoods of minor injury and serious injury/fatality sustained by its driver will increase by 1.5% and 4.1%, respectively. The results are consistent with engineering intuition and provide clear evidence for the speeding forbiddance in work zones.

In addition to vehicle speeding, speed limit also has a significant effect on driver injury severity. Specifically, compared to a speed limit of less than 30 mph (48 km/h), a speed limit of 30 mph−60 mph (48 km/h−96 km/h) or of over 60 mph (96 km/h) increases the likelihood of minor injury by 4.8% and 10.4, respectively, while increasing the likelihood of serious injury/fatality by 1.7% and 3.9%, respectively. The results are reasonable, because the vehicle speed is usually higher on roadways with higher speed limits. Zeng et al. [41] found that higher speed is linked to not only more severe driver injury of their own vehicle but also more severe driver injury of the other vehicle(s) in the same crash.

With regard to weather conditions, ‘cloudy’ has a significantly negative effect on the threshold. The likelihood of serious injury/fatality incurred by drivers involved in work zone crashes in cloudy weather is 1.4% higher than in clear weather (the reference category). The ‘rain/snow’ has a significantly negative effect on the latent injury propensity, which indicates that drivers are less likely to be severely injured in rain/snow weather. Its estimated marginal effects show that the likelihood of serious injury/fatality is 1.3% lower than in clear weather. While it is somewhat counterintuitive, some previous studies [36, 42] have found similar findings and concluded that drivers may be more cautious and slow down in inclement weather conditions (such as heavy rain, snow), which is helpful to minimize the crash injury severity.

The significantly positive effect of ‘dark’ on the latent injury propensity suggests that driver injury severity tends to be higher in a dark environment than in daylight (the reference category). Specifically, the likelihoods of resulting in minor injury or serious injury/fatality in a dark environment are 13.3% and 5.7% higher than in daylight, respectively. The results are probably due to the worse vision of drivers in a dark environment, which means less time for them to properly respond to any hazards [12].

Regarding the variables related to collision type, ‘rear end’ and ‘sideswipe’ have negative effects on the latent injury propensity, while ‘rear end’ and ‘angle’ have positive effects on the threshold. The marginal effects of these variables manifest that the likelihood of serious injury/fatality caused by a rear-end crash, angle crash and sideswipe collision are 5.1%, 1.9% and 3.7% lower than that caused by a fixed-object collision (the reference category), respectively. The results are generally in line with the existing findings [37, 43].

5. Conclusions

This research identified the factors contributing to driver injury severity in work zone crashes using the Bayesian hierarchical generalized ordered probit model, which is capable of accounting for the ordinal nature and within-correlation in the crash data. A three-year work zone crash dataset was collected from the CRSS for the empirical analysis and model demonstration.

The proposed Bayesian hierarchical generalized ordered probit model provided a substantially better performance than the counterpart generalized ordered probit model. The Bayesian parameter estimates revealed the significant with-correlation across the injury severity of drivers within the same work zone crash. These findings manifest the reasonableness of the advocated Bayesian hierarchical model.

The Bayesian estimation results showed that in work zone crashes: 1) drivers who are female, older or under the influence of alcohol are more likely to sustain severe injury; 2) automobiles, older vehicles and speeding have a increased severe injury risk for drivers; 3) driver injury is more severe occurring on roadways with a higher speed limit, in dark conditions and in cloudy weather; 4) driver injury in rear-end, angle and sideswipe collisions has a lower propensity to be severe than that in fixed-object crashes.

The findings from the study can provide practical implications in developing countermeasures aimed to mitigate the injury severity of work zone crashes. For example, targeted training and education efforts towards female and older drivers to regulate their behaviours when encountering an emergency is helpful to decrease their injury severity. Traffic police should improve the enforcement strategies against illegal driving behaviours, such as drunk driving and speeding near work zones. Vehicle manufacturers should put more efforts into enhancing the crashworthiness of passenger cars. Implementing some emerging technologies, such as intelligent connected vehicles and 5G communication, can improve work zone safety in adverse weather and lighting conditions.

Lastly, it is noteworthy to mention the limitations of the current study, which should be addressed in follow-up studies. First, the work zone characteristics (such as the type and location of work zones) are not included in the collected crash data. A more comprehensive dataset is needed to reveal their effects on the driver injury severity of work zone crashes. Second, exploring the differences in driver injury severity between single-vehicle and multi-vehicle crashes might also be an interesting topic for further investigation, especially for multi-vehicle crashes where the interaction across the driver-vehicle units is an important issue to be considered [41].

Acknowledgements

The work was supported by the Open Fund of the Key Laboratory of Highway Engineering of Ministry of Education (Changsha University of Science & Technology) (Grant No. kfj230301).

Conflict of interest statement

None declared.

References