category: literaturenote citekey: elbasyounyapplicationgeneralizedlinkfunctions2010 title: Application of generalized link functions in developing accident prediction models authors: "El-Basyouny, Karim; Sayed, Tarek" year: 2010 date: 2010-00-00 2010 doi: 10.1016/j.ssci.2009.09.007 publication: Safety Science url: "https://doi.org/10.1016/j.ssci.2009.09.007" zotero_key: FDHG26KS zotero_storage: 5EWGTQW7 collections: transpordi_planeerimine folder: 001_artiklid/Road Safety Indicators firstAuthor: "El-Basyouny, Karim"

status: converted

Contents lists available at ScienceDirect

Safety Science

Application of generalized link functions in developing accident prediction models

Karim El-Basyouny *, Tarek Sayed ¹

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

article info

Article history: Received 29 September 2008 Received in revised form 8 September 2009 Accepted 16 September 2009

Keywords: Canonical link Full Bayes estimation Generalized link Markov Chain Monte Carlo Poisson–gamma hierarchy

abstract

The Bayesian Poisson–Gamma hierarchy, leading to the negative binomial distribution, has been the standard practice in developing accident prediction models. To linearize the relationship connecting the mean of the negative binomial distribution to relevant covariates, a canonical log link has traditionally been used. Typically, little information is available regarding the choice of a particular link. To avoid link misspecification, it is proposed to nest the canonical log link model within a generalized link family and subsequently use the full Bayes method for parameter estimation, performance evaluation and inference. The proposed approach was applied to a sample of accident and traffic volume data corresponding to 99 intersections in the city of Edmonton, Alberta. The results showed that both the generalized link model and the traditional canonical link model provided adequate fit to the data. However, the Bayes factor provided a clear statistical support for the use of the generalized link approach. A procedure for link validation is also described. It allows the users (e.g., road authorities) to consider the changes in predicted accidents that will result if a generalized link is used instead of a canonical link. If a certain maximal change is tolerated, the canonical link can be used to analyze the data; otherwise the generalized link is worth the extra efforts and should be adopted. When compared with the traditional approach, the generalized link model was found to predict a lower number of accidents whenever there is a heavy traffic at the major approach, especially if combined with light flow on the minor approach. The paper concludes by identifying out areas for further research.

Crown Copyright -2009 Published by Elsevier Ltd. All rights reserved.

1. Introduction

The application of accident prediction models (APMs) in assessing the safety of a road entity (intersection, road segment, etc.) has become a standard practice among safety researchers and practitioners. An APM is a mathematical model that relates the expected accident frequency at a road entity to the various traffic and geometric characteristics of that entity. Accident prediction models have several road safety applications such as: evaluating the safety of various road facilities, identifying and ranking accident prone locations and evaluating the effectiveness of road safety improvement measures (Sawalha and Sayed, 2006).

Current practice utilizes the negative binomial (NB) model, which results from a Poisson–Gamma hierarchy, to develop APMs. The log function (canonical link) is usually used to connect the random (accidents) and systematic (traffic and geometric characteristics) components of the NB model to produce a linear predictor. Nevertheless, there is often no or little information available concerning the appropriate link to use in a given application. An improper choice of link may result in a link misspecification leading to substantial bias in the model coefficients and the mean response estimates (Czado and Munk, 2000; Czado and Raftery, 2006; and the references therein). To avoid this problem, link families have been considered by Pregibon (1980) and Czado (1992). Such families generalize the canonical link through the use of an extra parameter h > 0 characterizing the generalized link. The basic idea has been to embed the canonical log link of the NB model, which corresponds to the choice h = 1, in a wide parametric class of links and allow for a data selected choice of the link function.

This paper proposes that generalized links are considered while developing APMs. If the canonical link is validated against the generalized link, then it can be used, otherwise the generalized link is to be adopted. The proposed approach is applied to a sample of accident and traffic volume data corresponding to 99 intersections in the city of Edmonton, Alberta. The objective is to compare the generalized and canonical link functions with respect to accident prediction. The statistical software WinBUGS (Lunn et al., 2000), which provides MCMC sampling for a variety of Bayes models, was used in the analysis.

The rest of the paper is organized as follows: Section 2 describes previous work. A family of generalized links is introduced in Section 3, where the issues of link misspecification and validation are discussed. Section 4 describes the methodology used for

^* Corresponding author. 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.

estimating the models as well as goodness of fit measures. The data set and the model used in the analysis are described in Section 5. The results are presented in Section 6 and are discussed in Section 7. Section 8 presents a summary and conclusions.

2. Previous work

Some of the most commonly applied statistical models to accident data include: binomial, Poisson, Poisson–gamma (negative binomial), Zero-Inflated Poisson (ZIP), Zero-Inflated Negative Binomial (ZINB), and multinomial probability models (Lord et al., 2005). They showed that the Poisson and NB models serve as statistical approximations to the accident process, but since over-dispersion is exhibited in most accident data the APM models based on NB are usually preferred over Poisson regression models.

Let Y be the random variable that represents the accident frequency at a given location during a specific time period, and let y be a certain realization of Y. The mean of Y, denoted by K, is itself a random variable. For K = k, Y is Poisson distributed with parameter k. It is the usual practice to assume that the distribution of K can be described by a gamma probability density function with shape parameter j and scale parameter j/l. Then, the distribution of Y is negative binomial with a mean l and variance l + l² /j.

Several other distributions were proposed for developing APMs. For instance, assuming a normal distribution for K leads to the Poisson-lognormal model. Since the lognormal tails are asymptotically heavier than those of the Gamma distribution (Kim et al., 2002), the Poisson–lognormal is a good candidate for modeling accident rates in the presence of outliers. Lord et al. (2008) used the Conway–Maxwell–Poisson distribution to develop APMs. Such a distribution is a generalization of the Poisson that can cope with both over- and under-dispersion. They found that the generalized linear models based on the Conway–Maxwell–Poisson distribution perform as well as NB models in terms of goodness of fit and predictive performance.

Recently, researchers have tried to address a problem arising from the ''excess" zeroes in accident data, by applying zero-inflated probability models such as ZIP and ZINB. The models assume that a dual-state process ''safe" and ''unsafe" is responsible for generating the accident data. The general consensus indicated an improved fit to data compared to the regular Poisson and NB models. Lord et al. (2005) criticized the dual-state approach and argued that ''excess" zeroes arise from a combination of low exposure, high heterogeneity and sites categorized as high risk; inappropriate selection of time/space scales; under-reporting; and/or omitting important covariates. They recommended that a careful selection of the time/space scales for analysis, including an improved set of explanatory variables and/or unobserved heterogeneity effects in count regression models, or applying small area statistical methods (for data characterized by low exposure) represent the most defensible modeling approaches for datasets with a preponderance of zeros.

Developing APMs using data characterized by a low sample mean is problematic, especially if combined with a small sample size due to the expensive costs of collecting accident data and other relevant variables. The low mean problem has been studied in the traffic safety literature, see Maher and Summersgill (1996), who showed how such a problem affects the goodness of fit of APMs, Wood (2002), who proposed a method to test the fit of APMs developed using data exhibiting a low mean, and Lord (2006), who confirmed that a low sample mean combined with a small sample size can seriously affect the estimation of the dispersion parameter and undermine common analyses performed in highway safety.

The majority of APMs were developed using a NB model with a fixed dispersion parameter. Such an assumption was challenged and various dispersion parameter relationships were examined (El-Basyouny and Sayed, 2006). In particular, Mitra and Washing ton (2007) indicated that the modeling of the dispersion parameter depends on how the mean structure is modeled. It appears that extra-variation is a function of covariates when the mean structure is poorly specified and suffers from omitted variables. In contrast, when sufficient explanatory variables are used to model the mean extra-Poisson variation is not significantly related to these variables.

Let Xj, j = 1, ... , p, represent the independent variables (e.g., traffic flow and other site characteristics). The log function, which is the canonical link for the NB distribution, is then used to relate the mean l to the covariates

$$ln(\mu) = \beta_0 + \beta_1 X_1 + \dots + \beta_p X_p. \tag{1}$$

The functional form in Eq. (1) has been extensively reviewed (Miaou and Lord, 2003). The most common functional form used for developing APMs for road intersections is

$$ln(\mu) = ln(\alpha_0) + \alpha1 \, ln(V{\it mj}) + \alpha2 \, ln(V{\it mi}) + \sum_{j=1}^{\it m} \gamma_j X_j, \eqno(2)$$

where Vmj is the annual average daily traffic (AADT) on the major approach, Vmi is the AADT on the minor approach, the Xj are covariates representing certain traffic and geometric characteristics, and a0, a1, a2, c1, ... , c^p are model parameters. The corresponding form used for developing APMs for road segments is

$$\ln(\mu) = \ln(\alpha_0) + \alpha_1 \ln(L) + \alpha2 \ln(V) + \sum{i=1}^{m} \gamma_i X_i, \tag{3}$$

where L is segment length and V is the annual average daily traffic (AADT). Miaou and Lord (2003) examined the mathematical limitations of the various functional forms, particularly their properties at the boundaries. The authors demonstrated that, for a given data set, a large number of plausible functional forms with almost the same overall statistical goodness of fit are possible. Furthermore, they introduced an alternative class of logical formulations that may enable a richer interpretation of the data.

Other methods for developing APMs that were recently researched include: the application of the multivariate Poisson-lognormal (MVPLN) to model accident count data at different levels of severity (El-Basyouny and Sayed, 2009a), the use of a two-state Markov Switching model (Malyshkina et al., 2009) to analyze collision frequencies, and the use of models with random parameters (El-Basyouny and Sayed, 2009b).

3. Generalized links

As we have seen above, previous work in the traffic safety literature has investigated many aspects of APMs, except the link function. Canonical links in generalized linear models, such as the log link in NB regression, are used to linearize the relationship between the predicted value and the predictors. They are popular among practitioners as they provide maximum information and simple interpretation of the regression parameters (McCullagh and Nelder, 1989). Yet, a canonical link may not furnish the best fit available for a given data set.

It has been shown in the statistical literature that the use of an inappropriate link results in a link misspecification, which might lead to a substantial bias in the model estimates (Czado and Munk, 2000; Czado and Raftery, 2006; and the references therein). Noting that the choice of the link function is often made arbitrarily, recent research has: (i) permitted the data to estimate the link function by incorporating it as an unknown in the model (Mallick and Gelfand, 1994), (ii) described full Bayesian inference for generalized linear models where uncertainty exists about the structure of the linear predictor, the linear parameters and the link function (Ntzoufras et al., 2003), and (iii) used a non-parametric approach to estimate a generalized additive model with unknown link function (Horowitz and Mammen, 2007).

To guard against link misspecification, the canonical link is usually embedded into a wide class of links. One such useful family of link functions is given by Czado (1997), Czado and Raftery (2006)

$$\mu(\eta, \theta) = \exp{[(1 + \eta)^{\theta} - 1]/\theta}; \quad \eta \geqslant 0, \theta > 0,$$ (4)

where g = Xb, X = (1, X1, ... , Xp) and b ¼ ðb0; b1; - ; bpÞ 0 . In the generalized linear models terminology, Eq. (4) is actually an inverse link function corresponding to the following (forward) link function

$$g(\mu, \theta) = (\theta \ln(\mu) + 1)^{1/\theta} - 1; \quad \mu \ge 0, \theta > 0.$$ (5)

Regarding notation, the simplified notation lh will be used instead of l(g, h) unless it is needed to stress the dependence of l on g. Note also that as h increases away from 1, lh increases away from l¹ on one hand, while the scale of the regression parameters bⁱ is reduced on the other hand as a tradeoff in order to adequately fit the data. A similar tradeoff (albeit reversed) occurs as h decreases away from 1. Thereby, the notation bh will be used occasionally to highlight the fact that the regression parameters are on different scales for different link functions. Also, the notation g(b) will be used sometimes to emphasize the dependence of the linear predictor g on b.

Let M^h denote the NB APM based on g(l, h) and note that the canonical log link is the special case h = 1, corresponding to M1. Note also that M^h is a generalized linear model when h is known, while this is no longer true if h, b and j are jointly estimated by the data.

Eq. (4) defines a family of power transformations to impose linearity, i.e., to make the linear model more appropriate to the data in the sense that the transformation is approximately linear over the range of observations. To illustrate, let us consider the consequences of fitting a model with the canonical link g(l, 1) (=ln (l)) when the correct but unknown link function is indeed g(l, h ). Using a first order Taylor series expansion about the canonical link, we have the approximate equation

$$g(\mu, \theta^) \cong \ln(\mu) + (\theta^ - 1)[\ln(\mu) - \ln(\ln(\mu) + 1)(\ln(\mu) + 1)].$$ (6)

The correct link function g(l, h ) = g = Xb may thus be approximated by

$$ln(\mu) = X\beta + Z\delta,$$ (7)

with an obvious definition for Z and d = (h 1). The term Zd, which describes the local differences between the two link functions, represents the misspecification of the model using the wrong link function.

We propose to test M¹ against M^h so that the traditional model M¹ is used only if supported by the data, i.e., if the hypothesis h = 1 is plausible. Otherwise M^h should be adopted. Such a testing problem is discussed in the following subsection.

3.1. Link validation

The usual practice of treating model selection (M¹ vs. Mh) as a testing problem of a point null hypothesis h = 1 was criticized by many researchers (Czado and Munk, 2000). Alternatively, they proposed a hypothesis of the form

$$H_0: \theta \notin [1 - \theta_l, 1 + \theta_u]$$ vs. $H_1: \theta \in [1 - \theta_l, 1 + \theta_u],$ (8)

where (1 hl, 1+ hu) is a neighborhood of h within which the canonical link is validated. Eq. (8) represents the validation problem and its converse (H¹ vs. H0) represents the discrimination problem. The following p-value was then derived

$$p(\hat{\theta}, s, \theta_l, \theta_u) = \Phi{(\hat{\theta} - 1 - \theta_u)/s} - \Phi{-(\hat{\theta} - 1 + \theta_l)/s}, \tag{9}$$

where ^h and s denote the maximum likelihood estimate of h and its standard deviation. (For a full Bayes estimator, Eq. (9) can still be used to approximate the p-value due to the asymptotic normality of the posterior distribution.)

The bounds 1 h^l and 1 + h^u depend on the APM and the interests of the practitioner. To fix ideas, assume that data are available on yi and xi = (xi1, ... , xip) for a sample of n locations (i = 1, ... , n). In the NB APM context, an appropriate discrepancy measure between the two links g(l, h) and g(l, 1) can be defined as the maximal change in the predicted accidents when switching from the canonical link to the generalized link

$$\Delta(\theta) = \max_{i=1,\dots,n} |\mu(\etai(\beta\theta), \theta) - \mu(\eta_i(\beta_1), 1)|. \tag{10}$$

The next step is to determine a tolerance bound D⁰ such that testing

$$H_0: \Delta(\theta) > \Delta_0$$ vs. $H_1: \Delta(\theta) \leq \Delta_0$ , (11)

is equivalent to (8). Since the criterion D(h) is convex with a minimum at h = 1, the bounds 1 h^l and 1 + h^u are uniquely determined by solving h = D¹ (D0).

For a given data set, Czado and Munk (2000) indicated that the two-dimensional surface p(...,hl, hu) is reduced to the one-dimensional curve p(h.,., D0), where the level sets a = p give the asymptotic minimal bound D⁰ for which H⁰ can be rejected at level a as well as the maximal bound for which H¹ can be rejected by the discrimination test at level 1 a.

The performance of M¹ vs. M^h can be also assessed by way of simulation. Such an assessment will have to be based on criteria similar to (10), e.g., mean absolute deviation (difference) or mean squared deviation, since the regression parameters b¹ and bh are on different scales and cannot be directly compared.

4. Methodology

The joint maximum likelihood estimation of the parameters h, b and j of the NB APM based on the generalized link g(l, h) is outlined in Czado and Munk (2000). The Empirical Bayes (EB) estimates may then be used in subsequent applications of APMs. Alternatively, the full Bayes estimates can be obtained using the MCMC method. For a discussion of ''EB versus full Bayes", the reader is referred to (Miaou and Lord, 2003; Mitra and Washington, 2007).

Subsection 4.1 describes the full Bayes approach that was used for model specification, parameter estimation and performance evaluation in this paper, whereas Subsection 4.2 is devoted to goodness of fit measures.

4.1. Full Bayes estimation

The posterior distributions needed in the full Bayes approach were obtained using MCMC sampling of the NB APM in Eq. (1) by means of the statistical software WinBUGS.

MCMC methods were used to repeatedly sample from the joint posterior distribution of the parameters h, b and j. These techniques generate sequences (chains) of random points (h(t) , b(t) , j(t) t = 1, ... , T (iterations), whose distributions converge to the target posterior distributions. A sub sample (h(t) , b(t) , j(t) ), t = 1, ... , T⁰ < T, was used to monitor convergence and was then excluded as a burn-in sample. The remaining iterations (t = T⁰ + 1, ... , T) were used for parameter estimation, performance evaluation and inference.

Monitoring convergence is important because it ensures that the posterior distribution has been ''found" and thereby indicates when sampling of parameters should begin. A common methodology to check convergence is by tracking two or more parallel chains with diverse starting values to ensure full coverage of the sample space of the parameters. Convergence of multiple chains was assessed using the Gelman–Rubin convergence statistic as modified by Brooks and Gelman (1998). For a given parameter monitored over J chains, the Brooks–Gelman–Rubin (BGR) statistic is

$$BGR = I_P/I_U, (12)$$

where IP is the length of the 100(1 a)% interval for the parameter based on the pooled TJ samples and IU is the mean of the J withinchain interval lengths. A value under 1.2 of the BGR statistic indicates convergence. The BGR diagnostics are available in WinBUGS.

To obtain the full Bayes estimates the user is required to specify prior distributions for the parameters. In the absence of sufficient prior knowledge of the distributions for individual parameters on the risk factors, uninformative proper prior distributions are usually specified (Miaou and Lord, 2003; Mitra and Washington, 2007).

Under M1, it is usually assumed that j gamma (0.001, 0.001) and b^j N(0, 1002 ), j = 0, 1, ... , p. Such priors represent proper distributions that are flat; hence the term uninformative (also called diffused) proper priors.

Under Mh, the same uninformative prior was used for j. A truncated normal distribution is then used as a proper prior for h, namely, h N(1, r² ^h ), 0 < h < 1, as it seems natural to have the prior distribution of the link parameter centered at (or symmetric about) the canonical log link (Czado and Raftery, 2006). It can be verified that, the symmetry of these prior choices requires small standard deviations (rh 0.3, or smaller). The authors considered rh = 1, 2, 4 and noted that link values (h) larger than 1 are favored when these parameter values are used. Alternatively, Pareto priors with high density at small values of h and monotonically decreasing tails were also considered, as the generalized link function in Eq. (4) changes most rapidly for small values of h.

The following development is along the lines of Czado and Raf tery (2006), who advocated a reference proper prior distribution for bh. These priors vary with h and assume little prior information. Let zhⁱ = ln(lhi) + (yi lhi)/lhⁱ denote the adjusted dependent variable, xhⁱ = j lhi/(j + lhi) denote its weight and define the weighted summary statistics:

$$\bar{z}{\theta} = \sum{i=1}^{n} \omega{\theta i} z{\theta i} / \sum{i=1}^{n} \omega{\theta i}, s{\theta 0}^{2} = \sum{i=1}^{n} \omega{\theta i} (z{\theta i} - \bar{z}{\theta})^{2} / \sum{i=1}^{n} \omega{\theta i}, \quad \bar{x}{\theta j} = \sum{i=1}^{n} \omega{\theta i} x{i j} / \sum{i=1}^{n} \omega{\theta i}, s{\theta j}^{2} = \sum{i=1}^{n} \omega{\theta i} (x{i j} - \bar{x}{\theta j})^{2} / \sum{i=1}^{n} \omega{\theta i}, \quad j = 1, \dots, p.$$ (13)

For the (p + 1)-dimensional regression coefficient vector, including an intercept, it is assumed that bh Np+1(th, Rh), where the mean vector is given by th ¼ ðzh; 0; - ; 0Þ ⁰ and the covariance matrix is given by

$$\Sigma{\theta} = s{\theta 0}^{2} \begin{pmatrix} 1 + \sum{j=1}^{p} \bar{\chi}{\theta j}^{2} / s{\theta j}^{2} & -\bar{\chi}{\theta 1} / s{\theta 1}^{2} & \cdots & -\bar{\chi}{\theta p} / s{\theta p}^{2} \ -\bar{\chi}{\theta 1} / s{\theta 1}^{2} & 1 / s{\theta 1}^{2} & \cdots & 0 \ \cdots & \cdots & \cdots & \cdots \ -\bar{\chi}{\theta p} / s{\theta p}^{2} & 0 & \cdots & 1 / s_{\theta p}^{2} \end{pmatrix}.$$ (14)

4.2. Goodness of fit measures

Let D denote the un-standardized deviance of the negative binomial model. Spiegelhalter et al. (2002) proposed the deviance information criteria (DIC) as a measure of model complexity and fit

$$DIC = \bar{D} + p_D; \quad p_D = \bar{D} - \hat{D}, \tag{15}$$

where D is the posterior mean of D, D^ is the point estimate obtained by substituting the posterior means of (li, j) in D and pD is a measure of model complexity estimating the effective number of parameters. As a goodness of fit measure, DIC is a Bayesian generalization of Akaike's information criteria (AIC) that is penalizing larger parameter models. According to the WinBUGS' DIC guide (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 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.

The Bayes factor Bh¹ was also employed to assess the improvement in fit over a NB APM with a canonical link (M1) when a generalized link (Mh) is used. Following Czado and Raftery (2006), Bh¹ is approximated by

$$2 \ln(B_{\theta 1}) \approx D(M1) - D(M{\theta}) + (E_{\theta} - E_1),$$ (16)

where D(M1) and D(Mh) are obtained by substituting the corresponding posterior means of (li, j) in D. The computation of E^h in Eq. (16) requires the maximum likelihood estimates of the parameters and their covariance matrix in addition to the prior means and covariance matrices specified for j and bh. E¹ is the corresponding function under the canonical link. The details are given in (Cza do and Raftery, 2006). According to Raftery (1996), M¹ is supported if 2 ln(Bh1) is less than 0, otherwise M^h is supported; with the range 0–2.2 indicating weak support, 2.2–6 representing positive evidence, 6–10 demonstrating strong evidence and higher values signifying very strong evidence.

Other measures such as the model deviance G² and R² have been suggested (Mitra and Washington, 2007). Yet, the G² and R² measures proposed in their Eqs. (14) and (15) use a Poisson-distributed error structure where the mean is equal to the variance and hence need to be modified to account for the extra-variation (Wood, 2002). This was not pursued for the following reasons. According to the WinBUGS' DIC guide (Spiegelhalter et al., 2005), the posterior mean of the deviance (Dbar = G² ) should perhaps not really be considered as a measure of fit. The fit might better be measured by deviance evaluated at the parameter of interest (Dhat). As emphasized by Linde (2005), the posterior mean of the deviance (Dbar) is more a measure of model 'adequacy'. This is evident since the addition of a covariate might in certain cases increase the Dbar value. On the other hand, the R² measure of fit has been criticized in the literature and it has been shown to provide misleading results (Miaou et al., 1996).

5. Data and model description

A total of 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 Intersection Safety Camera (Sayed and de Leur, 2007). The data on accident frequencies and traffic volumes within the intersections were obtained from the city of Edmonton and covered the period from 1997 to 2003.

Average values for the traffic volumes (over the 7 years period) are used to build APMs for predicting the total (aggregate) number of accidents using

$$\eta = \ln(\alpha_0) + \alpha1 \ln(V{mj}) + \alpha2 \ln(V{mi}).$$ (17)

The aggregation is justified on several grounds. For instance, an aggregate APM was found to perform well compared with APMs developed to handle temporal correlation (Lord and Persaud, 2000). Also, the aggregation of accidents over a period of reasonable length helps to avoid confounding effects and such a phenomenon like regression-to-the-mean (Cheng and Washington, 2005). A statistical summary of the data is shown in Table 1.

The AADT functional form (17) is popular among practitioners and is often preferred over models that include several covariates because it can be easily re-calibrated (Lord et al., 2008; and the references therein). Further, it will be adopted for estimating safety

Table 1 Statistical summary of data set (99 intersections).

| | Minimum | Maximum | Mean | Standard
deviation | | |---------------------------------|---------|---------|--------|-----------------------|--| | AADT on minor
approach (Vmi) | 2809 | 29,128 | 16,043 | 6564 | | | AADT on major
approach (Vmj) | 5642 | 61,744 | 28,010 | 9705 | | | Total number of
accidents | 0 | 193 | 70 | 39 | |

Table 2 Model statistics.

| Canonical link | Estimate | Standard
deviation | | 95% Credible intervals | | | | |-----------------------|----------|-----------------------|----------------|------------------------|--|--|--| | (M1) | | | Lower
limit | Upper
limit | | | | | Intercept | 8.96 | 1.20 | 11.28 | 6.64 | | | | | AADT major | 0.88 | 0.13 | 0.64 | 1.12 | | | | | AADT minor | 0.44 | 0.09 | 0.25 | 0.63 | | | | | j | 7.02 | 1.18 | 4.95 | 9.58 | | | | | DIC | 924.87 | | | | | | | | Generalized link (Mh) | | | | | | | | | Intercept | 3.35 | 0.60 | 4.41 | 2.25 | | | | | AADT major | 0.34 | 0.06 | 0.22 | 0.46 | | | | | AADT minor | 0.29 | 0.04 | 0.18 | 0.34 | | | | | j | 6.76 | 1.16 | 4.78 | 9.28 | | | | | h | 1.48 | 0.14 | 1.25 | 1.81 | | | | | DIC | 926.96 | | | | | | | | | | | | | | | | | Bayes factor (Bh1) | | | | | | | | | 2ln (Bh1) | 15.11 | | | | | | |

performance in the forthcoming highway safety manual (HSM) (Hughes et al., 2005). It should be noted that AADT models may suffer from omitted variables bias since many non-flow related factors are known to affect accident frequency. Nevertheless, such models are well recognized in the traffic safety literature.

6. Results

This section compares the results of the canonical log link and those of the generalized link function. The results were evaluated in terms of model prediction and goodness of fit.

The posterior summaries in Table 2 were obtained using J = 2 chains, T = 10,000 iterations, T⁰ = 5000 burn-in sample in Win-BUGS. The BGR statistics in Eq. (12), for the parameters h, b and j were less than 1.2, indicating convergence.

Table 2 summarizes the model parameter estimates and their associated statistics using the canonical and generalized links, respectively. The table shows that the parameter estimates are significant as the 95% credible intervals are bounded away from zero. The coefficients of the major and minor flows are all positive, indicating an increase in the mean accident frequency with traffic flow. The shape parameter j of the negative binomial model is significant, demonstrating the presence of over-dispersion in the data. The credible interval for the parameter h does not straddle 1, signifying the need for a generalized link function.

The DIC statistics (Eq. (15)) were 924.87 and 926.96 under M¹ and Mh, respectively, giving a slight advantage to M¹ over Mh. In contrast, the Bayes factor calculation (Eq. (16)) produced

$$2\ln(B_{\theta 1}) = 920.7 - 923.7 + (-39.64 + 57.75) = 15.11,$$

providing very strong support for Mh.

In order to validate the canonical model M¹ against Mh, the maximal change discrepancy measure in Eq. (10) is adopted. PROC GENMOD in SAS was used, with the inverse and forward link functions in Eqs. (4) and (5), to compute D(h), h 2 [0.5, 2.5]. The maximal change and p-value curves are portrayed in Figs. 1 and 2, respectively.

The p-curve in Fig. 2 shows that using a generalized link will result in a maximal absolute difference of 2.4 in predicted accidents at a = 0.1 compared to a canonical link. A canonical link can be validated at a = 0.1, if one is willing to tolerate a change of 6.3 predicted accidents.

Fig. 1. Maximal change in predicted accidents when switching from M¹ to Mh.

Fig. 2. p-value curve for assessing the maximal change in predicted accidents.

7. Discussion

For the current data set, convergence was indicated with 10,000 iterations, of which 5000 were discarded as a burn-in sample. For larger data sets and/or the addition of covariates, a larger number of iterations might be needed to achieve convergence. The user might also consider alternative priors, e.g., gamma (0.1, 0.1) for j, N(0, 102 ), for bj, a truncated normal with mean 1 and larger standard deviation or a Pareto prior for h.

The two models M^h and M1, may predict very different values depending on the traffic flow data as depicted in Fig. 3, which displays the difference between the two means lh and l¹ at various combinations of traffic flows (10,000–60,000 for the major approach and 5000–30,000 for the minor approach). The response

Fig. 3. Response surface of differences in predicted accidents (canonical – generalized link functions) by major and minor flows.

surface in Fig. 3 indicates that the difference l¹ lh increases considerably with the major flow. In contrast, it decreases as the minor flow rises (with few exceptions in the right upper corner of the graph). Consequently, the difference l¹ lh is positive and large whenever there is a heavy traffic at the major approach, especially if combined with light traffic at the minor approach. On the other hand, the difference l¹ lh is negative and smallest whenever the opposite holds true.

For the current data set, the results of goodness of fit as measured by DIC seem to be similar under M^h and M1. However, as pointed out by O'Hara (2007), comparisons of DIC from different models may lead to similar results unless DIC is focused on the right part of the model, which is not a straightforward task. Alternatively, the Bayes factor can be used to assess the improvement in fit over M¹ when M^h is used. In this case the Bayes factor lends strong support to the generalized link model Mh.

On the other hand, link validation has indicated that an absolute difference of 2 predicted accidents is sufficient to discriminate between the two models, but the canonical link can be validated only if an absolute difference of 7 predicted accidents is tolerated. Whether or not a change of 7 predicted accidents is acceptable is entirely up to the practitioner (e.g., road authority). If a change of this magnitude is allowed, the canonical link can be used; otherwise the generalized link should be adopted.

8. Summary and conclusions

The generalized link was employed in association with the traditional negative binomial model to predict the number of accidents at 99 intersections using traffic flows at the major and minor approaches as covariates.

Using DIC as a goodness of fit measure, the NB APM adopting the generalized link was found to fit the data as well as the corresponding model based on the usual canonical log link. In contrast, the Bayes factor provided strong support to the generalized link over the canonical link.

A procedure for link validation was also described. It allows the practitioners (road authorities) to consider the changes in predicted accidents that will result if a generalized link is used instead of a canonical link. If a certain maximal change is tolerated, the canonical link can be used to analyze the data; otherwise the generalized link is worth the extra efforts and should be adopted.

The two models differ in terms of prediction. Compared with the NB APM based on the canonical log link, the model associated with the generalized link seems to produce lower predictions whenever the major flow is heavy, especially if combined with light flow on the minor approach. It also seems to predict a higher number of accidents when the opposite holds true. The canonical and generalized links may also lead to large differences in the identification and ranking of accident prone locations.

It should be emphasized that erroneous selection of one model over the other can lead to serious consequences. Thereby, it is crucial to assess the goodness of fit of the two models and to validate the generalized link as outlined above. Also, while fitting the generalized link may take longer (in terms of convergence), it is not markedly complicated (in terms of WinBUGS coding) than the canonical link. Finally, whether a specific accident difference (e.g., 7 or more or less) is worth the extra effort is completely up to the user.

It should be noted that these findings maybe be valid to this data set only. In order to generalize the comparison of the two models, it is perhaps necessary to involve other datasets and to carry on a large simulation study considering different data scenarios (different roadway elements, sample sizes, period of time, etc.) for locations with known safety status.

This study can be extended in many directions. For example, a log-normal prior could be used instead of the gamma prior leading to generalized links associated with Poisson log-normal models. Also, a generalized link could be developed to connect the dispersion parameter j to relevant covariates. Other extensions include applying the generalized link approach to other data sets, in order to evaluate its performance over a wide class of traffic environments, and conducting a simulation study as suggested above.

References

Brooks, S.P., Gelman, A., 1998. Alternative methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics 7, 434– 455.
Cheng, W., Washington, S.P., 2005. Experimental evaluation of hotspot identification methods. Accident Analysis and Prevention 37, 870–881.
Czado, C., Munk, A., 2000. Noncanonical links in generalized linear models when is the effort justified? Journal of Statistical Planning and Inference 87, 317–345.
Czado, C., Raftery, A.E., 2006. Choosing the link function and accounting for link uncertainty in generalized linear models using Bayes factors. Statistical Papers 47, 419–442.
Czado, C., 1992. On link selection in generalized linear models. In: Fahrmein, L., Francis, B., Gilchrist, R., Tutz, G. (Eds.), Advances in GLIM and Statistical Modeling, Lecture Notes in Statistics, vol. 78. Springer, New York, pp. 60–65.
Czado, C., 1997. On selecting parametric link transformation families in generalized linear models. Journal of Statistical Planning and Inference 61, 125–139.
El-Basyouny, K., Sayed, T., 2006. Comparison of two negative binomial regression techniques in developing accident prediction models. TRR: Journal of the TRB. 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.
Horowitz, J.L., Mammen, E., 2007. Rate-optimal estimation for a general class of nonparametric regression models with unknown link functions. Annals of Statistics 35 (6), 2589–2619.
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)). Washington, DC. http://onlinepubs.trb.org/onlinepubs/nchrp/nchrp\_w62.pdf.
Kim, H., Sun, D., Tsutakawa, R.K., 2002. Lognormal vs. gamma: extra variations. Biometrical Journal 44 (3), 305–323.
Linde, A., 2005. DIC in variable selection. Statistica Neerlandica 59 (1), 45–56.
Lord, D., Persaud, B., 2000. Accident prediction models with and without trend: application of the generalized estimating equation. Transportation Research Record 1717, 102–108.
Lord, D., 2006. Modeling motor vehicle crashes using Poisson–gamma models: examining the effects of low sample mean values and small sample size on the estimation of the fixed dispersion parameter. Accident Analysis and Prevention 38 (4), 751–766.
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., 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.
Maher, M.J., Summersgill, I., 1996. A comprehensive methodology for the fitting of predictive accident models. Accident Analysis and Prevention 28 (3), 281–296.
Mallick, B., Gelfand, A.E., 1994. Generalized linear models with unknown link functions. Biometrika 81 (2), 237–245. 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.
McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models, second ed. Chapman & Hall, London.
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.
Miaou, S.P., Lu, A., Lum, H.S., 1996. Pitfalls of using R-squared to evaluate goodness of fit of accident prediction models. Transportation Research Record 1542, 6–
Mitra, S., Washington, S., 2007. On the nature of over-dispersion in motor vehicle crash prediction models. Accident Analysis and Prevention 39, 459–468.
Ntzoufras, I., Dellaportas, P., Forster, J., 2003. Bayesian variable and link determination for generalised linear models. Journal of Statistical Planning and Inference 111, 165–180.
O'Hara, B., 2007. Focus on DIC. <http://www.deepthoughtsandsilliness.blogspot. com/2007/12/focus-on-dic.html>, (accessed May 2008).
Pregibon, D., 1980. Goodness of link tests for generalized models. Journal of the Royal Statistical Society C 29, 5–24.
Raftery, A.E., 1996. Approximate Bayes factors and accounting for model uncertainty in generalized linear models. Biometrika 83, 251–266.
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. In: Proceedings of the 2007 Annual meeting of the Transportation Research Board, Washington, DC.
Spiegelhalter, D., Thomas, A., Best, N., Lunn, D., 2005. WinBUGS user manual. MRC Biostatistics Unit, Cambridge. http://www.mrc-cam.ac.uk/bugs.
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.
Wood, G.R., 2002. Generalised linear accident models and goodness of fit testing. Accident Analysis and Prevention 34 (4), 417–427.

elbasyounyapplicationgeneralizedlinkfunctions2010.md 46 KB Riwayat Mentahan