Halaman utama Jelajahi Bantuan
Masuk
ardo
/
Obsidian
1
0
Fork 0
Berkas Masalah 0 Permintaan Tarik 0 Wiki
Pohon: 49b8810ad7
Ranting Tag
Obsidian
/
Zotero
/
001_artiklid
/
Road Safety Indicators
/
stipancicsurrogatesafetynetworkscreening2018.md

stipancicsurrogatesafetynetworkscreening2018.md 78 KB
Riwayat Mentahan

category: literaturenote citekey: stipancicsurrogatesafetynetworkscreening2018 title: "Surrogate safety and network screening: Modelling crash frequency using GPS travel data and latent Gaussian Spatial Models" authors: "Stipancic, Joshua; Miranda-Moreno, Luis; Saunier, Nicolas; Labbe, Aurelie" year: 2018 date: 2018-00-00 2018 doi: 10.1016/j.aap.2018.07.013 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2018.07.013" zotero_key: U8MBJ2EC zotero_storage: ZMFTRSVZ 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

Surrogate safety and network screening: Modelling crash frequency using GPS travel data and latent Gaussian Spatial Models

Joshua Stipancica,⁎ , Luis Miranda-Morenob , Nicolas Saunierc , Aurelie Labbed

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

ARTICLE INFO

Keywords: Crash frequency GPS Spatial Full Bayes Latent Gaussian Model INLA

ABSTRACT

Improving road safety requires accurate network screening methods to identify and prioritize sites in order to maximize the effectiveness of implemented countermeasures. In screening, hotspots are commonly identified using statistical models and ranking criteria derived from observed crash data. However, collision databases are subject to errors, omissions, and underreporting. More importantly, crash-based methods are reactive and require years of crash data. With the arrival of new technologies including Global Positioning System (GPS) trajectory data, proactive surrogate safety methods have gained popularity as an alternative approach for screening. GPS-enabled smartphones can collect reliable and spatio-temporally rich driving data from regular drivers using an inexpensive, simple, and user-friendly tool. However, few studies to date have analyzed large volumes of smartphone GPS data and considered surrogate-safety modelling techniques for network screening. The purpose of this paper is to propose a surrogate safety screening approach based on smartphone GPS data and a Full Bayesian modelling framework. After processing crash data and GPS data collected in Quebec City, Canada, several surrogate safety measures (SSMs), including vehicle manoeuvres (hard braking) and measures of traffic flow (congestion, average speed, and speed variation), were extracted. Then, spatial crash frequency models incorporating the extracted SSMs were proposed and validated. A Latent Gaussian Spatial Model was estimated using the Integrated Nested Laplace Approximation (INLA) technique. While the INLA Negative Binomial models outperformed alternative models, incorporating spatial correlations provided the greatest improvement in model fit. Relationships between SSMs and crash frequency established in previous studies were generally supported by the modelling results. For example, hard braking, congestion, and speed variation were all positively linked to crash counts at the intersection level. Network screening based on SSMs presents a substantial contribution to the field of road safety and works towards the elimination of crash data in evaluation and monitoring.

1. Introduction

Network screening is the low-cost examination of an entire road network to identify a smaller subgroup of sites for detailed investigation or site diagnosis (Hauer et al., 2002). This smaller subgroup (hotspots, blackspots, hazardous road locations, or high risk sites) is expected to include locations where design or operation "create an increased risk of unforeseeable accidents" (Algerholm and Lahrmann, 2012) with potential for crash reduction (Aguero-Valverde and Jovanis, 2009). Considering parties involved in road safety management have limited budgets, sites should be identified and prioritized to maximize the efficiency of implemented countermeasures, the specifics of which are determined in the diagnosis phase. Methods for network screening commonly use statistical (regression) models or safety performance functions and Bayesian statistics to estimate the expected number of crashes at each location of interest in the road network (Park and Sahaji, 2013). Through these techniques, the risk factors contributing to crash occurrence can also be uncovered (Chang and Wang, 2006). In addition to the expected number of crashes, various risk measures can be derived including the posterior probability of excess and posterior of

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.

ranks among others (Miranda-Moreno et al., 2007; Heydecker and Wu, 2001).

In network screening, safety must be objectively quantified. Most existing techniques use posterior ranking criteria based on Bayesian methods (Huang et al., 2009). Using statistical models, the relationships between attributes of traffic, geometry, environment, and driver (Abdel-Aty and Pande, 2005) and crash frequency and severity is first established (Lu, 2007). Then, Bayesian posterior analysis is used to quantify risk for each location. Though popular, methods derived from crash data are subject to errors in collision databases and are sensitive to underreporting (Kockelman and Kweon, 2002). Furthermore, as traffic collisions are relatively rare events (the low mean problem), long collection periods are required to accumulate sufficient crash data for analysis (Lee et al., 2006; Lord and Miranda-Moreno, 2008). For this reason, crash-based network screening cannot be carried out continuously, but is often performed periodically (for example, once every few years) so that crashes can accrue and databases can be updated. This highlights what is perhaps the most critical flaw of existing models. Crash-based methods are reactive, requiring crashes to occur before hazardous sites are identified and improvements are made (Algerholm and Lahrmann, 2012). For these reasons, alternative screening methods would be valuable for identifying high risk sites more quickly, more accurately, and with limited reliance on crash data (Cafiso and Di Graziano, 2011).

An alternative screening method requires an alternative data source from which crash risk can be constantly and systematically estimated throughout the network. Instrumented, or probe, vehicles that act "as moving sensors, continuously feeding information about traffic conditions" (El Faouzi et al., 2011) are an important source of such data. Vehicles instrumented with Global Position Systems (GPS) or other sensor types can assist in reducing dependence on crash databases by supporting the development of screening methods based on surrogate safety measures (SSMs) rather than collision statistics, a concept previously explored by Strauss et al. (2015) using instrumented bicycles. SSMs are any non-crash measures that are predictably related to crashes (Tarko et al., 2009). Although such methods rely on crash data for calibration, the application of the developed models to continuously monitor safety depends only on the input probe data that is continuously available. With the proliferation of GPS-enabled smartphones, which themselves contain many of the same sensors used for instrumenting vehicles, large volumes of reliable and spatio-temporally rich driving data can now be collected unobtrusively from regular drivers (Jun et al., 2007) in crashes, near crashes, and under normal conditions (Bagdadi, 2013; Wu and Jovanis, 2013). Smartphones are inexpensive and user-friendly, minimally impact behaviour, eliminate the need for external sensors (Eren et al., 2012; Johnson and Trivedi, 2011), take advantage of widespread technology, and exploit existing communication infrastructure (Herrera et al., 2010).

Despite advances in modelling techniques and data collection technologies, some areas of interest remain overlooked. Few studies have analyzed large volumes of GPS probe vehicle data collected from the smartphones of regular drivers across a large urban road network. Other than previous studies by the authors (Stipancic et al., 2018, 2017a, 2017b), even fewer have considered the link between GPS-derived SSMs and large volumes of historical crash data at the network level. Finally, few studies to date have considered advanced modelling techniques, including recent developments in Bayesian inference and spatial models, for screening in road networks using GPS-derived SSMs. The purpose of this paper is to propose a surrogate safety screening approach based on smartphone GPS data and a Full Bayesian modelling framework. Various SSMs, including vehicle manoeuvres and traffic flow measures, are first extracted. A Full Bayesian Spatial Latent Gaussian Model (LGM) is then developed with SSMs as covariates using the Integrated Nested Laplace Approximation (INLA) technique. Finally the model is validated on an independent data set.

2. Literature review

Bayesian techniques have been widely used within the field of transportation for modelling the likelihood of crashes and their frequency throughout road networks. Empirical Bayes (EB) models were popular in the 1980s and 1990s. In EB models, the probability of a crash is determined, in part, by using observed historical crash data (Jiang et al., 2014). Hauer (1992) described the EB process, noting that the safety of a site is described by both its characteristics and historical crash record, and presented applications of the model for estimating crashes in the US and Ontario, Canada. Mountain et al. (1996) applied the EB technique to a series of at grade crossings in the UK, showing an improvement over naïve regression models. Later, Miaou and Lord (2003) compared EB and Full Bayes (FB) techniques for the specific application of crash modelling, noting that EB estimates deviated from the FB estimates, and that those deviations could become significant for some data sets. FB techniques for complex problems (such as nonconjugate models including LGMs) typically determine the posterior distributions by first assuming a prior distribution and then iteratively computing and updating the posterior marginal using a Monte Carlo Markov Chain (MCMC) simulation. Many examples of complex FB models can be found in the existing literature, including models which incorporate random effects (Jiang et al., 2014; El-Basyouny and Sayed, 2009) and/or spatial correlations (Song et al., 2006; Miaou and Song, 2005). Of most interest to this study are the relatively rare examples of INLA used to approximate FB models of crash risk. For example, Hu et al. (2013) studied the patterns of highway crashes over time using a dynamic generalized linear model approximated using INLA, and found that safety was improving over time. Similarly, Serhiyenko et al. (2014) used INLA to approximate Bayesian models describing temporal trends in pedestrian crashes. Serhiyenko et al. (2016) also developed a multivariate Poisson Lognormal model estimated using INLA to model various types of crash counts.

Several studies have attempted to extract SSMs from probe vehicle data. Event-based techniques identify individual driver manoeuvres including steering, braking, or accelerating (Dingus et al., 2006). Fazeen et al. (2012) used smartphone accelerometer data to classify 'safe' accelerations and decelerations from 'unsafe' ones, though no evidence was provided demonstrating 'unsafe' behaviour led to increased risk. Jun et al. (2007) analyzed the relationship between temporal-spatial driving behaviour activity and likelihood of crash involvement, finding that drivers involved in crashes tended to travel longer distances at higher speeds and "engaged in hard deceleration events" more frequently. Algerholm and Lahrmann (2012) correlated jerk and crash occurrence both across drivers and across sites. Using GPS, accelerometer, radar, and self-reported collision data, Bagdadi (2013) proposed a jerk-based surrogate measure that correctly identified 86% of near misses. Traffic flow techniques rely on aggregate volume, speed, and density to measure risk (Yan et al., 2008). Though speed, flow, and variation in speed and flow have been suggested as potential SSMs in several studies (Abdel-Aty and Pande, 2005), traffic flow SSMs have rarely been studied using GPS data. Speed profiles from GPS devices were considered by Moreno and Garcia (2013) and Boonsiripant (2009).

Despite previous work on crash frequency modelling and extracting SSMs from GPS data, several shortcomings remain in the existing literature. Although many methods for extracting and analyzing SSMs have been proposed, few studies have extracted SSMs from instrumented vehicles, and very few have extracted such measures from smartphone GPS travel data alone. Few studies have used large amounts of data and validated the proposed SSMs using crash data. Instead, results are often compared to self-reported safety data or to near misses. More effort is needed to compare any proposed SSM with a reasonable amount of historical crash data to demonstrate that such a relationship exists. In terms of crash modelling, although more complex models continue to improve estimates of road traffic crashes, very few studies to date have incorporated SSMs into statistical models of crash frequency. Finally, although Bayesian modelling is a well-accepted approach for crash modelling, current estimation methods are computationally expensive and time consuming. This study takes advantage of recent advances in Bayesian inference, namely INLA, as a state-of-theart method to solve a complex problem in the field of road safety. This paper aims to combine SSMs based on vehicle manoeuvres and traffic flow extracted from GPS smartphone data and network screening models to minimize dependence on crash data in the network screening process.

3. Methodology

The study methodology consists of three steps, namely; i) collection, integration, and processing of the various data sources; ii) extraction of SSMs from smartphone GPS data; and iii) development, calibration, and validation of a crash frequency model. Additional details are provided in the following sections.

3.1. Data collection, integration, and processing

3.1.1. Smartphone GPS data

Collection and processing of smartphone GPS data is covered in detail by Stipancic et al. (2017a, 2017b), and the process is illustrated inFig. 1. Studied smartphone users are assumed to be representative of all drivers. While smartphone users are not truly representative, collecting data from the smartphones of regular drivers represents the least biased method currently feasible for collecting large volumes of GPS data, especially compared to taxis or other fleet vehicles.

For each logged trip, i, GPS travel data is returned as a series of observations (Oij), with j = … 0, ,ni where ni is the number of data points for trip i. These observations include the unique coordinate identifier (cij), datetime ( ) tij , latitude (xij), longitude (y )ij , altitude ( ) zij , and speed (vij). Time between consecutive observations is typically between 1 and 2 s. To reduce positional noise, the trip data is map matched to the OpenStreetMap (OSM) road network (OpenStreetMap, 2015) using TrackMatching, a commercially available web-based map matching service (Marchal, 2015), which returns a new latitude and longitude, ′ xij and ′ yij , corresponding to a specific OSM link ID, lij, and the source, sij, and destination nodes, dij. Next, a Savitzky-Golay filter is applied to reduce noise in the GPS measured speeds. This filter is "a weighted moving average–based filter" suitable for time series' with fixed and uniform intervals and with limited discontinuities in the data (Zaki et al., 2014). Optimal parameters for the filter were determined in a previous study (Stipancic et al., 2018) and were selected to maximize the correlation strength between crashes and the extracted SSMs. The filter is applied to speeds vij to yield filtered speeds ′ vij and the acceleration rate, aij, for every observation, yielding the final data structure below.

$$trip{i} = \begin{cases} O{i0} \ O{i1} \ \vdots \ O{ip} \ \vdots \ O{in{i}} \end{cases} = \begin{cases} i, \ c{i0}, \ t{i0}, \ x{i0}, \ y{i0}, \ z{i0}, \ v{i0} \ i, \ c{i1}, \ t{i1}, \ x{i1}, \ y{i1}, \ z{i1}, \ v{i1} \ \vdots \ i, \ c{ip}, \ t{ij}, \ x{ij}, \ y{ij}, \ z{ij}, \ v{ij} \ \vdots \ i, \ c{in{i}}, \ t{in{i}}, \ x{in{i}}, \ y{in{i}}, \ z{in{i}}, \ v{in{i}} \end{cases}$$

$$\Rightarrow \begin{cases} i, \ c{i0}, \ t{i0}, \ x{i1}', \ y{i0}', \ z{i0}, \ v{i0}', \ a{i0}, \ t{i0}, \ s{i0}, \ d{i0} \ \vdots \ i, \ c{i1}, \ t{i1}, \ x{i1}', \ y{i1}', \ z{i1}, \ v{i1}', \ a{i1}, \ t{i1}, \ s{i1}, \ d{i1} \ \vdots \ i, \ c{ij}, \ t{ij}, \ x{ij}', \ y{ij}', \ z{ij}, \ v{ij}', \ a{ij}, \ t{ij}, \ s{ij}, \ d{ij} \ \vdots \ i, \ c{in{i}}, \ t{in{i}}, \ x{in{i}}', \ y{in{i}}', \ z{in{i}}, \ v{in{i}}', \ v{in{i}}', \ a{in{i}}, \ t{in{i}}, \ s{in{i}}, \ d{in{i}} \end{cases}$$

3.1.2. Network map data

Ideally, each link in the network should connect adjacent intersections (Sioui and Morency, 2013). As the OSM road network is generated non-systematically by users, OSM links often connect multiple intersections, and should be redefined. Redefining the network requires the following steps, which can be completed in any GIS software environment. The result is further illustrated in Fig. 2.

1 Identify all nodes that represent an intersection in the road network.

Fig. 1. Collection and processing of smartphone-collected GPS data.

Fig. 2. Redefinition of OSM Links.

  • 2 Split the road network at the identified nodes.
  • 3 Rename each link according to its original ID and the nodes on either end of the link.
  • 4 Remap the GPS observations to the new network.

3.1.3. Collision data

The final step in data processing is assigning crash data to the links and intersections in the newly redefined road network to obtain crash counts at each link and intersection for model calibration and validation. Some collision reports may contain a latitude and longitude defining the location of the crash. However, for cases where coordinates are not provided, are provided inconsistently, or are inaccurate, a geocoding procedure may be used to determine the coordinates from text based fields, including address or intersection. This project made use of a geocoding procedure developed previously by Burns et al. (2014). Even with a latitude and longitude for each collision, crashes do not fall directly onto the network. In previous studies (Stipancic et al., 2018), buffers were used to assign crashes to links and intersections. However, using buffers allows for the over-counting of crashes, which can artificially create spatial autocorrelations. One option to solve this problem is to assign crashes to the nearest link or intersection in the network. However, this method allows for crashes to be mis-located, which can have serious implications in the modelling results. Instead, this study uses 50 m buffers around links and intersections, but does not allow link and intersection buffers to overlap. Although this method still allows for a small amount of overcounting (less than 20% at the link-level in the considered dataset), it improves the likelihood that crashes are assigned to the sites where they actually occurred. Though both methods have benefits and limitations, preliminary testing showed that either method provided substantially similar results.

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., 2018, 2017a, 2017b). The SSMs incorporated into the model are defined below.

3.2.1. Decelerations and accelerations

Deceleration is the most common evasive manoeuvre in urban areas (Algerholm and Lahrmann, 2012). Most studies focussed on deceleration have used accelerometers to calculate jerk as an SSM (Bagdadi, 2013; Algerholm and Lahrmann, 2012). Though the resolution of GPS data is too coarse for calculations of jerk, a simple deceleration threshold has been suggested to define hard braking events (HBEs) (Fazeen et al., 2012). Some have suggested that hard accelerations (HAEs) may also be good predictors of safety (Wu and Jovanis, 2013; Algerholm and Lahrmann, 2012). Each aij is compared to a braking threshold, amin, and acceleration threshold, amax , to determine which observations correspond to HBEs or HAEs. The status of Oij is determined using the following logic. For each series of consecutive negative (or conversely, positive) accelerations, ldec (lacc), the minimum (maximum) value is obtained. If this value is inferior (superior) to the threshold, min( ) ldec < amin (max( ) l a acc ma > x ), then consider that observation an HBE (HAE). This algorithm is illustrated in Fig. 3. Based on previous work (Stipancic et al., 2018), a threshold of −2 m/s2 was chosen for amin to maximize the correlation between events and crashes. The total number of HBEs was divided by the number of trips on each link or intersection. HAEs were omitted from the model because they are highly correlated with HBEs.

3.2.2. Congestion index

Though several congestion measures based on travel time have been proposed (Palen, 1997; Coifman, 2002; Falcocchio and Levinson, 2015), link travel time is dependent on position which is altered during map matching. Therefore, a congestion measure based on speed is preferred. Dias et al. (2009) proposed CI calculated as:

$$CI = \frac{\text{free flow speed} - \text{actual speed}}{\text{free flow speed}} \quad \text{if } CI > 0$$

$$= 0 \quad \text{if } CI \le 0 \quad (1)$$

This formulation yields CI values ranging from 0 (speed equal to the free flow speed) and 1 (speed is zero). As congestion is generally constrained to the AM and PM peak periods, speeds observed outside of these times can be used as estimates for the free flow speed. For each link, Lij, the free flow speed, FFSLij , is simply the average of all observed speeds falling on link Lij outside of the peak periods, defined as 6:00 to 10:00 AM and 3:00 to 7:00 PM. Next, CI for every observation can be computed according to:

$$CI{ij} = \frac{FFS{Lij} - v{ij}}{FFS{Lij}} \quad \text{if } FFS{Lij} > v'{ij}$$

$$= 0 \quad \text{otherwise}$$ (2)

where CIij is the congestion index for observation Oij, FFSLij is the free flow speed on link Lij, and *v*′ij is the filtered speed.

Finally, CI was calculated for each link by taking the average of each observed CIij:

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

Fig. 3. Algorithm for extracting vehicle manoeuvres from GPS trip data.

where N is the count of all observations on link Lij. As in previous work (Stipancic et al., 2017a, 2017b), filters were used to eliminate links with too few observations for a valid calculation (less than two trips with two observations per trip). Preliminary results indicated that CI during the PM peak had the strongest relationship with crash frequency and severity (Stipancic et al., 2017a, 2017b).

3.2.3. Average speed

Traditionally, travel speed was believed by many to be an indicator of crash risk at the link level (although this relationship is debated). Average speed is calculated for each link as the average of all the smoothed speeds falling on that link, or:

$$\overline{V}{L{ij}} = \frac{\sum{i} \sum{j} v'_{ij}}{N} \tag{4}$$

where N is the count of those observations on link Lij. The off-peak period was chosen to calculate this indicator to avoid collinearity with CI. In this case, VLij is exactly equal to FFSLij defined above as off-peak speed is the selected covariate (VL )ij . In practice, average speed can be computed for other time periods, so it is not strictly necessary that these measures be equal.

3.2.4. Speed uniformity

Although the magnitude of speed is widely believed to contribute to

crash occurrence, much of the existing literature supports that variation in speed may be a better predictor of risk. The coefficient of variation of speed (CVS) was used in this study. For each link, Lij, CVS is computed as:

$$CVS{Lij} = \frac{\sigma(v^i{ij})}{\overline{V}_{Lij}}$$ (5)

where σ ( ) *v*′ij is the standard deviation of all speeds on link Lij during the considered time period. CVS was found to be most strongly related to crash frequency and severity during the off-peak period (Stipancic et al., 2017a, 2017b).

3.3. Crash frequency model

3.3.1. Latent Gaussian Models

LGM are a subclass of structured additive models, in which the response variable yi for each subject i is assumed to follow a distribution from the exponential family (Normal, Poisson, or Binomial). In this case, yi represents the crash frequency (number of crashes) at link or intersection i and the mean of yi , noted μ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{kl} + \sum_{j=1}^{nf} f^{(j)}(u{ji}) + \varepsilon_i$$ (6)

where *β*0is the intercept, βk are the coefficients representing the linear effect of covariates zki, f ( )j are functions of covariates uji used to relax these linear relationships or introduce random effects, and εi is the unstructured (i.e. containing no spatial or serial structure) error component (Rue et al., 2009). LGM are an extremely flexible family of models because of the forms that f ( )j can take, such as introducing temporal or spatial dependence leading to dynamic or spatial models (Biangiardo and Cameletti, 2015).

From a likelihood perspective, LGM can be represented as a threestage hierarchical structure, beginning with the conditionally independent likelihood function:

$$\pi(yx, \theta) = \prod_{i=1}^{n} \pi(y_i | \eta_i(x), \theta)$$ (7)

where y is the response vector, x is the latent field, θ is the vector of hyperparameters, and η (x) i is the i th additive predictor. Next, the latent Gaussian field is formally defined with a mean m θ( ) and precision matrix Q− (θ) 1 conditioned on hyperparameters θ:

$$x|\theta \sim \mathcal{N}(m(\theta), Q^{-1}(\theta))$$ (8)

Finally, a prior distribution is assigned to the hyperparameters:

$$\theta \sim \pi(\theta)$$ (9)

Although the prior distribution of the latent field must be Gaussian by definition, the prior distributions of the hyperparameters are not subject to this constraint. For more details on LGMs, readers are referred to Rue et al. (2009).

3.3.2. Spatial correlation using the Besag proper model

Previous research has highlighted the need to account for spatial correlations, which can be incorporated to account for similarities of adjacent links or intersections (Biangiardo and Cameletti, 2015) due to the "effect of unknown confounders" (Latouche et al., 2007). Failure to account for spatial dependence may lead to model biases (Latouche et al., 2007; Biangiardo and Cameletti, 2015). 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). Assuming the outcome yi follows a Poisson distribution (note that the same model can be written using a Negative Binomial (NB) distribution), the model is written as:

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

where $$\log(\mu_i) = \beta0 + \sum{k=1}^{p} \betak z{ik} + u_i + v_i$$ (11)

where βk are regression coefficients for fixed effects, zik are covariates,

vi is a site-specific random effect modeled using an exchangeable correlation structure across sites, and ui is a site-specific spatially structured random effect. Several structures can be specified for u = … (, , ) u u 1 n . Here, a proper version of the conditional autoregressive (CAR) structure is chosen:

$$ui \left| u{-i}, \tau, d \sim \mathcal{N} \left( \frac{1}{d+ni} \sum{i \sim j} u_{j}, \frac{1}{\tau(d+n_i)} \right) \right|$$

$$(12)$$

where ni is the number of neighbours for link or intersection i, *u−i* 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).

3.3.3. Integrated nested Laplace approximation

The greatest challenge with Bayesian models is estimating the posterior distributions of the latent field. Traditionally, MCMC simulations have been used to iteratively compute and update the posterior marginal of the latent field. However, this process is time and resource intensive, requiring tens of thousands of iterations (Rue et al., 2009). This shortcoming has led to both new techniques for MCMC simulation (including Hamiltonian Monte Carlo (Radford, 2011)) and new developments in Bayesian inference. The INLA approach was proposed by Rue et al. (2009) to perform Bayesian approximations on LGMs. The INLA approach uses a combination of Laplace approximations and numerical integration to estimate the posterior marginal of the latent field, and offers accurate approximations with a significant reduction in computational time (Schrodle and Held, 2011). The authors demonstrate that this approach provides precise estimates in seconds or minutes where MCMC computations would take hours or days. Furthermore, the authors note that the approximation error in the INLA approach is nearly equal to the estimation error in typical MCMC methods (Rue et al., 2009). A package for the R programming language (R-INLA) has been developed to deploy INLA in a useable and stable format (Schrodle and Held, 2011).

3.3.4. Model calibration and validation

Calibration and validation of the crash prediction model was completed in several steps. Using the INLA approach, Poisson and NB Bayesian models were first estimated at both the link and intersection levels for a subset of the road network. These models use observed crashes as the outcome, y, the extracted SSMs, number of GPS trips (used as a proxy of exposure, or volume), and the roadway length and functional classification as independent variables. Spatial correlation was then incorporated in the best performing Poisson or NB model based on Deviance Information Criteria (DIC) to create a third model for both links and intersections. For the spatial model, this study utilizes Python scripts to generate graphs describing the network topology and defines neighbours as links or intersections which are immediately

Fig. 4. Illustration of a sites (black) and their neighbours (red) for links (a) and intersections (b) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Fig. 5. Smartphone application interfaces.

Fig. 6. Map of study location.

adjacent to one another. For links to be considered neighbours, they must share a common intersection (the links must intersect), while intersections must share a common link (be connected by a single link). This is illustrated further in Fig. 4. The spatial models are then validated on a separate subset of the network with similar geometric and built environment characteristics. 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. The models are compared based on goodness-of-fit variables including mean square error (MSE) and correlation coefficient (CORR).

4. Results

4.1. Data description

The GPS data utilized for this study was collected using the Mon Trajet smartphone application (City of Quebec, 2015) developed by Brisk Synergies (2015) for the City of Quebec, Canada. The application was installed voluntarily by drivers who anonymously logged trips using a simple interface, shown in Fig. 5. This study made use of a sample of data that 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 Quebec (MTQ) for the period between 2000 and 2010. The longest possible time period was chosen to ensure that as many sites as possible had observed crashes. In total, 14,278 minor injury, major injury, and fatal collisions involving at least one vehicle were identified during this time period. 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

A subset of the road network near Laval University, west of the city center as shown in Fig. 6, was used for model calibration and validation. This area was selected for its road density, diversity of roadway functional class, and average trip volumes. An area containing approximately 1000 links and intersections was selected, from which approximately 50% of the sites were selected randomly for calibration, with the remaining sites used for validation. At the link level, the calibration data set contained 71 motorway links, 237 arterial and collector links (grouped from the primary, secondary, and tertiary classifications within OSM), and 220 residential links, for 528 links in total. In terms of intersections, 82 were classified as motorway, 157 as arterial/collector, and 275 as residential, for 514 total intersections, classified by taking the highest classification for links adjacent to the intersection. The validation set contained 58 motorway links and 87 intersections, 288 arterial and collector links and 166 intersections, and 220 residential links and 291 intersections (566 links and 544 intersections total).

Descriptive statistics for all model variables are provided in Table 1. An average link experiences 0.58 crashes over the study period, through the maximum value is 35 crashes. 3.89 crashes occurred at an average intersection, up to 53 crashes at the most extreme site. The number of trips on each facility varies greatly, from a minimum of 2 to a maximum of several thousand. In this study, the number of GPS trips is used as a proxy for exposure (volume). It is assumed that the distribution of GPS trips is proportional to the total demand. Based on a small sample of freeways, arterials, and collectors (45 links) the correlation between average annual daily traffic (AADT) and observed GPS trips was 0.75. This result supports the assumption that GPS trip counts can be used as a proxy of traffic exposure. On average, each trip along a link will experience 0.11 HBEs, while for intersections, the number is 0.23. Congestion tends to be higher along links, while CVS is higher at intersections. Average travel speed is approximately 42 km/h. Correlations for the model variables are provided in Table 2. Note that both the

Table 1 Variables and Descriptive Statistics for the Calibration Data Set.

Mean Minimum Maximum Std. Dev.
Links (sample size = 528)
Crashes 0.58 0.00 35.00 3.15
ln(Trips) 4.38 0.69 7.61 1.22
ln(Length) 4.69 1.50 6.95 0.74
HBEs/Trip 0.11 0.00 0.82 0.13
Congestion Index 0.11 0.00 0.78 0.15
CVS 0.33 0.02 1.40 0.19
Average Speed (km/h) 41.26 7.13 107.42 18.83
Intersections (sample size = 514)
Crashes 3.89 0.00 53.00 7.55
ln(Trips) 5.31 1.79 8.06 1.10
HBEs/Trip 0.23 0.00 4.66 0.38
Congestion Index 0.12 0.00 0.77 0.14
CVS 0.32 0.02 0.97 0.15
Average Speed (km/h) 42.91 8.32 107.42 18.9

Table 2 Correlations between Model Variables at the Link and Intersection Level.

Trips HBEs /
Trip		 Congestion
Index		 CVS Average
Speed
Links
ln(Trips) 1.00 0.36 0.29 −0.13 0.56
HBEs / Trip – 1.00 0.20 0.24 −0.10
Congestion – – 1.00 −0.19 0.30
Index
CVS – – – 1.00 −0.62
Average Speed – – – – 1.00
Intersections
ln(Trips) 1.00 −0.29 0.23 −0.12 0.50
HBEs / Trip – 1.00 −0.08 0.07 −0.14
Congestion – – 1.00 −0.20 0.29
Index
CVS – – – 1.00 −0.64
Average Speed – – – – 1.00

number of trips and link length were converted using the natural log. For intersections, multiple functional forms relating exposure to crash frequency, such as those proposed by Miaou and Lord (2003), were tested. However, the natural log of the sum of all trips passing through the intersection was found to provide the best results.

4.3. Model calibration

At both the link and intersections levels, three separate models were estimated, resulting in six models total. First, the classic Poisson and NB Bayesian models are estimated using INLA and compared. Next, the Besag Spatial model is fitted using the best performing distribution (either the Poisson or NB) and compared by DIC to determine superiority. DIC was selected because it is less sensitive to extreme values. Goodness of fit measures, including MSE and CORR, are potentially sensitive to extreme values estimated by R-INLA for sites without neighbours.

4.3.1. INLA Poisson and negative binomial models

Link level results for the Poisson and NB models are presented in Table 3, where variables statistically significant at 95% confidence are bolded in red. The NB model outperformed the Poisson by DIC, though it had a poorer goodness-of-fit by mean-square-error (MSE). Goodnessof-fit is illustrated further in Fig. 7a and c, where fitted values are plotted against the number of observed crashes. Both models perform well, resulting in MSE below 9.0 and CORR above 0.40.

Considering the posterior means of the covariates, results generally supported the relationships between SSMs and crash frequency established in previous studies (Stipancic et al., 2018, 2017a, 2017b). In both models, the posterior mean for the number of trips is positive. This result mirrors expectations that increasing exposure (increasing traffic volumes) should increase the likelihood of a crash. Therefore, it is logical that sites with higher traffic volumes also have a higher number of crashes. One major departure from expectation is the negative posterior mean associated with HBEs observed in both models. In previous studies (Fazeen et al., 2012; Jun et al., 2007), HBEs were found to be positively correlated with crash frequency, while the posterior mean for HBE/Trip in both the Poisson and NB models is negative. Initially, it is unclear why the effect is negative in these most basic models.

For the traffic flow SSMs, the posterior means for both CI and CVS are positive (increasing congestion and speed variation are generally associated with increasing crash frequency). Again, it is logical and expected that sites with higher levels of congested and higher variation in travel speeds would also have a higher frequency of crashes. Not only does this result mirror intuition, but it supports results from previous study (Stipancic et al., 2017a, 2017b). For average speed, the mean is also positive (all else being equal, facilities with higher average speeds

Table 3 Link Model Results for Poisson and Negative Binomial Models (For interpretation of the references to colour in this table legend, the reader is referred to the web version of this article).

Note: Variables significant at 95% confidence are bolded in red.

tend to have more crashes). Although this result contradicts earlier research (Stipancic et al., 2017a, 2017b), it is logical that crash frequency would increase with speed. Longer links experience more crashes, and when it comes to functional classification, all else being equal, motorways experience fewer crashes compared to residential streets (negative posterior mean), while and arterials/collectors experience more (positive posterior mean). In both models, most of the proposed covariates were statistically significant at 95% confidence.

Table 4 contains the results for the intersection level Poisson and NB models. Model fit is illustrated in Fig. 7b and d. Although the CORR values are comparable to the link-level models, the intersection level models have a higher MSE (approximately 40 compared to less than 10 for link-level models). However, this is expected simply because crash counts at intersections are much higher on average. As with the link level model, the NB was the superior by DIC though it had a slightly higher MSE. Although this result may seem counterintuitive, considering results for sites with zero observed crashes (along the y-axis), the NB model estimates higher crash counts in both the link- and intersection-level models compared to Poisson, effectively increasing MSE. Furthermore, because DIC is penalized for number of parameters, and MSE is not, it is possible that the measures can give conflicting results.

The posterior mean for trips was, again, positive and statistically significant. However, at the intersection level, the posterior mean for HBEs was also positive, in accordance with previous work (Stipancic et al., 2018). It is expected that intersections with more braking events would also have more crashes (as HBEs are considered evasive manoeuvres), a result that is observed in the intersection-level modelling results. For CI and CVS, results were also similar, with both posterior means positive. Unlike the link level model, the posterior mean of average speed is negative, which matches expectations from earlier research (Stipancic et al., 2017a, 2017b). This result may seem counterintuitive, with Imprialou et al. (2016) attributing this negative relationship to link-aggregation bias. However, in this context, it is suggested that FFS captures other variables not considered, including geometry. In these models, the posterior mean for the motorway variable is negative (motorway intersections, including on- and off-ramps, are less likely than residential intersections to experience crashes), while for arterials and collectors, it is positive. Again, most of the proposed covariates were statistically significant at 95% confidence.

4.3.2. Spatial models

The NB model was the best performing model by DIC for both links

and intersections. Therefore, the spatial component was incorporated into the NB model formulation to yield spatial NB models, which are summarized in Table 5. By DIC, the spatial NB models were the best performing of all considered model types, and improved the goodnessof-fit compared to the non-spatial NB models, as shown Fig. 7e and f. The spatial model greatly improved the goodness-of-fit at the link level, while having virtually no effect for intersections. The reason for this may be that adjacent links truly are adjacent (are physically connected at intersections) and have much stronger spatial correlations than intersections, where neighbouring sites may be separated by several hundred metres or more. Not only were these models statistically superior while improving goodness-of-fit on the calibration data set, supporting previous work (Jiang et al., 2014) indicating the importance of including spatial correlations in crash frequency models, but the signs of the posterior means for covariates generally matched expectations from previous studies. Means for trips were consistently positive. Although HBEs/Trip had a negative posterior mean in the linklevel model, the mean was positive at the intersection level. CI and CVS were similarly positive, while average speed had a positive mean for links and a negative mean for intersections. However, these models also contain more non-significant variables than the non-spatial models (at 95% confidence). For example, trips, HBEs, congestion, and speed were non-significant in the link-level model, while HBEs, congestion, speed variation, and average speed are non-significant in the intersectionlevel model. It is important to note that these models are estimated on a relatively small sample of the whole road network, and expanding the model to the entire network will likely reveal the true impact and significance of these variables.

4.4. Model validation

The NB Spatial models developed above were validated on a separate data set randomly selected from the data sample near Laval University. The results of this validation procedure are summarized in Fig. 8a for the link-level model and Fig. 8b at the intersection level. As can been seen in Fig. 8b, the model makes accurate predictions for the intersection-level data set, and a CORR of 0.56 is quite promising. It is believed that the sample size used in this study was not sufficient to yield accurate link-level predictions. Still, prediction at the link level produces an MSE below 4.0. At the time of publication, R-INLA has no specific module for model validation. Prediction must be accomplished by setting null values for the validation sites during model calibration. This limitation presents challenges when validating spatial models, as

Fig. 7. Fitted values versus observed crashes for Poisson models for links (a) and intersections (b), NB models for links (c) and intersections (d) and NB Spatial models for links (e) and intersections (f).

some extreme values were inevitable for specific sites without neighbours for this specific model. Therefore, predictions for sites without neighbours are based on the fixed effects of the model only. Although both models are capable of describing the relationships between SSMs and crash frequency, the intersection-level model shows more promise

of predicting crash occurrence, though increasing the sample size used for model calibration is expected to further improve the accuracy of predicted values.

Table 4 Intersection Model Results for Poisson and Negative Binomial Models (For interpretation of the references to colour in this table legend, the reader is referred to the web version of this article).

Note: Variables significant at 95% confidence are bolded in red.

Table 5 Negative Binomial Spatial Model Results for Links and Intersections (For interpretation of the references to colour in this table legend, the reader is referred to the web version of this article).

Note: Variables significant at 95% confidence are bolded in red.

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

Fig. 9. Sites ranked by posterior expected number of crashes for links (a) and intersections (b).

4.5. Network screening

The following section describes how such a model could be practically applied to screen the road network for dangerous sites. For this work, the posterior expected number of crashes is used as a ranking criterion. Obviously, this approach is simplistic, but illustrates the most crash prone sites throughout the network as predicted by the model. The posterior expected number of crashes is computed using the mean of the covariates for each site combined with the spatial effect contributed by the BYM model. In the specific case of sites without neighbours, the number of crashes is predicted using only the fixed effects of the model (the mean of the covariates, with no spatial component). This addresses the issue of extreme crash counts predicted by INLA for solitary sites. The outcomes of the posterior analysis are illustrated in Fig. 9. These maps clearly show the locations of dangerous sites (sites with a high expected number of crashes) as predicted by the model, shown in red. For links, the most extreme sites are expected to experience 15 crashes in an 11 year period, while for intersections, 25 crashes is the expected maximum. Furthermore, the most dangerous sites tend to be along the major arterials within the study area or on freeway ramps and interchanges, as would be expected for a major city. Once such a model is calibrated for a given network, the results (the posterior number of crashes) could be used to estimate crashes based on surrogate indicators (when crash data is not available) or to supplement historical crash data in the network screening process (when crash data is available).

5. Conclusions

The purpose of this paper was to propose a surrogate safety screening approach based on smartphone GPS data and an FB modelling framework. SSMs were derived from GPS data and used as predictive variables along with roadway functional classes. Then, latent Gaussian Spatial Models for crash frequency were estimated using the INLA technique. The used data set was a large sample of GPS data collected in Quebec City. Extracted SSMs were related to vehicle manoeuvres (hard braking) and measures of traffic flow (congestion index, average speed, and speed variation). A subset of the road network near Laval University was used to calibrate the statistical models, which were estimated on 11 years of crash data. INLA NB models outperformed Poisson models at both the link and intersection level. In general, the relationships between SSMs and crash frequency established in previous studies were supported by the modelling results. The results indicate that crashes are more likely at links and intersections with more trips (as a proxy of exposure). Positive means for CI, CVS, and link length also indicate a positive relationship with crash frequency. The only disagreement between link- and intersection-level models was for HBEs and average speed. At the link-level, HBEs had a negative mean (positive at the intersection-level) and speed had a positive mean (negative at the intersection-level).

The greatest improvement in model fit was achieved through a spatial model. By DIC, the NB Spatial models were the best performing of all considered model types and improved the goodness-of-fit at the link level. This result supports previous work indicating the importance of including spatial correlations in crash frequency models. Some of the SSMs were observed to be non-significant in the spatial models. Despite this, the sign of the posterior means of the significant variables generally supported the results of previous research (the relationships between SSMs and crash frequency were confirmed). Although further work and larger data sets are required to improve model fit and prediction power, these results are promising, and support the implementation of a Spatial LGM based on SSMs and GPS data in network screening. When considering the predictive power of the model based on a separate validation data set, the intersection-level model performed well, providing a relatively high correlation between observed and predicted crashes. Although the link-level model struggled with prediction, it still provided a low MSE, and it is believed that prediction will be improved using a larger smartphone GPS data set, consisting of 3 weeks of continuous data covering the entire Quebec City road network.

One shortcoming of this work is that all traffic flow measures are calculated at the link level and must be aggregated for intersection-level analysis. Despite this, the fit of the intersection models are comparable to their link-level counterparts, even if more traffic flow measures are insignificant. Additionally, link-level data is not direction dependent, as direction of travel was not supported by the crash dataset. Incorporating direction of travel is theoretically possible and relatively simple, if it is supported by additional detail in the crash dataset. The small sample size likely contributes to poor predictive power and insignificance of variables. Clearly, more work is needed to improve predictions at both the link and intersection levels. Future work will aim to improve these models in several ways. First, the crash models will be expanded to the entire Quebec City road network, which will reveal the true effect and significance of the proposed covariates. Second, other model formulations will be tested in the future, to improve model flexibility with respect to the selected covariates and explore other spatial correlation structures. More complex ranking criteria, such the posterior probability of excess or the posterior of ranks (Miranda-Moreno et al., 2007; Heydecker and Wu, 2001), could also be considered. In addition to improving the crash frequency models presented herein, crash severity should also be considered, as collision frequency and severity are two independent dimensions of risk. Future work should also consider the impact of the crash assignment method. Geo-locating crashes on the road network based only on text from police reports is a well-known problem. Although several assignment algorithms were tested, more work is needed to determine the impact of assignment including model sensitivity to crash mis-location. Finally, other smartphone sensor data, including accelerometer, magnetometer, and gyroscope data, could be used to improve the accuracy of the SSMs. The ability to screen the network based on SSMs represents a substantial contribution to the field of road safety, and importantly works towards the elimination of crash data as a necessity in safety evaluation and monitoring.

Acknowledgement

Funding for this project was provided in part by the Natural Sciences and Engineering Research Council of Canada.

References