Identifying crash risk factors and high risk locations on an interstate network

Abstract

Highway safety improvement projects are identified by using either (i) a site-specific or (ii) a systemic approach. In the site-specific approach, locations for improvements are ranked according to different performance measures such as critical crash rate, expected crash rate or equivalent property damage only. Alternatively, in the systemic approach, roadway characteristics such as number of lanes, shoulder width, etc. are flagged as a ‘risk’ (or ‘preventative’) feature that increases (decreases) the risk of negative outcomes. Using the Highway Safety Information System database, we seek to merge the two approaches by, first, identifying roadway factors associated with an increased occurrence of car crashes (features we call ‘risk factors’) and, subsequently, identifying roadway segments with a higher crash risk. Specifically, we model the locations of crashes as a realization from a spatial point process. We then parameterize the associated intensity surface of this spatial point process as the sum of a regression on roadway characteristics and spatially correlated error terms. Thus, through the regression piece, we identify hazardous roadway features and through the spatially correlated error terms, we identify locations of high risk.

Keywords

Bayesian Gaussian process Highway Safety Information System Point process spatial statistics

1 Introduction

1.1 Problem statement and data

The goal of highway safety improvement projects is to minimize the risk to users (drivers and pedestrians). Traditionally, improvement projects are identified by using either (i) a site-specific or (ii) a systemic approach. In the site-specific approach, locations for improvements are ranked according to different performance metrics such as critical crash rate, expected crash rate or equivalent property damage only (EPDO; Persaud et al., 1999; Heydecker and Wu, 2001). Alternatively, in the systemic approach, roadway characteristics such as number of lanes, shoulder width, etc. are flagged as a ‘risk’ (or ‘preventative’) feature that increases (decreases) the risk of negative outcomes (Knapp et al., 2014; Lord et al., 2007). Safety improvements are then made across the whole road system to decrease crash risk. Typically, systemic analyses rely on aggregating data to crash counts across similar roadways. However, each road segment has its own unique attributes based on the location of the roadway (e.g., off-road distractions) that may impact its risk level. In this way, spatial location of crashes may act as a confounding factor in systemic analyses.

For this research, we seek to merge the two approaches by both identifying roadway factors associated with an increased occurrence of car crashes (features we call ‘risk factors’) and identifying roadway segments with a higher crash risk. In this way, we seek to inform possible systemic safety improvements to the roadway system while accounting for spatial location which facilitates locating regions where improvements are needed. By combining the two approaches (site-specific and systemic) into one model, we account for the uniqueness in each road as we estimate the effects of the various roadway characteristics. This allows us to estimate the amount of total risk due to location, as well as estimate the risk due to roadway construction.

To accomplish these goals, this analysis seeks to model data from the Highway Safety Information System (HSIS; www.hsisinfo.org); a multistate roadway database overseen by the University of North Carolina Highway Safety Research Center. Currently, states participating in the HSIS include California, Illinois, Maine, Minnesota, North Carolina, Ohio and Washington. However, because of the differences in the variables from state to state (definition and availability), this research focuses on data from five interstate highways in Washington: Interstate 5 (I-5), Interstate 90 (I-90), Interstate 82 (I-82), Interstate 405 (I-405) and Interstate 705 (I-705). Figure 1 shows the locations of these roads (I-705, a 1.5-mile-long road connecting I-5 to Tacoma is not labelled on the map due to its small size). While we choose to focus on interstate highways in Washington here, we note that similar analyses can be performed for other states and roadway types.

Figure 1

Source: Image obtained from geology.com

The HSIS data include roadway inventory information (number of lanes, type of pavement, speed limit, traffic volume, etc.) and crash information (crash type, vehicle type, number of vehicles involved, etc.) for crashes that occur between the years 2010–2012. Because roadway characteristics vary slightly from year to year, we used roadway inventory data from 2012 to keep the roadway definitions temporally consistent. Given the short time span considered here (three years), the 2012 roadway characteristics are likely representative of the state of the roadway. Table 1 displays the variables used in this study, along with the associated HSIS database where the information was obtained.

Table 1

Variables used in this study, and the associated HSIS database where the data was obtained

Variable Name	Description	Database
rd_inv	Road ID	Accident, road, curve
rte_nbr	Route number	Accident, road, curve
milepost	Milepost where crash occurred	Accident
lshl_typ	Type of left shoulder	Road
rshl_typ	Type of right shoulder	Road
med_type	Type of median	Road
terrain	Terrain of roadway	Road
spd_limt	Speed limit	Road
lshldwid	Left shoulder width (in feet)	Road
rshldwid	Right shoulder width (in feet)	Road
medwid	Median width (in feet)	Road
no_lanes	Number of lanes	Road
lanewid	Lane width (in feet)	Road
aadt	Average annual daily traffic	Road
deg_curv	Degree of curvature	Curve

The response variable in the HSIS data considered here is the location of a crash, which we denote by $s_{i} = (r_{i}, m_{i})^{'}$ for $i = 1, \dots, N$ , where $N =$ 33 589 is the total observed number of crashes on the five interstates, $r_{i}$ is the roadway (interstate) and $m_{i}$ is the mile marker. Of the five interstates, I-5 had the most crashes (21 089 crashes on 276.62 miles of highway), followed by I-90 (5 979 crashes on 297.51 miles of highway), I-405 (5 156 crashes on 30.3 miles of highway), I-82 (1 269 crashes on 132.57 miles of highway) and I-705 (96 crashes on 1.5 miles of highway). Figure 2 shows the density of crashes in our data along the interstates of interest. From these plots, we see that crash rates vary along the interstates, but tend to be highest near large cities. It is worth noting that both I-90 and I-405 have areas of relatively high crash occurrence which do not coincide with the locations of large cities.

Figure 2

Kernel density estimates of crashes along each interstate highway considered in our analysis. Tick marks along the horizontal axis indicate exact crash locations. Vertical lines indicate the approximate centres of large cities along the interstate

1.2 Modelling challenges and model overview

Statistically speaking, the goal here is to, first, regress the road factors from Table 1 onto the observed crash locations to determine which characteristics are associated with the high levels (density) of crashes. Second, after accounting for the roadway factors, we seek to identify road segments that have a high risk for reasons not attributable to the roadway factors in the dataset. Thus, by accomplishing these two goals, we will not only identify possible systemic improvements to the highway system (risk factors) but also specific locations that could benefit from safety improvements.

There are two main challenges associated with these goals. First, as mentioned previously, the response variable is the location of the crash (interstate and mile marker). While regression is a simple goal, standard regression techniques such as linear or generalized linear models (Faraway, 2014, 2016) are not applicable for this type of response. One approach would be to aggregate (count) the number of crashes along predefined road segments. In doing so, models for count data such as zero-inflated Poisson (ZIP), negative binomial models or Poisson models may then be appropriate (Lord et al., 2005, 2007; El-Basyouny et al., 2014; Chiou and Fu, 2013), which can be implemented using either a frequentist or Bayesian paradigm (Schultz et al., 2011; Sacchi et al., 2015). Concerns with aggregating crash counts along predefined roadway segments include, first, the defined segments may be arbitrarily selected, second, information at the sub-segment level is forfeit and, third, it converts an inherent spatially continuous process into a discrete one.

Rather than aggregrating as suggested previously, we opt to treat the locations $s_{i}$ for $i = 1, \dots, N$ as a realization from a non-homogenous spatial point process (Illian et al., 2008; Baddeley et al., 2015). A spatial point process is a random process over a spatial domain $D \subset ℝ^{d}$ , where a realization from this process consists of, first, the number of events $N$ and, second, a finite set of locations $s_{i}, i = 1, \dots, N$ within that space. The number $N$ and realized locations $s_{i}, \dots, s_{N}$ are governed by an intensity surface $Λ (s) : D \to ℝ^{+}$ in the following manner. For any subset $B \subseteq D$ , the number of events within $B$ is distributed as $P (δ)$ , where $P (δ)$ indicates a Poisson distribution with mean $δ = \int_{B} Λ (s) d s$ . Intuitively, for regions $B$ , where $Λ (s)$ is high, more crashes are expected to occur and vice versa. Conditional on the number of events $N$ , the locations $s_{1}, \dots, s_{N}$ are then drawn independently from the normalized intensity $Λ (s) / \int_{D} Λ (u) d u$ .

Our research approach is to model $Λ (s)$ as the sum of three components. The first component is a regression term on the roadway characteristics such that positive regression coefficients will indicate which road characteristics are associated with areas of high expected crash counts (satisfying our first goal). The second goal (identifying road segments with additional risk not associated with the roadway factors) will be accomplished by the second two components: spatially correlated and white noise error terms. These error terms will be interpretable as additional risk after accounting for traffic and the road characteristics. Hence, regions with high errors terms identify segments of interest.

A second challenge in analysing the HSIS data is that the locations of crashes on road segments within a roadway network are spatially correlated. In other words, two locations on an interstate that are close together should have more similar crash rates than segments that are far from each other. In this spatial point process framework used here, this challenge translates to the need to model a smooth intensity surface in space. Work by Barua et al. (2015), Ahmed et al. (2011), Guo et al. (2010), Aguero-Valverde and Jovanis (2008, 2010) and Barua et al. (2014) suggest that including spatial correlation in the statistical model provides superior fit while providing unbiased estimates of model parameters. While methods for incorporating spatial correlation (smoothness) are well established (El-Basyouny et al., 2014), traditional spatial methods generally utilize straight-line distance to estimate the correlation between two points. Straight-line distance is certainly inappropriate for roadway networks because these networks do not allow for straight-line travel between any two points, but rather require travellers to use existing roads to get from point A to point B. As a solution, we incorporate spatial correlation into our model based on driving distance, rather than Euclidean distance.

To summarize, the primary contribution of this article is to simultaneously perform a site-specific and systemic analysis of the interstate network in Washington shown in Figure 1. In doing so, this article proposes a spatial point process model with appropriate regression terms and spatially correlated error structure. The remainder of this article is outlined as follows. Section 2 details the point process model and reviews relevant methodology from the literature. Section 3 presents results and interpretations, as well as assessing the model fit. Section 4 details our conclusions and plans for future work.

2 Model

2.1 Model details

As in Section 1.2, let $s_{i} = (r_{i}, m_{i}) \in D$ for $i = 1, \dots, N$ denote the spatial location of the $i^{th}$ crash, where $r_{i} \in {I - 5, I - 90, I - 82, I - 405, I - 705}$ is the interstate, $m_{i}$ is the mile marker and $D$ is the spatial domain of interest (the space of interstate and mile markers). We assume that,

s_{1}, . . ., s_{N}, N \sim PP (Λ (s)),

(2.1)

where

PP (\cdot)

indicates a spatial point process and

Λ (s) : D \to ℝ^{+}

is an intensity surface. By the assumption in equation (2.1),

N \sim P (\int_{D} Λ (s) d s)

, where

P (μ)

denotes the Poisson distribution with mean

μ

, and

s_{i}

are independent and identically distributed as

Λ (s) / \int_{D} Λ (u) d u

. We emphasize that the model in (2.1) assumes each crash is independent (i.e., the presence of one crash does not affect the presence of another crash). Strictly speaking, this assumption may be violated in that the presence of one crash could create conditions that lead to a second crash (a process referred to as an attractive point process; see Goldstein et al., 2015). However, the HSIS data was thinned prior to analysis to ensure that each crash is unique.

In the point process framework, the intensity surface $Λ (s)$ is the quantity (parameter) of interest as this governs not only the locations of crashes but also the number of crashes. In essence, $Λ (s)$ represents a spatial ‘map’ over $D$ , where large $Λ (s)$ correspond to regions with a high rate of crashes and vice versa. Following Banerjee et al. (2014), the associated log-likelihood of the observed data $s_{1}, \dots, s_{N}$ and $N$ used to estimate $Λ (\cdot)$ is given by,

L (Λ) = - \int_{D} Λ (s) d s + \sum_{i = 1}^{N} \log (Λ (s_{i})) .

(2.2)

The use of equation (2.2) for statistical inference poses a few interesting challenges. First, recall that

Λ (s)

is an unknown (and infinite dimensional) parameter such that

\int_{D} Λ (s) d s

is an unknown integral. Solutions to calculating this integral include using a discrete Riemann approximation (Liang et al., 2008; Chakraborty et al., 2011; Waagepetersen, 2004) or flexible density functions (Kottas and Sansó, 2007; Heaton et al., 2015). In either solution, computation can be slow, so here we use a computationally simpler alternative based on a regionally constant specification for

Λ (s)

Let $R_{1}, \dots, R_{K}$ denote a partition of the spatial domain $D$ into $K$ road segments such that $R_{k_{1}} ⋂ R_{k_{2}} = \emptyset$ for $k_{1} \neq k_{2}$ . We parameterize,

Λ (s) = \sum_{k = 1}^{K} \frac{λ_{k}}{| R_{k} |} 1 {s \in R_{k}},

(2.3)

where

R_{k}

indicates the

k^{th}

road segment and

| R_{k} |

indicates the length of

R_{k}

. Under this parameterization notice that

λ_{k} = \int_{R_{k}} Λ (s) d s

is the expected number of crashes that occur in region

R_{k}

. Further notice that equation (2.2) becomes

\begin{array}{l} L (Λ) = - \sum_{k = 1}^{K} λ_{k} + \sum_{i = 1}^{N} \log (λ_{{k : s_{i} ϵ R_{k}}}) \\ = - \sum_{k = 1}^{K} [λ_{k} + N_{k} \log (λ_{k})], \end{array}

(2.4)

where

N_{k}

is the number of crashes that occur in region

k

and

\sum_{k = 1}^{K} N_{k} = N

. Notably, (2.4) is void of integrals and is much easier to use computationally.

To intuitively understand the parameterization in (2.3), consider Figure 3 which shows a one dimensional intensity surface (solid curve) along with two resolutions of regions (two values of $K$ ) indicated by the vertical dotted lines. The horizontal line segments indicate the maximum likelihood estimate of the constant height of the intensity surface (the $λ_{k}$ ) within each region (estimated independently). As can be seen from the Figure, as the number of regions $K$ increases, the discrete parameterization of $Λ (s)$ , indicated by the horizontal line segments, becomes closer and closer to the continuous $Λ (s)$ . Figure 3

Two discrete approximations of continuous intensity surface (solid curve) using two different values for K. Horizontal line segments indicate maximum likelihood estimates of the surface and dashed lines indicate region boundaries

Intuitively, more crashes are expected in areas with higher traffic. However, since our first goal (see Section 1.2) is to identify roadway characteristics associated with a higher rate of crashes and safety improvement projects are not able to alter traffic levels, we factor

λ_{k} = E_{k} μ_{k},

(2.5)

where

E_{k} > 0

, the expected number of crashes due only to traffic, and

μ_{k} > 0

is the relative risk. Specifially, we use internal standardization (Waller et al., 1997) to define

E_{k} = \frac{N}{\sum_{j = 1}^{K} {AADT}_{j}} {AADT}_{k},

(2.6)

where

{AADT}_{k}

is the Average Annual Daily Traffic (AADT) for segment

k

(the total volume of vehicle traffic over a road for a year divided by the number of days in a year) and

N

is the total number of observed crashes in our data. Inspecting (2.6), note that the first term is the number of crashes per AADT which when multiplied by

{AADT}_{k}

gives the expected number of crashes on segment

k

if traffic volume were the only factor that determined crash locations.

The $μ_{k}$ parameter in (2.5) represents increased or decreased risk of segment $k$ relative to the expected number of crashes given the traffic levels ( $E_{k})$ . Specifically, if $0 < μ_{k} < 1$ , then segment $k$ had fewer crashes than expected given traffic levels (lower risk). In contrast, if $μ_{k} > 1$ , segment $k$ had more crashes than expected given traffic levels (higher risk). Hence, performing inference for $μ_{k}$ provides information regarding high- and low-risk areas.

At this point, recall that our two goals are to, first, identify roadway features associated with high risk areas and, second, to identify ‘hot-spots’ where the crash risk is high after accounting for traffic and roadway features. These goals can be accomplished by regressing roadway features on $μ_{k}$ to determine which road factors are associated with high- or low-road risk. Let,

\log (μ_{k}) = x_{k}^{'} β + w_{k} + ε_{k},

(2.7)

where

x_{k} = (1, x_{k 1}, \dots, x_{kP})^{'}

is a vector of

P

roadway features taken from Table 1 with associated coefficients

β = (β_{0}, β_{1}, \dots, β_{P})^{'}

w_{k}

is a zero mean, spatially correlated random effect for region

k

that accounts for spatially correlated homogeneity between road segments and

ε_{k} \sim N (0, τ^{2})

are independent, uncorrelated, mean zero Gaussian errors that account for small-scale heterogeneity between roadway factors.

The coefficients $β_{1}, \dots, β_{P}$ in (2.7) determine which roadway factors are classified as ‘risk’ factors and which are classified as ‘preventative’. That is, if $β_{j} > 0$ , then the $j^{th}$ covariate can be viewed as a risk factor because increases in the covariate correspond to increases in the relative risk. Conversely, if $β_{j} < 0$ , then the $j^{th}$ covariate can be viewed as ‘preventative’ because increases in the covariate correspond to decreases in the relative risk.

The mean zero error terms $w_{k}$ and $ε_{k}$ represent additional risk beyond that explained by the roadway characteristics. In contrast to $ε_{k}$ , we model $w_{k}$ as additional risk attributable to similarities across roadway segments. Hence, we assume the spatial random effects $w = (w_{1}, \dots, w_{K})^{'}$ are correlated in space such that $w \sim N (0, σ^{2} M_{α})$ , where $M_{α} = {M (d (k, k^{'}) ∣ α)}_{k, k^{'} = 1}^{K}$ is the $K \times K$ matrix of Matern spatial correlations with decay $α$ and smoothness $ν = 1 / 2$ (Banerjee et al., 2014). Recall that, since travel within a roadway network can only be accomplished within the network's existing framework, $d (k, k^{'})$ indicates the driving distance between segments $k$ and $k^{'}$ , rather than Euclidean distance. The $α$ parameter in the correlation structure controls the extent of the spatial correlation (i.e., the rate at which the correlation decreases to zero as a function of distance). As an example, consider two points at a fixed distance of $d$ driving miles apart. With a fixed $d$ , as $α$ increases, the correlation decreases. In contrast, the smoothness parameter $ν$ controls the smoothness of the process as a function of distance. For example, if $ν$ is large, then the corresponding spatial process varies smoothly in space (i.e., is highly differentiable), whereas if $ν$ is small, then the process is not smooth (i.e., perhaps not differentiable). As is common practice, we choose to fix $ν = 1 / 2$ for this analysis.

Marginalizing equation (2.7) over $w \sim N (0, σ^{2} M_{α, ν})$ and $ε \sim N (0, τ^{2} I)$ , we have

\log (μ) \sim N (X β, Σ),

(2.8)

where

μ

is the vector of relative risks,

X

is a

K \times (P + 1)

matrix of roadway characteristics with associated coefficients,

β

and

Σ = σ^{2} M_{α, ν} + τ^{2} I

is the total variance of the spatial and non-spatial random effects. For emphasis, note that the model in (2.8) accomplishes our first research goal through

β

which quantifies the relationship between the road factors and the risk of the segment. The second goal is met when we estimate

w_{k}

, the spatial random effects. Large positive

w_{k}

’s indicate segments with a large amount of risk not attributable to traffic or the road factors in our model, and these segments can be flagged as ‘hot spots’ and in need of improvement.

2.2 Model fitting and parameter estimation

We opted to use the Bayesian paradigm (Gelman et al., 2014; Carlin and Louis, 2008) for parameter estimation, although we do note that parameter estimates can also be obtained via maximum likelihood. As prior distributions, we generally assumed vague priors, with a few exceptions which we detail here. First, we assumed $β \sim N (0, 10 I)$ , where a variance of 10 is quite vague considering these are effects on $\log (μ_{k})$ .

Theoretical work by Zhang (2004), Kaufman and Shaby (2013) show that estimation for $α$ can be difficult due to a flat likelihood surface. For this reason, we adopt a discrete uniform prior distribution for $α$ . Specifically, we assume that, $α \in {0.60, 0.30, 0.20, 0.15, 0.12}$ with equal prior probability. These $α$ values produced correlation structures which dropped to a correlation of 0.05 at a distance of 5, 10, 15, 20 and 25 miles, respectively.

In order to estimate $σ^{2}$ and $τ^{2}$ , we reparameterized $Σ$ to $Σ = γ^{2} (ω M_{α} + (1 - ω) I)$ such that $γ^{2} = σ^{2} + τ^{2}$ and $ω = σ^{2} / (σ^{2} + τ^{2})$ such that $γ^{2} > 0$ is the total variance and $ω \in [0, 1]$ is the proportion of the total variance attributable to spatial variation. Backtransforming, we have $σ^{2} = γ^{2} ω$ and $τ^{2} = γ^{2} (1 - ω)$ such that prior distributions for $γ^{2}$ and $ω$ induce a prior distribution for $σ^{2}$ and $τ^{2}$ . Because $0 \leq ω < 1$ , we again used a discrete uniform prior distribution on $ω$ for 20 values evenly spaced between 0 and 0.900. Additionally, we assumed that $γ^{2}$ a priori followed an inverse-gamma distribution with shape 2.1 and rate 1.

3 Results

Using the draws from the posterior distribution, Figure 4 compares the observed crash counts on each road segment to the expected crash counts (defined as the posterior mean of $λ_{k}$ in equation (2.5)) on these same segments. Because the association between the observed and expected counts is essentially linear (correlation of 0.99), our model fits the data well, thereby validating our results.

Figure 4

Comparison between actual crash counts for each 1-mile road segment in our data and expected crash counts for these same segments generated by our model. The extremely high correlation between the counts shows that our model fits our data.

3.1 Systemic results

Through our MCMC algorithm, we obtained draws from the posterior distributions for the effects of each of the roadway characteristics (i.e., $β$ in equation (2.7)). From these draws, we estimated the mean effect of each roadway factor, a central 95% credible interval and the probability of each parameter being greater than or less than zero (depending on which probability was higher), which are displayed in Table 2. The probabilities $\Pr (| β | > 0 ∣ data)$ provide information on how significant a road characteristic is.

Table 2

The mean and 95% credible interval for the effect of each roadway characteristic and the probability of the magnitude of the effect being greater than 0. The type of variable is either quantitative (Q) or categorical (C)

	Type	Posterior
Variable	Q/C	Mean	Lower	Upper	$\Pr (\| β \| > 0)$
Average curvature	Q	0.16	0.08	0.23	1.00
Left shoulder not asphalt	C	0.24	0.02	0.46	0.98
Rolling terrain	C	0.25	0.01	0.48	0.98
Mountainous terrain	C	0.42	0.01	0.83	0.98
Median not soil	C	−0.25	−0.47	−0.00	0.98
Lane width $>$ 12 ft	C	0.13	0.97	0.97	0.97
6 lanes	C	−0.15	−0.32	0.04	0.94
$\leq$ 60 MPH speed limit	C	0.14	−0.07	0.36	0.90
65 mph speed limit	C	0.29	−0.22	0.80	0.87
$>$ 6 ft left shoulder	C	0.07	−0.09	0.22	0.82
$\geq$ 10 ft right shoulder	C	−0.05	0.75	0.75	0.75
Right shoulder not asphalt	C	0.06	−0.14	0.26	0.74
$\geq$ 7 lanes	C	−0.05	−0.24	0.16	0.68
Median width	Q	−0.00	−0.00	0.00	0.67
Urban	C	−0.01	−0.18	0.17	0.54
4–6 ft left shoulder	C	−0.24	−5.33	4.87	0.54
2–4 ft left shoulder	C	−0.25	−5.26	4.75	0.54

Consider, first, the quantitative variable ‘average curvature’ as an example. From Table 2, as the average curvature of a road segment increases by one degree, we would expect the log relative risk of the segment to increase by 0.16 (0.08, 0.23), on average. In other words, segments with higher average curvature tend to have a higher relative risk, or have more crashes than we expect given traffic volume. Conversely, because the posterior mean of the effect of median width is −9.3 $\times 10^{- 5}$ (−0.000 5, 0.000 3), segments with wider medians tend to have a lower relative risk, on average (fewer crashes than expected given traffic volume). However, since the 95% credible interval for this effect contains zero, this effect is likely not significant.

The effects of the categorical variables speed limit, terrain type, median type, shoulder type, shoulder width, number of lanes, lane width and urban vs rural area are relative to the baseline categories which we chose based on the most common category observed in the data. Specifically, the baseline categories are a speed limit of 70 MPH, level terrain, soil median, asphalt shoulders, 0–2 ft left shoulder, smaller than 10 ft right shoulder, 5 or fewer lanes, 12 ft lane width and rural area type. As an example in how to interpret these effects, consider the effect of rolling terrain. Because the effect of rolling terrain is 0.25 (0.01, 0.48), road segments with rolling terrain tend to have a higher relative risk than those with level terrain. Conversely, because the effect of a median pavement type that is not soil is −0.25 (−0.47, −0.003), segments with non-soil medians (e.g., concrete or asphalt medians) tend to have a lower relative risk than those with soil medians.

Based on the results in Table 2, in terms of roadway construction, our recommendations would be to build straight, flat roads with asphalt left shoulders, medians of types other than soil (e.g., concrete medians), and lanes of 12 ft in width. Initially, it may seem surprising that speed limits lower than 70 MPH were estimated to increase the crash risk of road segment; however, we speculate that the high coefficients for the 65 and $\leq$ 60 MPH speed limits are misleading. That is, we hypothesize that the high speed limits themselves are not what is associated with lower risk, but rather other characteristics confounded with high speed limits which were not included in our data, such as pavement type, more specific gradient information or other environmental factors (trees, cliffs, risk of falling rocks). The standard lane width for road construction in the United States is 12 feet. According to the FHWA website, the only types of roadways constructed with widths larger than 12 ft are freeway ramps. This was verified for our data by checking road segments specified to have larger average lane widths than 12 ft against maps of the road segments. Therefore, the high risk factor for lane widths greater than 12 feet actually indicate that entrance and exit ramps are a potential risk factor for crashes.

3.2 Site-specific results

As discussed in the Introduction, we seek to merge the site-specific and systemic approaches to highway safety analysis. The coefficients discussed previously highlight potential systemic improvements to the roadways. However, as discussed in Section 2, we expect that not all variability in the relative risk is explained by roadway characteristics; hence, we included the spatially correlated and independent error terms into our model specification ( $w_{k}$ and $ε_{k}$ in equation (2.7)). Let ${\tilde{\in}}_{k} = w_{k} + \in_{k}$ denote the residual associated with road segment $k$ . These residuals, then, represent the increased (or decreased, if negative) risk associated with segment $k$ after having accounted for the roadway characteristics in our database. The posterior mean of these residuals for I-5 and I-90 are shown in Figure 5.

Figure 5

This plot shows additional risk (not attributable to traffic levels or the roadway factors in our data) identified by our model for I-5 (left) and I-90 (right). The solid curves represent a kernel smoother of the residuals to highlight patterns

Using the spatial residuals along I-5 in the left panel of Figure 5 as an example, we can see evidence of spatial clustering which, subsequently, justifies our use of spatial modelling. For example, notice that the residuals near Tacoma and Seattle are clustered above 0, indicating spatially related elevated risk for these areas that is not accounted for by our roadway characteristics. If not roadway characteristics, what contributes to the elevated risk in these areas? One possibility that we hypothesize accounts for these high residuals is driver distraction. Figure 6 displays a Google map image of miles 165–166 and 166–167 on I-5. Notice that there is a clear view of the Space Needle, a well-known Seattle landmark, along these road segments. This ‘off-road’ distraction, as well as other features in the major downtown area, could be a possible contributor to this added risk.

Figure 6

The view from mile 166 on I-5, obtained from Google Street View. The Space Needle, as well as other Seattle landmarks, is clearly visible from the road, which may distract drivers, possibly contributing to the increased risk in this area

Examining the right panel in Figure 5, there is a clear increase in risk for Seattle, a smaller increase near Spokane, and also a clear increase in additional risk from mile 50 to approximately mile 100. As displayed in Figure 7, a winter recreation resort, the Summit at Snoqualmie, is located at mile 51–53 and the Snoqualmie Pass, an occasionally hazardous mountain pass, is located at mile 54. The area is heavily forested and also includes other mountain recreation areas. We hypothesize that these factors may be associated with interstate entrance and exit ramps as well as further distractions for drivers. Furthermore, the highest peak in the residuals, within the mile 50 to mile 100 region of I-90, occurs at mile 86–87. This location is where I-90 crosses the Yakima River, another potential distraction to drivers.

Figure 7

Google Maps image of miles 50–54 on I-90, an area identified by our model as a region of high risk not attributable to the roadway factors in our data. A ski resort, the Summit at Snoqualmie, as well as other recreational areas, are located here suggesting that driver distraction and frequent interstate entrances and exits may contribute to risk

Finally, beyond systemic and site-specific approaches individually, our modelling strategy has the ability to merge these two approaches to estimate the overall risk $μ_{k}$ . Figure 8 displays the posterior means of $μ_{k}$ also for I-5 and I-90 and represent the cumulative crash risk along these roads from the measured roadway characteristics and underlying residuals.

Figure 8

The posterior mean of the relative risk $μ$ over I-5 (left) and I-90 (right). If $μ < 1$ , there were fewer crashes than expected due to traffic levels and if $μ > 1$ there were more crashes than expected

Looking at Figure 8, recall that $μ > 1$ suggests that there were more crashes on the segment than expected given traffic levels and that $μ < 1$ suggests that there were fewer crashes than expected given traffic levels. Comparing Figure 8 to Figure 5, often when a segment has a large spatial residual, it also has a large relative risk. For example, segments from mile 165–166 and 166–167 on I-5, which we looked at before as having large spatial residuals have $μ_{k}$ values of 4.39 (4.06, 4.74) and 4.10 (3.81, 4.39) respectively suggesting about 4 times more crashes occur than expected on each of these segments. Mile 52–53 on I-90 (by the Snoqualmie Resort) has a $μ_{k}$ of 3.03 (2.49, 3.65) and mile 86–87 (by the Yakima River) on I-90 has a $μ_{k}$ of 2.24 (1.76, 2.79). These $μ_{k}$ values identify specific segments worth looking into further. In addition, when we combine our knowledge of the $μ$ ’s and the spatial residuals, we can determine what proportion of the additional risk is due to roadway information in our data and what is due to spatial information.

4 Conclusion

In this article, we implemented a framework to use spatial statistical methods over a roadway network in order to identify factors associated with high crash counts (a systemic approach), as well as to pinpoint road segments of interest which had higher crash rates than expected (a site-specific approach). To do this, we assumed crash locations followed a non-homogeneous Poisson process. The intensity surface of the Poisson process was subsequently modelled as the sum of a regression piece (to estimate the impact of roadway features) as well as spatial and non-spatial residuals. Systemically, our results indicated that increasing the shoulder width, decreasing curvature and including a non-soil median (i.e., concrete barriers) decreases the overall crash risk. Likewise, our site-specific analysis suggests that distractions and frequent entrances and exits from the interstate contribute to increased crash risk. While our results are informative for highway safety improvement projects, we recommend that more research be done to determine the aspects of major metropolitan and recreational areas which increase the riskiness of roads.

We emphasize that our results here are associations rather than causal relationships. We do not have information on what caused these specific crashes (i.e., from police records). Rather, we merged roadway databases with crash databases to investigate trends in systemic features and roadway segments with crash locations. An interesting area of future research would be to investigate the use of causal methods for highway safety analyses and contrast our results with those of a causal analysis (Karwa et al., 2011).

For this analysis, we considered all crashes. However, we hypothesize that roadway features correlate differently for various types of crashes. For example, adding a concrete barrier could decrease the number of ‘run off the road’ crashes yet increase the number of ‘multiple vehicle crashes’ because cars, possibly, would veer into other lanes to avoid hitting the barrier. In the future, we plan to extend our methodology to account for the types of crashes in a multivariate framework. Additionally, this research only considered an interstate network. We are also interested in using this methodology to examine more complex roadway networks, possibly by including highways as well as interstates.

Another possible direction for future research would involve incorporating crash severities into the relative risk of road segments. For example, road segments with lower speed limits may observe a high number of total crashes but have a low average severity, while segments with higher speed limits may observe fewer total crashes but with a higher average severity. Incorporating the crash severities into the analysis would likely give a better idea of the true risk of the segments.

Since our dataset is quite large (observations on 33 589 crashes on 741 road segments), one last possibility we suggest for further research would be to use machine learning methods, such as unsupervised classification, to examine the relationships between the roadway factors.

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant Number DMS-1417856. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

References

Aguero-Valverde

Jovanis

(2008) Analysis of road crash frequency with spatial models. Transportation Research Record: Journal of the Transportation Research Board , 55–63.

Aguero-Valverde

Jovanis

(2010) Spatial correlation in multilevel crash frequency models: Effects of different neighboring structures. Transportation Research Record: Journal of the Transportation Research Board 21–32.

Ahmed

Huang

Abdel-Aty

Guevara

(2011) Exploring a Bayesian hierarchical approach for developing safety performance functions for a mountainous freeway. Accident Analysis & Prevention , 43, 1581–89.

Baddeley

Rubak

Turner

(2015) Spatial Point Patterns: Methodology and Applications with R . Boca Raton, FL: CRC Press.

Banerjee

Carlin

Gelfand

(2014) Hierarchical Modeling and Analysis for Spatial Data . CRC Press.

Barua

El-Basyouny

Islam

(2014) A full Bayesian multivariate count data model of collision severity with spatial correlation. Analytic Methods in Accident Research , 3, 28–43.

Barua

El-Basyouny

Islam

(2015) Effects of spatial correlation in random parameters collision count-data models. Analytic Methods in Accident Research , 5, 28–42.

Carlin

Louis

(2008) Bayesian Methods for Data Analysis . CRC Press.

Casella

George

(1992) Explaining the Gibbs sampler. The American Statistician 46, 167–74.

10.

Chakraborty

Gelfand

Wilson

Latimer

Silander

(2011) Point pattern modelling for degraded presence-only data over large regions. Journal of the Royal Statistical Society, Series C (Applied Statistics) , 60, 757–76.

11.

Chib

Greenberg

(1995) Understanding the Metropolis–Hastings algorithm. The American Statistician , 49, 327–35.

12.

Chiou

Y-C

(2013) Modeling crash frequency and severity using multinomial-generalized poisson model with error components. Accident Analysis & Prevention , 50, 73–82.

13.

El-Basyouny

Barua

Islam

(2014) Investigation of time and weather effects on crash types using full Bayesian multivariate Poisson lognormal models. Accident Analysis & Prevention , 73, 91–99.

14.

Faraway

(2014) Linear models with R . CRC Press.

15.

Faraway

(2016) Extending the Linear Model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models ., Vol. 124. CRC press.

16.

Gelman

Carlin

Stern

Rubin

(2014) Bayesian Data Analysis . Vol. 2. Boca Raton, FL: Chapman & Hall/CRC.

17.

Goldstein

Haran

Simeonov

Fricks

Chiaromonte

(2015) An attraction-repulsion point process model for respiratory syncytial virus infections. Biometrics , 71, 376–85.

18.

Guo

Wang

Abdel-Aty

(2010) Modeling signalized intersection safety with corridor-level spatial correlations. Accident Analysis & Prevention , 42, 84–92.

19.

Haario

Saksman

Tamminen

(2001) An adaptive metropolis algorithm. Bernoulli , 223–42.

20.

Heaton

Sain

Monaghan

Wilhelmi

Hayden

(2015) An analysis of an incomplete marked point pattern of heat-related 911 calls. Journal of the American Statistical Association , 110, 123–35.

21.

Heydecker

(2001) Identification of sites for road accident remedial work by Bayesian statistical methods: An example of uncertain inference. Advances in Engineering Software , 32, 859–69.

22.

Illian

Penttinen

Stoyan

(2008) Statistical Analysis and Modelling of Spatial Point Patterns . Vol. 70. John Wiley & Sons.

23.

Karwa

Slavkovic

Donnell

(2011). Causal inference in transportation safety studies: Comparison of potential outcomes and causal diagrams. The Annals of Applied Statistics , 1428–55.

24.

Kaufman

Shaby

(2013) The role of the range parameter for estimation and prediction in geostatistics. Biometrika , 100, 473–84.

25.

Knapp

Hallmark

Bou-Saab

(2014) Systemic safety improvement risk factor evaluation and countermeasure summary (Technical report).

26.

Kottas

Sanso

(2007) Bayesian mixture modeling for spatial poisson process intensities, with applications to extreme value analysis. Journal of Statistical Planning and Inference , 137, 3151–63.

27.

Liang

Carlin

Gelfand

(2008) Analysis of Minnesota colon and rectum cancer point patterns with spatial and nonspatial covariate information. The Annals of Applied Statistics , 3, 943.

28.

Lord

van Schalkwyk

Chrysler

Staplin

(2007) A strategy to reduce older driver injuries at intersections using more accommodating roundabout design practices. Accident Analysis & Prevention , 39, 427–32.

29.

Lord

Washington

Ivan

(2005) Poisson, Poisson-gamma and zero-inflated regression models of motor vehicle crashes: Balancing statistical fit and theory. Accident Analysis & Prevention , 37, 35–46.

30.

Lord

Washington

Ivan

(2007) Further notes on the application of zero-inated models in highway safety. Accident Analysis & Prevention , 39, 53–7.

31.

Persaud

Lyon

Nguyen

(1999) Empirical Bayes procedure for ranking sites for safety investigation by potential for safety improvement. Transportation Research Record: Journal of the Transportation Research Board , 7–12.

32.

Sacchi

Sayed

El-Basyouny

(2015) Multivariate full Bayesian hot spot identication and ranking: New technique. Transportation Research Record: Journal of the Transportation Research Board , 1–9.

33.

Schultz

Thurgood

Olsen

Reese

(2011) Analyzing raised median safety impacts using Bayesian methods. Transportation Research Record: Journal of the Transportation Research Board , 96–103.

34.

Spiegelhalter

Best

Carlin

Van Der Linde

(2002) Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 64, 583–639.

35.

Waagepetersen

(2004) Convergence of posteriors for discretized log Gaussian cox processes. Statistics & Probability Letters , 66, 229–35.

36.

Waller

Carlin

Xia

Gelfand

(1997) Hierarchical spatio-temporal mapping of disease rates. Journal of the American Statistical Association , 92, 607–17.

37.

Zhang

(2004) Inconsistent estimation and a symptotically equal interpolations in model-based geostatistics. Journal of the American Statistical Association , 99, 250–61.

Identifying crash risk factors and high risk locations on an interstate network

Abstract

Keywords

1 Introduction

1.1 Problem statement and data

Variables used in this study, and the associated HSIS database where the data was obtained

Kernel density estimates of crashes along each interstate highway considered in our analysis. Tick marks along the horizontal axis indicate exact crash locations. Vertical lines indicate the approximate centres of large cities along the interstate

2 Model

2.1 Model details

Two discrete approximations of continuous intensity surface (solid curve) using two different values for K. Horizontal line segments indicate maximum likelihood estimates of the surface and dashed lines indicate region boundaries

3 Results

Figure 4

Comparison between actual crash counts for each 1-mile road segment in our data and expected crash counts for these same segments generated by our model. The extremely high correlation between the counts shows that our model fits our data.

Table 2

The mean and 95% credible interval for the effect of each roadway characteristic and the probability of the magnitude of the effect being greater than 0. The type of variable is either quantitative (Q) or categorical (C)

Figure 5

This plot shows additional risk (not attributable to traffic levels or the roadway factors in our data) identified by our model for I-5 (left) and I-90 (right). The solid curves represent a kernel smoother of the residuals to highlight patterns

The view from mile 166 on I-5, obtained from Google Street View. The Space Needle, as well as other Seattle landmarks, is clearly visible from the road, which may distract drivers, possibly contributing to the increased risk in this area

The posterior mean of the relative risk μ over I-5 (left) and I-90 (right). If μ < 1 , there were fewer crashes than expected due to traffic levels and if μ > 1 there were more crashes than expected

Acknowledgments

References

The posterior mean of the relative risk $μ$ over I-5 (left) and I-90 (right). If $μ < 1$ , there were fewer crashes than expected due to traffic levels and if $μ > 1$ there were more crashes than expected