Applying a joint model of crash count and crash severity to identify road segments with high risk of fatal and serious injury crashes

is an important, unanswered question. While a joint econometric model of crash count and crash severity has the ﬂ exibility to account for the limitations mentioned previously, its ability to identify high-risk sites also needs to be examined. This study aims to ﬁ ll this research gap by employing the joint model for blackspot identi ﬁ cation. Using data from state-con- trolled roads in Queensland, Australia, a new risk score is developed based on predicted crash counts by severity, weighted by the cost ratio of severity levels. This weighted risk score is then used for identifying road segments with high risk of fatal and injury crashes. Results show that the joint model of crash count and crash severity has substantially improved prediction accuracy compared to the traditional count models. The correlation between crash counts of di ﬀ erent severity levels captures the unobserved heterogeneity caused by the extra-variation in total crash counts and moderates the parameters in the joint model. In comparison with the traditional ap- proaches, the proposed weighted risk score approach with the joint model of crash count and crash severity leads to the identi ﬁ cation of a higher number of fatal and serious injury crashes in the top ranked sites ﬂ agged for safety improvements.


Introduction
Crash count and crash severity have been widely used in transport safety as two indicators of crash risk along road segments (Washington et al., 2018). Crash count is a measure of the likelihood of crash occurrence, whereas crash severity is a measure of the societal impact of crashes and their harm to the society. Although federal road agencies and departments of transport aim to reduce crash counts across their road networks, safety improvement programs are primarily focused on preventing fatal and serious injury crashes, as the cost per person of fatal and serious injury crashes are substantially higher than minor injury and non-injury crashes (Harmon et al., 2018). In addition, the collective social cost of crashes in the society is substantially higher than the individual social costs (Tay, 2002). As a result, considering crash severity in conjunction with crash count is paramount for identification of high-risk sites (also referred to as blackspots or hotspots).
The several studies have incorporated crash severity into the identification of high-risk sites (which will be comprehensively reviewed in the next section) can be divided into three groups based on their methodological approaches: multivariate crash count models, equivalent property damage only (EPDO) models, and two-stage mixed models of crash count and crash severity (Washington et al., 2018). All of these approaches, perhaps not surprisingly, have some methodological/empirical limitations. The multivariate modelling of crash counts has been mostly achieved by assuming a Poisson lognormal distribution for crashes of different severity levels because the resulting count model has a hierarchical representation in the Bayesian platform (Shaon et al., 2019). Other distributions such as Poisson-Gamma (or Negative Binomial) make the multivariate crash count models computationally burdensome. More importantly, multivariate crash count models neglect the ordered nature of crash severity. The EPDO approach requires nonparametric models (such as quantile regression) because EPDO data are not distributed according to well-known statistical distributions. As such, additional complexity in crash severity data, for example unobserved heterogeneity, is not easily accounted for in this approach. The two-stage mixed modelling of crash count and crash severity again fails to account for the possible correlation between crash counts of different severity levels. These shortcomings render the existing approaches suitable for certain, 'well behaved' conditions only, and perhaps leave researchers to develop a rigorous methodology that can capture additional complexity combining crash frequency and injury severities.
Recently, an alternative approach was proposed by Yasmin and Eluru (2018) to explicitly incorporate the severity of crashes into crash count prediction. Unlike the two-stage mixed models, this approach is based on a joint (one-stage) modelling of crash count and crash severity, whereby a crash count model is jointly estimated with a crash severity model, and a correlation term is added to the overall model to account for the common factors between the two model components. The overall joint model has two important properties: 1) the parameters of the crash count model are influenced by the parameters of the crash severity model in a joint estimation process, and 2) the overall model parameters are moderated by the correlation between the two model components. Yasmin and Eluru (2018) showed that the joint model of crash count and crash severity is superior to traditional crash count models in terms of statistical fit. However, the performance of this model was not evaluated in identification of high-risk sites mainly because a blackspot identification criterion for this model does not exist to date, and thus, its potential for blackspot identification remains unexplored. This study aims to bridge this gap by investigating the applicability and the performance of the recently developed joint model of crash count and crash severity in identifying highway segments with high risk of fatal and injury crashes.

literature review
Multivariate modelling of crash counts is one of the most popular approaches in the literature to consider crash severity while identifying high-risk sites. Aguero-Valverde and Jovanis (2009) applied a full Bayesian multivariate Poisson lognormal regression to model crashes in Pennsylvania, United States and found a high correlation (from 0.47 to 0.97) between crash counts of different severity levels. They used the costs of different crash severity levels to convert the outcome of the multivariate model to expected crash cost and compared the expected crash cost with the observed crash cost to rank highway segments. They found that high-risk segments identified by the outcome of the multivariate model have consistently higher excess costs than those segments identified by the univariate model. The assumption of a lognormal distribution for the correlation between the crash counts of different severity levels in the multivariate setting resulted in crash counts following a Poisson lognormal distribution. This assumption, however, may not be entirely accurate as crash counts of different severity levels (fatal crashes with excess zero counts for example) may not be Poissonlognormally distributed (Aguero-Valverde, 2013). In addition, the multivariate model does not consider the ordinal nature of crash severity levels and thus may lead to incorrect inferences about the effects of explanatory variables and result in misidentification of blackspots and inefficient allocation of resources.
EPDO is another popular approach to incorporate crash severity into the selection of high-risk sites. Washington et al. (2014) converted the crash counts of different severity levels to an equivalent property damage only count by applying weights based on average cost ratios by severity. They applied a non-parametric quantile regression on EPDO data in order to estimate the effects of covariates on various quantiles of the population, rather than the population mean. Finally, they compared the outcome of the quantile regression model with the observed EPDO and ranked sites accordingly. They found that the high-risk sites identified by the EPDO approach places more emphasis on more harmful crashes than does the conventional count approach. As EPDO data are not distributed like any well-known distributions, non-parametric quantile regression models must be used in this approach. This requirement is a limitation of this EPDO approach as there is no likelihood function to be maximized in this approach and thus the complexity of crash data (e.g. unobserved heterogeneity) cannot be readily incorporated into the modelling process. Recently, some studies have used parametric regression models (e.g. lognormal hurdle and Tobit regression) to model EPDO rate as a function of exogenous covariates (Ma et al., 2016); however, the appropriateness of these distributional assumptions has not been tested.
Finally, two-stage mixed modelling of crash count and crash severity is another approach to consider crash severity for identifying high-risk sites. Miranda-Moreno et al. (2009) first introduced this approach and defined a total risk index as a multiplicative factor of crash count at a site and its expected consequences including expected number of fatal, serious, and minor injuries. They employed Bayesian framework to estimate a hierarchical Poisson model of crash count and a multinomial model of crash severity, independently. They constructed a total risk score per segment and compared it with a standard value established by decision makers and ranked highway segments accordingly. They found that the two-stage mixed model were able to identify high-risk highway segments based on crash count and crash severity and thus improved the effectiveness of allocating resources for safety improvements. In more recent studies, Wang et al. (2011) andStipancic et al. (2019) used a similar two-stage mixed modelling approach to identify high-risk sites based on crash count and crash severity and found that this approach is superior to the traditional count based approaches. However, the two-stage (independent) modelling of crash count and crash severity in all of these studies crucially ignores the possible correlation between crash counts of different severity levels. This ignorance may result in biased parameters and in turn may lead to incorrect predictions of crash count and crash severity.
Recently, a joint model of crash count and crash severity has been proposed by Yasmin and Eluru (2018) which is more flexible than the existing approaches in that it does not require the methodological/ empirical assumptions of the above approaches. In particular, this approach is not constrained to any distributional assumptions for the extra-variation in crash counts over and above that accounted for by the Poisson density. The joint model is estimated using maximum likelihood estimation methods and so it can accommodate the complexities associated with crash data, e.g., unobserved heterogeneities. The crash count and crash severity model components are jointly estimated (in one stage) and the parameters of the joint model are moderated by the correlation between crash count and crash severity. As such, the joint model of crash count and crash severity represents a promising alternative to the existing approaches, but importantly, requires further exploration regarding its ability for blackspot identification.

Methodology
This section presents the structure of the joint model of crash count and crash severity which was first introduced by Yasmin and Eluru (2018). To describe the model structure more elaborately, the crash count and crash severity model components are presented separately. The joint model of crash count and crash severity is then presented, followed by the selection criteria for identifying highway segments with high risk of fatal and injury crashes. Let Y it represent the total observed crash count on the ith highway segment (i = 1, 2, 3, …, N) and in the tth year (t = 1, 2, 3, …, T).

Crash count model
The transport safety literature has shown that Y it follows a negative binomial (NB) distribution with mean μ it and inverse dispersion parameter φ: Assuming an exponential function for the mean of the negative binomial distribution, the total predicted crash count (μ it ) on the ith segment and in the tth year can be expressed as a function of exogenous explanatory variables (Afghari et al., 2018a): where X it are other explanatory variables and β i are estimated regression parameters (including the intercept) and ε exp( ) it is a random error term, which follows a Gamma distribution with mean 1 and variance φ 1 . To account for unobserved heterogeneity, model parameters (β i ) are allowed to vary across highway segments (Anastasopoulos and Mannering, 2009). Note that parameters are fixed across time to account for multiple observations on a segment during different time periods (i.e. panel data setting). Such a model specification is referred to as grouped random parameters (Sarwar et al., 2017) in which model parameters are assumed to follow probabilistic distributions (e.g. Normal, Uniform, Triangular, etc.) across the observations within each group (or panel). The probability density function of the RPNB model is: where Γ (.) is the gamma function and f β ( ) is the density function of the model parameters. The log-likelihood function (LL) of the model is obtained by integrating the probability density function of the model over the entire set of random parameters, applying the logarithm transformation and summing it over observations to yield: To estimate this complex log-likelihood function, Maximum Simulated Likelihood Estimation is used where quasi random draws from Halton sequences are employed to simulate the densities of the random parameters (Bhat, 2001). It has been shown that this simulated maximum likelihood estimator is unbiased and consistent for a large number of draws (Munkin and Trivedi, 1999).

Crash severity model
To incorporate crash severity into the models, let s (s = 1, 2, 3, …, S) represent crash severity categories (e.g. property damage only crashes, minor injury crashes, serious injury crashes, fatal crashes, etc.) on highway segments. In ordered models, the actual proportion of total crashes by severity levels (Y sit ) is associated with an underlying latent variable (Y * it ). This latent variable is then mapped to the actual severity proportions by thresholds (τ ) and using the following linear function: where κ is the vector of parameters, X it is the vector of covariates and δ i is the random error term.
To estimate the latent propensity of crash severities, it is assumed that: where H (.) sit is the probability density function for the severity category s. Depending on the distributional assumption for the probability of error terms, H (.) sit can take standard normal or standard logistic probability density functions for the ordered probit or ordered logit models, respectively. The latter functional form is used in this study to construct an ordered logit model for crash severity. The probability of each crash severity category is then presented as: where φ (.) is the standard logistic cumulative probability density function. The corresponding quasi log-likelihood function is then expressed as: where w sit is the fraction (proportion between 0 and 1) of crashes in severity category s at road segment i and time period t, and the rest of notations are as previously stated. These fractions sum to unity over the ). This model is referred to as fractional split (Afghari et al., 2018b). Note that w sit takes binary values (0 or 1) in conventional logit models; one for the chosen alternative and zero for the non-chosen alternative. Maximum simulated likelihood approach is used to estimate this log-likelihood function.

Joint model of crash count and crash severity
To generate the correlation between the above mentioned ordered severity model and the previously described count model, a correlation term (υ it ) is now defined as: where γ i is the vector of parameters and m i is the vector of covariates capturing the observed correlation between crash count and crash severity components. This correlation term is then added to the propensity function of the severity model and the predicted mean of the count model: The log-likelihood function of the overall joint model of crash count and severity is then expressed as: where Ω represents the vector of random parameters in the joint model (β and κ) the and the rest of the notations are as previously stated.
Maximum simulated likelihood approach is used to estimate the joint econometric model.

Weighted risk score to combine crash count and crash severity
The estimated crash count and crash severity models can be used to derive the crash count predictions of different severity levels on road segments by multiplying the total crash count and the proportions of each severity level: A new weighted risk score (WRS) is now developed based on the cost ratio of crash severity levels: where WRS i is the predicted weighted risk score for segment i and cr s is the ratio of the cost of a crash in severity level s to the cost of a reference crash severity level. This weighted risk score is analogous to the property damage only equivalency factor in the EPDO approach (Washington et al., 2014) but the weights are applied to the predicted crash counts post estimation and thus the distribution of observed crash counts are intact. This weighted risk score can now be used in ranking highway segments based on their crash count and crash severity.
3.5. Selection criteria: potential for improvement and excess weighted risk score Two selection criteria are used in this study for identifying high-risk sites: potential for improvement (PFI) following the Empirical Bayesian (EB) method with the independent crash count model and excess weighted risk score for the joint model of crash count and crash severity. The EB approach has widely been used in the literature to identify high-risk sites (Afghari, 2019). The EB approach combines the predicted and observed number of crashes and accounts for the regression to the mean effect (Cheng and Washington, 2005;Washington et al., 2018). Many studies have shown that the EB approach is superior to alternative approaches in identifying high-risk sites (Montella, 2010). The EB estimator is a weighted sum of the predicted and the observed crash counts such that: where k is the over-dispersion parameter estimated during the SPF calibration process, and w 1 and w 2 are weights calculated based on the mean and variance of the SPF estimate (Persaud et al., 2010). Potential for improvement (PFI) has been used as the selection criteria for identifying high-risk sites based on their EB estimate (Washington et al., 2018). PFI is defined as the difference between the EB estimate and the predicted mean of crash count at a site: Highway segments are ranked according to decreasing PFI and segments with higher PFI are identified as high-risk sites.
Despite its appealing properties, the EB approach rests on distributional assumptions about crash occurrence (Poisson-Gamma or negative Binomial) which do not hold for the weighted risk score approach in this study. As such, excess weighted risk score is used as the selection criteria for ranking highway segments based on crash count and crash severity. The excess weighted risk score (EWRS) is defined as: where WRS it and WRS it are observed and predicted weighted risk scores. Segments are ranked according to the decreasing EWRS it and segments with higher EWRS it are identified as high-risk sites.

Empirical data
The data used to compare the performance of candidate models were collected for a random sample of highway segments and major arterial road sections (rural and urban) along the state-controlled roads in Queensland, Australia. The extent of the studied network is 1,477 km consisting of 521 road segments. The dataset includes crash data as well as roadway geometric and traffic characteristics data. Four years of crashes (from 2010 to 2013) were collected in three severity categories: fatal, serious injury and minor injury. The Queensland Department of Transport and Main Roads stopped collecting property damage only crash data in 2010 and so crashes in this category are not included in this study. Descriptive statistics of the crash data are presented in Table 1.
Roadway geometric and traffic operational characteristics include Annual Average Daily Traffic (AADT), percentage of heavy vehicle traffic, segment length, number of lanes, lane width, functional classification of the road (urban\rural), radius of horizontal curves, degree of horizontal curves, general terrain (vertical alignment), pavement seal conditions, speed limit, level of service, pavement roughness and rutting conditions. Furthermore, presence of shoulder, shoulder type (paved\unpaved), shoulder width, presence of shoulder marking, presence of median, median type (paved\unpaved), median width, and presence of median marking were extracted manually for these 521 segments and added to the dataset. Dummy values were assigned to the categorical variables for the functional classification of the road, presence of shoulder, shoulder type, presence of shoulder marking, presence of median, median type, presence of median marking, general terrain, speed limit, level of service and pavement seal conditions of road segments. The roadway geometric and traffic characteristics data were merged with crash data based on spatial coordinates of crashes and road segments. Table 2 presents descriptive statistics of explanatory variables used in this study.

Results and discussion
The crash count model (Eq. 3) and the crash severity model (Eq. 7) were first estimated to be used as the baseline for comparison purposes. The joint model of crash count and crash severity (Eq. 11) was then estimated and compared with the independent crash count and crash severity models. In all models, explanatory variables were tested for multicollinearity by computing the Pearson correlation coefficients, and the variables with unacceptably high (> 0.7) correlation coefficients were excluded from the models. In addition, AADT and segment length were used as the measures of exposure for both models. The relationship between these two variables and crash counts has been extensively studied in the road safety literature. While some research findings indicate that AADT and segment length have linear relationship with crash counts (Geedipally et al., 2009;Qin et al., 2005) and thus may be and that the non-linear relationship is warranted by estimating parameters for these two variables. In the dataset used for this study, the scatterplots of 'Crash counts versus AADT' and 'Crash counts versus segment length', as presented in Fig. 1, confirm the possible non-linear relationships between these two variables and crash counts. In addition, Anastasopoulos and Mannering (2009) have shown that estimating a parameter for segment length may reflect the boundary effect of road segmentation-crash counts may be clustered at the boundary of road segments because of a sudden change of roadway geometry. In accordance with the latter rationale, the logarithm of AADT and segment length are used as explanatory variables with estimable parameters (coefficient) in the models in this study to account for possible non-linear relationship between these two measures of exposure and crash count, and to account for potential boundary effects of road segmentation.

Baseline: independent models of crash count and crash severity
The grouped random parameters negative binomial crash count model and the fractional split ordered logit crash severity model were estimated separately using maximum simulated likelihood approach with 800 Halton draws. The required number of Halton draws was selected so that further increasing the number of draws does not change the estimates significantly. The estimated parameters of these models are presented in Table 3 and Table 4.
The results of the crash count model (Table 3) show that among all of the contributing factors to crash occurrence, ten factors are significant with at least 95 % confidence. The positive parameters of logarithm of AADT (0.286) and segment length (0.556) indicate that the likelihood of crash occurrence increases with increased exposure to crashes. However, the effect of AADT varies significantly across road segments with mean 0.286 and standard deviation 0.054. The positive parameters of number of lanes (0.283), presence of shoulder (0.063), level of service (0.087) and medium speed limit (0.598) indicate that these factors have increasing effects on total crash counts. The increasing effect of number of lanes might be due to aggregation of crashes by crash type in our study. Previous research has shown that the number of lanes have increasing effect on lane-changing related crashes along road segments (Venkataraman et al., 2014). In addition, the  increasing effect of presence of shoulder might be related to unsealed surface of the shoulder and other shoulder characteristics (e.g. edge drop). Past research has shown that unsealed shoulder is associated with increased crash counts (Cairney and McGann, 2000;Meuleners et al., 2011). On the contrary, percentage of heavy vehicle traffic (-0.222), rural functional classification (-0.243), median width (-0.095) and sealed pavement (-0.200) have decreasing effects on total crash counts. The varying parameter for percentage of heavy vehicle traffic (mean= -0.222 and standard deviation = 0.395) shows that this factor has heterogeneous effect across road segments. Finally, the inverse dispersion parameter indicates that total number of crashes are over dispersed and are correctly modelled by negative binomial specification.
The results of the crash severity model (Table 4) show that among all factors, only functional classification as the rural road and posted speed limit more than 100 km/h are significantly associated with crash severities. The positive parameters of these two variables are intuitive, indicating that these variables have increasing effect on the severity of crashes along road segments.

Joint model of crash count and crash severity
The joint model of crash count and crash severity was also estimated using 800 Halton draws in the maximum simulated likelihood estimation. Again, the number of Halton draws was selected to guarantee stability of the estimates. The estimated parameters of this joint model are presented in Table 5.
Results show that the joint model of crash count and crash severity consists of a different combination of explanatory variables compared to the crash count and crash severity models. Among the statistically significant explanatory variables, some are unique to each model component and some are common between the two model components (capturing their correlation).
Segment length, number of lanes, rural roads, presence of shoulder, median width and level of service are the statistically significant variables within the crash count model component with their parameter estimates slightly moderated than their counterparts in the independent count model. In addition, the terrain of the road is statistically significant with a negative parameter (-0.063) indicating that rolling and mountainous terrain is associated with less number of crashes and implying that drivers are more cautious in such circumstances. Sealed pavement, radius of horizontal curves and rural roads are the statistically significant variables within the crash severity model component. The negative parameter of sealed pavement indicates that sealed pavement is associated with decreased severity of crashes whereas the positive parameter of rural roads indicates that these roads are associated with increased severity of crashes. The positive parameter of radius of horizontal curves indicates that larger radius of curves (i.e. sharper curves) is associated with more injuries. This finding might suggest that drivers are more cautious along sharper curves. Numerous studies (e.g., Schneider et al., 2009;Anastasopoulos et al., 2012;Fitch and Hanowski, 2015;Oviedo-Trespalacios et al., 2018, 2019 have reported that that the complexity of the road geometry triggers riskcompensating behaviour among drivers and reduces crash risk. AADT, percentage of heavy vehicle traffic, medium speed limit and pavement roughness are the statistically significant explanatory variables that are common between the two model components and have plausible parameter estimates. The positive parameters of the logarithm of AADT (0.237), medium speed limit (0.393) and roughness (0.149) imply that these variables have increasing effects on crash count and crash severity. The negative parameter of the percentage of heavy vehicle traffic indicates that higher heavy vehicle traffic is associated with less number of crashes and lower severity of crashes.
An interesting finding from the joint model is that the dispersion parameter ( φ 1 ) of the count model component within the joint model is extremely small-albeit it is statistically significant-indicating that the over dispersion is very small. This finding implies that the unobserved heterogeneity resulted from the extra-variation in total number of crashes across sites mostly arises from ignoring the correlation between crash counts of different injury severity levels.

Comparison of goodness-of-Fit
The goodness-of-fit measures associated with the count model and  A.P. Afghari, et al. Accident Analysis and Prevention 144 (2020) 105615 the count component of the joint model shows that the latter model has lower MSPE and MAD (8.287 and 1.825, reduced from 43.230 and 1.990 respectively) and thus has improved statistical fit. In addition to MSPE and MAD, the cumulative residual plots are also plotted against increasing order to AADT to shed more light on the statistical fit of the two models. Cumulative residual plots are helpful tools in demonstrating a model fit with respect to its covariates and identifying potential and systematic bias e.g. over/under prediction (Hauer, 2015). A superior fit occurs when the plots oscillate close to zero. Excess oscillations above/under the zero axis, on the other hand, are a sign of under/over prediction. In addition, a less biased model has an approximately equal amount of positive and negative residuals. Fig. 2 presents the cumulative residual plots (adjusted to terminate at zero) for the two models and shows that the crash count component of the joint model oscillates substantially closer to zero, maintaining more balance between the positive and negative sides and staying closer within the 95 % boundaries of cumulative residuals.
The substantially improved goodness-of-fit and cumulative residuals plot of the count component of the joint model suggests that for this sample of data, at least, the joint model is preferred for predicting crash counts for the identification of high-risk sites.

Identification of high-risk sites
To identify high-risk sites across the network, the estimated models were used by three blackspot identification methods: (i) potential for improvement approach with a crash count model, (ii) excess weighted risk score approach with a two-stage mixed model of crash count and crash severity, and (iii) excess weighted risk score approach with a joint model of crash count and severity (Fig. 3).
More specifically, the crash count model was first used to determine the total predicted crash counts at individual road segments. These *Numbers inside brackets are the estimates for the standard deviations of random parameters. **Log-Likelihood of the count component within the joint model.  predicted crash counts were then used with the PFI approach to identify high-risk sites. The crash severity model was then used to determine the predicted probabilities of each severity level. The product of these probabilities and the previously mentioned predicted crash counts were then calculated to determine the predicted crash counts by severity (in a two-stage mixed model) and were used with the excess weighted risk score approach to identify high-risk sites. Finally, the joint model of crash count and crash severity were used to determine the predicted crash counts by severity while taking into account the correlation between different severity levels, and were used with the excess weighted risk score approach to identify high-risk sites. The weighted risk score for predicted crash counts of different severity levels (either from the two-stage mixed model or from the joint model) were calculated using the human capital costs (i.e. social costs) associated with Queensland crashes collected from Austroads Guide to Road Safety Part 8 (Austroads, 2015). These costs vary depending on the functional classification of the road and thus their weighted averages were used for the road network in this study. Table 6 shows the original crash costs and their weighted average across the network. Road segments were ranked by total decreasing potential for improvement (based on the crash count model) and total decreasing excess weighted risk score (based on the two-stage mixed model and the joint model of crash count and crash severity) for the period of 2010-2013.

Potential for improvement vs excess weighted risk score
To illuminate the importance of incorporating crash severity into blackspot identification, the PFI approach is first compared with the EWRS approach. Table 7 shows the top 20 road segments (out of 521 road segments) identified as high-risk sites by the two approaches. The results show that 7 out of 20 road segments (shaded cells in Table 7) are unique to the blackspot identification approach and are not commonly identified by the two approaches. The sites that are unique to the risk score approach have higher number of fatal crashes compared to the sites that are unique to the potential for improvement approach. As reported in Table 7, the ranking of sites by two approaches is different. As expected, the risk score approach puts more emphasis on fatal crashes in ranking sites, whereas the potential for improvement approach puts more emphasis on the total crashes. To obtain a more tangible understanding of the performance of these two approaches in identifying high-risk sites, the total number of crashes identified across the top 20 sites were also used for the comparison between the two approaches (Afghari et al., 2016).
at a site.
The results show that out of 1356 crashes identified by the potential for improvement approach, only 408 crashes are excess (30.0 %) whereas out of 1178 crashes identified by the excess weighted risk score approach, 618 are excess (52.4 %). This finding shows that the excess weighted risk score approach leads to the identification of more treatable crashes along the top ranked sites. The ranking of sites also paints a similar picture of model performance. The top ranked sites identified by the risk score approach have substantially higher excess crashes than the top ranked sites identified by the potential for improvement approach (for example, see sites of rank 1, rank 3, rank 4 and rank 5 in Table 8). Moreover, the excess weighted risk score approach provides additional information about excess crashes by severity levels which cannot be obtained by the other approach. This additional information indicate that the majority of the excess crashes (498 out of 618) are fatal and serious injury crashes.

Two-stage mixed model vs joint model
As illustrated previously, the PFI approach identifies higher total crashes in high rank sites while the EWRS approach identifies crashes with higher severities. However, this comparison may not provide a thorough understanding of the blackspot identification performance of the joint model within the EWRS approach because an approach that only looks into total crashes is not going to identify as many fatal and severe injury crashes compared to an approach that emphasizes severe crashes, regardless of model specification. As a result, the EWRS approach is selected as the blackspot identification approach and the performance of the joint model is now compared with that of a twostage (independent) mixed model of crash count and crash severity. Table 9 and Table 10 show the top 20 road segments identified by the EWRS approach using these two model specifications.
The results show that 4 out of 20 road segments (shaded cells in Table 9 and Table 10) are unique to the model specification and are not commonly identified by the two models. The sites that are unique to the joint model have higher number of fatal crashes compared to the sites that are unique to the two-stage mixed model. In addition, the ranking of sites by these two models is different. The two-stage mixed model identified 1257 crashes whereas the joint model identified 1178 crashes along the top 20 sites. The joint model identified 19 fatal crashes whereas the two-stage mixed model identified 14 fatal crashes along those sites. Similar to the PFI approach with the crash count model, the EWRS approach with the two-stage mixed model identified higher number of minor and serious injury crashes compared to the EWRS approach with the joint model. However, a closer look at 'excess crashes' revealed the true benefit of the joint model (Table 10). Out of 1257 crashes identified by the two-stage mixed model, only 586 crashes are excess (46.6 %), whereas out of 1178 crashes identified by the joint model, 618 are excess (52.4 %). This finding shows that the joint model leads to the identification of more treatable crashes along the top ranked sites. The ranking of sites also paints a similar picture of model performance.
The above findings together illustrate the true benefit of the joint model of crash count and crash severity when using within the excess weighted risk score approach in that it ultimately leads to the identification of higher number of treatable fatal and serious injury crashes. This benefit is the direct consequence of applying weights to the predicted crash counts by severity, and accounting for the correlation between the predicted crash counts of different severity which is in turn the consequence of employing the joint model of crash count and crash severity.

Conclusions
Total crash counts at transport locations have widely been used to identify high-risk sites across the network. However, the societal impact of crashes and their harm to the society is not accounted for in the existing approaches because crash severity has not been well incorporated into the traditional blackspot identification method. This study aimed to fill this gap by using the joint modelling approach to incorporate crash severity into crash count prediction and identify road segments with high risk of fatal and serious injury crashes.
Findings indicate that the incorporation of crash severity into the crash count model and capturing the correlation between crash counts of different severity levels improve the accuracy (i.e. statistical fit) of crash count predictions. In addition, accounting for such correlations influences the over dispersion of crash data and the statistical significance of explanatory variables. The joint model of crash count and crash severity provides unique information about the probability of Shaded: sites that are uniquely identified by each approach.
each severity level on road segments. From the pragmatic point of view and in terms of identifying highrisk sites, the joint model of crash count and crash severity enables the analyst to predict crash counts of different severity levels and provides more information about crashes compared to the traditional count model. Such information can be readily used to rank transport locations not only based on the total crashes but also based on the severe crashes such as fatal and serious injury crashes. However, it is important to note that determining which blackspot identification approach is more "useful" in practice may require additional considerations. For example, economic analysis (determining overall benefit and overall return on investment given a fixed budget) of the final selection of high-risk sites and the required countermeasures may change the superiority of one approach over another (Gross et al., 2016). Nonetheless, the weighted risk score approach presented in this study is not a replacement but may serve as a complement to the existing approaches because it provides additional information for identifying high-risk sites.

Table 10
Excess crashes along top 20 high-risk sites identified by the Excess Weighted Risk Score (EWRS) approach using two-stage mixed vs joint models of crash count and crash severity.

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.