category: literaturenote citekey: stipancicnetworkscreeninglargeurban2019 title: "Network screening for large urban road networks: Using GPS data and surrogate measures to model crash frequency and severity" authors: "Stipancic, Joshua; Miranda-Moreno, Luis; Saunier, Nicolas; Labbe, Aurélie" year: 2019 date: 2019-00-00 2019 doi: 10.1016/j.aap.2019.02.016 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2019.02.016" zotero_key: AFC39XZ5 zotero_storage: BABBFEMT collections: transpordi_planeerimine folder: 001_artiklid/Road Safety Indicators firstAuthor: "Stipancic, Joshua"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

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

Network screening for large urban road networks: Using GPS data and surrogate measures to model crash frequency and severity

Joshua Stipancica,⁎ , Luis Miranda-Morenob , Nicolas Saunierc , Aurélie Labbed

^a Department of Civil Engineering and Applied Mechanics, McGill University, Room 391, Macdonald Engineering Building, 817 Sherbrooke Street West, Montréal, Québec, H3A 0C3, Canada
^b Department of Civil Engineering and Applied Mechanics, McGill University, Room 268, Macdonald Engineering Building, 817 Sherbrooke Street West, Montréal, Québec, H3A 0C3, Canada
^c Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, C.P. 6079, succ. Centre-Ville, Montréal, Québec, H3C 3A7, Canada
^d Department of Decision Sciences, HEC Montréal, 3000 Chemin de la Côte-Sainte-Catherine, Montréal, Québec, H3T 2A7, Canada

ARTICLE INFO

Keywords: Network screening Bayesian models Crash modelling GPS data Surrogate safety

ABSTRACT

Crash frequency and injury severity are independent dimensions defining crash risk in road safety management and network screening. Traditional screening techniques model crashes using regression and historical crash data, making them intrinsically reactive. In response, surrogate measures of safety have become a popular proactive alternative. The purpose of this paper is to develop models for crash frequency and severity incorporating GPS-derived surrogate safety measures (SSMs) as predictive variables. SSMs based on vehicle manoeuvres and traffic flow were extracted from data collected in Quebec City. The mixed multivariate outcome is estimated using two models; a Full Bayes Spatial Negative Binomial model for crash frequency estimated using the Integrated Nested Laplace Approximation approach and a fractional Multinomial Logit model for crash severity. Model outcomes are combined to generate posterior expected crash frequency at each severity level and rank sites based on crash cost. The crash frequency model was accurate at the network scale, with the majority of proposed SSMs statistically significant at 95% confidence and the direction of their effect generally consistent with previous research. In the crash severity model, fewer variables were significant, yet the direction of the effect of all significant variables was again consistent with previous research. Correlations between rankings predicted by the mixed multivariate model and by the crash data were adequate for intersections (0.46) but were poorer for links (0.25). The ability to prioritize sites based on GPS data and SSMs rather than historical crash data represents a substantial contribution to the field of road safety.

1. Introduction

Crash frequency and injury severity are independent dimensions of road safety to consider when prioritizing sites for remediation. Network screening, the first task in road safety management, involves a low-cost examination of the entire road network to identify sites with potential for improvement (Hauer et al., 2002). The identified sites, termed hotspots, blackspots, or hazardous road locations may be areas where design or operation "create an increased risk of unforeseeable accidents" (Agerholm and Lahrmann, 2012). Once sites are identified, a detailed investigation, or diagnosis, is undertaken to identify contributing factors or potential causes of any safety issues and recommend efficient countermeasures. Most traditional screening techniques estimate crash frequency alone using regression techniques or other methods (Park and Sahaji, 2013; Anastasopoulos and Mannering, 2011) based on historical crash data and factors related to traffic exposure, geometry, and environment (Chang and Wang, 2006). As with Vision Zero, some have suggested that the most efficient way to improve safety is to reduce injuries rather than to reduce crashes (Chang and Wang, 2006). Therefore, the aim of safety management should be to reduce both frequency and severity of crashes (Aguero-Valverde and Jovanis, 2009). Common frequency models may be combined with injury severity models which are "conditional on the crash having occurred" (Anastasopoulos and Mannering, 2011), though these methods are data heavy, requiring detailed information for each distinct crash. Multivariate Bayesian models (Miaou and Song, 2005) have also been explored, though large-scale estimation is time-intensive.

Regardless of approach, most existing network screening techniques

E-mail addresses: joshua.stipancic@mail.mcgill.ca (J. Stipancic), luis.miranda-moreno@mcgill.ca (L. Miranda-Moreno), nicolas.saunier@polymtl.ca (N. Saunier), aurelie.labbe@hec.ca (A. Labbe).

^⁎ Corresponding author.

share critical drawbacks. By relying on ranking criteria derived from historical crash data (Huang et al., 2009) and aggregate exposure measures, most existing techniques are intrinsically reactive (Algerholm and Lahrmann, 2012), often requiring crashes to occur before hotspots can be identified (with one exception found in the Empirical Bayes technique). Crash-based methods require long collection periods to accumulate necessary data (Lee et al., 2006) and are subject to errors, omissions, and underreporting in crash databases (Kockelman and Kweon, 2002). Importantly, observed crashes themselves are not complete predictors of safety (crashes occurring today may not reveal where crashes will occur tomorrow). In response, surrogate safety measures (SSMs) have become a popular alternative to crash-based methods. Developed beginning in the 1960s, SSMs are any non-crash measures that are physically and predictably related to crashes (Tarko et al., 2009). Despite the promise of surrogate safety analysis (Tarko et al., 2009), most work to date has focussed on traffic conflicts and other interactions collected by trained field observers and, more recently, video sensors and computer vision techniques at individual sites. Notable exceptions include naturalistic driving studies from the SHRP 2 Safety Program in the United States (Hankey et al., 2016) and ITS field operational tests across Europe and the US (Benmimoun et al., 2013). Although these studies used probe vehicles, or mobile sensors, they relied on volunteers and/or test vehicles rather than regular drivers as demonstrated herein. Integrating SSMs into network screening requires a network-wide approach using a data source from which crash frequency and severity can be estimated for the entire road network. In the authors' previous work (Stipancic et al., 2018b), crash frequency models using SSMs were investigated, though severity was not considered.

This paper aims to address two outstanding issues related to screening large networks using SSMs. First, probe vehicles, which act "as moving sensors, continuously feeding information about traffic conditions" (El Faouzi et al., 2011), are perhaps the only data source feasible for large scale evaluation of SSMs across time and space. GPSenabled smartphones expand data collection possibilities to regular drivers (Jun et al., 2007), using a system that is inexpensive, userfriendly, and minimally impacts behaviour. Though SSMs based on vehicle manoeuvres (Stipancic et al., 2018a) and traffic flow (Stipancic et al., 2017a) have been successfully extracted from GPS smartphone data, they have not been incorporated into crash models. Second, though Bayesian methods are popular, estimation by Monte Carlo Markov Chain (MCMC) simulation can be computationally expensive for large networks and spatial models. Despite recent advances in Bayesian inference, the Integrated Nested Laplace Approximation (INLA) approach has been applied to road safety only very recently. The purpose of this paper is to propose a network screening approach based on SSMs derived from GPS data and using a mixed-multivariate model for crash frequency and severity. The specific objectives are to discuss the potential SSMs available from smartphone GPS data, calibrate a full Bayesian Spatial Latent Gaussian Model for crash frequency using the R-INLA program, model crash severity using a fractional Multinomial Logit Model (FMNL), and compare site rankings developed by the model and a traditional crash-based approach.

2. Literature review

Common techniques for crash frequency modelling include statistical count models such as Poisson (Mannering and Bhat, 2014), Negative Binomial (NB) (Lord et al., 2005), or Zero-inflated Poisson models (Mannering and Bhat, 2014). More advanced regression models incorporate random effects, multivariate outcomes, and hierarchical structures (Lord and Mannering, 2010). Conditional crash severity models require detailed environment, roadway, user, and vehicle information for crashes that have already occurred. Early binary models have evolved to include multiple discrete outcome models, whether unordered (multinomial and nested logit models) or ordered (ordered probit and logit models) (Yasmin and Eluru, 2013). Aggregate severity models use average geometric, traffic, and environmental data rather than detailed crash-specific data. Anastopoulos and Mannering (Anastasopoulos and Mannering, 2011) estimated proportions of crashes at each severity level using an aggregate model that was comparable to a conditional model in terms of identified hotspots (Anastasopoulos and Mannering, 2011).

In regression models, coefficients take fixed values. In Bayesian models, coefficients are defined by a probability distribution (Biangiardo and Cameletti, 2015). Complex Full Bayes (FB) models are typically estimated by assuming a prior distribution and iteratively updating using MCMC simulation (Rue et al., 2009). Examples can be found in El-Basyouny and Sayed (2009) and Aguero-Valverde and Jovanis (2008). Empircal Bayes (EB) models fix some parameters based on observed crash data (Jiang et al., 2014) instead of using hyperpriors. Though simpler to estimate, deviations between EB and FB estimates can be significant (Miaou and Lord, 2003). As with regression, Bayesian techniques have been extended to improve crash estimates, including the addition of temporal and spatial correlations which significantly improve accuracy (Quddus, 2008; Aguero-Valverde and Jovanis, 2008; Miaou and Lord, 2003). INLA has replaced MCMC for estimating univariate FB models in several studies such as Hu et al. (2013) and Serhiyenko et al. (2014). Crash severity can be integrated into FB models using either two-step or multivariate models. Wang et al. (2011) developed a two-step model, using an FB spatial model to estimate frequency and an unordered nominal response model to determine proportion by severity type. Multivariate Bayesian models were estimated by Miaou and Song (2005), Aguero-Valverde and Jovanis (2008), and Park and Lord (2007) to simultaneously estimate crash counts at several injury severity levels. Multivariate FB models may also be approximated using INLA (Serhiyenko et al., 2016).

SSMs have rarely been integrated into statistical models for network screening, though some studies have used conflicts to improve estimates of crash frequency using simulated (Ariza, 2011; Lorion, 2014) or observed data (Li et al., 2016). Other SSMs available from probe data, including individual driver manoeuvres of steering, braking, or accelerating (Ellison et al., 2013; Algerholm and Lahrmann, 2012; Bagdadi, 2019), and traffic flow SSMs of speed or speed variation (Moreno and Garcia, 2013; Boonsiripant, 2009) and congestion (Dias et al., 2009), have not been incorporated into screening models. When this paper was written, only exploratory work by Kluger (2017) could be found. Shortcomings in the existing literature generally fall at the intersections of crash modelling, Bayesian inference, and SSMs. First, although FB techniques are the most accurate and well-accepted approach for modelling road crashes, current MCMC simulation approaches are computationally expensive and time consuming. This study takes advantage of recent advances in Bayesian inference, namely INLA, as a state-of-the-art method to solve a complex problem in the field of road safety. Second, although crash models continue to improve, most existing approaches are crash-based. SSMs are under continuous development and can be extracted network-wide based on probe vehicle data. Besides earlier work by the authors (Stipancic et al., 2018b), almost no studies to date have incorporated SSMs into statistical crash models. This paper aims to address these gaps by developing a mixedmultivariate model for crash frequency and severity incorporating the INLA technique and SSMs extracted from the smartphones of regular drivers.

3. Methodology

The methodology for this study consists of several steps. First, data are collected and processed and SSMs are extracted. These steps have been covered in detail in previous studies (Stipancic et al., 2017b, 2018a; Stipancic et al., 2017a). The focus of this section is to describe the methods for modelling crash frequency and severity and the comparison of sites ranked using both the surrogate safety model and

Fig. 1. Methodological steps for crash modelling and site ranking.

traditional crash-based techniques. An overview of the methodological steps is presented in Fig. 1.

3.1. Data collection and processing

This project exploits three sources of data: smartphone GPS data, historical crash data, and digital map data. Each of these data sources must be processed before the development of the SSMs and crash model. A detailed description of data processing is provided by Stipancic et al. (2017b), but is covered briefly in the following paragraphs. Trips logged using a smartphone application are returned as a series of observations ( ) Oij for each user i , from j = 0 to the total the number of data points for trip i (typically one per second) containing information such as time, position, and speed ( ) vij . The raw GPS traces are processed using a commercially available map matching algorithm, TrackMatching (Marchal, 2015), which returns a latitude, longitude, and road link ID corrected according to the OpenStreetMap (OSM) road network (OpenStreetMap, 2015). In addition to positional noise, noise in the GPS measured speeds is reduced using a Savitzky-Golay filter, "a weighted moving average–based filter" suitable for time series with fixed and uniform intervals and limited discontinuities (Zaki et al., 2014), parameters for which were determined in a previous study (Stipancic et al., 2018a). The filter provides a filtered speed ( ) v ij and acceleration rate ( ) aij for each observation.

The use of second-by-second GPS data has the potential to introduce errors which must be carefully managed. As with all mobile sensor data, the GPS observations may not be well located, with errors in reported vehicle positions potentially biasing speed measurements and leading to false or missed detections of safety-critical events, though this is largely addressed through map matching and speed filtering. When several consecutive observations are missing (due to weak GPS signals), no safety-related information can be captured. However, unless missing observations are consistent and location-specific (tunnels, for example), then they should not systematically bias the results. One-second data frequency may not be fine enough to capture all desired events. For example, much of the existing literature uses accelerometer data (10 Hz) to measure acceleration. Although some events may be missed using GPS data alone, previous work has shown the strong correlation between actual crashes and GPS-measured SSMs (Stipancic et al., 2018a).

Map data is obtained from OSM to be consistent with the map matching results. The OSM road network is user-generated, leading to issues with the definition of some links. Ideally, each link in the network should connect adjacent intersections (Sioui and Morency, 2013). The network is redefined by splitting all links at each intersection and renaming the link according to the nodes at each end. The GPS data is similarly remapped according to the new network definition. As with GPS trips, historical crash data must be matched to the links and intersections of the network where they occur to obtain crash counts at each severity level for model estimation. Collision reports are first processed using a geocoding procedure developed previously by Burns et al. (2014) to obtain latitude and longitude from text-based fields such as address or intersection, with over 90% success.

Although the geo-location of historical crashes was largely successful, error may still be introduced when assigning crashes to specific links or intersections. In the existing literature, buffers are a popular method (i.e. creating a spatial buffer around each site and counting the contained crashes). However, this allows crashes to be counted at multiple sites, which introduces artificial autocorrelations in spatial modelling applications. In previous work (Stipancic et al., 2018b), several techniques were considered for avoiding this issue. In the end, crashes were assigned to links and intersections by using 50 m nonoverlapping buffers (Stipancic et al., 2018b). This technique was found to greatly reduce crash overcounting with only minor changes to modelling results compared to other techniques. However, it should be noted that this technique does permit some overcounting of crashes. More advanced crash assignment techniques will be considered in future work. Other potential sources of error include the map data itself, as OSM data may be unreliable, and results achieved using any two unique mapping sources may differ simply because of minor differences in the mapping data. Additionally, geolocating crashes is never 100% accurate as it relies on imperfect information present in crash databases. Care should be taken in the steps of geolocation, network definition, and crash assignment to ensure that errors are minimized relative to the sample size and that the errors do not introduce any systematic bias into the data set.

3.2. Extraction of surrogate safety measures

The extraction of SSMs from GPS smartphone data, and the strength of their relationships with crash frequency, is covered in detail in two previous studies (Stipancic et al., 2018a; Stipancic et al., 2017a). For measures to be considered valid, each site must contain two trips with at least two observations per trip. In terms of event-based measures, hard breaking events (HBEs) were considered for crash modelling, noting that braking is the most common evasive manoeuvre in urban areas (Algerholm and Lahrmann, 2012). Hard accelerating events (HAEs) were considered in previous studies but were omitted in the models due to high correlation with HBEs. For each observation, aij is compared to a braking threshold, determined previously to be −2 m/s2 (Stipancic et al., 2018a). Consecutive observations exceeding the threshold are considered a single HBE. The total number of HBEs for each site is obtained using the same buffers as for crashes, and then is normalized based on the total number of trips observed on each link or intersection to obtain a hard-braking rate.

As in previous work (Stipancic et al., 2017a), three traffic flow SSMs were also considered in the crash model. A measure of congestion, the congestion index ( ) CI as proposed by Dias et al. (2009) was calculated for each observation as

$$CI_{ij} = \frac{FFSL - v{ij}}{FFS_L}$$ if $FFSL > v'{ij}$
= 0 otherwise (1)

where FFSL is the free-flow speed for link L on which Oij falls, computed as the average speeds observed in the off-peak period (peak periods were defined as 6:00 to 10:00 AM and 3:00 to 7:00 PM). This yields CI values ranging from 0 (speed equal to the free flow speed) and 1 (speed is zero). Next, CIij is aggregated to each link using

$$CI{L} = \frac{\sum{i} \sum{j} CI{ij}}{N}$$ (2)

where N is the total number of observations on link L. The models use CI during the PM peak period, which was shown to have the strongest relationship with crash frequency and severity (Stipancic et al., 2017a). Average speed (*V*¯ ) is calculated for each link L as

$$\bar{V}{L} = \frac{\sum{l} \sum{j} v{lj}'}{N}$$ (3)

where N is the total number of observations on link L . The off-peak period was chosen for the average speed to avoid collinearity with CI (in this case, *V*¯ ^L is exactly equal to FFSL). In addition to the magnitude of speed, speed variation was considered using the coefficient of variation of speed ( ) CVS computed as

$$CVS{L} = \frac{\sigma(v{ij})}{\bar{V}_{L}}$$ (4)

where ( ) v ij is the standard deviation of all speeds used to compute *V*¯ L. CVS was found to be most strongly related to crash frequency and severity during the off-peak period (Stipancic et al., 2017a). All traffic flow SSMs were subsequently aggregated at the intersection level by taking the average of all intersecting links.

3.3. Modelling crash frequency

In Latent Gaussian Models (LGM), the response variable yi (in this case, total crash frequency) for each link or intersection i is assumed to follow a distribution from the exponential family (Normal, Poisson, or Binomial). This model can be written as a structured additive model, in which the mean of yi , noted µi , at site i is related to the predictors through a link function g () such that

$$g(\mu_i) = \eta_i = \beta0 + \sum{k=1}^{n_{\beta}} \betak z{ki} + \sum_{j=1}^{nf} f^{(j)}(u{ji}) + \varepsilon_i$$ (5)

where ⁰ is the intercept, ^k are the coefficients representing the linear effect of covariates zki (in this case, SSMs, trip counts, and roadway functional classification), f ( )^j are functions of covariates uji used to relax these linear relationships or introduce random effects, ⁱ is the unstructured error component, and ⁱ is the structured additive predictor (Rue et al., 2009). f ( )^j can take many forms, resulting in extreme flexibility incorporating aspects like spatial correlations critical for reducing biases in crash modelling (Latouche et al., 2007; Biangiardo and

Cameletti, 2015).

Previously (Stipancic et al., 2018b), an NB Spatial model was shown to accurately model crash frequency at the network scale. In this study, the spatial component is accounted for using a modified version of the Besag–York–Molliè (BYM) model presented by Besag et al. (1991). If outcome yi follows a Poisson distribution, the model is written as

$$y_i \sim Poisson(\mu_i)$$ (6)

where $$\log (\mu_i) = \beta0 + \sum{k=1}^{n_{\beta}} \betak z{ki} + u_i + v_i$$ (7)

where ^k are regression coefficients for fixed effects, zik are covariates, vi is a site-specific random effect modelled using an exchangeable correlation structure across sites (equal correlation between all pairs of considered sites), and ui is another site-specific random effect modelled as spatially structured (additional correlation between neighbouring sites). Several structures can be specified for u u u = … ( , , ) ¹ ⁿ . Here, a proper version of the conditional autoregressive (CAR) structure is chosen

$$u{l}|u{-l}, \tau, d^{-} \mathcal{N}\left(\frac{1}{d+n{l}} \sum{i = j} u{j}, \frac{1}{\tau(d+n{l})}\right)$$ (8)

where ni is the number of neighbours for link or intersection i, u ⁱ represents all members of u excluding the i th, i j ˜ indicates if links or intersections are neighbours, d > 0 controls the "properness", and > 0 is a scaling parameter (Biangiardo and Cameletti, 2015). Python scripts were used to generate graphs describing the network topology and define neighbours as links or intersections which are immediately adjacent to one another.

In order to reduce the time required to estimate LGM using MCMC approaches, this work proposes the implementation of the INLA technique for Bayesian inference, which has been shown to produce accurate approximations with a significant reduction in computational time (Schrodle and Held, 2011). The INLA approach was proposed by Rue et al. (2009) to perform Bayesian approximations on LGMs using a combination of Laplace approximations and numerical integration to estimate the posterior marginal of the latent field. Besides reducing potential computational time from the range of days to the range of hours, the approximation error in the INLA approach is nearly equal to the estimation error in typical MCMC methods (Rue et al., 2009). The necessity of incorporating a spatial component and its associated complexity makes INLA "particularly suitable in this context" (Biangiardo and Cameletti, 2015). A package for the programming language R has been developed to easily deploy INLA (R-INLA) (Schrodle and Held, 2011). Model performance was quantified using the deviance information criteria (DIC), mean square error (MSE), and the correlation coefficient (CORR).

3.4. Modelling crash severity

In order to integrate crash severity modelling into the above approach, the method for estimating a mixed multivariate model presented by Wang et al. (2011) is adopted. The mixed multivariate outcome is estimated using two models. First, crash counts, or frequency, are estimated using the presented Spatial LGM. Second, crash severity is integrated through a discrete choice model. At each site, the proportion of crashes at each severity level (fatal, major injury, and minor injury) are modelled using an FMNL estimated using quasi maximum likelihood. The probability for a crash at a given severity level m for link or intersection i is

$$Pl(m) = Pr\left(U{mi} \ge \suml U{li}\right) \tag{9}$$

$$Pl(m) = \frac{\exp(U{ml})}{\suml \exp(U{ll})}$$ (10)

where Umi is the utility associated with severity level m at site i , and Uli are the utilities for each severity level l . The utility for alternative m can be written as the sum of a deterministic component, Vmi, and a random error component, mi

$$U{mi} = V{mi} + \varepsilon_{mi} \tag{11}$$

where Vmi can be further decomposed as

$$V{ml} = \alpha{m0} + \sum{k=1}^{n{ci}} \alpha{mk} z{kl}$$ (12)

where ^m⁰ is the intercept, and mk are the coefficients of covariates zki (including SSMs, trip counts, and roadway functional classification) for severity level m. In the FMNL, rather than restricting the observed choice to 0 or 1, observations may vary as a fraction between 0 and 1, with the only constraint being that for any link and intersection, observations across all severity levels must sum to exactly 1. In this model, minor injury crash is taken as the reference level.

3.5. Site ranking

Once both models are estimated, the results must be combined so that sites can be ranked based on the estimated number of crashes at each severity level. Site ranking can be done in several ways by considering either frequency/occurrence of crashes or crash rates. These different criteria may be weighted by expected severities or consequences. To illustrate the developed methods and compare sites objectively, the method of Wang et al. (2011) is adopted. The authors propose the use of a decision parameter which combines the estimated crash counts systematically for each site under investigation Wang et al. (2011). The authors note that the choice of the decision parameter is context dependent, and can take several forms, including expected crash frequency, rate, or economic cost. This paper uses crash cost per vehicle-km to rank sites: the decision parameter for site i, i model , , is calculated first for the model as

$$\delta{l, model} = \frac{\sum{m} \mu{l} \cdot P{l}(m) \cdot C(m)}{t{l} \cdot l{l}}$$ (13)

where, as defined before, µi is the posterior mean expected number of crashes at site i and P mi( ) is the probability of a crash at severity level m at site i. Then, the product µ P m( ) ⁱ ⁱ yields the crash counts at each severity level and is multiplied by C m( ) , the relative cost of a crash at severity level m. In this case, C m( ) was chosen as the relative cost of a crash at each severity level according to Transport Canada (Economic Analysis Directorate of Transport Canada, 2008). Minor injury crashes were assigned cost of 1, major injury crashes a cost of 10, and fatalities a cost of 160. This is then divided by ti, the number of observed GPS traces at each site (used as a proxy of traffic exposure), and li, the length of the link (intersections are assigned a length of 1). A link or intersection with a higher value of ⁱ has a higher crash cost and is considered more hazardous and should therefore be given a higher priority in the network screening process. Next, the decision parameter is calculated for the crash data, i crash ,

$$\delta{l, crush} = \frac{\sum{m} Y{mi} \cdot C(m)}{l{l} \cdot l{i}} \text{ if } \sum{m} Y_{mi} > 0$$

$$= \frac{1}{100 \cdot l{i} \cdot l{i}} \text{ if } \sum{m} Y{mi} = 0$$ (14)

where Ymi is the number of observed crashes at severity level m and site i. This specification prevents a high number of duplicate rankings for sites with zero observed crashes, while ensuring sites with zero observed crashes are ranked below those with at least one observed crash. Recognizing that the crash-based rankings are somewhat naïve, they are presented for comparison only. As the model-based rankings rely on the posterior mean determined in the Bayesian analysis, the rankings are more robust, avoiding the regression-to-the-mean problem. The crash-based and modelled rankings are compared using Spearman's rank correlation. Hotspots identified using the crash data and the safety model are compared using the percent deviation, calculated as

% deviation= $$100 \cdot \left(1 - \frac{\kappa}{r}\right)$$ (15)

where is the number of locations found in the hotspot lists determined both by the crash-based and surrogate safety modelling methods and r is the number of hotspots considered (the list size) (Miranda-Moreno and Fu, 2006).

3.6. Model validation

While the calibrated models demonstrate the effect of each SSM on crash frequency and severity, the fit of the calibrated models overstates the predictive power of the model. In order to evaluate model predictions, the model is validated using a 10-fold cross validation. Sites are randomly split into 10 groups (or folds). The models are then estimated 10 times, with a different fold set aside for validation each time. Each calibrated model is then used to predict the number and severity of crashes at each site in the validation fold. Due to the way solitary sites (sites without neighbours) are handled by R-INLA, the crash counts at these sites are predicted using only the fixed effects of the model to prevent prediction of extreme mean crash counts. To demonstrate the predictive power, the sites are ranked again based on the predicted crash cost and compared to the crash-based site rankings using Spearman's rank correlation and percent deviation. Finally, maps are generated to observe the locations of the identified hotspots.

4. Results

4.1. Data description

GPS data was collected in Quebec City, Canada using the Mon Trajet smartphone application (City of Quebec, 2019). Originally developed for the City of Quebec by Brisk Synergies (2019), the application was installed voluntarily by drivers who anonymously logged their commute trips in the application, shown in Fig. 2. This study made use of a sample of the data, which contained over 4000 drivers and nearly 22,000 individual trips recorded between April 28 and May 18, 2014. Models were calibrated based on 11 years of crash data obtained from the Ministry of Transportation of Quebec (MTQ) for the period between 2000 and 2010. In total, 14,278 collisions involving at least one vehicle were identified during this time period. Although the pooling of crash data increases the number of observations for modelling, it precludes the ability of including temporal correlations. This will be explored using larger datasets in the future. Perhaps the greatest limitation of this work, and in fact most surrogate safety studies, is the fact that the temporal coverage of the SSM data and the crash data do not overlap. However, the assumption underlying the validity of surrogate safety methods is that the relationship between SSMs and safety should remain fairly stable, though more research is needed in this area.

4.2. Data exploration

There were sufficient collected data (at least two trips with two observations per trip) to calculate the proposed SSMs for 4623 links, of which 494 were classified as motorways, 2166 were classified as arterials and collectors (grouped from the primary, secondary, and tertiary classifications within OSM), and 1963 were classified as residential. In terms of intersections, 4429 had sufficient data for modelling. Intersections were classified by taking the highest classification for adjacent links, yielding 605 motorway intersections (mainly at the locations of access ramps), 1223 arterial/collector intersections, and 2601 residential intersections. Descriptive statistics for model variables are provided in Table 1. On average, all sites experience slightly more than one crash every two years, with the vast majority

Fig. 2. Smartphone application interfaces.

Table 1 Variables and Descriptive Statistics.

	Mean	Minimum	Maximum	Std. Dev.
Links
Crashes
Minor Injury	0.31	0.00	40.00	1.76
Major Injury	0.02	0.00	2.00	0.15
Fatal	2·10−3	0.00	1.00	0.05
Trips	128.71	2.00	2140.00	218.60
HBEs/Trip	0.08	0.00	2.58	0.14
Congestion Index	0.08	0.00	0.78	0.12
CVS	0.31	0.01	1.67	0.19
Average Speed (m/s)	11.32	1.05	30.46	5.26
Length (m)	149.96	4.49	3719.39	168.38
Intersections
Crashes
Minor Injury	2.70	0.00	51.00	5.60
Major Injury	0.15	0.00	4.00	0.48
Fatal	0.02	0.00	2.00	0.15
Trips	271.34	4.00	4282.00	417.87
HBEs/Trip	0.22	0.00	7.50	0.39
Congestion Index	0.09	0.00	0.78	0.11
CVS	0.31	0.01	1.40	0.16
Average Speed (m/s)	11.72	1.38	30.40	5.34

resulting in only minor injuries. The most crash prone links experience up to 43 crashes over the observation period, while up to 57 crashes occurred at the most extreme intersections. For an average link, about one in 12 trips experience an HBE, while for intersections, the number is one in five. Average CI was 0.08 and 0.09 for links and intersections respectively, average CVS was 0.31 for both, and average speed was approximately 42 km/h. Correlations between the considered SSM variables are provided in Table 2. In general, it would be beneficial to include more traditional geometric or signalization data, such as number of lanes/intersection legs and traffic signal presence. However, this information is not exhaustive or complete in the OSM database. More detailed geometric data could be included in future studies. However, as some geometric/signalization variables are correlated with functional class (the majority of urban arterial/collector intersections are signalized, for example), some of this information is captured by proxy.

4.3. Model calibration

4.3.1. Modelling crash frequency

The results of the Spatial NB model of crash frequency are presented in Table 3, and two promising observations are evident. First, most of the proposed SSMs are statistically significant at 95% confidence in both the link- and intersection-level models. Second, the direction of the effect of all variables (whether the posterior mean is positive or negative) is generally consistent with expectation and results from previous work. The posterior mean for the natural log of GPS trips, a proxy for exposure, is positive in both models. As exposure increases, crash frequency also increases. Previously, HBEs were shown to be positively correlated with crash frequency, and the posterior mean of HBEs/Trips is positive in the intersection model, providing more evidence for this positive correlation (Stipancic et al., 2018a). Additionally, this result supports several earlier studies finding positive relationships between braking and crashes (Bagdadi and Varhelyi, 2013; Jun et al., 2007; Agerholm and Lahrmann, 2012).

The signs for CVS and CI are also consistent with previous research by the authors (Stipancic et al., 2017a). Increased congestion increases crash frequency, a result demonstrated by several other authors as well (Dias et al., 2009; Golob et al., 2004; Shi and Abdel-Aty, 2015). The

Table 2 Correlations Between Considered SSM Variables.

	Total Crashes	Trips	HBEs/Trip	Congestion Index	Average Speed	CVS
Total Crashes	1.00	0.01	0.04	0.03	0.00	0.03
Trips		1.00	0.15	0.15	0.47	−0.12
HBEs/Trip			1.00	0.17	−0.07	0.22
Congestion Index				1.00	0.18	−0.11
Average Speed					1.00	−0.58
CVS						1.00

Table 3 Results of the Negative Binomial Spatial Model for Links and Intersections.

Explanatory variables	Links				Intersections
	mean	std dev	95% CI		mean	std dev	95% CI
Intercept	−15.83	0.87	−17.60	−14.22	−2.423	0.20	−2.811	−2.039
ln(Trips)	0.495	0.09	0.320	0.672	3.826	0.17	3.502	4.153
HBEs/Trip	−0.438	0.65	−1.735	0.806	1.131	0.10	0.940	1.328
Congestion Index	1.436	0.56	0.339	2.521	0.892	0.28	0.346	1.440
CVS	0.972	0.51	−0.044	1.978	0.905	0.26	0.402	1.408
Average Speed	−0.046	0.03	−0.106	0.014	−0.064	0.01	−0.085	−0.043
ln(Length)	2.110	0.12	1.881	2.354	N/A	N/A	N/A	N/A
Motorway	−3.037	0.47	−3.997	−2.135	−0.670	0.13	−0.921	−0.420
Arterial/Collector	0.848	0.19	0.485	1.220	0.449	0.06	0.323	0.576
Number of cases	4623				4429
DIC	3756.5				15920.7
MSE	2.5				19.3
CORR	0.63				0.67

Note: Variables significant at 95% confidence are bolded.

Fig. 3. Fitted values versus observed crashes for links (a) and intersections (b).

same is true for CVS, mirroring previously observed relationships between crash frequency and speed variation (Golob et al., 2004; Abdel-Aty and Pande, 2005; Moreno and Garcia, 2013). Sites with a higher average speed tend to have fewer crashes overall. This result seems counterintuitive, contradicting some earlier findings (Elvik, 2009), with Imprialou et al. (2016) attributing this negative relationship to linkaggregation bias. However, in this context, it is suggested that FFS captures other variables not considered, including geometry. In this case, links with a higher FFS are likely easier to navigate, and therefore fewer crashes should be expected. Increasing link length also increases crash frequency. All else being equal, motorways have fewer crashes than residential streets, while arterials and collectors tend to have a greater number.

Model fit is further demonstrated in Fig. 3. Model performance, as measured by MSE and CORR, is comparable between the models. Although the MSE is higher in the intersection level model, this is simply because the crash counts at intersections are higher overall (mean crash counts are 2.87 at intersections and only 0.33 on links). Correlations between observed crashes and fitted values are relatively high, which is a promising result. Based on the calibration data set, the selected covariates are able to explain 63% and 67% of the variation in the crash counts for links and intersections respectively. It is observed that the model tends to underestimate the number of expected crashes, as most of the data points fall below the diagonal line which represents ideal performance.

4.3.2. Modelling crash severity

Results of the FMNL model of crash severity are presented in Table 4. The trip and length variables are omitted as the volume/exposure and length are not typically associated with crash severity. Additionally, the arterial/collector variable had to be omitted from the link-level model in order to achieve convergence of the log-likelihood. Minor injury crashes were selected as the base case, so coefficients represent the change in utility for major injury and fatal crashes. Compared to the frequency model, fewer variables are significant at 95% confidence. However, considering those variables that are statistically significant, results are generally consistent with expectations based on previous work. For example, intersections with a higher number of HBEs/Trip are more likely to experience fatal collisions, confirming a positive correlation between braking and crash severity shown previously (Stipancic et al., 2018a). In contrast, a higher CI significantly reduces the chance of a fatal crash at links. This supports both previous research and intuition (Stipancic et al., 2017a). Although congestion is likely to increase crash frequency due to increased exposure, it is also likely to reduce severity through speed reduction. Both average speed and variation in speed are positively linked to fatal crashes.

4.3.3. Site ranking

Next, the results of both models are combined to rank sites based on the decision parameter, ⁱ, calculated using both crash data and the modelling results. Results are compared using Spearman's rho and by

Table 4 Results of the Fractional MNL for Links and Intersections.

	Links				Intersections
Explanatory variables	coeff	std err	95% CI		coeff	std err	95% CI
Major Injury
Intercept	−4.022	0.98	−5.939	−2.105	−2.835	0.42	−3.656	−2.014
HBEs/Trip	−0.504	2.23	−4.872	3.863	−0.234	0.18	−0.592	0.125
Congestion Index	−0.500	1.76	−3.941	2.942	−0.450	0.55	−1.531	0.631
CVS	0.124	1.17	−2.164	2.412	−0.322	0.59	−1.475	0.832
Average Speed	0.114	0.05	0.007	0.221	0.031	0.02	−0.017	0.079
Motorway	−1.543	0.84	−3.198	0.111	−0.890	0.36	−1.603	−0.178
Arterial/Collector	N/A	N/A	N/A	N/A	−0.054	0.14	−0.320	0.211
Fatal
Intercept	−11.17	2.81	−16.67	−5.662	−7.483	1.07	−9.573	−5.393
HBEs/Trip	−1.91	5.66	−13.00	9.18	0.512	0.14	0.245	0.779
Congestion Index	−9.393	4.15	−17.53	−1.253	−2.164	1.65	−5.408	1.080
CVS	5.729	2.38	1.072	10.385	1.938	0.88	0.213	3.663
Average Speed	0.366	0.16	0.047	0.685	0.197	0.06	0.079	0.314
Motorway	−1.723	2.43	−6.487	3.040	−1.555	0.76	−3.051	−0.058
Arterial/Collector	N/A	N/A	N/A	N/A	−0.862	0.34	−1.522	−0.201
Number of cases	453				2204
Log likelihood	−117.7				−574.4

Note: Variables significant at 95% confidence are bolded.

percent deviation. Fig. 4 compares the percent deviation for the crashbased and modelling approaches for hotspot lists of 100–1000, increasing in 100 site increments. Note, a lower percent deviation indicates superior performance. Percent deviation for links is around 50% regardless of the list size, while for intersections, percent deviation gradually decreases with increasing list size, from about 70% down to 40%. Based on the calibration data set, there is about 50% agreement between site rankings determined using the model and the crash data when considering the top 500 dangerous sites. Considering all links, Spearman's rho between the ranks based on crashes and the model is 0.47, while for intersections Spearman's rho was 0.63. The slightly better performance at the intersection level is again attributed to the larger sample size and overall higher observed crash counts which aid the calibration of the models. Additionally, model ranks are much less accurate for sites with few observed trips. Considering links with less than 100 observed trips, the average difference between model- and crash-based ranks was 1137, while for links with more than 500 observed trips, this number was 580.

The disagreement between the model and observed data can be ascribed to two sources. First, sites may be ranked higher by the crash data than by the model. This is generally not an issue for network screening if the site rankings are quite high (they would not be selected for remediation by either crash data or the developed model). However, there are several sites in both models with a high history of crashes which are not ranked highly by the surrogate model. For links, 42% of sites in the 90th percentile of observed crash frequency are not ranked in the top 1000 sites by the model (72 out of 170), while for intersections the number is 59% (256 out of 433). At these sites, the selected covariates are simply unable to capture the source of historical crashes. Additional variables or increased model flexibility may be able to capture the correlates of crashes at these few sites. Conversely, there are many sites which are expected to have a high crash cost based on the modelling results, which have zero (or very few) observed crashes (they are ranked much higher by the surrogate safety model than by the crash data). However, this does not necessarily mean that the model is performing poorly. These sites have the features of high-risk sites as determined by the model, even though they have no observed crashes. Therefore, these sites would be of interest because of their potential for future crashes. Although low-volume roads contribute disproportionally to percent deviation, they do not appear to bias the results (sites with less than 100 observed trips were equally likely to be ranked higher or lower by the model). More work is needed to determine which attributes contribute most to poor ranking. However, it is believed that an increased sample size (with more GPS observations per site) will aid in improving the percent deviation substantially.

Fig. 4. Percent deviation for hotspots generated by crash data and modelled using the calibration data.

Fig. 5. Predicted values versus observed crashes for links (a) and intersections (b).

4.4. Model validation

The final objective of this paper is to evaluate the model's prediction power using a 10-fold cross validation. The crash counts from the 10 folds predicted using the NB spatial model are compared to observed crashes in Fig. 5. For the link-level model, CORR decreased from 0.63 to 0.25 (a reduction of 0.39) while CORR for the intersection-level model decreased from 0.67 to 0.46 (a reduction of 0.21). Based on this figure, the prediction power is much stronger at the intersection level. Although the correlation at the link-level is still positive, the model struggles to predict crash counts. As before, the better performance for intersections is attributed to sample size. For sites without neighbours, crash counts are predicted using only the fixed effects of the model. This effectively reduces the MSE by eliminating extreme values otherwise predicted by R-INLA for solitary sites. Next, the predicted crash counts and severity proportions are again combined to calculate the decision parameter and rank the sites.

Fig. 6 shows the percent deviation between the list of hotspots generated by the crash data and safety model for various list sizes. When considering a small number of hotspots, the two methods deviate significantly from each other (approximately 85%). However, as the number of considered hotspots inceases, the deviation between the lists decreases (deviation is 55% for intersections and 65% for links when considering 1000 sites). Again, it should be noted that while some of this deviation is due to deficiencies in the model (truly dangerous sites are not captured by the model), some of this deviation is attributed to sites with potential for future crashes (sites that have not yet experienced a high number of crashes, but which are likely to in the future). As with the calibration data set, low-volume sites contribute disproprotionately to the deviation. A decrease in the goodness-of-fit is expected and observed for the predicted data. Spearman's rho for the link-level rankings decreased from 0.47 to 0.32 (a reduction of 0.15), while rho for the intersection-level rankings decreased from 0.63 to 0.45 (a reduction of 0.18). Although a 10-fold cross validation was used in this paper, in order to demonstrate the veracity of the SSMs of predictors of future crashes, future validation could aim to predict crashes for the year after the GPS data was collected.

Finally, the spatial distribution of the predicted rankings is illustrated in Fig. 7. Considering the better performing intersection-level model, there is a concentration of sites with a high predicted ranking near the downtown core of the city. This mirrors expectations of more dangerous sites being located in denser areas with higher traffic volumes. Second, the freeways leading to and from the city centre are generally receive a low ranking. This demonstrates the effect of controlling for exposure using the number of observed trips. Although freeways are expected to have a high number of crashes, because they also experience the highest volumes of all considered sites, the crash cost per vehicle-km (or the crash rate) is expected to be lower than other sites. However, considering the link-level results, the dangerous sites are dispersed throug the network. There is no observable spatial trend between model and crash-based rankings.

Fig. 6. Percent deviation for hotspots generated by crash data and modelled using the validation data.

Fig. 7. Sites ranked by the model from highest (red) to lowest (green) and difference between model and crash rankings, ranked higher by the model (blue) and higher by crashes (green) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

5. Conclusions

The purpose of this paper was to develop mixed-multivariate model capable of estimating the posterior of the expected crash frequency and severity at the link and intersection level across a large urban road network, using SSMs. GPS data and available geometric information are used to derive the contributing factors. This paper uses two models, combining a Bayesian Spatial LGM to model crash frequency and an FMNL model to estimate crash severity. From this, crash frequency, severity, and cost measures can be derived for use in network screening. A large quantity of GPS smartphone data was processed using both map-matching and speed filtering algorithms to reduce signal noise. Both event-based (HBEs) and traffic flow SSMs (congestion, speed, and speed variation) were extracted and integrated into the multivariate approach. A Bayesian Spatial NB model was estimated using the INLA technique and was shown to model crash frequency well at the network scale. Critically, most of the proposed SSMs were observed to be statistically significant at 95% confidence in both the link- and intersection-level models. Not only were these variables significant, but the direction of their effect was consistent with previous research. Namely, HBEs, congestion, and speed variation were all positively correlated with crash frequency, while average speed was negatively correlated. This result demonstrates the positive relationship between crash frequency and the proposed SSMs. Model fit, as measured by MSE and correlation coefficient, was comparable for both link- and intersection level models, with correlation exceeding 0.60 using the calibration data set.

In the second stage, crash severity was accounted for using a discrete choice framework estimated using quasi maximum likelihood. The probability for crashes at three distinct severity levels was estimated using an FMNL model. In this model, fewer variables were significant compared to the frequency model. Yet, the direction of the effect of all significant variables was again consistent with previous research or intuition. Namely, congestion tends to reduce the likelihood of fatal crashes, while increases in HBEs, speed, and speed variation are linked with an increase in crash severity. Each link and intersection was ranked according to the cash cost per vehicle-km, calculated first using the historical crash data and second using the modelling results. The ranked lists generated by the mixed multivariate model and the ranked lists based on crash data had a correlation of 0.47 for links and 0.63 for intersections. Percent deviation was approximately 50% when considering 500 hotspots. The model was finally validated using a 10-fold cross validation approach, and site rankings were again generated based on the predicted crash counts and severity proportions. For the validation data set, Spearman's rho was observed to decrease by 0.15 for links and 0.18 for intersections. This result shows that the predicted site rankings differ from the fitted site rankings by less than 20%. Although the intersection level model provided promising results, identifying many hotspots near the city centre, the link level model had poorer prediction power. Still, the selected covariates are able to explain 32% of the variation in crash cost for links and 45% of the variation for intersections. Again, not all of the of the discrepancy is considered a negative, as this approach is able to identify sites with a high potential for crashes (based on the selected SSMs) even if those sites have not historically experienced crashes.

Several areas have been identified for future work. First, additional model complexity, including spatial and/or temporal correlations and increased flexibility through interaction variables or more complex model structures, could help to further improve the accuracy of the model. Additional variables, such as geometric and signalization variables not captured in the OSM database, could also be considered by adopting a different mapping data source. Second, other SSMs, whether extracted from GPS or other in-vehicle sensors, could be incorporated as additional covariates in either model. Other measures extracted from probe vehicles may be able to capture other factors related to crash frequency or severity not considered in the existing methods. Third, considering the impact of the crash assignment method on modelling results is required. Geo-locating crashes on the road network based only on text from police reports is a well-known problem. Although several assignment algorithms were previously tested (Stipancic et al., 2017b), more work is needed to determine the impact of assignment including model sensitivity to of crash mis-location. Fourth, additional validation using larger datasets is expected to improve the accuracy of the results, especially for low-volume sites. Rather than relying on 10-fold cross validation, future studies could calibrate models based on historical data, and then predict crashes for the year following the GPS data collection. The greatest strength of this proposed approach is that, as GPS data is continually collected within a city, these rankings can be continuously updated. Although historical crash data is still useful and needed to calibrate these crash models, changes to site rankings can occur before additional crashes occur at the most dangerous sites, and this proactive approach has the potential to reduce road traffic crashes, injuries, and fatalities. If the developed techniques are proven reliable in other cities, prioritizing sites based on GPS data and SSMs rather than historical crash data could represent a substantial contribution to the field of road safety by allowing for a proactive approach to network screening.