category: literaturenote citekey: rezaeeprocessconvolutionmodelcrash2022 title: A process convolution model for crash count data on a network authors: "Rezaee, Hassan; Schmidt, Alexandra M.; Stipancic, Joshua; Labbe, Aurélie" year: 2022 date: 2022-00-00 2022 doi: 10.1016/j.aap.2022.106823 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2022.106823" zotero_key: UFNPF8Q6 zotero_storage: 5XKGBIIK collections: transpordi_planeerimine folder: 001_artiklid/Road Safety Indicators firstAuthor: "Rezaee, Hassan"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

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

A process convolution model for crash count data on a network

Hassan Rezaee a , Alexandra M. Schmidt b , Joshua Stipancic c , Aurélie Labbe a,∗

^a Department of Decision Sciences, HEC Montréal, Montréal, QC, Canada
^b Department of Epidemiology, Biostatistics and Occupational Health, McGill University, Montréal, QC, Canada
^c Data Lab, Intact Insurance, Montréal, QC, Canada

A R T I C L E I N F O

Keywords: Crash modelling Road network Spatial correlation Process convolution Path distance Autoregressive model

A B S T R A C T

Crash data observed on a road network often exhibit spatial correlation due to unobserved effects with inherent spatial correlation following the structure of the road network. It is important to model this spatial correlation while accounting for the road network structure. In this study, we introduce the network process convolution (NPC) model. In this model, the spatial correlation among crash data is captured by a Gaussian Process (GP) approximated through a kernel convolution approach. The GP's covariance function is based on path distance computed between a limited set of knots and crash data points on the road network. The proposed model offers a straightforward approach for predicting crash frequency at unobserved locations where covariates are available, and for interpolating the GP values anywhere on the network. Inference procedure is performed following the Bayesian paradigm and is implemented in R-INLA, which offers an estimation procedure that is very efficient compared to Markov Chain Monte Carlo sampling algorithms. We fitted our model to synthetic data and to crash data from Ottawa, Canada. We compared the proposed approach with a proper Conditional Autoregressive (pCAR) model, and with Poisson Regression (PR) and Negative Binomial (NB) models without latent effects. The results of the study indicated that although the pCAR model has comparable fitting performance, the NPC model outperforms pCAR when the main goal is to predict unobserved locations of interest. The proposed model also offers lower mean absolute error rates for cross validated crash counts, latent variable values, fixed-effect coefficients, as well as shorter interval scores for singletons. The NPC provides a natural way to account for the road network structure when considering the inclusion of spatially structured latent random effects in the modelling of crash data. It also offers an improved predictive capability for crash data on a road network.

1. Introduction

Despite worldwide efforts to improve road safety, road crashes are still among the leading causes of death and injuries. The modest decrease in the death rate per 100,000 population, from 18.8 in 2000 to 18.2 in 2016 (World Health Organization, 2018), indicates that road accident prevention has not yet reached its full potential. With the advent of new GPS and on-board diagnostics technologies for collecting telematics data, large spans of road networks can now be covered with features extracted from thousands of daily recorded GPS trace data. These features carry important information on congestion, average speeds or speed variations on a network level, and also about driving behaviours at the per-driver level. An important factor in modelling road network crashes is being able to exploit these large and inexpensive data sets and quantify their association with large quantities of crash data collected from the road network. Hence, a proactive, datadriven approach should be adopted to predict and potentially prevent road crashes by identifying major contributing factors and potential hotspots in time and space.

The main goal of road crash data analysis is to model the crashgenerating process so that one can (i) infer the association between crashes and environmental or driving behaviour variables, and (ii) predict the expected number of crashes at a given location on the road network. In the past two decades, several crash modelling approaches have been developed and applied to crash data, providing many new insights about the role of contextual variables in crash frequency or severity level, and potentially predicting crashes at target locations. Various techniques have also been proposed for the spatial analysis and modelling of road crash data; see Ziakopoulos and Yannis (2020) for a review. Among these tools, one can distinguish three main approaches to model spatial correlation: (i) neighbourhood-based approaches, (ii) machine learning approaches and (iii) distance-based approaches. The first approach relies heavily on an arbitrary and discretized definition of

E-mail address: aurelie.labbe@hec.ca (A. Labbe).

^∗ Corresponding author.

the neighbouring structure of the road network, and prediction at target unobserved locations is based mostly on simple interpolations. The second approach relies on an adequate tuning of the required hyperparameters during the training phase. Furthermore, deep learning models are developed for prediction purposes and cannot perform inference about the association between the covariates and road crashes. Unlike Bayesian methods, the deep learning methods do not offer a straightforward approach to assessing the uncertainty around the estimated parameters. Quantifying uncertainties around the estimated values of crash predictors and predictions for unobserved locations of interest is of paramount importance as it can directly affect decision- and policymaking in long-term road safety management and network screening. Finally, the third approach, used in this paper, assumes a distancebased spatial correlation structure. This approach is more appealing, especially if the distance used accounts for the road network structure. Such distances being non-Euclidean, they can result in models that do not guarantee the positive definiteness of the covariance function, a problem mostly ignored in the crash modelling literature.

Despite these steady methodological developments, several fundamental issues related to the way spatial correlation is accounted for in the model remain to be addressed. Crash data on adjacent spatial units often show spatial correlation because they share similar crash risk factors. This spatial correlation can be caused by observed effects such as traffic flow and road geometry or by unobserved effects such as environmental characteristics throughout the road network. For the former, a clear example is illustrated by Castro et al. (2013) for the case of crash frequencies around an intersection where the presence of an island can not only increase the probability of a crash occurring at that location, but can also ''have a spatial spillover effect on the injury propensity at proximally-located crash sites''. However, the most common source of spatial correlation in crash data is due to unobserved factors with inherent spatial correlations such as pavement irregularities, microclimate effects, and local sight-distance restrictions (Mannering and Bhat, 2014). This unobserved heterogeneity can complicate the modelling process, and ignoring spatial correlation from unobserved factors can result in models that lack sufficient predictive power and lead to a mis-specification of important contributing risk factors (Mannering and Bhat, 2014; Chiou et al., 2014).

The goal of this paper is to present an efficient and valid statistical crash model, along with an available implementation, which accounts for the road network structure. Spatial correlation is modelled to perform both inference and prediction tasks within a Bayesian framework, by taking advantage of the latest developments in the literature on spatial analysis. Building on the ideas proposed by Ver Hoef et al. (2006) and Okabe et al. (2009), an important part of this work is to account for the road network structure in the spatial random effect component of a generalized linear mixed model. We define the spatial random effect as a Gaussian Process (GP) approximated via the convolution of a moving average function (kernel) and an unstructured latent random process (Higdon, 2002). This work also builds partly on the kernel density estimation technique for road networks proposed by Okabe et al. (2009). This allows a class of flexible GPs whose covariance functions account for the network structure.

In the remainder of this article, we first provide detailed information of the existing studies in Section 2. In Section 3, we introduce the methodology where the distance metric, the statistical model including the likelihood function, and fixed and random effects are described in detail. A simulation study comparing our proposed approach with a proper CAR model and a Poisson and Negative Binomial regression is presented in Section 4, which also includes a case study using crash data from the city of Ottawa, Canada, along with an implementation of the model using the Integrated Nested Laplace Approximation (INLA) technique. Finally, Section 5 provides some concluding remarks and potential avenues for future research.

2. Literature review

Neighbourhood-based approaches model the spatial correlation based on a predefined discretized neighbouring structure. In such models, defining the neighbourhood is an essential step in (i) measuring the strength of the spatial relationships between observed data and (ii) quantifying the influence of neighbours on each other. One of the techniques used in crash data spatial modelling is the Geographically Weighted Regression (GWR) method (Fotheringham et al., 2003; Liu et al., 2017; Gomes et al., 2017; Hadayeghi et al., 2003), which can be considered an extension of the traditional regression model fitted over local neighbourhoods. In this model, the neighbourhood structure exerts a large influence on the smoothing properties of the model for which no clear understanding is available. Bayesian models using an autoregressive prior (Besag et al., 1991) form another popular class of methods in road crash spatial modelling. The two main variants of the autoregressive models are the SAR (simultaneous) and CAR (conditional) models. In the CAR model, the spatial correlation is described via a neighbourhood matrix, which is straightforward to define but not flexible enough to adequately model complex road networks with semi-linear structure. Note that autoregressive models were originally developed for disease mapping problems and hence are more suitable for areal data. The application of such models to road crash data depends on the definition of the neighbourhoods on the road network, which is subjective and depends on the practitioner to a large extent. This model has also been applied to crash data on a road network at both link and intersection levels by Stipancic et al. (2019), where both crash frequency and crash severity levels were modelled using Poisson and Negative Binomial likelihood functions. Aguero-Valverde (2014) used the CAR model in comparison with an exponential covariance model to analyse crash data on a link level in Centre County, Pennsylvania, USA. The effective spatial correlation range obtained from the exponential spatial correlation model indicated that a first-degree neighbourhood structure is sufficient to capture the spatial correlation existing among crash data. The definition of the neighbourhood structure in the CAR model is again shown to be dependent on the study area and the particular data set under investigation. Although neighbourhood-based approaches address crash data observed on both intersections and road links, segment-level analysis is a more common approach. For example, Aguero-Valverde and Jovanis (2008) applied the CAR framework to explore the effect of spatial correlation in models of road crash frequency at the segment level. By testing several segment neighbouring structures and their effect on model results, this study concluded that the simplest neighbourhood structure based on the first-order adjacency matrix performed best. Satria et al. (2021) also used the CAR model to evaluate the spatial correlation between observed severity of crash data at the segment level in Indonesia, highlighting the importance of incorporating spatial random effect in the model. Finally, two spatial modelling techniques including CAR and the multiple membership (MM) model were compared with the nonspatial Poisson–lognormal model in El-Basyouny and Sayed (2009). Both spatial models had significant estimates for the spatial correlation parameters where a considerable portion of the total variability was shown to come from the spatial random effect.

Machine learning and deep learning approaches have also been used for road crash spatial modelling; see e.g. Bao et al. (2019) and Effati et al. (2015). Such data-driven methods are interesting choices considering the massive quantities of contextual variables that are collected easily using a variety of tools such as on-board cameras, GPS dongles connected to the vehicle's on-board computers, or smartphone applications. The main drawback is that these models are mostly predictive, and hence one cannot make inferences about the significance of each variable on crashes as in a statistical model. Also, no standard guideline is currently available on how to tune hyperparameters during the training phase of these models. Cai et al. (2019) provides an application of the Convolutional Neural Network (CNN) using high-resolution data. In their study, the input layer of the CNN was constructed by splitting the grid containing the road network into small cells and adding different layers of data such as land use, roadway, and socio-demographic data. The proposed model was reported to perform better than ''traditional'' Negative Binomial and spatial Poisson lognormal models. There have been other attempts at analysing the spatial and temporal correlations of crash frequency data using deep learning techniques. In Ren et al. (2018), a deep recurrent neural network approach was used to predict traffic accidents in Beijing, considering both the spatial and temporal correlations. Although their results proved to be effective and can be applied to identify hazardous locations, it did not take into account the real road network structure. Moreover, these methods do not offer a straightforward approach for quantifying uncertainty around their predictions. Also, as pointed out in Ziakopoulos and Yannis (2020), these models require custom data sets that are collected specifically for the particular study at hand in order to exploit their full potential. Interested readers can refer to Abduljabbar et al. (2019) for a complete review on the application of deep learning in transportation research.

Capturing the spatial correlation over a road network using a path distance metric that takes into account the structure of the network is intuitively appealing. In this case, the spatial correlation is determined based on the path cost rather than a simple proximity in the Euclidean space. The difference between the two becomes more striking if we consider an extreme case with two crashes, one of which is located on a bridge and the other on the road passing beneath it. Computing the path distance over a road network is straightforward, and efforts have been made to introduce non-Euclidean distances to model data observed over a network, see e.g., Ver Hoef et al. (2018), Curriero (2006), Ver Hoef and Peterson (2010), Cressie and Majure (1997), Cressie et al. (2006). However, the use of such distance in models initially developed in a Euclidean setting can result in models that do not guarantee the positive definiteness of the resultant covariance matrix. For example, Gardner et al. (2003) used a Spherical covariance model developed for Euclidean distance to perform Kriging based on the stream distance. In the context of networks of trails, Ladle et al. (2016) used a network distance to perform ordinary Kriging estimations. In road safety, Shiode and Shiode (2011) and Selby and Kockelman (2013) used the network distance between sites on a road network instead of the Euclidean distance. However, these distances were directly used in Kriging to perform prediction at unsampled locations, which can result in invalid covariance functions. Euclidean distance was also replaced by shortest-path distance in Little et al. (1997), Rathbun (1998), and Jensen et al. (2006) which can also lead to invalid models. This issue has been discussed in detail by Curriero (2006) and Ver Hoef (2018). In general, using non-Euclidean distances along networks, e.g. roads, trails and streams, requires the use of one of the following strategies: (i) retain the Euclidean distance-based covariance model but modify the distance matrices such that they guarantee a positive definite covariance matrix; or (ii) develop new covariance functions that are compatible with both non-Euclidean and Euclidean distances. For the former, an idea proposed by Curriero (2006) is to use multidimensional scaling to map the non-Euclidean distance matrix onto a lower dimensional Euclidean space, while minimizing the enforced distortion. This method is associated with two potential problems. First, as noted in Ver Hoef et al. (2006), the addition of a new prediction point can alter the predictions at all other points. Second, the computational complexity of the multidimensional scaling for large data sets can be prohibitive. For the latter, Ver Hoef et al. (2006), Ver Hoef and Peterson (2010) have proposed covariance models based on moving average functions that can model the observed data and potentially predict at unsampled locations based purely on stream distances. This is the approach taken in this paper.

Fig. 1. Sample road network.

3. Method

3.1. General framework

In this paper, we consider crash data and covariates aggregated either at intersections or at road segments called links. The road network can be seen as a network = ( , ) where are the vertices, i.e. intersections, and L are the edges, i.e. road links. In the remainder of this paper, we use the terms links and intersections to denote the edges and vertices of the network, respectively. Using the standard network terminology, the degree of an intersection, (.), represents the number of connected links. We consider data observed at links or intersections over the network. The corresponding spatial locations are denoted ¹ , …, . If the observed data is measured on a link, corresponds (without loss of generality) to the middle point of the link. If the observed point is measured at an intersection, is the location of the intersection. The number of observed crashes at any location on is denoted (). This framework is illustrated in Fig. 1, by displaying a sample road network with 15 vertices (shown in black circles) and 10 crashes observed at locations (¹ , …, 10) among which eight are on links (represented by grey squares on the link midpoints) and two are at intersections (² and ⁹ ). The corresponding link is shown with an arrow in Fig. 1. Note that crash data are not necessarily available for all links and intersections. If a connecting path exists between spatial locations and ′ on , then we denote by *,′ the set of starting and ending locations, and , and all the intervening intersections (,*′) on the connecting path between and ′ on , i.e., *,*′ = { ∪ *,*′ ∪ ′}. For the spatial sites illustrated in Fig. 1, this would translate as ¹ *,*8 = {¹ , ⁵ , ⁸ , 12, ⁸ }, ¹ *,*9 = {¹ , ⁵ , ⁸ , ⁹ ² *,*9 = {² , ⁸ , ⁹ } and ¹ *,*2 = {¹ , ² }. All notations and symbols used in the paper are summarized in Table 1.

Let () be the number of crashes at location ∈ . We assume that conditional on (), () follows a Poisson distribution, such that

$$Y(s) \mid \lambda(s) \sim Poisson(\lambda(s)),$$ (1)

where () is the Poisson mean at spatial location and log( ()) is decomposed as the sum of fixed and latent (random) effects, that is,

$$\log \lambda(s) = \log l(s) + \beta_0 + X(s)\beta + Z(s), \tag{2}$$

where () is an offset given by the length of the link, ⁰ is an intercept, () = (¹ (), …, ()) is a -dimensional vector of known covariates at each spatial location ∈ , and ∈ R is a −dimensional vector of the coefficients associated with each (⋅). The latent component

Table 1 Notations and symbols used in the paper.

Notation	Description
𝑁 = (𝑉 , 𝐿)	Road network 𝑁 with vertices 𝑉 (intersections) and edges 𝐿 (road links)
𝑠	Any spatial locations on the road network 𝑁
𝑌 (𝑠𝑖 )	Number of crashes observed at location 𝑠𝑖 (𝑖 = 1,… , 𝑛)
𝑣𝑠,𝑠′	1 , 𝑣2 ′ Set of vertices (𝑣 , …) on the connecting path between locations 𝑠 and 𝑠 𝑠,𝑠′ 𝑠,𝑠′
𝛺𝑠,𝑠′	′ Set of starting and ending locations, 𝑠 and 𝑠 , and all the 𝑛𝑣 intervening intersections (𝑣𝑠,𝑠′)
	′ on 𝑁 i.e., 𝛺𝑠,𝑠′ = {𝑠 ∪ 𝑣𝑠,𝑠′ ∪ 𝑠 ′} on the connecting path between 𝑠 and 𝑠
𝑒 𝑒 𝑑 = 𝑑 (𝑠, 𝑢) 𝑠,𝑢	Euclidean distance between any two locations 𝑠 and 𝑢 on N
𝑠𝑝(𝑠, 𝑢) 𝑑	Shortest path distance between any two locations 𝑠 and 𝑢 on N

(.) is present to capture whatever remains left in the log-mean of the outcome, after adjusting for the covariates. We expect that this structure has some spatial structure, therefore, we assume that (.) follows a priori a zero mean GP with some valid covariance function. In classical spatial modelling, the covariance function defining a GP is constructed based on the Euclidean distance. See Banerjee et al. (2014) for examples of valid covariance functions based on the Euclidean distance among locations. Our main contribution lies in proposing a spatial random effect that accounts for the road network structure. We follow the Bayesian paradigm to make inferences about the parameters of the proposed model. The next subsection discusses in detail the specification of the distribution of (.).

3.1.1. Prior specification of the latent effect (.)

A major problem associated with a Gaussian process is known as the ''big n problem'', which occurs in the presence of a large number of observed data points. This is caused by the extremely expensive computation time required to invert the covariance matrix and to compute its determinant, leading to a computational complexity of ( ). There have been different proposals in the literature to handle high-dimensional spatial data while still keeping the properties of the underlying spatial process. See Heaton et al. (2018) for a review.

We propose an approximation to (.) based on the convolution approach originally proposed by Thiébaux and Pedder (1987). This technique approximates a stationary GP (.) as a function of a Brownian motion (.) and a moving average/kernel function (.) (Higdon et al., 1999). Due to the semi-linear structure of a road network, both (.) and (.) are defined here on the real line R; however, their extension to R² and R³ is straightforward (Higdon, 2002). A GP () for any location ∈ can be written as

$$Z(s) = \int{-\infty}^{\infty} k(d{s,u}^e) \psi(u) du, \qquad (3)$$

where () is an un-structured latent random process following a normal distribution (0, ² ), defined independently on spatial locations , ( , ) is a kernel function defined in R where, commonly, , = (, ) is the Euclidean distance between any two locations and on . This integral cannot be numerically solved for the majority of the kernel functions. An alternate approach is to approximate (3) through a discretization (Higdon, 2002). This is done by specifying a regular grid with, say, points (knots) over the region of interest, and assuming

$$Z(s) = \int{-\infty}^{\infty} k(d{s,u}^{e}) \psi(u) du \approx \sum{j=1}^{m} k(d{s,u{j}}^{e}) \psi(u{j}), \tag{4}$$

where (¹ )*,*… , () are latent process values at knot locations ¹ *,*… , . Using this lower dimensional representation, a computational gain can be achieved when ≪ . One great advantage of using the process convolution approach when modelling (.) is that it is non-parametric and hence does not require computing positive definite covariance matrices and estimating the covariance function parameters (Higdon, 2002). As a result and because the elements of are independent, the covariance function resulting from this construction is defined solely by the kernel (.), and is expressed as

$$c{s,s'} = c(d{s,s'}^e) = Cov(Z(s), Z(s')),$$

$$= \sigma{\psi}^2 \int{-\infty}^{\infty} k(d{s,u}^e)k(d{s',u}^e)du,$$

$$\approx \sigma{\psi}^2 \sum{j=1}^m k(d_{s,uj}^e)k(d{s',u_j}^e),$$ (5)

where again *,′ = (,* ′ ) is the Euclidean distance between spatial locations and ′ on . Next we describe how to accommodate the network structure on the description of the spatial random variable (.).

The Euclidean distance is not an appropriate metric to measure distance between events on a road network. Instead, we propose using the shortest path distance. The concept is similar to the stream distance on river networks as proposed by Ver Hoef and Peterson (2010). The shortest path on a road network for traffic-related data can be computed only if network links are accompanied with flow direction (traffic flow in our case). In the absence of flow direction data, one can simply consider an undirected graph. Considering the number of intersections on the shortest path trajectory between two locations and ′ as , the shortest path distance (, ′ ) with > 0 on a road network is computed as

$$d^{sp}(s,s') = d^{e}(s,v{s,s'}^{1}) + \sum{i=1}^{n{v}-1} d^{e}(v{s,s'}^{i},v{s,s'}^{i+1}) + d^{e}(v{s,s'}^{n_{v}},s'), \tag{6}$$

which is simply the sum of the Euclidean distances between all the spatial locations along the shortest path in *,′ . The shortest path distance between locations and ′ with = 0 is simply (,* ′ ) where *,*′ = ∞. For spatial sites illustrated in Fig. 1, (¹ , ⁸ ) = (1 , ⁵ ) + (5 , ⁸ ) + (8 , 12) + (12, ⁸ ) and (¹ , ⁴ ) = (1 , ⁵ ) + (5 , ⁴ ). The Dijkstra's algorithm (Dijkstra, 1959) implemented in the R package igraph (Csardi et al., 2006) is used to find the shortest path trajectory and the distance between any two locations on . This algorithm offers reasonable computation times for road networks defined on a metropolitan scale.

3.1.2. Kernel reweighting procedure to account for the network structure and ensure a valid covariance function

In the literature, there are different covariance functions based on the Euclidean distance (Banerjee et al., 2014). When an alternative metric is used, these functions do not necessarily lead to a valid covariance function, i.e. they are not necessarily positive definite. See Ver Hoef et al. (2006) or Curriero (2006) for a review. Therefore, care must be taken when choosing functions to describe the covariance function of a GP based on a metric other than the Euclidean function. In this study, we adapt our model's kernel function to ensure a valid spatial covariance function, based on the kernel convolution method proposed for networks by Ver Hoef et al. (2006). First, we note that the construction of (.) over the network using the kernel approach presented in Eq. (4) is defined on R. However, a road network is not a one-dimensional space as the links get split into branches. The kernel function used in the covariance function in Eq. (4) must be reweighted to account for the network structure (Ver Hoef et al., 2006; Ver Hoef and Peterson, 2010; Okabe et al., 2009). Failing to do so will make the

Fig. 2. [Left] Kernel function values without weighting. [Right] The same kernel function re-weighted according to Eq. (7).

random variable (.) nonstationary at the intersections, meaning that the total area under the kernel curve will be greater for points near an intersection than for locations located further away. This implies that the variance of (.) will be artificially inflated for those points.

Following (Okabe et al., 2009) and Ver Hoef et al. (2006), the weight given to ( ,) for a pair of data and knot locations (, ), respectively, and with > 0, is a simple function of the number of branches between and as

$$w{s,u} = \begin{cases} \frac{1}{\prod{i=1}^{nv} deg(v{s,s'}^i) - 1}, & \text{if } s \notin V, \ \frac{2}{deg(s)} \frac{1}{\prod_{i=1}^{nv} deg(v{s,s'}^i) - 1}, & \text{if } s \in V, \end{cases}$$ (7)

where ( *,*′ ) is the degree of the th intersection in *,*′ along the shortest path identified as , as defined in Section 3.1. For pairs of data and knot locations (, ) with = 0, , = 1. Eq. (7) states that the kernel function value must be split starting at the kernel centre if is located exactly on a vertex of degree (). The two components of Eq. (7) are consistent given that each spatial site located on a link can be seen as an intersection of degree 2, i.e., () = 2, ∀ ∉ .

The kernel weighting process is illustrated in Fig. 2. In this small example, a simple road network is considered, with two intersections ¹ and ² with degrees 4 and 3, respectively, one data point ¹ , and six knot locations ¹ , .. . ⁶ . We are interested in evaluating ( 1 , ), for = 1, …, 6. On the left hand side, the original unweighted kernel function values centred at ¹ are shown, similar to the one-dimensional case along the main line connecting ¹ and the locations. On the righthand side, the same road network configuration is shown with kernel function values weighted according to the structure of the road network at intersections ¹ and ² . The kernel function is naturally highest at the centre ¹ and decays in either direction with increasing shortest path distances, with splits at intersections where the kernel function values are reweighted according to Eq. (7). The weighted kernel values, ∗ (.), evaluated at pairs (¹ , ), for = 1*,*… , 6 are:

$$\begin{split} k^(d^{sp}_{3_1,u1}) &= \frac{1}{\sqrt{2}}k(d^{sp}{3_1,u_1}), & k^(d^{sp}_{3_1,u2}) = k(d^{sp}{3_1,u2}), \ k^*(d^{sp}{3_1,u3}) &= \frac{1}{\sqrt{3}}k(d^{sp}{3_1,u3}), & k^*(d^{sp}{3_1,u4}) &= \frac{1}{\sqrt{3}}k(d^{sp}{3_1,u4}), \ k^*(d^{sp}{3_1,u5}) &= \frac{1}{\sqrt{3}\times 2}k(d^{sp}{3_1,u5}), & k^*(d^{sp}{3_1,u6}) &= \frac{1}{\sqrt{3}\times 2}k(d^{sp}{3_1,u_6}). \end{split}$$

It should be noted that weighting the kernel function based on Eq. (7) corresponds to the ''equal-split kernel function'' proposed by Ok abe et al. (2009). As shown in Fig. 2, this technique produces discontinuities at the intersections of the network; to generate smoother transitions, Okabe et al. (2009) also proposed the ''equal-split continuous kernel function''. In our study, this approach did not improve the model estimation results so we have used Eq. (7) as the weighting strategy. Next, a valid covariance function ∗ can be obtained as shown in Ver Hoef et al. (2006) as

$$c{s,s'}^* = c^*(d{s,s'}^{sp}) = Cov^*(Z(s), Z(s')),$$

$$= \sigma{\psi}^{2} \int{-\infty}^{\infty} \sqrt{w{s,u} w{s',u}} k(d{s,u}^{sp}) k(d{s',u}^{sp}) du,$$

$$\approx \sigma{\psi}^{2} \sum{j=1}^{m} \sqrt{w{s,u{j}} w{s',u{j}}} k(d{s,u{j}}^{sp}) k(d{s',u{j}}^{sp}).$$ (8)

Hence, ∗ is defined as the × weighted kernel matrix assuming knots ¹ , … , , with ∗ , = ^√, ( , ) for = 1, …, and = 1, … , .

3.1.3. Choice of the kernel function

As we have seen above, the process convolution model depends on the choice of the kernel function (.). In theory, any type of kernel can be used, as long as it is unimodal, symmetric and continuous (Okabe et al., 2009). We investigated the influence of the kernel on the covariance function based on the following functions.

Triangular ∶ () = 1 − | | ℎ |

Epanechnikov ^∶ () = ³ 4 ( 1 − ( ℎ ) for < ℎ and 0 otherwise,

Quartic ^∶ () = ¹⁵ 16( 1 − ( ℎ )2

Different from the Triangular and Quartic kernels, the Epanechnikov kernel (Epanechnikov, 1969) assigns a zero density for values of greater than the kernel width ℎ. Fig. 3 plots the empirical covariance values calculated using the three tested kernel types and three ℎ values as a function of the path distance between 100 randomly selected locations scattered on the road network shown in Fig. 1. As the panels of Fig. 3 show, kernel width controls the scale of the spatial correlation in (.) in a similar fashion to the characteristic length-scale parameter in parametric statistical covariance functions. As expected, the variance near the origin ( = 0) increases with increasing kernel widths as shown in the columns of Fig. 3. In this figure, the thick black line represents the average of the covariance values for each distance .

Generally, two approaches are used to estimate kernel width: (i) cross validation based on certain comparison criteria, including model performance for individual in-sample or out of sample data locations, and (ii) estimating this parameter along with other model parameters in a fully Bayesian framework using a Markov Chain Monte Carlo (MCMC) algorithm (Higdon, 2006). The test results (not shown) indicated that the MCMC method can successfully recover the true value of the kernel width used to simulate spatial crash data. However, estimating the kernel width under an MCMC procedure results in a prohibitive computational time depending on the size of the available data. Appendix A provides some results on the sensitivity to kernel width of the proposed process convolution model. With regard to the kernel type, it is clear from the panels of Fig. 3 that the resulting covariance does not change considerably with the three tested kernel functions. This indicates that the impact of the kernel type on the covariance, and hence on the model results, is not substantial. The same results were obtained by Xie and Yan (2008) and Bíl et al. (2013).

Fig. 3. Empirical covariance values computed using three kernel functions (rows) and varying kernel widths ℎ (columns). The solid black line represents the average of the covariance values for each distance .

3.2. Inference procedure

As mentioned above, the inference procedure is performed following the Bayesian paradigm. After specifying the prior distribution of the latent process (.), we must assign the prior distribution for the other unknowns in the model. More specifically, the parameter vector is defined as = {⁰ , ×1, ×1, ² }, resulting in a (++ 2)-dimensional vector. Following Bayes' theorem, the posterior probability distribution of , (| ), is proportional to

$$\Pi(\Theta|Y) \propto \Pi(Y|\Theta) \times \Pi(\psi \mid \sigma{\psi}^2) \times \Pi(\sigma{\psi}^2) \times \Pi(\beta_0) \times \Pi(\beta),$$ (9)

where the likelihood function ( |) is assumed to follow a Poisson distribution as in Eq. (1) and the priors are defined as

$$\beta0 \sim N(0, \sigma{\beta}^2),$$

$\betaj \sim N(0, \sigma{\beta}^2) \text{ for } j = 1, \dots, p,$
$\psi | \sigma{\psi}^2 \sim N(0, \sigma{\psi}^2),$
$\sigma_{\psi}^2 \sim \text{Gamma } (a, b).$

We fix at a large value to let the data drive the inference procedure for the coefficients . The variance parameter for the latent process plays an important role in the model, and care must be taken when fixing the parameters and of this prior distribution. In this study, we have followed the procedure recommended by Blangia rdo and Cameletti (2015) (p. 132). The technique involves acquiring prior information on the range of variation of the variance parameter [, ], simulating random variance parameter values from a uniform distribution within the prespecified range , ∼ (, ) and finally calculating the Gamma distribution parameters according to

$$a = \frac{E(\sigma{\beta,sim}^2)}{Var(\sigma{\beta,sim}^2)}$$ and $b = \frac{a}{E(\sigma_{\beta,sim}^2)}$ .

The resultant posterior distribution in Eq. (9) does not have a closed form. One possible approach is to obtain samples from it through MCMC. However, these methods can be computationally intensive, especially, when one needs to make inference about a high-dimensional parameter vector, as is typical when modelling crash data. To approximate the resultant posterior distribution, we use the integrated nested Laplace approximation (INLA) (Rue et al., 2009). Specifically, the model was implemented in the R-INLA package. More details on the implementation of the proposed algorithm in INLA is provided in Appendix B. This well-established approach approximates Bayesian inference for Latent Gaussian models; see Lindgren and Rue (2015), Blangiardo and Cameletti (2015) for more details. Before running the proposed model on the crash data described next, we compared the estimates obtained from INLA with those obtained from a MCMC approach implemented in Stan (Stan Development Team, 2016) (see Appendix C for more details). For the tests conducted in this study, INLA provided comparable results with computation times at least one order of magnitude faster than the MCMC. Hence, all the results in this study are produced with INLA.

3.3. Prediction under the network process convolution model

The network process convolution model proposed here can also be used effectively when the goal is to predict values of at unobserved locations ′ , = 1, …, where the covariates are available. Given a set of observed crash data, = ( (¹ ), …, ( ))′ , and observed covariates on the network , we denote by and the set of unobserved measurements to be predicted and the available covariates. The complete set of data is then defined as

$$Y = \begin{pmatrix} Y{obs} \ Y{new} \end{pmatrix}, \quad X = \begin{pmatrix} X{obs} \ X{new} \end{pmatrix}.$$

Given the parameter vector = {⁰ , , , ² } including as the vector of (.) values at unobserved locations , the predictive distribution of the unobserved crash response is obtained as

$$\begin{split} \Pi(Y{new}|Y{obs}, \boldsymbol{X}{obs}, \boldsymbol{X}{new}) &= \int{\Theta{C}} \Pi(Y{new}, \Theta{C}|Y{obs}, \boldsymbol{X}{obs}, \boldsymbol{X}{new}) d\Theta{C}, \ &= \int{\Theta{C}} \Pi(Y{new}|\boldsymbol{X}{new}, \Theta{C}) \Pi(\Theta{C}|Y{obs}, \boldsymbol{X}{obs}) d\Theta_{C}. \end{split} \tag{10}$$

For the th target location, ′ , the value of the spatial random effect, ( ′ ), is computed as

$$Z(si') \approx \sum{j=1}^{m} k(d_{s_i'u_j}^{sp})\psi(u_j), \tag{11}$$

where is a realization from the posterior distribution of the latent process when fitting the model to and . The kernel function ( , ) based on the path distance is evaluated between number of target locations and knot locations, , on . As discussed both the posterior and the predictive distributions are approximated through the algorithms available from INLA. For further details refer to Blangiardo and Cameletti (2015).

4. Results

This section presents the analysis of artificial data, as well as the data from the city of Ottawa in Canada. A subset of the Ottawa road network is used for our tests with 318 crash data points (Fig. 8). On the same spatial locations where the crash data are observed there were also a set of covariates from anonymized raw GPS trace data collected from clients of Intact Insurance's usage-based insurance (UBI) program (Stipancic et al., 2021). The simulated and real data sets are described in detail later in Sections 4.2 and 4.3, respectively. In each case, we compare NPC with three different models with and without latent effects. The first model assumes a usual Poisson regression (PR) model that does not account for latent effects; therefore, (.) is not included in Eq. (2). In the second model, we consider that crash data follow a Negative Binomial distribution accounting for the overdispersion of data. This model is called NB hereafter. The third model assumes that the (.)'s follow a proper conditional autoregressive (AR) model (Gelfand, 2003), a priori. The conditional autoregressive (CAR) distribution (Besag, 1974) is widely used to accommodate latent spatial structures in this kind of problems. Commonly, the spatial structure is captured through a neighbourhood structure. The issue is that the CAR model yields an improper joint distribution of the latent effects (.). Although the prior distribution is improper, it can be shown that the resultant posterior distribution is proper (Banerjee et al., 2014) if all links are connected. On a road network, this assumption means that two links are connected if they share one intersection and two intersections are neighbours if they share a link. In practice, covariates on a road network might not be collected over all links and intersections, especially if such covariates are based on mobility data. This missing pattern creates unwanted gaps among the observed data, and isolated data points with no neighbours are called singletons. As the CAR distribution can result in an improper distribution when there are singletons in the observed data, it should not be used as a prior distribution for (.). An alternative is to assume a proper CAR distribution (Banerjee et al., 2014), which is defined as

$$Z(s_i)|Z(sj){j\neq i}, \taui^{-1} \sim \mathcal{N}\left(\alpha \sum{i=j} \frac{w_{ij}Z(sj)}{w{i+}}, \frac{\tau^2}{w_{i+}}\right), \quad i, j = 1, \dots, n,$$ (12)

where defines a weighting component, and ⁺ is the sum of the weights associated with link . It is common practice to assume a 0 − 1 neighbourhood structure, such that = 1 if links and are neighbours, and 0, otherwise. The parameter in Eq. (12) controls the strength of the spatial association, where = 0 indicates spatial independence.

Under a proper CAR prior for (.), the singletons are assumed independent. This is achieved simply by reordering the data in such a way that singletons form a separate block from the connected nodes where their covariance is a diagonal matrix described with a single variance parameter (Freni-Sterrantino et al., 2018; Ver Hoef et al., 2018). In our analyses below, we denote this as the pCAR model. The proposed model, described in Section 3, is denoted as NPC (Network Process Convolution).

4.1. Model comparison criteria

The model's goodness of fit is evaluated through the Deviance Information Criterion or proposed by Spiegelhalter et al. (2002). This measure is based on two components, namely, the likelihood of observing the data given the estimated parameters ( |) and the number of parameters used to fit the model in a way that higher likelihoods and a lower number of parameters both favour lower s and hence better model performance. The based on data set is computed as

$$\begin{aligned} DIC &= D(\overline{\Theta}) + 2p_D, \ p_D &= \overline{D(\Theta)} - D(\overline{\Theta}), \end{aligned}$$

where () = −2( |) is defined as the deviance measured for the posterior means of the parameter set . The component is the complexity penalty term and refers to the effective number of parameters in the model and is defined as the deviance statistic posterior mean. In addition, we have included the Widely Applicable Information Criterion (WAIC; (Watanabe and Opper, 2010)) for model assessment and selection in this study.

The predictive performance of the model is assessed via two measures: Interval Score (; (Gneiting and Raftery, 2007)) and the mean absolute error () for and (.). A mean absolute error is chosen over the mean squared error as it measures the error on the same scale as the inputs, and (.). Predictions are carried out in the form of a K-fold cross-validation leading to one predicted crash value per data location for , = 1*,*… , .

The interval score measure, , incorporates the sharpness of the posterior predictive distribution as well as the posterior means location with respect to the actual observed values. For an individual data point located on spatial location , the interval score is computed as

$$IS^{\alpha} = \frac{1}{n}\sum{i=1}^{n} \left( (\overline{b}{u,i} - \overline{b}{l,i}) + \frac{2}{\alpha}(\overline{b}{l,i} - Y_l) \mathcal{I}{Yl < \overline{b}{l,i}} + \frac{2}{\alpha}(Yl - \overline{b}{l,i}) \mathcal{I}{Yl > \overline{b}{u,i}} \right)$$

where , and , are the lower and upper bound values of the (1−)% credible interval for data located on spatial location , and is the indicator function. A lower score indicates a narrow prediction interval, where a penalty that depends on is added if the observation misses the interval. The mean absolute error () has been calculated for both and (.). is computed between true values and the corresponding posterior means at the prediction sites, . Similarly, is computed between (.) values and the corresponding posterior means at the prediction sites, . Note that the is computed only for simulated data where the true values of (.) are known. and are computed as

$$MAEY = \frac{1}{n} \sum{i=1}^{n} \left| Y_i - \overline{Y}_i \right|,$$

$$MAEZ = \frac{1}{n} \sum{i=1}^{n} \left| Z_i - \overline{Z}_i \right|$$

where is the predicted value of at spatial location , and likewise, is the predicted value of the process (.) at spatial location . To assess the performance of the model in capturing the true association between crash data and the covariates, we have also computed the bias, , between the true values and the estimated posterior means, which is computed as

$$B_{\beta_p} = \left| \beta_p - \overline{\beta}_p \right|$$

where represents the posterior mean of the th coefficient. Similar to , we are able to compute only in the simulation study as the true values of the parameters are known. In addition, we have considered the standard deviations associated with the posterior distributions of (noted ) as a measure of the significance of estimates produced by the models.

Among the comparison criteria described above, , , and are computed based on in-sample data, i.e. the model is fitted to the observed data and the criteria are computed, while , and are computed based on out of sample data, i.e. based on a Kfold cross validation procedure. In the simulation study of Section 4.2, we have used all the comparison criteria mentioned above, and for the test on real crash data (Section 4.3), we have used the , , and as comparison criteria between the different models.

4.2. Simulation study

In this simulation study, we considered the network shown in Fig. 11 with = 318 points. This network consists of a small subset of the Ottawa data set described in the real case study in Section 4.3. Data simulation for this study is challenging, particularly because of the lack of available approaches to simulate realistic spatial correlation over simulated crash data. As noted by Ver Hoef et al. (2018), the spatial correlation can be thought of as the effect of a missing spatially structured covariate. Following this comment, crash count data were generated from a Poisson distribution with log mean ⁰ + () ′, as in Eq. (2); that is, a spatial random effect component is not included in the generation of the data. For this experiment, we set ⁰ = −3 and four covariates were considered: the first three are simulated and do not have a spatial structure, and are sampled from (0, 1) with coefficients ¹ = 2, ² = 0*.5 and ³ = −0.*5, while the fourth covariate is based on the normalized average speed taken from the real Ottawa data set described in the real case study in Section 4.3. In this data set, speed ( ) was among the covariates that displayed a clear spatial structure over the road network, which can induce a spatial correlation among the simulated data as well. The coefficient for this covariate is 1, i.e., ⁴ = 1. When fitting the models to this artificial data, the covariate with spatial structure (speed) is removed and it is expected that the spatial latent structures of models pCAR and NPC will capture this structure, which is present in the generation of the data. To investigate the performance of the different models, one hundred artificial data sets were obtained using the set of fixed parameters and input covariates introduced above. For each sample, model estimation was carried out using the INLA algorithm. The NPC model was fitted with = 173 knots, distributed evenly on the road network. The pCAR model was fitted using a first-order neighbourhood structure, which left more than 16% of the crash data locations as singletons.

The performance of the proposed NPC model in comparison with the PR, NB and pCAR models was assessed via the measures mentioned in Section 4.1 including , , , , (computed at 0.05), , and ̂ . The true values of (.) are defined based on the speed covariate used as spatial random effect. As discussed in Section 4.1, a K-fold cross-validation (K = 5) approach was used to evaluate the predictive performance of the proposed model. Panels of Figure 4 summarize the distribution of and obtained across the 100 data sets. These measurements were computed considering three different aspects of the observations: ''Y all'' — these measurements were obtained considering all observed values; ''Y singletons'' — these measurements were computed considering only the singletons; and ''Y connected'' — and were computed based only on the set of observations that are connected. As the pCAR treats the singletons in a different way than the NPC, we find this comparison important. Overall, the panels of Fig. 4 suggest that the NPC provides smaller values of and than pCAR, NB and PR. If we focus on , NPC performs better for all the three scenarios shown, i.e., ''Y all'', ''Y singletons'' and ''Y connected''. For the singletons, pCAR results in a more concentrated distribution of the but with greater values than those obtained under NPC. When the latent effects are considered (bottom row of Fig. 4), NPC performs better than pCAR for all different scenarios considered, ''Z all'', ''Z singletons'' and ''Z connected''. In other words, the latent effect present in the NPC model better approximates the missing covariate when compared to the latent effect in the pCAR model.

Fig. 5 shows the scatter plots between the true values (based on speed) used to enforce spatial correlation among the simulated crash data and the posterior mean of the spatial random effect . The different panels of this figure show scatter plots for the connected nodes (top row) and singletons (bottom row) for the predicted values at unobserved locations. The predictions are made in a cross-validation mode and results presented are from a randomly selected simulated data set. Also, as the results show, the NPC model has better recovered ( values indicated) the true values of , particularly for the singletons. Predicted values from both the NPC and pCAR models have lower variance compared to the original values; however, the NPC model has produced results with variance closer to that of the true values. In the pCAR model, the values for the singletons are drawn from normal distributions independent from the other data locations which has resulted in greater uncertainty around predictions.

Panels of Fig. 6 show the box-plots of the measured between the true coefficients and their posterior mean estimates (top row) and of the respective posterior standard deviations (bottom row) of the coefficients. In general, NPC results in the smallest values of for ⁰ and ¹ , while for ² and ³ pCAR and NPC have similar bias values. The standard deviations of the coefficients are smallest under the PR model. This is probably related to the fact that this model has fewer parameters than the others; therefore, the marginal distributions are obtained under integrals of smaller dimensions than those under the pCAR and NPC models. When we compare between pCAR and NPC, the former results in greater posterior uncertainty of the coefficients compared to NPC.

Panels of Fig. 7 show the distribution of the interval score, under the same scenarios considered in Fig. 4, and the distribution of the and across the 100 samples. The results are based on the mean interval scores computed for the 2.5 and 97.5 quantiles, i.e., for = 0*.*05. NPC performs better under the ''IS all'' scenario, resulting in overall shortest for the singletons. This is expected as the NPC model captures the correlation between singletons and other links in the network, while pCAR assumes it to be independent, a priori, from the other links in the network. The panels based on and suggests that in terms of model fitting, pCAR and NPC provides similar results.

4.3. Analysis of crash data in Ottawa, Canada

We have used a small data set collected in the city of Ottawa, Canada. The data set used here covers the Ottawa neighbourhood of Centretown, as shown in Fig. 8, accounting for 318 crash vehicle (collisions) data locations where contextual variables from GPS trace data were readily available. Covariates were extracted from GPS trace data collected over one month (April 2016) by Intact Insurance as part of their UBI program. The GPS data were collected using dedicated dongles connected to the on-board diagnostics port of clients' vehicles with access to the vehicle speed sensor. The crash count data were obtained for a period of five years (2013 to 2018). The covariates include the number of trips per link (log scale), the length of the link in km (log scale), the average deceleration rate computed over all drivers passing through the link (identified as drate here), the congestion index (CI), the average speed (V), the coefficient of variation of speed (CVS),

Fig. 4. Box-plots of the distribution of (top row) and (bottom row) based on the 100 data sets obtained for the simulation study.

Fig. 5. Scatter plots between true values and the predicted posterior means of along with the 95% credible intervals. Scatter plots are displayed separately for singletons and connected data.

Fig. 6. Box-plots of the distribution of (top row) and ̂ (bottom row) based on the 100 data sets obtained for the simulation study.

and the road type. Road type was an indicator variable for primary (0) vs. highway (1). Readers are referred to Stipancic et al. (2018) for a complete description of the covariates and Stipancic et al. (2021) for an overview of the observed crash count data. Table 2 provides some descriptive statistics for the count crash data and the selected covariates. Over the observation period, the average crashes per data location is 0.8 on the sample road network, while the maximum number of observed crashes is 20. Note that the crash data are over-dispersed with a coefficient of variation of 2.45.

All four models were fitted to the Ottawa data set, using the INLA algorithms implemented in the R-INLA package, with the same parameters (neighbouring structure and number of knots) described in

Fig. 7. Box plots of the distribution of and based on the 100 simulated data sets. The labels on the subfigures pertaining to stand for the following: ''IS all'' was obtained considering all observed values; ''IS singletons'' — was computed considering only the singletons; and ''IS connected'' — was computed based only on the set of observations that are connected.

Fig. 8. Subset of the road network in Ottawa used in this study highlighted in orange (left). Observed crash data in the form of the bubble plot overlaid on the network (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2 Descriptive statistics of the outcome and covariates used in this study.

	Minimum	Mean	Max	Std. dev.
Crash data (Y)	0	0.84	20	2.06
ln(trips)	0.76	3.35	5.53	1.21
ln(length)	2.56	6.14	8.39	1.09
drate	0.00	0.037	0.31	0.053
CI	0.00	0.09	0.94	0.14
V	0.94	16.23	40.99	7.62
CVS	0.17	0.95	5.33	0.45
Road type	0.00	0.61	1.00	0.48

the simulation study. Table 3 provides the posterior summaries of the coefficients of the covariates under the different models. Among the covariates considered, only link length, CVS and road type result in 95% posterior credible intervals that do not include zero.

Table 4 presents the values of , , and obtained from the four models fitted to the Ottawa data set. values and were computed based on the fit and considering the entire data set (one model run) while and were computed through a K-fold cross-validation exercise after acquiring predicted values and the credible intervals for each data location. Clearly, pCAR results in the smallest value of and when compared with NPC, PR and NB. Therefore, in terms of model fitting, according to , pCAR performs better than NPC. However, in terms of prediction according to , NPC performs best, while in terms of , NPC performs second best. The reason why pCAR performs worse in terms of prediction is probably due to the fact that the latent effects at the singletons are considered to follow independent normal distributions a priori. In other words, the singletons do not borrow information from other locations, leading to greater uncertainty in the prediction.

To further visualize this issue regarding the higher posterior uncertainty at the singletons under the pCAR model, the panels of Fig. 9 show the posterior distribution of the latent effects (.) for three randomly selected singletons (53, ¹⁶⁷ and 148), and connected crash locations (83, ³⁷ and 135), where refers to the exact spatial location. The posterior distributions based on the pCAR model (solid lines) are associated with very large uncertainties around the singletons when compared to those obtained under the NPC model (dashed lines). Furthermore, for some locations, the resultant posterior distributions under the NPC model tend to be further away from 0 than their counterparts under the pCAR model. This suggests that the latent effect (.) under the NPC is better at capturing the local structure of the network, as it is able to borrow strength from other locations even if the location is a singleton. Under the pCAR model, the prior specification of a singleton is independent from the other locations, therefore, there is no borrowing of strength from other locations.

To evaluate the predictive performance of the proposed model, we have used a K-fold cross-validation (K = 3). Fig. 10 shows the scatter plots between the observed values at the holdout links versus the predicted values (posterior means), together with the limits of the 95% credible intervals (grey segments).

The panels of Fig. 11 show the posterior mean estimates of (.) for the set of data locations (left) and over the entire network (right). This interpolated map can be seen as a proxy measure for risk and be

Table 3 Posterior summaries (mean and 95% credible intervals in brackets) of the coefficients .

	PR	NB	pCAR	NPC
Intercept	−14.03 (−16.95;−11.23)	−14.09 (−17.46;−14.057)	−7.86 (−10.11;−5.66)	−8.17 (−10.31;−6.07)
ln(trips)	2.49 (2.16;2.83)	2.55 (2.14;3.01)	1.97 (1.66;2.29)	1.97 (1.69;2.26)
ln(length)	0.56 (0.29;0.84)	0.53 (0.20;0.86)	0.06 (−0.24;0.36)	0.11 (−0.18;0.39)
drate	0.15 (−2.36;2.64)	0.30 (−3.65;2.934)	−0.22 (−1.92;1.48)	−0.33 (−1.96;1.30)
CI	0.76 (−0.34;1.81)	1.02 (−0.24;2.28)	0.44 (−0.65;1.49)	0.34 (−0.66;1.32)
V	−0.11 (−0.69;0.36)	−0.18 (−0.91;0.40)	−0.60 (−1.32;−0.03)	−0.55 (−1.21;−0.01)
CVS	−0.06 (−0.10;−0.02)	−0.06 (−0.11;−0.023)	−0.10 (−0.15;−0.05)	−0.08 (−0.13;−0.04)
Highway/primary	0.63 (0.07;1.18)	0.82 (0.12;1.54)	1.40 (0.73;2.10)	1.07 (0.46;1.69)

Fig. 9. Marginal posterior distribution of (.) for singletons (top row) and connected data locations (bottom row) for pCAR (solid) and NPC (dashed) models.

Fig. 10. Scatter plots between observed data and the predicted posterior means along with the 95% credible intervals.

Table 4 , and values computed over Ottawa data set.

	PR	NB	pCAR	NPC
𝐷𝐼𝐶 (in sample)	488.67	469.45	456.78	470.81
𝑊 𝐴𝐼𝐶 (in sample)	495.37	469.98	455.51	485.45
𝑀𝐴𝐸𝑌 (out of sample)	0.68	0.86	0.66	0.60
𝐼𝑆 (out of sample)	0.39	0.36	0.18	0.24

used for site-was ranking purposes. It can also help highway authorities to better identify high-risk areas within the road network, where contextual variables are not available.

5. Discussion and conclusion

In this paper, we introduced a model for crash data that incorporates the structure of the road network using a process convolution approach to describe the latent random effect that is commonly included in the mean of a Poisson distribution. By modelling the complex spatial relationship between the number of crashes observed on a road network, the proposed model enables us to estimate crash intensity surfaces that can be used to identify hotspot regions within the road network. In terms of practical application, these measures can be used to rank and identify the high-risk locations in the network. In the public sector, road safety professionals aim to allocate their limited budgets to those sites in the network where improvements will produce significant crash reductions. By improving network screening (increasing the accuracy of identified crash-intensive locations within the road network), funds can be allocated more efficiently, yielding further reductions in both the monetary and social costs associated with road traffic crashes. In the private sector, the spatial distribution of traffic crashes could be used as a feature for improving the accuracy of risk pricing.

Different from other approaches available in the literature, e.g., the pCAR model, the proposed model provides spatial interpolation to unobserved locations of interest in a straightforward fashion. This is because the latent processes for observed and unobserved locations come from the same multivariate normal distribution, regardless of whether the locations are singletons or connected to other links. In both public and private practice, where data sparsity is often a concern, this contribution is especially important. In network screening, for

Fig. 11. [Left] Posterior mean of the latent effect of the NPC model at the observed links. [Right] Spatial interpolation of the spatial process (.) of the NPC model over the entire network.

example, either the independent variables (traffic data) or dependent variables (crash data) are unobserved (due to low traffic volumes or unreported crashes) at specific sites throughout the network. In this approach, even these unobserved locations could be candidates for safety remediation, giving practitioners a more complete view of risk throughout the network.

The proposed model incorporates the structure of the road network into crash modelling through the use of the path distance rather than the Euclidean distance. However, in practice the outcome of interest, e.g. crash frequency or crash severity injury level, is influenced not only by factors operating within the road network, i.e. on the road structure, but also by other factors that may operate in the Euclidean space. For example, an increase in crash frequency might be due to the deterioration of the road surface due to the geological properties of the surface sediments, which does not follow a network distance but rather a Euclidean distance, at least on the scale of a metropolitan city. As discussed by Ver Hoef et al. (2018) and others (Dale and Fortin, 2010; Peterson et al., 2013), one option is to consider a variance component model measuring the correlations in both spaces separately (Ver Hoef et al., 2018). It is, however, beyond the scope of this work to investigate this subject further.

The pCAR model that was fit to the data in this study is based on a first-order neighbourhood structure. The advantage of following this approach is that this neighbourhood structure results in a sparse precision matrix, reducing the computational burden of the inference procedure. Our goal in including the CAR model among the fitted models was to provide a comparison with what is commonly done in the literature. However, one can change the neighbourhood structure of the CAR model to one that uses the path distance. However, some care should be taken as it would be interesting to have a parameter controlling for the decay of the weights as the path distance increases. This will increase the complexity of the model as we would have an extra parameter to estimate in the neighbourhood structure of the CAR model. As crash data are commonly high-dimensional, this can have a great impact on the computation time needed to perform the inference procedure.

From the perspective of traffic stream and flow direction, other structure-related factors may come into play, such as the segment types that can contribute to the spatial pattern of crashes observed on a road network. We consider the proposed approach a first attempt to account for the road network structure when including latent spatial effects to describe the number of crashes as a function of covariates. We plan to extend this work to account for directional effects that lead to a non-stationary covariance structure of the latent spatial effect . A possibility is to borrow ideas from Ver Hoef and Peterson (2010), wherein they account for the direction of the (stream) flow when defining the latent effects of the model. An alternative is to borrow ideas from Schmidt and Rodriguez (2011) and Neto et al. (2014) and make use of covariate information in the covariance structure of the latent effect.

Typically, spatial models are limited in practice because of scalability; the rapidly increasing computation time associated with highdimensional spatial models make them difficult to implement at the scales necessary in practice. The proposed model can be easily scaled up to applications with a large number of crash data points, such as an entire metropolitan area (necessary for network screening) or an entire province or country (necessary for risk pricing). This is because the spatial process is approximated through a kernel convolution approach in a non-parametric fashion. To define spatial random effects (.) we do not compute covariance matrices; rather, we define it through a linear combination of the latent process values at a limited set of knots in the vicinity of the data location (depending on the chosen kernel width). This results in a considerable reduction of CPU time to define (.), and consequently to estimate/approximate the resultant posterior distribution using MCMC or INLA. Therefore, the CPU time used to fit the proposed model is largely determined by the number of knots used to define (.). Furthermore, the INLA implementation of the proposed model offers a considerable computational advantage over MCMC implementation. Appendix C provides a comparison between INLA and MCMC computation times for a varying number of knots.

Using path distances comes with both theoretical and empirical challenges. The theoretical challenge is related to the definition of a valid model based on the non-Euclidean shortest-path distances rather than the original Euclidean distances. Our aim in this study has been to tackle this problem by proposing the NPC model. The empirical challenge is the calculation of shortest-path distances and the associated weights between sets of locations over a road network. The calculation of the distances over a large network can be computationally demanding; however, due to the numerous splittings, the weights computed between locations far apart on a road network become close to zero, hence, their contribution to the overall kernel density would become negligible. As a result, the distances do not need to be computed between all locations. Also, the distance and weight matrices need to be computed only once, in advance, and be used in all the following modelling steps.

Lastly, given that the crash count data sets typically contain many zeros, other models such as zero-inflated models can be used. However, the zero-component of such models must also account for the spatial correlation present in the data. The proposed NPC model can be used to describe the spatial correlation in both zero and non-zero components. This is the subject of future research.

CRediT authorship contribution statement

Hassan Rezaee: Methodology, Investigation, Experiments design, Software - designing computer programs; implementation of the computer code and supporting algorithms; testing of existing code components, Validation, Writing – original draft. Alexandra M. Schmidt: Conceptualization, Methodology, Formal analysis – application of statistical techniques, Writing – review & editing, Supervision. Joshua Stipancic: Data curation, Formal analysis – synthesize study data, Validation, Writing – review & editing. Aurélie Labbe: Conceptualization, Methodology, Experiments design, Formal analysis – application of statistical techniques, Writing – review & editing, Supervision, Project administration, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The authors do not have permission to share data.

Acknowledgements

This research was funded by the FRQ-IVADO Chair in Data Science at HEC Montréal and Mitacs Elevation program (grant ID IT21937). The authors thank Intact Insurance for providing the data and collaborating in this study. We also thank the anonymous reviewers whose constructive comments helped improve the work.

Appendix A. Sensitivity to kernel width

To analyse the effect of kernel width (ℎ) on results obtained from the proposed model 100 data sets were generated from the base model. The data were produced using a set of fixed coefficients, kernel width (indicated as ℎ⁰ ) and latent process values at number of knot locations. The model parameters in = {⁰ , , , ² } were estimated in INLA using varying kernel widths ℎ forming ten *ℎ*∕ℎ⁰ ratios equally spaced between 0.25 to 4. Fig. A1 shows the effect of the kernel width ratio on results obtained for both in sample and out of sample (prediction mode) data locations. In sample refers to the training data used for modelling fitting and out of sample data refer to data held out (through a K-fold cross validation procedure) and used for testing/prediction only. On this figure the vertical grey lines indicate the 95% credible intervals. For each *ℎ*0∕ℎ ratio we have tested the algorithm for all the 100 data sets and the 2.5-th, 50-th and 97.5-th percentiles were computed. The predictions for out of sample data locations are made using a K-fold cross validation process. The panels of Fig. A1 unanimously show that using kernel width that is too small compared to the width used to generate the data can negatively influence the results by increasing the MAE's computed for crash values as well as the latent random effect (.). The results are, however, robust to kernel widths up to 4 times larger than the true width used to simulate the data. The same pattern can be seen for , , and shown in Figs. A2 and A3. We carried out a similar test for different kernel functions and its effect on model results which were found insignificant as supported by Fig. 3 of the main text.

Appendix B. INLA implementation

One of the main advantages of the proposed model is that it can be easily implemented in R-INLA. The R-INLA package proposes a wide variety of models for data likelihood and spatial random effects. The kernel convolution model of Eq. (4) is known as the z model, i.e. = . It should be noted that the notation (.) used for the design matrix in the INLA implementation corresponds to the matrix in our model. The INLA formulation corresponding to our NPC model is therefore

formula = crash ~ 1 + covar1 + covar2 + f(id, model = "z", Z = K),

assuming an intercept and two covariates named covar1 and covar2 in the model. The most important part of this implementation is the computation of the matrix that accounts for the structure of the road network.

The other two models (pCAR and PR) are also implemented in INLA. The implementation of the PR and pCAR model is straightforward. The former does not include a random effect, and the fixed effects are introduced using the formula statement, e.g.,

formula = crash ~ 1 + covar1 + covar2,

where 1 defines the intercept of the model. The pCAR implementation is similar to the PR with an additional spatial random effect. In this case the INLA formula will become

formula = crash ~ 1 + covar1 + covar2 + f(id, model

= "besagproper", graph = graph)

where the besagproper represents the pCAR prior on the random effects with the binary adjacency matrix given as graph. The R code of all the three models implemented in INLA are available on GitHub (https://github.com/aurelielabbe/Process-convolution-model for-crash-counts-data-on-a-network).

Appendix C. Comparison between INLA and MCMC

The posterior densities approximated with INLA were compared with the estimated densities from MCMC. We have generated 100 simulated data sets and estimated/approximated model parameters using both MCMC and INLA for each data set. Results are reported in Figs. B1, B2 and B3 comparing MAEs for and (.) for in sample and out of sample data locations (prediction mode) as well as the , and . Fig. B4 displays sample posterior distributions acquired from MCMC and INLA for coefficients and (.) values for four randomly selected locations. The results clearly show that INLA is equivalent to MCMC for parameter estimation in the current study. We also tested the case where crash data were statistically very rare on the road network, and the same results were obtained; although, both methods were shown to be less accurate in recovering the true values of the parameters.

The advantage of using INLA is considerable given that the computation burden of using INLA is orders of magnitude lower than MCMC. For example, the present study incorporated 318 data locations and the analysis was completed in less than 30 s (CPU time). The CPU time for a larger data set ( = 20,000) was about 4.7 h on a low-cost laptop equipped with an i7 2.40 GHz processor with 8 GB of memory. For comparison purposes, we ran the MCMC implementation of our model in Stan (Stan Development Team, 2016) for = 20,000 and the CPU time was on the order of days. Fig. B5 shows the CPU time comparison between INLA and MCMC implementation of NPC. In this figure, CPU times on the vertical axis are plotted on the logarithmic scale. For each on the X axis, model fitting was done three times based on different number of knots: = 0*.25 (dashed), = 0.5 (solid) and = 0.*75 (dotted). For each value on the X axis of this plot, the MCMC sampling was continued until at least an effective sample size of 500 was acquired. For this test, on average INLA is about 35 times faster than the MCMC while both implementations produced comparable results.

Fig. A1. Sensitivity of the model to kernel width ratio measured in the form the in sample and out of sample MAE rates for and (.).

Fig. A2. Sensitivity of the model to kernel width ratio measured in bias values . The vertical grey lines indicate the 95% confidence intervals.

Fig. A3. Sensitivity of the model to kernel width ratio measured in the and values. The vertical grey lines indicate the 95% confidence intervals.

Fig. B1. Box-plots of the distribution of based on the 100 data sets obtained for the simulation study for both the MCMC and INLC. measures are computed separately for each coefficient. Values are plotted on a logarithmic scale for better comparison.

Fig. B2. Box-plots of the distribution of MAE and MAE based on the 100 simulated data sets. The labels on top of the sub-figures stand for the following: ''in sample'': MAE was obtained considering observed values and the fitted values; ''out of sample'': MAE was computed considering the test set and the predicted values.

Fig. B3. Box plots of the distribution of and based on the 100 simulated data sets.

Fig. B4. Marginal posterior distribution of s (top row) and (.) at four data locations (bottom row) for MCMC (histograms) and INLA (density lines). The locations of the true values are indicated with the blue vertical lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. B5. CPU time comparison between INLA (black) and MCMC (grey) implementation of the proposed algorithm.

References

Abduljabbar, R., Dia, H., Liyanage, S., Bagloee, S.A., 2019. Applications of artificial intelligence in transport: An overview. Sustainability 11 (1), 189.
Aguero-Valverde, J., 2014. Direct spatial correlation in crash frequency models: Estimation of the effective range. J. Transp. Saf. Secur. 6 (1), 21–33. http://dx. doi.org/10.1080/19439962.2013.799108.
Aguero-Valverde, J., Jovanis, P.P., 2008. Analysis of road crash frequency with spatial models. Transp. Res. Rec. 2061 (1), 55–63.
Banerjee, S., Carlin, B.P., Gelfand, A.E., 2014. Hierarchical Modeling and Analysis for Spatial Data. CRC Press.
Bao, J., Liu, P., Ukkusuri, S.V., 2019. A spatiotemporal deep learning approach for citywide short-term crash risk prediction with multi-source data. Accid. Anal. Prev. 122, 239–254. http://dx.doi.org/10.1016/j.aap.2018.10.015.
Besag, J., 1974. Spatial interaction and the statistical analysis of lattice systems. J. R. Stat. Soc. Ser. B Stat. Methodol. 36 (2), 192–225. http://dx.doi.org/10.1111/j.
Besag, J., York, J., Mollié, A., 1991. Bayesian image restoration, with two applications in spatial statistics. Ann. Inst. Statist. Math. 43 (1), 1–20. http://dx.doi.org/10. 1007/bf00116466.
Bíl, M., Andrášik, R., Janoška, Z., 2013. Identification of hazardous road locations of traffic accidents by means of kernel density estimation and cluster significance evaluation. Accid. Anal. Prev. 55, 265–273. http://dx.doi.org/10.1016/j.aap.2013. 03.003.
Blangiardo, M., Cameletti, M., 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. John Wiley & Sons.
Cai, Q., Abdel-Aty, M., Lee, J., Huang, H., 2019. Integrating macro- and micro-level safety analyses: A Bayesian approach incorporating spatial interaction. Transportmetrica A: Transp. Sci. 15 (2), 285–306. http://dx.doi.org/10.1080/23249935. 2018.1471752.
Castro, M., Paleti, R., Bhat, C.R., 2013. A spatial generalized ordered response model to examine highway crash injury severity. Accid. Anal. Prev. 52, 188–203. http: //dx.doi.org/10.1016/j.aap.2012.12.009.
Chiou, Y.-C., Fu, C., Chih-Wei, H., 2014. Incorporating spatial dependence in simultaneously modeling crash frequency and severity. Anal. Methods Accid. Res. 2, 1–11. http://dx.doi.org/10.1016/j.amar.2013.12.001.
Cressie, N., Frey, J., Harch, B., Smith, M., 2006. Spatial prediction on a river network. J. Agric. Biol. Environ. Stat. 11 (2), 127–150. http://dx.doi.org/10.1198/ 108571106x110649.
Cressie, N., Majure, J.J., 1997. Spatio-temporal statistical modeling of livestock waste in streams. J. Agric. Biol. Environ. Stat. 2 (1), 24. http://dx.doi.org/10.2307/1400639. Csardi, G., Nepusz, T., et al., 2006. The igraph software package for complex network research. InterJournal, Complex Syst. 1695 (5), 1–9.
Curriero, F.C., 2006. On the use of non-Euclidean distance measures in geostatistics. Math. Geol. 38 (8), 907–926. http://dx.doi.org/10.1007/s11004-006-9055-7.
Dale, M., Fortin, M.-J., 2010. From graphs to spatial graphs. Annu. Rev. Ecol. Evol.
Syst. 41 (1), 21–38. http://dx.doi.org/10.1146/annurev-ecolsys-102209-144718. Dijkstra, E.W., 1959. A note on two problems in connexion with graphs. Numer. Math. 1 (1), 269–271. http://dx.doi.org/10.1007/bf01386390.
Effati, M., Thill, J.-C., Shabani, S., 2015. Geospatial and machine learning techniques for wicked social science problems: analysis of crash severity on a regional highway corridor. J. Geogr. Syst. 17 (2), 107–135. http://dx.doi.org/10.1007/s10109-015- 0210-x.
El-Basyouny, K., Sayed, T., 2009. Urban arterial accident prediction models with spatial effects. Transp. Res. Rec. 2102 (1), 27–33.
Epanechnikov, V., 1969. Nonparametric estimation of a multidimensional probability density. Teoriya Veroyatnostei I Ee Primeneniya 14 (1), 156–161.
Fotheringham, A.S., Brunsdon, C., Charlton, M., 2003. Geographically weighted regression: The analysis of spatially varying relationships. Geogr. Anal. 35 (3), 272–275. http://dx.doi.org/10.1111/j.1538-4632.2003.tb01114.x.
Freni-Sterrantino, A., Ventrucci, M., Rue, H., 2018. A note on intrinsic conditional autoregressive models for disconnected graphs. Spatial Spatio-Temporal Epidemiology 26, 25–34. http://dx.doi.org/10.1016/j.sste.2018.04.002.
Gardner, B., Sullivan, P.J., Arthur J. Lembo, J., 2003. Predicting stream temperatures: geostatistical model comparison using alternative distance metrics. Can. J. Fish. Aquat. Sci. 60 (3), 344–351. http://dx.doi.org/10.1139/f03-025.
Gelfand, A.E., 2003. Proper multivariate conditional autoregressive models for spatial data analysis. Biostatistics 4 (1), 11–15. http://dx.doi.org/10.1093/biostatistics/4. 1.11.
Gneiting, T., Raftery, A.E., 2007. Strictly proper scoring rules, prediction, and estimation. J. Amer. Statist. Assoc. 102 (477), 359–378. http://dx.doi.org/10.1198/
Gomes, M.J.T.L., Cunto, F., da Silva, A.R., 2017. Geographically weighted negative binomial regression applied to zonal level safety performance models. Accid. Anal. Prev. 106, 254–261. http://dx.doi.org/10.1016/j.aap.2017.06.011.
Hadayeghi, A., Shalaby, A.S., Persaud, B., 2003. Macrolevel accident prediction models for evaluating safety of urban transportation systems. Transp. Res. Rec.: J. Transp. Res. Board 1840 (1), 87–95. http://dx.doi.org/10.3141/1840-10.
Heaton, M.J., Datta, A., Finley, A.O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R.B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D.W., Sun, F., Zammit-Mangion, A., 2018. A case study competition among methods for analyzing large spatial data. J. Agric. Biol. Environ. Stat. 24 (3), 398–425. http://dx.doi.org/ 10.1007/s13253-018-00348-w.
Higdon, D., 2002. Space and space-time modeling using process convolutions. In: Quantitative Methods for Current Environmental Issues. Springer, pp. 37–56.
Higdon, D., 2006. A primer on space-time modeling from a Bayesian perspective. Monogr. Statist. Appl. Probab. 217–279. http://dx.doi.org/10.1201/ 9781420011050.ch6.
Higdon, D., Swall, J., Kern, J., 1999. Non-stationary spatial modeling. Bayesian Statistics 6. In: Proceedings of the Sixth Valencia International Meeting. pp. 761–768.
Jensen, O.P., Christman, M.C., Miller, T.J., 2006. Landscape-based geostatistics: A case study of the distribution of blue crab in chesapeake bay. Environmetrics 17 (6), 605–621. http://dx.doi.org/10.1002/env.767.
Ladle, A., Avgar, T., Wheatley, M., Boyce, M.S., 2016. Predictive modelling of ecological patterns along linear-feature networks. In: McCrea, R. (Ed.), Methods Ecol. Evol. 8 (3), 329–338. http://dx.doi.org/10.1111/2041-210x.12660.
Lindgren, F., Rue, H., 2015. Bayesian spatial modelling withR-INLA. J. Stat. Softw. 63 (19), http://dx.doi.org/10.18637/jss.v063.i19.
Little, L.S., Edwards, D., Porter, D.E., 1997. Kriging in estuaries: as the crow flies, or as the fish swims? J. Exp. Mar. Biol. Ecol. 213 (1), 1–11. http://dx.doi.org/10.1016/ s0022-0981(97)00006-3.
Liu, J., Khattak, A.J., Wali, B., 2017. Do safety performance functions used for predicting crash frequency vary across space? Applying geographically weighted regressions to account for spatial heterogeneity. Accid. Anal. Prev. 109, 132–142. http://dx.doi.org/10.1016/j.aap.2017.10.012.
Mannering, F.L., Bhat, C.R., 2014. Analytic methods in accident research: Methodological frontier and future directions. Anal. Methods Accid. Res. 1, 1–22. http: //dx.doi.org/10.1016/j.amar.2013.09.001.
Neto, J.H.V., Schmidt, A.M., Guttorp, P., 2014. Accounting for spatially varying directional effects in spatial covariance structures. J. R. Stat. Soc. Ser. C. Appl. Stat. 63 (1), 103–122.
Okabe, A., Satoh, T., Sugihara, K., 2009. A kernel density estimation method for networks, its computational method and a GIS-based tool. Int. J. Geogr. Inf. Sci. 23 (1), 7–32. http://dx.doi.org/10.1080/13658810802475491.
Peterson, E., Ver Hoef, J., Isaak, D., Falke, J., Fortin, M., Jordan, C., McNyset, K., Monestiez, P., Ruesch, A., Sengupta, A., et al., 2013. Stream networks in space:
concepts, models, and synthesis. Ecol. Lett. 16 (5), 707–719, https://doi.org/false. Rathbun, S.L., 1998. Spatial modelling in irregularly shaped regions: kriging estuaries. Environmetrics: Official J. Int. Environmetrics Soc. 9 (2), 109–129.
Ren, H., Song, Y., Wang, J., Hu, Y., Lei, J., 2018. A deep learning approach to the citywide traffic accident risk prediction. In: 2018 21st International Conference on Intelligent Transportation Systems. ITSC, IEEE, pp. 3346–3351.
Rue, H., Martino, S., Chopin, N., 2009. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J. R. Stat. Soc. Ser. B Stat. Methodol. 71 (2), 319–392. http://dx.doi.org/10.1111/j.1467- 9868.2008.00700.x.
Satria, R., Aguero-Valverde, J., Castro, M., 2021. Spatial analysis of road crash frequency using Bayesian models with integrated nested Laplace approximation (INLA). J. Transp. Saf. Secur. 13 (11), 1240–1262.
Schmidt, A.M., Rodriguez, M.A., 2011. Modelling multivariate counts varying continuously in space. Bayesian Stat. 9, 611–638.
Selby, B., Kockelman, K.M., 2013. Spatial prediction of traffic levels in unmeasured locations: applications of universal kriging and geographically weighted regression. J. Transp. Geogr. 29, 24–32. http://dx.doi.org/10.1016/j.jtrangeo.2012.12.009.
Shiode, N., Shiode, S., 2011. Street-level spatial interpolation using network-based IDW and ordinary kriging. Trans. GIS 15 (4), 457–477. http://dx.doi.org/10.1111/j. 1467-9671.2011.01278.x.
Spiegelhalter, D.J., Best, N.G., Carlin, B.P., van der Linde, A., 2002. Bayesian measures of model complexity and fit. J. R. Stat. Soc. Ser. B Stat. Methodol. 64 (4), 583–639. http://dx.doi.org/10.1111/1467-9868.00353.
Stan Development Team, 2016. Stan Modeling Language Users Guide and Reference Manual. Technical Report.
Stipancic, J., Miranda-Moreno, L., Saunier, N., Labbe, A., 2018. Surrogate safety and network screening: Modelling crash frequency using GPS travel data and latent Gaussian spatial models. Accid. Anal. Prev. 120, 174–187. http://dx.doi.org/10. 1016/j.aap.2018.07.013.
Stipancic, J., Miranda-Moreno, L., Saunier, N., Labbe, A., 2019. Network screening for large urban road networks: Using GPS data and surrogate measures to model crash frequency and severity. Accid. Anal. Prev. 125, 290–301. http://dx.doi.org/ 10.1016/j.aap.2019.02.016.
Stipancic, J., Racine, E.B., Labbe, A., Saunier, N., Miranda-Moreno, L., 2021. Massive GNSS data for road safety analysis: Comparing crash models for several Canadian cities and data sources. Accid. Anal. Prev. 159, 106232. http://dx.doi.org/10.1016/ j.aap.2021.106232.
Thiébaux, H., Pedder, M., 1987. Spatial Objective Analysis: With Applications in Atmospheric Science. Academic Press, San Diego.
Ver Hoef, J.M., 2018. Kriging models for linear networks and non-Euclidean distances: Cautions and solutions. Methods Ecol. Evol. 9 (6), 1600–1613. http://dx.doi.org/
Ver Hoef, J.M., Peterson, E.E., 2010. A moving average approach for spatial statistical models of stream networks. J. Amer. Statist. Assoc. 105 (489), 6–18. http://dx.doi. org/10.1198/jasa.2009.ap08248.
Ver Hoef, J.M., Peterson, E.E., Hooten, M.B., Hanks, E.M., Fortin, M.-J., 2018. Spatial autoregressive models for statistical inference from ecological data. Ecol. Monograph 88 (1), 36–59. http://dx.doi.org/10.1002/ecm.1283.
Ver Hoef, J.M., Peterson, E., Theobald, D., 2006. Spatial statistical models that use flow and stream distance. Environ. Ecol. Stat. 13 (4), 449–464. http://dx.doi.org/ 10.1007/s10651-006-0022-8.
Watanabe, S., Opper, M., 2010. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.. J. Mach. Learn. Res. 11 (12).
World Health Organization, 2018. Global Status Report on Road Safety 2018: Summary. tech. rep., World Health Organization.
Xie, Z., Yan, J., 2008. Kernel density estimation of traffic accidents in a network space. Comput. Environ. Urban Syst. 32 (5), 396–406. http://dx.doi.org/10.1016/ j.compenvurbsys.2008.05.001.
Ziakopoulos, A., Yannis, G., 2020. A review of spatial approaches in road safety. Accid. Anal. Prev. 135, 105323. http://dx.doi.org/10.1016/j.aap.2019.105323.

rezaeeprocessconvolutionmodelcrash2022.md 104 KB Cronologia Originale