Obsidian/Zotero/001_artiklid/Road Safety Indicators/cuiacceleratedhierarchicalbayesiancrash2021.md

---

category: literaturenote citekey: cuiacceleratedhierarchicalbayesiancrash2021 title: An accelerated hierarchical Bayesian crash frequency model with accommodation of spatiotemporal interactions authors: "Cui, Haipeng; Xie, Kun" year: 2021 date: 2021-00-00 2021 doi: 10.1016/j.aap.2021.106018 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2021.106018" zotero_key: JVZNVMFP zotero_storage: 4DVUH564 collections: transpordi_planeerimine folder: 001_artiklid/Road Safety Indicators firstAuthor: "Cui, Haipeng"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

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

An accelerated hierarchical Bayesian crash frequency model with accommodation of spatiotemporal interactions

Haipeng Cui ^a , Kun Xie b, *

• ^a Department of Civil and Environmental Engineering, National University of Singapore, 117576, Singapore
• ^b Department of Civil & Environmental Engineering, Old Dominion University, Norfolk, VA 23529, USA

ARTICLE INFO

Keywords: Crash frequency Spatiotemporal interaction Integrated nested Laplace approximation Bayesian approach

ABSTRACT

Although spatial and temporal correlations of crash observations have been well addressed in the literature, the interactions between them are rarely studied. This study proposes a Bayesian spatiotemporal interaction (BSTI) approach for crash frequency modeling with an integrated nested Laplace approximation (INLA) method to greatly expedite the Bayesian estimation process. Manhattan, which is the most densely populated urban area of New York City, is selected as the study area. Hexagons are used as the basic geographic units to capture crash, transportation, land use, and demo-economic data from 2013 to 2019. A series of Bayesian models with various spatiotemporal specifications are developed and compared. The BSTI model with Type II interaction, which assumes that the structured temporal random effect interacts with the unstructured spatial random effect is found to outperform the others in terms of goodness-of-fit and the ability to reduce the dependency of residuals. It is also found that the unobserved heterogeneity is mostly attributed to the spatial effects instead of temporal effects. In addition, the BSTI Type II model also yields the lowest predictive error when the last year's data are used as the test set. The proposed BSTI approach can potentially advance safety analytics by achieving high prediction accuracy and computational efficiency while maintaining its interpretability on the effects of contributing factors and the unobserved heterogeneity.

1. Introduction

Road traffic crashes result in approximately 1.35 million deaths and up to 50 million injuries each year globally (WHO, 2018). Crash frequency modeling can assist informed decision-making in road safety management such as hotspot identification (Cheng et al., 2017; Krueger et al., 2020) and before-after evaluation of safety countermeasures (Elvik et al., 2017; Tang et al., 2020). Traditional models such as Poisson models (Jones et al., 1991), Poisson-gamma models (Poch and Mannering, 1996) and Poisson-lognormal models (Lord and Miranda-Moreno, 2008) have been used to model crash frequencies, with the assumption that the crash observations are independent from each other. However, this assumption can be easily violated. Firstly, spatial correlation can exist in crash observations of neighboring sites. Risk factors that affect the likelihood of crashes like traffic demand, road design geometry and socioeconomic features tend to be similar among neighboring sites. Secondly, temporal correlation can exist in crash observations of consecutive time periods. For example, work zones that last for a few years would interrupt traffic operations and lead to more crashes in those years. Moreover, traffic demand, which is largely associated with population and economic activities, is generally temporally correlated and thus it could result in similar exposure to crashes in adjacent years. Thirdly, there could be spatiotemporal interaction (STI) of crash data collected in neighboring sites and consecutive time periods. STI is defined as the changing pattern of the spatial/temporal correlation of the data as time passes/location changes (Blangiardo and Cameletti, 2015). For instance, the spatial correlation among crash data in off-peak hours would be relatively weak, while would be intensified in peak hours due to queue propagation.

Despite the methodological advancements during the past few years in accounting for: (1) spatial correlation (Lord et al., 2000; Narayanamoorthy et al., 2013); (2) temporal correlation (Wang and Abdel-Aty, 2006; Wang et al., 2019); and (3) spatial and temporal correlations (Persaud et al., 2010; Ulak et al., 2018; Liu et al., 2019), addressing the interaction between spatial and temporal correlations is rarely studied. Additionally, numerous existing studies adopt Markov Chain Monte Carlo (MCMC) method to estimate the spatiotemporal models with complex specifications (e.g., to address the issues that

^* Corresponding author. E-mail address: kxie@odu.edu (K. Xie).

likelihood functions are not available in closed form), which is computationally expensive. The expedition of the estimation process is still a challenge (Lord and Mannering, 2010).

To address these challenges, we propose a Bayesian spatiotemporal interaction (BSTI) approach for crash frequency modeling at zone-level in Manhattan along with an integrated nested Laplace approximation (INLA) to greatly fasten the estimation process. The proposed BSTI model utilizes a hierarchical framework to jointly account for spatial random effect, temporal random effect and STI effect simultaneously. The BSTI model is estimated using an analytical approach, namely, the INLA approach, as a substitute of the conventional MCMC method to greatly expedite the estimation process. The remaining of the paper is organized as follows. Section 2 presents a literature review regarding state-of-the-art approaches in addressing spatial and/or temporal correlations among crash data. Section 3 presents a detailed description of the dataset utilized. Model development is thoroughly illustrated in Section 4. Section 5 presents the model results and relevant discussions. Finally, Section 6 concludes the paper.

2. Literature review

In this section, we present a review on the methodology used in addressing spatial correlation, temporal correlation, spatial and temporal correlations, and interactions between spatial and temporal correlations when modeling crash data. The strengths and weaknesses of each approach are reviewed. Besides, the efficiency issue in estimating Bayesian models is also discussed.

2.1. Spatial correlation

Spatial correlation in crash frequency modeling refers to that observed crash frequency in one location is correlated with those observed in the adjacent locations (Ziakopoulos and Yannis, 2020). To address spatial correlation, extended versions of GLM, for example the generalized estimating equations (GEE) and the multinomial generalized Poisson with error component (EMGP), have been proposed mostly in the early literature (Lord and Persaud, 2000; Abdel-Aty and Wang, 2006; Chiou et al., 2014). However, these methods cannot specify different spatial correlation matrix to address the variation between different groups (Chiou et al., 2014; Xie et al., 2014).

In the recent literature, two different workhorse characterization methods arise to account for spatial correlations, namely the spatial lag and the spatial error specifications (Narayanamoorthy et al., 2013). The simultaneous autoregressive (SAR) models (Whittle, 1954) and conditional autoregressive (CAR) models (Besag, 1974) are widely used to implement these two characterizations, respectively (Xie et al., 2014; Ziakopoulos and Yannis, 2020). SAR models assume the spatial correlation come from spatial error correlation effects (unobserved exogenous variables of one location affecting the dependent variable of its own and neighboring locations) and/or spatial lag effects (observed exogenous variables of one location affecting the dependent variable of its own and neighboring locations) (Narayanamoorthy et al., 2013; Cai et al., 2016; Ziakopoulos and Yannis, 2020). Some researchers have utilized the SAR to model crash frequency (Quddus, 2008; Narayanamoorthy et al., 2013). Compared to the SAR models, CAR models assume the spatial correlation comes from only the spatial error correlation effects (Narayanamoorthy et al., 2013). The SAR models are generally estimated using the maximum likelihood estimation method (Narayanamoorthy et al., 2013; Xie et al., 2015) and it is challenging to use the SAR models in a count-response setting (Wang and Kockelman, 2013). The CAR models adopt a Bayesian framework, which provides flexibility to accommodate count response variables (Xie et al., 2019). Therefore, the CAR models are more frequently used in crash frequency modeling (Ziakopoulos and Yannis, 2020). Quddus (2008) compared the two specifications of the SAR model, namely, spatial error specification and spatial lag specification with the CAR model and found that the CAR model was more appropriate in developing a relationship between area-wide traffic crashes and the contributing factors.

Besag et al. (1991) firstly implement a Bayesian model to account for spatial autocorrelation with a CAR prior. This implementation has been widely recognized and followed in the later literature regarding road safety analysis (El-Basyouny and Sayed, 2009; Xie et al., 2014; Cai et al., 2018; Zeng et al., 2020a, b). In addition, CAR models are usually developed under full Bayesian framework due to the above-mentioned advantage of addressing spatial correlation among crash count observations when compared to empirical Bayesian framework (Persaud et al., 2010; Xie et al., 2019). Besides, CAR models provide more flexibility in accommodating complex model structure and choosing crash frequency distributions (Lan et al., 2009; Xie et al., 2019). Although the abovementioned research has contributed greatly to the literature in the aspect of addressing spatial correlation among crash data, they have not accounted for temporal correlation.

2.2. Temporal correlation

It is commonly understood that researchers would collect crash data in multiple years to get sufficient crash observations for statistical modeling. To account for the temporal correlation of crash observations, GEE models (Wang et al., 2006) and linear temporal trend models (Andrey and Yagar, 1993) are first developed in the early literature. However, these models either lack flexibility on specifying the correlation matrix or lack of utilizing the information of explanatory variables (Xie et al., 2014). As an alternative, the autoregressive (AR) models can be considered for temporal correlation. Unlike CAR models previously discussed, the AR models usually specify first order or second order random walk priors on structured temporal correlation term (Ma et al., 2017; Liu and Sharma, 2018). It is worth noting that these AR models are mostly implemented under full Bayesian framework for better model flexibility and the incorporation of hierarchy. Furthermore, random-effects (RE) models (Zeng et al., 2018) and random-parameters (RP) models (Xiong et al., 2014) can be considered for addressing temporal correlation (Mannering et al., 2016; Mannering, 2018). RE models assume a common unobserved effect exist over the temporal units/groups with a predefined distribution and the unobserved effects are assumed to be uncorrelated with explanatory variables (Lord and Mannering, 2010). For crash frequency modeling, the RE models recomposed the Poisson parameter λ to address the temporal effects by incorporating an extra term (Zeng et al., 2018). As an extension to RE models, RP models allow a variable temporal effect across different temporal units/groups (Xiong et al., 2014). The abovementioned models are widely utilized in the domain of crash modeling. However, accounting for temporal correlations merely is not enough in modeling crash frequency observations (Chiou and Fu, 2015; Truong et al., 2016).

2.3. Spatial and temporal correlations

Accounting for spatial or temporal correlation solely would lead to biased inference when modeling crash frequencies (Lord and Mannering, 2010; Liu and Sharma, 2018; Zeng et al., 2019). To jointly account for the spatial and temporal correlations in crash frequency modeling, RE models (Truong et al., 2016; Ma et al., 2017) and RP models (Truong et al., 2016) can still be considered (Lord and Mannering, 2010). Different from addressing temporal/spatial correlation in RE models separately, the RE term introduced in a spatiotemporal RE model jointly address the spatial and temporal effect with a single predefined distribution, i.e. unstructured spatiotemporal correlations are assumed. However, a number of studies have revealed the two drawbacks of this specification: 1) accounting for spatial and temporal correlations with one RE term would neglect the possibility that the spatial and temporal correlations may have different structures and this would lead to a biased inference; 2) imposing unstructured effects to both spatial and temporal correlations would neglect the potential structural information

of the correlations and thus lead to biased inference (Mannering et al., 2016; Truong et al., 2016; Liu and Sharma, 2017; 2018). To tackle this problem, some researchers adopted Bayesian hierarchical models for flexible specifications of the spatiotemporal effect structures (Liu and Sharma, 2017; Cheng et al., 2018b; Liu and Sharma, 2018; Wen et al., 2019). For instance, Truong et al. (2016) compared the random-effect negative binomial (RENB) model with the random-effect Bayesian spatiotemporal model with CAR priors. The CAR term in the proposed model is considered as random effects. Results showed that the Bayesian model with CAR priors outperform the RENB model.

In some of the abovementioned research, the spatial and temporal effects are specified as structured effects and unstructured effects, respectively, and are incorporated in the Bayesian hierarchical model simultaneously (Liu and Sharma, 2018). For instance, Ma et al. (2017) proposed a series of space-time models under the Bayesian hierarchical framework with different assumptions on the spatial and temporal random effects, results showed that the spatiotemporal model with structured spatial effects (with CAR priors) and unstructured temporal effects outperformed the other models. Liu and Sharma (2017) compared Bayesian spatiotemporal models with other models including Bayesian spatiotemporal models omitting spatial or temporal effects, Poisson models and zero-inflated Poisson models. Results showed that the Bayesian spatiotemporal models outperformed the others. Some research utilizes multivariate conditional autoregressive (MVCAR) models to address the spatiotemporal correlations. Yang et al. (2019) utilizes MVCAR model to establish time-dependent safety performance functions. Crash frequencies from different time of day are considered as multivariate response variables of the model. A CAR prior is used to specify the spatial correlations. Results showed the superiority of the proposed model. Cui et al. (2019) extended the use of MVCAR model into the analysis of bus operation speed and obtained useful insights.

Despite the superior performance in accounting for spatiotemporal correlations in crash data, these abovementioned Bayesian spatiotemporal models have not considered the potential interaction effects between spatial and temporal trends, i.e., these models implicitly assume that the spatial and/or temporal trends are constant across all spatial/ temporal units.

2.4. Spatiotemporal interactions

There is a severe dearth of research that considered STI in the literature. Dong et al. (2016) modeled crash frequency using spatiotemporal model under Bayesian framework and considers the difference between the mean of each area and the overall mean of temporal trend as the interaction. They assumed varying temporal trend across different zones, and results proved that the proposed model outperformed the others. Cheng et al. (2018b) also modeled crash frequency using the spatiotemporal model considering interaction. They specified an unstructured interaction term to capture: 1) the overall spatiotemporal trend; and 2) the deviation from overall trend. Similar specifications can be found in Cheng et al. (2018a). The interaction term specified in Cheng et al. (2018b) and Dong et al. (2016) neglected the possible STI structures, such as: (1) correlated temporal trends in each location with no spatial correlations; (2) correlated spatial trends in each time unit with no temporal correlations; (3) correlated spatial and temporal trends at each spatial/temporal unit, i.e. the temporal correlation structure for each zone is correlated with the temporal pattern of the neighboring zones as well (Blangiardo and Cameletti, 2015). Despite their contributions, the models they proposed did not considered the essence of STI from a more general way.

In addition, most of the abovementioned research adopt the MCMC method to handle the estimation of the complex functional model forms (Ma and Kockelman, 2006; Xie et al., 2013; Zeng et al., 2018; Zhang et al., 2018; Wen et al., 2019; Zeng et al., 2019). When using an MCMC method for model estimation, a large amount of burn-in iterations is needed, and this process can be very time-expensive. The simulation time consumed by this method is still a barrier to complex model forms (Lord and Mannering, 2010). All the abovementioned studies that address the STI effects have used the MCMC method, which would be very time-consuming. Li et al. (2019) modeled alcohol/drug impaired-driving crash frequency using a Bayesian STI model. They used the MCMC method for model estimation, and a numerous of iterations is needed for model estimation.

With the model estimation efficiency being elevated, researchers can develop more complex model forms for better accounting unobserved heterogeneities, and the model calibration would also be more efficient. Besides, the expedited estimation process would widen the potential for online prediction of crash frequency. To tackle this problem, we utilize an INLA method (Blangiardo et al., 2013; Liu and Sharma, 2017) to greatly expedite the estimation process. Compared to Li et al. (2019), we explore the proposed BSTI model in multi aspects: (1) we firstly adopt the INLA method to greatly expedite the estimation of the BSTI model; (2) we evaluate the performance of the proposed model in addressing spatiotemporal correlations; and (3) we further examine the predictive performance of the proposed model.

In summary, existing research have provided insightful and comprehensive knowledge on accounting spatial and/or temporal correlations when modeling crash frequency. However, most of the studies that jointly account for spatiotemporal heterogeneity assume an overall spatiotemporal trend across all spatial/temporal units without considering the possibility of STIs. For the very limited number of studies that addressed the STIs, most of them assume a single interaction structure and does not examine the other possible structures thoroughly. Besides, the estimation efficiency is still an open issue for crash frequency models, especially for complex model specifications. To fill these research gaps, we propose a BSTI approach for crash frequency modeling. We contribute to the literature in the following aspects: (1) we explicitly account for various STI patterns; (2) we propose a BSTI model associated with a swift estimation method based on the INLA approach to model zonal crash frequency data; (3) we examine the potential predictive performance of the proposed model using test data; and (4) we study the effect of the STI terms by examining the residuals of the proposed models.

3. Data

In this study, the dataset is generated from the raw dataset collected in Manhattan during the years of 2013–2019 sourced by the New York State Department of Transportation. This raw dataset consists factors of three different aspects, namely transportation-related factors such as logged vehicle miles traveled (LogVMT), truck ratio and etc., land use factors for example the commercial ratio and etc., and demo-economic factors for instance population over 65 years old, unemployment rate and etc. For more information regarding this raw dataset, please refer to (Xie et al., 2019).

It is worth noting that the raw dataset is obtained using census tracts as spatial analyzing units. However, as found out in (Xie et al., 2019), the edges of the census tracts overlapped with major arterials and crashes occurring on those arterials are arbitrarily assigned to adjacent census tracts and this would lead to biased inferences. For instance, the whole central park zone is demarcated as a large square-shaped census tract of size 4000 m × 1000 m, and its two long edges overlap respectively with the 5th Avenue and the Central Park West Avenue, which are two major arterials in Manhattan with numerous crash observations. This demarcation is obviously problematic since it is not realistic for a park zone to have so many crash observations. To deal with this drawback, a natural idea is to use identical spatial units. It is found that if one wants to completely cover a plain using identical spatial units, hexagons provide the least total length of perimeter (Ball, 2016). For crashes located on the boundary of a zone, it is challenging to determine exactly which zone should those crashes be assigned to. To alleviate this issue, the hexagon that provides the least total length of boundaries was

Fig. 1. Hexagon-level data from 2013 to 2018 in Manhattan: (a) crash frequencies and (b) VMT.

selected as the spatial unit of analysis.

We use 228 hexagons covering Manhattan as spatial analyzing units for data preparation and crash data modeling, and each hexagon has an edge length of 351.9 m and an area of 0.322 km² . Therefore, an aggregation process is conducted to transform census tract-based data into hexagon-based data. During the aggregation process, the distributable factor (factor that can be assumed to be uniformly distributed over the zone, i.e., the factor value of a subzone is proportional to the size of the subzone, such as population) is summed according to the size of each subzone in the original census tract, while the non-distributable factor (factor that is assumed to be fix over the zone and any subzone has the same value as the census tract, such as median income) is averaged according to the size of each subzone in the hexagon. After the aggregation, the irregular and extremely small hexagons that located on the border of Manhattan with no LogVMT value are deleted for robust model inference. A total of 18 hexagons are deleted, left with 210 hexagons for safety modeling.

The hexagon-level data of the crash frequencies and VMT are presented in Fig. 1. It is found that for crash frequencies, the zones with similar colors tend to be clustered, showing the spatial correlations of the crash frequencies. Besides, the zones with high crash frequencies are likely to present high crash frequencies in the following years, which shows the potential temporal correlations of the crash frequencies. The hexagons on the boundaries of Manhattan present large VMT than the inside. This pattern is consistent with the reality that an expressway exist along the boundaries of Manhattan. It is also found that the zones with higher VMT also tend to have higher crash frequencies. This finding suggest that VMT has significant impact on crash frequencies.

The description and descriptive statistics of each factor are listed in Table 1.

4. Methods

To explore the characteristics of the STI and evaluate their effect in crash frequency modeling, 10 models under a 3-level Bayesian hierarchical framework were established in a progressive order from the Poisson model that does not address any spatiotemporal effect to the Type IV model that addresses the full STI effect for comparison. STI effect and model specification was elaborated (Section II.A and II.B). Secondly, the deduction of the proposed INLA method was demonstrated highlighting its advantage in saving computational cost against MCMC method while maintain the ability of statistical interpretation.

4.1. Models without addressing STI

Let *Y^{t i}*denote the crash counts at i th hexagon in t th year (i = 1, 2, …, n, n is the total number of hexagons, n = 210; t = 1, 2, …, T, T is the total number of years, T = 6). It is commonly assumed that *Y^{t i}*follows Poisson distribution with *λ^{t i}*as the mean. Therefore, on the first level, it is assumed that:

$$P(Y_i^t = y_i^t \middle| \lambda_i^t) = \frac{e^{-\lambda_i^t} (\lambda_i^t)^{y_i^t}}{y_i^{t!}}$$ (1)

where *λ^t i*is further decomposed into a combination of different components on the second level of the hierarchical framework.

4.1.1. Model 1: Poisson model (P model)

By decomposing *λ^t i*into merely intercept and covariates (fixed effects) without spatial or temporal effects, the baseline model is formulated as follows:

Table 1 Description and descriptive statistics of data (sample size: 210 hexagons).

Variable	Description	Mean	Standard
			Deviation
Crash	Count of crashes	203.18	172.82
Transportation
LogVMT	The logarithm of annual average daily vehicle miles traveled	4.54	0.58
	(vehm• ile)
Truck ratio	The average ratio of truck flow	0.05	0.03
	to total flow
Freeway length	Total length of freeway	0.31	0.67
	segments (km)
Avenue length	Total length of avenues (km)	0.87	0.91
Street length	Total length of streets (km)	1.98	3.32
Road density	Ratio of total road length to	11.34	4.77
	hexagon zone (km/km2 )
Intersection	Count of intersections in each	76.74	59.28
number	hexagon
Subway station	Count of subway stations in each	0.72	1.11
number	hexagon
Bus stop number	Count of bus stops in each hexagon	10.70	7.63
Land use
Commercial ratio	The ratio of commercial zone	0.25	0.30
	area to the whole area
Residential ratio	The ratio of residential zone	0.47	0.33
	area to the whole area
Mixed ratio	The ratio of mixed zone area to	0.10	0.21
	the whole area
Park ratio	The ratio of park area to the whole area	0.18	0.27
Demo-economic
Population	Total population (103 )	7.55	5.02
Population under	Population under 14 years (103 )	0.94	0.69
14
Population over	Population 65 years and over	1.00	0.73
65	(103 )
Median age	Median age of population	37.55	5.46
Median income	Median income per household (103 \$)	74.12	36.80
Household size	Number of people per household	2.18	0.38
Unemployment	Share of the labor force that is	0.09	0.04
rate	unemployed
Drive alone ratio	The ratio of commuters by driving alone	0.07	0.03
Carpool ratio	The ratio of commuters by carpooling	0.02	0.02
Public transit ratio	The ratio of commuters by public transit	0.59	0.10
Walk ratio	The ratio of commuters by walking	0.21	0.09

$$\log(\lambdai^t) = \beta^0 + \sum{p=1}^p \beta^p X_{pi}^t$$ (2)

where *β*0 is the intercept term, the X^t pi denotes the pth observed covariate for hexagon i in year t, P is the total number of explanatory variables, β^p is the regression coefficients of the pth covariate. On the third level, priors are assigned to each component of the decomposition of λ^{t i}, due to the lack of prior information, *β*0 and *β^p*are assumed to follow the Gaussian distribution (0105 ).

As previously mentioned, temporal correlations can exist among crash data. This kind of temporal correlations can be addressed by structured temporal effects, while for some temporal heterogeneity that only affects its local time period, unstructured temporal effects needs to be considered. Similarly, structured spatial effects can be used to address spatial correlations, while unstructured spatial effects can be used to account for local spatial disturbance. Therefore, we introduce four random effect terms into Eq. (2) to account for structured and unstructured spatial/temporal effects, respectively (Knorr-Held, 2000). We use u and v to denote structured and unstructured spatial effects, respectively, and use γ and ζ to denote structured and unstructured temporal effects, respectively.

In order to evaluate the effects of each of the four abovementioned random effects, we introduce them progressively as follows.

4.1.2. Model 2: Poisson model with structured temporal effect (PTstrc Model)

By introducing the structured temporal effect term γ into Eq. (2), the PTstrc model can be formulated as:

$$\log(\lambdai^t) = \beta^0 + \sum{p=1}^P \beta^p X_{pi}^t + \gamma_t \tag{3}$$

where *γt*denotes the structured temporal effect in year t. Multiple choices can be specified for modeling γt. For identifiability, we followed Clayton's rules, a first order random walk (RW1) model is specified to represent temporal correlation (Clayton, 1996), which assumes that the parameter of each year depends on the previous one, the RW1 is can be given as follow:

$$\gamma' | \gamma'^{-1} \sim Normal(\gamma'^{-1}, \sigma_{\gamma}^2)$$ (4)

where σ² *^γ*is the global variance of γ.

4.1.3. Model 3: Poisson model with structured and unstructured temporal effects (PT Model)

By introducing the unstructured temporal effect term ζ into in Eq. (3), the PT model can be formulated as:

$$\log(\lambdai^t) = \beta^0 + \sum{p=1}^P \beta^p X_{pi}^t + \gamma_t + \zeta_t$$ (5)

where *ζt*denotes the unstructured temporal effect in year t and is assumed to be i.i.d. among different years, an exchangeable normal prior with 0 mean and precision τζ is assigned to ζt:

$$\zetai \sim Normal(0, \frac{1}{\tau{\zeta}})$$ (6)

4.1.4. Model 4: Poisson model with structured spatial effect (PSstrc Model) By introducing the structured spatial effect term u into Eq. (2), the PSstrc model can be formulated as:

$$\log(\lambda{i}^{t}) = \beta^{0} + \sum{p=1}^{P} \beta^{p} X{pi}^{t} + u{i}$$ (7)

where ui denotes the structured spatial effect at zone i.

Prior specifications proposed by Besag et al. (1991) are commonly adopted in recent studies for its easiness to interpret and simplicity (Knorr-Held, 2000; Blangiardo et al., 2013; Zhang et al., 2018). Following Besag et al. (1991), the intrinsic conditional autoregressive (ICAR) prior was chosen for ui . In the Bayesian framework, the full conditional distribution of ui given u⁻ *ⁱ*is assumed to follow a normal distribution:

$$u^{i} | u^{-i} \sim Normal\left(\frac{1}{N{i}} \sum{l \neq i} u^{l}, \frac{\sigma{u}^{2}}{N{i}}\right)$$ (8)

where u⁻ *ⁱ*is the set of uj where j ∕= i, *Nⁱ*is the number of neighboring zones that are contiguous with the zone i (i.e., sharing boundaries). Noncontiguous zones are regarded as independent in this study. Such definition of zone dependency is widely adopted in the literature (Quddus, 2008; Ukkusuri et al., 2011; Xie et al., 2019; Cheng et al., 2020).

In practice, a zone with multiple neighbors will have a more precise estimate of ui than those rather isolated, thus for the i th zone, the vari-

Fig. 2. Graphical demonstration for different STI structures (circles represent independence, rounded rectangles represent dependence).

ance of ui is a global variance σ² *^u*divided by the number of its neighbors.

4.1.5. Model 5: Poisson model with structured and unstructured spatial effects (PS Model)

By introducing the unstructured spatial effect term v into Eq. (7), the PS model can be formulated as:

$$\log(\lambdai^t) = \beta^0 + \sum{p=1}^P \beta^p X_{pi}^t + u_i + v_i$$ (9)

where vi denotes the unstructured spatial effect at zone i and is considered independent and identically distributed (i.i.d.) among different locations, an normal prior with 0 mean and precision *τv*is assigned to vt:

$$v_t \sim Normal(0, \frac{1}{\tau_v})$$ (10)

4.1.6. Model 6: Spatiotemporal model without interaction (ST Model) To jointly account for spatial and temporal effects, we introduce u, v, γ and ζ into Eq. (2), the ST model can be formulated as:

$$\log(\lambda{i}^{t}) = \beta^{0} + \sum{p=1}^{P} \beta^{p} X{pi}^{t} + (u{i} + v{i}) + (\gamma{t} + \zeta_{t})$$ (11)

All the notations are introduced above. The variance and the precision of both spatial and temporal priors (σ^{2 u} , τv, τζ ) are assigned with a gamma(0.5, 0.0005).

4.2. Models addressing STI

Based on the independent space- and time-effects introduced in Eq. (11) space-time interaction can be easily deducted as a combination of one of the spatial effects and one of the temporal effects (Knorr-Held, 2000). Therefore, four types of interaction exist. This interaction explains the relations between the time trends in different locations. Based on model 6, a STI effect term δit is further embedded into the model as follows:

$$\log(\lambda{i}^{t}) = \beta{0}^{t} + \sum{p=1}^{P{t}} \beta{p}^{t} X{pi}^{t} + (u{i} + v{i}) + (\gamma{t} + \zeta{t}) + \delta_{it}$$ (12)

where δit is the interaction effect for zone i in year t. For computational convenience, the spatial/temporal components ui, vi, γt, *ζt*and δit are assumed to be multivariate Gaussian distributed with zero mean and precision matrix ηK, where η is a scalar to be estimated and K is a known structured matrix different for each component to accommodate different assumptions about the prior interrelationship within the component (Clayton, 1996; Knorr-Held, 2000).

K has non-diagonal elements equal to -1 for neighboring hexagons or years and diagonal elements equal to the number of adjacent hexagons or years. All other elements in K are zero (Knorr-Held, 2000). Following this rule, for example, 4 hexagons numbering from 1 to 4 are presented and only hexagon 1 and 4 are neighboring hexagons, the spatial adjacency matrix of this example K^u eg is as follows:

$$K^{u-eg} = \begin{bmatrix} 1 & 0 & 0 & -1 \ 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 \ -1 & 0 & 0 & 1 \end{bmatrix}$$ (13)

For structured temporal random effect *γt*with a first order random walk prior, the structured matrix of the above example K^γ eg is given below:

$$K_{\gamma-eg} = \begin{bmatrix} 1 & -1 & 0 & 0 \ -1 & 2 & -1 & 0 \ 0 & -1 & 2 & -1 \ 0 & 0 & -1 & 1 \end{bmatrix}$$ (14)

For unstructured spatial and temporal random effect vi and ζt, due to the i.i.d. assumptions, their corresponding structured matrices are identity matrices I.

Given that there are two components in both the spatial and temporal effect, it is straightforward to come up with four different types of interactions (Knorr-Held, 2000). The graphical demonstrations of these four types of interaction are shown in Fig. 2. Clayton (1996) suggests that *Kδ*be specified as the Kronecker product (denoted by " ⊗ ") of the structure matrices of the spatial and temporal effects that are assumed to interact. More details regarding this specification can be found in Clayton (1996).

4.2.1. Model 7: spatiotemporal model with type I interaction (type I model)

Type I interaction assumes that the two unstructured random effects vi and *ζt*interact (as shown in Fig. 2(a)). According to Clayton (1996), the structure matrix for i.i.d. random effects is identity matrix I, and the structured matrix K^I *^δ*is given below:

$$K{\delta}^{\mathrm{I}} = K{\zeta}^{\mathrm{I}} \otimes K_{\nu}^{\mathrm{I}} = I \otimes I = I \tag{15}$$

where the non-italic "I" denotes type I interaction. The prior distribution of δ under Gaussian family is given below:

$$p(\delta | \eta{\delta}) \propto \exp\left(-\frac{\eta{\delta}}{2} \sum{i=1}^{n} \sum{t=1}^{T} (\delta_{it})^{2}\right)$$ (16)

Type I interaction accounts merely unobserved covariates for each data in space and time, yet they do not present any spatial or temporal structure. This is also the difference between Type I Model and ST Model.

4.2.2. Model 8: spatiotemporal model with type II interaction (type II model)

Type II interaction assumes that the structured temporal random effect *γt*interacts with the unstructured spatial random effect vi (as shown in Fig. 2(b)). Since *γt*follows a random walk distribution, δi follows a random walk distribution:

$$\deltat | \delta{t-1} \sim Normal(\delta^{t-1}, \sigma_{\delta}^2)$$ (17)

The structured matrix *K*II *^δ*is given as:

$$K{\delta}^{\mathrm{II}} = K{\gamma}^{\mathrm{II}} \otimes K{\nu}^{\mathrm{II}} = K{\gamma}^{\mathrm{II}} \otimes I = K_{\gamma}^{\mathrm{II}}$$

$$\tag{18}$$

The prior of δ can be written as

$$p(\delta \middle| \eta{\delta}) \propto \exp \left( -\frac{\eta{\delta}}{2} \sum{i=1}^{n} \sum{t=1}^{T} \left( \delta{it} - \delta{i,t-1} \right)^{2} \right)$$ (19)

Type II interaction model indicates that each zone has a different temporal trend, and the time trends do not present any spatial structure.

4.2.3. Model 9: spatiotemporal model with type III interaction (type III model)

Type III interaction assumes that the unstructured temporal random effect *ζt*interacts with the structured spatial random effect ui (as shown in Fig. 2(c)). Since ui follows an intrinsic autoregression, *δi*follows an intrinsic autoregression and is defined as follows:

$$\delta^{i} \left| \delta^{-i} \sim Normal \left( \frac{1}{N{i}} \sum{j \in i, \gamma i \neq j} \delta^{i}{,i}, \frac{\sigma{\delta}^{2}}{N_{i}} \right) \right. \tag{20}$$

where σ² *^δ*is a global variance for δj, i ∼ j denotes j is a neighbor of hexagon i. The structured matrix *K*III *^δ*is given below:

Table 2 Summary of STI structures.

Model	Interaction term	Description
Model 7 (Type I)	δit ∼ N(0, σ2 it)	Interactions have no spatial or temporal structure and can be treated as unobserved covariates.
Model 8 (Type II)	δi,1:t ∼ RW(γt, σ2 t )	Interactions have temporal trends at each hexagon different from neighbors.
Model 9 (Type III)	δ1:i,t ∼ ICAR(ui, σ2 δ )	Interactions have different spatial trends at each year without temporal structure.
Model 10 (Type IV)	δ1:n,1:T ∼ ICAR(ui, σ2 δ ) ⊗ RW(γt, σ2 t )	Interactions have temporal-correlated spatial trends at each hexagon.

$$K{\delta}^{\mathrm{III}} = K{\zeta}^{\mathrm{III}} \otimes K{u}^{\mathrm{III}} = I \otimes K{u}^{\mathrm{III}} = K_{u}^{\mathrm{III}}$$ (21)

The prior of δ can be written as follows:

$$p(\delta | \eta{\delta}) \propto \exp\left(-\frac{\eta{\delta}}{2} \sum{i=1}^{T} \sum{(i,j) \in i \sim j} (\delta{ii} - \delta{ji})^{2}\right)$$ (22)

Type III model indicates that different spatial structure exists for each year, but spatial structures is independent across different years.

4.2.4. Model 10: spatiotemporal model with type IV interaction (type IV model)

The interaction between the structured spatial random effect ui and the structured temporal random effect *γt*would form the Type IV model (as shown in Fig. 2(d)). Since ui follows an intrinsic autoregression and γ^t follows a random walk distribution, δit is completely dependent over space and time and therefore cannot be parameterized as either of these two structured effects. δit can be considered now as a partial combination of Type II and III model. The structured matrix *K*IV *^δ*is given below:

$$K{\delta}^{\rm IV} = K{u}^{\rm IV} \otimes K_{\gamma}^{\rm IV} \tag{23}$$

The prior of δ can be written as follows:

$$p(\delta \middle| \eta{\delta}) \propto \exp \left( -\frac{\eta{\delta}}{2} \sum{t=2}^{T} \sum{(i,j) \in l \sim j} \left( \left( \delta{it} - \delta{jt} \right) - \left( \delta{i,t-1} - \delta{j,t-1} \right) \right)^{2} \right)$$

$$= \exp \left( -\frac{\eta{\delta}}{2} \sum{t=2}^{T} \sum{(i,j) \in l \sim j} \left( \delta{it} - \delta{jt} - \delta{i,t-1} - \delta_{j,t-1} \right)^{2} \right)$$ (24)

Note that Type IV interaction prior is a Markov random field where both spatial and temporal term enter in the full conditional form for δit (Knorr-Held, 2000). This prior assumes that spatial structure exists among zones for a certain year, and this spatial structure are also dependent across different years. Equivalently, one could emphasize that temporal structure exists among years for a certain zone, and this temporal structure are also dependent across neighboring zones. A summary of these four types of interactions is presented in Table 2.

4.3. Model assessment

The deviance information criterion (DIC) is widely used as a Bayesian measure of model performance and complexity (Spiegelhalter et al., 2002, 2005). Specifically, DIC is defined as follows:

$$DIC = \overline{D}(\theta) + p_D \tag{25}$$

where D(θ) is the posterior mean of the deviance of θ, θ is the vector of parameters required for the likelihood function (regarding the abovementioned models, θ contains β⁰, β^p, ui, vi, γt, ζt, δit), D(θ) can be considered as a measure of model goodness-of-fit, pD is the posterior mean of the deviance minus the deviance of the posterior means and indicates the complexity of models (Spiegelhalter et al., 2002).

Besides DIC, mean absolute error (MAE), mean absolute percentage

error (MAPE), and mean squared error (MSE) were used as additional model assessment measures. MAE, MAPE and MSE are expressed as follows:

$$MAE = \frac{1}{n} \sum_{i=1}^{n} |(y_i^{est} - y_i^{obs})|$$ (26)

$$MAPE = \frac{1}{n} \sum_{i=1}^{n} \frac{\left| \left( y_i^{ext} - y_i^{obs} \right) \right|}{y_i^{obs}}$$ (27)

$$MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i^{ext} - y_i^{obs})^2$$ (28)

where *yest ⁱ*is the estimated crash frequency for the i th hexagon, *yobs ⁱ*is the observed crash frequency for the i th hexagon; n is the total number of zones. Models associated with lower MAE, MAPE and MSE fit better to the data.

4.4. Bayesian approach using INLA approximation

4.4.1. Estimations in Bayesian framework

Bayesian method combines prior distributions with a likelihood function obtained to create posterior distributions as estimates. Comparing to the classical frequentist approach, which only returns a fixed unknown quantity of the estimator, Bayesian approach returns a probability distribution, where the quantity returned in the classical frequentist approach being the mean of the distribution (Blangiardo and Cameletti, 2015). This feature benefits us much more information on the parameter we are interested in, and detailed statistical inference could be generated after. Therefore, all the aforementioned models are estimated in the full Bayesian framework. The theoretical framework for Bayesian inference can be expressed as follows:

$$P(\mathbf{\theta}|\mathbf{y}) = \frac{p(\mathbf{\theta}) \cdot L(\mathbf{y}|\mathbf{\theta})}{p(\mathbf{y})}$$ (29)

where y is the vector of observed data, which is crash counts in this study; θ is the vector of parameters required for the likelihood function; P(θ|y) is the posterior distribution of θ given y; p(θ) is the prior distribution of θ; L(y|θ) is the likelihood function; p(y) is the marginal distribution of the observed data. Specifically, for i and t, the likelihood function for space-time interaction model is as follows:

$$L{i}^{t}(\mathbf{y}|\mathbf{\theta}) = \frac{1}{\sqrt{2\pi}\sigma^{u}} \exp\left(-\frac{\left(y{i}^{t} - \left(\beta{0}^{t} + \sum{p=1}^{P{t}} \beta{p}^{t} X{pi}^{t} + \left(u{i} + \lambda{i}\right) + \left(\gamma{r} + \zeta{r}\right) + \delta{ii}\right)\right)^{2}}{2\sigma_{ii}^{2}}\right)$$ (30)

4.4.2. MCMC method

In practice, posterior distribution is usually too complicated to be analytically calculated. To overcome this issue, MCMC method is widely used for solving Bayesian models (Gilks et al., 1995; Xie et al., 2019). The key technique of MCMC is Gibbs sampling (Geman and Geman, 1984). A burn-in phase is firstly required for the Markov chain to receive convergence. The posterior marginal distributions for the parameters of interest are then derived. Subsequently, each iteration draws a sample for each parameter of interest using Gibbs sampling for statistical inference (Lunn et al., 2000). The computational burden is tremendous, especially for complicated Bayesian models (Knorr-Held, 2000).

4.4.3. INLA method

As an alternative to MCMC method, the INLA method, which is a deterministic approach, was proposed by Rue et al. (2009) to solve Bayesian models in a tremendously faster way. For Bayesian inference, the posterior marginal distributions of the parameters are needed. However, if some variables are with non-Gaussian distributions, it is rather difficult to derive closed-form posterior marginals. For such models, the MCMC method is extensively utilized to generate the posterior marginal distributions. The MCMC method can be considered as a simulation process using the technique of Gibbs sampling and the Markov property (Geman and Geman, 1984). A large amount of iterations are needed for the Markov chain to converge, sampling based on the converged posterior marginal distributions can then be conducted for statistical inference. This process can be very time-consuming. To make the Bayesian inference more efficient, the INLA method has been proposed Rue et al. (2009). The INLA method is able to approximate any distributions using Gaussian distributions. The Laplace approximation technique is the core of the INLA method. Under the Bayesian framework, the posterior marginal distributions of interest can be expressed as follow,

$$p(\theta|y) = \int p(\theta, \pi|y) d\pi \tag{31}$$

where π denotes the hyperparameters and other variables are previously defined. Using the INLA method, no matter what the form of ^∫ p(θ, π|y) dπ is, it can be approximated by a Gaussian distribution. As such, a deterministic posterior marginal distribution can be derived.

Here we demonstrate the Laplace approximation in a mathematical manner. Detailed information can be found in Rue et al. (2009). Suppose we have a density function f(x) of a random variable X, and we are interested in the following integration:

$$\int f(x)dx = \int \exp(\log(f(x)))dx \tag{32}$$

where log(f(x)) can be represented by a Taylor series expansion evaluated in x = *x*0, which is given below:

$$\log(f(x)) \approx \log(f(x_0)) + (x - x0) \frac{\partial \log(f(x))}{\partial x} \bigg|{x = x_0} + \frac{(x - x0)^2}{2} \frac{\partial \log(f(x))}{\partial x^2} \bigg|{x = x_0}$$ (33)

by setting x⁰ = x^∗ = argmax*x*log(f(x)), then ∂log(f(x)) ∂x x=x^∗ = 0, and Eq. (33) turns into as follows:

$$\left.\log(f(x)) \approx \log(f(x^)) + \frac{(x-x^)^2}{2} \left.\frac{\partial \log(f(x))}{\partial x^2}\right|_{x=x^*} \tag{34}\right.$$

Thus Eq. (32) becomes as follows:

$$\int f(x)dx \approx \int \exp\left(\log(f(x^)) + \frac{(x - x^)^2}{2} \frac{\partial \log(f(x))}{\partial x^2}\Big|_{x = x^*}\right) dx$$

$$= \exp(\log(f(x^))) \int \exp\left(\frac{(x - x^)^2}{2} \frac{\partial \log(f(x))}{\partial x^2}\Big|_{x = x^*}\right) dx,$$ (35)

where the integration can be represented by the density of a Normal distribution if we set σ^2∗ = − 1*/*∂log(f(x)) ∂x² x=x^∗ , and Eq. (35) becomes as follows:

$$\int f(x)dx \approx \exp(\log(f(x^))) \int \exp\left(-\frac{(x-x^)^2}{2\sigma^{2*}}\right) dx$$ (36)

where inside the integration is a Normal distribution with mean equal to x^∗ and variance equal to σ^2∗, Normal(x^∗,σ²∗). Therefore, the integration can be expressed as follows:

$$\int_{n}^{m} f(x)dx \approx f(x^{})\sqrt{2\pi\sigma^{2}}(F(m) - F(n))$$ (37)

where F(⋅) is the cumulative density function of Normal(x^∗, σ²∗). So far, an arbitrary density function has successfully been

Table 3 Comparison of model performance using the training set.

Model						Estimation time
	DIC	DIC/n1	MAE	MAPE (%)	RMSE	INLA (s)	MCMC (s)
Poisson	46161.02	36.64	58.19	136.01	80.90	4	363
PTstrc	46069.57	36.56	58.06	135.87	80.82	4	416
PT	46059.61	36.56	56.47	134.07	79.36	5	492
PSstrc	12676.23	10.06	17.56	21.25	26.73	5	1023
PS	12676.22	10.06	17.56	21.25	26.73	7	1360
ST	12544.26	9.96	16.91	21.02	25.78	8	1894
Type I	10665.48	8.46	4.17	12.97	5.88	17	2194
Type II	10468.40	8.31	4.76	10.77	7.01	84	4319
Type III	10536.19	8.36	5.30	12.97	7.28	116	5069
Type IV	10415.03	8.27	5.91	10.85	8.59	191	7378

Note: the bold numbers denote the smallest numbers in the corresponding columns. 1 Standardized DIC (Miaou and Song, 2005).

transformed into an approximation of Normal distribution, which saves a tremendous amount of time compared to MCMC method (Blangiardo et al., 2013).

5. Results and discussion

All models introduced in the previous section are implemented in R using the R-INLA package (Rue et al., 2009; Blangiardo and Cameletti, 2015), which utilizes INLA approach for Bayesian model estimation. The experiments are conducted on a platform with an Intel 6-core i7 at 2.6 GHz processor and 16.0 GB RAM. Using the INLA method, the estimation of each of the proposed models takes less than 200 s. The time consumed for the MCMC method are also tested. For the MCMC method, two simulation chains are established with 100,000 iterations to run in each chain. The first 50,000 iterations are discarded as burn-ins. The thinning interval is set to 1. The MCMC method is implemented in WinBUGS (Spiegelhalter et al., 2005).

5.1. Model performance comparison

The model comparison results using multiple goodness of fit measures like DIC, MAE, MAPE and RMSE are presented in Table 3. The DIC value of the PTstrc model is almost the same as that of the PT model and is at least 92 less than the baseline model (Poisson model), showing that accounting for temporal correlations solely is not a suitable choice for modeling crash frequencies. Compared to temporal effects, accounting for spatial correlations among data (PSstrc model) would lead to a significant decrease in the DIC value (from 46,161 to 12,676). This finding proves the existence of a strong spatial correlation among data. By accounting for both spatial and temporal effects, the DIC value of the ST model decreases for another 132 compared to the PS model. This finding proves that addressing spatiotemporal correlations in crash frequency modeling would lead to better results than addressing them separately. Furthermore, when incorporating the STI term into the ST model, four different STI models can be obtained by specifying different type of interactions between the spatial and temporal effects, i.e., the Type I, Type II, Type III and Type IV models. It is found that the DIC values of all the STI models decrease at least 1878, where Type IV model presents the smallest DIC among all the proposed models. A Type IV model assumes that the structured temporal effects interact with the structured temporal effects, indicating that the temporal trend of one location is correlated with that of adjacent locations. The Type II model presents the second smallest DIC values of all the proposed models with 50 in DIC higher than the Type IV model. A Type II model assumes that each location has its own temporal trend independent from adjacent locations. Following Miaou and Song (2005), the standardized DIC is used as another metric for model comparison and is also presented in Table 3. The standardized DIC is denoted as DIC/n, where n = 210 × 6 = 1260 is the sample size of the data set. Based on the DIC/n results, the Type II model and the Type IV model present very similar DIC/n values (with a difference of 0.04). Therefore, more metrics should be used to select the best model.

Regarding other model measurements in Table 3, the Type I model presents the smallest MAE (4.17) and RMSE (5.88), while Type II model presents the smallest MAPE of 10.77%. It is worth noting that the model with the lowest DIC value (Type IV model) does not provide the smallest MAE, MAPE or RMSE. Therefore, further investigation on the performance of model fitting is needed.

Based on the model estimation time presented in Table 3, it is obvious that the INLA method can greatly expedite the estimation process. For the Type IV model, which is the most complex model in this study, the estimation time using INLA is less than 200 s. However, the estimation time needed for the MCMC method is over two hours.

5.2. Residual analysis

To quantitatively evaluate the effect of addressing spatial autocorrelation of variables of interest, the global univariate Moran's I tests on

Table 4 Results of Moran's I tests on crash frequencies and residuals of the proposed models.

	2013		2014		2015		2016		2017		2018
	I	p-value1	I	p-value	I	p-value	I	p-value	I	p-value	I	p-value
Crash	0.57	0.000	0.56	0.000	0.53	0.000	0.53	0.000	0.50	0.000	0.48	0.000
Poisson	0.20	0.000	0.23	0.000	0.20	0.000	0.19	0.000	0.13	0.002	0.18	0.000
PTstrc	0.20	0.000	0.24	0.000	0.20	0.000	0.20	0.000	0.14	0.002	0.18	0.000
PT	0.20	0.000	0.24	0.000	0.20	0.000	0.20	0.000	0.14	0.002	0.18	0.000
PSstrc	0.24	0.000	0.24	0.000	0.10	0.018	0.19	0.000	0.14	0.003	0.36	0.000
PS	0.24	0.000	0.24	0.000	0.10	0.018	0.19	0.000	0.14	0.003	0.36	0.000
ST	0.22	0.000	0.24	0.000	0.11	0.010	0.15	0.003	0.15	0.002	0.25	0.000
Type I	0.05	0.095	0.15	0.002	0.09	0.017	0.12	0.008	0.17	0.000	0.18	0.000
Type II	0.00	0.474	− 0.03	0.304	0.06	0.061	0.09	0.031	0.05	0.091	0.04	0.210
Type III	0.06	0.068	0.09	0.016	0.09	0.018	0.28	0.000	0.13	0.003	0.09	0.019
Type IV	− 0.21	0.000	− 0.14	0.001	− 0.18	0.000	− 0.16	0.000	− 0.12	0.003	− 0.21	0.000

¹ Pseudo p-value, the bold numbers denote insignificant spatial correlations, i.e. pseudo p-value ≥ 0.05.

Table 5 Results of autocorrelation tests on crash frequencies and residuals of the proposed models.

Model	Lag1	Lag2	Lag3	Lag4	Lag5	95% CIa
Crash	− 0.664	0.510	− 0.337	0.039	− 0.049	(− 0.816,
						0.816)
Poisson	− 0.564	0.129	− 0.086	− 0.029	0.049	(− 0.816,
						0.816)
PTstrc	− 0.582	0.140	− 0.083	− 0.013	0.037	(− 0.816,
						0.816)
PT	− 0.585	0.146	− 0.087	− 0.011	0.036	(− 0.816,
						0.816)
PSstrc	− 0.447	0.162	− 0.130	− 0.270	0.186	(− 0.816,
						0.816)
PS	− 0.447	0.162	− 0.130	− 0.270	0.186	(− 0.816,
						0.816)
ST	0.036	− 0.252	− 0.426	− 0.135	0.278	(− 0.816,
						0.816)
Type I	0.006	− 0.230	− 0.466	− 0.033	0.224	(− 0.816,
						0.816)
Type II	0.157	− 0.372	− 0.456	− 0.126	0.297	(− 0.816,
						0.816)
Type III	− 0.265	− 0.514	0.246	0.244	− 0.210	(− 0.816,
						0.816)
Type IV	0.245	− 0.510	− 0.485	0.011	0.239	(− 0.816,
						0.816)

^a CI denotes confidence interval.

crash counts and residuals of the proposed models are conducted. As mentioned previously, the spatial relationships of zones are defined based on adjacency, i.e., the neighbors of a zone are defined as the zones that share borders with itself. We apply the same spatial relationship for the Global Moran's I tests as in the model development. The Moran's I statistic is bounded between -1 and 1, where a value of 0 indicates no spatial autocorrelations. Positive Moran's I values denote positive autocorrelation, which means that sites with high/low crash counts are likely to be surrounded by sites with relatively high/low crash counts. Similarly, negative Moran's I values denote negative autocorrelation, which means that high/low crash counts are likely to be surrounded by sites with relatively low/high crash counts. Pseudo-p-value is used to test the significance of the spatial correlation, where a pseudo- p-value less than 0.05 suggests a significant spatial correlation of observations. For more details on Moran's I test, please refer to Xie et al. (2015) and Xie et al. (2019). The test results are summarized in Table 4.

For the crash frequency, a strong and significant positive spatial correlation pattern is found in each year. This pattern presents a weakening trend as the Moran's I statistic decreases from 0.57 in 2013 to 0.48 in 2018. For brevity, we refer the Poisson model, the PTstrc model, the PT model, the PSstrc model, the PS model and the ST model as noninteraction models. Regarding the residuals of the non-interaction models, we can see that the Moran's I statistics of these models decrease at least half of that of the crash frequency. This indicates that a great proportion of the spatial correlations can be addressed by these models. However, the pseudo p-values of these models are smaller than 0.05, which means that the residuals of these models still present significant spatial correlations. Therefore, these models are not suitable for modeling crash frequencies. The ST model presents smaller Moran's I values than the other non-interaction models, indicating the necessity of addressing both spatial and temporal correlations in crash frequency modeling. However, the pseudo p-values of the ST model are smaller than 0.05, which shows that the ST model is not a good choice for modeling crash frequencies neither. For the STI models, it is interesting to see that the residuals of the Type IV model present significant spatial correlations in each year. This connote that a more complex model does not guarantee a better performance in crash frequency modeling. Compared to the Type IV model, the Type I model and the Type III model present much smaller Moran's I values in each year. However, except for year 2013, the residuals for the latter two models in other years still

Table 6 Results of the multicollinearity test.

Variable	VIF
Transportation
LogVMT	1.270
Bus stop number	2.092
Truck ratio	1.428
Avenue length	2.271
Intersection number	2.045
Land use
Commercial ratio	3.117
Residential ratio	3.418
Mixed ratio	1.496
Demo-economic
Population	2.812
Drive alone ratio	1.291
Unemployment rate	1.239

present significant spatial correlations. This indicates the incapability of these models in addressing spatial correlations among crash frequency data. On the contrast, except for the year 2016, the Type II model presents insignificant and close-to-zero Moran's I values, which shows its superior ability in addressing spatial correlations among crash data.

For the time series of both the crash data and the residuals of the proposed models, the autocorrelation tests regarding their own lagged values are conducted, respectively. The results are summarized in Table 5. Surprisingly, the temporal correlations of the crash counts are not significant. A possible reason is that the temporal data is insufficient. The results presented in Table 5 are obtained by averaging each hexagon, and each hexagon has only 6 time periods. Therefore, for lag1, only 5 data points for each hexagon is utilized and less data points will be utilized as the degree of lag increases. Based on the above reason, we therefore focus on lag1 and lag2 cases for further analysis as they are more reliable than other cases. It can be found that the lag1 autocorrelation coefficient of the Type II model, Type I model and ST model are respectively 0.036, 0.006 and 0.157, which are much smaller than the coefficients of the crash data and other models. And the lag2 autocorrelation coefficients of these three models are also smaller than the coefficient of the crash model. This indicates that these three models can effectively address the temporal correlations among crash data.

Based on the above analysis, the Type II model outperforms other models (including the Type IV model) greatly in addressing spatial correlations and performs well in addressing temporal correlations. In addition, the Type II model presents the smallest MAPE in model estimation. Furthermore, the DIC value of the Type II model is only about 50 larger than the smallest DIC value presented by the Type IV model. Therefore, it is reasonable to choose the Type II model as the best model.

5.3. Variable interpretation

Different combinations of the variables are tested. We first construct a simple linear regression (LR) model without any random effect in the R platform using function regsubsets (Lumley and Lumley, 2013) to identify the best LR-based combination. This best LR-based combination serves as the start point for testing using the Type II model. Subsequently, 15 combinations are tested using the Type II model. Each variable can stay in the Type II model if the variable meets the following two requirements: (i) is statistically significant and (ii) has a correct positive/negative symbol. If a combination has at least one variable that violates all of these requirements, the combination is discarded. Based on the DIC values and the two requirements, the best tested combination is selected for variable interpretation. All the variables in the best tested combination are statistically significant and have correct positive/negative symbols.

The multicollinearity of the selected variables are measured using variance inflation factors (VIF) (Xie et al., 2017). A VIF greater than 5 indicates the existence of a multicollinearity problem (O'brien, 2007).

Table 7 Posterior parameter estimates of the STI models.

Variables	Type I model		Type II model		Type III model		Type IV model
	Mean	95% BCI	Mean	95% BCI	Mean	95% BCI	Mean	95% BCI
Intercept	2.230	(1.426, 3.037)	0.509	(0.106, 0.912)	0.651	(0.275, 1.027)	2.236	(− 34.549, 38.991)
Transportation
LogVMT	0.101	(0.005, 0.198)	0.095	(0.020, 0.171)	0.197	(0.130, 0.263)	0.152	(− 0.012, 0.318)
Bus stop number	0.037	(0.020, 0055)	0.040	(0.013, 0.069)	0.042	(0.028, 0.057)	0.084	(− 2.860, 3.027)
Truck ratio	0.119	(− 4.374, 4.606)	7.133	(0.647, 13.612)	3.968	(0.452, 7.518)	0.000	(− 61.870, 61.818)
Avenue length	0.268	(− 0.137, 0.673)	0.626	(0.025, 1.225)	0.351	(0.024, 0.679)	0.128	(− 46.493, 46.710)
Intersection number	0.002	(0.000, 0.005)	0.004	(0.001, 0.008)	0.003	(0.001, 0.005	0.004	(− 0.375, 0.382)
Land use
Commercial ratio	2.261	(1.522, 3.003)	1.243	(0.349, 2.137)	1.725	(1.190, 2.255)	0.543	(− 48.013, 49.059)
Residential ratio	1.932	(1.322, 2.544)	0.719	(0.066, 1.370)	1.420	(0.996, 1.840)	0.127	(− 45.926, 46.142)
Mixed ratio	1.619	(− 4.374, 4.606)	1.249	(0.095, 2.403)1	1.364	(0.737, 1.989	0.426	(− 50.756, 51.566)
Demo-economic
Population	0.013	(− 0.021, 0.046)	0.047	(0.005, 0.090)	0.024	(− 0.002 0.050)	0.045	(− 5.223, 5.308)
Drive alone ratio	− 1.554	(− 5.984, 2.868)	5.437	(0.268, 10.601)	1.393	(− 1.853, 4.660)	− 0.024	(− 61.880, 61.780)
Unemployment rate	− 3.122	(− 6.074, -0.181)	3.631	(0.471, 6.792)1	− 0.500	(− 2.744, 1.766)	− 0.027	(− 61.403, 61.297)

Note: Bold indicates significance. 1 Significant at 90% confidence level.

The results are presented in Table 6. It can be seen that all the covariates have VIF values smaller than 5, indicating that no multicollinearity problem exists among the covariates.

The Bayesian posterior estimate results for the Type I model, the Type II model, the Type III model and the Type IV model regarding the best tested combination are summarized in Table 7. In this paper, the 95% Bayesian Credible Interval (BCI) is used to test the significance of the covariates. Estimates are regarded as significant at the 95% confidence level if the BCIs do not cover 0 and vice versa. Except for mixed ratio and unemployment rate, which are statistically significant at the 90% confidence level, all other covariates included in the Type II model are statistically significant at the 95% confidence level. The mean values of the significant variables are denoted in bold.

In Table 7, compared to the Type II model, other models have variables that are insignificant and/or with incorrect symbols. For example, all the demo-economic variables in the Type I model, the Type III model and the Type IV model are insignificant. Furthermore, all variables in the Type IV model are insignificant. For the significant variables in the models excluding the Type IV model, almost all the variables present relatively large difference in model estimates as the structure of STI changes. These findings prove the strong impact of different STI on variable estimates, also the importance of fitting the correct STI to the data. On the other hand, only bus stop number presents similar estimates for the Type I model, the Type II model and the Type III model. The above findings also prove the superiority of the Type II model, its Bayesian posterior estimates are used for parameter interpretation.

Regarding the transportation parameters, vehicle miles traveled (VMT) is found to be positively correlated with the crash frequency, indicating that the zones with larger VMT are more likely to have more crashes. This finding is consistent with numerous previous studies (Abdel-Aty et al., 2013; Ma et al., 2017; Liu and Sharma, 2018; Li et al., 2019; Xie et al., 2019). The reason is intuitive that traveling more miles would increase the exposure to crash risk. The coefficients of the logarithmic VMT can be interpreted as: 1% increase in VMT would lead to a 0.095% increase in crash counts. Avenue length shows positive impact on crash counts, which is consistent with multiple previous studies (Ukkusuri et al., 2011; Abdel-Aty et al., 2013; Wang et al., 2013; Wang and Kockelman, 2013; Truong et al., 2016). One-unit increase in the avenue length is expected to raise crash counts by e0.626 – 1 = 87%. Besides, bus stop number is also found to be positively related to crash frequencies, which is consistent with the findings in Ukkusuri et al. (2011) and Lee et al. (2015). A possible reason is that more bus stops would lead to more conflicts between buses and other vehicles during the entrance and exit of bus stops. However, as discovered in Munira et al. (2020), the presence of bus stops at intersections may decrease the chances of incapacitating injury crashes while increase the chances of non-incapacitating injury crashes. As such, the influence of bus stops on crashes can indeed be complex and more efforts are worth devoting to study this part. Truck ratio is found to have positive impacts on crash counts, which is consistent with Xie et al. (2017) and Martin (2002). The explanation is straightforward, the power performance of trucks are much worse than most of other vehicles, leading to a bigger gap in speed variance and thus increase the risk of crashes, and their size are usually larger, it is more likely to block the horizon of drivers in adjacent vehicles, this would also increase the risk of crashes. We also find that the increase of intersection number could raise the possibility of crashes. This finding is consistent with multiple previous studies (Ukkusuri et al., 2011; Abdel-Aty et al., 2013; Truong et al., 2016; Zeng et al., 2019). Due to the merge of different traffic flows, intersections are one of the most hazardous locations in a road network (Guo et al., 2010; Zeng et al., 2019). More intersections would increase the possibility of conflicts between different traffic participants, and therefore may increase the crash risk. Besides, more signal changes may lead to the increase of rear-end crashes at signalized intersections (Wang and Abdel-Aty, 2006).

Land use-related factors are widely studied and recognized as influencing factors for crash frequency in the previous research (Ukkusuri et al., 2011; Wang and Kockelman, 2013; Xie et al., 2017). As presented in Table 7, the ratio of commercial area, the ratio of residential area and the ratio of mixed land use area have positive effect on crash frequencies. A possible reason is that regardless of the land use type, more land use would increase more transportation activities in the zone and therefore increase the risk of crash. Among all three land-use covariates, the ratio of residential area has the smallest posterior coefficient estimate. This finding is consistent with previous works (Pulugurtha and Sambhara, 2011; Xie et al., 2017). A possible reason is that the commercial zones attract more traffic than the residential zones. This would increase the exposure to conflicts and thus increase the crashes. Moreover, the commercial zones have more commercial-related activities than other zones. These commercial activities often involve package (un)loading using large vehicles like trucks, which can block the eyesight of other drivers. Consequently, these commercial-related activities would increase the chance of exposure to conflicts at commercial areas.

Previous studies have recognized the influence of demo-economic factors on crash occurrence including population (Ukkusuri et al., 2011; Wang and Kockelman, 2013; Lee et al., 2015; Xie et al., 2017; Cheng et al., 2018a) and unemployment rate (Liu and Sharma, 2017; Xie et al., 2017; Liu and Sharma, 2018). Population is found to be positively associated with crash counts. This is intuitive because more population implies more developed areas with more transportation activities, this could increase the likelihood of crashes. Another finding is that more

Table 8 Summary of random effects based on Type II.

Effects	Estimated	95% BCI
Spatial structured
μu	− 1.75E-08	(− 0.022, 0.022)
σ2 u	4.43E-05	(5.91E-04, 1.15E-05)
Spatial unstructured
μv	6.64E-06	(− 0.024, 0.024)
σ2 v	4.75E-05	(6.51E-04, 1.26E-05)
Temporal structured
μγ	2.28E-07	(− 0.269, 0.269)
σ2 γ	4.48E-05	(5.25E-04, 1.24E-05)
Temporal unstructured
μξ	1.10E-03	(− 0.031, 0.035)
σ2 ξ	4.64E-04	(2.16E-03, 1.71E-04)
STI
μδ	4.83E-01	(0.034, 0.994)1
σ2 δ	4.20E-03	(5.61E-03, 3.19E-03)

Note: μ denotes the estimated mean and *σ*2 denotes the estimated variance. 1 BCI at 90% confidence level.

commuters that drive alone to work could results in more crashes. As expected, unemployment rate has positive effects on crash counts. A possible explanation could be that unemployment could increase mental stress and result in aggressive driving behavior, which could increase crash risk (Liu and Sharma, 2018).

The abovementioned variable interpretations also provide practical implications for road safety improvement. For instance, it is recommended to develop effective traffic separation scheme for trucks and other vehicles to reduce truck ratios in areas with high traffic demand. In the planning stage, it is suggested to use greater intersection spacings to avoid creating more intersections. Besides, agencies may encourage carpooling by releasing incentives such as dedicated carpooling lanes.

5.4. Random effect interpretation

The estimated mean (μ) and variance (*σ*2) of the spatial, temporal and STI effects of the Type II model and their 95% BCI are summarized in Table 8. The variance of the STI effect accounts for 4.20E-03 / (4.43E− 05 + 4.75E− 05 + 4.48E− 05 + 4.64E− 04 + 4.20E− 03) = 87.5% of the total variance of all random effects and both its mean and variance are statistically significant. This result indicates that strong interactions between spatial effects and temporal effects exist among the study area and this interaction explains 87.5% of the variability of the

Table 9 Comparison of prediction performance using the test set.

Model	MAE	MAPE (%)	RMSE
Poisson	60.49	104.04	84.25
PTstrc	60.68	104.42	84.67
PT	60.71	104.47	84.73
PSstrc	35.56	34.95	54.26
PS	35.56	34.95	54.26
ST	32.92	33.91	49.63
Type I	36.46	34.17	55.7
Type II	28.51	22.51	43.58
Type III	164.7	91.86	209.46
Type IV	30.64	23.18	46.35

model. This proves the dominance of the STI effect and the necessity of addressing this effect in the model. The exponential of the mean of the STI effect is e0.0483 = 1.62 > 1, which indicates that the STI effect has a positive impact on crash frequency.

As shown in Table 8, compared to the STI effects, the spatial and temporal effects seem less "effective" in accounting for variance of the model. For spatial effects, the structured spatial effect (u) accounts for 48.3% of the total spatial effects (u + v), indicating that the structured spatial effects played a slightly less important role than unstructured spatial effect in addressing spatial correlations. For temporal effects, the structured temporal effect (γ) accounts for 8.8% of the total temporal effects (γ + ξ), indicating that the unstructured temporal effects played a more essential role than structured temporal effect in addressing temporal correlations.

The difference of the STI effects between every two consecutive years are presented in Fig. 3. A difference greater than zero (denoted in red) indicates an increase of the STI effects as the time evolves, while a difference smaller than zero (denoted in purple) indicates a decrease in the STI effects over time and zero difference in STI is denoted in white. Under this graphical representation, in general, the difference of the STI effect through time is decreasing as the colors of red and purple become lighter. This indicates the STI effect tends to be stable over time. The reason may be multifold, for example, as the development of the area being more and more saturated, the traffic characteristics in Manhattan become stable and this would lead to the decreasing of the STI effects. Furthermore, no zones are changing randomly between increasing and decreasing over time, this indicates that the STI effects are structured, which is consistent with the findings in Table 8.

Concretely, three major patterns can be recognized in Fig. 3: (1) a change from dark red to light red. For instance, the upper-east side of Manhattan. This indicates that the STI effects in this area is increasing

Fig. 3. Differenced STI effects between every two consecutive years.

Fig. 4. Predictive performance of the Type II model using the test set.

slower and slower through time, which also implies that the incrementing effect of the STI on crash frequency is vanishing. The results show that given this pattern, multizone-level countermeasures may be needed to mitigate the potential increase in crash frequency; (2) a change from purple to white. For instance, the lower-middle part of Manhattan. This shows that the STI effects in this area is decreasing slower and slower through time and tends to remain the same as previous years. Although the decreasing trend of the STI effects under this pattern is fading, the actual effect of the STI on decreasing the crash frequency is still growing. For target zones under this pattern, relevant authorities may focus more on alleviating the crash counts in adjacent zones to decrease the crash frequencies in the target zone; (3) a change from red to purple. For instance, the lower-east part of Manhattan. This pattern shows that the STI effects increases at first and then decreases and it may indicate that the neighboring areas have improved the safety performance and thus lead to the decrease of STI effects in the center area.

5.5. Predictive performance

The crash data in 2019 is considered as the test set in this study to examine the predictive performance of the proposed models, i.e., the information hidden in the crash frequency data in 2019 was not introduced into any of the proposed models. The results of the evaluation metrics are presented in Table 9.

As shown in Table 9, the Type II model presents the best predictive performance with smallest MAPE being 22.51% (the MAE and RMSE are also the smallest). This is consistent with the above analysis regarding the best selected model and further proves the superiority of the Type II model. Fig. 4 presents the zonal-level comparison between the observed crash frequency and predicted crash frequency in 2019. Fig. 4(c) presents the prediction error in a zonal level. In Fig. 4(c), the prediction error is randomly distributed, which again proves the capability of the proposed Type II model in dealing with future data.

Fig. 4(d) shows the scatter plot between the observed crash frequency and predicted crash frequency, and the solid line is the diagonal line. We can see that in overall, the points are very close to the diagonal line, indicating that the predicted crash frequency value is consistent with the observed crash frequency value. Specifically, for crash frequencies under 250, the points are distributed around the diagonal line, while for the crash frequencies over 250, the points tend to be distributed above the diagonal line.

Based on the above findings and the fact that the estimations of the proposed models take only between 4–200 s, the proposed Type II model shows great potential in crash frequency prediction.

6. Conclusions

This paper proposes a BSTI approach for crash frequency modeling at hexagon-level with an INLA method to greatly expedite the estimation process. Crash, transportation, land use and demo-economic data of Manhattan from 2013 to 2019 were collected for this study, where data from 2013 to 2018 are used for model estimation and the data in 2019 are used for testing. Ten different models are developed and compared. Results show that the proposed Type II model outperforms the others by assuming that the structured temporal random effect interacts with the unstructured spatial random effect. This proves the necessity and the superiority of addressing the STIs in crash frequency modeling. By conducting Moran's I test, residuals of the Type II model are found to be randomly distributed over space with insignificant spatial correlation. Results also show that VMT, bus stop number, truck ratio, avenue length, intersection number, ratio of commercial area, ratio of residential area, ratio of mixed land use area, population, ratio of commuters by driving alone and unemployment rate have positive impact on crash occurrence. The adopted INLA has greatly improved the efficiency in model estimation, the estimation time is in the range from 4 s to 200 s.

The predictive performance of ten models is further investigated and compared using the test data. The Type II model yields the lowest MAPE 22.51%, with randomly distributed prediction error over space. The scatter plot of the observed crash counts vs the predicted crash counts also indicates great predictive performance of the Type II model.

This research makes both methodological and empirical contributions. The study tackles the rarely addressed spatiotemporal interaction problem in crash frequency modeling by introducing a novel BSTI model. This model takes the advantage of the INLA method to greatly expedite the estimation process and thus provide the empirical superiority in prediction. The proposed BSTI approach can potentially advance safety analytics by achieving high prediction accuracy and computational efficiency while maintaining its interpretability on the effects of contributing factors and the unobserved heterogeneity.

In this study, the temporal effects of crash observations are insignificant likely due to the data limitation. In future research, more data of a longer study period will be collected to further compare the proposed models, especially for temporal correlations. Besides, crashes may be categorized according to severity to present deeper understanding regarding the STI effects between different crash types. Furthermore, besides first order/second order random walk priors and conditional autoregressive prior, other priors for the random effect terms can be utilized and compared. In addition, it should be noted that modeling results might be subject to the number of analysis units selected. Therefore, it would also be interesting to investigate how the results are sensitive to the number of analysis units. The advantage of INLA in computational efficiency can be utilized for scenarios with large analysis units.

CRediT authorship contribution statement

Haipeng Cui: Methodology, Software, Writing - review & editing, Validation, Data curation. Kun Xie: Conceptualization, Methodology, Software, Writing - review & editing.

Declaration of Competing Interest

The authors report no declarations of interest.

Acknowledgments

The contents of this paper present views of the authors who are responsible for the facts and accuracy of the data presented herein. The contents of the paper do not reflect the official views or policies of the agencies.

References

• Abdel-Aty, M., Wang, X., 2006. Crash estimation at signalized intersections along corridors: analyzing spatial effect and identifying significant factors. Transport. Res. Rec. (1953), 98–111.
• Abdel-Aty, M., Lee, J., Siddiqui, C., Choi, K., 2013. Geographical unit based analysis in the context of transportation safety planning. Transp. Res. Part A Policy Pract. 49, 62–75.
• Andrey, J., Yagar, S., 1993. A temporal analysis of rain-related crash risk. Accid. Anal. Prev. 25 (4), 465–472.
• Ball, P., 2016. Patterns in Nature: Why the Natural World Looks the Way It Does. University of Chicago Press.
• Besag, J., 1974. Spatial interaction and the statistical analysis of lattice systems. J. R. Stat. Soc. Ser. B 192–236.
• Besag, J., York, J., Molli´e, A., 1991. Bayesian image restoration, with two applications in spatial statistics. Ann. Inst. Stat. Math. 43 (1), 1–20.
• Blangiardo, M., Cameletti, M., 2015. Spatial and Spatio-temporal Bayesian Models With r-inla. John Wiley & Sons.
• Blangiardo, M., Cameletti, M., Baio, G., Rue, H., 2013. Spatial and spatio-temporal models with r-inla. Spat. Spatiotemporal Epidemiol. 4, 33–49.
• Cai, Q., Lee, J., Eluru, N., Abdel-Aty, M., 2016. Macro-level pedestrian and bicycle crash analysis: incorporating spatial spillover effects in dual state count models. Accid. Anal. Prev. 93, 14–22.
• Cai, Q., Abdel-Aty, M., Lee, J., Wang, L., Wang, X., 2018. Developing a grouped random parameters multivariate spatial model to explore zonal effects for segment and intersection crash modeling. Anal. Methods Accid. Res. 19, 1–15.
• Cheng, W., Gill, G.S., Dasu, R., Xie, M.Q., Jia, X.D., Zhou, J., 2017. Comparison of multivariate poisson lognormal spatial and temporal crash models to identify hot spots of intersections based on crash types. Accid. Anal. Prev. 99, 330–341.
• Cheng, W., Gill, G.S., Ensch, J.L., Kwong, J., Jia, X.D., 2018a. Multimodal crash frequency modeling: multivariate space-time models with alternate spatiotemporal interactions. Accid. Anal. Prev. 113, 159–170.
• Cheng, W., Gill, G.S., Zhang, Y.P., Cao, Z., 2018b. Bayesian spatiotemporal crash frequency models with mixture components for space-time interactions. Accid. Anal. Prev. 112, 84–93.
• Cheng, W., Gill, G.S., Zhang, Y.P., Vo, T., Wen, F., Li, Y.H., 2020. Exploring the modeling and site-ranking performance of bayesian spatiotemporal crash frequency models with mixture components. Accid. Anal. Prev. 135.
• Chiou, Y.C., Fu, C., 2015. Modeling crash frequency and severity with spatiotemporal dependence. Anal. Methods Accid. Res. 5-6, 43–58.
• 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.
• Clayton, D., 1996. Generalised linear mixed models. In: Gilks, W., Richardson, S., Spiegelhalter, D. (Eds.), Markov Chain Monte Carlo in Practice. Chapman and Hall, New York, pp. 275–301.
• Cui, H.P., Xie, K., Hu, B., Lin, H., Zhang, R., 2019. Analysis of bus speed using a multivariate conditional autoregressive model: contributing factors and spatiotemporal correlation. J. Transport. Eng. A-Syst. 145 (4).
• Dong, N., Huang, H., Lee, J., Gao, M., Abdel-Aty, M., 2016. Macroscopic hotspots identification: a bayesian spatio-temporal interaction approach. Accid. Anal. Prev. 92, 256–264.

• El-Basyouny, K., Sayed, T., 2009. Urban arterial accident prediction models with spatial effects. Transp. Res. Rec. (2102), 27–33.

• Elvik, R., Ulstein, H., Wifstad, K., Syrstad, R.S., Seeberg, A.R., Gulbrandsen, M.U., Welde, M., 2017. An empirical bayes before-after evaluation of road safety effects of a new motorway in norway. Accid. Anal. Prev. 108, 285–296.

• Geman, S., Geman, D., 1984. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell. PAMI-6 6, 721–741.

• Gilks, W.R., Richardson, S., Spiegelhalter, D., 1995. Markov Chain Monte Carlo in Practice. CRC press.

• Guo, F., Wang, X.S., Abdel-Aty, M.A., 2010. Modeling signalized intersection safety with corridor-level spatial correlations. Accid. Anal. Prev. 42 (1), 84–92.

• Jones, B., Janssen, L., Mannering, F., 1991. Analysis of the frequency and duration of freeway accidents in seattle. Accid. Anal. Prev. 23 (4), 239–255.

• Knorr-Held, L., 2000. Bayesian modelling of inseparable space-time variation in disease risk. Stat. Med. 19 (17–18), 2555–2567.

• Krueger, R., Bansal, P., Buddhavarapu, P., 2020. A new spatial count data model with bayesian additive regression trees for accident hot spot identification. Accid. Anal. Prev. 144, 105623.

• Lan, B., Persaud, B., Lyon, C., Bhim, R., 2009. Validation of a full bayes methodology for observational before–after road safety studies and application to evaluation of rural signal conversions. Accid. Anal. Prev. 41 (3), 574–580.

• Lee, J., Abdel-Aty, M., Choi, K., Huang, H., 2015. Multi-level hot zone identification for pedestrian safety. Accid. Anal. Prev. 76, 64–73.

• Li, Z.N., Chen, X.F., Ci, Y.S., Chen, C., Zhang, G.H., 2019. A hierarchical bayesian spatiotemporal random parameters approach for alcohol/drug impaired-driving crash frequency analysis. Anal. Methods Accid. Res. 21, 44–61.

• Liu, C.H., Sharma, A., 2017. Exploring spatio-temporal effects in traffic crash trend analysis. Anal. Methods Accid. Res. 16, 104–116.

• Liu, C.H., Sharma, A., 2018. Using the multivariate spatio-temporal bayesian model to analyze traffic crashes by severity. Anal. Methods Accid. Res. 17, 14–31.

• Liu, J., Hainen, A., Li, X., Nie, Q.F., Nambisan, S., 2019. Pedestrian injury severity in motor vehicle crashes: an integrated spatio-temporal modeling approach. Accid. Anal. Prev. 132, 12.

• Lord, D., Mannering, F., 2010. The statistical analysis of crash-frequency data: a review and assessment of methodological alternatives. Transp. Res. Part A Policy Pract. 44 (5), 291–305.

• Lord, D., Miranda-Moreno, L.F., 2008. Effects of low sample mean values and small sample size on the estimation of the fixed dispersion parameter of poisson-gamma models for modeling motor vehicle crashes: a bayesian perspective. Saf. Sci. 46 (5), 751–770.

• Lord, D., Persaud, B., 2000. Accident prediction models with and without trend: application of the generalized estimating equations procedure. Transport. Res. Rec. (1717), 102–108.

• Lord, D., Persaud, B.N., Trb, Trb, Trb, 2000. Accident Prediction Models With and Without Trend - Application of the Generalized Estimating Equations Procedure. Highway and Traffic Safety: Crash Data, Analysis Tools, and Statistical Methods: Safety and Human Performance. Transportation Research Board Natl Research Council, Washington, pp. 102–108.

• Lumley, T., Lumley, M.T., 2013. Package 'leaps'. Regression Subset Selection. Thomas Lumley Based on Fortran Code by Alan Miller. Available online: http://CRAN.R-pro ject.org/package=leaps (accessed on 18 March 2018)..

• Lunn, D.J., Thomas, A., Best, N., Spiegelhalter, D., 2000. Winbugs a bayesian modelling framework: concepts, structure, and extensibility. Stat. Comput. 10 (4), 325–337.

• Ma, J.M., Kockelman, K.M., 2006. Bayesian Multivariate Poisson Regression for Models of Injury Count, by Severity. Statistical Methods and Crash Prediction Modeling.

• Transportation Research Board Natl Research Council, Washington, pp. 24–34. Ma, X., Chen, S., Chen, F., 2017. Multivariate space-time modeling of crash frequencies by injury severity levels. Anal. Methods Accid. Res. 15, 29–40.

• Mannering, F., 2018. Temporal instability and the analysis of highway accident data. Anal. Methods Accid. Res. 17, 1–13.

• Mannering, F.L., Shankar, V., Bhat, C.R., 2016. Unobserved heterogeneity and the statistical analysis of highway accident data. Anal. Methods Accid. Res. 11, 1–16.

• Martin, J.-L., 2002. Relationship between crash rate and hourly traffic flow on interurban motorways. Accid. Anal. Prev. 34 (5), 619–629.

• Miaou, S.-P., Song, J.J., 2005. Bayesian ranking of sites for engineering safety improvements: decision parameter, treatability concept, statistical criterion, and spatial dependence. Accid. Anal. Prev. 37 (4), 699–720.

• Munira, S., Sener, I.N., Dai, B., 2020. A bayesian spatial poisson-lognormal model to examine pedestrian crash severity at signalized intersections. Accid. Anal. Prev. 144, 105679.

• Narayanamoorthy, S., Paleti, R., Bhat, C.R., 2013. On accommodating spatial dependence in bicycle and pedestrian injury counts by severity level. Transp. Res. Part B Methodol. 55, 245–264.

• O'brien, R.M., 2007. A caution regarding rules of thumb for variance inflation factors. Qual. Quant. 41 (5), 673–690.

• Persaud, B., Lan, B., Lyon, C., Bhim, R., 2010. Comparison of empirical bayes and full bayes approaches for before–after road safety evaluations. Accid. Anal. Prev. 42 (1), 38–43.

• Poch, M., Mannering, F., 1996. Negative binomial analysis of intersection-accident frequencies. J. Transp. Eng. 122 (2), 105–113.

• Pulugurtha, S.S., Sambhara, V.R., 2011. Pedestrian crash estimation models for signalized intersections. Accid. Anal. Prev. 43 (1), 439–446.

• Quddus, M.A., 2008. Modelling area-wide count outcomes with spatial correlation and heterogeneity: an analysis of london crash data. Accid. Anal. Prev. 40 (4), 1486–1497.

• Rue, H., Martino, S., Chopin, N., 2009. Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations. J. R. Stat. Soc. Series B Stat. Methodol. 71 (2), 319–392.

• 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. Series B Stat. Methodol. 64 (4), 583–639.

• Spiegelhalter, D., Thomas, A., Best, N., Lunn, D., 2005. Winbugs User Manual Mrc Biostatistics Unit. UK Cambridge.

• Tang, H., Gayah, V.V., Donnell, E.T., 2020. Crash modification factors for adaptive traffic signal control: an empirical bayes before-after study. Accid. Anal. Prev. 144, 105672.

• Truong, L.T., Kieu, L.M., Vu, T.A., 2016. Spatiotemporal and random parameter panel data models of traffic crash fatalities in vietnam. Accid. Anal. Prev. 94, 153–161.

• Ukkusuri, S., Hasan, S., Aziz, H.M.A., 2011. Random parameter model used to explain effects of built-environment characteristics on pedestrian crash frequency. Transp. Res. Rec. 2237 (1), 98–106.

• Ulak, M.B., Ozguven, E.E., Vanli, O.A., Dulebenets, M.A., Spainhour, L., 2018. Multivariate random parameter tobit modeling of crashes involving aging drivers, passengers, bicyclists, and pedestrians: spatiotemporal variations. Accid. Anal. Prev. 121, 1–13.

• Wang, X., Abdel-Aty, M., 2006. Temporal and spatial analyses of rear-end crashes at signalized intersections. Accid. Anal. Prev. 38 (6), 1137–1150.

• Wang, Y., Kockelman, K.M., 2013. A poisson-lognormal conditional-autoregressive model for multivariate spatial analysis of pedestrian crash counts across neighborhoods. Accid. Anal. Prev. 60, 71–84.

• Wang, X., Abdel-Aty, M., Brady, P.A., 2006. Crash estimation at signalized intersections: significant factors and temporal effect. Transp. Res. Rec. 1953 (1), 10–20.

• Wang, C., Quddus, M., Ison, S., 2013. A spatio-temporal analysis of the impact of congestion on traffic safety on major roads in the uk. Transportmetrica A-Transport Sci. 9 (2), 124–148.

• Wang, K., Bhowmik, T., Yasmin, S., Zhao, S.S., Eluru, N., Jackson, E., 2019. Multivariate copula temporal modeling of intersection crash consequence metrics: a joint estimation of injury severity, crash type, vehicle damage and driver error. Accid. Anal. Prev. 125, 188–197.

• Wen, H.Y., Zhang, X., Zeng, Q., Sze, N.N., 2019. Bayesian spatial-temporal model for the main and interaction effects of roadway and weather characteristics on freeway crash incidence. Accid. Anal. Prev. 132, 6.

• Whittle, P., 1954. On stationary processes in the plane. Biometrika 434–449.

• Who, 2018. Global Status Report on Road Safety 2018: Summary. World Health Organization.

• Xie, K., Wang, X.S., Huang, H.L., Chen, X.H., 2013. Corridor-level signalized intersection safety analysis in shanghai, china using bayesian hierarchical models. Accid. Anal. Prev. 50, 25–33.

• Xie, K., Wang, X., Ozbay, K., Yang, H., 2014. Crash frequency modeling for signalized intersections in a high-density urban road network. Anal. Methods Accid. Res. 2, 39–51.

• Xie, K., Ozbay, K., Yang, H., 2015. Spatial analysis of highway incident durations in the context of hurricane sandy. Accid. Anal. Prev. 74, 77–86.

• Xie, K., Ozbay, K., Kurkcu, A., Yang, H., 2017. Analysis of traffic crashes involving pedestrians using big data: investigation of contributing factors and identification of hotspots. Risk Anal. 37 (8), 1459–1476.

• Xie, K., Ozbay, K., Yang, H., 2019. A multivariate spatial approach to model crash counts by injury severity. Accid. Anal. Prev. 122, 189–198.

• Xiong, Y., Tobias, J.L., Mannering, F.L., 2014. The analysis of vehicle crash injuryseverity data: a markov switching approach with road-segment heterogeneity. Transp. Res. Part B Methodol. 67, 109–128.

• Yang, D., Xie, K., Ozbay, K., Yang, H., Budnick, N., 2019. Modeling of time-dependent safety performance using anonymized and aggregated smartphone-based dangerous driving event data. Accid. Anal. Prev. 132, 105286.

• Zeng, Q., Wen, H.Y., Huang, H.L., Pei, X., Wong, S.C., 2018. Incorporating temporal correlation into a multivariate random parameters tobit model for modeling crash rate by injury severity. Transportmetrica A-Transport Sci. 14 (3), 177–191.

• Zeng, Q., Guo, Q., Wong, S.C., Wen, H.Y., Huang, H.L., Pei, X., 2019. Jointly modeling area-level crash rates by severity: a bayesian multivariate random-parameters spatiotemporal tobit regression. Transportmetrica A-Transport Sci. 15 (2), 1867–1884.

• Zeng, Q., Wang, X., Wen, H., Yuan, Q., 2020a. An empirical investigation of the factors contributing to local-vehicle and non-local-vehicle crashes on freeway. J. Transp. Saf. Secur. 1–15.

• Zeng, Q., Wen, H., Wong, S., Huang, H., Guo, Q., Pei, X., 2020b. Spatial joint analysis for zonal daytime and nighttime crash frequencies using a bayesian bivariate conditional autoregressive model. J. Transp. Saf. Secur. 12 (4), 566–585.

• Zhang, S., Liu, X., Tang, J.J., Cheng, S.W., Qi, Y., Wang, Y.H., 2018. Spatio-temporal modeling of destination choice behavior through the bayesian hierarchical approach. Physica A 512, 537–551.

• Ziakopoulos, A., Yannis, G., 2020. A review of spatial approaches in road safety. Accid.