Assessing trends in Duane plots using robust fits

Reliability Engineering and System Safety 90 (2005) 106–113 www.elsevier.com/locate/ress

Assessing trends in Duane plots using robust fits* Luis A. Garcı´a-Escuderoa, Miguel A. Ferna´ndeza,*, Oscar Duqueb, Angel Zoritab a

Departamento de Estadı´stica e Investigacio´n Operativa, Facultad de Ciencias, Universidad de Valladolid, C/Prado de la Magdalena s/n, 47071 Valladolid, Spain b Departamento de Ingenierı´a Ele´ctrica, E.T.S.I.I., Universidad de Valladolid, Paseo del Cauce s/n, 47071Valladolid, Spain Received 21 January 2004; accepted 21 March 2005 Available online 11 May 2005

Abstract The Duane plot is a simple and widely used graphical technique in the analysis of repairable systems. The fitting of a straight line to points in that graph serves to determine the behavior of the system assuming a power law process. However, the classical least squares fitting suffer from lack of robustness and it is specially heavily affected by the more remote and perhaps least interesting points. In order to solve this drawback we propose a robust fit based on Least Median of Squares (LMS) or Least Trimmed Squares (LTS) regression estimators of [1]. This methodology has been applied in a project aimed to study the reliability, availability and maintainability of the aerial contact line of the Spanish railways. q 2005 Elsevier Ltd. All rights reserved. Keywords: Duane plots; Power law process; Robustness; LMS; LTS

1. Introduction This work is motivated by the study of a highly dispersed and complex system which contains a large number of different components (the system is described with more detail in Section 3). In this system, the repair time of the components may be considered negligible with respect the total time (time units are ‘months’ and few hours are needed to take the system up when a failure happens). We also assume that a ‘renewal’ (or perfect) repair does not take place (i.e. the system is not brought to a like new state after a repair because only the failed component is replaced or repaired). So, we consider ‘minimal’ repairs that leave the system in exactly the same conditions as just before the failure happens. These assumptions lead us to the nonhomogeneous Poisson process as an appropriate model for our problem. In this paper we are primary interested in the behavior and evolution of each of the components of *

Research partially supported by FIMALAC Project. * Corresponding author. Tel./fax: C34 983 423 013. E-mail addresses: [email protected] (L.A. Garcı´a-Escudero), [email protected] (M.A. Ferna´ndez), [email protected] (O. Duque), [email protected] (A. Zorita). 0951-8320/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ress.2005.03.013

the system, that is, to know if the failure rate of the component is improving or worsening throughout the study period. The Duane plot turns out to be a very convenient tool for this sort of study. The Duane plot [2] is a widely used graphical technique in the analysis of repairable systems. This graph is simply a scatter plot of points fðlog ti ; logðNðti Þ=ti ÞÞgniZ1

(1)

where N(ti) is the cumulative number of failures through time ti (which is assumed to be instantaneous). These plots should be roughly linear if a power law process governs the failure times of the repairable system. For the power law process the expected number of failures through time t is E½NðtÞ Z


t b q

(2)

with q a scale parameter and b a shape parameter describing the behavior of the system (if bO1 the system is deteriorating as time passes and if b!1 the system is improving). More precisely, the power law process is a nonhomogeneous Poisson process with intensity function lðtÞZ ðb=qÞðt=qÞbK1 . The process N(t) is also called a Weibull process in the sense that the previous intensity function has the same form than the hazard rate for

L.A. Garcı´a-Escudero et al. / Reliability Engineering and System Safety 90 (2005) 106–113

a Weibull distribution. If bZ1 we can model our system through a standard Poisson process with constant rate of failures. From (2), we have that log E½NðtÞ=t Z ðb K 1Þlog t K b log q: Consequently, this fact suggests the fitting of a straight line to points (1) and using the slope of this line as an estimator of bK1. Also an estimator of parameter q can be obtained from the intercept and the slope of this line. Apart from the easy graphical interpretation of this line, [3] show that the least squares estimator is almost as efficient as the maximum likelihood estimator of b. However, we believe that there are some drawbacks for the use of the least squares approach. First, if we analyze the expected values of points in a Duane plot they are not exactly linear. If Ti is the random variable reflecting the occurrence time of the i-th failure, then ðEðlog Ti Þ; Eðlog ðNðTi Þ=Ti ÞÞÞ   jðiÞ jðiÞ ; log i K log q K Z log q C ; i Z 1; .; n b b where j denotes the digamma function (i.e. the derivative of the Euler’s gamma function). The derivation of this fact appears in [4] where also appear graphs showing that these expected values are only approximately linear in the ‘longterm’ and that the first expected values are far from being exactly linear. For instance, Fig. 1 shows the corresponding graph when qZ1 and bZ1.5. The least squares fit fails to catch the final linear trend, yielding to a worse determination of 1-b. The LS estimated slope is equal to 0.3095. This value is further from the real value 0.5 than the value 0.4469 obtained using the robust LTS procedure described later in this paper. Moreover, from a sample viewpoint, we will see how the first points in the graph are likely to be influence points

Expected values of log(N(T) / T)

1.0 Robust Non robust

0.9 0.8

107

determining heavily the slope of the regression line. These initial points are somehow the least interesting ones because, in general, we are mainly interested in the behavior of the system near the end of the period of time (because it is closely related with the present time in which we are working or we are willing to make inferences). The reason why these points tend to be influence points of our regression is that, due to the ‘log’ transformation on the time scale, the points on the left of the Duane plot (very old data) are more scattered than the points in the right of the graph (closer data) in the sense that they are farther from the center of the transformed data. Therefore, these points are frequently outliers in ‘x’ coordinate, and it is well known that this kind of outliers easily constitute influence points in Regression Analysis. Notice that, given data points fðxi ; yi ÞgniZ1 for a simple regression, the leverage hii defined as
2 1 ðx K xÞ hii Z C Pn i n
2 jZ1 ðxj K xÞ points out which are the observations that potentially exert an undue influence on the least squares regression coefficients (see e.g. [5]). If we now consider xiZlog ti the first observations will be those with largest leverages as
2 they are more remote from the center of data (i.e. ðxi K xÞ is large). Fig. 2 is intended to make this point clear. Fig. 2 shows the estimation of b (‘slope’ plus 1) for a simulated data set from a power law process with bZ1.5 and qZ100. The 15 failures appeared at 94, 139, 188, 214, 269, 280, 320, 354, 379, 420, 427, 454, 495, 519 and 534 time units. Suppose that we decrease 40 units the lowest value 94 (Fig. 2a) or that we increase 40 units this value (Fig. 2b), and we do the same for the highest value 534 (Fig. 1c and d, respectively). These figures make apparent the change in the slope estimations, b’s. Notice that the change in the first observation implies a more dramatic change than the same variation in the last observation. In order to overcome this drawback we propose the use of a robust fit to points (1) which will not be affected by these initial influence points. The robust fit will be obtained through the consideration of the Rousseeuw’s LMS and LTS estimators ([1,6] and [7]). The choice of these two robust fits will be discussed in Section 2.

0.7 0.6

2. LMS and LTS estimators

0.5

In the linear regression model we assume that the i-th observation can be written as

0.4

logðNðti Þ=ti Þ Z b0 C b1 log ti C 3i ; 0.0

0.5

1.0

1.5

2.0

Expected values of log(T)

Fig. 1. Expected values for the points in a Duane plot when qZ1 and bZ1.5. A robust and a nonrobust fit also appear.

with 3i being the error term. Given fitted values b^ 0 and b^ 1 , the associated residual is defined as ri ðb^ 0 ; b^ 1 Þ Z logðNðti Þ=ti Þ K b^ 0 K b^ 1 log ti :

L.A. Garcı´a-Escudero et al. / Reliability Engineering and System Safety 90 (2005) 106–113

108

(a) Slope = 0.260

(b) Slope = 0.676 –3.6

log(N(t) / t)

log(N(t) / t)

–3.6

–3.8

–4.0

–4.0

–4.4

–4.8

–4.2 4.0

4.5

5.0

5.5

6.0

5.0

5.2

5.4

log(t)

5.6

5.8

6.0

6.2

log(t)

(c) Slope = 0.543

(d) Slope = 0.508 –3.6

log(N(t) / t)

log(N(t) / t)

–3.6

–4.0

–4.0

–4.4

–4.4 4.5

5.0

5.5

6.0

4.5

5.0

5.5

6.0

log(t)

log(t) Fig. 2. Slope estimations for different changes in data points.

From the fitted values, b^ 0 and b^ 1 , we construct the ^ 1C b^ 1 , and qZ ^ expðKb^ 0 =ð1C b^ 1 ÞÞ for the estimators bZ shape and the scale parameters of the power law process. The (classical) least squares estimator is based on the minimization of min b0 ;b1

n X

ri2 ðb0 ; b1 Þ;

iZ1

but the lack of robustness of this method is well known. For instance, the change of just a single observation could make b^ 1 arbitrarily large (in robustness language, the estimator has breakdown point equal to 0). Instead of considering this objective function we could obtain the estimators b^ 0 and b^ 1 as the solutions of the following expression: arg min med ri2 ðb0 ; b1 Þ b0 ;b1

i

It is not difficult to see that the fitted line may be geometrically seen as the middle of the narrowest string containing at least half of the data points in (1). So, the LMS estimator fits the majority of the data and it is not affected by the presence of few anomalous observations achieving a breakdown point equal to 1/2 (the perturbation of more than a half of the data points is needed for making the estimators of b0 and b1 completely useless). Another possibility is the use of the LTS estimator. For given values b0 and b1, let rð1Þ ðb0 ; b1 Þ% rð2Þ ðb0 ; b1 Þ% .%

rðnÞ ðb0 ; b1 Þ be the ordered residuals. The LTS estimators are now based on the minimization problem: arg min b0 ;b1

h X

2 rðiÞ ðb0 ; b1 Þ

iZ1

where h is a trimming factor which can be chosen by the analyst. Therefore the LTS regression estimate is the LS fit to those h observations with the smallest sum of squares. The breakdown point of the LTS is also a half when hZ ½ðnC 2C 1=2Þ where [z] denotes the integer part of z. An interesting property of these two robust methods is the ‘exact fit’ property which means that if the majority of the data points are perfectly aligned then the line fit is given by this alignment. The least square fit does not have this property. Fig. 3 shows the result of applying a robust fits to the same data sets leading to Fig. 2. We can see that these robust fits are not seriously affected by changes in few data values as they try to fit the majority of data. Specially, notice that changes in remote times are not necessarily more harmful than changes in recent periods of time. Finally, notice that the LTS fits are quite similar to those obtained by the LMS procedure. The statistical efficiency of the LTS is better than the LMS because it is asymptotically normal. Additionally, rather efficient algorithms have been recently proposed for its computation (see, [8]). So, in what follows, we will focus on the LTS estimation.

L.A. Garcı´a-Escudero et al. / Reliability Engineering and System Safety 90 (2005) 106–113

(b) LTS slope = 0.427

–3.6

–3.6

–3.8

–4.0

log(N(t) / t)

log(N(t) / t)

(a) LTS slope = 0.427

109

–4.0

–4.4

–4.8

–4.2 4.0

4.5

5.0

5.5

6.0

5.0

5.2

5.4

log(t)

5.6

5.8

6.0

6.2

log(t)

(c) LTS slope = 0.422

(d) LTS slope = 0.423 –3.6

log(N(t) / t)

log(N(t) / t)

–3.6

–4.0

–4.0

–4.4

–4.4 4.5

5.0

5.5

6.0

4.5

5.0

5.5

log(t)

6.0

log(t)

Fig. 3. Robust fits for the same data sets as in Fig. 2. The solid lines correspond to the LTS fits and the dashed lines to the LMS fits.

The maximum likelihood estimator of parameter b for the power law process is n : logðt n =ti Þ iZ1

b^ MLE Z Pn

n=20

Some references on maximum likelihood estimation (MLE) in power law processes are [9–13]. In order to compare the performance of the proposed procedure (LTS) with respect to the estimation based on least square (LS) fits and the MLE, we have simulated failures times from

n=40

n=100

2.5

2.5

2.5

2.0

2.0

2.0

1.5

1.5

1.5

1.0

1.0

1.0

0.5

0.5

0.5 MLE

LTS

LS

MLE

LTS

LS

MLE

LTS

LS

Fig. 4. Box-plots for the estimation of bZ1.5 using Maximum Likelihood (MLE), Least Trimmed Squares fit-based (LTS) and Least Squares fit-based (LS) estimators.

L.A. Garcı´a-Escudero et al. / Reliability Engineering and System Safety 90 (2005) 106–113

110 Table 1

Intercept) log(t)

Value

Standard error

t Value

Pr(Ojtj)

K6.5095 0.4661

0.1246 0.0212

K52.2378 22.0095

0.0000 0.0000

a power law process with qZ1 and bZ1.5. The resulting box-plots for the estimations of parameter b with these three methods based on 4000 power law process simulations are shown in Fig. 4. Different number of failures times nZ20, 40 and 100 have been considered. As it was expected, MLE provides the best results for large n due its asymptotic optimality. However, the LTS fit seems to provide a less biased estimation of parameter b especially for small sample sizes. The result of the classical LS fit was not good for any of the sample sizes considered and its use is not advisable unless we consider very high n’s. We would like to stress that, in these simulations, we are always sampling from the ‘ideal’ power law process without the presence of anomalous failure times. In a ‘nonideal’ case, the robust fit would be preferable as it will be less affected by the nonfitting observations. We try to justify this last claim in Section 3. The ability of doing (statistical) inference (tests, confidence intervals,.) based on these procedures is also important. This can be done using classical weighted least ^ with ri being squares with weights given by ui Z uðjri =sjÞ the residual from the LTS fit and u($) a decreasing function

(see [7]). For instance, we could give weights equal to 1 to observations with reasonably small residual (e.g. jri j% 2:5s^ with s^ a robust estimation of the scale) and 0 weights to points with large residuals and, then, apply weighted least squares. For the 15 failure times previously considered (without any modification), we can follow this approach obtaining. Table 1 which permits us to perform significance tests or confidence intervals. Finally, another advantage of these kind of estimators is that both LMS and LTS are implemented in standard statistical software as the latest versions of S-PLUS and SAS. The LMS has been also implemented in MATHEMATICA (see [14]). The weighted least squares regression is also implemented in most of the statistical packages.

3. Example The examples appearing in this section are taken from a project named FIMALAC. The aim of this project is to study the reliability, availability and maintenance of the aerial contact line of the Spanish railway company RENFE. Part of this project is the study of the reliability of the different components appearing in the contact line system. The procedures developed in this paper are designed to help to achieve one to the final aims of the project which is (b) LTS-based estimation of beta = 0.782

0.4

0.4

0.2

0.2

0.0

0.0

log(N / t)

log(N / t)

(a) LS-based estimation of beta = 0.871

–0.2

–0.2

–0.4

–0.4

–0.6

–0.6

1

2

3

4

1

3

2 log(t)

log(t) Fig. 5. A11 Isolators.

4

L.A. Garcı´a-Escudero et al. / Reliability Engineering and System Safety 90 (2005) 106–113

(a) LS-based estimation of beta = 1.093

(b) LTS-based estimation of beta = 0.559 0.0

–0.5

–0.5

–1.0

–1.0 log(N / t)

0.0

log(N / t)

111

–1.5

–1.5

–2.0

–2.0

–2.5

–2.5 2.5

3.0

3.5

4.0

4.5

2.5

3.0

3.5

4.0

4.5

log(t)

log(t) Fig. 6. Rt-51 spirelec isolators.

to apply RCM to this kind of large railway systems. Some very interesting recent references in this line are [15,16,17]. The data have been collected in the whole Spanish railway network from January 1990 until December 2000. The database stores for each of the failures when, where, how and why did the failure appear. There is a high level of detail which allows us to study each component in the system. In this section we show the Duane plots for three components of the contact line system, namely isolators A11 (Fig. 5), Rt-51 Spirelec (Fig. 6) and the feeder (Fig. 7). It is clear that fitting the data using LS (see Fig. 5a) may lead us to a very poor estimator of the b parameter due to the appearance of the influence points described in the previous sections. On the other hand the LTS proposed approach (see Fig. 5b) yields estimators which are not so influenced by the most ancient observations placing more importance to the latest ones and, therefore, giving us more accurately the present behavior of the components. Notice in Fig. 6 that the unadvised analyst could make an important mistake from the LS estimator of the slope. He might conclude that the failure rate of the component is increasing while it has been decreasing for some time. This could lead us to a wrong planning decision such as to modify the maintenance strategy thinking that is inappropriate when, precisely, it could turn out to be a suitable one. It could even influence the investment decisions, planning the replacement of the component.

The system analyzed had high complexity due to the consideration of many components along a large period of time with high heterogeneity. Although the methods developed in this paper can be used for the important task of estimating the parameters of a power law process, where appropriate, this was not the present aim in this project. At this early stage of the project the staff of RENFE was primarily interested in deciding if the present maintenance strategies for each of the components in the system were leading to a reduction of the frequency of failures. Of course, we do not argue that the power law process is the most appropriate model for modelling this data, but the use of these robust fits in the associated Duane plots was extremely useful for pointing out these trends component by component. Moreover, the system is composed of a great quantity of different components (over 50 types) and structured into a huge number of geographical zones (over 100) each of them having its own maintenance responsibilities. Thus, it was convenient to develop rather unsupervised method that allows evaluating and comparing the results of the maintenance strategies applied to each zone and type of component. Therefore, it was quite interesting to perform an hypothesis tests (in an automated fashion) like (

H0 : b% 1 H1 : bO 1

L.A. Garcı´a-Escudero et al. / Reliability Engineering and System Safety 90 (2005) 106–113

112

(a) LS-based estimation of beta = 0.802

(b) LTS-based estimation of beta = 1.497

0.5

0.5 log(N / t)

1.0

log(N / t)

1.0

0.0

0.0

–0.5

–0.5 0

1

3

2

4

5

0

1

3

2

4

5

log(t)

log(t) Fig. 7. Feeder.

that would allows us to detect a significative deterioration of any component. For doing this task, we could apply the one-side significance test from the reweighed least squares regression with weights coming from a LTS fit to the Duane plot points. We compare the results with the classical MLE-based tests that rejects the null hypothesis whenever b^ MLE O

2n c22ðnK1Þ;1Ka

where c2n;a is the value that leaves a probability a in the right tail of the c2 distribution with n degrees of freedom. Although, the results were quite similar in most of the cases, there were some cases where the results obtained with the LTS fit appear to be more convincing than the ones derived using the MLE. Just as an example, Fig. 7 shows the Duane plot corresponding to the feeder. The estimation of the b parameter using the LTS fit was equal to 1.497 leading to a p-value !0.0001 for rejecting H0 and detecting a worsening in that system. However, we obtained b^ MLE Z 1:1630 with a p-valueZ0.0602 (provoking the nondetection of the deterioration in that system by using the classical significance levels: 0.01 or 0.05). It is true that a careful revision of each Duane plot would surely reveal that feature, but it might be overlooked if a nonrobust unsupervised procedure is used.

4. Conclusions In this paper we have studied possible problems in the interpretation of Duane plots. Since the Duane plot is a very useful technique frequently used for exploratory purposes a mistake in assessing the slope of the plot can lead to take an erroneous maintenance decision. Our first contribution is to show why this sort of mistakes are likely to appear while the second one is to explain how those hazardous situations can be dealt with using robust fits that can be easily implemented in standard statistical software. This technique has proved to be very useful applied for the analysis of the maintenance strategies of a complex and extended system such as the aerial contact line of a National railway network, allowing to modify some of the current strategies and even speed up the replacement of components.

Acknowledgements The authors wish to thank Roberto Lo´pez de la Cruz and Manuel Carmona from RENFE for helpful comments. We also want to thank the associate editor and two anonymous referees for their careful reading of the paper which lead to this improved version.

L.A. Garcı´a-Escudero et al. / Reliability Engineering and System Safety 90 (2005) 106–113

References [1] Rousseeuw PJ. Least median of squares regression. J Am Stat Assoc 1984;79:871–80. [2] Duane JT. Learning curve approach to reliability monitoring. IEEE Trans Reliab 1964;R-21:563–6. [3] Molitor JT, Rigdon SE. A comparison of estimators of the shape parameter of the power law process. Proceedings of the section on physical and engineering sciences. Alexandria, VA: American Statistical Association; 1993 pp. 143–147. [4] Rigdon SE, Basu AP. Statistical methods for the reliability of repairable systems. New York: Wiley; 2000. [5] R.H. Myers, Classical and Modern Regression with Applications, PWS-KENT; 1990. [6] Rousseeuw PJ. Multivariate estimation with high breakdown point. In: Grossmann W, Pflug G, Vincze I, Wertz W, editors. Mathematical statistics and applications. Dordrecht: Reidel Publishing Company; 1985. p. 283–97. [7] Rousseeuw PJ, Leroy AM. Robust regression and outlier detection. New York: Wiley; 1987. [8] Rousseeuw PJ, Van Driessen K. An algorithm for positive-breakdown methods based on concentration steps. In: Gaul W, Opitz O, Schader M, editors. Data analysis: scientific modelling and practical application. New York: Springer; 2000. p. 335–46. [9] Finkelstein JM. Confidence bounds on the parameters of a Weibull process. Technometrics 1976;18:115–7.

113

[10] Crow LH. Confidence interval procedures for the Weibull process with applications to reliability growth. Technometrics 1982; 24:67–70. [11] Bain LJ, Engelhardt M, Wright FT. Tests for an increasing trend in the intensity of a Poisson process: Apower study. J Am Stat Assoc 1985; 80:419–22. [12] Guida M, Calabria R, Pulcini G. Bayes inference for a nonhomogeneous poisson process with power intensity law. IEEE Trans Reliab 1989;38:603–9. [13] MIL-HDBK-781. Handbook for reliability test methods, plans, and environments for engineering, development qualification, and production. Washington, DC: United States Government Printing Office; 1997. [14] Shaw WT, Tigg J. Applied mathematica. Reading, MA: AddisonWesley; 1994. [15] Carretero J, Pe´rez JM, Garcı´a-Carballeira F, Caldero´n A, Ferna´ndez J, Garcı´a JD, Lozano A, Cardona L, Cotaina N, Prete P. Applying RCM in large scale systems: a case study with railway networks. Reliability Engineering and System Safety 2003;82:257–73. [16] Pedregal DJ, Garcı´a FP, Schmid F. RCM2 predictive maintenance of railway systems based on unobserved components models. Reliability Engineering and System Safety 2004;83:103–10. [17] Garcı´a Ma´rquez FP, Engelhardt F, Schmidt FT, Conde Collado J. A reliability centered approach to remote condition monitoring. A railway points case study. Reliability Engineering and System Safety 2003;80:33–40.