category: literaturenote citekey: chengbayesianspatiotemporalcrashfrequency2018 title: Bayesian spatiotemporal crash frequency models with mixture components for space-time interactions authors: "Cheng, Wen; Gill, Gurdiljot Singh; Zhang, Yongping; Cao, Zhong" year: 2018 date: 2018-00-00 2018 doi: 10.1016/j.aap.2017.12.020 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2017.12.020" zotero_key: 3MGF5U8S zotero_storage: G9VD3HH2 collections: transpordi_planeerimine folder: 001_artiklid/Road Safety Indicators firstAuthor: "Cheng, Wen"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

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

Bayesian spatiotemporal crash frequency models with mixture components for space-time interactions

Wen Chenga,⁎ , Gurdiljot Singh Gilla , Yongping Zhanga , Zhong Caob

^aDepartment of Civil Engineering, California State Polytechnic University, Pomona 3801, W. Temple Ave., Pomona, CA 91768, United States
^b Engineering Technology Research Center for Mathematical Educational Software, Guangzhou University, Waihuanxi Road No.230, Guangzhou Higher Education Mega Center, China

ARTICLE INFO

Keywords: Mixture Spatiotemporal Time adjacency Interaction Predictive accuracy Cross validation

ABSTRACT

The traffic safety research has developed spatiotemporal models to explore the variations in the spatial pattern of crash risk over time. Many studies observed notable benefits associated with the inclusion of spatial and temporal correlation and their interactions. However, the safety literature lacks sufficient research for the comparison of different temporal treatments and their interaction with spatial component. This study developed four spatiotemporal models with varying complexity due to the different temporal treatments such as (I) linear time trend; (II) quadratic time trend; (III) Autoregressive-1 (AR-1); and (IV) time adjacency. Moreover, the study introduced a flexible two-component mixture for the space-time interaction which allows greater flexibility compared to the traditional linear space-time interaction. The mixture component allows the accommodation of global space-time interaction as well as the departures from the overall spatial and temporal risk patterns. This study performed a comprehensive assessment of mixture models based on the diverse criteria pertaining to goodness-of-fit, cross-validation and evaluation based on in-sample data for predictive accuracy of crash estimates.

The assessment of model performance in terms of goodness-of-fit clearly established the superiority of the time-adjacency specification which was evidently more complex due to the addition of information borrowed from neighboring years, but this addition of parameters allowed significant advantage at posterior deviance which subsequently benefited overall fit to crash data. The Base models were also developed to study the comparison between the proposed mixture and traditional space-time components for each temporal model. The mixture models consistently outperformed the corresponding Base models due to the advantages of much lower deviance.

For cross-validation comparison of predictive accuracy, linear time trend model was adjudged the best as it recorded the highest value of log pseudo marginal likelihood (LPML). Four other evaluation criteria were considered for typical validation using the same data for model development. Under each criterion, observed crash counts were compared with three types of data containing Bayesian estimated, normal predicted, and model replicated ones. The linear model again performed the best in most scenarios except one case of using model replicated data and two cases involving prediction without including random effects. These phenomena indicated the mediocre performance of linear trend when random effects were excluded for evaluation. This might be due to the flexible mixture space-time interaction which can efficiently absorb the residual variability escaping from the predictable part of the model. The comparison of Base and mixture models in terms of prediction accuracy further bolstered the superiority of the mixture models as the mixture ones generated more precise estimated crash counts across all four models, suggesting that the advantages associated with mixture component at model fit were transferable to prediction accuracy. Finally, the residual analysis demonstrated the consistently superior performance of random effect models which validates the importance of incorporating the correlation structures to account for unobserved heterogeneity.

E-mail addresses: wcheng@cpp.edu (W. Cheng), gurdiljotg@cpp.edu (G.S. Gill), yongpingz@cpp.edu (Y. Zhang), zhongc@gzhu.edu.cn (Z. Cao).

^⁎ Corresponding author.

1. Introduction

Roadway crashes have caused an immense burden on society with respect to emotional and financial losses. Researchers are entrusted to develop analytic approaches to gain a better understanding of the causal factors for crash occurrence and develop more accurate crash prediction models to formulate road-safety policies and engineering solutions for mitigation of crashes. However, the accuracy of inferences drawn from the statistical analysis of crashes is highly dependent on the robustness of crash data and an array of potential influential factors such as roadway geometric (number of lanes, lane width, radius of horizontal curve, etc.), traffic flow (vehicle density, volume, real-time speed, speed deviations, etc.), environment (lighting, weather), driver characteristics and mental state (gender, response time, age, etc.), among others. Unfortunately, the crash related data collected by safety agencies may be inadequate or unavailable for detailed investigation (Lord and Mannering, 2010). Hence, the researchers managed to handle this issue to study the significant factors by virtually enhancing the quantity of dataset by disaggregation over some geographical space (micro or macro level) and some specified time period (e.g., division of five year accumulated crash data into five individual subsets). The crash-frequency data are obtained in the form of non-negative integers allowing the application of count-based regression models.

These regression models (or, crash prediction models) have been used in research and practice for determination of influential factors, planning purposes, or site ranking. Models of varying complexity have been employed, ranging from very basic to sophisticated. The traditional approach to analyzing roadway crashes employed generalized linear models (McCullagh and Nelder, 1989;Zeger and Karim, 1992) to establish a linear relationship between explanatory variables and logtransformed outcomes such as crash frequencies of different severities or vehicle modes. This allowed for clear interpretation of inferences drawn from model estimates. To handle over-dispersion commonly associated with crash data, over-dispersed generalized linear models such as Poisson mixtures (e.g., negative binomial or Poisson-gamma, Poisson-lognormal, etc.) were introduced (Persaud, 1994; Hauer, 1997; Milton and Mannering, 1998; Karlaftis and Tarko, 1998). These models may not fully incorporate the unobserved heterogeneity as in case of count-data models, the overdispersion may be attributed to various factors, such as the grouping of data over space (segments, neighborhood, cities, regions, etc.), unaccounted temporal correlation, and model miss-specification (Gourieroux and Visser, 1997;Poortema, 1999;Colin Cameron and Trivedi, 1998).

Many studies have been proposed to address the unobserved heterogeneity shared by roadway entities in close proximity. These studies have focused on various spatial units including micro-levels such as intersections (Wang and Abdel-Aty, 2006; Cheng et al., 2017a), segments (Aguero-Valverde and Jovanis 2008), and macro-levels like corridors (Abdel-Aty and Wang, 2006; Guo et al. 2010), census tracts (MacNab, 2004), traffic analysis zones (Washington et al. 2010), counties (Gill et al., 2017), and so on. The significance of incorporating spatial correlations was highlighted by some studies (Guo et al., 2010; Abdel-Aty & Wang, 2006) with the consistently better performance of the spatial models over those accounting for heterogeneity random effect only.

In addition to spatial correlations, the disaggregation of crash data over specified time periods also leads to temporal correlation as those datasets may share unobserved effects which remain constant over time. To remove the potential bias of estimated model parameters, some researchers addressed the serial correlation in crash data by employing different temporal treatments such as linear and/or quadratic trend (Andrey and Yagar, 1993; Hay and Pettitt, 2001), autoregressive correlation structure with a time step of one year (lag-1) (Huang et al., 2009; Wang et al., 2013; Cheng et al., 2018), fixed-over-time and independent-over-time random effects (Aguero-Valverde, 2013; Jiang et al., 2014), and time-varying model coefficients and intercepts (Cheng et al., 2017b), and so on.; These models revealed that accounting for temporal correlations significantly improved the capability of the models to fit the crash data.

Building on the advantages of spatial and temporal correlation structures to address the issue of unobserved heterogeneity, some studies incorporated temporal dimension for spatial models as the crash analysis is not curbed to a single time period. Wang and Abdel-Aty (2006) investigated the rear-end crashes at intersections while employing the generalized estimating equations with the negative binomial link function to account for the temporal and spatial correlation among the 3-year longitudinal crash data. The study by Blazquez and Celis (2013) did a spatial and temporal analysis of child pedestrian crashes occurring during a period of nine years. The spatial autocorrelation analysis indicated that the responsibility of pedestrians is the major contributing factor for the generation of child pedestrian crashes with a tendency to cluster in space and time. Also, spatial clustering distribution of crashes in terms of time of the day was also observed. While these studies treated spatial and temporal correlations independently, some studies noted that vehicle crashes tend to cluster both spatially and temporally, hence space-time interaction specification was employed at different spatial scales. Aguero-Valverde and Jovanis (2006) proposed a spatiotemporal model with a linear timespace interaction term to study the fatal and injury crashes in Pennsylvania. The authors noticed that spatial correlation, time trend, and space–time interactions were significant in the proposed county-level Bayesian crash models. They also recommended such model should be extended to road segment and intersection-level crash models, where spatial correlation is likely to be even more pronounced. Subsequently, Plug et al. (2011) explored the variation of spatial distribution of single vehicle crashes (SVCs) according to different time periods (time of day and day of the week) by employing visualization technologies. The results showed significant differences in spatiotemporal patterns of SVCs for various crash causes.

The literature review illustrated notable benefits at various fronts associated with inclusion of spatial and temporal correlations and their interactions. However, in comparison with other types of models, a very limited body of research dedicated to spatiotemporal models exists in the field of traffic safety. Moreover, most of the current limited spatialtemporal models in the field assume a linear temporal trend and linear space-time interaction which may be seen as a restrictive assumption (Lawson and Clark, 2002; Lawson et al., 2003). For example, as discussed previously, the temporal random effects may take on nonlinear shape or have autocorrelation with previous crash counts. In addition, the temporal trend might have a non-linear change across the spatial units. To add to the current literature with more spatiotemporal models, the authors aimed to develop four alternative spatial-temporal models which employ different temporal treatments with the varying complexity of random effects: (I) linear time trend; (II) quadratic time trend; (III) Autoregressive-1 (AR-1); and (IV) time adjacency. Furthermore, instead of using linear space-time interaction, the authors borrowed a flexible two-component-mixture interaction from one previous disease-mapping study (Abellan et al., 2008). Such mixture can easily capture the space-time trends that depart from the predictable patterns of overall temporal and spatial risk surface as it allows the smoothness as well as discontinuities in the space-time variations within the roadway entities. The interested readers can be referred to this study for more details of the mixture components. The study results demonstrated a number of benefits associated with the proposed mixture model. However, the performance of the mixture model in traffic safety area is unknown and is therefore worth studying. In addition, in order for a comprehensive comparison of the predictive accuracy of the four models, five evaluation criteria were utilized which include log pseudo marginal likelihoods (LPML), mean square predictive error (MSPE), mean absolute deviation (MAD), residual sum of squares (RSS), and total rank difference (TRD). LPML assesses the predictive capability using cross-validation while the rest utilize the same dataset for model development. Three different approaches (model replicated, model predicted, and Bayesian estimated) were used to generate the crash dataset for evaluation as an effort to quantify the prediction performance and understand the benefits of concerned models from different perspectives. Finally, a residual analysis was also performed to observe the impact of the inclusion of random effects on the prediction accuracy.

2. Methodology

This section is divided into two main subsections: Model Development and Evaluation Criteria. The former subsection first describes the formulation of Poisson lognormal model which was employed for the estimation of crash counts; then introduces the Base model specification which is usually employed in conventional safety studies focused on space-time interaction and discusses the advantage of Mixture model over the Base; and finally, four different types of temporal treatments are introduced which represent the temporal specifications employed in different studies. In the second subsection (Evaluation Criteria), different criteria are discussed which were used for assessment of model goodness-of-fit and comparison of models' predictive accuracy based on both out-of-sample (and, the data used for cross validation) and in-sample crash data (or, the same data used for model development).

2.1. Model development

2.1.1. Model specification

The Full Bayesian (FB) framework was employed to estimate the crash frequency of vehicles at the segment level. To account for the unobserved heterogeneity from different perspectives, the FB approach is preferable due to its flexibility and effectiveness to incorporate complex correlations with a hierarchical structure of data. FB approach provides a posterior distribution of parameters from Markov-chain Monte Carlo (MCMC) simulation which samples the variables as random, unlike the point estimates generated by the traditional approach of maximum likelihood estimation. This approach has been widely used for crash prediction models due to the multilevel and correlated nature of data (Lord and Mannering, 2010). In order to account for the overdispersion usually associated with crash data, the present study chose Poisson-Lognormal model given its capability of better handling of low sample mean and small sample size due to the heavier tails associated with Lognormal distribution (Lord and Miranda-Moreno, 2008). This model assumes that crash count (yi) at a given segment i, obeys Poisson distribution, while the corresponding observation specific error term εⁱ follows a normal distribution:

$$y_i | \lambda_i \sim Poisson(\lambda_i)$$ (1)

Where λⁱ is the expected Poisson count for location i, and can be expressed as follows:

$$\ln(\lambda{i}) = \beta{0} + X{i}^{'}\beta + \varepsilon{i} \tag{2}$$

Where β⁰ is the intercept, X' is the matrix of risk factors, β is the vector of model parameters, εⁱ is the independent random effect which captures the extra-Poisson heterogeneity among locations. εⁱ can be assumed with the following noninformative normal priors:

$$\varepsilon_i \sim Normal(0, \tau^2)$$ (3)

Where τ ² is the variance of the normal distribution for εi. The inverse of τ ² is called precision and it can be modeled using the following gamma prior with prior mean equal to one and its prior variance large (equal to one thousand), representing high uncertainty or prior ignorance:

$$\tau^{2-1} \sim Gamma(0.001, 0.001)$$ (4)

The model presented in Eq. (2) serves as the basic formulation for

crash prediction in most safety studies. The current study also employs this model as the basis for development of spatiotemporal models.

2.1.2. Two groups of models

This section first introduces the Base model specification which serves to represent the assignment of spatial, temporal, and space-time interactions by conventional traffic safety studies. The limitation associated with such specification is discussed and then the mixture model specification is introduced.

2.1.2.1. Base models. As previously mentioned, the spatial random effects were incorporated over the model represented in Eq. (2). The subsequent model which accounts for the spatial correlation among neighboring segments is obtained with the following formulation:

$$\ln(\lambda{i}) = \beta{0} + X{i}'\beta + \varepsilon{i} + u_{i}$$ (5)

Where ui is the spatially structured random effect which follows the CAR (conditional autoregressive) formulation to incorporate the spatial correlation among crashes occurring at neighboring segments.

$$u_i \sim N\left(\frac{\overline{u}_i, 1}{\tau_i}\right) \overline{u}i = \frac{\sum{i \neq k} uk w{ik}}{\sum{i \neq k} w{ik}} \tau_i = \frac{\tauc}{\sum{i \neq k} w_{ik}}$$ (6)

Where wik is the proximity weight matrix and τc is the precision parameter in the CAR prior. As evident from the above equation, estimation of the risk in any site is conditional on risks in neighboring locations. Subscripts i and k refer to a segment and its neighbor, respectively, and k ε Ni, where Ni represents the set of neighbors of segment i. Besides the identification of neighbors, the assigned weights also affect the risk estimation. In the past studies (Aguero-Valverde and Jovanis, 2010; Xu and Huang, 2015), weight structures including various adjacency-based, distance-based models, and semi-parametric geographically weighted, and so on, have been explored for incorporating the spatial correlations at different levels of roadway entities. The current study employs the adjacency-based structure to explore the spatial correlation, which is commonly used for a microlevel entity of segments. The weights for the adjacency matrix are given by (wik) = 1 if segments i, k are adjacent, and 0 otherwise. It is noteworthy that this adjacency-based weight matrix gives a dichotomous output for weight with only two responses, zero for nonneighbors and one for neighbors. Hence, the proximity of immediate neighbors establishes a spatial dependency while the impact of more distant neighbors is assumed to be negligible.

If temporal treatment also needs to be considered, then Eq. (5) can be easily modified and extended to account for serial correlations and associated space-time interaction. The resultant formulation can be expressed in Eq. (7).

$$\ln(\lambda_{it}) = \beta_0 + X_i'\beta + f(t) + ui + v{it}$$ (7)

Where f t( ) is the temporal treatment, and vit is the space-time interaction parameter. vit can take different forms such as linear (Aguero-Valverde and Jovanis, 2006) and time-varying spatial correlation (Lawson et al., 2003). However, the existing space-time interaction models have a limitation in the sense that they assume all the sites under focus experience a global space-time pattern where the error (i.e., the part not accounted by isolated spatial and temporal terms) falls into one distribution common to all sites. This may be interpreted as the inefficiency to distinguish between the stable and unstable underlying risk pattern, where the stable (or predictable) patterns are represented by the overall time trend (f (t)) and spatial crash risk (ui), and unstable (or unpredictable) patterns are the atypical departures from the stable pattern. This study introduces the mixture component for accommodating the large fluctuations of the space-time interaction which allows the identification of unstable locations which do not conform to the global crash risk pattern of space and time variability.

2.1.2.2. Mixture component. The mixture component for space-time interaction was introduced by Abellan et al., (2008) to "uncover the full space-time profile of risks" for an epidemiologic study. The concerned study proposed a 2-component mixture which combines variabilities of different levels to represent the space-time interaction. The accommodation of different levels of variabilities is required to distinguish the stable and unstable pattern. Specifically, vit from Eq. (7) takes the following form:

$$v_{it} \sim p \text{Normal}(0, \tau_1^2) + (1-p) \text{Normal}(0, \tau_2^2)$$ (8)

The prior for p is uniform on [0, 1] through a Dirichlet distribution. This varied prior distribution renders the flexibility to accommodate discreet and continuous space-time variations in crash data, where the discreet variations refer to the outstanding or outlier sites in the plan of space-time interaction which deviate from the continuous global pattern of space-time variations. The first component captures the small residual noise which is generally adequately captured by the overall time trend and spatial component. The second component serves the purpose of capturing the "true" departure from the global pattern by accommodating the large fluctuations which indicate the higher risk relative to other sites. Half-normal hyper-prior distributions for the standard deviations τk, (k = 1, 2), are specified which reflect that τ¹ (corresponding to "stable" risk pattern) has to be small to effect shrinkage so as to capture the noise, whereas the prior for τ² (corresponding to "unstable" risk pattern) allows a large range of values for this parameter to capture the large fluctuations. It is worth mentioning that other positive distributions available in WINBUGS such as gamma, exponential can also be tested and compared with the half-normal distribution used hereby. The different range of variance is shown in the following equations:

$$\tau1 \sim Normal(0, 0.01). I{(0, +\infty)}$$ (9)

$$\tau2 \sim Normal(0, 0.01). I{(0,+\infty)} + Normal(0, 100). I_{(0,+\infty)}$$ (10)

Based on both simulation data and real case study, the authors revealed a number of benefits of the mixture model in terms of interpretation and potential for detection of localized excess for diseases. As a result, the present study borrowed such structure for the space-time interaction and tested its capability in addressing the prediction accuracy for the traffic crash data. Additionally, another characteristic of the current research lies in the different ways of handling temporal treatment (i.e., f t( ) in Eq. (7)).

2.1.3. Four types of temporal treatments

The following subsections present the corresponding model formulations in order of complexity, from the linear time trend to more sophisticated ones. These models were developed using both the Base and Mixture specifications. It is noteworthy that all models vary on the basis of the assignment of different temporal treatments over the Poisson lognormal models which incorporate the spatial correlation and space-time interaction.

2.1.3.1. Model 1: Linear time trend. In this model, a linear trend is introduced where time is regarded as a potential influential covariate, and the model estimates the coefficient for the trend. Eq. (7) assumes the following form:

$$\ln(\lambda_{it}) = \beta_0 + X_i'\beta + ui + v{it} + \beta_t * T$$ (11)

Where βt is the scalar parameter for linear yearly trend T (T = 1 to 6 in the present study). It is expected that the passage of time will impact the crash risk on the segments which would be captured by the model.

2.1.3.2. Model 2: Quadratic time trend. In this model, a non-linear impact of time is considered. Eq. (7) takes the following form:

$$\ln(\lambda_{it}) = \beta_0 + X_i'\beta + ui + v{it} + \betat * T + \beta{t2} * T^2$$ (12)

Where βt² is the coefficient for the quadratic trend term. It is expected that this model will allow more subtle approach compared to the linear trend as it allows flexibility to fit the crash data by virtue of its additional quadratic trend.

2.1.3.3. Model 3: Autoregressive-1 (AR-1). This model accounts for the autoregressive safety effect by specifying the distribution of εit as a lag-1 dependence on errors, where lag-1 means that the time is varying yearly. It incorporates the weighted sum of the past one year of values together with a random term. It is important to note that, in addition to AR-1, there are a large number of models involving autoregressive error process (Miaou and Song, 2005), such as higher order AR model AR(p), autoregressive–moving-average (ARMA), Integer valued autoregressive (INAR), Autoregressive Conditional Heteroscedastic (ARCH), and Generalized Autoregressive Conditional Heteroscedastic (GARCH), and so on. The authors opted for the popular AR-1 specification to represent such models based on the assumptions of stationarity restriction. Under this model, Eq. (7) assumes the following form:

$$\ln(\lambda_{it}) = \beta_0 + Xi^{'}\beta + \varepsilon{it} + ui + v{it}$$ (13)

The weighted sum is fixed, and the random terms change at every time step following the same distribution, which means this model is homoscedastic. The distributions are given by:

$$\varepsilon{it} \sim normal\left(0, \frac{\sigma{it}^2}{(1-\gamma^2)}\right)$$ (14)

$$\varepsilon{it} \sim normal(\gamma \varepsilon{i,t-1}, \sigma_{it}^2) \text{ for } t > 1$$ (15)

Where γ is the autocorrelation coefficient with the range of −1 < γ < 1. This model addresses the potential correlation between successive time periods and is expected to deliver a more precise estimation of the model parameters.

2.1.3.4. Model 4: Time adjacency. In this model, a CAR formulation is employed for the temporal aspect of the data, similar to the spatial aspect discussed previously. In this case, Eq. (7) assumes the following form:

$$\ln(\lambda_{it}) = \beta_0 + X_i'\beta + ui + v{it} + \alpha_t$$ (16)

Where αt is the temporal correlation term which follows the CAR formulation given in Eq. (6). The six-year time period of this study assumes strong influence of the adjacent years on crash risk, similar to the dependence of adjacent segments in case of spatial correlation. Compared to the AR-1 model which captures the impact of one year prior to the concerned time, this model also accounts for the potential impact of the subsequent year. The similar formation of such temporal treatment can be found in a previous study (Abellan et al., 2008).

2.2. Evaluation criteria

This study assessed the model performance by employing different evaluation criteria which were found in the past research. The authors intend to evaluate the model performance from different perspectives such as goodness-of-fit (while accounting for complexity), and out-ofsample cross validation and in-sample predictive accuracy for crash counts. The brief background and analytics of each criterion are provided as follows.

2.2.1. DIC and its components

The Deviance Information Criterion (DIC) (Spiegelhalter et al., 2003) served as the measure for assessment of goodness-of-fit for the models. It is a penalized fit criterion which can be calculated from the samples generated by an MCMC simulation. Usually, the estimation of DIC can be conducted through the following expression.

$$\begin{cases} DIC = \overline{D} + P_D \ D(\theta) = -2 \cdot \log(p(y\theta)) \ \overline{D} = E^{\theta}[D(\theta)], P_D = \overline{D} - D(\overline{\theta}) \end{cases}$$ (17)

Where, y is the data; θ are the unknown parameters of the model; p yθ ( ) is the likelihood function; is the expectation of ^θ; ⁻^D is the posterior mean deviance, or the expectation of the deviance, D θ( ); ⁻^D measures how well the model fits the data; the smaller the ⁻D, the better the fit. PD represents the effective number of parameters. Usually, smaller values of DIC are preferred. Some general guidelines are also suggested by Lunn et al. (2012): models with DIC values within 5 of the 'best' model are comparable and also strongly supported, values within 5 and 10, weakly supported, and models with a DIC greater than 10 are substantially inferior.

2.2.2. Comparison of models based on out-of-sample cross-validatory predictive density

This study employed the conditional predictive ordinate (CPO) (Gelfand, 1996) for cross validation based on predictive densities. To prevent potential bias due to random selection of data division into subsets (one subset for development and another for checking), the current study calculated the CPO by implementing a CV-1 (leave-oneout) distribution which removed the selection bias by employing a continuous approach of selecting all data points, except one, for model development and the left out data point to verify the prediction accuracy of the calibrated model. Such implementation also has benefits over traditional modeling evaluations which are prone to overestimation due to double usage of data (once during model development and then again for model checking). Under the MCMC framework, the estimate of CPO for each observation i can be calculated as:

$$CPO = \left(\frac{1}{T} \sum_{t=1}^{T} \frac{1}{f(Y_t | \beta^{(t)})}\right)^{-1}$$ (18)

Where Yi is the ith observation (i = 1, 2, 3,. . ., n) for all segments and β is the vector of estimated model parameters. This harmonic mean of density (CPO) may be extended to calculate the goodness-of-fit of models by computing the product of CPOs over all observations, which is known as the pseudo marginal likelihood. For computational convenience, the log pseudo marginal likelihoods (LPML) is calculated:

$$LPML = \sum_{i=1}^{n} \log(CPO_i)$$ (19)

2.2.3. Evaluation criteria for predictive accuracy based on the same dataset In this study, to obtain a more comprehensive comparison of four alternate models, the traditional evaluation criteria such as MSPE, MAD, RSS, and TRD were also calculated to assess the crash prediction accuracy of the models while treating the observed crash counts as the reference. These criteria are expected to quantify the variability in the prediction performance of models, which will help understand the benefits of concerned models from different perspectives. The evaluation was conducted by employing three different approaches, essentially governed by the selection of datasets for comparison with the observed crash counts. In the first case, the Bayesian estimated crash frequency was utilized which was obtained from Eqs. (11)–(13), and 16 for corresponding models. The difference of Bayesian estimated and observed crash counts is referred as Bayesian residuals (Lawson et al., 2003). In the second case, the model predicted crash frequency was compared with the observed counts. To obtain the model predicted

$$ln(PD_{it}) = \beta_0 + X_i'\beta \tag{20}$$

counts, Eq. (7) may be modified as follows:

Where PDit is regarded as the model predicted crash count. The comparison of model predicted and observed crash counts is referred as normal residuals. This predicted count differs from the Bayesian estimated as no random effect terms are incorporated, hence the comparison of Bayesian and normal residuals is expected to quantify the variability in predictive performance of models without the impact of spatial and temporal correlation structures. In the third case, the model replicated crash counts are utilized for comparison with observed counts. The replicated counts were generated from the posterior predictive distribution by introducing a simple step in MCMC sampler as follows:

$$y(rep)_i | \lambda(rep)_i \sim Poisson(\lambda_i)$$ (21)

Where λ(rep) refers to the expected crash counts estimated based on the assumption that adopted model truly fits the observed crash data. The closeness of λ(rep) and λi signifies that the model was able to precisely replicate the crash counts based on the observed ones and hence may be employed for accurate representation of crash risk in the future (Lawson et al., 2003). For the sake of simplicity in the description of evaluation criteria in the subsequent paragraph, the word "estimated" will be used for Bayesian estimated, model predicted, and model replicated crash counts.

As previously stated, four evaluation criteria were employed to assess prediction accuracy of the four competing models from different perspectives based on three cases discussed above, namely: mean square predictive error (MSPE), mean absolute deviation (MAD), residual sum of squares (RSS), and total rank difference (TRD). Specifically, the MSPE was calculated as follows:

$$MSPE = \frac{1}{n} \sum_{i=1}^{n} (Y_i - O_i)^2$$ (22)

Where Yi is the estimated crash frequency and Oi is the observed crash count for entity i by a model during the same time period. The smaller value of MSPE for a model indicates that the discrepancy between the estimated crash counts for the sites is relatively lower which demonstrates the accuracy of the model at prediction. The second criterion that was employed for assessment of prediction accuracy is MAD, which is defined as:

$$MAD = \frac{1}{n} \sum_{i=1}^{n} |Y_i - O_i|$$ (23)

Where terms are as defined previously. The third criterion, Chi-squared RSS (Residual Sum of Squares) is defined as:

$$RSS = \sum{i=1}^{n} \frac{(O{i} - Y{i})^{2}}{Y{i}}$$ (24)

Finally, TRD was employed as a more sensitive criterion to account for the rank deviations based on the observed and estimated crash counts (Cheng and Washington, 2008). The rank difference is calculated by using the following equation.

$$TRD = \sum_{i=1}^{n} |R(i_o) - R(i_e)|$$ (25)

Where R( ) io is the observed data rank at entity i and R(ie) is the rank based on estimated crash counts of each of the three approaches for the same time period. A particular model is considered more reliable if smaller TRD value is observed, which signifies that the specific model assigns rankings close to the observed crash counts.

3. Data description

The data used for this study were provided by HSIS (Highway Safety Information System) which collected the data in the form of different raw files from California TASAS (Traffic Accident Surveillance and Analysis System) and NHTSA (National Highway Traffic Safety Administration). The data from each year had four files pertaining to different types of data linked with Road, Occupancy, Vehicle, and Crash characteristics. The information extracted from these four files had crash number along with other factors like geometric (lane width,

Table 1 Summary Statistics of Variables for segments.

Variable	Mean	SD	Minimum	Maximum
Total vehicle crash counts per year	0.34	0.81	0	8
Lane width 9 and 10 ft (yes or no)	0.33	0.47	0	1
Lane width 11 ft (yes or no)	0.10	0.30	0	1
Lane width 12 ft (yes or no)	0.47	0.50	0	1
Lane width greater than 12 ft (yes or no)	Ref
Speed limit 25–40 mph (yes or no)	0.247	0.431	0	1
Roadway surface width (ft)	23.35	3.22	20	48
Speed limit 45 mph (yes or no)	0.616	0.486	0	1
Speed limit 50–65 mph (yes or no)	Ref
Left Shoulder width (ft)	2.824	2.163	0	12
Right Shoulder width (ft)	2.800	2.270	0	12
AADT	3,149.80	3,818	610	21,241
Plain cement concrete Bridge Deck (yes or no)	0.12	0.32	0	1
Asphalt concrete 7 in thick (yes or no)	0.34	0.47	0	1
Asphalt Concrete (Oiled Earth-Gravel) (yes or no)	Ref

¹ SD refers to standard deviation.

number of lanes, median type), traffic (Average Annual Daily Traffic), design speed, and so on. The crash data for 279 segments in one rural roadway section of California (State Route 23) was considered for model development over a period of six years (2007–2012). The 11 mile roadway section runs from the City of Moorpark at the south to the City of Fillmore at the north. The roadway has 2-lane in each direction and mainly follows within a mountainous area. This relatively long period was required to fully incorporate the temporal trends and correlation structures for different treatments of time in four alternate models. The summary statistics of crash counts and roadway pertinent covariates are shown in Table 1.

4. Results

4.1. Modeling results

The models were developed in the statistical software WinBUGS (Lunn et al., 2000) to generate MCMC samples for Bayesian posterior inferences. For the model calibration, two chains of 13,500 iterations were set up. Convergence was ensured by visual inspection of chains and observing the desired threshold condition of MC errors to be lower than 5% of the standard deviation of parameters. After ensuring the convergence, first 9000 samples were discarded as adaptation and burn-in, and rest of the samples were used to draw parameter estimates. The model complexity was observed to have a direct impact on model running time, and this criterion may influence the selection of models by safety agencies in terms of the trade-off between efficiency and precision. AR-1 (Model 3) and time adjacency (Model 4) took the highest running time for the same number of iterations, which was almost 7% higher than the most efficient model (Model 1). In case of the independent variables incorporated for model development, it is important to note that indicator variables (categorical) were used for lane width, speed limit, and surface type, with the baseline or reference as lane width more than 12 ft, speed limit 50–65 mph, and Asphalt Concrete (Oiled Earth-Gravel), respectively. The variables selected for model development were checked for multi-collinearity and correlation issues and only the satisfying variables were incorporated as an effort to arrive at potentially less biased inferences drawn from relatively more precise estimates. The right shoulder width was therefore excluded from model development due to its high correlation with the left should width, so was roadway surface width due to high collinearity with lane width. As shown in Table 2, all four models identified the vehicle volume (AADT) to be statistically significant along with Speed limit of 25–40 mph for two models and width of left shoulder and Plain cement concrete Bridge Deck type for Model 4. The parameter estimates illustrate the robustness of the developed models while employing different specifications of temporal aspect of crash data. In specific, the larger AADT would lead to more crash occurrences, while the narrower left shoulder width tends to increase crash probabilities. Interestingly, in comparison with the reference speed limit (50–65 mph), the segments with smaller speed limits (such as 25–40 mph, 45 mph) are inclined to generate vehicle crashes. The potential explanation might be due to worse driving conditions (e.g., sharp curves, limited visions) associated with the segments of lower speed limits in the mountainous area. It should be noted that the focus of the current study is on the model comparison based on predictive accuracy rather than the ability to filter influential factors, hence the explanation of results is focused in that direction hereafter.

4.2. Goodness-of-fit

This study employed DIC and its components for assessment of goodness-of-fit of alternate models. DIC was selected for the model comparison as it is a penalized criterion which acts as a trade-off between model complexity and model fit which are represented by the ^effective number of parameters (Pd) and posterior deviance (−D). It was expected that mixture models would generate better fit for the crash data due to the addition of more parameters compared to the Base models. However, the advantage of model fit was expected to be coupled with higher complexity. As shown in Table 3, the mixture models consistently exhibited higher complexity for each model, as demonstrated by the variation in values of Pd. The highest difference was depicted by Model 1 where the Pd value seems to indicate that effective number of parameters employed by the mixture model may be three times the corresponding Base model (Pd for mixture was 174 while for Base it was 54). The least discrepancy between the complexity of Base and mixture was observed for Model 4, which also demonstrated the only case of comparable DIC (difference of only 1 point). The understandable higher complexity associated with the mixture models was compensated by the significantly lower (relative to Base) posterior deviance across all models which subsequently resulted in the consistently lower DIC values. These findings clearly demonstrate that the addition of mixture component for spatiotemporal models provides the advantage of superiority at fitting the observed crash data. The significant differences in the DIC values of mixture and Base models, which were the largest for the simplest Model 1 and slowly reduced which progression of model complexity (Model 4), clearly conveys the flexibility of mixture components. In terms of comparison between different temporal treatments, Model 4 demonstrated the best overall fit with the lowest DIC value of 2,043, which was 84 points lower than Model 3 (2,127). Even though the Pd value for Model 4 indicates the higher complexity but the significant superiority for posterior deviance generated the advantage at overall fit. This may be explained by the fact that AR-1 employs the lag-1 technique to borrow information from the past year which may not be as efficient as consideration of temporal correlation among all years due to the time matrix employed by Model 4. Based solely on the goodness-of-fit criteria, the time-matrix specification for addressing the temporal aspect performed the best, followed by linear trend and quadratic, while the AR-1 specification demonstrated the worst performance.

4.3. Predictive accuracy

4.3.1. Cross-validation

The cross-validation was adopted as one approach to compute the conditional predictive ordinate (CPO) and eventually calculate log pseudo marginal likelihood (LPML) for comparison of model fit. The higher value of LPML reflects relatively superior model fit property. The difference between the LPML values of two concerned models is

² Ref signifies the variables which act as reference for dummy variables.

Table 2 Statistics of estimated model parameters.

Model	Variable Coefficient	Mean	SD	MC error	2.50%	97.50%
Model 1	Intercept	−6.734	0.5839	0.0473	−7.888	−5.662
	Lane width 9 and 10 ft	0.0903	0.2444	0.0148	−0.3561	0.6136
	Lane width 11 ft	0.0096	0.2546	0.0138	−0.4578	0.5429
	Lane width 12 ft	−0.1744	0.2232	0.0114	−0.5953	0.2797
	Speed limit 25–40 mph	0.4311	0.1939	0.0094	0.0611	0.8273
	Speed limit 45 mph	0.1452	0.1838	0.0093	−0.1954	0.5268
	Left Shoulder width	−0.0242	0.0259	0.0011	−0.0769	0.0248
	AADT	0.8724	0.0699	0.0057	0.7392	1.0200
	Plain cement concrete Bridge Deck	−0.1862	0.3201	0.0172	−0.8171	0.4323
	Asphalt concrete 7 in thick	−0.0032	0.2107	0.0110	−0.4280	0.4028
	Linear trend	−0.0309	0.0288	0.0011	−0.0895	0.0237
Model 2	Intercept	−6.6830	0.7520	0.0627	−8.1280	−5.4420
	Lane width 9 and 10 ft	0.0981	0.2497	0.0144	−0.3734	0.5923
	Lane width 11 ft	0.0069	0.2603	0.0140	−0.4874	0.5036
	Lane width 12 ft	−0.2159	0.2360	0.0126	−0.6474	0.2825
	Speed limit 25–40 mph	0.3429	0.1853	0.0067	−0.0148	0.7016
	Speed limit 45 mph	0.0710	0.1757	0.0087	−0.2740	0.4136
	Left Shoulder width	−0.0217	0.0254	0.0011	−0.0727	0.0276
	AADT	0.8486	0.0760	0.0063	0.7064	0.9906
	Plain cement concrete Bridge Deck	−0.1082	0.3083	0.0159	−0.7176	0.4813
	Asphalt concrete 7 in thick	0.0157	0.2045	0.0104	−0.3860	0.4203
	Linear trend term	0.1472	0.1426	0.0114	−0.1145	0.4175
	Quadratic trend term	−0.0254	0.0201	0.0016	−0.0639	0.0118
Model 3	Intercept	−7.1500	0.6263	0.0519	−8.1440	−5.7500
	Lane width 9 and 10 ft	0.1459	0.2327	0.0136	−0.2932	0.6201
	Lane width 11 ft	0.1094	0.2459	0.0132	−0.3539	0.6080
	Lane width 12 ft	−0.0883	0.2137	0.0108	−0.5102	0.3393
	Speed limit 25–40 mph	0.3721	0.1820	0.0082	0.0167	0.7380
	Speed limit 45 mph	0.1264	0.1687	0.0079	−0.1999	0.4653
	Left Shoulder width	−0.0336	0.0249	0.0010	−0.0826	0.0159
	AADT	0.9138	0.0736	0.0061	0.7486	1.0340
	Plain cement concrete Bridge Deck	−0.2443	0.3006	0.0144	−0.8394	0.3468
	Asphalt concrete 7 in thick	0.0138	0.1963	0.0090	−0.3756	0.3844
	Gamma	0.54	0.26	0.034	0.03	0.89
	sigma	68670	41260	5425	14750	166600
Model 4	Intercept	−6.6440	1.6520	0.1419	−9.4350	−3.5490
	Lane width 9 and 10 ft	0.5564	0.3060	0.0203	0.0103	1.1950
	Lane width 11 ft	0.2971	0.3034	0.0173	−0.2715	0.9211
	Lane width 12 ft	0.5269	0.3245	0.0211	−0.0917	1.1900
	Speed limit 25–40 mph	0.5439	0.3546	0.0252	−0.2143	1.2060
	Speed limit 45 mph	0.1020	0.3401	0.0248	−0.6385	0.7376
	Left Shoulder width	−0.0803	0.0350	0.0019	−0.1516	−0.0140
	AADT	0.8001	0.2074	0.0178	0.3885	1.1320
	Plain cement concrete Bridge Deck	−0.7369	0.3883	0.0203	−1.5300	−0.0213
	Asphalt concrete 7 in thick	−0.0322	0.2465	0.0140	−0.5198	0.4533

	Variance for αt	66080	118900	7592	33.23	405900

Notes: Bold text refers to statistically significant variables at 95% confidence level.

Table 3 Goodness-of-fit of alternate models.

Temporal Types	Model Types	−D	Pd	DIC
Model 1 (linear)	Base	2121	54	2175
	Mixture	1891	174	2065
Model 2 (quadratic)	Base	2040	94	2134
	Mixture	1888	184	2072
Model 3 (AR)	Base	2122	52	2175
	Mixture	1996	131	2127
Model 4 (time matrix)	Base	1953	90	2044
	Mixture	1883	160	2043

referred to as the log pseudo Bayes factor (LPBF) (Basu and Chib, 2003), where a LPBF greater than 5 reflects the superiority of the model of interest (Ntzoufras, 2009). As shown in Table 4, Model 1 was observed to have the best performance with the highest value of LPML (−862.86), which indicates that consideration of time as a linear covariate offers the best fit with crash data coupled with the benefit of least computational effort. This superior performance by Model 1 was followed by Model 4 (−876.08) and Model 2 (−886.05). The worst performance was exhibited by Model 3 which employed the AR-1 specification with an LPML value of −918.06, which had a very significant LPBF of 55.2 when measured against the best performing Model 1. It was expected that Model 3 would borrow information from the prior one year for estimation of crash risk for the concerned time period but the results seem to demonstrate that such advantage, when used along with 2-component mixture for space-time interaction, may not be transferred to cross-validated predictive performance (CPO), which served as the basis for computation of LPML. In terms of comparison between mixture and Base models, the mixture models consistently demonstrated superior prediction accuracy across all models. The largest discrepancy of 127 points of LPML was observed for Model 1 where the Base model has (−989) while mixture model was superior with (−862). These findings corroborate the previously observed trends in case of model fit, hence indicating that the advantages of mixture models at model fit are transferable to superior performance at prediction of crash counts.

Table 4 Evaluation results of predictive performance of alternate models.

Evaluation Criteria	Comparison Dataset	Model 1		Model 2		Model 3		Model 4
		Mixture	Base	Mixture	Base	Mixture	Base	Mixture	Base
LPML	NA	−862.86	−989.2	−886.65	−932.3	−918.06	−1007	−876.08	−929.45
MAD	Replicated vs. Observed	0.32	0.37	0.33	0.35	0.33	0.37	0.34	0.35
	Predicted vs. Observed	0.39	0.43	0.37	0.39	0.38	0.40	0.38	0.39
	Bayesian vs. Observed	0.32	0.37	0.33	0.35	0.33	0.37	0.33	0.35
MSPE	Replicated vs. Observed	0.26	0.41	0.29	0.38	0.28	0.41	0.31	0.37
	Predicted vs. Observed	0.46	0.52	0.46	0.50	0.47	0.50	0.49	0.51
	Bayesian vs. Observed	0.25	0.40	0.28	0.38	0.29	0.41	0.28	0.37
RSS	Replicated vs. Observed	814.32	4530.6	814.02	1777.71	813.15	7057.45	813.74	1766.99
	Predicted vs. Observed	7334.60	4912.7	9043.77	5685.39	7898.23	6092.34	7056.49	6689.79
	Bayesian vs. Observed	813.74	4548.9	847.92	1769.50	4467.75	5259.25	912.78	1670.22
TRD	Replicated vs. Observed	546652	591886	554405	579092	556374	598923	567821	579691
	Predicted vs. Observed	618411	622451	618623	627074	618563	623069	621886	636765
	Bayesian vs. Observed	546128	591852	554587	579114	559146	598809	555804	578820

1. Bold text represents the best performance for that category.
2. Bold text in Italics represents the best performance for the Base models among that category.
3. To prevent the infinity value obtained for RSS, a small value of 10^−4 was added to the denominator in Eq. (24) for RSS.

4.3.2. Validation based on same data for model development

Apart from the out-of-sample cross-validation evaluation, four criteria (MAD, MSPE, RSS, and TRD) were employed to evaluate the alternate models for predictive accuracy based on the same dataset (insample data) that was utilized for model development. As shown in Table 4, the best performing model at cross-validation predictive capability (Model 1) was also advantageous for nine out of twelve cases, except two cases of RSS and one case of MAD. This phenomenon indicates the high correlation between cross evaluation and the typical evaluations. Comparatively speaking, the other three models had similar performance, but all performed much worse than Model 1 with each of them claiming the leading position in just one or two situations. This finding is somewhat counterintuitive as Model 1 actually has the simplest temporal treatment. The potential reason might be due to the flexible mixture spatiotemporal interaction which can efficiently absorb the residual variability escaping from the predictable part of the model (Abellan et al., 2008). Such fact might help explain why Model 1 even performed worst when normal residual (i.e., Predicted vs. Observed) was used for MAD, while it performs best under all four criteria when random effects were included for Bayesian residual (i.e., Bayesian vs. Observed). In other words, the powerful mixture components might offset the mediocre performance associated with the less sophisticated temporal treatment such as linear one. Overall, the trends discussed above warrant the sensitive consideration of inferences drawn from the evaluation criteria which are employed for assessment of prediction capabilities as the advantages of certain model specifications may not be fully exploited for varied crash datasets. For example, the possible justification for inferior performance of Model 4 may be attributed to the inclusion of a small sample size with respect to the number of years of crash data which eventually resulted in a relatively small matrix (6 × 6) for the temporal CAR specification.

In terms of comparison between the Base and mixture models, the previously discussed goodness-of-fit and predictive accuracy based on the out-of-sample data clearly established the consistent superiority of mixture models. The authors were also interested in seeing if such advantage also extends to the prediction accuracy based on the in-sample crash data. Hence, the Base models were also evaluated for prediction accuracy and the results are summarized in Table 4. The dominance of the mixture models is evident from the observed results as they consistently generated more precisive crash estimates across all temporal models as well as three different datasets.

The authors were also interested in quantifying the influence of the inclusion of random effects on predictive accuracy. The impact of incorporation of random effects to account for spatial and temporal correlation structures, and the flexibility of 2-component mixture may be observed by comparing the Bayesian residuals (Bayesian vs Observed) and normal residuals (Predicted vs Observed), which reflect the predictive performance of random effect models and normal models, respectively. Table 4 demonstrates the significant superiority of sophisticated models which employed random effects as the Bayesian residuals exhibited consistently more accurate predictions of crash counts for all competing models as determined through the four evaluation criteria. The Bayesian residuals exhibited lower value for each evaluation criterion, though the differences were subtle or dramatic depending on the nature of criterion. In the case of MAD, the maximum difference observed was 0.7 for Model 1 while MSPE demonstrated a maximum difference of almost double points in case of Model 1 where normal residuals (0.46) recorded a 0.21 points difference from Bayesian residuals (0.25). A similar difference of 0.21 points was also observed for Model 5 where the Bayesian and normal residuals were relatively high compared to Model 1 and Model 2. The disparity of predictive performance was notably more in the case of TRD, but very significant deviations were recorded for RSS where the Bayesian residuals exhibited eight or nine times lesser values. This consistent superior trend of random effect models reinforces the advantages of incorporation of spatially structured and unstructured correlations to account for the unobserved heterogeneity which benefits the model accuracy at the prediction of crash counts.

5. Conclusions and recommendations

This study introduced the mixture models for space-time interaction in the safety research. The models with the mixture component allow the accommodation of global space-time patterns, which are referred as "stable", while also capturing the atypical departures from the stable pattern which may not be captured by the conventional spatiotemporal models in safety research. This was possible by incorporating two components for the space-time interaction term where one component allowed the accommodation of noise to convey stable pattern while the second component captured the larger fluctuations. The mixture component was explored by collaboration with four different temporal treatments and the models were evaluated based on the goodness-of-fit, and in-sample validation and out-of-sample cross-validation for predictive accuracy. Full Bayesian framework was employed to build Poisson lognormal models where conditional autoregressive specification was employed to incorporate spatially structured random effects using adjacency-based weight matrix, and the space-time interaction term was allowed the flexibility to capture the discrete and continuous space-time variations of the crash data that may have escaped from the predictable patterns of overall temporal and spatial risk surface. This flexible framework served as the base for development of four models which differed on the treatment of temporal aspect: (I) linear time trend; (II) quadratic time trend; (III) Autoregressive-1 (AR-1); and (IV) time adjacency. Moreover, five evaluation criteria were employed for assessment of the predictive capability of alternate models from different perspectives. The Base models were also developed to serve as a reference for demonstrating the advantages of mixture models associated with model fit and prediction capabilities.

Using six-year crash data from various roadway segments of one state route in California, the modeling results demonstrated the robustness of the proposed modeling framework with the vehicle volume being identified as statistically significant for all four models. Plus, Speed limit of 25–40 mph, width of left shoulder and Plain cement concrete Bridge Deck type were also flagged out as highly influential factors of crash occurrence for some of the models. The penalized criterion of DIC was employed for assessment of goodness-of-fit while accounting for the model complexity. Understandably, the mixture models depicted larger complexity associated with the incorporation of the higher number of parameters due to the two-component specification. However, due to the remarkable reduction in posterior deviance, the mixture models subsequently demonstrated lower (compared to corresponding Base models) values of DIC across all four models. The discrepancy of DIC values among Base and mixture models was observed to be largest for the simplest temporal model (linear trend) and experienced constant reduction as the complexity of temporal treatment increased, with the least discrepancy for the time matrix model. The assessment of model fit clearly established the advantages associated with the mixture models. Considering the temporal treatment, the model with time matrix was the best as it seemed to borrow strength from the neighboring years for fitting the crash data.

For the comprehensive comparison of the predictive accuracy of model-estimated crash counts, this study employed an array of evaluation measures based on cross validation and in-sample validation. In the former case, leave-one-our cross-validation (CV-1) technique was utilized to avoid the potential bias of results due to different splits of data. The resultant value of log pseudo marginal likelihood (LPML) illustrated that Model 1 (i.e., the model with the linear time trend treatment) was the best model, whose log pseudo Bayes factor (LPBF) values against competing ones ranged from 4.2 to 55.2. Four different evaluation criteria, including MAD, MSPE, RSS, and TRD, were used for the typical in-sample validation. Under each criterion, the corresponding statistics were obtained by comparing the observed crash counts with three different types of data, respectively, which contain Bayesian estimated counts (or, the predication with random effects), the normal predicted counts (the prediction without random effects), and the model replicated ones. Out of the 12 scenarios, Model 1 again claimed the first place in 9 of them, indicating the somewhat high correlation between cross validation and typical validation results. The remaining three models demonstrated similar performance even though mixed rankings were observed under different conditions. Further examination of results illustrated that Model 1 was placed 2nd and last for two cases when the prediction without random effects was used for evaluation, while its superior performance is consistent across all situations involving prediction with random effects. Such phenomenon indicated that the least sophisticated temporal treatment, linear one, exhibited merely mediocre performance when only the predictable part of the models was considered. However, such limitation was masked when the 2-component-mixtrure spatiotemporal interaction was employed. Indeed, the inclusion of such mixture rendered Model 1 the leading positions in most cases, which might attest to the high efficiency of the mixture components in capturing the unobserved heterogeneity escaping from the risk explained by the covariates and temporal treatment. Moreover, the residual analysis revealed the consistent superior trend of random effect models which bolsters the incorporation of spatially structured and unstructured correlations to account for the unobserved heterogeneity which benefits the model accuracy at the prediction of crash counts. In terms of comparison between Base and mixture models on the basis of predictive accuracy for both in-sample and out-of-sample data, the mixture models demonstrated consistent superiority across all models and crash datasets. This finding supplements the dominance of mixture models observed for model fit and establishes them to be beneficial from different perspectives required in studies of traffic safety (superior fit to crash data and higher precision for crash estimation).

The present study conducts a comprehensive evaluation of four spatiotemporal models with different temporal treatments. However, the results observed in this paper require some caveats. First, the roadway segment crash data are subject to a low-sample mean issue, even though the chosen Poisson-Lognormal model could better handle low sample mean and small sample size as mentioned previously. The data collected from larger spatial units such as traffic analysis zone, census tract, city, or county, are recommended to verify the results obtained hereby. Second, the mixture components are employed for the space-time interactions. The relative performances of models might change when other forms of interactions are used which include linear, quadratic, and time-varying spatial correlation, and so on. Third, only six years of data are used for performance evaluation of various models. It is, therefore, worth exploring longer time periods which might lead to different findings, especially for time adjacency model which is subject to limited sample size in terms of the number of years in the study. Fourth, the comparison is done on the basis of univariate models. The analysis is also recommended for the multivariate settings. Finally, the promising results demonstrated by the mixture models at goodness-offit and predictive accuracy warrant the need for exploring their advantages at another important use of crash prediction models, i.e. site ranking or hot spot identification. Since the mixture component allows the detection of higher fluctuations (variability) in the space-time profile of sites under focus, hence it will be beneficial at identification of sites which deviate from the global space-time pattern observed for all the sites. Such high-risk sites may remain hidden in the traditional spatiotemporal models. The practical applications at the segment level could be the detection of sites with underlying risk which might not be incorporated directly from the data. The possible examples could be some short-term change in traffic behavior such as a musical concert which attracted more activity pertaining to pedestrians, the impact of sun glare during specific months on a section of roadway which may increase crash risk, the construction work for repair of pavement or new construction which may not be evident in data, the issues with line of sight for drivers due to the presence of some obstruction such as a tree, and so on. The mixture model may help highlight such sites with underlying crash risk, which may not be evident in traditional models.

Acknowledgements

Support for this study was provided by 2017 UCCONNECT Research Program (Grant number: 65A0529-052.3). This research was supported, in part, by the Highway Safety Information System (HSIS) through the contribution of data. Opinions expressed in this paper are those of the authors and may not necessarily reflect the views and policies of HSIS. The authors are thankful to the anonymous reviewers for their constructive comments which helped raise the quality of this paper.

References

Abdel-Aty, M., Wang, X., 2006. Crash estimation at signalized intersections along corridors: analyzing spatial effect and identifying significant factors. Transp. Res. Rec. J. Transp. Res. Board 1953, 98–111.

Abellan, J.J., Richardson, S., Best, N., 2008. Use of space-time models to investigate the stability of patterns of disease. Environ. Health Perspect. 116 (8), 1111.
Aguero-Valverde, J., Jovanis, P., 2006. Spatial analysis of fatal and injury crashes in Pennsylvania. Accid. Anal. Prev. 38 (3), 618–625.
Aguero-Valverde, J., Jovanis, P., 2008. Analysis of road crash frequency with spatial models. Transp. Res. Rec. J. Transp. Res. Board 2061, 55–63.
Aguero-Valverde, J., 2013. Full bayes poisson gamma, poisson lognormal, and zero inflated random effects models: comparing the precision of crash frequency estimates. Accid. Anal. Prev. 50, 289–297.
Andrey, J., Yagar, S., 1993. A temporal analysis of rain-related crash risk. Accid. Anal. Prev. 25 (4), 465–472.
Aguero-Valverde, J., Jovanis, P., 2010. Spatial correlation in multilevel crash frequency models: effects of different neighboring structures. Transp. Res. Rec. J. Transp. Res. Board 2165, 21–32.
Basu, S., Chib, S., 2003. Marginal likelihood and Bayes factors for dirichlet process mixture models. J. Am. Stat. Assoc. 98 (461), 224–235.
Blazquez, C.A., Celis, M.S., 2013. A spatial and temporal analysis of child pedestrian
crashes in Santiago, Chile. Accid. Anal. Prev. 50, 304–311. Cheng, W., Washington, S., 2008. New criteria for evaluating methods of identifying hot spots. Transp. Res. Rec. J. Transp. Res. Board 2083, 76–85.
Cheng, W., Gill, G.S., Dasu, R., Xie, M., Jia, X., Zhou, J., 2017a. Comparison of multi variate 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., Sakrani, T., Dasu, M., Zhou, J., 2017b. Predicting motorcycle crash injury severity using weather data and alternative Bayesian multivariate crash fre-
quency models. Accid. Anal. Prev. 108, 172–180. Cheng, W., Gill, G.S., Choi, S., Zhou, J., Jia, X., Xie, M., 2018. Comparative evaluation of temporal correlation treatment in crash frequency modelling. Transportmetrica A: Transp. Sci. 1–40. http://dx.doi.org/10.1080/23249935.2017.1418458. (justaccepted).
Colin Cameron, A., Trivedi, P.K., 1998. Regression Analysis of Count Data 30 ECONOMETRIC SOCIETY MONOGRAPHS.
Gelfand, A., 1996. Model determination using sampling-based methods. In: Gilks, W., Richardson, S., Spiegelhalter, D. (Eds.), Markov Chain Monte Carlo in Practice. Chapman & Hall, Suffolk.
Gill, G., Cheng, W., Xie, M., Vo, T., Jia, X., Zhou, J., 2017. Evaluating the influence of neighboring structures on spatial crash frequency modeling and site ranking performance. Transp. Res. Rec. J. Transp. Res. Board 2659. http://dx.doi.org/10.3141/ 2659-13.
Gourieroux, C.A., Visser, M., 1997. A count data model with unobserved heterogeneity. J. Econometrics 79 (2), 247–268.
Guo, F., Wang, X., Abdel-Aty, M.A., 2010. Modeling signalized intersection safety with
corridor-level spatial correlations. Accid. Anal. Prev. 42 (1), 84–92. Hauer, E., 1997. Observational Before-After Studies in Road Safety. Elsevier Science Ltd.,

Oxford, United Kingdom.

Hay, J.L., Pettitt, A.N., 2001. Bayesian analysis of a time series of counts with covariates: an application to the control of an infectious disease. Biostatistics 2 (4), 433–444. Huang, H., Chin, H., Haque, M., 2009. Empirical evaluation of alternative approaches in

chengbayesianspatiotemporalcrashfrequency2018.md 82 KB Historie Surový