Applied Science, Faculty of
Engineering, School of (Okanagan)
DSpace
UBCO
Faghihi, Seyed Farhad
Doctor of Philosophy - PhD
University of British Columbia
Many sets of road traffic data are a collection of ordered measures over the course of time. The key characteristic of time series data is the autocorrelation inside the observed points. Time series analysis is a technique to find the patterns inside time series data, while considering the autocorrelation. The existing time series statistical models in the literature suffer from several assumptions such as stationarity, normality, or homoscedasticity. In this research, a two-stage methodology of Discrete Wavelet Transform (DWT) coupled with Dynamic Artificial Neural Network (DANN) was proposed for predicting a given collision time series. The DANN model enabled the mapping of nonlinear, nonstationary relationship between the future values of collision series and its previous values or, moreover, the previous values of other time series. The DWT facilitated the learning process for limited, noisy observations by decomposing a time series into multiple subseries, which represented specific features of data. The relative performance of the proposed methodology was compared to the conventional AutoRegressive Integrated Moving Average (ARIMA) model for road safety studies. The variables used to build and validate the models were monthly time series data recorded from January 2003 to June 2016 in the City of Kelowna, BC. 54,000 DANN models were trained in total to determine the optimum network configurations for the subseries. The findings of this research showed that the combined DWT and DANN methodology outperformed the ARIMA models in terms of several indices. In addition, the proposed methodology was applied to evaluate the safety phases implemented on Springfield Road, Kelowna at various times. The results demonstrated that the methodology was also capable of using for the intervention analysis, where the traditional time series modelling approaches cannot be executed. The comprehensive study of this thesis is expected to have important implications for both theory and practice in transportation engineering and road safety.
APPLICATION OF DYNAMIC ARTIFICIAL NEURAL NETWORK TO ROAD COLLISION PREDICTION MODELLING by Seyed Farhad Faghihi B.A.Sc., In Civil Engineering, Azad University, 2008 M.A.Sc., In Transportation Engineering, Azad University, 2012 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE COLLEGE OF GRADUATE STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan) September 2018 © Farhad Faghihi, 2018 The following individuals certify that they have read, and recommend to the College of Graduate Studies for acceptance, a thesis/dissertation entitled: APPLICATION OF DYNAMIC ARTIFICIAL NEURAL NETWORK TO ROAD COLLISION PREDICTION MODELLING submitted by Seyed Farhad Faghihi in partial fulfillment of the requirements of the degree of Doctor of Philosophy . Dr. Gordon Lovegrove, School of Engineering Supervisor Dr. Abbas Milani, School of Engineering Supervisory Committee Member Dr. Ahmed Idris, School of Engineering Supervisory Committee Member Dr. Mohammad H. Zarifi, School of Engineering University Examiner Dr. Lina Kattan, School of Engineering, University of Calgary External Examiner ii Abstract Many sets of road traffic data are a collection of ordered measures over the course of time. The key characteristic of time series data is the autocorrelation inside the observed points. Time series analysis is a technique to find the patterns inside time series data, while considering the autocorrelation. The existing time series statistical models in the literature suffer from several assumptions such as stationarity, normality, or homoscedasticity. In this research, a two-stage methodology of Discrete Wavelet Transform (DWT) coupled with Dynamic Artificial Neural Network (DANN) was proposed for predicting a given collision time series. The DANN model enabled the mapping of nonlinear, nonstationary relationship between the future values of collision series and its previous values or, moreover, the previous values of other time series. The DWT facilitated the learning process for limited, noisy observations by decomposing a time series into multiple subseries, which represented specific features of data. The relative performance of the proposed methodology was compared to the conventional AutoRegressive Integrated Moving Average (ARIMA) model for road safety studies. The variables used to build and validate the models were monthly time series data recorded from January 2003 to June 2016 in the City of Kelowna, BC. 54,000 DANN models were trained in total to determine the optimum network configurations for the subseries. The findings of this research showed that the combined DWT and DANN methodology outperformed the ARIMA models in terms of several indices. In addition, the proposed methodology was applied to evaluate the safety phases implemented on Springfield Road, Kelowna at various times. The results demonstrated that the methodology was also capable of using for the intervention analysis, where the traditional time series modelling approaches iii cannot be executed. The comprehensive study of this thesis is expected to have important implications for both theory and practice in transportation engineering and road safety. iv Lay Summary The number of accidents during equal periods, for example months, is a measure of road safety on a road or within a region. This research attempted to develop a method, which is inspired by human brain function that better predicts the monthly number of accidents in the City of Kelowna, BC. The same method was also used to assess how changes implemented on Springfield Road, Kelowna affected the number of accidents, and thus how the changes influenced the transport safety of this road. v Preface Portions of this thesis have been published as part of a poster presentation and conference proceeding respectively as follows:  Faghihi, F. & Lovegrove G. (2018). Application of Dynamic (Time Series) Artificial Neural Network Model to Develop Collision Prediction Models: Case Study in City of Kelowna, Canada. The 97th Transportation Research Board (TRB) Annual Meeting, January 7th to 11th.  Faghihi et al. (2016). Effectiveness of Applying ‘5E’ Approach Integrated with ITS and Speed Sensitive Traffic Signals to Support Complete Streets–Case Study Springfield Road Kelowna, BC. Canadian Institute of Transportation Engineer Conference. Kelowna, BC. vi Table of Contents Abstract .......................................................................................................................................... ii Lay Summary ............................................................................................................................... iv Preface .............................................................................................................................................v Table of Contents ......................................................................................................................... vi List of Tables ............................................................................................................................... xii List of Figures ............................................................................................................................ xvii List of Symbols .............................................................................................................................xx List of Abbreviations .................................................................................................................xxv Chapter 1: Introduction ................................................................................................................1 1.1 Overview ...................................................................................................................... 1 1.2 Road safety issues ........................................................................................................ 1 1.3 Problem statement ........................................................................................................ 2 1.4 Research objectives ...................................................................................................... 4 1.5 Thesis framework ......................................................................................................... 4 Chapter 2: Literature Review .......................................................................................................6 2.1 Overview ...................................................................................................................... 6 2.2 Road safety studies....................................................................................................... 6 2.3 Introduction to time series analysis .............................................................................. 7 2.4 Overview of time series features .................................................................................. 8 2.5 Time series statistical models ...................................................................................... 9 2.5.1 Preliminary time series studies in road safety ............................................................ 13 vii 2.5.2 Time domain models .................................................................................................. 14 2.5.2.1 Deterministic time series model........................................................................ 14 2.5.2.2 Autoregressive integrated moving average model ............................................ 16 2.5.2.3 State space model .............................................................................................. 17 2.5.2.3.1 Advantages and disadvantages of state space models ................................. 20 2.5.2.4 Vector autoregressive model ............................................................................. 21 2.5.3 Frequency domain models ......................................................................................... 23 2.5.4 Nonlinear time series models ..................................................................................... 24 2.6 History of artificial neural networks .......................................................................... 25 2.6.1 Artificial neural network in road safety analysis ....................................................... 28 2.7 Summary .................................................................................................................... 31 Chapter 3: Autoregressive Integrated Moving Average Modelling ........................................32 3.1 Overview .................................................................................................................... 32 3.2 Introduction to ARIMA model .................................................................................. 32 3.3 General form of ARIMA models ............................................................................... 33 3.4 Autocorrelation and partial autocorrelation function ................................................. 37 3.5 Summary .................................................................................................................... 38 Chapter 4: Artificial Neural Network Modelling .....................................................................39 4.1 Overview .................................................................................................................... 39 4.2 Biological neural network .......................................................................................... 39 4.3 Artificial neuron ......................................................................................................... 40 4.4 Transfer function ........................................................................................................ 42 4.5 One-layer artificial neural network ............................................................................ 44 viii 4.6 Artificial neural network layers ................................................................................. 45 4.7 Artificial neural network architecture ........................................................................ 46 4.8 Artificial neural network learning algorithm ............................................................. 47 4.8.1 Backpropagation ........................................................................................................ 48 4.8.2 Artificial neural network parameters update .............................................................. 51 4.9 Dataset division .......................................................................................................... 52 4.10 Artificial neural network generalization .................................................................... 52 4.11 Summary .................................................................................................................... 56 Chapter 5: Model Development ..................................................................................................57 5.1 Overview .................................................................................................................... 57 5.2 Methodology of ARIMA model ................................................................................ 57 5.3 Methodology of ARIMAX modelling ....................................................................... 63 5.4 Methodology of dynamic artificial neural network modelling .................................. 65 5.4.1 Pre-training steps........................................................................................................ 67 5.4.1.1 Data collection .................................................................................................. 67 5.4.1.2 Data preprocessing ............................................................................................ 67 5.4.1.2.1 Outlier detection .......................................................................................... 67 5.4.1.2.2 Wavelet transform ....................................................................................... 68 5.4.1.2.3 Data normalization ...................................................................................... 71 5.4.1.3 Choice of network architecture ......................................................................... 72 5.4.1.4 Choice of network architecture specifics .......................................................... 73 5.4.1.5 Choice of optimum lag of inputs ...................................................................... 73 5.4.2 Training steps ............................................................................................................. 75 ix 5.4.2.1 Choice of learning algorithm ............................................................................ 75 5.4.2.2 Network parameters initialization ..................................................................... 75 5.4.2.3 Stopping criteria ................................................................................................ 76 5.4.3 Post-training steps ...................................................................................................... 77 5.4.3.1 Network output denormalization ...................................................................... 77 5.4.3.2 Choice of network performance measure ......................................................... 77 5.5 Summary .................................................................................................................... 79 Chapter 6: Results and Discussion .............................................................................................80 6.1 Overview .................................................................................................................... 80 6.2 Assemble datasets ...................................................................................................... 80 6.3 Case study one: City of Kelowna ............................................................................... 81 6.3.1 ARIMA models .......................................................................................................... 89 6.3.1.1 Estimation of ARIMA models .......................................................................... 95 6.3.1.2 Diagnostic statistical tests of ARIMA models .................................................. 97 6.3.1.3 Predicting with ARIMA models ..................................................................... 101 6.3.2 ARIMAX models ..................................................................................................... 104 6.3.2.1 Estimation of ARIMAX models ..................................................................... 107 6.3.2.2 Diagnostic statistical tests of ARIMAX models ............................................. 108 6.3.2.3 Predicting with ARIMAX models .................................................................. 112 6.3.3 Dynamic artificial neural network models ............................................................... 115 6.3.3.1 Preprocess data................................................................................................ 117 6.3.3.2 Choice of network architecture ....................................................................... 120 6.3.3.3 Choice of learning algorithm .......................................................................... 122 x 6.3.3.4 Network generalization ................................................................................... 122 6.3.3.5 Stopping criteria .............................................................................................. 122 6.3.3.6 Choice of optimum input lags and network specifics ..................................... 123 6.3.3.7 NAR networks ................................................................................................ 123 6.3.3.8 NARX networks.............................................................................................. 131 6.3.4 Comparison of ARIMA and DANN models............................................................ 137 6.4 Intervention analysis, Case study two: Springfield Road Corridor ......................... 142 6.5 Summary .................................................................................................................. 157 Chapter 7: Conclusion ...............................................................................................................158 7.1 Summary of the thesis research ............................................................................... 158 7.2 Contributions of this research .................................................................................. 160 7.2.1 Contribution 1 – Methodology to develop dynamic artificial neural networks ....... 160 7.2.2 Contribution 2 – Apply dynamic artificial neural networks .................................... 161 7.2.3 Contribution 3 – Case study application of dynamic artificial neural networks to evaluate safety interventions .................................................................................... 162 7.3 Limitations of this research ...................................................................................... 163 7.3.1 Application of sensitivity analysis ........................................................................... 163 7.3.2 Application of popular generalization techniques ................................................... 163 7.3.3 Application of unsupervised learning for clustering analysis .................................. 164 7.4 Future research ......................................................................................................... 164 Appendices ..................................................................................................................................178 Appendix A Average number of collisions for each month ................................................... 178 Appendix B Dynamic artificial neural network modelling code ............................................ 179 xi Appendix C Preliminary results of DANN modelling ............................................................ 183 Appendix D Performance of artificial neural networks .......................................................... 187 xii List of Tables Table 2.1 Summary of time series statistical modelling approaches in road safety studies ........ 12 Table 4.1 Transfer functions ........................................................................................................ 42 Table 5.1 Attribute of autocorrelation and partial autocorrelation function for non-seasonal ARIMA model ........................................................................................................... 60 Table 5.2 Attribute of autocorrelation and partial autocorrelation function for seasonal ARIMA model .......................................................................................................................... 60 Table 6.1 Summary statistics for monthly number of collisions in the City of Kelowna ............ 82 Table 6.2 Summary statistics for explanatory variables .............................................................. 83 Table 6.3 Spearman’s rank correlation between collision series data and candidate explanatory variables ..................................................................................................................... 89 Table 6.4 Outliers for the monthly collision series in Kelowna .................................................. 90 Table 6.5 Tentative ARIMA models for monthly collision series in Kelowna ........................... 95 Table 6.6 Dynamic estimated coefficients of ARIMA models for monthly collision series in Kelowna ..................................................................................................................... 96 Table 6.7 Ljung-Box statistic of selected ARIMA models........................................................ 101 Table 6.8 Performance of ARIMA models for monthly collision series in Kelowna ............... 102 Table 6.9 Explanatory and dynamic estimated coefficients of ARIMAX models for monthly collision series in Kelowna ...................................................................................... 107 Table 6.10 Ljung-Box statistic of selected ARIMAX models ................................................... 112 Table 6.11 Performance of ARIMAX models for monthly collision series in Kelowna .......... 112 xiii Table 6.12 Mean squared error and standard deviation of NAR networks corresponding to monthly number of TOTAL collisions in Kelowna ................................................. 125 Table 6.13 Mean squared error and standard deviation of NAR networks corresponding to monthly number of PDO collisions in Kelowna ...................................................... 126 Table 6.14 Mean squared error and standard deviation of NAR networks corresponding to monthly number of KSI collisions in Kelowna ....................................................... 127 Table 6.15 Optimal NAR network configurations of approximated (A2) and detailed (D1 and D2) subseries corresponding to each collsion type in Kelowna .............................. 128 Table 6.16 Performance of NAR networks corresponding to each collsion type in Kelowna .. 129 Table 6.17 Mean squared error and standard deviation of NARX networks corresponding to monthly number of TOTAL collisions in Kelowna ................................................. 132 Table 6.18 Mean squared error and standard deviation of NARX networks corresponding to monthly number of PDO collisions in Kelowna ...................................................... 133 Table 6.19 Mean squared error and standard deviation of NARX networks corresponding to monthly number of KSI collisions in Kelowna ....................................................... 134 Table 6.20 Optimal NARX network configurations of approximated (A2) and detailed (D1 and D2) subseries corresponding to each collsion type in Kelowna .............................. 135 Table 6.21 Performance of NARX networks corresponding to each collsion type in Kelowna 135 Table 6.22 A comparison of models’ performances on testing set ............................................ 140 Table 6.23 Summary of safety measures implemented in each phase ....................................... 143 Table 6.24 Outliers of TOTAL collision series for the Springfield Road Corridor corresponding to the training data.................................................................................................... 146 xiv Table 6.25 Mean squared error and standard deviation of NAR networks corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase one .................................................................................................................................. 148 Table 6.26 Mean squared error and standard deviation of NAR networks corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase two .................................................................................................................................. 149 Table 6.27 Optimal NAR network configurations of the approximated and detailed subseries corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase one .................................................................................................. 149 Table 6.28 Optimal NAR network configurations of approximated and detailed subseries corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase two .................................................................................................. 149 Table 6.29 Performance of NAR network corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase one ........................................... 150 Table 6.30 Performance of NAR network corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase two ........................................... 150 Table 6.31 Descriptive statistics of intervention analyses ......................................................... 154 Table 6.32 Seasonal effectiveness in the number of TOTAL collisions due to the intervention(s) .................................................................................................................................. 156 Table A.1 Average number of collisions for each month over 13.5 years for different collision series in the City of Kelowna ................................................................................... 178 Table C.1 Coefficient of determination of NAR networks for TOTAL collisions in Kelowna 184 Table C.2 Coefficient of determination of NAR networks for PDO collisions in Kelowna ..... 185 xv Table C.3 Coefficient of determination of NAR networks for KSI collisions in Kelowna ....... 186 Table D.1 Mean squared error of NAR networks for the A2 subseries for TOTAL collisions . 188 Table D.2 Mean squared error of NAR networks for the D1 subseries for TOTAL collisions . 189 Table D.3 Mean squared error of NAR networks for the D2 subseries for TOTAL collisions . 190 Table D.4 Mean squared error of NAR networks for the A2 subseries for PDO collisions ...... 191 Table D.5 Mean squared error of NAR networks for the D1 subseries for PDO collisions ...... 192 Table D.6 Mean squared error of NAR networks for the D2 subseries for PDO collisions ...... 193 Table D.7 Mean squared error of NAR networks for the A2 subseries for KSI collisions ........ 194 Table D.8 Mean squared error of NAR networks for the D1 subseries for KSI collisions ........ 195 Table D.9 Mean squared error of NAR networks for the D2 subseries for KSI collisions ........ 196 Table D.10 Mean squared error of NARX networks for the A2 subseries for TOTAL collisions .................................................................................................................................. 197 Table D.11 Mean squared error of NARX networks for the D1 subseries for TOTAL collisions .................................................................................................................................. 197 Table D.12 Mean squared error of NARX networks for the D2 subseries for TOTAL collisions .................................................................................................................................. 198 Table D.13 Mean squared error of NARX networks for the A2 subseries for PDO collisions . 198 Table D.14 Mean squared error of NARX networks for the D1 subseries for PDO collisions . 199 Table D.15 Mean squared error of NARX networks for the D2 subseries for PDO collisions . 199 Table D.16 Mean squared error of NARX networks for the A2 subseries for KSI collisions ... 200 Table D.17 Mean squared error of NARX networks for the D1 subseries for KSI collisions ... 200 Table D.18 Mean squared error of NARX networks for the D2 subseries for KSI collisions ... 201 xvi Table D.19 Mean squared error of NAR networks for the A2 subseries for TOTAL collision series before installing phase one in Springfield Road Corridor ............................. 201 Table D.20 Mean squared error of NAR networks for the D1 subseries for TOTAL collision series before installing phase one in Springfield Road Corridor ............................. 202 Table D.21 Mean squared error of NAR networks for the D2 subseries for TOTAL collision series before installing phase one in Springfield Road Corridor ............................. 202 Table D.22 Mean squared error of NAR networks for the A2 subseries for TOTAL collision series before installing phase two in Springfield Road Corridor ............................. 203 Table D.23 Mean squared error of NAR networks for the D1 subseries for TOTAL collision series before installing phase two in Springfield Road Corridor ............................. 203 Table D.24 Mean squared error of NAR networks for the D2 subseries for TOTAL collision series before installing phase two in Springfield Road Corridor ............................. 204 xvii List of Figures Figure 1.1 Monthly number of collisions in the City of Kelowna ................................................. 3 Figure 2.1 Linear threshold neuron .............................................................................................. 26 Figure 4.1 An artificial neuron..................................................................................................... 41 Figure 4.2 Transfer function graphs............................................................................................. 43 Figure 4.3 A layer of neurons ...................................................................................................... 45 Figure 4.4 Architecture of feedforward network with one hidden layer ..................................... 46 Figure 4.5 Architecture of recurrent network .............................................................................. 47 Figure 5.1 Flowchart of building ARIMA model ........................................................................ 58 Figure 5.2 Flowchart of training dynamic artificial neural network (DANN) ............................ 66 Figure 5.3 Wavelet decomposition tree ....................................................................................... 71 Figure 6.1 Time series plots of monthly number of collisions in the City of Kelowna............... 85 Figure 6.2 Seasonal subseries plot of monthly collision series in the City of Kelowna .............. 87 Figure 6.3 Collision series division for developing ARIMA approach in Kelowna ................... 90 Figure 6.4 Autocorrelation functions (ACFs) of monthly collision series for the training set .... 91 Figure 6.5 Autocorrelation functions (ACFs) of differenced collision series ............................. 93 Figure 6.6 Partial autocorrelation functions (PACFs) of differenced collision series ................. 94 Figure 6.7 Autocorrelation functions (ACFs) of residuals for selected ARIMA models ............ 98 Figure 6.8 Standardized residuals of selected ARIMA models ................................................... 99 Figure 6.9 Normal probability plots of selected ARIMA models ............................................. 100 Figure 6.10 Estimations of ARIMA models for monthly collision series in Kelowna.............. 103 xviii Figure 6.11 Cross-correlation functions (CCFs) of the monthly employment (EMP) series against collision series ............................................................................................. 105 Figure 6.12 Cross-correlation functions (CCFs) of pre-whitened monthly employment (EMP) series against filtered collision series ....................................................................... 106 Figure 6.13 Autocorrelation functions (ACFs) of residuals for selected ARIMAX models ..... 109 Figure 6.14 Standardized residuals of selected ARIMAX models ............................................ 110 Figure 6.15 Normal probability plots of selected ARIMAX models ......................................... 111 Figure 6.16 Estimations of ARIMAX models for monthly collision series in Kelowna ........... 114 Figure 6.17 Collision series division for developing DANN approach in Kelowna ................. 115 Figure 6.18 Flowchart of dynamic artificial neural network (DANN) coupled with discrete wavelet transform (DWT) ........................................................................................ 118 Figure 6.19 Approximated (A2) and detailed (D1 and D2) subseries obtained by using a two-level discrete wavelet transform (DWT) of monthly collision series in Kelowna ........... 119 Figure 6.20 Approximated (A2) and detailed (D1 and D2) subseries obtained by using a two-level discrete wavelet transform (DWT) of monthly employment (EMP) number in the Province of British Columbia............................................................................. 120 Figure 6.21 Dynamic artificial neural networks (DANNs) ....................................................... 121 Figure 6.22 Tapped delay line ................................................................................................... 121 Figure 6.23 Estimations of NAR networks corresponding to monthly collision series in Kelowna .................................................................................................................................. 130 Figure 6.24 Estimations of NARX networks corresponding to collision series in Kelowna .... 136 Figure 6.25 A comparison of the observed monthly number of collisions versus predicted ones over the testing set.................................................................................................... 138 xix Figure 6.26 Mean absolute percentage error of the observed values of each collision series versus predicted ones for different quarterly predicting horizons ........................... 141 Figure 6.27 Approximated (A2) and detailed (D1 and D2) subseries obtained by using a two-level discrete wavelet transform (DWT) of TOTAL collision series in the Springfield Road Corridor .......................................................................................................... 147 Figure 6.28 Estimations of two NAR networks corresponding to TOTAL collision series in the Springfield Road Corridor ....................................................................................... 151 Figure 6.29 Intervention plots for installing safety management phases in the Springfield Road Corridor. The dashed line represents the time of phase two implementation .......... 152 Figure 6.30 Time series plots of the difference between predicted and observed values by two NAR networks.......................................................................................................... 153 Figure 6.31 Available monthly average daily traffic (ADT) for Springfield Road Corridor .... 156 xx List of Symbols ∂E∂wijL partial derivative of error function with respect to each weight ∂fL(sjL)∂sjL derivative of transfer function to summation A scale (a.k.a. dilation) factor Aj and Bj amplitudes α and β regularization parameters a neuron output α learning rate ai denormalized network output ain original network output B shift (a.k.a. translation or position) factor B back-shift operator β0 intercept β1 trend β2 level change between regions 1 and 2, and region 3 β3 level change between region 1 and, regions 2 and 3 β4 trend change between region 1 and, regions 2 and 3 β6 null hypothesis of: H0 = l0 – l1 = 0 β7 null hypothesis of: H0 = δ0 – δ1 = 0 βjt and λt coefficients c constant xxi CWT(a, b) wavelet coefficients γ vector of indicators corresponding to the seasonal effect γi seasonal indicator variable γt seasonal component y̅ mean value of time series yt collision time series ŷt predicted value of yt γxy sample cross-covariance Δt time interval δ vector of coefficients corresponding to the time trend δ0 trend change in treatment series δ1 trend change in comparison series εt, ξt, and ζt irregular components EMPt-1 number of employments at month t-1 et sum of discrepancies between target and corresponding network output F(xk) first derivative of F(x) estimated at point xk ▽d difference operator S D seasonal difference operator fmin minimum range of nonlinear transfer function fmax maximum range of nonlinear transfer function Η matrix of coefficients corresponding to the explanatory variable ΘQ(BS) seasonal moving average operator θ1, θ2, …, θq coefficients of the moving average orders xxii θq(B) moving average operator int integer K window half-width KSIt number of KSI collisions at month t L number of layer L number neurons in hidden layer L0 level change in treatment series L1 level change in comparison series µt trend M number of coefficients mt median of data window N number of network weights N1 number of inputs NL number of neurons in layer L n number of data points PDOt number of PDO collisions at month t p maximum autoregressive order ø1, ø2, …, øp coefficients of the autoregressive orders øp(B) autoregressive operator Q number of training sets R number of inputs connected to hidden layer rc ratio of collisions in after to before intervention period on comparison road ri ratio of collisions in after to before intervention period on treated road xxiii S matrix of coefficients corresponding to the seasonal effect σω standard deviation of the residuals σt standard deviation T threshold value T number of points in a seasonal period T length of time series xt TOTALt number of TOTAL collisions at month t t time t target output tmin minimum of target set tmax maximum of target set ut explanatory variable V number of validating samples ʋt slope Фt vector of coefficients corresponding to the autoregressors ΦP(BS) seasonal autoregressive operator W weights vector WT transpose of the weights vector ‖wq‖ Euclidean length of weights in row q wijL(k + 1) updated weight gained at epoch k X neuron inputs vector x1 time x2 1 if in regions 1 or 2 and -1 if in region 3 xxiv x3 1 if in region 1 and -1 if in regions 2 or 3 x2 1 if in treatment series and -1 if not x3 1 if it is pre-intervention and -1 if not xjt jth explanatory variable at time t xt pre and post-intervention period xin normalized value of xi xmin minimum of training set xmax maximum of training set ψj psi-weight ψ(t) mother wavelet. ωt intervention variable zjt explanatory variables xxv List of Abbreviations A Approximations AADT Average Annual Daily Traffic ACF AutoCorrelation Function ADF Augmented Dickey-Fully ADL Actively Driver Licenses AI Artificial Intelligence AIC Akaike Information Criteria AIV Actively Insured Vehicles ANN Artificial Neural Network AR AutoRegressive ARCH AutoRegressive Conditionally Heteroscedastic ARIMA AutoRegressive Interacted Moving Average ARIMAX AutoRegressive Interacted Moving Average with eXplanatory variables Bagging Bootstrap AggreGatING BIC Bayesian Information Criteria BLUE Best Linear Unbiased Estimator CBD Central Business District CCF Cross-Correlation Function CWT Continuous Wavelet Transform D Details DANN Dynamic Artificial Neural Network xxvi DWT Discrete Wavelet Transform EB Empirical Bayes EMP EMPloyment GLM Generalized Linear Model ICBC Insurance Corporation of British Columbia iid Independently and Identically Distributed INAR INteger-valued AutoRegressive MA Moving Average MAE Mean Absolute Error MAPE Mean Absolute Percentage Error MAXT Mean mAXimum Temperature MINT Mean mINimum Temperature ML Machine Learning MPE Mean Percentage Error MSE Mean Squared Error MT Mean Temperature NAR Nonlinear AutoRegressive NARX Nonlinear AutoRegressive with eXternal inputs NDP Number of Days with Precipitation OLS Ordinary Least Squares PACF Partial AutoCorrelation Function RMSE Root Mean Squared Error SAR Seasonal AutoRegressive xxvii SARIMA Seasonal Autoregressive Interacted Moving Average SD Standard Deviation SMA Seasonal Moving Average SPF Safety Performance Function SSTS Speed Sensitive Traffic Signals TDL Tapped Delay Line TP Total Precipitation UEMP UnEMPloyment VAR Vector AutoRegressive VKT Vehicle Kilometres Traveled VMS Variable Message Sign 1 Chapter 1: Introduction 1.1 Overview This chapter sets out the layout of this thesis. Section 1.2 reports global and Canadian status of road transport safety and its socioeconomic impacts. Section 1.3 presents an overview of relevant literature on road safety analysis, followed by the motivation of this research in the form of a thesis problem statement. Section 1.4 introduces three overarching objectives of the research on sustainable transport safety. Section 1.5 concludes this chapter with the thesis framework. 1.2 Road safety issues Injures due to road traffic collisions are major public health concern on a global scale. It was estimated that there are over 1.2 million people killed and roughly 50 million injured or disabled each year on the world’s roads. Road traffic collisions are the leading cause of death among young adults of 15 to 29 years old, and across all age groups, the 9th leading cause of death and predicted to rise to the 7th by 2030. The lethality of road traffic collisions is not only a human tragedy, but also a significant economic burden, averaging 3% of GDP worldwide [1]. Hence, the United Nations announced the decade of action for road safety 2011–2020 [2]. Although there has been a substantial decrease of collision fatalities in recent years among industrialized countries, the stress that road collisions place on society is still a serious challenge [3], [4]. For example, approximately 120,000 people were injured and 1,900 people lost their lives annually for the years 2009 to 2013 in Canada [5]. In the same period, 260,000 collision claims on 2 average were reported each year to the Insurance Corporation of British Columbia (ICBC) of which 52,000 included casualties [6]. The above issues, as well as the limitations on capital budget and time, demonstrate a strong need for a reliable and practical technique in the field of road safety to study and monitor progress in reaching safety targets. It is also necessary for road authorities and decision makers to have insight into road safety status in the future. 1.3 Problem statement Many road safety data are aggregated and studied over time, a technique that is known as time series analysis. The main feature of such data is the correlation between an observation and its previous values in time. This causes that the standard error of the regression coefficients may be underestimated or overestimated when linear or generalized linear models (GLMs) are built. Several statistical approaches have been proposed in the literature to address the time dependencies within the observed values. Those approaches have primarily focused on the prediction of linear and stationary data. However, in a real-world collision time series, there is a great deal of nonlinearity and nonstationarity, especially when the mean is rather small. For example, Figure 1.1 illustrates an overview of the collected monthly time series in this research (to be discussed in detail in Chapter 6) for three collision types in the City of Kelowna. All collision series depict an overall upward trend and inconstant mean over time, which suggest that the time series are nonstationary. In addition, the trend of two collision series (i.e. property damage only collisions and killed or seriously injured) are stochastic, as shown by the dashed lines in Figure 1.1. 3 Therefore, the traditional statistical models are unable to capture the nonlinear and nonstationary patterns inside the given collision series. Figure 1.1 Monthly number of collisions in the City of Kelowna (January 2003 to June 2016). The dashed lines represent the local average. An alternative approach is to use artificial intelligence (AI) models that have been successfully applied for modelling in various fields to cope with the assumptions of the traditional time series statistical approaches. This also improves the predictive performance of the model. The motivation of this thesis was to predict the number of collisions using the most popular AI model, the artificial neural network (hereafter denoted as ANN or neural network), as a substitute for the traditional time series models in road safety studies. Two case studies were defined for this thesis: (1) the entire City of Kelowna, where the ANN’s performance can be evaluated with the autoregressive integrated moving average (ARIMA) model 0200400600800100012002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of collisionsTime (year)Total collisions Property damage only collisions Killed or seriously injured collisions4 and (2) the Springfield Road Corridor safety study, where traditional time series statistical approaches are inapplicable because the normality assumption is not met. 1.4 Research objectives Building a dynamic ANN (DANN) as a nonlinear time series model is the main goal of this research. The objective of this research was defined to be three-fold: 1) Develop a methodology to conduct DANN modelling to predict collision series in road safety analysis, 2) Assemble the required datasets to apply a DANN approach and evaluate its performance, and 3) Address the City of Kelowna’s need for improved evaluation tools, by exploring a methodology to evaluate the impact of each safety intervention implemented on Springfield Road Corridor at various times. 1.5 Thesis framework The remainder of the thesis is organized as follows. Chapter 2 discusses the main properties of the time series and reviews the statistical approaches conducted on analyzing time series data in road safety. The preliminary developments of ANNs and limited studies of this model in road safety studies are also presented in this chapter. Chapter 3 is devoted exclusively to ARIMA modelling, including components and assumptions. Chapter 4 explains the fundamentals of ANN. Chapter 5 sets out the methodologies to design a univariate and multivariate ARIMA and DANN models. Chapter 6 presents the results of applying the DANN technique to the Kelowna dataset and associated models, including a discussion on their validity compared to traditional evaluation tools. Chapter 7 summarizes each chapter, provides the main findings and contributions of the thesis, 5 points out the limitations of the current research, and outlines directions for future research to refine the DANN modelling. 6 Chapter 2: Literature Review 2.1 Overview This chapter aims to review the main properties of time series data and relevant time series analysis in the literature to provide for a better picture of the problem and criticize the existing models. Subsequently, the history of ANNs and studies that developed this approach in road safety are investigated. A summary and general comments are given at the end of the chapter. Thus, this chapter is structured as follows. Section 2.2 explains types of road safety studies. Sections 2.3 and 2.4 present the concept of time series analysis and features of time series data, respectively. Section 2.5 provides statistical methods that were applied in road safety studies to regress collision time series. Section 2.6 presents a brief history of ANNs, followed by an investigation of road safety studies using this approach. Section Error! Reference source not found. presents a conclusion for the chapter. 2.2 Road safety studies The road safety of a geographical area (e.g. city, neighbourhood, road segment, or intersection) is defined as the expected number of collisions (by total or type) in a specified period [7]. Hakim et al. [8] classify road safety studies into three categories, according to the type of data being analyzed: 1) Cross-sectional studies: predicting how safe a location is through transferring knowledge about similar locations to that location, 2) Time series studies: predicting the future safety characteristics of a specific location by transferring knowledge across time, and 7 3) Longitudinal studies: containing time series observations that belong to several locations to deal with spatial and time correlations. This research focused on investigating data with time series analysis since the interest was to predict the collision counts of a location over time and evaluate the safety phases implemented in one location at several points in time. 2.3 Introduction to time series analysis Time series data is a chronological sequence of observations repeated through time [9]. Observations are most often collected at equidistant time intervals (Δt), denoted as a set of discrete values of {y1, y2, …, yt, yt+1, …, yn}, where t and n denote the time and number of observed samples, respectively. Many road traffic datasets appear as time series such as monthly aggregated collision counts or annual vehicle kilometres traveled (VKT). Unlike cross-sectional data, the order of observations in time series is important and changing the order would vary the meaning of the data. Since measured or observed values are often corrupted by noise in engineering practice, time series data generally includes a deterministic component and a stochastic component resulting from the errors that cause fluctuations around deterministic values. The analysis of a time series is primarily aimed at identifying and extracting the internal structure of the time series data to understand its dynamic meaningful patterns [10]-[12]. 8 The generalized term of time series analysis in road safety can be categorized into three groups: (1) predicting road safety trend development in the future, (2) describing the evolution of road safety over time, and (3) explaining the effect of selected factors (e.g. VKT, unemployment rate, or weather situation) or specific events (e.g. fuel shortage, economic crisis, or safety measures) on road safety through time [13]. Next section will outline the major features of time series data. 2.4 Overview of time series features The major features of time series data include seasonality, trend, cycle, irregular component, and autocorrelation. Seasonality or seasonal variation is a repeating and predictable fluctuating pattern that occurs within a fixed period such as one year. Seasonality can be observed in quarterly, monthly, weekly, daily, etc. depending on the nature of time series. For example, collision count is often higher in wintertime whereas lower in summertime due to weather conditions and daylight duration. Trend is an unrepeatable regular movement without time effect that represents an underlying level in a time series. Trend is steadier and its direction rarely changes when compared to seasonality. For example, collision series can indicate an upward trend due to a rising number of road trips or downward due to enhancement of standards for vehicle and road design. Cycle is a pattern that is longer than fixed period such as long-run climatic fluctuations due to flood or El Nino; however, cyclic behaviour cannot be found in traffic data. An irregular component is a random part of a time series, which is unpredictable. In other words, the residuals of a time series after modelling the patterns are referred to as an irregular component [10]. An intrinsic feature of time series data is time dependency between sampling points. Such dependency between an observation with its own time lags (i.e. calendar interval between 9 observations) is called autocorrelation (a.k.a. serial correlation or dependence) [10]. The autocorrelation phenomenon arises due to several reasons. First, the adjacent observations are influenced by the same conditions. This causes high values tending to be followed by high values and vice versa. Second, the error terms of the model are autocorrelated when a key input variable(s) is not considered and its successive values are correlated [14]. The autocorrelation between nearby points in time severely restricts the applicability of most conventional statistical approaches because those analyses are established on an independence assumption [11]. Based on the Gauss–Markov theorem, the ordinary least squares (OLS) method provides the best linear unbiased estimator (a.k.a. BLUE) for regression coefficients if the error terms are random, i.e. independently and identically distributed (iid). However, this assumption is violated when discussing a time series. In summary, autocorrelation is a distinguishing characteristic of time series from other types of data. The proposed time series analysis techniques in the literature will be presented in the following section. 2.5 Time series statistical models Time series statistical modelling is a valuable tool to identify the inherent fluctuation of the observed values for better understanding the underlying behaviour of the system. The statistical modelling of time series may take place in two general approaches [11]: 1) Time domain approach, which is based on an autocovariance function. In this approach, autocorrelation among adjacent observations may be addressed in two ways. First, deterministic time indices are added as explanatory and indicator variables to classical regressions that respectively represent trend and seasonality of the series (known as deterministic time series regression). Trend can be characterized by linear (Equation 2.1), 10 exponential (Equation 2.2), or polynomial (Equation 2.3) relation of a time index to the given series [15]: yt = αt + β + εt (2.1) yt = exp(αt + β + εt) (2.2) yt = αt + βt2 + γ + εt (2.3) where:  yt denotes time series,  t denotes the time,  α and β denote the unknown coefficients, and  εt denotes the error term. Seasonal variation can be modelled by adding S indicator variables corresponding to the effect of every season with no intercept (Equation 2.4), or S – 1 indicator variables corresponding to the seasons of 2 to S – 1 with an intercept corresponding to the first season (Equation 2.5) [15]: yt = αt + ∑γiSi=1+ εt (2.4) yt = αt + β + ∑γiS−1i=1+ εt (2.5) where:  γi denotes the seasonal indicator variables and the overall sum of seasonal indicators in one seasonal span should be always one. 11 Second, the current value of the dependent variable is related to its own past values (known as univariate dynamic regression) and, moreover, the past values of other series as explanatory variables (known as multivariate dynamic regression). The main advantage of dynamic regression over deterministic time series regression is that the features of a time series (e.g. trend and seasonality) can be adjusted for a dataset by the model rather than some fixed function of time in classical regression [16]. 2) Frequency domain approach, which is based on spectral density function. The major features of the series are related to harmonic terms (i.e. composite of sine and cosine waves) of varying periods and amplitudes. Both approaches yield the same results for long-term observed series while the time domain approach has better predictive performance in a short series like road traffic data [11]. It should be noted that although long collision series data is available in most cases, their reliability decreases especially in regard to nonfatal collisions as they go farther back in time [17]. Time series statistical modelling methods proposed in the literature, and their advantages and disadvantages are listed in Table 2.1. The time domain approach is more common in road safety analysis due to its simplicity and ease of application for the general practitioners. The preliminary time series studies in road safety and various models used in each approach are reviewed in the following sections. 12 Table 2.1 Summary of time series statistical modelling approaches in road safety studies Reference Model Advantages Disadvantages [18] Deterministic time series model  It is similar to GLMs.  It is easy to develop for general practitioners.  Trend and seasonality are considered deterministic.  Linearity assumption is a major drawback. [19] ARIMA model  It is easily interpreted.  It can be calculated by various software.  Prediction is straight forward.  Highly accurate predictions in short-term horizon.  Linearity and stationarity assumptions are major drawbacks.  It falls off in predicting long-term horizon.  It requires visual analysis.  There is no specific instruction for step-by-step model identification.  It is unstable in terms of changes in observations.  It might not be unique and can be equivalent to another model with different orders. [13] State space model  Stationarity of the data is not required.  Time series components can be explicitly modelled.  Time series components can be modelled either deterministic or stochastic.  It is very adaptive in terms of adding or removing components.  It can easily handle missing data.  Linearity assumption is a major drawback.  Some preliminary models should be developed based on user’s experience.  The software needs to be programmed by its users. [20] Vector autoregressive model  The stationarity assumption of time series data is not required.  A combination of outputs is modelled simultaneously.  There should be linear relationship among time series.  Trend and seasonality are considered deterministic. [21] Spectral analysis  Unlike time domain approach, it shows how a series changes over various frequencies.  Stationarity assumption is a major drawback.  Trend and seasonality are considered deterministic. 13 2.5.1 Preliminary time series studies in road safety By the 1970s, road safety studies had gained further attention with the gradual increase of motorization and road collision fatalities. The earliest attempts to describe the substantial issues of time series and intervention analysis in road safety were done by Glass [22], and Campbell and Ross [23], which were later characterized comprehensively by Glass et al. [24]: 1) History: This refers to the possibility that external factors, rather than the countermeasure, may change the road safety of before and after implementation. For example, there may be a dry year with little precipitation after implementation, consequently the number of collisions would decrease, even if the countermeasure had not been installed. 2) Maturation: This refers to a downward or upward collision trend over time, rather than change attributed to a safety countermeasure. For example, a reduction in collisions may not only be explained by not only the effect of a specific countermeasure, but also by road, vehicle, or drivers’ attribute enhancement over many years. Maturation is distinguished from history as it represents time trend rather than discrete events. 3) Testing: This refers to pre-intervention changes without implementing the countermeasure. For example, the collisions may reduce because drivers are informed by publicizing and consequently more cautious. 4) Instrumentation: This refers to changes in measuring or the type of recording instrument used. For example, the definition of property damage collision and the severity of collisions would change by the police or insurance companies. 5) Instability: This refers to change in observations due to numerous change-producing events. The instability is greater when the population is smaller. For example, instability in collision counts of a road segment is greater than that of a city. 14 2.5.2 Time domain models Four time-domain models were found in the literature: (1) deterministic time series model, (2) autoregressive integrated moving average (ARIMA) model, (3) state space model, and (4) vector autoregressive (VAR) model, which are explained in the following sections. 2.5.2.1 Deterministic time series model Before systematic time series analysis, the correlated residuals in classical linear or generalized linear regression models of a collision series were addressed by adding time indices that represent trend and seasonality. Ellingstad and Westra [18] investigated the effect of two interventions (1972 speed limit reduction in Nebraska State and 1973 energy shortage) on a seasonally adjusted property damage collision series. For this purpose, three change parameters were measured in this study, as shown in Equation 2.6: (1) time trend, which indicates the linear growth or decline over the entire collision series and deals with autocorrelation, (2) level change, which indicates the sudden shift in the collision mean right after each intervention, and (3) trend change, which indicates gradual modification of collision slope from before to after intervention. The change parameters were fitted by indicator variables in a linear regression. yt = β0 + β1x1 + β2x2 + β3x3 + β4x1x3 (2.6) where:  x1 denotes the explanatory time variable,  x2 denotes 1 if in regions 1 (before the first intervention) or 2 (between two interventions) and -1 if in region 3 (after the second intervention),  x3 denotes 1 if in region 1 and -1 if in regions 2 or 3,  β0 denotes the intercept, 15  β1 denotes the trend,  β2 denotes the level change between regions 1 and 2, and region 3,  β3 denotes the level change between regions 1 and 2, and region 3, and  β4 denotes the trend change between region 1 and, regions 2 and 3. Then, t-test was applied to evaluate if the mean is dissimilar in different regions of time series (before, between, and after two interventions). In addition, a comparison group, referring to a unit group that is not affected by the interventions of interest, was used to detect the effect of external factors given by Equation 2.7. However, assumptions of normality and constant variance were ignored in both linear regressions as the authors mentioned. yt = β0 + β1x1 + β2x2 + β3x3 + β4x1x2 + β5x1x3 + β6x2x3 + β7x1x2x3 (2.7) where:  x2 denotes 1 if in treatment series and -1 if not,  x3 denotes 1 if it is pre-intervention and -1 if not,  β6 tests the null hypothesis of: H0 = L0 – L1 = 0, where L0 and L1 are the level change in treatment series and comparison series, respectively, and  β7 tests the null hypothesis of: H0 = δ0 – δ1 = 0, where δ0 and δ1 are the trend change in treatment series and comparison series, respectively. In addition to the basic regression models, some applications of systematic time series models took place in road safety analysis. Three dedicated time domain models were proposed within the literature: (1) ARIMA model, (2) state space model, and (3) VAR model, which will be introduced in the next. 16 2.5.2.2 Autoregressive integrated moving average model One of the most popular systematic time series statistical approach in road safety analysis is the ARIMA model. This model is also referred to as the Box-Jenkins model as it was popularized by Box and Jenkins in 1976 for the first time [25]. The ARIMA model is a general form and has some subclasses that will be discussed in detail in Chapter 3. The first attempt on using an ARIMA model in road safety was done by Wiorkowski and Heckard [16] in 1977, who modelled six traffic indices for the State of Texas on a monthly basis: (1) vehicle miles in units of 100 million vehicle miles, (2) fatality rate, (3) fatal collision rate, (4) injury rate, (5) injury collision rate, and (6) severity index (fatalities or injuries). Next, the effect of two interventions (i.e. fuel shortage and speed limit reduction) were evaluated in this study. Despite the widespread use of ARIMA modelling, many researches document its drawbacks. For example, the founders of this technique pointed out that there is no specific formulation available for step-by-step model identification [26]. Thus, the researchers recommended using graphical methods and experienced judgements. However, visual analysis of graphs is quite difficult for inexpert users to find the best orders of differencing, autoregressive, and moving average terms. Moreover, Commandeur et al. [27] argued that to make data stationary, the trend and seasonal should be filtered out before modelling. Scott [28] compared ARIMA with a regression model built on an assumption of uncorrelated residuals. The author claimed that the ARIMA model cannot provide better results when the noise is very high and consequently significant autocorrelation is not found among residuals. The performance of the ARIMA model is also weak when trying to predict long horizons, especially if the data being analysed is nonlinear [29]. Furthermore, the ARIMA model usually performs well when the mean of aggregated count series 17 is rather high [30]. Quddus [31] suggested the application of the integer-valued autoregressive (INAR) Poison model instead. Although the use of INAR Poisson model outperforms the ARIMA for a time series exhibiting low count values, this technique suffers from equidispersed and stationary assumptions that are often not satisfied in real applications. 2.5.2.3 State space model State space models decompose a time series into several possible unobserved dynamic properties, which are referred to as states, to develop explicit models for both additive and multiplicative linear regressions. That is why this model is also known as structural time series analysis or unobserved components model. The basic components of this model may include trend (level and slope), seasonal, and irregular components to predict or describe time series. The components can also be extended by adding explanatory and intervention variables to show the underlying development of the series, as given by Equation 2.8 [13], [32], [33]: yt = µt + γt+∑βjtxjtkj=1+ λtωt + εt (2.8) where:  µt denotes the trend,  γt denotes the seasonal component,  xjt denotes the jth explanatory variable at time t,  ωt denotes the intervention variable,  βjt and λt denote the coefficients, and  εt denotes the irregular component and should be Gaussian white noise. 18 The corresponding multiplicative representation model is given by Equation 2.9: yt = µt γt∑βjtxjtkj=1 λtωt εt (2.9) The main characteristic of state space modelling is that a parameter can be related to the same parameter at the previous point in time by a dynamic linear regression in order to capture the correlation among observations [17]. The merit of state space models over classical regression models is the flexibility of its components to have either deterministic or stochastic forms. In other words, state space models can be viewed as a general form of classical regression models in which the time coefficients can be dynamic (i.e. change over time). For example, the state space model in the simplest form, which is called a local level model, can be given by Equation 2.10 and Equation 2.11 [13]: yt = αt+ʋt + εt (2.10) αt+1 = αt + ξt (2.11) where:  αt denotes the level at time t,  ʋt denotes the slope, and  εt and ξt, denote the irregular components. This model represents a random walk that means the observed value of time t is more like the observation at t+1 than other observations. In this model, the slope is treated deterministically in Equation 2.10, which is called the observation or measurement equation, and level is treated 19 stochastically in Equation 2.11, which is called state or system equation. The slope can also be considered stochastic like in the following set of equations, which is called local trend model: yt = αt+εt (2.12) αt+1 = αt + ʋt + ξt (2.13) ʋt+1 = ʋt + ζt (2.14) where:  ζt denotes the irregular component. The idea of state space time series analysis was introduced by Kalman [34] and the technique for estimation, analysis, and prediction is called Kalman filtering. The state space model was first presented in road safety by Harvey and Durbin [32] in 1986 to investigate the introduction of compulsory wearing of a seat belt in the United Kingdom. The same methodology was later used to compare road safety progress in ten European countries [3]. In this study, the Smeed’s model representing the rate of fatalities as a function of population and number of registered vehicles was extended by a stochastic trend as an indicator of the rate progress in road safety. The calculation process of state space models was met with a variety of computational bottlenecks when first developed. Gradually the method gained popularity in analyzing road safety problems as statistical software improved. For example, Scuffham and Langley [35] investigated the quarterly number of traffic fatalities in New Zealand. The trend behaved stochastically with two state equations representing the level and slope. Also, the seasonal component was modelled stochastically with the sum of individual season in a trigonometric form. However, the effect of explanatory variables (including the unemployment rate, GDP, proportion of motorcycles, 20 proportion of young males in the population, alcohol consumption, and open road speed limit) and intervention variables (including the 1973 and 1979 oil crises, and seat belt legislation) were considered deterministic. Nghiema, Commandeur, and Connelly [36] recently examined the effect of several road safety interventions intended to reduce impaired driving on traffic fatalities using state space model to control the unobserved components. Furthermore, multivariate state space time series analysis can be adopted to model a risk system in road safety, which is known as latent risk time-series model. Risk here refers to as [37]: (1) exposure of an event, (2) probability that the event occurs, and (3) expected severity of the event. In other words, state space models allow for formulating three factors of exposure, risk, and severity simultaneously. For example, Bijleveld et al. [37] presented a multivariate state space model to investigate the risk of moped and bicycle users that might be killed or hospitalized in a collision in the Netherlands. The same methodology was later employed on the Smeed’s model to allow for a nonlinear relationship between variables in order to analyze the risk of fatal collisions and exposure, defined here as kilometers driven by motor vehicles [38]. In another study, multivariate state space models were used for the analysis and prediction of macroscopic road safety trends in Greece [39]. 2.5.2.3.1 Advantages and disadvantages of state space models In addition to the capability of modelling each component separately in order to observe the data in more detail, state space modelling has other advantages. A state space model is a general approach to predict nonstationary time series such that trend and seasonal components can be modelled explicitly. However, those components are first removed by difference filtering in 21 ARIMA models [40]. State space models are very adaptive for different problems by adding or removing components. Due to the development of statistical packages, the complicated state space models can be calculated relatively easy. The prediction computation of state space modelling is smooth as the smoothing method is used for fitting through the Kalman filter [41]. However, the fundamental tenet of state space modelling is the linear relationship between components. Moreover, it is rather challenging to develop a state space model because some preliminary structural models, which are built based on prior knowledge, need to be conducted first to investigate how each component behaves either deterministically or stochastically. It contrasts with the ARIMA models that have clearer instruction. It should be noted that any state space model can be represented by an ARIMA model with the same predictors but with different coefficients. In fact, they are just two different systematic classes of model fitting for identifying coefficients [27]. For example, it can be demonstrated that the structural models presented in Equation 2.10 and Equation 2.11 are equivalent to ARIMA(0, 1, 1). 2.5.2.4 Vector autoregressive model The VAR model is a general form of a univariate autoregressive model that allows the inclusion of multiple time series in the same model (a.k.a. simultaneous dynamic equations or multivariate transfer function models). In other words, the VAR model is a multivariate time series analysis that regresses each dependent variable on its own lags and on the lags of other dependent variables. The objective of VAR models is not only to predict the future according to previous values, like univariate time series analysis, but also to obtain insight into the dynamic structure of a system [42]. For example, Michalaki et al. [43] applied the VAR model to enable simultaneous modelling 22 of two monthly aggregated collision series of the emergency lane (called hard-shoulder lane in this study) and main traffic lanes, while also examining the possible impact that they had on each other. Serhiyenkoa et al. [20] used the VAR model to predict and uncover temporal trends of compositional pedestrian collisions for several severity levels, as shown in Equation 2.15: (1) fatal, (2) severe injury, (3) evident minor injury, and (4) non-evident possible injury and property damage only. The value of each level is calculated as the proportion of total monthly pedestrian collisions. Then compositional time series is modelled as a vector autoregressor, time trend, and seasonal matrix, with vehicle miles travelled as the explanatory variable. yt = ∑ Φkxt−k + ηut + γt + Sδ + ωtpk=1 (2.15) where:  yt denotes the vector of injury severity levels at time t,  p denotes the maximum autoregressive order,  Фt denotes the vector of coefficients corresponding to the autoregressors,  η denotes the matrix of coefficients corresponding to the explanatory variable,  ut denotes the explanatory variable,  δ denotes the vector of coefficients corresponding to the time trend,  S denotes the matrix of coefficients corresponding to the seasonal effect,  γ denotes the vector of indicators corresponding to the seasonal effect, and  ωt denotes the vector of white noises. 23 2.5.3 Frequency domain models An alternative way to express the periodic behavior of a stationary time series is a sum of sines and cosines at the set of possible frequencies ωj = j/n where j = 1, 2, …, n/2: yt = ∑(Aj cos(2πωjt) + Bj sin(2πωjt))n/2j=1 (2.16) where:  n denotes the number of observations (assume that n is even) and  Aj and Bj denote the amplitudes representing seasonal variation. Aj and Bj are random independent variables with a mean of zero, which can be estimated as a classical regression. Time domain approach is referred to as a regression of the present time on past values while frequency domain approach is the regression of the present time on a linear combination of sinusoidal functions [11]. Meyerhoff [44] conducted this approach in 1978 to reveal the repetitive patterns of fatal collisions in the United States using the Fourier series. The underlying frequencies were identified by a power spectrum. Then, the impact of daylight saving time, during March and April of 1974, was compared to the same time in 1973, which time was not adjusted, and investigated by two comparisons: (1) transition data: for the period immediately before and after daylight saving time, and (2) year-to-year (filtered) data: for the same periods in two or more different years, one with intervention and the other without it. Novoa et al. [45] also investigated the effectiveness of road safety prioritisation in the year 2004 on a number of road traffic injuries in Spain using an interrupted time series design. In this study, a quasi-Poisson regression model was developed to 24 address the presence of overdispersion. Linear trend and seasonality were handled by considering time index and sinusoidal functions in the model, respectively: ln(yt) = β0 + β1t + ∑ [β2k sin(2kπtT) + β3kcos(2kπtT)]6k=1+ β4xt + β5xtt (2.17) +∑(β6jzjt) + εtj where:  yt denotes the monthly number of injured people,  T denotes the number of points in a seasonal period (i.e. T = 12),  xt denotes the pre and post-intervention period (i.e. xt = 0 for pre-intervention period and xt = 1 for post-intervention), and  zjt denotes the explanatory variables. 2.5.4 Nonlinear time series models Road safety researchers have rarely focused on developing nonlinear alternative approaches to the classical linear time series statistical models. For example, Yannis et al. [46] introduced an autoregressive version of the Smeed’s model to predict the fatality rate of sixteen European countries. Four models were evaluated in this study: (1) Smeed’s original model, (2) extended version of the Smeed’s model by adding the first autoregressive orders of dependent and independent variables (i.e. ratio of number of fatalities to the number of vehicles and number of vehicles to the population, respectively), (3) log-transformed Smeed’s model, and (4) log-transformed autoregressive Smeed’s model. The results revealed that the original and log-transformed Smeed’s model are unable to capture the autocorrelations inside observations, 25 whereas the autoregressive models generally demonstrate satisfactory results with exceptions in only a few countries. In addition, both autoregressive models yield a more satisfactory performance on the average of sixteen countries compared to the base models. Despite the nonlinear relationship among the variables, the nature of this model is still a linear regression because the OLS method was used to estimate the unknown coefficients. The dependencies of higher autoregressive orders, on the other hand, were ignored, and this may explain why the residuals of some autoregressive models were autocorrelated. Given all the above weaknesses of time series statistical models developed to date, some researchers have begun searching for an artificial intelligence (AI) model like ANN that holds potential to address more of the nonlinear and irregular relationship. The preliminary structures of ANNs will be outlined in the next section to obtain a fundamental understanding of how they evolve into today’s networks, but the details of ANN models will be discussed in Chapter 4. 2.6 History of artificial neural networks A good definition of ANN was given by Haykin [47] as a combination of computing units understanding the system through a statistical learning process and keeping the knowledge in its connections. Erguo and Jinshou [48] also described ANN as a form of machine learning (ML) that is considered to be a nonlinear data modelling tool instead of linear correlative regression between output variable and input variable(s). The first modern view of the artificial neurons was introduced by McCulloch, a neuroscientist, and Pitts, a logician, [49] in 1943, in terms of one basic principle: The output of an artificial neuron (a 26 in Figure 2.1) can have two states, on or off. The output is one when the weighted sum of neuron inputs is equal to or greater than a threshold value, otherwise the output would be zero: a = {1 if WTX ≥ T0 otherwise (2.18) where:  a denotes the neuron output,  W = (w1, w2, …, wn) denotes the weights vector,  WT denotes the transpose of the weights vector,  X = (x1, x2, …, xn) denotes the neuron inputs vector, and  T denotes the threshold value. Figure 2.1 Linear threshold neuron The McCulloch-Pitts network is a highly-simplified model that can only generate a binary output. Also, due to the lack of learning algorithm, the threshold and weight values are fixed, which was not as sophisticated as today’s ANN structures. T a xn x1 x2 wn w1 w2 27 The Hebbian rule was proposed in 1949 as the first learning law to train neural networks, although Hebb [50] used this theory to explain a conditional learning mechanism for biological neurons in the brain. The Hebbian learning states that: “when an axon of cell A is near-enough to excite cell B and when A repeatedly and persistently takes part in firing of B, then some growth process or metabolic change takes place in one or both these cells such that the efficiency of cell A is increased.” The mathematical interpretation of the Hebbian rule is that the weight is increased when the input (xi) and output (a) of a neuron have both positive or negative values, and vice versa: winew = wiold + xia (2.19) The Hebbian rule is not commonly used today; however, all ANN learning algorithms work by changing weights of the network. In addition, the Hebbian rule variations such as unsupervised Hebbian rule or delta rule are still very popular in training ANNs. The first true application of ANNs began by introducing the perceptron network by Rosenblatt [51] in 1958. Further pioneer work was done by Widrow and Hoff [52] by presenting an adaptive perceptron-like network and a least mean square learning algorithm that is more practical than perceptron learning. ANN gained the general interest of researchers by the early success of these two networks in pattern recognition. However, they could only solve limited problems because both are only linear classifiers. Also, the perceptron network can converge to the correct weighting sets if the weights exist. The inherent limitations of the perceptron and least mean square algorithm, which is widely criticized by the Minsky and Papert’s book [53], ceased further study on ANNs for over a decade. 28 With the arrival of the 1980s and development of powerful personal computers, neural networks gained attention once again among researchers. One of the major contributions toward the rebirth of ANNs was a backpropagation algorithm for training multilayer perceptron networks. The process of this algorithm was first explained in the doctoral thesis of Werbos [54] in 1974. It was also rediscovered independently in other studies like [55], but not publicized until 1986 by Rumelhart and McClelland [56]. Backpropagation is a generalization of the least mean square rule by allowing modellers to calculate the derivatives of errors with respect to the weights in hidden layers. In fact, least mean square is an explicit derivative of the squared error in a single-layer network. This prevents nonlinear discriminating capability of networks with hidden layers. Currently, a multilayer network with a backpropagation learning algorithm is the most conventional ANN model since it can be used as a universal function approximator [57]. 2.6.1 Artificial neural network in road safety analysis ANN has been successfully applied for solving problems in a wide range of research areas. The most important applications of ANNs are function approximation (a.k.a. nonlinear regression), pattern recognition, clustering, and prediction (a.k.a. time series analysis or dynamic modelling) [57]. However, it has not been deeply studied in correlation to road safety and is a new analysis technique in this field. Moreover, some limited works that have been conducted in this field have focused on modelling cross-sectional data and not time series ones. For example, Chang [58] used ANN to find the relationship between traffic collisions and risk factors, including highway geometric variables, traffic characteristics, and environmental features in various segments of a highway. The predictive performance of the ANN models was compared with the most widely used collision prediction models, i.e. Poisson and negative binomial regression models. The results 29 indicated similar performances for both approaches. Nevertheless, he concluded that ANN could be an alternative approach to overcome the underlying constraint of GLMs, i.e. the pre-defined relationship between the dependent and independent variables. Moreover, the interaction between independent variables might inflate the estimated coefficients in GLMs, whereas it is not a critical concern in ANN modelling. Chiou [59] proposed an ANN-based expert system for collision appraisal. In this study, five degrees of responsibilities (i.e. full, major, equal, minor, and no responsibility) of all parties involved in two motor-vehicle collision cases were modelled based on police investigation reports. A discrimination analysis model was then calibrated for comparison. The study found that the ANN model had a higher correlation coefficient and correctness rate (i.e. the proportion of correctly predicted parties) in both the training and validating set. Pande and Abdel-Aty [60] explored the applicability of ANNs as a classifier in a traffic management system. The authors used the probabilistic neural network to identify the conditions prone to rear-end collisions on the freeway. In another study, the ANN model evaluated the safety level of a traffic collision scene. A collision scene has the potential of any secondary collisions due to drivers’ dangerous behaviours like lane changing. A five-grade evaluation index (from better to worse condition) was developed to prevent any secondary collisions around a collision scene [61]. Luca [62] compared ANN and calibrated multivariate analysis model presented in the Highway Safety Manual 2009. The findings demonstrated that ANN predictive performance is better in terms of total sum of residuals. A series of studies investigated ANN as the modelling approach in the analysis of collision-related injury severity in the event of a collision. Abdelwahab and Abdel-Aty [63] modelled the driver injury severity with crash factors related to the driver, vehicle, roadway, and environmental 30 characteristics. The study field was restricted to two-vehicle collisions that occurred in signalized intersections. The collisions were further classified into one of three injury-based severity categories. A multilayer neural network was compared to the performance of the ordered logit model and proved to have a better classification performance. The power of ANN in the analysis of collision injury severity was confirmed by several subsequent studies. For example, Delen et al. [64] presented a series of neural network models to associate the risk of injury to crash factors and then further conducting a sensitivity analysis on trained neural network models. To retrieve the information that is required to find the true relationship between crash-related factors and injury severity, the five-class prediction was transferred into a set of binary prediction models. The result from this study verified the strong performance of ANNs in modelling injury severity levels. Liu et al. [65] used ANN to predict the injury occurrence based on demographic and environmental factors involving in the collision. The findings revealed that the optimum network can correctly predict approximately 69% of the injured and non-injured cases. A sensitivity analysis was also conducted to determine the importance of each predictor in classifying the injury occurrence of two motor-vehicle collisions. In another study, the performance of ANN was compared to the ordered logit model in predicting the collision injury severity (i.e. injury, possible injury, no-capacitating injury, and capacitating injury or fatality) of two-vehicle crashes [66]. The results were validated by the findings of another study to model four severity classes by sixteen attributes collected at the time of each collision [67]. In this study, data was first clustered to improve the predictive performance, then modelled by the ANN classifier. Sameen et al. [68] designed a recurrent neural network to predict the injury severity of traffic collision based on 1130 collision records. The proposed model was compared with multilayer perceptron and Bayesian logistic regression models. The accuracy of the recurrent neural network was higher than two comparison 31 models. The above studies only applied ANNs as the function approximation to predict the severity of a collision when it occurs; however, the core of this research is aimed at designing an ANN model to predict the number of collisions in future. 2.7 Summary This chapter reviewed the most common statistical regression models in the literature. Compared to these traditional models, ANN models have been far less researched in the field of collision prediction modelling and model development. Historically, due to the high complexity and black-box nature of ANN models, researchers mainly opted for simpler approaches among the wide range of statistical regression techniques available. However, it is well-known that ANN could be extremely efficient in predictive models defined in complex engineering systems. With the current improvement in data acquisition systems in engineering fields, a large dataset analyses are becoming more common needs of people making managerial decisions. This has made sophisticated artificial intelligence (AI) models more popular due to their significant ability in dealing with large and complex datasets. Therefore, not only the field of transportation, but all engineering fields would be left with no choice other than to address emerging data analysis and technology needs. There is a growing necessity to insert these AI models to interpret the trend or futures states of certain parameters in transportation systems. This would enable practitioners to be better prepared for the design aspects of a system, or more importantly as in this thesis, they could add valuable information regarding an intervention of a change in current transportation systems being utilized by users. The prime objective of this thesis is defined around a comparative assessment between ARIMA model as a conventional statistical time series model in road safety and ANN models to point out the conclusive remarks of using ANN in this field. 32 Chapter 3: Autoregressive Integrated Moving Average Modelling 3.1 Overview This Chapter outlines the fundamentals of the ARIMA modelling for time series analysis. Section 3.2 discusses the subclasses and assumptions of ARIMA models. Section 3.3 presents the symbols of ARIMA models in general. Section 3.4 introduces the concepts of autocorrelation and partial autocorrelation functions. Section Error! Reference source not found. presents a conclusion for the chapter. 3.2 Introduction to ARIMA model Generally, the ARIMA model includes autoregressive (AR) terms, moving average (MA) terms, and differencing operations. The model also has several subclasses for each unique case, including:  AR for when a model only contains autoregressive terms,  MA for when a model only contains moving average terms,  ARMA for when a model contains both autoregressive and moving average terms,  SARIMA for when a model contains seasonal terms, and  ARIMAX for when a model contains other time series as the explanatory variables. The assumptions of the ARIMA model are as follows: 1) The error terms are representative of Gaussian white noise (ωt~N(0, σω2 )), meaning that the residuals of the model are independent and identically distributed (iid) with a mean of zero and homoscedastic variance. 33 2) The time series is stationary. A time series is considered to be (weakly or second-order) stationary if the following properties are satisfied [26]: 2.1) The mean is constant over time, 2.2) The variance is constant over time, and 2.3) The autocorrelation of each lag is constant over time. If the data exhibits nonstationary behaviour, three methods can be used to transform the data to be stationary data in order for the ARIMA model to be applied: 1) Differencing to account for non-stationarity resulting from a trend in the series. A first difference is applied to remove a linear trend and second difference is applied to remove a quadratic trend in the series, 2) Seasonal differencing to account for non-stationarity resulting from seasonality, and 3) Logarithm transformation to account for non-stationarity resulting from heteroscedastic variance over time. It should be noted that over-differencing may add unnecessary dependency to the model. For example, yt = ωt is white noise while yt = ωt − ωt−1 is a first order MA model. 3.3 General form of ARIMA models As mentioned previously, the ARIMA model may include different terms. The AR term for the variable yt is defined as the time lagged value. The autoregressive model with an order of p, indicating the most extreme lag of yt that is used as an independent variable to predict values at time t and denoted by AR(p), for a series with zero mean is given by: 34 yt = ø1yt−1 + ø2yt−2 + ⋯+ øpyt−p + ωt (3.1) where:  ø1, ø2, …, øp denote the coefficients of the AR orders and  ωt denotes the error term that should be Gaussian white noise. If the mean is nonzero, yt will be replaced by yt-µ and the AR(p) model is given by: yt = c + ø1yt−1 + ø2yt−2 + ⋯+ øpyt−p + ωt (3.2) where:  µ denotes the mean and  c denotes the constant of the model and is given by: c = µ(1 − ø1 − ⋯− øp) (3.3) The AR(p) model is also written by an autoregressive operator: ∅(B)yt = ωt (3.4) ∅p(B) = 1 − ø1B − ø2B2 − ⋯− øpBp (3.5) where:  Øp(B) denotes the autoregressive operator and  B denotes the back-shift operator, e.g. Bpyt = yt−p. The MA terms are the past errors of a time series model. The moving average model with the order of q, denoted by MA(q), is given by: yt = µ + θ1ωt−1 + θ2ωt−2 + ⋯+ θqωt−q + ωt (3.6) 35 where:  θ1, θ2, …, θq denote the coefficients of the MA orders. The MA(q) model is also written by a moving average operator: yt = µ + θq(B)ωt (3.7) θq(B) = 1 + θ1B + θ2B2 + ⋯+ θqBq (3.8) where:  θq(B) denotes the moving average operator. The ARMA model for a stationary series is given by: yt = ø1yt−1 + ø2yt−2 + ⋯+ øpyt−p + ωt + θ1ωt−1 + θ2ωt−2 + ⋯+ θqωt−q (3.9) The model can also be written in terms of the operators: ∅p(B)yt = θq(B)ωt (3.10) In most real cases, it can be assumed that a time series contains two components: a nonstationary d-th order polynomial trend and stationary series. Differencing with respect to the d-th order removes the nonstationary components of the series, so the remaining component will be stationary and the ARIMA model can be applied. A generalization of the ARIMA model can be written by: ∅p(B)dyt = θq(B)ωt (3.11) where:  ▽d denotes the difference operator, which is defined by: dyt = (1 − B)dyt (3.12) 36 If the mean is nonzero, the ARIMA model will be: ∅p(B)d(yt − µ) = θq(B)ωt (3.13) Seasonal behaviour would produce a series that is nonstationary because the average of certain distinct times inside of the seasonal span would be dissimilar to the averages of different times. Therefore, it is appropriate to consider seasonal terms in the ARIMA model. Seasonal span (S) is defined as the number of time points until the seasonal pattern repeats, e.g. S = 12 for monthly data or S = 4 for quarterly data. The seasonal form of the ARIMA model (SARIMA) is denoted by ARIMA(p,d,q) × (P,D,Q)S where p represents the non-seasonal AR order, d represents the non-seasonal differencing, q represents the non-seasonal MA order, P represents the seasonal AR (SAR) order, D represents the seasonal differencing, and Q represents the seasonal MA (SMA) order. This form of ARIMA model can be exhibited in the following Equation 3.14: ФP(BS)∅p(B)SDd(yt − µ) = ΘQ(BS)θq(B)ωt (3.14) where:  ΦP(BS) denotes the SAR operator, which is defined by: ΦP(BS) = 1 − Φ1BS − Φ2B2S − ⋯− ΦPBPS (3.15)  ΘQ(BS) denotes the SMA operator, which is defined by: ΘQ(BS) = 1 + Θ1BS + Θ2B2S + ⋯+ ΘQBQS (3.16)  SD denotes the seasonal difference operator, which is defined by: SD = (1 − BS)D (3.17) 37 3.4 Autocorrelation and partial autocorrelation function Autocovariance demonstrates the linear dependency between two observed points belonging to the same series at different times. The sample autocovariance at time lag k is given by: γ(k) =1n∑(yt−k − y)(yt − y̅)n−kt=1 (3.18) where:  n denotes the sample size and  y̅ denotes the mean value of time series. Autocorrelation measures the autocovariance of a time series by giving a value between [-1, 1]. The sample autocorrelation at lag k is equal to the sample autocovariance at lag k divided by the variance (σ2), which is given by: ρ(k) =γ(k)γ(0)=γ(k)σ2=∑ (yt−k − y̅)(yt − y̅)n−kt=1∑ (yt − y̅)2nt=1 (3.19) Autocorrelation function (ACF, a.k.a. correlogram in old references) is the plot of autocorrelation values within a time series at various lags versus time. ACF is an aid to time series analysis to check for: (1) stationarity, (2) seasonality, and (3) randomness in the residuals of a time series model. It can be proven that if xt represents Gaussian white noise and the number of samples is large, the standard deviation (SD) of the ACF is equal to σρ(k) = 1 √n⁄ [26]. According to the empirical rule (a.k.a. three-sigma rule), if the error of a model is Gaussian white noise, 95% of the ACF values should lie within the interval of two standard deviations (2σρ(k) = ±2 √n⁄ ). Hence, the autocorrelation for any lag is significant if its value extends outside of this interval. 38 Partial autocorrelation of a time series is a conditional correlation between yt and yt-k, if the relationship between intermediate lags, i.e. {yt-k+1, yt-k+2, …, yt-1}, is taken to consideration. For example, the partial autocorrelation between xt and xt-3 can be calculated from the autocorrelation between the residuals (ε1 and ε2) of the following two regressions: yt = c + ø1yt−1 + ø2yt−2 + ε1 (3.20) yt−3 = c + ø1yt−1 + ø2yt−2 + ε2 (3.21) The first order of a partial autocorrelation is equal to the first order of an autocorrelation as there is no lag in between. The sample partial autocorrelation for k ≥ 2 is denoted by: φ(k) = ρ(yt, yt−k|yt−k+1, yt−k+2, … , yt−1) (3.22) Partial autocorrelation function (PACF) coupled with ACF are used to determine the orders of AR and MA for a stationary time series. Theoretically, for an AR(p) model, the first p lags of the PACF are nonzero and the rest are zero. For a MA(q) model, the first q lags of the ACF are nonzero and the rest are zero [11]. However, the practical ACF and PACF patterns do not precisely match the theoretical patterns due to sampling error. 3.5 Summary In this chapter, the form and functionality of ARIMA models were described, based on an assimilation of the latest model development techniques found in the review of research literature. In Chapter 4, the basics of ANN modelling approach will be presented in order to develop a methodology for road safety applications. 39 Chapter 4: Artificial Neural Network Modelling 4.1 Overview This chapter outlines the fundamental concepts of ANNs. Section 4.2 describes similarities between biological and artificial neural networks. Sections 4.3 to 4.7 introduce the computational model of artificial neurons and how they are interconnected to form a network. Sections 4.8 to 4.10 explain the learning procedure and generalizing the network output. Section Error! Reference source not found. presents a conclusion for the chapter. 4.2 Biological neural network A biological neural network (e.g. the human brain) is a complex and intricate system, consisting of billions of interacting neuron cells. Each neuron cell is composed of a cell body with protruding notches, called dendrites, and a connection like wire called the axon. Neurons interact with one another via the dendrites (input) and the synaptic regions (output) on the tip of the axon, through process involving electric triggers. The electric triggers signal neurotransmitters to diffuse through the synaptic membrane and across the synaptic region to fuse with the dendrites of the corresponding neuron. This process is repeated from many other neurons, resulting in the interconnected neural network. The neurotransmitters are not equally weighted; some are more important than others and some act as a blockade to messages. These fundamental differences are based on the chemical makeup of the neurons and the substances surrounding and interacting with the neurons [69]. 40 The basis of this process is essential to ANNs. An artificial neuron receives multiple unequally weighted messages, comparable to the dendrites receiving neurotransmitters via diffusion in the biological neural network. The artificial neuron then processes and distributes the varying output, which is similar to the neurotransmitters travelling through the axon and pre-synaptic region. Artificial networks do not have the power of biological networks; however, they have many practical applications. 4.3 Artificial neuron ANN is composed of a collection of interconnected nodes (i.e. artificial neurons) that are a simple abstraction of biological neurons as the mechanism of biological neurons is not fully understood. Artificial neurons are referred to as computing units or processing elements in some references [70]. An artificial neuron (hereafter called neuron) is a processing element that is fed by a signal from each input (xi) with its strength represented by a weight (wi) in Figure 4.1. The weights and inputs can be respectively represented by the following matrices: W = [w1 w2 … wn] (4.1) X = [x1x2…xn] (4.2) 41 Figure 4.1 An artificial neuron Each input is multiplied by its corresponding weight. Also, an input with the value of one that operates similarly to an intercept in statistical regressions is multiplied by a bias (b, a.k.a. offset). In other words, a bias is equivalent to a weight on an extra input that always has a value of one. A neuron can be developed without a bias but networks with biases usually provide more satisfactory results. These bias and weighted inputs are added in summation (), then it is modified by a transfer function (f) that makes a scalar output (a): s = w1x1 + w2x2 + ⋯+ wnxn + b = ∑wixi + bni=1 (4.3) a = f(WX + b) (4.4)  f Transfer function s a Inputs xn x1 1 b w1 w2 wn x2 Summation Weights Output 42 The transfer function is fixed and selected before training the network. However, weights and biases are adjustable parameters of the neuron and updated by a learning algorithm until predefined goals are met [57]. 4.4 Transfer function A transfer (a.k.a. activation) function limits the amplitude of the neuron output. It can be selected by the modeler according to the problem that the network will solve. A transfer function can be a linear or nonlinear function of the weighted sum of inputs (s in Figure 4.1). The popular transfer functions and their graphs are respectively presented in Table 4.1 and Figure 4.2. Table 4.1 Transfer functions [3] Name Equation Derivative with respect to s Range Binary threshold f(s) = {0 s < θ1 s ≥ θ f(s)́ = 0 {0, 1} Linear f(s) = s f(s)́ = 1 (-∞, ∞) Rectified linear f(s) = {0 s < 0s s ≥ 0 f(s)́ = {0 s < 01 s ≥ 0 [0, ∞) Log-sigmoid f(s) =11 + e−s f(s)́ = f(s)(1 − f(s)) (0, 1) Hyperbolic tangent sigmoid f(s) = tanh(s) =es − e−ses + e−s f(s)́ = 1 − f(s)2 (-1, 1) 43 (a) Binary threshold function (b) Linear function (c) Rectified linear function (d) Log-sigmoid function (e) Hyperbolic tangent sigmoid function Figure 4.2 Transfer function graphs A binary threshold function sets the output as one if the weighted sum exceeds the threshold value, otherwise the output is zero. This function is appropriate for an operating classifier. A linear transfer function sends out the same value as it receives, which indicates that s = a. Rectified linear function sends out the same value if the weighted sum exceeds zero, otherwise the output is zero. A log-sigmoid function (a.k.a. Fermi function) gives a continuous output in the interval (0, 1) and is a smooth and bounded function of the total input. A hyperbolic tangent sigmoid function provides a symmetrical output value with a range of (-1, 1). Both sigmoid transfer functions are 44 used most often in practice, especially in any backpropagation learning algorithm as they are continuously differentiable and have simple derivative equations shown in Table 4.1. There are some transfer functions that are not the function of a neuron summation, which are out of the scope of this research. For example, a softmax function forces the outputs of all neurons in a layer to sum to one. Therefore, the network output can be a probability distribution across mutually discrete alternatives instead of across real values. 4.5 One-layer artificial neural network Multi-neurons often work in a parallel system that comprises a network layer. Figure 4.3 illustrates how every input variable is connected to all neurons in a layer. The calculation for each neuron is similar to the procedure presented in Section 4.3. The only differences are the definition of the weight matrix that is given by Equation 4.5 and the notation, where the first and second subscripts represent the input and neuron number associating together, respectively. For example, w1,2 indicates the weight connecting the first input to the second hidden neuron. W =[ w1,1 w1,2 … w1,qw2,1 w2,2 … w2,q............wn,1 wn,2 … wn,q] (4.5) 45 Figure 4.3 A layer of neurons 4.6 Artificial neural network layers ANN is organized into three types of layers: input, hidden, and output layer. The input layer contains no neurons; however, it is considered as a separate layer for simplicity. From a statistical viewpoint, the input and output layers respectively belong to independent and dependent variables. The major processing is done via the hidden layer(s). A network could have one or multiple hidden layers. ANNs with no hidden layers have limited applications, whereas adding hidden layer(s) transforms them into universal function approximators [57]. The number of inputs and outputs are obvious for most systems; however, the optimal number of hidden layers and neurons depend on the complexity of the problem that the network is positioned to solve. It should be noted that the optimum lag time of the inputs should also be determined in the DANN.  f Transfer function s1 a1 x1 x2 xn Summation b1  f s2 a2 b2  f sq aq bq Inputs Outputs w1,1 w2,1 wn,q Weights w1,2 46 4.7 Artificial neural network architecture ANN architecture defines the way neurons are connected. There are three general types of network architectures [57]: (1) feedforward network (Figure 4.4): The information flows in one direction through the input layer to the hidden layers until it reaches the output neuron(s), (2) recurrent (a.k.a. feedback) network (Figure 4.5): This network adds backward connections to a feedforward network so that the output of a layer is fed back to the input, and (3) competitive network: This network uses both feedforward and feedback networks and is explicitly applied for binary pattern recognition. This network architecture is not illustrated here as it is not used in this research. Each type of ANN architecture has several unique types (e.g. time delay and convolutional network are special types of feedforward network). Figure 4.4 Architecture of feedforward network with one hidden layer. Superscripts represent the layer number.  f 1 s11 a11 x1 x2 xn b11  f 1 a12 b12  f 1 s1q a1q1 b1qw11,1 w1n,q Input Layer Output Layer Hidden Layer  f 2 s21 a21 b21  f 2 s22 a22 b22  f 2 s2q a2q2 b2qw21,2w2n,q 47 Figure 4.5 Architecture of recurrent network. TDL represents tapped delay line, which will be explained in Section 6.3.3.2. The architecture of L-layer ANN is given by: N1 − N2 − ⋯− NL (4.6) where:  N1 denotes the number of inputs and  NL denotes the number of neurons in layer L. 4.8 Artificial neural network learning algorithm A network learning (a.k.a. training) algorithm is defined as the process of updating (a.k.a. optimizing) the adjustable parameters of a network (i.e. weights and biases) by minimizing the performance (a.k.a. cost) function of the network to estimate the rules that map the inputs to outputs [57]. Weights and biases should be initialised randomly before training the network, which will be explained in Section 5.4.2.2. A learning algorithm is broadly divided into three categories: (1) supervised learning: to predict an output when the desired output (hereafter called target) is presented for each input matrix, (2) x1(t)  f s at b w1 wi TDL x2(t) xi(t) wi+1 TDL TDL 48 reinforcement learning: to select an action to maximize the future rewards gained by a sequence of inputs, instead of presenting the correct output to the network, and (3) unsupervised learning: to discover a good internal representation of the inputs or to cluster them. One major aim of unsupervised learning is to create an internal representation of the input that will be used for the subsequent supervised learning. The learning type of time series analysis is unclear. On one hand, time series analysis is considered to be supervised learning since the output is provided to the network. On the other hand, it is unsupervised learning because a different output is not used and the network only finds the patterns inside the input (i.e. the best past and present values to predict the future). In other words, a supervised learning method is used while there is no need for a separate target dataset. There are two groups of supervised learning: (1) classification where the output is a class label (e.g. zero or one) and (2) regression where the target output is a real number or a matrix of real numbers. 4.8.1 Backpropagation Backpropagation is the learning algorithm to train a multilayer ANN. Perceptron and least mean square learning rules are designed to train single-layer networks; however, each layer might have different transfer functions and a different number of neurons in multilayer networks that prevent the error to be written explicitly as a function of the weights. Backpropagation is a generalized form of least mean square learning that allows the calculation of the implicit derivatives using the chain rule. Backpropagation is an approximate gradient (a.k.a. steepest) descent procedure that minimizes the squared error as the performance function. The term “error” here represents the discrepancy between each target and corresponding network output. It should be noted that the error is divided by two in order to remove the power when the derivatives are computed: 49 E =12(t − a)2 (4.7) where:  t denotes the target output and  a denotes the network output. Gradient descent is an iterative optimization technique that minimizes a function by: F(xk+1) = F(xk) − αF(xk) (4.8) where:  α denotes the learning rate (α < 1) and  F(xk) denotes the first derivative of F(x) estimated at point xk. Similarly, the updated weights in ANN are gained by: wijL(k + 1) = wijL(k) − α∂E∂wijL (4.9) where:  L denotes the number of the layer,  wijL(k + 1) denotes the updated weight gained at epoch k (i.e. one forward pass and one backward pass in ANN), and  ∂E∂wijL denotes the partial derivative of error function with respect to each weight. As the errors are the implicit function of the weights, the partial derivative of the error function to the weights are calculated by the chain rule: 50 ∂E∂wijL =∂E∂siL ∂siL∂wijL = ciL ajL−1 (4.10) ciL =∂fL(sjL)∂sjL wijL+1ciL+1 (4.11) where:  ∂fL(sjL)∂sjL denotes the derivative of the transfer function to the summation. It can be seen why differentiable transfer functions are applied for backpropagation. Equation 4.9 can be replaced by Equation 4.10: wijL(k + 1) = wijL(k) − αciLajL−1 (4.12) ciL is obtained in the last layer by: ciL = −∂fL(sjL)∂sjL (t − aiL) (4.13) The updated biases are also gained by replacing (∂siL∂bijL = 1) in Equation 4.10: bijL(k + 1) = bijL(k) − αciL (4.14) The reader is referred to [71] for the fundamental concept behind the backpropagation algorithm. Other types of optimization algorithms can also be used in backpropagation to minimize the network performance function such as Newton’s method and conjugate gradient, which are developed based on the second order of the Taylor series expansion. These algorithms are 51 appropriate for training networks with substantial number of parameters as they provide faster training than gradient descent. However, there is no guarantee for them to converge to the global minimum (i.e. xk is a global minimum point if F(xk) is smaller than F(xk + Δx) for all Δx ≠ 0) [57]. Neural network training is an iterative process where the network outputs are calculated through feedforward and the weights and biases are subsequently updated through backpropagation. The training process is terminated when one of the stopping criteria is met, which will be discussed in Section 5.4.2.3. 4.8.2 Artificial neural network parameters update The network parameters can be updated in three ways: 1) Stochastic gradient descent (a.k.a. incremental, sequential, or online learning): The parameters are updated after presenting every training sample to the network as shown in Equations 4.12 and 4.14. The term “stochastic” refers to the fact that it represents a rough estimate of the overall gradient. 2) Batch gradient descent (a.k.a. offline learning): The parameters are updated after applying the entire training samples to the network. The gradients (Jacobians in the Newton’s method and conjugate gradient) are calculated for each sample and averaged over all samples to estimate a more accurate gradient. Hence, Equations 4.7, 4.12, and 4.14 are respectively replaced by [57]: E =12Q∑(tq − aq)2Qq=1 (4.15) 52 wijL(k + 1) = wijL(k) −αQ∑ ciL(q)ajL−1(q)Qq=1 (4.16) bijL(k + 1) = bijL(k) −αQ∑ ciL(q)Qq=1 (4.17) where:  Q denotes the number of training sets. 3) Mini-batch gradient descent: It is similar to online learning while a subset of training samples is used each time to update the parameters. 4.9 Dataset division Datasets are widely divided into three subsets in neural networks: training, validating, and testing set. When datasets are time series, they should be partitioned into contiguous segments. A training set typically constitutes 70% of the original dataset, and validating and testing set, each 15% [57]. The training set is for estimating gradients, and updating weights and biases after each learning epoch. A validating set is not used in learning algorithm; however, it does prevent the network from overfitting by executing an early stop, which will be explained in Section 4.10. A testing set gives an unbiased estimate of the fitted model performance in predicting. Hence, it must not be used to train the network. In other words, training and validating sets are for network training, whereas a testing set demonstrates the performance accuracy of the trained network [72]. 4.10 Artificial neural network generalization ANN generalization is defined as how well the network performance is when new data is simulated by the trained network [57]. Three types of error might occur when the new data is predicted by 53 the network [73]: (1) overfitting (a.k.a. overlearning) meaning that the network indicates low training error but a high generalization error (i.e. error on the new data). This occurs when the model overestimates the capacity of the training set and is unlikely to perform as well as trained data over the new data, (2) underfitting meaning the network shows both high training and generalization error, and (3) extrapolation: The network fails to approximate the function for the new data because the training set does not represent the features of the new data, not due to overfitting. References are often discussed in overfitting because extrapolation error is impossible to prevent in ANNs. Conversely, underfitting is easy to detect and the solution is to attempt to use an alternate network configuration. Overfitting occurs because each training dataset includes the regularities in the mapping from input to output as well as sampling error. There are many ways to prevent overfitting and improve the generalization of ANNs. The most widely used ones are to: 1) Collect more data samples: Overfitting regularly happens when the amount of training data is limited. 2) Combine multiple neural networks: Averaging the predictions of many models is an effective way to tradeoff between fitting the true regularities and overfitting sampling error in the data. Moreover, a neural network might converge to a local minimum point on the error surface with an initial condition instead of converging to the global minimum with another initial condition. If the individual networks provide different results, averaging over the testing sets of the combined networks yield more satisfactory results than all the individual networks. The most frequently used methods to combine multiple neural networks are: 54 2.1) Try different initial weights and biases: Each network is trained different times with various initial values of the parameters. 2.2) Bootstrap aggregating (bagging): This method introduced by Breiman [74] to improve classification by averaging various networks that are trained on different subsets of the training data. 2.3) k-fold cross-validation (a.k.a. rotation estimation): Like bagging, the aim of this method is to combine the results of different networks to derive more accurate and precise network performance. The advantage of this method over bagging is that each observation is used only once for validating set. In this method, the testing set is first put aside from dataset. Then, the rest is divided into k-equal sized folds (usually 10). Each time the network is trained, the kth fold is considered to be the validating set and the remaining k-1 folds are used as the training set. Finally, the results from k trained networks are averaged to produce a single estimation of testing set [75]. 3) Use a network that has the right capacity: This method aims to find a simplest network that is complicated enough to obtain the right regularities considering the limited and noisy available data. There are generally five approaches to limit the capacity of the network. The first three approaches are related to the network configuration by constraining the number of parameters (equivalently, the number of hidden layers and neurons) and the last two control the magnitude of parameters [57]: 3.1) Growing: Start with the simplest network with no neurons, then add more neurons until the network starts to overfit. 3.2) Pruning: Remove the neurons one by one from a large network including many neurons until the network does an adequate job. 55 3.3) Global searches: Search among all possible networks to find the simplest network with satisfactory performance. 3.4) Early stopping: The complexity of the network increases after each learning iteration. This approach prevents overfitting by stopping the iteration when the error on the validating set grows for several consecutive iterations. Then, the weights and biases for the iteration that exhibits the best performance (i.e. lowest error) will be used as the trained network. In other words, the error of the validating set rises when the network starts to overfit. The assumption of this approach is to find the true regularities before the network finds the spurious regularities belonging to the sampling error of a specific training set. 3.5) Regularization: This approach adds the sum of the squared network weights to the error on the validating set to penalize the model complexity: E = β∑(ti − ai)2Vi=1+ α∑wj2nj=1 (4.18) where:  V denotes the number of validating samples,  n denotes the number of network weights, and  α and β denote the regularization parameters and ratio α β⁄ controls the complexity of the ANN. A combination of these methods is usually used to prevent overfitting and testing set represents the generalization capability of the model in future. An inferior method is to try a group of candidate networks and see which gives the best performance on the testing set because there is 56 no guarantee that the chosen network will have the same performance on another testing set. Hence, a testing set should be never utilized to train networks, design network architectures, or select a network among some network alternatives. 4.11 Summary In this chapter, the fundamentals of ANNs were presented. The ANN is a complex computational process, and as such, this chapter aims to provide clarity on the computational model of ANNs and how the components interconnect to form a network, the process of learning that the network undergoes. The principles provided in this chapter will be applied to develop a methodology for predicting a collision series in the subsequent chapter. 57 Chapter 5: Model Development 5.1 Overview This chapter presents the model development methodology. Section 5.2 describes the methodology of designing an ARIMA model. Section 5.3 describes the methodology of developing an ARIMA with explanatory variables (ARIMAX) model. Section 5.4 explains the methodology of training a DANN model. Section Error! Reference source not found. presents a conclusion for the chapter. 5.2 Methodology of ARIMA model Figure 5.1 presents the procedure for building the ARIMA model, including four general steps: model identification (steps 1 to 3), model estimation (step 4), model diagnostic checking (steps 5 and 6), and predicting (step 7). 58 Figure 5.1 Flowchart of building ARIMA model Step 1) Plot time series: Time series plot of the observed data to reveal if the data contains: Select the best accepted model Check if time series is stationary using the plots Yes Yes No Plot time series data Plot ACF of the raw data Identify the orders of (S)AR and (S)MA using ACF and PACF of the stationary data Identify the order of (seasonal) differencing Estimate the model coefficients If the series is stationary Make predictions If not satisfied If the residuals are white noise 59 1.1) Any repeating seasonal pattern that depends on the nature of the series (i.e. whether the data was aggregated hourly, daily, weekly, monthly, or quarterly): If the series indicates seasonal behaviour, the SARIMA model must be developed. It should be noted that the yearly aggregated time series does not need seasonal terms. 1.2) Any consistent downward or upward trend over time: For a linear or quadratic trend, the series would need one or two differences, respectively. 1.3) Inconstant variance over time: Logarithm transformation would eliminate the increasing variance through time. Other models, like the autoregressive conditionally heteroscedastic (ARCH) model, should be used for a series with volatile variance. 1.4) Any obvious outliers: The outliers should be examined from the series. 1.5) Any abrupt changes: Abrupt changes may indicate a need for more investigation. Step 2) Plot the ACF of the raw data: ACF is very helpful in identifying the non-seasonal or seasonal differencing order. When the time series displays seasonal behaviour, with S seasonal span, seasonal differencing is necessary if there is a slow decay in the ACF at lags kS (k = 1, 2, …). Non-seasonal differencing is also needed if the ACF stays close to 1 or slowly decays over many lags. Also, some statistical tests may reveal whether or not a series is stationary when it is not obvious from the plots. The most popular test is the Augmented Dickey-Fully (ADF) test, which determines that the null hypothesis will be satisfied if the data is nonstationary and has a unit root in the AR model, whereas the alternative hypothesis will be satisfied if the data is stationary. Step 3) Plot the ACF and PACF of stationary series: 60 The ACF and PACF of the differenced (if necessary) series are used together to determine the tentative model(s). Table 5.1 and Table 5.2 describe the behaviour of the ACF and PACF for different non-seasonal and seasonal ARIMA models, respectively. It should be noted that in some cases, it is challenging to determine the orders of (S)AR and (S)MA. Table 5.1 Attribute of autocorrelation and partial autocorrelation function for non-seasonal ARIMA model [26] Plot AR(p) MA(q) ARMA(p,q) Autocorrelation function Tails off Cuts off at lag q Tails off Partial autocorrelation function Cuts off at lag p Tails off Tails off Table 5.2 Attribute of autocorrelation and partial autocorrelation function for seasonal ARIMA model [26] Plot SAR(P) SMA(Q) SARMA(P,Q) Autocorrelation function Tails off at lags kS Cuts off at lag QS Tails off Partial autocorrelation function Cuts off at lag PS Tails off at lags kS Tails off Step 4) Estimate the coefficients: The coefficients of the tentative model(s) are estimated by the maximum likelihood method. Step 5) Diagnostic checking: 5.1) The significance of the coefficients: The constant, AR, or MA term should be removed from the model if its coefficient is insignificant. 5.2) Analysis of the residuals: The residuals of the model should be considered to be Gaussian white noise. To check the accuracy of this assumption, static tests should be performed. First, 61 the standardized residuals should represent no pattern with a mean of zero and homeostatic variance. The standardized residual for an observation at time t is calculated by: Standardized residual at time t =yt − ŷtσω (5.1) where:  ŷt is the predicted value of yt by the model and  σω is the standard deviation of the residuals. Second, the ACF of the residuals belonging to a good fitted model falls approximately within the given confidence limits. Due to sampling error, a significant lag at an unusual point may occur, which can be ignored. Third, the overall autocorrelation of the residuals for a certain number of lags (h) should be insignificant. The Ljung-Box (a.k.a. modified Pierce-Box) test is a diagnostic statistic test used to examine if the accumulated autocorrelations differ from zero. The residuals of the model for the first h lags are independently distributed when the p-value of this test is not significant and obtained by [76]: Q(H) = n(n + 2)∑ρ(k)2n − khk=1 (5.2) where:  n denotes the number of observations left after differencing. 5.3) Normal distribution: The residuals should be normally distributed. Different model(s) should be selected in the event that any assumption is rejected. Step 6) Compare models: 62 Provided that more than one tentative model yields satisfactory results, the model with the fewest number of coefficients is preferred. If the models have an equal number of coefficients, the smaller performance measure values indicate the better model: (1) mean squared error (MSE), (2) Akaike information criteria (AIC), and (3) Bayesian information criteria (BIC). MSE measures the overall performance of the model while AIC and BIC consider not only goodness of fit but also complexity of the model. Their formulas are as given by: MSE =1n∑(yt − ŷt)2 (5.3) AIC = log (∑(yt − ŷt)2n) +n + 2mn (5.4) BIC = log (∑(yt − ŷt)2n) +mlog(n)n (5.5) where:  m denotes the number of coefficients. Step 7) Predicting: The first future value of a time series is calculated from the developed ARIMA model using the last observed values. The next future value is obtained by using the first predicted value. This pattern of value determination continues. The standard error of the predicted error with any lead time (l) is [26]: σ(l) = σω2 (∑ ψj2l−1j=0)12⁄ (5.6) where: 63  ψj denotes the psi-weight and obtained by converting the ARIMA model to an infinite order of the MA model: xt − µ = (1 + ψ1B + ψ2B2 + ⋯)ωt = ∑ψjωt−j∞j=0 (5.7) When normally distributed errors are assumed, the predicted value with a 95% probability limit will lie within: yt+l = ŷt+l ± 1.96σ(l) (5.8) If the lower limit of any prediction interval is negative, which is impossible for collision counts, the lower limits assign zero for those points in practice. It should be mentioned that as the prediction goes farther in time, the standard error increases and the prediction values converge to the mean and remain constant for a stationary series. 5.3 Methodology of ARIMAX modelling An autoregressive integrated moving average with explanatory variables (ARIMAX) model is an extended version of the ARIMA model including explanatory variables that explains the relationship among several time series [77]. The ARIMAX model is similar to a linear lagged regression; however, it allows for the advantages of the ARIMA model in the presence of autocorrelation among residuals to still be utilized. Cross-correlation function (CCF) shows the sample cross-correlations between yt and various lags of xt and can be employed to identify lag times of xt that might be useful to predict the output (yt) in ARIMAX models. The sample cross-correlation between two series of xt and yt at lag time k is given by: 64 ρxy(k) =γxy(k)√γx(0)γy(0) (5.9) where:  γxy denotes the sample cross-covariance between two series of xt and yt at lag time k and is calculated by: γxy(k) =1n∑(xt−k − x̅)(yt − y̅)n−kt=1 (5.10) The CCF may produce a misleading or unclear interpretation of the relationship between the input and output due to common trends that exist between two series over time. To make the CCF easily identifiable, the following pre-whitening procedure should be developed: Step 1) Transform the input series to white noise by determining the residual series of a possibly approximate ARIMA model. For simplicity, only an AR model with required differencing is considered. Step 2) Filter the output series by calculating the differences between the observed values of yt and its estimated values using the developed ARIMA model with the same coefficients from step 1. Step 3) Identify the possible lagged terms from the CCF of the residuals of xt from step 1 and filtered values of yt from step 2. Step 4) Develop a lagged regression with the original series and calculate the residuals. Step 5) Identify the possible ARIMA model for the residuals obtained from step 4. Step 6) Develop an ARIMAX model using the possible lagged terms from step 3 and ARIMA terms from step 4. 65 Step 7) If any explanatory variable coefficients are insignificant, they should be excluded from the model and steps 4 to 6 should be repeated without the insignificant coefficients. Step 8) If any diagnostic tests corresponding to the ARIMA model are not satisfied, steps 5 and 6 should be repeated. Step 9) Predict the future values using the fitted ARIMAX model and Equation 5.8. 5.4 Methodology of dynamic artificial neural network modelling Chapter 4 has focused on ANN architecture and learning process, with an emphasis on the fundamental concept behind ANNs. In this section, theoretical and practical aspects are combined to explain the design procedure of DANN modelling. The entire DANN design is an iterative process and stays until satisfactory performance is produced. Figure 5.2 illustrates the flowchart of ANN design process inclusively for dynamic networks and predictive models, which can be divided into three pre-training, training, and post-training levels. The following sections of this chapter will discuss the steps of each level. 66 Figure 5.2 Flowchart of training dynamic artificial neural network (DANN)Pre-training level Training level Post-training level Network performance is satisfactory End Yes Collect data Preprocess data Select network architecture Select optimum lag of inputs Select optimum number of hidden layers and neurons Select learning algorithm Initialize network parameters Select stopping criteria Train network Denormalize network output Analyze network performance Try different network No 67 5.4.1 Pre-training steps Several steps are required before training the network in order to speed up the learning process and increase the network performance. The pre-training steps consist of data collection, data preprocessing, choice of network architecture, choice of network specifics (i.e. optimum number of hidden layers and neurons), and choice of optimum lag of inputs. 5.4.1.1 Data collection The datasets should be sufficiently sampled before the network is trained. However, there is no control over the data collection process in most cases related to road safety analysis. That is, collisions are rare events, and as such, collision data is often scarce, which requires the entire set to be used as opposed to just a sample of the collision data set. It is, therefore, not discussed further. 5.4.1.2 Data preprocessing Data preprocessing refers to any type of preliminary processing technique performed on raw data to ensure that it is understandable training set for the model. Three data preprocessing methods of outlier detection, wavelet transform, and normalization were used in this research. 5.4.1.2.1 Outlier detection An outlier is an unusual or non-repetitive observation that influences time series data and may have a substantial effect on the accuracy prediction of the model [78]. Hampel filtering is a common method to detect outliers in time series analysis. This filter is based on the moving data window: 68 WtK = {xt−K, xt−K+1, … , xt, … , xt+K−1, xt+K} (5.11) where:  K denotes the window half-width and is a positive integer. The Hampel filter was only applied on the training data. As no points can be removed from time series, this filter replaced the central value in the data window with the median of the data window if it was far enough from the median. The filter response can be written by [79]: yt = {xt |xt − mt| ≤ Tσtmt |xt − mt| > Tσt (5.12) where:  yt denotes the filtered time series,  mt denotes the median of the data window: mt = median{xt−K, xt−K+1, … , xt, … , xt+K−1, xt+K} (5.13)  σt denotes the standard deviation, which is calculated by the scaled median absolute deviation: σt = 1.4826 median(|xt−K − mt|, … , |xt+K − mt|) (5.14)  T denotes the threshold parameter, which is usually set to three. 5.4.1.2.2 Wavelet transform Many studies have demonstrated that a network trained on a time series preprocessed by wavelet transform analysis into its subseries produces more satisfactory results than a network on raw data [80]. Wavelet is a waveform with an average of zero and a limited duration. Wavelet transform decomposes the original time series into scaled (i.e. stretched or compressed) and shifted (i.e. 69 delayed or hastened) versions of the original function called mother wavelet to extract various features of the series. In other words, wavelet transform offers a time-scale representation of the series [81]. Primarily, wavelet transforms are divided into two groups: continuous or discrete. Continuous wavelet transform1 (CWT) is the sum of the shifted and scaled versions of the mother wavelet over time [82]: CWT(a, b) =1√a∫ xt+∞−∞ψ(t − ba)dt (5.15) where:  CWT(a, b) denotes the wavelet coefficients, which are a function of scale and position,  a denotes the scale (a.k.a. dilation) factor (a ≠ 0),  b denotes the shift (a.k.a. translation or position) factor, and  ψ(t) denotes the mother wavelet. CWT is computationally impossible because of the inability to analyze an infinite number of scaled and shifted mother wavelets, as scales and positions are real numbers (a, b ∈ ℝ). In contrast, discrete wavelet transform (DWT) is practical to analyze as the wavelet factors are always a power of two known as dyadic scales and positions (a = 2j and b = ka, for (j, k) ∈ ℤ) [83]: 1 CWT is inherently discrete because time series data is measured at discrete time intervals. The CWT is continuous in terms of scale and shift, which can obtain any numbers [81]. 70 DWT(j, n) =1√2j∑ xtψ(t − k2j2j)T−1t=0 (5.16) where:  t denotes the discrete time and  T denotes the length of time series xt. Mallat [84] proposed a practical algorithm to implement DWT by passing the original time series through two (low-pass and high-pass) filters to obtain a series of approximations (A) and details (D), respectively. The approximations contain low-frequency components, while the details contain high-frequency components of the original series. The decomposition process can be iterated by breaking down the approximations into higher resolution levels, which is called the wavelet decomposition tree (Figure 5.3). The suitable number of decomposition levels can be selected by [85]: L = int[log(n)] (5.17) where:  L denotes the decomposition level,  n denotes the number of data points, and  int stands for integer number. 71 Figure 5.3 Wavelet decomposition tree. A and D represent the approximations and details, respectively. The subscript represents the level of decomposition. The original time series can be reconstructed from the components without loss of information, for example after Nth decomposition level: Original time series = AN + D1 + D2 + ⋯+ DN (5.18) where:  The subscript denotes the level of approximations or details. 5.4.1.2.3 Data normalization Data normalization is the process of transforming the input matrix and target output to have a unit variance before training the network to prevent [86]: (1) overriding the inputs with large values and (2) early saturation of hidden neurons. Saturation occurs when a neuron virtually reaches a boundary of the nonlinear transfer function and the training process becomes very slow. Generally, data is normalized based on the range of a nonlinear transfer function used in the hidden layer by: xin =(fmax − fmin)  (xi − xmin)xmax − xmin+ fmin (5.19) where: Original time series D1 A1 D2 A2 D3 A3 ... 72  xin denotes the normalized value of xi,  xmin and xmax respectively denote the minimum and maximum of training set, and  fmin and fmax respectively denote the minimum and maximum of the range of nonlinear transfer function. Specifically, the normalized equations for the log-sigmoid and tangent hyperbolic sigmoid functions are respectively: xin =xi − xminxmax − xmin (5.20) xin = 2xi − xminxmax − xmin− 1 (5.21) After the network is trained, the reverse normalization is necessary to convert the network output to the original unit of the target output. It should be noted that validating and testing data are normalized using summary statistics of the training data. Otherwise, the network results will be biased. 5.4.1.3 Choice of network architecture The network architecture is selected based on the type of problem that the network is supposed to solve, e.g. fitting, pattern recognition, prediction, or clustering. According to the data type, ANN is classified into either a: (1) static network where the network output is calculated directly from the concurrent values of the input vector or (2) dynamic network where the output is regressed on the previous values of the output, which is known as a nonlinear autoregressive (NAR) network or previous values of explanatory inputs as well, which is known as a nonlinear autoregressive 73 with external inputs (NARX) network. The standard neural network architecture for predicting a time series is a dynamic network. 5.4.1.4 Choice of network architecture specifics After the network architecture has been selected, the optimum number of hidden layers and neurons should be determined. Since no systematic way has been discovered to select the optimal number of hidden layers and neurons, a trial-and-error method is typically used. Thus, it is usually unknown before network training. 5.4.1.5 Choice of optimum lag of inputs The choice of optimum lag of inputs is a controversial issue in parametric, semi-parametric, and especially data-driven time series modelling approaches. On the one hand, the redundant dimension of the lags in the input vector would cause the network training process to be longer and more complicated. On the other hand, if the input vector contains irrelevant inputs, the network will yield an unsatisfactory performance [87]. Algorithms to find the optimum lag time of inputs are viewed to be either a model-based (known as wrapper) or model-free (known as filter) algorithm. Wrapper algorithms are considered as part of network architecture optimisation. Wrappers search through all or a subset of possible input lag vectors and select the one that provides the best generalization performance on a validating set. Any network performance measures can correspond to find the optimal lag of inputs vector but the most conventional criterion is MSE [88]. Wrappers can be implemented in three ways: (1) feedforward selection, (2) backward elimination, and (3) global search. A feedforward selection 74 starts by training the network with a single lag then the lags are added one by one. The algorithm stops when more lag is failed to improve the network performance. A backward elimination is the reverse of feedforward algorithm. ANN is trained with a substantial number of lag inputs. The relative unimportant lag inputs are removed one at a time until eliminating any lags fails to improve the network performance. Compared to feedforward selection, this method is more inefficient because a large network may need to be developed and it may be difficult to find the relevant input lags among many candidates. A global search evaluates many combinations together to select the optimal lag inputs vector. Generally, global optimization is classified into two major groups of heuristic and exhaustive search. Heuristic search is suitable when the number of available input variables is large. An example of this strategy is genetic algorithm. Exhaustive search is feasible for a relatively small number of available input series and assesses all possible input lags and selects according to the predefined optimally measure. It is the only strategy that is promised to select the best lag inputs for an ANN [89]. Filter algorithms measure the dependency between input lags and the desired output [87], [88]. Cross-correlation and autocorrelation were applied to identify the appropriate lags [90]. However, linear-based correlation techniques used in traditional time series analyses to select input variables fail to reveal the nonlinearity of a dynamic system. Drezga and Rahman [91] proposed a phase-space embedding method for input variable selection [92]. This method is developed based on the Taken’s theorem that an observed time series can be reconstructed into a series of lag-time vectors: xt = f{xt, xt−τ, … , xt−(m−1)τ} where τ and m denote the lag time and embedding dimension, respectively. 75 5.4.2 Training steps After preparing the data and selecting network architecture, the network can be trained. The network training process has various steps, including choice of learning algorithm, network parameters initialization, and choice of stopping criteria for training process. 5.4.2.1 Choice of learning algorithm Jacobian-based optimization algorithms (e.g. Newton’s method or conjugate gradient) are appropriate for training networks with a substantial number of parameters (i.e. over a thousand) as they provide faster training than gradient descent. However, there is no guarantee for them to converge to the global minimum (i.e. xk is a global minimum point if F(xk) is smaller than F(xk + Δx) for all Δx ≠ 0) [57]. The learning algorithm can be performed in either a stochastic or batch mode. Stochastic learning is faster especially when the number of training samples are very large. However, full-batch mode produces smaller error in terms of residuals. Also, the performance function might oscillate and never converge to a minimum point in stochastic mode whereas full-batch is more stable. The stochastic or mini-batch learning is preferable to speed up the network training when the training set is vey large and redundant. 5.4.2.2 Network parameters initialization Weights and biases should be initialized to start training the network. The initial weights of two neurons in the same layer cannot be equal because both neurons will have the same gradients and consequently learn the same features. Moreover, if the initial weights are set with large 76 magnitudes, the training process of the neurons that have transfer functions with finite boundaries like sigmoid families will become very slow as the neurons become saturated. Nguyen and Widrow [93] developed an algorithm to initialize the weights and biases of the hidden neurons to speed up the training process in backpropagation. The Nguyen-Widrow method is the most widely used initialization algorithm in ANNs and its assumption is that each element of the input vector is scaled in the range of [-1, 1]. In this method, the magnitude of each column of the weight matrix given by Equation 4.5 are randomly set to have the Euclidean length of: ‖wq‖ = 0.7L1 R⁄ (5.22) where:  ‖wq‖ denotes the Euclidean length of the weights in row q,  L denotes the number neurons in the hidden layer, and  R denotes the number of inputs connected to the hidden layer. The magnitude of bias corresponding to hidden neuron q are also randomly set in the range of −‖wq‖ and ‖wq‖. 5.4.2.3 Stopping criteria Some criteria should be defined to stop training as the performance error of ANNs never converges to exactly zero for most cases. The training procedure usually stops when one of the following criteria is met: 1) Minimum level of training error, 2) Maximum number of epochs, 3) Maximum amount of training time, 77 4) Minimum gradient of the performance function, 5) Minimum reduction in performance function, or 6) After consecutive number of iterations that the validating performance is failed to improve if the early stopping method is implemented. 5.4.3 Post-training steps The promising success of a trained network should be tested before using it. The post-training methods for DANNs include: network output denormalization and choice of network performance measure. 5.4.3.1 Network output denormalization The original network output can be denormalized as follows [57]: ai = (ain + 1) (tmax − tmin2) + tmin (5.23) where:  ai denotes the denormalized network output,  ain denotes the original network output, and  tmin and tmax denote the minimum and maximum of target sets, respectively. 5.4.3.2 Choice of network performance measure Various performance measures have been proposed in the ANN literature to compare the predictive accuracy of different trained ANNs. The most popular performance measures and their properties are as follows [94]: 78 1) Mean absolute error (MAE) measures the absolute average deviation of predicted values from targets. It provides a scale of the total error because the errors with opposite signs do not cancel each other out. MAE is given by: MAE =1n∑|et|nt=1 (5.24) where:  n denotes the number of samples and  et denotes the sum of discrepancies between target and corresponding network output at time t and given by: et = tt − at (5.25) 2) Mean absolute percentage error (MAPE) measures the percentage of average absolute error. Similarity to MAE, it does not show the direction of the error. The only difference with MAE is that MAPE provides the percentage of the error. MAPE is given by: MAPE =1n∑|ettt|nt=1× 100 (5.26) 3) Mean percentage error (MPE) measures the percentage of the average error. In contrast to MAPE, it shows the direction of the error. The zero value does not indicate a perfect prediction because the errors with opposite signs cancel each other out. MPE is given by: MPE =1n∑(ettt) × 100nt=1 (5.27) 4) Mean squared error (MSE) measures the average squared deviation of the prediction error and is the most widely adopted criterion. It does not show the direction of the error but represents the scale of the total error as the positive and negative error values do not cancel each other 79 out. In contrast to the previous performance measures, it panelizes extreme errors. MSE is given by: MSE =1n∑et2nt=1 (5.28) 5) Root mean squared error (RMSE) is the square root of calculated MSE. RMSE is given by: RMSE = √1n∑et2nt=1 (5.29) The preferred model is the one that corresponds to the minimum performance measure. It is better to use more than one performance measure to judge the networks, as every measure has its own properties. The above performance measures compare the relative performance of different ANNs while the correlation of determination (denoted as R2) can be used to evaluate the performance of ANNs individually. The R2 value measures the degree of correlation among the observations and predicted values by the model, which is given by [95]: R2 = 1 −∑ et2nt=1∑ tt2nt=1 −1n∑ at2nt=1 (5.30) 5.5 Summary In this chapter, the procedure to develop the ARIMA and DANN models was solidified. The ideas presented in this chapter will guide the assembly of the univariate and multivariate ARIMA and DANN models with regard to road collisions in the subsequent chapter. 80 Chapter 6: Results and Discussion 6.1 Overview This chapter presents and discusses results based on the proposed methodology given in Chapter 5. Section 6.2 outlines data collection resources and summary statistics of the collision series for reporting a quantitative description. Sections 6.3.1 and 6.3.2 present the ARIMA and ARIMAX models, and their diagnostic checking of the residuals. Section 6.3.3 provides the results of DANN models. Section 6.3.4 compares the performance of the models yielded in previous sections. Section Error! Reference source not found. evaluates the impact of safety management phases implemented on Springfield Road, Kelowna by analyzing collision counts. Section Error! Reference source not found. summarizes the findings of this chapter. 6.2 Assemble datasets Two distinct case studies were examined in this research: (1) City of Kelowna, where the performance of DANN can be compared with ARIMA approach (presenting in Sections 6.3.1 to 6.3.4) and (2) Springfield Road, where the traditional time series analyses in road safety studies were not applicable because the collision series violated the normal distribution assumption (presenting in Section Error! Reference source not found.). All datasets were time series that comprised 162 monthly aggregated measures covering the period between January 2003 and June 2016. Five types of datasets (collision claims, exposure data, weather data, economic conditions, and safety interventions) were assembled from various sources. As summarized in Table 6.1, the collision claims and exposure data were received from ICBC. Collision claims included the date, time, location, and severity type of each collision. Two variables of number of actively insured 81 vehicles and active driver licenses were considered as two alternative proxies of exposure data. Five weather indices of mean minimum, mean and mean maximum temperature, total precipitation, and number of days with precipitation were collected from the website of Environment and Climate Change Canada (https://www.canada.ca/en/environment-climate-change.html). To represent economic conditions, monthly employment and unemployment number in the Province of British Columbia were obtained from the website of Statistics Canada (http://www.statcan.gc.ca/eng/start). The type and implementation date of safety interventions on Springfield Road were acquired from the City of Kelowna. 6.3 Case study one: City of Kelowna Three collision types were the dependent variables of this case study, including: (1) total (TOTAL), (2) property damage only (PDO), and (3) killed or seriously injured (KSI) collisions. As depicted in Table 6.1, for all series, the highest collision counts occurred in 2015, excluding 2016 that was averaged over the first six months of the year. However, the year 2008 had the highest standard deviation (SD). One probable reason is a substantial unexpected number of collisions took place in December 2008. Table 6.2 summarizes descriptive statistics for the candidate explanatory variables. 82 Table 6.1 Summary statistics for monthly number of collisions in the City of Kelowna Year Collision type TOTAL PDO KSI Mean* Standard deviation* Mean Standard deviation Mean Standard deviation 2003 524 72 394 56 130 21 2004 578 57 447 48 131 15 2005 628 65 486 51 142 23 2006 683 118 538 99 145 26 2007 726 93 578 73 148 26 2008 751 141 599 117 152 29 2009 718 51 572 45 146 15 2010 671 72 520 53 151 29 2011 639 61 499 45 140 20 2012 653 70 506 56 147 21 2013 687 78 535 64 153 22 2014 718 82 563 64 156 23 2015 795 85 637 68 157 24 2016** 830 60 672 51 158 20 Total 686 79 539 64 147 23 *Mean and standard deviation are calculated over 12 months for a given year. **These values include 6 months of the year (i.e. from January to June). 83 Table 6.2 Summary statistics for explanatory variables. Standard deviation is in italics under the mean. Year Exposure data Weather data Economic conditions AIV* ADL* MINT* MT* MAXT* TP* NDP* EMP* UEMP* 2003 71982 77592 23.3 9.9 -3.4 23.7 9.5 1552.0 173.5 2454 237 12.5 8.8 7.6 14.9 4.6 49.3 16.4 2004 74888 78791 22.3 9.2 -4.1 39.5 12.0 1603.3 157.6 2981 564 11.2 8.3 8.8 24.4 4.1 50.6 20.0 2005 77979 80783 21.9 7.8 -6.1 27.7 19.8 1659.2 130.0 2983 925 10.1 8.1 9.8 14.2 6.6 42.5 23.1 2006 81051 83236 23.5 2.4 -4.4 28.4 10.8 1702.4 107.2 3123 882 11.3 0.3 9.6 14.3 4.4 46.6 11.8 2007 85044 85789 22.4 2.2 -4.1 33.1 10.1 1756.7 98.2 3543 868 11.0 0.0 8.9 24.9 4.4 44.0 7.4 2008 87807 88241 22.9 2.3 -5.7 23.6 9.9 1786.3 107.6 3661 721 11.5 0.0 10.4 14.7 4.4 45.6 6.8 2009 88264 89852 21.7 2.5 -6.6 21.8 7.3 1712.0 183.1 3594 454 12.6 0.0 10.8 10.5 3.3 34.2 12.2 2010 89778 91373 22.1 2.6 -6.0 18.7 8.9 1723.4 182.1 3688 558 9.9 0.0 8.5 11.9 3.1 38.5 15.4 2011 90574 92425 20.8 2.7 -5.6 19.6 9.1 1724.8 181.2 3691 411 9.6 0.0 7.9 14.9 3.9 43.6 21.3 2012 91886 93648 22.0 2.8 -5.3 29.1 9.6 1769.5 166.1 4025 490 9.8 0.0 7.6 24.1 5.5 48.4 10.8 2013 94077 95213 21.3 2.9 -5.3 32.0 9.9 1775.8 159.7 4120 941 11.9 0.0 8.4 24.6 3.9 49.4 9.6 2014 96676 97333 22.6 8.5 -6.0 26.7 9.3 1781.9 147.0 4202 567 10.2 9.0 10.3 11.1 3.0 46.5 13.6 2015 100113 99325 23.7 9.6 -4.9 22.7 6.9 1824.0 151.4 4259 1031 11.2 8.3 7.8 10.2 1.9 51.3 6.5 2016** 101618 101101 22.1 8.2 -6.8 21.6 6.8 1836.1 159.2 5445 358 10.5 7.2 8.2 16.8 4.4 49.7 17.3 Total 87572 89720 22.3 5.1 -5.2 26.6 10.2 1725.1 149.9 9104 6820 10.6 6.0 8.7 17.7 5.1 88.1 31.7 *AIV: Actively insured vehicles; ADL: Active driver licenses; MINT: Mean minimum temperature; MT: Mean temperature; MAXT: Mean maximum temperature; TP: Total precipitation; NDP: Number of days with precipitation; EMP: Employment number ( 1,000); UEMP: Unemployment number ( 1,000). **These values include 6 months of the year (i.e. from January to June). 84 Figure 6.1 shows time series plots of the collision series in the City of Kelowna. Regarding the TOTAL and PDO collisions, the mean was not constant over time. An upward trend was observed between 2003 and 2008, followed by a noticeable decrease between 2009 and 2011, and again rising until the end of the period. The reduction in collision counts between 2008 and 2011 may have been due to the economic crisis that hit in 2008. The KSI collisions illustrated similar behavior, but the slope changed slightly. One possible outlier can be seen at December 2008 for all collision series. 85 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.1 Time series plots of monthly number of collisions in the City of Kelowna (January 2003 to June 2016). One possible outlier can be seen at December 2008 for all collision series. 3004005006007008009001000110012002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of TOTAL collisions20030040050060070080090010002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of PDO collisions70901101301501701902102302003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of KSI collisionsTime (year)86 Figure 6.2 shows the seasonal subseries plot2 of each collision series from 2003 to 2016. The thin lines indicate the trend of each month. The thick horizontal lines represent the mean number of collisions for each month over 13.5 years. The dashed line represents the mean over all observed points. The means for each month displayed a strong seasonal pattern. The TOTAL and PDO collision series reached annual peaks in the months of July and August, possibly due to the high number of tourists in Kelowna, and, in December, possibly due to the heavy snow and high VKT before the holidays. In contrast, a drop in collisions was observed from February to April. The KSI collision series reached an annual minimum in February, rose through September or October, and then began dropping until February. The average number of collisions for each month is presented in Appendix A. 2 A boxplot may be preferable in the case of a large number of observations. 87 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.2 Seasonal subseries plot of monthly collision series in the City of Kelowna (January 2003 to June 2016) 3004005006007008009001000110012000 1 2 3 4 5 6 7 8 9 10 11 12Number of TOTAL collisionsOverall mean20030040050060070080090010000 1 2 3 4 5 6 7 8 9 10 11 12Number of PDO collisionsOverall mean70901101301501701902102300 1 2 3 4 5 6 7 8 9 10 11 12Number of KSI collisionsMonthOverall mean88 Table 6.3 illustrates the Spearman’s rank correlation between the collision series and candidate explanatory variables introduced in Section 6.2. Monthly number of employments (EMP) was the only variable providing high correlation, thus it was selected to feed the multivariate dynamic models. The statistical association between EMP and collisions could be due to the correlation between employment level and economic activity. Generally, economic activity is associated with higher traffic volumes overall. Also, it is obvious that number of collisions can be dominated by traffic volumes in general. The influence of economic conditions on road crashes was studied thoroughly in [96]. Other explanatory variables were taken as the input to the ARIMAX models, but the preliminary results demonstrated that their coefficients were insignificant. In addition, feeding these variables did not improve the performance of the NARX networks. The reason might be that the exposure variables had a continuous growing trend over the study period while the collision series showed a downward trend between 2009 and 2011. As a result, the exposures were unable to characterise that period in the models. Moreover, the sinusoidal behaviour of the weather variables did not represent the seasonal pattern of the collision series because the collision count was high both in summer months and December, which was dissimilar to the weather variables. 89 Table 6.3 Spearman’s rank correlation between collision series data and candidate explanatory variables Collision type Type of dataset Exposure data Weather data Economic conditions AIV* ADL* MINT* MT* MAXT* TP* NDP* EMP* UEMP* TOTAL 0.59 0.46 0.20 -0.11 0.13 0.07 -0.19 0.74 -0.26 PDO 0.58 0.46 0.12 0.15 -0.14 0.09 0.06 0.72 -0.25 KSI 0.46 0.29 0.35 0.10 0.27 0.12 -0.07 0.60 -0.20 *AIV: Actively insured vehicles; ADL: Active driver licenses; MINT: Mean minimum temperature; MT: Mean temperature; MAXT: Mean maximum temperature; TP: Total precipitation; NDP: Number of days with precipitation; EMP: Employment number ( 1,000); UEMP: Unemployment number ( 1,000). 6.3.1 ARIMA models Since the ARIMA approach does not require a validating set, 85% of the collision series (January 2003 to June 2014) were applied for training and 15% (July 2014 to June 2016) for testing the models (Figure 6.3). Before developing the models, the outliers of the training samples were detected using the Hampel filter, presented in Section 5.4.1.2.1, as any sample data point differed from the median more than three standard deviations. The outliers were subsequently replaced by the corresponding median since data points cannot be removed from a time series. The outliers found for each collision type are shown in Table 6.4. 90 Figure 6.3 Collision series division for developing ARIMA approach in Kelowna Table 6.4 Outliers for the monthly collision series in Kelowna (January 2003 to June 2014) Collision type Number of outliers Date TOTAL 1 December 2008 PDO 2 February 2003 December 2008 KSI 2 February 2003 December 2008 The ARIMA models were constructed using the ‘astsa’ package in the RStudio 3.3.3 software program [97]. Figure 6.4 shows the autocorrelation functions (ACFs) of the collision series for the traininig set. Two horizontal dashed lines represent the bounds for the 95% statistical significance. Any lags that lied within the limits are considered unautocorrelated. The ACF of the TOTAL and PDO collision series gradually tapered toward zero over many lags suggesting non-seasonal differencing. The seasonal lags (i.e. k = 12 and 24) decayed slowly in all series indicating a seasonal differencing. 0200400600800100012002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of collisionsTime (year)Total collisions Property damage only collisions Killed or seriously injured collisions85%Training set15%Testing set91 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.4 Autocorrelation functions (ACFs) of monthly collision series for the training set. The horizontal axis represents the size of the lag between observed points of a collision series (in months). The vertical axis represents the autocorrelation value. The horizontal dashed lines represent the bounds for the 95% statistical significance. Any lags that lie within the limits are considered independent. -0.200.20.40.60.81 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.200.20.40.60.81 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.4-0.200.20.40.61 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24AutocorrelationLag (month)92 The plotted ACFs and PACFs of the differenced series (i.e. and12 for the TOTAL and PDO collision series, and12 for the KSI) are shown in Figure 6.5 and Figure 6.6, respectively. There were some statistically significant spikes at early lags of both the ACF and PACF of the TOTAL and PDO collision series followed by the non-significant values for other lags, which suggested the AR, MA, or both terms in the model. Seasonal lags also decayed slowly as the lag increases in both the ACF and PACF, which indicated the SAR, SMA, or both terms in the model. The ACF showed a single spike at lag 12 for the KSI collision series. This was accompanied by a tapering pattern over seasonal lags in the PACF suggesting the SMA(1) term. The tentative models for each collision type are presented in Table 6.5. Since the interpretation of ACF and PACF were not clear for the TOTAL and PDO collisions, there were several possible models. 93 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.5 Autocorrelation functions (ACFs) of differenced collision series. The horizontal axis represents the size of the lag between observed points of a differenced collision series (in months). The vertical axis represents the autocorrelation value. The horizontal dashed lines represent the bounds for the 95% statistical significance. Any lags that lie within the limits are considered independent. -0.6-0.4-0.200.20.41 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.4-0.200.20.41 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.6-0.4-0.200.21 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24AutocorrelationLag (month)94 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.6 Partial autocorrelation functions (PACFs) of differenced collision series. The horizontal axis represents the size of the lag between observed points of a differenced collision series (in months). The vertical axis represents the autocorrelation value. The horizontal dashed lines represent the bounds for the 95% statistical significance. Any lags that lie within the limits are considered independent. -0.6-0.4-0.200.20.41 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Partial autocorrelation-0.6-0.4-0.200.20.41 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Partial autocorrelation-0.6-0.4-0.200.20.41 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Partial autocorrelationLag (month)95 Table 6.5 Tentative ARIMA models for monthly collision series in Kelowna TOTAL collisions PDO collisions KSI collisions ARIMA(2,1,0) × (2,1,0)12* ARIMA(0,1,2) × (2,1,0)12 ARIMA(1,1,1) × (2,1,0)12 ARIMA(2,1,0) × (0,1,2)12 ARIMA(2,1,0) × (1,1,1)12 ARIMA(0,1,2) × (0,1,2)12 ARIMA(0,1,2) × (1,1,1)12 ARIMA(1,1,1) × (0,1,2)12 ARIMA(1,1,1) × (1,1,1)12 ARIMA(2,1,0) × (2,1,0)12 ARIMA(0,1,2) × (2,1,0)12 ARIMA(1,1,1) × (2,1,0)12 ARIMA(2,1,0) × (0,1,2)12 ARIMA(2,1,0) × (1,1,1)12 ARIMA(0,1,2) × (0,1,2)12 ARIMA(0,1,2) × (1,1,1)12 ARIMA(1,1,1) × (0,1,2)12 ARIMA(1,1,1) × (1,1,1)12 ARIMA(0,0,0) × (0,1,1)12 *ARIMA(p,d,q) × (P,D,Q)S where p represents the non-seasonal autoregressive order, d represents the non-seasonal differencing, q represents the non-seasonal moving order, P represents the seasonal autoregressive order, D represents the seasonal differencing, and Q represents the seasonal moving average order. 6.3.1.1 Estimation of ARIMA models The only candidate model with significant coefficients corresponding to the TOTAL and PDO collision series was ARIMA(2,1,0) × (2,1,0)12. It should be noted that the constant in the KSI collision model was insignificant, thus it was excluded from the model. The results of dynamic estimated coefficients of the ARIMA models are presented in Table 6.6. The coefficients (a.k.a. parameters or estimators) providing a p-value < 0.05 were considered significantly different from zero. The standard error measures how precisely the model estimates the unknown coefficient. The t-value is computed by dividing the estimated coefficient by the standard error. The decreasing absolute magnitude from ø1 to ø2 or Ф1 to Ф2 is confirmed by Figure 6.5 because a reduction was observed from the autocorrelation of lag one to lag two. The negative magnitudes of (S)AR coefficients indicate that the previous values of monthly TOTAL and PDO collisions have a negative impact on the predicted values. The positive magnitude of SMA coefficient in the KSI collision model shows the positive effect of previous errors on the predicted values. The 96 coefficients had opposite signs in the literature, e.g. [27]. The reason may lie in the fact that unlike Kelowna, most collision series have shown a downward trend in general over the last two decades. Table 6.6 Dynamic estimated coefficients of ARIMA models for monthly collision series in Kelowna Collision type Model type Coefficient Estimate Standard error T-value P-value TOTAL ARIMA(2,1,0) × (2,1,0)12 ø1* -0.506 0.083 -6.135 0.000 ø2* -0.448 0.081 -5.558 0.000 Ф12* -0.608 0.082 -7.376 0.000 Ф24* -0.495 0.084 -5.928 0.000 PDO ARIMA(2,1,0) × (2,1,0)12 ø1 -0.521 0.083 -6.314 0.000 ø2 -0.433 0.082 -5.260 0.000 Ф12 -0.583 0.081 -7.129 0.000 Ф24 -0.541 0.083 -6.541 0.000 KSI ARIMA(0,0,0) × (0,1,1)12 Θ12* 0.799 0.098 7.944 0.000 *ø1 denotes the coefficient of AR(1); ø2 denotes the coefficient of AR(2); Ф12 denotes the coefficient of SAR(1); Ф24 denotes the coefficient of SAR(2); Θ12 denotes the coefficient of SMA(1). The fitted ARIMA models belonging to the TOTAL, PDO, and KSI collision series are respectively given by Equation 6.1, Equation 6.2, and Equation 6.3. The TOTAL and PDO collisions are linear functions of their own lags while KSI is a linear function of the regression error. The TOTAL and PDO collision models have similar parameters, meaning that they have a roughly parallel trend and dynamic variation. (1 + 0.608B12 + 0.495B24)(1 + 0.506B + 0.448B2) 12 ln(TOTALt) = ωt (6.1) (1 + 0.583B12 + 0. 541B24)(1 + 0.521B + 0.433B2) 12 ln(PDOt) = ωt (6.2)  12 ln(KSIt) = (1 + 0.799B12)ωt (6.3) where: 97  B, B2, B12, and B24 respectively denote the back-shift operators with the orders of 1, 2, 12, and 24 months,  ▽ and 12respectively denote the difference operators of 1 and 12 months,  TOTALt denotes the number of TOTAL collisions at month t,  PDOt denotes the number of PDO collisions at month t,  KSIt denotes the number of KSI collisions at month t, and  ωt is the error term. 6.3.1.2 Diagnostic statistical tests of ARIMA models This section investigates whether the residuals of the ARIMA models follow the assumption that the error terms are the “Gaussian white noise”, which was discussed in Section 3.2. The plotted ACFs of the residuals are statistically insignificant except for a few lags that can be ignored, which represent that the residuals are not autocorrelated (Figure 6.7). The standardized residuals illustrate a homoscedastic variance with a mean of zero (Figure 6.8). The normal probability plots show that the points fall roughly along the straight line, which justifies the normal distribution of the residuals (Figure 6.9). The p-values of the Ljung-Box test do not reject the null hypothesis that the accumulated sample autocorrelations are insignificant (p-value > 0.05), except for the PDO series that is marginally significant (Table 6.7). Consequently, the residuals are well behaved for all ARIMA models and represent Gaussian white noise. 98 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.7 Autocorrelation functions (ACFs) of residuals for selected ARIMA models. The horizontal axis represents the size of the lag between observed points of a collision series (in months). The vertical axis represents the autocorrelation value. The dashed lines represent the bounds for the 95% statistical significance. Any lags that lie within the limits are considered independent. -0.200.21 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.200.20.41 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.4-0.200.21 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24AutocorrelationLag (month)99 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.8 Standardized residuals of selected ARIMA models -3-2-101232003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014-4-3-2-101232003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014-3-2-1012342003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014Time (year)100 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.9 Normal probability plots of selected ARIMA models. The horizontal axis represents the ordered residual values. The vertical axis represents the normal order statistic medians. -3-2-101234-100 -50 0 50 100Normal order static medians-3-2-101234-100 -50 0 50 100 150Normal order static medians-3-2-101234-40 -20 0 20 40 60 80Normal order static mediansOrdered residual values101 Table 6.7 Ljung-Box statistic of selected ARIMA models Collision type P-value TOTAL 0.217 > 0.05 PDO 0.047 < 0.05 KSI 0.280 > 0.05 6.3.1.3 Predicting with ARIMA models In the previous sections, the ARIMA models were developed for each collision series. The model assumptions were also tested and satisfied. The next step is to predict the collision series for the testing set through the fitted ARIMA models. The performance of training and testing sets of ARIMA models is presented in Table 6.8. Commandeur et al. [27] performed a comprehensive study and developed multiple ARIMA models for different regions and collision types. The R2 values fluctuated from 0.38 to 0.95 among the models. All R2 values of the ARIMA models in this research are within this range. MAE represents the magnitude of the overall error and the PDO collision model shows greater a MAE value than the TOTAL collision model for the testing set. Although the average of TOTAL collisions is larger than the average for the PDO collisions, it can be concluded that the ARIMA model provides better predictions for the TOTAL collision series. All ARIMA models provide overestimating prediction on the training set (for negative MPEs) while underestimating for the testing set (for positive MPEs). For the TOTAL and PDO collision models, the performance of the testing set is lower than the training set; however, the KSI collision model performs similarly on both sets. This means that the TOTAL and PDO collision models are overfitted, unlike the KSI collision model. 102 Table 6.8 Performance of ARIMA models for monthly collision series in Kelowna Collision type Training set Testing set R2 MAE* MAPE MPE MSE* R2 MAE MAPE MPE MSE TOTAL 0.64 40.56 5.95% -0.59% 2824 0.54 45.50 5.79% 4.25% 3638 PDO 0.63 35.34 6.62% -0.41% 2041 0.45 49.29 7.87% 6.19% 3889 KSI 0.41 13.58 9.63% 4.26% 313 0.41 16.58 10.69% 9.07% 447 *MAE has the same unit as the data, i.e. number of collisions, and MSE has the squared unit of the data, i.e. (number of collisions)2. The observed and predicted number of collisions over time are illustrated in Figure 6.10. The uncertainty with regard to the 95% confidence interval is constant on the training set of the models. The uncertainty of the testing set increases considerably in the TOTAL and PDO collision models as it goes further in time, even though it is constant for the KSI collision models. One likely reason may be the remarkable growth of the TOTAL and PDO collisions, unlike the KSI collisions, over the testing period.103 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.10 Estimations of ARIMA models for monthly collision series in Kelowna300600900120015002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of TOTAL collisionsObserved Predicted (training set) Predicted (testing set) Lower confidence limits Upper confidence limits200400600800100012002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of PDO collisions50901301702102502003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of KSI collisionsTime (year)104 6.3.2 ARIMAX models The first step in building the ARIMAX models was to identify the time lags of the EMP series that were the best predictors of the collision series. The cross-correlation functions (CCFs) between the collision series and lagged values of the EMP are shown in Figure 6.11. The CCFs provide an unclear relationship between these two series and indicate a need for the “pre-whitening strategy”, which was discussed in Section 5.3, as an aid in interpreting the CCF. The CCF between pre-whitened EMP (i.e. residuals after a non-seasonal and seasonal differencing) and filtered collision series corresponding to the differencing operators are presented in Figure 6.12. As the EMP variable leads the collision series, the negative values are only considered. There are clear spikes at k = {-1, -7, -12} for the TOTAL collision series, k = {-1, -2, -3, -7} for the PDO series, and k = {-1, -3, -4, -6, -8} for the KSI series in Figure 6.12. 105 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.11 Cross-correlation functions (CCFs) of the monthly employment (EMP) series against collision series. The horizontal line represents the size of the lag between observed points of employment and collision series (in months). The vertical line represents the cross-correlation value. The dashed lines represent the bounds for the 95% statistical significance. Any lags that lie within the limits are considered independent. -0.200.20.40.60.8-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11Cross-correlation-0.200.20.40.60.8-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11Cross-correlation-0.200.20.40.6-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11Cross-correlationLag (month)106 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.12 Cross-correlation functions (CCFs) of pre-whitened monthly employment (EMP) series against filtered collision series. The horizontal axis represents the size of the lag between observed points of employment and collision series (in months). The vertical axis represents the cross-correlation value. The dashed lines represent the bounds for the 95% statistical significance. Any lags that lie within the limits are considered independent. -0.4-0.200.20.4-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11Cross-correlation-0.4-0.200.20.4-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11Cross-correlation-0.4-0.200.20.4-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11Cross-correlationLag (month)107 6.3.2.1 Estimation of ARIMAX models The results of explanatory and dynamic estimated coefficients of the ARIMAX models are presented in Table 6.9. EMPt-1 provided a significant coefficient for all models. Therefore, other predictors were removed from the models. Unlike the ARIMA models, the constants of ARIMAX models were significant. The reason is that implementing two differencing operators made the mean of the series zero in the ARIMA models. The positive coefficients of the EMP variable suggest a direct association with the number of collisions. Table 6.9 Explanatory and dynamic estimated coefficients of ARIMAX models for monthly collision series in Kelowna Collision type Coefficient Estimate Standard error T-value P-value TOTAL c* -9.242 2.092 -4.418 0.000 ø1* 0.338 0.082 4.146 0.000 Ф12* 0.479 0.076 6.267 0.000 β1* 2.112 0.281 7.510 0.000 PDO c -7.834 2.226 -3.519 0.001 ø1 0.329 0.083 3.990 0.000 Ф12 0.467 0.078 5.949 0.000 β1 1.890 0.299 6.314 0.000 KSI c -8.360 2.555 -3.272 0.001 ø1 0.212 0.086 2.471 0.015 Ф12 0.348 0.083 4.187 0.000 β1 1.790 0.343 5.212 0.000 *c denotes the constant of the model; ø1 denotes the coefficient of the AR(1); Ф12 denotes the coefficient of the SAR(1); β1 denotes the coefficient of the EMPt-1 variable. The fitted ARIMAX models corresponding to the TOTAL, PDO, and KSI collision series are respectively given by Equation 6.4, Equation 6.5, and Equation 6.6. Interestingly, the ARIMA 108 structure of all ARIMAX models are the same. The association between economic conditions and the number of collisions of the next month is consistent with the literature [96]. (1 − 0.479B12)(1 − 0.338B) ln(TOTALt) = −9.242 + 2.112 ln(EMPt−1) + ωt (6.4) (1 − 0.467B12)(1 − 0.329B) ln(PDOt) = −7.834 + 1.89 ln(EMPt−1) + ωt (6.5) (1 − 0.348B12)(1 − 0.212B) ln(KSIt) = −8.36 + 1.79 ln(EMPt−1) + ωt (6.6) where:  EMPt-1 denotes the number of employment at month t-1 and  ωt is the error term. 6.3.2.2 Diagnostic statistical tests of ARIMAX models The assumptions of the ARIMAX models were tested using a procedure similar to the ARIMA models. The ACFs of the residuals were independent except for some negligible lags (Figure 6.13). The standardized residuals were well scattered around zero (Figure 6.14). The points formed approximately a straight line in the normal probability plots (Figure 6.15). 109 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.13 Autocorrelation functions (ACFs) of residuals for selected ARIMAX models. The horizontal axis represents the size of the lag between observed points of employment (EMP) and collision series (in months). The vertical axis represents the cross-correlation value. The dashed lines represent the bounds for the 95% statistical significance. Any lags that lie within the limits are considered independent. -0.20.00.21 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.20.00.20.41 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Autocorrelation-0.20.00.21 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24AutocorrelationLag (month)110 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.14 Standardized residuals of selected ARIMAX models -3-2-101232003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014-3-2-1012342003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014-3-2-101232003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014Time (year)111 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.15 Normal probability plots of selected ARIMAX models. The horizontal axis represents the ordered observed values. The vertical axis represents the normal order statistic medians. -3-2-101234-150 -100 -50 0 50 100 150Normal order static medians-3-2-101234-150 -100 -50 0 50 100 150Normal order static medians-3-2-101234-60 -40 -20 0 20 40 60Normal order static mediansOrdered residual values112 The p-value of the Ljung-Box test corresponding to the TOTAL and KSI collision models was acceptable, but the p-value of the PDO model rejected the overall randomness of the residuals (Table 6.10). Consequently, the ARIMAX models satisfied all diagnostic tests except the Ljung-Box test of the PDO collision model. Table 6.10 Ljung-Box statistic of selected ARIMAX models Collision type P-value TOTAL 0.246 > 0.05 PDO 0.035 < 0.05 KSI 0.554 > 0.05 6.3.2.3 Predicting with ARIMAX models The ARIMAX models were constructed and the assumptions were checked. In this section, the future values of collision series were predicted using the fitted models. Table 6.11 gives the models’ performance on the training and testing sets by means of different indices. Table 6.11 Performance of ARIMAX models for monthly collision series in Kelowna Collision type Training set Testing set R2 MAE MAPE MPE MSE R2 MAE MAPE MPE MSE TOTAL 0.67 40.15 5.98% 0.67% 2743 0.37 62.38 8.09% 7.26% 6226 PDO 0.63 35.18 6.78% 0.88% 2074 0.33 66.04 10.85% 10.32% 6499 KSI 0.38 13.57 9.27% 1.04% 292 0.34 14.50 8.81% 0.99% 288 The literature revealed dissimilar results when adding explanatory variables to the ARIMA models. For example, Commandeur et al. [27] showed that this action improved the R2 between 1 to 4% except for one model. However, it dramatically enhanced the R2 by 45% in another study 113 [98]. When comparing Table 6.8 and Table 6.11, adding an EMP variable improves the R2 of TOTAL collisions by 4.7%, does not affect the PDO, and worsens the KSI by 7.3% on the training set. The KSI collision model presents roughly similar predictive performance for both the training and testing sets; however, the TOTAL and PDO collision models provide poor performance on the testing set when compared to the training set. This implies that the models are overfitted. Generally, all models yield underestimated predictions due to the positive MPE values. High MSE values of the TOTAL and PDO collision models indicate extreme individual errors. Figure 6.16 shows the observed and predicted monthly number of collisions over time using the ARIMAX models. The uncertainty of all models is constant through time. 114 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.16 Estimations of ARIMAX models for monthly collision series in Kelowna30040050060070080090010002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of TOTAL collisionsObserved Predicted (training set) Predicted (testing set) Lower confidence limits Upper confidence limits2003004005006007008002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of PDO collisions60901201501802102402003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of KSI collisionsTime (year)115 6.3.3 Dynamic artificial neural network models The collision series were divided into three contiguous segments; 70% of the data (January 2003 to June 2012) were applied for training, 15% (July 2013 to June 2014) for validating, and 15% (July 2014 to June 2016) for testing the networks (Figure 6.17). The DANNs were trained based on the methodology outlined in Figure 5.2, which is repeated below for convenience, and used the MATLAB R2016b software program [99]. A part of the code is presented in Appendix B. The pre-training, training, and post-training steps are provided together since in practice they are hardly separable. First, the results of the nonlinear autoregressive (NAR) networks are presented followed up by the nonlinear autoregressive with extremal inputs (NARX) networks. Figure 6.17 Collision series division for developing DANN approach in Kelowna 0200400600800100012002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of collisionsTime (year)Total collisions Property damage only collisions Killed or seriously injured collisions70%Training set15%Testing set15%Validating set116 Figure 5.2 Flowchart of training dynamic artificial neural network (DANN)Pre-training level Training level Post-training level Network performance is satisfactory End Yes Collect data Preprocess data Select network architecture Select optimum lag inputs Select optimum number of hidden layers and neurons Select learning algorithm Initialize network parameters Select stopping criteria Train network De-normalize network output Analyze network Try different network No 117 6.3.3.1 Preprocess data Appendix C presents the performance of the validating sets corresponding to the networks that were trained on the raw collision series in City of Kelowna. As can be seen, the R2 values were deficient. Thus, three data preprocessing techniques were performed on the raw data to prepare for feeding the networks: 1) Outlier detection: The data points referred to as outliers in Table 6.4 were also applied here. 2) Wavelet transform: The optimal wavelet resolution level was selected based on the amount of training data for the given time series using Equation 5.17: L = int[log(138)] = 2 (6.7) For a two-level discrete wavelet transform (DWT), the original time series was decomposed into approximations of A1 and details of D1. A1 was then decomposed into A2 and D2. As depicted in Figure 6.18, the subseries of A2, D1, and D2 were taken as inputs to the networks to predict the future patterns of each subseries. At the end, the predicted values of subseries were summed to form the predictions of the original collisions series. 118 Figure 6.18 Flowchart of dynamic artificial neural network (DANN) coupled with discrete wavelet transform (DWT) The A2, D1, and D2 subseries of training data belonging to each collision type are illustrated in Figure 6.19. The same methodology was applied for the NARX networks; the external input was first decomposed into the A2, D1, and D2 subseries (Figure 6.20) and they were subsequently taken as the input to the corresponding network. As can be seen, the approximations contain the trend structure and low-frequency features of the original series; however, the details represent discontinuities, singularities and high-frequency features. Original time series (training samples) D1 DANN2 DANN1 Predict the original time series by adding up the predicted values of the subseries DANN3 D2 A2 Wavelet decomposition D෡1 D෡2 A෡2 119 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.19 Approximated (A2) and detailed (D1 and D2) subseries obtained by using a two-level discrete wavelet transform (DWT) of monthly collision series in Kelowna 4006008002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014A2-700702003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D1-900902003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D23505006502003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014A2-600602003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D1-600602003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D21001401802003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014A2-300302003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D1-250252003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D2Time (year)120 Figure 6.20 Approximated (A2) and detailed (D1 and D2) subseries obtained by using a two-level discrete wavelet transform (DWT) of monthly employment (EMP) number in the Province of British Columbia 3) Data normalization: The input and target data were scaled between [-1, 1] according to Equation 5.21 because hyperbolic tangent sigmoid was selected as the transfer functions in the hidden layer(s). Also, a linear function was applied in the output layer to generate the predicted values. 6.3.3.2 Choice of network architecture As the target was available for training a dynamic network, series-parallel (a.k.a. open-loop) feedforward network architecture was first constructed using one tapped delay line (TDL) for NAR networks and two for NARX networks (Figure 6.21a). After the network was trained, the feedforward network was replaced by parallel (a.k.a. closed-loop) recurrent network architecture to simulate multi-step ahead prediction (Figure 6.21b). It should be noted that the real future values of the external input (i.e. monthly number of employments) were fed to the NARX networks to predict the testing set since its projection was out of the scope of this research. -250252003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D2-300302003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D1-250252003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014D2Time (year)121 (a) Series-parallel network architecture (b) Parallel network architecture Figure 6.21 Dynamic artificial neural networks (DANNs) As depicted in Figure 6.22, a TDL has p outputs, including the delays (i.e. lags) from 1 to t-p time steps of the input signal. Figure 6.22 Tapped delay line D yt yt-1 D D yt-p Feedforward network TDL TDL xt ?̂?𝑡 yt xt Recurrent network TDL TDL ?̂?𝑡 122 6.3.3.3 Choice of learning algorithm Gradient descent optimization was selected for training the networks since it is fast enough for a network having up to a few hundred parameters. Network weights and biases were updated using batch mode because it is more stable and fast enough to train a network with roughly 100 training samples. 6.3.3.4 Network generalization Global searches and early stopping that were discussed in Section 4.10 were applied to control overfitting and to find the right capacity. Moreover, the performance of each network configuration was averaged over several training runs with different initial conditions. It should be noted that other network generalization methods were not suitable for this research because data samples were limited, and the order of observations cannot change in time series data. 6.3.3.5 Stopping criteria The network training stops when any of the following criteria occurs: 1) zero for the minimum level of error, 2) 1,000 for the maximum number of epochs, 3) infinite for the maximum amount of training time, 4) 10-6 for the minimum gradient of the performance function, 5) zero for the minimum reduction in performance index (i.e. MSE), and 6) six times for the consecutive number of iterations that the performance index is failed to improve on the validating set. 123 For most networks in this research, the early stop (number 6) occurred to stop the training process, then the maximum number of epochs (number 5). 6.3.3.6 Choice of optimum input lags and network specifics To find the optimal number of input lags, the model-based exhaustive search technique, which was discussed in Section 5.4.1.5, was selected because this method can be employed for choosing optimal lag order and network specifics simultaneously. Multiple networks were trained with different numbers of hidden layers, ranging from one to four and neurons ranging from one to ten through one to twelve lag orders, for a total of 4  10  12 = 480 network configurations. An ANN model might generate dissimilar results each time the network is trained with a specific configuration with different initial conditions. To provide better generalization results, the literature recommended five to ten times training runs using different initial conditions [100]. Each potential network configuration was therefore retrained ten times with a separate set of initial weights and biases for a total of 480  10 = 4800 number of networks. Then, the network performance was averaged over ten runs. The following presents the performance of the trained NAR and NARX networks in Sections 6.3.3.7 and 6.3.3.8, respectively. 6.3.3.7 NAR networks Table D.1 to Table D.9 in Appendix D present the performance (in terms of MSE of validating set) of 480 network configurations corresponding to three subseries for each collision type that were derived from the DWT. Networks with three or four hidden layers are neglected because they yielded higher MSE values when compared to one or two hidden layers for all subseries, which represented the networks start to be overfitted. The best performance among the networks with an 124 equal number of hidden layers (one or two) and neurons but a different number of input lags are highlighted in these tables. As shown in Table D.1, for the approximated subseries (A2) of TOTAL collision series, networks with the lag order of eleven and twelve provided the lowest MSE in most configurations with one or two hidden layers, except for the network with two hidden layers and seven neurons. Accordingly, the optimum number of input lags was eleven or twelve. The same procedure was applied for the rest of subseries. To evaluate the precision of the network performance, standard deviation (SD) of every configuration among ten runs was measured for the networks with the chosen lag orders (e.g. eleven and twelve for the A2 subseries of TOTAL collision series). The MSE of the chosen lag orders and their corresponding SD (shown in italic) present together in Table 6.12 to Table 6.14 for a better comparison. The configuration providing both the lowest MSE and SD was selected as the optimal network (i.e. number of input lags, hidden layers, and hidden neurons). For example, for the A2 subseries of TOTAL collision series, networks having one hidden layer and two neurons with eleven or twelve lags demonstrated the lowest MSE. However, the SD of the network with two lags was significantly lower. Consequently, the network with one hidden layer and two neurons and the input vector of {Tt-1, Tt-2, …, Tt-12} was selected. It should be noted that the MSEs and SDs of the networks corresponding to the KSI collision series were very similar (Table 6.14). Hence, the simplest network was selected. The selected network configurations were highlighted in the tables. 125 Table 6.12 Mean squared error and standard deviation of NAR networks corresponding to monthly number of TOTAL collisions in Kelowna. Standard deviation is in italics under the mean squared error. Subseries Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 1 11 55 53 56 56 62 54 64 58 67 89 2 6 9 12 11 16 18 14 18 60 12 57 54 59 54 55 56 64 68 84 91 1 3 2 11 11 6 21 15 28 33 2 11 156 426 54 70 98 178 147 153 112 217 1 166 8 18 67 132 86 109 42 123 12 157 55 60 431 103 118 150 315 171 159 3 7 9 145 79 54 72 388 227 54 D1 1 3 215 131 122 116 110 86 111 90 103 91 194 19 30 57 53 18 36 17 26 46 7 75 89 91 106 164 112 145 139 138 154 3 19 19 23 90 21 62 44 41 55 2 3 216 139 168 123 93 104 126 97 127 147 194 24 116 66 33 26 38 28 45 74 7 264 227 175 107 122 278 180 262 181 191 271 269 216 30 38 237 63 230 84 41 D2 1 9 42 42 36 44 42 46 48 39 42 43 3 5 11 7 13 14 9 14 10 18 10 79 32 35 32 34 38 41 44 40 60 144 3 5 7 8 8 9 6 9 49 12 31 34 31 36 32 41 47 42 49 47 1 4 11 8 7 10 16 13 16 17 2 9 408 113 155 168 204 48 57 62 63 56 391 193 227 273 263 15 10 27 38 22 10 408 113 155 168 204 48 57 62 63 56 377 314 42 122 283 109 176 14 83 92 12 451 333 69 65 80 60 105 186 163 142 421 392 90 51 113 23 84 223 221 117 126 Table 6.13 Mean squared error and standard deviation of NAR networks corresponding to monthly number of PDO collisions in Kelowna. Standard deviation is in italics under the mean squared error. Subseries Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 1 10 65 61 59 59 67 64 78 78 84 76 3 11 10 18 10 19 24 17 21 28 11 52 50 53 52 52 56 57 66 71 82 4 4 7 10 7 7 12 11 15 19 12 54 52 55 48 54 62 74 63 80 79 3 5 10 13 10 10 20 13 28 18 2 10 61 62 77 74 68 93 100 180 134 144 7 3 19 32 20 28 25 129 70 54 11 56 58 58 58 74 160 90 103 129 165 1 7 8 14 17 199 25 19 67 56 12 54 54 67 77 76 145 114 202 155 186 3 5 34 28 13 126 38 178 58 71 D1 1 4 117 106 102 99 93 90 84 90 84 88 17 15 12 18 16 22 16 10 22 42 5 94 89 90 92 89 92 88 95 92 89 38 37 25 37 31 49 33 35 35 41 2 4 255 146 192 120 98 89 99 111 96 130 243 138 181 36 15 23 32 76 35 115 5 283 168 91 120 145 99 122 106 150 105 267 131 176 35 27 27 35 31 40 28 D2 1 10 74 28 64 21 20 21 20 57 16 25 110 3 145 5 6 6 9 132 5 8 12 31 28 23 48 20 33 70 22 34 47 1 4 6 103 6 48 132 7 41 76 2 10 356 107 57 40 36 107 209 63 107 72 290 211 110 38 32 183 250 81 144 74 12 288 147 146 96 45 55 66 108 80 74 353 215 258 170 52 79 93 143 67 77 127 Table 6.14 Mean squared error and standard deviation of NAR networks corresponding to monthly number of KSI collisions in Kelowna. Standard deviation is in italics under the mean squared error. Subseries Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 1 10 2 2 1 1 2 1 1 1 2 2 0 1 0 1 1 1 0 0 1 1 11 1 1 1 1 1 1 1 2 1 2 0 0 1 0 1 1 1 1 1 1 12 1 1 1 1 1 1 1 2 2 2 0 0 0 0 1 1 0 1 1 0 2 10 2 2 1 2 5 3 7 3 2 3 2 0 0 0 9 4 8 1 1 1 11 1 1 1 3 2 3 3 4 3 11 0 1 1 4 1 2 3 5 1 11 12 1 3 2 2 3 3 2 5 15 4 0 3 0 1 6 2 1 6 20 2 D1 1 9 9 6 6 6 6 8 7 8 7 9 8 1 1 1 2 3 2 3 1 3 10 7 10 6 8 8 7 8 7 9 9 1 1 1 1 2 2 2 2 1 2 2 9 32 14 9 9 11 8 9 14 13 9 16 12 5 6 9 3 2 9 9 4 10 13 13 7 11 12 12 15 12 12 13 13 10 2 10 8 6 8 3 10 4 D2 1 10 6 5 5 5 5 5 5 6 6 7 1 1 2 2 1 1 2 2 3 0 11 9 5 5 5 5 6 5 5 12 7 1 1 2 3 2 1 3 15 2 9 12 5 5 4 4 5 5 8 9 11 7 1 1 1 1 1 1 9 12 19 3 2 10 24 15 10 23 18 8 14 12 13 11 23 8 23 25 4 8 4 14 4 28 11 13 8 6 8 14 10 14 9 17 10 23 7 1 3 13 8 20 2 20 3 12 12 12 20 11 15 10 8 21 15 13 23 20 26 14 16 4 3 21 16 4 128 The optimal network configuration for each subseries corresponding to each collision type is presented in Table 6.15. Since each subseries represents distinctive features, they have different configurations. Table 6.15 Optimal NAR network configurations of approximated (A2) and detailed (D1 and D2) subseries corresponding to each collsion type in Kelowna Collision type Subseries Optimal network configuration TOTAL A2 12-2-1 D1 7-1-1 D2 12-1-1 PDO A2 11-2-1 D1 4-7-1 D2 10-9-1 KSI A2 10-3-1 D1 9-2-1 D2 12-3-1 Table 6.16 shows the models’ performance by means of different indices. Compared to the validating set, the error of the testing sets is higher, especially for the KSI collision network. Despite overestimated predictions for the validating set (negative MPE value), the TOTAL collision network underestimates the testing set (positive MPE value). However, the other two NAR networks underestimate both the training and validating sets. The MSE values of the training and testing sets are considerably larger than validating set. This suggests that there are some extreme induvial prediction errors in the training and testing sets. 129 Table 6.16 Performance of NAR networks corresponding to each collsion type in Kelowna Collision type Data set Performance measure R2 MAE MAPE MPE MSE TOTAL Training 0.87 25.18 3.74% -0.06% 949 Validating 0.96 10.14 1.45% -0.12% 135 Testing 0.88 19.21 2.44% 0.75% 708 PDO Training 0.92 15.81 3.00% -0.14% 400 Validating 0.95 8.75 1.63% 0.15% 111 Testing 0.92 14.30 2.30% 1.40% 388 KSI Training 0.97 3.17 2.20% 0.17% 16 Validating 0.97 2.46 1.63% 0.22% 11 Testing 0.84 7.16 4.44% 1.13% 72 Figure 6.23 shows the observed versus the predicted monthly number of collisions over time using the NAR networks. The NAR networks performed well because the predicted values are roughly the same as observations except for some seasonal peaks. 130 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.23 Estimations of NAR networks corresponding to monthly collision series in Kelowna 30040050060070080090010002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of TOTAL collisionsObserved Predicted (training set) Predicted (validating set) Predicted (testing set)3004005006007008002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of PDO collisions60901201501802102402003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of KSI collisionsTime (year)131 6.3.3.8 NARX networks As seen in Table D.1 to Table D.9 in Appendix D, the NAR networks with three or four hidden layers perform poorly when compared with networks with one or two hidden layers (i.e. higher MSE values). In addition, networks with two hidden layers yield unstable performance when compared with networks with one hidden layer (i.e. higher SD values). The NARX networks were thus constructed with only one hidden layer. The first lag of EMP series was only fed to the networks because they reveal satisfactory performance in validating set. Table D.10 to Table D.18 in Appendix D provide the MSE of various network configurations on validating set corresponding to each collision type and DWT subseries. The best performance among the networks with an equal number of hidden neurons but a different number of input lags are highlighted. Table 6.17 to Table 6.19 depict the MSEs of the chosen lags and their corresponding standard deviations (SDs, shown in italic). The configuration showing both the lowest MSE and SD is selected as the optimal network. For example, for the A2 subseries of TOTAL collision series, networks having three hidden neurons with eleven dependent lags demonstrate the lowest MSE and SD (Table 6.17). It should be noted that the SD of the networks with two hidden neurons and eleven or twelve lags is slightly lower, but their MSE values are considerably higher. Consequently, the network with three hidden neurons and the input vector of {EMPt-1, Tt-1, Tt-2, …, Tt-12} was selected. The same procedure was employed for the rest of subseries. 132 Table 6.17 Mean squared error and standard deviation of NARX networks corresponding to monthly number of TOTAL collisions in Kelowna. Standard deviation is in italics under the mean squared error. Subseries Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 11 41 46 41 52 48 73 60 138 88 64 10 5 6 19 15 34 14 210 53 20 12 177 45 50 48 47 47 47 73 189 81 425 5 14 12 18 10 17 28 345 34 D1 3 152 160 135 134 96 119 129 112 105 139 4 23 25 48 18 37 83 33 28 51 6 213 192 194 186 181 174 196 698 162 206 15 23 30 22 59 62 52 1659 73 53 11 149 94 103 126 126 158 180 183 253 248 174 7 16 48 34 53 66 68 146 103 D2 10 35 37 36 39 43 44 46 51 47 58 1 3 5 7 8 14 5 14 12 23 11 87 38 40 46 38 48 45 76 51 68 162 3 9 9 9 10 11 34 20 32 12 30 38 38 40 40 41 52 55 44 57 1 9 7 9 17 5 20 29 14 24 133 Table 6.18 Mean squared error and standard deviation of NARX networks corresponding to monthly number of PDO collisions in Kelowna. Standard deviation is in italics under the mean squared error. Subseries Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 11 77 44 44 51 51 48 65 46 57 60 0.1 0.2 0.3 0.4 0.3 0.3 0.6 0.5 0.5 0.7 12 173 46 43 47 57 51 74 61 77 147 44.4 0.3 0.3 0.2 0.4 0.4 0.5 3.3 1.2 0.6 D1 5 125 83 84 85 100 106 106 114 108 133 113 8 11 24 16 19 19 31 32 33 12 116 76 83 84 155 92 110 131 173 174 124 15 13 23 192 21 31 35 120 157 D2 10 44 41 41 40 35 30 81 48 70 43 0 6 11 15 7 5 115 15 103 19 11 44 39 42 74 44 41 42 41 53 54 0 8 12 113 17 15 19 17 13 20 12 70 34 34 42 43 43 74 53 52 53 96 5 10 11 8 18 95 21 24 18 134 Table 6.19 Mean squared error and standard deviation of NARX networks corresponding to monthly number of KSI collisions in Kelowna. Standard deviation is in italics under the mean squared error. Subseries Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 10 1 1 1 1 1 1 2 1 1 3 0.3 0.2 0.4 0.5 0.3 0.2 0.7 0.3 0.4 5.3 11 1 1 1 1 1 1 1 1 1 1 0.3 0.2 0.3 0.4 0.3 0.3 0.6 0.5 0.5 0.7 12 20 1 1 1 1 1 2 2 2 1 44.4 0.3 0.3 0.2 0.4 0.4 0.5 3.3 1.2 0.6 D1 7 10 10 10 10 10 12 12 11 15 17 0.4 0.7 1.3 2.1 2.4 2.4 2.4 4.6 4.4 5.7 8 10 10 9 9 12 10 12 14 16 16 0.2 2.3 1.6 2.1 2.4 2.8 3.8 4.5 6.0 3.5 D2 11 6 6 5 5 5 7 6 7 11 8 0.2 0.6 0.9 1.4 1.3 2.4 2.6 3.2 8.2 3.2 12 4 4 4 5 6 5 6 6 10 8 0.4 0.6 0.6 1.1 1.5 0.8 2.4 1.9 11.0 2.3 Table 6.20 presents the optimal network configurations. Table 6.21 displays the performance of summing the predicted values of three subseries for each collision type on training, validating, and testing set. The performance of the validating and testing sets are similar for the TOTAL and PDO collision networks while the KSI collision network is rather overfitted (i.e. showing higher error on the testing set). 135 Table 6.20 Optimal NARX network configurations of approximated (A2) and detailed (D1 and D2) subseries corresponding to each collsion type in Kelowna Collision type Subseries Optimal lag order of dependent data Optimal network configuration TOTAL A2 11 12-3-1 D1 11 12-2-1 D2 12 13-1-1 PDO A2 12 13-3-1 D1 12 13-2-1 D2 12 13-2-1 KSI A2 10 11-2-1 D1 8 9-1-1 D2 12 13-1-1 Table 6.21 Performance of NARX networks corresponding to each collsion type in Kelowna Collision type Data set Performance measure R2 MAE MAPE MPE MSE TOTAL Training 0.99 8.34 1.23% -0.10% 109.72 Validating 0.96 9.87 1.45% -0.02% 156.49 Testing 0.98 9.52 1.18% -0.18% 140.22 PDO Training 0.99 6.98 1.32% 0.01% 79.15 Validating 0.95 10.68 2.00% 0.41% 157.05 Testing 0.98 8.90 1.38% 0.15% 109.32 KSI Training 0.96 3.69 2.55% 0.17% 23.01 Validating 0.96 3.00 2.03% -0.20% 16.43 Testing 0.91 5.34 3.26% 0.27% 42.10 Figure 6.24 compares the observed versus the NARX network simulated data for every collision series. Unlike the NAR networks, the NARX networks accurately predict the seasonal peaks. Intuitively, this suggests that similar seasonal variation of the external input series with the collision series may assist the networks in training. 136 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.24 Estimations of NARX networks corresponding to collision series in Kelowna 30040050060070080090010002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of TOTAL collisionsObserved Predicted (training set) Predicted (validating set) Predicted (testing set)3004005006007008002003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of PDO collisions60901201501802102402003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016Number of KSI collisionsTime (year)137 6.3.4 Comparison of ARIMA and DANN models This section compares the performance of the ARIMA and DANN approaches in predicting collision series. Figure 6.1 illustrates the observed versus the predicted values of the models for the 24-month testing period. In general, the predicted values of the NAR(X) networks are much closer to the observations compared with ARIMA(X) models. The ARIMA model underestimates all collision series and the ARIMAX model behaves similarly for the TOTAL and PDO collision series. 138 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.25 A comparison of the observed monthly number of collisions versus predicted ones over the testing set (July 2014 to June 2016) 60070080090010001 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Number of TOTAL collisionsObserved ARIMA model ARIMAX model NAR net. NARX net.4705406106807501 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Number of PDO collisions901201501802101 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24Number of KSI collisionsTime (month)139 Table 6.22 presents the models’ performance on the testing set in terms of different indices. The NAR(X) networks provide lower errors than the ARIMA(X) models in all collision series. It is interesting that adding the EMP as an explanatory variable to the ARIMA models of the TOTAL and PDO collision series reduced the quality of the model and makes the KSI slightly better. However, adding EMP variable improved the predictability of the NAR networks. The R2 values of the KSI collision models are generally lower than the other two collision types, which suggest that unpredictable, unique conditions may also affect this type. The negative value of MPE corresponding to the NARX network of the TOTAL collisions means that the predicted values are higher than the targets in TOTAL, while other models underestimate the collision series due to positive MPEs. The substantial difference between MAPE and MPE (around 8%) of the ARIMAX model belonging to the KSI collision series explains larger underestimation. The MSE values of the ARIMA(X) models are drastically greater than neural networks indicating large individual errors in the predictions. 140 Table 6.22 A comparison of models’ performances on testing set (July 2014 to June 2016) Collision type Model Performance measure R2 MAE MAPE MPE MSE TOTAL ARIMA 0.54 45.50 5.79% 4.25% 3638 ARIMAX 0.37 62.38 8.09% 7.26% 6226 NAR network 0.88 19.21 2.44% 0.75% 708 NARX network 0.98 9.52 1.18% -0.18% 140 PDO ARIMA 0.45 49.29 7.87% 6.19% 3889 ARIMAX 0.33 66.04 10.85% 10.32% 6499 NAR network 0.92 14.30 2.30% 1.40% 388 NARX network 0.98 8.90 1.38% 0.15% 109 KSI ARIMA 0.41 16.58 10.69% 9.07% 447 ARIMAX 0.34 14.50 8.81% 0.99% 288 NAR network 0.84 7.16 4.44% 1.13% 72 NARX network 0.91 5.34 3.26% 0.27% 42 The MAPE provides a relative magnitude to the targets as a percentage (Equation 4.34), thus it is used as an indicator to measure models’ performance over time. Figure 6.26 depicts the predictive accuracy of the models in terms of MAPE versus eight quarterly horizons across two years. The prediction of the ARIMA(X) models are fairly accurate in the short-term horizon (i.e. first quarter), but the accuracy broadly degrades as the predicting horizon increases. The performance of NAR(X) networks are roughly constant except for the first four quarters of the KSI collision series. 141 (a) TOTAL collisions (b) PDO collisions (c) KSI collisions Figure 6.26 Mean absolute percentage error of the observed values of each collision series versus predicted ones for different quarterly predicting horizons (July 2014 to June 2016) 024681 2 3 4 5 6 7 8MAPE (%)ARIMA model ARIMAX model NAR net. NARX net.02468101 2 3 4 5 6 7 8MAPE (%)0246810121 2 3 4 5 6 7 8MAPE (%)Predicting horizon (quarter)142 6.4 Intervention analysis, Case study two: Springfield Road Corridor The second case study of this research was the Springfield Road Corridor starting from Ziprick Road to Rutland Road with a length of 2.67 km. This road has been classified as an east-west arterial by the City of Kelowna. The road corridor has two lanes and a bike lane in each direction. The posted speed limit is 50 km/hr along the corridor and the traffic volume is low. For example, the average annual daily traffic (AADT) was approximately 17,000 vehicles per day at the east end in 2015. When the corridor was expanded from two to four lanes in 2003, it became a safety hot spot as demonstrated by the ICBC collision reports and resident safety-issues complaints. The corridor represents a mixture of traffic functions. Not only does it connect the Rutland community to the CBD as an arterial, but also provides access to a large number of driveways in the neighborhood as a local road. The low traffic volume and arterial layout may lead drivers to drive faster, on average, than the posted speed limit. In order to address the safety and speeding issue by modifying driving behaviour, the City of Kelowna decided to use a holistic multidisciplinary road safety management system called the ‘3E’ (Engineering, Enforcement, and Education) approach, which was implemented in three phases. The activation date and type of safety measures for each phase are summarized in Table 6.23. 143 Table 6.23 Summary of safety measures implemented in each phase Phase number Implementation date Safety measure One  February 2010 - Present  Variable message signs Two  December 2012 - Present  New variable message signs  Red 50 pavement markings  Traffic signal optimization  City web page  Media campaign  Police speed enforcement Three  April 21 to May 31st of 2015  February 1 to March 24th of 2016  Speed sensitive traffic signals Phase one was installed in February 2010 including two variable message signs (VMSs, a.k.a. speed reader boards or driver feedback signs), which were located west of Gerstmar Road for eastbound traffic and east of Gerstmar Road for westbound traffic. The signs display the speed of the passing vehicle and flashes “SLOW DOWN” when the vehicle’s speed is 5 km/hr over the speed limit. In December 2012, phase two was implemented including: four new VMSs, four red 50 pavement markings, and traffic signal optimization (green wave). In addition, a city web page and media campaign were launched, coupled with an increase in police speed enforcement. The new VMSs are located at west of Ziprick Road and Tamarack Road for eastbound traffic and west of Hollywood Road and Tamarack Road for westbound traffic. The locations of four red 50 pavement markings are east of Ziprick Road and Quigley Road for eastbound traffic and west of Hollywood Road and Graham Road for westbound traffic. Speed sensitive traffic signals (SSTSs) were applied in phase three to stop the vehicles travelling over 70 km/hr by the red signal at the next intersection. SSTSs were launched in ghost mode from the 10th to the 20th of April 2015 and an all-red feature was activated for a brief time from April 21 to May 31st of 2015 and again from 144 February 1 to March 24th of 2016. The time of activation was from 10:00pm to 6:00am to prevent traffic interruption during the day. Since phase three was implemented for such a short time, it was premature to confirm the effectiveness. Phase three was thus removed from this research and should be part of a future research follow-up. The study design that is used to evaluate the effect of any changes in road safety is referred to as a before-after observational safety study [7]. In this regard, only two general scientific approaches have been found in the literature: Empirical Bayes (EB) and intervention analysis. In the EB approach, a safety performance function (SPF) is first developed for similar, but untreated, road sites (known as comparison group) to consider the effect of confounding factors in the treated road site(s). Then the expected number of collisions for the treated road site(s) is estimated using the calibrated SPF. Several types of cross-sectional data can be employed when using SPFs such as traffic, geometric, and weather variables. Traffic volume has always been a pivotal factor in SPF; however, traffic data for this road has been collected in an ad hoc informal manner, relying on the availability of the City staff and budget. For example, data is not available for 2007 and 2009. Continuous traffic volume has been collected via the VMSs since the installation of phase one (February 2010). However, data is missing for various periods because it was corrupt, or the devices were broken. In addition, the new VMSs have been collecting traffic data only for one lane in one direction. Therefore, the accurate traffic volume cannot be attained for Springfield Road, or for other similar arterials in Kelowna. On the other hand, the safety phases were implemented at separate times. Consequently, the EB approach is inapplicable for this case. 145 An alternative approach is intervention analysis or event study, which is referred to as the evaluation of the impact of an external event or special circumstance on the behaviour of a time series when the timing of the event is known [26]. In order to complete intervention analysis and evaluate the impact of safety improvement phases one and two, two separate DANN models were developed using distinctive training sets using the following steps: 1) First DANN model: Develop the first DANN model for the time series by means of the data before the first intervention point (January 2003 to January 2010), 2) Use the trained network from step one to predict the values of what number of collisions should have been had if no safety improvements taken place. These predictions will then be compared with what was observed after the implementation of first intervention, 3) Measure the difference between the observed and predicted values after the first intervention point to estimate the impact of the first safety management phase (December 2010 to November 2012). The differences are also measured after the second intervention point to determine the effectiveness of phase one and two together (December 2012 to June 2016), 4) Second DANN model: Apply step one using the time series data from before the second intervention point (January 2003 to June 2016) in order to predict the values for after the second safety management phase was installed, and 5) Measure the difference between the observed and predicted values obtained from step four to evaluate the impact of second safety management phase (December 2012 to June 2016). Since the average collision count on this road was low, the TOTAL collision series was only considered for the intervention analysis. The same methodology presented in Section 6.3.3 was used to train the NAR networks corresponding to the time series data before installing phase one, 146 including 85 points between January 2003 and January 2010, and also phase two, including 119 points between January 2003 and November 2012. The outliers found for the observed points corresponding to the training data of two NAR networks were the same (Table 6.24). Table 6.24 Outliers of TOTAL collision series for the Springfield Road Corridor corresponding to the training data (i.e. January 2003 to December 2008 for before installing phase one and January 2003 to May 2011 for before installing phase two) Number of outliers Date 2 January 2005 February 2006 The optimal DWT resolution level was two for both time series, as shown in Equation 6.8 and Equation 6.9. The A2, D1, and D2 subseries for the training data are illustrated in Figure 6.27. L = int[log(85)] = 2 (6.8) L = int[log(119)] = 2 (6.9) 147 (a) Pre-phase one (b) Pre-phase two Figure 6.27 Approximated (A2) and detailed (D1 and D2) subseries obtained by using a two-level discrete wavelet transform (DWT) of TOTAL collision series in the Springfield Road Corridor. Figure (a) represents the training data before installing phase one (January 2003 to December 2008). Figure (b) represents the training data before installing phase two (January 2003 to May 2011). The MSE values corresponding to the validating sets of each subseries are presented in Table D.19 to Table D.21 for the NAR networks trained on the data points before installing phase one, and Table D.22 to Table D.24 for the NAR networks trained on the data points before installing phase 0362003 2004 2005 2006 2007 2008A2-5052003 2004 2005 2006 2007 2008D1-2.502.52003 2004 2005 2006 2007 2008D21472003 2004 2005 2006 2007 2008 2009 2010 2011A2-5052003 2004 2005 2006 2007 2008 2009 2010 2011D1-2.502.52003 2004 2005 2006 2007 2008 2009 2010 2011D2Time (year)148 two. The MSE and SD values of the lag orders representing the best performances are shown in Table 6.25 and Table 6.26 for before phase one and before phase two, respectively. Standard deviations are shown in italics under the mean squared error. Table 6.25 Mean squared error and standard deviation of NAR networks corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase one (January 2003 to January 2010). Standard deviation is in italics under the mean squared error. Subseries Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 6 0.029 0.057 0.105 0.066 0.359 0.185 0.639 0.225 0.260 0.475 0.00 0.06 0.17 0.04 0.34 0.18 1.28 0.20 0.15 0.35 12 0.077 0.025 0.072 0.275 0.264 0.267 0.210 0.647 0.479 0.739 0.19 0.01 0.07 0.22 0.24 0.29 0.16 0.73 0.64 0.41 D1 4 1.37 1.89 1.91 2.63 2.13 3.02 2.99 6.82 3.41 3.79 0.01 1.13 0.44 0.75 0.58 1.33 0.86 7.18 1.26 1.94 5 1.11 1.67 2.17 2.17 4.61 3.12 3.62 3.88 5.39 6.15 0.05 0.92 0.88 0.45 3.72 1.46 1.08 1.89 3.77 4.88 D2 4 0.17 0.15 0.14 0.15 0.17 0.14 0.16 0.16 0.24 0.23 0.00 0.02 0.02 0.03 0.06 0.03 0.05 0.04 0.08 0.11 6 0.11 0.11 0.11 0.12 0.14 0.15 0.17 0.17 0.21 0.17 0.01 0.02 0.02 0.03 0.07 0.05 0.08 0.08 0.07 0.07 149 Table 6.26 Mean squared error and standard deviation of NAR networks corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase two (January 2003 to November 2012). Standard deviation is in italics under the mean squared error. Subseries Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 A2 2 0.024 0.019 0.020 0.020 0.018 0.020 0.020 0.020 0.021 0.022 .0015 .0026 .0018 .0034 .0019 .0026 .0040 .0024 .0043 .0038 D1 5 0.19 0.20 0.21 0.20 0.26 0.24 0.26 0.34 0.35 0.33 0.02 0.04 0.05 0.04 0.07 0.05 0.08 0.06 0.10 0.07 D2 8 0.154 0.154 0.132 0.152 0.159 0.163 0.227 0.185 0.235 0.221 0.011 0.026 0.027 0.055 0.118 0.045 0.041 0.074 0.100 0.103 The optimal network configuration for each subseries providing the lowest MSE and SD is depicted in Table 6.27 and Table 6.28 for before phase one and before phase two, respectively. Table 6.27 Optimal NAR network configurations of the approximated and detailed subseries corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase one (January 2003 to January 2010) Subseries Optimal network configuration A2 12-2-1 D1 5-1-1 D2 6-1-1 Table 6.28 Optimal NAR network configurations of approximated and detailed subseries corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase two (January 2003 to November 2012) Subseries Optimal network configuration A2 2-5-1 D1 5-1-1 D2 8-3-1 150 The performance of the sum of the predicted values of three NAR networks (i.e. A2, D1, and D2) on the training, validating, and testing set are presented in Table 6.29 and Table 6.30 for before phase one and before phase two, respectively. The errors of the testing set are higher than validating set of the NAR network corresponding to before installing phase one (Table 6.29). This indicates that the networks are slightly overfitted; however, the errors of validating and testing sets belonging to before installation of phase two are similar (Table 6.30). Table 6.29 Performance of NAR network corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase one (January 2003 to January 2010) Data set Performance measure R2 MAE MAPE MPE MSE Training 0.88 0.26 7.04% -0.18% 0.26 Validating 0.92 0.69 16.19% 12.86% 1.15 Testing 0.83 0.54 6.96% -3.30% 0.85 Table 6.30 Performance of NAR network corresponding to TOTAL collision series in the Springfield Road Corridor before installing phase two (January 2003 to November 2012) Data set Performance measure R2 MAE MAPE MPE MSE Training 0.87 0.44 11.33% 3.68% 0.52 Validating 0.91 0.39 7.47% 2.67% 0.39 Testing 0.90 0.28 6.21% -1.69% 0.28 Figure 6.28 compares the observed versus predicted values of the NAR networks. The predicted values virtually fall on the observations, which suggest that the DANN models performed well in general. However, the networks illustrate poor performance for the seasonal peaks. 151 (a) Pre-phase one (b) Pre-phase two Figure 6.28 Estimations of two NAR networks corresponding to TOTAL collision series in the Springfield Road Corridor. Figure (a) represents the estimations of pre-phase one (January 2003 to January 2010). Figure (b) represents the estimations of pre-phase two (January 2003 to November 2012). The observed versus predicted values of two NAR networks after the installation of phase one and phase two are given in Figure 6.29. The differences between the predicted and observed series (i.e. predicted minus observed values) are shown in Figure 6.30. 0246810122003 2004 2005 2006 2007 2008 2009 2010Number of TOTAL collisionsObserved Predicted (training set) Predicted (validating set) Predicted (testing set)0246810122003 2004 2005 2006 2007 2008 2009 2010 2011 2012Number of TOTAL collisionsTime (year)152 (a) Post-phase one (b) Post-phase two Figure 6.29 Intervention plots for installing safety management phases in the Springfield Road Corridor. The dashed line represents the time of phase two implementation. Figure (a) represents the observed versus predicted values by NAR network trained on data points before phase one. Figure (b) represents the observed versus predicted values by NAR network trained on data points before phase two. 024681012Feb/2010 Feb/2011 Feb/2012 Feb/2013 Feb/2014 Feb/2015 Feb/2016Number of TOTAL collisionsObserved PredictedPost-phase twoPost-phase one024681012Dec/2012 Dec/2013 Dec/2014 Dec/2015Number of TOTAL collisionsTime (month)153 (a) Post-phase one (b) Post-phase two Figure 6.30 Time series plots of the difference between predicted and observed values by two NAR networks. The dashed line represents the time of phase two implementation. Figure (a) represents the difference between predicted and observed values by the NAR network trained on data points before installing phase one. Figure (b) represents the difference between predicted and observed values by the NAR network trained on data points before installing phase two. While it is not immediately obvious from the above figures whether any benefit to the interventions has occurred, Table 6.31 outlines the descriptive statistics of the intervention analyses. The total and mean values represent the accumulated and average observed or predicted monthly number of -6-4-202468Feb/2010 Feb/2011 Feb/2012 Feb/2013 Feb/2014 Feb/2015 Feb/2016Difference in number of TOTAL collisionsPost-phase twoPost-phase one-10-8-6-4-20246Dec/2012 Dec/2013 Dec/2014 Dec/2015Difference in number of TOTAL collisionsTime (year)154 TOTAL collisions during the corresponding data points. The effectiveness shows the percentage change between the total number of observed and predicted TOTAL collisions (observed−predictedobserved × 100). Negative percentage means that the collisions were reduced after the intervention(s), and vise versa. The ranks column represents the number of months where the difference between the predicted and observed TOTAL collision counts were positive, negative, or tie. Table 6.31 Descriptive statistics of intervention analyses. The total and mean columns represent the accumulated and average observed or predicted monthly number of TOTAL collisions during the corresponding data points. The effectiveness column represents the percentage change in the number of collisions due to the intervention(s). The ranks column represents the number of months where the difference between the predicted and observed values are positive, negative, or tie. Intervention Number of data points Data type Total Mean Effectiveness Ranks Negative Positive Tie Phase one 34 observed 118 3.47 -16% 10 16 8 predicted 140 4.12 Phases one and two 43 observed 154 3.65 -12% 16 20 7 predicted 175 4.12 Phase two 43 observed 154 3.65 +10% 20 18 5 predicted 140 3.26 The impact of the safety phases are as follows:  Phase one: The implementation of the first safety management phase saved 16% of collisions on average during 34 months before the installation of the second phase (i.e. February 2010 to November 2012). The mean of monthly number of collisions reduced 155 from 4.12 to 3.47 (i.e. 0.65 collisions per month). The collisions decreased in 16 months while it increased in 10 months. It did not change in 8 months.  Phase one and two: The effect of phase one and two together yielded a 12% decrease on average for 43 months (i.e. December 2012 to June 2016). The mean of monthly number of collisions reduced from 4.12 to 3.65 (i.e. 0.47 collisions per month). The collisions declined in 20 months but increased in 16 months.  Phase two: The collision counts raised by 10% on average during 34 months after the installation of the second phase. The mean of monthly number of collisions increased from 3.26 to 3.65. It appears that the installation of phase one (two VMSs) was successful in collision reduction, but phase two approach was failed to decrease. Table 6.32 presents the seasonal effects of the intervention(s). Interestingly, the installation of phase one, phase two, and phase one and two together were substantially effective in spring and summer seasons; however, phase one, and phase one and two together failed to influence, or had a small influence, on collision counts in fall and winter. Moreover, phase two increased collision counts in fall and winter. This suggests that: 1) VMSs, green wave signs, and enforcement might be more visible in spring and summer, 2) Speeding might not be the main cause of collisions in fall and winter, or 3) Seasonal or general trend of traffic volume might have a significant effect. However, traffic volume was not available for every month and excluded from the explanatory variable of the DANN models. The scatter points in Figure 6.31 show the available monthly traffic volume that was collected from VMSs. 156 Table 6.32 Seasonal effectiveness in the number of TOTAL collisions due to the intervention(s). Negative percentage means that the collision count was reduced after intervention(s), and vise versa. Intervention Season Winter Spring Summer Fall Phase one 0% -7% -78% 0% Phases one and two 4% -20% -33% 2% Phase two 47% -7% -22% 43% Figure 6.31 Available monthly average daily traffic (ADT) for Springfield Road Corridor. The dashed lines represent the implementation date of the safety phases. Further investigation is therefore needed to analyze the reason of collision events on this road corridor so appropriate safety improvements can be recommended. One alternative is to adopt the complete street approach. This improvement will ensure that the local residents are able to have a safe, convenient access to their homes, and the road users are able to continue using this road in a safe and effective manner. 15,00017,50020,00022,50025,00027,50030,00032,50035,000Aug/2009 Jan/2011 May/2012 Oct/2013 Feb/2015 Jun/2016Traffic volume (ADT)Time (month)Phase one Phase two157 6.5 Summary In this chapter, the ARIMA and DANN modelling approaches were developed for three types of collision series in the City of Kelowna. Then, the results were compared and demonstrated that DANN approach provided substantially higher predictive performance than the ARIMA approach. In addition, the NARX network outperforms the NAR network for all collision types. Finally, the effect of safety management phases that were implemented on the Springfield Road Corridor were investigated by the DANN model. The results showed that in terms of monthly number of TOTAL collisions, phase one improved the safety of this road by 16% while phase two reduced the safety by 10%. 158 Chapter 7: Conclusion 7.1 Summary of the thesis research One important purpose of road safety studies is to monitor, describe, evaluate, and predict road safety over time. Time series analysis techniques take the autocorrelations between observations into considerations and provide reliable inferences. While time series models do this quite well, research to date has found that they have a common drawback. Namely, these models place restrictions on the residuals such as normality, stationarity, or homoscedasticity. Such assumptions may be violated in road safety studies due to dynamic fluctuation and the complex, nonlinear behavior of collision time series data. As a result, a more powerful model is needed to cope with those assumptions. This research set out to test whether artificial neural networks (ANNs) might hold the key to overcoming these shortcomings in the existing statistical collision series analysis techniques used to date. The preliminary results of the trained networks showed that despite the advantages of ANNs, there was a challenging issue that should be considered. ANNs required large training samples to learn features of the collision series. In this thesis research, a two-stage methodology was proposed to deal with the limited observations of collision data. First, a wavelet transform allowed a breakdown of a given collision series into multiple subseries that represented unique features of the dataset. Then, dynamic artificial neural networks (DANNs), where the output depends on the previous inputs or outputs, were used to model each subseries separately. This action enabled the networks to learn the features of limited noisy data points. Consequently, the sum of simulated values by the networks provided a highly accurate prediction on the testing set. 159 For the specific case of this thesis, three types of collision series were collected for the City of Kelowna, BC; each was decomposed into three subseries. In total, 54,000 DANNs were trained to find the best predictions for the subseries. Next, the predicted values of subseries were accumulated to yield the projection of the original collision series. Finally, the predictive performance of the proposed approach was compared to the autoregressive integrated moving average (ARIMA) model, as it is a conventional time series statistical analysis tool for road safety studies. This process is summarized in Figure 6.18, which is repeated below for convenience. Figure 6.18 Flowchart of dynamic artificial neural network (DANN) coupled with discrete wavelet transform (DWT) The results of models’ performance for the three collision types are summarized in Table 6.22, which is repeated below for convenience. Original time series (training samples) D1 DANN2 DANN1 Predict the original time series by adding up the predicted values of the subseries DANN3 D2 A2 Wavelet decomposition D෡1 D෡2 A෡2 160 Table 6.22 A comparison of models’ performances on testing set (July 2014 to June 2016) Collision type Model Performance measure R2 MAE MAPE MPE MSE TOTAL ARIMA 0.54 45.50 5.79% 4.25% 3638 ARIMAX 0.37 62.38 8.09% 7.26% 6226 NAR network 0.88 19.21 2.44% 0.75% 708 NARX network 0.98 9.52 1.18% -0.18% 140 PDO ARIMA 0.45 49.29 7.87% 6.19% 3889 ARIMAX 0.33 66.04 10.85% 10.32% 6499 NAR network 0.92 14.30 2.30% 1.40% 388 NARX network 0.98 8.90 1.38% 0.15% 109 KSI ARIMA 0.41 16.58 10.69% 9.07% 447 ARIMAX 0.34 14.50 8.81% 0.99% 288 NAR network 0.84 7.16 4.44% 1.13% 72 NARX network 0.91 5.34 3.26% 0.27% 42 These results have demonstrated that the coupled DANN and wavelet transform outperformed the traditional method. Moreover, the proposed methodology was used to train 3,600 networks in order to evaluate the safety phases installed on Springfield Road, Kelowna. The results demonstrated the applicability of this approach for the intervention analysis. 7.2 Contributions of this research This research has made three major contributions in the field of road safety time series analysis: (1) methodology to develop DANNs, (2) apply DANNs, and (3) case study application of DANNs to evaluate safety interventions. 7.2.1 Contribution 1 – Methodology to develop dynamic artificial neural networks Traditional dynamic models, such as the ARIMA model, state space model, and vector autoregressive (VAR) model, have been widely applied to collision predicting. The above models 161 yield good prediction results when the collision series is linear. However, no earlier studies have focused on exploring the nonlinear time series models. This research is the first to present a comprehensive methodology for developing a DANN in this field to achieve the first objective of this thesis, i.e. to develop a methodology to conduct DANN modelling to predict collision series in road safety analysis. On the other hand, the earlier limited road safety analyses using a static network, discussed in Section 2.6.1, have not deeply studied ANN. This thesis has set out a methodology to undertake the profound investigation of the network pre-training steps (involving preprocessing methods, optimum lag inputs selection, network architecture selection, and optimum number of hidden layers and neurons selection), training steps (involving learning algorithm selection, network parameters initialization, and stopping criteria selection for training process), and post-training steps (involving output denormalization and network performance measures selection). Thus, the proposed methodology can be easily applied to other studies that are suitable to use either static or dynamic ANNs. 7.2.2 Contribution 2 – Apply dynamic artificial neural networks Multiple data types were collected from various sources to meet the second objective of this thesis, i.e. to assemble the required datasets to apply a DANN approach and evaluate its performance. Dynamic models were built for both univariate and multivariate conditions to investigate the effect of the explanatory variable on the models’ performance. The findings of this research showed that: 1) In the case of decomposing a given collision series by the discrete wavelet transform (DWT), the DANN approach can be simply applicable for collision prediction models, 162 2) ARIMA models have shown dissimilar behaviors towards different collision series, discussed in Section 6.3.1.3, which depend on the dynamic behavior of the time series. However, the proposed methodology filled this knowledge gap and made better predictions where the ARIMA approach produced low performance, and 3) DANNs efficiently mapped the association between the dependent and independent variables, but the ARIMA with explanatory variable (ARIMAX) models were unable to enhance the predictive performance of the corresponding ARIMA models. 7.2.3 Contribution 3 – Case study application of dynamic artificial neural networks to evaluate safety interventions The traditional time series statistical analyses cannot detect the nonlinear and irregular behaviour of the micro-level case studies where the mean is low. This research attempted to fill this gap in the literature and explore the third objective of this thesis, i.e. to address the City of Kelowna’s need for improved evaluation tools. The proposed methodology was developed to evaluate the impact of safety interventions implemented on Springfield Road Corridor at various times. The DANNs were trained quite well on a short and noisy data frame of this case study. Consequently, the results of this research demonstrated the capability of the proposed approach to be used not only in macro-level time series analysis of road safety development, but also for micro-level. In addition, this approach can be straightforwardly implemented by the City of Kelowna to: 1) Investigate any safety interventions for other case studies without imposing any limitations from traditional techniques, 2) Move being reactive to proactive in terms of safety management strategy, and 163 3) Provide predictions for future trend of other road transport factors, such as vehicular traffic volume, and bus transit and parking space demand. 7.3 Limitations of this research This research may have encountered three limitations by not applying the following techniques: (1) sensitivity analysis, (2) popular generalization methods, and (3) unsupervised learning for clustering analysis. 7.3.1 Application of sensitivity analysis ANN as a black-box mechanism does not provide any insights regarding the structure of the function. A sensitivity analysis, either derivative or variance based, can extract understandable knowledge from the effect of each input vector on the model output. Nevertheless, nine explanatory variables, which defined as input datasets for the nonlinear autoregressive with external inputs (NARX) networks, illustrated low correlation with the collision series except one variable, which was the employment number (EMP) in the province of British Columbia. As a result, there were not two or more external inputs in the networks to compare their effects. 7.3.2 Application of popular generalization techniques The commonly used generalization methods to improve the predictive accuracy and stability of the learning algorithms in machine learning, such as bagging or k-fold cross-validation, often divide a given training dataset into multiple subsets. The network is trained for each new training set separately. The outputs are combined to generate better results. However, the order of observed measures is fixed in time series data; therefore, the training samples cannot be partitioned and fed 164 to different networks. On the other hand, the appropriate generalization methods for time series data, like walk-in forward, divide the dataset into multiple overlapping subsets. The training set moves forward for each training run. Nevertheless, the limited observations restricted the applicability of these methods in this research. 7.3.3 Application of unsupervised learning for clustering analysis Unsupervised learning can be used to model the underlying features of the data in order to have a better understanding about it. For example, clustering algorithms, such as k-means or hierarchical method, organize similar data points into a group. A short data frame in this research prevented from performing an unsupervised network to find data with similar features in conjunction with dynamic networks as an alternative to the wavelet transform. The reason lies in the fact that each cluster would contain a few data points that are insufficient to train a network. 7.4 Future research ANNs have a noteworthy contribution to the advancement of transportation engineering and road safety. The following recommendations are proposed as potential future research directions: 1) There is a significant overlap between conventional statistical models and machine learning (ML) algorithms used in predictive analytics in the world of data science. However, the current complexity of datasets due to digitization of engineering systems and easier methods of data collection, made ML algorithms extremely popular. The term “machine learning” is often misused by scholars as a technique that can do wonders or magics to a dataset since it does not require significant statistical assumptions prior to the learning process of algorithms. Considering the literature of engineering fields, one could notice how these techniques found 165 their ways as means of predictive analytics in the complex system. Although these methods can perform effectively, one can compare the strengths and weaknesses of different ML algorithms in their basic general concepts such as input design, over or under-fitting, emerging data mining, and deep learning concepts. It would be interesting to apply different so-called artificial intelligence-based ML algorithms such as support vector machines or genetic programming for a comparative assessment with ANNs used in this thesis. Also, this research deployed a wavelet decomposition technique as a strong data mining to improve the performance of networks. Other data mining approaches that are used in the field of ML are clustering, dimensionality reduction, anomaly detection, or reinforcement learning. However, many of these proposed methods are applicable to large datasets, which was the primary limitation of this thesis. 2) The behavior of a collision series may be affected by several events at various temporal resolutions like weekend effect, but collisions are usually aggregated and analyzed over a fixed time resolution like monthly in this thesis. Instead, the collision events can be characterized by multiple models representing separate series with different temporal resolutions. Some techniques such as a minimum description length procedure can extract the right resolutions and decompose a given collision series into its underlying components. 3) Confidence intervals quantify the uncertainty of the predicted values and show the promising precision of the estimations. The assumptions of traditional confidence interval methods are frequently violated by the ML algorithms. To improve the model results, the uncertainty level of the predictions should be quantified. There is no generalized method to estimate the confidence interval of ANN models, but some methods like delta, Bayesian, mean-variance 166 estimation, or the bootstrap resampling method were proposed in the literature. Hence, comprehensive research is needed to investigate and compare these methods. 4) Some tasks for developing DANNs should be done manually like selecting the optimal network configuration. As a result, the success of the models profoundly relies on the modeler’s judgement and expertise. Moreover, due to the complexity and required statistical knowledge of these tasks, it is often beyond the capability of inexperienced practitioners. An automated tool to design the DANN models and optimize their performance can help inexpert developers. It also reduces human errors and efforts. It should be noted that some techniques have been proposed to automate the ML design, but they are still at their early stages. 5) In this research, the multivariate dynamic models were developed using only one explanatory variable since other candidate variables didn’t characterize the dynamic patterns of the collision series. Long short-term memory networks can be tested as an alternative to model multiple input variables. This network is a special type of recurrent neural network that can not only learn long-term dependencies, but also stores short-memories. The network decides what to remember and therefore new data will not affect valuable information. 6) Develop a model using data points between interventions one and two to better explain the DANN results for post-intervention two. 7) Measure traffic volume and add to the multivariate modeling approaches as an explanatory variable to represent the real exposure data. 167 Bibliography [1] World Health Organization, Global Status Report on Road Safety 2015. Geneva, Switzerland: World Health Organization, 2015. [2] United Nations General Assembly, "Improving global road safety," Geneva, Switzerland, Tech. Rep. A/RES/64/255, 2010. [3] S. Lassarre, "Analysis of progress in road safety in ten European countries," Accident Analysis and Prevention, vol. 33, (6), pp. 743-751, 2001. DOI: 10.1016/S0001-4575(00)00088-9. [4] R. Elvik, "The stability of long-term trends in the number of traffic fatalities in a sample of highly motorised countries," Accident Analysis and Prevention, vol. 42, (1), pp. 245-260, 2010. DOI: 10.1016/j.aap.2009.08.002. [5] Transport Canada, "Canadian motor vehicle traffic collision statistics 2014," Tech. Rep. T45-3E-PDF, 2014. [6] Insurance Corporation of British Columbia. "Quick Statistics", 2016. [online]. Available: http://www.icbc.com/about-icbc/newsroom/Documents/quick-statistics.pdf. [7] E. Hauer, Observational before-After Studies in Road Safety: Estimating the Effect of Highway and Traffic Engineering Measures on Road Safety. Oxford, England: Pergamon Press, 1997. [8] S. Hakim et al, "A critical review of macro models for road accidents," Accident Analysis and Prevention, vol. 23, (5), pp. 379-400, 1991. DOI: 10.1016/0001-4575(91)90058-D. [9] M. Pickup, Introduction to Time Series Analysis (Quantitative Applications in the Social Sciences). Los Angeles, California: SAGE, 2015. [10] P. S. P. Cowpertwait and A. V. Metcalfe, Introductory Time Series with R. 2009. 168 [11] R. H. Shumway and D. S. Stoffer, Time Series Analysis and its Applications: With R Examples. (3rd ed.) 2011. [12] A. K. Palit and D. Popovic, Computational Intelligence in Time Series Forecasting: Theory and Engineering Applications. (1st ed.) 2005. DOI: 10.1007/1-84628-184-9. [13] J. J. F. Commandeur, S. J. Koopman and MyiLibrary, An Introduction to State Space Time Series Analysis. 2007. [14] S. Chatterjee and J. S. Simonoff, Handbook of Regression Analysis. (1st ed.) 2013. [15] C. B. Antoniou et al, "Multilevel modelling and time series analysis in traffic safety research: methodology," 2007. [16] J. J. Wiorkowski and R. F. Heckard, "The use of time series analysis and intervention analysis to assess the effects of external factors on traffic indices: a case study of the effects of the speed limit reduction and energy crisis in the State of Texas," Accident Analysis and Prevention, vol. 9, (4), pp. 229-247, 1977. [17] D. Karlis and E. Hermans, "Time series models for road safety accident prevention," Steunpunt Mobiliteit & Openbare Werken – Spoor Verkeersveiligheid, Tech. Rep. RA-MOW-2011-020, 2012. [18] V. S. Ellingstad and D. P. Westra, "Evaluation methodologies for traffic safety programs," Hum. Factors, vol. 18, (4), pp. 313-325, 1976. [19] S. Hess and J. Polak, "Effects of Speed Limit Enforcement Cameras on Accident Rates," Transportation Research Record: Journal of the Transportation Research Board, vol. 1830, pp. 25-33, 2003. 169 [20] V. Serhiyenko et al, "Dynamic compositional modeling of pedestrian crash counts on urban roads in Connecticut," Accident Analysis and Prevention, vol. 64, pp. 78-85, 2014. DOI: 10.1016/j.aap.2013.11.006. [21] A. Tolón-Becerra, X. Lastra-Bravo and I. Flores-Parra, "National and regional analysis of road accidents in Spain," Traffic Injury Prevention, vol. 14, (5), pp. 486-495, 2013. [22] G. V. Glass, "Analysis of Data on the Connecticut Speeding Crackdown as a Time-Series Quasi-Experiment," Law & Society Review, vol. 3, (1), pp. 55-76, 1968. [23] D. T. Campbell and H. L. Ross, "The Connecticut Crackdown on Speeding: Time-Series Data in Quasi-Experimental Analysis," Law & Society Review, vol. 3, (1), pp. 33-53, 1968. [24] G. V. Glass, V. L. Willson and J. M. Gottman, Design and Analysis of Time-Series Experiments. Boulder, CO: Colorado Associated University Press, 1975. [25] G. M. J. G.E.P. Box, Time Series Analysis: Forecasting and Control. San Francisco, CA: Holden Day, 1976. [26] G. E. P. Box et al, Time Series Analysis: Forecasting and Control. (4th ed.) 2008. [27] J. Commandeur et al, "On statistical inference in time series analysis of the evolution of road safety," Accident Analysis and Prevention, vol. 60, pp. 424-434, 2013. DOI: 10.1016/j.aap.2012.11.006. [28] P. P. Scott, "Modelling time-series of British road accident data," Accident Analysis and Prevention, vol. 18, (2), pp. 109-117, 1986. DOI: 10.1016/0001-4575(86)90055-2. [29] M. Ghiassi, H. Saidane and D. K. Zimbra, "A dynamic artificial neural network model for forecasting time series events," Int. J. Forecast., vol. 21, (2), pp. 341-362, 2005. DOI: 10.1016/j.ijforecast.2004.10.008. 170 [30] P. Michalaki et al, "A time-series analysis of motorway collisions in England considering road infrastructure, socio-demographics, traffic and weather characteristics," Journal of Transport & Health, vol. 3, (1), pp. 9-20, 2016. [31] M. Quddus, "Time series count data models: An empirical application to traffic accidents," Accident Analysis and Prevention, vol. 40, (5), pp. 1732-1741, 2008. DOI: 10.1016/j.aap.2008.06.011. [32] A. C. Harvey and J. Durbin, "The Effects of Seat Belt Legislation on British Road Casualties: A Case Study in Structural Time Series Modelling," Journal of the Royal Statistical Society (Serie A), vol. 149, (3), pp. 187-227, 1986. [33] A. C. Harvey and C. Fernandes, "Time series models for insurance claims," Journal of the Institute of Actuaries, vol. 116, (3), pp. 513-528, 1989. [34] R. E. Kalman, "A New Approach to Linear Filtering and Prediction Problems," Journal of Basic Engineering, vol. 82, (1), pp. 35, 1960. [35] P. A. Scuffham and J. D. Langley, "A model of traffic crashes in New Zealand," Accident Analysis and Prevention, vol. 34, (5), pp. 673-687, 2002. DOI: 10.1016/S0001-4575(01)00067-7. [36] S. Nghiem, J. J. F. Commandeur and L. B. Connelly, "Determinants of road traffic safety: New evidence from Australia using state-space analysis," Accident Analysis and Prevention, vol. 94, (September), pp. 65-72, 2016. [37] F. Bijleveld et al, "Model-Based Measurement of Latent Risk in Time Series with Applications," Journal of the Royal Statistical Society (Serie A), vol. 171, (1), pp. 265-277, 2008. DOI: 10.1111/j.1467-985X.2007.00496.x. 171 [38] F. Bijleveld et al, "Multivariate Non-Linear Time Series Modelling of Exposure and Risk in Road Safety Research," Journal of the Royal Statistical Society, vol. 59, (1), pp. 145-161, 2010. DOI: 10.1111/j.1467-9876.2009.00690.x. [39] C. Antoniou and G. Yannis, "State-space based analysis and forecasting of macroscopic road safety trends in Greece," Accident Analysis Prevention, vol. 60, pp. 268-276, 2013. DOI: 10.1016/j.aap.2013.02.039. [40] J. Commandeur and R. Bergel, "Equivalence between ARIMA and state space models," in Multilevel Modelling and Time Series Analysis in Traffic Research: Methodology, E. Dupont and H. Martensen, Eds. 2007, pp. 337-400. [41] D. KARLIS and E. HERMANS, "Time series models for road safety accident prediction," Steunpunt Mobiliteit & Openbare Werken – Spoor Verkeersveiligheid, Tech. Rep. RA-MOW-2011-020, 2012. [42] H. Lütkepohl, New Introduction to Multiple Time Series Analysis. Berlin Heidelburg: Springer-Verlag, 2005. [43] P. Michalaki et al, "A time-series analysis of motorway collisions in England considering road infrastructure, socio-demographics, traffic and weather characteristics," Journal of Transport & Health, vol. 3, (1), pp. 9-20, 2016. [44] N. J. Meyerhoff, "The influence of daylight saving time on motor vehicle fatal traffic accidents," Accident Analysis and Prevention, vol. 10, (3), pp. 207-221, 1978. [45] A. M. Novoa et al, "Road safety in the political agenda: the impact on road traffic injuries," Journal of Epidemiology and Community Health, vol. 65, (3), pp. 218-225, 2010. 172 [46] G. Yannis, C. Antoniou and E. Papadimitriou, "Autoregressive nonlinear time-series modeling of traffic fatalities in Europe," European Transport Research Review, vol. 3, (3), pp. 113-127, 2011. DOI: 10.1007/s12544-011-0055-4. [47] S. S. Haykin, Neural Networks: A Comprehensive Foundation. New Jersey: Prentice Hall, 1994. [48] L. Erguo and Y. Jinshou, "An input-training neural network based nonlinear principal component analysis approach for fault diagnosis," in Proceedings of the 4th World Congress on Intelligent Control and Automation, 2002. [49] W. S. McCulloch and W. Pitts, "A logical calculus of the ideas immanent in nervous activity," Bulletin of Mathematical Biophysics., vol. 5, (4), pp. 115-133, 1943. [50] D. O. Hebb, The Organization of Behavior: A Neuropsychological Theory. New York: Wiley, 1949. [51] F. Rosenblatt, "The perceptron: A probabilistic model for information storage and organization in the brain," Psychol. Rev., vol. 65, (6), pp. 386-408, 1958. [52] B. Widrow and M. E. Hoff, "Adaptive switching circuits," Adaptive Switching Circuits, (1960 IRE WESCON Convention Record, New York: IRE Part 4), pp. 96-104, 1960. [53] M. Minsky and S. Papert, Perceptrons; an Introduction to Computational Geometry. Cambridge, MA: MIT Press, 1969. [54] P. J. Werbos, "Beyond regression: New tools for prediction and analysis in the behavioral sciences," Doctoral Dissertation, Applied Mathematics, Harvard University, Cambridge, MA, 1974. 173 [55] D. B. Parker, "Learning-logic: Casting the cortex of the human brain in silicon," Center for Computational Research in Economics and Management Science, MIT, Cambridge, MA, Tech. Rep. TR-47, 1985. [56] D. E. Rumelhart and J. L. McClelland, Parallel Distributed Processing: Explorations in the Microstructure of Cognition. Cambridge, MA: MIT Press, 1986. [57] H. B. Demuth et al, Neural Network Design. Boulder, CO: Campus Pub. Service, Colorado University Bookstore, 2014. [58] L. Chang, "Analysis of freeway accident frequencies: Negative binomial regression versus artificial neural network," Safety Science, vol. 43, (8), pp. 541-557, 2005. DOI: 10.1016/j.ssci.2005.04.004. [59] Y. Chiou, "An artificial neural network-based expert system for the appraisal of two-car crash accidents," Accident Analysis and Prevention, vol. 38, (4), pp. 777-785, 2006. [60] A. Pande and M. Abdel-Aty, "A computing approach using probabilistic neural networks for instantaneous appraisal of rear-end crash risk," Computer-Aided Civil and Infrastructure Engineering, vol. 23, (7), pp. 549-559, 2008. [61] Z. Wenhui, L. Shengqin and W. Biao, "Safety evaluation of traffic accident scene based on artificial neural network," in 2009. [62] M. De Luca, "A comparison between prediction power of artificial neural networks and multivariate analysis in road safety management," Transport, pp. 1-7, 2015. DOI: 10.3846/16484142.2014.995702. [63] H. Abdelwahab and M. Abdel-Aty, "Development of artificial neural network models to predict driver injury severity in traffic accidents at signalized intersections," Transportation Research Record, (1746), pp. 6-13, 2001. 174 [64] D. Delen, R. Sharda and M. Bessonov, "Identifying significant predictors of injury severity in traffic accidents using a series of artificial neural networks," Accident Analysis and Prevention, vol. 38, (3), pp. 434-444, 2006. DOI: 10.1016/j.aap.2005.06.024. [65] D. Liu et al, "Classifying Injuries Occurrence in Motor Vehicles Collisions Using Artificial Neural Network," Proceedings of the Human Factors and Ergonomics Society Annual Meeting, vol. 55, (1), pp. 1808-1812, 2011. [66] Q. Zeng et al, "Developing an optimized artificial neural network to predict traffic crash injury severity," in CICTP 2014: Safe, Smart, and Sustainable Multimodal Transportation SystemsAnonymous 2014, pp. 2396-2407. [67] S. Alkheder, M. Taamneh and S. Taamneh, "Severity Prediction of Traffic Accident Using an Artificial Neural Network: Traffic Accident Severity Prediction Using Artificial Neural Network," Journal of Forecasting, vol. 36, (1), pp. 100-108, 2017. [68] M. Sameen and B. Pradhan, "Severity Prediction of Traffic Accidents with Recurrent Neural Networks," Applied Sciences, vol. 7, (6), pp. 476, 2017. [69] D. Graupe, Principles of Artificial Neural Networks. (2nd ed.) Hackensack, Singapore: World Scientific, 2007. [70] R. Rojas, Neural Networks: A Systematic Introduction. 1996. [71] D. E. Rumelhart, G. E. Hinton and R. J. Williams, "Learning representations by back-propagating errors," Nature, vol. 323, (6088), pp. 533-536, 1986. [72] M. Qi and G. P. Zhang, "An investigation of model selection criteria for neural network time series forecasting," European Journal of Operational Research, vol. 132, (3), pp. 666-680, 2001. DOI: 10.1016/S0377-2217(00)00171-5. 175 [73] S. Xu and L. Chen, "A novel approach for determining the optimal number of hidden layer neurons for FNN's and its application in data mining," in 5th International Conference on Information Technology and Applications (ICITA 2008), 2008, pp. 683-686. [74] L. Breiman, "Bagging predictors," Department of Statistics, University of California, Berkeley, CA, Tech. Rep. 421, 1994. [75] R. Kohavi, "A study of cross-validation and boostrap for accuracy estimation and model selection," in Proceedings of the Fourteenth International Joint Conference on Artificial Intelligence, San Mateo, CA: Morgan Kaufmann, 1995, pp. 1137–1143. [76] G. M. Ljung and G. E. P. Box, "On a Measure of Lack of Fit in Time Series Models," Biometrika, vol. 65, (2), pp. 297, 1978. [77] C. Kongcharoen and T. Kruangpradit, "Autoregressive integrated moving average with explanatory variable (ARIMAX) model for Thailand export," in 33rd International Symposium on Forecasting, South Korea, 2013, pp. 1-8. [78] C. Chen and L. Liu, "Joint Estimation of Model Parameters and Outlier Effects in Time Series," Journal of the American Statistical Association, vol. 88, (421), pp. 284-297, 1993. [79] R. K. Pearson et al, "Generalized Hampel Filters," EURASIP Journal on Advances in Signal Processing, vol. 2016, (1), pp. 1-18, 2016. [80] J. Adamowski and H. F. Chan, "A wavelet neural network conjunction model for groundwater level forecasting," Journal of Hydrology, vol. 407, (1), pp. 28-40, 2011. [81] M. Misiti et al, Wavelet Toolbox for use with MATLAB. Natick, MA: The MathWorks Inc., 1996. [82] S. G. Mallat, A Wavelet Tour of Signal Processing. (2nd ed.) San Diego: Academic Press, 1999. 176 [83] P. Areekul et al, "Neural-wavelet Approach for Short Term Price Forecasting in Deregulated Power Market," Journal of International Council on Electrical Engineering, vol. 1, (3), pp. 331-338, 2011. [84] S. G. Mallat, "A theory for multiresolution signal decomposition: the wavelet representation," IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, (7), pp. 674-693, 1989. [85] V. Nourani, M. T. Alami and M. H. Aminfar, "A combined neural-wavelet model for prediction of Ligvanchai watershed precipitation," Eng Appl Artif Intell, vol. 22, (3), pp. 466-472, 2009. [86] I. A. Basheer and M. Hajmeer, "Artificial neural networks: fundamentals, computing, design, and application," J. Microbiol. Methods, vol. 43, (1), pp. 3-31, 2000. [87] R. J. May et al, "Non-linear variable selection for artificial neural networks using partial mutual information," Environmental Modelling & Software, vol. 23, (10), pp. 1312-1326, 200. [88] R. May, G. Dandy and H. Maier, "Review of input variable selection methods for artificial neural networks," in Artificial Neural Networks - Methodological Advances and Biomedical Applications, K. Suzuki, Ed. InTech, Rijeka, Croatia, 2011, pp. 19-44. [89] B. V. Bonnlander and A. S. Weigend, "Selecting input variables using mutual information and nonparametric density estimation," in Proceedings of the 1994 International Symposium on Artificial Neural Networks, 1994, pp. 42-50. [90] K. P. Sudheer, A. K. Gosain and K. S. Ramasastri, "A data-driven algorithm for constructing artificial neural network rainfall-runoff models," Hydrol. Process., vol. 16, (6), pp. 1325-1330, 2002. [91] F. Takens, "Detecting strange attractors in turbulence," Lecture Notes in Mathematics, vol. 898, (1), pp. 366-381, 1981. 177 [92] I. Drezga and S. Rahman, "Input variable selection for ANN-based short-term load forecasting," IEEE Trans. Power Syst., vol. 13, (4), pp. 1238-1244, 1998. [93] D. Nguyen and B. Widrow, "Improving the learning speed of 2-layer neural networks by choosing initial values of the adaptive weights," in Proceedings of the IJCNN, 1990, DOI: 10.1109/IJCNN.1990.137819. [94] R. Adhikari and R. K. Agrawal, An Introductory Study on Time Series Modeling and Forecasting. Germany: Lambert Academic Publishing, 2013. [95] P. D. Sreekanth et al, "Forecasting groundwater level using artificial neural networks," Curr. Sci., vol. 96, (7), pp. 933-939, 2009. [96] A. C. Wagenaar, "Effects of macroeconomic conditions on the incidence of motor vehicle accidents," Accident Analysis and Prevention, vol. 16, (3), pp. 191-205, 1984. [97] D. S. Stoffer, "astsa: Applied statistical time series analysis. RStudio package version 3.3.3." 2014, 2017. [98] C. C. Ihueze and U. O. Onwurah, "Road traffic accidents prediction modelling: An analysis of Anambra State, Nigeria," Accident Analysis and Prevention, vol. 112, pp. 21-29, 2018. [99] Anonymous "MATLAB and Statistics Toolbox Release R2016b," 2016. [100] L. Hamm, B. W. Brorsen and M. T. Hagan, "Comparison of Stochastic Global Optimization Methods to Estimate Neural Network Weights," Neural Process Letters, vol. 26, (3), pp. 145-158, 2007. 178 Appendices Appendix A Average number of collisions for each month Table A.1 presents the average number of collisions for each month over 13.5 years for different collision series in the City of Kelowna. Table A.1 Average number of collisions for each month over 13.5 years for different collision series in the City of Kelowna Collision type Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sept. Oct. Nov. Dec. TOTAL 680 570 582 622 702 696 754 719 686 687 696 789 PDO 547 456 462 488 551 546 598 558 524 525 537 630 KSI 133 114 121 134 151 150 155 161 162 161 160 159 179 Appendix B Dynamic artificial neural network modelling code This appendix presents a part of the DANN modelling code developed in the MATLAB R2016b software program. clear, clc % DEFINE COLLISION SERIES collisions = [dataset]; % A=[maxNumDelays numHiddenNeurons numTrainingNetworks] A = [12 10 10]; % CONVERT NUMERIC ARRAY TO CELL ARRAY T = num2cell(collisions); % NUMBER OF DELAYS j = 0; for delay = 1 : A(1) % NUMBER OF HIDDEN NEURONS for numHiddenNeurons = 1 : A(2) % NUMBER OF TRAINING NETWORKS for run = 1 : A(3) % DEFINE NETWORK ARCHITECTURE net = narnet(1:delay,numHiddenNeurons); %narnet(inputDelays,numHiddenNeurons) % DEFINE TRANSFER FUNCTIONS FOR HIDDEN layer AND OUTPUT LAYER net.layers{1}.transferFcn = 'tansig'; net.layers{2}.transferFcn = 'purelin'; 180 % DIVIDE DATA INTO THREE CONTIGUOUS BLOCKS net.divideFcn = 'divideind'; % SET UP DIVISION OF TRAINING, VALIDATING, AND TESTING SETS net.divideParam.trainInd = 1:114-delay; net.divideParam.valInd = 114-delay+1:138-delay; % PREPARE INPUT AND TARGET TIME SERIES DATA FOR NETWORK TRAINING [Xs,Xi,Ai,Ts] = preparets(net,{},{},T); %[inputs,inputStates,layerStates,targets] = preparets(net,{},{},target) % TRAIN NETWORK [net,tr] = train(net,Xs,Ts,Xi,Ai); % VIEW NETWORK ARCHITECTURE view(net) % CALCULATE NETWORK OUTPUT Y = net(Xs,Xi,Ai); % DEFINE TRAINING, VALIDATING, AND TESTING TARGET SETS j = j + 1; targetTrain = cell2mat(Ts(1:length(tr.trainInd))); targetTrainMatrix(j,1:length(tr.trainInd)) = targetTrain; targetVal = cell2mat(Ts(length(tr.trainInd)+1:length(tr.trainInd)+length(tr.valInd))); targetValMatrix(j,1:length(tr.valInd)) = targetVal; targetTest = cell2mat(Ts(length(tr.trainInd)+length(tr.valInd)+1:end)); targetTestMatrix(j,1:length(tr.testInd)) = targetTest; 181 % DEFINE TRAINING AND VALIDATING PREDICTED SETS predictTrain = cell2mat(Y(1:length(tr.trainInd))); predictTrainMatrix(j,1:length(tr.trainInd)) = predictTrain; predictVal = cell2mat(Y(length(tr.trainInd)+1:length(tr.trainInd)+length(tr.valInd))); predictValMatrix(j,1:length(tr.valInd)) = predictVal; % CALCULATE ERRORS OF TRAINING and VALIDATING SETS errorTrain = targetTrain - predictTrain; errorVal = targetVal - predictVal; % CORRELATION OF DETERMINATION (R2) OF TRAINING AND VALIDATING SETS R = corrcoef(predictTrain,targetTrain); performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+1) = (R(2))^2; R = corrcoef(predictVal,targetVal); performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+2) = (R(2))^2; % MEAN ABSOLUTE ERROR (MAE)OF TRAINING and VALIDATING SETS performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+4) = mae(targetTrain,predictTrain); performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+5) = mae(targetVal,predictVal); % MEAN ABSOLUTE PERCENTAGE ERROR (MAPE) OF TRAINING and VALIDATING SETS performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+7) = (mean(abs(errorTrain)/abs(targetTrain)))*100; 182 performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+8) = (mean(abs(errorVal)/abs(targetVal)))*100; % MEAN PERCENTAGE ERROR (MPE) OF TRAINING and VALIDATING SETS performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+10) = (mean(errorTrain/targetTrain))*100; performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+11) = (mean(errorVal/targetVal))*100; % MEAN SQUARED ERROR (MSE) OF TRAINING AND VALIDATING SETS performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+13) = immse(targetTrain,predictTrain); performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+14) = immse(targetVal,predictVal); % ROOT MEAN SQUARED ERROR (RMSE) OF TRAINING AND VALIDATING SETS performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+16) = sqrt(immse(targetTrain,predictTrain)); performance(delay*A(3)-A(3)+run,numHiddenNeurons*18-18+17) = sqrt(immse(targetVal,predictVal)); 183 Appendix C Preliminary results of DANN modelling This appendix presents the R2 values corresponding to the validating sets of the networks that were trained on the raw data of collision series. 184 Table C.1 Coefficient of determination of NAR networks for TOTAL collisions in Kelowna Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 0.13 0.14 0.12 0.13 0.14 0.12 0.11 0.12 0.15 0.15 2 0.13 0.13 0.09 0.16 0.14 0.10 0.12 0.17 0.14 0.11 3 0.10 0.08 0.17 0.16 0.07 0.07 0.12 0.13 0.20 0.19 4 0.28 0.07 0.07 0.09 0.10 0.23 0.20 0.21 0.16 0.11 5 0.09 0.08 0.09 0.17 0.18 0.11 0.14 0.26 0.24 0.28 6 0.05 0.13 0.21 0.19 0.35 0.18 0.16 0.19 0.21 0.04 7 0.08 0.20 0.08 0.28 0.30 0.07 0.05 0.13 0.05 0.20 8 0.07 0.34 0.09 0.32 0.39 0.07 0.37 0.28 0.23 0.23 9 0.07 0.04 0.01 0.07 0.28 0.33 0.35 0.33 0.38 0.15 10 0.20 0.37 0.32 0.28 0.38 0.40 0.17 0.31 0.41 0.43 11 0.64 0.42 0.40 0.27 0.49 0.44 0.46 0.29 0.43 0.47 12 0.43 0.53 0.42 0.44 0.24 0.44 0.38 0.36 0.23 0.52 2 1 0.15 0.15 0.08 0.12 0.10 0.16 0.10 0.09 0.15 0.15 2 0.11 0.14 0.17 0.11 0.19 0.11 0.20 0.12 0.04 0.14 3 0.10 0.13 0.15 0.08 0.06 0.09 0.17 0.15 0.11 0.12 4 0.03 0.10 0.05 0.05 0.22 0.27 0.18 0.26 0.20 0.18 5 0.09 0.10 0.09 0.10 0.20 0.16 0.13 0.21 0.28 0.20 6 0.05 0.09 0.18 0.21 0.06 0.32 0.34 0.27 0.37 0.22 7 0.07 0.01 0.20 0.16 0.26 0.21 0.23 0.19 0.35 0.26 8 0.00 0.00 0.37 0.07 0.00 0.20 0.10 0.18 0.07 0.35 9 0.07 0.10 0.26 0.11 0.39 0.15 0.11 0.07 0.36 0.27 10 0.01 0.17 0.00 0.17 0.19 0.22 0.30 0.35 0.49 0.26 11 0.10 0.37 0.50 0.40 0.32 0.37 0.04 0.40 0.26 0.54 12 0.36 0.00 0.32 0.32 0.22 0.40 0.47 0.41 0.18 0.31 3 1 0.13 0.03 0.12 0.18 0.13 0.13 0.13 0.13 0.13 0.11 2 0.00 0.11 0.17 0.06 0.12 0.16 0.20 0.17 0.03 0.21 3 0.04 0.05 0.02 0.01 0.17 0.20 0.15 0.08 0.15 0.26 4 0.03 0.05 0.12 0.23 0.05 0.14 0.29 0.11 0.17 0.17 5 0.13 0.11 0.20 0.02 0.15 0.22 0.23 0.18 0.26 0.17 6 0.08 0.02 0.06 0.10 0.22 0.17 0.25 0.35 0.12 0.25 7 0.00 0.02 0.06 0.16 0.19 0.29 0.07 0.16 0.23 0.16 8 0.00 0.14 0.24 0.40 0.28 0.01 0.25 0.24 0.00 0.27 9 0.06 0.02 0.01 0.37 0.05 0.03 0.33 0.09 0.09 0.16 10 0.05 0.13 0.00 0.34 0.07 0.46 0.48 0.24 0.27 0.12 11 0.00 0.52 0.27 0.39 0.00 0.48 0.61 0.33 0.48 0.28 12 0.40 0.53 0.00 0.55 0.46 0.23 0.62 0.33 0.47 0.51 4 1 0.12 0.02 0.18 0.14 0.11 0.13 0.12 0.16 0.16 0.14 2 0.14 0.16 0.12 0.14 0.14 0.14 0.15 0.15 0.20 0.01 3 0.08 0.00 0.17 0.13 0.06 0.18 0.09 0.14 0.02 0.00 4 0.04 0.12 0.04 0.09 0.11 0.03 0.00 0.07 0.27 0.14 5 0.10 0.01 0.15 0.20 0.32 0.00 0.16 0.27 0.12 0.14 6 0.00 0.15 0.15 0.00 0.24 0.08 0.25 0.14 0.31 0.09 7 0.06 0.04 0.05 0.03 0.00 0.05 0.30 0.17 0.00 0.25 8 0.00 0.00 0.02 0.04 0.06 0.20 0.00 0.00 0.02 0.07 9 0.02 0.02 0.01 0.04 0.16 0.08 0.16 0.16 0.20 0.25 10 0.00 0.00 0.10 0.09 0.28 0.38 0.13 0.43 0.40 0.19 11 0.58 0.12 0.67 0.19 0.47 0.21 0.13 0.06 0.65 0.47 12 0.00 0.41 0.54 0.00 0.47 0.05 0.30 0.42 0.18 0.33 185 Table C.2 Coefficient of determination of NAR networks for PDO collisions in Kelowna Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 0.06 0.04 0.08 0.09 0.08 0.07 0.10 0.10 0.07 0.09 2 0.06 0.07 0.05 0.10 0.09 0.12 0.11 0.08 0.07 0.12 3 0.08 0.11 0.09 0.04 0.01 0.03 0.05 0.08 0.08 0.10 4 0.02 0.00 0.01 0.05 0.01 0.05 0.00 0.00 0.00 0.04 5 0.05 0.07 0.11 0.00 0.02 0.12 0.02 0.04 0.22 0.00 6 0.06 0.07 0.11 0.03 0.11 0.26 0.17 0.08 0.08 0.29 7 0.00 0.10 0.21 0.09 0.07 0.06 0.06 0.22 0.26 0.01 8 0.18 0.17 0.13 0.01 0.14 0.03 0.05 0.09 0.00 0.06 9 0.04 0.04 0.12 0.34 0.08 0.15 0.32 0.01 0.14 0.15 10 0.11 0.11 0.25 0.19 0.13 0.25 0.32 0.07 0.32 0.04 11 0.05 0.15 0.19 0.11 0.29 0.40 0.22 0.16 0.20 0.33 12 0.25 0.21 0.27 0.28 0.24 0.33 0.25 0.19 0.24 0.23 2 1 0.06 0.03 0.08 0.07 0.06 0.06 0.07 0.05 0.04 0.09 2 0.06 0.04 0.09 0.01 0.05 0.13 0.23 0.18 0.18 0.11 3 0.06 0.01 0.01 0.04 0.05 0.07 0.11 0.02 0.04 0.11 4 0.03 0.04 0.01 0.02 0.00 0.03 0.32 0.04 0.11 0.08 5 0.09 0.12 0.10 0.14 0.18 0.05 0.06 0.01 0.02 0.00 6 0.07 0.01 0.01 0.11 0.20 0.29 0.02 0.03 0.08 0.14 7 0.05 0.12 0.08 0.09 0.03 0.02 0.04 0.07 0.22 0.22 8 0.10 0.18 0.07 0.11 0.09 0.07 0.14 0.03 0.21 0.15 9 0.02 0.01 0.01 0.00 0.00 0.02 0.03 0.14 0.12 0.01 10 0.01 0.26 0.05 0.22 0.03 0.14 0.20 0.06 0.20 0.08 11 0.00 0.33 0.12 0.33 0.09 0.35 0.21 0.47 0.00 0.21 12 0.00 0.08 0.31 0.31 0.21 0.23 0.38 0.49 0.19 0.19 3 1 0.05 0.03 0.00 0.07 0.08 0.03 0.07 0.09 0.07 0.03 2 0.08 0.00 0.02 0.00 0.08 0.24 0.19 0.02 0.18 0.22 3 0.01 0.06 0.07 0.04 0.12 0.02 0.03 0.03 0.05 0.21 4 0.13 0.05 0.00 0.01 0.03 0.04 0.01 0.18 0.00 0.04 5 0.03 0.11 0.00 0.04 0.02 0.03 0.02 0.06 0.28 0.04 6 0.00 0.10 0.07 0.00 0.00 0.32 0.09 0.03 0.11 0.16 7 0.00 0.03 0.00 0.00 0.14 0.09 0.11 0.01 0.24 0.04 8 0.00 0.00 0.10 0.12 0.02 0.07 0.06 0.08 0.00 0.22 9 0.01 0.04 0.20 0.09 0.06 0.01 0.10 0.04 0.13 0.00 10 0.00 0.02 0.01 0.00 0.05 0.00 0.36 0.32 0.16 0.00 11 0.26 0.03 0.62 0.27 0.05 0.44 0.30 0.03 0.15 0.36 12 0.23 0.21 0.35 0.26 0.26 0.11 0.21 0.07 0.33 0.35 4 1 0.06 0.06 0.04 0.00 0.08 0.07 0.10 0.07 0.05 0.06 2 0.00 0.09 0.07 0.07 0.09 0.17 0.04 0.16 0.06 0.20 3 0.01 0.01 0.09 0.07 0.03 0.00 0.03 0.07 0.01 0.01 4 0.00 0.05 0.08 0.04 0.00 0.03 0.00 0.02 0.20 0.09 5 0.00 0.08 0.17 0.13 0.03 0.15 0.02 0.01 0.01 0.00 6 0.00 0.05 0.00 0.03 0.29 0.02 0.05 0.10 0.10 0.05 7 0.00 0.14 0.04 0.04 0.15 0.04 0.01 0.10 0.01 0.30 8 0.15 0.02 0.05 0.00 0.23 0.24 0.15 0.03 0.11 0.16 9 0.04 0.07 0.20 0.02 0.01 0.05 0.32 0.00 0.12 0.21 10 0.00 0.00 0.06 0.00 0.04 0.00 0.10 0.35 0.21 0.01 11 0.00 0.13 0.15 0.13 0.40 0.17 0.28 0.00 0.31 0.05 12 0.05 0.00 0.17 0.22 0.20 0.28 0.24 0.08 0.18 0.00 186 Table C.3 Coefficient of determination of NAR networks for KSI collisions in Kelowna Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 0.37 0.41 0.42 0.42 0.46 0.39 0.33 0.56 0.45 0.34 2 0.43 0.39 0.38 0.39 0.30 0.33 0.21 0.29 0.15 0.34 3 0.49 0.52 0.50 0.43 0.25 0.45 0.43 0.33 0.39 0.40 4 0.54 0.54 0.53 0.43 0.47 0.42 0.58 0.50 0.48 0.41 5 0.34 0.25 0.51 0.44 0.53 0.58 0.48 0.35 0.40 0.48 6 0.55 0.58 0.36 0.40 0.54 0.26 0.29 0.38 0.48 0.24 7 0.51 0.49 0.37 0.55 0.44 0.41 0.21 0.41 0.47 0.29 8 0.54 0.34 0.44 0.52 0.35 0.55 0.18 0.49 0.43 0.64 9 0.47 0.45 0.39 0.51 0.54 0.43 0.27 0.47 0.01 0.22 10 0.61 0.49 0.57 0.43 0.58 0.36 0.43 0.43 0.33 0.62 11 0.56 0.59 0.55 0.47 0.46 0.34 0.54 0.44 0.53 0.33 12 0.44 0.54 0.46 0.52 0.64 0.33 0.42 0.35 0.42 0.29 2 1 0.39 0.43 0.44 0.40 0.58 0.46 0.27 0.60 0.38 0.49 2 0.00 0.39 0.31 0.00 0.23 0.21 0.34 0.32 0.36 0.40 3 0.50 0.46 0.48 0.48 0.37 0.48 0.49 0.23 0.25 0.40 4 0.54 0.43 0.51 0.49 0.48 0.48 0.39 0.59 0.30 0.57 5 0.33 0.35 0.16 0.53 0.45 0.50 0.56 0.37 0.35 0.45 6 0.00 0.38 0.49 0.55 0.36 0.44 0.42 0.22 0.38 0.49 7 0.00 0.29 0.56 0.49 0.29 0.24 0.44 0.42 0.30 0.34 8 0.54 0.09 0.57 0.35 0.33 0.28 0.42 0.38 0.33 0.30 9 0.39 0.05 0.40 0.41 0.56 0.23 0.44 0.59 0.19 0.37 10 0.52 0.55 0.50 0.26 0.40 0.36 0.42 0.35 0.48 0.43 11 0.53 0.40 0.63 0.50 0.58 0.29 0.46 0.55 0.48 0.24 12 0.47 0.63 0.55 0.44 0.10 0.47 0.50 0.50 0.72 0.29 3 1 0.42 0.31 0.17 0.40 0.45 0.32 0.39 0.33 0.45 0.39 2 0.00 0.26 0.34 0.34 0.42 0.32 0.27 0.17 0.37 0.28 3 0.04 0.02 0.41 0.34 0.37 0.23 0.40 0.41 0.36 0.45 4 0.08 0.51 0.58 0.00 0.60 0.41 0.51 0.18 0.56 0.54 5 0.55 0.18 0.04 0.54 0.44 0.35 0.54 0.51 0.36 0.26 6 0.50 0.32 0.52 0.37 0.32 0.37 0.39 0.44 0.50 0.37 7 0.16 0.52 0.59 0.41 0.42 0.34 0.31 0.23 0.48 0.61 8 0.32 0.44 0.55 0.58 0.51 0.49 0.09 0.29 0.36 0.44 9 0.52 0.69 0.47 0.62 0.53 0.00 0.19 0.42 0.38 0.43 10 0.53 0.04 0.63 0.42 0.36 0.48 0.50 0.31 0.50 0.67 11 0.03 0.04 0.41 0.54 0.50 0.44 0.42 0.56 0.42 0.56 12 0.46 0.04 0.53 0.47 0.56 0.28 0.36 0.13 0.39 0.32 4 1 0.41 0.40 0.54 0.37 0.35 0.40 0.39 0.29 0.28 0.47 2 0.00 0.15 0.01 0.36 0.31 0.28 0.21 0.41 0.20 0.38 3 0.02 0.01 0.33 0.39 0.48 0.34 0.23 0.35 0.36 0.35 4 0.61 0.48 0.42 0.45 0.49 0.03 0.38 0.37 0.43 0.38 5 0.30 0.57 0.10 0.51 0.34 0.50 0.62 0.54 0.34 0.39 6 0.00 0.52 0.02 0.31 0.19 0.30 0.44 0.60 0.35 0.60 7 0.54 0.19 0.09 0.11 0.46 0.22 0.48 0.21 0.19 0.30 8 0.56 0.43 0.48 0.59 0.47 0.01 0.41 0.47 0.04 0.58 9 0.34 0.04 0.46 0.60 0.15 0.39 0.43 0.24 0.43 0.52 10 0.00 0.62 0.47 0.09 0.36 0.34 0.46 0.20 0.41 0.32 11 0.37 0.63 0.32 0.54 0.45 0.24 0.02 0.14 0.44 0.33 12 0.20 0.71 0.56 0.09 0.60 0.15 0.67 0.26 0.09 0.35 187 Appendix D Performance of artificial neural networks The network performances obtained from an exhaustive search are presented here. Each table represents the MSE values of the validating set of the DANNs corresponding to each subseries for every collision type. The best performance among the networks with an equal number of hidden layers (only one and two) and neurons, but a different number of input lag orders were highlighted. 188 Table D.1 Mean squared error of NAR networks for the A2 subseries for TOTAL collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 1006 1018 1020 1023 1021 1027 1018 1022 1020 1037 2 304 361 322 314 326 329 323 329 337 342 3 278 289 246 256 244 257 261 239 273 249 4 279 302 276 251 235 256 241 250 249 239 5 225 221 231 225 245 256 244 192 215 239 6 194 192 187 150 139 186 204 175 174 187 7 148 163 131 109 144 152 170 158 156 175 8 120 108 113 124 112 116 126 112 113 189 9 111 112 119 97 103 133 132 161 117 147 10 66 60 67 70 68 60 75 80 81 91 11 55 53 56 56 62 54 64 58 67 89 12 57 54 59 54 51 50 64 68 84 91 2 1 1001 1516 1166 1049 1181 1047 1075 1133 1077 1125 2 2239 331 328 336 327 3586 361 358 344 334 3 469 580 567 289 279 242 258 247 242 254 4 285 280 253 257 277 222 249 262 255 240 5 239 234 190 399 246 223 264 281 326 281 6 216 220 209 186 238 220 222 228 295 236 7 143 350 127 181 151 195 201 217 246 273 8 121 139 159 157 129 139 173 172 146 195 9 238 382 140 146 190 185 161 269 236 208 10 866 96 70 75 186 195 103 233 134 159 11 156 426 54 70 98 178 147 153 112 217 12 157 55 60 431 103 118 150 315 171 159 3 1 1007 1320 1043 1200 1127 1026 1108 1088 1179 1198 2 893 1115 321 660 379 483 351 335 370 415 3 1669 986 537 582 687 248 268 279 302 286 4 1223 1512 689 743 246 444 243 513 248 304 5 1203 1213 1289 714 659 911 624 459 307 381 6 642 1064 214 710 730 260 331 346 306 498 7 819 452 485 495 394 393 435 445 811 541 8 2317 1619 289 144 265 607 503 396 307 400 9 1244 1563 1104 324 506 477 218 468 249 439 10 1533 430 151 339 760 366 502 713 438 346 11 1068 907 1044 296 255 683 572 463 360 304 12 1534 806 672 502 387 545 151 333 452 401 4 1 1017 1576 1525 1015 1386 1135 1103 1068 1158 1178 2 1026 1304 677 628 829 402 411 512 550 427 3 1896 1517 272 659 285 863 261 388 603 280 4 1872 1504 1199 1430 558 348 341 558 636 550 5 2429 1766 700 369 1317 513 330 744 504 699 6 2161 1523 475 1131 738 1315 587 387 466 1053 7 475 2147 510 232 642 1204 529 732 659 2705 8 2044 641 1369 1198 1053 193 702 379 942 549 9 2268 2282 855 193 250 408 792 1087 441 337 10 2109 969 1103 541 289 298 366 1304 768 537 11 2261 2093 553 1048 606 438 783 889 535 366 12 1729 1679 940 453 246 642 1182 680 597 914 189 Table D.2 Mean squared error of NAR networks for the D1 subseries for TOTAL collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 406 445 408 405 401 382 397 374 402 392 2 364 363 378 332 339 351 355 350 336 332 3 215 131 122 116 110 86 111 90 103 91 4 133 123 115 111 115 122 115 113 129 128 5 127 113 131 107 135 128 117 127 139 143 6 142 105 104 110 137 112 148 136 129 173 7 75 89 91 106 164 112 145 139 138 154 8 120 124 119 131 133 136 149 173 159 165 9 115 107 121 107 143 144 130 195 160 186 10 110 111 118 109 123 151 160 171 201 181 11 195 92 97 107 118 128 121 139 232 200 12 237 101 119 113 94 160 120 140 193 191 2 1 433 440 452 407 366 390 369 371 362 381 2 475 451 465 360 408 353 371 366 374 379 3 216 139 168 123 93 104 126 97 127 147 4 199 153 191 142 156 167 166 194 173 159 5 455 322 170 183 128 135 132 151 178 246 6 295 364 174 117 128 161 166 234 236 213 7 264 227 175 107 122 278 180 262 181 191 8 250 139 252 142 143 158 194 228 222 234 9 224 220 300 100 169 193 270 282 232 250 10 105 252 275 219 283 222 250 248 268 224 11 286 271 195 219 160 295 213 323 408 316 12 318 216 172 137 303 270 302 240 243 270 3 1 406 519 476 403 391 356 383 356 408 383 2 521 456 393 399 363 369 391 366 346 358 3 285 258 228 109 158 115 208 123 130 161 4 200 122 193 136 235 178 142 284 183 192 5 302 319 125 274 189 170 154 206 344 184 6 163 105 318 208 160 175 208 254 201 236 7 457 318 170 194 199 225 197 263 225 347 8 294 259 214 161 142 160 293 312 301 326 9 348 142 325 229 191 213 330 267 345 317 10 353 230 197 268 279 216 251 460 461 379 11 362 230 319 183 383 304 334 389 411 505 12 355 375 345 243 314 300 351 423 518 408 4 1 639 447 414 455 419 410 404 394 400 419 2 495 514 471 431 377 1099 397 383 401 365 3 461 461 239 144 118 119 208 210 182 177 4 387 286 203 146 233 278 179 182 290 235 5 576 256 281 149 131 174 276 207 276 245 6 323 284 427 326 251 239 351 300 357 315 7 213 259 131 310 168 211 288 318 373 458 8 303 288 453 190 244 287 453 319 479 489 9 230 339 373 213 385 417 389 485 521 626 10 486 314 237 396 422 416 389 405 412 568 11 373 218 317 320 271 369 525 502 364 612 12 342 178 560 218 372 474 482 524 505 367 190 Table D.3 Mean squared error of NAR networks for the D2 subseries for TOTAL collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 808 809 789 769 768 763 718 677 713 570 2 191 182 153 141 120 114 107 114 116 113 3 191 176 163 143 126 141 136 111 110 107 4 143 124 120 114 117 104 91 105 95 85 5 135 126 111 107 109 94 97 79 79 77 6 106 90 72 76 73 63 64 57 51 56 7 93 91 82 83 90 87 79 93 145 88 8 60 57 51 55 52 63 75 55 89 58 9 42 42 36 44 42 46 48 39 42 43 10 79 32 35 32 34 38 41 44 40 60 11 32 37 37 36 44 42 56 48 50 55 12 31 34 31 36 32 41 47 42 49 47 2 1 840 816 805 754 702 681 718 676 657 510 2 485 203 172 128 114 112 116 118 106 106 3 470 177 264 141 113 98 101 86 96 89 4 221 157 119 101 104 84 87 94 93 89 5 214 155 116 103 77 78 97 72 79 78 6 253 192 77 76 77 129 86 61 78 59 7 314 210 76 144 83 161 111 74 130 92 8 396 160 167 69 133 79 84 68 76 98 9 408 113 155 168 204 48 57 62 63 56 10 382 182 51 93 211 87 135 51 86 112 11 287 195 136 81 63 115 222 145 119 128 12 451 333 69 65 80 60 105 186 163 142 3 1 832 795 799 694 749 673 631 574 610 635 2 608 501 226 182 110 107 110 125 130 133 3 486 168 212 284 115 107 110 115 87 83 4 379 282 113 176 176 92 93 127 99 132 5 214 304 179 133 169 159 85 105 83 86 6 343 253 180 97 162 56 120 107 90 82 7 590 193 98 86 134 190 115 77 68 135 8 409 215 213 178 68 229 168 94 90 167 9 552 169 209 55 132 82 93 82 98 138 10 221 184 102 91 126 203 114 242 101 626 11 522 207 212 54 178 91 188 166 149 140 12 457 277 199 50 155 205 177 332 224 173 4 1 838 834 806 759 694 699 568 634 693 535 2 538 397 164 135 199 194 135 131 117 117 3 486 316 403 140 177 125 107 114 122 123 4 482 399 337 202 128 233 88 117 100 99 5 455 110 90 341 155 210 84 120 72 94 6 519 240 140 337 91 197 134 88 137 156 7 139 481 299 264 230 166 92 253 98 119 8 566 189 203 207 161 198 158 101 156 220 9 355 355 275 197 605 202 271 192 206 221 10 102 240 293 115 218 186 137 130 227 235 11 300 343 242 193 204 222 284 308 220 325 12 669 223 301 152 263 443 198 1575 235 199 191 Table D.4 Mean squared error of NAR networks for the A2 subseries for PDO collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 560 553 559 565 569 561 564 567 568 572 2 198 209 209 200 222 235 244 234 220 228 3 167 158 154 163 158 139 148 150 162 161 4 179 183 150 151 128 174 166 154 169 162 5 144 154 126 123 126 164 150 175 161 165 6 126 143 123 133 116 116 121 149 167 116 7 94 87 74 97 98 88 107 116 140 110 8 81 82 81 93 77 96 85 105 104 103 9 77 76 74 84 86 79 97 87 97 113 10 65 61 59 59 67 64 78 78 84 76 11 52 50 53 52 52 56 57 66 71 82 12 54 52 55 48 54 62 74 63 80 79 2 1 565 820 568 738 575 563 585 572 593 565 2 202 209 203 225 224 232 237 252 269 259 3 178 155 182 163 143 166 159 163 170 174 4 175 191 159 157 160 167 178 150 177 173 5 147 156 164 155 215 169 193 160 190 169 6 378 125 125 154 137 157 166 161 204 161 7 122 115 106 122 97 151 135 148 125 184 8 83 82 87 139 109 109 169 135 144 163 9 80 80 88 99 124 143 138 189 219 208 10 61 62 77 74 68 93 100 180 134 144 11 56 58 58 58 74 160 90 103 129 165 12 54 54 67 77 76 145 114 202 155 186 3 1 711 927 754 636 572 578 576 586 684 581 2 878 1122 698 223 557 229 282 245 293 339 3 1008 617 366 970 190 172 226 167 174 165 4 646 497 586 498 228 189 199 175 203 179 5 270 813 426 986 136 339 158 262 401 246 6 262 739 480 202 154 396 171 268 412 219 7 779 819 395 269 113 2307 214 202 369 429 8 1087 917 342 289 224 153 654 374 176 361 9 1242 399 230 569 303 320 196 494 358 410 10 691 1063 454 662 220 516 128 489 488 364 11 1012 609 125 456 145 187 198 342 284 462 12 753 517 462 101 223 333 328 697 365 278 4 1 899 1013 746 648 591 567 577 580 586 585 2 870 531 1056 309 240 707 373 238 261 271 3 1161 670 588 250 168 356 177 180 295 171 4 846 606 455 207 308 183 184 237 233 182 5 151 689 718 614 164 206 255 248 218 226 6 749 938 145 484 292 406 213 186 210 292 7 1133 467 638 440 251 182 158 1451 379 284 8 495 637 816 430 152 129 328 223 220 400 9 1209 652 427 261 199 456 480 218 522 361 10 769 690 539 343 182 409 295 177 202 184 11 646 568 532 409 185 276 290 292 376 313 12 57 961 1177 258 166 457 138 447 205 315 192 Table D.5 Mean squared error of NAR networks for the D1 subseries for PDO collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 358 368 353 366 369 371 364 391 392 383 2 375 348 339 327 345 325 355 339 369 350 3 187 159 141 145 139 140 137 132 138 163 4 117 106 102 99 93 90 84 90 84 88 5 94 89 90 92 89 92 88 95 92 89 6 90 91 96 103 109 111 112 109 119 133 7 76 98 88 99 90 106 107 106 126 136 8 115 103 108 107 107 113 124 110 141 131 9 99 90 92 98 109 111 110 120 112 125 10 136 101 98 99 132 134 118 108 195 129 11 82 136 99 105 94 96 102 126 128 126 12 82 97 99 103 94 132 117 100 139 104 2 1 476 387 364 381 388 392 387 410 399 421 2 423 407 354 358 373 387 385 386 380 396 3 416 225 177 152 139 125 151 142 147 162 4 255 146 192 120 98 89 99 111 96 130 5 283 168 91 120 145 99 122 106 150 105 6 281 150 129 225 113 116 159 140 125 169 7 235 146 156 187 127 120 159 155 197 145 8 390 211 122 101 124 164 146 206 148 201 9 283 168 96 137 175 171 182 269 221 235 10 273 195 121 246 224 187 206 335 204 269 11 368 160 100 133 231 148 241 185 213 293 12 288 172 196 147 196 217 222 243 223 253 3 1 371 368 456 416 379 411 394 399 407 480 2 492 402 363 338 377 380 390 398 383 424 3 302 368 153 187 184 173 223 161 171 168 4 285 172 265 123 143 183 111 103 139 149 5 274 171 182 91 98 121 123 210 132 192 6 353 217 99 191 153 131 203 115 154 200 7 319 258 238 292 213 155 207 206 195 377 8 323 297 161 124 249 178 155 196 242 453 9 222 188 177 188 256 267 276 233 302 324 10 553 218 140 228 259 220 182 288 312 369 11 431 187 225 200 185 373 368 416 459 322 12 225 95 290 307 145 238 272 315 318 328 4 1 435 421 402 407 377 390 406 426 418 414 2 516 451 459 472 406 402 361 394 445 436 3 320 259 322 152 190 233 192 272 208 129 4 266 215 261 136 257 169 193 103 159 189 5 379 161 208 236 197 146 146 172 223 200 6 402 359 140 173 167 278 175 281 250 297 7 345 300 251 253 173 224 184 330 274 222 8 299 385 175 143 168 156 291 196 283 344 9 402 365 261 228 254 269 363 356 448 317 10 398 270 264 263 273 332 361 529 370 295 11 376 162 325 252 252 403 371 492 510 426 12 338 130 347 289 356 363 349 277 383 423 193 Table D.6 Mean squared error of NAR networks for the D2 subseries for PDO collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 742 742 741 742 744 759 753 748 777 810 2 196 173 158 145 122 113 105 95 97 94 3 187 163 129 127 78 76 51 48 50 45 4 129 106 90 72 60 60 56 54 45 55 5 123 96 83 69 58 58 51 54 50 42 6 131 80 61 55 52 37 41 36 35 37 7 95 69 59 56 50 52 51 56 46 46 8 54 49 45 36 37 35 35 30 33 33 9 32 28 25 25 24 26 25 39 28 24 10 74 28 64 21 20 21 20 57 16 25 11 82 56 23 23 19 24 21 25 19 26 12 31 28 23 48 20 33 70 22 34 47 2 1 766 769 745 752 759 771 787 778 800 837 2 281 305 236 110 180 95 90 99 101 99 3 345 333 90 124 61 52 58 67 55 73 4 196 148 90 63 62 56 51 70 49 74 5 234 85 123 93 43 42 42 75 41 93 6 270 278 112 58 53 109 43 74 39 99 7 334 165 122 78 93 107 46 95 42 99 8 211 209 96 147 69 64 69 73 40 82 9 317 171 67 123 81 63 71 70 71 75 10 356 107 57 40 36 107 209 63 107 72 11 317 79 67 121 48 87 32 95 67 75 12 288 147 146 96 45 55 66 108 80 74 3 1 769 774 752 749 752 773 799 823 801 817 2 446 310 241 104 108 96 102 94 113 110 3 422 132 111 147 85 81 394 81 75 86 4 319 107 86 83 86 84 88 84 85 134 5 197 389 184 153 122 85 83 87 105 110 6 256 257 137 66 86 89 86 88 99 80 7 501 243 145 62 98 126 92 72 94 69 8 177 185 80 81 85 193 113 74 82 81 9 412 75 166 106 135 80 70 138 102 71 10 210 228 107 125 75 149 124 112 148 121 11 339 109 334 86 86 148 96 86 82 100 12 254 293 97 82 117 166 109 135 190 116 4 1 773 782 770 748 781 778 772 785 808 838 2 432 274 228 108 92 100 93 116 115 103 3 474 274 390 145 210 75 76 71 81 52 4 277 610 222 314 56 115 56 72 103 102 5 354 281 224 85 68 73 83 137 71 232 6 271 261 291 164 61 80 156 54 87 89 7 250 151 224 52 67 100 79 103 72 128 8 296 199 291 68 42 94 167 160 154 122 9 299 112 225 115 55 94 78 107 184 147 10 446 83 212 89 97 137 146 120 88 176 11 648 264 147 75 176 101 116 126 272 111 12 271 276 322 48 124 176 125 153 194 122 194 Table D.7 Mean squared error of NAR networks for the A2 subseries for KSI collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 87 87 87 87 85 88 82 87 86 84 2 17 15 14 13 12 12 11 11 10 10 3 16 14 11 11 9 10 7 8 7 7 4 16 13 12 8 8 5 7 5 4 4 5 14 11 9 6 6 6 4 4 4 4 6 13 8 7 6 4 3 4 4 3 4 7 10 6 6 5 5 3 3 3 3 3 8 8 4 2 2 3 2 1 2 2 2 9 6 4 3 3 2 1 1 2 2 2 10 2 2 1 1 2 1 1 1 2 2 11 1 1 1 1 1 1 1 2 1 2 12 1 1 1 1 1 1 1 2 2 2 2 1 88 103 86 85 87 84 86 88 88 87 2 15 14 13 12 10 11 9 9 9 9 3 16 84 15 9 12 7 6 7 7 5 4 115 13 10 9 6 6 8 6 5 10 5 13 10 9 6 6 6 5 5 4 6 6 14 58 7 5 5 3 3 3 4 5 7 14 15 6 8 5 3 15 6 5 4 8 94 6 5 3 5 2 5 5 2 4 9 134 4 24 4 5 3 6 7 5 4 10 2 2 1 2 5 3 7 3 2 3 11 1 1 1 3 2 3 3 4 3 11 12 1 3 2 2 3 3 2 5 15 4 3 1 133 147 109 111 84 84 92 84 92 90 2 98 98 69 10 11 26 12 8 8 8 3 40 64 73 10 31 8 6 7 7 6 4 92 57 41 8 13 12 12 8 10 5 5 76 34 8 9 6 9 10 20 7 9 6 123 79 7 20 7 18 21 13 10 8 7 68 67 38 22 26 18 14 11 11 10 8 132 51 35 19 12 8 7 6 18 15 9 87 30 39 39 4 34 5 8 10 13 10 98 53 25 5 8 57 15 19 7 9 11 86 18 78 27 13 16 14 14 12 13 12 70 50 27 12 16 12 18 5 8 25 4 1 132 175 133 89 87 86 85 83 87 88 2 52 128 42 87 49 9 11 10 10 11 3 152 76 36 46 15 10 8 14 7 6 4 51 97 82 10 7 24 6 22 8 26 5 127 130 47 16 56 14 21 18 11 17 6 47 8 40 35 300 29 16 12 32 22 7 78 74 89 16 11 21 38 13 28 28 8 68 62 74 20 6 12 24 37 20 28 9 115 103 36 43 19 3 17 24 30 15 10 62 40 66 57 40 33 19 30 27 22 11 134 125 39 53 59 15 41 18 20 20 12 145 132 104 29 3 12 23 9 18 23 195 Table D.8 Mean squared error of NAR networks for the D1 subseries for KSI collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 27 26 29 26 31 31 33 34 39 34 2 20 20 18 19 19 20 19 21 21 22 3 19 18 18 18 18 17 16 17 20 18 4 13 13 12 12 12 12 12 12 13 12 5 15 13 12 12 13 12 12 12 14 13 6 9 10 10 10 10 11 11 11 11 11 7 13 10 9 11 11 12 10 13 16 13 8 9 10 10 9 10 9 11 12 11 13 9 9 6 6 6 6 8 7 8 7 9 10 7 10 6 8 8 7 8 7 9 9 11 11 7 7 7 7 8 11 13 12 8 12 9 7 7 16 7 8 10 11 15 12 2 1 36 28 29 28 33 33 35 36 36 40 2 24 22 20 20 22 21 23 23 25 24 3 23 25 17 18 17 16 19 19 19 19 4 30 21 13 15 14 12 12 12 15 13 5 24 16 13 13 13 15 12 15 13 18 6 22 16 16 20 11 12 16 16 17 16 7 28 10 13 12 19 15 22 17 18 15 8 18 11 11 14 17 16 21 18 20 20 9 32 14 9 9 11 8 9 14 13 9 10 13 13 7 11 12 12 15 12 12 13 11 18 11 8 10 11 12 11 18 26 16 12 20 15 11 12 16 16 18 14 20 20 3 1 33 35 31 31 31 33 33 34 36 38 2 23 28 23 24 22 24 22 24 29 26 3 29 37 24 26 18 21 20 22 27 25 4 15 19 19 20 15 22 13 19 18 17 5 27 22 19 14 16 20 19 20 21 14 6 23 18 12 15 16 20 18 25 16 21 7 35 25 15 15 25 16 18 21 26 20 8 26 24 16 14 19 15 25 26 30 23 9 22 14 11 15 13 14 15 16 18 16 10 27 15 15 21 11 17 20 19 15 22 11 16 11 17 12 11 14 16 19 24 11 12 27 12 25 24 19 23 24 29 32 27 4 1 33 32 32 32 32 33 34 36 36 38 2 30 26 23 23 23 25 26 26 25 27 3 26 28 42 26 21 19 18 23 25 24 4 39 35 26 20 17 18 19 18 26 19 5 26 27 26 18 24 19 23 16 14 20 6 24 20 26 22 18 18 25 19 22 25 7 29 27 18 14 32 20 26 26 29 34 8 46 14 22 14 19 27 37 27 40 28 9 25 17 10 18 18 16 19 17 12 26 10 27 13 24 14 13 21 14 36 26 27 11 17 28 19 18 17 28 21 28 28 32 12 17 25 27 20 21 24 35 25 30 35 196 Table D.9 Mean squared error of NAR networks for the D2 subseries for KSI collisions. The best performance among the networks were highlighted. Number of hidden layers Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1 71 71 71 72 72 74 79 82 79 75 2 23 22 21 19 19 17 18 19 18 16 3 23 20 20 17 19 15 15 14 15 14 4 15 15 11 11 11 10 11 9 9 9 5 14 13 12 11 10 9 11 10 11 10 6 10 9 8 8 9 9 8 8 8 9 7 12 9 8 8 8 8 8 8 8 10 8 8 8 7 7 8 7 7 8 8 9 9 13 7 7 7 6 7 7 7 9 9 10 6 5 5 5 5 5 5 6 6 7 11 9 5 5 5 5 6 5 5 12 7 12 5 5 4 4 5 5 8 9 11 7 2 1 73 71 74 73 71 74 80 88 92 78 2 47 20 25 20 21 22 20 19 19 22 3 26 32 18 18 18 15 15 11 12 12 4 28 40 11 14 13 12 10 11 13 12 5 33 26 21 13 17 11 12 11 12 11 6 30 15 8 13 10 17 18 18 11 13 7 22 9 15 12 9 14 13 14 11 11 8 10 9 12 12 10 24 17 12 11 13 9 20 21 14 8 12 12 10 10 13 12 10 24 15 10 23 18 8 14 12 13 11 11 13 8 6 8 14 10 14 9 17 10 12 12 12 20 11 15 10 8 21 15 13 3 1 74 78 70 71 73 73 79 86 83 87 2 21 27 22 20 21 25 19 24 18 20 3 39 36 25 26 17 16 15 17 18 15 4 28 26 12 15 15 12 22 15 15 14 5 26 26 17 14 15 17 12 15 14 19 6 16 17 17 19 29 17 13 13 18 15 7 18 17 17 15 20 20 13 17 13 15 8 9 15 22 20 12 9 19 22 17 21 9 29 12 22 11 21 21 24 16 24 18 10 29 10 24 10 12 17 22 16 25 22 11 16 10 7 20 13 16 16 20 11 29 12 35 16 19 19 18 9 27 25 25 18 4 1 77 75 76 76 77 73 77 79 83 79 2 26 50 27 35 20 24 23 23 23 23 3 27 21 22 23 20 17 27 17 22 18 4 48 32 26 19 21 20 16 19 19 13 5 33 24 18 23 24 17 21 28 20 23 6 23 15 22 19 15 16 31 16 16 27 7 25 21 26 18 22 18 26 22 21 34 8 29 18 18 19 25 20 23 25 23 17 9 29 33 25 35 19 30 21 53 29 31 10 21 21 21 15 31 29 29 30 26 24 11 22 116 29 16 19 19 25 25 19 21 12 30 16 21 20 20 27 17 21 25 22 197 Table D.10 Mean squared error of NARX networks for the A2 subseries for TOTAL collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 895 928 940 1000 963 961 980 958 987 950 2 342 320 296 298 314 273 306 313 305 275 3 283 289 241 245 244 250 244 389 268 254 4 435 257 262 227 210 235 217 242 380 248 5 305 216 221 218 201 211 236 239 265 229 6 213 192 194 186 181 174 196 698 162 206 7 139 144 125 141 131 176 145 178 139 192 8 111 116 91 78 101 104 126 141 124 135 9 283 98 116 110 130 125 118 169 165 165 10 46 49 49 60 66 72 79 76 69 87 11 41 46 41 52 48 73 60 138 88 64 12 177 45 50 48 47 47 47 73 189 81 Table D.11 Mean squared error of NARX networks for the D1 subseries for TOTAL collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 483 431 471 432 459 488 490 474 498 508 2 347 367 352 348 358 384 365 352 372 350 3 152 160 135 134 96 119 129 112 105 139 4 154 118 110 116 131 120 136 140 176 143 5 166 113 110 112 127 123 142 148 136 154 6 91 95 130 108 127 136 167 210 172 197 7 135 95 113 108 125 127 156 158 150 156 8 132 136 151 142 165 166 164 192 168 199 9 122 120 123 143 150 165 177 163 232 262 10 117 113 147 129 136 155 156 222 294 192 11 149 94 103 126 126 158 180 183 253 248 12 124 127 103 113 132 218 153 230 243 184 198 Table D.12 Mean squared error of NARX networks for the D2 subseries for TOTAL collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 900 896 957 915 940 848 966 987 950 1160 2 235 267 225 207 173 187 260 193 204 219 3 259 191 201 199 207 176 211 177 187 195 4 153 151 152 154 152 150 172 164 175 179 5 148 133 134 124 122 162 139 151 118 186 6 166 117 110 112 116 99 127 112 142 116 7 106 102 92 103 101 89 102 127 143 130 8 77 75 68 136 77 85 69 89 89 100 9 78 75 67 73 78 74 94 102 84 102 10 35 37 36 39 43 44 46 51 47 58 11 87 38 40 46 38 48 45 76 51 68 12 30 38 38 40 40 41 52 55 44 57 Table D.13 Mean squared error of NARX networks for the A2 subseries for PDO collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 524 535 552 543 535 559 550 571 572 593 2 200 202 199 193 180 213 208 222 231 204 3 176 174 163 156 140 190 162 261 172 205 4 365 164 154 156 131 172 174 162 237 164 5 148 134 148 125 138 162 139 175 183 185 6 120 126 142 119 124 140 125 144 157 141 7 189 75 84 98 98 130 119 104 310 139 8 445 82 74 91 97 103 105 102 202 106 9 76 81 76 78 102 208 106 92 147 158 10 351 272 54 58 56 75 56 60 75 144 11 77 44 44 51 51 48 65 46 57 60 12 173 46 43 47 57 51 74 61 77 147 199 Table D.14 Mean squared error of NARX networks for the D1 subseries for PDO collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 408 371 394 399 382 392 387 410 388 400 2 351 330 323 344 337 343 330 324 329 349 3 179 174 152 152 139 159 144 151 141 168 4 106 106 93 103 100 106 117 105 127 135 5 125 83 84 85 100 106 106 114 108 133 6 89 102 129 108 122 98 117 129 678 145 7 74 83 87 90 100 115 119 123 113 144 8 101 99 108 111 115 123 128 135 121 188 9 105 190 118 99 165 143 175 124 157 180 10 101 106 106 111 112 110 124 126 151 196 11 79 167 84 155 108 115 154 140 132 132 12 116 76 83 84 155 92 110 131 173 174 Table D.15 Mean squared error of NARX networks for the D2 subseries for PDO collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 823 796 776 778 751 769 795 799 814 789 2 215 178 179 157 187 153 162 129 123 141 3 182 194 145 131 144 127 447 107 150 115 4 144 112 100 122 99 129 105 81 110 181 5 200 103 102 118 106 105 92 81 83 117 6 157 98 91 79 82 80 75 71 100 80 7 170 91 80 88 80 80 117 82 91 78 8 69 59 59 59 56 82 59 65 62 107 9 70 67 54 55 66 52 63 67 70 72 10 44 41 41 40 35 30 81 48 70 43 11 44 39 42 74 44 41 42 41 53 54 12 70 34 34 42 43 43 74 53 52 53 200 Table D.16 Mean squared error of NARX networks for the A2 subseries for KSI collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 74 68 62 57 53 50 51 53 54 52 2 15 13 12 12 11 10 10 11 10 10 3 27 14 10 10 9 7 9 8 9 9 4 14 12 11 9 8 8 7 6 7 6 5 13 10 7 6 6 6 5 5 4 4 6 29 7 4 5 4 3 3 3 4 4 7 26 6 5 3 4 3 4 3 3 4 8 7 4 4 2 2 2 2 3 3 3 9 5 3 3 2 2 2 4 2 3 5 10 1 1 1 1 1 1 2 1 1 3 11 1 1 1 1 1 1 1 1 1 1 12 20 1 1 1 1 1 2 2 2 1 Table D.17 Mean squared error of NARX networks for the D1 subseries for KSI collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 27 28 28 30 33 32 32 34 35 41 2 19 18 19 19 19 19 20 20 20 21 3 18 18 18 18 19 18 18 18 18 17 4 13 13 13 12 12 12 13 14 13 14 5 13 12 13 13 13 14 13 13 15 16 6 13 10 12 12 11 12 13 13 12 11 7 10 10 10 10 10 12 12 11 15 17 8 10 10 9 9 12 10 12 14 16 16 9 13 11 11 11 12 12 14 14 15 18 10 10 12 10 10 11 11 13 14 14 18 11 10 10 12 10 12 14 18 13 20 16 12 11 12 13 13 15 21 19 24 24 19 201 Table D.18 Mean squared error of NARX networks for the D2 subseries for KSI collisions. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 64 66 64 68 68 78 83 81 82 84 2 21 24 23 24 26 27 27 29 35 29 3 21 22 19 20 21 22 27 30 24 29 4 14 13 13 13 14 14 14 18 15 16 5 12 13 12 11 16 16 15 13 17 19 6 7 8 7 8 8 8 9 9 10 12 7 7 8 7 8 8 8 10 10 10 11 8 8 7 7 7 6 8 7 7 7 8 9 9 7 6 7 7 10 7 8 8 8 10 5 6 5 5 6 7 8 11 8 8 11 6 6 5 5 5 7 6 7 11 8 12 4 4 4 5 6 5 6 6 10 8 Table D.19 Mean squared error of NAR networks for the A2 subseries for TOTAL collision series before installing phase one in Springfield Road Corridor. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 0.289 0.491 0.351 0.634 0.412 0.461 1.643 1.976 0.821 1.612 2 0.172 0.201 0.230 0.193 0.218 0.257 0.283 0.838 0.504 0.356 3 0.125 0.117 0.139 0.258 0.240 0.331 0.433 0.350 0.490 0.379 4 0.101 0.155 0.112 0.398 0.540 0.187 0.345 0.760 0.641 0.303 5 0.161 0.101 0.119 0.286 0.312 0.239 0.274 0.710 0.342 0.350 6 0.029 0.057 0.105 0.066 0.359 0.185 0.639 0.225 0.260 0.475 7 0.019 0.057 0.075 0.107 0.181 0.297 0.246 0.550 0.394 0.355 8 0.018 0.033 0.082 0.147 0.137 0.322 0.499 0.391 0.385 0.372 9 0.117 0.268 0.081 0.149 0.194 0.342 0.311 0.506 0.487 0.414 10 0.114 0.129 0.094 0.330 0.176 0.483 0.231 0.619 0.330 0.478 11 0.142 0.099 0.185 0.162 0.221 0.252 0.230 0.566 0.451 0.787 12 0.077 0.025 0.072 0.275 0.264 0.267 0.210 0.647 0.479 0.739 202 Table D.20 Mean squared error of NAR networks for the D1 subseries for TOTAL collision series before installing phase one in Springfield Road Corridor. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1.48 2.89 3.85 4.63 5.41 6.21 4.03 5.05 7.91 5.67 2 2.26 2.72 2.28 2.71 3.29 3.80 4.56 6.89 3.53 6.70 3 1.39 2.91 2.35 2.25 2.34 3.21 3.34 4.65 4.72 4.51 4 1.37 1.89 1.91 2.63 2.13 3.02 2.99 6.82 3.41 3.79 5 1.11 1.67 2.17 2.17 4.61 3.12 3.62 3.88 5.39 6.15 6 1.47 5.56 2.22 2.64 2.75 3.38 8.34 4.29 4.22 6.17 7 1.59 1.92 2.79 3.05 3.76 4.24 3.26 4.96 4.70 5.12 8 1.48 2.12 2.36 2.60 3.27 3.32 3.80 3.83 5.71 3.42 9 1.82 1.87 2.71 2.27 2.28 3.20 3.79 5.28 3.47 4.99 10 1.72 2.22 2.66 3.03 3.27 3.85 3.77 4.83 4.51 5.45 11 1.62 2.34 2.57 2.25 4.04 3.37 5.53 3.97 3.82 8.55 12 1.58 2.16 3.31 2.39 2.43 5.62 3.02 4.69 4.06 6.61 Table D.21 Mean squared error of NAR networks for the D2 subseries for TOTAL collision series before installing phase one in Springfield Road Corridor. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1.06 1.02 1.06 1.07 1.12 1.08 1.21 1.23 1.33 1.28 2 0.41 0.38 0.34 0.39 0.36 0.44 0.59 0.36 0.37 0.38 3 0.36 0.28 0.29 0.22 0.27 0.28 0.28 0.22 0.38 0.32 4 0.17 0.15 0.14 0.15 0.17 0.14 0.16 0.16 0.24 0.23 5 0.17 0.16 0.16 0.16 0.17 0.17 0.21 0.19 0.23 0.20 6 0.11 0.11 0.11 0.12 0.14 0.15 0.17 0.17 0.21 0.17 7 0.12 0.11 0.12 0.12 0.14 0.16 0.25 0.22 0.24 0.33 8 0.16 0.18 0.15 0.13 0.17 0.16 0.21 0.19 0.25 0.28 9 0.18 0.16 0.13 0.14 0.18 0.26 0.31 0.39 0.34 0.40 10 0.16 0.13 0.13 0.15 0.24 0.25 0.29 0.35 0.35 0.43 11 0.16 0.13 0.11 0.27 0.27 0.31 0.47 0.61 0.49 0.53 12 0.10 0.17 0.18 0.23 0.29 0.45 0.48 0.61 0.52 0.90 203 Table D.22 Mean squared error of NAR networks for the A2 subseries for TOTAL collision series before installing phase two in Springfield Road Corridor. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 0.048 0.308 0.049 0.058 0.060 0.051 0.056 0.076 0.069 0.063 2 0.024 0.019 0.020 0.020 0.018 0.020 0.020 0.020 0.021 0.022 3 0.021 0.021 0.021 0.021 0.019 0.021 0.025 0.031 0.021 0.024 4 0.031 0.033 0.027 0.030 0.038 0.041 0.043 0.044 0.050 0.046 5 0.033 0.040 0.047 0.046 0.054 0.062 0.066 0.060 0.076 0.074 6 0.033 0.031 0.033 0.032 0.039 0.039 0.041 0.044 0.046 0.042 7 0.027 0.026 0.026 0.026 0.030 0.031 0.029 0.029 0.044 0.039 8 0.027 0.024 0.027 0.026 0.030 0.030 0.027 0.034 0.051 0.042 9 0.059 0.026 0.029 0.033 0.025 0.030 0.036 0.031 0.044 0.052 10 0.022 0.022 0.021 0.020 0.025 0.026 0.031 0.037 0.036 0.040 11 0.020 0.022 0.025 0.024 0.030 0.033 0.034 0.043 0.039 0.043 12 0.024 0.024 0.029 0.025 0.022 0.042 0.039 0.085 0.064 0.062 Table D.23 Mean squared error of NAR networks for the D1 subseries for TOTAL collision series before installing phase two in Springfield Road Corridor. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 0.51 0.51 0.48 0.55 0.51 0.49 0.50 0.50 0.48 0.49 2 0.30 0.31 0.31 0.32 0.33 0.34 0.35 0.36 0.35 0.34 3 0.30 0.32 0.29 0.30 0.31 0.35 0.35 0.32 0.35 0.41 4 0.25 0.27 0.25 0.26 0.35 0.37 0.38 0.38 0.36 0.41 5 0.19 0.20 0.21 0.20 0.26 0.24 0.26 0.34 0.35 0.33 6 0.27 0.35 0.30 0.44 0.45 0.46 0.46 0.50 0.49 0.55 7 0.26 0.34 0.33 0.39 0.37 0.41 0.44 0.51 0.60 0.70 8 0.26 0.30 0.29 0.40 0.44 0.44 0.50 0.58 0.75 0.73 9 0.37 0.27 0.37 0.35 0.45 0.43 0.62 0.57 0.66 0.56 10 0.35 0.34 0.33 0.33 0.46 0.37 0.53 0.46 0.59 0.63 11 0.21 0.26 0.26 0.33 0.35 0.44 0.47 0.53 0.52 0.68 12 0.22 0.23 0.36 0.36 0.45 0.40 0.52 0.53 0.55 0.59 204 Table D.24 Mean squared error of NAR networks for the D2 subseries for TOTAL collision series before installing phase two in Springfield Road Corridor. The best performance among the networks were highlighted. Number of lags Number of hidden neurons 1 2 3 4 5 6 7 8 9 10 1 1.721 1.646 1.616 1.587 1.502 1.595 1.559 1.651 1.500 1.621 2 0.545 0.590 0.634 0.596 0.613 0.545 0.602 0.582 0.708 0.696 3 0.478 0.486 0.428 0.408 0.448 0.536 0.342 0.413 0.292 0.318 4 0.281 0.266 0.225 0.243 0.264 0.295 0.241 0.275 0.243 0.310 5 0.290 0.246 0.294 0.313 0.212 0.256 0.247 0.274 0.320 0.247 6 0.200 0.194 0.177 0.183 0.161 0.164 0.182 0.223 0.292 0.236 7 0.203 0.352 0.165 0.185 0.179 0.212 0.210 0.215 0.243 0.280 8 0.154 0.154 0.132 0.152 0.159 0.163 0.227 0.185 0.235 0.221 9 0.159 0.162 0.150 0.153 0.212 0.182 0.177 0.230 0.246 0.294 10 0.139 0.162 0.161 0.172 0.161 0.188 0.195 0.428 0.241 0.298 11 0.133 0.167 0.152 0.166 0.207 0.196 0.175 0.189 0.266 0.305 12 0.114 0.175 0.149 0.167 0.198 0.202 0.223 0.256 0.258 0.348
Thesis/Dissertation
eng
Civil Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Application of dynamic artificial neural network to road collision prediction modelling
Text