Measuring direct and indirect treatment effects using safety performance intervention functions

Safety Science 50 (2012) 1125–1132

Contents lists available at SciVerse ScienceDirect

Safety Science journal homepage: www.elsevier.com/locate/ssci

Measuring direct and indirect treatment effects using safety performance intervention functions Karim El-Basyouny a,⇑, Tarek Sayed b,1 a b

City of Edmonton Assistant Professor of Urban Traffic Safety, Department of Civil and Environmental Engineering, University of Alberta, Edmonton, AB, Canada T6G 2W2 Department of Civil Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z4

a r t i c l e

i n f o

Article history: Received 2 February 2011 Received in revised form 24 June 2011 Accepted 7 November 2011 Available online 27 December 2011 Keywords: Intervention models Full Bayes estimation Poisson–lognormal regression Random parameters models Treatment effectiveness index

a b s t r a c t Recent traffic safety research has advocated the use of the Full Bayes (FB) approach to conduct controlled before–after safety evaluations. The FB approach was shown to offer both methodological and data advantages. However, there is still a lack of complete understanding of the role that the model parameters play in the formulation of treatment effectiveness. This paper offers a novel approach to compute components related to direct and indirect treatment effects under a linear intervention model for the Safety Performance Function (SPF). The isolation of a component corresponding to the direct treatment effects enables the analyst to assess the effectiveness of the countermeasures apart from local (siterelated) environmental factors. Two case studies are used to demonstrate the approach. The first case study is based on the results published from a previous safety evaluation in Iowa and the second uses a new data set from British Columbia based on the Insurance Corporation of British Columbia Road Improvement Program. The various direct and indirect components have yielded valuable insight into the effects of the model parameters on treatment impacts. Moreover, the calculation of the treatment effectiveness indices have been simplified by providing straightforward equations in terms of model parameters for the computation of their components without resorting to additional algorithms. More importantly, the results can have a significant impact on the economic evaluation of safety programs and countermeasures by allowing the calculation of collision modification factors (CMFs) that vary with time. Ó 2011 Published by Elsevier Ltd.

1. Introduction There are several methods to estimate the safety effect of changes in roadway characteristics such as Randomized Controlled Trials (RCTs), quasi-RCTs, and Controlled Before–After (CBA) studies. Although CBA studies may not be the best way for evaluating road safety interventions, they are effective and very common in road safety applications (e.g., Aeron-Thomas et al., 2005; Wilson et al., 2006; and the references therein). A CBA study involves a design where there is contemporaneous data collection before-and-after the intervention and an appropriate control site or activity. The safety evaluations of CBA studies involve a time series intervention analysis, as there is a clearly defined point in time when the intervention occurred (countermeasures implemented). These longitudinal evaluations base their results on actual changes that have occurred over a period of time extending from the before condition to the after condition. The ⇑ Corresponding author. Tel.: +1 780 492 9564 E-mail addresses: [email protected] (K. El-Basyouny), tsayed@ civil.ubc.ca (T. Sayed). 1 Tel.: +1 604 822 4379. 0925-7535/$ - see front matter Ó 2011 Published by Elsevier Ltd. doi:10.1016/j.ssci.2011.11.008

outcome of these evaluation studies is used to quantify the economic benefits of road safety improvement programs and to produce collision modification factors for estimating the safety implications of traffic operations and highway design decisions. It is therefore imperative that the evaluation methodology is rigorous such that the results are robust and can withstand technical scrutiny (Sayed et al., 2004). Traditionally, the Empirical Bayes (EB) methodology is used to determine the effectiveness of safety treatments (Hauer, 1997; Hauer et al., 2002; Sayed et al., 2004; Persaud and Lyon, 2007). A typical EB CBA study requires the collection of three distinct sets of data: (i) treatment sites, (ii) comparison sites, and (iii) reference sites. The comparison group is used to correct for time trend effects and includes sites that have not been treated but they experience similar traffic and environmental conditions. Therefore, the geographic proximity and comparability to the treatment sites should be the main criterion for comparison sites selection. The treatment effectiveness index is used to capture the change of collisions in the comparison group to the change of collisions in the treatment group. The reference group is used to correct for the regressionto-the-mean artifact. Usually, the reference group includes a larger number of sites that are similar to the treatment sites and is used

1126

K. El-Basyouny, T. Sayed / Safety Science 50 (2012) 1125–1132

to develop a Safety Performance Function (SPF). The EB approach is used to refine the estimate of the expected number of collisions at a location by combining the observed number of collisions (at the location) with the predicted number of collisions from the SPF. On the other hand, Full Bayes (FB) estimation is a wellestablished statistical methodology, and there is a much greater statistical literature available to guide its correct use in practice. Yet, its application in CBA studies has been proposed in the road safety literature only recently (Carriquiry and Pawlovich, 2005; Pawlovich et al., 2006; Aul and Davis, 2006; Li et al., 2008; Park et al., 2009; Lan et al., 2009; El-Basyouny and Sayed, 2010, 2011). A number of studies have compared the use of EB and FB approaches in traffic safety applications (Huang et al., 2009; Lan et al., 2009; Persaud et al., 2009). The results of such comparisons have outlined the several advantages of the FB technique, which can be categorized into methodological and data advantages. In terms of methodological advantages, the FB approach is appealing for many reasons that include having the ability to account for all uncertainty in the data, to provide more detailed inference, and to allow inference at more than one level for hierarchical models, among others. In terms of data requirements, the FB approach efficiently integrates the estimation of the SPF and treatment effects in a single step, whereas these are separate tasks in the EB method. To benefit from the additional advantages of the FB approach, several researchers have proposed the use of intervention models in the context of a CBA safety evaluation (Pawlovich et al., 2006; Li et al., 2008). For example, Li et al. (2008) considered various intervention forms and hierarchical models (Poisson–Gamma or Poisson–Lognormal) to conduct a safety evaluation. The forms/ models were developed to deal with both immediate and gradual treatment impacts while accounting for countermeasure implementation, time effects, traffic volumes as well as the effects of other covariates representing various site characteristics. Park et al. (2009) extended the intervention model to include multivariate dependent variables with multiple regression links and proposed an algorithm for the computation of the treatment effectiveness index to determine the efficacy of the countermeasure. To gain new insights into how collision counts are influenced by covariates and to account for heterogeneity, the use of multivariate intervention models with random parameters among matched treatment-comparison pairs were advocated by El-Basyouny and Sayed (2011). It should be noted that in all of the above studies, linear slopes were assumed to represent the time and treatment effects across the treated and comparison sites. Although the results of the above mentioned studies have demonstrated the superiority of the FB CBA safety approach, there is still a lack of complete understanding of the role that the model parameters play in the estimation of the treatment effectiveness. This is partially due to the perceived notion that the FB methodology is complex and that a particular algorithm (as proposed by Park et al., 2009) is required to compute the treatment effectiveness index. This paper offers a novel approach to compute the various components of the index under a linear intervention model for the SPF. The various components are related to direct and indirect treatment effects. The direct effects are further decomposed into long term trend changes and overall mean level changes whereas the indirect effects are imposed on the predicted collisions through traffic volumes and other site characteristics that vary with time (e.g., weather conditions). The importance of isolating a component corresponding to the direct treatment effects cannot be over emphasized as it enables analysts to assess the effectiveness of the countermeasures apart from local (site-related) environmental factors. These specific aspects of treatment effectiveness are important for traffic safety applications and help overcome the imperfection of the rather ‘‘blind’’ treatment effectiveness index.

The objectives of this study are threefold: (i) to provide a better understanding of the role of model parameters in estimating treatment effectiveness, (ii) to simplify the calculation of the treatment effectiveness indices by providing straightforward equations in terms of model parameters for the computation of their components without resorting to additional algorithms, and (iii) to provide tools for investigating the safety changes in treated sites due to the implementation of countermeasures. The approach is demonstrated using two case studies. The first case study is based on the results published from the Iowa study conducted by Li et al. (2008). The second case study uses a new data set from British Columbia based on the Insurance Corporation of British Columbia (ICBC) 2009 Road Improvement Program (RIP). It is hoped that the approach described in this paper will help advance the use of the FB CBA safety analysis in order to benefit from the many advantages offered by the technique. 2. Methodology 2.1. The models Consider a CBA study where it is assumed that the data are available for a reasonable period of time before and after the intervention. Let Yit denote the collision count observed at site i (i = 1, 2, . . . , n) during year t (t = 1, 2, . . . , m). To introduce the linear intervention models, let Ti denote the treatment indicator (equals 1 for treated sites, zero for comparison sites), t0i denote the intervention year for the ith treated site and its matching comparison group, Iit denote the time indicator (equals 1 in the after period, zero in the before period) and V1it and V2it denote the annual average daily traffic (AADT) at the major and minor approaches, respectively. Further, let (X3i, ... , XJi) denote additional covariates representing geometric and environmental site characteristics. 2.1.1. The Poisson–Lognormal Intervention (PLNI) models For the PLNI model (M1a), it is assumed that the Yit are independently distributed as

Y it jhit  Poissonðhit Þ;

ð1Þ

lnðhit Þ ¼ lnðlit Þ þ ei ;

ð2Þ

lnðlit Þ ¼ a0 þ a1 T i þ a2 t þ a3 ðt  t 0i ÞIit þ a4 T i t þ a5 T i ðt  t 0i ÞIit þ b1 lnðV 1it Þ þ b2 lnðV 2it Þ þ b3 X 3i þ    þ bJ X Ji ;

ð3aÞ

ei  Nð0; r2e Þ;

ð4Þ

where (a0, ... , bJ) are the regression coefficients and re represents the extra-Poisson variation. In Eq. (3a), the parameter a1 represents the countermeasures effect (the difference in log collision count between treated and comparison sites), a2 represents a linear time trend, a3 represents the slope due to the intervention, a4 and a5 allow for different time trends and different intervention slopes across the treated and comparison sites, whereas the remaining regression coefficients represent the effects of traffic volume and other traffic-related site characteristics on log collision counts. Since the changes at the treatment sites might not be gradual and a sudden drop (or increase) in collision counts are likely to occur at those sites which received the intervention, Li et al. (2008) proposed an additional parameter in Eq. (3a) in order to accommodate such an effect leading to the model M1b 2

lnðlit Þ ¼ a0 þ a1 T i þ a2 t þ a3 ðt  t 0i ÞIit þ a4 T i t þ a5 T i ðt  t 0i ÞIit þ a6 T i Iit þ b1 lnðV 1it Þ þ b2 lnðV 2it Þ þ b3 X 3i þ    þ bJ X Ji : ð3bÞ

K. El-Basyouny, T. Sayed / Safety Science 50 (2012) 1125–1132

2.1.2. The PLNI models with random parameters among pairs Since the matched comparison sites were selected to be as similar to treatment sites as possible, this may induce a correlation in collision count between sites within comparison-treatment pairs. To account for this correlation, suppose that the ith site belongs to the pair p(i) e {1, 2, ... , NC}, where NC denotes the number of comparison groups. The variation due to the comparison-treatment pairing can be represented by allowing the regression coefficients in Eq. (3a) to vary randomly from one pair to another. This leads to the model M2a

lnðlit Þ ¼ apðiÞ;0 þ apðiÞ;1 T i þ apðiÞ;2 t þ apðiÞ;3 ðt  t 0i ÞIit þ apðiÞ;4 T i t þ apðiÞ;5 T i ðt  t 0i ÞIit þ bpðiÞ;1 lnðV 1it Þ þ bpðiÞ;2 lnðV 2it Þ þ bpðiÞ;3 X 3i þ    þ bpðiÞ;J X Ji ;

ð5aÞ

and





apðiÞ;j  Nðaj ; r2a;j Þ; j ¼ 0; 1; . . . ; 5; bpðiÞ;j  N bj ; r2b;j ; j ¼ 1; . . . ; J:

ð6Þ

It should be noted that several alternative distributions were considered in (6) but the normal distribution was found to provide the best statistical fit (Li et al., 2008; Milton et al., 2008; Anastasopoulos and Mannering, 2009). In practice, the random parameter ^ a;j ap(i),j (bp(i),j) is used in (5a) whenever the posterior estimate r ^ b;j ) is significantly greater than 0, otherwise the parameter aj (r (bj) is fixed across pairs (Anastasopoulos and Mannering, 2009). The random parameters model (5a) enables meaningful inference both at the pair-level and the overall-level. The PLNI model with random parameters among pairs and a jump (denoted M2b) is given by

lnðlit Þ ¼ apðiÞ;0 þ apðiÞ;1 T i þ apðiÞ;2 t þ apðiÞ;3 ðt  t 0i ÞIit þ apðiÞ;4 T i t þ apðiÞ;5 T i ðt  t 0i ÞIit þ apðiÞ;6 T i Iit þ bpðiÞ;1 lnðV 1it Þ þ bpðiÞ;2 lnðV 2it Þ þ bpðiÞ;3 X 3i þ ::: þ bpðiÞ;J X Ji ;

misleading just to report the model with the lowest DIC. Posterior predictive checking can be used for assessing the goodness-of-fit for various competing models (Gelman et al., 2004; Li et al., 2008). 3. Measuring treatment effectiveness Let lTBi and lTAi denote the predicted collision counts for the ith treated site averaged over appropriate years during the before and after periods, respectively, and let lCBi and lCAi denote the corresponding quantities for the matching comparison group where the predicted collision counts are averaged over appropriate sites (all sites in the matching comparison group) and years. Following Park et al. (2009), the ratio lCAi/lCBi can be used to adjust the prediction for general trends between the before and after periods at the ith treated site. Thus, the predicted crashes in the after period for the ith treated site had the countermeasures not been applied is given by pTAi = lTBi(lCAi/lCBi). The index of effectiveness of the countermeasures at the ith treated site is given by the ratio lTAi/pTAi, which reduces to hi ¼ lTAi lCBi =lTBi lCAi , or

lnðhi Þ ¼ lnðlTAi Þ þ lnðlCBi Þ  lnðlTBi Þ  lnðlCAi Þ:

lnðh Þ ¼ ð1=nÞ

n X

  ln hi :

The overall treatment effect is calculated from (h⁄  1), while the overall percentage of reduction in predicted collision counts is given by (1  h⁄)  100. In the same way, let LTAi, LTBi, LCAi and LCBi denote the averages of ln (lit) corresponding to lTAi, lTBi, lCAi and lCBi, respectively. Hence, by averaging the predicted collisions on the log-scale, a corresponding set of indices can be obtained from

ð9Þ

and

lnðhÞ ¼ ð1=nÞ

To obtain the full Bayes estimates of the unknown parameters, it is required to specify prior distributions for the hyper parameters. Such prior distributions were recently discussed in the safety literature (Li et al., 2008; Ma et al., 2008; Aguero-Valverde and Jovanis, 2009; Park et al., 2009; El-Basyouny and Sayed, 2010, 2011). The most commonly used priors are diffused normal distributions (with zero mean and large variance) for the regression parameters and Gamma(e, e) or Gamma(1, e) for the precision (inverse variance) parameters, where e is a small number (e.g., 0.01 or 0.001). In the current study, the priors N(0, 1002) and Gamma(0.001, 0.001) were used for the regression and precision parameters, respectively. The posterior distributions needed in the full Bayes approach were sampled using the Markov Chain Monte Carlo (MCMC) techniques available in WinBUGS (Lunn et al., 2000). The BGR statistics (Brooks and Gelman, 1998), ratios of the Monte Carlo errors relative to the standard deviations of the estimates and trace plots for all model parameters were monitored for convergence. The Deviance Information Criteria (DIC) was used for model comparisons (Spiegelhalter et al., 2002). As a goodness of fit measure, the DIC is a Bayesian generalization of Akaike’s information criterion (AIC) that penalizes larger parameter models. According to the DIC guidelines (Spiegelhalter et al., 2005), 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

ð8Þ

i¼1

where apðiÞ;6  Nða6 ; r2a;6 Þ. 2.2. Full Bayes estimation

ð7Þ

The overall index can be computed from

lnðhi Þ ¼ LTAi þ LCBi  LTBi  LCAi : ð5bÞ

1127

n X

lnðhi Þ:

ð10Þ

i¼1

Using Eqs. (9) and (10), the overall treatment effect is calculated from (h  1), while the overall percentage of reduction in predicted collision counts is given by (1  h)  100. Although the logarithmic function ln(x) is not linear over the range 0 < x < 1, it is approximately linear for x > 2 and the linear approximation is indeed excellent for larger values of x. Thus, it is conjectured that h would be fairly close to h⁄ in most of the CBA traffic safety studies carried out in practice unless the observed collision counts are fairly small. To gain some insight into the effects of the regression parameters on the index of treatment effectiveness, consider the model M1b in (3b) and note that the regression parameters (a0, a1, a2, a3, b3, ... , bJ) cancel out of (9). Further, for the ith treated site and t > t0i, let s(i) = t  t0i denote the number of post-treatment years. To simplify the notation, the subscript i is dropped and the notation s = 1, 2, ... is used instead, but it is important to keep in mind that s vary with i. Let his denote the treatment effectiveness index for the ith treated site after s post-intervention years. It is easily verified (see Appendix A) that

his ¼ C 1 ði; sÞC 2 ði; sÞC 3 ði; sÞC 4 ði; sÞ;

ð11Þ

C 1 ði; sÞ ¼ expf0:5a4 ðs þ t0i þ 1Þg;

ð12aÞ

C 2 ði; sÞ ¼ expfa6 þ 0:5a5 ðs þ 1Þg;

ð12bÞ

C 3 ði; sÞ ¼ ½hðV 1is Þb1 ;

ð12cÞ

1128

K. El-Basyouny, T. Sayed / Safety Science 50 (2012) 1125–1132

C 4 ði; sÞ ¼ ½hðV 2is Þb2 ;

ð12dÞ

where h(V1is) and h(V2is) are indices defined in terms of the major and minor traffic volumes the same way his is defined in terms of the predicted collision counts. Eqs. (12a)–(12d) hold also for the model M1a in (3a) and for the random parameters PLNI models M2a and M2b at the overall-level of inference. Yet, the parameter a6, which appears in C2(i, s), is set to zero for the models without a jump. It should also be noted that the smaller are these components, the larger are the reductions in predicted collisions and the more effective are the countermeasures (treatments). For the models M1a and M2a, the two components C1(i, s) and C2(i, s) involve the difference in time trends and intervention slopes between treated and comparison sites. If the time trend for the ith treated site is smaller than that for the matching comparison group, then a4 would be negative and the index will be reduced. In contrast, the index would be increased for positive values of a4. On the other hand, if the decrease in collision counts accelerates after treatment implementation then a5 would be negative and the index will be reduced and vice versa. For the two models M1b and M2b, the component C2(i, s) involves the additional parameter a6 representing a possible post-treatment sudden drop/increase in collisions. Typically, a sudden post-treatment drop in collisions is likely to reduce the index and vice versa. The two components C3(i, s) and C4(i, s) involve the respective indices for the major and minor traffic volumes raised to the fractional powers b1 and b2, which typically satisfy the conditions 0 < b1 < 1 and 0 < b2 < 1. Thus, both components would be closer to 1 than the two traffic volume indices themselves. In fact, unless the traffic volumes are subject to significant annual fluctuations, both components are expected to be near 1. Since C1(i, s) and C2(i, s) measure different aspects of the treatment effect, they can be combined to obtain a component that is inversely related to the direct treatment impact

C T ði; sÞ ¼ expfa þ bsg; a ¼ a6 þ 0:5a5 þ 0:5a4 ðt0i þ 1Þ;

b ¼ 0:5ða5 þ a4 Þ: ð13aÞ

According to (13a), CT(i, 1) = exp {a + b} and CT(i, s) = exp {b}CT(i, s  1) for s > 1. Subsequently, this combined component would either decline or grow according to whether b < 0 or b > 0. Further, the percent annual change, [oCT(i, s)/os]/CT(i, s) = b, is constant during the period following the intervention. Similarly, the components C3(i, s) and C4(i, s) can be combined to obtain

C V ði; sÞ ¼ ½hðV 1is Þb1 ½hðV 2is Þb2 :

ð13bÞ

The combined component in (13b) is inversely related to the indirect (through traffic volumes) treatment impact under the PLNI models. 4. Iowa case study Li et al. (2008) used a dataset including 28 road segments in the State of Iowa, 14 of which were converted from four through lanes to three lanes (two through lanes and a center turning lane) sometime within the period 1982–2004. The other 14 sites in the study were selected to act as comparison or control sites. The number of collisions per month was recorded at the sites over segments of different lengths (0.2–2.5 miles) between January 1982 and December 2004. Monthly traffic volumes at each of the sites were estimated by the IA-DOT using average daily traffic (ADT). There is considerable collision information for the period preceding the intervention at all sites. On the other hand, collision information during the period following the intervention is rather limited at several of the study sites. The median number of months for the

before period is 215 (with a mean of 214.6), whereas the median number of months for the after period is 60 (with a mean of 60.4). A variety of intervention models were considered where road length and traffic volume were used as offset variables in the regression equations, which included three additional variables representing quarterly seasonal effects. The best fits were provided by hierarchical PLNI models that allow different values of each regression coefficient for each site leading to models similar to (5a) and (6) where p(i) is replaced by i. Thus, the components of the treatment effectiveness index are given by

his ¼ C 1 ði; sÞC 2 ði; sÞC V ði; sÞ;

s ¼ 1; 2; . . . :

where CV(i, s) = h(Vis) and Vis represents the traffic volume on the ith road segment during month s following the intervention. Consider a typical (median) treated site i with 215 and 60 months before and after the intervention, respectively, and denote the two hierarchical models by H1 (without a jump) and H2 ^ 4 ¼ 0:001; a ^ 5 ¼ 0:020Þ were ob(with a jump). The estimates ða ^ 4 ¼ 0:000; a ^ 5 ¼ 0:010; tained under H1, while the estimates ða a^ 6 ¼ 0:355Þ were obtained under H2. Using these estimates along with t0i = 216, the component C1(i, s) representing the treatment by ^ 4 ðs þ 217Þg, whereas time interaction was computed from expf0:5a the component C2(i, s) representing the treatment by intervention ^ 6 þ 0:5a ^ 5 ðs þ 1Þg. date interaction was computed from expfa Fig. 1 shows a plot of the treatment components C1(i, s), C2(i, s), and CT(i, s) on the y-axis versus the number of post treatment months (s) on the x-axis. For H1, where the estimate of the treatment by time interaction is positive, Fig. 1 shows that C1(i, s) has increased gradually from 1.115 (11.5% increase) one month following the intervention to 1.149 (14.9% increase) after 60 months. However, the negative estimate of the treatment by intervention date interaction caused C2(i, s) to drop from 0.980 (2% reduction) at the first post-treatment month to 0.543 (45.7% reduction) after 60 post-treatment months. It is readily seen that the decline in C2(i, s) has counterbalanced the increase in C1(i, s) resulting in a combined component ranging from 1.093 (9.3% increase) right after the intervention to 0.624 (37.6% reduction) at December 2004. For H2, the combined component reduces to C2(i, s) due to the zero estimate of the treatment by time interaction. The negative estimates of a5 and a6 caused the combined component to drop from 0.694 (30.6% reduction) to 0.517 (48.3% reduction) over the 60 post-treatment months. Fig. 2 shows a plot of the combined treatment components CT(i, s) on the y-axis versus the number of post treatment months (s) on the x-axis, for the models without (H1) and with (H2) a jump. Fig. 2 shows that the combined treatment curve for the model without a jump is located above that with a jump, albeit steeper. In fact, the percentages of annual decrease were 0.0095 and 0.0050, under H1 and H2, respectively. Hence, the large difference in the treatment components between the two hierarchical PLNI models diminishes with the evolution of time following the intervention. Li et al. (2008) used posterior predictive checks to assess the goodness-of-fit of H1 and H2 and concluded that ‘‘deciding in favor of one of the two hierarchical models in this study is more a matter of suitability to the application itself than to differences between models in term of their fit or predictive ability’’ and that ‘‘from a statistical point of view, either one of the two models appears to be equally appropriate’’. In fact, for the evolution of the treatment effects over the time periods following the intervention and using the statistical tools leading to Fig. 2, it appears that H2 (w/jump) is a more reasonable model than H1 (w/o jump), since the latter has a relatively steeper slope in order to make up for its failure to reflect the sudden drop in collision rates immediately after the intervention.

K. El-Basyouny, T. Sayed / Safety Science 50 (2012) 1125–1132

1129

Fig. 1. Monthly treatment components (y-axis) for a median site under H1 (w/o a jump) vs. number of post treatment months (x-axis), Li et al. (2008).

Fig. 2. Monthly combined treatment components (y-axis) for a median site under H1 and H2, vs. number of post treatment months (x-axis), Li et al. (2008).

^4 ; a ^5 ; a ^ 6 Þ along with the Note that the parameters’ estimates ða intervention dates (t0i) are the only information needed to compute the treatment components. However, the intervention dates for the 14 treated sites were not published in Li et al. (2008) and only the completion year was given in Pawlovich et al. (2006). Had the exact intervention dates (by month and year) been available, it would have been possible to compute the treatment components for i = 1, 2, ... , 14 and s = 1, 2, .... Furthermore, if the index is set to hi ¼ hi;276t0i , the overall index could have been computed from (10). Under H2, the reduction in collision rate over the study period was estimated as 39.1% for all sites (Li et al., 2008). 5. British Columbia case study In 1989 the Insurance Corporation of British Columbia (ICBC) started a Road Improvement Program (RIP). The program’s main

goal is to reduce the frequency and severity of collisions, thereby reducing deaths, injuries and insurance claim costs. Under this program a group of NT = 18 intersections (sites) in the Greater Vancouver area were selected for general intersection improvements within a controlled before–after experiment. The treated intersections were matched with NC = 16 comparison groups comprising 98 intersections. The comparison groups range in size from 1 to 10 intersections. Thus, there are n = 116 intersections in the sample. The countermeasures were implemented during the years 2004, 2005 or 2006. Data on collision injury counts along with the corresponding traffic volumes were collected for the years 2001–2008, i.e., m = 8. Hence, the total sample size is n ⁄ m = 928 observations. The average collision injury counts per year before and after the treatment are summarized in Table 1 for the treated sites. The highest annual collision injury frequencies were observed for Site

1130

K. El-Basyouny, T. Sayed / Safety Science 50 (2012) 1125–1132

Table 1 Annual collision injury frequency for treated sites. Average counts (per year) ID

Treatment year

Before

After

% Change

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

2004 2004 2004 2004 2005 2005 2005 2005 2005 2005 2005 2005 2005 2006 2006 2006 2006 2006

12.00 35.00 80.67 32.33 7.75 95.75 57.50 31.50 28.50 126.25 5.75 20.75 36.75 12.20 6.00 4.40 7.20 7.00

13.25 26.75 54.00 14.50 4.67 85.00 53.00 21.33 22.00 144.33 4.67 14.00 31.67 14.00 2.00 4.50 4.00 7.50

10 24 33 55 40 11 8 32 23 14 19 33 14 15 67 2 44 7

10 followed by sites 6, 3, and 7. Of the 18 treated sites, thirteen sites have experienced reductions (negative changes) in collision injury frequency ranging from 8–67%. In contrast, the remaining five sites have experienced increases ranging from 2–15%. The above preliminary analysis shows that the treatments may have an impact on the safety of these sites (intersections). A rigorous before–after analysis of the effects of different parameters on collision injury counts is conducted below. The posterior estimates (means) for the parameters of the PLNI models under consideration were obtained via two chains with 10,000 iterations 5000 of which were excluded as a burn-in sample using WinBUGS. The major and minor AADT statistics were standardized (by subtracting the mean) to speed convergence and alleviate autocorrelations for the regression parameters. The BGR statistics were below 1.2, the ratios of the Monte Carlo errors relative to the standard deviations of the estimates were below 0.05 and trace plots for all model parameters indicated convergence. The autocorrelation functions died out reasonable well, particularly for the intervention parameters. The regular PLNI models M1a and M1b were outperformed by the random parameters models M2a and M2b with significant reductions in DIC of 130 and 140, respectively. Also, the estimate ^ 2e ¼ 0:391 under M1a was reduced of the extra-Poisson variance r to 0.320 (18% reduction) under M2a. A similar reduction of 15% was obtained for M2b relative to M1b. The random parameters PLNI ‘‘jump’’ model M2b was ruled out as the reduction in DIC with respect to M2a was not significant. In addition, the 95% credible interval for the sudden jump parameter ^ 6 Þ was (0.191, 0.051), which includes zero. ða A posterior predictive approach (Gelman et al., 2004; Li et al., 2008) was used to assess the goodness-of-fit (adequacy) of the PLNI model M2a. Such procedures involve generating replicates under the postulated model and comparing the distribution of a certain discrepancy measure such as the chi-square statistic to the value of chi-square obtained using observed data. A model does not fit the data if the observed value of chi-square is far from the predictive distribution, i.e., the discrepancy cannot reasonably be explained by chance if the p-values are close to 0 or 1 (Gelman et al., 2004). Under M2a, the observed value of chi-square was 337.9 with an associated p-value of 0.123 indicating that the model seems to perform well in terms of accommodating the variation in accident frequency across pairs. Eqs. (7) and (9) gave close estimates for hi and hi. In fact, the difference was within 0.025 at the site level under all models

considered. For the overall level, the difference in the estimates of h⁄ and h using Eqs. (8) and (10) was within 0.001 under all models considered. For the present data, the posterior estimates were ^ 4 ¼ 0:031; a ^ 5 ¼ 0:058Þ under the PLNI model with random ða parameters among matching pairs (M2a). Thus, for those sites where the countermeasures were implemented in 2004 (t0i = 4), the treatment component (13a) simplifies to CT(i, s) = exp {0.1065–0.0445s} = 0.956CT(i, s  1), which is depicted in Fig. 3 along with C1(i, s) and C2(i, s). Fig. 3 shows a plot of the treatment components C1(i, s), C2(i, s), and CT(i, s) on the y-axis versus the number of post treatment years (s) on the x-axis. Fig. 3 shows that the treatment by time component C1(i, s) has decreased gradually from 0.911 in 2005 to 0.870 in 2008. Correspondingly, the treatment by intervention date component C2(i, s) has decreased from 0.944 to 0.865. As a result, the combined treatment component showed a decline from 0.860 (14% reduction) at the end of 2005 to 0.752 (24.8% reduction) at the end of 2008. Further, the percent annual decrease for the combined treatment component is 0.0445. This behavior is clearly manifested in Fig. 3 as CT(i, s) is decreasing steadily. The analogous curves for those sites where the countermeasures were implemented in 2005 and 2006 (t0i = 5, 6) are quite similar. It should be noted that the curves represented by C1(i, s), C2(i, s) and CT(i, s) are linear in s only on the log-scale (although they look almost linear in Figs. 1–3). Also, for a positive estimate of the treatment by time trend interaction and a negative estimate of the treatment by intervention date interaction, the combined treatment was diluted and its curve is located between C1(i, s) and C2(i, s) (Fig. 1). In contrast, when negative estimates were obtained for both components, the combined treatment was compounded and its curve is located below the two other curves (Fig. 3). The ranges of the indices for the major and minor traffic volumes were (0.984, 1.009) and (0.965, 1.035), respectively, and the posterior estimates of the regression coefficients associated with ^1 ¼ 0:290; b ^2 ¼ 0:440Þ. Substituting the two traffic volumes were ðb these results in (12c) and (12d) led to the ranges (0.995, 1.003) and (0.984, 1.015) for C3(i, s) and C4(i, s), respectively, which are closer to 1 as discussed above. Furthermore, the combined traffic volume component CV(i, s) ranged from 0.983 to 1.017. Under the random parameters PLNI model M2a, the estimate of the overall index was 0.787 ± 0.041 with the 95% credible interval (0.709, 0.868) indicating a significant overall reduction in predicted collision injury counts of 21.3%.

6. Summary and future research This paper presented a novel approach to compute various direct and indirect components of the treatment effects. The proposed approach yields valuable insights into the effects of the model parameters on treatment impacts. The analysis revealed that the parameters (a4, a5, a6), which are related to the intervention time and date interactions, have a direct influence on treatment effectiveness. For instance, if the time trend for the ith treated site is smaller than that for the matching comparison group, then a4 would be negative and the treatment effects will be increased and vice versa. On the other hand, if the decrease in collision counts accelerates after the implementation of the countermeasures then a5 would be negative and the treatment effects will be increased and vice versa. In contrast, the parameters (b1, b2), which are related to the effects of major and minor traffic volumes, represent the indirect treatment impact. In fact, it has been argued that unless the traffic volumes are subject to significant annual fluctuations, their effect on treatment effectiveness would be marginal.

1131

K. El-Basyouny, T. Sayed / Safety Science 50 (2012) 1125–1132

Fig. 3. Yearly treatment components (y-axis) under the random parameters PLNI model vs. number of post treatment years (x-axis), (completion year 2004).

It should be noted that the standard errors and credible intervals for the proposed components can be easily obtained within the FB-MCMC framework to account for prediction errors. The proposed approach provides also a simple and intuitive way of computing the overall treatment effectiveness index without resorting to additional algorithms. Given that the model specification is unchanged, then the computation of the overall index re^1 ; b ^2 Þ along with the ^4 ; a ^5 ; a ^6 ; b quires the regression coefficients ða intervention dates (t0i) and the two indices for the major and minor traffic volumes. Moreover, the analyst would be able to compute the different components of the index such as the difference in time trends and intervention slopes between treated and comparison sites and plot their change over time. This has been demonstrated using the information provided in the Iowa study by Li et al. (2008). The treatment effectiveness indices were computed and plotted for two hierarchical models: H1 (without a jump) and H2 (with a jump). Under H1, the combined treatment component ranged from 1.093 (9.3% increase) right after the intervention to 0.624 (37.6% reduction) at December 2004. Under H2, the combined treatment component ranged from 0.694 (30.6% reduction) to 0.517 (48.3% reduction) over the 60 post-treatment months. Furthermore, the results were substantiated using a new data set from British Columbia. Four PLNI models were considered for conducting a before–after evaluation of collision injury counts at certain intersections improved by the installation of countermeasures under the ICBC’s RIP. The results showed that the PLNI models allowing for random parameters among matching comparisontreatment pairs have significantly outperformed the regular PLNI models, while reducing the estimates of the extra-Poisson variance. Also, the random parameters PLNI model without a jump was preferred for analyzing the current data set, since the PLNI model with a sudden jump parameter did not lead to a significant reduction in DIC. The estimate of the overall index was 0.787 ± 0.041 implying a significant 21.3% reduction in predicted collision injury counts. The components associated with direct and indirect (manifested through traffic volumes) treatment effects, were computed for each year following the intervention. For instance, the treatment component has declined steadily from .860 to 0.752 during the period 2005–2008 for the sites whose treatment was completed in 2004.

It should be noted that although the proposed approach has been applied to CBA studies, it can be used equally well to analyze the results of Randomized Controlled Trials (RCTs) and quasi-RCTs. Overall, the results suggest that incorporating such design features as matched (yoked) comparison groups in SPFs leads to superior models. Moreover, the results of this paper can have a significant impact on the economic evaluation of safety programs and countermeasures. Current evaluation techniques assume constant countermeasure effectiveness or a collision modification factor that does not change with time. This seems unrealistic. Given the way future benefits are discounted to present values, the results of using a collision modification factor that changes with time can affect the cost effectiveness of safety countermeasures significantly. However, the linear intervention model furnishes only some restricted treatment profiles. As a result, a non-linear intervention model affording a larger variety of treatment profiles is being currently developed and the results will be submitted for publication in a forthcoming paper.

Appendix A Let the superscripts T and C, denote treated and comparison sites and let the superscripts A and B denote the after and before periods, respectively, and note that

ln



ln



ln



ln

















CB lCB ¼ a0 þ a2 t B þ c1 ln V CB it 1it þ c2 ln V 2it ;















lTA ¼ a0 þ a1 þ a6 þ ða2 þ a4 ÞtA þ ða3 þ a5 Þðt A  t 0i Þ þ c1 it     TA  ln V TA 1it þ c2 ln V 2it :





TB lTB ¼ a0 þ a1 þ ða2 þ a4 Þt B þ c1 ln V TB it 1it þ c2 ln V 2it ;

Hence, we have that

ln



CA lCA ¼ a0 þ a2 t A þ a3 ðt A  t 0i Þ þ c1 ln V CA it 1it þ c2 ln V 2it ;















lTA þ ln lCB  ln lTB  ln lCA it it it it A B A ¼ a4 ðt  t Þ þ a5 ðt  t0i Þ þ a6 þ c1 E1it þ c2 E2it ;

1132

K. El-Basyouny, T. Sayed / Safety Science 50 (2012) 1125–1132

where



E1it ¼ ln V TA 1it



      TB CA þ ln V CB 1it  ln V 1it  ln V 1it ;

        CB TB CA E2it ¼ ln V TA 2it þ ln V 2it  ln V 2it  ln V 2it : The average of tB, which runs from 1 to (t0i  1), is 0.5t0i, whereas the average of tA is t0i + 0.5(s + 1), s = 1, 2, .... Thus, the averages of (tA  tB) and (tA  t0i) are 0.5(t0i + s + 1) and 0.5(s + 1). Substituting the averages of (tA  tB), (tA  t0i), E1it and E2it in Eq. (9) leads to Eqs. (11)–(12). References Aeron-Thomas, A., Hess, S., 2005. Red-light cameras for the prevention of road traffic crashes. Cochrane Database of Systematic Reviews 2005, Issue 2. Art. No.: CD003862. Aguero-Valverde, J., Jovanis, P.P., 2009. Bayesian multivariate Poisson log-normal models for crash severity modeling and site ranking. Transportation Research Record 2136, 82–91. 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. Aul, N., Davis, G., 2006. Use of propensity score matching method and hybrid Bayesian method to estimate crash modification factors of signal installation. Transportation Research Record 1950, 17–23. Brooks, S.P., Gelman, A., 1998. Alternative methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics 7, 434– 455. Carriquiry, A., Pawlovich, M.D., 2005. From empirical Bayes to full Bayes: methods for analyzing traffic safety data. . El-Basyouny, K., Sayed, T., 2010. A Full Bayes Approach to before–after safety evaluation with matched comparisons: a case study of stop-sign in-fill program. In: Presented at the 89th Annual Meeting of the Transportation Research Board, in print, Transportation Research Record, Washington, D.C. El-Basyouny, K., Sayed, T., 2011. A multivariate intervention model with random parameters among matched pairs for before–after safety evaluation. Accident Analysis and Prevention 43 (1), 87–94. Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 2004. Bayesian Data Analysis, 2nd ed. Chapman & Hall, London.

Hauer, E., 1997. Observational Before–After Studies in Road Safety: Estimating the Effect of Highway and Traffic Engineering Measures on Road Safety. Elsevier Science Ltd.. Hauer, E., Harwood, D.W., Council, F.M., Griffith, M.S., 2002. The empirical Bayes method for estimating safety: a tutorial. Transportation Research Record 1784, 126–131. Huang, H., Chin, H.C., Haque, M.M., 2009. Empirical evaluation of alternative approaches in identifying crash hot spots: naive ranking, empirical Bayes, and full Bayes methods. Transportation Research Record 2103, 32–41. Lan, B., Persaud, B., Lyon, C., Bhim, R., 2009. Validation of a full Bayes methodology for observational before–after road safety studies and application to evaluation of rural signal conversions. Accidents Analysis and Prevention 41 (3), 574–580. Li, W., Carriquiry, A., Pawlovich, M., Welch, T., 2008. The choice of statistical models in road safety countermeasures effectiveness studies in Iowa. Accident Analysis and Prevention 40 (4), 1531–1542. 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. Milton, J., Shankar, V., Mannering, F., 2008. Highway accident severities and the mixed logit model: an exploratory empirical analysis. Accident Analysis and Prevention 40 (1), 260–266. Park, E.S., Park, J., Lomax, T.J., 2009. A fully Bayesian multivariate approach to before–after safety evaluation. In: Presented at the 89th Annual Meeting of the Transportation Research Board, Washington, D.C. Pawlovich, M., Li, W., Carriquiry, A., Welch, T., 2006. Iowa’s experience with road diet measures: use of Bayesian approach to assess impacts on crash frequencies and crash rates. Transportation Research Record 1953, 163–171. Persaud, B., Lyon, C., 2007. Empirical Bayes before–after safety studies: lessons learned from two decades of experience and future directions. Accident Analysis and Prevention 39 (3), 546–555. Persaud, B., Lan, B., Lyon, C., Bhim, R., 2009. Comparison of empirical Bayes and full Bayes approaches for before–after road safety evaluations. In: Presented at the 89th Annual Meeting of the Transportation Research Board, Washington, D.C. Sayed, T., deLeur, P., Sawalha, Z., 2004. Evaluating the Insurance Corporation of British Columbia Road Safety Improvement Program. Transportation Research Record 1865, 57–63. 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. Spiegelhalter, D., Thomas, A., Best, N., Lunn, D., 2005. WinBUGS User Manual. MRC Biostatistics Unit, Cambridge. . Wilson, C., Willis, C., Hendrikz, J.K., Bellamy, N., 2006. Speed enforcement detection devices for preventing road traffic injuries. Cochrane Database of Systematic Reviews 2006, Issue 2. Art. No.: CD004607.