Application of the Cox Regression Model for Analysis of Railway Safety Performance
Preprint – Released for Publication p. 1 of 12
Application of the Cox Regression Model for analysis of Railway Safety Performance
Hendrik Schäbe (TÜV Rheinland) & Jens Braband (Siemens Mobility)
Abstract
The assessment of in-service safety performance is an important task, not only in railways. For example it is important to identify deviations early, in particular possible deterioration of safety performance, so that corrective actions can be applied early. On the other hand the assessment should be fair and objective and rely on sound and proven statistical methods. A popular means for this task is trend analysis. This paper defines a model for trend analysis and compares different approaches, e. g. classical and Bayes approaches, on real data. The examples show that in particular for small sample sizes, e. g. when railway operators shall be assessed, the Bayesian prior may influence the results significantly.
Introduction
The analysis of railway safety performance has been of interest for many years. On a nationwide level already so called Common Safety Targets and Common Safety Indicators have been proposed (EU (2009)). Especially, it is important to analyse trends and to derive judgement on whether the performance is improving, deteriorating or unchanged. This is performed today by the European Union Agency for Railways mainly based on weighted averages or moving averages (European Union Agency for Railways (2018)). Due to the normally sparse data, in many cases no judgement can be made, since the small sample size does not allow to provide a statistically significant result. In this paper we analyse numbers of severe railway accidents as reported by Evans (2020).Note that this topic is of particular relevance as the European Railway Agency is drafting a legal text on the assessment of safety levels and safety performance of railway operators in the EU. But in the current draft (European Union Agency for Railways (2020)) the term safety level is used for the statistical safety performance and is defined currently in a very narrow sense: “‘safety level’ means the weighted sum of occurrences of eligible events … corresponding to a given volume of operation … and normalised by this volume of operation…”. Currently also no statistical procedure for the evaluation of these safety levels has been proposed so that this paper tries to fill the gap. For this purpose already several proposals have been made e. g. Evans (2020) or Andrasik (2020), which are taken into account and compared.
The model
Railway accidents, in particular severe accidents with fatalities, are rare events and as such the number of severe accidents, say K, follow a Poisson distribution, see Braband and Schäbe (2014): P(K=k) = 𝜆 𝑘 𝑘! exp (−𝜆) . (1) If the accidents relate to n different time intervals, which are numbered by the index i, this easily yields P(Ki=ki) = 𝜆 𝑖𝑘𝑖 𝑘 𝑖 ! exp (−𝜆 𝑖 ) . (2) For different time intervals, we have here assumed different intensities i, since we may not assume a priori, that the intensity with which the accidents occur is always the same. Preprint – Released for Publication p. 2 of 12 Often the i are expected to form a decreasing sequence of values. In order to model this effect, we use the Cox regression model using the time as explanatory variable, see Cox and Oakes (1984). Then we have: i = ti). (3) Practically (3) means that we assume n equidistant time intervals, where ti is the value describing the center of the i-th time interval. Moreover, ∏ 𝜆 𝑖𝑘𝑖 𝑘 𝑖 ! exp (−𝜆 𝑖 ) 𝑛𝑖=1 . (4) This yields the maximum likelihood estimators, which are obtained by the following equations. 𝜆̂ = ∑ 𝑘 𝑖𝑛𝑖=1 ∑ exp (−𝛽̂𝑡 𝑖 ) 𝑛𝑖=1 (5) and ∑ 𝑘 𝑖 𝑡 𝑖𝑛𝑖=1 ∑ exp(−𝛽̂𝑡 𝑖 ) = ∑ 𝑘 𝑖 ∑ 𝑡 𝑖𝑛𝑖=1𝑛𝑖=1 exp (−𝛽̂𝑡 𝑖 ) 𝑛𝑖=1 (6) Note that the estimator for must be obtained by numerical solution of equation (6). In the next step is larger than 0 or not. If this is the case, one could conclude a decreasing behavior of the severe accidents. The fact that a maximum likelihood estimator is asymptotically normally distributed with variance 1/I can be used. Here I is the Fisher information, computed by I =- 𝜕 𝜕𝛽 ln(𝐿𝑖𝑘) = ∑ 𝜆 𝑡 𝑖2 exp (−𝛽𝑡 𝑖 ) 𝑛1 . (7) Then, a confidence interval for b is [ 𝛽̂ -z ; 𝛽̂ +z ]. (8) with = √𝐼 (9) and z is the 1- Quantile of the Normal distribution. Then, the confidence interval has a coverage of 1-2* . The result (8) can also be used for a statistical test for the Hypothesis: there is no decrease in the number of severe accidents. If UConf = 𝛽̂ - () is larger than zero, the hypothesis would would be declined and improvement, i.e. a significant decrease of severe accidents over time would be concluded. Preprint – Released for Publication p. 3 of 12 It has to be noted that the maximum likelihood estimators are asymptotically efficient, i.e. have the smallest spread among all estimators if the sample size is large enough.
Examples
Using the data of Evans (2020) we have analysed the data of several countries. The figures below show the counts and the fitted curve. Also the last two points of the counts represent the counts for 1980-2019 and 1990-2019 for a 5 years period – for comparison with the other data. We have used data from the entire European union plus Switzerland plus Norway and also the date from several large countries as Germany, France, Italy and the United Kingdom. As an example for a small country we have used Estonia. Here it becomes evident, that small sample sizes are a real problem. Figure 1 Number of accidents with fatalities (counts and model) of the 28 EU countries + Switzerland + Norway Since UConf takes the value of 0,043 with first kind error of 10% we conclude a significant decrease. EU28+NO+CH
Exp Theoret
Preprint – Released for Publication p. 4 of 12 Figure 2 Number of accidents with fatalities (counts and model) of Germany Since UConf takes the value of 0,046 with first kind error of 10% we conclude a significant decrease. Figure 3 Number of accidents with fatalities (counts and model) of France Since UConf takes the value of 0,0517 with first kind error of 10% we conclude a significant decrease. DE Exp Theoret FR Exp Theoret
Preprint – Released for Publication p. 5 of 12
Figure 4 Number of accidents with fatalities (counts and model) of the United Kingdom Since UConf takes the value of 0,0757 with first kind error of 10% we conclude a significant decrease.
Figure 5 Number of accidents with fatalities (counts and model) of Italy UK Exp Theoret IT Exp Theoret
Preprint – Released for Publication p. 6 of 12 Since UConf takes the value of -0,00403 with first kind error of 10% we cannot conclude a significant decrease, although visually a decreasing tendency can be seen. This, however, is not significant and can be caused by random influences.
Figure 6 Number of accidents with fatalities (counts and model) of Estonia Since UConf takes the value of -0,38 with first kind error of 10% we cannot conclude a significant decrease. This, however is a result of the small number of accidents (only two). A general decreasing trend can be seen. EE Exp Theoret
Preprint – Released for Publication p. 7 of 12 Figure 7 Number of accidents with fatalities (counts and model) of Romania Since OConf takes the value of -0,00094 with first kind error of 10% we cannot conclude a significant increase, which also visually can be seen The results how, that the model allows a statistical analysis and a significance test if enough data are available. We have analysed the data collected from 1980 until 2019, which is a long time. This allowed us to derive results in many cases, but not in all.
Bayesian approach
Frequently, expert or prior information shall be used. In order to avoid unnecessary complexity, we use a conjugate prior ( ) in the same form as the likelihood function (4) with (3), i.,e. ln ( ( )) ~ ∑ [𝑎 𝑖𝑛𝑖=1 (ln(𝜆 ) − 𝛽𝜏 𝑖 ) − 𝜆 exp (−𝛽𝜏 𝑖 )]. (11) This prior has the effect, as if in the time intervals with central points i ai severe accidents would have been occurred. If all ai would be 0, we would have a non-informative prior. Due to the use of the conjugate prior, for ti = i a Bayes estimator for the maximum of the posterior density is equivalent to the maximum likelihood estimator with data qni+(1-q)ai using (5) and (6). Here qi is the weight that we give to the prior information. An approximate HPD interval (high probability interval) of the posterior density can be obtained from (8) using data ni+ai. Here, we have approximated the posterior by a normal distribution around it maximal value. One has to note that within this model the information from the prior distribution and from the sample are combined and RO Exp Theoret
Preprint – Released for Publication p. 8 of 12 that leads to a seemingly higher sample size and therefore, to larger values of
0. With one exception- if a non-informative prior is used. If, one is only interested in the parameter , for deriving a judgement on an increase or decrease of the number of fatalities on this would not play a role. Statistics with different time interval lengths
The model can easily be adopted for the case of time intervals with different length. Then, the formulae (2) – (6) would be rewritten in the following manner The Poisson distribution is kept unchanged: P(Ki=ki) = 𝜆 𝑖𝑘𝑖 𝑘 𝑖 ! exp (−𝜆 𝑖 ) . (12) The different intensities i now also depend on the lengths of the time intervals, say Ti. Then (3) is changed into: i = ti). (13) That means, the intensity is proportional to the length of the interval and to the well-known regression factor and the
0, which is now a parameter per time unit, which is different from the model in (1) – (4). The time ti lies in the center of the i-th interval having length Ti. Now, the log likelihood function is l = ln ( lik) = ∑ 𝑘 𝑖 ( 𝑛𝑖=1 ln 𝜆 − 𝛽𝑡 𝑖 + ln(𝑇 𝑖 )) − 𝑙𝑛𝑘 𝑖 ! − 𝜆 𝑇 𝑖 exp (−𝛽𝑡 𝑖 ) (14) This yields again the maximum likelihood estimators, which are obtained by the following equations. 𝜆̂ = ∑ 𝑘 𝑖𝑛𝑖=1 ∑ 𝑇 𝑖 exp (−𝛽̂𝑡 𝑖 ) 𝑛𝑖=1 (15) and ∑ 𝑘 𝑖 𝑡 𝑖𝑛𝑖=1 ∑ 𝑇 𝑖 exp(−𝛽̂𝑡 𝑖 ) = ∑ 𝑘 𝑖 ∑ 𝑡 𝑖𝑛𝑖=1𝑛𝑖=1 𝑇 𝑖 exp (−𝛽̂𝑡 𝑖 ) 𝑛𝑖=1 (16) Formula (7) becomes I =- 𝜕 𝜕𝛽 ln(𝐿𝑖𝑘) = ∑ 𝜆 𝑇 𝑖 𝑡 𝑖2 exp (−𝛽𝑡 𝑖 ) 𝑛1 . (17) For simplicity assume now only two time intervals with T1 = T and T2 = (1- ) T so that T1/T2= (−) and T= T1+T2 and k = k1+k2. Then, t1 = T/2 and t2 = (1+ )T/2 𝜆̂ = 𝑘α Texp(− 𝛽̂αT2 )+(1−α)Texp (− 𝛽̂(1+𝛼)𝑇2 ) (18) Preprint – Released for Publication p. 9 of 12 [k1 T+k2(1+ )T] [ α Texp (− 𝛽̂αT2 ) + (1 − α) Texp (− 𝛽̂(1+𝛼)𝑇2 )] = k [ α 𝑇 exp (− 𝛽̂αT2 ) + (1 − 𝛼 ) 𝑇 exp (− 𝛽̂(1+𝛼)𝑇2 )] (19) I = T/2) + (1+ )(1-
2) exp(- (1+ )T/2)}/4 (20) Using a Bayesian prior distribution with the same and T but parameters a1 and a2, (18), (19) and (20) keep valid, only with k1 replaced by k1 + a1 and k2 replaced by k2+a2. Again, one need to note that for a statistical inference on the parameter l0 the effect of the Bayes estimator here is that of an increased number of severe accidents. Two sample comparison – an example In this section we will provide an example for a comparison of the safety behavior of national railway systems in 2010-2014 compared with 2005-2009behavior three cases: Greece, Estonia, European Union. We will use the data from Evans (2020). We will use the Bayesian methods given in the last section. Two different priors will be studied. The first one is a non-informative one with a1=a2=0. In the second case, we use a prior proposed by Andrasik. Andrasik used a Binomial distribution to model the fatalities in the both time intervals. For intervals of equal length he arrives at parameters a1 =a2=2 in our notation. Formulae (18) – (20) give in our special case with = ½ and T = 10 years 𝜆̂ = 𝛽̂T4 )+exp (− )) (21) (k1+3 k2) [ exp (− 𝛽̂T4 ) + exp (− )] = k [ exp (− 𝛽̂T4 ) + 3 exp (− )] (22) This equation can be simplified to give 𝛽̂ = − 𝑙𝑛 ( 𝑧−13−𝑧 ) with z = (k1+3k2)/k = 1 + 2k1/k (23) The information becomes I = T/4) + 9 exp(-3 T/4)} / 32 (24) As a result we get – Released for Publication p. 10 of 12 . The result is given in the following table. Classical Test Bayes Test 2005-2009 2010-2014 Infor-mation I Beta (esti-mator) lower confidence level upper confidence level Infor-mation I Beta (esti-mator) lower confidence level upper confidence level EU28+NO+CH 35 34 431,25 5,798E-03 -7,341E-02 8,500E-02 456,25 5,480E-03 -7,153E-02 8,249E-02 Greece (EL) 2 1 18,75 1,386E-01 -2,412E-01 5,185E-01 43,75 5,754E-02 -1,911E-01 3,062E-01 France (FR) 1 2 18,75 -1,386E-01 -5,185E-01 2,412E-01 43,75 -5,754E-02 -3,062E-01 1,911E-01
Table 1 Result of Classical and Bayesian two sample tests
The Bayesian confidence intervals are HPD (High Probability Distribution) intervals. The intervals have been computed with coverage of 90%. It can be seen that with the present data the classical and the Bayesian point estimator of show the same tendency regarding improvement (positive b) or deterioration (negative ). In none of the three cases the change is significant. Furthermore it can be seen that the Fisher information of the estimator for is larger for the Bayes estimator than for the classical one. This is especially evident for small sample sizes. That means that a Bayes estimator can influence the result and this influence is stronger the more information is introduced into the posteriori distribution via the prior distribution. In the following table we have provided an example with the country “Ex1”. Here, the classical statistics would see a significant change (deterioration), whereas the Bayes statistics would see no significant change. A general tendency can be summarised: In Bayesian statistics the prior distribution can influence the result, the more, the higher the information contained in the prior is compared with the sample. In these cases, an inappropriately chosen prior might corrupt the result. This is a general result and does not hold true only for a specific statistical model, see e.g. Zellner (1971) chapter 2.11 or Schäbe (1993). Preprint – Released for Publication p. 11 of 12
Classical Test Bayes Test 2005-2009 2010-2014 Infor-mation I Beta (esti-mator) lower confidence level upper confidence level Infor-mation I Beta (esti-mator) lower confidence level upper confidence level Ex1 84 100 1150 -3,487E-02 -8,337E-02 1,363E-02 1175 -3,413E-02 -8,211E-02 1,386E-02
Table 2 Result with a hypothetical country
Conclusion
In this paper we have shown how railway safety behavior can be modeled with the help of a Poisson distribution connected with the Cox model. We have set up the model for different situations, with equally spaced time intervals, time intervals of different length and developed estimators and tests. This has been done for classical statistics as well as for Bayesian statistics. We have also considered, how a Bayesian prior can influence the result, especially with small sample sizes. This becomes even more important if not only the safety performance of countries shall be assessed but that of single operators (European Union Agency for Railways (2020)). Nevertheless, statistics on railway safety behavior should be based on characteristics giving larger sample size as incident statistics rather than on accident statistics, see Braband and Schäbe (2013a).
References
Andrasik, R., (2020) Evaluation of safety targets based on Bayesian inference Braband, J. Schäbe, H., (2013a) Assessment of national reference values for railway safety: a statistical treatment, Journal of Risk and Reliability, 227(4) 405 – – Steenbergen et al. (Eds) 2014 Taylor & Francis Group, London, ISBN 978-1-138-00123-7), pp. 1297-1302. Cox D.R., Oakes, D. (1984), Analysis of Survival Date, London, New York, Chapman and Hall European Union Agency for Railways: Report on Railway Safety and Interoperability in the EU, 2018 European Union Agency for Railways: COMMISSION DELEGATED REGULATION (EU) establishing common safety methods for assessing the safety level and the safety performance of railway operators at national and Union level, Draft, July 2020 EU: COMMISSION DECISION of 5 June 2009 on the adoption of a common safety method for assessment of achievement of safety targets, as referred to in Article 6 of Directive 2004/49/EC of the European Parliament and of the Council (2009/460/EC), 2009 Preprint – Released for Publication p. 12 of 12 Evans, A. W., (2020),