---
category: literaturenote
citekey: fucomparisonthresholddeterminationmethods2021
title: Comparison of threshold determination methods for the deceleration rate to avoid a crash (DRAC)-based crash estimation
authors: "Fu, Chuanyun; Sayed, Tarek"
year: 2021
date: 2021-00-00 2021
doi: 10.1016/j.aap.2021.106051
publication: "Accident Analysis & Prevention"
url: "https://doi.org/10.1016/j.aap.2021.106051"
zotero_key: XPJ7HC4U
zotero_storage: A845NAKI
collections: transpordi_planeerimine
folder: 001_artiklid/Road Safety Indicators
firstAuthor: "Fu, Chuanyun"
status: converted
---

Contents lists available at [ScienceDirect](www.sciencedirect.com/science/journal/00014575)
# Accident Analysis and Prevention
journal homepage: [www.elsevier.com/locate/aap](https://www.elsevier.com/locate/aap)

# Comparison of threshold determination methods for the deceleration rate to avoid a crash (DRAC)-based crash estimation

Chuanyun Fu a,b, \*, Tarek Sayed b
- a *School of Transportation and Logistics, Southwest Jiaotong University, China* b *Department of Civil Engineering, The University of British Columbia, Canada*
- ARTICLE INFO
*Keywords:* Traffic conflicts Crash estimation DRAC MADR Bayesian hierarchical framework Extreme value theory
## ABSTRACT
There is an increased interest in the use of traffic conflicts as a surrogate safety measure and several traffic conflict indicators have been developed. One of these indicators is the deceleration rate to avoid a crash (DRAC). Generally, the greater the DRAC value, the higher the crash risk and a crash would occur when the DRAC exceeds the maximum available deceleration rate (MADR). It is noted that the MADR varies considerably for individual vehicles and depends on many factors such as the pavement conditions, vehicle weight, tire, and the braking system. Previous studies usually either set a specific value for the MADR or randomly sample values from a truncated normal distribution of MADR. However, little is known about which threshold determination approach is better. Therefore, this study aims to compare the threshold determination approaches for DRAC-based crash estimation applying Bayesian hierarchical extreme value modeling. Using traffic conflict and crash data collected from four signalized intersections in the city of Surrey, several Bayesian hierarchical models are developed for five specific values of MADR and values from two truncated normal distributions of MADR. The crash frequencies estimated from these models were compared with observed crashes. The results show that, in terms of DRACbased crash estimation accuracy, the truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68) of MADR outperforms other determination methods of MADR. Moreover, in terms of DRAC-based crash estimation accuracy and precision, the use of multisite Bayesian hierarchical models outperforms the at-site models. The truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68) of MADR is therefore recommended for DRAC-based crash estimation.
#### **1. Introduction**
The problems associated with the use of crash data in road safety analysis are well-documented [\(Sayed and Zein, 1999\)](#page-8-0). There are quality and availability problems associated with crash data. As well, the use of crash data in safety analysis is a reactive approach. Crashes need to be observed over a prolonged time period to identify safety problems and to correct them. Therefore, the development of surrogate safety measures such as traffic conflicts for use in safety analysis has attracted considerable research interest. Compared to traffic crashes, traffic conflicts occur more frequently, can be easily observed, and provide better insight into the collision failure mechanism ([Sayed and Zein, 1999](#page-8-0)).
Over the past decades, a number of traffic conflict indicators have been established and used in traffic safety analysis. Examples of most widely used conflict indicators include temporal and spatial proximity indicators and deceleration-based indicators. Several researchers have explicitly recognized the deceleration rate to avoid a crash (DRAC) as a useful surrogate safety measure [\(Archer, 2005](#page-8-0)). The DRAC is the rate at which a following vehicle must decelerate to avoid a collision with the leading vehicle ([Zheng and Sayed, 2019a](#page-9-0)). It is defined as the differential speed between a following vehicle and its leading vehicle divided by their closing time ([Cooper and Ferguson, 1976;](#page-8-0) [Gettman and Head,](#page-8-0) [2003\)](#page-8-0). Compared to the other surrogate safety measures, such as time to collision (TTC) and the post encroachment time (PET), the DRAC has some advantages, since it could reflect the evasive action (e.g., required deceleration) of the following vehicle to come to a timely stop or reach the corresponding speed of the leading vehicle and thereby avoid a collision ([Cunto and Saccomanno, 2008](#page-8-0)).
The greater the DRAC value the higher the crash risk, and a crash would occur if the DRAC exceeds the maximum available deceleration rate (MADR). Therefore, the crash risk based on DRAC can be calculated given a MADR. Since the MADR is unique for individual vehicles and
\* Corresponding author at: School of Transportation and Logistics, Southwest Jiaotong University, China. *E-mail addresses:* [fcy1128@mail.ubc.ca](mailto:fcy1128@mail.ubc.ca) (C. Fu), [tsayed@civil.ubc.ca](mailto:tsayed@civil.ubc.ca) (T. Sayed).
depends on pavement conditions, vehicle weight, tire and braking system, and other factors, it is usually assumed to follow a truncated normal distribution N(8.45, 1.42 ), with an upper limit of 12.68 m/s2 and a lower limit of 4.23 m/s2 ([Cunto and Saccomanno, 2008](#page-8-0)). Several studies have investigated the crash risk based on DRAC by mainly using the extreme value distribution (EVT) and the truncated normal distribution. The former treats the conflict extremes characterized by DRAC as a random variable, while the latter treats the MADR as a random variable. Additionally, a threshold needs to be developed to identify traffic conflict extremes for use in the establishment of the EVT model based on DRAC.
Previous studies used different specific values for MADR to conduct DRAC-based traffic conflict safety evaluation. A value of 8 m/s2 for MADR was considered as a conservative value that all vehicles could reach in one previous study [\(Ceunynck, 2017\)](#page-8-0). [Wang et al. \(2019\)](#page-8-0) employed the MADR values of 8 m/s2 and 12 m/s2 for crash estimation from the bivariate extreme value models, respectively. In the study by ([Zheng et al., 2019b](#page-9-0)), the mean value (i.e., 8.45 m/s2 ) of the truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68) of MADR was selected as the threshold for crash occurrence, which was applied to develop the bivariate extreme value models for crash estimation. It is noted that a specific value of MADR may lead to DRAC-based crash underestimation or overestimation at some sites. Therefore, [Zheng and Sayed \(2019a\)](#page-9-0) used the values of DRAC sampled from the fitted extreme value distribution and the values of MADR sampled from a given truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68) to calculate the DRAC-based crash risk. As well, [Fu et al. \(2020\)](#page-8-0) determined the values of MADR for DRAC-based crash estimation by sampling from the truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68).
In addition to the extreme value distribution, the truncated normal distribution is also used to estimate the crash risk based on DRAC. In other words, the MADR is regarded as a random variable and is assumed to follow a truncated normal distribution. The crash risk can be calculated as the probability of the MADR being less than the DRAC. [Meng](#page-8-0) [and Weng \(2011\)](#page-8-0) evaluated the rear-end crash risk measured by DRAC at work zones. They assumed that the MADR of cars and trucks followed truncated normal distributions N(8.45, 1.42 )I(1.23, 12.68) and N(6.82, 1.42 )I(0.60, 10.05), respectively. This approach was also applied to investigate the work zone rear-end crash risk based on DRAC for different vehicle-following patterns ([Weng et al., 2014](#page-8-0)), assess the performance of two microscopic safety indices (i.e., DRAC and TTC) on safety evaluation for expressways [\(Qu et al., 2014](#page-8-0)), and explore the rear-end crash risk based on DRAC and TTC in work zone merging areas ([Weng et al., 2015](#page-8-0)). With the updated data of the braking performance of various automobiles, a new truncated normal distribution of MADR N (9.7, 1.32 )I(4.2, 12.7) was obtained ([Wang and Stamatiadis, 2013](#page-8-0), [2014\)](#page-8-0). The two studies employed the distribution of MADR to calculate the probability of MADR being less than the required braking rate to avoid a crash conditional on the reaction time being less than the TTC.
The review of previous studies using DRAC in road safety analysis shows that limited research has been undertaken on selecting which threshold determination approach for DRAC-based traffic conflict analysis (i.e., setting specific values or distributions of MADR) would yield the best results in terms of accuracy and precision.
Therefore, the objective of this study is to compare the threshold determination approaches for DRAC-based crash estimation in terms of accuracy and precision. For the comparison, a Bayesian hierarchical extreme value model that combines traffic conflict extremes from different sites, incorporates the impact of various covariates, and considers the unobserved site-specific heterogeneity is developed.
The structure of this paper is organized as follows: Section 2 presents the data used in this study. Section 3 provides details of the proposed Bayesian hierarchical extreme value model, including the generalized extreme value model, the Bayesian hierarchical structure, the model estimation method, the model choice, and the crash estimation approach. The model estimation results, model comparisons, and variables interpretation are discussed in Section [4](#page-4-0). The last section, Section [5](#page-8-0), presents the conclusions of this study and the recommended future work.
## **2. Data description**
The data used for this study include video data and crash data. The video data were collected from four signalized intersections in the city of Surrey, British Columbia, Canada. They are the intersections of 72 Ave & 128 St, 72 Ave & 132 St, 64 Ave & King George Blvd, and Fraser Hwy & 168 St (see [Fig. 1](#page-2-0)).
Video data were collected by cameras that focus on the approaches where most of rear-end conflicts occur. Computer vision techniques were applied to identify signal cycles, extract vehicle trajectories, and generate space-time diagram of vehicles from the video footages. The detailed video processing procedures are presented in ([Essa and Sayed,](#page-8-0) [2018,](#page-8-0) [2019](#page-8-0)). The DRAC was then automatically calculated from the extracted vehicle trajectories. It is noted that the conflict indicator is computed continuously (i.e., frame by frame and one second = 30 frames) over time for an event, and the maximum value is selected for DRAC. For each pair of consecutive vehicles on the through lane, the DRAC could be acquired. Additionally, three traffic parameters at the signal cycle level were extracted from the space-time diagram, including: traffic volume (V, vehicles/lane/cycle) which is the number of vehicles per cycle on the studied lane; shock wave area (A, km s) which is a triangle area where it is the state II of classical four traffic flow states theory ([May, 1990\)](#page-8-0), as illustrated in [Fig. 2](#page-2-0); and platoon ratio (P) which is defined as the proportion of all vehicles arriving during green multiplied by the ratio of the signal cycle length to the effective green time of the subject movement [\(Highway Capacity Manual, 2000\)](#page-8-0).
In this study, only DRAC was utilized to characterize traffic conflicts. It is calculated as ([Cooper and Ferguson, 1976](#page-8-0); [Gettman and Head,](#page-8-0) [2003\)](#page-8-0):
$$DRAC_{F,I} = \frac{\left(v_{F,I} - v_{L,I}\right)^2}{2\left(x_{L,I} - x_{F,I} - D_L\right)} \tag{1}$$
where, *vF,t*and *vL,t*are speeds of following vehicle and leading vehicle at time *t*, respectively; *xF,t*and *xL,t*are positions of following vehicle and leading vehicle at time *t*, respectively; *DL* is the length of leading vehicle.
It is noted that there may be zero or even no values of DRAC for a pair of consecutive vehicles. Therefore, only events with DRAC *>* 0 m/s2 were considered as candidate traffic conflicts for the following analysis. A total of 2202 events with DRAC values were obtained from 416 signal cycles, which is on average 5.3 events per cycle. The summary statistics of traffic conflict and crash data are shown in [Table 1.](#page-3-0)
Three years of crash data (2013–2015) at the four signalized intersections in the city of Surrey, British Columbia was used. The data, provided by the Insurance Corporation of British Columbia, included the information such as location, date, time, crash type, and crash severity, among others. Since the crash data is used to assess the performance of the proposed model, only rear-end crashes that occurred in the through lanes at the target approaches were selected. Moreover, since rear-end conflicts were observed in the daytime period, only rear-end crashes happened during the daytime were kept. Eventually, a total of 43 rearend crashes remained, as presented in [Table 1.](#page-3-0)
#### **3. Model development**
# *3.1. Generalized extreme value model*
In this study, the generalized extreme value (GEV) distribution is fitted based on the block maxima (BM) extremes. The observations are aggregated into fixed blocks, and the maxima of each block is treated as extreme ([Coles, 2001;](#page-8-0) [Zheng et al., 2018](#page-9-0)). Suppose that *X*1*, X*2*,* ⋯*, Xn* is a sequence of independent random variables having a common distribution function and *Mn* = max{*X*1*, X*2*,* ⋯*, Xn*} represents the block maxima.

**Fig. 1.** Study intersections and corresponding approaches.

**Fig. 2.** Classical four different traffic flow states at signalized intersections.
If there exist sequences of constants {*an >* 0} and {*bn*} such that
$$\Pr\{(M_n - b_n)/a_n \le z\} = G(z) \text{ as } n \to \infty$$
(2)
where, *G* is a non-degenerate distribution function. It is a member of GEV family [\(Coles, 2001\)](#page-8-0), which is given by:
$$G(z) = exp\left\{ -\left[1 + \xi \frac{(z - \mu)}{\sigma}\right]^{-1/\xi} \right\}$$
(3)
defined on the set {*z* : 1 + *ξ*(*z* − *μ*)*/σ >* 0}, where *μ* is the location parameter, − ∞ *< μ <* ∞; *σ* is the scale parameter, *σ >* 0; *ξ* is the shape parameter, − ∞ *< ξ <* ∞.
# *3.2. Bayesian hierarchical structure*
A Bayesian hierarchical structure is developed for the univariate extreme value model, which has three layers: data layer, process layer, and prior layer [\(Zheng and Sayed, 2019b; Zheng et al., 2019a;](#page-9-0) [Fu et al.,](#page-8-0) [2020;](#page-8-0) [Zheng and Sayed, 2020](#page-9-0)).
## *3.2.1. Data layer*
A GEV distribution given by Eq. (2) is employed to model the data at the first layer of hierarchy, which is
$$G(z_{ij} < z; \mu_i, \sigma_i, \xi_i) = exp \left\{ -\left[1 + \xi_i \frac{(z - \mu_i)}{\sigma_i}\right]^{-1/\xi_i} \right\}$$
(4)
where, *zij* is the *j* th recorded block maxima at site *i*, with *i* = 1*,* 2*,* ⋯*, m* and *j* = 1*,*2*,*⋯*,ni*; (*μi , σi, ξi*) represent the parameters at site *i*. Given that *zij* corresponds to the block maxima, it can be depicted by the GEV distribution whose parameters are dependent on the site. To ensure a positive scale parameter, the *σ* is reparameterized as *ψ* = log*σ*. Therefore, the Eq. (4)
can be rewritten as:
$$G(z_{ij} < z; \mu_i, \sigma_i, \xi_i) = exp \left\{ -\left[1 + \xi_i \frac{(z - \mu_i)}{\exp(\psi_i)}\right]^{-1/\xi_i} \right\}$$
(5)
**Table 1** Summary statistics of data.
| Data information | | 72 Ave &
128 St | 72 Ave &
132 St | 64 Ave &
King
George
Blvd | Fraser
Hwy &
168 St |
|------------------------------------------|---------------------------------|----------------------------------------------|-----------------------------------|------------------------------------|--------------------------------------|
| Video was recorded in | | March 28,
2012
(10− 12 am;
2− 4 pm) | April 3,
2012
(9− 12
am) | June 10,
2015 (1− 5
pm) | June 10,
2015 (9
am − 1
pm) |
| Number of cycles | | 130 | 72 | 103 | 111 |
| | | 928 | 303 | 436 | 535 |
| Event counts
DRAC value
(m/s2
) | Mean
Min
Max | 2.05
0.00
7.10 | 1.81
0.16
6.38 | 2.41
0.01
5.98 | 2.60
0.06
6.70 |
| V (vehicle | Mean | 13.52 | 13.43 | 11.62 | 13.71 |
| count per | Min | 2 | 7 | 4 | 3 |
| cycle) | Max | 28 | 25 | 18 | 23 |
| A (km⋅s) | Mean | 1.38 | 0.96 | 1.19 | 1.10 |
| | Min | 0.00 | 0.00 | 0.00 | 0.00 |
| | Max | 4.45 | 2.90 | 3.03 | 5.66 |
| P | Mean | 1.10 | 1.07 | 1.40 | 1.12 |
| | Min | 0 | 0.14 | 0.00 | 0.32 |
| | Max | 1.88 | 1.72 | 2.11 | 1.75 |
| Observed
crashes | 2013
2014
2015
Average | 8
3
3
4.7 | 0
1
0
0.3 | 9
7
2
6.0 | 4
3
3
3.3 |
#### *3.2.2. Process layer*
A latent Gaussian process underlying extremes is characterized through developing a structure that connects the parameters in the data layer with certain covariates. The identity link function is used to build the connection. It is noted that since it is hard to precisely estimate the shape parameter *ξ* of GEV distribution, it is not possible to model *ξ* as a smooth function of certain covariates. Hence, the shape parameter *ξ* is hypothesized to be unknown but not influenced by covariates during the process [\(Coles, 2001](#page-8-0); [Cooley et al., 2006\)](#page-8-0). As such, the GEV distribution parameter can be expressed as:
$$\begin{cases} \mu_i = \beta_{\mu 0} + \beta_{\mu} \mathbf{X} \\ \psi_i = \beta_{\mu 0} + \beta_{\psi} \mathbf{X} \\ \xi_i = \beta_{20} \end{cases}$$
$$\tag{6}$$
where, *βμ*0, *βψ*0, and *βξ*0 are intercepts; **X** is the vector of covariates; **β***μ* and **β***ψ*are the vectors of estimated parameters of covariates. In fact, there usually seems to be additional heterogeneity that is not directly addressed by covariates. In order to capture the additional unobserved heterogeneity, random error terms are also introduced as:
$$\begin{cases} \mu_{i} = \beta_{\mu 0} + \beta_{\mu} \mathbf{X} + \varepsilon_{\mu i} \\ \psi_{i} = \beta_{\psi 0} + \beta_{\psi} \mathbf{X} + \varepsilon_{\psi i} \\ \xi_{i} = \beta_{20} + \varepsilon_{\xi i} \end{cases}$$
(7)
where, *εμi*, *εψi*, and *εξi* are random effects reflecting between-site variation, which are consistent for extremes in the same site but vary between different sites. Therefore, the model for the process layer can be expressed in the following form of a random intercept:
$$\begin{cases} \mu_{i} = \beta_{\mu i} + \beta_{\mu} \mathbf{X} \\ \psi_{i} = \beta_{\psi i} + \beta_{\psi} \mathbf{X} \\ \xi_{i} = \beta_{\xi i} \end{cases}$$
(8)
where, *βμi*, *βψi*, and *βξi*are site-specific random intercepts.
#### *3.2.3. Prior layer*
In this layer, priors are assigned to site-specific random intercepts and the estimated parameters of covariates for characterizing the latent process. It is hypothesized that each parameter in this layer is independent to the others. Due to no available prior information regarding how the parameters of GEV distribution are associated with covariates, uninformative priors for these parameters, i.e., {*βμi*, *βψi*, **β***μ*, **β***ψ* }, are adopted. They are hypothesized to be normally distributed with zero mean and large variance, i.e., *N*(0*,*106). It is worth noting that the GEV model is sensitive to its shape parameter, so an improper defined prior would lead to a convergence issue during the model estimation. According to the previous road safety studies such as [\(Songchitruksa and](#page-8-0) [Tarko, 2006\)](#page-8-0) and ([Zheng et al., 2014](#page-9-0)), the estimated shape parameters are usually within the range of (-1.0, 1.0). Therefore, the priors for shape parameters are assumed to follow a uniform distribution, i.e., Unif(-1.0, 1.0).
## *3.3. Model estimation method*
In the above-mentioned three-layer hierarchical model, according to the Bayes theorem, the posterior inference for the univariate Bayesian hierarchical extreme value model parameters **θ** given the data **Ω** is
$$p(\mathbf{\theta}|\mathbf{\Omega}) \propto \underbrace{p_1(\mathbf{\Omega}|\mathbf{\theta}_1) p_2(\mathbf{\theta}_1|\mathbf{\theta}_2) p_3(\mathbf{\theta}_2)}_{Data} \underbrace{p_{rocess} p_{rior}}_{Process}$$
(9)
where, *p*(**θ**|**Ω**) is the posterior distribution of the parameter **θ**; *ps* is the density function associated with layer *s* (*s* = 1*,* 2*,* 3) of the hierarchical model and depends on parameters **θ***v*(*v* = 1*,*2).
The density associated with the data layer is written as:
$$p_{1}(\mathbf{Z}|\boldsymbol{\theta}_{1}) = \prod_{i=1}^{m} \prod_{j=1}^{n_{1}} \left(\frac{1}{\exp\left(\boldsymbol{y}_{i}\right)}\right) \exp \left\{ -\left[1 + \xi_{i} \frac{\left(z_{ij} - \mu_{i}\right)}{\exp\left(\boldsymbol{y}_{i}\right)}\right]^{-1/\xi_{i}} \right\} \left[1 + \xi_{i} \frac{\left(z_{ij} - \mu_{i}\right)}{\exp\left(\boldsymbol{y}_{i}\right)}\right]^{-1-1/\xi_{i}}$$
$$(10)$$
where, **θ1** = [**μ***,* **ψ***,* **ξ**] *T*; *m* is the number of study sites (i.e., intersections); *ni* is the number of blocks at site *i*. According to the GEV parameters specified in Eq. (8), the density related to the process layer is obtained as:
$$\begin{split} p_{2}(\boldsymbol{\theta}_{1}|\boldsymbol{\theta}_{2}) &= \frac{1}{\delta_{\mu}\sqrt{2\pi}} \exp\left[-\frac{(\mu-\mu_{i})^{2}}{2\delta_{\mu}^{2}}\right] \times \frac{1}{\delta_{\psi}\sqrt{2\pi}} \exp\left[-\frac{(\psi-\psi_{i})^{2}}{2\delta_{\psi}^{2}}\right] \\ &\times \frac{1}{\delta_{\xi}\sqrt{2\pi}} \exp\left[-\frac{(\xi-\xi_{i})^{2}}{2\delta_{\xi}^{2}}\right] \end{split} \tag{11}$$
where, **θ2** = [**βμ***,* **βψ***,* **βξ**] *T*. With the priors specified, the density associated with the prior layer is given by
$$p_3(\mathbf{\theta}_2) = p_{\beta_{\mu}}(\beta_{\mu}) \times p_{\beta_{\psi}}(\beta_{\psi}) \times p_{\beta_{\xi}}(\beta_{\xi})$$
(12)
The posterior distributions of model parameters can be inferred by virtue of Markov Chain Monte Carlo (MCMC) simulation using Gibbs Sampling.
#### *3.4. Model choice*
When developing the univariate Bayesian hierarchical extreme value model with relating the parameters of GEV distribution to several covariates, there will be many models to choose from. In the Bayesian framework, the Deviance Information Criterion (DIC) is commonly adopted to compare models. The basic principle of DIC is parsimony, obtaining the simplest model that explains as much of the variation in the data as possible ([Spiegelhalter et al., 2002\)](#page-8-0). The DIC is computed as follows:
$$DIC = \overline{D} + n_D \tag{13}$$
where, *D* is the posterior mean deviance that measures the model fitting and *nD* is the effective number of parameters in the model. Generally, a smaller DIC value indicates a better model. It is also suggested that a difference greater than 10 can definitely rule out the model with the higher DIC value, a difference ranging from 5 to 10 is considered substantial, and a difference less than 5 makes the models competitive ([El-Basyouny and Sayed, 2012b; 2012a](#page-8-0)).
## *3.5. Crash estimation method*
When the GEV (*μ, σ, ξ*) is fitted based on the data from individual sites, it is called the at-site model. As well, it is called the multisite model when the GEV (*μ, σ, ξ*) is fitted based on the data from several sites simultaneously. In this study, a traffic conflict is characterized by DRAC. A crash would occur if the DRAC exceeds the MADR [\(Cunto and Sac](#page-8-0)[comanno, 2008](#page-8-0)). As such, with the fitted GEV distribution of maximum DRAC, the risk of crash can be calculated as:
$$R_{i} = \Pr\{Z_{i} \ge MADR\} = 1 - G_{i}(MADR)$$
$$= 1 - \exp\left\{-\left[1 + \xi_{i}^{MADR} - \mu_{i}\right]^{-1/\xi_{i}}\right\}$$
(14)
where, *Ri* is the risk of crash for block *i*; *Xi* is the maximum DRAC in block *i*; *Gi*(⋅) is the fitted univariate Bayesian hierarchical extreme value distribution for block *i*. Five specific values (i.e., 7 m/s2 , 8 m/s2 , 8.5 m/ s 2 , 9.7 m/s2 , and 12 m/s2 ) of MADR are employed to calculate the risk of crash, respectively. As well, two truncated normal distributions (i.e., N (8.45, 1.42 )I(4.23, 12.68) and N(9.7, 1.32 )I(4.23, 12.68)) of MADR are used to generate the sampled values of MADR for crash risk calculation.
#### **4. Results and discussions**
#### *4.1. Model estimation results*
With the maximum values of DRAC from different signal cycles of four signalized intersections, the multisite Bayesian hierarchical extreme value models were developed. The estimations of the proposed models were conducted using WinBUGS ([Lunn et al., 2000](#page-8-0)). Two separate chains for each model parameter with diverse initial values were set to run 50,000 iterations and the first 20,000 iterations were discarded as burn-in samples. The posterior estimates were acquired from the remaining 30,000 iterations. The model convergence was checked by calculating the Gelman-Rubin statistic of two chains for each parameter and visually examining the trace plots of parameter chains ([Gelman and](#page-8-0)
[Rubin, 1992](#page-8-0); [El-Basyouny and Sayed, 2009; Zheng and Ismail, 2017\)](#page-8-0).
Four multisite Bayesian hierarchical models were established: BHM\_GEV(0,0,0), BHM\_GEV(1,0,0), BHM\_GEV(0,1,0), and BHM\_GEV (1,1,0). All the four models have site-specific random error terms for the location parameter, scale parameter, and shape parameter. BHM\_GEV (0,0,0) is the multisite stationary Bayesian hierarchical model without covariates. BHM\_GEV(1,0,0) and BHM\_GEV(0,1,0) are the multisite nonstationary Bayesian hierarchical models with the three covariates added to the location parameter and scale parameter separately. BHM\_GEV (1,1,0) is the multisite non-stationary Bayesian hierarchical model with covariates added to both location parameter and scale parameter simultaneously, and only significant covariates are included in this model. The estimation results of four multisite Bayesian hierarchical models are presented in Table 2.
As shown in Table 2, the DIC values are 1562.089, 1394.052, 1543.146, and 1384.523 for the multisite BHM\_GEV(0,0,0), multisite BHM\_GEV(1,0,0), multisite BHM\_GEV(0,1,0), and multisite BHM\_GEV (1,1,0), respectively. In general, the stationary model has a greater DIC value than those of non-stationary models. The smallest DIC difference between stationary model and non-stationary models is 18.943, which shows that the non-stationary models considerably outperform the stationary model. This also confirms that the introduced covariates in the non-stationary models explain significant variations in the observed data. Moreover, the multisite BHM\_GEV(1,1,0) has the smallest DIC value, which indicates that it is the best fitted model and better than other models in terms of goodness of fit.
## *4.2. Covariates interpretation*
As discussed in [\(Zheng et al., 2019a](#page-9-0)), it is not straightforward to interpret the influences of covariates on safety, since they directly impact the location parameter and scale parameter of the GEV distribution. The location and scale parameters determine the shape convexity of GEV distribution, which includes the part of the distribution where most crash events are located. Given a specific MADR, an illustration of GEV distributions of DRAC with different location parameters and scale parameters is presented in [Fig. 3](#page-5-0). Given the fixed scale parameter and shape parameter, the GEV distribution with a smaller location parameter is farther away from the crash risk area, implying higher safety level. Similarly, given the fixed location parameter and
**Table 2** Estimation results of multisite Bayesian hierarchical models.
| | | BHM_GEV(0,0,0) | | BHM_GEV(1,0,0) | | BHM_GEV(0,1,0) | | BHM_GEV(1,1,0) | |
|-----------------|--------|----------------|--------|----------------|--------|----------------|--------|----------------|--------|
| Model parameter | | Mean | S.D. | Mean | S.D. | Mean | S.D. | Mean | S.D. |
| | βμ0[1] | 2.9090** | 0.1612 | 1.5410** | 0.3118 | 2.5930** | 0.1714 | 1.6250** | 0.2769 |
| μ | βμ0[2] | 2.0830** | 0.1330 | 0.9711** | 0.3025 | 2.0100** | 0.1287 | 1.0950** | 0.2840 |
| | βμ0[3] | 2.7340** | 0.1910 | 1.7710** | 0.3607 | 2.4670** | 0.1963 | 1.6880** | 0.3364 |
| | βμ0[4] | 2.9670** | 0.1978 | 1.6910** | 0.3202 | 2.7480** | 0.1962 | 1.7450** | 0.2965 |
| | βμ(V) | – | – | 0.0999** | 0.0182 | – | – | 0.0825** | 0.0202 |
| | βμ(A) | – | – | 0.6016** | 0.0878 | – | – | 0.7287** | 0.1088 |
| | βμ(P) | – | – | − 0.6242** | 0.2183 | – | – | − 0.5974** | 0.2060 |
| | βΨ0[1] | 0.4562** | 0.0794 | 0.2203** | 0.0694 | 0.1955 | 0.1457 | − 0.0259 | 0.0999 |
| ψ = | βΨ0[2] | 0.0150 | 0.0972 | 0.0362 | 0.0900 | − 0.2311 | 0.1580 | − 0.1767 | 0.1128 |
| log(σ) | βΨ0[3] | 0.5244** | 0.0929 | 0.2391** | 0.0833 | 0.3245* | 0.1634 | 0.0179** | 0.1092 |
| | βΨ0[4] | 0.5208** | 0.0966 | 0.2511** | 0.0866 | 0.2907** | 0.1509 | 0.0651** | 0.1023 |
| | βΨ(V) | – | – | – | – | 0.0167** | 0.0091 | – | – |
| | βΨ(A) | – | – | – | – | 0.1309** | 0.0497 | 0.1619** | 0.0524 |
| | βΨ(P) | – | – | – | – | − 0.0964 | 0.0926 | – | – |
| | βξ0[1] | − 0.2369** | 0.0893 | − 0.2131** | 0.0521 | − 0.3069** | 0.0688 | − 0.1450** | 0.0586 |
| | βξ0[2] | 0.0318 | 0.0851 | − 0.1363* | 0.0690 | − 0.0458 | 0.0833 | − 0.1611** | 0.0754 |
| ξ | βξ0[3] | − 0.4058** | 0.1055 | − 0.2970** | 0.0812 | − 0.4181* | 0.0838 | − 0.2333** | 0.0767 |
| | βξ0[4] | − 0.2262** | 0.1275 | − 0.1135** | 0.1005 | − 0.3026** | 0.0954 | − 0.1304 | 0.0946 |
| D | | 1550.000 | | 1380.000 | | 1530.000 | | 1370.000 | |
| nD | | 12.089 | | 14.052 | | 13.146 | | 14.523 | |
| DIC | | 1562.089 | | 1394.052 | | 1543.146 | | 1384.523 | |
Note: "–" means the covariate is not included; S.D. is the standard deviation; "\*" and "\*\*" mean the covariate is significant at 0.1 and 0.05 level, respectively; "V" is the traffic volume; "A" is the shock wave area; and "P" is the platoon ratio.

**Fig. 3.** Illustration of GEV distributions with different location and scale parameters.
shape parameter, the GEV distribution with a smaller scale parameter has a skinny shape, and its tail area within the crash risk area becomes narrower, indicating higher safety level as well.
According to the results from the best fitted model, multisite BHM\_GEV(1,1,0), all the three covariates are significantly related to the location parameter of GEV distribution, while only one covariate is significantly associated with the scale parameter of GEV distribution. Specifically, the traffic volume and shock wave area are positively related to the location parameter of GEV distribution, whereas the platoon ratio is negatively associated with the location parameter of GEV distribution. Moreover, the shock wave area is positively related to the scale parameter of GEV distribution. Both the greater traffic volume and shock wave area will increase the location parameter, which yields more expected number of crashes. However, a greater platoon ratio will decrease the location parameter, and thus produce less expected number of crashes. Similarly, a greater shock wave area will increase the scale
parameter, which generates more expected number of crashes.
#### *4.3. Comparison of threshold determination methods for DRAC*
The best fitted model, multisite BHM\_GEV(1,1,0), was used to predict the annual number of crashes for each intersection, where *T* = 365 × 12 = 4380, and t = 4, 3, 4, and 4 for intersections of 72 Ave & 128 St, 72 Ave & 132 St, 64 Ave & King George Blvd, and Fraser Hwy & 168 St, respectively. Five specific values (i.e., 12 m/s2 , 9.7 m/s2 , 8.5 m/s2 , 8 m/s2 , and 7 m/s2 ) and two truncated normal distributions (i.e., N(9.7, 1.32 )I(4.23,12.68) and N(8.45, 1.42 )I(4.23,12.68)) for MADR were used. In addition, the uncertainties of estimated crashes were determined directly in WinBUGS during the model posterior inference via different quantiles, and this is also an advantage of Bayesian inference that can incorporate all uncertainties in the data and model parameters ([El-Basyouny and Sayed, 2012b;](#page-8-0) [Zheng and Sayed, 2020\)](#page-9-0). The crash estimation results from the multisite BHM\_GEV(1,1,0) are shown in Table 3. In the table, the 95% confidence intervals (95% C.I.) of observed crashes are also presented, given that the observed crashes at each intersection are just one observation of the underlying process ([Zheng and Sayed, 2020](#page-9-0)). By hypothesizing that the observed crashes follow the Poisson distribution, the 95% C.I.s were computed as { *λ* : 1 2*nχ*2 2*y*0*,*0*.*975 ≤ *λ* ≤ 1 2*nχ*2 2(*y*0+1)*,*0*.*025} , where *n* = 3 is the number of years for crash data collection, and *y*0 is the total number of crashes in the 3-year period for an intersection [\(Songchitruksa and Tarko, 2006\)](#page-8-0).
As shown in Table 3, there are three intersections where the estimated crashes from multisite BHM\_GEV(1,1,0) using the MADR of 12 m/ s 2 are smaller than the observed crashes and outside the 95% confidence intervals of the observed crashes. Similar results are obtained from multisite BHM\_GEV(1,1,0) using the MADR of 9.7 m/s2 . There are two intersections where the estimated crashes from multisite BHM\_GEV (1,1,0) using the MADR of 8.5 m/s2 are smaller than the observed crashes and outside the 95% confidence intervals of the observed crashes, whereas the estimated crashes of other two intersections are within the 95% confidence intervals of the observed crashes. The multisite BHM\_GEV(1,1,0) using the MADR of 8 m/s2 generates crash underestimation at two intersections with the estimated crashes outside the 95% confidence intervals of the observed crashes, while produces crash overestimation at one intersection with the estimated crashes
**Table 3** Crash estimation results of multisite Bayesian hierarchical models with different thresholds.
outside the 95% confidence intervals of the observed crashes. Although the multisite BHM\_GEV(1,1,0) using the MADR of 7 m/s2 yields the estimated crashes inside the 95% confidence intervals of the observed crashes at two intersections, it produces serious crash overestimations at other two intersections. The estimated crashes from multisite BHM\_GEV (1,1,0) using the values of MADR sampled from a truncated normal distribution N(9.7, 1.32 )I(4.23, 12.68) are similar to those from multisite BHM\_GEV(1,1,0) using a specific value of 8.5 m/s2 for MADR. However, the estimated crashes from multisite BHM\_GEV(1,1,0) using the values of MADR sampled from a truncated normal distribution N(8.45, 1.42 )I (4.23, 12.68) are all within the 95% confidence intervals of the observed crashes.
The results of crash estimation indicate that in terms of accuracy, the multisite BHM\_GEV(1,1,0) using the sampled values for MADR is better than the one using a specific value for MADR. This is due to the MADR varying with the vehicle, pavement, and environmental conditions at different sites. A larger specific value of MADR would lead to crash underestimation at some sites, while a smaller specific value of MADR would generate crash overestimation at these sites. By using the multisite BHM\_GEV(1,1,0), the comparison of estimated crashes among different MADR values is presented in Fig. 4. In terms of crash estimation accuracy, among 16 different specific values of MADR, the best ones are 7.5 m/s2 , 8.4/8.5/8.6 m/s2 , 6.8 m/s2 , and 8.7 m/s2 for the four intersections, respectively. According to the 95% confidence interval of observed crashes, the ranges of appropriate MADR are 7.2–7.8 m/s2 , 7.3–8.6 m/s2 , 6.6–7.1 m/s2 , 8.2–9.5 m/s2 for the four intersections, respectively. Therefore, it is hard to find an appropriate specific value or range of MADR for DRAC-based crash estimation at all these intersections.
Additionally, the multisite BHM\_GEV(1,1,0) using the values of MADR sampled from a truncated normal distribution N(8.45, 1.42 )I (4.23, 12.68) outperforms the one with the sampled values of MADR from a truncated normal distribution N(9.7, 1.32 )I(4.23, 12.68). Crash underestimation at two intersections from the multisite BHM\_GEV (1,1,0) using the values of MADR sampled from a truncated normal distribution N(9.7, 1.32 )I(4.23, 12.68) may be a result of the inappropriate upper and lower limits of such a truncated normal distribution. Therefore, in terms of crash estimation accuracy, it is recommended to determine the value of MADR by sampling from the truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68).
# *4.4. Comparison of multisite and at-site Bayesian hierarchical models*
The at-site Bayesian hierarchical models were developed for each site by following the same procedure as the multisite Bayesian hierarchical models. All the three covariates were also considered in the model development. Those at-site Bayesian hierarchical models with the smallest DIC values were recognized as the best fitted models, and their estimation results are shown in [Table 4.](#page-7-0) It can be seen that the sign of significant covariates of the at-site models is the same as the multisite BHM\_GEV(1,1,0), but the significance of covariates are varying across the sites.
With the values of MADR sampled from the truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68), the best fitted at-site models were employed to predict the annual crash frequencies, which were then compared to the observed crashes and the estimated crashes from the

**Fig. 4.** Comparison of estimated crashes among different MADR values (m/s2 ).
**Table 4** Estimation results of the best at-site Bayesian hierarchical models.
| | 72 Ave & 128 St | | 72 Ave & 132 St | | 64 Ave & King George Blvd | | Fraser Hwy & 168 St | | |
|-----------------|-----------------|----------------|-----------------|----------------|---------------------------|----------------|---------------------|----------------|--------|
| Model parameter | | BHM_GEV(1,1,0) | | BHM_GEV(1,1,0) | | BHM_GEV(1,0,0) | | BHM_GEV(1,1,0) | |
| | | Mean | S.D. | Mean | S.D. | Mean | S.D. | Mean | S.D. |
| μ | βμ0 | 1.5760** | 0.4281 | 1.3260** | 0.1477 | 1.2050** | 0.6346 | 0.6576** | 0.4479 |
| | βμ(V) | 0.1327** | 0.0324 | – | – | 0.0989** | 0.0453 | 0.0902** | 0.0354 |
| | βμ(A) | 0.3402** | 0.1334 | 1.0380** | 0.1929 | 1.0350** | 0.1997 | 1.0670** | 0.2185 |
| | βμ(P) | − 0.7599** | 0.3658 | – | – | − 0.5886** | 0.3608 | – | – |
| log(σ) | βΨ0 | − 0.1577** | 0.2235 | − 0.5869** | 0.1421 | 0.1947** | 0.0815 | 0.0088 | 0.1445 |
| | βΨ(V) | 0.0237* | 0.0141 | – | – | – | – | – | – |
| | βΨ(A) | – | – | 0.5125** | 0.1159 | – | – | 0.2181** | 0.1110 |
| | βΨ(P) | – | – | – | – | – | – | – | – |
| ξ | βξ0 | − 0.1373** | 0.0673 | − 0.1592* | 0.0797 | − 0.2525** | 0.0748 | − 0.1391 | 0.0946 |
| D | | 430.000 | | 200.000 | | 330.000 | | 390.000 | |
| nD | | 5.868 | | 4.451 | | 5.153 | | 6.082 | |
| DIC | | 435.868 | | 204.451 | | 335.153 | | 396.082 | |
Note: "–" means the covariate is not included; S.D. is the standard deviation; "\*" and "\*\*" mean the covariate is significant at 0.1 and 0.05 level, respectively; "V" is the traffic volume; "A" is the shock wave area; and "P" is the platoon ratio.
multisite BHM\_GEV(1,1,0), respectively. The results of the comparison are presented in Table 5. The estimated crashes of four intersections from the multisite BHM\_GEV(1,1,0) are all within the 95% confidence intervals of the observed crashes, while the estimated crashes of only one intersection from the at-site model are within the 95% confidence intervals of the observed crashes. Moreover, the estimated crashes of three intersections from the multisite BHM\_GEV(1,1,0) are closer to the observed crashes than those from the at-site models. Additionally, the 95% confidence intervals of estimated crashes from the multisite BHM\_GEV(1,1,0) are narrower than those from the at-site Bayesian hierarchical models. This indicates that as compared to the at-site Bayesian hierarchical models, the multisite BHM\_GEV(1,1,0) distinctly improves crash estimation in terms of accuracy and precision. Such significant improvement is probably due to synthesizing the strengths across different sites for model estimation.
However, it is noted that the confidence intervals of estimated crashes from the multisite BHM\_GEV(1,1,0) are still much wider than the 95% Poisson intervals of observed crashes. As discussed in one previous study [\(Zheng et al., 2019a\)](#page-9-0), the large uncertainty in crash estimation based on the multisite Bayesian hierarchical model is mainly a result of the paucity of data. On one hand, the short-term observation period of three-four hours may not completely reflect the daily/weekly/seasonal variation of traffic in one year. On the other hand, any uncertainty in the estimated risk of crashes during the traffic conflict observation period will be magnified by at least *T/t* = 365 × 12*/*4 = 1095 times. From this point of view, the precision issue can be handled with a larger data sample collected from a longer observation period.
**Table 5** Crash estimate comparison between multisite model and at-site models.
| Site | Observed | crashes per year | per year | Multisite model
estimated crashes | At-site models
estimated crashes
per year | | |
|------------------------------|----------|------------------|----------|--------------------------------------|-------------------------------------------------|-------------|--|
| | Mean | 95%
C.I. | Mean | 95% C.I. | Mean | 95% C.I. | |
| 72 Ave & 128
St | 4.7 | [2.6,
7.8] | 6.7 | [0.0,505.5] | 8.3 | [0.0,647.6] | |
| 72 Ave & 132
St | 0.3 | [0.0,
1.9] | 1.9 | [0.0,174.9] | 1.1 | [0.0,176.2] | |
| 64 Ave & King
George Blvd | 6.0 | [3.6,
9.5] | 4.4 | [0.0,440.1] | 2.7 | [0.0,492.5] | |
| Fraser Hwy &
168 St | 3.3 | [1.6,
6.1] | 5.6 | [0.0,479.6] | 13.2 | [0.0,969.6] | |
## *4.5. Comparison of extreme value distribution and truncated normal distribution*
Given that the MADR follows a truncated normal distribution N (8.45, 1.402 ), with a lower limit of 4.23 m/s2 and an upper limit of 12.68 m/s2 ([Cunto and Saccomanno, 2008\)](#page-8-0), the risk of crash can also be calculated as:
$$R_{i} = \Pr(MADR \le DRAC_{i}) = \frac{\Phi(\frac{DRAC_{i} - \mu_{MADR}}{\sigma_{MADR}}) - \Phi(\frac{C_{MADR} - \mu_{MADR}}{\sigma_{MADR}})}{\Phi(\frac{U_{MADR} - \mu_{MADR}}{\sigma_{MADR}}) - \Phi(\frac{C_{MADR} - \mu_{MADR}}{\sigma_{MADR}})}$$
(15)
where, *μMADR* and *σMADR* are the mean and standard deviation of MADR distribution (here refers to a truncated normal distribution); *LMADR* and *UMADR* are the lower and upper limits of MADR; *Φ*(⋅) is the cumulative distribution function of the standard normal distribution.
Based on Eq. (15), the estimated crashes for each site were obtained and compared to the observed crashes and the estimated crashes from the multisite BHM\_GEV(1,1,0) using the values of MADR sampled from the truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68), respectively. The comparison results are shown in Table 6. The estimated crashes from the Bayesian hierarchical extreme value distribution are closer to the observed crashes than those from the truncated normal distribution. This confirms that the Bayesian hierarchical extreme value distribution is more suitable to fit the traffic conflict extremes than the truncated normal distribution.
**Table 6** Crash estimate comparison between Bayesian hierarchical extreme value and truncated normal distributions.
| Site | per year | Observed crashes | | Estimated crashes
from extreme value
distribution per year | Estimated
crashes from
truncated
normal
distribution per
year | |
|------------------------------|----------|------------------|------|------------------------------------------------------------------|------------------------------------------------------------------------------|-------------|
| | Mean | 95% C.
I. | Mean | 95% C.I. | Mean | 95%
C.I. |
| 72 Ave & 128 St | 4.7 | [2.6,
7.8] | 6.7 | [0.0,505.5] | 9.3 | – |
| 72 Ave & 132 St | 0.3 | [0.0,
1.9] | 1.9 | [0.0,174.9] | 4.3 | – |
| 64 Ave & King
George Blvd | 6.0 | [3.6,
9.5] | 4.4 | [0.0,440.1] | 3.8 | – |
| Fraser Hwy & 168
St | 3.3 | [1.6,
6.1] | 5.6 | [0.0,479.6] | 14.7 | – |
## **5. Conclusions**
This study compared the threshold determination methods (i.e., setting specific values or distributions of MADR) for DRAC-based crash estimation using Bayesian hierarchical extreme value modeling. The crash data and traffic conflicts defined by DRAC were collected from four signalized intersections in the city of Surrey, British Columbia, Canada. Four multisite Bayesian hierarchical models were developed, including a stationary model and three non-stationary models. For comparison purposes, the at-site non-stationary models for four sites were also developed separately. As well, the estimated crashes from a truncated normal distribution of MADR for four sites were calculated. The estimated crashes from the best fitted multisite Bayesian hierarchical model with five specific values (i.e., 7 m/s2 , 8 m/s2 , 8.5 m/s2 , 9.7 m/s2 , and 12 m/s2 ) of MADR and two truncated normal distributions (i. e., N(8.45, 1.42 )I(4.23, 12.68) and N(9.7, 1.32 )I(4.23, 12.68)) of MADR were compared to the observed crashes, respectively. With the best threshold determination approach, the estimated crashes from the best fitted multisite Bayesian hierarchical model were compared to those from the best fitted at-site non-stationary models and the truncated normal distribution, respectively.
The results show that, in terms of DRAC-based crash estimation accuracy, (i) the truncated normal distribution N(8.45, 1.42 )I(4.23, 12.68) of MADR outperforms both the truncated normal distribution N(9.7, 1.32 )I(4.23, 12.68) of MADR and the specific values (i.e., 7 m/s2 , 8 m/s2 , 8.5 m/s2 , 9.7 m/s2 , and 12 m/s2 ) of MADR; (ii) the Bayesian hierarchical extreme value distribution outperforms the truncated normal distribution. In terms of DRAC-based crash estimation accuracy and precision, the multisite Bayesian hierarchical model is better than the at-site Bayesian hierarchical models. Moreover, finding specific appropriate values of MADR to use for DRAC-based crash estimation at all intersections is difficult. Additionally, it is also found that during a signal cycle, more vehicles and greater shock wave area will generate higher crash risks, whereas a greater platoon ratio will produce a lower crash risk.
In conclusion, this study highlights the difference in DRAC-based crash estimation with different thresholds. It also demonstrates the usefulness and advantages of the Bayesian hierarchical extreme value modeling in crash estimation compared to the truncated normal distribution. However, some limitations in this study merit further investigations. Firstly, the traffic conflicts are observed during a short period and from only four signalized intersections. Data from more sites and longer observation periods can be used to confirm the results of this study. Furthermore, only rear-end traffic conflicts and crashes are used. Other types of traffic conflicts and crashes need to be investigated. Secondly, as to a truncated normal distribution of MADR, its parameters (e.g., lower limit, upper limits, mean, and standard error) need to be further investigated, since the MADR is subject to various vehicles, road, and environmental conditions. Third, other distributions of MADR can be investigated. Fourth, to avoid a collision, the driver of the following vehicle may change lane instead of braking heavily. The evasion maneuvers mainly dependent on the severity level of the conflict. Similar to the deceleration to safety time (Hupfer, 1997), the conflict level can be determined based on DRAC which can be investigated in future studies.
# **CRediT authorship contribution statement**
**Chuanyun Fu:** Conceptualization, Methodology, Investigation, Validation, Writing - original draft. **Tarek Sayed:** Conceptualization, Methodology, Data curation, Supervision, Writing - review & editing.
# **Declaration of Competing Interest**
None.
## **Acknowledgments**
This study was jointly supported by the China Scholarship Council (CSC) and the National Natural Science Foundation of China (Grant No. 71801182).
#### **References**
- [Archer, J., 2005. Indicators for Traffic Safety Assessment and Prediction and Their](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0005) [Application in Micro-simulation Modelling: A Study of Urban and Suburban](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0005) [Intersections. Doctoral Dissertation. Department of Infrastructure, Royal Institute of](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0005) [Technology, Sweden](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0005).
- [Ceunynck, T.D., 2017. Defining and Applying Surrogate Safety Measures and](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0010) [Behavioural Indicators Through Site-based Observations. Doctoral Dissertation.](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0010) [Department of Technology and Society, Lund University, Sweden.](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0010)
- [Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. Springer,](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0015) [London, UK](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0015).
- Cooley, D., Naveau, P., Jomelli, V., Rabatel, A., Grancher, D., 2006. A Bayesian hierarchical extreme value model for lichenometry. Environmentrics 17 (6), 555–574. [https://doi.org/10.1002/env.764.](https://doi.org/10.1002/env.764)
- [Cooper, D.F., Ferguson, N., 1976. Traffic studies at T-junctions](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0025) a conflict simulation [model. Traffic Eng. Control 17, 306](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0025)–309.
- Cunto, F., Saccomanno, F.F., 2008. Calibration and validation of simulated vehicle safety performance at signalized intersections. Accid. Anal. Prev. 40 (3), 1171–1179. [https://doi.org/10.1016/j.aap.2008.01.003.](https://doi.org/10.1016/j.aap.2008.01.003)
- El-Basyouny, K., Sayed, T., 2009. Urban arterial accident prediction models with spatial
- effects. Transp. Res. Rec. 2102, 27–33. [https://doi.org/10.3141/2102-04.](https://doi.org/10.3141/2102-04) El-Basyouny, K., Sayed, T., 2012a. Measuring direct and indirect treatment effects using safety performance intervention functions. Saf. Sci. 50 (4), 1125–1132. [https://doi.](https://doi.org/10.1016/j.ssci.2011.11.008)
- [org/10.1016/j.ssci.2011.11.008.](https://doi.org/10.1016/j.ssci.2011.11.008) El-Basyouny, K., Sayed, T., 2012b. Measuring safety treatment effects using full Bayes non-linear safety performance intervention functions. Accid. Anal. Prev. 45,
- 152–163. . Essa, M., Sayed, T., 2018. Traffic conflict models to evaluate the safety of signalized intersections at the cycle level. Transp. Res. Part C 89, 289–302. [https://doi.org/](https://doi.org/10.1016/j.trc.2018.02.014)
- [10.1016/j.trc.2018.02.014.](https://doi.org/10.1016/j.trc.2018.02.014) Essa, M., Sayed, T., 2019. Full Bayesian conflict-based models for real time safety evaluation of signalized intersections. Accid. Anal. Prev. 129, 367–381. [https://doi.](https://doi.org/10.1016/j.aap.2018.09.017)
- [org/10.1016/j.aap.2018.09.017.](https://doi.org/10.1016/j.aap.2018.09.017) Fu, C., Sayed, T., Zheng, L., 2020. Multivariate Bayesian hierarchical modeling of the non-stationary traffic conflict extremes for crash estimation. Anal. Methods Accid. Res. 28, 100135 .
- Gelman, A., Rubin, D.B., 1992. Inference from iterative simulation using multiple
- sequences. Stat. Sci. 7 (4), 457–472. . [Gettman, D., Head, L., 2003. Surrogate Safety Measures From Traffic Simulation Models,](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0070) [Final Report. U.S. Department of Transportation, Washington, D. C.](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0070)
- [Highway Capacity Manual, 2000. Transportation Research Board, Washington, D.C.](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0075). [Hupfer, C., 1997. In: Deceleration to Safety Time \(DST\) - a Useful Figure to Evaluate](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0080) [Traffic Safety? The 10th ICTCT Workshop. Lund, Sweden.](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0080)
- Lunn, D.J., Thomas, A., Best, N., Spiegelhalter, D., 2000. WinBUGS a Bayesian modelling framework: concepts, structure, and extensibility. Stat. Comput. 10, 325–337. [https://doi.org/10.1023/A:1008929526011.](https://doi.org/10.1023/A:1008929526011)
- [May, A.D., 1990. Traffic Flow Fundamentals. Prentice Hall, New Jersey, USA.](http://refhub.elsevier.com/S0001-4575(21)00082-8/sbref0090)
- Meng, Q., Weng, J., 2011. Evaluation of rear-end crash risk at work zone using work zone traffic data. Accid. Anal. Prev. 43 (4), 1291–1300. [https://doi.org/10.1016/j.](https://doi.org/10.1016/j.aap.2011.01.011) [aap.2011.01.011](https://doi.org/10.1016/j.aap.2011.01.011).
- Qu, X., Kuang, Y., Oh, E., Jin, S., 2014. Safety evaluation for expressways: a comparative study for macroscopic and microscopic indicators. Traffic Inj. Prev. 15 (1), 89–93. [https://doi.org/10.1080/15389588.2013.782400.](https://doi.org/10.1080/15389588.2013.782400)
- Sayed, T., Zein, S., 1999. Traffic conflict standards for intersections. Transp. Plan. Technol. 22 (4), 309–323..
- Songchitruksa, P., Tarko, A.P., 2006. The extreme value theory approach to safety estimation. Accid. Anal. Prev. 38 (4), 811–822. [https://doi.org/10.1016/j.](https://doi.org/10.1016/j.aap.2006.02.003)
- Spiegelhalter, D.J., Best, N.G., Carlin, B.P., Linde, A.V.D., 2002. Bayesian measures of model complexity and fit. J. R. Stat. Soc. Stat. Methodol. Ser. B 64, 583–639. [https://](https://doi.org/10.1111/1467-9868.00353) [doi.org/10.1111/1467-9868.00353.](https://doi.org/10.1111/1467-9868.00353)
- Wang, C., Stamatiadis, N., 2013. Surrogate safety measure for simulation-based conflict
metric. Accid. Anal. Prev. 71, 82–92..
- study. Transp. Res. Rec. 2386 (1), 72–80. . Wang, C., Stamatiadis, N., 2014. Evaluation of a simulation-based surrogate safety
- Wang, C., Xu, C., Dai, Y., 2019. A crash prediction method based on bivariate extreme value theory and video-based vehicle trajectory data. Accid. Anal. Prev. 123, 365–373. .
- Weng, J., Meng, Q., Yan, X., 2014. Analysis of work zone rear-end crash risk for different vehicle-following patterns. Accid. Anal. Prev. 72, 449–457. [https://doi.org/](https://doi.org/10.1016/j.aap.2014.08.003) [10.1016/j.aap.2014.08.003](https://doi.org/10.1016/j.aap.2014.08.003).
- Weng, J., Xue, S., Yang, Y., Yan, X., Qu, X., 2015. In-depth analysis of drivers' merging behavior and rear-end crash risks in work zone merging areas. Accid. Anal. Prev. 77, 51–61. .
- Zheng, L., Ismail, K., 2017. A generalized exponential link function to map a conflict indicator into severity index within safety continuum framework. Accid. Anal. Prev. 102, 23–30. [https://doi.org/10.1016/j.aap.2017.02.013.](https://doi.org/10.1016/j.aap.2017.02.013)
- Zheng, L., Sayed, T., 2019a. Comparison of traffic conflict indicators for crash estimation using peak over threshold approach. Transp. Res. Rec. 2673 (5), 493–502. [https://](https://doi.org/10.1177/0361198119841556) [doi.org/10.1177/0361198119841556.](https://doi.org/10.1177/0361198119841556)
- Zheng, L., Sayed, T., 2019b. Bayesian hierarchical modeling of traffic conflict extremes for crash estimation: a non-stationary peak over threshold approach. Anal. Methods Accid. Res. 24, 100106 .
- Zheng, L., Sayed, T., 2020. A bivariate Bayesian hierarchical extreme value model for traffic conflict-based crash estimation. Anal. Methods Accid. Res. 25, 100111 [https://doi.org/10.1016/j.amar.2020.100111.](https://doi.org/10.1016/j.amar.2020.100111)
- Zheng, L., Ismail, K., Meng, X., 2014. Freeway safety estimation using extreme value theory approaches: a comparative study. Accid. Anal. Prev. 62, 32–41. [https://doi.](https://doi.org/10.1016/j.aap.2013.09.006) [org/10.1016/j.aap.2013.09.006.](https://doi.org/10.1016/j.aap.2013.09.006)
- Zheng, L., Sayed, T., Tageldin, A., 2018. Before-after safety analysis using extreme value theory: a case of left-turn bay extension. Accid. Anal. Prev. 121, 258–267. [https://](https://doi.org/10.1016/j.aap.2018.09.023) [doi.org/10.1016/j.aap.2018.09.023.](https://doi.org/10.1016/j.aap.2018.09.023)
- Zheng, L., Sayed, T., Essa, M., 2019a. Bayesian hierarchical modeling of the nonstationary traffic conflict extremes for crash estimation. Anal. Methods Accid. Res. 23, 100100 [https://doi.org/10.1016/j.amar.2019.100100.](https://doi.org/10.1016/j.amar.2019.100100)
- Zheng, L., Sayed, T., Essa, M., 2019b. Validating the bivariate extreme value modeling approach for road safety estimation with different traffic conflict indicators. Accid. Anal. Prev. 123, 314–323. [https://doi.org/10.1016/j.aap.2018.12.007.](https://doi.org/10.1016/j.aap.2018.12.007)