NEW TECHNIQUES FOR DEVELOPING SAFETY PERFORMANCE FUNCTIONS by Karim El-Basyouny M.A.Sc., University of British Columbia, 2006 B.Sc., United Arab Emirates University, 2003 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) February 2011 © Karim El-Basyouny, 2011 ABSTRACT While motorized travel provides many benefits, it can also do serious harm in the form of roadrelated collisions. The problem affects millions of human lives and costs billions of dollars in economic and social impacts every year. The problem could be addressed thorough several approaches with engineering initiatives being recognized as the most sustainable and cost effective. However, the success of the engineering approaches in reducing collision occurrences hinges upon the existence of reliable methods that provide accurate estimates of road safety. Currently, Safety Performance Functions (SPFs) are considered by many as the main tool in estimating the safety levels associated with different road entities. Therefore, the research in this thesis focuses on addressing key issues related to the development of SPFs for i) collision data analysis and ii) collision intervention analysis. Some of the key issues addressed include: 1) adding spatial effects to SPFs thereby recognizing the evident spatial nature of road collisions, 2) fitting hierarchical models to allow inference to be made on more than one level, 3) recognizing the multivariate nature of collisions as most data are available by collision type or severity and modeling the data as such, 4) identifying and accounting for outliers in the development of SPFs, 5) developing a novel evaluation methodology to estimate the effectiveness of safety countermeasures when subject to data limitations, and 6) compare different tools for investigating the safety change in treated sites due to the implementation of safety countermeasures. The applications of the various models have been demonstrated using several collision datasets and/or safety programs. The results provide strong evidence for (i) incorporating spatial effects in SPFs, (ii) clustering road segments or intersections into homogeneous groups (e.g., corridors, zones, districts, municipalities, etc.) and incorporating random cluster parameters in SPFs, (iii) developing robust multivariate models with multiple covariates for modeling collisions by severity and/or type concurrently, and (iv) the effectiveness of the proposed full Bayes safety assessment methods that account for several theoretical and ii practical issues concurrently. In addition to the improvement in goodness of fit, the proposed models have also improved inference and precision of expected collision frequency. iii PREFACE a. Articles published in refereed journals 1. El-Basyouny, K. and Sayed, T. (2009) Collision Prediction Models using Multivariate Poisson-lognormal Regression. Accident Analysis and Prevention, 41, 820-828. 2. El Basyouny, K. and Sayed, T. (2009) Accident Prediction Models with Random Corridor Parameters. Accident Analysis and Prevention, 41, 1118-1123. 3. El-Basyouny, K. and Sayed, T. (2009) Urban Arterial Accident Prediction Models with Spatial Effects. Transportation Research Record: Journal of the Transportation Research Board, 2102, 27-33. 4. El-Basyouny, K. and Sayed, T. (2010). A Method to Account for Outliers in the Development of Accident Prediction Models. Accident Analysis and Prevention, 42, 1266– 1272. 5. El-Basyouny, K. and Sayed, T. (2010). A Full Bayes Approach to Before-After Safety Evaluation with Matched Comparisons: A Case Study of Stop-Sign In-Fill Program. Transportation Research Record: Journal of the Transportation Research Board, 2148, 1–8. 6. El-Basyouny, K. and Sayed, T. (2011). A Multivariate Intervention Model with Random Parameters among Matched Pairs for Before-After Safety Evaluation. Accident Analysis and Prevention, 43, 87–94. b. Articles submitted to refereed journals 1. El-Basyouny, K. and Sayed, T. (2011). Full Bayes Before-After Safety Evaluations: Insight into Odds Ratio Calculation and Decomposition. Under Review. 2. El-Basyouny, K. and Sayed, T. (2011). A Non-Linear Intervention Model for Full Bayes Before-After Safety Evaluations. Under Review. iv TABLE OF CONTENTS ABSTRACT ................................................................................................................................. ii PREFACE .................................................................................................................................. iv TABLE OF CONTENTS .............................................................................................................. v LIST OF TABLES....................................................................................................................... xi LIST OF FIGURES ................................................................................................................... xii ACKNOWLEDGMENTS........................................................................................................... xiii DEDICATION ........................................................................................................................... xiv 1 INTRODUCTION ................................................................................................................ 1 1.1 Background ................................................................................................................ 1 1.1.1 The Road Safety Problem ................................................................................. 1 1.1.2 Engineering Approaches for Managing the Road Safety Problem ..................... 3 1.2 Research Objectives .................................................................................................. 6 1.3 Research Issues and Goals ....................................................................................... 8 Part (I): Safety Performance Functions (SPFs) for Collision Data Analysis ......................... 8 Part (II): Safety Performance Functions (SPFs) for Collision Intervention Analysis ............11 1.4 2 Thesis Structure ........................................................................................................14 LITERATURE REVIEW .....................................................................................................15 2.1 Background ...............................................................................................................15 2.2 Collision Modeling .....................................................................................................15 2.2.1 Poisson Regression Models .............................................................................17 2.2.2 Poisson-Gamma (Negative Binomial) Regression Models................................18 2.2.3 Poisson-Lognormal Regression Models ...........................................................19 2.2.4 Zero-Inflated Regression Models ......................................................................20 2.3 Enhanced Collision Modeling ....................................................................................21 v 2.3.1 Variable Variance Models.................................................................................22 2.3.2 Random Effects Models ...................................................................................22 2.3.3 Random Parameters Models ............................................................................23 2.3.4 Multivariate Models ..........................................................................................25 2.3.5 Finite Mixture/Markov Switching Models ...........................................................27 2.4 Safety Performance Functions (SPFs) ......................................................................28 2.4.1 Safety Performance Functions for Collision Data Analysis ...............................28 2.4.2 Safety Performance Functions for Collision Intervention Analysis ....................32 2.5 Parameter Estimation Methods .................................................................................35 2.5.1 Empirical Bayes (EB) Method ...........................................................................36 2.5.2 Full Bayes (FB) Method ....................................................................................38 2.6 Summary ..................................................................................................................42 PART I: SAFETY PERFORMANCE FUNCTIONS FOR DATA ANALYSIS ...............................43 3 4 SPATIAL MODELS ...........................................................................................................43 3.1 Background ...............................................................................................................43 3.2 The Models ...............................................................................................................46 3.2.1 Spatial Poisson Models ....................................................................................46 3.2.2 Gaussian Conditional Auto-Regressive (CAR) Models .....................................46 3.2.3 Multiple Membership (MM) Models ...................................................................47 3.2.4 Extended MM (EMM) Models ...........................................................................47 3.3 Spatial 'Intraclass' Correlation ...................................................................................48 3.4 Prior Distributions ......................................................................................................48 3.5 The Data ...................................................................................................................49 3.6 Results and Discussion .............................................................................................53 3.7 Summary ..................................................................................................................56 RANDOM PARAMETERS MODELS .................................................................................58 vi 5 4.1 Background ...............................................................................................................58 4.2 The Models ...............................................................................................................59 4.2.1 Poisson-Lognormal with Random Corridor Effects (PLN-VC) ...........................59 4.2.2 Poisson-Lognormal with Random Parameters among Pairs (PLN-RP).............60 4.3 Prior Distributions ......................................................................................................61 4.4 The Data ...................................................................................................................61 4.5 Results and Discussion .............................................................................................62 4.6 Summary ..................................................................................................................66 MULTIVARIATE MODELS ................................................................................................67 5.1 Background ...............................................................................................................67 5.2 The Multivariate Poisson-Lognormal Model ..............................................................70 5.3 Multivariate Identification ..........................................................................................70 5.4 Prior Distributions ......................................................................................................72 5.5 Model Comparisons ..................................................................................................73 5.6 The Data ..................................................................................................................73 5.7 Results and Discussion .............................................................................................75 5.7.1 Model Results ..................................................................................................75 5.7.2 Omitted Variable Bias .......................................................................................79 5.7.3 Multivariate Identification of Hazardous Locations ............................................80 5.8 6 Summary ..................................................................................................................82 OUTLIER RESISTANCE MODELS ...................................................................................84 6.1 Background ...............................................................................................................84 6.2 The Methodology ......................................................................................................86 6.2.1 The Scale-Mixture Model..................................................................................86 6.2.2 The Two-Component Mixture ...........................................................................86 6.2.3 Detecting Outliers .............................................................................................87 vii 6.2.4 Down-Weighting Outliers ..................................................................................88 6.3 Prior Distributions ......................................................................................................89 6.4 The Data ...................................................................................................................89 6.5 Results and Discussion .............................................................................................90 6.5.1 Outlier Detection and Weights ..........................................................................90 6.5.2 Parameter Estimates ........................................................................................93 6.5.3 Results of the Scale and Two-Component Mixture Models ...............................94 6.5.4 Results of the Outlier Rejecting Approach ........................................................94 6.6 Summary ..................................................................................................................95 PART II: SAFETY PERFORMANCE FUNCTIONS FOR INTERVENTION ANALYSIS..............97 7 TWO-WAY INTERACTION MODELS................................................................................97 7.1 Background ...............................................................................................................97 7.2 The Models ...............................................................................................................99 7.2.1 The Poisson Lognormal Model (PLN) .............................................................100 7.2.2 The Poisson Lognormal Model with a Random Zone Effects (PLN-VC) .........100 7.2.3 The Poisson Lognormal Model with Random Parameters among Zones (PLN- RP)……………….............................................................................................................101 7.2.4 8 The Poisson Model with Random Parameters among Zones (Poisson-RP) ...101 7.3 The Data .................................................................................................................101 7.4 Results and Discussion ...........................................................................................103 7.5 Summary ................................................................................................................108 LINEAR INTERVENTION MODELS ................................................................................110 8.1 Background .............................................................................................................110 8.2. The Models .............................................................................................................111 8.2.1 The Poisson-Lognormal Intervention (PLNI) Models ......................................111 8.2.2 The PLNI Models with Random Parameters among Matched Pairs................112 viii 8.3 Odds Ratio Decomposition......................................................................................113 8.4 Case Studies...........................................................................................................115 8.4.1 Iowa Study .....................................................................................................115 8.4.2 British Columbia Study ...................................................................................118 8.5. 9 Summary ................................................................................................................122 MULTIVARIATE INTERVENTION MODELS ...................................................................125 9.1 Background .............................................................................................................125 9.2 The Models .............................................................................................................126 9.2.1 The Multivariate Poisson-Lognormal Intervention (MVPLNI) Model ................126 9.2.2 The MVPLNI Model with Random Parameters among Pairs (MVPLNI-RP) ....126 9.3 Prior Distributions ....................................................................................................127 9.4 The Data .................................................................................................................127 9.5 The Odds Ratios .....................................................................................................129 9.6 Results and Discussion ...........................................................................................130 9.6.1 Model Estimates .............................................................................................130 9.6.2. The Site-Level and Treatment-level Odds Ratios ...........................................133 9.7 10 Summary ................................................................................................................135 NON-LINEAR INTERVENTION MODELS: THE KOYCK MODEL ................................138 10.1 Background .............................................................................................................138 10.2 The Koyck Model ....................................................................................................139 10.3 Prior Distributions ....................................................................................................141 10.4 Decomposing the Odds Ratio under the Koyck Model ............................................141 10.5 The Data .................................................................................................................145 10.6 Results and Discussion ...........................................................................................146 10.6.1 Treatment Impact under PLNI Models ............................................................147 10.6.2 Treatment Impact under Koyck Models ..........................................................149 ix 10.7 11 Summary ................................................................................................................151 CONCLUSIONS, CONTRIBUTIONS & FUTURE RESEARCH .....................................153 11.1 Summary and Conclusions .....................................................................................153 11.2 Research Contributions ...........................................................................................158 11.3 Future Research .....................................................................................................159 11.3.1 Network Safety Performance Functions .........................................................159 11.3.2 Implications for Collision Modification Factor Estimations and Benefit-Cost Analysis………….............................................................................................................161 11.3.3 Improvements to Safety Performance Functions ............................................162 REFERENCES .......................................................................................................................165 APPENDICES .........................................................................................................................181 Appendix A: Mathematical Justification for Model Comparison............................................181 Appendix B: Odds Ratio Decomposition under Linear Intervention Model ...........................183 Appendix C: Odds Ratio Decomposition under Nonlinear 'Koyck' Intervention Model .........184 x LIST OF TABLES Table 3.1 Model covariates and their corresponding description ........................................50 Table 3.2 Statistical summary of data set (n = 281 road segments) ...................................51 Table 3.3 Parameter estimates for the PLN, CAR, MM and EMM models ..........................53 Table 4.1 Model covariates and their corresponding description ........................................61 Table 4.2 Statistical summary of data set (n = 392 road segments) ...................................61 Table 4.3 Parameter estimates and 95% credible intervals for all models ..........................63 Table 5.1 Statistical summary of data set (99 intersections) ...............................................74 Table 5.2 Univariate models’ statistics ...............................................................................76 Table 5.3 Multivariate models’ statistics .............................................................................76 Table 5.4 Number of hot spots selected under different PLN models .................................80 Table 6.1 Parameter estimates and their associated standard deviations ..........................93 Table 7.1 Zones selected for evaluation ...........................................................................103 Table 7.2 Parameter estimates, standard deviations and 95% credible intervals .............106 Table 7.3 Zone-level odds ratios for Poisson-RP ............................................................107 Table 8.1 Annual injury collisions frequency for treated sites ...........................................119 Table 9.1 Selected treatment sites ...................................................................................128 Table 9.2 Annual collision frequency by severity level for treated sites.............................129 Table 9.3 Parameter estimates and standard errors.........................................................131 Table 9.4 Site-level odds ratios under MVPLNI-RP ..........................................................134 Table 9.5 Treatment-level odds ratios under MVPLNI-RP ................................................135 Table 10.1 Yearly collision injury rates for treated sites ......................................................146 xi LIST OF FIGURES Figure 3.1 GVRD study area. ..............................................................................................52 Figure 4.1 Regression coefficients effects ...........................................................................65 Figure 5.1 Wireframe plots ..................................................................................................75 Figure 5.2 Standard deviations (sd) of expected collision frequencies by severity and model ...........................................................................................................................78 Figure 5.3 95% credible ellipses and rectangle for (ln( λ PDO ),ln( λ I + F )) , using the data for Site 56 ................................................................................................................81 Figure 6.1 Scatter plot of Mahalanobis distance by site. ......................................................91 Figure 6.2 Scatter plots of scaling weights by site ...............................................................92 Figure 7.1 Zones selected for evaluation in the city of Vancouver. ....................................103 Figure 8.1 Monthly treatment components for a median site under H1 (w/o a jump), Li et al. (2008)...............................................................................................................116 Figure 8.2 Monthly combined treatment components for a median site under H1 & H2, Li et al. (2008). .........................................................................................................118 Figure 8.3 Yearly treatment components under the PLNC model (completion year 2004). 121 Figure 10.1 Novelty effects per year under the Koyck model. ..............................................142 Figure 10.2 Yearly treatment components by δ at ψ= 0.4, under the Koyck model..............144 Figure 10.3 Yearly treatment components by ψat δ = 0.3, under the Koyck model..............144 Figure 10.4 Yearly treatment components under the PLNC model (completion year 2002). 148 Figure 10.5 Yearly treatment components under the PLNC and Koyck models (completion year 2002). .......................................................................................................149 Figure 11.1 Richmond network boundaries .........................................................................161 xii ACKNOWLEDGMENTS I wish to thank my dissertation supervisor Professor Tarek Sayed, Civil Engineering Department, University of British Columbia, for his guidance, encouragement, patience, and financial support. Indeed, he has been a tremendous help for the past six years. The many hours of discussions, in which he showed his enthusiasm and positive attitude towards science in particular and life in general, kept me on the right track. I also owe sincere thanks to my supervisory committee for their advice, support, and guidance during the course of my graduate studies. I would also like to acknowledge the help of my professors at UBC, who impacted my graduate study through their coursework. Special thanks are due to my professors at the Civil Engineering Department, UAE University, for fostering me during my undergraduate study and for their encouragement to pursue my postgraduate degrees. In addition, I wish to thank my colleagues and friends in Canada and Egypt for their support and friendship. Finally, but foremost, I am deeply grateful for my parents, sister and wife for their unconditional support. Their prayers, help, and love continue to inspire me to achieve greater things in life. xiii DEDICATION (change color to white) To my parents, sister and lovely wife xiv 1 INTRODUCTION This chapter provides a general introduction to the thesis and is divided into three parts. The first part presents background information that is necessary to understand the significance of the research problem. The second part discusses the research objectives, issues and goals. The chapter concludes by describing the research structure which provides a roadmap for the different components of the thesis. 1.1 Background This section presents the magnitude of the road safety problem from both a social and economic perspective. Moreover, the section summaries some of the approaches utilized to manage the problem and the most common techniques employed. 1.1.1 The Road Safety Problem The significant impact of motorized vehicles on expediting the transportation of goods and people and on the economic and social development of many countries cannot be denied. However, while motorized travel provides many benefits, it can also cause serious harm in the form of road-related collisions. The problem is common to both industrialized and developing countries as it affects millions of human lives and costs billions of dollars in economic and social impacts every year. As a result, road-related injuries and fatalities have been recognized as a major health hazard around the globe. The increase in number of road deaths and injuries, particularly in developing countries, has forced many to label the road safety problem as the ‘neglected epidemic’ (WHO, 2004). It is estimated that every year 1.3 million fatalities occur due to road collisions worldwide (WHO, 2004). Meanwhile, millions of people sustain injuries with some suffering permanent disabilities. The number of fatalities in road traffic collisions for most countries ranges from 50 to 150 per million population (based on 2004 statistics). Injuries from traffic related incidents are estimated at approximately 50 million per year that is the combined population of five of the 1 world’s largest cities. In economic terms, the cost of road collisions injuries is estimated to range from 1% to 2% of Gross National Product (GNP) worldwide. The direct economic costs of global road collisions have been estimated at a staggering US$ 518 billion annually (Jacobs et al., 2000). The social and economic burden of road collisions in North America is also enormous. The Center for Disease Control has ranked transportation collisions as the third leading cause of death in the United States (after cancer and heart disease) each year from 1991 to 2000 (CDC, 2003). Collision costs in the United States are estimated at US$ 164 billion in 2005 dollars (Cambridge Systematics Inc., 2008). In Canada, vehicle collisions are a serious health and safety issue. Collisions were the seventh leading cause of Potential Years of Life Lost (National Cancer Institute of Canada, 2001) and they are the leading cause of death in Canadian children (American Academy of Pediatrics, 2002; Howard, 2002). Moreover, collision costs were estimated at CDN$ 62.7 billion each year representing approximately 4.9% of Canada's 2004 GDP (Vodden et al., 2007). Therefore, it is not surprising that many of the Canadian provinces suffer from road safety problems. British Columbia (B.C.) has a population of approximately 4.3 million people with 3.7 million licensed vehicles. Because of its size and varying topography, B.C. requires thousands of kilometers of provincial highways to connect its communities. In 2006, there were 370 fatalities, approximately 20,000 injuries, and 30,000 property damage only collisions in B.C. On average a fatality occurred every 21.2 hours and 79 people were reported injured each day in 2006 (Traffic Collision Statistics, 2006 annual report). Historical data show that the number of fatalities has been decreasing gradually since the late 80’s. However, the total number of collisions has continued to increase. The Insurance Corporation of British Columbia estimates that the direct annual cost of collision to the province exceeds two billion dollars. These statistics demonstrate that the road safety problem has a far reaching social and economic effect on the global community, affecting the lives and welfare of millions of people 2 around the world. If current (unsafe) trends continue unchanged, the number of people killed and injured on the world’s roads is projected to increase by more than 65% between 2000 and 2020 (WHO, 2004). Conversely, if fatality rates per vehicle were reduced by 30% by 2020, more than 2.5 million lives could be saved and 200 million injuries avoided. Consequently, the importance of reducing the social and economic costs of road collisions cannot be overstated. There are several approaches to address the road safety problem with engineering initiatives being recognized as the most sustainable and cost effective. The next section discusses some of the approaches to manage the problem with specific emphasis on engineering. 1.1.2 Engineering Approaches for Managing the Road Safety Problem At the highest level of abstraction the road system is usually represented by three components: the driver, the vehicle, and the road (Sayed et al. 1995). The dynamic and complex nature of the road system can be interpreted through the interactions and relationships that exist among the road system components. Road safety and other system performance measures (delay, congestion etc.) are regarded as an emergent property of the road system. These measures are generated by inappropriate or unsafe interactions between or within the components of the road system (Archer, 2001). Therefore, the first step in determining the safety of a particular road entity is to relate collisions to the failure in one or more combination of the three road system components. Several studies have shown that driver error contributes to approximately 90-95% of all collisions on the road. The findings may lead people to deduce that safety initiatives focusing on the “driver” should significantly reduce collision occurrence. According to this argument, information and publicity should form the backbone of road injury prevention. However, human driving behavior is not only governed by individual knowledge and skills but also by the surrounding environment. Indirect influences, such as the design and layout of the road, the type and nature of the vehicle, and traffic laws and their enforcement – or lack thereof – will affect the driving behavior. For this reason, the use of information and 3 publicity on their own was generally shown to be unsuccessful in reducing road collisions. Moreover, the effectiveness of driver-based initiatives has been proven most effective in the short-term but unsustainable in the long run (Sayed et al., 2010). Accordingly, many road authorities shifted their focus to improving the safety of the roads through the implementation of road-based safety initiatives. These initiatives have been found to be sustainable and effective in reducing, not only the frequency of collisions, but also their severity level. Broadly, these initiatives are divided into two categories: reactive initiatives (in response to existing problems) and proactive initiatives (taken to prevent the emergence of problems) (de Leur, 2001; Lovegrove, 2006). The reactive approach attempts to improve existing unsafe locations in order to reduce frequency and severity of collisions. The approach is an efficient way to improve road safety performance for existing road infrastructure. However, the approach has a major drawback since significant collision history must exist before any road improvements are implemented (de Leur, 2001; Lovegrove, 2006). An example of a reactive approach would be the identification of collision prone locations or “hot-spots” program. This program identifies locations that are considered hazardous or prone for exhibiting a higher than normal number of collisions. In order to identify and address “hot-spots”, historical collision data are required making this approach reactive in delivery. In contrast, the proactive approach attempts to prevent unsafe road conditions from occurring by implementing road safety measures early in the traffic planning and design process. Rather than waiting for road safety problems to occur, the proactive approach attempts to evaluate the safety throughout each stage of the planning process, thus minimizing the safety risks. Proactive approaches to safety advocate the development of tools that improve safety decision-making in the planning, design and implementation phases (de Leur, 2001). The success of the reactive and proactive approaches in reducing collision occurrences hinges upon the existence of consistent methods that provide reliable estimates of road safety. 4 Currently, Safety Performance Functions (SPFs) are considered by many as the main tool in estimating the safety levels associated with different road entities. SPFs are mathematical models that are statistically developed to conduct collision data analysis by attempting to explain collision occurrence on various road facilities as a function of the traffic and geometric characteristics of these facilities. More recently, SPFs have been proposed to conduct collision intervention analysis by relating the collision occurrence on various road facilities as a function of time, treatment and interaction effects. These SPFs acknowledges that safety treatment (intervention) effects do not occur instantaneously but are spread over future time periods and are used to capture the effectiveness of safety interventions. Under the reactive approach, SPFs are an invaluable tool that can be used to help road safety improvement programs (RSIPs) succeed in achieving their goals. The success of RSIPs relies on their abilities to accurately estimate the true safety level of any road location. SPFs are used to identify locations that are truly hazardous, since not every high collision frequency location is a hazardous location “hot-spot”. SPFs are also valuable in the diagnosis and remedy stages of RSIPs. For example, successful remedy of detected hazardous locations “hot-spots” requires the ability to estimate the change in safety level caused by changes in various road characteristics, which could be easily carried out using a SPF. More importantly, SPFs can assist RSIPs in quantifying the effectiveness of an implemented safety countermeasure in terms of collision reduction and economic benefits by analyzing the results of before-after evaluations or by conducting cross-sectional evaluations. SPFs are not just valuable to the success of the reactive approach of dealing with road safety problems; they are of vital importance to the success of the proactive approach. The primary objective of the proactive approach is to ensure that road safety is an explicit priority of land use and transportation planning policies. This objective cannot be effectively achieved without the existence of models that can be used to assess the impact of these policies on road safety. For example, SPFs are used to produce collision modification factors for estimating the 5 safety implications of future decisions related to changes in traffic operations and highway design. However, there are many issues related to collision modeling that either deserve further attention or have not been dealt with adequately in the road safety literature. The research issues presented in this thesis are concerned with a number of concepts, theories and methods related to developing models for road safety analysis. 1.2 Research Objectives The success of both the reactive/proactive approaches depends on an accurate estimation of the safety level associated with existing or planned road facilities. Important decisions are made and considerable resources are allocated based on the results of the studies that evaluate these approaches. However, the results are only as reliable as the safety estimation technique utilized. SPFs are considered by many as the main tool in estimating the safety levels associated with different road entities (intersection or road segments). Consequently, every effort should be made to ensure the development of sound and consistent SPFs. The underlying principle for a more effective safety evaluation strategy is to develop predictive models that are capable of accounting for the different variations (i.e., temporal, spatial, etc.) and properties (i.e., multivariate dependent and independent variables, errors-inpredictors, etc.) that exist in collision data. Many of the research issues investigated in this thesis have been overlooked in the past due to limitations in the modeling tools. However, recent advances in statistical modeling techniques have facilitated the fitting of more advanced models. The use of these models is justified as long as there are improvements made to the model’s predictive power (model inference) and fit to the data (goodness-of-fit). It is implicitly assumed that the model’s hypothesis should be valid on a theoretical basis and should conform to engineering judgment. 6 The importance of improving the predictive abilities of safety models cannot be overstated. Safety models are used primarily to identify and quantify road safety problems and evaluate the efficiency of road safety countermeasures. Consequently, improving SPFs results in refining the identification and ranking of collision prone locations or “hotspots” and provides more reliable estimates for location specific predictions and before-after safety evaluation. The research in this thesis has two objectives. The first objective focuses on addressing key issues related to the development of SPFs for collision data analysis. The goal is to improve the predictive power and fit of existing models. Some of the key modeling issues addressed include: 1) adding spatial effects to SPFs thereby recognizing the evident spatial nature of road collisions, 2) fitting hierarchical models to allow inference to be made on more than one level, 3) recognizing the multivariate nature of collisions as most data are available by collision type or severity and modeling the data as such, and 4) identifying and accounting for outliers in the development of SPFs. The second objective focuses on developing SPFs for collision intervention analysis of the results of Controlled Before-After (CBA) studies. The goal is to introduce effective methods for road safety assessment and to provide a better understanding of the factors that influence the frequency and severity of collisions in order to propose more successful safety-related countermeasures. Some of the key issues addressed include: 1) developing a novel evaluation methodology to estimate the effectiveness of safety countermeasures when subject to data limitations, 2) providing a better understanding of the role of model parameters and their effects on the overall odds ratio estimation and consequently safety treatment effectiveness, 3) simplifying the calculation of the odds ratio and subsequently the safety treatment effectiveness by providing straightforward equations in terms of model parameters, 4) comparing different tools for investigating the safety change in treated sites due to the implementation of safety countermeasures, and 5) investigating different modeling formulations for the before-after safety 7 evaluations based on the topics discussed in the first part of the thesis (i.e. random effects, spatial effects, multivariate models). The next section highlights some of the research issues and goals. 1.3 Research Issues and Goals Part (I): Safety Performance Functions (SPFs) for Collision Data Analysis 1. Spatial Models To develop SPFs, it is useful to classify the likely effects on the distribution of collisions into three separate categories: (i) Poisson variation; (ii) heterogeneity (extra-variation), due to withinsite effects reflecting its individual characteristics; and (iii) spatial effects (spatial correlation), since neighboring sites typically have similar environmental and geographic characteristics and thereby form a cluster that has similar collision occurrence. Most of the literature related to the development of SPFs accounts for the first two types of variation; namely: Poisson variation and Poisson extra-variation. Recently, the need to include the spatial effects in the development of SPFs has been gaining attention in the literature. It has been demonstrated that spatial dependence can be a surrogate for unknown and relevant covariates, thereby improving model estimation. Research has shown that ignoring spatial dependence leads to underestimation of variability. Despite the evident spatial nature of road collisions very few efforts have attempted to develop SPFs with spatial effects. The first goal of this research is to investigate the consequences of including spatial effects in SPFs. Two types of spatial modeling techniques, namely, the Gaussian conditional autoregressive (CAR) and the multiple membership (MM) models, will be compared to the traditional Poisson log-normal (PLN) model. A variation of the MM model (Extended MM or EMM) will be also investigated to study the effect of clustering segments within the same corridor on spatial correlation. All models are compared in terms of their inference and goodness-of-fit. 8 2. Random Parameters Models In contrast to traditional SPFs, which fit one regression model to the dataset, the random parameters approach develops different regression models for each individual site in the data set. The approach is particularly appealing when used to explicitly account for the variations of the effect of variables across observations. Constraining the parameters to be constant when they actually vary across observations could lead to inconsistent and biased parameter estimates. Given the potential heterogeneity in collision-frequency data, the random-parameters approach may be useful in several cases. However, using the random parameters approach to develop different regression models for each individual site in the data set may seem to be a case of over-fitting, which might affect the model’s generality. To benefit from the advantages of using the random parameters approach and avoid over-fitting, it is suggested to cluster the road entities into rather homogeneous groups (e.g., districts, municipalities, zones, corridors, etc.) and fit a different regression curve for each group rather than for each individual site. The second goal of this research is to investigate the consequences of incorporating spatial cluster effects in SPFs using random parameters models. Two models that account for the cluster variation through the mean and variance of an extended Poisson log-normal (PLN) model are compared to the traditional PLN model. Random parameters PLN models are proposed to account for the cluster variation through the mean by fitting a different regression curve for each cluster. This model focuses on explaining part of the extra-variation through improvements in the mean function. Alternatively, a spatial PLN model is proposed to account for the cluster variation through the variance by using an additional variance component. The model focuses on explaining part of the extra-variation through improvements in the variance function. All models are compared in terms of their inference and goodness-of-fit. 3. Multivariate Models Data on the number of collisions at a particular site are usually classified by severity (e.g., fatal, major injury, minor injury or property damage only), by the number of vehicles involved (e.g., 9 single or multiple), and/or by the type of collision (e.g., angle, head-on, rear-end, sideswipe or pedestrian-involved), etc. In spite of the multivariate nature of such data sets, they have been mostly analyzed by modeling each category separately, without taking into account the correlations that probably exist among the different levels. These correlations may be caused by omitted variables, which can influence collision occurrence at all levels of classification, or from ignoring shared information in unobserved error terms. Such univariate treatment of correlated counts as independent can lead to imprecise analysis of road safety. The third goal of this research is to investigate the nature of the correlations across collision severity levels and their influence on safety analyses. A bivariate Poisson log-normal (PLN) model is compared with the independent (separate) univariate PLN models. Moreover, a new multivariate hotspots identification technique, which generalizes the univariate posterior probability of excess that has been commonly applied in the literature, is proposed. All models are compared in terms of their goodness-of-fit, inference, identification of hotspots and the precision of expected collision frequency. 4. Mixture Models Collision 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. Typically, there are three approaches for analyzing data with potential outliers: (i) the “outlier ignoring” approach, (ii) the “outlier rejecting” approach, and (iii) the “outlier accommodating” approach. The problem with the outlier ignoring approach is that while including some non-influential outliers may be tolerable, influential outliers would significantly impact the model fit leading to faulty inference. On the other hand, the main problem with the popular outlier rejecting approach is its failure to reflect the uncertainty in the exclusion process, which results in too small standard errors that might lead to flawed inferences. In contrast, the third approach “outlier accommodating” uses robust methods to down-weight potential outliers. Such robust models have been proposed in 10 the statistical literature to down-weight the outlying observations rather than rejecting them. Few studies have discussed the effect of outliers and how to handle them while developing SPFs. The fourth goal of this research is to compare the three approaches for analyzing collision data with potential outliers. Two outlier resistance mixture models are introduced based on the Multivariate Poisson Lognormal (MVPLN) approach presented earlier. All models are compared in terms of their goodness-of-fit, inference, and precision of expected collision frequency. Part (II): Safety Performance Functions (SPFs) for Collision Intervention Analysis 1. Two-way Interaction Models In many circumstances data limitations might restrict the application of Bayesian analysis to conduct Controlled Before-After (CBA) safety studies. These data limitations are usually attributed to i) the lack of sufficient data to develop the SPFs that are necessary to conduct an Empirical Bayes (EB) analysis or ii) the availability of collision data in aggregated (either in space or time) forms that limits the use of a Full Bayes (FB) analysis. The fifth goal of this research is to develop a FB approach to conduct a CBA evaluation when subject to data restrictions. A two-way interaction SPF is proposed based on the PoissonLognormal (PLN) hierarchy with covariates representing the time, treatment (intervention) and interaction effects. To account for the correlation in collision counts between sites within comparison-treatment pairs, the FB two-way interaction SPF is extended to allow for the variation among the comparison-treatment pairs through the mean (using a random parameters model) and variance function (using a random effects model). All models are compared in terms of their goodness-of-fit and inference. 2. Linear Intervention Models When collision data is available in disaggregate (in time and space) form, a Full Bayes (FB) linear intervention SPF has been proposed to conduct a CBA safety evaluation. The approach acknowledges that the treatment (intervention) effects do not occur instantaneously but are spread over future time periods. Using dynamic regression models the approach identifies the 11 lagged effects of the intervention to describe their response. The literature has identified two types of linear intervention SPFs based on the Poisson-Lognormal hierarchy. Such models are denoted PLNI. The two models differ in the way the treatment impact is specified, the first deals with an immediate treatment impact (PLNI with jump) while the second model assumes a gradual treatment impact (PLNI without jump). In all forms, linear slopes were implicitly assumed to represent the time and treatment effects across the treated and comparison sites. However, there is still a lack of complete understanding of the role that the model parameters play in the estimation of the odds ratio (OR) and subsequently on treatment effectiveness. This is partially due to the perceived notion that the FB methodology is complex and that a particular algorithm is required to compute the OR. The sixth goal of this research is to offer a novel approach to decompose the odds ratio under the linear intervention form into basic components related to direct and indirect treatment effects, where the latter are imposed on the predicted collisions through traffic volumes and other site characteristics that vary with time (e.g., weather conditions). Also, the consequences of the linearity assumption on these components, and subsequently on the OR, are investigated. Another objective is to provide straightforward equations in terms of model parameters for the computation of the OR and its associated components without resorting to additional algorithms. The decomposition and computation of the OR are illustrated using two case studies. 3. Multivariate Linear Intervention Models The incorporation of random parameters models and multivariate models into the development of SPFs were discussed in part one of the thesis. However, no studies have attempted to incorporate both approaches in intervention analysis. The seventh goal of this research is to investigate the consequence of developing intervention SPFs using random parameters and multivariate models. Both approaches are combined to assess the factors that influence the frequency and severity of collisions. Two models are 12 considered for conducting CBA evaluations. The first model uses a Multivariate Poisson Lognormal Intervention (MVPLNI) model for collision counts by severity-level with covariates representing traffic volume, time, treatment (intervention) and interaction effects. The second model extends MVPLNI by using random parameters to account for the comparison-treatment pairing. All models are compared in terms of their goodness-of-fit and inference. 4. Non-Linear Intervention Models: The Koyck Model As indicated above, the use of a Poisson-Lognormal Intervention (PLNI) model in the context of a CBA safety evaluation was investigated in the literature. The linear intervention approach uses dynamic regression models to identify the lagged effects of the intervention in order to describe their response. The approach acknowledges the fact that the treatment (intervention) effects do not occur instantaneously but are rather spread over future time periods. In all of its forms, linear slopes were tacitly assumed to represent the time and treatment effects across the treated and comparison sites. It is important to note that the linear slopes of the intervention model can only furnish some restricted treatment profiles which can impact the estimation of Collision Modification Factors (CMFs) and the economic evaluation of safety programs and countermeasures. The eighth goal of this research is to challenge the linear slope assumption that is implicitly assumed in the PLNI SPF by proposing an alternative formulation to model the lagged treatment effects that are distributed non-linearly over time. Two models are considered for conducting CBA evaluations. The first model uses a Poisson Lognormal Intervention (PLNI) model whereas the second model is based on the Poisson Lognormal Koyck model. The Koyck model is an alternative dynamic regression form involving a first-order autoregressive (AR1) model that is based on infinite distributed lags. The Koyck model affords a rich family of forms that can accommodate various profiles for the time and treatment effects. The two models are extended to include random parameters to account for the variation among the comparisontreatment pairs. All models are compared in terms of their goodness-of-fit and inference. 13 1.4 Thesis Structure The thesis is divided into eleven chapters. Excluding the introduction, literature review, and concluding chapter, the thesis core elements are divided into two parts. The first part of the thesis focuses on developing Safety Performance Functions (SPFs) for collision data analysis. This part encompasses chapters 3 to 6. The second part of the thesis focuses on developing SPFs to analyze the results of CBA evaluations. Chapters 7 to 10 address issues related to this topic. The first chapter provides an introduction to the thesis by including background information on the road safety problem, an overview of the research objectives, issues and goals, and a description of the thesis structure. Chapter two provides a review of the road safety literature pertaining to collision modeling for purposes of data analysis and program evaluations. However, due to the diversity of topics being handled in this thesis each individual chapter has a more detailed literature review on its topic of discussion. Chapter three investigates the consequence of adding spatial effects to SPFs. Chapter four studies the effects of fitting models with random parameters to account for variation across spatial clusters. Chapter five addresses the multivariate nature of collisions and modeling the data as such. Chapter six focuses on identifying and accounting for outliers in the development of SPFs. Chapter seven proposes a novel SPF to estimate the effectiveness of safety countermeasures when subject to data limitations. Chapter eight highlights the most common SPFs forms to conduct CBA evaluations in the literature and extends them to include models with random parameters among matched comparison-treatment pairs. Chapter nine extends the SPFs developed in chapter eight by including multivariate models to assess the factors that influence the frequency and severity of collisions. Chapter ten compares linear and non-linear SPFs, to conduct CBA safety evaluations, highlighting the difference/similarities and advantages/disadvantages of both approaches. Finally, chapter eleven discusses the research conclusions, contributions and suggestions for future research. 14 2 LITERATURE REVIEW This section presents an overview of different collision modeling issues. The ultimate goal is to bring the reader up to date with the current body of literature by developing a well-structured literature review. However, due to the diversity of topics being handled in this thesis each individual chapter has a more detailed literature review. 2.1 Background Given the magnitude and consequences of the road safety problem, researchers have focused their attention on gaining a better understanding of the factors that affect the probability of collisions. It is hoped that identification of such factors would enable safety analysts to better predict the likelihood of collisions and provide directions for policies and countermeasures to reduce the frequency and severity of collisions. However, the absence of detailed traffic and collision data has been partially responsible for our lack of identifying the cause and effect relationships regarding collision probabilities. To address this problem, researchers framed their analytical approaches to study the factors that affect the frequency of collisions in both geographical space (intersection or road segment) and specified time period (day, month, year, etc.) (Lord & Mannering, 2010). Typically, this is achieved through the use of Safety Performance Functions (SPFs). This chapter provides an overview of the collision modeling methods as well as some of the new modeling techniques that are currently considered the state-of-the-art in safety analysis. Key issues related to the development/calibration of SPFs are also introduced. The information presented in this section will provide the basis for all the work in the thesis. 2.2 Collision Modeling Safety analysts need to analyze collision data to estimate the level of safety at different road entities (segments and intersections) in order to identify hazardous (unsafe) locations and to evaluate the effectiveness of road safety countermeasures. In this context, the term ‘safety’ 15 reflects the ‘true’ underlying collision frequency at a location. There are a variety of methods available to analyze collision data and determine the ‘true’ safety levels. Broadly, these methods are classified into two categories: conventional analysis and probabilistic analysis. Conventional methods assume that observed collisions at a specific site is considered to be an unbiased estimate of the true level of safety at a site. Therefore, observations of collision occurrence at a site can be used to provide an estimate of the mean collision frequency, which can be compared with appropriate standard values to indicate the level of safety. However, the uncertainty associated with the use of crude estimates is a major drawback. Probabilistic methods are distinguished from ‘conventional’ methods in that any parameter in a problem (such as the true collision frequency at a location) is regarded as a random variable with a probability distribution. Probabilistic methods account for the stochastic effects in collision data and recognize that collisions are rare random events. In contrast to the conventional method, the mean collision frequency is never known but is only estimated through the use of safety-based models. This has lead several researchers to discard the use of conventional analysis methods because they suffer from a large statistical uncertainty and for their failure to acknowledge the effects of site-specific attributes (such as traffic and geometric characteristics) on the collision occurrence at a particular site (Sawalha and Sayed, 2006). Within the probabilistic approach, there are several regression techniques to model road collisions. The model development and subsequently the results are strongly affected by the choice of the regression technique used. The earlier models were developed using ordinary or normal linear regression. These models assume a normal error structure for the response variable, a constant variance for the residuals, and the existence of a linear relationship between the response and explanatory variables. Such models were criticized by several researchers (Jovanis and Chang, 1986; Hauer et al., 1988; Miaou and Lum, 1993) indicating that road collisions are discrete, nonnegative and rare events that could not be adequately modeled using conventional linear regression. 16 Such limitations led to the use of Generalized Linear Regression Models. Some of the statistical programs specifically designed for fitting generalized linear models are GLIM (Francis et al., 1993), GENSTAT (Lane et al., 1988), and GENMOD (SAS Institute Inc., 1987). These computer packages can be used for modeling data that follow a wide range of probability distributions such as the Normal, Poisson, binomial, negative binomial, gamma, and many others. Furthermore, the computer packages allow the flexibility of using several non-linear model forms that can be converted into linear forms through the use of several built-in link functions. The Poisson regression model was first adopted to account for the discrete, nonnegative and random nature of collision data. However, due to its limitations and introduction of new numerical procedures, alternative models have been suggested. The next sections briefly discuss some of the known regression models, which will be referenced and used throughout this thesis. 2.2.1 Poisson Regression Models The Poisson distribution is commonly used to model discrete, nonnegative and random count data. Therefore, it was one of the first distributions used to model collisions. Let Y i denote the number of collisions at site i (i =1,…,n) and assume that collisions at the n sites are independent and that Y i | θ i ~ Poisson (θ i ) . (1) where θi is the Poisson parameter. The probability of a site i having yi collisions is given by Pr{Y i = y i | θ i} = e−θ i θ iy i / y i ! (2) The Poisson parameter θi is commonly specified as an exponential function of sitespecific attributes such as exposure, traffic and geometric characteristics (Miaou and Lord, 2003) usually expressed as θ i = exp( X i ' β ) (3) 17 where X i ' is a row vector of covariates representing site-specific attributes and β is a vector of regression parameters to be estimated from the data. In the Poisson regression model, the mean and variance of the count variable are constrained to be equal such that: E (Y i ) =Var (Y i ) =θ i (4) It has been shown (Kulmala, 1995; Cameron and Trivedi, 1998; Winkelmann, 2003) that most collision data are likely to be over-dispersed (the variance is greater than the mean) which makes the application of the simple Poisson regression problematic. Poisson models cannot handle over-dispersion. Assuming a Poisson distribution in the presence of over-dispersed data can underestimate the standard errors of the regression coefficients, which can lead to inflated values of t-test thereby affecting the significance level of the model regression coefficients. This will in turn lead to an incorrect selection of covariates resulting in faulty inference as well as affecting the goodness-of-fit. To overcome some of the problems associated with the simple Poisson regression the Conway-Maxwell-Poisson distribution was proposed to model collision data (Lord et al., 2007). The Conway-Maxwell-Poisson distribution is a generalization of the Poisson distribution and was found to handle both under-dispersed and over-dispersed data sets. Despite that advantage the model has seen limited use since its inception. 2.2.2 Poisson-Gamma (Negative Binomial) Regression Models To overcome some of the problems associated with the Poisson regression several researchers proposed the use of the Poisson-Gamma (PG) hierarchy leading to the Negative Binomial (NB) regression model. The PG or NB regression model is an extension of the Poisson model developed to account for the over-dispersion that commonly exists in collision data. To address over-dispersion for unobserved or unmeasured heterogeneity, it is assumed that θ i = µ i exp(u i ) , µ i = exp( X i ' β ) , (5) 18 where the term exp( u i ) represents a multiplicative random effect. The negative binomial (Poisson-Gamma) model is obtained by the assumption exp(u i ) | κ ~ Gamma(κ , κ ) , (6) where κ is the inverse dispersion parameter. The dispersion (or over-dispersion) parameter is usually referred to as α = 1 / κ . The probability density function of the PG or NB model is given by κ Γ ( yi + κ ) κ Pr{Y i = yi | µ i , κ } = yi !Γ (κ ) κ + µ i y µi i κ +µ . i (7) Under the PG or NB model the mean and variance are given by E (Y i ) = µ i , Var (Y i ) = µ i + µ i2 / κ . (8) When collision data are not over-dispersed, α will go to zero and the Poisson-Gamma model will reduce to the Poisson model. It is therefore clear that the selection between the two models is dependent upon the value of α. The Poisson-Gamma or Negative Binomial regression model has been widely applied in the road safety literature. Arguably, the main reason for the extensive use of this model is that it is simple to compute since the gamma distribution is a conjugate prior to Poisson leading to a gamma posterior distribution which simplifies the posterior analysis considerably. 2.2.3 Poisson-Lognormal Regression Models Recently, several researchers have proposed the use of the Poisson-Lognormal (PLN) model as an alternative to the Poisson-Gamma model to model collision data (Miaou et al., 2003; Lord and Miranda-Moreno, 2008; Aquero-Valverde and Jovanis, 2008). The PLN model, which is commonly used to address over-dispersion for unobserved or unmeasured heterogeneity, is obtained by the assumption exp(u i ) | σ u2 ~ Lognormal (0, σ u2 ) or 2 2 u i | σ u ~ Normal (0, σ u ) (9) 19 where σ u2 denotes the extra Poisson variance. Note that in the PLN model exp(ui) follow a Lognormal distribution whereas in the Poisson-Gamma model exp(ui) followed a gamma distribution. The PLN is a good candidate for modeling collision occurrence in the presence of outliers, since its tails are known to be asymptotically heavier than those of the Gamma distribution (Kim et al., 2002; Lord and Miranda-Moreno, 2008). Other empirical evidence and advantages of the PLN model are discussed in Winkelmann (2003). Under the PLN model E (Y i ) = µ i exp(0.5σ u2 ) , Var (Y i ) = E (Y i ) + [ E (Y i )]2 (exp(σ u2 ) − 1) . (10) A major limitation for the application of the Poisson-Lognormal (PLN) regression model is that it is marginal distribution does not have a closed form as the Poisson-Gamma model. As a result, the PLN has been less popular since it demands more computational effort and few statistical programs are readily available for their calibration. However, enormous progress has been recently achieved in numerical methods combined with the availability of free and easy to use software, which permitted the application of the PLN model to analyze collision data. The PLN regression model offers more flexibility than the Poisson-Gamma model and is relied on extensively in this thesis. 2.2.4 Zero-Inflated Regression Models Collision data can have a high proportion of observed zero counts, particularly when analyzing collisions occurrence in rural areas, and this is known to cause over-dispersion. The issue is problematic when the number of observed zero counts exceed the zero counts that can be tolerated by the simple Poisson regression model. If a large number of sites have reported zero collisions for a given time period of observation this will produce a collision frequency distribution with heavy proportions of zeros. Such data sets are characterized by having a low sample mean problem (Miranda-Moreno, 2006). Several studies tried to address the problem arising from the “excess” zeroes in collision data, by applying zero-inflated probability models such as zero-inflated Poisson (ZIP) and Zero- 20 inflated Negative Binomial (ZINB). The models assume that a dual-state process “safe” and “unsafe” is responsible for generating the collision 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. The authors 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. Analyzing data characterized by a low sample mean is problematic, especially if combined with a small sample size due to the expensive costs of collecting collision data and other relevant variables. The low mean problem has been studied in the road safety literature, see Maher and Summersgill (1996), who showed how such a problem affects the goodness-offit of SPFs; Wood (2002), who proposed a method to test the fit of SPFs 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. Several solutions were proposed to address this problem; some of these solutions will be discussed in the next sections. 2.3 Enhanced Collision Modeling All the regression models described above can be extended to include any of the modeling techniques discussed in the forthcoming sections. 21 2.3.1 Variable Variance Models In the traditional Poisson-Gamma and Poisson-Lognormal models, the dispersion parameter is held fixed. Such an assumption was challenged and various dispersion parameter relationships were examined (Heydecker and Wu, 2001; Miaou and Lord, 2003; Miranda-Moreno et al, 2005; El-Basyouny and Sayed, 2006). In essence, the dispersion parameter α (in Poisson-Gamma) or σ u2 (in Poisson-Lognormal) could be allowed to vary as a function of site-specific attributes similar to the mean function α i = exp(W i ' δ ) or σ ui2 = exp(W i ' δ ) , (11) where W i is a vector of covariates representing site-specific attributes that may or may not include some of the covariates that belong to X i and δ is a vector of regression parameters to be estimated from the data. Research on this topic has concluded that modeling the dispersion parameter depends on how the mean structure is specified. 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, the extra-Poisson variation is not significantly related to these variables (Mitra and Washington, 2007). 2.3.2 Random Effects Models The literature classifies the likely effects on the distribution of collisions into three separate categories: (i) Poisson variation; (ii) heterogeneity (extra-variation), due to within-site effects reflecting its individual characteristics; and (iii) random effects (spatial or temporal correlation). The first two types of variation; namely: Poisson variation and Poisson extra-variation were addressed in previous sections above. Temporal correlation arises from the fact that collision data are usually available per unit time (day, month, year, three years etc.) that generates multiple observations for each site. These observations are likely to be correlated over time because many of the unobserved 22 effects associated with a specific site will remain the same with time. Similarly, a spatial correlation might exist since neighboring sites typically have similar environmental and geographic characteristics and thereby form a cluster that has similar collision occurrence (Mountain et al., 1998; Shankar et al., 1998; Lord and Persaud, 2000; Ulfarsson and Shankar 2003; Guo et al., 2010). Spatial and/or temporal correlation can be accounted for through modifying the mean function. Temporal effects can be incorporated by adjusted the Poisson hierarchical model as θ i = µ i exp(u i ) exp(t i ) . (12) The temporal component ti accounts for possible time trends due to socioeconomic changes, traffic operation modifications, weather variations, etc. (Miranda-Moreno, 2006). Alternatively, spatial effects can be incorporated by adjusted the Poisson hierarchical model as θ i = µ i exp(u i ) exp(si ) . (13) The spatial component si suggests that sites that are closer to each other are likely to have common features affecting their collision occurrence. As noted by Miaou and Lord (2003), random variations across sites may be structured spatially due to the complexity of the traffic interaction around locations. Different formulations exist to capture both the spatial and temporal effects that exist in collision data. As a result, this issue is explored further in the thesis. 2.3.3 Random Parameters Models In the variable variance models, the focus was on studying the structure of the variance function. However, Lord and Park (2008) emphasized that transportation safety analysts should concentrate on the structure of the mean function. Ideally, a good structure along with a proper selection of the covariates for the mean function would simplify the structure of the variance 23 function or at least significantly minimize the magnitude of the variance (Miaou and Song, 2005; Mitra and Washington, 2007). A recent approach for modeling the mean function advocates the use of random parameters (Anastasopoulos and Mannering, 2009). In all traditional models only one regression equation was fit to the dataset. Using a random parameters model develops different regression equations for individual sites. This type of modeling technique has been considered by several researchers for its added flexibility and intuitive appeal (Shankar et al., 1998; Li et al., 2008; Milton et al., 2008). The model can be viewed as an extension of the random effects models since instead of varying only the intercept of a model, random parameters models allow each estimated parameter in the model to vary across each individual observation in the dataset. This model focuses on explaining part of the extra-variation through improvements in the mean function by accounting for the unobserved heterogeneity from one site to another. To be specific, note that Equation (5) for µ i can be written explicitly as follows ln(µ i ) = β 0 + β 1 X i1 + ...+ β p X ip , (14) where the X ij are covariates representing certain traffic and geometric characteristics, while the β ij are model parameters. The random parameters models can be represented by allowing all regression coefficients (not just the intercept) to vary randomly from one observation to another leading to ln(µ i ) = β i ,0 + β i ,1 X i1 + ...+ β i , p X ip , (15) where β i , j ~ N ( β j , σ 2j ) , j = 0,1, ..., p . (16) Several alternative distributions were considered in (16) but the normal distribution was found to provide the best statistical fit (Milton et al., 2008; Anastasopoulos and Mannering, 2009; Li et al., 24 2008). In practice, a random parameter β i, j is used in (15) whenever the posterior estimate σˆ j is significantly greater than 0, otherwise the parameter β j is fixed across segments (Anastasopoulos and Mannering, 2009). Random parameters models usually provide a statistically significant better fit than a model with traditional fixed parameters due to the fact that each observation having its own parameters. In certain cases, the n sites under consideration belong to K mutually exclusive clusters (e.g., districts, municipalities, zones, corridors, etc.). Suppose that the ith site belongs to cluster c( i ) ∈ { 1,2 ,..., K } . As a compromise between fitting the single Equation (14) and the n Equations (15), the cluster random parameters model can be represented by ln(µ i ) = β c (i ),0 + β c (i ),1 X i1 + ...+ β c (i ), p X ip , (17) where β c (i ), j ~ N ( β j , σ 2j ) , j = 0,1, ..., p . (18) The added flexibility provided by the cluster random parameters modeling allows for a more detailed inference by fitting a different regression curve for each cluster. This will allow for cluster- and overall-inference to be made. This topic will be addressed in more detail in subsequent sections. The applicability and practicality of the cluster random parameters models is demonstrated throughout this thesis. 2.3.4 Multivariate Models Collision data is typically available by severity (e.g., fatal, major injury, minor injury or property damage only), by the number of vehicles involved (e.g., single or multiple), and/or by the type of collision (e.g., angle, head-on, rear-end, sideswipe or pedestrian-involved), etc. However, most collisions have been mostly analyzed by modeling each category separately, without taking into account the correlations that probably exist among the different levels. These correlations may be caused by omitted variables, which can influence collision occurrence at all levels of 25 classification, or from ignoring shared information in unobserved error terms. Such univariate treatment of correlated counts as independent can lead to imprecise road safety analysis. Recently, multivariate SPFs have been proposed in the safety literature (Tunaru, 2002; Ma and Kockelman, 2006; Brijs et al., 2007; Park and Lord, 2007; Ma et al., 2008; Ye et al., 2009; Aguero-Valverde and Jovanis, 2009). These multivariate extensions are based on either multivariate Poisson (MVP) models (Tsionas, 2001; Karlis; 2003; Karlis and Meligkotsidou, 2005) or multivariate Poisson log-normal (MVPLN) models (Chib and Winklmann, 2001). The MVPLN regression is preferred to the MVP approach for the analysis of multivariate collision count data because (i) it accounts for over-dispersion (extra Poisson variation), which is often observed in collision data; and (ii) it allows for a fairly general correlation structure. For a set of data on road collisions at n locations, where the collisions at each location are classified into K categories, define the vector yi = ( yi 1 y i 2 ... yiK ) , where y ik denote ' the number of collisions at the ith location in category k . It is assumed that the y i are independently distributed and that the Poisson distribution of y ik , given θ ik , is Pr{Y ik = y ik θ ik} = e −θ ik θ iky ik yik ! , i = 1,2,..., n , k = 1,2,..., K . (19) To model extra variation, assume further that ln(θ ik ) = ln(µ ik ) + ε ik , ln(µ ik ) = β k 0 + β k1 X i1 + ... + β kp X ip , (20) and the ε ik denote multivariate normal errors distributed as ε i ~ N K (0, Σ ) , (21) where ε i1 ε i2 εi = , ... ε iK σ 11 σ 12 σ 21 σ 22 Σ = ... ... σ K1 σ K 2 ... σ 1K ... σ 2 K . ... ... ... σ KK 26 The importance of addressing the multivariate nature of collisions data is investigated in more details in this thesis particularly its consequence in terms of identifying hazardous ‘hotspots’ locations. 2.3.5 Finite Mixture/Markov Switching Models Finite mixture models assume that the data are generated from several distributions that are mixed together. Such an assumption implies that individual observations are generated from unknown number of subgroups. Similarly, Markov Switching models assume that a number of underlying distributions is used to generate the data and that each observation can switch among these distributions over time. The application of finite mixture models was studied by Park and Lord (2009) and Park et al. (2009). Markov Switching models were investigated by Malyshkina et al. (2009) and Malyshkina and Mannering (2010). Both finite mixture and Markov switching models offer considerable potential for providing important new insights into the analysis of collision data. Moreover, collision 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 (outliers) can have significant impact on the estimates of the regression parameters. Generally, it is recommended that outliers be identified, according to some statistical measures, and then removed from the data set, e.g., see Kulmala (1995) who used the leverage statistic and Sawalha and Sayed (2006) who used Cook’s distance for outlier detection. Few studies have discussed the effect of outliers and how to handle them in model development. Finite mixture 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 twocomponent mixture of gamma or normal priors (depending on the model's hierarchy), where it is assumed that most of the observations come from a basic distribution, whereas the remaining 27 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. The application of finite mixture models to identify and account for outliers in collision data is investigated in this thesis. 2.4 Safety Performance Functions (SPFs) In the Poisson hierarchal models the Poisson parameter θi is commonly specified as an exponential function of site-specific attributes multiplied by some random effect. Recall that µ i = exp( X i ' β ) . Such link functions are known as Safety Performance Functions (SPFs). For a long time the link functions were also known as Collision Prediction Models (CPMs). However, the term ‘CPMs’ seem to indicate that the models are used only for prediction purposes, which is a common misconception, since some models are developed for "prediction" while others are developed for "analysis". In collision data analysis one may study existing data to make fundamental inferences versus a study that seeks simply to predict. Not every model has the practical-oriented goal of prediction; some have the fundamental goal of providing new insights into the process which may not always lead to the best predicting model. Consequently, the term CPMs was dropped for the more general term SPFs. The preceding sections outlined the different modeling methods as well as some of the emerging modeling techniques in safety analysis. The next sections describe the specifications of SPFs when used for data and intervention analysis. 2.4.1 Safety Performance Functions for Collision Data Analysis Several approaches exist for estimating safety ranging from simply using collision rates to SPFs. Several researchers (Hauer, 1995) have shown that the relationship between collision frequency and exposure is frequently nonlinear which indicates that collision rates are not 28 appropriate representatives of safety. This finding has led most safety researchers to discard the use of collision rates as a measure of road safety and currently, safety performance functions constitute the primary tools for estimating road safety (Sawalha and Sayed, 2006). For the purpose of collision data analysis, a SPF is a mathematical model that is developed statistically to relate the expected collision frequency at a road location to its traffic and geometric characteristics. Collision data are usually available by collision severity, collision type, or other collision characteristics such as cause or occurrence time. Due to the lack of annual-based information on traffic flows, geometry modifications and other possible explanatory factors, collision data analysis is often limited to a cross-sectional data analysis. As a result, traffic variations at each location over time are usually unobserved. This could be attributed to the prohibitive costs of collecting relevant variables periodically (Miranda-Moreno, 2006). The analysis of collision data is done primarily to gain a better understanding of the factors that affect the likelihood of a road collision (Abdel-Aty et al., 2009; Haleem and AbdelAty, 2010). As mentioned before, the absence of detailed driving data has been partially responsible for our inability to identify the cause and effect relationships. To address this problem several researchers framed their analytical approaches to study the factors that affect the frequency of collisions in both geographical space (intersection or road segment) and specified time period (day, month, year) through the use of SPFs. A systematic statistical analysis of SPFs is typically used to develop collision frequency models that can explain the observed and random variations of collisions across sites. The mathematical form used for the SPFs should generally satisfy two conditions. First, it must yield logical results meaning that it should not lead to the prediction of a negative number of collisions. Second, the form should ensure a prediction of zero collision frequency for zero values of the exposure variables (Sawalha and Sayed, 2006). The most common functional form used for developing SPFs for road intersections is 29 ln(µ ) = ln(α 0) + α 1 ln(V mj ) + α 2 ln(V mi ) + ∑ Jj =1 γ j X j , (22) which is a special case of (14), were 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 α0, α1, α2, γ1, …, γJ are model parameters. For AADTbased models, Wang et al. (2006) explored different functional forms for modeling collisions at signalized intersections. The study recommended forms which were found to be a more representative indicator of traffic intensity at the intersections. Alternatively, the functional form (22), relating the frequency of collisions to the product of traffic flows entering the intersection, is the most commonly applied and investigated in the literature and it dates as far back as the 1950s (Tanner, 1953; Kulmala, 1995; Nicholson and Turner, 1996; Miaou, 1996; Miaou and Lord, 2003). Moreover, it is has been shown (Hauer et al., 1988) that a model that utilizes the product of the traffic flows, provides a better representation of the relationships between collisions and the traffic flows at intersections. This relationship was coagulated in a later publication where Hauer and Bamfo (1997) described two tools to determine the function that links the dependent variable to the explanatory variables. The corresponding form used for developing SPFs for road segments is ln(µ ) = ln(α 0) + α 1 ln( L) + α 2 ln(V ) + ∑ Jj=1 γ j X j , (23) which is a special case of (14), 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. 30 Canonical links in generalized linear models, such as the log link in NB and PLN regressions, are used to linearize the relationship between the predicted value and the predictors. They are popular among practitioners since 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). To guard against link misspecification, a Bayesian semi-parametric procedure was employed to estimate the link function, substantially reducing the risk of a mis-specified model (Shively et al., 2009). An alternative parametric approach is to embed the canonical (log) link into a wide class (generalized family) of links and use the full Bayes method for parameter estimation, performance evaluation and inference (El-Basyouny and Sayed, 2010a). Although the explanatory variables used in SPFs are supposed to be measured without errors, this may not necessarily hold true for some of these explanatory variables that are measured with uncertainty. Such errors-in-predictors can attenuate the predictors’ effects and increase dispersion (Chesher, 1991). For simple linear models measurement errors (ME) result in regression estimates that are biased downward, but the effect of unacknowledged ME in nonlinear models apparently resists simple characterization (Carroll et al., 1995). Traffic volume (or flow counts) is an example of predictors that are measured with uncertainty, which if not resolved can bias the regression estimates (Maher and Summersgill, 1996; Davis, 2000). This variable is perhaps the most commonly used in SPFs for both intersection and road segments. To account for ME, traffic volume time series data were used to develop SPFs where the underlying idea is to adjust the volume data for temporal correlation and measurement errors in order to improve the estimates of traffic flow (El-Basyouny and Sayed, 2010b). 31 2.4.2 Safety Performance Functions for Collision Intervention Analysis Time series and cross-section studies are two techniques that are frequently used to estimate the safety effect of specific countermeasures or location characteristics. The most common method to estimate the effectiveness of safety initiatives is a time series analysis, which is often referred to as “before-after” analysis. This approach attempts to measure the change in safety over time due to the implementation of a safety initiative. A controlled before-and-after (CBA) study is generally perceived by many to be an effective way to estimate the safety effect of changes in roadway characteristics. Alternatively, a cross-sectional analysis is sometimes used in the evaluation of road safety initiatives. A cross-sectional study compares the expected collision frequencies of a group of locations having a specific component of interest to the expected collision frequency of a group of locations with similar characteristics, yet they lack the presence of this specific component. Any differences in collision frequency between the two groups are attributed to the change in conditions, representing the safety effect. Cross-sectional studies are considered inferior to time series analysis (observational before-after studies) since no actual change has taken place. As a result, a cross-sectional analysis approach is usually used to develop an initial estimate of the safety effects. On the other hand, before-after safety evaluations base their results on actual changes that have occurred in one data set over a period of time extending from the before condition to the after condition (Ozbay et al., 2009). For the purpose of collision intervention analysis, a SPF is a mathematical model that is developed statistically to relate the collision occurrence on various road facilities as a function of time, treatment and interaction effects that acknowledges that safety treatment (intervention) effects do not occur instantaneously but are spread over future time periods. Several researchers have proposed the use of SPFs based on the intervention methodology to conduct CBA safety evaluations (Pawlovich et al., 2006; Li et al., 2008). In such studies, collision data is 32 usually available before and after the intervention for both a group of treated sites and a group of comparison sites. The comparison group of sites is used to correct for time trend effects, including the confounding factors of history and maturation. Comparison sites serve as a control group to minimize the unintended influence of other variables in the system. If the proposed countermeasure is actually a causal factor of improved safety, then logic dictates that the safety improvement should manifest itself more significantly in the treatment group than in the control group. To this end, the comparison sites are chosen based on their proximity to the treated sites so that they experience similar traffic and environmental conditions, but are not subject to any treatment. Care must be exercised in selecting the comparison group to ensure that the time periods for the treatment and comparison sites are similar so that the factors influencing safety are similar between the two groups (Sayed and de Leur, 2007). Li et al. (2008) considered various intervention forms and hierarchical models (PoissonGamma 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 intervention models to include multivariate dependent variables with multiple regression links and proposed an algorithm for the computation of the Odds Ratio (OR) to determine the effectiveness of the countermeasures. It should be noted that in all studies, linear slopes were assumed to represent the time and treatment effects across the treated and comparison sites. “Let Y it and µ it denote the observed and expected (predicted) collision counts at site i (i = 1, 2, …, n) during year t (t = 1,2, …, m).” 33 It is assumed that the data are available for a reasonable period of time before and after the intervention. The most common functional form to capture the gradual treatment effect using an intervention SPF is given by ln(µ it ) = α 0 + α 1 T i + α 2 t + α 3 (t − t 0i ) I it + α 4 T i t + α 5 T i (t − t 0i ) I it + γ 1 ln(V 1it ) + γ 2 ln(V 2it ) + γ 3 X 3i + ... + γ J X Ji , (24) where T i is a treatment indicator (equals 1 for treated sites, zero for comparison sites), t 0 i is the intervention year for the ith treated site and its matching comparison group, I it is a time indicator (equals 1 in the after period, zero in the before period), V 1it and V 2 it denote the annual average daily traffic (AADT) at the major and minor approaches, respectively, ( X 3i , ..., X Ji ) are additional covariates representing geometric and environmental site characteristics, (α 0 , ..., γ J ) are the regression coefficients. In Equation (24), the parameter α 1 represents the countermeasures effect (the difference in log collision count between treated and comparison sites), α 2 represents a linear time trend, α 3 represents the slope due to the intervention, α 4 and α 5 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 Equation (24) in order to accommodate such an effect leading to ln(µ it ) = α 0 + α 1 T i + α 2 t + α 3 (t − t 0i ) I it + α 4 T i t + α 5 T i (t − t 0i ) I it + α 6 T i I it , + γ 1 ln(V 1it ) + γ 2 ln(V 2it ) + γ 3 X 3i + ... + γ J X Ji . (25) 34 To determine treatment effectiveness, the Odds Ratio (OR) is used to assess the relative collision risk in a before-after framework involving treatment-comparison matching pairs. Let µ TBi and µ TAi denote the predicted collision counts for the ith treated site averaged over appropriate years during the before-after periods, respectively, and let µ CBi and µ CAi 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. The odds ratio for the ith treated site is given by ORi = µTAi µ CBi µTBi µ CAi . (26) Alternatively, following Park et al. (10), the ratio µ CA µ CB can be used to adjust the prediction for general trends between the before-after periods. Thus, the predicted collisions in the after period for the treatment group had the countermeasures not been applied is given by π TA = µ TB ( µ CA µ CB) . The index of effectiveness of the countermeasures is given by the ratio µ TA / π TA , which reduces to the odds ratio (26). The overall odds ratio can be computed from ln(OR) = (1 / n)∑in=1 ln(ORi ) . (27) The overall treatment effect is calculated from ( OR − 1 ) , while the overall percentage of reduction in predicted collision counts is given by ( 1 − OR ) × 100. 2.5 Parameter Estimation Methods There are several methods to calibrate the regression parameters of the SPFs. Perhaps the two most common estimation methods are the Empirical Bayes (EB) method using the Maximum Likelihood estimation technique and the fully Bayesian method. For many years, the preferred estimation technique was the EB method. The two most likely reasons for the popularity of the EB method are i) the technique uses the maximum likelihood estimation method which is 35 preprogrammed in many of the available statistical software and ii) the simplicity of conducting the posterior analysis since closed-form functions are available for most common distributions. However, this advantage is also a limitation since the EB method cannot be used when the likelihood function is difficult to characterize. On the other hand, recent years have witnessed enormous progress in numerical methods which, combined with the availability of free and easy to use software, have permitted the implementation of the full Bayes approach. This section discusses these two most common estimation methods. 2.5.1 Empirical Bayes (EB) Method Hauer (1997) developed a standard form of Bayesian statistical technique that is now widely applied to the analysis of road collision data (Heydecker and Wu, 2001). Bayesian methods treat unknown parameters such as predicted collision frequency as being random variables themselves having their own probability distributions. Once a parameter’s probability distribution has been determined, the parameter is estimated by the mean of its distribution. The probability distribution of an unknown parameter is obtained in two stages. The first stage involves determining a prior distribution before any data become available. The second stage involves implementing Bayes theorem to convert the prior distribution into a posterior distribution after the data become available. In this light, Bayes’ theorem provides a means to update the prior distribution in light of site-specific observations of collision occurrence resulting in a Bayesian posterior distribution. Bayesian methods account for the stochastic effects in collision data and enable analysts to combine evidence from observations of stochastic systems with that from other sources (such as models that describe standard operation) to provide estimates that are as accurate as possible given the information used (Heydecker and Wu, 2001). In the context of road safety, the Bayes' theorem is given by: p (θ | y, λ ) = p ( y | θ ) p (θ | λ ) / p ( y | λ ) , (28) where 36 • y is a random variable that represents the collision frequency at a given location during a specific time period; • p ( y | θ ) is the conditional probability of y given θ ; • θ is a random variable representing the mean collision frequency; • p (θ | λ ) is the prior probability of θ given ‘the second level’ parameter λ . It is "prior" in the sense that it does not take into account any information about y; • p ( y | λ ) is the conditional probability of y given λ and acts as a normalizing constant. • p (θ | y, λ ) is the conditional probability of θ given y and λ . It is also called the posterior probability because it is derived from or depends upon the specified value of y; Intuitively, Bayes' theorem in this form describes the way in which one's beliefs about ' θ ' are updated by having observed 'y'. In the Empirical Bayes approach, the regression parameters of the SPF are calibrated using the maximum likelihood technique or any other classical estimation technique using collision data (Hauer et al., 2002; Miaou and Lord, 2003; Miranda-Moreno, 2006). Since the second level parameter λ is unknown, it is typically estimated by maximizing the conditional distribution p ( y | λ ) with respect to λ (Carlin and Louis, 1996). A major drawback of the EB approach is that the estimate of λ , denoted by λˆ , is based on the data and once estimated, is taken as known and fixed. Thus, the inference about the first level parameter θ is based on the posterior distribution p (θ | y , λˆ ) and is conditional on the fixed estimate λˆ . This results in overconfident estimation for θ because the uncertainty in λ is not incorporated in the inference. However, the approach is widely recognized and applied in the literature for its ability to account for the regression-to-the-mean (RTM) bias and other confounding factors. The RTM is a statistical phenomena that refers to the tendency of extreme events to be followed by less extreme events even if no change has occurred in the underlying mechanism that generates the 37 process. Consequently, a site that experiences a high number of collisions during a certain period may experience few collisions during a later period without having been subject to any safety improvement measures (Sayed and de Leur, 2007). The EB approach is widely established and has been used extensively since its inception with the method of maximum likelihood being the most widely used to calibrate the regression models’ parameters. Maximum likelihood estimation is available in all commercial statistical packages, which also provide several measures that can be used to assess the goodness-of-fit of the developed SPF. Two commonly used measures are the Pearson χ2 statistic and the Scaled Deviance, which is defined as the likelihood ratio test statistic measuring twice the difference between the log likelihoods of the studied model and the full (saturated) model where each location has its own collision parameter. Both the scaled deviance and the Pearson χ2 have exact χ2 distributions for normal linear models, but are asymptotically χ2 distributed with k degrees of freedom for other distributions of the exponential family, where k is the sample size minus the number of estimated parameters (Sawalha and Sayed, 2006). 2.5.2 Full Bayes (FB) Method Sometimes, the classical estimation methods cannot be used or cannot be applied in a straightforward manner to estimate the parameters of a regression model. In such cases, the Full Bayes (FB) approach can be adopted for model calibration using Markov chain Monte Carlo (MCMC) simulation methods (Gilks et al., 1996). The main difference between the Empirical Bayes and Full Bayes approaches is in the way the prior parameters are handled. Recall in the EB approach the parameters are estimated using the maximum likelihood technique or any similar estimation techniques. In the FB approach the second level parameters are assigned prior distributions having third level parameters (hyper-parameters) based on prior idea or knowledge of the phenomenon being investigated. As a result, the FB approach overcomes one of the main shortcomings of the EB approach and is generally perceived as a more rigorous approach to collision modeling. 38 To account for all parameter uncertainties, the FB assumes that the second level parameters λ are not known a priori (even as point estimates) but are instead generated by a hyper-prior distribution p (λ | φ ) where typically φ is assigned a non-informative prior (Gelman et al., 2004). Inference about θ is then based on p (θ | y ) , which is obtained by integrating p (θ | y, λ ) with respect to the distribution of λ. Therefore, the FB approach is more flexible because it accounts for parameter uncertainties as well as spatial and/or temporal correlations, while modeling geographical and/or time dependence. To obtain the full Bayes estimates it is required to specify prior distributions for the parameters. Prior distributions are meant to reflect to some extent prior knowledge about the parameters of interest. If such prior information is available, it should be used to formulate the so-called informative priors. The specification of informative priors for generalized linear models was dealt with by Bedrick et al. (1996), who considered conditional means priors as well as data augmentation priors of the same form as the likelihood and showed that such priors result in tractable posteriors. A good discussion on the elicitation of priors in collision data analysis can be found in Schluter et al. (1997). In contrast, uninformative (vague) priors are usually used to reflect the lack of prior information (Miaou and Lord, 2003; Qin et al., 2005; Miaou and Song, 2005; Mitra and Washington, 2007; Lord et al., 2008; Li et al., 2008). It should be noted that introducing vague priors on all unknown parameters can be risky under some conditions such as the combination of low mean and small sample size (MirandaMoreno and Lord, 2007). In such cases, better results can be obtained by using semiinformative priors with a small mean and variance for the dispersion parameter. This is especially true for the negative binomial models. The PLN models seem to be more robust as the inverse-gamma prior assumed for σ u2 is a conjugate distribution of the normal distribution. The selection of priors is not entirely an analytical decision. Analysts should seek informative (or at least semi-informative) priors reflecting their previous knowledge about the 39 subject matter. If such information is absent, uninformative (vague) priors are usually used to reflect the lack of information. Posterior inference of the parameters in hierarchical Bayes models often involves working with high-dimensional integrals. This problem can be solved by generating a large number of samples from the posterior distribution using MCMC techniques. The MCMC method provides a unified framework, within which model identification and specification, parameter estimation, performance evaluation, inference, and communication of complex models can be conducted in a consistent and coherent manner (Miaou and Lord, 2003). From these samples, the posterior quantities of interest are computed for the model parameters. In fact, the rediscovery of MCMC methods by statisticians in the last 15 years have revolutionized the entire statistical field (Besag et al., 1995; Gilks et al., 1996; Robert and Casella, 1999; Roberts and Rosenthal, 1998). With the current computing power, sampling the posterior distributions using MCMC algorithms as needed in the full Bayes methods is easily achieved. Existing statistical programming software packages, such as WinBUGS (Lunn et al., 2000) and MLwiN (Yang et al., 1999) provide Gibbs and other MCMC sampling techniques which are needed to develop hierarchical Bayes models. The MCMC methods are used to repeatedly sample from the joint posterior distribution. The techniques generate sequences (chains) of 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 parameter 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 BrooksGelman-Rubin (BGR) statistic (Brooks and Gelman, 1998). A value under 1.2 of the BGR 40 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. Spiegelhalter et al. (2002) proposed the Deviance Information Criteria (DIC) as a measure of model complexity and fit. Let D denote the un-standardized deviance of the postulated model then DIC = D + p D ; p D = D − Dˆ , (29) ˆ is the point estimate obtained by substituting the where D is the posterior mean of D, D posterior means of the model’s parameters in D and p D 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 penalizes larger parameter models. According to Spiegelhalter (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. While the DIC is used for model comparisons, a posterior predictive approach (Gelman et al., 1996; Stern and Cressie, 2000; Li et al. 2008) can be used to assess the goodness-of-fit (adequacy) of the model with the lowest DIC. 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; the discrepancy cannot reasonably be explained by chance if the p-values are close to 0 or 1 (Gelman et al., 1996). 41 The replicates are best obtained simultaneously with model estimation in WinBUGS in order to account for all uncertainties in model parameters as reflected by the estimated distributions. The chi-square statistic is computed from χ 2 = ∑in=1 [yi − E (Y i )] / Var (Y i ) , 2 (30) where the yi denotes either the observed or the replicated collision frequencies. Several studies have compared the Full Bayes technique to the widely established empirical Bayes approach (Lan et al., 2009; Persaud et al., 2009) with the conclusion that the FB method allows additional flexibility in the development of SPFs which are capable of accounting for different variations (i.e., temporal, spatial, etc.) and properties (i.e., multivariate dependent and independent variables, errors-in-predictors, etc.) that exist in collision data. To benefit from the additional advantages of the FB approach, all the SPFs developed in this thesis were based on the full Bayesian estimation technique. 2.6 Summary This chapter provided an overview of the different collision modeling techniques that are commonly applied in the literature. Moreover, recent advances in collision modeling have been highlighted. These advances have facilitated the fitting of more complex SPFs that are used for both collision data analysis and collision intervention analysis. The development of these SPFs including: estimation and calibration, testing for goodness-of-fit, formulation, among many other things have been presented. However, there are many issues related to collision modeling that either deserve further attention or have not been dealt with adequately in the road safety literature. The next two parts of this thesis provides a more in-depth analysis of these issues. 42 PART I: SAFETY PERFORMANCE FUNCTIONS FOR DATA ANALYSIS 3 SPATIAL MODELS This chapter investigates the inclusion of spatial effects in the development of SPFs. 3.1 Background Despite the evident spatial nature of road collisions, relatively limited research has been conducted in road safety to account for spatial correlation. Several efforts focused on determining the practical consequences of not accounting for spatial effects in modeling accidents. Levine et al. (1995a) described spatial patterns in Honolulu motor vehicle collisions. They developed spatial software tools for describing the degree of spatial concentration. In a follow-up project (Levine et al., 1995b), the spatial relationship between activities which generate trips and motor vehicle collisions was examined and applied to the City and County of Honolulu. A spatial lag model was developed which examines the zonal relationship of motor vehicle collisions to population, employment and road characteristics. They argued that this method focuses attention on characteristics of neighborhoods and areas, and not just on the road system. Black and Thomas (1998) proposed a method to assess the extent to which the value of a variable on a given segment of a network influences values of that variable on contiguous segments. The method used network autocorrelation analysis, a network variant of spatial autocorrelation analysis, and Moran's I statistic to make the assessment. The study concluded that there was a positive and significant spatial correlation in the data. Nicholson (1999) advocated the need for spatial distributions of collisions to be analyzed. The author argued that current practice for assessing spatial distributions of collisions is insufficient. Certain statistical analysis techniques for spatial data (including quadrant, nearest-neighbor methods and K-function) were assessed and their implementations were 43 discussed. The nearest neighbor method was found to be the most powerful and robust technique. Miaou et al. (2003) indicated that mapping transforms spatial data into a visual form thereby enhancing the ability of users to observe, conceptualize, validate, and communicate information. The authors used a hierarchical Bayes model to build model-based risk maps for area-based road collisions. The research used county-level vehicle collision records and roadway data from Texas to illustrate the method. Conditional auto-regressive model (CAR) was used to model spatial correlation. The authors indicated that the models were useful tools for spatially correlated multivariate data to analyze all response variables simultaneously. MacNab (2004) used Bayesian spatial and ecological regression models to analyze small-area variation in the frequency and severity of collisions. The study serves to demonstrate how Bayesian modeling techniques could be implemented to assess potential risk factors measured at group (area) level. The study used hospital data for 83 local health areas in British Columbia (BC), Canada. The random spatial effects were included in the model and assumed to have a CAR distribution. The research proposed a unified modeling framework that enables thorough investigations into associations between injury rates and regional characteristics, residual variation and spatial autocorrelation. Miaou and Song (2005) investigated the impact of spatial effects on the overall model goodness of fit and site ranking performances. The research extended a univariate CAR spatial model to a multivariate setting following a Bayesian framework. The authors concluded that the inclusion of a spatial component in the collision prediction model significantly improved the overall goodness of fit performance of the model and affected the ranking results for the two data sets that were examined. Aguero-Valverde and Jovanis (2006) compared full Bayes (FB) hierarchical models (with spatial and temporal effects and space–time interactions) to traditional negative binomial estimates of annual county-level collision frequency. The study utilized injury and fatal collision 44 data for Pennsylvania between the years 1996-2000. Spatial correlation, time trend, and space– time interactions were found to be significant in the FB injury collision models. County-level FB models revealed the existence of spatial correlation in collision data and provided a mechanism to quantify, and reduce the effect of this correlation. The authors concluded by advocating the inclusion of spatial correlation in road segment and intersection-level collision models, where spatial correlation is likely to be even more pronounced. Wang and Abdel-Aty (2006) used the generalized estimating equations technique with the negative binomial link function for temporal and spatial analyses of rear-end collisions at signalized intersections. The intersections were grouped into clusters based on their spatial location, where those within a cluster were considered correlated while intersections from different clusters were considered independent. Three different correlation structures were examined where the correlation decreases as the gap between intersections increases. Their results demonstrated high spatial correlations between the intersections for rear-end collisions. Some intersection related variables were identified as significantly influencing rear-end collision occurrences at signalized intersections. Also, the introduction of spatial correlation has resulted in noticeable changes in the estimates of several regression coefficients, indicating possible model misspecification if the correlation is not accounted for. In a more recent publication, Aguero-Valverde and Jovanis (2008) explored the effect of spatial correlation in models of road collision frequency at the segment level. A FB hierarchical approach was used with CAR effects for the spatial correlation terms. The study used collision, traffic and roadway data from a rural county in Pennsylvania. The results indicated the importance of including spatial correlation in road collision models. The models with spatial correlation showed significantly better fit to the data than PLN. More importantly was the potential of spatial correlation to reduce the bias associated with model misspecification, as shown by the change on the estimate of the AADT coefficient. 45 Thus, the studies found in the literature dealing with spatial effects in collision modeling indicate the importance of including the spatial effects. To investigate the impact of incorporating spatial effects in SPFs, the Gaussian conditional auto-regressive (CAR), (Besag et al., 1991), and the multiple membership (MM), (Goldstein, 1995; Goldstein et al., 1998; Langford et al., 1999), models are compared to the traditional Poisson-Lognormal (PLN). The two spatial models pool information from neighboring sites to improve model estimation. It was argued in the literature (Lord and Persaud, 2000) that the aggregation of collisions over multiple years can be a useful tool for overcoming temporal correlation (at least partially). By the same logic, one is led to believe that clustering segments within the same corridor may reduce spatial correlation. To examine this potential, an additional variance component reflecting the variation among different corridors was incorporated in the MM model, which seems to be a logical choice for such an extension as each site is a member of multiple larger areas. The models were estimated in a Full Bayes (FB) context via MCMC simulation and were compared in terms of their goodness of fit and inference using a dataset composed of urban arterials in the city Vancouver, British Columbia. 3.2 The Models 3.2.1 Spatial Poisson Models The spatial NB (Poisson-Gamma) models are given by equations (1), (13), (14) and (6), whereas the spatial PLN (Poisson-Lognormal) models are given by equations (1), (13), (14) and (9). 3.2.2 Gaussian Conditional Auto-Regressive (CAR) Models The Gaussian CAR models (Besag et al., 1991) are among the most commonly used for modeling spatial effects. Recall Equation (13) and let ni , C ( i ) and s −i represent the number of 46 neighbors of site i, the set of neighbors of site i and the set of all spatial effects except s i , respectively. The Gaussian CAR model is given by 2 si | s −i ~ Normal ( s i , σ s / ni ) , s i = ∑ j∈C ( i ) s j / ni (31) Equation (31) is based on an adjacency-based proximity measure: wij = 1 if sites i and j are neighboring sites and wij = 0 otherwise. The conditional mean is the mean of adjacent special effects, while the conditional variance is inversely proportional to the number of neighbors. 3.2.3 Multiple Membership (MM) Models An alternative way to account for spatial correlation is to use MM models (Goldstein, 1995; Goldstein et al., 1998; Langford et al., 1999) where each site is considered a member of a higher level unit which contains its nearest neighbors. For a first-order spatial autocorrelation model, let u i and v i represent the random effects of site i and its effects on its neighbors, respectively. The spatial effect of site i is given by si = ∑ j∈C ( i ) v j / ni . To model a correlation between u i and vi , it is assumed that (u i , vi ) ′ ~ N (0, Σ ) , σ u2 ρ σ uσ v . Σ = σ v2 ρ σ uσ v (32) The MM models allow a more direct interpretation of the model parameters, as σ u2 and σ v2 represent marginal variances (as opposed to the conditional variance used in CAR models). Furthermore, the MM models provide an additional parameter ρ measuring the strength of the association between the unstructured and structured (spatial) effects. 3.2.4 Extended MM (EMM) Models Typically, the n sites under consideration belong to K mutually exclusive corridors. In such cases, an additional variance component can be included in the MM model to allow for the 47 possibility that different corridors have different collisions risks because traffic, geometric and environmental conditions vary across corridors. Suppose that the ith site belongs to corridor c( i ) ∈ { 1,2 ,..., K } . The extended MM model is given by θ i = µ i exp(u i ) exp(si ) exp(wc ( i )) , (33) where wc ( i ) ~ Normal (0, σ c2) and σ c2 denotes the additional variance component representing the variation among different corridors. 3.3 Spatial 'Intraclass' Correlation The impact of spatial correlation is assessed by computing the proportion of total variation that is due to spatial variation ψ s = var( s ) /(var( s ) + σ u2 ) , (34) where var( s ) represents the marginal variance of s . For CAR models, var( s ) can be estimated directly from the posterior distribution of s . For MM models, var( s ) can be estimated via var( s ) = σ v2 / nh , where n h represents the harmonic mean of the ni (Langford et al., 1999). For EMM, Equation (34) becomes ψ s = var(s ) /(var(s ) + σ u2 + σ c2) , while the corresponding proportion due to corridor variation is given by ψ c = σ c2 /(var(s ) + σ u2 + σ c2) . 3.4 (35) Prior Distributions The most commonly used priors are diffused normal distributions (with zero mean and large variance) for the regression parameters, Gamma (ε , ε ) or Gamma (1, ε ) for k (under the negative binomial model) and σ u−2 (under the PLN model), where ε is a small number, e.g., 0.01 or 0.001. An uninformative inverse-gamma prior is assumed for σ c2 as well. 48 For the Gaussian CAR model, it was shown (Breslow and Clayton, 1993) that Equation (31) is equivalent to f ( s1,..., s n | σ 2s ) ∝ σ −s n exp{− ∑ ni s i ( s i − s i ) / 2 σ 2s} , which provides the correct likelihood function for σ 2s . The term contributed by each site is li = ni s i ( s i − s i ) . A proper prior Gamma(1 + ∑ li / 2,1 + n / 2) is assumed for σ −s 2 . It was also shown (Bernardinelli et al., 1995) that the marginal variance σ 2sm of the spatial effects is approximated by σ 2sm ≈ 2 σ 2s / n , where n = ∑ ni / n . Aguero-Valverde and Jovanis (2008) used a “fair” prior to specify the relationship between the unstructured and structured spatial effects. However, in this study it is assumed that σ u2 = c 2 σ 2sm , where c is a discrete random variable having a discrete uniform prior over the 19 values {0.1, 0.2, … 0.9, 1, 2, …, 10} (Congdon, 2006). For the MM models, a Wishart ( P, r ) prior is assumed for Σ −1 , where P and r represent the prior guess at the order of magnitude of the precision matrix Σ −1 and the degrees of freedom, respectively. To represent vague prior knowledge, r is usually chosen to be as small as possible, i.e., r = rank (Σ ) = 2 (Spiegelhalter et al., 1996; Tunaru, 2002; Spiegelhalter et al. 2005) whereas P is the 2 x 2 identity matrix. 3.5 The Data A total of n=281 urban road segments in the city of Vancouver were investigated for the purpose of developing SPFs relating the safety of these road segments to their traffic and geometric characteristics (Sawalha and Sayed, 2001). The data were obtained from the urban city of Vancouver and covered the period from 1994 to 1996. A geographic map was obtained from which the neighboring segments and the number of neighbors (ni) were determined (Figure 3.1). Various neighboring structures were considered by Aguero-Valverde and Jovanis (2008). Based on their results as well as others reported in the literature, only first-order neighbors were considered to define the neighboring structure. First-order neighbors included all segments having direct connection with the segment under consideration. Also, it turned out that the 281 49 segments belong to K=37 corridors. Table 3.1 shows a list of the explanatory variables with their corresponding abbreviations and units. Table 3.1 Model covariates and their corresponding description Covariates Section length Annual average daily traffic Unsignalized intersection density Acronym L AADT UNSD Symbol L V X1 Undivided cross section IUND X2 Business land use IBUS X3 Between signal number of lanes Number of collisions NL Y X4 Y Units Kilometers Vehicles per day Intersections per kilometer Indicator variable yes = 1, no = 0 Indicator variable yes = 1, no = 0 The road segment length ( Li ) was incorporated as an offset in the following functional form, which captures the nonlinear relationship between the mean number of collisions ( µ i ) and the traffic flow volume ( V i ) (Sawalha and Sayed, 2001; Qin et al., 2004; Aguero-Valverde and Jovanis, 2008) ln(µ i ) = ln( Li ) + ln(α 0) + α 1 ln(V i ) + ∑ Jj =1 γ j X ij , (36) where X ij are covariates representing certain traffic and geometric characteristics, and α 0 , α 1 , γ 1 ,..., γ J are model parameters. Note that other functional forms are available in the literature (Miaou and Lord, 2003; Qin et al., 2004). However, since the focus of this research component is to demonstrate the effects of including spatial effects on the development of SPFs, only the above functional form is investigated. The form allows for no collisions if the exposure is zero and it has been investigated extensively in the literature since the early 1950s (Lord, 2000). Average values for the traffic volumes (over the 3 years period) were used to build SPFs for predicting the total (aggregate) number of collisions using Equation (36) with J=4. The aggregation is justified on several grounds. For instance, an aggregate SPF was found to 50 perform well compared with SPFs developed to handle temporal correlation (Lord and Persaud, 2000). Also, the aggregation of collisions 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 3.2. Table 3.2 Variable Statistical summary of data set (n = 281 road segments) Minimum Maximum L 0.11 V 4236 X1 0 X2 0 X3 0 X4 2 a ni 1 b Corridor size 1 a The harmonic mean is 4.24 ( n h ). b The number of road segments within the corridor. 3.61 62931 21.16 1 1 7 7 15 Mean Standard Deviation 0.79 26763 6.94 0.78 0.31 4.00 4.86 7.60 0.43 13794 3.29 0.42 0.46 1.51 1.32 4.14 To investigate the outcomes of incorporating spatial effects in SPFs, the CAR and MM models were compared to PLN. Furthermore, EMM was also fitted to assess the effects of clustering segments into corridors on spatial correlation. 51 Figure 3.1 GVRD study area. 52 3.6 Results and Discussion The posterior summaries in Table 3 were obtained from WinBUGS using two chains with 40000 iterations each, 20000 of which were excluded as a burn-in sample. 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. Table 3.3 summarizes the parameter estimates and their associated statistics for the PLN, CAR, MM and EMM models. The table shows that the parameter estimates were significant as the 95% credible intervals were bounded away from zero (except for IUND under CAR and MM). The regression coefficients were all positive, indicating that factors such as AADT, unsignalized intersections, undivided cross sections, business land use and number of lanes are positively associated with the number of collisions. Table 3.3 Parameter estimates for the PLN, CAR, MM and EMM models PLN CAR MM EMM Est. Sd Est. Sd Est. Sd Est. Sd Intercept V X1 X2 X3 X4 -2.854 0.548 0.078 0.260 0.312 0.148 0.821 0.087 0.011 0.084 0.083 0.034 -1.247 0.414 0.071 0.141 0.305 0.119 0.981 0.105 0.011 0.096 0.093 0.040 -1.696 0.337 0.083 0.123 0.349 0.135 1.255 0.136 0.015 0.122 0.137 0.053 -6.338 0.514 0.066 0.238 0.394 0.105 0.958 0.099 0.011 0.104 0.096 0.045 σ u2 σ v2 ρ 0.292 0.029 0.044 0.009 2.220 0.325 0.164 0.026 17.460 1.921 0.363 0.126 0.965 0.008 0.337 0.135 17.600 4.350 0.005 0.002 0.985 0.004 c σ σ c2 ψs ψc 2 s DIC 0.307 0.026 1.116 0.106 0.876 2101 2093 0.024 0.650 2105 0.030 2091 53 The estimates of σ u2 were significant under all four models, demonstrating the presence of over-dispersion in the data. The proportions of total variability due to spatial variation were 0.876, 0.650 and 0.005 for CAR, MM and EMM, respectively. Under the MM model, the estimate of ρ was 0.965, indicating that road segments with high extra-Poisson variation tend to have high spatial variation as well and vice-versa.. Under EMM, the contribution of corridors to total variability was 0.985 and the estimate of ρ was 0.337, barely significant. The DIC statistics were 2101, 2093, 2105 and 2091, under PLN, CAR, MM and EMM, respectively. The best fit model was EMM followed by CAR. The incorporation of the spatial correlation affected the estimation of parameters. The indicator variable for the presence of undivided cross-sections (IUND) was significant under PLN. After the inclusion of the spatial effects the variable became not significant. The spatial effects associated with the CAR and MM models seem to be related to the presence of undivided cross-sections. Such spatial multicollinearity appears to negate the need to include IUND in the SPF. Also, the regression coefficient of traffic volume (AADT) was reduced under the spatial models; more so under MM and CAR than under EMM, but the respective 95% credible intervals still overlapped. Table 3.3 reveals a significant negative impact of traffic volume on the safety of urban road segments in the city of Vancouver. The regression coefficient of AADT was 0.548 under PLN, which is smaller (albeit not significantly) than the estimate obtained by Aguero-Valverde and Jovanis (2008) for the rural segments in Central Pennsylvania. Urban segments differ from rural highways in the way segmentation can take place. On urban roads, signalized intersections are natural delineators of road sections whereas on rural roads the definition of segmentation is somehow ambiguous (Sawalha and Sayed 2001). The effect of AADT on the number of collisions was smaller under EMM and CAR (it was reduced further under MM) indicating a potential of reducing the bias associated with model 54 misspecification, due to the omission of spatial variables. While this result agrees with that of Aguero-Valverde and Jovanis (2008), the evidence of misspecification in the PLN model was not conclusive in both studies as the respective 95% credible intervals overlapped one another. The safety of Vancouver urban roads was also negatively impacted by business land use, the presence of undivided cross sections, the number of lanes between signals and the density of unsignalized intersections, as these covariates were positively associated with the number of collisions. These findings conform to those of the literature as reported by Shankar et al. (1997) and Park and Lord (2007). It should be noted that the main purpose of this research component is to show the importance of including spatial effects in Safety Performance Functions. A more thorough discussion on the affects of including different variables on the safety of urban arterials is found in Sawalha and Sayed (2001). The estimates of the parameters representing unstructured and structured (spatial) effects were significant under all four models. Furthermore, a high portion of the total variability (significantly larger than 50%) was explained by the spatial correlation under CAR and MM, but not under EMM where the corridor effects were overwhelming (ψ c = 0.985 ). There were considerable changes in the estimates of σ u2 , σ v2 and ρ when σ c2 was added to MM. The overlooked corridor variation in MM was reflected partly in σ u2 and mostly in σ v2 giving rise to a large spatial contribution (ψ s = 0.650 ). Under EMM, the estimates of σ u2 and σ v2 have been drastically decreased and the contribution due to spatial variation has subsequently become trivial (ψ s = 0.005 ), but the spatial effects were significant nevertheless. Further, the estimate of ρ , which measures now the partial association between the unstructured and structured spatial effects adjusted for corridor variation, was a modest 0.337. 55 According to the DIC guidelines, EMM provided the best fit. It has outperformed PLN and MM, while fitting the data a little better than CAR. CAR was the second best since it has outperformed PLN and MM. It should be mentioned that CAR was compared to the Poisson-Gamma (negative binomial) model using Equation (6) instead of Equation (9) and the results were found to be similar. A proper ICAR model (Miaou and Song, 2005; Congdon, 2006) was also fit and found to produce similar results to the CAR model in Equation (31). 3.7 Summary In this chapter, the traditional PLN was compared with two spatial models (CAR & MM) and an extended MM model (EMM) using a data set for urban arterials in the city of Vancouver, BC. The CAR model has (i) outperformed MM and PLN; (ii) explained a high proportion of spatial variability; and (iii) reduced the potential bias associated with model misspecification, due to the omission of spatial variables. The EMM model provided the best fit. It showed considerable improvement over MM and PLN, while fitting the data a little better than CAR. Although the spatial effects remained significant, they have been alleviated to a great extent by accounting for corridor variation. The results provide some strong evidence for incorporating spatial effects in SPFs. Apart from the improvement in goodness of fit, the inclusion of such effects could explain enough variation that some of the model covariates would be rendered non-significant and thereby the inference would be affected as well. Also, the current results reinforce the findings of other studies in the literature which suggested that spatial models have the potential for reducing the bias, in estimating the regression coefficients, resulting from the omission of spatial variables. However, the evidence was inconclusive and further research involving larger datasets and/or other types of road segments is needed. 56 It might be argued that the significance of spatial effects is simply an artifact of omitting important variables and/or inefficient determination of homogenous road segments. Thus, with appropriate definition/selection of road segments and corridors along with proper selection of (and collection of data on) pertinent covariates, the spatial effects would diminish or go away completely. While this may be possible in an ideal world, it is unlikely to happen in most applications because road segments are mainly determined by geographic/administrative considerations and data availability is usually limited. For instance, data on such characteristics as pavement surface, roadway alignment, speed limits, percentage of trucks, weather conditions, visibility, wind, etc., are rarely observed or measured despite their relevance to collision modeling. Accounting for spatial variation in the development of SPFs helps in overcoming such less than ideal circumstances. It would be interesting to investigate other techniques to account for the random spatial effects. While the CAR distribution is the most commonly used in disease mapping and road collisions, other techniques are available in the literature that warrants attention. Simultaneous autoregressive (SAR) and moving average (MA) are among the many alternative techniques used to model spatial correlation (Congdon, 2006). This research could be extended to compare other forms of incorporating the spatial effects. A distance-based weighting scheme could be utilized instead of neighbor proximity. A further extension to this work would be the development of a safety performance function that accounts for the spatial correlations between road segments and intersections. 57 4 RANDOM PARAMETERS MODELS This chapter investigates the use of random parameters in the development of SPFs. 4.1 Background In contrast to the traditional SPF, which fits one regression model to the dataset, the random parameters approach develops different regression models for each individual site in the data set. Several researchers have demonstrated the applicability of the random-parameters approach to explicitly account for the variations in the effects of variables across observations (Shankar et al., 1998; Milton et al., 2008; Li et al., 2008; Anastasopoulos and Mannering, 2009). Constraining the parameters to be constant when they actually vary across observations could lead to inconsistent and bias parameter estimates. Given the potential heterogeneity in collisionfrequency data, the random-parameters approaches may be appropriate in some cases. However, using the random parameters approach to develop a different regression model for each individual site in the data set may seem to be a case of over-fitting, which might affect the model’s generality. It is usually tempting to include many variables/parameters in a model in an effort to make it fit the observed data as closely as desired. However, the question that must be asked is whether a model having a perfect or extremely close fit to a particular set of observations can produce appropriate predictions when applied to a different set of locations. To answer this question, it should be noted that the random parameters approach permits both location-level inference and overall-level inference. Thus, while the overall-level inference might be not as much vulnerable, the location-level regression equations can be unstable and may perform poorly when applied to a new sample drawn from the same population. For the study of new locations, using more variables/parameters is not always better, and it may actually be worse (Sawalha and Sayed, 2006). 58 Therefore, to benefit from the advantages of using the random parameters approach and avoid over-fitting, it is proposed to cluster the road entities into rather homogeneous groups (e.g., districts, municipalities, zones, etc.) and fit a different regression curve for each group rather than for each individual site. Clustering has been show to improve model fit and predictive performance (Huang and Abdel-Aty, 2010). In this chapter, 392 road segments were clustered into 58 corridors. To investigate the consequences of incorporating corridor effects in SPFs, two models that account for the corridor variation through the mean and variance of an extended PLN model are compared to the traditional PLN model. Random parameters PLN models were used to account for the corridor variation through the mean by fitting a different regression curve for each corridor. This model focuses on explaining part of the extra-variation through improvements in the mean function. Alternatively, a spatial PLN model was used to account for the corridor variation through the variance by using an additional variance component. This spatial PLN model is somewhat similar to the EMM model proposed in chapter 3, which has been found to alleviate both spatial and extra Poisson variations. The model focuses on explaining part of the extra-variation through improvements in the variance function. The models were estimated in a Full Bayes (FB) context via Markov MCMC simulation and were compared in terms of their goodness of fit and inference using a data set composed of urban arterials in the city Vancouver, British Columbia. 4.2 The Models 4.2.1 Poisson-Lognormal with Random Corridor Effects (PLN-VC) Typically, the n segments under consideration belong to K mutually exclusive corridors. In such cases, an additional variance component can be included in the model to allow for the possibility that different corridors have different collisions risks because traffic, geometric and 59 environmental conditions vary across corridors. Suppose that the ith segment belongs to corridor c( i ) ∈ { 1,2 ,..., K } . This leads to the extended model PLN-VC based on (1), (14), (9) and θ i = µ i exp(ui ) exp(wc (i )) , (37) where wc ( i ) ~ N (0, σ c2) and σ c2 denotes the additional variance component representing the variation among different corridors. Under PLN-VC E (Y i ) = µ i exp(0.5[σ u2 + σ c2]) , Var( Y i ) = E( Y i ) + [ E( Y i )] 2 (exp( σ u2 + σ c2 ) − 1 ) . Let µ i+ = µ i exp(wc ( i )) and note that (38) ln(µ i+ ) = β c ( i ), 0 + β 1 X i1 + ...+ β p X ip , where β c ( i ),0 = β 0 + wc (i ) . The PLN-VC model is equivalent to a PLN with a random intercept that varies across corridors. 4.2.2 Poisson-Lognormal with Random Parameters among Pairs (PLN-RP) Alternatively, the corridor variation can be represented by allowing all regression coefficients (not just the intercept) to vary randomly from one corridor to another leading to the second extended model PLN-RP based on (1), (5), (17), (18), and (9). Under PLN-RP E (Y i ) = µ *i exp(0.5σ i2) , Var (Y i ) = E (Y i ) + [ E (Y i )]2 (exp(σ i2) − 1) . (39) where ln(µ *i ) = β 0 + β 1 X i1 + ...+ β p X ip , σ i2 =σ 02 + σ u2 + X i21σ 12 + ... + X ip2 σ 2p . (40) Although PLN is nested within PLN-VC, which in turn is nested within PLN-RP, the models have quite different characteristics. The PLN and PLN-VC models have the same form (14) for their mean functions and have fixed dispersion parameters σ u2 and σ u2 + σ c2 , respectively. In contrast, PLN-RP has a variable dispersion function (40), involving the covariates X ij , and permits both corridor-level inference, using (17), and overall-level inference, using (40). 60 4.3 Prior Distributions A diffused normal distribution (with zero mean and variance = 1002 ) was used for the regression parameters, whereas a Gamma (0.001, 0.001) was used for σ u−2 , σ c−2 and σ −j 2 . 4.4 The Data A total of n=392 urban road segments in the city of Vancouver were investigated for the purpose of developing Safety Performance Functions relating the safety of these road segments to their traffic and geometric characteristics (Sawalha and Sayed, 2001). The data were obtained from the urban city of Vancouver and covered the period from 1994 to 1996. The 392 segments belong to K=58 corridors. Table 4.1 shows a list of the explanatory variables with their corresponding abbreviations and units. Table 4.1 Model covariates and their corresponding description Covariates Segment length Annual average daily traffic Crosswalks density Acronym L AADT CROD Symbol L V X3 Business land use IBUS X4 Unsignalized intersection density Between signal number of lanes Number of collisions UNSD NL Y X5 X6 Y Table 4.2 Units Kilometers (km) Vehicles per day (veh/day) Crosswalks per km Indicator variable yes = 1, no = 0 Intersections per km Collisions per 3-years Statistical summary of data set (n = 392 road segments) Maximum Mean Standard Deviation L 0.11 3.61 V 4232 62931 X3 0 10 X4 0 1 X5 0 21.16 X6 2 7 Y 1 311 a Corridor size 1 17 a The number of road segments within the corridor. 0.82 24412 1.94 5.90 3.89 51.45 6.76 0.40 13151 2.11 3.40 1.36 50.12 3.83 Variable Minimum 61 In this application, Equations (14) and (17) are used with p=6, X1=ln(L), X2=ln(V) and X3X6 as described in Table 4.1. Although other functional forms are available in the literature (Miaou and Lord, 2003; Qin et al., 2004), only the above functional form is investigated to demonstrate the effect of including corridor variation on the development of SPFs. A statistical summary of the data is shown in Table 4.2. 4.5 Results and Discussion The posterior summaries in Table 4.3 were obtained via WinBUGS using two chains with 20000 iterations each, 10000 of which were excluded as a burn-in sample. 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. Table 4.3 summarizes the parameter estimates and their 95% credible intervals for the PLN, PLN-VC and PLN-RP models. The table shows that the parameter estimates are significant as the 95% credible intervals were bounded away from zero, except for X6 (NL) under PLN-RP. Apart from the intercepts, the regression coefficients are all positive, indicating that factors such as segment length, AADT, crosswalks density, business land use, unsignalized intersection density and number of lanes are positively associated with the number of collisions. The incorporation of corridor variation affected the estimation of parameters. In particular, the variable X6 representing the number of lanes (NL) is significant under PLN and PLN-VC, but is not significant under PLN-RP. In general, the regression coefficients have changed under the three models but the respective 95% credible intervals still overlapped. The estimates of σ u2 are significant under all three models, demonstrating the presence of over-dispersion in the data. However, the PLN estimate of 0.388 has been reduced by 29% and 44% under PLN-VC and PLN-RP, respectively. As expected, accounting for corridor variation reduces the estimates of extra Poisson variation. 62 The DIC statistics were 2840, 2832 and 2829 under the PLN, PLN-VC and PLN-RP models, respectively. Thus, the two extended PLN models fitted the data better than PLN, with the random corridor parameters PLN-RP model providing a significantly better fit than PLN, according to the DIC guidelines. Under PLN-RP, the observed value of chi-square was 290.8. To assess goodness-of-fit, the distribution of the chi-square discrepancy measure in replicated datasets was generated. The observed value of chi-square was located near the center of the replicated distribution, with an associated p-value of 0.485. As a result, PLN-RP seems to perform well in terms of accommodating the variation in collision frequency across segments. In contrast, the observed value of chi-square under PLN-VC was larger, 347.9, with a smaller p-value of 0.156. Table 4.3 Parameter estimates and 95% credible intervals for all models PLN PLN-VC PLN-RP -3.177 (-4.395, -2.055) -2.777 (-4.062, -1.496) -3.090 (-4.681, -1.587) 0.683 (0.542, 0.836) 0.711 (0.555, 0.871) 0.712 (0.558, 0.875) 0.561 (0.435, 0.698) 0.539 (0.400, 0.681) 0.571 (0.413, 0.738) 0.115 (0.075, 0.149) 0.112 (0.084, 0.143) 0.123 (0.081, 0.163) 0.249 (0.107, 0.397) 0.290 (0.113, 0.456) 0.286 (0.105, 0.483) 0.096 (0.075, 0.115) 0.081 (0.060, 0.102) 0.098 (0.066, 0.133) 0.118 (0.051, 0.186) 0.089 (0.023, 0.163) 0.063 (-0.019, 0.151) σ σ c2 σ0 σ1 σ2 σ3 σ4 σ5 σ6 0.338 (0.286, 0.400) 0.240 (0.199, 0.288) 0.188 (0.150, 0.233) DIC 2840 β0 β1 β2 β3 β4 β5 β6 2 u 0.123 (0.064, 0.211) 2832 2.294 (1.129, 4.404) 0.163 (0.067, 0.308) 0.214 (0.081, 0.436) 0.074 (0.047, 0.116) 0.224 (0.078, 0.445) 0.077 (0.049, 0.115) 0.135 (0.069, 0.224) 2829 63 The results of PLN-RP in Table 4.3 indicate that random parameters were associated with the intercept as well as all six covariates. Since the PLN-RP intercept is normally distributed with mean −3.090 and standard deviation 2.294, then it is expected to be less than the PLN fixed intercept of -3.177 on 48.5% of the segments and greater than the PLN intercept on 51.5% of the segments. This variability is capturing the unobserved heterogeneity across corridors. The road segment length also resulted in a random parameter that is normally distributed, with a mean 0.712 and standard deviation 0.163. Thus, for almost all corridors, collision counts are expected to increase with segment length although by varying magnitudes. A similar result was obtained for AADT, where collision counts are expected to increase with AADT, for the vast majority (99.6%) of the corridors, and are expected to decrease for only a small proportion (0.4%) of the corridors. As noted by Anastasopoulos and Mannering (2009), this AADT finding is likely picking up a complex interaction among traffic volume, driver behavior and collision frequency. It may be capturing, among other factors, the response and adaptation of drivers to various levels of traffic volume. The estimates of the distributional parameters for CROD, IBUS and UNSD indicate that collision frequency is expected to increase with crosswalks density, business land use and the density of unsignalized intersections on the majority of the corridors; the percentages were 95%, 89.9% and 89.9%, respectively. On the remaining corridors, collision frequency is expected to decrease; reflecting heterogeneity across corridors. As mentioned above, the number of lanes was found to be significant under PLN and PLN-VC, but not under PLN-RP. This finding maybe due to the fact that the corridor effects associated with traffic volumes, crosswalks density and business land use seem to be related to the number of lanes. Such multicollinearity appears to negate the need to include NL in PLNRP. 64 In addition to the overall estimates of β = ( β 0 , β 1 ,..., β 6 ) , provided in Table 4.3, there were 58 estimates of β corresponding to the 58 corridors. Such estimates can be used for corridor-level inference using (17). The nature of these regression coefficients is illustrated in Figure 4.1. Figure 4.1 Random regression coefficients (a) Scatter plot of β5 vs. average UNSD per corridor. (b) Bar chart of β 6 by corridor. Number of Lanes betw een Signals (NL) Unsignalized Intersections Density (UNSD) 0.30 0.25 0.25 0.20 0.15 0.10 0.15 β6 β5 0.20 0.05 0.10 0.00 0.05 -0.05 1 0.00 -0.10 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Average UNSD per Corridor (a) 5 9 13 17 21 25 29 33 37 41 45 49 53 57 -0.15 Corridor (b) For UNSD, the average values were computed for each of the 58 corridors and the 58 regression coefficients β 5 are plotted against these averages in Figure 4.1(a), which shows that: (i) the corridors’ estimates of β 5 varied between 0.02 and 0.25; (ii) the corridors varied in terms of UNSD as the average per corridor ranged from 0.6 to 13.8; and (iii) there is a tendency for the corridors’ regression coefficients to decrease as the average UNSD increases. The tradeoff suggested by (iii) is logical in terms of goodness-of-fit. The same pattern has occurred for β 1 ,..., β 4 . However, while the correlation coefficient between β 5 and average UNSD per corridor was -0.44, which is significant at p=0.001, the corresponding correlations for the other regression coefficients were weaker and not significant. Figure 4.1(b) displays the regression coefficients β 6 of the number of lanes by corridor. Not only there were 8 corridors with negative estimates, but also there were considerable variability as the coefficients ranged from -0.120 to 0.215. 65 4.6 Summary In this chapter, the traditional PLN model was compared with two extended PLN models using a data set for 392 urban arterials in the city of Vancouver, BC, that were clustered into 58 corridors. To assess the effects of incorporating corridor variation with alternate specifications, two models that account for corridor variation through the variance function (PLN-VC) and the mean function (PLN-RP) were considered. A variety of covariates representing segment length, AADT, crosswalks density, business land use, unsignalized intersection density and the number of lanes between signals were found to significantly influence collision frequencies. The PLN-RP model provided the best fit, representing a considerable improvement over PLN, while fitting the data a little better than PLN-VC. Under PLN-RP, the covariates resulted in random parameters and thereby their effects on collision frequency were found to vary significantly across corridors. The results provide strong evidence for the benefit of clustering road segments into homogeneous groups (e.g., corridors) and incorporating random corridor parameters in SPFs. Apart from the improvement in the goodness of fit, such an approach can be used to gain new insights into how the covariates affect collision frequencies, and to account for heterogeneity due to unobserved road geometrics, traffic characteristics, environmental factors and driver behavior. The inclusion of corridor effects in the mean function explained enough variation that some of the model covariates would be rendered non-significant and thereby the inference would be affected. This research could be extended to compare other ways of clustering road segments, e.g., districts, municipalities, zones, etc. 66 5 MULTIVARIATE MODELS This chapter investigates the use of multivariate dependent variables with multiple regression links in the development of SPFs. 5.1 Background Multivariate extensions of the simple Poisson distribution to account for the correlations among different count processes (such as counts of different collision severities) have been proposed in the literature. For instance, a model that assumes the same nonnegative covariance term for all pairs has been considered by Tsionas (2001) and Karlis (2003). Ma and Kockelman (2006) adapted a multivariate Poisson (MVP) regression approach to assess the effects of different covariates on collision counts at different severity levels. The assumption of equal covariance terms has been relaxed by Karlis and Meligkotsidou (2005), but the restriction to nonnegative covariance terms was still maintained. Brijs et al. (2007) used a similar non-regression model for ranking hazardous sites in Belgium. Multivariate zero-inflated Poisson (MVZIP) models were also considered (Li et al., 1999). Since the MVP regression models do not allow for over-dispersion, a multivariate negative binomial (MVNB) model was used to investigate the simultaneity of fatality and injury collision outcomes (Ladron de Guevara and Washington, 2004). A major drawback of the MVNB models is its inability to tolerate negative correlations. To address the shortcomings of MVP and MVNB, Chib and Winkelmann (2001) developed the MVPLN regression approach. The MVPLN model does not only account for overdispersion, but also it has a fairly general correlation structure allowing for (i) different covariance terms; and (ii) the possibility of negative correlations, which cannot be ruled out unquestionably for all classifications. Recently, applications of the MVPLN to model collision count data at different levels of severity were undertaken. Tunaru (2002) used an MVPLN non-regression model implemented 67 via WinBUGS for analyzing multiple response count data (classified by severity and the number of vehicles involved) taking into account complex correlation structures. Abdel-Aty et al (2006) analyzed multiple binary categorizations of collisions to identify the factors associated with their frequencies using seemingly unrelated negative binomial (SUNB) regression. These categorizations included different classifications including propertydamage-only and injury collisions. SUNB estimation proved to increase the efficiency of the crash frequency models by accounting for the disturbance correlation, reducing the standard errors, and providing better model fit. Park and Lord (2007) used MATLAB codes, tailored according to the MCMC algorithm of Chib and Winkelmann (2001), for jointly modeling collision frequency by severity using data from 451 three-leg un-signalized intersections in California obtained through the Highway Safety Information System (HSIS). Their results showed promise toward obtaining more accurate estimates by accounting for correlations in the multivariate collision counts and over-dispersion. Ma et al. (2008) coded a Gibbs sampler and two Metroplois–Hastings algorithms (Gilks et al., 1996) in the R language for the prediction of collision counts by severity using data collected from Washington State through the HSIS on 7773 rural two-lane roadways in the Puget Sound region. Their results indicated that the MVPLN approach offered better predictions than those from univariate Poisson and negative binomial models. The results indicated that there were statistically significant correlations between (as well as over-dispersion in) collision counts at different levels of injury severity. Several recommendations for highway safety treatments and design policies were suggested, e.g., they concluded that wide lanes and shoulders are key elements for reducing collision frequencies, as are longer vertical curves. Ye et al. (2009) developed a simultaneous equations model of collision frequencies by collision type for a sample of 165 rural intersections in 38 counties of Georgia. The model was formulated using multivariate Poisson regression with multivariate normal heterogeneity. A simulation-based maximum likelihood approach was used to estimate the parameters. The 68 model estimation results support the notion of the presence of significant common unobserved factors across collision types, although the impact of these factors on parameter estimates was found to be rather modest. Aguero-Valverde and Jovanis (2009) used Full Bayes MVPLN models implemented via WinBUGS to estimate the expected collision frequency for different collision severity levels using data on 6353 rural two-lane segments in central Pennsylvania. They found high correlations among collision severities, with highest values between contiguous severity levels. In terms of goodness-of-fit, the multivariate model outperformed the independent univariate PLN models and provided improved precision for collision frequency estimates. The multivariate estimates were combined with cost data from PennDOT to rank sites for safety improvements via expected collision cost and excess expected cost per segment. The results showed that the multivariate-based top ranked segments have consistently higher costs and excess costs than the univariate estimates as a result of higher multivariate estimates of fatalities and major injuries. These higher estimates resulted in different rankings for the multivariate and independent models. To corroborate the findings reported in the literature regarding the nature of the correlations across collision severity levels and their influence on safety analyses, a bivariate PLN model is compared with the independent (separate) univariate PLN models. Moreover, a new multivariate hotspots identification technique, which generalizes the univariate posterior probability of excess that has been commonly proposed and applied in the literature, will be presented. The models were estimated in a Full Bayes (FB) context via MCMC simulation and were compared in terms of their goodness-of-fit, inference, identification of hotspots and precision of expected collision frequency, using a data set composed of signalized intersections in the city of Edmonton, Alberta. 69 5.2 The Multivariate Poisson-Lognormal Model For a set of data on road collisions at n locations, where the collisions at each location are classified into K categories, define the vector yi = ( y i1 yi 2 ... yiK ) , where yik denote the ' number of collisions at the ith location in category k . Recall equations (19)-(21). In Equation (20) let X and β k denote the matrix of covariates and the vector of regression coefficients, respectively, and let β denote the set {β 1 , β 2 ,..., β K} . Thus, given ( X , β , Σ ) , the θ i are independently distributed as f (θ i | X , β , Σ ) = { }, exp − 0.5(θ *i − µ *i ) ′ Σ −1 (θ *i − µ *i ) ( 2π ) K /2 (∏ kK=1θ ik ) Σ 1/ 2 (41) which is a K-dimensional log-normal distribution, where ' θ i = (θ i1 θ i 2 ... θ iK ) , ' θ *i = (ln(θ i1) ln(θ i 2) ... ln(θ iK ) ) , ' µ i = (µ i 1 µ i 2 ... µ iK ) , ' µ*i = (ln( µ i 1 ) ln( µ i 2 ) ... ln( µ iK )) . It should be noted that the univariate PLN models are obtained as the special case where Σ is a diagonal matrix. 5.3 Multivariate Identification The treatment in this section is not meant to account for the extensive literature on the selection of hazardous locations, but rather to illustrate the generalization to the multivariate setting of some of the techniques that were developed for univariate analysis. Several hazardous location identification techniques exist in the literature. These include the probability of worst site, posterior distribution of the ranks of the collisions costs or expected number of collisions, or through expected collision cost and excess expected cost (Tunaru, 2002; Brijs et al., 2007; Aguero-Valverde and Jovanis, 2009). However, assigning costs to different types of injury can be controversial for a variety of reasons (Brijs et al., 2007). These reasons include ethical arguments (e.g., can we assign a 70 cost to a human life?) and economic arguments (what are the quantities that must be measured to estimate the cost of a seriously injured person?). To avoid such difficulties, an alternative approach of selecting hot spots is adopted which is based on the posterior probability of excess (Higle and Witkowski, 1988; Sayed and Abdelwahab, 1997; Heydecker and Wu, 2001). The Bayesian posterior probability that a site has an excessive mean collision frequency provides an indication of sites at which the future collision record is expected to be worse than usual. In the univariate analysis, this probability is expressed as: ∞ 2 ∫µ 0 f (θ i | y , X , β , σ ) d θ i > 1 − δ , i = 1,2,..., n , (42) where σ 2 , µ 0 and δ denote the extra-Poisson variation, an upper limit of the “acceptable” mean number of collisions and a threshold value between 0 and 1, respectively. In practice, µ 0 is typically specified as either the median or mean of the prior distribution, whereas δ has been arbitrarily assumed to equal 0.05 (Higle and Witkowski, 1988; Sayed and Abdelwahab, 1997). Under the multivariate Poisson-lognormal model, the criterion (42) generalizes to ∞ ∞ ∞ ∫µ10 ∫µ 20 ...∫µ K 0 f (θ i | y , X , β , Σ ) d θ i1 d θ i 2 ... d θ iK > 1 − δ , i = 1,2 ,..., n , (43) where µ k 0 denotes an upper limit of the “acceptable” mean number of collisions in category k. It should be noted that the univariate Equation (42) is the special case K=1. The evaluation of the multiple integral in (43) is rather complicated. Thereby, the following approximation, due to Clayton and Kaldor (1987), is proposed to simplify the computation of (43), f (θ i | y, X , β , Σ ) ≈ exp{− 0.5(θ *i − mi )′ S i−1 (θ *i − mi )} (2π ) K / 2 (∏kK=1θ ik ) S i 1/ 2 , (44) which is a K-dimensional log-normal distribution with ( y i1 + 0.5 ) ln ( y i1 + 0.5 ) − 0.5 −1 * ( yi 2 + 0.5 ) ln ( y i 2 + 0.5 ) − 0.5 mi = S i Σ µ i + , ... ( y + 0.5 ) ln ( y + 0.5 ) − 0.5 iK iK (45) 71 and yi 1 + 0.5 −1 0 Si = ∑ + ... 0 −1 0 yi 2 + 0.5 ... . ... ... ... 0 ... yiK + 0.5 0 ... 0 (46) The approximation (44) can be used in the Bayesian probability of excess criterion (43) to compare the univariate and multivariate approaches of selecting hot spots. ˆ i denote the Bayesian estimate of µ i , where the For univariate PLN models (K=1) let µ * * ˆi , subscript k has been dropped to simplify notation. Further, let µ 0 denote the average of µ * * i = 1,2 ,..., n . Then, the ith intersection would be selected as a hazardous location whenever Φ ( µ *i 0) < δ , (47) where Φ denotes the univariate standard normal distribution function and µ *i 0 = ( µ *0 − mi ) / S i . The standardized multivariate normal distribution function can be used in a similar procedure to select hot spots under MVPLN. Thus, using a similar notation, the ith intersection would be selected as a hazardous location whenever * * * Φ K ( µ i10 , µ i 20 ,..., µ iK 0 , Rˆ ) < δ , (48) where Φ K denotes the multivariate standard normal distribution function, µ *ik 0 = ( µ *k 0 − mik ) / S ikk , k = 1, 2, …, K, ˆ denotes the correlation matrix corresponding to Σˆ , the Bayesian estimate of the and R covariance matrix Σ . 5.4 Prior Distributions Diffused normal distributions (with zero mean and large variance) were used for the regression parameters and a Wishart ( P , r ) prior was used for Σ −1 , where P and r ≥ K represent the prior 72 guess at the order of magnitude of the precision matrix Σ −1 and the degrees of freedom, respectively. Choosing r = K as the degrees of freedom corresponds to vague prior knowledge (Spiegelhalter et al., 1996; Tunaru, 2002; Spiegelhalter et al., 2005). 5.5 Model Comparisons To show that DIC is additive under independent models and priors, let f ( y | θ ) and f ( y ) denote the conditional and marginal distributions of y, where θ denote the vector of parameters associated with y. Then, DIC = D + p , where p = D − D( θ ) , θ = E [ θ | y ] and D = E [ D( θ ) | y ] are the posterior means of θ and the Bayesian deviance D( θ ) = −2 ln{ f ( y | θ )} + 2 ln{ f ( y )} . (49) For K collision categories, let y and θ be partitioned as ( y1 ,..., y K ) and ( θ 1 ,...,θ K ) . Define DIC k = D k + p k , pk = D k − Dk ( θ k ) , D k = E [ Dk ( θ k ) | y k ] , θ k = E [ θ k | yk ] and Dk ( θ k ) = −2 ln{ f ( y k | θ k )} + 2 ln{ f ( y k )} . Under independent models and priors, we have that f ( y | θ ) = ∏ kK=1 f ( y k | θ k ) and f ( y ) = ∏kK=1 f ( y k ) . These multiplicative conditional and marginal distributions of y translate additively in the Bayesian deviance (49) leading to DIC = ∑ kK=1 DIC k . (50) A mathematical justification is provided in Appendix A. 5.6 The Data A total of n=99 signalized intersections in the city of Edmonton were investigated for the purpose of developing SPFs 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 collision frequencies and traffic volumes within the intersections were obtained from the city of Edmonton. The traffic and collision data were provided for a 3-year 73 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 term “collision” as defined by Edmonton’s Transportation Department, refers to reportable onstreet collisions that do not occur on private property, include at least one motor vehicle, and result in injury, at least $1,000 in property damage, or both. The database does not include nonvehicular collisions (e.g., cyclist hitting a pedestrian) and includes collisions that are forwarded to the department by the police service within a specified period. The term “intersection collisions” include collisions that occur inside and 10 m past the legally defined limits of the outer crosswalk lines of intersecting roads and include right-turn cut-offs. The collisions 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) are used to build SPFs for predicting aggregate number of collisions ln( µ ik ) = β k 0 + β k 1 X i 1 + β k 2 X i 2 , (51) where X i 1 = ln( V i 1 ) , V i 1 = AADT at the major approach, and X i 2 = ln( V i 2 ) , V i 2 = AADT at the minor approach. As mentioned in chapter 3 the aggregation is justifiable. A statistical summary of the data is shown in Table 5.1. Table 5.1 Statistical summary of data set (99 intersections) Minimum Maximum Mean Standard Deviation AADT on Major approach 5,642 61,744 28,010 9,705 AADT on minor approach 2,809 29,128 16,043 6,564 Number of PDO collisions 0 126 43.13 24.96 Number of I+F collisions 0 78 26.99 15.59 74 To explore the potential relationships between PDO/I+F and major & minor traffic volumes, Figure 5.1 displays the 3-dimensional relationships via wireframe plots. Figure 5.1 Wireframe plots (a) Wireframe plot of ln(PDO) vs. ln(V1) and ln(V2). (b) Wireframe plot of ln(I+F) vs. ln(V1) and ln(V2). (a) (b) Figure 5.1 lends support to the AADT functional form (51), as the relationship between each of PDO & I+F and major & minor traffic volumes, is effectively linear on the logarithmic scale, except for a very few cases where there were no collisions. The AADT functional form (51) 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). Further, it will be adopted for estimating safety performance in the forthcoming highway safety manual (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 collision frequency. Nevertheless, such models are well recognized in the road safety literature. 5.7 Results and Discussion 5.7.1 Model Results The posterior summaries in Tables 5.2 and 5.3 were obtained via two chains with 20,000 iterations 10,000 of which were excluded as a burn-in sample using WinBUGS. A Wishart prior with an identity scale matrix and two degrees of freedom was adopted (Chib and Winklemann, 2001; Congdon, 2006). Examination of the BGR statistics, ratios of the Monte Carlo errors 75 relative to the standard deviations of the estimates and trace plots for all model parameters indicated convergence. Table 5.2 Univariate models’ statistics Univariate Poisson-Lognormal PDO Intercept Ln(V1) (major AADT) Ln(V2) (minor AADT) σ 11 DIC I+F Intercept Ln(V1) (major AADT) Ln(V2) (minor AADT) σ 22 DIC Table 5.3 Estimate Standard Deviation 95% Credible Intervals Lower Limit Upper Limit -10.590 0.927 0.496 0.361 730.3 1.282 0.131 0.097 0.036 -13.090 0.667 0.307 0.297 -7.950 1.188 0.686 0.437 -9.267 0.754 0.493 0.415 684.3 1.447 0.148 0.105 0.043 -12.180 0.480 0.292 0.339 -6.440 1.042 0.706 0.507 Multivariate models’ statistics Multivariate Poisson-Lognormal 95% Credible Intervals Lower Limit Upper Limit Estimate Standard Deviation -10.680 0.902 0.531 0.163 1.347 0.138 0.104 0.030 -13.420 0.643 0.320 0.113 -8.250 1.179 0.723 0.230 -9.346 0.742 0.513 0.217 1.565 0.159 0.121 0.042 -12.410 0.435 0.276 0.145 -6.240 1.059 0.761 0.312 0.143 0.029 0.092 0.208 ρ =σ 12 / σ 11 σ 22 0.758 0.056 0.634 0.852 DIC 1373 PDO Intercept Ln(V1) (major AADT) Ln(V2) (minor AADT) σ 11 I+F Intercept Ln(V1) (major AADT) Ln(V2) (minor AADT) σ 22 Covariance σ 12 Correlation 76 The results of Tables 5.2 and 5.3 show that: (i) the parameter estimates are significant, as the 95% credible intervals are bounded away from zero; (ii) the regression coefficients are positive, indicating an increase in the mean collision 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 collisions than on I+F collisions. Table 5.2 and 5.3 reveal small differences in the estimates of the regression parameters between the univariate and multivariate models. However, the estimates of the extra Poisson variation parameters ( σ 11 and σ 22 ) were considerably smaller under MVPLN. In fact, the MVPLN 95% credible intervals were entirely shifted to the left of those obtained under the univariate PLN models, for both PDO and I+F. For i = 1,2 ,..., n and k = 1,2 ,..., K , the mean and variance of Y ik are given by E (Y ik ) = µ ik exp(0.5σ kk ) , (52) and Var (Y ik ) = E (Y ik ) + [ E (Y ik )]2 (exp(σ kk ) − 1) . (53) Since the second term in Equation (53) dominates the first, we have that Var (Y ik ) ≅ [ E (Y ik )]2 (exp(σ kk ) − 1) , leading to the ratio Var (Y ik | MVPLN ) / Var (Y ik | PLN ) ≅ (exp(σˆ kk ) − 1) /(exp(σ~ kk ) − 1) , (54) where σˆ kk and σ~ kk are the posterior estimates of σ kk under MVPLN and PLN, respectively. As precision is inversely proportional to the variance of expected collision frequency, it is estimated that MVPLN is more than twice as precise as the independent univariate PLN models; the actual precision ratios were computed via Equation (54) as 2.5 and 2.1 for PDO and I+F, respectively. The precision of the expected collision frequencies can be estimated also using 77 the standard deviations via Equation (53). Figure 5.2 displays these standard deviations by severity and model. Figure 5.2 reveals that (i) the standard deviations of the multivariate model are uniformly smaller than those of the PLN models, (ii) the relationship between the two sets of standard deviations appears to be linear for the two severity levels, and (iii) the standard deviations of the multivariate model were decreased 41% and 36% for PDO and I+F, respectively. With a larger data set for rural two-lane segments in central Pennsylvania, Aguero-Valverde and Jovanis (2009) obtained reductions in the standard deviations of 20.5% (PDO), 10.5% (minor injuries), 21.5% (moderate injuries), 48.1% (major injuries) and 40.7% (fatalities). Figure 5.2 Standard deviations (sd) of expected collision frequencies by severity and model 76 72 68 64 60 56 y = sd(MVPLN) 52 48 44 y = 0.59x 40 36 32 y = 0.64x 28 PDO 24 I+F 20 Linear 16 12 8 4 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 x = sd(PLN) The improvement in precision is due mainly to the fact that MVPLN accounts for the correlation between the latent variables θ i1 (PDO) and θ i 2 (I+F). This correlation ( ρ ) was estimated at 0.758, which is highly significant. In essence, higher PDO collisions are associated 78 with higher I+F collisions, as the collision likelihood for both types is likely to rise due to the same deficiencies in roadway design and/or other unobserved factors. The DIC statistics were 730.3, 684.3 and 1373 under the PDO PLN, I+F PLN and MVPLN, respectively. Thus, the MVPLN model provides a superior fit over the two univariate models as its (multivariate) DIC is much less than the sum of their (univariate) DICs; a very significant drop off of 41.6. 5.7.2 Omitted Variable Bias As noted previously, AADT-only models may suffer from omitted variables bias since the unobserved heterogeneity from other factors known to influence collision frequency (such as number of lanes, signal-control timing, speed limits, etc.) ends up in the correlation structure, thus affecting the estimated correlation. To study the magnitude of this effect on the correlation, the literature was consulted where several independent variables were used to develop multivariate SPFs. Park & Lord (2007) reported a number of correlations between PDO and four severity levels. Upon combining these 4 severity levels in one category (I+F) and using the covariance matrix, the correlation between PDO and I+F was estimated at 0.8608. A similar calculation using the results of Aguero-Velverde & Jovanis (2009) produced a correlation estimate of 0.58. Ma et al. (2008) didn’t report the variances associated with the various severity levels. However, assuming comparable precisions, their results produced a rough correlation estimate of 0.44. These results indicate that the current estimate of 0.758 is in the upper range of the estimates reported in the literature. Although the inclusion of additional relevant covariates is expected to reduce the correlation, it is conjectured to remain significant with crucial inferential impact. 79 5.7.3 Multivariate Identification of Hazardous Locations Hazardous locations were identified via Equations (47) and (48) using δ = 0.10, 0.05 and 0.01. For K=2, the standardized bivariate normal distribution function PROBBNRM (SAS, 2002) was used in the implementation of Equation (48). The PROBBNRM function required three input arguments ( µ 10 − mi 1 ) / * * S i 11 , ( µ 20 − mi 2 ) / S i 22 , and ρˆ = 0.758 . The outcomes appear in Table 5.4, where the results of the univariate Total PLN model were obtained through direct simulation of the sum PDO+I+F. Table 5.4 Number of hot spots selected under different PLN models Univariate PLN δ = 0.10 PDO I+F Total Hazardous? No Yes No Yes No Yes No 57 1 58 0 58 0 Yes 11 30 12 29 6 35 Multivariate PLN δ = 0.05 No 61 1 62 0 62 0 Yes 8 29 9 28 5 32 δ = 0.01 No 69 0 69 0 69 0 Yes 8 22 10 20 2 28 According to Table 5.4, the numbers of intersections identified as hot spots were 31, 29, 35 and 41 under PDO PLN, I+F PLN, Total PLN and MVPLN, respectively, for δ = 0.10 . The corresponding numbers were 30, 28, 32 and 37, for δ = 0.05 , and 22, 20, 28 and 30 for δ = 0.01 . As expected the numbers of hot spots was reduced as the threshold value decreased. However, several locations identified as hazardous by the MVPLN model were missed by the univariate models regardless of the choice of threshold value. It should also be noted that none of the intersections that were not identified as hazardous by MVPLN were identified by any of the univariate models except for one location under the PDO PLN model for δ = 0.10 and 0.05. The results of Table 5.4 demonstrate that some hazardous locations can be overlooked if the analysis was restricted to the univariate models. These results are consistent with those of 80 the cost ranking approach of Aguero-Valverde and Jovanis (2009) in that the MVPLN and independent PLN models yield different identification / ranking of hot spots. The advantage of MVPLN over separate univariate analyses is illustrated in Figure 5.3, where posterior information is also contrasted to prior knowledge. The figure displays two credible ellipses and a credible rectangle for (ln(θ PDO), ln(θ I + F ) ) , with 95% confidence. Figure 5.3 was prepared using the data for Site 56, which was identified as hazardous by MVPLN and I+F PLN but neither by PDO PLN nor by Total PLN. Similar patterns were observed for other identified sites. Figure 5.3 95% credible ellipses and rectangle for (ln( θ PDO ), ln( θ I + F )) , using the data for Site 56 The larger ellipse corresponds to the prior distribution in Equation (41). The area inside that ellipse corresponds to 95% probable values for (ln( θ PDO ),ln( θ I + F )) before the sample data 81 were collected. The smaller ellipse corresponds to the posterior distribution in Equation (44). The area inside that ellipse corresponds to 95% probable values for (ln( θ PDO ),ln( θ I + F )) after incorporating the sample data into the analysis. It is clear that the sample information has enhanced prior knowledge and resulted in a more focused “posterior” ellipse that is a (small) subset of the “prior” ellipse. The rectangle was constructed via the Bonferroni method using the two credible intervals obtained from the separate univariate models, PDO PLN and I+F PLN. The shaded area inside that rectangle corresponds to 95% probable values for (ln( θ PDO ),ln( θ I + F )) obtained by combining the two separate univariate posterior intervals. The posterior ellipse is nearly contained in the rectangle indicating that some values that are improbable under the MVPLN posterior analysis would be probable under the separate univariate analyses and vice versa. Also, some of the probable values under the rectangle seem to contradict prior knowledge as they lie outside the prior ellipse. 5.8 Summary In this chapter, a MVPLN regression was used to jointly analyze a sample of collision counts classified by two severity levels (PDO and I+F) at 99 urban intersections. To illustrate the importance of the multivariate technique, it was compared with the independent univariate PLN models with respect to model inference, goodness-of-fit, identification of hazardous locations and precision of expected collision frequency. The estimates of the extra Poisson variation parameters were considerably smaller under MVPLN. As precision is inversely proportional to the variance of expected collision frequency, it is estimated that the MVPLN model is more than twice as precise as the univariate PLN models. The improvement in precision is due mainly to the fact that MVPLN accounts for the correlation between the latent variables representing (PDO) and (I+F). This correlation ( ρ ) was estimated at 0.758, which is highly significant. This correlation was identified by other 82 researchers in the literature. This indicates that higher PDO collisions are associated with higher I+F collisions, as the collision likelihood for both types is likely to rise due to similar deficiencies in roadway design and/or other unobserved factors. In terms of the goodness-of-fit, the MVPLN provided a superior fit over the two univariate models as its (multivariate) DIC is much less than the sum of their (univariate) DICs. The differences between the DICs show a significant drop off of 41.6. In terms of model application, the research component introduces a new multivariate technique for the identification of hazardous locations. The new technique generalizes the univariate posterior probability of excess that has been commonly proposed and applied in the literature to fit the multivariate relationship between the latent variables. The results showed that all hazardous locations identified by the univariate models were identified by the multivariate model (except one location under PDO PLN). The results also demonstrated that some hazardous locations could be overlooked if the analysis was restricted to the univariate models. The multivariate nature of collision data is both logical and intuitive. Recent improvements in statistical tools allow for a more precise analysis which should be used to assist transportation engineers to better understand the relationships between the different modeling variables and techniques. 83 6 OUTLIER RESISTANCE MODELS This chapter compares three approaches for analyzing collision data with potential outliers when developing SPFs. The proposed techniques are based on the MVPLN regression approach that has been discussed in chapter 5. Robust multivariate techniques based on mixture models, which are capable of accommodating outliers, are proposed. 6.1 Background Collision 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., see Kulmala (1995) who used the leverage statistic and Sawalha and Sayed (2006) who used Cook’s distance for outlier detection. Typically, there are three approaches for analyzing data with potential outliers: (i) the “outlier ignoring” approach, (ii) the “outlier rejecting” approach, and (iii) the “outlier accommodating” approach. Proponents of the outlier 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). 84 On the other hand, the main problem with the popular outlier 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 chapter, 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 tdistribution 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. The authors 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. Two outlier resistance mixture models are proposed within a multivariate context and compared to the Multivariate Poisson Lognormal (MVPLN) approach presented in chapter 5. It is proposed that these mixture models are considered while developing SPFs as they provide robust inference. The models were estimated in a Full Bayes (FB) context via MCMC simulation and were compared in terms of their goodness-of-fit, inference, and precision of expected collision frequency, using a data set composed of signalized intersections in the city of Edmonton, Alberta. 85 6.2 The Methodology The three approaches discussed above “outlier ignoring”, “outlier rejection”, and “outlier accommodating” will be compared. The Multivariate Poisson-Lognormal (MVPLN), presented in chapter 5, will be used to construct models incorporating all data points (outlier ignoring approach). For the outlier 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 outlier 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. 6.2.1 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 (21) by a K-variate Student t distribution with zero mean, scale matrix Σ and ν degrees of freedom ε i ~ t K ( 0 , Σ ,ν ) . (55) However, the Student t-distribution in (55) can be derived by mixing the multivariate normal in (21) with a scaling variable ω i (Lange et al., 1989) so that ε i ~ N K ( 0 ,Σ / ω i ) , ω i ~ Gamma(ν / 2 ,ν / 2 ) . (56) For a univariate PLN model, where K = 1, Σ is replaced by σ 2 in (55) and (56). 6.2.2 The Two-Component Mixture A robust two-component mixture model is proposed, where (21) is replaced by ε i ~ π N K ( 0 ,Σ / δ ) + ( 1 − π ) N K ( 0 ,Σ ) . (57) 86 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 , 0.10 ) and δ ∈( 0.1, 0.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 σ 2 in (57). Alternatively, for a univariate NB model, (57) is replaced by the gamma mixture θ i ~ π Gamma ( δκ , δκ / µ i ) + ( 1 − π ) Gamma ( κ , κ / µ i ) , (58) where κ represents the inverse dispersion parameter of the NB model. 6.2.3 Detecting Outliers For univariate normal regression models, 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 outlier detection in NB regression. 87 For multivariate normal models, diagnostics measures were developed based on the Mahalanobis distance (Lange et al., 1989). Under (21) the Mahalanobis distance is given by * * 2 * * −1 d i = ( θ i − µ i )′ Σ ( θ i − µ i ) , 2 2 di ~ χ K (59) where θ *i = (ln( θ i1 ) ln( θ i 2 ) ... ln( θ iK )) and µ i = (ln( µ i 1 ) ln( µ i 2 ) ... ln( µ iK )) . Hence, ' ' * it is proposed that values of the Mahalanobis distance ( d i2 ) exceeding the critical value of the χ 2 -distribution with K degrees of freedom and an appropriately selected significance level α be identified as potential outliers. 6.2.4 Down-Weighting Outliers The Mahalanobis distance is closely related to the posterior mean of the scaling variable ω i in (56). In fact, going along the lines of Little (1988), it can be verified that E( ω i | y i ,Θ ) = (ν + K ) /(ν + d i2 ) , (60) where Θ denotes the set of all parameters under consideration. Thus, under the scale-mixture model, potential outliers are associated with small weights, particularly those ω i 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 (57) represent a discrete special case of those in (56), where ω i = δ with probability π , and ω i = 1 with probability 1 − π . Following Little (1988), it can be verified that 1 − π + π δ 1+ K / 2 e( 1−δ )d i / 2 E( ω i | yi ,Θ ) = . 2 1 − π + π δ K / 2 e( 1−δ )d i / 2 2 (61) However, as noted by Little (1988), the distributions of the posterior weights in (60) and (61) 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. 88 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. 6.3 Prior Distributions 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 . The following priors were used β kj ~ N ( 0 , 1002 ) , Σ −1 ~ Wishart ( P , 2 ) , where P = I the 2 x 2 identity matrix, ν ~ Exponentia l( γ ) and γ ~ Uniform ( 0.02 , 0.50 ) . 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 ) . The results showed that the models with δ = 0.1 provided better fits than those with δ = 0.2 (DIC = 1355 vs. 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. 6.4 The Data The same data in chapter 5 is used to compare the three approaches to analyze data with potential outliers. To recap a sample of n = 99 signalized intersections in the city of Edmonton 89 were investigated for the purpose of developing SPFs relating the safety of urban 4-leg intersections to their traffic flows. The collisions 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 a 3 years period) were used to build SPFs for predicting aggregate number of collisions using the link function (51). A detailed description of the data is given in Section 5.6. Note that Equation (51) is based on a simple AADT model that is 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. 6.5 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. 6.5.1 Outlier 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 90 the sites. Hence, the proposed mixture models were used to identify and down-weight potential outliers. The posterior estimates of the Mahalanobis distances in (59) are displayed in Figure 6.1 (left) for MVPLN along with a dashed line indicating the critical value of χ ( K = 2 , α = 0.01 ) = 2 9.210. Three outliers are easily identified, namely, those corresponding to sites 54, 4 and 5. These are the three highest points with d i2 exceeding the χ -critical value. It should be noted 2 that this list of outliers is different from the preliminary list obtained by the quick visual inspection. 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 Figure 6.1 (right) for the scale-mixture model. Under (56) the 2 2 Mahalanobis distance is given by d~ i = ω i d i2 and d~ i / K ~ F ( K ,ν ) . Using the estimate νˆ = 4.564 from Table 6.1, the dashed line in Figure 6.2 (right) indicates twice the critical value of the F-distribution, namely, c = 2 * F ( K = 2 ,νˆ = 4.564 , α = 0.15 ) = 6.328. It is seen that all 2 Mahalanobis distances d~i fall below c, even at such a high significance level as α = 0.15 . The scaling variable ω i was effective in down-weighting the outlying observations. Figure 6.1 Scatter plot of Mahalanobis distance by site. 2 (a) MVPLN model (dashed line corresponds to χ ( 2 , 0.01 ) = 9.210 ). (b) Scale-mixture model (dashed line corresponds to 2 * F ( 2 , 4.564 , 0.15 ) = 6.328 ). 16 6.5 4, 15.280 15 6.0 54, 14.360 14 5.5 13 5.0 12 4.5 10 Mahalanobis Distance (D) Mahalanobis Distance (D) 11 5, 9.911 9 8 7 6 4.0 3.5 3.0 2.5 5 2.0 4 1.5 3 1.0 2 0.5 1 0.0 0 0 5 10 15 20 25 30 35 40 45 50 Site (a) 55 60 65 70 75 80 85 90 95 100 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 Site (b) 91 Figure 6.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 ω i 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). Figure 6.2 Scatter plots of scaling weights by site. (a) Scatter plot of scaling weights by site for the scale-mixture model. (b) Scatter plot of scaling weights by site for the two-component mixture model. 1.4 1.0 1.3 0.9 1.2 0.8 1.1 1.0 0.7 0.9 0.6 Weights ω 0.8 0.7 0.5 0.6 0.4 0.5 0.3 0.4 0.3 0.2 5, 0.279 5, 0.138 0.2 0.1 4, 0.136 54, 0.100 4, 0.100 0.1 54, 0.053 0.0 0.0 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 Site (a) 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 Site (b) Since the χ - and F-critical values are readily available, the use of Mahalanobis 2 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 collisions 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. 92 6.5.2 Parameter Estimates Table 6.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 collision categories, the current application involves only two severity levels leading to bivariate distributions. The results in the second column of Table 6.1 (for MVPLN) are the same as those given in chapter 5. Such results correspond to the approach that ignores potential outliers. Table 6.1 PDO Intercept Ln(V1) Ln(V2) σ 11 I+F Intercept Ln(V1) Ln(V2)) σ 22 Parameter estimates and their associated standard deviations MVPLN 2-Component Mixture ScaleMixture Outlier Rejecting MVPLN -10.680 ± 1.347 0.902 ± 0.138 0.531 ± 0.104 0.163 ± 0.030 -11.320 ± 1.229 0.871 ± 0.118 0.633 ± 0.096 0.100 ± 0.019 -11.070 ± 1.201 0.908 ± 0.123 0.569 ± 0.101 0.097 ± 0.022 -11.560 ± 1.101 0.853 ± 0.113 0.676 ± 0.089 0.090 ± 0.017 -9.346 ± 1.565 0.742 ± 0.159 0.513 ± 0.121 0.217 ± 0.042 -9.793 ± 1.446 0.721 ± 0.145 0.585 ± 0.121 0.142 ± 0.028 ± 1.491 ± 0.153 ± 0.125 ± 0.032 -9.978 ± 1.315 0.685 ± 0.133 0.643 ± 0.111 0.134 ± 0.026 0.143 ± 0.029 0.078 ± 0.018 0.074 ± 0.021 0.065 ± 0.016 0.758 ± 0.056 0.651 ± 0.076 0.592 ± 0.086 1373 1355 0.644 ± 0.083 4.564 ± 1.737 1355 -9.582 0.750 0.535 0.135 Covariance σ 12 Correlation ρ ν DIC 1322 The results of Table 6.1 show that: (i) the parameter estimates are significant; (ii) the regression coefficients are positive, indicating an increase in the mean collision 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 collisions than on I+F collisions. 93 6.5.3 Results of the Scale and Two-Component Mixture Models Table 6.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 ( σ 11 and σ 22 ) were considerably smaller under the mixture models compared with MVPLN. Note that: Var ( Y ik | Mixture Model ) / Var ( Y ik | MVPLN ) ≅ (exp( σˆ kk / ωˆ i ) − 1 ) /(exp( σ~ kk ) − 1 ) , (62) 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-component 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 for 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. 6.5.4 Results of the Outlier 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 6.1. Although the estimates of the regression parameters under the outlier rejecting MVPLN are not significantly different from those of the 94 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 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. 6.6 Summary In this chapter, a scale-mixture model and a two-component mixture (contaminated normal) model were proposed as robust alternatives to MVPLN regression for the joint analysis of a sample of collision counts classified by two severity levels (PDO and I+F) at 99 urban intersections. 95 For such a multivariate 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. 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 outlier ignoring approach, and (ii) avoids the risk of underestimating this uncertainty, which results from using the outlier 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 collisions are associated with higher I+F collisions, as the collision likelihood for both types is likely to rise due to similar deficiencies in roadway design and/or other unobserved factors. The degree of “success” of the outlier accommodating approach will (to some extent) depend 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 datasets is required to confirm the research findings. Future topics of research might also include adapting the diagnostic measures that have been developed for univariate normal regression models to GLM settings as well as developing similar measures for multivariate models. 96 PART II: SAFETY PERFORMANCE FUNCTIONS FOR INTERVENTION ANALYSIS 7 TWO-WAY INTERACTION MODELS This chapter advocates the use of a two-way interaction model to conduct a full Bayes (FB) Controlled Before-After (CBA) evaluation when subject to data constraints or limitations. 7.1 Background In chapter 2 a variety of methods to analyze collision data and to determine the effectiveness of the applied countermeasures were discussed. It was concluded that Bayesian methods are frequently preferred to conventional methods for their ability to treat unknown parameters such as predicted collision frequency as random variables having their own probability distributions. Examples of Bayesian evaluation techniques include the Empirical Bayes (EB) and fully Bayes (FB) CBA studies (Sayed et al., 2004; Pawlovich et al., 2006; Persaud and Lyon, 2007). A typical EB CBA study requires the collection of data for three distinct sets of data: i) treatment sites, ii) comparison sites, and iii) reference sites. The treatment group is composed of all sites with an implemented safety countermeasure. The comparison group is used to correct for time trend effects and includes sites that have not been treated but exhibit similar traffic and environmental conditions. The reference group is used to correct for the regressionto-the-mean (RTM) artifact. Usually, the reference group includes a larger number of sites that are similar to the treatment sites and is used to develop a 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. Alternatively, the FB approach has been proposed in the road safety literature to conduct CBA studies. The FB approach is appealing for several reasons, which can be categorized into methodological and data advantages. In terms of methodological advantages, the FB approach has the ability to account for all uncertainty in the data, to provide more detailed inference, and 97 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 thereby negating the need for a reference group. 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). However, in many circumstances data limitations might restrict the application of an EB or a FB intervention approaches to conduct CBA safety studies. These data limitations are usually attributed to i) the lack of sufficient data to develop the SPF that is necessary to conduct an EB analysis or ii) the availability of collision data in aggregated (either in space or time) forms that limits the use of a FB intervention model. In the light of these two data restrictions, a novel FB two-way interaction model is proposed. The approach is demonstrated using a data set from the city of Vancouver, British Columbia that is subject to certain data restrictions. In 1998, the Insurance Corporation of British Columbia (ICBC) initiated the Stop-Sign Infill (SSIF) plan. The SSIF program funds the conversion from uncontrolled residential intersections to two-way stop controlled intersections in an alternating pattern. It is required to evaluate the effectiveness of this program. The data available for the analysis consist of the aggregate number of claims (based on collisions) ( Y i ) in zone i during the observation period ( t i ), i = 1,2 ,..., n , where n denotes the total number of observations (sample size). As there were a treatment group and another matched yoked comparison group for each zone and for both groups there were data during the before-after periods, the sample size is n = 4 × N z , where N z = 22 denotes the number of zones. Due to the above data restrictions the FB intervention approach (to be discussed in chapter 8) cannot be used to determine the effectiveness of the SSIF plan. Therefore, to benefit from the additional advantages of the FB approach a two-way interaction model is proposed 98 based on the Poisson-Lognormal hierarchy with covariates representing the time, treatment (intervention) and interaction effects (reflecting the cross product between the time and treatment effects). Since the selected intersections were aggregated into 22 zones, the FB twoway interaction model was extended to account for variation among zones through the mean (using a random parameters model) and variance function (using a random effects model). All models were estimated in a Full Bayes context via MCMC simulation and were compared in terms of their goodness-of-fit and inference. 7.2 The Models Four models are developed where the data are available in a before-after design involving matched yoked comparison groups. The first model is a Poisson-Lognormal (PLN) regression model with covariates representing the time, treatment (intervention) and interaction effects (reflecting the cross product between the time and treatment effects). A variance component is added in the second model (PLN-V) to account for the variation among zones. In the third model (PLN-RP), random parameters are incorporated to account for the between-zone variation. The last model (Poisson-RP) is a Poisson regression model similar to (PLN-RP). Each of the four models has its own assumptions.. However, there are three main assumptions that are made throughout the analysis. First, since the number of sites varies from one zone to another, zone results are bound to be influenced by the inherent variability in zone size. Thereby, it is assumed that the number of collisions at each site follows a Poisson distribution. Consequently, the aggregate number of collisions would also be a Poisson random variable, whose variance, which is equal to its mean, is going to be proportional to zone size. Second, since the selection of the sites for treatment was not based on their collision history, it is assumed that the regression to the mean (RTM) effect will be minimal. Third, since the treated areas are all within residential neighborhoods, it is assumed that there was no significant 99 increase in traffic volume. Further, knowing that the population is growing in the Greater Vancouver area of British Colombia, it is assumed that the traffic volumes will increase with increased population, which would result in under-estimation rather than an over-estimation of treatment effectiveness. 7.2.1 The Poisson Lognormal Model (PLN) The PLN model is used to address over-dispersion for unobserved or unmeasured heterogeneity. According to this model, it is assumed that Y i | θ i ~ Poisson ( t i × θ i ) , (63) where θ i , the collision rate (average collisions per year) at zone i, is modeled as θ i = µ i exp( ui ) , (64) ln( µ i ) = β 0 + β 1 X i 1 + β 2 X i 2 + β 3 X i 3 , (65) 2 u i ~ N ( 0 ,σ u ) , (66) X i 1 = 1 , in the after period, X i 1 = 0 , in the before period, X i 2 = 1 , for treated sites, X i 2 = 0 , for comparison sites, X i 3 = X i 1 × X i 2 , and σ u2 denotes the extra Poisson variance. Thus, X i 1 is a time indicator, X i 2 is a treatment indicator, and X i 3 is a two-way interaction term allowing different treatment effects across time. Under PLN, it is easily shown that the odds ratio is given by OR = exp( β 3 ) . (67) The treatment effect is calculated as (OR − 1), while the percentage of reduction in predicted collision frequency is given by (1 − OR) × 100. 7.2.2 The Poisson Lognormal Model with a Random Zone Effects (PLN-VC) To account for the between-zone variation, suppose that the ith observation belongs to zone z( i ) ∈ { 1,2 ,..., N z } . Under this model it is assumed that θ i = µ i exp( ui ) exp( δ z( i ) ) , δ z ( i ) ~ N ( 0 ,σ δ2 ) , (68) 100 where σ δ2 denotes the variance component representing the variation among different zones. It should be noted that each zone involves observations from treated and comparison sites before-after the intervention. Under PLN-VC, the odds ratio is given by Equation (67). + + Let µ i = µ i exp( δ z ( i ) ) and note that ln( µ i ) = β z ( i ),0 + β 1 X i 1 + β 2 X i 2 + β 3 X i 3 , where β z ( i ),0 = β 0 + δ z ( i ) . The PLN-VC model is thus equivalent to a PLN with a random intercept that varies across comparison-treatment pairs. 7.2.3 The Poisson Lognormal Model with Random Parameters among Zones (PLN-RP) Alternatively, the zone variation can be represented by allowing the regression coefficients in Equation (65) to vary randomly from one zone to another ln( µ i ) = β z ( i ),0 + β z ( i ),1 X i1 + β z ( i ),2 X i 2 + β z ( i ),3 X i 3 , (69) and β z( i ), j ~ N ( β j ,σ 2j ) , j = 0 ,1,2 ,3 . (70) The PLN-RP model permits the computation of the zone-level odds ratios ORk = exp( β k 3 ) , k = 1,2 ,..., N z , (71) as well as the overall odds ratio from Equation (67). 7.2.4 The Poisson Model with Random Parameters among Zones (Poisson-RP) According to this model, Y i | µ i ~ Poisson( t i × µ i ) , where µ i is given by equations (69) and (70). Under Poisson-RP, the zone-level odds ratios are given by Equation (71), while the overall odds ratio is given by Equation (67). 7.3 The Data The Insurance Corporation of British Columbia (ICBC) Stop-Sign In-Fill (SSIF) program involved the installation of stop signs at approximately 1450 intersections in Vancouver and Burnaby between 1998 and 2003. The SSIF program is a system wide improvement program. It differs 101 from a traditional “black spot” improvement program where the selection of treatment sites is based on selecting sites with unusually high frequency of collisions. The collision data set used in this study is based on ICBC insurance claims that have been proven to be very useful in road safety analysis (Sayed et al., 2004; de Leur and Sayed, 2001). As such, the before-after claims (collisions) data were obtained and used in this analysis. Unfortunately, traffic volume data at the treated locations were not available. Municipalities typically only gather traffic volume data at residential neighborhoods by request of residents. Due to the time-intensive process of generating claims data, it was decided to evaluate only the intersections improved in the city of Vancouver. The sample was further reduced to a smaller number of intersections based on the following criteria: i) an accurate stop sign installation date must be available for each intersection; and ii) the completion date must be in the years 1999, 2000, 2001, or 2002. These dates provided adequate time for a proper beforeafter safety evaluation. The city of Vancouver is split into zones as an organized method to complete the stop sign in-fill throughout the city. The size and number of intersections in each zone were dependent on the major surrounding roads. A data set of 22 zones (Figure 7.1) containing a total of 380 intersections treated by the SSIF program were used in the evaluation process (Table 7.1). Within the 22 traffic zones, there were 133 untreated intersections. These 133 intersections were used to form various matching comparison groups to correct for time trend effects, including the confounding factors of history and maturation. Thus, for each zone, the group of treatment sites was matched with a yoked comparison group of intersections satisfying two conditions: i) they experience similar traffic and environmental conditions to those of the treated sites, but were not subject to any treatment; and ii) they are in close proximity to the treated sites. 102 Figure 7.1 Table 7.1 Zone # 1 2 3 4 5 6 7 8 9 10 11 7.4 Zones selected for evaluation in the city of Vancouver. Zones selected for evaluation Total # of SSIF Program Intersections Treated 29 12 28 24 29 21 13 13 18 12 13 8 25 25 18 6 12 7 27 20 25 13 Zone # 12 13 14 15 16 17 18 19 20 21 22 Total Total # of Intersections 46 23 32 41 5 25 23 17 18 20 26 513 SSIF Program Treated 40 23 26 35 4 17 19 11 12 15 17 380 Results and Discussion The posterior summaries in Tables 7.2 and 7.3 were obtained via two chains with 20000 iterations each, 10000 of which were excluded as a burn-in sample, using WinBUGS with a diffused N ( 0 , 100 2 ) prior for the regression parameters and an uninformative inverse- 103 Gamma ( 0.01, 0.01 ) prior for the dispersion parameters. 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. Table 7.2 summarizes the parameter estimates and their 95% credible intervals for the four models. Albeit the interaction effect under PLN was negative, it is not significant as the 95% credible interval for β 3 includes 0, implying a non-significant reduction in predicted collision frequency. However, accounting for the zone variation in PLN-VC produced a significant (negative) interaction effect. Thus, although PLN-VC provided a similar fit to PLN (the DIC values were 522.5 and 525.9, respectively), the incorporation of the random zone effects has resulted in a significant OR whose 95% credible interval ranged from 0.252 to 0.905 (< 1). However, the random parameters PLN-RP model has fitted the data much better than PLN and PLN-VC with significant reductions in DIC of 45.2 and 41.8, respectively. The results of PLN-RP indicate that the parameter estimates are significant as the 95% credible intervals were bounded away from zero and that random parameters were associated with the intercept as well as all three covariates. The estimates, standard deviations and the associated 95% credible intervals for σ u2 were 0.797 ± 0.170 (0.525, 1.184), 0.604 ± 0.149 (0.373, 0.947) and 0.015 ± 0.010 (0.004, 0.041) under PLN, PLN-VC and PLN-RP, respectively. Hence, the estimate of 0.797 under PLN has been reduced by 24.2% and 98.1% under PLN-VC and PLN-RP, respectively. For PLN and PLN-VC, the standard deviations are rather small relative to the parameter estimates and the lower limits of the credible intervals are noticeably greater than zero demonstrating the presence of over-dispersion in the data. In contrast, under PLN-RP the standard deviation of 0.010 is relatively large in comparison to the parameter estimate of 0.015 indicating the nonsignificance of σˆ u2 . This is confirmed by the small credible interval (0.004, 0.041), which will never include zero for a variance estimate since the conditional posterior has always a positive 104 domain. Thus, allowing the PLN to have random parameters (among zones) seems to reduce the need to account for over-dispersion as the estimate of the extra Poisson variance became fairly small. This led to the Poisson-RP model, which appears to provide comparable results to PLN-RP with respect to goodness-of-fit and parameter estimates. Hence, for this data set, the simpler random parameters (among zones) Poisson-RP model is preferred over PLN-RP. The above findings confirm the conclusions of Miaou and Song (2005) and Aul & Davis (2006). The former authors have suggested that careful modeling of the random effects may reduce the unexplained over-dispersion typically observed in collision data, whereas the latter authors have argued that over-dispersion may to a large extent be an artifact of mixing fundamentally heterogeneous entities and that it seems prudent to avoid dependence on methods that require the observation of over-dispersion before analysis can proceed. The overall odds ratios were 0.533, 0.495, 0.472 and 0.489 leading to reductions in predicted collision frequency of 44.7%, 50.5%, 52.8% and 51.1%, under PLN, PLN-VC, PLN-RP and Poisson-RP, respectively. To determine the significance of these reductions the 95% credible intervals in Table 7.2 are examined. For PLN, the 95% credible interval for the OR straddles 1 indicating a non significant reduction in predicted collision counts. Moreover, the 95% credible interval for PLN-VC is too wide to be useful. In contrast, PLN-RP and Poisson-RP produced rather compact 95% credible intervals for OR entailing significant reductions ranging from 38.8% to 64.6% for PLN-RP and from 36.8% to 62.3% for Poisson-RP. For the overall OR, the traditional technique produced an estimate of 0.552 with a standard deviation of 0.152 (Sayed et al., 2006). These estimates are in close agreement with those of PLN and PLN-VC. Although not as close, the estimate of the traditional overall OR is compatible with PLN-RP and Poisson-RP. The small difference between the traditional and FB results can be attributed to the way sites were selected which does not produce high RTM effect (the SSIF is a system wide improvement program where sites are not selected because of high collision frequency). 105 Table 7.2 β0 β1 β2 β3 σ u2 Parameter estimates, standard deviations and 95% credible intervals PLN 0.796 ± 0.175 (0.458, 1.172) -0.219 ± 0.275 (-0.766, 0.317) 0.732 ± 0.250 (0.232, 1.166) -0.655 ± 0.358 (-1.375, 0.004) 0.797 ± 0.170 (0.525, 1.184) σ δ2 PLN-VC 0.777 ± 0.203 (0.341, 1.207) -0.190 ± 0.233 (-0.651, 0.247) 0.833 ± 0.226 (0.421, 1.329) -0.754 ± 0.314 (-1.377, -0.100) 0.604 ± 0.149 (0.373, 0.947) 0.226 ± 0.162 (0.015, 0.624) σ0 σ1 σ2 σ3 OR DIC 0.553 ± 0.194 (0.253, 1.004) 525.9 0.495 ± 0.160 (0.252, 0.905) 522.5 PLN-RP 0.716 ± 0.267 (0.163, 1.224) -0.210 ± 0.100 (-0.402, -0.018) 0.848 ± 0.324 (0.212, 1.493) -0.762 ± 0.142 (-1.039, -0.491) 0.015 ± 0.010 (0.004, 0.041) Poisson-RP 0.728 ± 0.269 (0.177, 1.243) -0.241 ± 0.099 (-0.442, -0.057) 0.838 ± 0.332 (0.192, 1.511) -0.725 ± 0.133 (-0.976, -0.459) 1.158 ± 0.230 (0.790, 1.677) 0.132 ± 0.057 (0.056, 0.271) 1.417 ± 0.261 (0.996, 2.009) 0.172 ± 0.085 (0.067, 0.387) 0.472 ± 0.067 (0.354, 0.612) 480.7 1.169 ± 0.235 (0.794, 1.707) 0.141 ± 0.063 (0.058, 0.298) 1.444 ± 0.271 (1.016, 2.069) 0.172 ± 0.081 (0.065, 0.373) 0.489 ± 0.066 (0.377, 0.632) 478.3 In addition to the overall estimate of β = ( β 0 , β 1 , β 2 , β 3 ) , provided in Table 7.2, there were 22 estimates of β under PLN-RP & Poisson-RP corresponding to the 22 zones (pairs). Such estimates were used in Equation (71) to compute the zone-level odds ratios. The results under Poisson-RP are reported in Table 7.3. According to Poisson-RP, the zone-level odds ratios ranged between 0.443 and 0.556 and there were significant reductions in predicted collision frequency in all zones, at the 5% level, as the upper bounds of all 22 credible intervals are less than 1. In contrast, according to the traditional approach, the zone-level odds ratios ranged between 0.260 and 1.333 and there were 16 zones exhibiting reductions, of which only 4 were significant; two zones exhibiting increase albeit non-significant; and 4 zones where zero collision counts were observed and 106 thereby their odds ratios could not be computed (Sayed et al., 2006). Furthermore, of the 18 odds ratios that could be computed by the traditional technique, only 10 were covered by the 95% credible intervals produced by the FB approach. Table 7.3 Zone 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 Overall Zone-level odds ratios for Poisson-RP OR 0.527 0.478 0.477 0.463 0.480 0.503 0.466 0.488 0.495 0.498 0.501 0.505 0.443 0.506 0.538 0.500 0.493 0.499 0.491 0.519 0.556 0.482 0.489 s. d. 0.128 0.095 0.099 0.101 0.101 0.115 0.095 0.111 0.111 0.102 0.111 0.096 0.089 0.112 0.116 0.122 0.104 0.101 0.101 0.117 0.142 0.102 0.066 95% C. I. (0.339, 0.830) (0.316, 0.686) (0.305, 0.697) (0.279, 0.679) (0.310, 0.705) (0.323, 0.762) (0.294, 0.669) (0.303, 0.735) (0.313, 0.751) (0.329, 0.723) (0.325, 0.764) (0.348, 0.720) (0.275, 0.632) (0.329, 0.781) (0.362, 0.814) (0.313, 0.782) (0.320, 0.729) (0.334, 0.720) (0.324, 0.718) (0.340, 0.785) (0.367, 0.925) (0.309, 0.710) (0.377, 0.632) The above results support the conclusions of Carriquiry & Pawlovich (2005) who noted that since scarce data would result in an unreliable safety estimates for individual sites; it becomes important to ‘borrow’ information from similar sites to complement the insufficient data. For PLN-RP and Poisson-RP, this is achieved via the FB approach, which combines information from all available exchangeable (or similar) sites in order to better estimate site-level odds ratios. The combined estimate typically takes on the form of a weighted average, where the weights depend on the relative amount of information on the site of interest. 107 Table 7.3 is also useful for retrospective analyses. For instance, if one discovers later on that certain zones have similar characteristics, one could group similar zones and assess the impact of applying the SSIF program given zone characteristics using the results of Table 7.3 and Equation (13) in Sayed et al. (2006). 7.5 Summary In this chapter, four models were proposed to conduct a CBA evaluation, within a full Bayes context, of the collision frequencies in the zones (neighborhoods) improved by the installation of stop signs under the ICBC’s SSIF plan. The first Model is a PLN regression model with covariates representing the time, treatment (intervention) and interaction effects. A variance component is added in PLN-VC to account for the variation among zones. In PLN-RP, random parameters among zones are incorporated to deal with the between-zone variation. The last model Poisson-RP is a Poisson regression model with random parameters among zones. The negative interaction effects were statistically significant, implying significant reduction in predicted collision frequency under all models except PLN. Recognizing the variation among zones has significantly improved the fit while reducing the need to account for over-dispersion. Thus, for this data set, the simpler random parameters Poisson model PoissonRP is preferred over PLN-RP. For the overall OR, there were close agreement in terms of primary and secondary magnitudes (means and standard deviations) between the estimates of the traditional approach and those of PLN and PLN-VC. Although not as close, the estimate of traditional overall OR is compatible with PLN-RP and Poisson-RP. It seems that the random selection of sites, which reduces the RTM effect, is the reason that the traditional approach is giving relatively similar overall-level results. At the zone-level, Poisson-RP (and PLN-RP) produced different OR estimates from those of the traditional technique. Fewer significant estimates were obtained from the latter and 108 many estimates were not covered by the 95% credible intervals of the FB approach, which combines information from all available exchangeable (or similar) sites in order to better estimate site-level odds ratios. Consequently, the standard deviations of the FB odds ratios were much smaller under PLN-RP and Poisson-RP than those of the traditional technique, which was also true for the overall OR. The results of this research provide some strong evidence for the benefit of incorporating such design features as matched yoked comparison groups as well as clustering (zoning) in SPFs. Apart from the improvement in the goodness of fit, such extended models can be used to account for heterogeneity due to unobserved road geometrics, traffic characteristics, environmental factors and driver behavior. As well, the random parameters models allow the computation of odds ratios at both the zone-level and the overall-level. 109 8 LINEAR INTERVENTION MODELS This chapter illustrates the use of an intervention model to conduct a full Bayes CBA evaluation. Moreover, the odds ratio is decomposed under the linear intervention SPF into basic components related to direct and indirect treatment effects expressed as a function of model parameters. 8.1 Background When collision data is available in disaggregate (in time and space) form, several researchers have proposed the use of a Full Bayes (FB) intervention model to conduct a CBA safety evaluation (Pawlovich et al., 2006; Li et al., 2008). The approach acknowledges that the treatment (intervention) effects do not occur instantaneously but are spread over future time periods. Using dynamic regression models the approach identifies the lagged effects of the intervention to describe their response. In this chapter, four different intervention models based on the Poisson-Lognormal hierarchy are compared. The first and second models were developed to deal with the immediate (PLNI with jump) and gradual (PLNI without jump) treatment impacts, respectively. The third and fourth models extend the PLNI models with and without a jump by including random parameters to account for the correlation between sites within comparison-treatment pairs. In all forms, linear slopes were tacitly assumed to represent the time and treatment effects across the treated and comparison sites. Moreover, the odds ratio (OR) is decomposed under the linear intervention SPF into basic components related to direct and indirect treatment effects, where the latter 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 an OR component corresponding to the direct treatment effects cannot be over emphasized as it enables the analyst to assess the effectiveness of the countermeasures apart from local (site-related) environmental factors. Also, 110 the consequences of the linearity assumption on these components, and subsequently on the OR, are investigated. Another objective is to provide straightforward equations in terms of model parameters for the computation of the OR and its associated components without resorting to additional algorithms such as the one proposed by Park et al. (2009). The objectives of this chapter are threefold: i) to provide a better understanding of the role of model parameters and their effects on the overall odds ratio estimation, ii) to simplify the calculation of the odds ratio by providing straightforward equations in terms of model parameters for the computation of the OR and its associated components without resorting to additional algorithms, and iii) to provide tools for investigating the safety change 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 ICBC’s Road Improvement Program (RIP), where the full Bayes approach is utilized to determine the effectiveness of certain countermeasures that were implemented in 18 treated intersections in the Greater Vancouver area. 8.2. The Models Let Y it denote the collision count observed at site i (i = 1, 2, …, n) during year t (t = 1,2, …, m). It is assumed that the data are available for a reasonable period of time before and after the intervention. 8.2.1 The Poisson-Lognormal Intervention (PLNI) Models The gradual (without a jump) PLNI model is given by (1), (5), (24) and (9), with the subscript i replaced by it where appropriate. Equations (5) and (24) show that it is possible to fit a different model to each site by including a site-level random effect ε i among the covariates, so that the resulting linear regression model involves different intercepts across sites. 111 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 Equation (25) in order to accommodate such an effect leading to the jump PLNI model. 8.2.2 The PLNI Models with Random Parameters among Matched 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 comparisontreatment pairs. To account for this correlation, suppose that the ith site belongs to the pair p (i ) ∈ {1,2,..., N C} , where N C denotes the number of comparison groups. The variation due to the comparison-treatment pairing can be represented by allowing the regression coefficients in Equation (24) to vary randomly from one pair to another. This leads to the gradual (without a jump) PLNI-RP model with random parameters among matched pairs ln( µ it ) = α p( i ),0 + α p( i ),1 T i + α p( i ),2 t + α p( i ),3 ( t − t 0 i ) I it + α p( i ),4 T i t + α p( i ),5 T i ( t − t 0 i ) I it + γ p ( i ),1 ln( V 1it ) + γ p( i ),2 ln( V 2 it ) + γ p( i ),3 X 3i + ... + γ p ( i ),J X Ji , (72) and 2 α p( i ), j ~ N ( α j , σ α2 , j ) , j = 0 , 1, ..., 5 , γ p( i ), j ~ N ( γ j , σ β , j ) , j = 1, ..., J . (73) The jump PLNI-RP model with random parameters among pairs is given by ln( µ it ) = α p( i ),0 + α p( i ),1 T i + α p( i ),2 t + α p( i ),3 ( t − t 0 i ) I it + α p( i ),4 T i t + α p( i ),5 T i ( t − t 0 i ) I it + α p( i ),6 T i I it + γ p( i ),1 ln( V 1it ) + γ p( i ),2 ln( V 2 it ) + γ p( i ),3 X 3 i + ... + γ p( i ),J X Ji , (74) where α p ( i ),6 ~ N ( α 6 , σ α2 ,6 ) . The random parameters PLNI-RP models enable meaningful inference both at the pairlevel and the overall-level. 112 8.3 Odds Ratio Decomposition Recall Equation (26) and let LTAi , LTBi , LCAi and LCBi denote the averages of ln(µ it ) corresponding to µ TAi , µ TBi , µ CAi and µ CBi , respectively. Hence, by averaging the predicted collisions on the log-scale, the odds ratio can be alternatively obtained from ln( ORi ) = LTAi + LCBi − LTBi − LCAi . (75) To gain some insight into the effects of the regression parameters on the odds ratio, consider the jump PLNI model and note that the regression parameters ( α 0 ,α 1 ,α 2 ,α 3 ,γ 3 , ...,γ J ) cancel out of (75). Further, for the ith treated site and t > t 0i , let s (i ) = t − t 0i 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 ORis denote the odds ratio for the ith treated site after s post-intervention years. It can be verified (see Appendix B) that OR is = C 1 ( i , s ) C 2 ( i , s ) C 3 ( i , s ) C 4 ( i , s ) , (76) C 1 ( i , s ) = exp{ 0.5 α 4 ( s + t 0 i + 1 )} , (77a) C 2 ( i , s ) = exp{ α 6 + 0.5α 5 ( s + 1 )} , (77b) γ C 3 ( i , s ) = [ OR( V 1is )] 1 , (77c) γ C 4 ( i , s ) = [ OR( V 2 is )] 2 , (77d) where OR( V 1is ) and OR (V 2is ) are odds ratios defined in terms of the major and minor traffic volumes the same way ORis is defined in terms of the predicted collision counts. Equations (77a)-(77d) hold also for the gradual PLNI model and for the random parameters PLNI models at the overall-level of inference. Yet, the parameter α 6 , which appears in C 2 (i, s ) , is set to zero for the models without a jump. It should also be noted that the smaller 113 these components, the larger the reductions in predicted collisions and the more effective are the countermeasures (treatments). For the gradual models, the two components C1 (i, s ) and C 2 (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 α 4 would be negative and the odds ratio will be reduced. In contrast, the odds ratio would be increased for positive values of α 4 . On the other hand, if the decrease in collision counts accelerates after treatment implementation then α 5 would be negative and the odds ratio will be reduced and vice versa. For the jump models, the component C 2 (i, s ) involves the additional parameter α 6 representing a possible post-treatment sudden drop/increase in collisions. Typically, a sudden post-treatment drop in collisions is likely to reduce the odds ratio and vice versa. The two components C 3 (i, s ) and C 4 (i, s ) involve the odds ratios of the major and minor traffic volumes raised to the powers γ 1 and γ 2 , which typically satisfy the conditions 0 < γ 1 ,γ 2 < 1 . Thus, both components would be closer to 1 than the two traffic volume odds ratios 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 C 2 (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 ) = exp{a + bs} , a = α 6 + 0.5 α 5 + 0.5 α 4 (t 0i + 1) , b = 0.5(α 5 + α 4) . (78a) According to (78a), C T (i,1) = exp{a + b} and C T ( i , s ) = exp{ b }. C T ( 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, [∂ C T (i, s ) / ∂s ] / C T (i, s ) = b , is constant during the period following the intervention. 114 Similarly, the components C 3 (i, s ) and C 4 (i, s ) can be combined to obtain γ γ C V ( i , s ) = [ OR( V 1is )] 1 [ OR( V 2 is )] 2 . (78b) The combined component in (78b) is inversely related to the indirect (through traffic volumes) treatment impact under the PLNI models. 8.4 Case Studies 8.4.1 Iowa 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). The authors considered a variety of intervention models where road length and traffic volume were used as offset variables in the regression equations, which included three additional variables representing quarterly seasonal effects on collisions. 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 (72) and (74) where p(i) is replaced by i. Thus, the odds ratios are decomposed as follows OR is = C 1 ( i , s ) C 2 ( i , s ) C V ( i , s ) , s = 1,2, ... . 115 where C V (i, s ) = OR (V is ) and V is 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 (with a jump). The estimates (αˆ 4 = 0.001, αˆ 5 = −0.020) were obtained under H1, while the estimates (αˆ 4 = 0.000, αˆ 5 = −0.010, αˆ 6 = −0.355) were obtained under H2. Using these estimates along with t 0 i = 216 , the component C1 (i, s ) representing the treatment by time interaction was computed from exp{0.5 αˆ 4 ( s + 217 )} , whereas the component C 2 (i, s ) representing the treatment by intervention date interaction was computed from exp{αˆ 6 + 0.5 αˆ 5 ( s + 1)} . Figure 8.1 Monthly treatment components for a median site under H1 (w/o a jump), Li et al. (2008). 116 For H1, where the estimate of the treatment by time interaction is positive, Figure 8.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 C 2 (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 C 2 (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 C 2 (i, s ) due to the zero estimate of the treatment by time interaction. The negative estimates of α 5 and α 6 caused the combined component to drop from 0.694 (30.6% reduction) to 0.517 (48.3% reduction) over the 60 posttreatment months. Figure 8.2 shows that the combined treatment curve for the model without a jump is above that with a jump, albeit steeper. In fact, they decrease at constant monthly rates of -0.0095 and -0.0050, respectively. Hence, the large difference in the treatment components between the two hierarchical PLNC 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 & 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 for odds ratio decomposition leading to Figure 8.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 collisions immediately after the intervention. 117 Note that the parameters’ estimates (αˆ 4 , αˆ 5 , αˆ 6) along with the intervention dates ( t 0 i ) 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 odds ratio is set to OR i = OR i , 276 −t 0 i , the overall odds ratio for collisions could have been computed from (27). Under H2, the reduction in collision rate over the study period was estimated as 39.1% for all sites (Li et al., 2008). Figure 8.2 Monthly combined treatment components for a median site under H1 & H2, Li et al. (2008). 8.4.2 British Columbia 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, 118 thereby reducing deaths, injuries and insurance claim costs. Under this program a group of N T = 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 N C = 16 comparison groups comprising 98 intersections. The comparison groups range in size from 1 to 10 intersections. The sample size is n = 116 intersections. 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. The average injury counts per year before and after the treatment are summarized in Table 8.1 for the treated sites. The highest annual frequencies were observed for Site 10 followed by sites 6, 3, and 7. Of the 18 treated sites, thirteen sites have experienced reductions (negative changes) in collision injury ranging from 8-67%. In contrast, the remaining five sites have experienced increases ranging from 2-15%. Table 8.1 ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Annual injury collisions frequency for treated sites Treatment Year 2004 2004 2004 2004 2005 2005 2005 2005 2005 2005 2005 2005 2005 2006 2006 2006 2006 2006 Average counts (per year) Before After 12.00 13.25 35.00 26.75 80.67 54.00 32.33 14.50 7.75 4.67 95.75 85.00 57.50 53.00 31.50 21.33 28.50 22.00 126.25 144.33 5.75 4.67 20.75 14.00 36.75 31.67 12.20 14.00 6.00 2.00 4.40 4.50 7.20 4.00 7.00 7.50 %Change 10% -24% -33% -55% -40% -11% -8% -32% -23% 14% -19% -33% -14% 15% -67% 2% -44% 7% 119 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 for the parameters of the PLNI models under consideration were obtained via two chains with 10000 iterations 5000 of which were excluded as a burn-in sample using WinBUGS with a diffused N ( 0 , 100 2 ) prior for the regression parameters and an uninformative inverse- Gamma ( 0.01, 0.01 ) prior for the dispersion 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 gradual and jump PLNI models were outperformed by the corresponding random parameters PLNI-RP models with significant reductions in DIC of 130 (5652-5522) and 140 (5656-5516), respectively. Also, the estimate of the extra-Poisson variance σˆ ε2 = 0.391 under the gradual PLNI was reduced to 0.320 (18% reduction) under the gradual PLNI-RP. A similar reduction of 15% was obtained for the jump PLNI-RP relative to the jump PLNI. The random parameters jump PLNI-RP model was ruled out as the estimate of the sudden jump parameter (αˆ 6) was not significant. In fact, the estimate ( ± standard error) was − 0.071 ± 0.061 . Also, equations (26) and (75) gave close estimates for the odds ratios. In fact, the difference was within 0.025 and 0.001 at the site level and the overall level, respectively, under all models considered. For the present data, the posterior estimates were (αˆ 4 = −0.031, αˆ 5 = −0.058) under the gradual PLNI-RP model with random parameters among matched pairs. Thus, for those sites where the countermeasures were implemented in 2004 (t 0 i = 4) , the combined treatment 120 component (78a) simplifies to C T (i, s ) = exp{−0.1065 − 0.0445 s} = 0.956 C T (i, s − 1) , which is depicted in Figure 8.3 along with C1 (i, s ) and C 2 (i, s ) . Figure 8.3 Yearly treatment components under the gradual PLNI-RP model with random parameters among matched pairs (completion year 2004). Figure 8.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 C 2 (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 combined treatment component has been declining at a constant rate of -0.0445. This behavior is clearly manifested in Figure 8.3 as C T (i, s ) is decreasing steadily. The analogous curves for those sites where the countermeasures were implemented in 2005 and 2006 (t 0i = 5, 6) are quite similar. 121 It should be noted that the curves represented by C1 (i, s ) , C 2 (i, s ) and C T (i, s ) are linear in s only on the log-scale (although they look almost linear in Figures 8.1-8.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 between C1 (i, s ) and C 2 (i, s ) (Figure 8.1). In contrast, when negative estimates were obtained for both components, the combined treatment was compounded and its curve is below the two other curves (Figure 8.3). The ranges of the odds ratios 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 the two traffic volumes were ( γˆ 1 = 0.290 , γˆ 2 = 0.440 ) . Substituting these results in (77c) and (77d) led to the ranges (0.995, 1.003) and (0.984, 1.015) for C 3 (i, s ) and C 4 (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 gradual PLNI-RP model, the estimate of the overall odds ratio was 0.787 ± 0.041 with the 95% credible interval (0.709, 0.868) indicating a significant overall reduction of 21.3%. 8.5. Summary This chapter presents a novel approach to decompose the odds ratio into several components corresponding to direct and indirect treatment impacts. The decomposition of the odds ratio has yielded valuable insights in the effects of the model parameters on treatment impacts. The analysis revealed that the parameters (α 4 , α 5 , α 6) , which are related to the intervention time and date interactions, have a direct influence on the odds ratio and treatment effectiveness. For instance, if the time trend for the ith treated site is smaller than that for the matching comparison group, then α 4 representing treatment by time interaction would be negative and the odds ratio 122 will be reduced and vice versa. Alternatively, if the decrease in collision counts accelerates after treatment implementation then α 5 representing treatment by intervention date interaction would be negative and the odds ratio will be reduced and vice versa. In contrast, the parameters ( γ 1 ,γ 2 ) , which are related to the effects of major and minor traffic volumes, represent the indirect treatment impact. In fact, the decomposition results have shown that unless the traffic volumes are subject to significant annual fluctuations, their effect on the odds ratio and treatment effectiveness would be marginal. The decomposition of the odds ratio also provided a simple and intuitive way of computing the overall odds ratio without resorting to additional algorithms. Given that the model specification is unchanged, then the computation of the overall odds ratio requires the regression coefficients ( αˆ 4 ,αˆ 5 ,αˆ 6 , γˆ 1 ,γˆ 2 ) along with the intervention dates ( t 0 i ) and the two odds ratios for the major and minor traffic volumes. Moreover, the analyst would be able to compute the different components of the odds ratio 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 odds ratio and treatment effectiveness were computed and plotted for the two proposed 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 CBA 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 matched comparison-treatment pairs have significantly outperformed the regular PLNI models, while 123 reducing the estimates of the extra-Poisson variance. Also, since the sudden jump parameter was not significant, the random parameters PLNI model without a jump was preferred for analyzing the current data set. The estimate of the overall odds ratio was 0.787 ± 0.041 implying a significant 21.3% reduction. The odds ratios for each year following the intervention were decomposed into basic components associated with direct and indirect treatment effects manifested through the traffic volumes. 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. Overall, the results suggest that incorporating such design features as matched (yoked) comparison groups in SPFs leads to superior models. Moreover, the results shown in this chapter can have a significant impact on the economic evaluation of safety programs and countermeasures. Current evaluation techniques assume constant countermeasure effectiveness or Collision Modification Factors (CMFs) that do not change with time. This seems unrealistic. Given the way future benefits are discounted to present values, the results of using CMFs that change with time can affect the cost effectiveness of safety countermeasures significantly. Since the linear intervention model furnishes only some restricted treatment profiles, a non-linear intervention model affording a larger variety of treatment profiles is discussed in chapter 10. 124 9 MULTIVARIATE INTERVENTION MODELS This chapter illustrates the use of a multivariate intervention model to conduct a full Bayes CBA safety evaluation. 9.1 Background The incorporation of random parameters models and multivariate models into the analysis of SPFs were discussed in chapters 4 and 5, respectively. Moreover, Park et al. (2009) proposed a Multivariate Poisson-Lognormal Intervention (MVPLNI) model to account for the correlation among severity levels of collision counts. The use of random parameters models to account for the correlation between sites within comparison-treatment pairs was proposed in chapter 8. In this chapter, both approaches are combined to assess the factors that influence the frequency and severity of collisions in order to provide more effective safety-related countermeasures. Two models are considered for conducting a CBA safety evaluation. The first model uses a MVPLNI model for collision counts by severity-level with covariates representing traffic volume, time, treatment (intervention) and interaction effects. The second model extends MVPLNI by using random parameters to account for the comparison-treatment pairing. The objectives of this chapter are twofold: i) to compare the results of the MVPLNI with the extend model including random parameters in terms of goodness-of-fit and inference, and ii) to provide tools for investigating the safety change, in a multivariate context, in treated sites due to the implementation of certain types of countermeasures. The approach is demonstrated using a data set from British Columbia based on the ICBC’s RIP, where the full Bayes approach is utilized to determine the effectiveness of certain countermeasures that were implemented in 20 treated intersections in the Greater Vancouver area. 125 9.2 The Models 9.2.1 The Multivariate Poisson-Lognormal Intervention (MVPLNI) Model For a collision count of severity level k (k = 1, 2, …, K), let Y itk denote the observed count at site i during year t. Under the MVPLNI model, equations (1), (5) and (24) remain the same except for adding the superscript k to indicate the collision severity level. However, in order to account for the correlations among collision counts of different severity levels at site i, it is now assumed that ε i = ( ε 1i , ε i2 , ..., ε iK )' ~ N K ( 0 ,Σ ) , (79) where Σ is a covariance matrix (the diagonal element σ kk represents the variance of ε ik , whereas the off-diagonal element σ jk represents the covariance of ε ij and ε ik ). The multivariate model in equations (1), (5), (24) and (79) is similar to the model given in Park et al. (2009) for road segments, who noticed that it can account not only for correlations among different collision severity levels but also for correlations among repeated observations (collisions of the same severity level observed over different time periods at each site). It should be noted also that this multivariate model incorporates the correlations among K different severity levels into estimation, which was not possible in previous univariate intervention models. 9.2.2 The MVPLNI Model with Random Parameters among Pairs (MVPLNI-RP) Under the MVPLNI model with random parameters among matching comparison-treatment pairs, equations (1), (5), (72) and (73) remain the same except for adding the superscript k to indicate the collision severity level. Additionally, the MVPLNI-RP model comprises Equation (79) to account for the correlations among collision counts of different severity levels. 126 Further, the random parameters model MVPLNI-RP permits a pair-level inference, using Equation (72) and an overall-level inference using Equation (24). Note that for sudden (with jump) MVPLNI models, equations (24) and (72) are replaced by (25) and (74). 9.3 Prior Distributions The following uninformative priors were used: β j ~ N ( 0 , 1002 ) , σ −j 2 ~ Gamma ( 0.001, 0.001 ) and Σ −1 ~ Wishart ( P , 2 ) , where P = I the 2 x 2 identity matrix. Yet, to examine the sensitivity of the results to the choice P = I (no correlation), the alternative P = ( Pi , j ) , where Pii = 0.10 and Pij = −0.05 , i ≠ j , (larger variances and an intermediate correlation of 0.5) is also considered. 9.4 The Data Under the ICBC’s RIP a number of road safety initiatives including general signal visibility improvements, left turn phase improvement and left turn lane installation have been funded in the Greater Vancouver area. Improving signal visibility includes using such treatments as larger signal heads and reflective backboards. Thus, a group of N T = 20 intersections (sites) were selected for treatment within a before-after experiment. Each treated intersection was matched with a comparison group of intersections satisfying the recommended similarity and proximity conditions. The treated intersections were matched with N C = 16 comparison groups comprising 98 intersections. The comparison groups range in size from 1 to 10 intersections. The sample size is n = 118 intersections. The countermeasures were implemented during the years 2004, 2005 or 2006. The treated sites are described in Table 9.1. Data on collision counts by severity levels along with the corresponding traffic volumes were collected for the years 2001-2008. This period provides adequate time for a proper beforeafter safety evaluation. The K = 2 severity levels were injuries and fatalities (I+F) and property damage only (PDO). 127 Table 9.1 Selected treatment sites Site City 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Vancouver North Vancouver City Coquitlam Coquitlam Port Coquitlam District of North Vancouver Vancouver Vancouver Vancouver Vancouver Vancouver Vancouver Coquitlam Coquitlam Vancouver Vancouver Vancouver District of North Vancouver Coquitlam Coquitlam Treat. Group 1 1 2 3 1 1 2 2 2 2 2 2 1 3 3 1 1 1 1 3 Description Improving signal visibility Improving signal visibility Left Turn Phase Improvement Left Turn Lane Installation Improving signal visibility Improving signal visibility Left Turn Phase Improvement Left Turn Phase Improvement Left Turn Phase Improvement Left Turn Phase Improvement Left Turn Phase Improvement Left Turn Phase Improvement Improving signal visibility Left Turn Lane Installation Left Turn Lane Installation Improving signal visibility Improving signal visibility Improving signal visibility Improving signal visibility Left Turn Lane Installation Treat. Year 2004 2004 2004 2004 2004 2005 2005 2005 2005 2005 2005 2005 2005 2005 2005 2006 2006 2006 2006 2006 The average collision counts per year before and after the treatment are summarized in Table 9.2 by severity level for the treated sites. For both severity levels, the highest annual collision frequency was observed for Site 5 followed by sites 11, 12, 7, 3 and 8. Also, except for site 1, the I+F annual collision frequency have been reduced after treatment. In fact, the observed reductions (negative changes) were -27%, -27% and -40% for treatments 1, 2, and 3, respectively. The overall reduction in the I+F annual collision frequency was -28%. In contrast, there were positive changes (increases) in the PDO annual collision frequency for eight of the twenty treated sites. The overall reduction in the PDO annual collision frequency was -6%. 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 treatments and parameters on the severity of collision counts is conducted using a random parameters MVPLNI model as formulated above. 128 Table 9.2 Annual collision frequency by severity level for treated sites Site Treat. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 1 2 3 1 1 2 2 2 2 2 2 1 3 3 1 1 1 1 3 9.5 Average I+F collisions per year Before After %Change 4.3 5.5 27 13.7 10.3 -25 35.7 19.3 -46 14.7 6.5 -56 99.3 72.5 -27 2.5 2.0 -20 38.8 23.7 -39 26.0 18.3 -29 12.8 10.3 -19 12.3 7.0 -43 42.8 35.3 -17 39.8 38.3 -4 2.3 1.3 -41 11.3 5.3 -53 12.3 11.3 -7 4.8 3.0 -38 1.6 0.0 -100 2.6 1.5 -42 3.0 1.5 -50 3.0 1.5 -50 Average PDO collisions per year Before After %Change 7.7 7.8 1 21.3 16.5 -23 45.0 34.8 -23 17.7 8.0 -55 147.7 104.8 -29 5.3 2.7 -49 57.0 61.3 8 31.5 34.7 10 18.8 11.0 -41 16.3 15.0 -8 83.5 109.0 31 58.0 72.0 24 3.5 3.3 -5 9.5 8.7 -9 24.5 20.3 -17 7.4 11.0 49 4.4 2.0 -55 1.8 3.0 67 4.2 2.5 -40 4.0 6.0 50 The Odds Ratios For s = 1, 2, 3 , let S ( s ) denote the subset of intersections that were subject to treatment s and let N s denote the number of intersections in this subset. For instance, note that S (1) = {1, 2, 5, 6, 13, 16, 17, 18, 19} , N 1 = 9 and so forth. Then, the treatments’ odds ratios can be computed from ln( OR s ) = ( 1 / N s )∑i∈S ( s ) ln( OR i ) . (80) The treatment effect is calculated from ( OR s − 1 ) , while the percentage of reduction in predicted collision counts, also called the Collision Reduction Factor (CRF), is given by ( 1 − OR s ) × 100. 129 For the random parameters models, two estimates of the overall odds ratio are available for analysis. The first estimate makes use of Equation (27), while the second estimate is obtained from OR = µ TA µ CB µ TB µ CA , (81) where the predicted collision counts are computed from Equation (24) and then averaged over appropriate sites and years. 9.6 Results and Discussion 9.6.1 Model Estimates The posterior summaries in Table 9.3 were obtained via two chains with 20000 iterations 10000 of which were excluded as a burn-in sample using WinBUGS. The major and minor AADT statistics were centered (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.065 and trace plots for all model parameters indicated convergence. The autocorrelation functions died out reasonable well, particularly for the intervention parameters. The alternative choices for the matrix P of the Wishart prior led to virtually the same posterior estimates of the regression parameters with differences in the third decimal place. Table 9.3 summarizes the parameter estimates and their standard errors for MVPLNI and MVPLNI-RP. For I+F (k=1) the estimates of the time trend ( αˆ 2 ) , intervention slope ( αˆ 3 ) and their differences ( αˆ 4 and αˆ 5 ) were negative under both models, except the non-significant intervention slope difference ( αˆ 5 ) under MVPLNI. Yet, these estimates are not significantly different from zero, except for (i) the intervention slope ( αˆ 3 ) indicating a significant decline in I+F collision counts subsequent to the intervention under both models, and (ii) the time trend 130 difference ( αˆ 4 ) indicating a significantly larger decline rate in I+F collision counts for treated sites than that for comparison sites under MVPLNI. The results for PDO (k=2) show that the estimates of the time trend ( αˆ 2 ) , intervention slope ( αˆ 3 ) and their differences ( αˆ 4 and αˆ 5 ) are not significantly different from zero under both models, except for the intervention slope difference ( αˆ 5 ) under MVPLNI. The estimate αˆ 5 = −0.106 indicates a significantly larger decline rate in PDO collision counts for treated sites than that for comparison sites subsequent to the intervention. Table 9.3 Parameter estimates and standard errors MVPLNI α0 σ α0 α1 σ α1 α2 σα2 α3 σα3 α4 σα4 α5 σα5 γ1 σγ1 γ2 σγ2 σ kk σ 12 ρ OR DIC MVPLNI-RP I+F (k=1) PDO (k=2) I+F (k=1) PDO (k=2) -6.782± 0.981 -7.051± 0.934 -5.164±1.517 2.592±1.114 -6.602±1.260 2.300±0.744 0.445± 0.180 0.117± 0.160 0.410±0.162 0.170±0.083 0.152±0.172 0.167±0.076 -0.012±0.009 -0.019±0.006 -0.011±0.019 0.058±0.014 -0.010±0.018 0.059±0.013 -0.047±0.018 0.020±0.012 -0.073±0.033 0.085±0.027 0.010±0.033 0.101±0.029 -0.045±0.020 0.032±0.014 -0.058±0.033 0.087±0.024 -0.004±0.032 0.088±0.023 0.010±0.039 -0.106±0.027 -0.016±0.060 0.113±0.041 -0.079±0.053 0.112±0.042 0.467±0.097 0.449±0.090 0.323±0.144 0.218±0.118 0.397±0.109 0.165±0.075 0.456±0.054 0.566±0.051 0.443±0.072 0.133±0.059 0.578±0.075 0.170±0.070 0.440±0.081 0.422±0.075 0.366±0.071 0.363±0.065 0.394±0.073 0.914±0.021 0.824±0.037 0.923±0.031 9750 0.330±0.062 0.906±0.024 0.771±0.054 0.846±0.047 9652 131 For I+F the estimates and standard deviations for σ 11 were 0.440 ± 0.081 and 0.366 ± 0.071 under the two models. Hence, the estimate of 0.440 under MVPLNI has been reduced by 17% under MVPLNI-RP. For both models, the standard deviations are small relative to the parameter estimates and the lower limits of the credible intervals are noticeably greater than zero demonstrating the presence of over-dispersion in the I+F collision counts. For PDO the estimates and standard deviations for σ 22 were 0.422 ± 0.075 and 0.363 ± 0.065 under the two models. Hence, the estimate of 0.422 under MVPLNI has been reduced by 14% under MVPLNIRP. These estimates demonstrate the presence of over-dispersion in the PDO collision counts as well. As stated earlier, these findings confirm the conclusions of Miaou and Song (2005) and Aul and Davis (2006). The former authors have suggested that careful modeling of the random effects may reduce the unexplained over-dispersion typically observed in collision data, whereas the latter authors have argued that over-dispersion may to a large extent be an artifact of mixing fundamentally heterogeneous entities and that it seems prudent to avoid dependence on methods that require the observation of over-dispersion before analysis can proceed. The two multivariate models account for the correlation between the latent variables θ 1it (I+F) and θ it2 (PDO). The estimate of this correlation ( ρ ) was 0.91, which is quite high and is obviously highly significant. In essence, higher PDO collisions are associated with higher injury collisions, as the collision likelihood for both levels is likely to rise due to the same deficiencies in roadway design and/or other unobserved factors. But equations (24) and (72) may suffer from omitted variables bias since the unobserved heterogeneity from other factors known to influence collision counts ends up in the correlation structure, thus affecting the estimated correlation. However, it has been shown earlier that the current estimate of 0.91 is in the upper range. Although the inclusion of additional relevant covariates is expected to reduce the correlation, it is conjectured to remain significant with crucial inferential impact. 132 The estimates ( σˆ j ) were significant under MVPLNI-RP for all covariates including the intercept suggesting that the incorporation of random parameters among pairs in the MVPLNI model is warranted. Further, according to the DIC guidelines (Spiegelhalter et al., 2005), the random parameters model MVPLNI-RP has fitted the data much better than the regular MVPLNI model with a significant reduction in DIC of 98. Under MVPLNI-RP, the Pearson’s residuals were examined and no anomalies were detected. The posterior estimates of the observed chi-square statistics were 537.4 and 475.4 for I+F and PDO, respectively. The associated p-values estimated from the distributions of the chisquare discrepancy measure in replicated datasets were 0.369 and 0.128, respectively. Accordingly, MVPLNI-RP seems to perform well in terms of accommodating the variation in collision frequency across intersections. The estimates of the I+F overall odds ratio were 0.824 ± 0.037 and 0.771 ± 0.054, implying significant reductions in predicted I+F collision counts of 18% and 23% under MVPLNI and MVPLNI-RP, respectively. In contrast, the estimates of the PDO overall odds ratio were 0.923 ± 0.031 and 0.846 ± 0.047, implying reductions in predicted PDO collision counts of 8% and 15% under MVPLNI and MVPLNI-RP, respectively. Recall that two estimates of the overall OR can be obtained under MVPLNI-RP using either Equation (27) or Equation (81). The former estimates are reported in Table 9.3, while the latter estimates were 0.754 ± 0.109 and 0.848 ± 0.117, for I+F and PDO, respectively. Although the point estimates agree closely, the standard errors of the second approach are just about double those of the first approach. 9.6.2. The Site-Level and Treatment-level Odds Ratios In addition to the overall estimate of ( α 0 ,α 1 , ...,α 5 ,γ 1 ,γ 2 ) , provided in Table 9.3, there were additional 16 sets of estimates under MVPLNI-RP corresponding to the 16 treatmentcomparison pairs. Such estimates were used in equations (26) to compute the site-level odds 133 ratios. The results are reported in Tables 9.4, which reveals that the site-level odds ratios were less precisely estimated than the overall OR. Table 9.4 Site-level odds ratios under MVPLNI-RP I+F a b Site OR SE 95% C. I. OR SE. 1 0.994 0.248 (0.591, 1.565) 0.992 0.220 2 0.768 0.140 (0.531, 1.084) 0.695 0.105 a 3 0.746 0.084 (0.585, 0.917) 0.622 0.062 0.525 0.095 4 0.650 0.124 (0.437, 0.926)a 5 0.862 0.098 (0.680, 1.062) 0.733 0.069 6 0.806 0.255 (0.404, 1.393) 0.736 0.201 7 0.844 0.092 (0.676, 1.044) 1.037 0.088 8 0.765 0.084 (0.615, 0.946)a 0.897 0.078 0.909 0.079 9 0.772 0.084 (0.621, 0.954)a 10 0.681 0.136 (0.447, 0.976)a 0.916 0.152 11 1.084 0.125 (0.861, 1.346) 1.364 0.106 12 1.185 0.125 (0.965, 1.448)b 1.144 0.090 13 0.766 0.288 (0.352, 1.455) 0.767 0.238 14 0.650 0.146 (0.406, 0.970)a 0.741 0.159 15 1.116 0.118 (0.907, 1.367)b 1.052 0.085 16 0.708 0.197 (0.391, 1.161) 1.376 0.310 17 0.500 0.198 (0.189, 0.951)a 0.626 0.189 18 0.813 0.276 (0.398, 1.475) 1.185 0.391 a 19 0.565 0.185 (0.268, 0.982) 0.551 0.169 20 0.898 0.268 (0.479, 1.520) 1.034 0.260 Overall 0.771 0.054 (0.669, 0.884)a 0.846 0.047 The OR is significant at the 0.05 level of significance. The 95% confidence interval is not overlapping that for the overall OR. PDO 95% C. I. (0.639, 1.494) (0.513, 0.928)a 0.510, 0.747)a,b (0.361, 0.731)a,b (0.611, 0.876)a (0.402, 1.187) (0.875, 1.222) (0.754, 1.063) (0.765, 1.076) (0.652, 1.253) (1.169, 1.579)b (0.976, 1.333)b (0.406, 1.325) (0.484, 1.105) (0.896, 1.227) (0.862, 2.081) (0.323, 1.059) (0.604, 2.123) (0.278, 0.950)a (0.610, 1.631) (0.757, 0.944)a Table 9.4 shows also that the site-level I+F odds ratios ranged between 0.500 and 1.185 and that there were 17 sites exhibiting reductions in predicted I+F collision counts. However, these reductions were significant for only 8 sites, as the upper bounds of the corresponding 95% credible intervals are less than 1. In contrast, the site-level PDO odds ratios ranged between 0.525 and 1.376 and there were 13 sites exhibiting reductions, of which only 5 were significant. Accordingly, the effectiveness of the treatment seems to vary by severity level from one location to another. Furthermore, the 95% credible intervals for sites 12 and 15 didn’t overlap that of the I+F overall OR, whereas the 95% credible intervals for sites 3, 4, 11 and 12 didn’t overlap that of the 134 PDO overall OR. Thus, these sites need to be investigated as some other factor(s) may be interacting with the implemented treatment. Equation (80) was used to compute the treatment-level odds ratios that are reported in Table 9.5, which reveals that there were reductions in both I+F and PDO collision counts for all three treatments (all odds ratios are less than 1). However, the reductions in PDO collision counts were not significant for left turn phase improvement. Table 9.5 Treatment-level odds ratios under MVPLNI-RP I+F Treatment OR SE 95% C. I. Improving signal 0.713 0.086 (0.554, 0.892)a visibility Left Turn Phase 0.847 0.054 (0.746, 0.958) a Improvement Left Turn Lane 0.794 0.089 (0.635, 0.982)a Installation Overall 0.771 0.054 (0.669, 0.884)a a The OR is significant at the 0.05 level of significance. PDO OR SE. 95% C. I. 0.794 0.077 (0.655, 0.956) a 0.956 0.048 (0.869, 1.054) 0.796 0.079 (0.650, 0.961) a 0.779 0.096 (0.757, 0.944)a The collision reduction factors (CRFs) were (29%, 21%), (15%, 4%) and (21%, 20%) for improving signal visibility, left turn phase improvement and left turn lane installation, respectively, where the first percentage corresponds to I+F collisions and the second to PDO. Thus, improving signal visibility has resulted in the largest CRFs, followed by the installation of a left turn lane. Left turn phase improvement has resulted in a significant CRF for I+F collisions, but not for PDO collisions. 9.7 Summary In this chapter, two models were considered for conducting CBA evaluation of the collision counts at the intersections improved by the installation of countermeasures under the ICBC’s RIP. The first model is a multivariate Poisson-lognormal intervention (MVPLNI) model for collision counts by severity-level with covariates representing traffic volume, time, treatment (intervention) and interaction effects. In the second model (MVPLNI-RP), random parameters 135 for comparison-treatment pairs are incorporated to deal with the between-pair variation. The results of the two models were compared using the full Bayes approach, within the framework of a before-after setup with matched (yoked) comparison groups. The results show that the variation among pairs has significantly improved the fit while reducing the estimates of the extra-Poisson variation. Further, the effects of the covariates on collision counts were found to vary significantly across treatment-comparison pairs justifying the use of random parameters. Thus, for this data set, the random parameters model MVPLNI-RP is preferred over MVPLNI. The correlation between the two severity levels (I+F and PDO) of collision counts was highly significant indicating that higher PDO collisions are associated with higher I+F collisions, as the collision likelihood for both levels is likely to rise due to similar deficiencies in roadway design and/or other unobserved factors. This correlation was identified by other researchers in the literature. According to MVPLNI-RP, the estimates of the overall odds ratio were 0.77 and 0.85 implying significant reductions in predicted collision counts of 23% and 15% for I+F and PDO, respectively. Therefore, the countermeasures seem to be more effective in reducing I+F collision counts more than PDO collision counts. It should be also emphasized that the overall odds ratio was more precisely estimated as an average of the site-level odds ratios using Equation (27) than collectively as formulated in Equation (81). Although the majority of the site-level odds ratios exhibited reductions in both I+F and PDO predicted collision counts, only some of these reductions were significant. Also, the effectiveness of the treatment seems to vary by severity level from one location to another. In addition, the odds ratios for few sites were not consistent with the overall odds ratio and deserve further investigation. 136 For both severity levels, the highest collision reduction factors were those of improving signal visibility followed by left turn lane installation. For left turn phase improvement the CRF was significant for I+F but not for PDO. The results suggest that incorporating such design features as matched (yoked) comparison groups in SPFs leads to superior models. Apart from the improvement in the goodness of fit, such extended models can be used to account for heterogeneity due to unobserved road geometrics, traffic characteristics, environmental factors and driver behavior. 137 10 NON-LINEAR INTERVENTION MODELS: THE KOYCK MODEL This chapter illustrates the use of the Koyck model as an alternative non-linear approach to conduct a full Bayes CBA safety evaluation. 10.1 Background In chapter 8, the use of Poisson-Lognormal intervention (PLNI) model in the context of a CBA safety evaluation was investigated. The approach uses dynamic regression models to identify the lagged effects of the intervention in order to describe their response and acknowledges the fact that the treatment (intervention) effects do not occur instantaneously but are rather spread over future time periods. In all SPFs considered, linear slopes were tacitly assumed to represent the time and treatment effects across the treated and comparison sites. Moreover, the odds ratio was decomposed under the linear intervention SPF into basic components related to direct and indirect treatment effects. The consequences of the linearity assumption on these components and subsequently on the OR were highlighted. It was noted that the linear slopes of the intervention model furnishes only some restricted treatment profiles which can impact the estimation of CMFs and the economic evaluation of safety programs and countermeasures. Alternatively, the Koyck model (Koyck, 1954) could be used to represent the lagged treatment effects that are distributed over time. The Koyck model is an alternative dynamic regression form involving a first-order autoregressive (AR1) model that is based on distributed lags (Judge et al., 1988; Pankratz, 1991). The Koyck model affords a rich family of forms (over the parameter space) that can accommodate various profiles for the time and treatment effects. A byproduct of the Koyck models is its ability to isolate an additional component of OR corresponding to the time (novelty) effects. The novelty phenomenon was first explored by Persaud (1986) who examined the decline in safety that is thought to accommodate a change in traffic control (two-stop, all stop, signalized etc.). Analyzing such a phenomenon provides valuable insights on whether the overall gain in safety more than compensates for the short- 138 term confusion caused by the introduction of the new safety countermeasure. Albeit its importance the novelty phenomenon has received little to no attention in the road safety literature and therefore it is explored in this chapter. The primary objective of this chapter is to demonstrate the effectiveness of the Koyck model to conduct a FB CBA safety evaluation. To achieve this objectives the following issues will be addressed: i) challenging the linear slopes assumption commonly used to represent the time and treatment effects across the treated and comparison sites statistically and logically, ii) providing a non-linear alternative to represent the lagged treatment effects that are distributed over time, iii) comparing the results of the non-linear Koyck model to the conventional linear intervention model, and iv) decomposing the odds ratio under the Koyck model to understand the role of the model parameters and their effects on the OR as well as providing straightforward equations for the computation of the OR and its components. The approach is demonstrated using a data set from British Columbia based on the ICBC’s RIP. 10.2 The Koyck Model In all intervention models studied in the literature, the treatment effects were modeled via linear trend and intervention slopes (with a possible jump). Alternatively, the consequences of the intervention can be modeled using distributed lags along with a first-order autoregressive (AR1) model as a proxy for the time effects (Judge et al., 1988; Pankratz, 1991). Recall that the occurrence of the intervention is indicated by the step function I it . The regression equation for the rational distributed lag model is given by ln( µ it ) = α 0 + α 1 T i + [ ω /( 1 − δ B )] I it + [ ω* /( 1 − δ B )] T i I it + γ 1 ln( V 1it ) + γ 2 ln( V 2it ) + γ 3 X 3i + ... + γ J X Ji + ν t , (82a) where B denotes the backshift operator ( B Z t = Z t −1 ) , δ < 1 and ν t satisfies the following stationary AR1 equation ν t = φ ν t −1 + et , φ < 1, 2 et ~ N (0, σ ν ) , t = 2 , 3 , ..., m . (83). 139 −1 Consider the expansion ( 1−δ B ) I it = I it + δ I i ,t −1 + δ 2 I i ,t −2 + ... , and note that the rational distributed lag model depicts an everlasting treatment effect as ln(µ it ) is tacitly assumed to be a function of the infinite distributed lags ( I it , I i ,t −1 , I i ,t −2 , ...) . The parsimonious model (82a) is known as the Koyck model (Koyck, 1954) in which the lag weights ω δ k and ω* δ k decline geometrically for k = 0 , 1, 2 , ... . Consequently, the earlier years following the intervention are more heavily weighted than distant years. It should also be noted that although the weights never reach zero, they will eventually become negligible. The two parameters ω and ω* are impact multipliers, whereas δ is a decay parameter controlling the rate at which the weights decline. Rewriting Equation (82a) as ln(µ it ) = C it + ν t , the AR1 Equation (83) implies that ν t = φ [ln(µ i ,t −1) − C i ,t −1] + et . Substituting this last expression in (82a) leads to ln( µ it ) = ( 1 − φ ) α 0 + ( 1 − φ )α 1 T i + [ ω /( 1 − δ B )] I *it + [ ω* /( 1 − δ B )] T i I *it + γ 1 X 1it + γ 2 X 2it + ( 1 − φ ) γ 3 X 3i + ... + ( 1 −φ )γ J X Ji + φ ln( µ i ,t −1 ) + et , where * I it = I it − φ I i ,t −1 , and X jit = ln( V jit ) − φ ln( V ji ,t −1 ) , (82b) j = 1, 2 . Finally, applying the operator (1 − δ B ) to both sides of (82b) yields ln(µ it ) = (1 − φ )(1 − δ ) α 0 + (1 − φ )(1 − δ ) α 1 T i + ω I *it + ω * T i I *it + γ 1 X *1it + γ 2 X *2 it + ( 1 − φ )( 1−δ )γ 3 X 3i + ... + ( 1−φ )( 1−δ )γ J X Ji + ( φ + δ ) ln( µ i ,t −1 ) − φ δ ln( µ i ,t − 2 ) + et , where X *jit = X jit (82c) − δ X ji ,t −1 , j = 1, 2 . Equation (82c) holds for t = 3 , 4 , ..., m . The regression models for the first two years are given by ln( µ i 2 ) = ( 1 − φ )α 0 + ( 1 − φ )α 1 T i + γ 1 [ ln( V 1i 2 ) − φ ln( V 1i 1 )] + γ 2 [ ln( V 2 i 2 ) − φ ln( V 2 i1 )] + ( 1 − φ ) γ 3 X 3i + ... + ( 1−φ )γ J X Ji + φ ln( µ i 1 ) + e2 , and 140 ln( µ i1 ) = α 0 + α 1 T i + γ 1 ln( V 1i 1 ) + γ 2 ln( V 2i 1 ) + γ 3 X 3i + ... + γ J X Ji + ν 1 , where ν 1 ~ N ( 0 ,σ ν2 /( 1 − φ )) . To derive the variance of ν 1 , the AR1 Equation (83) implies 2 that var(ν t ) = φ var(ν t −1 ) + σ ν2 . For φ < 1 (stationary AR1), var(ν t ) = σ ν2 /( 1 − φ ) , for all t. 2 2 The Koyck model’s version with random parameters among matched comparisontreatment pairs is denoted Koyck-RP. It is important to check the appropriateness of such models for a given dataset by monitoring the posterior probabilities of the stationary conditions ( δˆ ≤ 1) and ( φˆ ≤ 1) . Note also that the Koyck model is intrinsically non-linear (in the parameters). 10.3 Prior Distributions The uninformative priors N (0, 1002) and Gamma (0.001, 0.001) were used for the regression and precision parameters, respectively, while the uniform prior U(-1,1) was used for φ and δ . The uniform prior U(-1,1) is usually used for the decay parameter δ and the AR1 parameter φ in line with stationarity. The posterior probability of non-stationarity ( φ ≥ 1) should be monitored and a N (0, τ ) prior can be used if stationarity is not imposed, where τ is small, e.g., 1 or 0.5 (Congdon, 2006). 10.4 Decomposing the Odds Ratio under the Koyck Model Consider the Koyck model in Equation (82b) and note that the regression parameters ( α 0 ,α 1 ,ω ,γ 3 , ...,γ J ) cancel out of (75). Thus, for the ith treated site, the odds ratio after s postintervention years ( s = 1,2, ... ) can be decomposed (see Appendix C) as follows OR is = K 1 ( i , s ) K 2 ( i , s ) K 3 ( i , s ) K 4 ( i , s ) , (84) where φ s −1 φ K 1 ( i , s ) = [ ORi ,s −1 ] = [ ORi 1 ] , (85a) 141 * s * K 2 ( i , s ) = exp{ c + d ( 1 − δ ) / s } , c = ω ( 1 − φ ) /( 1 − δ ) , d = ω ( φ − δ ) / ( 1−δ ) , (85b) 2 φγ γ K 3 ( i , s ) = [ OR( V 1is )] 1 / [ OR( V 1i ,s −1 )] 1 , γ K 4 ( i , s ) = [ OR( V 2 is )] 2 / [ OR( V 2i ,s −1 )] φγ2 . (85c) (85d) Equations (85a)-(85d) which were obtained after lengthy algebraic manipulations, hold also for the Koyck-RP model (with the random parameters among pairs) at the overall-level of inference. The component K 1 ( i , s ) corresponds to the time (novelty) effects. After the first year of intervention, the subsequent odds ratios would either grow or decline exponentially at a rate of φ according to whether OR i1 < 1 or OR i1 > 1 . In both cases, K 1 ( i , s ) converges to 1 (since φ < 1 ). To illustrate, the time component is portrayed in Figure 10.1 where ORi1 = 0.80 . For φ > 0 , K 1 ( i , s ) converges to 1 from below and the convergence is slower for larger values of φ . Alternatively, for φ < 0 , K 1 ( i , s ) oscillates around 1 and the fluctuations are larger for smaller values of φ . Figure 10.1 Yearly novelty effects under the Koyck model. 142 The treatment component (85b) describes a non-linear function of s involving the impact multiplier ω * along with the AR1 parameter φ and the decay parameter δ . In the long-run K 2 ( i , s ) converges to exp{c} , which corresponds to the everlasting (permanent) treatment impact. To understand its behavior, K 2 ( i , s ) is shown in Figures 10.2 and 10.3 for various combinations of φ and δ , where ω* = −0.223 corresponding to K 2 ( i ,1 ) = 0.80 . Figure 10.2 shows the curves representing δ = −0.5 , − 0.2 , 0.2 , 0.5 for φ = 0.4 . Thus, starting with K 2 ( i ,1 ) = 0.80 , K 2 ( i , s ) decreases as δ increases for s > 1 . Thus, the curves for δ = 0.2 , 0.5 show smaller values of K 2 (i, s ) compared with those for δ = − 0.5 , − 0.2 . Also, while K 2 ( i , s ) increases temporarily for negative and small positive values of δ , the opposite holds true for δ = 0.5 . In all cases, K 2 ( i , s ) stabilizes within few years following the intervention. Figure 10.3 shows the curves representing φ = −0.4 , − 0.2 , 0.2 , 0.4 for δ = 0.3 . Thus, starting with K 2 ( i ,1 ) = 0.80 , K 2 ( i , s ) increases with φ for s > 1 . Yet, the curves for φ = − 0.4 , − 0.2 , show rather small values of K 2 ( i , s ) compared with those for φ = 0.2 , 0.4 . Also, while K 2 ( i , s ) declines temporarily for negative and small positive values of φ , the opposite holds true for φ = 0.4 . In all cases, K 2 ( i , s ) stabilizes within few years following the intervention, but the convergence is a little slower for φ < 0 . 143 Figure 10.2 Yearly treatment effects by δ at φ = 0.4, under the Koyck model. Figure 10.3 Yearly treatment effects by φ at δ = 0.3, under the Koyck model. 144 The two components K 3 ( i , s ) and K 4 ( i , s ) represent the effects of the major and minor traffic volumes, respectively. For both components, the numerator is the current traffic volume odds ratio raised to a fractional power ( γ 1 or γ 2 ) and thereby would be close to 1. Yet, the denominator would be even closer to 1 as the power of the previous year’s odds ratio is much smaller ( φ γ 1 or φ γ 2 ). Thus, unless the traffic volumes are subject to significant annual fluctuations, both components are expected to be near 1. The components K 3 ( i , s ) and K 4 ( i , s ) can be combined to obtain γ γ K V ( i , s ) = [ OR( V 1is )] 1 [ OR( V 2 is )] 2 /{ [ OR( V 1i ,s −1 )] φγ1 [ OR( V 2 i ,s −1 )] φγ2 }. (86) The combined component in (86) is inversely related to the indirect (through traffic volumes) treatment impact under the Koyck model. 10.5 The Data Under the ICBC’s RIP a group of N T = 25 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 N C = 14 comparison groups comprising 91 intersections. The comparison groups range in size from 3 to 10 intersections. The sample size is n = 116 intersections. The countermeasures were implemented during the years 2002 or 2003. Data on collision injury counts along with the corresponding traffic volumes were collected for the years 1999-2008. The average collision injury counts per year before and after the treatment are summarized in Table 10.1 for the treated sites. The highest yearly injury rates were observed for site 15 followed by sites 13, and 6. Of the 25 treated sites, about two thirds (16 sites) have experienced reductions (negative changes) in collision injury rates ranging from 9-70%. In contrast, the remaining 9 sites have experienced increases ranging from 4-47%. This preliminary analysis shows that the treatments may have an impact on the safety of these sites 145 (intersections). A rigorous before-after analysis of the effects of different treatments and parameters on collision injury counts is conducted in the following sections. Table 10.1 ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 Yearly collision injury rates for treated sites Treatment Year 2002 2003 2003 2003 2002 2002 2002 2003 2003 2002 2002 2002 2003 2003 2003 2003 2003 2003 2003 2003 2002 2003 2002 2003 2003 Average counts (per year) Before After 14.0 12.7 18.5 23.6 12.5 3.8 17.0 11.6 14.0 10.0 26.3 15.0 8.3 10.2 3.8 3.4 9.3 8.0 20.3 16.8 12.0 14.7 9.7 12.7 29.3 32.4 10.8 15.8 98.0 85.0 8.8 5.4 10.5 7.4 11.5 6.6 6.3 4.4 4.0 4.6 13.7 10.0 21.5 19.4 19.0 19.7 6.0 5.0 4.3 6.2 %Change -10 28 -70 -32 -29 -43 22 -9 -14 -17 22 31 11 47 -13 -38 -30 -43 -30 15 -27 -10 4 -17 46 10.6 Results and Discussion The posterior estimates of the parameters of the PLNI and Koyck models were obtained via two chains with 10000 iterations 5000 of which were excluded as a burn-in sample using WinBUGS. The BGR statistics, the ratios of the Monte Carlo errors relative to the standard deviations of the estimates and trace plots for all model parameters indicated convergence. The PLNI “jump” models were ruled out as the estimates of the sudden jump parameter (αˆ 6) were not significant and the differences in DIC were negligible. 146 The gradual PLNI-RP model (with random parameters among comparison-treatment pairs) has outperformed the regular PLNI model with a significant reduction in DIC of 161 (68586697). Also the estimate of the extra-Poisson variance σˆ ε2 = 0.335 under PLNI has been reduced to 0.137 (59% reduction) under PLNI-RP. Similarly, the random parameters model Koyck-RP has outperformed the Koyck model with a significant reduction in DIC of 167 (6851-6684) and 57% reduction in the estimate of the extra-Poisson variance. Yet, the Koyck-RP model provided the best fit with a DIC of 6684 compared to 6697 for PLNI-RP; a significant reduction of 13 according to the DIC guidelines. But there is more to the comparison of the two models than goodness-of-fit as they provide quite different structures for the treatment impact on collisions as discussed above and shown below. Also, equations (26) and (75) gave close estimates for the odds ratios. In fact, the difference was within 0.023 and 0.015 at the site level under the random parameters PLNI and Koyck models, respectively. At the overall level, the corresponding differences were even smaller (within 0.003). 10.6.1 Treatment Impact under PLNI Models For the present data, the posterior estimates were ( αˆ 4 = 0.031, αˆ 5 = −0.068 ) under the gradual PLNI-RP model with random parameters among matched pairs. Thus, for those sites where the countermeasures were implemented in 2002 ( t 0 i = 4 ) , the treatment component (78a) simplifies to C T ( i , s ) = exp{ 0.0435 − 0.0185 s } = 0.982 C T ( i , s − 1 ) , which is depicted in Figure 10.4 along with C 1 ( i , s ) and C 2 ( i , s ) . Figure 10.4 shows that the treatment by time component C1 (i, s ) has increased gradually from 1.097 in 2003 to 1.186 in 2008. Correspondingly, the treatment by intervention date component C 2 (i, s ) has decreased from 0.934 to 0.788. As a result, the combined treatment component showed a decline from 1.025 (2.5% increase) at the end of 2003 to 0.935 147 (6.5% reduction) at the end of 2008. The analogous curves for those sites where the countermeasures were implemented in 2003 (t 0i = 5) are quite similar. Figure 10.4 Yearly treatment components under the PLNI-RP model (completion year 2002). The ranges of the odds ratios for the major and minor traffic volumes were 0.83-1.16 and 0.86-1.26, respectively, and the posterior estimates of the regression coefficients associated with the two traffic volumes were ( γˆ 1 = 0.233, γˆ 2 = 0.409 ) . Substituting these results in (77c) and (77d) led to the ranges 0.96-1.03 and 0.94-1.10 for C 3 ( i , s ) and C 4 ( i , s ) , respectively, which are closer to 1 as discussed earlier. Furthermore, the combined traffic volume component C V ( i , s ) ranged from 0.92 to 1.11. Under the gradual PLNI-RP model, the estimate of the overall odds ratio was 0.935 ± 0.042 with the 95% credible interval (0.856, 1.018) indicating that the apparent overall reduction of (6.5%) is not significant. 148 10.6.2 Treatment Impact under Koyck Models For the present data, the posterior estimates were ( ω ˆ * = −0.076 , δˆ = 0.462 , φˆ = 0.009 ) under the Koyck-RP model with random parameters among matched pairs. For the sites whose treatment completion year is 2002, ORi1 was 0.916 and Equation (85a) yields OR i 2 = 0.999 . Hence, the initial reduction of 8.4% during 2003 faded away by 2004. The corresponding novelty effects for the 2003 completion year sites were also short-lived. Equation (85b) simplifies to K 2 (i, s ) = exp{−0.140 + 0.119 (1− 0.462 ) / s} , which is s depicted in Figure 10.5 along with the corresponding treatment component for the PLNI-RP model. Figure 10.5 Yearly treatment components under the PLNI-RP and Koyck-RP models (completion year 2002). 149 It is apparent from Figure 10.5 that K 2 (i, s ) has declined from 0.927 one year after the intervention to 0.89 within few years, before eventually stabilizing at 0.87. The immediate treatment impact after the intervention has thus resulted in a reduction of 7.3%. The reduction increased gradually during the following years as the treatment impact continued to build up. In the long-run, the reduction stabilizes at 12.3%, which represents the everlasting reduction due to the implementation of the countermeasures. Figure 10.5 shows also that the K 2 ( i , s ) curve lies below that for C T ( i , s ) , indicating higher direct treatment reductions under the Koyck-RP model compared to the PLNI-RP model. Further, since δˆ = 0.462 > 0 , the percent annual change varies with the number of postintervention years as follows: [ ∂ K 2 ( i , s ) / ∂s ] / K 2 ( i , s ) = d {[ 1 − ln( δ )s ] δ s − 1 } / s 2 = 0.119 [( 1 + 0.772 s ) 0.462 − 1 ] / s 2 . s Thus, while the treatment component has been declining at an increasing rate (ranging from -0.022 to 0) under the Koyck-RP model, it was decreasing at a constant rate (b = 0.0185) under the PLNI-RP model. This behavior is clearly manifested in Figure 10.5 as K 2 ( i , s ) is asymptotic (approaching a horizontal asymptote), while C T ( i , s ) is decreasing steadily. Since the estimate φˆ = 0.009 is not only fairly small but also is non-significant, the denominators of equations (85c), (85d) and (86) are nearly 1. The posterior estimates ( γˆ 1 = 0.247 , γˆ 2 = 0.418 ) under the Koyck-RP model were close to those under PLNI-RP. Thus, the ranges for K 3 ( i , s ) , K 4 ( i , s ) and K V ( i , s ) are approximately the same as those given previously for C 3 ( i , s ) , C 4 ( i , s ) and C V ( i , s ) . Under the Koyck-RP model, the estimate of the overall odds ratio was 0.877 ± 0.040 with the 95% credible interval (0.799, 0.957) indicating a significant reduction of 12.3%. 150 10.7 Summary In this chapter, four models were considered for conducting a CBA evaluation of collision counts at certain intersections improved by the installation of countermeasures under the ICBC’s RIP. The results showed that the PLNI and Koyck models allowing for random parameters among matched comparison-treatment pairs has significantly outperformed the corresponding regular models, while reducing the estimates of the extra-Poisson variance. Also, the random parameters Koyck-RP model was found to fit the current dataset considerably better than PLNIRP. The estimate of the overall odds ratio was 0.935 ± 0.042 implying a non-significant 6.5% reduction under the PLNI-RP model. In contrast, the corresponding estimate was 0.877 ± 0.040 indicating a significant 12.3% reduction under the Koyck-RP model. These odds ratios were decomposed into basic components associated with novelty effects, direct treatment effects and indirect treatment effects manifested through the traffic volumes. For the sites whose treatment was completed in 2002, the treatment component of the odds ratio has declined steadily from 1.025 to 0.935 during the after period under the PLNI-RP model. On the other hand, the corresponding treatment component has declined at an increasing annual rate under the Koyck-RP model. In fact, the immediate treatment impact resulted in a 7.3% reduction. This reduction has increased gradually in the short-term as the treatment impact continued to build up. In the long-run, the reduction stabilizes at 12.3% corresponding to the permanent treatment effect. In general, it was shown that the percent change in the treatment component is constant under PLNI, but varies with the number of post-treatment periods under the Koyck model. Thus, while the treatment component allows a limited number of treatment profiles under PLNI, a variety of profiles are possible under the Koyck model. 151 Unfortunately, the novelty effects were short-lived as the initial 8.4% reduction has evaporated within two years due to the negligible estimate of φ . Yet, the novelty effects might last longer for other datasets yielding higher estimates of the AR1 parameter. The results suggest that using a Koyck dynamic regression model and incorporating such design features as matched (yoked) comparison groups in SPFs leads to superior models. Apart from the improvement in the goodness of fit, such flexible models afford a variety of treatment profiles, besides accounting for heterogeneity due to unobserved road geometrics, traffic characteristics, environmental factors and driver behavior. 152 11 CONCLUSIONS, CONTRIBUTIONS & FUTURE RESEARCH This chapter summarizes the main conclusions and contributions of the thesis. The chapter concludes by highlighting areas for future research. 11.1 Summary and Conclusions The first main objective of this research was to address key issues related to the development of safety performance functions which have become the primary tool for evaluating road safety. The goal was is to improve the predictive power and fit of existing models since important decisions are made and considerable resources are allocated on the basis of the results of these models. The first issue investigated the inclusion of spatial effects in SPFs. Two types of spatial modeling techniques, namely, the Gaussian conditional auto-regressive (CAR) and the multiple membership (MM) models, were compared to the traditional Poisson log-normal (PLN) model. A variation of the MM model (Extended MM or EMM) was also investigated to study the effect of clustering segments within the same corridor on spatial correlation. The results provided strong evidence for incorporating spatial effects in SPFs. Apart from the improvement in goodness of fit, the inclusion of such effects explained enough variation that some of the model covariates would be rendered non-significant and thereby affecting the model’s inference. Furthermore, the results reinforce the findings of other studies in the literature which suggested that spatial models have the potential for reducing the bias, in estimating the regression coefficients, resulting from the omission of spatial variables. The second issue investigated the consequences of incorporating cluster effects in SPFs using random cluster effects and parameters. Two models that account for the cluster variation through the mean and variance of an extended Poisson log-normal (PLN) model were compared to the traditional PLN model. Random parameters PLN models were proposed to account for the cluster variation through the mean by fitting a different regression curve for each 153 cluster. Alternatively, a spatial PLN model was proposed to account for the cluster variation using an additional variance component. The results provided strong evidence for the advantages of clustering sites into rather homogeneous groups and incorporating random cluster parameters in SPFs. Apart from the improvement in the goodness of fit, such an approach can be used to gain new insights into how the covariates affect collision frequencies. As well the approach can account for heterogeneity due to unobserved road geometrics, traffic characteristics, environmental factors and driver behavior. The results have shown that the inclusion of random parameters (among clusters) in the mean function explained enough variation that some of the model covariates were rendered non-significant and thereby affecting the model inference. The third issue investigated the consequence of developing multivariate SPFs with multiple regressions to account for the nature of the correlations across collision severity levels and their influence on safety analyses. A bivariate Poisson log-normal (PLN) model, used to jointly analyze a sample of collision counts classified by two severity levels (PDO and I+F), was compared with the independent (separate) univariate PLN models. The results provided strong evidence for the benefits of developing multivariate SPFs, particularly in terms of improved precision and goodness-of-fit. The improvement in precision was due mainly to the fact that MVPLN accounts for the correlation between the latent variables underlying PDO and I+F. The results indicate that higher PDO collisions are associated with higher I+F collisions, as the collision likelihood for both types is likely to rise due to similar deficiencies in roadway design and/or other unobserved factors. This correlation was identified by other researchers in the literature. Also, the two models were compared in terms of the identification of hazardous “hotspots” locations. A new technique was developed to generalize the univariate posterior probability of excess that has been commonly proposed and applied in the literature to fit the multivariate relationship between latent variables. The results showed that almost all hazardous locations identified by the univariate models were identified by multivariate model. The results 154 also demonstrated that some hazardous locations could be overlooked if the analysis was restricted to the univariate models. The fourth issue was to investigate the difference between three approaches for analyzing collision data with potential outliers: i) outlier ignoring, ii) outlier rejecting, and iii) outlier accommodating. Two outliers accommodating approaches (a scale-mixture model and a two-component mixture/contaminated normal model) were proposed as robust alternatives to the Multivariate Poisson Lognormal (MVPLN) regression for the joint analysis of a sample of collision counts classified by two severity levels (PDO and I+F). The results provided strong evidence for the advantages of adopting an “outlier accommodating” approach when developing SPFs. The results indicated that the mixture models were successful in accommodating potential outliers while improving the model precision. This improved precision would (i) reduce the rather high uncertainty associated with the predicted values, if one adopts the outlier ignoring approach, and (ii) avoids the risk of underestimating this uncertainty, which results from using the outlier rejecting approach. A second main objective of this research was to develop SPFs for collision intervention analysis within the Controlled Before-After (CBA) framework. The outcome of these evaluations 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. The fifth issue investigated the development of a new FB approach to conduct CBA safety evaluations subject to data restrictions (aggregation). A two-way interaction model was proposed based on the Poisson-Lognormal hierarchy with covariates representing time, treatment (intervention) and interaction effects. Because the sites are clustered into homogeneous groups and the data are aggregated over these clusters, the FB two-way interaction model was extended to allow for cluster variation using either random cluster effects or parameters. The results of the analysis provided strong evidence for incorporating such 155 design features as clustering sites into homogeneous groups in CBA evaluations. In addition to the improvement in the goodness of fit, such extended models can be used to account for heterogeneity due to unobserved road geometrics, traffic characteristics, environmental factors and driver behavior. As well, the random parameters models allow the computation of odds ratios and subsequently the treatment effectiveness at both the cluster-level and the overalllevel. The sixth issue investigated the development of linear intervention models and the role that their parameters have in the computation of the odds ratios and subsequently on treatment effectiveness when conducting a FB CBA evaluation. As a result a method was developed to decompose the odds ratios for linear intervention SPFs. The decomposition of the odds ratio yielded valuable insights in the effects of the model parameters on treatment impacts. The analysis revealed that the certain parameters such as the intervention time and date interactions have a direct influence on the odds ratio and treatment effectiveness. In contrast, the parameters related to the effects of major and minor traffic volumes, represent the indirect treatment impact. In fact, the decomposition results have shown that unless the traffic volumes are subject to significant annual fluctuations, their effect on the odds ratio and treatment effectiveness would be marginal. The decomposition of the odds ratio also provided a simple and intuitive way of computing the odds ratio without resorting to additional algorithms. Given that the model specification is unchanged, then the computation of the odds ratios requires only the estimates of regression coefficients associated with the intervention time and date interactions as well as the two odds ratios for the major and minor traffic volumes. The seventh issue investigated the consequence of developing a FB CBA safety evaluation using random parameters and multivariate models. Both approaches were combined to assess the factors that influence the frequency and severity of collisions. Two models were considered for conducting the CBA evaluation. The first model uses a Multivariate Poisson Lognormal Intervention (MVPLNI) model for collision counts by severity-level with covariates 156 representing traffic volume, time, treatment (intervention) and interaction effects. The second model MVPLNI-RP extends MVPLNI by using random parameters to account for the comparison-treatment pairing. The results show that the variation among pairs has significantly improved the fit while reducing the estimates of the extra-Poisson variation. Further, the effects of the covariates on collision counts were found to vary significantly across treatmentcomparison pairs justifying the use of random parameters. The correlation between the two severity levels (I+F and PDO) of collision counts was highly significant indicating that higher PDO collisions are associated with higher I+F collisions, as the collision likelihood for both levels is likely to rise due to similar deficiencies in roadway design and/or other unobserved factors. The eighth issue investigated the validity of linear intervention models in conducting FB CBA safety evaluations. Two models were considered; the first of which is a Poisson Lognormal Intervention (PLNI) model. The second model is based on the non-linear Koyck model, which affords a rich family of SPFs that can accommodate various profiles for the time and treatment effects. The two models were extended to include random parameters to account for the variation among the comparison-treatment pairs. The results showed that the PLNI and Koyck models allowing for random parameters among matched comparison-treatment pairs has significantly outperformed the corresponding regular models, while reducing the estimates of the extra-Poisson variance. Also, the random parameters Koyck model was found to fit the current dataset considerably better than PLNI. The odds ratio was decomposed under both models into basic components associated with direct and indirect treatment effects manifested through the traffic volumes. An additional advantage of using the Koyck model is its ability to isolate another component related to the novelty effects of the intervention. This is an important phenomenon that captures how effective the countermeasure was at reducing collisions at the start of the implementation and for how long has the effectiveness (due to novelty of the countermeasure) lasted. In general, it was shown that the percent change in the treatment component is constant 157 under the linear PLNI model, but varies with the number of post-treatment periods under the non-linear Koyck model. Thus, while the treatment component allows a limited number of treatment profiles under PLNI, a variety of profiles are possible under the Koyck model. The results suggest that using a Koyck dynamic regression model and incorporating such design features as matched (yoked) comparison groups in SPFs leads to superior models. Apart from the improvement in the goodness of fit, such flexible models afford a variety of treatment profiles, besides accounting for heterogeneity due to unobserved road geometrics, traffic characteristics, environmental factors and driver behavior. 11.2 Research Contributions The following are the main contributions of this research: 1. The development of SPFs that account for spatial correlation in collision data using random effects models; 2. The development of SPFs that include random cluster effects and parameters in the variance and mean functions, respectively; 3. The development of multivariate SPFs that account for the correlations between different collision types, severity, etc., that commonly exist in collision data; 4. The development of a new hotspot multivariate identification technique that generalizes the univariate posterior probability of excess (that has been commonly proposed and applied in the literature) to fit the multivariate relationship between the latent variables underlying the various severity levels. 5. The development of techniques for the multivariate identification of potential outliers in collision data; 6. The development of multivariate mixture SPFs to accommodate potential outliers in collision data; 158 7. The development of fully Bayesian tools for investigating the safety change in treated sites due to the implementation of safety countermeasures. 8. The development of a fully Bayesian CBA evaluation approach when the data are subject to constraints; 9. Demonstrating the applicability of using linear intervention models with random parameters among matched comparison-treatment pairs to assess the effectiveness of safety countermeasures; 10. Demonstrating the applicability of using multivariate linear intervention models to account for the correlations between different collision types, severity levels, etc., while assessing the effectiveness of safety countermeasures; 11. Demonstrating the applicability of using random parameters and multivariate linear intervention models to assess the effectiveness of safety countermeasures; 12. The development of a novel non-linear technique to conduct CBA safety evaluations that can afford a variety of treatment profiles; 13. The development of a novel methodology to decompose the OR into basic components related to direct and indirect treatment effects; 14. The development of a simplified approach to calculate the OR by providing straightforward equations in terms of model parameters for the computation of the OR and its associated components without resorting to additional algorithms. 11.3 Future Research 11.3.1 Network Safety Performance Functions The development of SPFs has been restricted to modeling intersections or road segments separately. With the exception of few studies (Abdel-Aty and Wang, 2006; Wang et al., 2006; Das et al., 2008), this practice ignores the relationship between the safety levels at intersections and neighboring road segments or vice versa. The spatial models developed in this thesis could 159 be used to assess the safety of a network composed of both intersections and road segments. Preliminary analysis on this topic has provided some valuable insight on the relationship between collisions on intersections and road segments in a network. Two spatial SPFs were developed and linked to model collisions on intersection and road segments for the City of Richmond, British Columbia. Due to data limitations only traffic volumes were used as covariate for the intersection models whereas traffic volumes, segment lengths, and unsignalized intersection densities were used to develop SPFs for road segments. The network consisted of 72 road segments and 43 intersections, the network boundaries are shown in Figure 11.1. The network models were compared to spatial models and the independent models in terms of goodness-of-fit, inference and identification of hotspots. The results were very promising since the network models provided a better fit than the spatial models. In terms of hotspot identification, the independent and spatial models identified locations that were somewhat dispersed in the network. However, the identification based on the network-based models was more consistent by flagging three road segments and their in-between intersections making the identification more concentrated. Moreover, the road segments and intersections that were identified using the network models were on one of the major arterials in the City of Richmond. Unfortunately, some of the network model parameters were not significant. This represented a major drawback since the models were already developed using the minimum number of explanatory variables and therefore they were subject to omitted variable bias. It is obvious that for these network models to be applied a more comprehensive data set is needed and more research is required to determine the feasibility and application of network-based models. 160 Figure 11.1 Richmond network boundaries Courtesy of googlemaps.com 11.3.2 Implications for Collision Modification Factor Estimations and Benefit-Cost Analysis The second part of the thesis focused on developing a novel CBA evaluation methodology in a full Bayes setting. This research component could be extended to conduct an economic analysis and compute different monetary measures (i.e., internal rate of return, benefit-cost ratio or net present value) to estimate the cost effectiveness of different safety treatments or countermeasures. Economic analysis of countermeasures is important from a retrospective viewpoint. For example, after the implementation of any safety-based countermeasures road safety agencies need to ensure that they have achieved a proper return on their investment. Given the advantages offered by using a fully Bayes CBA framework, the approaches used in part two of the thesis could be extended to compute the several economic measures while 161 accounting for all sources of uncertainty including those related to discount rate, treatment service life, costs, etc. Moreover, the current practice associates the improvements of a particular safety countermeasure to a single Collision Modification Factor (CMF) with its standard error. However, the research on the linear (PLNI) and non-linear (Koyck) intervention models has shown that treatment impact varies with time. The implications of having a variable CMF on estimating the economic benefits of safety countermeasures needs to be investigated. 11.3.3 Improvements to Safety Performance Functions Driver Behavior Measures: A safety performance function is typically used to relate the occurrence of collisions to the different traffic and geometric variables at an intersection or road segment. It will be of interest to investigate the inclusion of human factors in the development of SPFs. Traditionally, it is assumed that driver behavior is indirectly accounted for by including geometric variables which in turn has an impact on collision occurrence. Dietze et al. (2005), for example, considered developing SPFs that included infrastructure impacts and driver behavior by grouping road entities that had similar geometric properties under the presumption that similar geometric properties cause similar driving behavior and therefore the effects on road safety have to be similar. The research by Dietze et al. (2005) provides a means to account for human factors but additional research is warranted on alternative approaches to account for driver behavior. For example, additional parameters can be derived from the knowledge of human perception, information processing and decision making and integrated in the development of SPFs. The topic is worth investigating to enhance the predictability of SPFs and provide insight into the driver behavior factors that influence collision occurrence. Exposure Measure: Traffic volume (or flow counts) is perhaps the most common predictor used in SPFs for both intersection and road segments. Ideally, traffic volume should be measured (e.g., counted) for the entire period over which collisions are aggregated. Yet, typically, a single- or 2-day sample count is obtained and the average of that sample is then 162 used as an estimate of the traffic volume. Improved estimates of traffic volume can be obtained from repeated observations of traffic volume or from validation data. However, both alternatives are expensive since carrying out the volume counts represents one of the major costs of any such safety study. Moreover, traffic volume (or flow counts) is an example of predictors that are measured with uncertainty, which if not resolved can bias the regression estimates (Maher and Summersgill, 1996; Davis, 2000). Given the modeling advantages of the full Bayes approach, the degree of that bias has to be assessed. The underlying idea would be to explore the effects of adjusting the volume data for temporal correlation and measurement errors in order to improve the estimates of traffic flow. Surrogate Safety Measures: The majority of road safety studies use collision records as the primary source of data. However, the reliability and usefulness of the collision data are often considered suspect. If collision data are faulty, caused either by a reduction in reporting levels or from the deterioration in quality, then the ability of road safety professionals to engineer solutions to address problems may be severely jeopardized (de Leur and Sayed, 2001). Due to the lack or poor quality of collision data several other surrogate measures were proposed as surrogate safety measures. Recently, the application of Traffic Conflict Technique (TCT) in assessing the safety of a road entity has been gaining attention among safety researchers and practitioners. Several studies have demonstrated the applicability of extracting conflict data by using i) field observers (William, 1972; Zegeer and Dean, 1978; Crowe, 1990; Sayed and Zein, 1999), ii) simulation models (Sayed et al. 1994; Persaud and Mucsi 1995; Huang and Pant 1999; Archer, 2001), and iii) video-camera (Ismail et al., 2009a,b; 2010a,b) to assess the safety of a particular road entity. Regardless of the method by which traffic conflicts are recorded, the approach offers traffic safety analysts a unique opportunity to observe a large number of unsafe vehicle interactions (potential collisions) in order to pinpoint the failure mechanism that leads to such an unsafe behavior. It has been argued in the literature that traditional collision-based analysis fails in enabling the identification of cause and effect relationships with regard to 163 collision probabilities (Lord and Mannering, 2010). Alternatively, the TCT offers a unique way to gain a better understanding of the factors that affect the probability of collisions thereby being able to better predict the likelihood of collisions, while providing directions for policies and countermeasures that target the reduction of collisions. As a result, the relationship between collisions and conflicts must first be established to use traffic conflicts as surrogates to collisions for safety analysis. Although traffic conflicts are more frequent than road collisions and are of marginal social cost; they provide more insight into the failure mechanism that leads to road collisions. Using some of the techniques developed in this thesis, future research can focus on (i) establishing a relationship between collisions and conflicts and (ii) developing a statistical methodology to model conflicts and relate it to the different traffic, road geometry, and driver behavior characteristics. 164 REFERENCES Abdel-Aty, M., and Wang, X., 2006. Crash Estimation at Signalized Intersections Along Corridors: Analyzing Spatial Effect and Identifying Significant Factors. Transportation Research Record: Journal of the Transportation Research Board. 1953, 98-111. Abdel-Aty, M., Pemmanaboina, R., and Hsia, L., 2006. Assessing Crash Occurrence on Urban Freeways by Applying a System of Interrelated Equations. Transportation Research Record: Journal of the Transportation Research Board. 1953, 1-9. Abdel-Aty, M., Devarasetty, P.C., and Pande, A., 2009. Safety evaluation of multilane arterials in Florida. Accident Analysis & Prevention, 41 (4), 777-788. Aguero-Valverde, J. and Jovanis, P.P., 2006. Spatial analysis of fatal and injury crashes in Pennsylvania. Accident Analysis and Prevention. 38, 618–625. Aguero-Valverde, J. and Jovanis, P.P., 2008. Traffic Analysis of Road Crash Frequency Using Spatial Models. 87th Annual Meeting of the Transportation Research Board, Washington, D.C. Aguero-Valverde, J. and Jovanis, P.P., 2009. Bayesian multivariate Poisson log-normal models for crash severity modeling and site ranking. Transportation Research Record. 2136, 82-91. Albert, J.H., 1999. Criticism of a hierarchical model using Bayes factors. Statistics in Medicine. 18, 287-305. American Academy of Pediatrics; Committee on Injury and Poison Prevention, 2002. Selecting and using the most appropriate car safety seats for growing children: Guidelines for counseling parents. Pediatrics. 109, 550-3. Anastasopoulos, P.Ch. and Mannering, F., 2009. A note on modeling vehicle accident frequencies with random-parameters count models. Accident Analysis and Prevention. 41 (1), 153–159. 165 Archer, J., 2001. Traffic Conflict Technique Historical to current State-of-the-Art. Effektmodeller för vägtrafikanläggningar (EMV), TRITA-INFRA 02-010, ISSN 1651-0216, ISRN KTH/INFRA--02/010-SE, Stockholm, September. Aul, N. and 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. Bedrick, E.J., Christensen, R. and Johnson, W., 1996. A new perspective on priors for generalized linear models. Journal of the American Statistical Association. 91, 14501460. Berkane, M. and Bentler, P.M., 1988. Estimation of contamination parameters and identification of outliers in multivariate data. Sociological Methods and Research. 17, 55-64. Bernardinelli, L., Clayton D. and Montomoli, C., 1995. Bayesian estimates of disease maps: How important are priors. Statistics in Medicine. 14, 2411-2431. Besag, J., Green, P., Higdon, D. and Mengersen, K., 1995. Bayesian Computation and Stochastic Systems (with Discussion). Statistical Science. 10, 3-66. Besag, J., York, J. and Mollié, A., 1991. Bayesian Image Restoration with Two Applications in Spatial Statistics. Annals of the Institute of Statistical Mathematics. 43, 1-75. Black, W.R. and Thomas, I., 1998. Accidents on Belgium’s motorways: a network autocorrelation analysis. Journal of Transport Geography. 6 (1), 23-31. Breslow, N.E. and Clayton, D.G., 1993. Approximate Inference in Generalized Linear Mixed Models. Journal of the American Statistical Association. 88 (421), 9-25. Brijs, T., Karlis, D., Van den Bossche, F. and Wets, G., 2007. A Bayesian model for ranking hazardous road sites. Journal of the Royal Statistical Society: Series A (Statistics in Society). 170 (4), 1001-1017. Brooks, S.P. and Gelman, A., 1998. Alternative methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics. 7, 434-455. 166 Cambridge Systematics Inc., 2008. Crashes vs. Congestion: What’s the Cost to Society? Presented to the Automobile Association of America, March. Cameron, A.C. and Trivedi, P.K., 1998. Regression Analysis of Count Data. Econometric Society Monographs 30, New York: Cambridge University Press. Carlin, B.P. and Louis, T.A., 1996. Bayes and Empirical Bayes Methods for Data Analysis. London: Chapman and Hall/CRC, 15. Carriquiry, A. and Pawlovich, M.D., 2005. From empirical Bayes to full Bayes: methods for analyzing traffic safety data. http://www.dot.state.ia.us/crashanalysis/pdfs/ eb_fb_comparison_whitepaper_october2004.pdf. Carroll, R., Rupppert, D. and Sefanski, L., 1995. Measurement Error in Nonlinear Models, Chapman-Hall. Cheng, W. and Washington, S.P., 2005. Experimental evaluation of hotspot identification methods. Accident Analysis and Prevention. 37, 870-881. Chesher, A., 1991. The effect of measurement error. Biometrika 78, 451–462. Chib, S. and Winkelmann, R., 2001. Markov chain Monte Carlo analysis of correlated count data. Journal of Business and Economic Statistics. 19, 428–435. Clayton, D. and Kaldor, J., 1987. Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics. 43 (3), 671-681. Congdon, P., 2006. Bayesian Statistical Modeling. John Wiley & Sons, New York, 2nd edition. Crowe, C.E., 1990. Traffic Conflict Values for three-leg, un-signalized intersections. Transportation Research Record. 1287, 185-194. Czado, C. and 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. and 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. 167 Das, A., Pande, A., Abdel-Aty, M., and Santos, J.B., 2008. Characteristics of Urban Arterial Crashes Relative to Proximity to Intersections and Injury Severity. Transportation Research Record: Journal of the Transportation Research Board. 2083, 137-144. Davis, G.A., 2000. Estimating traffic accident rates while accounting for traffic-volume estimation error: A Gibbs sampling approach. Transportation Research Record. 1717, 94–101. de Leur, P. and Sayed, T., 2001. Using claims prediction model for road safety evaluation. Canadian Journal Civil Engineering. 28(5), 804–812. de Leur, P., 2001. Improved approaches to manage road safety infrastructure. Ph.D. Dissertation. Department of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada. Dietze, M., Ebersbach, D., Lippold, C., Mallschutzke, K., Gatti, G. and Wieczynski, A., 2005. Safety Performance Function. RIPCORD-ISEREST. Final Report, 506184. El-Basyouny, K. and Sayed, T. 2006. Comparison of Two Negative Binomial Regression Techniques in Developing Accident Prediction Models. Transportation Research Record. 1950, 9-16. El-Basyouny, K. and Sayed, T., 2010a. Application of generalized link functions in developing accident prediction models. Safety Science. 48, 410–416 El-Basyouny, K. and Sayed, T., 2010b. Safety Performance Functions with measurement errors in traffic volume. Safety Science. 48, 1339–1344 Francis, B., Green, M. and Payne, C., 1993. The GLIM System: Release 4 Manual. Clarendon Press, Oxford. Gelman, A., Carlin, J.B., Stern, H.S. and Rubin, D.B., 2004. Bayesian Data Analysis. 2nd Edition. Chapman & Hall, London. Gelman, A., Meng, X. and Stern, H., 1996. Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica. 6, 733-807. 168 Geweke, J., 1993. Bayesian treatment of the independent Student-t linear model. Journal of Applied Econometrics. 8, S19-S40. Gilks, W.R., Richardson, S., and Spiegelhalter, D.J., 1996. Markov Chain Monte Carlo in Practice. London: Chapman and Hall. Goldstein, H., 1995. Multilevel Statistical Models. London: Arnold. Goldstein, H., Rasbash, J., Plewis, I., Draper, D., Browne, W., Yang, M., Woodhouse, G. and Healy, H., 1998. A User’s Guide to MLwiN. London: Institute of Education. Guo, F., Wang, X., and Abdel-Aty, M., 2010. Modeling signalized intersection safety with corridor-level spatial correlations. Accident Analysis & Prevention, 42(1), 84-92. Haleem, K. and Abdel-Aty, M. 2009. Examining traffic crash injury severity at unsignalized intersections. Journal of Safety Research, 41(4), 347-357. Hauer, E. and Bamfo, J., 1997. Two tools for finding what function links the dependent variable to the explanatory variables. Proceedings of the ICTCT 1997 Conference. Lund, Sweden. Hauer, E., 1995. On Exposure and Accident Rate. Traffic Engineering and Control. 36 (3), 134138 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., Kononov, J., Allery, B. and Griffith, M.S., 2002. Screening the Road Network for Sites with Promise. Transportation Research Record. 1784, 27–32. Hauer, E., Ng, J.C.N. and Lovell, J., 1988. Estimation of Safety at Signalized Intersections. Transportation Research Record. 1185, 48-61. Heydecker, B.G. and 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. 169 Higle, J. and Witkowski, J.M. 1988. Bayesian identification of hazardous sites. Transportation Research Record. 1185, 24-35. Howard, A.W., 2002. Automobile restraints for children: A review for clinicians. CMAJ. 167, 76973. Huang, H. and Abdel-Aty, M., 2010. Multilevel data and Bayesian analysis in traffic safety. Accident Analysis & Prevention, 42(6), 1556-1565. Huang, P. and Pant, P., 1996. Simulation neural-network model for evaluating dilemma zone problems at high-speed signalized intersections. Transportation Research Record. 1456. Hughes, W., Eccles, K., Harwood, D., Potts, I. and 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, D.C. Available from http://onlinepubs.trb.org/onlinepubs/nchrp/nchrp_w62.pdf Ismail, K., Sayed, T. and Saunier, N., 2009b. Automated pedestrian safety analysis using video data in the context of scramble phase intersections. In Annual Conference of the Transportation Association of Canada, Vancouver, B.C. Ismail, K., Sayed, T. and Saunier, N., 2010a. Automated safety analysis using video sensors: Technology and case studies. In Canadian Multidisciplinary Road Safety Conference, Niagara Falls, Ontario. Ismail, K., Sayed, T. and Saunier, N., 2010b. Automated Analysis of Pedestrian-vehicle Conflicts: A Context For Before-and-after Studies. In Transportation Research Board Annual Meeting Compendium of Papers, Washington D.C. 10-3739. Accepted for publication in Transportation Research Record: Journal of the Transportation Research Board. Ismail, K., Sayed, T., Saunier, N. and Lim, C., 2009a. Automated Analysis of Pedestrian-Vehicle Conflicts using Video Data. Transportation Research Record. 2140, 44-54. 170 Jacobs, G., Aeron-Thomas, A. and Astrop, A., 2000. Estimating global road fatalities. Crowthorne, Transport Research Laboratory, (TRL Report, No. 445). Jovanis, P.P. and Chang, H.L., 1986. Modeling the Relationship of Accidents to Miles Traveled. Transportation Research Record. 1068, 42-51. Judge, G.G., Hill, R.C., Griffiths, W.E. Lutkepohl, H. and Lee, T.C., 1988. Introduction to the Theory and Practice of Econometrics, 2nd Edition. Wiley: New York. Karlis, D. and Meligkotsidou, L., 2005. Multivariate Poisson regression with covariance structure. Statistics and Computing. 15, 255–265. Karlis, D., 2003. An EM algorithm for multivariate Poisson distribution and related models. Journal of Applied Statistics. 30, 63-77. Kim, H., Sun, D. and Tsutakawa, R.K., 2002. Lognormal vs. gamma: extra variations. Biometrical Journal. 44 (3), 305–323. Koyck, L.M., 1954. Distributed Lags and Investment Analysis. North-Holland: Amesterdam. Kulmala, R., 1995. Safety at Rural Three-and Four-arm Junctions. Development of Accident Prediction Models. Technical Research Centre of Finland, VTT 233, Espoo. Ladron de Guevara, F. and Washington, S., 2004. Forecasting crashes at the planning level. A simultaneous negative binomial crash model applied in Tucson, Arizona. Transportation Research Record. 1897, 191–199. Lan, B., Persaud, B., Lyon, C. and 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. Lane, P.W., Galwey, N.W. and Alvey, N.G., 1988. GENSTAT 5: an Introduction. Oxford University Press, Oxford. Lange, K.L., Little, R.J.A. and Taylor, J.M.G., 1989. Robust statistical modeling using the t distribution. Journal of the American Statistical Society. 84, 881-896. 171 Langford, I.H., Leyland, A.H., Rasbash, J. and Goldstein, H., 1999. Multilevel modeling of the geographical distributions of diseases. Applied Statistics. 48 (2), 253-268. Levine, N., Kim, K.E. and Nitz L.H., 1995b. Spatial analysis of Honolulu motor vehicle crashes. II. Zonal generators. Accident Analysis and Prevention. 27 (5), 675–685. Levine, N., Kim, K.E. and Nitz, L.H., 1995a. Spatial analysis of Honolulu motor vehicle crashes. I. Spatial patterns. Accident Analysis and Prevention. 27 (5), 663–674. Li, C.C., Lu, J.C., Park, J., Kim, K., Brinkley, P.A. and Peterson, J.P., 1999. Multivariate zeroinflated Poisson models and their applications. Technometrics. 41 (1), 29–38. Li, W., Carriquiry, A., Pawlovich, M. and Welch, T., 2008. The choice of statistical models in road safety countermeasures effectiveness studies in Iowa. Accident Analysis and Prevention. 40 (4), 1531–1542. 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. and Mannering, F., 2010. The statistical analysis of crash-frequency data: A review and assessment of methodological alternatives. Transportation Research Part A. 44, 291– 305. Lord, D. and Miranda-Moreno, L.F., 2008. Effects of low sample mean values and small sample sizes on the parameter estimation of hierarchical Poisson models for motor vehicle crashes: a Bayesian perspective. Safety Science. 46, 751-770. Lord, D. and 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. and Persaud, B., 2000. Accident prediction models with and without trend: application of the generalized estimating equation. Transportation Research Record, 1717. 02108. 172 Lord, D., 2000. The prediction of accidents on digital networks: characteristics and issues related to the application of accident prediction models. Ph.D. Dissertation. Department of Civil Engineering, University of Toronto, Toronto, Ontario, Canada. 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., and Geedipally, S.R., 2008. Application of the Conway-MaxwellPoisson generalized linear model for analyzing motor vehicle crashes. Accident Analysis and Prevention. 40, 1123-1134. Lord, D., Washington, S.P. and 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, 35-46. Lord, D., Washington, S.P. and Ivan, J.N., 2007. Further notes on the application of zero inflated models in highway safety. Accident Analysis and Prevention. 39 (1), 53–57. Lovegrove, G.R., 2006. Community-Based, Macro-level Collision Prediction Models. Ph.D. Dissertation. Department of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada. Lunn, D.J., Thomas, A., Best, N. and Spiegelhalter, D., 2000. WinBUGS - a Bayesian modeling framework: concepts, structure, and extensibility. Statistics and Computing, 10, 325337. Ma, J. and Kockelman, K.M., 2006. Bayesian multivariate Poisson regression for models of injury count by severity. Transportation Research Record. 1950, 24–34. Ma, J., Kockelman, K.M. and 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. 173 MacNab, Y.C., 2004. Bayesian spatial and ecological models for small-area accident and injury analysis. Accident Analysis and Prevention. 36 (6), 1028–1091. Maher, M.J. and Summersgill, I., 1996. A comprehensive methodology for the fitting of predictive accident models. Accident Analysis and Prevention. 28 (3), 281-296. Malyshkina, N. and Mannering, F., 2010. Zero-state Markov switching count-data models: An empirical assessment. Accident Analysis and Prevention. 42(1), 122-130. Malyshkina, N.V., Mannering, F.L. and 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. and Nelder, J.A., 1989. Generalized Linear Models, 2nd Edition. Chapman & Hall, London. Miaou, S. and Lum, H., 1993. Modeling Vehicle Accident and Highway Geometric Design Relationships. Accident Analysis and Prevention. 25(6), 689-709. Miaou, S. and Song, J., 2005. Bayesian ranking of sites for engineering safety improvements: Decision parameter, treatability concept, statistical criterion, and spatial dependence. Accident Analysis and Prevention. 37 (4), 699–720. Miaou, S., 1996. Measuring the goodness-of-fit of accident prediction models. Federal Highway Administration, Publication No. FHWA-RD-96-040, Virginia, U.S. Miaou, S., Song, J.J. and Mallick, B.K., 2003. Roadway traffic crash mapping: a space–time modeling approach. Journal of Transportation and Statistics. 6(1), 33–57. Miaou, S.P. and 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. 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. 174 Miranda-Moreno, L.F. and Lord, D., 2007. Evaluation of alternative hyper-priors for Bayesian road safety analysis. Presented at the 87th Annual Meeting of the Transportation Research Board, Washington, D.C. Miranda-Moreno, L.F., 2006. Statistical models and methods for identifying hazardous locations for safety improvements. Ph.D. Dissertation. Department of Civil Engineering, University of Waterloo, Waterloo, Ontario, Canada. Miranda-Moreno, L.F., Fu, L., Saccomanno, F.F. and Labbe, A., 2005. Alternative Risk Models for Ranking Locations for Safety Improvement. Transportation Research Record. 1908, 1-8. Mitra, S. and Washington, S., 2007. On the nature of over-dispersion in motor vehicle crash prediction models. Accident Analysis and Prevention. 39, 459-468. Mountain, L., Maher, M. and Fawaz, B. 1998. The influence of trend on estimates of accident at junctions. Accident Analysis and Prevention. 30 (5). National Cancer Institute of Canada. 2001. Canadian Cancer Statistics. Ottawa: National Institute of Cancer. Nicholson, A. 1999. Analysis of spatial distributions of accidents, Safety Science. 31, 71-91. Nicholson, A. and Turner, S., 1996. Estimating accidents in a road network. Proceedings of Roads 96 Conference, Part 5, New Zealand, 53-66. Ozbay, K., Yanmaz-Tuzel, O., Ukkusuri, S.V., and Bartin, B. (2009) Safety Comparison of Roadway Design Elements on Urban Collectors with Access, Report Number: FHWANJ-2009-008 2009. Pankratz, A., 1991. Forecasting with Dynamic Regression Models. Wiley: New York. Park, B.J. and Lord, D., 2009. Application of finite mixture models for vehicle crash data analysis. Accident Analysis and Prevention. 41, 683-691. 175 Park, E.S. and Lord, D., 2007. Multivariate Poisson–lognormal models for jointly modeling crash frequency by severity. Transportation Research Record. 2019, 1–6. Park, E.S., Park, J. and Lomax, T.J., 2009 A fully Bayesian multivariate approach to before-after safety evaluation. Presented at the 89th Annual Meeting of the Transportation Research Board, Washington, D.C. Pawlovich, M.D., Li, W., Carriquiry, A. and 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. and 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. and Mucsi, K., 1995. Microscopic accident potential models for two-lane rural roads. Transportation Research Record. 1485, 134-139. Persaud, B., Lan, B., Lyon, C. and Bhim, R., 2009. Comparison of empirical Bayes and full Bayes approaches for before-after road safety evaluations. Presented at the 89th Annual Meeting of the Transportation Research Board, Washington, D.C. Persaud, B.N., 1986. Safety migration, the influence of traffic volumes, and other issues in evaluating safety effectiveness: some findings on conversion of intersections to multiway stop control. Transportation Research Record. 1068, 108–114. Qin, X., Ivan, J.N. and Ravishankar, N., 2004. Selecting exposure measures in crash rate prediction for two-lane highway segments. Accident Analysis and Prevention. 36 (2), 183–191. Qin, X., Ivan, J.N., Ravishanker, N. and Liu, J., 2005. Hierarchical Bayesian estimation of safety performance functions for two-lane highways using Markov Chain Monte Carlo modeling. Journal of Transportation Engineering. 131(5), 345-351. 176 Relles, D.A., and Rogers, W.H., 1977. Statisticians are fairly robust estimators of location. Journal of the American Statistical Society. 72, 107-111. Robert, C.P. and Casella, G., 1999. Monte Carlo Statistical Methods. New York: SpringerVerlag. Roberts, G.O. and Rosenthal, J.S., 1998. Markov chain Monte Carlo: Some Practical Implications of Theoretical Results. Canadian Journal of Statistics. 26, 5-31. SAS Institute Inc., Version 9.1, PROC MIXED, PROC GENMOD, Cary, USA, 2002-2003. Sawalha, Z. and Sayed, T., 2001. Evaluating safety of urban arterial roadways. Journal of Transportation Engineering. 127 (2), 151-158. Sawalha, Z. and Sayed, T., 2006. Traffic accidents modeling: some statistical issues. Canadian Journal of Civil Engineering. 33 (9), 1115-1124. Sayed, T. and Abdelwahab, W., 1997. Using accident correctability to identify accident prone locations. Journal of Transportation Engineering, ASCE. 123 (2), 107-113. Sayed, T. and de Leur, P., 2007. Evaluation of Edmonton’s intersections safety camera programs. Transportation Research Record. 2009, 37-45. Sayed, T. and Zein, S., 1999. Traffic Conflict Standards for Intersections. Transportation Planning and Technology. 22, 309-323. Sayed, T. et al., 1994. Simulation of Traffic Conflicts at Unsignalized Intersections with TSCSim. Accident Analysis and Prevention. 26(5), 593-607. Sayed, T., Abdelwahab, W. and Navin, F., 1995. Identifying Accident-Prone Locations Using Fuzzy Pattern Recognition. Journal of Transportation Engineering. 121 (4), 352-358. Sayed, T., deLeur, P. and Sawalha, Z., 2004. Evaluating the Insurance Corporation of British Columbia Road Safety Improvement Program. Transportation Research Record. 1865, 57-63. Sayed, T., El-Basyouny, K. and Pump, J., 2006. Safety evaluation of stop sign in-fill program. Transportation Research Record. 1953, 201-210. 177 Sayed, T., Saunier, N., Lovegrove, G., and de Leur, P., 2010. Advances in Proactive Road Safety Planning. Proceedings of the 20th Canadian Multidisciplinary Road Safety Conference. Niagara Falls, Ontario, June 6-9. Schluter, P.J., Deely, J.J. and Nicholson, A.J., 1997. Ranking and selecting motor vehicle accident sites by using a hierarchical Bayesian model. The Statistician. 46 (3), 293– 316. Shankar V.N., Albin, R.B., Milton, J.C. and Mannering, F.L., 1998. Evaluating median cross-over likelihoods with clustered accident counts: An empirical inquiry using the random effects negative binomial model. Transportation Research Record. 1635, 44-48. Shankar, V., Milton, J., and Mannering, F.L., 1997. Modeling accident frequency as zero-altered probability processes: an empirical inquiry. Accident Analysis and Prevention. 29 (6), 829–837. Shively, T.S., Kockelman, K. and Damien. P., 2010. A Bayesian semi-parametric model to estimate relationships between crash counts and roadway characteristics. Transportation Research Part B: Methodological. 44(5), 699-715. Spiegelhalter, D., Thomas, A., Best, N. and Lunn, D., 2005. WinBUGS User Manual. MRC Biostatistics Unit, Cambridge. Available from http://www.mrc-cam.ac.uk/bugs. Spiegelhalter, D.J., Best, N.G., Carlin, B.P. and Van der Linde, A., 2002. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society B. 64, 1–34. Spiegelhalter, D.J., Thomas, A. and Best, N.G., 1996. Computation on Bayesian graphical models. In J.M. Bernardo, J.O. Berger, A.P. Dawid, and A.F.M. Smith, editors, Bayesian Statistics, 5, 407–425, Oxford University Press, Oxford. Stern, H.S. and Cressie, N., 2000. Posterior predictive model checks for disease mapping models. Statistics in Medicine. 19, 2377-2397. Tanner, J.C., 1953. Accidents at rural three-way junctions. J. Inst. Highway Eng. 2(11), 56–67. Traffic Collision Statistics, 2006. Annual Report on British Columbia traffic collision. 178 Tsionas, E.G., 2001. Bayesian multivariate Poisson regression. Communications in StatisticsTheory and Methods. 30, 243–255. Tunaru, R., 2002. Hierarchical Bayesian models for multiple count data. Austrian Journal of Statistics. 31 (3), 221-229. U.S. Department of Health and Human Services, Centers for Disease Control, National Center for Health Statistics, National Vital Statistics Reports: Deaths, 1991–2000 issues. Ulfarsson, G.F. and Shankar, V.N., 2003. Accident Count Model Based on Multiyear CrossSectional Roadway Data with Serial Correlation. Transportation Research Record. 1840, 193-197. Vodden, K., Smith, D., Eaton, F. and Mayhew, D., 2007. Analysis and estimation of the social cost of motor vehicle collisions in Ontario: final report. Presented to Ministry of Transportation, August. Wang, X. and Abdel-Aty, M., 2006. Temporal and spatial analysis of rear-end crashes at signalized intersections. Accident Analysis and Prevention. 38 (6), 1137-1150. Wang, X., Abdel-Aty, M., and Brady, P.A., 2006. Crash Estimation at Signalized Intersections: Significant Factors and Temporal Effect. Transportation Research Record: Journal of the Transportation Research Board. 1953, 10-20. William, T.M., 1972. An Evaluation of Traffic Conflict Technique. Hwy. Res. Rec. 384, 1-8. Winkelmann, R., 2003. Econometric Analysis of Count Data. Springer, Germany. Wood, G.R., 2002. Generalized linear accident models and goodness of fit testing. Accident Analysis and Prevention. 34(4), 417–427. World Health Organization, 2004. World report on road traffic injury prevention. Geneva. Yang, M., Rasbash, J., Goldstein, H. and Barbosa, M., 1999. MLwiN Macros for Advanced Multilevel Modeling. Version 2.0. Multilevel Models Project, Institute of Education, University of London, (http://www.ioe.ac.uk/multilevel/), December. 179 Ye, X., Pendyala, R.M., Washington, S.P., Konduri, K. and Oh, J., 2009. A simultaneous equations model of crash frequency by collision type for rural intersections. Safety Science. 47 (3). 443–452. Zegeer, C.V. and Deen, R.C., 1978. Traffic conflict as a diagnostic tool in highway safety. Transportation Research Record. 667, 48–55. 180 APPENDICES Appendix A: Mathematical Justification for Model Comparison Let f ( y | θ ) and f ( y ) denote the conditional and marginal distributions of y, where θ denote the vector of parameters associated with y. Note that the posterior distribution is given by f (θ | y ) = f ( y | θ ) f (θ ) / f ( y |) , where f (θ ) is the prior for ࢲ. The Bayesian deviance is given by (Speiglharter et al., 2002) D( θ ) = −2 ln{ f ( y | θ )} + 2 ln{ f ( y )} . (49) ∞ Let θ = E[θ | y ] = ∫0 θf (θ | y )dθ posterior mean of ࢲ D = E[ D(θ ) | y ] = ∫0∞ D(θ ) f (θ | y)dθ posterior mean of D(ࢲ) Define p = D − D( θ ) where D( ) is obtained by substituting in (49) Then DIC = D + p For K collision categories, let y and θ be partitioned as ( y1 ,..., y K ) and ( θ 1 ,...,θ K ) . Define DIC k = D k + p k , pk = D k − Dk ( θ k ) , D k = E [ Dk ( θ k ) | y k ] , θ k = E [ θ k | yk ] and Dk ( θ k ) = −2 ln{ f ( y k | θ k )} + 2 ln{ f ( y k )} . Under independent models and priors, we have that f ( y | θ ) = ∏ kK=1 f ( y k | θ k ) and f ( y ) = ∏kK=1 f ( y k ) . Hence, we have D (θ ) = −2 ln{∏ kK=1 f ( y k | θ k )} + 2 ln{∏kK=1 f ( y k )} = ∑ kK=1 {− 2 ln f ( y k | θ k ) + 2 ln f ( y k )} = ∑ kK=1 D k (θ k ) D = ∑kK=1 ∫0 {− 2 ln f ( y k | θ k + 2 ln f ( y k )} f (θ | y)dθ = ∑kK=1 ∫0 Dk (θ k ) f (θ | y )dθ ∞ ∞ This implies that 181 D = ∑ kK=1 E ( D k (θ k ) | y k ) = ∑ kK=1 D k Also D (θ ) = −2 ln{ f ( y | θ )} + 2 ln{ f ( y )} = −2 ln ∏ kK=1 f ( y k | θ k ) + 2 ln ∏kK=1 f ( y k ) This implies that D (θ ) = ∑ kK=1{−2 ln f ( y k | θ k ) + 2 ln f ( y k )} = ∑ kK=1 D k (θ k ) Thus p = D − D(θ ) = ∑kK=1 D k − ∑kK=1 Dk (θ k ) = ∑kK=1[ D k − Dk (θ k )] = ∑kK=1 p k And DIC = D + p = ∑kK=1 D k + ∑kK=1 p k = ∑kK=1 ( D k + p k ) = ∑kK=1 DIC k Thus, it can be seen that the multiplicative conditional and marginal distributions of y translate additively in the Bayesian deviance (49) leading to DIC = ∑ kK=1 DIC k . (50) 182 Appendix B: Odds Ratio Decomposition under Linear Intervention Model 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. Recall Equation (25) and note that B CB CB ln( µ CB it ) = α 0 + α 2 t + γ 1 ln( V 1it ) + γ 2 ln( V 2 it ) , A A CA CA ln( µ CA it ) = α 0 + α 2 t + α 3 ( t − t 0 i ) + γ 1 ln( V 1it ) + γ 2 ln( V 2 it ) , B TB TB ln( µ TB it ) = α 0 + α 1 + ( α 2 + α 4 ) t + γ 1 ln( V 1it ) + γ 2 ln( V 2 it ) , A A TA TA ln( µ TA it ) = α 0 + α 1 + α 6 + ( α 2 + α 4 ) t + ( α 3 + α 5 )( t − t 0 i ) + γ 1 ln( V 1it ) + γ 2 ln( V 2 it ) . Hence, we have that ln(µ it ) − ln(µ it ) − ln(µ it ) + ln(µ it ) = α 4 (t A − t B ) + α 5 (t A − t 0i ) + α 6 + γ 1 E1it + γ 2 E 2it , TA TB CA CB where TA TB CA CB E 1it = ln( V 1it ) − ln( V 1it ) − ln( V 1it ) + ln( V 1it ) , TA TB CA CB E 2 it = ln( V 2it ) − ln( V 2 it ) − ln( V 2 it ) + ln( V 2it ) . The average of t B , which runs from 1 to t 0i − 1 , is 0.5 t 0 i , whereas the average of t A , which depends on s is t 0 i + 0.5( s + 1 ) , s = 1, 2 ,... . Thus, the averages of ( t A − t B ) and ( t A − t 0 i ) are 0.5( t 0 i + s + 1 ) and 0.5( s + 1 ) . Substituting the averages of ( t A − t B ) , ( t A − t 0 i ) , E 1it and E 2 it in Equation (75) leads to the decomposition in Equation (76). 183 Appendix C: Odds Ratio Decomposition under Nonlinear 'Koyck' Intervention Model Recalling Equation (82b) and going along the lines of Appendix B, we obtain −1 TB CA CB * A A * * A ln( µ TA it ) − ln( µ it ) − ln( µ it ) + ln( µ it ) = ω ( 1−δ B ) ( I it − φ I i ,t −1 ) + γ 1 E 1it + γ 2 E 2 it + φ F i ,t −1 , where E*1it = E 1it − φ E 1i ,t −1 , E*2 it = E 2it − φ E 2i ,t −1 and F i ,t −1 = ln( µ i ,t −1 ) − ln( µ i ,t −1 ) − ln( µ i ,t −1 ) + ln( µ i ,t −1 ) . TA A TB CA CB Note that −1 ω* ( 1−δ B ) ( I it − φ I iA,t −1 ) = ω* ( 1+δ B +δ 2 B2 + ...)( I it − φ I iA,t −1 ) , A A = ω* I itA + ω* ( δ − φ ) I iA,t −1 + ω* δ ( δ − φ ) I iA,t − 2 + ω* δ 2 ( δ − φ ) I iA,t −3 + ..., = ω* , s = 1, = ω* + ω* ( δ − φ ) , s = 2, = ω* + ω* ( δ − φ ) + ω* δ ( δ − φ ) , s =3, = ω* + ω* ( δ − φ ) + ω* δ ( δ − φ ) + ω* δ 2 ( δ − φ ) , s = 4 , etc. The average of ω* ( 1−δ B )−1( I it − φ I iA,t −1 ) , which depends on s is given by A ω* + ω* ( δ − φ ) [( s − 1 ) + ( s − 2 )δ + ( s − 3 ) 2 + ... + 2 s −3 + s −2 ] , δ δ δ s * =ω + ω* ( δ − φ ) s ∑ks −=11 ( s − k )δ k −1 =ω* + Since ∑ks −=11 δ k −1 = ∑ s −1 k =1 kδ k −1 ω* ( δ − φ ) {s s s = 1, 2 ,... . ∑sk−=11 δ k −1 − ∑sk−=11 k δ k −1}. δ ( 1 − δ s −1 ) 1 − δ s −1 ∂ s −1 k , ∑ks −=11 δ k = and ∑ks −=11 kδ k −1 = ∑k =1 δ , then 1−δ 1−δ ∂δ 1−δs s δ s −1 A = − . Therefore, the average of ω* ( 1−δ B )−1( I it − φ I iA,t −1 ) is 2 ( 1 −δ ) 1 − δ 184 ω + * ω* ( δ − φ ) s( 1 − δ s −1 − 1 − δ s + s δ s −1 = ω* ( 1 − φ ) + ω* ( φ − δ ) 1 − δ s , 2 2 s 1−δ s ( 1 −δ ) 1 − δ ( 1 −δ ) 1−δ s = 1, 2 ,... . Substituting the averages of ω* ( 1−δ B )−1( I it − φ I iA,t −1 ) , E*1it , E*2 it and F iA,t −1 in Equation A (75) leads to the decomposition in Equation (84). 185
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- New techniques for developing safety performance functions
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
New techniques for developing safety performance functions El-Basyouny, Karim 2011
pdf
Page Metadata
Item Metadata
Title | New techniques for developing safety performance functions |
Creator |
El-Basyouny, Karim |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | While motorized travel provides many benefits, it can also do serious harm in the form of road-related collisions. The problem affects millions of human lives and costs billions of dollars in economic and social impacts every year. The problem could be addressed thorough several approaches with engineering initiatives being recognized as the most sustainable and cost effective. However, the success of the engineering approaches in reducing collision occurrences hinges upon the existence of reliable methods that provide accurate estimates of road safety. Currently, Safety Performance Functions (SPFs) are considered by many as the main tool in estimating the safety levels associated with different road entities. Therefore, the research in this thesis focuses on addressing key issues related to the development of SPFs for i) collision data analysis and ii) collision intervention analysis. Some of the key issues addressed include: 1) adding spatial effects to SPFs thereby recognizing the evident spatial nature of road collisions, 2) fitting hierarchical models to allow inference to be made on more than one level, 3) recognizing the multivariate nature of collisions as most data are available by collision type or severity and modeling the data as such, 4) identifying and accounting for outliers in the development of SPFs, 5) developing a novel evaluation methodology to estimate the effectiveness of safety countermeasures when subject to data limitations, and 6) compare different tools for investigating the safety change in treated sites due to the implementation of safety countermeasures. The applications of the various models have been demonstrated using several collision datasets and/or safety programs. The results provide strong evidence for (i) incorporating spatial effects in SPFs, (ii) clustering road segments or intersections into homogeneous groups (e.g., corridors, zones, districts, municipalities, etc.) and incorporating random cluster parameters in SPFs, (iii) developing robust multivariate models with multiple covariates for modeling collisions by severity and/or type concurrently, and (iv) the effectiveness of the proposed full Bayes safety assessment methods that account for several theoretical and practical issues concurrently. In addition to the improvement in goodness of fit, the proposed models have also improved inference and precision of expected collision frequency. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-14 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0063049 |
URI | http://hdl.handle.net/2429/31254 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_spring_elbasyouny_karim.pdf [ 3.44MB ]
- Metadata
- JSON: 24-1.0063049.json
- JSON-LD: 24-1.0063049-ld.json
- RDF/XML (Pretty): 24-1.0063049-rdf.xml
- RDF/JSON: 24-1.0063049-rdf.json
- Turtle: 24-1.0063049-turtle.txt
- N-Triples: 24-1.0063049-rdf-ntriples.txt
- Original Record: 24-1.0063049-source.json
- Full Text
- 24-1.0063049-fulltext.txt
- Citation
- 24-1.0063049.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0063049/manifest