Integrating Spatial and Temporal Approaches for Explaining Bicycle Crashes in High-Risk Areas in Antwerp ( Belgium )

The majority of bicycle crash studies aim at determining risk factors and estimating crash risks by employing statistics. Accordingly, the goal of this paper is to evaluate bicycle–motor vehicle crashes by using spatial and temporal approaches to statistical data. The spatial approach (a weighted kernel density estimation approach) preliminarily estimates crash risks at the macro level, thereby avoiding the expensive work of collecting traffic counts; meanwhile, the temporal approach (negative binomial regression approach) focuses on crash data that occurred on urban arterials and includes traffic exposure at the micro level. The crash risk and risk factors of arterial roads associated with bicycle facilities and road environments were assessed using a database built from field surveys and five government agencies. This study analysed 4120 geocoded bicycle crashes in the city of Antwerp (CA, Belgium). The data sets covered five years (2014 to 2018), including all bicycle–motorized vehicle (BMV) crashes from police reports. Urban arterials were highlighted as high-risk areas through the spatial approach. This was as expected given that, due to heavy traffic and limited road space, bicycle facilities on arterial roads face many design problems. Through spatial and temporal approaches, the environmental characteristics of bicycle crashes on arterial roads were analysed at the micro level. Finally, this paper provides an insight that can be used by both the geography and transport fields to improve cycling safety on urban arterial roads.


Introduction
As a key component of sustainable transportation systems, cycling has been actively promoted in cities throughout the world [1,2]; however, bicycle-related crashes have been associated with increasing numbers of fatalities and injuries [3][4][5][6] and the risk of crashes prevents people from using bicycles [7]. Compared with driving, cyclists have a higher probability of injuries in traffic accidents [8]. Unfortunately, bicycle crash risks are unclear [7] because current risk estimates mainly depend on general exposure (such as population or census data) [9,10] or insufficient exposure (such as traffic exposure or bicycle exposure) [11]. Despite the fact that bicycle crashes tend to cluster within a spatial

Research Goals
In order to determine future measures to improve cycling safety, this paper aims to: 1.
Reveal the spatial and temporal risk patterns of bicycle crashes (where? and when?) from a regional level to a locational level (a macro scale to a micro scale). To this end, a two-stage workflow (spatial and temporal approaches) is created for exploring bicycle crashes. Through the spatial approach, urban arterials are determined to have the highest bicycle-motorized vehicle BMV crash densities (see Figure 6a and Figure 7).

2.
Explain how arterial infrastructure affects bicycle crashes in the city of Antwerp (CA, Belgium) by examining possible risk factors.

Bicycle-Motorized Vehicle (BMV) Crash Studies
Bicycle-motorized vehicle crashes raised concerns in the 1970s when the annual fatality rate among children surpassed 500 in the Netherlands [14]. Since the 1990s, much research has focused on road crashes related to bicycles. Road crashes arise from the interaction between human factors (such as driving behaviours), environmental factors (such as traffic exposure, road facilities), and vehicle-related factors (e.g., driving speed) [15][16][17]. These factors account for 57%, 34%, and 13% of all bicycle crashes, respectively (Figure 1a) [18]. However, recent reviews on BMV crashes have identified that these factors account for 59.3%, 57.6%, and 15.3%, respectively [19]. As can be seen from Figure 1b, the influence of environmental factors in crashes involving bicycles and motorised vehicles may be comparable to that of human factors.
Sustainability 2019, 11, x FOR PEER REVIEW 2 of 28 area [12,13], and are likely to be over-dispersed during a time period, the majority of risk studies provide no explanation for these spatio-temporal aspects of bicycle crashes.

Research Gap
Bicycle crashes are mainly studied from the perspectives of transportation and geography. From the transportation perspective, most studies neglect the regional impact of bicycle safety on the macro-scale level. This holds especially true for crash analyses of specific locations, e.g., intersections or other crash-prone locations on the micro-scale level. Meanwhile, from a geographical perspective, many studies ignore the locational influence of safety attributes on the micro-scale level. This especially applies to crash analyses of geographic areas, such as spatial autocorrelation or many different kinds of cluster analyses on the macro-scale (regional) level. Moreover, from an epidemiological point of view, bicycle crashes are complex spatio-temporal phenomena with many contributing risk factors (e.g., weather conditions, motorized and non-motorized traffic volumes, road facilities, road traffic controls and driving behaviours) varying over space and time. Therefore, significant gaps are present in current bicycle crash studies.

Research Goals
In order to determine future measures to improve cycling safety, this paper aims to: 1. Reveal the spatial and temporal risk patterns of bicycle crashes (where? and when?) from a regional level to a locational level (a macro scale to a micro scale). To this end, a two-stage workflow (spatial and temporal approaches) is created for exploring bicycle crashes. Through the spatial approach, urban arterials are determined to have the highest bicycle-motorized vehicle BMV crash densities (see Figure 6a and Figure 7).
2. Explain how arterial infrastructure affects bicycle crashes in the city of Antwerp (CA, Belgium) by examining possible risk factors.

Bicycle-Motorized Vehicle (BMV) Crash Studies
Bicycle-motorized vehicle crashes raised concerns in the 1970s when the annual fatality rate among children surpassed 500 in the Netherlands [14]. Since the 1990s, much research has focused on road crashes related to bicycles. Road crashes arise from the interaction between human factors (such as driving behaviours), environmental factors (such as traffic exposure, road facilities), and vehicle-related factors (e.g., driving speed) [15][16][17]. These factors account for 57%, 34%, and 13% of all bicycle crashes, respectively (Figure 1a) [18]. However, recent reviews on BMV crashes have identified that these factors account for 59.3%, 57.6%, and 15.3%, respectively [19]. As can be seen from Figure 1b, the influence of environmental factors in crashes involving bicycles and motorised vehicles may be comparable to that of human factors.  [20,21]. A number of previous studies have explored the relationship between bicycle crashes and the physical environment [8,22,23]. Additionally, despite limitations due to issues of privacy, a few studies have explored the influences of human and vehicle factors on bicycle collisions [24][25][26].

BMV Crash Patterns
BMV crash patterns are fundamentally space-and time-related because the occurrence of a crash directly involves its specific location and time. Since bicycle crashes are mainly assessed through the perspectives of geography (which is related to spatial methodologies) and transportation (which is related to temporal methodologies), it is essential to review the current literature relating to these approaches.

Spatial Patterns of BMV Crashes
Many risk evaluations have focused on spatial crash analysis. Bíl [27] and Chen [28] studied the crash risk of road sections, and Hels and Orozova-Bekkevold [29] and Harris et al. [30] investigated the risk-influencing factors of road roundabouts and intersections, respectively. Compared with road sections, the majority of studies found an increased possibility of BMV crashes at intersections [4,28,[31][32][33]; however, the likelihood of severe injuries from BMV crashes on the road sections is more likely to increase [4]. In addition, Vandenbulcke et al. [34] indicated that distinguishing cities (built-up areas) from rural areas is important because bicycle crashes are closely associated with urban structures [21,35,36].
With regard to the spatial scale, crash risks can be assessed from a regional level (macro scale) [37][38][39][40][41][42] to a locational level (micro scale) [43,44]. Generally speaking, in the regional (macro) level, BMV risks can be analysed over certain geographic zones in order to understand the effect of environmental factors on accident occurrences, thereby improving traffic safety in the whole region [45]. On the other hand, in the locational (micro) scale, BMV risks can be observed on certain specific road entities (e.g., ramps, curved road sections) in order to determine explanatory factors contributing to accident events and thereby facilitating constructive countermeasures to reduce crashes [20]. However, few BMV crash studies so far have estimated the risk of certain specific locations (e.g., urban arterials, tunnels and bridges) [9]. The spatial approaches proposed in this paper provide generalisability, meaning the same countermeasures may be applied to areas with similar infrastructure characteristics [46].

Temporal Patterns of BMV Crashes
Although less bicycle risk estimations emphasise temporal approaches based on probability modelling due to inadequate time-related counting data (exposure), these approaches focus on the study of contributing factors that affect the frequency of BMV crashes occurring within a certain period [47,48]. A potential alternative approach is to collect exposure data at each location [49]. However, when there are high rates of bicycle collisions, such collected data are increasingly expensive and labour-intensive to obtain [12]. On the other hand, the definition of temporal approaches may refer to only time-related contributing factors (rather than time-related methodologies), where these estimations focus on bicycle risks influenced by traffic exposure, travel time, specific seasons, days, peak hours and road surface conditions related to the weather [20]. Dozza [50] examined cycling risk at hourly, daily, weekly, monthly and seasonal scales and discovered that BMV and single-bicycle crashes had significant risk differences in the dark and during weekend afternoons, peak hours and July. Kim [51] and Yan [52] both indicated that serious bicycle crashes occurred during peak hours, bad weather influencing driver visibility and in the absence of street lights at night. Chen [53] and Prati [21] showed that transient poor road surface conditions (e.g., following rain) may significantly increase the riding risk. Chen and Fuller (2014) [54] found that the probability of a bike crash increased by 2.13 times at night. In addition, early studies indicated that for a given location, rainfall was related to a 70% higher crash risk [55], and it might increase the severity of BMV crashes as well [56].
Note that the temporal analysis of bicycle crashes not only refers to time-related contributing risk factors, but also refers to a BMV crash methodology with exposure counting data (see Section 2.3).

Methods without Exposure Counting Data
From a methodological perspective, BMV crash estimates may be separated into two groups: exposure-based and non-exposure-based methods. So far, due to the lack of qualitative traffic counts within a location scale, few bicycle risk estimates have been investigated [9,24,53]. Non-exposure-based methods are viewed as preliminary steps to examine bicycle risks because they explore factors that may cause bicycle crashes with the use of descriptive statistics (such as medical records [57], police and questionnaire data [58]) or general statistics (such as odds ratios [59], logistic models [53,60], census data combined with spatial analyses [61]).
Of these spatial analyses, methods without exposure commonly analyse BMV crashes on a macro scale, such as traffic analysis zones [37,39,41,46], and are implemented in geographic information systems (GIS) spatial modelling, such as kernel density estimations [7,20,[62][63][64]. Since methods with exposure counting data are quite data-intensive, methods without exposure reduce the amount of required data and are important for crash locations with inadequate exposure.

Methods with Exposure Counting Data
Unlike the methods discussed above, exposure-related approaches using related statistical probabilities (such as Bayesian distributions [21], negative binomial distributions [41,43], Poisson distributions [46]), pay more attention to explaining the relationships between influencing factors and bicycle crashes. In addition, because "exposure to risk" is strongly related to average daily bicycle traffic (ADB) and average daily traffic (ADT) [11,65,66], the modelling is labour-intensive and time-consuming. However, exposure-based methods may have a better ability to explain their contributing factors [10,20]. Moreover, by combining these methods with spatial network reference units, such as hotspot analysis [67], cluster analysis [68], or density estimations [67], exposure-based methods may bring more accurate results.
Methods with exposure are strongly related to time because a crash probability model can be developed by controlling bicycle and motorized-vehicle exposure within a specified period to gain better understanding of risk factors [20,47]. However, exposure data are usually unknown [63], the role of exposure can be replaced or explained using (bicycle) lane kilometres [31,69], total number of trips [46], travel behaviour of cycling [49], (bicycle) commuting flows [63] or peak hour flows [70].

Risk Factors Associated with Cycling Environments
Considering the risk impact of cycling environments, the majority of studies have indicated that bicycle crashes are affected by road environments, road traffic controls, and road engineering facilities [8,20]. Particularly, some cycling environments may have an increased risk for bicycle crash frequency but a decreased risk for bicycle crash severity (such as highly-urbanized areas [44,71], peak hours [7], straight roadways [53], signalised intersections). The observed phenomenon may be entirely different, which means frequency risks may be lower but the severity of risks may be higher (e.g., curved road sections [51]). However, in most of the cases, crash frequencies may be consistent with their severity [20]. Overall, previous studies showed that heavy traffic (during peak hours), an increased number of lanes and a higher speed limit may lead to higher crash probabilities. However, up to now, a large number of studies have examined only risk factors related to the crash frequency [71] or severity [21,53]. For example, lower speed limits significantly lower bicycle and pedestrian crashes on urban road sections and junctions [4,59]; the turning movement of bicycles and motorised vehicles and urban road networks are aggravating and mitigating factors for bicycle crashes, respectively [4]; or increased bicycle exposure may reduce the severity of bicycle crashes [36,72,73].
Finally, based on the overall background of bicycle crashes, Section 5 collected 33 potential influencing factors with their detailed descriptions. The study assumed that other factors related to infrastructure (such as primary roads or secondary roads) may affect the risk of having a bicycle crash. Indeed, in the first stage of this study (the spatial approach), urban arterials are identified as high-risk areas, which may be attributed to complex traffic situations in road environments, including amongst others: higher traffic volume, more lanes and wider road widths [7,19,71]. However, prior bicycle crash studies have given little specific information about such risks, how these risks affect BMV crashes, and how they are affected by their BMV crash severities [7,20]. Therefore, the purpose of this research is to explore the patterns of bicycle crashes spatially and temporally, and to reveal how road environments influence bicycle crash risks. A further purpose is to explicate particularities in a case study (urban arterials of city of Antwerp) and reflect on existing BMV risk studies.

Methodology: Integrating Spatial and Temporal Approaches
For the spatial and temporal analysis of bicycle crashes, a two-stage workflow ( Figure 2) integrated geocoded data (such as X and Y coordinates, road sections, road junctions, among others), visual data (such as crash figures, aerial photos), traffic data and infrastructure data into the spatio-temporal analysis (for details, see Section 5). In stage 1, spatial references of analysis combined nine reference units with census data (the population at risk) for analysis. Within each reference unit, BMV crashes with corresponding severity levels were precisely located. Stage 1 included census data because an increased population may contribute to more BMV crashes. Using kernel density estimation as a spatial approach, urban arterials were found to have the highest density of bicycle crashes on a macro scale. The displayed crash patterns were further investigated based on road environments and statistical data by including more information on road networks (e.g., road sections, road junctions, annual updated transportation facilities, etc.). In stage 1, all bicycle crash data were implemented in ArcGIS 10.6.1 (A. 10.6.1; Esri: Redlands, CA, USA, 2019), and in stage 2, only urban arterial crashes were implemented in both ArcGIS 10.6.1 and Limdep 11 (L. 11; Econometric Software, Inc.: Plainview, NY, USA, 2016). Therefore, certain types of bicycle crashes, such as crashes on arterials, can be easier identified and further investigated.
Stage 2 defined the traffic volume data (along with bicycle count data) as temporal references for BMV crashes, and delineated potential influencing factors of arterial crashes. These factors varied over time (such as the timing when a bicycle facility was built) and remained independent through a correlation coefficient test. Negative binomial modelling was used as a temporal approach to estimate the risk of BMV crashes on a microscopic scale. Finally, risk factors associated with bicycle infrastructure were assessed by adopting a maximum likelihood estimation (MLE) and all arterial crash data were put in the Limdep environment.

Stage 1: Spatial Approach
To detect the spatial pattern of crash events, a weighted kernel density estimation (wKDE) [74,75] served as an initial step in the risk calculation to estimate crash density (Equation (1)). The wKDE placed a surface on each crash point, calculating the distance from a reference position to that point (di); wKDE was generally derived using: where the kernel function (K) was the intensity of the crash event i, influenced by population and distance; h was the search bandwidth (radius); n was the frequency of observed crashes; D(x,y) was the density at the position with x and y coordinates; W i was the weighted value of the crash event i, which considered the different severity of bicycle crashes [76]. Stage 1 replaced the fixed search bandwidth (h) with an adaptive search bandwidth P(u,v) in order to exclude uneven spatial distribution caused by the population at risk. Local population density was presented by the P function and was placed at the centre position (u,v) for reference. KDE was first applied to traffic crashes in 2008 [13,77]. The wKDE is the extension of KDE and has the advantage regarding determining the density of crash risks more accurately than KDE [74,75]. Second, such a spatial density approach makes a visual comparison analysis (Figure 7), providing risk class homogeneity for the entire studied area and enabling identification of high-density areas, such as arterial roads in this study.

Stage 1: Spatial Approach
To detect the spatial pattern of crash events, a weighted kernel density estimation (wKDE) [74,75] served as an initial step in the risk calculation to estimate crash density (Equation 1). The wKDE placed a surface on each crash point, calculating the distance from a reference position to that point (di); wKDE was generally derived using: Of note, this study proposes to use planar KDE rather than network KDE [13,78] for density estimations because: (1) By using an adaptive search radius, KDE may easily exclude the uneven spatial distribution caused by the population at risk. an area unit rather than that on a linear unit. However, a crash may be more properly recognised as a spatial point carrying a spread of risk (e.g., crashes caused by a discrepancy in elevation in the mountain areas [20]). (3) With a re-weighting function, wKDE can be adjusted to control for the severity of different crash points (e.g., slightly injured, severely injured, fatal crashes). However, network KDE does not distinguish crashes by injury severity, which means all bicycle crashes are treated the same. (4) Network KDE could possibly underestimate crash density. For example, KDE finds eight bicycle crashes within the search bandwidth, but network KDE only finds three bicycle crashes at the same studied area (Figure 3b) because network KDE does not consider a prolonging distance caused by road curvature or driving manoeuvres.
determining the density of crash risks more accurately than KDE [74,75]. Second, such a spatial density approach makes a visual comparison analysis (Figure 7), providing risk class homogeneity for the entire studied area and enabling identification of high-density areas, such as arterial roads in this study.
Of note, this study proposes to use planar KDE rather than network KDE [13,78] for density estimations because: (1) By using an adaptive search radius, KDE may easily exclude the uneven spatial distribution caused by the population at risk. (2) KDE ( Figure 3a) calculates the density on an area unit rather than that on a linear unit. However, a crash may be more properly recognised as a spatial point carrying a spread of risk (e.g., crashes caused by a discrepancy in elevation in the mountain areas [20]). (3) With a re-weighting function, wKDE can be adjusted to control for the severity of different crash points (e.g., slightly injured, severely injured, fatal crashes). However, network KDE does not distinguish crashes by injury severity, which means all bicycle crashes are treated the same. (4) Network KDE could possibly underestimate crash density. For example, KDE finds eight bicycle crashes within the search bandwidth, but network KDE only finds three bicycle crashes at the same studied area (Figure 3b) because network KDE does not consider a prolonging distance caused by road curvature or driving manoeuvres.
(a) (b) Figure 3. The different estimates between (planar) KDE (a) and network KDE (b) for the same crash points. To estimate the density value, KDE uses the whole two-dimensional space based on Euclidean distances and finds eight crash points within a search radius h, whereas network KDE only finds three crash points within the same radius in the network space based on linear distances.

Stage 2: Temporal Approach
The temporal approach estimated the crash risks with a suitable risk model. Temporal construction sites or bad road environments may accumulate more BMV crashes. With the help of aggregation by time, the characteristics of BMV crashes can be statistically described and visually detected (stage 2 of Figure 2). A Poisson distribution [47,79] was used for crash probability modelling by controlling motorized vehicle and bicycle exposure within an observed time. This modelling assumed that BMV crashes were random and obeyed a binomial distribution [71] in the observed period. According to this particular distribution, the probability of crashes and the influencing factors related to the probability of crashes can be estimated.
However, comparing the traffic flows (V i ) and bicycle risk (P i ), P i had a very small value because the motorized vehicle flows were usually much greater than the frequency of bicycle crashes. Therefore, a Poisson distribution was suitable to explain the binomial distribution of bicycle crashes [80]. The Poisson distribution with random, discrete, and non-negative characteristics was commonly applied to crash estimates; however, this distribution required that its mathematical expectation and variance were equal. In most situations, it was difficult to reach this constraint because the crash data were over-dispersed. To relax this constraint, inserting an independent error term ε i into the Poisson distribution allowed for the variance and mathematical expectation of the Poisson distribution to be unequal [61,81,82].
Moreover, it was assumed that, e ε i obeyed a gamma distribution, giving a variance of δ, and θ = 1/δ, thus a negative binomial (NB) distribution may be more suitable to explicate the prior Poisson distribution of bicycle crashes because it allowed the NB and Poisson distributions' variances to be different and their expected values to remain the same. The negative binomial modelling was given by (the detailed formula derivation can be found in recent research conducted by Wang et al. (2019) [20]): where i was the road junction or section index (see Figure 4a,b), n i was the number of BMV crashes under specific V i traffic flows, V i was the traffic flow of location i, P i was the BMV risk under traffic flow (V i ) and P(n i |ε i ) was the likelihood of n i bicycle crashes occurring under the inserted error term ε i . If θ was statistically significant, the NB distribution would be used in the BMV crash model. Otherwise, the Poisson distribution would be adopted. Equation (3) shows that the BMV crash risk (P i ) was determined by the bicycle flow F i of location i and a series of influencing risk factors (X i ), and β was a vector of coefficients of X i . To avoid the possible over-estimation of bicycle risk in stage 1 (Figure 6a), and to reflect the high cycling population at certain locations, F i was collected by this study according to field surveys on arterial roads. By adding the value of P i from Equation (3) into Equation (2), the final formula is shown in Equation (4) for the probability P(n i ) of having n i crashes.
This temporal approach hadd three advantages: (1) The crash risk approached 0 when there were few bicycle flows at location i. (2) β was made up of a series of coefficients. If the value of β was positive, it had an aggravating impact on the crash risk; otherwise, it had a mitigating impact on risk when the value of β was negative. (3) Unlike Wang et al. [20], historical crash severity index (CSI) was viewed as a contributing risk factor [20,41] integrated into the temporal modelling (rather than two separated models). Conducting the modelling this way estimated the crash risk on arterial roads because stage 1 had found that the severity of arterial crashes was in line with its crash frequency, which meant that on urban arterials, a location was found with an increased number of crash frequency, accompanied with increased crash injuries and fatalities. Determining and dealing with the road features of these collision-prone positions may greatly lower the risk of crash frequencies, severe injuries and fatalities.
To find the contribution of the crash risk, each influencing risk factor should be independent. Using Pearson's correlation examination, a pair of risk factors with a correlation ≥0.700 would not be together into the model [61,83], only the one with favourable explanatory abilities was included in the temporal modelling [84]. A total of 33 vectors of explanatory variables X i were later selected for the temporal modelling. These explanatory risk factors (X i ) were chosen based on crash types and their environmental features, which may have an aggravating or mitigating influence on the crash risk of urban arterials. Finally, using maximum likelihood estimation (MLE) [47], Equation (5) evaluated the unknown vector of coefficient β:  To find the contribution of the crash risk, each influencing risk factor should be independent. Using Pearson's correlation examination, a pair of risk factors with a correlation ≥0.700 would not be together into the model [61,83], only the one with favourable explanatory abilities was included in the temporal modelling [84]. A total of 33 vectors of explanatory variables Xi were later selected for the temporal modelling. These explanatory risk factors (Xi) were chosen based on crash types and their environmental features, which may have an aggravating or mitigating influence on the crash risk of urban arterials. Finally, using maximum likelihood estimation (MLE) [47], Equation (5) evaluated the unknown vector of coefficient β:

Case Study
The presented case study considered both statistical population and exposure counting data. BMV crashes were spatially and temporally combined. The proposed spatio-temporal approaches were applied to a case study of the city of Antwerp (Belgium) where a five-year BMV crash database was analysed, from 2014 to 2018. Before Section 6 presents the results, Sections 4 and 5 describe the data, studied area, and the detailed information regarding contributing risk factors.

Study Area
The city of Antwerp (CA) had approximately 521,000 residents and the population density was approximately 2500 residents per square kilometre. In the CA, the total length of the road network related to bicycle flow was 6269 kilometres. Around 44% of these roads (2750 kilometres) had bicycle infrastructure in the range of shared uni-directional bicycle lanes to separated bi-directional bicycle lanes, where 498 kilometres were exclusive cycle lanes ( Figure 5). Sustainability 2019, 11, x FOR PEER REVIEW 11 of 28

Analytical Environments
By exchanging file formats and combining three analysis software programs, the spatiotemporal approaches, introduced in Section 3, were built. ArcGIS 10.6.1 performed all spatial analysis, while Limdep 11 and SPSS 24 (S. IBM: USA, 2016) accomplished all temporal and statistical analysis. Manual input of police reports and address mapping successfully geocoded 99% of the bicycle crashes. By using the adaptive bandwidth [87], ArcGIS 10.6.1 performed kernel density estimations, obtaining high-density areas of BMV crashes.

Influencing Risk Factors
This study's purpose was to understand the spatial and temporal dynamics of road environments and engineering facilities. Other influencing risk factors, like the gender or age, could also provide important information for BMV crash prevention or intervention. However, it was not In addition, in the city of Antwerp, the average family owns 2.2 bicycles and 36.5% of the traffic modal share is from cycling traffic [85]. The current amount may be higher because within each studied year, the total amount of BMV crashes has steadily grown [86].

Data Attributes and Sources
Based on the GIS road networks and reported bicycle crashes from the police, the study developed a database. Figure 2 of Section 3 demonstrates the procedure of data collection. First, stage 0 conducted a complete literature review on the influencing factors. Bicycle crashes were then geocoded on the road network of the GIS (stage 1) and their road environment characteristics were determined from police reports, field surveys and aerial photos (stage 2). Finally, the database with influencing risk factors, information attributes and bicycle crashes was applied for the estimation of bicycle-motorized vehicle crash risks in the city of Antwerp. All information in the database was collected from (1) (6) field surveys conducted specifically for the study (see Table 1, column 3).
The following attributes were contained as a crash event's information: date, time, location, the weather, vehicle type, crash severity, crash type, lighting, road surface condition of the road and risk influencing factors (Table 1, column 1). Table 1, column 2, fully describes these risk factors and categorises them into three classifications: road environments, road traffic controls and road traffic facilities. Although the database did not contain information about liability, the scene of a bicycle crash was easier to reconstruct by using crash scene photos, standard crash report forms, official police reports and crash site figures depicted by trained police personnel. Bicycle crashes involving at least one motorised vehicle were considered in this study, resulting in 4162 reported crashes in this spatio-temporal database from January 2014 to December 2018. However, due to the invalidation of some crash locations, 42 reports were excluded. Finally, 4120 crashes (stage 1) and 1415 arterial crashes (stage 2) were the objects of analysis in this study.

Influencing Risk Factors
This study's purpose was to understand the spatial and temporal dynamics of road environments and engineering facilities. Other influencing risk factors, like the gender or age, could also provide important information for BMV crash prevention or intervention. However, it was not permissible to obtain data on human factors that might lead to identification of individuals in this research due to the privacy concerns [24]. Table 1 not only shows all of the risk factors, data sources, and detailed descriptions utilised in this study, but also mentions the relevant studies of these risk factors: (1) City of Antwerp (CA) had tram tracks and icy roads; (2) bicycle and traffic exposure data were acquired from Flemish Traffic Center (FTC) or from field surveys with the traffic cameras; (3) the study considered two traffic-calming measures: residential zones (30 and 20 km/h), and dynamic traffic-calming zones (car-free zones) during certain hours; and (4) the study also contained traffic signals, and turning movement of bicycles and motorised vehicles, since these were likely to be influencing factors leading to BMV crashes.

Results
This section concisely describes the spatial and temporal pattern of bicycle crashes. For the spatial scale, bicycle crashes tended to aggregate on arterial roads (see Figure 6c-f). A total of 38.42% and 35.68% of arterial crashes were caused by rear-end conflicts and side crashes with overtaking manoeuvres, respectively. Around 10% of BMV crashes were situated in the zone hotspots (see Figure 7). Within these zone hotspots, 78.52% of crashes were associated with urban arterials and 50.26% with road junctions. Similarly, throughout the city of Antwerp (CA), 21.70% of arterial crashes were in the red zones (zone hotspots), whereas only 3.01% of secondary road crashes were in these red zones. Arterial crashes were unevenly distributed and were disproportionately related to the length of the network. Figure 7 shows approximately 33.61% of bicycle crashes were concentrated in 28.79% of the total network length (i.e., the network length of arterials).

The Category and Number (%) of BMV Crash Occurrences
Category = 0-3 (see Table 1  For the temporal scale, arterial crashes were more likely to be clustered in summer (28.03%), on weekdays (82.99%) and during peak hours (46.41%, from 7 to 9 am, and from 4 to 6 pm). Descriptive statistics demonstrate that crash frequencies matched the Poisson distribution by following independent, scarce, and random events that occurred within the studied year, and were more consistent with the NB distribution because of their over-dispersed, non-negative features ( Table 2). The frequency's expected value was 5.16 but its variance was 303.78, unequal to the mean value, which preliminarily proved the NB modelling was more suitable for the explanation of influencing factors in the study. To understand crash patterns on arterial roads, Tables 2 and 3 categorised influencing factors into continuous and ordinal/nominal variables. Since some paired influencing factors were highly correlated (the correlation coefficient ≥ 0.7), each paired dependent variable was first examined for their explanatory ability regarding BMV crash modelling, and the ones with a less favourable explanation (e.g., light condition "night-time" was replaced with "night-time with lighting" and "night-time without lighting") were then excluded from the model.

The Category and Number (%) of BMV Crash Occurrences
Category = 0-3 (See Table 1 After running the estimation procedure in stage 2, the result of the likelihood-ratio test for the estimation of dispersion parameter (α) shows the significance was 0.5363 with a p-value = 0.0001. This meant that crash data were over-dispersed and NB modelling at this statistical significance was superior to and thus replaced Poisson modelling. This result reconfirmed that the NB modelling was more suitable for evaluating bicycle crashes in the city of Antwerp. By removing certain variables (e.g., light conditions) to solve correlations, 32 independent influencing factors remained in the final risk model. To identify significant influencing risk factors, NB modelling was first performed in the temporal workflow. Second, to enhance the modelling performance, the same NB modelling was performed using only influencing factors with significant levels. Figure 8 shows the estimated coefficient and their p-values for these risk factors. The coefficient of these factors indicated the comparative risk level of the overall NB modelling. A positive coefficient shows that Xi was a possible aggravating factor related to bicycle crashes, while a negative value suggested a mitigating factor, corresponding to Equation (3) (see Methodology). For example, a higher risk of BMV crashes was found in locations during the morning or afternoon peak hours. Oppositely, a lower risk was found in locations with main cycling routes or with signalised facilities. The results show that more than two thirds of the influencing risk factors were significant in the temporal model (Figure 8a-c), thus giving a suitable classification under the model.

Spatial Dimension
In stage 1, the 4120 reported crashes appeared at 1125 different locations of the CA. Within the studied areas, the crashes were strongly aggregated, as shown in Figure 6. During the entire survey period, the areas with the highest bicycle crash density occurred in districts 2, 4, 5, 7 and 8. These areas are located in the city centre, enclosed by the ring road R1 from the east, south and west, and have a relatively large number of bicycle facilities (see Figure 5) and urban road networks (see the top right corner of Figure 6). Arterial roads were identified as having the highest density of bicycle crashes in the city of Antwerp. A total of 85% of all arterial crashes were situated in these five (2, 4, 5, 7 and 8) districts. Descriptive statistics at this stage indicate that other crash features (e.g., road conditions, road engineering facilities, road traffic controls) in these districts did not significantly differ from the rest of the city of Antwerp.
The spatial aggregation of bicycle crashes was consistent with the study by Lovelace et al. [63] and Loidl et al. [24] who found significant spatial aggregations of bicycle crashes on the macro (regional) scale. The former also used KDE for the spatial approach. However, through a temporal approach, it becomes more essential for this case study to pay further attention to arterial crashes on the micro (locational) scale.

Temporal Dimension
Stage 2 focused on the arterial crash analysis. The 1415 reported arterial crashes were located at 274 different locations. The temporal pattern of BMV crashes was revealed by influencing risk factors (e.g., the relation of traffic flows and road environments [108], dynamic road traffic controls [46,53,93]). The results are discussed based on the existing bicycle crash literature.

Road Environments
The results reveal that traffic exposure may provide a significant contribution to the risk of bicycle crashes on urban arterials. During peak hours in the morning and afternoon, increased traffic may result in an increase in bicycle crashes. However, the degree of their influence is not huge, expressed by the coefficients of 0.1199 and 0.1910, respectively. This may be due to the fact that road users typically slow down under higher traffic flows during peak periods, thereby reducing the impact of urban traffic flows. This result is also confirmed by previous studies [3,71,109,110]. However, the results reveal that bicycle flows were not significant in aggravating the risk of bicycle crashes in the city of Antwerp. This may be attributable to the concept of "safety in numbers" [36,111] and is especially suitable for reported crashes (major and fatal injuries) [110,111]. Such a result may be caused by the effect of "risk compensation". A possible explanation for this phenomenon is that motorists may drive more carefully when they see groups of cyclists, implying that group cycling may change driving behaviour [7,34].
The study has also demonstrated that injury severity was positively related to the risk of BMV crashes, expressed by the coefficient of 0.8424. This means that an urban arterial section with higher crash frequencies may be considered more dangerous because it also had more serious bicycle injuries [20,21,112]. In addition, while bicycle crash risks at arterials were rarely affected by changes in the season, day or hour, darkness was a common influencing factor for reducing bicycle crash risks, which is consistent with previous findings: night-time may contribute to a decline in cycling, especially during wintertime [113], and sunshine is identified as an influence on bicycle ridership [114]. Although the risk of bicycle crashes was much higher in daylight than in darkness, our study has still indicated one out of five severe crashes were related to darkness, implying that improved lighting conditions or wearing visible clothing may reduce the value of the crash severity index (CSI) [69,115], and then may indirectly lower the risk of BMV crashes.

Road Engineering Facilities
Road engineering facilities have a significant impact on the risk of BMV crashes. This risk may be attributed to the geometric design of the road environment. Road sections (i.e., the absence of intersections) may significantly lower the risk of bicycle crashes, as expressed by the coefficient of −0.4423 (p-value < 0.01). This result concurs with recent studies that show the majority of BMV crashes happened at road junctions [7,12,31,116]. The absence of lighting systems at night influenced the risk of BMV crashes, as shown by the coefficient of 0.1976 (p-value < 0.05). This is consistent with the existing literature, indicating that inadequate road lighting is an important risk factor [69,115] and sufficient lighting may reduce bicycle crashes by 58% [117]; road lighting may improve the visibility of cyclists, thus reducing the risk of BMV crashes.
In addition, road categories and land use may also influence crash risks, as observed by Kaplan and Giacomo Prato [35], as well as Vandenbulcke et al. [7]. Regarding road categories, due to mixed traffic with heterogeneous vehicular velocities, urban arterials intersecting with secondary roads and residential areas may increase cycling risk, as shown by the coefficients 0.8298 and 0.4022, respectively. Secondary roads were twice as risky as residential areas because of more complex traffic situations. Implementing protected or divided bicycle facilities may greatly reduce the crash risk. Residential areas connected to arterial roads were related to increased bicycle crash risks. This may be because in narrow neighbourhood streets, the conflict between cyclists and motorists is more likely and the majority of neighbourhood streets lack separate bicycle paths. This conclusion is consistent with findings by Chen et al. (2018) [118]. Therefore, the implementation of one-way roads with bicycle paths may be an effective countermeasure to reduce risk. Second, crash risks may also be influenced by land use features [67]. In Antwerp, many roads within the central business areas were constructed earlier than the average urban road and have heavy traffic during rush hours. However, arterial roads located in the central business district (CBD) had lower crash risks. This may represent the effectiveness of municipal efforts to improve bicycle safety. For example, most arterial roads located in the CBD have been monitored with strict traffic law enforcement, thus enhancing bicycle safety [71]. Another possible explanation is that areas with dense road networks are mostly CBDs where bicycle facilities are well-designed with separated lanes, thus the conflicts between cyclists and motorists do not increase with road density [118].
The crash risk may also be influenced by lane types. In the CA, most lanes on arterial roads were correlated to the heterogeneity of traffic speed and were riskier than lanes with only fast or slow transport modes. The results demonstrate that tram routes and bus routes on the mixed lane may increase crash risks (coefficients = 0.4450 and 0.6624, respectively), while one-way bicycle paths and main cycling routes seemed to decrease crash risks, as shown by the coefficients −0.4976 and −0.1075, respectively. Previous studies also confirmed that tram track or bus route crashes on mixed lanes were significantly higher than arterial roads without bicycle infrastructure [7,23,43,119], and crash risks on the road with bicycle paths, compared to the one without, are reduced by about 25% [120,121]. Additionally, the results show that a longer distance from the bus stop [116] may have a lower possibility of having crashes, as shown by the coefficient of −0.1781. Therefore, physically separated bicycle paths near bus stops or tram tracks may greatly reduce this type of BMV crash [7]. However, an exclusive lane for two-way cyclists significantly raises BMV crash risks on arterial roads (coefficient = 0.4262), mainly because motorists do not see cyclists coming from right/left directions (two directions) [43,103,122]. Moreover, when two-way dedicated bicycle lanes became one-way (or the reverse), there was an increased risk because a sudden change required more reaction time for motorists and cyclists to respond. Marking bicycle crossings with coloured pavement at intersections, or providing physically divided and protected bicycle facilities, may greatly reduce the risk of BMV crashes, thereby improving road safety.
Illegal overtaking manoeuvres [123] have also been seen as a contributory factor for the risk of BMV crashes, as shown by the coefficient 0.2791. For example, on arterial roads, black spots were situated where there are overtaking-prohibited lanes (with marked lines). The results indicate that 72% of the conflicts at these locations were inappropriate lateral collisions, while few were frontal and rear collisions, emphasising the importance of maintaining adequate lateral clearance between bicycles and motorised vehicles (see References [98,124,125]). Physical lane boundaries or divisional facilities (e.g., channelization) may prevent vehicles from lane changes, effectively reducing BMV crash risks on arterial roads (Wang et al., 2019) [20]. Proper mitigating countermeasures for arterial roads (e.g., the expansion of the sidewalk for the purpose of cycling, the reallocation of bicycle paths on the sidewalk, the implementation of segregated bicycle paths or the narrowed width of traffic lanes to prevent motorised vehicles from inappropriate overtaking manoeuvres) may significantly lessen the risk of BMV crashes.
The dimension of road engineering facilities [98,123] may also have a significant effect on the risk of BMV crashes. The risk may be attributed to the increased size of arterial road junctions. In the CA, large-scale junctions often face far more complex traffic situations, such as high volumes of traffic, complex traffic compositions and speed differences between motorised vehicles and bicycles, thereby raising the risk of BMV crashes on arterial roads. Additionally, a road section with a longer length may raise the risk of BMV crashes, as shown by the coefficient 0.1279. In other words, increased length of road sections may result in increased activities of cycling, leading to a raised risk of cycling-related crashes and injuries [12,20,126].

Road Traffic Controls
Finally, the risk of BMV crashes may be associated with road traffic controls. This study has demonstrated that roads with higher speed limits may raise the risk of bicycle crashes. This study has further shown that roughly 75% of bicycle crashes took place on urban arterials with speed limits ranging from 50 to 70 km/h. Therefore, at some crash-prone locations, vehicles on arterial roads might be advised to limit their speed to below 50 km/h [4], and cyclists are advised to have a minimum sufficient sight distance of 38 meters [127,128]. In addition, when the speed limit of arterial roads exceeds 50 km/h, shoulder curbs or curb lanes should separate higher traffic flows from cyclists (of note, the Danish Cycling Embassy recommends a threshold of 60 km/h) [129].
The existence of signalling facilities may reduce the risk of BMV crashes. These facilities reduce the speed of vehicles, significantly decreasing the risk of conflicts between bicycle and motorised vehicles [130], as seen with a coefficient of −0.3017. The results are similar to those of Sweden [105] and Canada [30], suggesting that the presence of bicycle lanes, in conjunction with traffic signals, may notably reduce the risk of bicycle injuries. From the classification of these signals (Table 1, column 3), traffic signals prioritising pedestrians and bicyclists are advised as a measure to prevent BMV conflicts and improve road safety.
An increased cycle length of traffic signals (i.e., more than 120s to complete a cycle of green, amber, and red indications for both bicycle and motorised vehicle phases together) may result in a high risk of BMV crashes, as seen with the coefficient of 0.1669. The longer time of green indications may induce aggressive driving behaviour (e.g., violation of road markings, failure to make way, inappropriate lane-changing, and the increase of travelling speed [20]). On the other hand, the increased cycle length may entail a longer duration of red indications for cyclists. Bicycle queues may enlarge, occupying the whole lane (i.e., "spread effects" [107]), creating more conflicts between bicycles and motorised vehicles. Therefore, it is recommended that bicycle signal phases be added [131] or to appropriately shorten the cycle length of traffic signals to prevent erratic driving behaviour and potential conflicts, thereby lessening the risk of BMV crashes and injuries [106].
One-way roads [51] were correlated with a decreased risk of bicycle crashes (coefficient = −0.1177). One-way roads are expected to be safer than two-way roads because traffic situations are less complex, enabling motorists to more easily notice cyclists [51,132]. One-way arterial roads may also result in fewer crashes involving bicyclists because there are fewer turning movements. The results also demonstrated that motorist manoeuvres [133] induce a high risk of bicycle crashes (coefficient = 0.2388). More than half of the crashes (55.6%) occurred when cyclists were riding in a straight line and drivers were turning, similar to prior data [4]. One probable explanation is that motorists turning right may pay more attention to motorised vehicles or bicycles coming from the left, thus failing to see cyclists from the right until it is too late to react. However, the crash risk may still be decreased with proper interventions at arterial junctions with high turning manoeuvres (e.g., adopting bicycle (green) signal phases [116], implementing coloured bicycle crossings, eliminating junction bus stops or using junction refuges), which may greatly reduce crashes with left-turning and right-turning vehicles [134].

Meaning of the Two-Stage Workflow
A two-stage workflow focuses on bicycle and motorised vehicle crashes and combines advantages from recent studies conducted in geography and transportation fields. By using a two-stage strategy to assess bicycle crash risks, a spatio-temporal workflow opens new research directions for the analysis of traffic crashes (i.e., models aiming at estimating the risk of BMV crashes from a macro scale (a region) to a micro scale (a location/road network)). Compared to conventional methods for the analysis of bicycle safety (such as a crash frequency model), the adoption of this two-stage strategy has a number of methodological advantages: (1) In stage 1, the evaluation of bicycle crash risks had already been made possible despite inadequate bicycle traffic exposure; (2) both census/population data (stage 1) and exposure data (stage 2) were included in the models, improving the accuracy of estimates; (3) in the second stage, the detailed data collection of each crash point avoided potential errors arising from the arbitrary aggregation of point data (crash points) in the first stage, thus reducing the risk of point data aggregation; (4) the characteristics of bicycle crashes could be better understood through visual and statistical analysis from a macro-to a micro-level; (5) stage 1 avoided the expensive work of collecting counting data, thus stage 2 minimised labour-intensive and time-consuming analyses, providing conclusions about the influence of different environment conditions and road facilities; (6) in comparison with the traditional crash black-spot approach (Figure 8b), the prediction of potential crash risks could be provided for locations where bicycle crashes were unknown or underreported; (7) finally, stage 2 could resolve missing values as the scope had reduced from the regional scale to the local scale and could be achieved via field investigation and manual inputs from original reports (e.g., crash scene figures).

General Conclusions
In the case study, BMV crashes may be explained by a series of spatial and temporal phenomena. This study utilised the two-stage workflow, aiming to better understand "when" and "where" BMV crashes appeared from a regional (macro) scale to a locational (micro) scale. This study supports the use of the two-stage strategy because regional studies are not suitable for locational risk assessments [24]; however, locational studies may overlook the overall tendency of crashes on the regional scale. Although the study results are specifically associated with the presented case study, general conclusions may be drawn.
(1) The two-stage workflow may capture the patterns of BMV crashes in the city, thereby measures can be suggested to reduce bicycle crashes and crash risks. (2) This strategy may also be applied to other disciplines (or other cities) and makes analysing point events possible over space and time because stage 1 may be viewed as an initial step for the identification of hot-spot areas, where these risky areas may be further addressed carefully in stage 2. (3) Up to now, few studies have focused on the pattern of BMV crash risks from a macro-to a micro-scale, mainly due to a lack of available data. Through stage 1 to stage 2, the study makes the investigation of traffic and bicycle flows possible, thus understanding the influence of traffic controls, engineering facilities and road environments on BMV crashes. (4) Utilising spatio-temporal approaches to assess crash risk is more effective than utilising conventional black-spot approaches (Figure 6b). Spatio-temporal approaches enable the potential crash risk of each location in the entire region to be calculated such that the hazardous areas can be further addressed according to their different road characteristics (e.g., road pavement materials, type of lane, and traffic volume, etc.). (5) Finally, the spatio-temporal approach incorporates the construction year of road engineering facilities [7], dynamic road traffic controls and directional traffic volumes, leading to a more precise analysis of existing environment conditions on the arterial road.

Limitations and Further Research
It should be mentioned that this study has several limitations. First, the collection of data is labour-intensive and time-consuming. As can be seen from Table 1 and its footnote, although stage 2 has included 33 influencing risk factors, stage 1 has made only 15 risk variables available for the statistical and spatial analysis (4120 observed crashes). Second, this study has excluded some influencing risk factors associated with privacy concerns (like vehicle and human-related factors). Finally, the study aimed to understand the influence of bicycle crash risks by evaluating the influencing factors of road environmental conditions, road traffic controls, and road facilities. Further studies might include more influencing risk factors for road safety analysis.