category: literaturenote citekey: elbasyounymethodaccountoutliersdevelopment2010 title: A method to account for outliers in the development of safety performance functions authors: "El-Basyouny, Karim; Sayed, Tarek" year: 2010 date: 2010-00-00 2010 doi: 10.1016/j.aap.2010.02.003 publication: "Accident Analysis & Prevention" url: "https://doi.org/10.1016/j.aap.2010.02.003" zotero_key: 8EM87T8K zotero_storage: HXCQI8F5 collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "El-Basyouny, Karim"

status: converted

Contents lists available at ScienceDirect

Accident Analysis and Prevention

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

A method to account for outliers in the development of safety performance functions

Karim El-Basyouny∗, Tarek Sayed¹

Dept. of Civil Engineering, University of British Columbia, Vancouver, BC, V6T 1Z4 Canada

article info

Article history: Received 19 October 2009 Received in revised form 3 February 2010 Accepted 5 February 2010

Keywords: Contaminated normal models Discrete mixtures Full Bayes estimation Markov Chain Monte Carlo Multivariate Poisson lognormal regression Outlier detection Scale-mixture models Student t models

abstract

Accident data sets can include some unusual data points that are not typical of the rest of the data. The presence of these data points (usually termed outliers) can have a significant impact on the estimates of the parameters of safety performance functions (SPFs). Few studies have considered outliers analysis in the development of SPFs. In these studies, the practice has been to identify and then exclude outliers from further analysis. This paper introduces alternative mixture models based on the multivariate Poisson lognormal (MVPLN) regression. The proposed approach presents outlier resistance modeling techniques that provide robust safety inferences by down-weighting the outlying observations rather than rejecting them. The first proposed model is a scale-mixture model that is obtained by replacing the normal distribution in the Poisson-lognormal hierarchy by the Student t distribution, which has heavier tails. The second model is a two-component mixture (contaminated normal model) where it is assumed that most of the observations come from a basic distribution, whereas the remaining few outliers arise from an alternative distribution that has a larger variance. The results indicate that the estimates of the extra-Poisson variation parameters were considerably smaller under the mixture models leading to higher precision. Also, both mixture models have identified the same set of outliers. In terms of goodness-offit, both mixture models have outperformed the MVPLN. The outlier rejecting MVPLN model provided a superior fit in terms of a much smaller DIC and standard deviations for the parameter estimates. However, this approach tends to underestimate uncertainty by producing too small standard deviations for the parameter estimates, which may lead to incorrect conclusions. It is recommended that the proposed outlier resistance modeling techniques be used unless the exclusion of the outlying observations can be justified because of data related reasons (e.g., data collection errors).

1. Introduction

The application of safety performance functions (SPFs) in assessing the safety of a road entity (intersection, road segment, etc.) has become a standard practice among safety researchers and practitioners (Kulmala, 1995; Miaou and Lord, 2003; Sawalha and Sayed, 2001, 2006).

Recently, several advances have been made in the development of SPFs. These include: (i) the application of themultivariate Poisson lognormal (MVPLN) to model accident count data at different levels of severity (Tunaru, 2002; Park and Lord, 2007; Ma et al., 2008; Ye et al., 2009; Aguero-Valverde and Jovanis, 2009; El-Basyouny and Sayed, 2009a); (ii) the use of random parameters (Anastasopoulos and Mannering, 2009; El-Basyouny and Sayed, 2009b); (iii) the use of finite mixtures of Poisson or negative binomial (NB) regression models for capturing heterogeneity in crash count models (Park and Lord, 2009); (iv) the application of different modeling techniques to overcome the excessive zeroes observed in accident data (Shankar et al., 1997; Kumara and Chin, 2003; Qin et al., 2004; Lord et al., 2005), see also Malyshkina et al. (2009) for the use of a two-state Markov Switching model to analyze accident frequencies; and (v) the use of variable dispersion parameters to model accident count data (Heydecker and Wu, 2001; Miaou and Lord, 2003; Miranda-Moreno et al., 2005; El-Basyouny and Sayed, 2006; Mitra and Washington, 2007; Lord and Park, 2008).

Traditionally, SPFs are developed using the Poisson-gamma hierarchy, which leads to the NB regression model, or the Poissonlognormal (PLN) hierarchy. However, the successful application of these hierarchical models requires careful examination of various related aspects concerning the random effects (are the data overdispersed?), exchangeability (are there outlying observations?), model building (which covariates to include?), the choice of hyper

^∗ Corresponding author at: Dept. of Civil Engineering, University of British Columbia, 2002 - 6250 Applied Science Lane, Vancouver, BC, V6T 1Z4 Canada. Tel.: +1 604 716 4470.

E-mail addresses: basyouny@civil.ubc.ca (K. El-Basyouny), tsayed@civil.ubc.ca (T. Sayed).

¹ Tel.: +1 604 822 4379.

parameters, and the parametric form of the hierarchical model. All of these issues have been extensively researched in the traffic safety literature, except the treatment of outliers, which has received limited attention.

Accident data sets can include some unusual data points (observations) that are not typical of the rest of the data. The presence of these data points (usually termed outliers) can have significant impact on the estimates of the parameters of SPFs. Few studies have discussed the effect of outliers and how to handle them in developing SPFs. Generally, it is recommended that outliers be identified, according to some statistical measures, and then removed from the data set, e.g., seeKulmala (1995) who used the leverage statistic and Sawalha and Sayed (2006) who used Cook's distance for outliers' detection.

Typically, there are three approaches for analyzing data with potential outliers: (i) the "outliers ignoring" approach, (ii) the "outliers rejecting" approach, and (iii) the "outliers accommodating" approach.

Proponents of the outliers ignoring approach argue that (i) the distinction between "true" outliers (to be deleted) and "false" outliers (to be kept) is usually blurred in practice, (ii) every time an outlier is removed and the model is re-run new outliers are often discovered, leading to a vicious loop, (iii) if few outliers were identified that represent a small percentage of the sample size (e.g., less than 5%) it is still acceptable to include them, especially if the analyst is not certain about whether or not they are outliers. However, it should be noted that while including some non-influential outliers maybe tolerable, influential outliers would significantly impact the model fit leading to faulty inference (Sawalha and Sayed, 2006).

On the other hand, the main problem with the popular outliers rejecting approach is its failure to reflect the uncertainty in the exclusion process, which results in too small standard errors that might lead to faulty inferences (Lange et al., 1989).

The third approach "outlier accommodating", which is being advocated in this paper, uses robust methods to down-weight potential outliers. The Robust models have been proposed in the statistical literature to down-weight the outlying observations rather than rejecting them (Little, 1988; Lange et al., 1989; Geweke, 1993; Albert, 1999). In particular, a Student t distribution is proposed as an alternative to the normal in the Poisson-lognormal hierarchy leading to a scale-mixture model. Another variant is the two-component mixture of gamma or normal priors (depending on the hierarchy), where it is assumed that most of the observations come from a basic distribution, whereas the remaining few outliers arise from an alternative distribution that has a larger variance. Relles and Rogers (1977) compared the performance of the outlier rejecting approach with those of several robust estimators of location. They found that the outlier rejecting approach was at least 20% less efficient than the best estimate in all situations studied. Consequently, they recommended that popular robust estimators should be adopted.

Thus, it is recommended to fit the MVPLN (or PLN) model first and examine the residuals visually for potential outliers. A thorough investigation of such outliers should then be conducted. If there were errors due to data collection and recording, which cannot be corrected, or if the respective sites are genuinely different from the rest of the sites in the database, then these outliers should be excluded from further analysis. Otherwise, in cases where the status of the outliers cannot be determined, the proposed mixture methods can be used to minimize the potential biases caused by these outliers.

This paper introduces two outlier resistance mixture models that are based on the MVPLN approach. It is proposed that these mixture models are considered while developing SPFs as they provide robust inference. The proposed approach is applied to a sample of accident and traffic volume data corresponding to 99 signalized intersections in the city of Edmonton, Alberta, where the accidents were classified into property damage only (PDO) and injuries and fatalities (I + F). The objective is to compare the proposed models in terms of their goodness-of-fit and inference. The statistical software WinBUGS (Lunn et al., 2000), which provides MCMC sampling for a variety of Bayes models, was used in the analysis.

2. Methodology

The three approaches discussed above "outlier ignoring", "outlier rejection", and "outlier accommodating" will be compared. The multivariate Poisson lognormal (MVPLN) will be used to construct models incorporating all data points (outliers ignoring approach). For the outliers' rejection approach, the first step is to detect the outliers and then exclude them from further analysis. The outlier detection is accomplished by using the Mahalanobis distance as discussed below. The outliers accommodating approach utilizes outlier resistance modeling techniques that provide robust safety inferences by down-weighting the outlying observations rather than rejecting them. Two outlier resistance techniques are presented, namely: scale-mixture and two-component models. Prior specification, parameter estimation techniques, and model comparisons within a full Bayes context are also discussed.

2.1. The multivariate Poisson lognormal (MVPLN) model

The MVPLN model presented here is the same as that developed previously by El-Basyouny and Sayed (2009a). Suppose that data are available on road accidents at n sites, where the accidents at each site are classified into K categories. Let yik denote the number of accidents at location i in category k and define yi ⁼ (yi1, yi2, ..., yiK ) . It is assumed that the yi are independently distributed and that:

$$f(y{ik}|\lambda{ik}) = \frac{\lambda{ik}^{y{ik}}e^{-\lambda{ik}}}{y{ik}!}, \quad i = 1, 2, ..., n, \quad k = 1, 2, ..., K,$$ (1)

where ik denotes the parameter of the Poisson distribution. To accommodate the over-dispersion that is usually exhibited in most traffic data sets, it is further assumed that ln(ik) = ln(ik) + εik, where ik is linked to a group of covariates (Xij) representing relevant traffic and road characteristics by

$$\ln(\mu{ik}) = \beta{k0} + \beta{k1}X{i1} + \dots + \beta{kJ}X{iJ}, \qquad (2)$$

and εi ⁼ (εi1, εi2, ..., εiK ) is distributed as:

$$\varepsilon_i \sim N_K(0, \Sigma).$$ (3)

For a univariate PLN model, where y represents just one random variable, e.g., total accidents, the subscript k is dropped as K = 1 and ˙ is replaced by ² in (3).

2.2. The scale-mixture model

The Student t density, which has heavier tails than the normal, provides a robust alternative to the normal in case outliers are suspected in the data. Thus, it is proposed to replace the multivariate normal distribution in (3) by a K-variate Student t distribution with zero mean, scale matrix ˙ and degrees of freedom:

$$\varepsilon_i \sim t_K(0, \Sigma, \nu).$$ (4)

However, the Student t distribution in (4) can be derived by mixing the multivariate normal in (3) with a scaling variable ωⁱ (Lange et al., 1989) so that

$$\varepsilon_i \sim N_K(0, \Sigma/\omega_i), \quad \omega_i \sim Gamma(\nu/2, \nu/2).$$ (5)

For a univariate PLN model, where K = 1, ˙ is replaced by ² in (4) and (5).

2.3. The two-component mixture

A robust two-component mixture model is proposed, where (3) is replaced by

$$\varepsilon_i \sim \pi N_K(0, \Sigma/\delta) + (1-\pi)N_K(0, \Sigma).$$ (6)

The first component, which represents the outlying cases, is variance inflated by the parameter ^ı ^∈ (0, 1), while ^∈ (0, 1) represents the prior probability of an outlier.

The estimation of the parameters and ı is discussed by Berkane and Bentler (1988). However, in practice, small values are usually specified for and ^ı. The choices ^∈ (0.05, ⁰.10) and ^ı ^∈ (0.1, ⁰.2) yield suitable models when the data are contaminated with a small fraction of outliers (Little, 1988; Lange et al., 1989; Albert, 1999). This contaminated normal model can be effective in the presence of extreme outliers.

For a univariate PLN model, where K = 1, ˙ is replaced by ² in (6). Alternatively, for a univariate NB model, (6) is replaced by the gamma mixture:

$$\lambda_i \sim \pi$$ Gamma $(\delta \kappa, \, \delta \kappa / \mu_i) + (1 - \pi)$ Gamma $(\kappa, \, \kappa / \mu_i),$ (7)

where represents the inverse dispersion parameter of the NB model.

2.4. Detecting outliers

For univariate normal regressionmodels, a large number of diagnostic measures have been developed in the statistical literature for detecting outliers and for measuring influence on parameter and precision estimates and on fitted and predicted values. Most of these measures are available in statistical packages such as SAS (SAS, 2002–2003; influence diagnostics for PROC MIXED).

On the contrary, few diagnostics have been developed for generalized linear models (GLMs) and fewer are available in the widely used statistical packages. For instance, an experimental statement has been recently developed by SAS (the ASSESS statement in PROC GENMOD) to assess the adequacy of the covariates' functional forms as well as that of the link based on cumulative residuals. To compute any of the other measures in a GLM setting, the analyst has to write/use specialized codes. For example, Sawalha and Sayed (2006) developed an algorithm to adapt Cook's distance for outliers' detection in NB regression.

For multivariate normal models, diagnostics measures were developed based on the Mahalanobis distance (Lange et al., 1989). Under (3) the Mahalanobis distance is given by

$$d_i^2 = (\lambda_i^* - \mu_i^)' \Sigma^{-1} (\lambda_i^ - \mu_i^*), \quad d_i^2 \sim \chi_K^2,$$ (8)

where -∗ ⁱ ⁼ (ln(ⁱ1), ln(ⁱ2), ..., ln(iK )) and ^∗ ⁱ ⁼ (ln(i1), ln(i2), ..., ln(iK )). Hence, it is proposed that values of the Mahalanobis distance (d² ⁱ ) exceeding the critical value of the 2-distribution with K degrees of freedom and an appropriately selected significance level ˛ be identified as potential outliers.

2.5. Down-weighting outliers

The Mahalanobis distance is closely related to the posterior mean of the scaling variable ωⁱ in (5). In fact, going along the lines of Little (1988), it can be verified that

$$E(\omega_i|y_i,\Theta) = \frac{\nu + K}{\nu + d_i^2},\tag{9}$$

where denotes the set of all parameters under consideration. Thus, under the scale-mixture model, potential outliers are associated with small weights, particularly those ωⁱ that are significantly under 1, with the degree of down-weighting being enhanced for smaller values of (Lange et al., 1989; Congdon, 2006).

It should be noted that the scaling weights in (6) represent a discrete special case of those in (5), where ωⁱ = ı with probability -, and ^ωⁱ = 1 with probability 1 ⁻ -. Following Little (1988), it can be verified that

$$E(\omega_i|y_i,\Theta) = \frac{1 - \pi + \pi \delta^{1+K/2} e^{(1-\delta)d_i^2/2}}{1 - \pi + \pi \delta^{K/2} e^{(1-\delta)d_i^2/2}}.$$ (10)

However, as noted by Little (1988), the distributions of the posterior weights in (9) and (10) are quite different since the scale-mixture model produces relatively smoothly declining weights, while the two-component mixture (contaminated normal) model tends to concentrate the low weights in a few outlying cases.

It is important to note that (i) under the outlier ignoring approach outlying observations are assigned the same weights as those of other points in the data set, (ii) under the outlier rejection approach outlying observations are assigned zero weights, and (iii) under the robust mixture models outlying observations are assigned smaller weights, where the weight of each observation is inversely proportional to its Mahalanobis distance. Thus, the larger is the distance, the less significant is the weight.

2.6. Prior distributions

To obtain the full Bayes estimates of the unknown parameters, it is required to specify prior distributions for the hyper parameters ˇkj, ˙ and .

The prior distributions for ˇkj and ˙, are discussed in El-Basyouny and Sayed (2009a). The most commonly used priors are diffused normal distributions (with zero mean and large variance) for the regression parameters and a Wishart(P, r) prior for ˙−1, where P and r ≥ K represent the prior guess at the order of magnitude of the precision matrix ˙−¹ and the degrees of freedom, respectively. Choosing r = K as the degrees of freedom corresponds to vague prior knowledge.

For univariate PLN and NB models, the uninformative priors Gamma(ε, ε) or Gamma(1, ε) are usually used for −² and , where ε is a small number, e.g., 0.01 or 0.001.

Geweke (1993) noted that uninformative prior distributions for can be troublesome and proposed an exponential distribution of the form ^f() = exp(−) as a prior for . Since ^E() = ⁻1, then can be assumed to follow a uniform distribution over the interval (0.02, 0.50), which corresponds to the range = 2 to = 50 for effective normality (Congdon, 2006). Yet for small samples, Lange et al. (1989) advocated the use of a preset value = 4.

As for the parameters of the two-component mixture model, a limited sensitivity analysis was conducted comprising the four combinations: ( = 0.1, ı = 0.1), ( = 0.1, ı = 0.2), ( = 0.05, ı = 0.1) and ( = 0.05, ı = 0.2). Using the present data set, the results showed that the models with ı = 0.1 provided better fits than those with ı = 0.2 (DIC = 1355 versus 1359). For ı = 0.1, there were small differences in the estimates of the regression and variance–covariance parameters (up to 4.5%) between the model with = 0.1 and that with = 0.05. Hence, the combination ( = 0.1, ı = 0.1) was used throughout the study.

In addition, the following priors were used ^ˇkj ^∼ ^N(0, 1002), ˙ ^∼Wishart(I, 2), where ^I is the 2 ^× 2 identity matrix, ^∼ Exponential() and ^∼ Uniform(0.02, 0.50).

2.7. Full Bayes estimation

The posterior distributions needed in the full Bayes approach were sampled using the MCMC techniques available in WinBUGS (Lunn et al., 2000). The techniques generate sequences (chains) of

Fig. 1. Scatter plot of Mahalanobis distance by site. Left: MVPLN model (dashed line corresponds to 2(2, 0.01) = 9.210). Right: scale-mixture model (dashed line corresponds to 2*F(2, 4.564, 0.15) = 6.328).

random points, whose distributions converge to the target posterior distributions. A sub-sample is used to monitor convergence and then excluded as a burn-in sample. The remaining iterations are used for parameter estimation, performance evaluation and inference.

Monitoring convergence is important because it ensures that the posterior distribution has been "found". Thereby indicating when parameters sampling should begin. To check convergence, two or more parallel chains with diverse starting values are tracked to ensure full coverage of the sample space. Convergence of multiple chains is assessed using the Brooks–Gelman–Rubin (BGR) statistic (Brooks and Gelman, 1998). A value under 1.2 of the BGR statistic indicates convergence. Convergence is also assessed by visual inspection of the MCMC trace plots for the model parameters as well as by monitoring the ratios of the Monte Carlo errors relative to the respective standard deviations of the estimates; as a rule of thumb these ratios should be less than 0.05.

2.8. Model comparisons

The Deviance Information Criteria (DIC) is used for model comparisons (Spiegelhalter et al., 2002). As a goodness-of-fit measure, DIC is a Bayesian generalization of Akaike's Information Criteria (AIC) that penalizes larger parameter models. According to Spiegelhalter et al. (2005), it is difficult to determine what would constitute an important difference in DIC. Very roughly, differences of more than 10 might definitely rule out the model with the higher DIC. Differences between 5 and 10 are considered substantial. If the difference in DIC is less than 5, and the models make very different inferences, then it could be misleading just to report the model with the lowest DIC.

3. Data and model description

A sample of n = 99 signalized intersections in the city of Edmonton were investigated for the purpose of developing accident prediction models relating the safety of urban 4-leg intersections to their traffic flows. This data set was used as a reference group for the evaluation of the City of Edmonton's Intersection Safety Camera program (Sayed and de Leur, 2007). The data on accident frequencies and traffic volumes within the intersections were obtained from the city of Edmonton. The traffic and accident data were provided for a 3-year period before the implementation of the program and were selected to ensure that no changes in either the intersection layout or traffic volumes occurred over the course of the 3 years. The accidents were classified into K = 2 categories: property damage only (PDO), k = 1, and injuries and fatalities (I + F), k = 2.

Average values for the traffic volumes (over the 3 years period) were used to build SPFs for predicting aggregate number of accidents using the link function:

$$\ln(\mu{ik}) = \beta{k0} + \beta{k1} X{i1} + \beta{k2} X{i2}, \tag{11}$$

where ⁱ¹ = ln(Vi1), Vi¹ is the annual average daily traffic (AADT) at the major approach, and ⁱ² = ln(Vi2), Vi² is the AADT at the minor approach. A detailed description of the data is given in El-Basyouny and Sayed (2009a).

Note that Eq. (11) is based on a simple AADT model. AADT models are increasingly popular among practitioners, e.g., they are often preferred over models that include several covariates because they can be easily re-calibrated (Lord et al., 2008) or transferred. Further, such models are adopted for estimating safety performance in the forthcoming highway safety manual (Hughes et al., 2005). However, AADT models are susceptible to omitted variables bias. Although, analysts should always seek to use correctly specified models, sometimes data limitations force the omission of important relevant covariates. Therefore, if an analyst is planning to use a simple "volume" model (owing to data limitations or otherwise), then it is crucial to check and account for potential outliers as they contradict the exchangeability assumption.

4. Results and discussion

The posterior summaries for the fitted models were obtained via two chains with 20,000 iterations 10,000 of which were excluded as a burn-in sample using WinBUGS. For all fitted models, examination of the BGR statistics, ratios of the Monte Carlo errors relative to the standard deviations of the estimates and trace plots for all model parameters indicated convergence.

4.1. Outliers' detection and weights

After fitting the MVPLN model, a visual inspection of diagnosis plots of the residuals versus the predicted values for I + F, PDO and the total (I + F + PDO) suggests that the data may have several outliers corresponding to sites 4, 5, 35, 38 and 51. Checking the data collection and recording process showed no errors in the respective data points and further investigation of site characteristics revealed no genuine differences between these particular sites and the rest of the sites. Hence, the proposedmixturemodels were used to identify and down-weight potential outliers.

The posterior estimates of the Mahalanobis distances in (8) are displayed in Fig. 1(left) for MVPLN along with a dashed line indicating the critical value of 2(K = 2, ˛= 0.01) = 9.210. Three outliers are easily identified, namely, those corresponding to sites 54, 4 and 5. These are the three highest points with d² ⁱ exceeding the 2-critical value. It should be noted that this list of outliers is different from the preliminary list obtained by the quick visual inspection.

Table 1 Parameter estimates and their associated standard deviations.

	MVPLN	Two-component mixture	Scale mixture	Outlier rejecting MVPLN
PDO
Intercept	−10.680 ± 1.347	−11.320 ± 1.229	−11.070 ± 1.201	−11.560 ± 1.101
Ln(V1)	0.902 ± 0.138	0.871 ± 0.118	0.908 ± 0.123	0.853 ± 0.113
Ln(V2)	0.531 ± 0.104	0.633 ± 0.096	0.569 ± 0.101	0.676 ± 0.089
11	0.163 ± 0.030	0.100 ± 0.019	0.097 ± 0.022	0.090 ± 0.017
I + F
Intercept	−9.346 ± 1.565	−9.793 ± 1.446	−9.582 ± 1.491	−9.978 ± 1.315
Ln(V1)	0.742 ± 0.159	0.721 ± 0.145	0.750 ± 0.153	0.685 ± 0.133
Ln(V2)	0.513 ± 0.121	0.585 ± 0.121	0.535 ± 0.125	0.643 ± 0.111
22	0.217 ± 0.042	0.142 ± 0.028	0.135 ± 0.032	0.134 ± 0.026
Covariance
12	0.143 ± 0.029	0.078 ± 0.018	0.074 ± 0.021	0.065 ± 0.016
Correlation
	0.758 ± 0.056	0.651 ± 0.076	0.644 ± 0.083	0.592 ± 0.086
			4.564 ± 1.737
DIC	1373	1355	1355	1322

To confirm that the proposed mixture models represent robust alternatives to MVPLN, which are capable of accommodating outliers, the posterior estimates of the Mahalanobis distances are displayed in Fig. 1(right) for the scale-mixture model. Under (5) the Mahalanobis distance is given by d˜2 ⁱ ⁼ ωid² ⁱ and ^d˜2 ⁱ /K∼F(K, ). Using the estimate ˆ ⁼ ⁴.564 from Table 1, the dashed line in Fig. 2 (right) indicates twice the critical value of the F-distribution, namely, ^c ⁼ ² ^∗ ^F(^K ⁼ ², ^ˆ ⁼ ⁴.564, ˛ ⁼ ⁰.15) ⁼ ⁶.328. It is seen that all Mahalanobis distances d˜2 ⁱ fall below c, even at such a high significance level as ˛= 0.15. The scaling variable ωⁱ was effective in down-weighting the outlying observations.

Fig. 2 depicts two scatter plots of the scaling weights by site for the two mixture models. The three outliers are easily identified as they are the three lowest points in both plots with ωⁱ far under 1. Note also that the scale-mixture produced relatively smoothly declining weights (left), while the two-component mixture (contaminated normal) model concentrated the low weights in a few outlying cases (right).

Since the 2- and F-critical values are readily available, the use of Mahalanobis distances for identifying outliers is preferable to the use of scaling weights as there are no specific guidelines about the cut-off value of the latter.

Consider the outliers (sites 54, 4 and 5). These sites have observed accidents that are either much lower or higher than expected. Important missing variables such as pavement friction, angles of the approaches, signal timing, average approach speeds, etc. would likely reduce their Mahalanobis distances. However, lacking such data, the robust mixture models have down-weighted these outliers to reduce their impact on model estimates and predictions.

4.2. Parameter estimates

Table 1 summarizes the parameter estimates and their associated standard deviations using the proposed mixture models as well as the MVPLN. It should be noted that while these models can handle more than two accident categories, the current application involves only two severity levels leading to bivariate distributions. The results in the second column of Table 1 (for MVPLN) are the same as those given in El-Basyouny and Sayed (2009a). Such results correspond to the approach that ignores potential outliers.

The results of Table 1 show that: (i) the parameter estimates are significant; (ii) the regression coefficients are positive, indicating an increase in the mean accident frequency with traffic volumes; (iii) the impact of traffic volumes at the major approach is larger than that at the minor approach; and (iv) traffic volumes have a larger effect on PDO accidents than on I + F accidents.

4.3. Results of the two mixture models

Table 1 reveals small differences in the estimates of the regression parameters between the various multivariate models. However, the estimates of the extra-Poisson variation parameters (¹¹ and 22) were considerably smaller under the mixture models compared with MVPLN. Using Eq. (18) in El-Basyouny and Sayed (2009a), it is seen that

$$\frac{Var(Y{ik}|\text{Mixture Model})}{Var(Y{ik}|\text{MVPLN})} \cong \frac{\exp(\hat{\sigma}_{kk}/\hat{\omega}i) - 1}{\exp(\hat{\sigma}{kk}) - 1},$$ (12)

where ˆkk and ˆωi are the posterior estimates under the appropriate mixture model and ˜kk is the posterior estimate under MVPLN. As precision is inversely proportional to the variance, it is estimated that the two-mixture model is about 1.7 times as much precise as MVPLN (disregarding the three potential outliers). The actual precision ratios were 1.74 and 1.67 for PDO and I + F, respectively. For both categories under the scale-mixture model, the corresponding ratios ranged from 0.8 to 2.3 and exceeded 1.7 in the majority of sites (≥90%). The improvement in precision is due mainly to the down-weighting of outlying data.

The estimate of 0.758 for the correlation between PDO and I + F under MVPLN has been slightly reduced under the mixture models, but remains highly significant nevertheless.

The estimate of 4.564 for under the scale-mixture model supports the conjecture of Lange et al. (1989)that the preset value = 4 can be safely used in small samples.

The DIC statistics were 1373, 1355 and 1355 under MVPLN, the two-component mixture and the scale mixture, respectively. Thus, both mixture models have significantly outperformed MVPLN.

4.4. Results of the outliers rejecting approach

The results of applying MVPLN to the remaining 96 observations, after excluding the 3 potential outliers, are given in the last column of Table 1. Although the estimates of the regression parameters under the outlier rejecting MVPLN are not significantly different from those of the other models, the estimates of the variance–covariance parameters were considerably smaller compared with MVPLN.

Naturally, the exclusion of outliers has resulted in a substantially better fit where DIC was significantly reduced to 1322. However, data points (observations) cannot be eliminated just to obtain a

Fig. 2. Left: scatter plot of scaling weights by site for the scale-mixture model. Right: scatter plot of scaling weights by site for the two-component mixture model.

better fit. As pointed out by Lange et al. (1989), the outlier rejecting approach tends to underestimate uncertainty by producing too small standard deviations for the parameter estimates, which may lead to incorrect conclusions. For the present data set, the standard deviations resulting from the outlier rejecting MVPLN were much smaller relative to MVPLN. In fact, the standard deviations of the regression parameters were up to 18% smaller than MVPLN under the outlier rejecting MVPLN compared with a 10% decrease for the mixture models. The situation is more evident for the variance–covariance parameters where the reductions in standard deviations were up to 46% and 30% under the outlier rejecting MVPLN and the mixture models, respectively. Thus, while the two mixture models have yielded moderately smaller standard deviations of the parameter estimates (relative to MVPLN), the corresponding reductions in the standard deviations under the outlier rejecting approach were extreme. Hence, if there were no data collection error reasons to justify the elimination of outliers, the outlier rejecting approach would underestimate the uncertainty associated with the data set and risks reaching erroneous conclusions.

5. Conclusions

Mixture models, namely, a scale-mixture model and a twocomponent mixture (contaminated normal) model, were proposed as robust alternatives to MVPLN regression for the joint analysis of a sample of accident counts classified by two severity levels (PDO and I + F) at 99 urban intersections.

For such a multivariate GLM setting, the Mahalanobis distance was used to identify potential outliers. It was also shown that the proposed robust models were successful in accommodating these outliers. In terms of goodness-of-fit, both mixture models have outperformed MVPLN. Furthermore, while fitting the mixture models may take a little longer (in terms of convergence), it is not markedly complicated (in terms of WinBUGS coding) than MVPLN.

The results revealed small differences in the estimates of the regression parameters between the various multivariate models. Yet, the estimates of the extra-Poisson variation parameters were considerably smaller under the mixture models compared with MVPLN. As precision is inversely proportional to the variance, it is estimated that the mixture models are at least 1.7 times as much precise as MVPLN due to the down-weighting of outlying data. This improved precision would (i) reduce the rather high uncertainty associated with the predicted values, if one adopts the outliers ignoring approach and (ii) avoids the risk of underestimating this uncertainty, which results from using the outliers rejecting approach.

Although the estimate of the correlation between PDO and I + F under MVPLN has been slightly reduced under the mixture models, it remains highly significant indicating that higher PDO rates are associated with higher I + F rates, as the accident likelihood for both types is likely to rise due to similar deficiencies in roadway design and/or other unobserved factors.

In terms of model application, the paper introduced robust multivariate techniques that are capable of accommodating outliers. The proposed techniques are based on the MVPLN regression approach that has been gaining attention recently in the traffic safety literature. The results showed that potential outliers could be overlooked if the analysis is restricted to MVPLN.

The potential outliers should be carefully investigated. It is recommended that the proposed outlier resistance modeling techniques be used unless the exclusion of the outlying observations can be justified because of data related reasons (data collection errors).

The degree of "success" of the outlier accommodating approach will (to some extent) depends on the data set analyzed. The extent of the omitted variables bias will likely differ from one data set to the next. Although the present results conform to those in the literature, further research with different data sets is required to confirm the paper's findings. Future topics of research might also include adapting the diagnosticmeasures that have been developed for univariate normal regression models to GLM settings as well as developing similar measures for multivariate models.

References

Aguero-Valverde, J., Jovanis, P.P., 2009. Bayesian multivariate Poisson log-normal models for crash severity modeling and site ranking. In: Presented at the 88th Annual Meeting of the Transportation Research Board.

Albert, J.H., 1999. Criticism of a hierarchical model using Bayes factors. Statistics in Medicine 18, 287–305.

Anastasopoulos, P.Ch., Mannering, F., 2009. A note on modeling vehicle accident frequencies with random-parameters count models. Accident Analysis and Prevention 41 (1), 153–159.

Berkane, M., Bentler, P.M., 1988. Estimation of contamination parameters and indentification of outliers in multivariate data. Sociological Methods and research 17, 55–64.

Brooks, S.P., Gelman, A., 1998. Alternative methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics 7, 434–455.

Congdon, P., 2006. Bayesian Statistical Modeling, 2nd ed. Wiley, New York.

El-Basyouny, K., Sayed, T., 2006. Comparison of two negative binomial regression techniques in developing accident prediction models. Transportation Research Record 1950, 9–16.

El-Basyouny, K., Sayed, T., 2009a. Collision prediction models using multivariate Poisson-lognormal regression. Accidents Analysis and Prevention 41 (4), 820–828.

El-Basyouny, K., Sayed, T., 2009b. Accident prediction models with random corridor parameters. Accidents Analysis and Prevention 41 (5), 1118–1123.

Geweke, J., 1993. Bayesian treatment of the independent Student-t linear model. Journal of Applied Econometrics 8, S19–S40.

Heydecker, B.G., Wu, J., 2001. Identification of sites for accident remedial work by Bayesian statistical methods: an example of uncertain inference. Advances in Engineering Software 32, 859–869.

Hughes, W., Eccles, K., Harwood, D., Potts, I., Hauer, E., 2005. Development of a Highway Safety Manual. Appendix C: Highway Safety Manual Prototype Chapter: Two-Lane Highways. NCHRP Web Document 62 (Project 17-18(4)), Washing-

ton, DC, Available from http://onlinepubs.trb.org/onlinepubs/nchrp/nchrp w62. pdf.
Kulmala, R., 1995. Safety at Rural Three-and Four-arm Junctions: Development of Accident Prediction Models. Technical Research Centre of Finland, VTT 233, Espoo.
Kumara, S.S.P., Chin, H.C., 2003.Modeling accident occurrence at signalized tee intersections with special emphasis on excess zeros. Traffic Injury Prevention 3 (4), 53–57.
Lange, K.L., Little, R.J.A., Taylor, J.M.G., 1989. Robust statistical modeling using the t distribution. Journal of the American Statistical Society 84, 881–896.
Little, R.J.A., 1988. Robust estimation of the mean and covariance matrix from data with missing values. Applied Statistics 37, 23–39.
Lord, D., Guikema, S.D., Geedipally, S.R., 2008. Application of the Conway–Maxwell–Poisson generalized linear model for analyzing motor vehicle crashes. Accident Analysis and Prevention 40, 1123–1134.
Lord, D., Park, Y.-J., 2008. Investigating the effects of the fixed and varying dispersion parameters of Poisson-gamma models on empirical Bayes estimates. Accident Analysis and Prevention 40 (4), 1441–1457.
Lord, D., Washington, S.P., Ivan, J.N., 2005. Poisson, Poisson-gamma and zero inflated regression models of motor vehicle crashes: balancing statistical fit and theory. Accident Analysis and Prevention 37 (1), 35– 46.
Lunn, D.J., Thomas, A., Best, N., Spiegelhalter, D., 2000. WinBUGS—a Bayesian modelling framework: concepts, structure, and extensibility. Statistics and Computing 10, 325–337.
Ma, J., Kockelman, K.M., Damien, P., 2008. A multivariate Poisson-lognormal regression model for prediction of crash counts by severity, using Bayesian methods. Accident Analysis and Prevention 40, 964–975.
Malyshkina, N.V.,Mannering, F.L., Tarko, A.P., 2009.Markov switching negative binomial models: an application to vehicle accident frequencies. Accident Analysis and Prevention 41 (2), 217–226.
Miaou, S.P., Lord, D., 2003. Modeling traffic crash-flow relationships for intersections: dispersion parameter, functional form, and Bayes versus empirical Bayes. Transportation Research Record 1840, 31–40.
Miranda-Moreno, L.F., Fu, L., Saccomanno, F.F., Labbe, A., 2005. Alternative risk models for ranking locations for safety improvement. Transportation Research Record 1908, 1–8.
Mitra, S., Washington, S., 2007. On the nature of over-dispersion in motor vehicle crash prediction models. Accident Analysis and Prevention 39, 459–468.
Park, B.-J., Lord, D., 2009. Application of finite mixture models for vehicle crash data analysis. Accident Analysis and Prevention 41, 683–691.
Park, E.S., Lord, D., 2007.Multivariate Poisson-lognormalmodels for jointlymodeling crash frequency by severity. Transportation Research Record 2019, 1–6.
Qin, X., Ivan, J.N., Ravishankar, N., 2004. Selecting exposure measures in crash rate prediction for two-lane highway segments. Accident Analysis and Prevention 36 (2), 183–191.
Relles, D.A., Rogers, W.H., 1977. Statisticians are fairly robust estimators of location. Journal of the American Statistical Society 72, 107–111.
SAS Institute Inc., 2002–2003. Version 9.1, PROC MIXED, PROC GENMOD. SAS Institute Inc., Cary, USA.
Sawalha, Z., Sayed, T., 2001. Evaluating safety of urban arterial roadways. Journal of Transportation Engineering 127 (2), 151–158.
Sawalha, Z., Sayed, T., 2006. Traffic accidents modeling: some statistical issues. Canadian Journal of Civil Engineering 33 (9), 1115–1124.
Sayed, T., de Leur, P., 2007. Evaluation of Edmonton's intersections safety camera programs. Transportation Research Record 2009, 37–45.
Shankar, V., Milton, J., Mannering, F.L., 1997. Modeling accident frequency as zero-altered probability processes: an empirical inquiry. Accident Analysis and
Spiegelhalter, D.J., Best, N.G., Carlin, B.P., Van der Linde, A., 2002. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society B 64, 1–34.

Prevention 29 (6), 829–837.

Spiegelhalter, D., Thomas, A., Best, N., Lunn, D., 2005. WinBUGS User Manual. MRC
Biostatistics Unit, Cambridge, Available from http://www.mrc-cam.ac.uk/bugs. Tunaru, R., 2002. Hierarchical Bayesian models for multiple count data. Austrian Journal of Statistics 31 (3), 221–229.
Ye, X., Pendyala, R.M., Washington, S.P., Konduri, K., Oh, J., 2009. A simultaneous equations model of crash frequency by collision type for rural intersections. Safety Science 47 (3), 443–452.

elbasyounymethodaccountoutliersdevelopment2010.md 45 KB Verlauf Originalformat