category: literaturenote citekey: xusensitivityanalysiscontextregional2014 title: "Sensitivity analysis in the context of regional safety modeling: Identifying and assessing the modifiable areal unit problem" authors: "Xu, Pengpeng; Huang, Helai; Dong, Ni; Abdel-Aty, Mohamed" year: 2014 date: 2014-00-00 2014 doi: 10.1016/j.aap.2014.02.012 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2014.02.012" zotero_key: H3UZ9AKB zotero_storage: AVX9KDGC collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Xu, Pengpeng"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

jour nal home page: www.elsevier.com/locate/aap

Sensitivity analysis in the context of regional safety modeling: Identifying and assessing the modifiable areal unit problem

Pengpeng Xua, Helai Huang ^b,∗, Ni Donga, Mohamed Abdel-Aty^c

^a Urban Transport Research Center, School of Traffic and Transportation Engineering, Central South University, Changsha, Hunan 410075, PR China
^b School of Traffic and Transportation Engineering, Central South University, Changsha, Hunan 410075, PR China
^c Department of Civil, Environmental & Construction Engineering, University of Central Florida, Orlando, FL, USA

a r t i c l e i n f o

Article history: Received 26 July 2013 Received in revised form 31 January 2014 Accepted 19 February 2014 Available online 10 April 2014

Keywords: Modifiable areal unit problem Sensitivity analysis Regional safety prediction model REDCAP

a b s t r a c t

A wide array of spatial units has been explored in current regional safety analysis. Since traffic crashes exhibit extreme spatiotemporal heterogeneity which has rarely been a consideration in partitioning these zoning systems, research based on these areal units may be subjected to the modifiable areal unit problem (MAUP).

This study attempted to conduct a sensitivity analysis to quantitatively investigate the MAUP effect in the context of regional safety modeling. The emerging regionalization method-RECDAP (regionalization with dynamically constrained agglomerative clustering and partitioning) was employed to aggregate 738 traffic analysis zones in the county of Hillsborough to 14 zoning schemes at an incremental stepsize of 50 zones based on spatial homogeneity of crash risk. At each level of aggregation, a Bayesian Poisson lognormal model and a Bayesian spatial model were calibrated to explain observed variations in total/severe crash counts given a number of zone-level factors.

Results revealed that as the number of zones increases, the spatial autocorrelation of crash data increases. The Bayesian spatial model outperforms the Bayesian Poisson-lognormal model in accurately accounting for spatial autocorrelation effects, unbiased parameter estimates, and model performance, especially in cases with higher disaggregated levels. Zoning schemes with higher number of zones tend to have increasing number of significant variables, more stable coefficient estimation, smaller standard error, whereas worse model performance. The variables of population density and median household income show consistently significant effects on crash risk and are robust to variation in data aggregation. The MAUP effects may be significantly reduced if we just maintain at about 50% of the original number of zones (350 or larger). The present study highlights MAUP that is generally ignored by transportation safety analysts, and provides insights into the nature of parameter sensitivity to data aggregation in the context of regional safety modeling.

1. Introduction

Regional crash prediction model to explain observed crosssectional variations in crash counts using macro-structural covariates becomes a fairly routine component in traffic safety research. Three essential incentives ensure its rapid development. First, many factors affecting crashes operate at a spatial scale (e.g., trip distribution and generation, land use pattern, and various demographical characteristics). It is desired to relate safety with zonal level factors. Then, there is a need for state agencies to regularly monitor region level safety and provide incentives to reduce the number of traffic casualties in a region's safety program. Therefore, a reliable assessment of safety is indispensable by estimating the aggregating crash potentials associated with the target road network on different spatial scales. Furthermore, road safety is increasingly considered as a necessary component in transportation planning (FHWA, 2005). Regional safety prediction models have been suggested as means of incorporating safety considerations into long term transportation planning (Washington et al., 2006).

To response these needs, the last decade witnessed fast growing scope of scientific research in the context of regional safety analysis. Various spatial units are utilized in previous research. However, most of the employed zoning schemes have specific usage and may not be appropriate for regional safety analysis, as traffic crashes

^∗ Corresponding author. Tel.: +86 731 82656631.

E-mail addresses: huanghelai@csu.edu.cn, huanghelai@hotmail.com (H. Huang), mabdel@mail.ucf.edu (M. Abdel-Aty).

exhibit extreme spatiotemporal heterogeneity which has rarely been a consideration in partitioning zonal systems (Huang et al., 2013). Therefore, research based on these modifiable areal units is probably subject to problems related to the scale and zonal effects.

Up till now, there is a little comprehensive investigation assessing how scale and zonal effects influence statistical results of safety modeling. The assessment of MAUP effects becomes a pressing issue as it may lead to unreliable and inaccurate safety estimation (Wang et al., 2012a; Huang et al., 2013).

The objective of this study is to quantitatively investigate the impacts of specific aggregation configurations on the variations of regional safety modeling. The emerging regionalization method-REDCAP (regionalization with dynamically constrained agglomerative clustering and partitioning, Guo, 2008; Guo and Wang, 2011) is employed to produce a new set of aggregated areal units based on spatial homogeneity of safety related factors. Through this sensitivity analysis, it is possible to explore the scope of statistical results that can be produced for various zoning systems,to determine which safety factors are sensitive to the variation in scale or areal unit definition,to what extent, and to identify which level of data aggregation minimizes the impact of the MAUP. It is expected to highlight the MAUP that is generally ignored by transportation safety analysts, and provide insights into the nature of parameter sensitivity to data aggregation in the context of multivariate safety modeling analysis.

2. Literature review

2.1. Regional crash prediction model

Regional crash prediction models have been widely applied to explain observed variations in crash counts associated with areawide factors. Different area-wide characteristics were considered, including road characteristics such as intersections density (e.g., Huang et al., 2010), road length with different speed limit (e.g., Abdel-Aty et al., 2011; Siddiqui et al., 2012), road length with different functional classifications (e.g., Quddus, 2008; Hadayeghi et al., 2010; Huang et al., 2010), junctions and roundabouts (e.g., Noland and Quddus, 2004; Quddus, 2008); traffic patterns such as traffic flow and traffic speed (Noland and Quddus, 2005; Quddus, 2008; Hadayeghi et al., 2010); trip generation and distribution (Abdel-Aty et al., 2011); environment conditions such as total precipitation/snowfall, and number of rainy/snowy days per year (Aguero-Valverde and Jovanis, 2006); land use (e.g., Hadayeghi et al., 2010; Siddiqui et al., 2012; Pulugurha et al., 2013); and socioeconomic factors such as population density (e.g., Hadayeghi et al., 2003; Huang et al., 2010; Siddiqui et al., 2012), age cohorts (e.g.,Aguero-Valverde and Jovanis, 2006; Quddus, 2008; Hadayeghi et al., 2010; Huang et al., 2010), household incomes (e.g., Huang et al., 2010; Siddiqui et al., 2012) and employment (Quddus, 2008; Hadayeghi et al., 2010; Huang et al., 2010).

As crash data are typically collected with reference to location dimension, two problems arise (LeSage, 1999): (1) spatial correlation exists among the observations, and (2) spatial heterogeneity occurs in the relationships that are modeled. Traditional crash prediction models such as Poisson lognormal model and negative binomial model (i.e., Poisson gamma model) have largely ignored the issue of spatial correlation (spatial dependency) in traffic crash data, which would be misleading as this cannot reflect the true underlying data generating processes. Recently, Bayesian hierarchical models have been widely employed to address the issue of unmeasured spatial correlation in regional safety analysis (Miaou et al., 2003; Aguero-Valverde and Jovanis, 2006; Quddus, 2008; Huang et al., 2010; Siddiqui et al., 2012).

A wide array of spatial units have been explored in previous regional safety studies, such as regions (Washington et al., 1999), counties (Miaou et al., 2003; Aguero-Valverde and Jovanis, 2006; Huang et al., 2010), districts (i.e., a level of sub-national division used for the purposes of local government; Haynes et al., 2007), traffic analysis zones (TAZs, Hadayeghi et al., 2003; Guevara et al., 2004; Hadayeghi et al., 2006, 2010; Abdel-Aty et al., 2011; Siddiqui et al., 2012; Pulugurha et al., 2013), wards (i.e., the primary unit of electoral geography for civil parishes; Noland and Quddus, 2004; Quddus, 2008), enumeration districts (a geographic area assigned to an individual census taker, or enumerator, usually representing a specific portion of a city or county; Noland and Quddus, 2005), and grid-based structure (MacNab, 2004). However, most of above established zoning schemes have specific usage and may not be the optimal zone systems for macro-level safety analysis. As a result, road safety research based on these modifiable areal units is probably subject to modifiable areal unit problem (MAUP).

2.2. Sensitivity analysis for modifiable areal unit problem

MAUP refers to situations that when the boundary of zones used in a spatial analysis changes, the statistical inference and interpretation derived from the zones are also different (Openshaw, 1984). Typically, it can be decomposed into two components: the scale and the zoning effects (Openshaw, 1984; Fotheringham and Wong (1991)). While the scale effect describes the occurrence of statistical result variations using data that are aggregated at different levels, zoning effect refers to the variability introduced by different zoning configurations.

Although the MAUP research has a substantive history,the solutions to MAUP are not straightforward as different geographical phenomena subjects to varied spatial and temporal distributions. Fotheringham et al. (2000) suggested that a sensitivity analysis which examines the parameters estimates using various spatial and scale configurations is desired before presenting analysis results to decision makers. This suggestion of sensitivity analysis is usually implemented by aggregating the smallest available base units into large regions and alternate zones. Based on the sensitivity analysis, the impacts of specific configurations on the variations of statistical analysis results can then be quantitatively investigated.

By considering the MAUP on a line, Thomas (1996) empirically investigated the effect of systematic variation of road segment length on the statistical description of crash counts. Three distinct groups of roadway segments were found: (1) for small segments about 100 m, crash counts followed almost Poisson distribution, (2) for medium size segments (e.g., 300–2000 m), crash counts had an intermediate empirical distribution, and (3) for lager segments (e.g., more than 2000 m), crash counts were almost normally distributed. The empirical work of Thomas (1996) demonstrated that the definition of road segments may also affect the statistical results due to MAUP. However, the study was limited to a univariate exploratory sensitivity analysis and the author did not further investigate how this influence could be inserted into safety modeling.

Recently,Abdel-Aty et al.(2013)investigated the effects of zonal variation (e.g., TAZs, block groups, and census tracts) on the performance of Bayesian Poisson lognormal models with three different dependent variables (e.g., total crashes, severe crash and pedestrian crashes, respectively). Inconsistent parameter estimates were found across the three geographic units, which demonstrated the existence of MAUP effects in regional safety analysis. However, the research was limited to three regions, probably too few to fully assess the impacts of spatial aggregation.

In a more extended research, Fotheringham and Wong (1991) performed linear regression to comprehensively investigate the impact of scale and zoning effects on the stability of parameter estimates. 871 block groups were aggregated to 800, 400, 200, 100, 50 and 25 areal units according to a random procedure with a contiguity constraint. At each level of aggregation, 20 different zoning systems were created. They found that the changes were essentially complex and unpredictable. For example, an increase of 0.1 in the proportion of elderly would create a decrease in the predicted mean family income of only \$308 based on 800 zones. When the data were aggregated to only 25 zones, the results suggested that the same increase would produce a decrease of \$2654. Empirical work of Fotheringham and Wong (1991) provided important theoretical guide on the sensitivity of parameter estimates in multivariate analysis. However, as Openshaw and Rao (1995) claimed, due to the absence of effective zone design tools,the principal value of the above research is restricted to providing rational confirmation of the complexity of MAUP without offering much prospect of a solution.

Unlike conventional statistical perspective, Openshaw (1984) claimed that model building and calibration activities were not independent of the zonal design. By generating a random sample of 261 alternative twenty-two-zone systems and 87 forty-twozone systems from the set of 72 basic spatial units, Openshaw (1977) created optimal zoning systems for a spatial interaction model which ensured a least squares fit with its parameters fixed at arbitrary values. This idea subsequently expanded into an automated zoning procedure (AZP, Openshaw, 1977; Openshaw and Rao, 1995; Martin, 2003). AZP is a heuristic-based optimization approach which starts with an initial zonal system and iteratively refines the solution by reassigning objects to neighboring regions. The advantage ofAZP is thatit allows interaction betweenmodeling exercises and zoning generations since model performance can be selected as one of the optimization goals. However, this methodology is open to criticisms as: researchers are designing spatial units to achieve a satisfactory statistical result (Swift et al., 2008); what constitutes optimal in terms of multivariate analysis is likely to be rather subjective as a zoning system that is optimal for one variable might not be optimal for another (Fotheringham and Wong, 1991); it requires high computational cost, thus inability to process large datasets (Guo, 2008).

Instead of utilizing any aggregation procedure in accordance with above studies, Zhang and Kukadia (2005) employed two distinct scale and zonal schemes (i.e., conventional definition of the census geography and grid-based structures) to diagnose how MAUP influenced the effect of urban form indicators on travel mode choice in Boston. With respect to the grid scheme, estimated coefficients showed consistentincreasing trends as the level of data aggregation increased. In contrast, results with data aggregation at the census geographic unit were quite unstable and perplexing.

Meanwhile, as Anselin (2001) stated, with regards to MAUP, unless there is a homogeneous spatial process underlying the data, any aggregation will tend to be misleading. Yannis et al. (2007) employed a multilevel model in examining the regional effect of enforcement on road accidents in Greece. When using ad hoc geographic clustering scheme, they found that there were some inconsistency in spatial pattern which may be due to the MAUP. Alternatively, they developed a mathematical clustering scheme (k-mean algorithm) based on spatial homogeneity in transport characteristics and road safety parameters, and then a consistent result was found. The results implied thatthe statistical results may be more reliable if the spatial units used are more homogeneous. Nevertheless,the clusters scheme developed by Yannis et al.(2007) derived clusters based on attribute similarity only,then the clusters were divided or merged to form regions in the geographic space. While this type of approach is useful for identifying dependence of multivariate patterns, two drawbacks arise: (1) the clusters in the attribution space are segmented into pieces in the geographical space; (2) the final number of regions is unpredictable and difficult to control.

Regionalization has long been a common approach to aggregate univariate or multivariate data into homogeneous and spatially contiguous regions in geography (Cliff et al., 1975). Traditional regionalization methods place the first priority on attribute similarity within areas, and most are implemented manually or semi-automatically (Guo and Wang, 2011). Advancements in Geographical Information Systems (GIS) have enabled researchers to develop innovative, automated regionalization methods. Among them, REDCAP is potential and has been successfully employed in geography and ecology (Kupfer et al., 2012; Wang et al., 2012b). The method is relatively innovative as compared with the random procedure adopted by Fotheringham and Wong (1991), REDCAP is capable of creating homogeneous regions; with regards to mathematical clustering scheme developed by Yannis et al. (2007), it satisfies the spatial contiguity constraint among derived clusters and the final number of regions can be predetermined; it is efficient to regionalize large datasets that have 10,000 or more spatial objects (Guo, 2008), which outperforms the AZP used by Openshaw (1977); with a regionalization toolkit, it provides an integrated environment for the entire process of data loading, algorithm configuration, result visualization, interactive exploration and result saving. Based on above analysis, REDCAP was adopted to generate homogeneous zoning schemes.

With the aforementioned issues and potential solutions revealed by prior research, this study attempts to conduct a sensitivity analysis to quantitatively investigate the MAUP effect in the context of regional safety modeling. The analysis is based on a 3-year data set from the county of Hillsborough, Florida.

3. Data

3.1. Data collection

Data sets were collected from Hillsborough County, Florida in the year 2005–2007. Hillsborough County is a county located in the west of central Florida with a population of 1,230,000 in 2010. It is the largest county in the Tampa-St. Petersburg-Clearwater Metropolitan Area, with Tampa forming the region's hub.

The smallest available base units were suggested to conduct a sensitivity analysis (Fotheringham et al., 2000). According to the U.S. Census Bureau, census blocks are the smallest geographical unit. However, due to the data unavailability, block groups (BGs) were selected as a substitute in previous research (Fotheringham and Wong, 1991). BGs are collection of census blocks and the lowest level for which census data are made available (Abdel-Aty et al., 2013). Meanwhile, TAZs are the only traffic related zone system and have been widely adopted for crash analysis (Guevara et al., 2004; Hadayeghi et al., 2006, 2010; Abdel-Aty et al., 2011; Siddiqui et al., 2012; Pulugurha et al., 2013). Compared with BGs, TAZs are thought to have better homogeneity as they are special areas delineated by state, MPOs or local transportation officials particularly for tabulating traffic-related data (U.S. Census Bureau, 2001). Hillsborough County contains 796 BGs and 738 TAZs in total. Median areas of BGs and TAZs are 0.366 and 0.370 square mile, respectively, which show that their sizes are quite comparable. Therefore, TAZs were finally adopted as the base units, and the following data sets were aggregated at the TAZ-level.

Crash data were obtained from the Florida Department of Transportation (FDOT) Crash Analysis Reporting System. During the period of 2005–2007, a total of 57,694 crashes occurred, among which4,854 (8.41%) were severe crashes (i.e.,fatal andsevere injury crashes).

Table 1 Summary of variable and descriptive statistics.

Variables	Definition	Mean	S.D.	Min.	Max.
Predictor variable
Total cra	Total number of crashes per TAZ	78.18	72.89	0.00	481
Severe cra	Total number of fatal and severe injury crashes per TAZ	6.58	7.02	0.00	47
Exposure variable
DVMT	Daily vehicle miles traveled (in thousands)	95.07	110.24	0.06	788.77
Explanatory variables
DVMT	Daily vehicle miles traveled (in thousands)	95.07	110.24	0.06	788.77
Inter density	Number of intersections/road length	3.17	5.61	1.00	66.12
Seglen25	Percent of road segment length with 25 mph speed limit	0.72	0.21	0.00	1.00
Seglen45	Percent of road segment length with 45 mph speed limit	0.02	0.05	0.00	0.44
Seglen55 65	Percent of road segment length with 55 65 mph speed limit	0.05	0.10	0.00	0.83
POP density	Population/area (per acre)	3.76	3.38	0.00	19.01
MHINC	Median household income (in thousands)	40.14	20.24	0.00	115.66

Road and traffic-related data were collected mainly from two sources: FDOT's Roadway Characteristics Inventory and GIS maps with Hillsborough. These include total roadway length, road segment length per TAZ with 25/35/45/55/65 mph speed limit, intersections and daily vehicle miles traveled (DVMT).

Moreover, a number of demographic and socioeconomic factors were also investigated, which were downloaded from the U.S. Census. These data include geographical area of each TAZ, population, and incomes.

3.2. Data preparation

To appropriately reflect the effects of various factors on crash occurrence, the raw data were carefully processed and rescaled. In this study, numbers of total crashes and severe crashes were selectedas thedependent variables. DVMT wasutilizedas the exposure variable. After controlling for the exposure variable, a number of factors potentially affecting crash risk could be constructed.

Several road and traffic related variables were coded in the analysis. To avoid the problem of variables redundancy, road segments with posted speed limit were divided into four groups: segments with speed limit 25/35/45/55 65 mph. The proportion of road segments with various speed limits was considered potential factors influencing crash risk. It was found thatthe following variables:the percent of segments with speed limit 25 and the percent of segments with speed limit 35 had some kind of collinearity (i.e., the value of Pearson product-moment correlation is 0.795). Considering the fact that the countermeasures of speed limit 25 mph was widely adopted (the percent of road segmentlength with 25 mph in Hillsborough is 0.76, while the same for 35 mph is 0.14),the percent of segments with speed limit 35 was finally omitted.

Intersections are generally associated withmore traffic conflicts, and as such the intersection density (i.e., the number of intersection/road length) was investigated to understand the variation of crash risk across different geographical units. Even though the amount of travel (DVMT) was utilized in the models as an offset term, DVMT was also included as an explanatory variable to account for possible non-linearity between crash frequency and DVMT (Aguero-Valverde and Jovanis, 2006).

With regards to the socio-demographical variables, population density was the first important factor considered as demonstrated in prior studies (Guevara et al., 2004; Huang et al., 2010; Siddiqui et al., 2012). To reflect area deprivation level, yearly median income per household was employed as a surrogate.

Table 1 provides descriptive statistics of both predictor and explanatory variables considered in the analysis. This includes mean, standard deviation (S.D.), minimum and maximum values of the variables.

4. Methodology

4.1. Regionalization method

REDCAP refers to a family of regionalization methods, termed "regionalization with dynamic constrained agglomerative clustering and partitioning". The goal of REDCAP is to construct a set of homogeneous regions by aggregating contiguous small areas of similar attribute values (e.g., crash risk). To achieve this goal, RED-CAP constructs a cluster hierarchy based on attribute similarities among small areas and then partitions the spatially contiguous cluster tree while explicitly optimizing an overall homogeneity measure (Guo, 2008). The homogeneity measure is the total sum of squared deviations (SSD) (Everitt, 2002), which is defined as:

SSD = $$\sum{r=1}^{k} \sum{i=1}^{nr} \sum{j=1}^{d} (x_{ij} - \bar{x}_j)^2$$

where k is the number of regions, nr is the number of objects in region r, d is the number of variables considered, xij is the value for the jth attribute of the ith object, x¯^j is the mean value of the jth attribute for all objects in region r.

As shown in Fig. 1, REDCAP consists oftwo phases:(1) clustering data with contiguity constraints to produce a spatially contiguous tree; (2) partitioning the tree to generate regions while optimizing the homogeneity measure (i.e., SSD). The color shade of each polygon represents its attribute value and similar color represents similar values. Two polygons are considered contiguous if they share a segment of border. In the first step, as shown in Fig. 1(A), two adjacent and most similar areas are aggregated to form the first cluster; two adjacent and most similar clusters are grouped together to form a high-level cluster; and so on until the whole study area is one cluster. In the second step, as shown in Fig. 1(B), REDCAP partitions the tree by removing edge 11, 15 to create two regions while optimizing the SSD measure. In other word, the two regions are generated in such a way that the total within-region homogeneity is maximized. The partitioning process continues until the desired number of regions is reached.

There are six candidate methods in RECADP which are based on three agglomerative clustering approaches, including the single linkage (SLK), average linkage (ALK), and the complete linkage (CLK), each of which is constrained with spatial contiguity in two different ways (i.e., the first-order constraining and the full-order constraining). As studied by Guo (2008), the Full-Order-CLK outperforms in terms of maximum homogeneity, minimum internal variations, best regional size balance, and retaining original data distribution. Therefore, the Full-Order-CLK is used as the default regionalization method in this research to create homogeneous zoning schemes for this sensitivity analysis.

Fig. 1. An example to illustrate REDCAP: (A) a spatially contiguous tree constructed with a contiguity-hierarchical clustering method; (B) partitioning the tree by removing edge 11, 15 to create two regions that optimize the SSD measure. Source: Guo and Wang, 2011.

Considering the fact that traffic crashes has rarely been a consideration in partitioning the TAZs, crash risk is selected as homogeneity criteria, and the 738 TAZs in the area of Hillsborough are aggregated to 14 zonal schemes from 50 zones to 700 zones. The data for each zonal scheme are derived from the constituent TAZs of the regions, which are automatically generated by REDCAP. It should be noted that in the case of household income, the initial value reported is a median and each aggregated value is formed by calculating the mean of the medians of the constituent TAZs.

4.2. Model configuration

4.2.1. Bayesian Poisson lognormal model

Poisson lognormal model is comparable with negative binomial model while provides more flexibility, thus has been suggested as a substitute for the crash count data (Lord and Mannering, 2010). Poisson lognormal models in a Bayesian framework are developed to reflect the effects of various zonal level factors on crash risk in this study.

Let Yi denotes the number of total/severe crashes in areal unit i, EVi is the exposure variable (i.e., DVMT) of areal unit i, Xi is the covariate vector with parameter coefficient vector, a Poisson lognormal can be specified as follows:

$$Y_i \sim Poisson(\lambda_i)$$

$$\ln(\lambda_i) = \ln(EV_i) + \beta_0 + \mathbf{X_i}\mathbf{\beta} + \theta_i$$

where ⁱ is the parameter of Poisson model, ˇ⁰ is the intercept term, is a vector of coefficient estimates, ⁱ is a random effect to account for unstructured over-dispersion, which is specified via an ordinary, exchangeable normal prior:

$$\theta_i \sim N\left(\frac{0,1}{\tau_h}\right)$$

where ^h is the precision parameter, which follows a prior gamma (0.5, 0.0005).

4.2.2. Bayesian spatial model

Although the Bayesian lognormal model presented above is capable of capturing unstructured heterogeneity, it cannot account for possible spatial dependency of traffic crashes among adjacent zones. For this purpose, Bayesian spatial models were proposed (Besag, 1974) and have been widely applied in regional safety modeling.

The Bayesian spatial model is presented below:

$$Y_i \sim Poisson(\lambda_i)$$

$$\ln(\lambda_i) = \ln(EV_i) + \beta_0 + \mathbf{X_i}\boldsymbol{\beta} + \theta_i + \phi_i$$

where ⁱ is the correlation heterogeneity or spatial correlation, and all other terms are as previously defined.

For the spatial correlation term i, the conditional autoregressive (CAR) prior proposed by Besag (1974) is used:

$$\phi_i \sim N\left(\bar{\phi}_i, \frac{1}{\tau_i}\right) \text{ with } \bar{\phi}i = \frac{\sum{i \neq j} \phij \omega{ij}}{\sum{i \neq j} \omega{ij}} \text{ and } \tau_i = \frac{\tauc}{\sum{i \neq j} \omega_{ij}}$$

in which ωij is binary entries of proximity matrix, and if i and j are adjacent (i.e., a shared border), ωij = 1, otherwise, ωij = 0. ^c is the precision parameter, which is also assumed as a prior gamma (0.5,0.0005).

The proportion of variability in the random effects that is due to spatial autocorrelation is calculated as follows:

$$\alpha = \frac{sd(\phi)}{sd(\phi) + sd(\theta)}$$

in which, sd is the empirical marginal standard deviation function. Finally, four Bayesian models are developed for each of the 15

zonal schemes (i.e., 14 zone partitioning schemes and 1 original base unit), and as such a total of 60 models (i.e., 4 × 15) are estimated. The elaborated model specification is illustrated in Table 2.

4.2.3. Model assessment

In response to compare above model performance, Deviance Information Criterion (DIC) is calculated. DIC is based on the posterior distribution of the deviance statistic and can be seen as a

Table 2 Model specifications.

Model	1T-BP	2T-BS	3S-BP	4S-BS
Yi	Total cra	Total cra	Severe cra	Severe cra
EVi	DVMT	DVMT	DVMT	DVMT
Xi	DVMT	DVMT	DVMT	DVMT
	Inter density	Inter density	Inter density	Inter density
	Seglen25	Seglen25	Seglen25	Seglen25
	Seglen45	Seglen45	Seglen45	Seglen45
	Seglen55 65	Seglen55 65	Seglen55 65	Seglen55 65
	POP density	POP density	POP density	POP density
	MHINC	MHINC	MHINC	MHINC

Note: 1T-BP: Bayesian Poisson lognormal model with total crash frequency as dependent variable; 2T-BS: Bayesian spatial model with total crash frequency as dependent variable; 3S-BP: Bayesian Poisson lognormal model with severe crash frequency as dependent variable; 4S-BS: Bayesian spatial model with severe crash frequency as dependent variable; and all other abbreviations are defined as in Table 1.

generalization of the Akaike Information Criterion. Specifically, the DIC is defined as:

$$DIC = D(\bar{\theta}) + 2p_D = \bar{D} + p_D$$

in which D( ¯) is the deviance evaluated at ¯, the posterior means of the parameters of interest, pD is the effective number of parameters in the model, and D¯ is the posterior mean of the deviance statistic D(). Models with lower DIC are preferred.

Another two common measures, Mean Absolute Deviance (MAD) and Mean Squared Predictive Error (MSPE) are also be used to test goodness of fit of models. The MAD and MSPE are expressed as:

$$\mathsf{MAD} = \frac{1}{n} \sum_{i=1}^{n} \left| y_i^{pred} - y_i^{obs} \right|$$

$$MSPE = \frac{1}{n} \sum_{i=1}^{n} (y_i^{pred} - y_i^{obs})^2$$

where yobs ⁱ is the observed number of total/severe crashes in areal unit ⁱ, and ^ypred ⁱ is the expected total/severe crash number obtained by above regional crash prediction models. Models associated with lower MAD and MSPE fit better to the data.

4.3. Standard-difference-in-means test

In order to shed light on the sensitivity of the parameter estimates to variations in data aggregation, it would be useful to determine whether the differences in parameter estimates are statistically significant. The standard-difference-in-means test (Fotheringham and Wong, 1991), shown in the following equation, provides a tool to test whether the coefficients obtained from one geographical aggregation are different in statistical terms:

$$t = \frac{\hat{\beta}{ij} - \hat{\beta}{ik}}{\text{SE} \left| \hat{\beta}{ij} - \hat{\beta}{ik} \right|}$$

where ˇˆij and ˇˆik are estimated coefficients ofindependent variable i at the jth and kth aggregated level, respectively. SE is the standard error.

5. Results and discussion

Based on above data and methodology availability, we present the detailed results for Moran's I univariate evaluation and parameter estimates multivariate evaluation in the following section.

5.1. Moran's I univariate evaluation

Previous studies (Aguero-Valverde and Jovanis, 2006; Quddus, 2008; Huang et al., 2010; Siddiqui et al., 2012) revealed that spatial correlation of traffic crashes existed widely across adjacent spatial zones. Traditional crash prediction models such as Bayesian Poisson lognormal models have largely ignored the issue of spatial correlation in crash data. If residuals in a Bayesian Poisson lognormal are spatially autocorrelated across the study area, the underlying model assumptions are violated and biased estimates may be produced (Anselin, 2001). To explore the slope of spatial autocorrelation to variations in data aggregation, Moran's I evaluation is conducted for each of the 15 zoning schemes based on eight variables, including number of total/severe crashes, total/severe crash risk, and residuals for four Bayesian models. The variations of Z-scores of Moran's I are plotted in Fig. 2.

As shown in Fig. 2(A), when the number of zones is small enough (i.e., 50/100 zones), there appeared no strong spatial correlation among crash data. As the total number of zones increases, the spatial autocorrelation of traffic crashes is statistically significant and increases in a general trend. This may be due to that at higher aggregation, objects with similar values of crash risk are gradually aggregated together, and the spatial units across the study area become less homogeneous, as a result, the spatial correlation of the variables would be less statistically significant. Moreover, it can also justify the need to incorporate a spatial error term into a crash model, especially when the number of zones reaches 150 or larger.

In order to check if after controlling the independent variables, there is still any spatial correlation left, a Moran's I evaluation is conducted for model residuals ofthe four Bayesian models. As illustrated in Fig. 2(B), with regards to the Bayesian Poisson lognormal models, when the number of zones is less than 150, the spatial autocorrelation among the residuals are not significant at a 95% confidence level (i.e., the Z-scores are less than 1.96). When the number of zones gradually increases, the spatial autocorrelation of the residuals is statistically significant and shows a tendency to increase. In contrast, the residuals of Bayesian spatial models demonstrate no significantly spatial autocorrelation regardless of the scale definitions. The results are unsurprising as the Bayesian Poisson lognormal model is essentially a non-spatial model which accounts for unobserved heterogeneity only, while the Bayesian spatial model is capable of taking into account unobserved heterogeneity and spatial autocorrelation effects simultaneously. This finding demonstrates the Bayesian spatial model outperforms in accurately accounting for spatial autocorrelation effects and unbiased parameter estimates, supporting its widely use in existing safety analysis (Miaou et al., 2003; Aguero-Valverde and Jovanis, 2006; Quddus, 2008; Huang et al., 2010; Siddiqui et al., 2012).

Fig. 2. Z-scores of Moran's I for: (A) the dependent variables: total/severe crashes, and total/severe crash risk; (B) residuals of four Bayesian models: 1T-BP; 2T-BS; 3S-BP; 4S-BS.

5.2. Parameter estimates multivariate evaluation

The model estimation results are summarized in Tables 3 and 4. The analyses place emphasis on the sensitivity of parameter estimates, the effects of spatial correlation, and model fitting to variation in data aggregation.

If the MUAP did not exist, there would be relatively little variation in modeling results and certainly no systemic variations. From Table 3, results reveal obvious variations in parameter estimation under different zoning schemes, which indicates the extensive presence of MAUP in crash analysis. For example, in the case of DVMT estimates in Model 2T-BS, when the model is calibrated on 738 TAZs, the coefficient is −0.34, which is statistically different from zero at the 95% confidence interval. When the data are aggregated to only 50 zones, the mean estimate is only −0.13 and not statistically significant. This fact sets alarm on previous regional crash prediction models without taking the MAUP issues into account, from which biased model inferences may have been induced latently.

To determine whether the variations depicted above were statistically significant, the standard difference-in-means tests were performed. Considering the fact that the calculated results are numerous and the variation tendencies of regression parameters between total crashes models and severe crashes models are similar, we only present the test results for total crashes models (i.e., Model 1T-BP and Model 2T-BS). The specific results are illustrated in Fig. 3, which depict the mean of the parameter estimates with 1.5 standard errors at either side. Apparently, if the lower 1.5 standard error symbol of one parameter estimate is above the upper 1.5 standard error symbol of another, the two are very likely to be significantly different at the 95% confidence.

As shown in Fig. 3, several general observations are worth noting. First, with the increase of number of zones, all parameter estimates except for DVMT tend to be stable and the disparity of parameter values between different zoning schemes is slight, especially when the number of zones reaches 350 or lager. This result indicates that the parameter estimates are more sensitive to variation when the number of zones is small, implying the MAUP effect is more serious with higher aggregation level. Second, all the variables exhibit the same trend that as the data become increasingly aggregated, the standard errors of the parameter estimates increase, which is expected as the standard error depends

Table 3 Models results with total crashes frequency as the dependent variable.

Zones	50	100	150	200	250	300	350	400	450	500	550	600	650	700	738
Model 1T-BP
Intercept	−6.81	−6.88	−6.96	−6.96	−6.91	−6.93	−6.94	−6.94	−6.98	−6.97	−6.96	−6.95	−6.94	−6.95	−6.95
DVMT	−0.13	−0.08	−0.10	−0.12	−0.12	−0.13	−0.19	-0.19	−0.21	−0.22	−0.23	−0.25	−0.28	−0.32	−0.34
Inter density	0.47	0.31	0.26	0.18	0.18	0.15	0.14	0.14	0.13	0.13	0.13	0.13	0.13	0.12	0.13
Seglen25	0.10	0.08	0.07	0.08	0.09	0.10	0.10	0.08	0.09	0.08	0.07	0.06	0.05	0.03	0.04
Seglen45	0.02	0.06	0.14	0.11	0.12	0.06	0.07	0.03	0.04	0.04	0.03	0.04	0.05	0.07	0.08
Seglen55 65	−0.29	−0.28	−0.24	−0.19	−0.17	−0.11	−0.07	−0.07	−0.07	−0.07	−0.08	−0.08	−0.08	−0.05	−0.03
POP density	0.62	0.52	0.43	0.40	0.41	0.37	0.35	0.32	0.30	0.29	0.28	0.28	0.27	0.28	0.26
MHINC	−0.35	−0.26	−0.25	−0.29	−0.24	−0.20	−0.22	−0.21	−0.25	−0.24	−0.24	−0.22	−0.20	−0.19	−0.19
sd()	1.42	1.40	1.23	1.14	1.04	1.00	0.94	0.93	0.89	0.86	0.83	0.81	0.78	0.75	0.73
DIC	456	864	1282	1676	2076	2473	2860	3244	3601	3989	4356	4724	5087	5432	5705
MAD	2.06	2.38	2.50	2.67	2.91	3.03	3.22	3.29	3.44	3.59	3.71	3.86	3.98	4.18	4.30
MSPE	7.07	8.24	9.99	11.73	14.26	15.57	17.59	18.40	19.84	21.52	23.09	24.85	26.73	29.08	30.60
Model 2T-BS
Intercept	−6.82	−6.89	−6.96	−6.96	−6.92	−6.93	−6.94	−6.94	−6.98	−6.97	−6.96	−6.95	−6.94	−6.95	−6.95
DVMT	−0.10	−0.09	−0.12	−0.12	−0.13	−0.13	−0.18	−0.19	−0.21	−0.22	−0.22	−0.24	−0.26	−0.32	−0.34
Inter density	0.22	0.16	0.20	0.11	0.09	0.08	0.07	0.07	0.06	0.06	0.06	0.06	0.07	0.07	0.09
Seglen25	0.09	0.10	0.08	0.09	0.11	0.12	0.12	0.11	0.10	0.09	0.07	0.06	0.05	0.03	0.05
Seglen45	0.09	0.03	0.12	0.10	0.10	0.05	0.06	0.03	0.03	0.03	0.03	0.04	0.05	0.08	0.09
Seglen55 65	−0.42	−0.40	−0.26	−0.21	−0.19	−0.13	−0.11	−0.10	−0.10	−0.10	−0.11	−0.10	−0.11	−0.07	−0.05
POP density	0.48	0.41	0.39	0.37	0.37	0.31	0.29	0.26	0.24	0.23	0.22	0.21	0.21	0.21	0.19
MHINC	−0.48	−0.41	−0.29	−0.32	−0.28	−0.24	−0.29	−0.21	−0.25	−0.22	−0.20	−0.16	−0.14	−0.13	−0.12
˛	0.30	0.32	0.33	0.30	0.31	0.30	0.30	0.32	0.31	0.32	0.34	0.35	0.35	0.37	0.39
sd()	0.56	0.58	0.54	0.45	0.44	0.40	0.38	0.40	0.37	0.37	0.38	0.39	0.38	0.38	0.40
sd()	1.29	1.25	1.08	1.04	0.94	0.91	0.87	0.84	0.81	0.78	0.75	0.72	0.69	0.65	0.63
DIC	457	862	1282	1675	2075	2472	2859	3242	3599	3987	4353	4721	5084	5427	5700
MAD	1.88	2.16	2.43	2.63	2.84	2.99	3.16	3.20	3.34	3.45	3.54	3.68	3.79	3.98	4.07
MSPE	6.67	7.45	9.75	11.64	13.98	15.26	17.24	17.74	19.24	20.74	21.88	23.36	25.05	26.96	28.07

Note: The titled bold numbers mean statistical significance at 90% significance level; while the bold ones mean statistical significance at 95% significance level.

Table 4 results withsevere crashes frequency as the dependent variable.

Mo del 3S -BP Int −9 .55 −9 .54 −9 .54 −9 .56 −9 .51 −9 .52 −9 .54 −9 .53 −9 .57 −9 .56 −9 .56 −9 .55 −9 .54 −9 ept erc DV MT −0 .10 −0 .03 −0 .08 −0 .09 −0 .09 −0 .10 −0 .14 −0 .15 −0 .16 −0 .16 −0 .17 −0 .19 −0 .21 −0 Int d ity −0 .13 −0 .11 −0 .08 −0 .18 −0 .18 −0 .21 −0 .19 −0 .18 −0 .20 −0 .19 −0 .21 −0 .20 −0 .19 −0 er ens Seg len 25 0.1 4 0.1 0 0.0 9 0.0 8 0.1 5 0.1 5 0.1 4 0.1 3 0.1 3 0.1 3 0.1 1 0.1 1 0.0 8 0.0 len Seg 45 −0 .08 −0 .09 0.0 0 −0 .04 −0 .04 −0 .06 −0 .04 −0 .06 −0 .03 −0 .03 −0 .04 −0 .01 −0 .01 0.0 Seg len 55 65 −0 .21 −0 .19 −0 .15 −0 .11 −0 .10 −0 .09 −0 .07 −0 .07 −0 .08 −0 .07 −0 .09 −0 .09 −0 .09 −0 PO d P ity 0.4 8 0.4 2 0.3 5 0.3 1 0.3 2 0.2 5 0.2 2 0.2 1 0.1 9 0.1 7 0.1 7 0.1 6 0.1 5 0.1 ens MH INC −0 .51 −0 .47 −0 .35 −0 .35 −0 .30 −0 .25 −0 .25 −0 .22 −0 .24 −0 .22 −0 .22 −0 .20 −0 .18 −0 sd( ) 1.2 8 6 1.0 3 0.9 6 0.9 2 0.8 8 0.8 6 0.8 6 0.8 6 0.8 3 0.8 0.7 9 0.7 1.1 1.1 1 1 DIC 330 588 879 112 7 138 4 164 3 188 6 214 0 233 7 257 8 279 5 30 10 322 6 343 1 MA D 2.2 2 2.5 6 2.7 4 2.8 9 3.1 0 3.2 4 3.3 2 3.3 7 3.3 7 3.4 0 3.4 4 3.5 0 3.5 0 3.6 MS PE 7.8 7 9.5 9 11. 60 13. 09 14. 98 16. 64 17. 03 17. 70 17. 69 18. 01 18. 58 19. 30 19. 74 21. Mo del 4S -BS Int −9 .55 −9 .58 −9 .58 −9 .59 −9 .55 −9 .56 −9 .57 −9 .55 −9 .59 −9 .58 −9 .58 −9 .57 −9 .59 −9 ept erc DV MT −0 .10 −0 .03 −0 .08 −0 .10 −0 .10 −0 .10 −0 .15 −0 .16 −0 .17 −0 .18 −0 .18 −0 .21 −0 .24 −0 Int d ity −0 .04 −0 .03 −0 .02 −0 .09 −0 .07 −0 .03 −0 .04 −0 .03 −0 .04 −0 .02 −0 .04 −0 .02 −0 .02 −0 er ens Seg len 25 0.1 8 0.1 5 0.0 9 0.1 2 0.1 9 0.2 0 0.1 8 0.1 6 0.1 5 0.1 5 0.1 2 0.1 2 0.0 9 0.0 len Seg 45 −0 .14 −0 .11 −0 .06 −0 .08 −0 .04 −0 .03 −0 .01 −0 .03 0.0 0 0.0 1 0.0 0 0.0 3 0.0 4 0.0 Seg len 65 −0 .26 −0 .23 −0 .14 −0 .07 −0 .06 −0 .06 −0 .05 −0 .06 −0 .07 −0 .06 −0 .08 −0 .06 −0 .07 −0 55 PO d P ity 0.5 3 0.4 9 0.4 0 0.3 3 0.3 4 0.2 8 0.2 6 0.2 2 0.2 0 0.1 7 0.1 7 0.1 7 0.1 6 0.1 ens MH INC −0 .54 −0 .48 −0 .35 −0 .39 −0 .35 −0 .31 −0 .27 −0 .26 −0 .27 −0 .24 −0 .22 −0 .17 −0 .15 −0 0.4 6 0.4 8 0.4 9 0.5 0 0.5 1 0.5 1 0.4 9 0.4 8 0.4 8 0.4 8 0.4 9 0.4 8 0.4 7 0.5 ˛ sd( ) 0.9 6 0.8 4 0.8 0 0.7 4 0.7 2 0.6 8 0.6 6 0.6 1 0.6 0 0.5 9 0.6 0 0.5 7 0.5 4 0.5 sd( ) 1.1 2 0.9 2 0.8 2 0.7 3 0.6 8 0.6 6 0.6 8 0.6 6 0.6 5 0.6 4 0.6 2 0.6 1 0.6 1 0.5 DIC 330 587 878 112 5 137 9 163 3 187 6 212 6 232 3 256 2 277 4 298 5 319 9 340 0 MA D 2.0 4 2.4 5 2.6 0 2.7 7 2.9 6 3.0 7 3.2 0 3.2 2 3.2 2 3.2 6 3.2 9 3.3 6 3.3 8 3.4	738	700

	−9 .51	.52
	−0 .28	.27
	−0 .17	.19
	0.0 6	5
	0.0 1	2
	−0 .05	.07
	0.1 6	7
	−0 .16	.17
	0.7 5	6
	359 7
	3.7 4	6
	22. 62	40

	−9 .53	.54
	−0 .28	.26
	−0 .05	.04
	0.0 7	6
	0.0 6	5
	−0 .01	.03
	0.1 4	6
	−0 .11	.14
	0.5 4	0
	0.6 0	6
	0.5 1	6
	355 7
	3.5 0	5
MS PE 7.4 6 9.1 3 10. 30 11. 75 13. 33 14. 35 15. 47 15. 81 15. 79 16. 22 16. 46 17. 19 17. 72 18.	19. 07	48

Note: The titled bold numbers meanstatistical significance at 90% significance level; while the bold ones mean statistical significance at 95% significance level.

Fig. 3. Mean and 1.5 standard error intervals of regression parameter estimates for model 1T-BP and model 2T-BS.

in part on the number of observations used in the model calibration. Third, with regards to the statistical significance, more variables are significant with the increased number of zones, indicating that the safety interpretation power by those factors also increases. It may be due to the fact that the aggregate models have relatively smaller sample sizes, resulting in larger standard errors, thus smaller number of significant variables. Similar findings could also be found in previous research (Fotheringham and Wong, 1991; Wang et al., 2012a). Forth, concerning the magnitudes of parameter estimates, the coefficients associated with the variables of Inter density, Seglen55 65, POP density, and MHINC systematically decrease as the number of zones increases, whereas the DVMT coefficients display a reverse trend. For the other two sets of parameters (i.e., Seglen25 and Seglen45), there is relatively little systematic variation. Fifth, it can also be observed that: only the parameter estimates of DVMT at 700/738-unit level are significantly different from the ones at 300-unit or less, implying that

compared with other parameter estimates, DVMT is relatively more sensitive to zonal definition when data were measured at higher aggregation; the variables of POP density and MHINC have consistent significant effects on crash risk and are robust to changes in data aggregation.

Based onthe above analysis, it canbe concluded that wemay significantly reduce the MAUP effects for all the explanatory variables employed in this study if we just maintain the number of zones at 350 or larger (i.e., about 50% of the original number of zones).

In the model specification, we used multiple zone-wide explanatory variables, i.e., DVMT, intersection density, percent of road length with different speed limits, population density and median household income.

DVMT is found to be statistically significant with negative effect on total/severe crash risk, especially in the models with higher number of zones (i.e., 350 or larger). The results suggest that the total/severe crash risk decrease with DVMT. Similar results were also reported in the study of Aguero-Valverde and Jovanis (2006).

Intersections are well-known as hot spots in the road network. Given equal road length, more intersections were revealed to be associated withhigher crashfrequency inprevious research(Huang et al., 2010; Abdel-Aty et al., 2011; Siddiqui et al., 2012). In our study, intersection density is observed to be significant in most Bayesian Poisson lognormal models, whereas statistically insignificant in most Bayesian spatial models. As mentioned above, the residuals in the Bayesian lognormal models were found to be spatially correlated across the study area, which may result in biased estimates of the parameter. This may be the reason why Bayesian Poisson lognormal models show inconsistent results in the significance of intersection density.

Regarding the percent of roadway with different speed limits, two variables are found significant for the total crashes models. The percent of road length with speed limit 45 mph was found positively associated at the 700/738-unit level. While road length with higher posted speed limit(i.e., 55 and 65 mph) have a negative effect when the number of zones reaches 100 or larger in Model 2T-BS. For the severe crash models, the percent of road length with speed limit 25 mph have a significantly positive relationship when the number of zones ranges from 250 to 650 in Model 4S-BS. The results indicate that roadways with higher speed limits have lower crash risk. The association between speed and crash risk has been well recognized in previous road safety analysis (Aarts and Schagen, 2006).

With respect to the socio-economic factors, population density is positively related to the risks of total/severe crashes, suggesting more residents in an area have more activities that could result in more traffic crashes. Similar findings were also reported in prior studies (Guevara et al., 2004; Siddiqui et al., 2012), though Noland and Quddus (2004) reported an opposite result: the lower population density experiences relatively more casualties in England. In addition, median household income was found to have negative effects on total/severe risk, implying that regions with higher median household income are relatively safer. It may be due to that individuals in relatively affluent areas are likely more educated and seek to avoid risk activities and generally own cars with better safety performance (Huang et al., 2010).

Concerning the effect of spatial autocorrelation, the variations contributed by spatial correlations in all models are always statistically significant, indicating that the effect of spatial clustering is substantial in the crash data. It is also interesting to find the proportions of spatial effects do not exhibit the same tendency as the dependent variables illustrated in Fig. 2(A). This may be due to that the spatial correlation of total/severe crashes is explained by the spatial correlation of explanatory variables specified in the models with constant rate, regardless of the scale definitions for data aggregation.

In order to compare model performance, DIC, MAD and MSPE were calculated. As shown in Tables 3 and 4, it can be seen that: (1) for the models with the same response variable at identical aggregate level, the Bayesian spatial models outperformed Bayesian Poisson lognormal models with lower DIC, MAD and MSPE; (2) the DIC values increase with the increase of number of zones. Given the theoretical nature of DIC, the result is not surprising as the increase of number of zones (i.e.,the increase of sample sizes) leads to significantincrease ofthe effective number of parameters, which is an importation component of DIC; (3) zoning schemes with smaller number of zones have lower MAD and MSPE, indicating that the model performance increases with respect to aggregation. This finding is similar to previous studies, such as Fotheringham and Wong (1991), Zhang and Kukadia (2005), and Wang et al. (2012a).

6. Conclusions

Toward a proactive approach, regional safety prediction models are increasingly considered as a vital component in transportation safety planning (TSP) to forecast crash potentials (measures) for alternative transportation planning schemes. Once aggregation is conducted,the issue of MAUP must be explicitly evaluated.And this sensitivity analysis would be an effective approach for analysts to quantitatively measure the scope and magnitude of the MUAP by reporting results at different levels of data aggregation and with different areal unit definitions.

Based on a 3-year data set from the county of Hillsborough, Florida, this study conducts a sensitivity analysis to quantitatively investigate the MAUP effect in the context of regional safety modeling. Results revealed that: as the number of zones increases, the spatial autocorrelation of crash data increase in a general trend; the Bayesian spatial model outperforms the Bayesian lognormal model in accurately accounting for spatial autocorrelation effects, unbiased parameter estimates, and model performance, especially in the case with higher disaggregate level; zoning schemes with higher number of zones (i.e., 350 or larger) tend to have increasing number of significant variables, more stable coefficient estimation, smaller standard error, whereas worse model performance; the variables of population density and median household income show consistent significant effects on crash risk and are robust to variation in data aggregation; we may significantly reduce the MAUP effects in this study if we just maintain at about 50% of the original number of zones (350 or larger).

Admittedly, TAZs are now the only traffic related zone system and superior in easily being integrated with the transportation planning process, thus having been widely utilized in regional safety analysis. However, considering the fact that TAZs are not specially delineated for traffic crash analysis but for long range transportation plans, it might not be the optimal zone systems for regional safety analysis. Therefore, future studies are desired to explore some adjustments to current TAZs structure from the perspective of safety analysis.

Acknowledgements

This work was jointly supported by:(1) Natural Science Foundation of China (Nos. 71371192, 51108465); (2) FHWA Broad Agency Announcement DTFH61-10-R-00013, Transformational Changes and Revolutionary Advances for Transportation Planning: Developing Advanced Spatial Analysis Tools in Transportation Safety Planning: Examining Planning Zones and Trip Data; (3) Technology Project of Ministry of Transport of the P.R. China (No. 20113187851460).