USING SPATIAL INTERPOLATION TO DETERMINE IMPACTS OF SNOWFALL ON TRAFFIC CRASHES By Gentjan Heqimi                A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Civil EngineeringMaster of Science 2016    ABSTRACT USING SPATIAL INTERPOLATION TO DETERMINE IMPACTS OF SNOWFALL ON TRAFFIC CRASHES By Gentjan Heqimi --   iii This thesis is dedicated to my family. Thank you for your support throughout all these years.  iv ACKNOWLEDGMENTS  I would like to express my sincere gratitude to my advisor Dr. Timothy J. Gates at Michigan State University for the continuous support and insight throughout not only this master thesis, but also throughout my entire experience . This was truly a very productive and extremely gratifying experience. More than I would have ever expected.                   v TABLE OF CONTENTS  LIST OF TABLES ......................................................................................................................... vi LIST OF FIGURES ...................................................................................................................... vii KEY TO ABBREVIATIONS ...................................................................................................... viii INTRODUCTION ...........................................................................................................................1 CHAPTER 1: LITERATURE REVIEW .........................................................................................3 1.1 General Snow Effects ....................................................................................................3 1.2 Weather Data Spatial Interpolation ................................................................................6 1.3 Crash Modeling ..............................................................................................................9 CHAPTER 2: OBJECTIVE & STUDY FRAMEWORK .............................................................11 CHAPTER 3: DATA COLLECTION ...........................................................................................15 3.1 Crash, Traffic, and Roadway Data ...............................................................................15 3.2 Snowfall Data...............................................................................................................19 3.2.1 Spatial Interpolation Implementation & Optimization .................................21 3.2.2 Spatial Interpolation Output ..........................................................................25 CHAPTER 4: STATISTICAL ANALYSIS ..................................................................................29 4.1 Statistical Methods .......................................................................................................29 4.1.1 Model 1  Natural Log of Snowfall ..............................................................29 4.1.2 Model 2  Categorized Snowfall ..................................................................30 4.2 Model 1  Natural Log of Snowfall .............................................................................31 4.2.1 Midblock Crashes .........................................................................................31 4.2.2 Interchange Crashes ......................................................................................33 4.2.3 Model 1 Results Discussion ..........................................................................35 4.3 Model 2  Categorized Snowfall .................................................................................36 4.3.1 Midblock Crashes .........................................................................................36 4.3.2 Interchange Crashes ......................................................................................38 4.3.3 Model 2 Results Discussion ..........................................................................40 4.4 Midblock Snowfall Crash Curves  ...............................................................................41 CHAPTER 5: SPATIAL PATTERN ANALYSIS  .......................................................................49 5.1 Crashes & Kernel Density Function ............................................................................49 5.2 Case Study 1  Interstate 75.........................................................................................52 5.3 Case Study 2  Interstate 94.........................................................................................54 CHAPTER 6: CONCLUSIONS ....................................................................................................58 REFERENCES ..............................................................................................................................62  vi LIST OF TABLES  TABLE 1  Characteristics of  the Study Network .........................................................................13 TABLE 2  2004-2014 Winter Crash Count Descriptive Statistics of Study Network...................16 TABLE 3  2004-2014 AADT, commercial AADT, Length, and Horizontal Curvature Descriptive Statistics of Study Network ...........................................................................................................17 TABLE 4  Ordinary Kriging Target Prediction Errors ..................................................................23 TABLE 5  Weather Stations and Prediction Errors for Ordinary Kriging Outputs .......................24 TABLE 6  Network Mean Predicted Snowfall (inches) for Winter Months 2004-2014 ...............27 TABLE 7  Model 1 (Natural Log of Snowfall) - Midblock Crashes NB Regression Results .......32 TABLE 8  Model 1 (Natural Log of Snowfall) - Interchange Crashes NB Regression Results ...34 TABLE 9  Model 1 Snowfall Elasticities by Crash Category .......................................................35 TABLE 10  Model 2 (Categorized Snowfall) - Midblock Crashes NB Regression Results .........36 TABLE 11  Model 2 (Categorized Snowfall) - Interchange Crashes NB Regression Results ......38 TABLE 12  Model 2 Snowfall Elasticities by Crash Category .....................................................40           vii LIST OF FIGURES  FIGURE 1  Michigan Freeway Study Network .............................................................................14 FIGURE 2  2004-2014 Network Average AADT Distribution .....................................................18 FIGURE 3  Weather Station Locations for a Typical Winter Month, January 2010 .....................20 FIGURE 4  Typical Ordinary Kriging Output for Snowfall, January 2010 ..................................26 FIGURE 5  Segment Snowfall Assignment Algorithm .................................................................27 FIGURE 6  Network Average Annual Snowfall during January, February, and December 2004-2014................................................................................................................................................28 FIGURE 7  Effect of Annual Snowfall on Winter Midblock Crashes Combined for Tangent and Curved Segments ...........................................................................................................................42 FIGURE 8  Effect of Annual Snowfall on Winter Midblock Crashes for Tangent and Curved Segments ........................................................................................................................................43 FIGURE 9  Effect of Annual  Snowfall on Non-Truck/Bus Winter Midblock Crashes for Tangent and Curved Segments ....................................................................................................................44 FIGURE 10  Effect of Annual Snowfall on Truck/Bus Winter Midblock Crashes for Tangent and Curved Segments ...........................................................................................................................45 FIGURE 11  Effect of Annual Snowfall on Injury Winter Midblock Crashes for Tangent and Curved Segments ...........................................................................................................................46 FIGURE 12  Effect of Annual Snowfall on PDO Winter Midblock Crashes for Tangent and Curved Segments ...........................................................................................................................47 FIGURE 13  Simplified Visual Representation of Kernel Density Process ..................................50 FIGURE 14  2004-2014 Network Winter Midblock Crash Rate & Average Annual Snowfall....51 FIGURE 15  2004-2014 I-75 Winter Midblock Crash Rate & Average Annual Snowfall ...........52 FIGURE 16  I-75 Kernel Density Distribution ..............................................................................53 FIGURE 17  2004-2014 I-94 Winter Midblock Crash Rate & Average Annual Snowfall ...........55 FIGURE 18  I-94 Kernel Density Distribution ..............................................................................56  viii KEY TO ABBREVIATIONS  AADT  Annual Average Daily Traffic BL  Business Loop BR  Business Route DOT  Department of Transportation GIS  Geographic Information Systems HSM  Highway Safety Manual MDOT  Michigan Department of Transportation MGF  Michigan Geographic Framework MISH  Meteorological Interpolation based on Surface Homogenized Data Basis MPH  Miles per Hour MSP  Michigan State Police NB  Negative Binomial NFC  National Functional Classification NOAA  National Oceanic and Atmospheric Administration PDO  Property Damage Only PR  Physical Road PRISM Parameter Elevation Regressions on Independent Slopes Model IDW  Inverse Distance Weighting ITS  Intelligent Transportation Systems US  United States   1 INTRODUCTION  Inclement weather can influence crash occurrence in a transportation system. Nearly one-quarter of all crashes in the United States (US) are a result of adverse weather patterns (FWHA, 2016). Severe weather conditions can impact vehicle performance, driver behavior, and infrastructure through reduced visibility, pavement friction performance, vehicle stability and maneuverability (Pisano et al., 2008; Liu, 2013; Leard and Roth, 2015). Anecdotal knowledge suggests that these effects are more pronounced during winter weather. Indeed, existing literature on adverse weather effects on crashes indicates that the  While several studies have evaluated the effects of snowfall on crashes and in particular on crash severity outcomes, there has been little consideration on the effects of snowfall on additional types of crashes as well as differentiations between the various facility types of the network on which these crashes occur. Accordingly, an investigation into the effects of snowfall on different crash categories and locational scenarios can provide revealing information that could prove useful to transportation agencies in implementing mechanisms to minimize crashes due to inclement weather conditions. These effects are particularly pertinent in regions like the Great Lake basin of the US where snow events occur frequently during winter months due to lake effect climactic conditions. Although drivers in these regions are presumed accustomed to driving in snowy weather, crash occurrence continues to be high (Andreescu and Frost, 1998). More recently, severe winter events have resulted in dramatic accidents such as the events of January 2015 in southwest Michigan, where extreme whiteout conditions led to a 193 car pile-up along Interstate 94 (I-94) resulting in loss of life, multiple injuries, and significant economic cost (Sell, 2016).   2 Based on this premise, this study aims to investigate the effects of annual snowfall on freeway crashes during winter periods in the State of Michigan for interchange and non-interchange (i.e. midblock) segment crashes based on spatial analytical methods and appropriate statistical techniques. These two locational circumstances are analyzed separately to account for the unique road, operational, and behavioral elements which defines and distinguishes them. Several winter crash categories are considered during the analysis. These include all winter crashes, crashes involving a truck or bus, crashes not involving a truck or bus, injury crashes (including fatalities), and Property Damage Only (PDO) crashes. The time period of the analysis consists of the winter months of January, February, and December for the 11-year period of 2004 to 2014. The study network consists of the entire freeway network of the State of Michigan, where freeways are defined as National Functional Class (NFC) equal to either Principal Arterial Interstate or Principal Arterial-Other Freeways. Cumulatively this represents approximately 1,955 miles of roadway throughout Michigan. The selected roadway network provides variability in snowfall along the network with snowfall typically peaking along the western coastline of the state and in the Upper Peninsula due to lake-effect snow impacting the region (Andresen et al., 2012). Additional Michigan historical weather conditions for the 1981 to 2010 winter months of January, February, and December indicate spatial variability in the mean temperature as well, with the northern half of the state experiencing mean temperatures in the 15-25 oF range and the southern half experiencing temperatures in the 25-30 oF range (Andresen et al., 2012). Precipitation, while varying in value between the northwest and southeast region of the State, occurs less often during winter months due to characteristic sub-freezing temperatures experienced during this time (Andresen et al., 2012).  3 CHAPTER 1: LITERATURE REVIEW  1.1 General Snow Effects Severe winter weather conditions can have direct and indirect effects upon the transportation system. These effects can impact vehicle performance, driver behavior, and the infrastructure through reduced visibility, pavement friction performance, roadway operations (i.e. travel speeds, roadway capacity, delay), vehicle stability and maneuverability, thus increasing the risk of accidents (Pisano et al., 2008; Liu, 2013; Leard and Roth, 2015).  In regions like the Great Lakes basin of the US, snowfall is often the primary weather factor impacting traffic operations and safety. Prior research has investigated the effects of snowfall on crashes utilizing various levels of detail and across several roadway types. It is commonly accepted that overall, snowfall results in an increase in the crash occurrence rate (Qiu and Nixon, 2008). These effects are observed through traditional crash modeling means such as Negative Binomial models specified by the Highway Safety Manual (HSM) (AASHTO, 2010), and through more unconventional means incorporating spatial statistical methods (Khan et al., 2008).   While overall crash risks are expected to increase during snow conditions, effects of snowfall on different vehicle types (i.e. commercial vehicles, non-commercial vehicles) and crash outcomes (i.e. fatal, injury, PDO)  may vary. With regards to crash outcomes, the literature appears to be in general agreement that snowfall results in an increase in PDO crashes (Eisenberg, 2004; Eisenberg and Warner, 2005; Blionis, 2013). These findings are not surprising given that, throughout snow conditions, accidents can occur at low speeds due to poor pavement friction performance and decreased vehicle control capabilities.    4 Comparatively, findings on injury crashes and fatal crashes are less decisive. For example Fridstrom et al. (1995) found a negative relationship between snow and injury crashes in a study on various locations in the Nordic region. By contrast Eisenberg (2004) found that snowfall results in increased numbers of non-fatal crashes. These findings are corroborated by additional studies which found that snowfall exhibits a positive significant impact on injury crashes, albeit at lower levels than PDO crashes (Eisenberg and Warner, 2005; Bilionis, 2013; Liu, 2013). The difference between non-fatal injury and PDO crashes may be attributed to heightened driver alertness and lower vehicle speeds to adjust to the adverse weather conditions (Eisenberg, 2004; Bilionis, 2013). Additionally Brorsson et al. (1988) suggested that lower levels of injuries involving single motor vehicles could be a result of snow walls developed along the roadway which help vehicles decelerate prior to collision. This effect would also translate in decreasing numbers of non-fatal injury crashes with increasing snow depth and consistent sub-zero temperatures.   The impact of snowfall on fatal crashes is similarly inconclusive. For instance Eisenberg and Warner (2005), in a study on the effects of snowfall on US crash rates for 1975-2000, found that snowfall does not result in an increase in fatal crashes. These findings are also supported in Eisenberg (2004) for monthly snowfall values. Comparatively, Perry and Symons (1991) found fatal crash rates to increase for snowy days in the UK. Eisenberg (2004) also found an increase in fatal crashes when assessing daily snowfall effects on crashes in the US. The author however attributes this positive relationship to a lag effect on precipitation (i.e. crashes decrease when it snows every day and increase significantly when the time since the last snowfall increases). This lag effect is illustrated in additional studies (Eisenberg and Warner, 2005; El-Basyouny et al., 2014; Seeherman and Liu, 2015) and indicates that snow events which follow a dry season tend  5 to be more dangerous, a result likely related to drivers being unaccustomed to sudden changes in driving conditions.  Parallel variances can also be found in crashes involving commercial vehicles and non-commercial vehicles whether due to differences in trip characteristics or mechanistic characteristics. For example, drivers perceive severe snow events as dangerous and may avoid unnecessary trips. Commercial vehicle trips however are commonly business related and thus less flexible in route choice (Pisano et al., 2008). Truck performance is also more vulnerable to decreasing visibility as it can affect stopping distance more significantly than smaller and lighter vehicles, among other physical characteristics attributed to weight and truck performance (Pisano et al., 2008). Effects of snowfall on crashes are additionally shown to exhibit a non-linear relationship with crash occurrences. Eisenberg (2004) found crashes to peak at moderate snowfall levels and decline at higher ranges due to potential reductions in travel as weather severity increases. Additional studies show similar relationships where larger amounts of snowfall do not necessarily translate to higher crash occurrences (Khattak and Knapp, 2001), or they exhibit decreasing marginal impacts with increasing snowfall amounts (Seeherman and Liu, 2015). In circumstances where crashes continue to increase at larger snowfall levels, the data may be displaying effects of underlying events such as higher frequencies of snowstorms and/or storms of greater intensity.  These particular weather patterns are shown to result in higher crash frequencies (Khattak and Knapp, 2001).  A secondary noteworthy dimension towards the effects of snowfall on crashes is the facility in which these crashes occur. Two particular facility types can be defined in a freeway setting based on their unique characteristics of geometry, operational behavior and driver  6 behavior: the interchange and non-interchange (i.e. midblock) segment area of the freeway. Existing literature suggests that interchange and non-interchange segments should be treated separately in crash modeling as crash rates in the interchange sphere of influence are much higher than their non-interchange counterpart (Kiattikomol, 2005). Given their distinctive characteristics, the effects of snowfall on interchange versus non-interchange segment crashes may vary as well. 1.2 Weather Data Spatial Interpolation Since historical snowfall data is a key component of this study, the estimation of historical snowfall amounts along the study network must be reliable and reflect with a degree of accuracy likely conditions experienced by drivers along each point on the network for that particular time period. Presently there is no practical method to explicitly measure such data continuously along an expansive freeway system. To overcome this barrier snowfall data are commonly estimated from weather stations which are randomly distributed throughout the US (Bostan et al., 2012). Each of these stations captures and reports weather conditions on that particular location over pre-established time periods. Because distances between these stations may vary in magnitudes of less than a mile to several miles, values reported by one station do not provide the level of accuracy required to estimate snowfall when applied to specific locations along a freeway. Additionally, weather patterns can be localized in small geographic regions that experience relatively different snowfall amounts among each other or their surroundings due to lake-effect climatic conditions, wind patterns, or terrain. Thus it is naive to average these values across a specific region (Leard and Roth, 2015).  To overcome these shortfalls weather data prediction models must incorporate the density and locations of weather stations to obtain accurate weather values on desired localities (Ashraf  7 et al., 1997). Spatial interpolation is the procedure utilized in incorporating these variables to estimate weather related data at specific locations. Sluiter (2009) states that spatial interpolation methods can be grouped into three categories: deterministic, probabilistic, and other; where deterministic methods produce continuous surfaces based on specific geometric characteristics of existing observations (i.e. nearest neighborhood, triangulation, Inverse Distance Weighting (IDW), splines, linear regression); probabilistic methods produce continuous surfaces based on statistical theory (i.e. Ordinary Kriging, Simple Kriging, Universal Kriging); and other consists of a combination of deterministic and probabilistic methods (i.e. Parameter Elevation Regressions on Independent Slopes Model (PRISM), Meteorological Interpolation based on Surface Homogenized Data Basis (MISH)).  Traditionally deterministic methods have an extended history of use in predicting meteorological data. In recent trends however, probabilistic methods have become a more preferred approach in predicting weather values since they provide statistical reliability, consider the spatial correlation between the observations, and allow for the inclusion of secondary explanatory variables (i.e. elevation) in improving the estimation on unknown locations (Mair and Fares, 2011). While more demanding to implement, several studies have found that probabilistic methods deliver superior estimates than their deterministic counterparts. For example Mair and Fares (2011) found Ordinary Kriging and Simple Kriging with local means outperform Thiesen polygons, IDW, and linear regression in rainfall estimations in a mountainous region in Hawaii. Tabios and Salas (1985) found comparable findings in their review of 30 years of annual precipitation in Region II of North Central US, with Kriging and other optimization interpolation techniques outperforming deterministic approaches (i.e. Thiesen polygons, IDW, Lagrange approach). Ashraf et al. (1997) in their review of two years of daily  8 climate data in three states in the US found that Kriging produces the lowest root mean square interpolation error compared to inverse distance methods. Goovaerts (2000) found probabilistic methods (i.e. Ordinary Kriging, Co-Kriging, Simple Kriging with local means) outperform Thiessen polygons, IDW, and linear regression in rainfall estimations even when elevation data is not incorporated. When elevation is included in the process, results produce more reliable values. Similar findings are reflected in Huang et al. (2015) in their analysis of daily snow depth, with the authors noting that Ordinary Kriging produces the best results when elevation is not correlated with predicted snow depth. In contrast, Mair and Fares (2011) did not find any accuracy improvement in Ordinary Kriging results when secondary data (i.e. elevation) are highly correlated to the predictor, thus indicating that different and more advanced types of Kriging techniques may be more appropriate when secondary variables are considered (i.e. Co-Kriging, Universal Kriging).  Despite the increasing widespread use of the probabilistic family of Kriging, several authors underscore the importance of using proper care throughout the estimation process as these techniques are fundamentally linear optimization processes (Lanciani and Salvati, 2008; Oliver and Webster, 2015). In its most basic form, Kriging can be defined as a linear and optimization spatial interpolation method used in estimating unknown values from existing observations where near sample points receive more weight than the ones further away (Oliver and Webster, 2015; Sluiter, 2016). Weights are modeled through semivariograms which identify the spatial correlation between observations (Oliver and Webster, 2015; Sluiter, 2016). The semivariogram is the basic framework of Kriging which defines how space fits the data. Thus if the fit is erroneous, Kriging predictions will not be accurate. (Ali, 2013; Oliver and Webster, 2015). To avoid biased Kriging outputs, Oliver and Webster (2015) recommend using an  9 adequate number of data points, selecting the proper semivariogram form (i.e. lag interval, bin width), transforming data to achieve a near normal distribution when significant skewness is present, careful examination of outliers to identify erroneous observations, and de-trending data when significant trends are evident, depending on the Kriging method used. 1.3 Crash Modeling Crash frequency is the primary indicator of safety of a roadway, where crash frequency is defined as the number of crashes occurring in a site over a period of one year (AASHTO, 2010). Crash modeling entails estimating the average expected crash frequency of a particular site given its specific geometric, operational, and local conditions over a pre-defined time period (AASHTO, 2010).  Typically, such traffic crashes are assumed to be random occurrences in time and space (AASHTO, 2010; Bilionis, 2013; Zou et al., 2015). Since crashes represent counts of specific events, they have been originally presumed to follow a Poisson distribution where the mean is equal to the variance of the data (Kim et al., 2010; Lord and Mannering, 2010; Bilionis, 2013; Seeherman and Liu, 2015). In these cases crash modeling has been conducted via a Poisson regression model (Lord and Mannering, 2010). The equivalency of the mean and the variance however is not always achieved, as researchers have found that crashes often exhibit extra dimensions of variance or over-dispersion. To account for this inequality where the variance exceeds the mean, crash modeling has been most often conducted by employing a Negative Binomial regression model (AASHTO, 2010; Kim et al., 2010; Lord and Mannering, 2010; Bilionis, 2013; Seeherman and Liu, 2015). The Negative Binomial model is an extension of the Poisson model which accounts for the over-dispersion in the data by including a gamma distributed error term with mean 1 and variance  (Lord and Mannering, 2010). While the Negative Binomial model has consistently been the most often used method for crash modeling,  10 it falls short when the data is under-dispersed or the sample size is too small (Lord and Mannering, 2010). To account for some of these shortfalls additional alternatives to the Negative Binomial model have been used like the Poisson-Lognormal model and Conway-Maxwell-Poisson model (Lord and Mannering, 2010). There are also cases however where there may be time periods or sections of a study site that exhibit a high frequency of no crashes. These could occur when the time period or the dataset utilized is relatively small. In these cases a zero-inflated regression model may be most appropriate to account for the recurring zero variables (Lord and Mannering, 2010; Bilionis, 2013).                 11 CHAPTER 2: OBJECTIVE & STUDY FRAMEWORK  The objective of this study is to investigate the effects of annual snowfall on freeway crashes during winter periods in the State of Michigan. Separate analyses are performed for interchange and non-interchange segment (i.e. midblock) crashes due to their differing geometric, operational, and behavioral characteristics. These two primary crash classifications are based on the annual Michigan State Police (MSP) statewide crash database codebook and can be defined as follows:   Midblock crashes  Traffic crashes not associated with an interchange or intersection (typically occurring between two interchanges)  Interchange crashes  Traffic crashes associated with an interchange (typically occurring between ramp termini) The specific crash categories analyzed include:  All winter crashes  All crashes occurring on the study network during January, February, and December between 2004 and 2014.  Truck/bus winter crashes  Crashes involving at least one truck or bus during January, February, and December between 2004 and 2014.  Non-truck/bus winter crashes  Crashes involving no trucks or buses during January, February, and December between 2004 and 2014.  Injury winter crashes  Crashes where at least one injury (or fatality) was reported during January, February, and December between 2004 and 2014.  PDO winter crashes  Crashes involving only property damage during January, February, and December between 2004 and 2014  12 The inclusion of the listed crash categories as well as their differentiation between interchange and non-interchange (i.e. midblock) crashes can contribute and potentially fill a void in the existing literature with regards to the impacts of snowfall on these crash scenarios.   The effects of snowfall on winter-season crashes are assessed in terms of the estimated annual winter snowfall occurring along the study network. The relationship between snowfall and crashes is further examined by categorizing annual snowfall in terms of its quartile intervals. The time period of the analysis includes the months of January, February, and December during the 11-year period of 2004 and 2014, where the selected months represent the winter period in Michigan which historically experiences the most snowfall amounts. Accordingly data used in this study correspond exclusively to this time period.  The study network consists of the entire limited access freeway network in the State of Michigan (i.e. NFC equal to Principal Arterial-Interstate or NFC equal to Principal Arterial-Other Freeways). This system contains approximately 1954.9 one directional miles of freeway which are primarily concentrated in the southern half of Michigan and decrease towards the more remote northern parts of the state. The composition of this network contains a mix of Interstate, Interstate Business Loop (BL), US, US Business Route (BR), State, and Connector routes. The dominant speed limit on these routes is 70 mph (60 mph for trucks); with a relatively small group of segments comprising an estimated five percent of the total mileage of the network having posted speed limits of less than 70 mph. The latter occurs only within selected major urban areas of the state.  The study network is additionally comprised of 2,398 discrete freeway segments of lengths varying from 0.034 to 9.35 miles. These segments are based on the 2014 Michigan Department of Transportation (MDOT) sufficiency database and represent the base framework  13 of the study. All data used in this study are binned on these segments prior to any crash analysis. Characteristics of these freeway segments are presented in Table 1, while the freeway network is illustrated in its entirety in Figure 1. TABLE 1  Characteristics of the Study Network Route Total Miles* Mean Number of Lanes Typical Speed Limit Mean Shoulder Width Mean AADT* Mean CMV AADT* Interstate 1241.2 2.6 70 mph 9.5 30,658 3,117 US 564.3 2.1 70 mph 9.0 15,599 1,385 Michigan 119.1 2.6 55/70 mph 8.7 35,742 1,231 Interstate - BL 3.6 2.0 55/70 mph 8.8 16,293 367 US - BR 18.5 1.9 55/65/70 mph 8.0 8,933 699 Connector 8.2 1.9 70 mph 9.1 8,569 421 All Routes 1954.9 2.5 70 mph 9.3 27,137 2,476 Lane width = 12 ft in nearly all segments       *Includes both directions        14  FIGURE 1  Michigan Freeway Study Network    15 CHAPTER 3: DATA COLLECTION  The data used in this study includes traffic crashes, snowfall, directional Annual Average Daily Traffic (AADT) volumes, directional commercial AADT volumes, horizontal curvature, and segment length for each of the 2,398 freeway segments during the 11-year period of 2004 to 2014. With the exception of snowfall which represents the primary explanatory variable of interest, the other four variables are included based on their notable role in predicting crashes and their common use in crash modeling. Specifically AADT or commercial AADT represent the primary exposure variable for the corresponding crash category, length acts as a normalizing explanatory variable due to the varying lengths of the segments, and horizontal curvature is included based on the relationship between road alignment and crash occurrence. These variables are used in prior studies investigating the impact of weather on crashes and are shown to be significant factors in crash occurrence under these conditions (Khattak and Knapp, 2001; Ahmed et al., 2011; Ahmed et al., 2012; Bilionis, 2013). Rainfall was also included in the initial investigation phase of the study, however it exhibited a high degree of correlation with snowfall thus was removed from the analysis to avoid multicollinearity bias.  3.1 Crash, Traffic, and Roadway Data The crash data is obtained from the annual MSP statewide crash database for the months of January, February, and December for the 11-year period of 2004 to 2014. The crash categories are isolated for both midblock and interchange scenarios based on the coded values provided in the MSP database, and are matched to the corresponding freeway segments as annual winter-season totals based on their Physical Road (PR) identification number and mile point. The PR system is a linear referencing system which distinctively identifies roadway events (i.e. crashes,  16 segments) s transportation network. Descriptive statistics for these crashes are provided in Table 2.  TABLE 2  2004-2014 Winter Crash Count Descriptive Statistics of Study Network Category 2004-2014 Winter Crash Totals  Winter Crashes per Segment per Year Count Percent  Min Max Mean St Dev All Crashes 125,665 100%  0 61 4.76 5.13 Non-Truck/Bus 114,824 91%  0 49 4.35 4.76 Truck/Bus 10,841 9%  0 18 0.41 0.84 Injury 24,462 19%  0 17 0.93 1.38 PDO 101,203 81%  0 50 3.83 4.20 Cumulative  100%      All Crashes 49,944 100%  0 43 1.89 3.02 Non-Truck/Bus 45,377 91%  0 39 1.72 2.76 Truck/Bus 4,567 9%  0 15 0.17 0.55 Injury 9,640 19%  0 15 0.37 0.82 PDO 40,304 81%  0 40 1.53 2.50 Midblock  40%*      All Crashes 74,141 100%  0 44 2.81 3.88 Non-Truck/Bus 67,977 92%  0 39 2.58 3.62 Truck/Bus 6,164 8%  0 7 0.23 0.59 Injury 14,543 20%  0 15 0.55 1.07 PDO 59,598 80%  0 34 2.26 3.15 Interchange  59%*      1% of cumulative crashes are non-traffic coded crashes *Relative to the cumulative 2004-2014 winter crash population  Overall, a total of 125,665 crashes occurred on the study network for the 2004-2014 winter-season period. Interchange crashes comprise the majority of the crashes with approximately 59% of the total amount; while approximately 40% represent midblock crashes. The difference or approximately 1% are crashes coded as non-traffic and are not included in this study. With respect to crash types, the distribution is nearly equivalent between the three crash categories (i.e. all winter, midblock, interchange) and ranges in the 91-92% for non-truck/bus crashes vs 8-9% for truck/bus crashes, and 80-81% for PDO crashes vs 20-21% for injury crashes.   17 With respect to traffic volumes, directional AADT and directional commercial AADT values are obtained from the 2014 MDOT sufficiency database for all freeway segments for each year of the 2004 to 2014 time period. These values represent the average AADT (or commercial AADT) experienced on each of the freeway segment for each year of the study period.  Horizontal curvature is extracted from a comparable database developed by calculating the radii of curved segments from the 2014 Michigan Geographic Framework (MGF) shapefile and adapting it to the 2014 MDOT sufficiency database segments. This database provides horizontal curvature information binned in various design speed of curve formats by assuming a 7% maximum superelevation, which corresponds to the maximum superelevation allowed on etwork. In this study, this variable is presented as the fraction of the segment with a horizontal curve with a design speed less than 85 mph and is assumed temporally constant.  Lastly, the length of each of the 2,398 freeway segments is calculated based on their geographic length in miles using the ArcGIS for Desktop software. Table 3 presents descriptive statistics for these four variables, while Figure 2 illustrates the overall average AADT distribution throughout the study network.  TABLE 3 2004-2014 AADT, commercial AADT, Length, and Horizontal Curvature Descriptive Statistics of Study Network  Per Segment per Year Category Min Max Mean St Dev AADT 100 107,000 27,144 21,472 Commercial AADT 25 8,846 2,477 1,771 Segment Length (miles) 0.03 9.35 1.62 1.38 Fraction of Horizontal Curvature with design speed less than 85 mph 0.00 1.00 0.11 0.24  18  FIGURE 2  2004-2014 Network Average AADT Distribution The average AADT distribution along the network indicates that the highest traffic volumes are experienced southeast of the state in the Metro Detroit region and tend to dissipate outwards of this area. Significant volumes are also experienced in regions surrounding primary    19 cities such as Grand Rapids, Lansing, Kalamazoo, Flint, and Saginaw. Comparatively, areas experiencing the lowest volumes are located primarily in the northern half of the state given the  3.2 Snowfall Data To assess historical snowfall conditions along the study network, monthly snowfall data is extracted from the National Oceanic and Atmospheric Administration (NOAA) climate data center for January, February, and December of the 11-year period of 2004 to 2014. The climate data center provides weather data as captured and reported by weather stations across the US for their respective location and time period selected. These stations are randomly located with varying distances between each other ranging from less than a mile to several miles. A similar pattern is observed for weather stations surrounding the study network (Figure 3).  Consequently, the random distribution of the stations necessitates the use of spatial interpolation methods to predict likely historical snowfall values throughout the network with a degree of accuracy. While a multitude of spatial interpolation methods exist, the probabilistic method of Ordinary Kriging is employed due to its consistent and superior performance in estimating precipitation versus deterministic methods such as IDW, nearest neighborhood, splines, and linear regression (Tabios and Salas, 1985; Ashraf et al, 1997; Goovaerts, 2000; Mair and Fares, 2011; Huang et al. 2015), as well as due to its prevalent use in spatial interpolation as it requires a minimal number of assumptions to be satisfied (Oliver and Webster, 2015).   20  FIGURE 3  Weather Station Locations for a Typical Winter Month, January 2010  21 3.2.1 Spatial Interpolation Implementation & Optimization Ordinary Kriging is one form of Kriging used to predict values from known observations, where the mean is assumed unknown, constant, and is estimated from the local neighborhood (Oliver and Webster, 2014). Since Ordinary Kriging is a local neighborhood estimator, nearby values closest to the subject location receive more weight than the ones further away. The weights are identified through a spatial function called a semivariogram which assesses the spatial correlation of the existing observations and the extent of this correlation (Oliver and Webster, 2015). In its most basic form, the semivariogram relates graphically the semivariance of the observations to the distances (lag) between these observations (ArcGIS 9: Using ArcGIS Geostatistical Analyst, 2003). The primary goal then in this process is to identify the best fit for this structure by selecting the optimal semivariogram model, where the fit is sensitive to the sample size of the observations, lag interval and size (i.e. grouping of observations based on distance to help identify spatial correlation patterns), data distribution where a normal distribution is preferred, outliers, and trends in the data (Oliver and Webster, 2015).  While the semivariogram identifies the weights used in the estimation process, the predicted values of any unknown location in the Ordinary Kriging process are defined in the following general form (Huang et al., 2015; Oliver and Webster, 2015):   (1) Where,  = predicted monthly snowfall at unknown location  = observed monthly snowfall at known location   = estimated weight for , where   = number of weather stations  22 Typically at least 50 observations points are required to produce accurate predictions (Holdaway, 1996), although 100 to 150 observations may be preferred (Oliver and Webster, 2015). Ordinary Kriging in this study is implemented via the ArcGIS geostatistical extension. Predictions are made for each individual winter month of the 2004-2014 time period by using all stations reporting snowfall data in Michigan as well as stations located within approximately 50 miles of its borders. This additional extent is deemed appropriate through observational examinations of the region as it provides an adequate number of weather stations to use for interpolation even in the more remote endpoints of the study network. Outliers are also examined and removed where applicable. In this context, outliers are defined as those stations which report monthly snowfall values of 0 inches, are surrounded by at least two nearby stations reporting significant monthly snowfall values, and repeat this pattern for at least one more month in the study time period. Lastly, data is further de-trended where appropriate and normalized if possible to facilitate optimized Ordinary Kriging outputs. To further optimize the models and resulting prediction accuracy, output errors are inspected and minimized via the cross-validation process.  The cross-validation process is a leave-one-out algorithm which withholds an observation point, makes a prediction about that point, and compares and validates the two values by outputting various error measurements. This procedure is then looped for the remaining observations (Ali, 2013; Laaha et al, 2013; Laaha et al., 2014; Oliver and Webster, 2015). The error output contains the following measurement errors which assist in optimizing and selecting the ideal prediction model (ArcGIS 9: Using ArcGIS Geostatistical Analyst, 2003):   23     Where,  = predicted monthly snowfall at location   = observed monthly snowfall at location   = predicted standard error at location  (measure of uncertainty between the true and predicted value at location , represented by the squared root of the prediction variance at location )  = number of weather stations Target prediction errors needed to obtain ideal prediction values are listed in Table 4 (ArcGIS 9: Using ArcGIS Geostatistical Analyst, 2003; Oliver and Webster, 2015). These values are used as reference throughout the Ordinary Kriging optimization process for each monthly output case.  TABLE 4  Ordinary Kriging Target Prediction Errors Mean Error Root-Mean-Square Error Average Standard Error Mean Standardized Error Root-Mean-Square Standardized Error very small (0) very small (0), within range of Average Standard Error very small  (0), within range of Root-Mean-Square Error very small (0) 1   24 The number of weather stations used in each monthly output in the Ordinary Kriging process along with the five corresponding prediction errors are listed in Table 5.   TABLE 5  Weather Stations and Prediction Errors for Ordinary Kriging Outputs Year Month No. of  Stations Error Mean RMS Mean Std RMS Std Avg Std 2004 January 211 -0.020 9.139 -0.003 1.211 7.327  February 212 0.056 5.248 0.007 1.177 4.420  December 214 0.018 7.239 -0.004 1.146 6.209 2005 January 209 -0.084 5.983 -0.012 0.989 6.013   February 209 0.017 1.058 0.017 1.101 0.959   December 218 0.007 7.327 0.003 1.148 6.234 2006 January 212 0.158 3.889 0.005 0.995 5.548  February 215 -0.109 7.218 -0.020 1.243 5.723  December 229 -0.028 4.761 -0.008 1.158 3.944 2007 January 227 -0.016 7.138 -0.004 1.240 5.599   February 227 -0.015 7.063 -0.001 1.177 5.764   December 224 0.009 6.474 0.004 1.198 5.340 2008 January 225 -0.058 8.586 -0.006 1.185 7.126  February 225 -0.058 6.845 -0.007 1.174 5.744  December 241 0.069 10.801 0.004 1.084 9.754 2009 January 239 -0.051 8.601 -0.003 1.116 7.415   February 246 0.020 6.322 0.003 1.197 5.153   December 249 -0.053 8.734 -0.004 1.144 7.508 2010 January 253 -0.045 6.845 -0.004 1.063 6.377  February 255 -0.051 6.359 -0.005 1.098 6.263  December 246 0.009 7.658 -0.003 1.087 6.985 2011 January 265 0.019 11.375 0.003 1.130 10.068   February 261 -0.107 8.289 -0.012 1.005 8.247   December 268 -0.001 3.676 0.000 1.027 3.566 2012 January 264 0.016 8.932 0.002 1.119 7.955  February 269 -0.056 5.082 -0.010 1.029 4.943  December 275 -0.011 4.777 -0.001 1.053 4.548 2013 January 261 -0.052 8.111 -0.004 1.121 7.204   February 260 0.016 10.470 0.002 1.045 10.048   December 262 -0.124 10.812 -0.010 1.063 10.162 2014 January 254 -0.030 14.446 0.000 1.016 14.277  February 255 -0.058 7.876 -0.007 0.981 8.060   December 259 0.042 4.587 0.010 1.180 3.827 Average 241 -0.017 7.325 -0.002 1.112 6.615   25 In all cases the minimum number of stations reporting data is greater than 200, thus meeting the observation sample recommended to produce accurate predictions. Similarly, while fluctuations exist among the individual monthly Ordinary Kriging outputs, obtained average prediction errors are within range of the target values for all five error measurements. These fluctuations are largest in instances where weather stations in proximity of one another report values of significant differences (i.e. 40 inch snowfall on Station 1 vs 5 inch snowfall on Station 2) and cannot be identified as outliers. 3.2.2 Spatial Interpolation Output The output of each monthly Ordinary Kriging model consists of a continuous raster surface covering the entire study network, the cells of which correspond to predicted monthly snowfall values in inches. Following each raster output, predicted values are than assigned to the endpoints of each roadway freeway segment and averaged for that segment in Geographical Information Systems (GIS) space. To ensure uniformity and integrity throughout the monthly snowfall segment assignment, any segment greater than 0.25 miles is split in 0.25 mile intervals. Predicted values for these segments are than assigned to the endpoints of the 0.25 mile breaks and averaged to obtain the final predicted monthly snowfall along that segment. Figure 4 illustrates a typical Ordinary Kriging output for monthly snowfall; Figure 5 illustrates the developed GIS algorithm for segment snowfall assignment, while Table 6 presents the mean predicted snowfall in inches for each 2004-2014 winter month for all freeway segments (i.e. network).     26  FIGURE 4  Typical Ordinary Kriging Output for Snowfall, January 2010   27  FIGURE 5  Segment Snowfall Assignment Algorithm TABLE 6  Network Mean Predicted Snowfall (inches) for Winter Months 2004-2014   2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 Avg January 26.3 22.2 6.8 12.0 19.6 23.4 9.8 16.3 12.0 8.3 25.3 16.5 February 6.7 0.0 8.6 14.3 23.9 9.3 18.8 21.8 7.8 17.4 16.8 13.2 December 12.8 20.2 5.9 16.2 32.6 13.0 10.3 3.3 6.4 14.3 1.0 12.4 Monthly Avg 15.3 14.1 7.1 14.1 25.4 15.3 13.0 13.8 8.7 13.3 14.4 14.0 Annual Total 45.8 42.4 21.3 42.4 76.1 45.8 38.9 41.5 26.1 40.0 43.1 42.1  Lastly, snowfall values are converted into annual totals for use with the crash model. Figure 6 below presents the average 2004-2014 annual winter snowfall distribution along the network. The snowfall distribution indicates that historically the western part of mainland Michigan and the Upper Peninsula experience more snowfall than the rest of the network. Snowfall tends to decline moving towards the eastern portion of the network in mainland Michigan. This variation could be explained by the Great Lakes lake-effect and is consistent with other documented historical snowfall trends for this region (Andresen et al., 2012).   28  FIGURE 6  Network Average Annual Snowfall during January, February, and December 2004-2014      29 CHAPTER 4: STATISTICAL ANALYSIS  4.1 Statistical Methods A series of regression models are conducted to investigate snow effects on crash types and crash outcomes for midblock and interchange crashes along the network. Since all of the crash categories in the dataset reflect a Poisson distribution with over-dispersion (i.e. variance > mean), the Negative Binomial regression model is employed for each case. The basic form of this model can be expressed as:   (7) Where,  = expected number of crashes of segment i  = explanatory variables  = regression coefficient  = gamma distributed error term with mean 1 and variance  ( = dispersion paramater Two Negative Binomial models are developed to assess the effects of snowfall on crashes. The first model (i.e. Model 1) aims to assess the overall effects of snowfall on crashes and represents the primary crash model in this study.  The second model (i.e. Model 2) aims to complement the first model by further investigating the relationship between snowfall and crashes. In this context, Model 1 represents the more preferred approach in evaluating these effects, with the two models differing solely based on the format of the snowfall variable. 4.1.1 Model 1  Natural Log of Snowfall The first model aims to assess the overall effects of snow on crashes using the natural log of snowfall amounts. Explanatory variables include annual winter snowfall, AADT or  30 commercial AADT (where commercial AADT is applicable for truck/bus crashes only as it is a more representative exposure variable for these types of vehicles), length, and horizontal curvature. The variables of snowfall, AADT or commercial AADT, and length are included in the model in natural log form, while horizontal curvature is expressed as a fraction. Thus Model 1 can be stated as:   (8) Where,  = expected number of crashes of segment i per year    = fraction of segment i with horizontal curve w/ design speed < 85 mph  = length of segment i in miles   = annual average daily traffic of segment i for year OR commercial annual average daily traffic of segment i for year (applicable for truck/bus crashes only)  = annual winter snowfall of segment i in inches   = model 1 regression coefficients The inclusion of the snowfall variable in natural log form further facilitates the interpretation of the regression coefficients, where the latter represents the percent change in crashes associated with a one-percent change in the snowfall variable.  4.1.2 Model 2  Categorized Snowfall The second model aims to investigate more closely the relationship between snowfall and crashes. This is achieved by breaking down annual snowfall amounts into three categories, or dummy variables, based on its four quartiles, where:  1st and 2nd annual snowfall quartile: 5.2 in   X < 37.4 in  3rd annual snowfall quartile: 37.4 in   X < 50.2 in  4th annual snowfall quartile: 50.2 in   X < 157.1 in  31 AADT or commercial AADT, and length are retained in the model in natural log form, horizontal curvature is retained as a fraction, while snowfall is expressed as a binary (0/1) variable for each of the three categories. The equation for Model 2 then becomes:    (9) Where,  = expected number of crashes of segment i per year  = fraction of segment i with horizontal curve w/ design speed < 85 mph  = length of segment i in miles   = annual average daily traffic of segment i for year OR commercial annual average daily traffic of segment i for year (applicable for truck/bus crashes only)  = categorical variable (0/1) of segment i experiencing annual snowfall amounts in inches which fall in the 3rd or 4th quartile range. Reference category is the bottom two quartiles.  = model 2 regression coefficients The interpretation of the categorical snowfall variable coefficients in this case can be facilitated by taking their exponent, subtracting one, and multiplying by 100, where the obtained value than represents the percent change in crashes relative to the reference point (i.e. 1st and 2nd annual snowfall quartiles).  4.2 Model 1  Natural Log of Snowfall 4.2.1 Midblock Crashes The Model 1 results for midblock crashes are provided in Table 7. These results indicate that the effects of snowfall on winter crashes differ in magnitude between crash types and crash severity outcomes. Each obtained equation along with the variables of length, AADT or commercial AADT, horizontal curvature, and snow are statistically significant at the 0.01  32 significance level.  Furthermore, all variable coefficients are positive indicating increasing crash occurrence with increasing values of the explanatory variables.  TABLE 7  Model 1 (Natural Log of Snowfall) - Midblock Crashes NB Regression Results Variable B Std Error Lower CI Upper CI Wald Chi-Square df p-value All Winter Crashes Intercept -8.429 0.146 -8.715 -8.143 3329.615 1 <0.001 Ln Length 1.210 0.012 1.187 1.234 10012.338 1 <0.001 Ln AADT 0.667 0.011 0.645 0.689 3533.089 1 <0.001 Horiz. Curve < 85 mph 0.497 0.037 0.425 0.569 181.405 1 <0.001 Ln Snow 0.501 0.017 0.468 0.534 876.531 1 <0.001 Negative Binomial 0.781 0.015 0.753 0.811 na na na Non-Truck/Bus Crashes Intercept -8.464 0.149 -8.755 -8.173 3240.817 1 <0.001 Ln Length 1.213 0.012 1.189 1.237 9637.573 1 <0.001 Ln AADT 0.661 0.011 0.639 0.683 3384.573 1 <0.001 Horiz. Curve < 85 mph 0.533 0.038 0.459 0.606 201.408 1 <0.001 Ln Snow 0.498 0.017 0.464 0.532 827.817 1 <0.001 Negative Binomial 0.771 0.015 0.742 0.802 na na na Truck/Bus Crashes Intercept -13.235 0.276 -13.776 -12.694 2298.446 1 <0.001 Ln Length 1.220 0.028 1.166 1.274 1962.894 1 <0.001 Ln Commercial AADT 1.119 0.027 1.067 1.171 1776.413 1 <0.001 Horiz. Curve < 85 mph 0.563 0.092 0.382 0.744 37.180 1 <0.001 Ln Snow 0.596 0.037 0.525 0.668 266.987 1 <0.001 Negative Binomial 0.864 0.061 0.752 0.992 na na na Injury Crashes Intercept -10.943 0.245 -11.423 -10.462 1991.654 1 <0.001 Ln Length 1.245 0.021 1.204 1.285 3639.899 1 <0.001 Ln AADT 0.765 0.019 0.729 0.802 1679.932 1 <0.001 Horiz. Curve < 85 mph 0.580 0.061 0.461 0.700 91.174 1 <0.001 Ln Snow 0.463 0.028 0.409 0.518 276.768 1 <0.001 Negative Binomial 0.788 0.036 0.720 0.863 na na na PDO Crashes Intercept -8.604 0.153 -8.903 -8.305 3174.904 1 <0.001 Ln Length 1.213 0.013 1.188 1.238 9129.151 1 <0.001 Ln AADT 0.655 0.012 0.632 0.678 3138.357 1 <0.001 Horiz. Curve < 85 mph 0.467 0.039 0.391 0.543 145.291 1 <0.001 Ln Snow 0.522 0.018 0.488 0.557 869.662 1 <0.001 Negative Binomial 0.753 0.016 0.723 0.785 na na na na = not applicable; CI = Confidence Interval; df = degrees of freedom           33 Overall, winter crashes increase by 0.5% for each one-percent increase in annual snowfall. Snow effects are most pronounced for those winter crashes involving a truck or bus, which show the most elastic relationship where each one-percent increase in snowfall results in a 0.6% percent increase in truck/bus crashes. Conversely, the effect of snow on crashes not involving a truck or bus experiences a lower elastic relationship with an increase in approximately 0.5% for each additional one-percent increase in snowfall. Among crash severity outcomes, injury crashes show a less pronounced effect than PDO crashes. While injury crashes increase by 0.46% for each additional one-percent increase in snowfall, PDO crashes increase by 0.52% for each additional one-percent increase in snowfall.  4.2.2 Interchange Crashes The Model 1 results for interchange crashes are provided in Table 8. Similar to the Model 1 midblock crash results, these findings illustrate that snowfall has an impact on and varies among the crash categories analyzed, albeit less emphasized in magnitude as opposed to midblock crashes. In all cases the attained equations along with the variables of length, AADT or commercial AADT, horizontal curvature, and snow are significant at the 0.01 significance level. Akin to the Model 1 midblock crash results, all variable coefficients are positive indicating increasing crash occurrence with increasing values of the explanatory variables.  Overall, winter interchange crashes increase by 0.4% for each one-percent increase in annual snowfall. Unlike the midblock crash results, crashes involving a truck or buss are less susceptible to snow with a 0.26% increase in crashes for each one-percent increase in annual snowfall. Comparatively, the effects of snow on non-truck/bus crashes increase by 0.4% for each one-percent increase in annual snowfall. Among crash severity outcomes, results are relatively similar to midblock crashes where injury crashes are less susceptible to snow as opposed to PDO  34 crashes, and experience a 0.37% increase for each one-percent increase in snowfall; while PDO crashes experience a 0.41% increase for each one-percent increase in snowfall.  TABLE 8  Model 1 (Natural Log of Snowfall) - Interchange Crashes NB Regression Results Variable B Std Error Lower CI Upper CI Wald Chi-Square df p-value All Winter Crashes Intercept -11.984 0.125 -12.229 -11.738 9165.422 1 <0.001 Ln Length 0.388 0.010 0.369 0.407 1621.871 1 <0.001 Ln AADT 1.124 0.010 1.105 1.143 13379.684 1 <0.001 Horiz. Curve < 85 mph 0.481 0.028 0.427 0.535 305.285 1 <0.001 Ln Snow 0.396 0.014 0.368 0.424 759.827 1 <0.001 Negative Binomial 0.631 0.011 0.610 0.652 na na na Non-Truck/Bus Crashes Intercept -12.050 0.128 -12.301 -11.798 8822.444 1 <0.001 Ln Length 0.388 0.010 0.369 0.408 1552.878 1 <0.001 Ln AADT 1.123 0.010 1.103 1.142 12749.754 1 <0.001 Horiz. Curve < 85 mph 0.517 0.028 0.462 0.572 340.363 1 <0.001 Ln Snow 0.393 0.015 0.364 0.422 708.554 1 <0.001 Negative Binomial 0.640 0.011 0.619 0.663 na na na Truck/Bus Crashes Intercept -12.076 0.243 -12.553 -11.600 2464.885 1 <0.001 Ln Length 0.220 0.021 0.179 0.261 110.795 1 <0.001 Ln Commercial AADT 1.224 0.025 1.175 1.272 2444.429 1 <0.001 Horiz. Curve < 85 mph 0.485 0.061 0.366 0.604 63.473 1 <0.001 Ln Snow 0.255 0.031 0.195 0.315 69.650 1 <0.001 Negative Binomial 0.866 0.052 0.770 0.974 na na na Injury Crashes Intercept -15.402 0.216 -15.824 -14.979 5096.870 1 <0.001 Ln Length 0.404 0.016 0.373 0.434 667.505 1 <0.001 Ln AADT 1.305 0.017 1.272 1.338 6011.224 1 <0.001 Horiz. Curve < 85 mph 0.499 0.040 0.421 0.577 156.899 1 <0.001 Ln Snow 0.368 0.023 0.322 0.413 252.326 1 <0.001 Negative Binomial 0.587 0.024 0.543 0.636 na na na PDO Crashes Intercept -11.935 0.130 -12.190 -11.680 8430.201 1 <0.001 Ln Length 0.385 0.010 0.366 0.405 1499.181 1 <0.001 Ln AADT 1.096 0.010 1.076 1.115 11805.188 1 <0.001 Horiz. Curve < 85 mph 0.440 0.028 0.384 0.495 241.655 1 <0.001 Ln Snow 0.408 0.015 0.378 0.437 748.062 1 <0.001 Negative Binomial 0.609 0.011 0.587 0.632 na na na na = not applicable; CI = Confidence Interval;  df = degrees of freedom                   35 4.2.3 Model 1 Results Discussion The Model 1 results on the effects of snowfall on all of the crash categories are summarized in terms of their elasticity in Table 9.  TABLE 9  Model 1 Snowfall Elasticities by Crash Category Crash Category Midblock  Interchange Percent Increase in Annual Crashes Associated with 1-pct. Increase in Annual Snowfall p-value  Percent Increase in Annual Crashes Associated with 1-pct. Increase in Annual Snowfall p-value All Crashes 0.501 <0.001  0.396 <0.001 Non-Truck/Bus 0.498 <0.001  0.393 <0.001 Truck/Bus 0.596 <0.001  0.255 <0.001 Injury 0.463 <0.001  0.368 <0.001 PDO 0.522 <0.001  0.408 <0.001  These results, in particular for midblock crashes, are generally in line with prior studies on the effects of snowfall on crashes (Qiu and Nixon, 2008), and are particularly intuitive for crashes involving trucks or buses. For instance, non-commercial trips are flexible in departure times, route changes, or route cancelation during severe snow events; while trips taken by truck or buses are generally business oriented thus less flexible on time and route choice. Truck performance is also more susceptible to decreasing visibility during adverse snow weather due to greater stopping distances required for larger and heavier vehicles (Pisano et al., 2008).  Existing literature also supports results for injury and PDO crashes with snowfall having a larger influence on PDO crashes and less, albeit positive, on injury crashes (Eisenberg, 2004; Eisenberg and Warner, 2005; Bilionis, 2013; Liu, 2013). These effects can be attributed to heightened driver alertness and lower speeds exhibited during snow events thus reducing the likelihood of high speed collisions and resulting injuries. Snow walls developed along the roadway can also be a cause for lower injury crash frequencies, and would be mostly applicable  36 on segments experiencing the highest amounts of snowfall supplemented by continuous sub-freezing temperatures.  In comparison, the effects of snowfall on interchange crashes, while not as emphasized in magnitude as midblock crashes, are not particularly surprising and appear to be in concert with Kiattikomol (2005) findings that interchange and non-interchange segments should be treated separately in crash modeling. While snowfall has a positive significant impact on all crash categories analyzed, its magnitude is likely suppressed due to the various effects in play in an interchange setting such as over-emphasized roadway geometry, sudden deceleration and acceleration of vehicles, and additional traffic flow patterns and driver behavior characteristics unique to these locations. These factors are not captured in this model but are potentially reflected in the lower magnitude and variability of the snowfall coefficient.  4.3 Model 2  Categorized Snowfall 4.3.1 Midblock Crashes The Model 2 results for midblock crashes are provided in Table 10. These results indicate that the impact of snowfall on winter crashes varies among the different snowfall intervals. In all cases the attained equations along with the variable of length, AADT or commercial AADT, horizontal curvature, and the snowfall intervals in the 3rd and 4th quartile ranges are significant at the 0.01 significance level.  TABLE 10  Model 2 (Categorized Snowfall) - Midblock Crashes NB Regression Results Variable B Std Error Lower CI Upper CI Wald Chi-Square df p-value Exp (B) All Winter Crashes Intercept -6.671 0.117 -6.900 -6.442 3251.847 1 <0.001 na Ln Length 1.215 0.012 1.191 1.239 9997.124 1 <0.001 na Ln AADT 0.655 0.011 0.633 0.678 3363.958 1 <0.001 na Horiz. Curve < 85 mph 0.499 0.037 0.427 0.572 181.420 1 <0.001 1.648  37          Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.171 0.019 0.134 0.208 80.850 1 <0.001 1.186 Snow - 4th Quartile 0.527 0.020 0.489 0.566 715.332 1 <0.001 1.694 Negative Binomial 0.797 0.015 0.768 0.827 na na na na Non-Truck/Bus Crashes Intercept -6.698 0.119 -6.931 -6.465 3183.369 1 <0.001 na Ln Length 1.218 0.012 1.193 1.242 9614.877 1 <0.001 na Ln AADT 0.649 0.011 0.626 0.671 3207.715 1 <0.001 na Horiz. Curve < 85 mph 0.534 0.038 0.460 0.608 200.237 1 <0.001 1.705 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.167 0.019 0.129 0.205 74.951 1 <0.001 1.182 Snow - 4th Quartile 0.514 0.020 0.475 0.554 657.344 1 <0.001 1.672 Negative Binomial 0.789 0.015 0.759 0.820 na na na na Truck/Bus Crashes Intercept -11.13 0.220 -11.563 -10.700 2555.502 1 <0.001 na Ln Length 1.225 0.028 1.171 1.279 1966.418 1 <0.001 na Ln Commercial AADT 1.098 0.026 1.046 1.150 1721.704 1 <0.001 na Horiz. Curve < 85 mph 0.553 0.093 0.372 0.735 35.676 1 <0.001 1.739 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.264 0.042 0.181 0.348 38.701 1 <0.001 1.303 Snow - 4th Quartile 0.618 0.041 0.538 0.698 228.126 1 <0.001 1.856 Negative Binomial 0.898 0.062 0.784 1.028 na na na na Injury Crashes Intercept -9.307 0.197 -9.693 -8.920 2225.145 1 <0.001 na Ln Length 1.249 0.021 1.208 1.289 3638.146 1 <0.001 na Ln AADT 0.754 0.019 0.717 0.790 1613.762 1 <0.001 na Horiz. Curve < 85 mph 0.581 0.061 0.462 0.701 91.000 1 <0.001 1.789 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.165 0.030 0.105 0.224 29.526 1 <0.001 1.179 Snow - 4th Quartile 0.481 0.031 0.420 0.543 233.479 1 <0.001 1.618 Negative Binomial 0.805 0.037 0.736 0.881 na na na na PDO Crashes Intercept -6.761 0.122 -7.001 -6.521 3057.707 1 <0.001 na Ln Length 1.218 0.013 1.193 1.243 9105.181 1 <0.001 na Ln AADT 0.643 0.012 0.619 0.666 2970.205 1 <0.001 na Horiz. Curve < 85 mph 0.469 0.039 0.393 0.546 145.287 1 <0.001 1.599 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.178 0.020 0.139 0.216 80.828 1 <0.001 1.194 Snow - 4th Quartile 0.540 0.020 0.500 0.580 702.204 1 <0.001 1.717 Negative Binomial 0.772 0.016 0.741 0.804 na na na na na = not applicable; CI = Confidence Interval; df = degrees of freedom            38 Overall, snowfall effects show monotonic increases when moving to higher quartiles or segments experiencing increasing annual snowfall amounts. And while the incremental differences in the regression coefficients between each of the snowfall quartiles are relatively consistent across all crash categories, annual snowfall amounts in the 4th quartile are associated with the greatest percent increase in crashes relative to the reference point (1st and 2nd quartile). Similar to the Model 1 midblock crash results, crashes involving a truck or bus are the most susceptible to increasing snowfall amounts. Translated in practical terms, while the 3rd snowfall quartile displays 30% greater crashes compared to the reference point, the 4th quartile range displays 86% greater crashes compared to the reference point. Comparatively, non-truck/bus crashes display approximately 18% and 67% greater crashes for the 3rd and 4th quartile intervals relative to the reference point. Similar effects are observed for injury and PDO crashes, where akin to Model 1 results, PDO crashes are more susceptible to higher snowfall amounts as opposed to injury crashes.  4.3.2 Interchange Crashes Lastly, the Model 2 results for interchange crashes are presented in Table 11. Similar to the Model 2 midblock crash findings, these results indicate that the effects of snowfall on winter crashes vary among the different snowfall intervals for this particular crash scenario. Furthermore, all attained equations along with all of the explanatory variables are statistically significant at the 0.01 significance level.  TABLE 11  Model 2 (Categorized Snowfall) - Interchange Crashes NB Regression Results Variable B Std Error Lower CI Upper CI Wald Chi-Square df p-value Exp (B) All Winter Crashes Intercept -10.47 0.101 -10.672 -10.275 10672.373 1 <0.001 na Ln Length 0.389 0.010 0.370 0.408 1619.511 1 <0.001 na  39          Ln AADT 1.103 0.010 1.084 1.122 12779.120 1 <0.001 na Horiz. Curve < 85 mph 0.482 0.028 0.428 0.537 302.460 1 <0.001 1.620 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.184 0.016 0.153 0.215 134.706 1 <0.001 1.202 Snow - 4th Quartile 0.385 0.018 0.350 0.419 479.367 1 <0.001 1.469 Negative Binomial 0.644 0.011 0.623 0.666 na na na na Non-Truck/Bus Crashes Intercept -10.54 0.104 -10.741 -10.334 10305.667 1 <0.001 na Ln Length 0.390 0.010 0.370 0.409 1549.626 1 <0.001 na Ln AADT 1.101 0.010 1.081 1.120 12158.944 1 <0.001 na Horiz. Curve < 85 mph 0.518 0.028 0.463 0.574 336.442 1 <0.001 1.679 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.181 0.016 0.149 0.213 124.614 1 <0.001 1.198 Snow - 4th Quartile 0.373 0.018 0.338 0.409 429.839 1 <0.001 1.453 Negative Binomial 0.655 0.011 0.633 0.678 na na na na Truck/Bus Crashes Intercept -11.11 0.201 -11.502 -10.713 3039.152 1 <0.001 na Ln Length 0.227 0.021 0.186 0.268 117.659 1 <0.001 na Ln Commercial AADT 1.204 0.025 1.156 1.252 2393.418 1 <0.001 na Horiz. Curve < 85 mph 0.467 0.061 0.348 0.587 58.615 1 <0.001 1.596 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.227 0.034 0.160 0.295 43.830 1 <0.001 1.255 Snow - 4th Quartile 0.214 0.038 0.140 0.288 31.903 1 <0.001 1.239 Negative Binomial 0.874 0.052 0.778 0.982 na na na na Injury Crashes Intercept -14.06 0.179 -14.410 -13.709 6168.780 1 <0.001 na Ln Length 0.406 0.016 0.375 0.437 673.368 1 <0.001 na Ln AADT 1.291 0.017 1.257 1.324 5819.473 1 <0.001 na Horiz. Curve < 85 mph 0.499 0.040 0.421 0.577 156.377 1 <0.001 1.647 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.160 0.024 0.113 0.206 44.412 1 <0.001 1.173 Snow - 4th Quartile 0.390 0.029 0.333 0.446 184.053 1 <0.001 1.476 Negative Binomial 0.593 0.024 0.548 0.642 na na na na PDO Crashes Intercept -10.38 0.105 -10.587 -10.174 9700.057 1 <0.001 na Ln Length 0.387 0.010 0.367 0.407 1495.482 1 <0.001 na Ln AADT 1.074 0.010 1.054 1.094 11236.257 1 <0.001 na Horiz. Curve < 85 mph 0.441 0.028 0.385 0.497 239.595 1 <0.001 1.554 Snow - 1st & 2nd Quartile base na na na na na na na Snow - 3rd Quartile 0.188 0.016 0.156 0.220 132.440 1 <0.001 1.207 Snow - 4th Quartile 0.393 0.018 0.358 0.429 470.284 1 <0.001 1.482 Negative Binomial 0.624 0.012 0.602 0.648 na na na na na = not applicable; CI = Confidence Interval; df = degrees of freedom                    40 Overall, the Model 2 results for interchange crashes reflect combined findings for the Model 1 interchange and Model 2 midblock crashes. Akin to Model 1 interchange crash results, the snowfall coefficients in this model are lower in magnitude for each snowfall category. Analogous to Model 2 midblock crashes, effects generally display monotonic increases when moving to higher quartiles, with the incremental differences between each of the snowfall categories being reasonably consistent across the various types of crashes considered in the model. With the exception of crashes involving a truck or bus, which show similar crash increases for the 3rd and 4th quartile range, annual snowfall amounts in the 4th quartile are generally associated with the greatest percent increase in crashes relative to the reference point. 4.3.3 Model 2 Results Discussion The Model 2 results on the effects of snowfall on all of the crash categories are summarized in terms of their elasticity in Table 12. TABLE 12  Model 2 Snowfall Elasticities by Crash Category Crash Category Snowfall Interval  in Inches Midblock  Interchange Percent Increase in Annual Crashes Associated with Increase in Annual Snowfall Over Base Category p-value  Percent Increase in Annual Crashes Associated with Increase in Annual Snowfall Over Base Category p-value All Crashes 37.4 - 50.2 18.6 <0.001  20.2 <0.001  50.2 - 157.1 69.4 <0.001  46.9 <0.001 Non-Truck/Bus 37.4 - 50.2 18.2 <0.001  19.8 <0.001  50.2 - 157.1 67.2 <0.001  45.3 <0.001 Truck/Bus 37.4 - 50.2  30.3 <0.001  25.5 <0.001  50.2 - 157.1 85.6 <0.001  23.9 <0.001 Injury 37.4 - 50.2 17.9 <0.001  17.3 <0.001  50.2 - 157.1 61.8 <0.001  47.6 <0.001 PDO 37.4 - 50.2 19.4 <0.001  20.7 <0.001  50.2 - 157.1 71.7 <0.001  48.2 <0.001 Base Category for Annual Snowfall = 5.2 inches  37.4 inches (Quartiles 1 and 2) Quartile 3 = 37.4 inches  50.2 inches; Quartile 4 = 50.2 inches  157.1 inches   41 Overall, these findings are surprising as prior studies have suggested that the effects of snowfall on crashes peak at mid ranges and decline or have a lower marginal effect at the higher snowfall levels (Khattak and Knapp, 2001; Eisenberg, 2004; Seeherman and Liu, 2015). One possible explanation for these differences is that the annual snowfall data represents an aggregation of weather patterns and is displaying the effects of underlying events such as higher frequencies of snowstorms and/or storms of greater intensity. This would certainly be plausible given that in nearly all cases the effects of snowfall are highest for those segments experiencing the largest amount of snow.  4.4 Midblock Snowfall Crash Curves In order to present a secondary avenue for the examination of snowfall effects on crashes and introduce practical implications, Model 1 midblock Negative Binomial crash findings are translated in graphical form for each of the crash types analyzed. To establish the graphical relationship, crashes are predicted as a function of AADT or commercial AADT, length, horizontal curvature, and snowfall, all of which are found significant  in the statistical models at the 0.01 significance level. For each graph, predicted crashes are presented as crashes per mile by fixing length to 1 (Y-axis) and plotting it against the AADT or commercial AADT network data range (X-axis), where similar to the Negative Binomial models commercial AADT is only applicable for crashes involving a truck or bus. Crash curves are established for six snowfall 20 inch intervals reflecting values in the network annual snowfall data range. Two graphs are produced for each crash type, one for tangent segments (i.e. fraction of segment with horizontal curve with design speed < 85 mph = 0) and one for curved segments using the mean value of the horizontal curvature network data as baseline (i.e. fraction of segment with horizontal curve with design speed < 85 mph = 0.1). Additional baseline values applicable to these graphical  42 relationships include: Speed Limit = 70 mph / 60 mph, Average Number of Lanes = 2.5, Average Lane Width = 12 ft, Average Shoulder Width = 9.3 ft.  Lastly, since interchange crashes are subject to additional factors unique to their design which are not fully captured by the presented models, they are not given the same graphical treatment as midblock crashes.  Figure 7 presents an overlay of selected crash curves for all winter midblock crashes for both tangent and curved segments to allow for the direct comparison of the two curvature scenarios, while Figure 8 through 12 illustrate the developed crash curves for each midblock crash type analyzed.   FIGURE 7  Effect of Annual Snowfall on Winter Midblock Crashes Combined for Tangent and Curved Segments  43   FIGURE 8  Effect of Annual Snowfall on Winter Midblock Crashes for Tangent and Curved Segments  44   FIGURE 9  Effect of Annual Snowfall on Non-Truck/Bus Winter Midblock Crashes for Tangent and Curved Segments  45   FIGURE 10  Effect of Annual Snowfall on Truck/Bus Winter Midblock Crashes for Tangent and Curved Segments  46   FIGURE 11  Effect of Annual Snowfall on Injury Winter Midblock Crashes for Tangent and Curved Segments  47   FIGURE 12  Effect of Annual Snowfall on PDO Winter Midblock Crashes for Tangent and Curved Segments  48 Similar to the negative binomial results, predicted crashes for each crash type increase with increasing values of AADT or commercial AADT, and snowfall. Crash occurrence is furthermore larger in the presence of horizontal curvature. These results are overall consistent for each crash category, with the crash curves for each snowfall interval illustrating similar graphical patterns among them. Crashes involving a truck or bus represent the only exception to the latter where commercial AADT has a larger impact on these types of crashes with a parameter coefficient greater than one, thus leading to an exponential growth pattern for the provided crash curves.                 49 CHAPTER 5: SPATIAL PATTERN ANALYSIS  5.1 Crashes & Kernel Density Function While statistical analysis demonstrates that annual snowfall has a statistically significant positive effect on all crash types and locational scenarios analyzed, an examination of the spatial dispersion pattern of these crashes can provide practical insights which may assist agencies in identifying potential candidate segment areas for countermeasure implementation. Spatial patterns for freeway crashes in this study are examined using the kernel density function. This method is used at length  in crash analysis (Xie and Yan, 2008; Anderson, 2009; Blazquez and Celis, 2012).  Kernel density is a spatial pattern identifier which calculates the density of data points in a circular region. The circular area over which density is determined is called the kernel (Rushton and Tiwari, 2009). Within the kernel, density peaks at the location of the data point and becomes zero at the kernel boundary. When kernel density is applied over counts (i.e. individual crashes in XY spaces), the total sum value under the kernel is 1. When counts are assigned a secondary value (i.e. snowfall amount), the kernel assumes this value within its perimeter (How Kernel Density works, 2016). Additionally, kernels are drawn for each cell in the geographic extent of the dataset thus resulting in a smooth raster output also known as spatial filtering which helps to identify potential patterns (Rushton and Tiwari, 2009; How Kernel Density works, 2016). This process can be utilized for both points (i.e. crashes) and lines (i.e. segments) in GIS space. A simplified visual representation of the kernel density process is presented in Figure 13.    50  FIGURE 13  Simplified Visual Representation of Kernel Density Process Kernel density on crash analysis is typically applied to the count distribution of crashes to identify potential patterns. This process however is not applicable for this particular study since crashes along the subject freeway network are overly concentrated in the Metro Detroit region due to elevated AADT volumes and a relatively high concentration of interchanges, thus producing unusable patterns. To overcome this challenge, crashes are applied to the freeway segments in crash rate form for midblock scenarios. Crash rate in this case is defined as:   (10) Where,  = 2004-2014 winter midblock crash rate for segment i  = 2004-2014 winter midblock crashes for segment i  = 2004-2014 average AADT for segment i  = number of years of data  = length in miles for segment i The applied transformation eliminates the problem of high crash densities in high AADT regions since crashes are presented as a function of AADT and length. These two variables are also included as explanatory variables in the Negative Binomial models and are found significant  51 at the 0.01 significance level for all crash types analyzed. The additional omission of interchange crashes further facilitates the pattern identification process since only midblock crashes are considered. These are correspondingly shown to be more susceptible to snowfall effects for all crash categories. Not surprisingly the relationship between average annual snowfall and winter midblock crash rates for the network has a linear correlation with a Pearson coefficient of 0.327 (p-value <0.001) (Figure 14). The linear relationship indicates that spatial identification of segments most prone to snowfall is feasible under this structure.  FIGURE 14  2004-2014 Network Winter Midblock Crash Rate & Average Annual Snowfall Following the crash rate transformation, the kernel density function is performed on all winter midblock crash rates for two case studies, Interstate 75 (I-75) and Interstate 94 (I-94). These two freeway corridors are selected based on their lateral and longitudinal spread over the state of Michigan, thus providing a relatively complete coverage of snowfall distribution for both north-south and east-eest directions.  To provide comparable metrics, the kernel density function is also applied to the AADT volumes and annual snowfall amounts averaged for the 11 year  52 period of 2004-2014 for each freeway segment corresponding to the two corridors. Kernel outputs are then stacked against each other to facilitate the identification of crash rate patterns. The two case studies are presented below.  5.2 Case Study 1  Interstate 75 Case study 1 is represented by I-75. This freeway corridor runs from the Ohio border to the south up to the Canadian border to the north. It is the longest freeway in Michigan and covers a distance of approximately 397 miles per direction. The lateral spread provides an ideal distribution in assessing crash patterns between the northern and southern parts of the state. Similar to the network, I-75 exhibits a linear relationship between average annual snowfall and winter midblock crash rates with a Pearson coefficient of 0.423 (p-value <0.001) (Figure 15). The spatial distribution of this relationship is presented below in Figure 16. This image compares in sequential order the AADT, snowfall, and winter midblock crash rate kernel density distribution along the entire interstate.  FIGURE 15  2004-2014 I-75 Winter Midblock Crash Rate & Average Annual Snowfall  53  FIGURE 16  I-75 Kernel Density Distribution The kernel density output indicates very clear and noticeable crash rate patterns for I-75. While the segments with the highest AADT volumes are concentrated in the Metro Detroit area,  54 average annual snowfalls and winter midblock crash rates for this region are relatively low. Conversely the northern section of I-75 between the city of Gaylord and the Canadian border has the lowest AADT volumes, highest average annual snowfall, and highest winter midblock crash rate concentration. The density for the three variables peaks in these corresponding directions for the freeway segments located between the Mackinac Bridge and the Canadian border. This particular distribution indicates that crash rates in this region are more susceptible to snowfall, potentially due to a higher frequency of snowstorms and/or storms of greater intensity experienced during winter periods.  Accordingly countermeasures installed in these segments can have the most impact on winter crash rates. Examples of potential countermeasures may include environmental sensor stations installed near the most susceptible segments to improve winter maintenance, or linking environmental sensors to ITS devices like flashing weather alert signs or variable message signs to alert drivers of adverse weather conditions when passing specific locations along the corridor. On a macro level scale, the attained crash rate patterns can be used to assist agencies in countermeasure implementation segment prioritization.  5.3 Case Study 2  Interstate 94 While case study 1 presents a lateral cross-section of the state, case study 2 or I-94 presents a longitudinal one. This freeway corridor runs from the Canadian border to the east down to the Indiana border to the west for approximately 276 miles in length per direction. Similar to the I-75 corridor, I-94 provides an ideal longitudinal distribution in assessing crash patterns between the eastern and western portions of the state. Comparably to both the network and the I-75 corridor, I-94 exhibits a linear relationship between average annual snowfall and winter midblock crash rates with a Pearson coefficients of 0.397 (p-value <0.001) (Figure 17).  55 The spatial distribution of this relationship is presented below in Figure 18. This image lists in sequential order the AADT, snowfall, and winter midblock crash rate kernel density distribution along the entire interstate.   FIGURE 17  2004-2014 I-94 Winter Midblock Crash Rate & Average Annual Snowfall           56  FIGURE 18  I-94 Kernel Density Distribution Similar to the I-75 case study, the kernel density function for the I-94 corridor presents clear and interpretable patterns. Comparable to I-75, the Metro Detroit region is characterized by the concentration of the highest AADT volumes, and relatively lower average annual snowfall amounts and winter midblock crash rates. Whereas the western section of the corridor in proximity of the coastline (i.e. Van Buren and Berrien County) displays relatively lower AADT  57 volumes, and the highest average annual snowfall amounts and winter midblock crash rate concentration for the corridor. Akin to the northern parts of I-75, this freeway section can benefit from winter weather related countermeasures such as winter weather maintenance sensors and/or ITS signs to alert drivers of adverse winter weather conditions.                      58 CHAPTER 6: CONCLUSIONS  This study investigated the effects of annual snowfall on freeway crashes during winter periods in the State of Michigan. Two principal Negative Binomial regression models are established to assess these effects. Explanatory variables for Model 1 include AADT or commercial AADT, segment length, horizontal curvature, and snowfall. Explanatory variables for Model 2 include all four initial variables from Model 1 and three categorical snowfall intervals based on its quartile distribution. Crash categories considered in the analysis include total winter, truck/bus, non-truck/bus, fatal and injury, and PDO crashes. Separate analyses are additionally performed for each model for interchange versus non-interchange segment crashes to account for their unique physical, operational and behavioral characteristics. This inclusion of the listed crash categories as well as their differentiation between interchange and non-interchange crashes provides an additional perspective not previously covered in the existing literature with regards to the impacts of snowfall on crashes.  The Negative Binomial regression results indicate that annual snowfall has a statistically significant positive correlation with winter crashes for all of the crash categories analyzed. The snowfall effects are stronger for non-interchange segment crashes (i.e. midblock) as opposed to interchange segment crashes. Among crash categories, crashes involving a truck or bus experience the most elastic relationship with snowfall for midblock scenarios. Among crash outcomes, PDO crashes are shown to be more strongly affected by snow conditions as opposed to injury and fatal crashes. The effects of snowfall on crashes are further exacerbated along those segments experiencing the largest amount of annual snowfall (i.e. top 25th percentile of annual snowfall totals). The models further indicate that the risk of crash occurrence increases with  59 increasing AADT or commercial AADT values and along segments characterized by the presence of horizontal curvature. Overall these findings suggest that drivers are at elevated risks of crash occurrence on freeways experiencing greater annual snowfall totals. The risk is dependent on vehicle type, and is highest for PDO crashes as opposed to injury or fatal crashes. These outcomes can be attributed to differing trip and mechanistic characteristics between commercial and non-commercial vehicles, and heightened driver alertness and lower speeds exhibited during severe snow events (Eisenberg, 2004; Pisano et al., 2008; Bilionis, 2013). Although the risk of crash occurrence is exacerbated along non-interchange segment areas, interchange crashes continue to occur at a higher rate thus indicating that there are other factors at play in such settings not captured by the models. Similarly, while crash risks are highest for the segments experiencing the greatest amounts of snowfall, the findings are inconsistent with prior studies (Khattak and Knapp, 2001; Eisenberg, 2004; Seeherman and Liu, 2015) which have suggested that the effects of snowfall on crashes peak at mid ranges and decline or have a lower marginal effect at the higher snowfall levels, thus likely reflecting underlying events such as higher frequencies of snowstorms and/or storms of greater intensity.  Nonetheless, the areas which experience the highest midblock winter crash rates correspond to those segments experiencing the highest amount of snowfall totals. For both case studies investigated (i. e. I-75 and I-94), the implemented kernel density function displays robust patterns in its identification of the most susceptible freeway segments. These patterns, which can be replicate on other freeways with similar geographic and geometric characteristics, provide practical insights which can assist agencies in countermeasure implementation. Examples of potential countermeasures include environmental sensor stations installed near the most  60 susceptible segments to improve winter maintenance, and linking such sensors to ITS devices to alert drivers of adverse weather conditions when passing particular locations along a freeway corridor.  In conclusion, while these findings are applicable to the freeway network in the State of Michigan, the methodology can be replicated and applied on other regions experiencing significant snowfall. These methods are particularly applicable for annual, monthly, and weekly weather data, as long as such data are available from weather stations located in proximity to the highway in question. Caution should be used when applying these methods to daily, and/or hourly weather as the number of weather stations reporting information at this level of detail is generally low, thus leading to inaccurate results. Likewise, while the Ordinary Kriging method performs satisfactory on regions with a relatively flat terrain (i.e. Michigan), elevation data should be incorporated into the spatial interpolation procedures when predicting weather data on regions characterized by a relatively rugged and in particular mountainous terrain. The inclusion of additional explanatory variables such as vertical curvature and other freeway geometric and operational characteristics could further strengthen the regression models.  Lastly, analysis on a micro-scale (i. e. daily snowfall) could help provide additional detail regarding the winter weather crash causal factors. To minimize data bias due to a limited weather station sample, only corridors surrounded by a sufficient number of weather stations should be considered. Potential variables for such microscopic forms of analysis can include an examination of daily snowfall totals in combination with the number of days with snowfall occurring along a segment. These two variables when combined together reflect more closely the frequency and intensity of snowstorms occurring on a particular roadway section. The additional inclusion of weather related factors such as average temperature and wind speeds on those days  61 experiencing snowfall could help further explain some of the conditions that drivers may be experiencing on the road. Additional non-weather related variables that can strengthen these microscopic models can include applicable geometric and operational roadway characteristics, as well roadside terrain characteristics which represent the natural and man-mad barriers located along the subject segments. The latter may be pertinent in explaining scenarios where snowfall   62 REFERENCES  63 REFERENCES  Ahmed, M. M., H. Huang, M. Abdel-Aty, and B. Guevara. Exploring a Bayesian hierarchical approach for developing safety performance functions for a mountainous freeway. Accident Analysis and Prevention, Vol. 43, 2011, pp. 1581-1589. Ahmed, M. M., M. Abdel-Aty, and R. Yu. Assessment of Interaction of Crash Occurrence, Mountainous Freeway Geometry, Real-Time Weather, and Traffic Data. Transportation Research Record, No. 2280, 2012, pp. 51-59. Ali, H. Z. Cross-Validation of Elevation Data Interpolation. Journal of Al Rafidain University College, No. 32, 2013, pp. 113-123. American Association of State Highway and Transportation Officials (AASHTO). Highway Safety Manual. AASHTO, Washington D.C., 2010. Aderson, T. K. Kernel density estimation and K-means clustering to profile road accident hotspots. Accident Analysis & Prevention, Vol. 41, No. 3, 2009, pp. 359-364. Andersson, A. K., and L. Chapman. The impact of climate change on winter road maintenance and traffic accidents in West Midlands, UK. Accident Analysis & Prevention, Vol. 43, No. 1, 2011, pp. 284-289. Andresen, J., S. Hilberg, and K. Kunkel. Historical Climate and Climate Trends in the Midwestern USA. Prepared for the US National Climate Assessment Midwest Technical Report, Great Lakes Integrated Sciences and Assessments Center, Ann Arbor, MI, 2012. Andrescu, M., and D. Frost. Weather and traffic accidents in Montreal, Canada. Climate Research, Vol. 9, 1998, pp. 225-230. ArcGIS 9: Using ArcGIS Geostatistical Analyst. Esri, Redland, CA, 2003. Ashraf, M., J. C. Lofits, and K. G. Hubbard. Application of geostatistics to evaluate partial weather station networks. Agricultural and Forest Meteorology, Vol. 84, No. 3, 1997, pp. 255-271.  Bilionis, D. Interaction effects of prevailing weather conditions and crash characteristics on crash severity: A case study of two corridors in Iowa. MS thesis, Iowa State University, Ames 2013. http://lib.dr.iastate.edu/cgi/viewcontent.cgi?article=4189&context=etd. Accessed Feb. 24, 2016. Blazquez, C. A., and M. S. Celis. A spatial and temporal analysis of child pedestrian crashes in Santiago, Chile. Accident Analysis and Prevention, Vol. 50, 2013, pp. 304-311. Brorsson, B., J. Ifver, and H. Rydgren. Injuries from single-vehicle crashes and snow depth. Accident Analysis & Prevention, Vol. 20, No. 5, 1988, pp. 367-377.  64 Bostan, P. A., G. B. M. Heuvelink, and S. Z. Akyurek. Comparison of regression and kriging techniques for mapping the average annual precipitation of Turkey. International Journal of Applied Earth Observation and Geoinformation, Vol. 19, 2012, pp. 115-126. Eisenberg, D. The mixed effects of precipitation on traffic crashes. Accident Analysis and Prevention, Vol. 36, No. 4, 2004, pp. 637-647. Eisenberg, D. and K. E. Warner. Effects of snowfalls on motor vehicle collisions, injuries, and fatalities. American Journal of Public Health, Vol 95, 2005, pp. 120-214. El-Basyouny, K., S. Barua, and M. T. Islam. Investigation of time and weather effects on crash types using full Bayesian multivariate Poisson lognormal models. Accident Analysis and Prevention, Vol. 73, 2004, pp. 91-99. Fridstrom, L., J. Liver, S. Ingebrigsten, R. Kulmala, and L. Thomsen. Measuring the contribution of randomness, exposure, weather, and daylight to the variation in road accident counts. Accident Analysis and Prevention, Vol. 27, No. 1, 1995, pp. 1-20. Goovaerts, P. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. Journal of Hydrology, Vol. 228, No. 1, 2000, pp. 113-129.  Holdaway, M. R. Spatial modeling and interpolation of monthly temperature using kriging. Climate Research, Vol. 6, No. 3, 1996, pp. 215-225. How Kernel Density works. Esri, Redlands, CA, 2016. http://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/how-kernel-density-works.htm. Accessed Oct. 20, 2016. Huang, C. L., H. W. Wang, and J. L. Hou. Estimating spatial distribution of daily snow depth with kriging methods: combination of MODIS snow cover area data and ground-based observations. The Cryosphere Discussions, Vol. 9, No. 5, 2015, pp. 4997-5020.  Khan, G., X. Qin, and D. A. Noyce. Spatial Analysis of Weather Crash Patterns. Journal of Transportation Engineering, Vol. 134, No. 5, 2008, pp. 191-202. Khattak, A., and K. Knapp. Snow Event Effects on Interstate Highway Crashes. Journal of Cold Regions Engineering, Vol. 15, No. 4, 2001, pp. 219-229.  Kiattikomol, V. Freeway Crash Prediction Models for Long-Range Urban Transportation Planning. PhD dissertation, University of Tennessee, Knoxville, 2005. http://trace.tennessee.edu/cgi/viewcontent.cgi?article=3602&context=utk_graddiss. Accessed Oct. 5, 2016. Kim, T. Y., K. W. Kim, and B. H. Park. Accident Models of Trumpet Interchange S-type Ramps Using by Poisson, Negative Binomial Regression and ZAM. KSCE Journal of Civil Engineering, Vol. 15, No. 3, 2010, pp. 545-551. Laaha, G., J. O. Skoien, F. Nobilis, and G. Bloschl. Spatial Prediction of Stream Temperatures Using Top-Kriging with an External Drift. Environmental Modeling & Assessment, Vol.  65 18, No. 6, 2013, pp. 671-683.  Laaha, G., J. O. Skoien, and G. Bloschl. Spatial prediction on river networks: comparison of top-kriging with regional regression. Hydrological Processes, Vol. 28, No. 2, 2014, pp. 315-324. Lanciani, A., and M. Salvati. Spatial interpolation of surface weather observations in alpine meteorological services. FORALPS Technical Report 2, Universita degli Studi di Trento, Trento, Italy, 2008. Leard, B., and B. Roth. Weather, Traffic Accidents and Climate Change. SSRN Electronic Journal. Resources for the Future Discussion Paper, 2015, pp. 15-19.  Liu, Y. Weather Impact on Road Accident Severity in Maryland. MS thesis, University of Maryland, College Park  2013. http://drum.lib.umd.edu/bitstream/1903/14263/1/Liu_umd_0117N_14019.pdf. Accessed Feb.24, 2016. Lord, A. and F. Mannering. The statistical analysis of crash-frequency data: A review and assessment of methodological alternatives. Transportation Research Part A: Policy Practice, Vol. 44, No. 5, 2010, pp. 291-305. Mair, A. and A. Fares. Comparison of Rainfall Interpolation Methods in a Mountainous Region of a Tropical Island. Journal of Hydrological Engineering, Vol. 16, No. 4, 2011, pp. 371-383. Oliver, M. A., and R. Webster. Basic Steps in Geostatistics: The Variogram and Kriging. Springer, Cham, Switzerland, 2015.  Perry, A. H., and L. J. Symons. Highway Meteorology. E & FN Spon, Swansea, Wales, 1991. Pisano, P., L. Goodwin, and M. Rossetti. U.S. highway crashes in adverse road weather conditions. Presented at the 24th Conference on Information and Processing Systems during the 88th annual meeting of the American Meteorological Society, New Orleans, LA, 2008. Qiu, L., and W. A. Nixon. Effects of Adverse Weather on Traffic Crashes: Systematic Review and Meta-Analysis. Transportation Research Record, No. 2055, 2008, pp 139-146. Rushton, G., and C. Tiwari. Spatial Filtering/Kernel Density Estimation. International Encyclopedia of Human Geography, Vol. 10, 2009, pp. 359-364. Seeherman, J., and Y. Liu. Effects of extraordinary snowfall on traffic safety. Accident Analysis and Prevention, Vol. 81, 2015, pp. 194-203. Sell, S. One Year Later: 193 Vehicle Pile-up on I-94. Detroit Free Press, January 2016. http://www.freep.com/story/news/local/michigan/2016/01/08/i94-car-pileup-snow-storm/78531362/. Accessed Oct. 19, 2016.  66 Sluiter, R. Interpolation methods for climate data  literature review. KNMI intern rapport Royal Netherlands Meteorological Institute, De Bilt, Netherlands, 2009. Tabios, G. Q., and J. D. Salas, A comparative analysis of techniques for spatial interpolation of precipitation. Journal of the American Water Resources Association, Vol. 21, No. 3, 1985, pp. 365-380. U.S. DOT Federal Highway Administration (FHWA), How Do Weather Events Impact Roads? May 2016. http://www.ops.fhwa.dot.gov/weather/q1_roadimpact.htm, Accessed Jul. 28, 2016. Xie, Z., and J. Yan. Kernel Density Estimation of Traffic Accidents in a Network Space. Computers, Environment, and Urban Systems, Vol. 32, No. 5, 2008, pp. 396-406. Zou, Y., W. Lingtao, and D. Lord. Modeling over-dispersed crash data with long tail: Examining the accuracy of the dispersion parameter in Negative Binomial models. Analytic Methods in Accident Research, Vol. 5-6, 2015, pp. 1-16.