category: literaturenote citekey: brizredoninvestigationconsequencesmodifiableareal2019 title: "Investigation of the consequences of the modifiable areal unit problem in macroscopic traffic safety analysis: A case study accounting for scale and zoning" authors: "Briz-Redón, Álvaro; Martínez-Ruiz, Francisco; Montes, Francisco" year: 2019 date: 2019-00-00 2019 doi: 10.1016/j.aap.2019.105276 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2019.105276" zotero_key: TQRRFFJI zotero_storage: 3G2H855K collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Briz-Redón, Álvaro"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

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

Investigation of the consequences of the modifiable areal unit problem in macroscopic traffic safety analysis: A case study accounting for scale and zoning

Álvaro Briz-Redóna,⁎ , Francisco Martínez-Ruizb , Francisco Montesa

^a Statistics and Operations Research, University of València, C/ Dr. Moliner, 50, Burjassot 46100, Spain
^b Statistics Office, City Council of València, Spain

ARTICLE INFO

Keywords: Modifiable areal unit problem Traffic safety analysis Traffic analysis zones Conditional autoregressive model Geographically weighted regression

ABSTRACT

Traffic safety analysis at the macroscopic level usually relies on previously defined areal traffic analysis zones (TAZs) that are used as the units of investigation. Hence, statistical inference is made on the basis of such units, implying that the consideration of a certain TAZ configuration may influence the results and conclusions achieved. Regarding this, the modifiable areal unit problem (MAUP) is a well-known issue in the field of spatial statistics, which refers to the effects that arise in statistical properties and estimations when there is a change in areal units of analysis.

In this paper, the consequences of MAUP have been investigated through a dataset of traffic crashes that occurred in Valencia within the years 2014 and 2015 and two common statistical models: a conditional autoregressive model and a geographically weighted regression. In the absence of an established TAZ scheme for the city, four classes of basic spatial units (BSUs) were considered: census tracts, hexagonal units and two types with construction based on the structure of main roads and intersections of the city. Each of these BSU types was specified at different levels of spatial aggregation. The main research objective was to investigate the final effects that changes in BSU type and scale have on model parameter estimations, but also the specific alterations that MAUP causes to data in terms of the distributional characteristics of the response, multicollinearity among the covariates and covariates' spatial autocorrelation.

The results showed the presence and severity of MAUP for the dataset and area that were analysed. Although effects from scale variations were more moderate, changing the BSU type affected the results severely. The joint use of hexagonal units and a conditional autoregressive model achieved the best performance among all the possibilities explored, but the choice of a proper BSU unit should rely on more factors. Despite MAUP effects, educational centres showed a consistent (and negative) association with traffic crashes, a fact possibly related to their distribution across the whole city. Other covariates revealed a positive correlation with crash counts, but these findings were more uncertain given the discrepancies found at different scales and zonings.

1. Introduction

1.1. The modifiable areal unit problem

Traffic safety analysis at the macroscopic level requires the definition of a basic spatial unit (BSU) for performing the analysis. Hence, the whole area of investigation needs to be covered by BSUs that allow researchers to analyse the incidence and causality of traffic crashes across it. The definition of BSUs can be done both manually, through the advice of experts of the field, or automatically on the basis of an algorithm specifically designed for BSU delineation.

The choice of a certain BSU over an area of interest is closely related to the well-known modifiable areal unit problem (MAUP). MAUP refers to the effects that carries the change from a collection of BSUs to another with regard to statistical inference and interpretation (Openshaw, 1984). In a seminal paper, Openshaw (1977) presented the two main factors that need to be addressed for the delineation of an area into BSUs: scale and zoning. Scale, or aggregation level, refers to the number of zones in which the whole area of study is subdivided for performing the analysis. Hence, given a scale, zoning is the way the BSUs are joined forming the zones of analysis while preserving the specified scale. Openshaw (1977) proved that the election of the zones has an effect on

E-mail address: alvaro.briz@uv.es (Á. Briz-Redón).

^⁎ Corresponding author.

spatial interaction models, in terms of fitting and parameter estimates. For this reason, he proposed a methodology in order to find the zoning subdivision of the area of analysis that optimizes model performance. More specifically, Openshaw also studied the consequences of MAUP on linear regression (Openshaw, 1978) and correlation coefficients (Openshaw, 1979), although he recognized that there are high difficulties for assessing the problem theoretically, leaving simulation studies as the main tool available for its approach.

Years later, Fotheringham and Wong (1991) extended the examination of MAUP to multivariate statistical analysis, considering multiple linear regression and multiple logistic regression within the context of two classical administrative divisions: block groups and census tracts. The aggregation of both divisions at different scales allowed observing that MAUP was capable of creating a severe instability in parameter estimates when the zoning or, more remarkably, the scale were modified. These authors found that the interpretation of some of the variables included in the models could be dramatically altered due to MAUP, as changes in the signs of the parameter estimates were appreciated. In addition, goodness of fit was observed to grow monotonically as aggregation level got increased.

Despite the fact that the study of the MAUP and its consequences in statistical inference have mainly been of descriptive or exploratory nature, recent works are trying to fill this gap by providing more accurate measurements of MAUP effects. Remarkably, Duque et al. (2018) have proposed a nonparametric test, S-maup, that measures the sensitivity to MAUP of a spatially intensive variable. Therefore, the S-maup test can be used to determine the level of aggregation at which MAUP effects do not impact the statistical analysis severely. As a drawback, the S-maup test lacks a theoretical definition. Indeed, an extensive simulation procedure was implemented by the authors in order to be able to supply critical values for different levels of scale and autocorrelation for the spatial variable. In a more observational work, Lee et al. (2018) have investigated the effects of MAUP in means, variances and Moran coefficients, considering several scales and levels of autocorrelation for the variables involved. They have concluded that MAUP effects are not strong on means, unless a very high spatial autocorrelation is present, and that higher levels of aggregation tend to decrease the variance.

1.2. TAZ delineation and MAUP effects in traffic safety analysis

The convenient delineation of an area into traffic analysis zones (TAZs), which behave as BSUs from the perspective of the present research, requires several considerations to be made, although the guidelines suggested in literature are usually varied and even contradictory. In a pioneering work, O'Neill (1991) proposed six criteria for the delineation of TAZs, including the homogeneity of socioeconomic characteristics, population and trip attraction levels, the minimization of intrazonal trips and the employment of physical, historical or administrative boundaries. Martínez et al. (2009) made use of a mobility survey available for the city of Lisbon (Portugal) in order to design TAZs fitting the following four criteria: boundaries are set over roads presenting a low trip generation density, intra-TAZ trips are minimized, TAZs with a very low or large number of trips are avoided and homogeneity within a TAZ is pursued as much as possible. Dong et al. (2015) used a K-means algorithm to classify a predefined set of cell areas according to primary features (traffic volume, hourly inflow, outflow and incremental flow) and optimizing features (peak and valley values for the primary features).

Efforts have also been made in order to account for the boundary effect in TAZ delineation. Siddiqui and Abdel-Aty (2012) proposed the distinction between boundary and interior pedestrian crashes considering a buffer of 100 ft from TAZs boundaries. Covariate information was weighted in the case of boundary crashes depending on the length of shared boundary between contiguous TAZs. Then, the specification of two analogous models for both types of crashes allowed the detection of differential effects for some of the covariates included in the study. Furthermore, there exists some simple methods that allow the allocation of traffic crashes occurred near TAZ boundaries, including halfand-half ratio, one-to-one ratio and ratio of exposure. Very recently, Zhai et al. (2018a) proposed a novel model-based iterative method for assigning crashes located close to boundaries. This was proven to produce better predictions at the BSU level than the other boundary assignation methods and to increase the number of significant covariates detected.

Several authors have investigated MAUP in the context of traffic safety analysis, which are now briefly discussed. First, Thomas (1996) noted that changes in scale may alter the probability distribution that best fits the nature of the available crash counts. Indeed, Thomas (1996) worked at the road segment level in order to infer three length thresholds that would require a distinct modelling strategy for crash counts: a Poisson distribution for very short segments (less than 1 hm), an intermediate empirical distribution for middle segments and a normal distribution for long segments (more than 20 hm).

In the last decade, however, most of the research studies related to MAUP were settled in the context of areal units of analysis. For instance, Lee et al. (2014) took the more than 1000 TAZs already defined for several counties in Central Florida (USA) and combined them into new subdivisions of the space containing from 100 to 1000 TAZs (in intervals of 100 TAZs). Specifically, total crash rates available for the period of study were employed by the regionalization algorithm for the obtention of sets of homogeneous zones containing the different number of TAZs specified. The Brown-Forsythe test was applied in order to check how the changes in TAZs affected the variance of crash rates. A moderate value of this test with a certain TAZ system represents an optimal situation, which means that the scale at which the TAZs are defined is suitable for detecting both local and global variation. Thus, the definition of 500–700 TAZs was found optimal for the data analysed in Lee et al. (2014). Xu et al. (2014) tested different TAZ schemes including from 50 to 738 units of analysis. They suggested the use of 350 or more TAZs (for their case study) in order to reduce MAUP effects because for this scale they found a superior number of significant covariates and more stable coefficient estimations. Similarly, Ukkusuri et al. (2012), used ZIP codes and census tracts as TAZs, determined that a finer aggregation level (census tracts in their study) was more suitable for data modelling as it enables a higher data variability and greater explanatory power. Abdel-Aty et al. (2013) modelled crash counts occurred at two American counties at the level of TAZs, block groups and census tracts, considering total, severe and pedestrian crashes. Relevant differences were found in terms of the number of significant variables that were yielded by models based on different spatial units and, more specifically, in the type of factor (roadway related vs. commute related) providing more significant variables. Amoh-Gyimah et al. (2017) investigated several spatial units that may be used as TAZs, including statistical area levels, postal areas, state electoral divisions, grid cells and Thiessen polygons developed around the centroids of Melbourne Integrated Transport Model. These authors made use of several statistical models to provide a more complete perspective of the effects that the choice of TAZs can lead to. The presence of MAUP was evident as they observed that a reduction in the number of zones produced an increase in the number of significant variables. Furthermore, they concluded that the selection of the modelling technique is another important factor that may reduce MAUP. Indeed, geographically weighted Poisson regression appeared to be less affected by MAUP than random parameter negative binomial. Finally, these authors suggest the use of Thiessen-based and grid cells for prediction purposes, according to their results. Zhai et al. (2018b) applied a multivariate Poisson lognormal model with multivariate conditional auto-regressive prior on block groups, census tracts, zip codes and predefined BSUs for the analysis of traffic crashes occurred in a county of Florida (USA). Important variations were found in relation to coefficient sign, magnitude and significance, and the larger units showed a superior forecasting performance. In addition, the detection of high-crash locations revealed some unexpected situations, as certain zones that were shared by all the BSU configurations showed a completely opposite behaviour depending on the underlying BSU-dependent model being used for such assessment.

Finally, in a review paper regarding the effects of MAUP in traffic safety analysis, Xu et al. (2018) proposed four potential solutions: avoiding data aggregation, considering the spatial variation of the covariates employed for data modelling (an issue that is usually skipped), defining an optimal zoning system for the analysis and conducting sensitivity analyses in order to check for MAUP presence and magnitude, regardless of the strategies undertaken for attempting its reduction.

In this paper, we carry out a complete investigation of MAUP effects from a dataset of traffic crashes occurred in Valencia (Spain). Whereas some related papers in the field have only focused on scale (Lee et al., 2014; Xu et al., 2014), most of them have tested different zonings without controlling the scale factor explicitly. This fact makes it challenging to determine if MAUP effects are a consequence of scale, zoning or the interaction between the two. Our paper tries to fill this gap with a simultaneous investigation of several BSUs and aggregation levels that allow the distinction between scale and zoning effects, in seeking to provide a more complete depiction of the phenomenon. Two modelling approaches, conditional autoregressive models and geographically weighted regressions have been used for this objective, following the choices of similar papers. Furthermore, we have specifically investigated how the changes in scale or zoning affect several questions involved in any macroscopic statistical modelling. These include the spatial autocorrelation of the covariates, multicollinearity among covariates and the basic distributional characteristics of the response variable. The investigation of MAUP usually focuses on the changes that finally arise in the estimation of model parameters after a switch of scale or zoning, but the changes in the underlying characteristics of the data being modelled are frequently overlooked. We also try to provide more insights on this issue.

2. Data

2.1. Crash dataset and road structure characteristics

A total of 18,037 traffic crashes that took place in the city of Valencia (Spain) during the years 2014 and 2015 were analysed (Fig. 1a). Geographical coordinates for each of these crashes and information regarding the date and hour of occurrence were provided by the Local Police of Valencia. The available coordinates were used to locate the crashes on a spatial representation of the road network of the city (linear network), as a guarantee of accuracy. This road network has a length of 840.3 km (with a diameter of almost 11.6 km) and contains 6110 road intersections. Arterial roads of Valencia, which were employed to define BSUs, extend up to 168.3 km and are also displayed in Fig. 1a.

2.2. Covariate definition

Several covariates were constructed to explain the incidence of traffic crashes among BSUs for the years of study, which were classified into environmental, network-related and socioeconomic. Environmental covariates included the consideration of different services (public or private) that are located along the road network which are known to influence the dynamics of traffic flow and in consequence are likely to affect crash rates. The services selected were schools (from preschool to high school level), bars/restaurants, hotels, private companies (mainly financial, legal or insurance), bus and tram stops.

Network-related covariates were precisely derived from the information provided by the road network structure, and included nonpedestrian road length, which was considered as an exposure, average road betweenness and number of road intersections (involving any road type, main or not). Betweenness is a measure of network connectivity and was computed for each road segment of the network according to the next formula (Freeman, 1977):

$$Be = \sum{i \sim i} \frac{\sigma{ij}(e)}{\sigma{ij}}$$

where i and j are vertex of the network that are connected by a path (i ∼ j), σij the number of shortest paths between i and j and σij(e) the number of shortest paths that connect i and j while passing through the edge e of the network.

Finally, socioeconomic and demographic information was introduced through the percentages of population in the range 16–24 and over 65 years, and also with the average power of cars (in hp), which clearly correlates with economic status.

It is of need to highlight that the data that was used to construct this set of covariates for different zonal schemes and levels of aggregation was available in point-referenced format. Hence, it was possible to aggregate the data at any desired level of aggregation or zoning system. In addition, the availability of the digitized version of the road network of the city allowed the computation of the betweenness or the number of road intersections. All these steps were carried out through specific GIS functions available in the R programming language (R Core Team, 2018) (Table 1).

2.3. BSU definitions

In the absence of an established TAZ configuration for the city of Valencia (which is probably the most used areal unit in the field of traffic safety analysis), several possibilities were explored for the

Fig. 1. Points representing the locations of traffic crashes that occurred in Valencia during the years 2014 and 2015 (a) and time series (displayed by hour and weekday) of traffic crashes observed in Valencia during the same period (b).

Table 1 Description of the covariates defined for the analysis and basic statistics of these covariates for the four BSU configurations tested (in their original configuration, prior to aggregation/regionalization).

| Type | Variable | | BSU configuration | | | | | | | | | |---------------------------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------|--------------------------------------------------------|--------------------------------------------------------|---------------------------------------------------------|---------------------------------------------------------|---------------------------------------------------------|--------------------------------------------------------|---------------------------------------------------------|--|--| | | | CTs | | TMs | | | TIs | | | | | | | | | SD | Mean | SD | Mean | SD | Mean | SD | | | | Crashes | No. of traffic crashes (CRASH) | 31.74 | 36.11 | 31.29 | 30.50 | 47.52 | 46.52 | 34.43 | 35.86 | | | | Environmental | No. of undergraduate educational centres (EDU)
No. of bars/restaurants (BAR)
No. of companies (COMP)
No. of hotel rooms (HOT)
No. of parking zones (PARK)
No. of bus stops (BUS)
No. of tram stops (TRAM) | 0.62
8.49
29.54
14.34
0.16
1.67
0.05 | 0.87
9.35
49.66
71.58
0.47
2.18
0.30 | 0.61
8.37
29.13
14.14
0.16
1.64
0.05 | 1.03
10.56
34.28
53.21
0.44
1.74
0.26 | 0.93
12.71
44.23
21.47
0.24
2.50
0.08 | 1.29
13.97
49.12
66.85
0.62
2.39
0.31 | 0.68
9.31
32.41
15.76
0.18
1.82
0.06 | 1.05
13.54
60.72
63.04
0.52
1.77
0.25 | | | | Traffic-related | Average betweenness (BETW)
Intersection density per road km (INT) | 518.97
8.71 | 1013.62
6.33 | 786.17
8.61 | 2020.76
6.42 | 760.32
9.01 | 1703.99
6.78 | 384.92
6.96 | 853.54
5.51 | | | | Socioeconomic/demographic | % of young (16–24 years) population (YP)
% of old (≥65 years) population (OP)
Average horsepower of cars (HP) | | 2.12
5.97
0.72 | 9.06
23.27
11.78 | 4.43
10.12
2.88 | 9.41
24.54
12.09 | 3.49
11.11
2.13 | 7.96
21.74
10.67 | 5.47
17.23
4.22 | | |

investigation of MAUP effects. The use of census tracts (CTs) of Valencia and a grid of hexagonal BSUs (HEXAs) are two easy-to-implement options that were tested. Particularly, CTs have been investigated in many previous studies (Wier et al., 2009; Abdel-Aty et al., 2013; Cai et al., 2017). Regarding the use of hexagons, these have been recommended over square grids in related literature on traffic safety and MAUP given its more compact shape (Loidl et al., 2016). The scarcity of road network at some areas in the North of Valencia led to join some of the hexagonal units that were covering them, but these modifications were minimal in relation to the whole hexagonal grid.

Furthermore, two more specific BSU schemes were delineated on the basis of two capital elements of any urban traffic network: main roads (segments) and intersections between main roads (points). It is known that these two road entities absorb a substantial percentage of traffic crashes, being the case of road intersections especially treated in literature (Miaou and Lord, 2003; Huang et al., 2017; Lee et al., 2017). Thiessen polygons (also known as Voronoi or Dirichlet polygons) were constructed around points located along main roads of the city and exactly at main intersections, generating two BSU types that hereinafter are referred to as TMs (Thiessen polygons based on main roads) and TIs (Thiessen polygons based on intersections between main roads). Given a collection of locations in a planar space, the Thiessen polygon built from one of these locations, P, contains all the points of the space that are closer to P than to any of the other locations established. Hence, each of the Thiessen polygons defined as a BSU was associated with a particular point along the main road structure (in-between a main road) or to a main intersection of the city. The use of TMs, TIs or HEXAs clearly alleviates the uncertainties derived from crashes located near BSU boundaries, which may have a strong effect in the case of CTs given the historical tendency of defining administrative divisions along main roads, where many crashes occur (Table 2).

Then, Fig. 2 includes the four types of BSU configurations that were

Table 2 Percentage of crashes located near BSU boundaries considering five distance thresholds (5, 10 and 20 m) and mean distance from crashes to BSU boundaries for the four BSU configurations employed in the analysis.

BSU configuration	< 5 m (%)	< 10 m (%)	< 20 m (%)	Mean distance (m)
CTs	35.36	45.29	59.56	27.99
TMs	8.75	18.32	30.23	48.77
TIs	6.65	11.54	22.73	64.78
HEXAs	5.83	11.91	22.40	52.19

defined over the region of study, which provide enough evidence of how the change of the system alters substantially the spatial distribution of traffic crashes across the city. For instance, the central district of Valencia includes several CTs and HEXAs where the crash rate belongs to the highest quintile (Fig. 2a and d), but this effect clearly reduces when the TMs and TIs are considered (Fig. 2b and c).

The number of CTs in Valencia at the beginning of the year 2015 (566) served as a guide in order to define the other three BSU systems in a way they presented a comparable scale (similar number of BSUs). In the case of TIs, the initial scale was conditioned by the number of main intersections in Valencia, rendering it impossible the implementation of a finer BSU scheme of this nature, than that presented in Fig. 2c. Thus, the four baseline configurations in Fig. 2 composed of 566 CTs, 574 TMs, 378 TIs and 515 HEXAs were chosen to analyse the MAUP effect in the modelling of traffic crash counts for the available dataset.

3. Methodology

3.1. Software

The R programming language (3.5.1 version, R Development Core Team, Vienna, Austria) (R Core Team, 2018) was used to obtain all the results presented in this work. The R packages ClustGeo (Chavent et al., 2017b), ggplot2 (Wickham, 2016), INLA (Rue et al., 2009; Martins et al., 2013; Lindgren and Rue, 2015), rgeos (Bivand and Rundel, 2018), spatstat (Baddeley et al., 2015), spded (Bivand and Piras, 2015), spgwr (Bivand and Yu, 2017) and SpNetPrep (Briz-Redón, 2019) were specifically required for performing the analysis.

3.2. Regionalization algorithm

The term regionalization was defined by Guo (2008) as the process of aggregating a set of spatial entities into a reduced number of regions in a way that a predefined objective function is optimized. There are several important regionalization algorithms, including SKATER (Assunção et al., 2006), REDCAP (Guo, 2008) and ClustGeo (Chavent et al., 2017b). In this paper the latter was chosen, which is implemented in the R package ClustGeo (Chavent et al., 2017a). The next paragraphs contain a brief description of how this method works and how it was used.

Given a number of clusters, K, to be formed and two matrices, D⁰ and D1, that represent the homogeneity and physical distances (respectively) between the spatial units available before regionalization, the ClustGeo algorithm relies on the minimization of a measure called

Fig. 2. Crash counts (CRASH) in Valencia considering a BSU configuration composed of CTs (a), TMs (b), TIs (c) and HEXAs (d). Districts of Valencia are overlayed (thicker lines, in black) for better readability and comparison.

mixed within-cluster inertia, defined as the sum of the mixed inertias of all of the clusters established. The mixed inertia of a cluster, Ck ^α, follows the next expression (Chavent et al., 2017b):

$$I_{\alpha}(Ck^{\alpha}) = (1 - \alpha) \sum{i \in Ck^{\alpha}} \sum{j \in C_k^{\alpha}} \frac{w_i w_j}{2\muk^{\alpha}} d{0,ij}^2 + \alpha \sum_{i \in Ck^{\alpha}} \sum{j \in C_k^{\alpha}} \frac{w_i w_j}{2\muk^{\alpha}} d{1,ij}^2$$

where α ∈ [0, 1] is a parameter that controls the importance that the homogeneity and physical distances (represented by D⁰ and D1) have in the clustering procedure, k ∈ {1, …, K} is the index for the cluster, wi is the weight of spatial unit ⁱ, ⁼ ^∑ ^∈ μk ^w ^α i C ⁱ ^k ^α and d0,ij (resp. d1,ij) is the normalized dissimilarity between spatial units i and j in D⁰ (resp. D1).

In our paper, the total number of crash counts registered per BSU during the period 2014–2015 at four time slots (23h–7h, 7h–14h, 14h–20h and 20h–23h, which were selected according to the daily trends observable in Fig. 1b) and at the weekends were used to define the dissimilarity matrix D0. Regarding D1, this matrix was constructed from the Euclidean distances between the centroids of the BSUs. Furthermore, the weights (wi) were set equal for all BSUs and a value of α = 0.1 was chosen, giving much importance to the spatial distances between the BSUs during the aggregation procedure (the investigation of the optimal value of α suggested by ClustGeo led to this choice).

It needs to be remarked that there is a technical difference between SKATER and REDCAP algorithms and the method implemented in ClustGeo. Indeed, the choice of K in ClustGeo does not represent the number of contiguous and homogeneous regions that are created, but the number of homogeneous regions (according to the variables provided to the algorithm) that need to be regrouped later in order to fully satisfy the contiguity constraints. Hence, the input K in ClustGeo is a lower bound of the number of BSUs that are generated at the end of the process, although both values barely differ. The use of several values of K, from 100 to 500 in intervals of 100, allowed MAUP to be investigated in the present study with five different levels of spatial aggregation, which are denoted by AG100, AG200, AG300, AG400 and AG500 within the rest of the paper.

3.3. Crash counts modelling

3.3.1. Conditional autoregressive model

The modelling of crash counts at the macroscopic level requires the consideration of data overdispersion. Two common choices in the field of traffic safety analysis to address this issue are negative binomial (also known as Poisson–Gamma) and Poisson lognormal probability distributions (Lord and Mannering, 2010). Both have their own advantages and disadvantages, Poisson lognormal being more recommended in cases of high overdispersion (particularly skewed distribution of the counts), whereas negative binomial has been suggested to be more suitable for moderately overdispersed counts, and also for counts with a large number of zeros (Khazraee et al., 2018; Shirazi and Lord, 2018). In the context of this study, it is hard to choose between one distribution or the other, as the change in scale or zoning alters the statistical properties that are involved in our decision. Anyhow, as the crash counts available under the different combinations of aggregation level and BSU type were overall only moderately overdispersed, we decided to select the negative binomial modelling approach.

Therefore, a conditional autoregressive (CAR) model with negative binomial (NB) response was chosen to fit the crash counts recorded for each BSU. The use of a CAR structure for modelling crash counts is a usual strategy in traffic safety analysis to account for spatial heterogeneity (Quddus, 2008; Huang et al., 2010).

If Y ∼ NB(μ, ψ) (basic NB distribution of mean μ and shape ψ) then it holds that E(Y) = μ, V ( ) Y μ = + ^μ and P

= = + − −++ ( ) ( )( ) ( ) Y x x ψ ψ μ ψ ψ μ μ ψ x 1 ¹ . Then, assuming a NB distribution for the response (crash counts) the following spatial model was implemented:

$$Y_i \sim NB(\mu_i, \psi)$$

$$\log(\mu_i) = \log(E_i) + \beta0 + \sum{m=1}^{p} \betam X{im} + \phi_i$$ (1)

where Yi is the number of crashes observed at BSU i, μⁱ and ψ are, respectively, the mean risk (for BSU i) and overdispersion (1/ψ) values for the NB distribution, the natural logarithm acts as a link function for μi, Ei (exposure at BSU i) is the length of non-pedestrian road at BSU i which acts as an offset of the equation, Xim represents the value of the mth covariate at BSU i, β^m is the coefficient that controls the effect of the mth covariate and ϕⁱ represents a spatial effect for BSU i. Regarding the selection of the exposure, the unavailability of vehicle miles travelled data (traffic volume) for non-main roads of Valencia left nonpedestrian road length as the natural choice, a possibility already considered in previous research studies (Qin et al., 2004; Imprialou et al., 2016).

The spatial effect in Eq. (1) was modelled using the following CAR structure (Besag, 1974; Besag et al., 1991):

$$\phi_i | \phij, j \neq i \sim N \left( \sum{j=1}^n w_{ij} \phi_j, \tau_i^{-1} \right)$$

where *w*ij is an indicator parameter that is 1 if BSUs i and j are contiguous and 0 otherwise, and τⁱ is a precision parameter that varies with BSU i.

3.3.2. Geographically weighted regression

Geographically weighted regression (GWR) is a form of linear regression that captures the spatial heterogeneity present in the data by allowing model parameters to vary locally (Brunsdon et al., 1996; Fotheringham et al., 2002; Nakaya et al., 2005). GWR has already been used in traffic safety analysis (Hadayeghi et al., 2010; Matkan et al., 2011; Xu and Huang, 2015; Gomes et al., 2017), including some analyses from the perspective of MAUP effects (Amoh-Gyimah et al., 2017).

The mathematical expression that corresponds to the GWR model is the following:

$$\log(\mu_i) = \log(E_i) + \beta_0(BSUi) + \sum{m=1}^{p} \beta_m(BSUi)X{lim}$$ (2)

where μi, Ei and Xim are as in Eq. (1). The main feature of GWR is the consideration of local regression parameters (in contrast to global parameters of Eq. (1)) which are denoted by βm(BSUi) in Eq. (2). As in Eq. (1), a NB distribution was used in the definition of the model to consider overdispersion. A modification of GWR called semiparametric GWR consisting in the combination of fixed and spatially-varying effects for the covariates involved in the model has also been used in traffic safety analysis (Xu and Huang, 2015; Amoh-Gyimah et al., 2017). However, we decided to stay with the classical version of the GWR model in order to provide a more unified framework for the comparison of the set of models obtained for each aggregation level and zoning, which is the main purpose of this work.

Hence, a GWR model behaves similarly to a generalized linear model (GLM), although for the former the parameters that compose the model are estimated locally, at each BSU, depending on the crash counts and covariate values at the surrounding areal units. The influence that BSU i produces on another BSU j (denoted as *w*ij) was controlled by the following Gaussian kernel function:

$$w{ij} = e^{-0.5(d{ij}^2/\sigma^2)}$$

where dij is the Euclidean distance between BSU i and BSU j (between their centroids) and σ is the fixed bandwidth employed by the kernel function, which represents the level of influence that the rest of BSUs have on a given BSU with regard to model fitting (a higher value for σ means that model parameters are estimated on the basis of a wider zone around each BSU). Several other kernel functions are available instead of the Gaussian (bisquare, for instance), but this choice is usually not responsible of strong effects on the results (Silverman, 2018).

Regarding the bandwidth, a value of σ = 2 km was chosen in this study for all the BSUs and aggregation levels being considered. Other authors opted for the choice of a specific optimal bandwidth for each BSU and aggregation level (Amoh-Gyimah et al., 2017), but here a fixed value was used in order to avoid the presence of a source of variation other than scale or zoning, which are the focus of the analysis. The value of 2 km was chosen because it was close to the optimal values that were observed for the different BSUs and aggregation levels tested.

3.4. Assessment of model performance

The goodness of fit of the CAR models was assessed through Bayesian deviance information criterion (DIC) (Spiegelhalter et al., 2002). Similarly, AIC was used for GWR models. In addition, several measurements of model performance typically used in other traffic safety analysis papers on the MAUP for model comparison were considered for both CAR and GWR models: mean absolute deviation (MAD) (Lee et al., 2014; Xu et al., 2014; Amoh-Gyimah et al., 2017; Zhai et al., 2018b), sum of absolute deviation (SAD) (Lee et al., 2014; Xu et al., 2018; Zhai et al., 2018b) and percent mean absolute deviation (PMAD) (Lee et al., 2014; Xu et al., 2018).

The formulas for MAD, SAD and PMAD are the following:

$$MAD = \frac{1}{n} \sum_{i=1}^{n} |e_i|$$

$$SAD = \sum_{i=1}^{n} |e_i|$$

$$PMAD = \frac{\sum{i=1}^{n} |e{i}|}{\sum{i=1}^{n} |y{i}|}$$

where ⁿ is the number of BSUs and eyy = − ^ˆ ⁱ i i represents the difference between the number of crashes observed at BSU i (yi) and the number fitted by the model (*y*ˆ i).

From the perspective of interpreting the results, a lower value of any of the aforementioned statistics (DIC, AIC, MAD, SAD or PMAD) indicates a better fit to the available data.

3.5. Statistical tools for covariate exploration

Several statistical tools were used to explore the covariates provided to the models at different scales and zonings, and hence provide more instruments to analyse their sensitivity to the MAUP. This section includes a brief description of these tools.

Average nearest neighbour index (NNI) measures the level of clustering/dispersion of a point pattern. Hence, it is suitable for the exploration of a covariate constructed from a pattern of points located across the area of investigation (EDU, BAR, COMP, HOT, PARK, BUS and TRAM among the set of covariates used in the present research). The definition of NNI is the following (Clark and Evans, 1954; Cressie, 1993):

$$NNI = \frac{(1/P) \sum{i=1}^{P} D{NN}(i)}{(1/2) \sqrt{A/P}}$$

where P is the number of points that form the pattern, A is the area of the whole space where the pattern lies and DNN(i) is the distance from point i to its nearest neighbour (the closest point to i in the pattern). The NNI represents a ratio between the average nearest-neighbour distance observed for the pattern and the value that would be expected under the hypothesis of random spatial distribution. A NNI lower than 1 indicates that the pattern is clustered, whereas a value higher than 1 is a sign of the dispersion of the pattern.

Moran's I (Moran, 1950a,b) was computed for every combination of BSU, aggregation level and covariate available. Moran's I is defined as follows:

$$I = \frac{\sum{i=1}^{n} \sum{j \in N_{\text{dir}}(i)} (1/n_i)(x_i - \bar{x})(xj - \bar{x})}{\sum{i=1}^{n} (x_i - \bar{x})^2}$$

where xi is the value of a covariate at BSU i and *x*¯ the mean value of the covariate across all spatial units available. Moran's I behaves as a spatial autocorrelation coefficient for areal-based data. Under the hypothesis of no spatial autocorrelation, it holds that E(I) = −1/(n − 1), where n is the number of spatial units in each case. A higher Moran's I value indicates a higher tendency of the covariate to show strongly associated values for neighbouring BSUs.

Multicollinearity among the covariates considered was investigated through the variance inflation factor (VIF) (Fox, 1991), which is defined as follows for a given covariate or predictor, Xj:

$$VIF_j = \frac{1}{1 - R_j^2}$$

where Rj ² is the R² found when regressing all other covariates onto Xj (Miles, 2014). A higher value of VIF suggests that the covariate is more susceptible to lead to multicollinearity issues. A value of VIF not greater than 10 is usually interpreted in literature as a sign of no severe multicollinearity (Miles, 2014).

4. Results and discussion

Tables 5 8 include the parameters estimated for the CAR models considering each BSU configuration and aggregation level. On the other hand, Figs. 6–9 display the densities (distributions) of the local parameter estimates obtained from the GWR models for each BSU type, aggregation level and covariate. The distributions of these local parameters were scaled (divided by their standard deviation) to facilitate the graphical comparison.

The main conclusion that yields from all these results is that MAUP effects have heavily affected the macroscopic traffic safety analysis performed. In the present section, an initial subsection gives insights into the varying nature of the data (response and covariates) as one changes scale or zoning. The subsequent subsections contain a description of how the variations in aggregation level (scale) or BSU type (zoning) changed model parameter estimates and a performance comparison of all the models implemented. After that, the associations derived from all the models and BSU configurations tested are globally evaluated in pursuit of more solid conclusions (despite the consequences of MAUP) that allow the identification of some factors that correlate with more traffic crashes.

4.1. Effects of MAUP on input data

Before analysing scale and zoning effects on parameter estimations considering both CAR and GWR models, an investigation of the consequences of MAUP on the input data that is afterwards modelled (response and covariates) was performed.

The distributional characteristics of the response (crash counts) for different BSUs and aggregation levels are shown in Table 4, from which it yields that overdispersion (through the coefficient of variation) and kurtosis reduced as the level of aggregation was increased. Furthermore, it can also be observed in Table 4 that the percentage of zeros was almost negligible for most of the combinations between a BSU and an aggregation level, excluding the HEXAs.

The magnitude of the spatial autocorrelation shown by the response variable and the covariates also suffers from scale and zoning changes. Fig. 3 displays all Moran's I indexes computed for each combination of a BSU and aggregation level, which suggests that it is hard to predict the level of spatial autocorrelation after a change of scale. Indeed, some covariates tend to be more spatially autocorrelated as aggregation increases (the number of parking zones, for TMs), whereas other covariates show the opposite behaviour (the number of companies, for HEXAs). In addition, CTs and HEXAs show overall higher levels of spatial autocorrelation for the covariates than TMs and TIs. This is remarkable for traffic crashes, which display a particularly high spatial autocorrelation in the case of HEXAs, rather than in the other three BSUs considered for investigation. The spatial autocorrelation of some of the covariates is more deeply discussed in the following subsections.

Finally, the computation of variance inflation factors (VIF) leads to the conclusion that multicollinearity issues were not present for the distinct dataset analysed (at different scales and zonings), as VIF factors were always below 10 (Fig. 4). However, it is important to appreciate that VIF consistently increased as the level of aggregation increased, specially for some covariates such as the number of bars/restaurants (BAR), the number of companies (COMP) and intersection density (INT). Hence, our analysis suggests that the level of aggregation should not be excessively increased in order to avoid multicollinearity among the covariates.

4.2. Parameter variations across aggregation levels

Given the moderate level of significance achieved by the set of covariates, an 80% credibility level was also considered along with the most usual 90% level for the estimations yielded by the CAR models. Fig. 5 displays a graphical summary of the significance achieved by all the covariates involved in the analysis. According to Fig. 5, the hypothesis of obtaining more significant variables and the consequent higher interpretation power at lower levels of aggregation suggested by Xu et al. (2014) was true in the case of CTs, but it was not clear, at all, for the rest of BSU configurations. Anyhow, this question could have been better addressed in the presence of a higher number of significant covariates.

The level of aggregation applied to each BSU configuration through the regionalization algorithm produced moderate-to-severe effects in parameter estimates for the CAR models. Hence, although the parameter estimates evolved moderately with changes in the scale (Tables 5 –8), several covariates were only significant at some of the aggregation levels tested. However, some of the covariates did not seem affected by MAUP and remained significant with each aggregation level (old population percentage and average horsepower for CTs, and number of educational centres and average horsepower for both TMs and TIs). Despite not being significant for the most aggregated scheme considered (AG100), the positive association between traffic crashes and the number of bus stops for HEXAs was also consistent. Remarkably, none of the covariates that were found significant (at 80% of credibility) experimented a change of effect (from positive to negative, or vice versa) after a shift in the scale. This is a positive result, since it indicates that MAUP effects were not the strongest possible across aggregation levels.

Regarding the GWR models, Figs. 6–9 show that local parameter distributions may vary acutely after some changes on the aggregation level. It is hard to assess if the distribution of the local parameters tends to be more concentrated (leptokurtic) or flat (platykurtic) as the level of aggregation increases/decreases, as this seems strongly dependent on the covariate and the BSU. Furthermore, the contradictory presence of local parameters of opposite signs that takes place for most of the covariates is a well-known issue that often arises in GWR models (Hadayeghi et al., 2010; Xu and Huang, 2015; Amoh-Gyimah et al., 2017). Fig. 10 shows the behaviour of the local parameter estimates obtained from the GWR models through the signs of 5th, 10th, 20th, 80th, 90th and 95th percentiles. Hence, a negative value for 80th, 90th or 95th percentiles indicates a high agreement among the local GWR coefficients and a negative association of the covariate with crash

Fig. 3. Moran's I values for the response (CRASH) and the covariates for each BSU type and aggregation level tested.

Fig. 4. Assessment of multicollinearity among the covariates through variance inflation factor (VIF) for each BSU type and aggregation level tested.

counts. Analogously, a positive 5th, 10th or 20th percentiles means the same but for a positive association.

In contrast to CAR models, for which no covariate showed a significant change of effect after a variation in scale, the GWR models experimented this issue for the covariates representing the number of companies (COMP), parking zones (PARK) and intersection density (INT) when using CTs (considering the percentile-based criteria that has been employed for assessing the association between crash counts and covariates with the GWR models). This lack of coherence was specifically serious for the number of companies, a covariate that showed a strong autocorrelation according to Moran's I at this BSU system (Fig. 3). The exploration of crash counts and the number of companies (COMP) at AG100, AG300 and AG500 (Fig. 12) unveils some singular patterns around the city centre (some surrounding Districts are highlighted in blue in Fig. 12). Thus, whereas at AG100 most of the BSUs presenting a high number of companies were located within these Districts, the use of a more disaggregated configuration provided greater variation in the number of companies across the whole city,

Fig. 5. Summary of the results obtained for the CAR models, considering the four types of BSUs and the levels of aggregation that were applied.

with many more BSUs in the periphery of Valencia presenting high values. With regard to intersection density (INT), this covariate also presented high Moran's I values (Fig. 3), but visual inspection was far less clear than in the case of the number of companies, becoming challenging to figure out how the effect of intersections changed from AG300 to AG400 and again from AG400 to AG500 (for CTs and the GWR models). On the other hand, it is remarkable how the number of educational centres (EDU), which presents the more coherent behaviour across BSUs and aggregation levels in the case of the GWR models, is one of the covariates that showed a lower range of values for Moran's I statistic. Similarly, this result also agrees with that provided by the NNI (Table 3), as some of the covariates based on point patterns lying over the city presenting a high level of clustering (companies, NNI = 0.26) display a more sensitive to MAUP behaviour than other that, albeit clustered, show a more regular pattern (educational centres, NNI = 0.86). A similar level of consistence to that found for the

Fig. 6. Combined graph showing the distributions of local parameter estimates, for the covariates used in the GWR models (in rows) and each level of spatial aggregation (in columns) tested for the CTs.

Fig. 7. Combined graph showing the distributions of local parameter estimates, for the covariates used in the GWR models (in rows) and each level of spatial aggregation (in columns) tested for the TMs.

educational centres was also obtained for the number of tram stops (TRAM) in the case of the GWR models, although this covariate resulted non-significant for almost all combinations of aggregation level and BSU for the CAR models. The NNI of the point pattern formed by the tram stops was 1.30 (Table 3), clearly indicating the dispersed configuration of these stops across Valencia.

4.3. Parameter variations across BSU types

Contrary to scale, MAUP effects from zoning variations were by far more severe in this case study. Indeed, some covariates showed a significant and opposite effect depending on the BSU system being considered, including old population percentage, the number of bus stops and intersection density (Fig. 5). The cases of both the number of bus

Fig. 8. Combined graph showing the distributions of local parameter estimates, for the covariates used in the GWR models (in rows) and each level of spatial aggregation (in columns) tested for the TIs.

Fig. 9. Combined graph showing the distributions of local parameter estimates, for the covariates used in the GWR models (in rows) and each level of spatial aggregation (in columns) tested for the HEXAs.

stops and intersection density were specifically related to the differential behaviour exhibited by the highest aggregation level, AG100. This aggregation level was clearly the least coherent among all the levels tested, possibly indicating its unsuitability to capture some microvariations present in the data. On the other hand, the percentage of old population (OP) appeared as a highly sensitive-to-MAUP covariate, standing out from all the ones supplied to the models. Thus, whereas this covariate showed a clear positive association with traffic crashes

Table 3 Nearest neighbour indexes (NNI) and p-value associated with each index.

	EDU	BAR	COMP	HOT	PARK	BUS	TRAM
NNI	0.86	0.51	0.26	0.35	1.08	0.70	1.30
p-value	0.00*	0.00*	0.00*	0.00*	0.15	0.00*	0.00*

* An asterisk indicates the statistical significance of the index (p < 0.05), which indicates clustering (NNI < 1) or dispersion (NNI > 1).

Fig. 10. Summary of the results obtained for the GWR models, considering the four types of BSUs and the levels of aggregation that were applied.

Table 4 Basic distributional properties of the response variable (crash counts) corresponding to each BSU type and aggregation level, where CV means coefficient of variation (standard deviation to mean ratio).

BSU	AG	CV	Kurtosis	Zeros (%)
CT	AG100	0.71	1.77	0.00
CT	AG200	0.87	3.41	0.00
CT	AG300	0.88	4.71	0.00
CT	AG400	0.92	10.86	0.25
CT	AG500	1.04	16.07	0.40
TM	AG100	0.54	1.76	0.00
TM	AG200	0.72	0.59	0.00
TM	AG300	0.76	2.11	0.33
TM	AG400	0.79	7.07	0.50
TM	AG500	0.88	12.03	0.60
TI	AG100	0.69	1.65	0.00
TI	AG200	0.77	1.49	0.00
TI	AG300	0.81	6.77	0.00
HEXA	AG100	0.85	2.32	2.68
HEXA	AG200	0.83	−0.15	6.50
HEXA	AG300	0.88	1.63	7.00
HEXA	AG400	0.84	3.75	5.75
HEXA	AG500	1.01	4.70	10.00

for CTs and, to a lesser extent, for TIs, it associated with a decrease in crash counts with HEXA units at the three most aggregated levels. This contradictory result was investigated through the cartographic representation of crash counts (CRASH) and the old population covariate for CTs and HEXAs at AG100, AG200 and AG300 (Fig. 13). In all the maps available in Fig. 13, the border of a census tract located in the South of Valencia (which is the largest of the city) is highlighted in blue. This census tract constitutes a wide area of low population density and a high percentage of residents with 65 or more years of age (OP). In addition, the area is not dense in road network, which naturally reduces the number of traffic crashes. Hence, the use of HEXAs led to a covering of this census tract with several hexagons of very high percentage of old population and very low number of crashes, which surely affected the estimation of the parameter related to this covariate. On the other hand, the use of CTs summarizes this part of the city in only one area presenting a high percentage of old population and moderate value of crash counts, which can barely alter model estimations. In conclusion, the use of covariates that depend on possibly sparse population should be used with special care, as the choice of the wrong BSU type in such cases could lead to artefactual associations between the covariate and crash rates observed.

The differences in local parameter estimates for the GWR models across the four BSU types is rather evident (Fig. 10). Several covariates presented a controversial behaviour with a strong dependence to the BSU system. These discrepancies are more obvious than with the CAR models given the higher number of covariates that are highlighted with the specified percentile criteria. Leaving apart the AG100 aggregation level because it has overall produced more disparate results (reducing its reliability), there were still some covariates showing inconsistent associations with crash counts, confirming the consequences of MAUP in this case study. Furthermore, a global high level of coincidence between two distributions of local parameter estimates derived from GWR (at two different scales and/or zonings) does not guarantee, at all, that the local estimates vary similarly across space, which is clear in view of the examples shown in Fig. 11.

4.4. Model performance comparisons

Table 9 provides information with regard to model fitting for all the BSU systems and aggregation levels employed in the study. It is appreciable that CAR models performance improved gradually (decrease in DIC) as the aggregation increased (reaching the minimums at AG100 for all the BSU systems). Although this is an issue already pointed out by Fotheringham and Wong (1991), that does not always represent a real improvement in model quality and interpretation, in our case may be also the consequence of a weak multicollinearity among the covariates at AG100 (specially for HEXAs).

On the other hand, for any fixed scale with the exception of AG100, HEXAs appeared as the optimal choice for the CAR models. Hence, the use of hexagonal grids would be a reasonable recommendation, although care must be taken with areas of low population density if population-related covariates are being used, as shown in the previous subsection. For the latter, CTs or other administrative division should be more convenient.

One positive conclusion is the substantial level of agreement shown by the CAR and GWR models for each BSU system and level of aggregation, as it can be observed from the comparison of Figs. 5 and 10 . However, GWR models showed superior values for MAD, SAD and PMAD for most of the combinations of BSU types and aggregation levels, an opposite result to that found by other authors (Xu and Huang, 2015; Amoh-Gyimah et al., 2017). The use of semiparametric GWR models or an adaptive version of their kernel's bandwidth may have led to more close performance results for some combinations of scale and zoning, but this possibility was discarded to guarantee a fair comparison between models in the context of our analysis, which is more focused on model parameter estimations rather than on model performances.

Table 5 Estimates with standard deviation (SD) for the parameters involved in the CAR model, considering CTs as BSUs.

Covariate	AG100		AG200		AG300		AG400		AG500
	Est.	SD	Est.	SD	Est.	SD	Est.	SD	Est.	SD
Intercept	−13.2477*	2.6119	−11.4619*	1.8909	−8.8597*	1.6753	−8.0273*	1.1308	−8.1816*	1.0403
EDU	−0.0327	0.0261	−0.0661*	0.0260	−0.0797*	0.0285	−0.0778*	0.0303	−0.0656*	0.0308
BAR	0.0025	0.0027	0.0030	0.0027	0.0021	0.0033	−0.0017	0.0035	0.0004	0.0042
COMP	−0.0005	0.0006	−0.0005	0.0006	0.0000	0.0009	0.0002	0.0009	0.0002	0.0009
HOT	0.0002	0.0004	0.0002	0.0003	0.0003	0.0004	0.0001	0.0004	0.0001	0.0004
PARK	−0.0230	0.0657	−0.0970	0.0607	−0.0187	0.0665	0.0102	0.0659	−0.0249	0.0661
BUS	0.0006	0.0134	0.0224*	0.0121	0.0186	0.0141	0.0156	0.0150	0.0222	0.0155
TRAM	−0.0436	0.0949	0.0424	0.0763	−0.0065	0.0924	−0.0711	0.0906	−0.0586	0.0911
BETW	0.0001	0.0001	0.0001	0.0001	0.0001	0.0000	0.0001	0.0000	0.0001*	0.0000
INT	−0.0130	0.0187	−0.0053	0.0122	−0.0010	0.0093	0.0095	0.0076	0.0051	0.0065
YP	−0.0265	0.0523	−0.0019	0.0313	0.0260	0.0245	0.0122	0.0206	0.0170	0.0176
OP	0.0322*	0.0151	0.0156	0.0096	0.0161*	0.0085	0.0120*	0.0072	0.0131*	0.0063
HP	0.7517*	0.2087	0.6004*	0.1487	0.3605*	0.1355	0.3111*	0.0921	0.3169*	0.0842
ψ	3.0831*	0.4127	8.8277*	2.9744	6.5329*	1.8981	5.4961*	1.3313	5.2406*	1.0798

^* An asterisk indicates the significance of a parameter with a 90% credibility.

Table 6 Estimates with standard deviation (SD) for the parameters involved in the CAR model, considering TMs as BSUs.

Covariate	AG100 AG200			AG300		AG400			AG500
	Est.	SD	Est.	SD	Est.	SD	Est.	SD	Est.	SD
Intercept	−8.3199*	1.2364	−7.7404*	1.0021	−7.2087*	0.7263	−6.7423*	0.6376	−6.0589*	0.5407
EDU	−0.0475*	0.0270	−0.0515*	0.0279	−0.0732*	0.0316	−0.0963*	0.0325	−0.0967*	0.0324
BAR	0.0019	0.0026	0.0029	0.0028	−0.0002	0.0034	−0.0016	0.0035	−0.0031	0.0038
COMP	−0.0002	0.0006	0.0001	0.0009	−0.0002	0.0010	0.0002	0.0011	−0.0001	0.0013
HOT	−0.0001	0.0004	0.0007	0.0005	0.0009	0.0006	0.0008	0.0006	0.0010*	0.0006
PARK	−0.0126	0.0594	−0.0397	0.0625	−0.0185	0.0800	−0.0870	0.0799	−0.0357	0.0788
BUS	−0.0380*	0.0118	−0.0010	0.0140	−0.0039	0.0163	−0.0086	0.0193	0.0079	0.0208
TRAM	−0.1556*	0.0792	0.0021	0.0988	−0.0631	0.1141	−0.1625	0.1170	−0.1163	0.1138
BETW	0.0000	0.0001	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
INT	−0.0033	0.0164	0.0081	0.0107	0.0115	0.0098	0.0000	0.0075	−0.0006	0.0066
YP	0.0082	0.0398	0.0365	0.0268	0.0042	0.0195	0.0084	0.0165	0.0061	0.0143
OP	0.0339*	0.0122	−0.0039	0.0083	0.0026	0.0059	0.0038	0.0052	0.0021	0.0047
HP	0.3635*	0.0993	0.3028*	0.0755	0.2849*	0.0568	0.2565*	0.0501	0.1996*	0.0428
ψ	3.8533*	0.5287	8.4472*	1.7932	3.1448*	0.5569	2.7796*	0.3599	2.9332*	0.3634

^* An asterisk indicates the significance of a parameter with a 90% credibility.

Table 7 Estimates with standard deviation (SD) for the parameters involved in the CAR model, considering TIs as BSUs.

| Covariate | AG100 | | AG200 | | AG300 | | | |-----------|------------|--------|----------|--------|----------|--------|--| | | Est.
SD | | Est. | SD | Est. | SD | | | Intercept | −6.7825* | 1.5880 | −6.8832* | 0.9019 | −6.2381* | 0.6777 | | | EDU | −0.0399 | 0.0281 | −0.0633* | 0.0330 | −0.0580* | 0.0298 | | | BAR | 0.0028 | 0.0028 | 0.0022 | 0.0034 | 0.0039 | 0.0035 | | | COMP | 0.0007 | 0.0007 | 0.0012 | 0.0009 | −0.0004 | 0.0009 | | | HOT | −0.0002 | 0.0005 | 0.0002 | 0.0006 | 0.0003 | 0.0005 | | | PARK | −0.0133 | 0.0632 | −0.0286 | 0.0724 | −0.0143 | 0.0638 | | | BUS | −0.0182 | 0.0128 | −0.0123 | 0.0167 | 0.0002 | 0.0161 | | | TRAM | −0.0962 | 0.0754 | −0.2074* | 0.1079 | −0.1248 | 0.1088 | | | BETW | −0.0001 | 0.0001 | 0.0000 | 0.0001 | 0.0000 | 0.0000 | | | INT | −0.0126 | 0.0178 | −0.0052 | 0.0125 | −0.0126 | 0.0085 | | | YP | −0.0452 | 0.0485 | 0.0329* | 0.0193 | 0.0226 | 0.0142 | | | OP | 0.0207 | 0.0145 | 0.0138* | 0.0074 | 0.0010 | 0.0043 | | | HP | 0.2806* | 0.1284 | 0.2234* | 0.0604 | 0.1994* | 0.0471 | | | ψ | 3.2962* | 0.4602 | 2.3784* | 0.2386 | 7.0513* | 2.1572 | |

^* An asterisk indicates the significance of a parameter with a 90% credibility.

4.5. Parameter interpretations

Despite MAUP effects, several associations between crash counts and some of the covariates were observed at several combinations of aggregation level and BSU type, deserving a deeper analysis. In particular, the number of undergraduate educational centres showed a consistent negative correlation with crash counts, whereas average horsepower of the cars in the BSU generally associated with more traffic crashes. Other covariates, such as the number of bus stops in the BSU, the average betweenness of its road segments or the percentage of population with 65 or more years living in the BSU also suggested the presence of a positive relationship, but uncertainties from MAUP issues were stronger for them. The knowledge of the city being analysed arouses the suspicion that some of these correlations could be related to a hidden (not included in the models) factor as it is the distance to the city centre of Valencia (case of average horsepower and old population), but this question would require a specific research.

5. Conclusions

This paper is, to the best of our knowledge, the first one that provides a simultaneous investigation of scale and zoning effects regarding the modifiable areal unit problem in the context of traffic safety analysis. Furthermore, another capital objective was to specifically assess how a change in the aggregation level or BSU type may affect the basic characteristics of both the response variable being considered (crash counts) and the set of covariates included in the models. The consequences of MAUP for the data analysed were notorious from the perspective of both scale and (specially) zoning alterations. Some of the effects of MAUP were understandable from visual inspection of the data, as shown through some exemplifications, but it is really tough

Table 8 Estimates with standard deviation (SD) for the parameters involved in the CAR model, considering HEXAs as BSUs.

| Covariate | AG100 | | AG200 | | | AG300 | | AG400 | | AG500 | | |-----------|----------|--------|------------|-----------|-----------|----------|----------|----------|-----------|----------|--| | | Est. | SD | Est. | SD | Est. | SD | Est. | SD | Est. | SD | | | Intercept | −8.8001* | 1.3441 | −4.8880* | 0.7877 | −5.2212* | 0.6348 | −4.3829* | 0.4955 | −5.1725* | 0.3578 | | | EDU | −0.0457 | 0.0396 | −0.0470 | 0.0299 | −0.0759* | 0.0311 | −0.0814* | 0.0339 | −0.0435 | 0.0334 | | | BAR | 0.0001 | 0.0035 | −0.0026 | 0.0030 | −0.0010 | 0.0037 | −0.0008 | 0.0038 | 0.0012 | 0.0036 | | | COMP | −0.0011 | 0.0010 | −0.0002 | 0.0009 | −0.0001 | 0.0009 | 0.0011 | 0.0011 | 0.0011 | 0.0010 | | | HOT | 0.0004 | 0.0006 | 0.0004 | 0.0007 | 0.0002 | 0.0006 | 0.0006 | 0.0006 | 0.0005 | 0.0006 | | | PARK | 0.1677 | 0.1075 | 0.0371 | 0.0781 | 0.0741 | 0.0796 | 0.0406 | 0.0805 | 0.0068 | 0.0766 | | | BUS | 0.0267 | 0.0231 | 0.0816* | 0.0204 | 0.0902* | 0.0212 | 0.0934* | 0.0229 | 0.1336* | 0.0222 | | | TRAM | −0.1348 | 0.1352 | 0.1095 | 0.1633 | −0.0518 | 0.1356 | −0.0057 | 0.1233 | 0.0244 | 0.1290 | | | BETW | 0.0001 | 0.0002 | 0.0001 | 0.0001 | 0.0001 | 0.0001 | 0.0001 | 0.0001 | 0.0001 | 0.0000 | | | INT | 0.0908* | 0.0409 | 0.0151 | 0.0199 | 0.0097 | 0.0153 | 0.0021 | 0.0102 | −0.0012 | 0.0090 | | | YP | 0.0455 | 0.0397 | 0.0036 | 0.0261 | 0.0079 | 0.0126 | 0.0175 | 0.0121 | 0.0021 | 0.0095 | | | OP | −0.0199* | 0.0103 | −0.0127* | 0.0075 | −0.0096* | 0.0056 | −0.0038 | 0.0047 | −0.0033 | 0.0037 | | | HP | 0.3438* | 0.0971 | 0.0395 | 0.0520 | 0.0703 | 0.0464 | −0.0055 | 0.0347 | 0.0480* | 0.0274 | | | ψ | 1.4053* | 0.1862 | 1133.6054* | 9922.7680 | 205.1350* | 579.7774 | 79.1317* | 100.2495 | 103.7747* | 162.6996 | |

^* An asterisk indicates the significance of a parameter with a 90% credibility.

Fig. 11. Estimated local parameters for the GWR model considering CTs, TMs, TIs and HEXAs for EDU (a–d), BAR (e–h) and BUS (i–l) at AG300. Districts of Valencia are overlayed (thicker lines, in black) for better readability and comparison.

Fig. 12. Crash counts (a–c) and COMP values (d–f) for CTs at AG100, AG300 and AG500 (in order of appearance at each row, from more to less aggregated). Some central Districts of Valencia are highlighted in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

sometimes to explain certain model parameter disagreements that arise from a change in aggregation level or BSU configuration.

The comparative analysis yielded that CAR models using hexagonal gridded units (HEXAs) were the best choice according to the performance measurements adopted. The employment of BSU types based on the road network being analysed (TMs and TIs), which naturally avoided boundary effects (although HEXAs even outperformed them in this aspect), did not lead to better model performances. Anyhow, model performance measures should not be the only instrument to select one combination of scale and zoning over others. Indeed, the use of CAR models and HEXAs also unveiled controversial behaviours for some model parameter estimates. These were found to be a consequence of the fact of using a BSU type (hexagonal unit) that may not be the best one to represent demographic characteristics of the area of

Fig. 13. Crash counts for CTs (a–c) and HEXAs (d–f) and OP values for CTs (g–i) and HEXAs (j–l) at AG100, AG200 and AG300 (in order of appearance at each row, from more to less aggregated). A CT in the South of Valencia is highlighted in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

investigation. Thus, even though the results for census tracts were more modest in terms of model performance, this kind of administrative unit is possibly the most suitable one to seek more robust conclusions if several demographic covariates are present.

The analysis of the MAUP presented in this paper has also emphasized how the changes in scale or zoning alter the typology of the response and predictor variables that are eventually provided as the input to a statistical model. Specifically, higher aggregation levels associated with a reduction in overdispersion and kurtosis. This result suggests that the choice of a modelling approach once a change in scale or zoning has been produced should be well addressed, implying the reconsideration or even rejection of a previously selected approach. Furthermore, higher levels of spatial aggregation yielded an overall increase of variance inflation factors, a sign of multicollinearity risk that leads to the conclusion that an excessive aggregation of the data should be avoided, or at least properly checked. The spatial nature of some of the covariates has also provided some clues on their sensitivity towards MAUP. Indeed, covariates having low levels of spatial autocorrelation or generated from point patterns not extremely clustered have displayed a more coherent behaviour among scales and zonings. However, this issue requires a deeper investigation.

Moreover, some limitations of this paper deserve some comment. First, it is worth noting that selecting a proper exposure measure is essential to avoid bias in parameter estimations. Indeed, the lack of consideration of a exposure measure could have a greater impact on statistical estimations than a variation in scale or zoning. Due to the unavailability of traffic volume, we used non-pedestrian road length as a proxy for exposure, but we may be missing some important

Table 9 Models performance in terms of DIC and AIC (for CAR and GWR models, respectively), MAD, SAD and PMAD for all the BSU configurations and levels of aggregation tested.

| AG | BSU | CAR | | | | GWR | | | | | |-------|------|---------|--------|----------|-------|---------|-------|----------|-------|--| | | | DIC | MAD | SAD | PMAD | AIC | MAD | SAD | PMAD | | | AG100 | CT | 1317.37 | 60.55 | 6963.82 | 38.77 | 1348.66 | 54.29 | 6243.25 | 35.21 | | | AG100 | TM | 1311.01 | 62.20 | 7029.14 | 39.13 | 1276.00 | 51.04 | 5767.66 | 32.53 | | | AG100 | TI | 1243.94 | 72.99 | 7737.23 | 43.07 | 1241.80 | 55.19 | 5849.80 | 32.99 | | | AG100 | HEXA | 1306.55 | 102.84 | 11517.79 | 64.95 | 1330.17 | 56.50 | 6327.53 | 35.68 | | | AG200 | CT | 1920.69 | 18.34 | 3667.84 | 20.42 | 2185.88 | 37.25 | 7449.06 | 42.01 | | | AG200 | TM | 1881.92 | 10.35 | 2101.65 | 11.70 | 2156.57 | 37.89 | 7691.15 | 43.37 | | | AG200 | TI | 2116.80 | 42.06 | 8411.20 | 46.83 | 2169.72 | 34.07 | 6813.67 | 38.43 | | | AG200 | HEXA | 1586.44 | 2.56 | 512.40 | 2.89 | 2164.42 | 37.73 | 7546.31 | 42.56 | | | AG300 | CT | 2694.35 | 14.33 | 4299.53 | 23.94 | 3073.02 | 29.43 | 8828.77 | 49.79 | | | AG300 | TM | 2906.43 | 22.33 | 6698.53 | 37.29 | 3017.90 | 31.04 | 9313.18 | 52.52 | | | AG300 | TI | 2739.21 | 11.64 | 3490.79 | 19.43 | 3068.26 | 26.49 | 7947.48 | 44.82 | | | AG300 | HEXA | 2311.54 | 3.57 | 1071.61 | 6.04 | 3050.20 | 27.52 | 8256.20 | 46.56 | | | AG400 | CT | 3455.81 | 11.98 | 4792.86 | 26.68 | 3918.14 | 23.66 | 9463.93 | 53.37 | | | AG400 | TM | 3700.36 | 20.20 | 8078.93 | 44.98 | 3868.05 | 25.05 | 10021.12 | 56.51 | | | AG400 | HEXA | 2981.30 | 2.99 | 1194.31 | 6.74 | 3903.94 | 22.20 | 8879.82 | 50.08 | | | AG500 | CT | 4103.36 | 10.73 | 5367.13 | 29.88 | 4772.90 | 19.95 | 9976.61 | 56.26 | | | AG500 | TM | 4386.43 | 15.87 | 7934.65 | 44.17 | 4724.74 | 21.20 | 10597.96 | 59.77 | | | AG500 | HEXA | 3275.34 | 1.99 | 997.13 | 5.62 | 4768.70 | 17.18 | 8590.50 | 48.45 | |

information. Second, the choice of certain homogeneity criteria to carry out the regionalization process may be another factor, other than scale and zoning, that affects model fitting. In this paper, homogeneity criteria were solely based on crash counts. Other factors, such as land use and socio-economic characteristics, should also be considered in future studies.

To conclude, this paper has provided more evidence regarding the complications that the MAUP can create in the context of a spatial traffic safety analysis. The performance of sensitivity analyses suggested by Xu et al. (2018) considering model estimates for several scales or zonings (or both) seems unavoidable, but this kind of analysis should also include the investigation of the "intermediate" factors that affect statistical inference such as the modelling approach, the multicollinearity shown by the covariates and their spatial autocorrelation. The consideration of all of these factors should help researchers to achieve firmer conclusions, although one cannot forget that it is still likely that the MAUP will never be solved (Manley, 2014).