A Review of Statistical Failure Time Models with Application of a Discrete Hazard Based Model to 1Cr1Mo-0.25V Steel for Turbine Rotors and Shafts

Producing predictions of the probabilistic risks of operating materials for given lengths of time at stated operating conditions requires the assimilation of existing deterministic creep life prediction models (that only predict the average failure time) with statistical models that capture the random component of creep. To date, these approaches have rarely been combined to achieve this objective. The first half of this paper therefore provides a summary review of some statistical models to help bridge the gap between these two approaches. The second half of the paper illustrates one possible assimilation using 1Cr1Mo-0.25V steel. The Wilshire equation for creep life prediction is integrated into a discrete hazard based statistical model—the former being chosen because of its novelty and proven capability in accurately predicting average failure times and the latter being chosen because of its flexibility in modelling the failure time distribution. Using this model it was found that, for example, if this material had been in operation for around 15 years at 823 K and 130 MPa, the chances of failure in the next year is around 35%. However, if this material had been in operation for around 25 years, the chance of failure in the next year rises dramatically to around 80%.


Introduction
The prediction of long-term creep properties from short timescale experiments is rated as the most important challenge to the UK Energy Sector in a recent UK Energy Materials Review [1]. Creep strain (ε) is a function not only of stress (τ) and absolute temperature (T), but also of time (t) After an initial strain on loading, a decaying creep rate ( . ε) during the primary stage of creep is followed by an accelerating strain during the tertiary stage. A minimum creep rate ( . ε m ) occurs at the boundary of these two stages. As such, Equation (1a) is often represented in differential form . ε = f 2 (τ, T, t) (1b) When it comes to extrapolating from short term accelerated test data, three very broad approaches can be identified. Whole creep curve methods work by relating the whole creep curve to the test conditions under which that creep curve was obtained. A single creep curve at steady uniaxial stress τ and absolute temperature T can be modelled using a general functional form ε = η t, Ψ 1 , Ψ 2 , . . . . , Ψ j , . . . . , Ψ q (2a) where η is some non-linear function and the Ψ j are numerical parameters. At any one test condition, the Ψ j parameters are constant, but they do vary systematically with the test conditions. It is this fact that enables creep curve predictions to be made Ψ j = g j (τ, T, b j,1 , b j,2 , . . . , b j,k , . . . ., b j,p ) where g j are non-linear functions, and b j,k are additional numerical parameters that can be estimated using a suitable estimation technique. However, the form of the η and g j functions are not known and consequently the literature contains many representations of these including, for example, the Theta methodology proposed by Evans and Wilshire [2] Whilst Evans [3] derived Equation (3a) from creep deformation mechanism theory, the all important extrapolation function given by Equation (3b) is mainly empirical in nature. There are many other approaches in the literature including those by McVetty [4], Garofalo [5], Ion et al. [6], Prager [7], Othman and Hayhurst [8], Kachanov [9] and Rabotnov [10].
Secondly, parametric techniques work by relating a measured point on the creep curve to the test conditions under which that measurement was made. This point is typically the minimum creep rate or the time to rupture. Around the minimum creep rate there remains a considerable period of time where . ε remains more or less constant, so that Equation (1b) reduces to with, following Monkman and Grant [11], where here t represent the time at which failure occur. Because the functional form of f 3 is not known the literature again contains many different representations of Equations (4a) and (4b). For example, Dorn [12] and Larson and Miller [13], both assumed that at a constant stress where Q c is the activation energy in J/mol. These two approaches then diverge with the incorporation of stress with Dorn suggesting the parameter C 0 varies with stress, whilst the Larson and Miller model has Q c varying with stress-but in both cases the form of the stress function was empirically specified (typically involving the use of polynomials in stress or the log of stress). The literature contains many other variations including Manson and Haferd [14], Manson and Muraldihan [15], Manson and Brown [16] and Trunin et al. [17]. Unfortunately, all these parametric models suffer from parameter instability with respect to stress and temperature making reliable long term life predictions from accelerated short term testing impossible-as empirically illustrated by Abdallah et al. [18]. These empirical models are now quite old and, despite their known short comings, are still extensively used for safe life estimation. The hyperbolic tangent method [19] and the Wilshire Equation [20] can be seen as the most recent types of parametric model, with the later having the form where R is the universal gas constant, τ* = ln(−ln(τ/τ TS )) with τ TS being the tensile strength. Unlike the above models, a raft of recent publications [21][22][23][24][25][26] on a wide range of high temperature materials have demonstrated the parameters of this model (b 0 , b 1 , and Q c ) are stable and so reliable long term predictions have been made from this model using short term data of no more than 5000 h duration. Finally, there are computational/numerical approaches that often incorporate detailed deformation mechanism into finite element code to obtain creep property predictions. Many of these numerical models are based on remaining life assessment with abridged accelerated testing. In some of these approaches, for example [27,28], uniaxial test specimens are cut from removed components that have been in service for over a long time (typically over 100,000 h) and re-tested under laboratory conditions for short times until failure occurs by accelerating the temperatures (but using the in service stress). Such testing yields the remaining life or useable residual life (around 60% of remaining life). The above parametric models are then used to predict/extrapolate these residual lives to operating temperatures. Often numerical models are used to extrapolate such abridged short-term testing. An alternative to this destructive approach is the non-destructive disc test, where small discs are taken from in service components without destroying their integrity. Again numerical models can be built for this disc test, e.g., Evans [29] or parametric procedures can be used for extrapolative purposes.
What all the above studies have in common however is that they are all deterministic in nature. As the level of stress increases the time to failure diminishes and the primary component of the creep curve becomes less pronounced. These variations in creep properties are governed by (as yet not fully understood) physical laws that can be used to determine creep properties at any test condition. These physical laws are embedded into the above mechanistic models that can then be used to explain variations in creep properties as a function of test conditions alone.
However, creep is not just a deterministic process that is predetermined by physical laws. It is also a random process. If the time function, f 1 , could be quantified, it could then be used to predict the strain at any time. Unfortunately, the nature of creep is such that this function could not then be used to predict the strain for another specimen tested under exactly the same conditions. This random component of creep is in turn very large. For example, in the NIMS [30] 1Cr-1Mo-0.25V steel database used in this paper, the time to failure at 773 K and 373 MPa varies between the limits of 125 h and 1360 h depending on the batch. A lot of the random variation seen in creep databases of this nature are down to variations in chemical composition and heat treatments experienced by the different batches of the test material. However, even when such factors are removed, failure times are still highly stochastic in nature. This was illustrated by Evans [31], who tested 15 specimens of Ti-6.2.4.6 at different temperatures and stresses. These test specimens were all cut from the same batch and tested within the same laboratory on the same make of calibrated uniaxial test apparatus. The results are reproduced in Figure 1 where it can be seen that the random variation is great enough to encompass variations induced by changes in stress.
This random component of creep is not just important because it is large in size, but also because if predictions are to be made for more than just the average time at failure, then this random component must be modelled with the same degree of vigour as displayed by all the above deterministic models. However, they are not. We can summarise the above as stating creep life has both a deterministic and a random component where u is the random component and f 4 is some function that describes how the random component is distributed.
The first aim of this paper is to address this shortcoming by providing a detailed, although by no means complete (as this is a very large subject area), review of statistical failure time models that describe different ways that f 4 (u) can be specified, so as to provide materials scientists with a framework for further developing their deterministic creep models, such as the Wilshire Equation, so that they become cable of providing predictions that have levels of confidence attached to them. It should be noted that these statistical models say nothing about the deterministic models reviewed above and do not imply that one model is any better at prediction than another. However, when a model of the random component is combined with the above deterministic components, they become more capable of predicting both the systematic variation with test conditions and the observed random variation at each test condition. Without this, the deterministic models can do no more than accurately predict the average safe life and not the safe life corresponding to say a 1% chance of failure. This review is provided in Section 3 of this paper and where appropriate illustrated using data on 1Cr-1Mo-0.25V steel.
Materials 2017, 10, 1190 4 of 30 random variation at each test condition. Without this, the deterministic models can do no more than accurately predict the average safe life and not the safe life corresponding to say a 1% chance of failure. This review is provided in Section 3 of this paper and where appropriate illustrated using data on 1Cr-1Mo-0.25V steel. The second aim of this paper is to illustrate one of the many ways that these deterministic and random component models can be combined. It is impossible in one paper to do an illustration for all of the above deterministic models reviewed above, and so the Wilshire model is selected for this purpose. This approach is selected because it has been shown to outperform the others in terms of accurately predicting the average time to failure beyond 100,000 h using very short term data (less than 5000 h). The Wilshire equation is combined with a discrete hazard based model for the random component. A hazard based model was chosen because it offers extra flexibility on distributional shape and form compared to other approaches as discussed in the review section below. Further, a discrete version of the hazard function is used because it helps empirically quantify the form of the hazard function (and because it has never been used within the context of creep failure beforewhereas other approaches have [32]).
The review and illustration of combining deterministic and random creep models illustrated using the NIMS database on 1Cr-1Mo-0.25V steel [30]. This is carried out in Section 4, and conclusions are then drawn in Section 5.

The Data
This present study features forged 1Cr-1Mo-0.25V steel for turbine rotors and shafts. For multiple batches of this bainitic product, both the creep and creep fracture properties have been documented comprehensively by the National Institute for Materials Science (NIMS), Japan [30]. NIMS creep data sheet No. 9B includes information on nine batches of as tempered 1Cr-1Mo-0.25V steel. Each batch of material had both a different chemical composition and a different thermal and processing history-details of which can be found in creep data sheet No. 9B. Specimens for the tensile and creep rupture tests were taken radially from the ring shaped samples which were removed from the turbine rotors. Each test specimen had a diameter of 10 mm with a gauge length The second aim of this paper is to illustrate one of the many ways that these deterministic and random component models can be combined. It is impossible in one paper to do an illustration for all of the above deterministic models reviewed above, and so the Wilshire model is selected for this purpose. This approach is selected because it has been shown to outperform the others in terms of accurately predicting the average time to failure beyond 100,000 h using very short term data (less than 5000 h). The Wilshire equation is combined with a discrete hazard based model for the random component. A hazard based model was chosen because it offers extra flexibility on distributional shape and form compared to other approaches as discussed in the review section below. Further, a discrete version of the hazard function is used because it helps empirically quantify the form of the hazard function (and because it has never been used within the context of creep failure before-whereas other approaches have [32]).
The review and illustration of combining deterministic and random creep models illustrated using the NIMS database on 1Cr-1Mo-0.25V steel [30]. This is carried out in Section 4, and conclusions are then drawn in Section 5.

The Data
This present study features forged 1Cr-1Mo-0.25V steel for turbine rotors and shafts. For multiple batches of this bainitic product, both the creep and creep fracture properties have been documented comprehensively by the National Institute for Materials Science (NIMS), Japan [30]. NIMS creep data sheet No. 9B includes information on nine batches of as tempered 1Cr-1Mo-0.25V steel. Each batch of material had both a different chemical composition and a different thermal and processing history-details of which can be found in creep data sheet No. 9B. Specimens for the tensile and creep rupture tests were taken radially from the ring shaped samples which were removed from the turbine rotors. Each test specimen had a diameter of 10 mm with a gauge length of 50 mm. These specimens were tested at constant load over a wide range of conditions: 47-333 MPa and 723-923 K. In addition to failure time (t) measurements, values of the 0.2% proof stress (τ Y ) and the ultimate tensile strength (τ TS ) determined from high strain rate (~10 −3 s −1 ) tensile tests carried out at the creep temperatures for each batch of steel investigated were also reported.
The review section below (Section 3) is illustrated using all batches of data, whilst the discrete hazard based model outlined in Section 4 is illustrated using only a single batch of materials from this database-VaA.

A Statistical Description of Continuous Failure at Fixed Test Conditions
Due to batch to batch variations in chemistry and heat treatment and within batch variations in microstructure, creep failure times for a high temperature material (or indeed any given material), even under fixed test conditions, are stochastic in nature. Therefor such failure times need to be described through a random variable T. In reality, T can take on a large and continuous number of different values at a given test condition, t 1 , t 2 , . . . , t n , with 0 ≤ t 1 ≤ t 2 ≤ . . . ≤ t n . As such it cannot be known with certainty when failure will occur and so failure must be expressed using the survivor function which gives the probability of surviving beyond a certain length of time, S(t) The probability of failing at or before a given length of time is then given by F(t) = 1 − S(t). The probability of failure in a very small increment of time, ∆t, is then f(t) = ∆F(t)/∆t. F(t) is often referred to as the cumulative distribution function (cdf) and its derivative, f(t), the probability density function (pdf).
The probability of failure can also be expressed through the hazard function. This function gives the rate of failure at time t, given the specimen survives up to time t where |T ≥ t reads given that T is greater than or equal to t. As such, the hazard rate is a conditional probability of failure. A conditional probability is defined as P(A|B) = P(A and B)/P(B), where in terms of the hazard rate event B is the probability of surviving a length of time t and so equals S(t). Event A and B would then be the probability of failing in the small increment of time ∆t beyond t, which is the pdf at time t. Thus If follows from these definitions that the hazard function can also be found from the survivor function using and the cumulative (or integral) hazard function is given by Approaches to estimating the survivor function generally fall under three headings: parametric, non-parametric and semi-parametric. The assumption behind the parametric approach is that the form of the survivor function can be captured through a small number of parameters. For example, if failure times at a fixed test condition are normally distributed, then the survivor function is fully defined through two parameters-the mean and the standard deviation. In contrast, the non-parametric approach is model (parameter) free and as such makes no assumptions about how failure times are distributed. The semi-parametric approach combines these two approaches, for example, by specifying a base line hazard function at a particular test condition non-parametrically and then using a few parameters to model how this baseline function changes with the test conditions.

Non-Parametric Estimation
The starting point for many non-parametric techniques is to partition time into j = 1 to k equal intervals, with k being as large as practically possible. If n equals the number of specimens placed on test at the same test condition and d j the number of specimens failing during the kth interval, then Kaplan and Meier [33] proposed the following estimator of the survivor function (for uncensored data) that has as its basis the binomial distribution where d j is the number of failures in time interval j. This estimator is also referred to as the product-limit estimator as originally these authors justified this estimator based on its properties when k tended to infinity or as the time interval tended to zero. Nelson [34] and Aalen [35] proposed the following non-parametric estimator of the cumulative hazard functionΛ where r j is the total number of specimens at risk (or not yet failed) just prior to time t i . The Fleming-Harrington [36] estimator of the survivor function is, from Equations (7e) and (8b), The above are of course estimates (designated by the hat symbol) of the survivor function computed in the above ways, but from a population or very large sample. The standard deviation of the above estimators provides a way to quantify the possible size of the difference between the true or population survivor function-S(t)-and that calculated from a small sample or a randomly selected sub set of the population. The standard deviation of these estimators are in turn estimated bŷ In large samples, these estimates are unbiased and the Nelson Aalen estimator is then also approximately normally distributed.
As an illustration, Figure 2 compares these non-parametric estimators of the survivor function for the 1Cr-1Mo-0.25V specimens in the NIMS dataset tested at 823 K and 294 MPa. At low times to failure the above two estimators provide very similar values for the survivor function, but these estimators start to diverge at around 250 h-with the Fleming-Harrington estimator exceeding the Kaplan-Meier estimator. The Nelson-Aalen estimator of the cumulative hazard function is shown on the right hand side vertical axis. The errors bars associated with the estimated cumulative hazard function, which are made equal to one standard deviation, are also shown. As can be seen, the standard error increases quite dramatically with the time to failure, making the estimates at high survival probabilities quite unreliable in a sample this small. where, like T and Y, Z is a random (but standardised) variable, z ∈ (−∞, ∞). To prevent the occurrence of a degenerate distribution for large values of k1 and/or k2, the following re-parameterisation is used , σ = b/δ and where W is therefore another standardised random variable defined as W = δZ. By specifying a very general distribution for Z, it is possible to identify many of the familiar failure time distributions used in failure time analysis. Prentice [37] and Kalbfleisch and Prentice [38] for example defined the probability density function for Z as where Z is said to be distributed as the logarithm of an F random variable with 2k1 and 2k2 degrees of freedom. T is described as following a four parameter generalised F distribution, T ~ GENF(μ, σ, k1, k2). Г(k) is the gamma function at k. The Appendix A to this paper also shows that the pdf of this generalised gamma distribution can be re-parameterised as a function of time

Parametric Estimation
Any distribution defined for t ∈ (0, ∞) can be used to specify parametric survivor and hazard rate functions. A good transformation for visualising many commonly used parametric distributions is the log transformation of failure time, Y = ln(T), with y ∈ (−∞, ∞). Then, a whole family of distributions for Y opens up by introducing location (via parameter µ) and scale (via parameter b) changes of the form where, like T and Y, Z is a random (but standardised) variable, z ∈ (−∞, ∞). To prevent the occurrence of a degenerate distribution for large values of k 1 and/or k 2 , the following re-parameterisation is used where δ = (k 1 k 2 /(k 1 + k 2 ), σ = b/δ and where W is therefore another standardised random variable defined as W = δZ. By specifying a very general distribution for Z, it is possible to identify many of the familiar failure time distributions used in failure time analysis. Prentice [37] and Kalbfleisch and Prentice [38] for example defined the probability density function for Z as where Z is said to be distributed as the logarithm of an F random variable with 2k 1 and 2k 2 degrees of freedom. T is described as following a four parameter generalised F distribution, T~GENF(µ, σ, is the gamma function at k. The Appendix A to this paper also shows that the pdf of this generalised gamma distribution can be re-parameterised as a function of time where λ = exp(−µ) and β = 1/(δσ) = 1/b. Except, under some restricted values for k 1 and k 2 , there is no closed form expression for the survivor and hazard functions, but they are related to the incomplete beta function and Appendix A shows how this can be computed using percentiles from the F distribution. Equation (9d) is however degenerate when k 1 = k 2 = ∞ and then a different specification of the pdf must be used (see Appendix A).
Particular values for these parameters define important sub families within the GENF family and these sub families are summarised in Figure 3. It can be seen that some of these distributions are commonly used within engineering. When k 2 = ∞, failure times have a Generalised Gamma distribution, T~GENG(µ, σ, k 1 ). There are then three well known two parameter distributions within this Generalised Gamma family. T is gamma distributed, T~GAM(µ, σ, k 1 ), when k 2 = ∞ and σ = 1. T is log normally distributed, T~LOGNOR(µ, σ), when k 2 = k 1 = ∞; and T is Weibull distributed, T WEIB(µ, σ), when k 2 = ∞ and k 1 = 1. In turn, the Weibull distribution collapses to the exponential distribution when k 2 = ∞, k 1 = 1 and σ = 1. The family, T~BURR(µ, σ, k 1 ), is obtained when either k 1 = 1 (Burr III) or k 2 = 1 (Burr XII). Then when k 2 = k 1 = 1, T has a log-logistic distribution, T LOGLOGIS(µ, σ), and when k 2 = k 1 = 1 = σ the log-logistic distribution collapses to the logistic distribution, T~LOGIS(µ). The form and characteristics of all these special cases are further described in the Appendix A.
where λ = exp(−μ) and β = 1/(δσ) = 1/b. Except, under some restricted values for k1 and k2, there is no closed form expression for the survivor and hazard functions, but they are related to the incomplete beta function and Appendix A shows how this can be computed using percentiles from the F distribution. Equation (9d) is however degenerate when k1 = k2 = ∞ and then a different specification of the pdf must be used (see Appendix A).
Particular values for these parameters define important sub families within the GENF family and these sub families are summarised in Figure 3. It can be seen that some of these distributions are commonly used within engineering. When k2 = ∞, failure times have a Generalised Gamma distribution, T ~ GENG(μ, σ, k1). There are then three well known two parameter distributions within this Generalised Gamma family. T is gamma distributed, T ~ GAM(μ, σ, k1), when k2 = ∞ and σ = 1. T is log normally distributed, T ~ LOGNOR(μ, σ), when k2 = k1 = ∞; and T is Weibull distributed, T ~ WEIB(μ, σ), when k2 = ∞ and k1 = 1. In turn, the Weibull distribution collapses to the exponential distribution when k2 = ∞, k1 = 1 and σ = 1. The family, T ~ BURR(μ, σ, k1), is obtained when either k1 = 1 (Burr III) or k2 = 1 (Burr XII). Then when k2 = k1 = 1, T has a log-logistic distribution, T ~ LOGLOGIS(μ, σ), and when k2 = k1 = 1 = σ the log-logistic distribution collapses to the logistic distribution, T ~ LOGIS(μ). The form and characteristics of all these special cases are further described in the Appendix A. Evans [32] has shown how the parameters of these distributions can be estimated using maximum likelihood procedures. An alternative semi-parametric approach is to use the least-square procedure in conjunction with a probability plot. The procedure here is to linearise a plot of t against S(t) by finding suitable transformations of S(t) and possibly t. A least squares best fit line to the data on such a plot then yields estimates of the parameters μ (given by the intercept of the best fit line) and σ (the slope of the best fit line). However, as seen in Figure 1, the non-parametric estimate ) ( i km t S is a step function increasing by an amount 1/n at each recoded failure time. Plotting at the bottom (top) of the steps would lead to the best fit line being above (below) the plotted points and so Evans [32] has shown how the parameters of these distributions can be estimated using maximum likelihood procedures. An alternative semi-parametric approach is to use the least-square procedure in conjunction with a probability plot. The procedure here is to linearise a plot of t against S(t) by finding suitable transformations of S(t) and possibly t. A least squares best fit line to the data on such a plot then yields estimates of the parameters µ (given by the intercept of the best fit line) and σ (the slope of the best fit line). However, as seen in Figure 1, the non-parametric estimateŜ km (t i ) is a step function increasing by an amount 1/n at each recoded failure time. Plotting at the bottom (top) of the steps would lead to the best fit line being above (below) the plotted points and so lead to a bias in the resulting parameter estimates. A reasonable compromise plotting position is the mid-point of the jump where i indexes the ordered failure times (i = 1 for the smallest failure time, i = 2 for the next smallest all the way up to n for the largest failure time), with t 1 being the smallest failure time up to t n the largest failure time. From the Appendix A to this paper, the log of the pth percentile for t is given by where w k1,k2,p is the pth quantile of an F distribution with (2k 1 , 2k 2 ) degrees of freedom. Percentiles of the F distribution are tabulated at the back of many well know engineering statistical text books (it can also be found in Excel using the FINV function). Usingp i in Equation (4a) for p in Equation (10b) allows values for w k1,k2,p to be computed. Thus, when w k1,k2,p i is plotted against the ordered values for ln(t), ln(t i ), the data points will reveal scatter around a linear line provided the data have a generalised F distribution with given values for k 1 and k 2 .
The generality of Equation (10b) is clearly seen by considering the special case of k 2 = ∞ and k 1 = 1, which is the Weibull distribution, whose survivor function is shown in the Appendix A of this paper to be This can be linearised as Then, replacing S(t) with the parametric estimatorp i and t by its ordered value t i gives Equations (11c) and (10b) imply that w k1,k2,p collapses to ln{−ln[S(t)]} when k 2 = ∞ and k 1 = 1.
As an illustration, Figure 4a is a ln(t i )/ ln{− ln[p i ]} plot for the ten 1Cr-1Mo-0.25V specimens tested at 823 K and 294 MPa. The best fit line obtained using the least squares technique is also shown. The slope of this best fit line is σ = 0.2938, with intercept µ = 5.8385. This implies β = 1/σ = 3.4 and λ = exp(−5.8385) = 0.0029. However, with a coefficient of determination (R 2 ) of just 85%, the Weibull distribution is unlikely to be the best description of this sample of failure times. This R 2 value was computed over the range p = 0 to 2 and q = 0 to 1 (both in increments of As such, this range covered all the distributions shown in Figure 3. It was found that R 2 was maximised when p = q = 0, i.e., when k 1 = k 2 = ∞. It therefore appears that, within the generalised F distribution family, it is the log normal distribution that best describes the specimens tested at 823 K and 294 MPa. This is consistent with the findings by Evans [32]. Figure 4b plots ln(t i ) against w k1,k2,p i when k 1 = k 2 = ∞, so that the variable on the horizontal axis is essentially a standard normal variate. The R 2 value is 93% and so much higher than the Weibull case. The slope of this best fit line is σ = 0.3788 and can be interpreted as an estimate of the standard deviation in log times to failure. The intercept is µ = 5.6768 and can be interpreted as an estimate for the mean of the log times to failure at the stated test conditions. The survivor function associated with his normal distribution (using these parameter estimate) is shown in Figure 2. It tends to be lie above the parametric estimators at intermediate failure times, but below it as the higher failure times.

A Statistical Description of Continuous Failure Times at Varying Test Conditions
There are a number of approaches to extending the above concepts to the case of varying test conditions.

Accelerated Failure Time Models (AFT)
In this type of model, μ in Equation (9b) is made a function of the test conditions where x1 to xm are separate variables describing the test condition (for example, x1 may be stress, x2 temperature etc.) and r is an un-specified function (its form being best suggested by creep theory). x is a 1 by m matrix containing the m test variables that describe the test conditions for each of the N specimen placed on test. A commonly used specification for r(x) is where b1 to bm are parameters that require estimation. As the name suggests, this approach has an accelerated life interpretation. In this formulation, the error term σW is seen as a base or reference distribution that applies when x1 = x2 = … = xm = 0. This base distribution can be translated to a time scale by defining T0 = exp{σW}. The probability that a test specimen will survive time t, S0(t), is then In this accelerated model, T is distributed as and so the test conditions act multiplicatively on survival times. Therefore, the probability that a test specimen with test conditions x will be survive time t is

A Statistical Description of Continuous Failure Times at Varying Test Conditions
There are a number of approaches to extending the above concepts to the case of varying test conditions.

Accelerated Failure Time Models (AFT)
In this type of model, µ in Equation (9b) is made a function of the test conditions where x 1 to x m are separate variables describing the test condition (for example, x 1 may be stress, x 2 temperature etc.) and r is an un-specified function (its form being best suggested by creep theory). x is a 1 by m matrix containing the m test variables that describe the test conditions for each of the N specimen placed on test. A commonly used specification for r(x) is where b 1 to b m are parameters that require estimation. As the name suggests, this approach has an accelerated life interpretation. In this formulation, the error term σW is seen as a base or reference distribution that applies when x 1 = x 2 = . . . = x m = 0. This base distribution can be translated to a time scale by defining T 0 = exp{σW}. The probability that a test specimen will survive time t, S 0 (t), is then S 0 (t) = Pr{T 0 > t} = Pr{W > ln(t)/σ} In this accelerated model, T is distributed as and so the test conditions act multiplicatively on survival times. Therefore, the probability that a test specimen with test conditions x will be survive time t is S(t,x) = Pr{T > t|x} = Pr{T 0 e r(x) > t} = Pr{T 0 > te r(x) } = S 0 (te b1x1+ . . . +bmxm ) (12c) Thus, the probability that a specimen with test conditions x will survive time t is the same as the probability that a base test specimen will be alive at time texp{r(x)}. This can be interpreted as time passing more rapidly by a factor exp{r(x)}-for example, twice as fast or half as fast. (A good analogy here is the use by humans of pet years to describe the age of their pets in relation to their life). Consider for example a multiplier of two for a specimen with test condition x. In terms of survival, this means that the probability that the specimen would be alive at any given time is the same as the probability that a base specimen would be alive at twice the length of time. In terms of risk, this model implies that an engineering component is exposed at any service life to double the risk of a base component that has been in service for twice as long.
The importance of Equation (12c) for this paper is that Evans [32] has shown, when using an AFT model, that whilst a generalised F distribution explained the shape of the failure time distribution at most test conditions for 1Cr-1Mo-0.25V steel, none of the distributions contained as special cases within the generalised F distribution adequately explained the shape of the actual failure time distributions at the remaining test conditions. This failure is explained by Equation (6c) as it shows that the survivor function should have the same form at all test conditions (namely that form identified for specimens tested at the base conditions)-unless time is stretched too much. However, in hazard based models, to be discussed below, the survivor function at a particular test conditions can differ markedly from that identified at the base test conditions and so offers extra flexibility over AFT models.

Proportional Odds Models
Another approach assumes that the effect of the test conditions is to increase or decrease the odds of failure by a given duration by a proportionate amount: where S 0 (t,x) is a baseline survivor function, taken from a suitable distribution, and exp{b 1 x 1 + . . . + b m x m } is a multiplier reflecting the proportionate increase in the odds associated with test condition values x. Taking natural logs, gives so the test conditions effects are linear in the logit scale. A somewhat more general version of the proportional odds model is known as the relational logit model. The idea is to allow the log-odds of failing in a given population to be a linear function of the log-odds in a reference or baseline population, so that logit(1 − S(t)) = α + θlogit(1 − S 0 (t)) (13c) The proportional odds model is the special case where θ = 1 (and where the constant α depends on the test conditions).
As an example, consider a proportional odds model with a log-logistic baseline. The corresponding survival function and the odds of failure are Multiplying the odds by exp(b 1 x 1 + . . . + b m x m ) yields another log-logistic model. However, this is not true of other distributions: if the baseline survivor function is Weibull then this baseline multiplied by the odds of failing is not a Weibull survivor function.

Proportional Hazard Models (PH)
The PH model of Cox [39] has a baseline hazard function h 0 (t) that shows how the hazard rate increases with time when this linear combination of test conditions equals unity The log hazard function is then additive Obviously, the cumulative hazards would follow the same relationship, as can be seen by integrating both sides of the previous equation. Exponentiating minus the integrated hazard, we find the survivor functions to be S(t,x) = S 0 (t) exp(b1x1+b2x2+ . . . +bmxm) (14c) so the survivor function for test conditions x is the baseline survivor raised to a power that is dependent upon the test condition. If a test specimen is exposed to twice the risk of a reference specimen at every point in time, then the probability that the specimen will be alive at any given time is the square of the probability that the reference or base specimen would be alive at the same time. In this PH model, a simple relationship in terms of hazards translates into a more complex relationship in terms of survival functions. Choosing a different parametric form for the baseline hazard, leads to a different model in the proportional hazards family. Apart from when the baseline hazard function corresponds to that of the Weibull hazard function, the hazard function at all other test conditions will be different in form from the baseline hazard function. A possible limitation of the PH model is seen in Equation (14a), which implies that hazard functions associated with different test conditions are always constant multiples of one another-hence the name "proportional" hazards. One way to relax this proportionality assumption is to allow the test variables to interact with time or equivalently to allow b 1 , b 2 etc. to be time dependent. Then, If, for example, the base line hazard function corresponds to the log normal distribution (so under test conditions r(x) = 1 the underlying failure time distribution is log normal), the underlying failure time distribution will not be log normal at any other test condition (i.e., when r(x) = 1). This is a major advantage of building a failure time model around the hazard function rather than around the pdf or f(t).

Modelling Discrete Failure Times at Varying Test Conditions
Another major issue with hazard based models is to do with the identification of the baseline hazard function, h o (t). Without having many repeat tests carried out at a single test condition, it is difficult to accurately identify its functional form. One solution to this problem is to create a discrete failure time dataset from the original continuous one, i.e., split the continuous failure time data up into small but equally sized time spans. By doing so, it is possible to calculate a piecewise hazard function for each interval of time, which over all time intervals allows the shape of the base line hazard function to be identified. Springer and Willett [40] provide a good review of this approach. This is the approach taken in Section 4.

Creating Discrete Data from Continuous Data
The first step required in building a discrete hazard function is to create the specimen-specimens dataset from the continuous failure time database. Here time is partitioned into k equal intervals I j = (a j−1 to a j ), j = 1 to k and with k being as large as practically possible and a being a point in time.
As an illustration of how this is done, consider batch VaA of the NIMS database, where i = 1, N = 43 specimens are tested, with each specimen receiving a different stress-temperature test combination. If x 1 in Equation (12a) represent stress, then in batch VaA of NIMS this was varied from 412 MPa to 47 MPa and if x 2 represents temperature this varied over the range 723 K to 948 K. The smallest recorded failure time was 338,760 s and the largest was 407,844,720 s. Many creep prediction models work with the log time to failure and so in natural log units these failure time limits corresponded to 12.73 and 19.83. The researcher then needs to decide upon how many discretised time intervals to work with. For example, the NIMS data could discretised into 15 equal (log) time intervals, respectively giving (log) time intervals of width 0.5. In this example, k = 15 and a j−1 − a j = 0.5 with a 0 = 12.5 and a k = 20). The first interval this NIMS dataset is therefore 12.5-13.0 and corresponds to j = 1, the second is 13.0-13.5 and corresponding to j = 2 all the way up to the interval 19.5-20 and corresponding to j = 15. The specimens-specimen data are then generated by creating a binary variable, v, for each time interval. Thus, the binary variable equals 0 in time interval I t if the specimen does not fail in that interval and 1 if it does. This binary variable is created for each specimen in the test matrix. Table 1 illustrates the start of the creation of this specimens-specimen format by considering just the firsts two NIMS specimen in batch VaA using a j−1 − a j = 0.5. The first specimen was tested at x 1 = 412 MPa and x 2 = 723 K. It failed at 16.36 log seconds. The second was tested at x 1 = 373 MPa and x 2 = 723 K and it failed at 17.84 log seconds. Continuing the process shown in Table 2, creates M values for v, where M = kN = 358 for this NIMS batch.
In Equation (14b), x is a Nk by 2 matrix where each column contains the i different stress and temperature combinations that each specimen was tested at and x i is the ith first row of the matrix x. However, when the data are discretised in this way, h(ij|x) in Equation (14a) is not directly observable. Instead, there is the binary variable v ij that equals zero when the specimen is un-failed in time interval a j−1 − a j or 1 if it fails in that time interval. Therefore, what is required is a non-linear function that maps between v ij = 0 and v ij = 1 to give the hazard rate between 0 and 1 for each specimen in each time interval. Two commonly used functions that achieve this are the logistic so Equation (16d) collapses to Equation (16b) when α = 1 and tends to Equation (16c) as α tends to ∞.
Plotting estimated values of b 0,j against time allows the shape of the baseline hazard function to be observed, which in turn could lead to a simplification of Equation (16d). For example, if such a plot reveals a straight line, the k b o,j D j terms can be replaced by a Non-proportionality can then be accommodated by allowing the x variables to interact with time t in the following way

Estimation
These discrete hazard models are essentially binary response regression models and so the unknown parameters can be estimated by maximising the log likelihood function given by where the binary response variable v ij and h(t|x) is given by Equation (16d), Equation (16e) or Equation (16f). Hence, Equation (17) will be maximised yielding some given set of values for the unknown parameters on the right hand side of Equation (16d), Equation (16e) or Equation (16f). This function can be maximised for various values of α, with the chosen value for α corresponding to the largest of these maximised log likelihoods. Direct maximisation can be carried out using standard non-linear optimisation algorithms [41], or alternatively Equation (17) can be maximised indirectly by using the iteratively reweighted least squares algorithm of McCullagh and Nelder [42].

Assessing Model Adequacy
There are a number of ways to assess the adequacy of these discrete hazard based models. The discrete hazard models given by either Equation (16d), Equation (16e) or Equation (16f) produce a predicted hazard rate or probability of failure within the time intervals a t−1 to a t , rather than a specific time at which failure occurs. This makes it a little more complicated to assess whether the model is capable of accurately predicting the observed failure times. To obtain predicted times to failure (more specifically the time interval in which failure is predicted to occur), a cut-off point c is needed, such that if the predicted hazard rate exceeds c in time interval a t−1 to a t then failure is predicted to have occurred in that time zone. Otherwise, the prediction is that the specimen will not fail in that time interval. The usual value for c is 0.5, and once chosen a simple classification table, such as that in Table 2, can be constructed.
In Table 2, there are M 2 time intervals where a specimen failed. The model correctly predicted correctly m 4 of these, but incorrectly predicted m 3 of these intervals, that is m 3 specimens failed in times zones different from those the model predicted them to fail in (where M 2 = m 3 + m 4 ). Similarly, there are M 1 time zones where v ij actually equalled zero (M 1 time zones not containing failed specimens) and of these, the model predicted correctly m 1 of these time zones, but incorrectly predicted m 2 of these zones (i.e., there were m 2 times zones where a specimen survived, but the model predicted failures to occur). M can be found from summing either M 1 and M 2 or M 3 and M 4 .
The models success rate in predicting the time zones where specimens will not fail is given by m 1 /M 1 . This can be taken to be the probability of the model correctly detecting a false signal and is called the models specificity. The models success rate in predicting the time zones where specimens will fail is given by m 4 /M 2 . This can be taken to be the probability of the model correctly detecting a true signal and is called the models sensitivity. The models over-all rate of correct classification is then given by (m 1 + m 4 )/M. Given the way in which a specimens-specimen dataset is constructed, there are many more values of v ij = 0 compared to v ij = 1, it is common to observe sensitivity values well below specificity values. Thus a well specified model will have high values for both sensitivity and specificity. However, the sensitivity and specificity depend in part on the chosen value for the cut-off point c and c = 0.5 may not be the optimal value as far as failure time prediction is concerned. One solution to this is to construct a classification table for a range of cut-off points to see how the discrete hazard model works as a classifier of when failure will occur. However, a neater way of doing this is given by the area under the receiver operating characteristic (ROC) curve. The ROC emerges on a graph that plots the sensitivity against (1-specificity) associated with all possible values for c (c = 0 to 1). The resulting area under the ROC curve lies between zero and unity and measures the ability of the hazard model to discriminate between time zones where a specimen will fail and zones where it will not. Hosmer and Lemenshow [43] suggest the following benchmarks for this area: An ROC of 0.5 provides no discrimination implying the model performance no better than tossing a coin to decide if a specimen fails in a given time interval. An ROC between 0.7 and 0.8 gives acceptable discrimination, whilst a ROC between 0.8 and 0.9 gives excellent discrimination. Finally, a ROC above 0.9 gives outstanding discrimination.
One way to determine an optimal value for c is to choose that value that yields the best failure time predictions. For example, if c inserted for v in say Equation (16e), then the variable t on the right hand side of the equation becomes the predicted time zone at which a specimen fails,t This can be solved forĵ Thus,ĵ is a prediction of the time interval where a specimen fails. For example, ifĵ = 2.5, failure is predicted to occur in the time interval I 2 = a 2 − a 1 , and the actual predicted time is then equal to (a 2 + a 1 )/2 (or another example, actual failure time = 0.25a 2 + 0.75a 1 ifĵ = 2.75). Notice the j = 2 time interval is, from Table 1, the 12.5-13 logged seconds and its mid point is therefore 12.75 logged seconds (or 96 h).
A plot can then be made of actual failure time against this predicted time. If a best fit line is put through the data on such a plot, the optimal value for c can be taken to be the one that produces a best fit line closest to the 45 • line on such a plot. As an alternative, c can be chosen to minimise the mean squared error defined as Using c = 0.5 in Equation (18b) can also be interpreted as yielding a median predicted time to failure, whilst using c = 0.05 produces a time interval prediction such that there is only a 5% chance of failure occurring in that or an early time interval. Likewise using c = 0.95 produces a time interval prediction such that there is a 95% chance of failure occurring in that or an early time interval. These then come together to define a 90% confidence interval for the time interval where failure will occur.

Application of Discrete Hazard Function to Batch
where x 1 = τ* and x 2 = 1000/RT. However, h(ij|x) is not observable. Instead there is the binary variable v ij that equals zero when the specimen is un-failed in time interval a j−1 − a j or 1 if it fails in that time interval. Therefore, what is required is a non-linear function that maps between v ij = 0 and v ij = 1 to give the hazard rate between 0 and 1 for each specimen in each time interval. Section 3.3.2 outlined a general form of such a function allowing Equation (19a) to be written as If the b 0,j parameters trace out a linear time trend, then Equation (19b) can be written as As will be seen below, this plausibility of such a simplification can be assessed by plotting out the estimated b 0j values. Non-proportionality can then be accommodated by allowing τ and 1000/RT to interact with time in the following way In addition, it is known (see Wilshire [21] and Evans [32]) that for this material the relationship between t and τ* changes at some critical value for τ* (i.e., b 1 changes value at this point) and that the activation energy changes at around 823 K (i.e., b 2 changes value at this point). Thus, Equation (19b) can be written as where B 1 = 0 when τ* ≤ τ crit and unity otherwise. Similarly, B 2 = 0 when T ≤ 823 and unity otherwise. The reason for doing this is that it allows stress and temperature to have a different effect on the hazard rate either side if 823 K and τ crit . The explanation provided by Wilshire is that as τ crit is close in value to the yield stress, dislocation movement is confined to grain boundaries below the yield stress so that the activation energy falls below that for self diffusion through the crystal-so causing the value for b 2 to change. Thus, when τ* ≤ τ crit , B 2 = 0, and so the effect of temperature on the hazard rate is determined by the value for b 2 . However, when τ* > τ crit , B 2 = 1, and so the effect of temperature on the hazard rate is determined by the value for b 2 + b 4. The role of stress is also different either side of τ crit (changing from b 2 to b 2 + b 3 ). Often, the values for b 0,j reveal a linear trend or some well defined non-linear trend such as a polynomial, exponential or power law trend. For example, if a linear trend is revealed by a plot of the b 0,t values against t, then Equation (19e) takes the form In Equation (19e,f), the effect of changing test conditions is to shift in a parallel fashion the log baseline hazard function, but it is also possible to allow the slope of the baseline hazard function to depend on x 1 and/or x 2 . For example, if b 0 depends on x 1 then Equation

Model Given by Equation (19e)
A specimens-specimen dataset was created for batch VaA using k = 15 and I j = a j−1 − a j = 0.5 with a 0 = 12.5 and a k = 20. This dataset consisted of M = 358 observations on v. Using these data, the parameters of Equation (19e) were estimated for a range of values for α and τ crit . The values for α and τ crit that produced the highest ROC were α = 1 and τ crit = 0.1. This α value suggests the logistic discrete hazard model is preferable to the log-log discrete hazard model, whilst the τ crit value is a little higher than that identified by Evans [32] and Wilshire [21] but is broadly similar in value. Table 3 shows the results of applying the McCullagh and Nelder algorithm to this data. The values for x 1 and x 2 were normalised to be zero at 823 K and 294 MPa so that the resulting estimated value for b 0,t give the baseline log hazard rates that corresponds to this test condition. x 4 642.0808 −4.43 *** Parameters estimates using the iteratively reweighted least squares technique of McCullagh and Nelder [42]. Student t-values test the null hypothesis that the true parameter values equal zero. ** identifies statistically significant variables at the 5% significance level, and *** identifies statistically significant variables at the 1% significance level. These levels of significance are based on the student t-statistic that has a student t distribution.
The student t-values associated with x 1 and x 2 in Table 3 reveals that both τ* and the reciprocal of temperature are statistically significant variables so yielding support for the Wilshire methodology, whilst the last two rows show a discontinuity in the Wilshire model at τ* = 0.1 and at a temperature of 823 K. These estimates are not comparable in value to those in Equation (5b) because the latter show the impact of τ* and T on failure times directly, rather than the log hazard rate as is so in the former case. Reading across row one of Table 3, the estimated log hazard rate for time interval j = 1 (12.5-13.0 log seconds) when temperature is 823 K and stress is 294 MPa is −4.7526, implying a hazard rate of exp(−4.7526) = 0.0086. This hazard rate corresponds to the log time interval of 12.5 to 13.0, which in turn corresponds to a time interval in hours of 75 to 123. Associating this hazard rate with the mid point of this time interval gives a hazard rate of 0.0086 at time 100 h. Proceeding in the same way for this next two rows of Table 3 gives a cumulative hazard rate of 0.189 at 163 h and 0.795 at 268 h. Recall that Figure 2 show the non-parametric estimate of the cumulative hazard rate associated with the test conditions of this estimated baseline hazard rate. Reading of this graph at 270 h shows the non-parametric estimate of the cumulative hazard rate is 0.85, showing this model is consistent with this non-parametric estimator. Table 4 shows the classification table for this model using c = 0.5 revealing that sensitivity equals 68% and specificity equals 98% with an overall correct failure time zone prediction rate of 94%. Figure 5 shows the ROC for this model and the area under this curve is 0.972, which according to Hosmer and Lemenshow, makes this model outstanding in its ability to discriminate between failure and non-failure in the 15 time zones.   The MSE is minimised at c = 0.53 with a value of 0.0526. Figure 6 then plots the actual times to failure (hours) for specimens in batch VaA, together with the failure time zones predicted by the model using c = 0.53. The error bars shown around the predictions reflect the fact that this model predicts the time interval in which failure occurs and the width of this time interval is 0.5 log hours. The models predictions are taken to be the mid points of these error bars. With only a few exceptions (for example, at 773 K with 373 MPa and 823 K with 157 MPa), the model predicts the time interval at which failure actually occurs correctly. The worst prediction comes at 873 K and 47 MPa-but this point was also poorly predicted in the original Wilshire [21] paper as well. Figure 7 plots the piece-wise log hazard rates associated with each time zone (the b 0t in Table 3) against the numbered time zone and this plot reveals a well-defined linear trend. The trend is very strong with an R 2 value of 98.4% with no obvious deviation from linearity. This suggests that it should be possible to create a more parsimonious version of Equation (19e), without affecting the predictive ability of the simpler model, by replacing the fifteen b 0t parameters with a linear trend containing just two parameters-a and b 0 . Figure 7 implies that a = 6.57 and b 0 = 1.86. This parsimonious version is estimated in the following subsection. Figure 5. The receiver operating characteristic (ROC) from model given by Equation (13e).

Model Given by Equation (19f)
When Equation (19g) was estimated, the t value associated with parameter b5 suggested that b5 was insignificantly different from zero at the 1% significance level so that the data were not supportive of the slope of the log base hard function estimated in the previous sub section changing at some critical value for the normalised stress. Table 5 shows the estimated parameters of Equation (19f) obtained using the McCullagh and Nelder algorithm (again the values for α and τcrit that produced the highest ROC were α = 1 and τcrit = 0.1).

Model Given by Equation (19f)
When Equation (19g) was estimated, the t value associated with parameter b 5 suggested that b 5 was insignificantly different from zero at the 1% significance level so that the data were not supportive of the slope of the log base hard function estimated in the previous sub section changing at some critical value for the normalised stress. Table 5 shows the estimated parameters of Equation (19f) obtained using the McCullagh and Nelder algorithm (again the values for α and τ crit that produced the highest ROC were α = 1 and τ crit = 0.1).  [42]. Student t-values test the null hypothesis that the true parameter values equal zero. *** identifies statistically significant variables at the 1% significance level. These levels of significance are based on the student t-statistic that has a student t distribution.
The estimated values for b 1 to b 4 in Table 5 are consistent with those shown in Table 3 and, again, the student t-values associated with these parameters reveal they are significantly different from zero at either the 1% or 5% significance level. Further, the values for a and b 0 are not that dissimilar from the values sown in Figure 6. These parameter estimates can be used to calculate the hazard rate at the base or reference test conditions of temperature = 823 K and stress = 294 MPa. For example, in time zone j = 1 (corresponding to the log time interval of 12.5 to 13 or 75 to 123 h) the hazard rate is predicted to be exp(−5.8351 + 2.0736 × 1) = 0.0232. Table 6 shows the classification table for this model when c = 0.5 revealing that sensitivity equals 69.8% and (1-specificity) equals 96.8% with an overall correct failure time zone prediction rate of 95%. Figure 8 shows the ROC for this model and the area under this curve is 0.951, which according to Hosmer and Lemenshow, makes this model outstanding in its ability to discriminate between failure and non-failure in the 15 time zones. Further, this value is only slightly below that from the previous model, showing that the use of a simple linear time trend instead of a piece-wise hazard function has not resulted in a significant reduction in predictive ability.  Table 6 shows the classification table for this model when c = 0.5 revealing that sensitivity equals 69.8% and (1-specificity) equals 96.8% with an overall correct failure time zone prediction rate of 95%. Figure 8 shows the ROC for this model and the area under this curve is 0.951, which according to Hosmer and Lemenshow, makes this model outstanding in its ability to discriminate between failure and non-failure in the 15 time zones. Further, this value is only slightly below that from the previous model, showing that the use of a simple linear time trend instead of a piece-wise hazard function has not resulted in a significant reduction in predictive ability.  The MSE is minimised at c = 0.53 and Figure 9 then plots the actual times to failure for specimens in batch VaA against the failure time predicted by the model using c = 0.53. As a continuous base hazard function is now used instead of the piece-wise hazard function, an actual failure time, rather than an interval time, prediction can be made. Figure 9 plots the actual failure times against the predicted failure times (n natural logs). The best fit line on this plot is very close to the ideal outcome associated with the 45 degree line-corresponding to a situation where the model predicts each  The MSE is minimised at c = 0.53 and Figure 9 then plots the actual times to failure for specimens in batch VaA against the failure time predicted by the model using c = 0.53. As a continuous base hazard function is now used instead of the piece-wise hazard function, an actual failure time, rather than an interval time, prediction can be made. Figure 9 plots the actual failure times against the predicted failure times (n natural logs). The best fit line on this plot is very close to the ideal outcome associated with the 45 degree line-corresponding to a situation where the model predicts each failure time perfectly. Hence, the bias in prediction remains small in this more parsimonious model.  Figure 10 shows a different representation of these predictions-where stress is shown on the vertical axis and times to failure on the horizontal. This time the error bars show a 50% prediction band based on using c = 0.25 and c = 0.75 in Equation (18b). Again, and with only a few exceptions, the actual failure times fall within the models 50% prediction bands. Figure 11 illustrates how this model can be used to predict the hazard rates associated with in service life under various conditions-in this illustration operating at 130 MPa and 823 K. It can be seen that the hazard rate remains very close to zero up to 50,000 h or about 8 years of in service use. The risk of failure then starts to rise quite dramatically. For example, if this material had been in operation for around 15 years, the chances of it failing in the next year is around 35%, but, if this material had been in operation for around 25 years, the chances of it failing in the next three years rises dramatically to around 80%.

Conclusions
This paper has provided a summary review of some statistical failure time models with the aim of aiding the assimilation of such models with existing predictive models for creep life. This will enable an enrichment of prediction to be achieved with a move away from predicting failure times on the average towards predicting the safe life associated with a minimum chance of failure This was followed by an illustration of one possible assimilation, namely-the deterministic Wilshire equation and the statistical discrete hazard model. This statistical model was chosen because it provided the capability of estimating failure probabilities in future time intervals for materials that have already been in service for various lengths of time.
Estimation of this model revealed that at a fixed test condition, the log of the probability of failure in the next time interval (given survival up to then) is a linear function of time. This log base

Conclusions
This paper has provided a summary review of some statistical failure time models with the aim of aiding the assimilation of such models with existing predictive models for creep life. This will enable an enrichment of prediction to be achieved with a move away from predicting failure times on the average towards predicting the safe life associated with a minimum chance of failure This was followed by an illustration of one possible assimilation, namely-the deterministic Wilshire equation and the statistical discrete hazard model. This statistical model was chosen because it provided the capability of estimating failure probabilities in future time intervals for materials that have already been in service for various lengths of time.
Estimation of this model revealed that at a fixed test condition, the log of the probability of failure in the next time interval (given survival up to then) is a linear function of time. This log base line hazard function then shifted in a parallel fashion with the well-known Wilshire variables of τ* and the reciprocal of temperature. Like in the original Wilshire methodology, this shifting nature of the base hazard function was different above and below 823 K and τ* = 0.1. The model was shown to produce outstanding discrimination with respect to which time interval a specimen would fail in. Finally, and as an illustration of the output this model was capable of producing, it was found that if this material had been in operation for around 15 years at 823 K and 130 MPa, the chances of it failing in the next year is around 35%. However, if this material had been in operation for around 25 years under this condition, the chances of it failing in the next year rose dramatically to around 80%.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. The Generalised F Family of Failure Time Distributions
Using the change of variables technique, and substituting Z = (Y − µ)/b into the pdf of Equation (9c) in Section 3.1.2 yields f(y) = f(z) ∂z ∂y (k 1 /k 2 ) k 1 {Γ(k 1 )Γ(k 2 )}/Γ(k 1 +k 2 ) (λt) k 1 β 1+{k 1 /k 2 }(λt) k 1 β (k 1 +k 2 ) = (k 1 /k 2 ) k 1 λ β {Γ(k 1 )Γ(k 2 )}/Γ(k 1 +k 2 ) (λt) β k 1 −1 This equation is not valid when k 1 = k 2 = ∞, which corresponds to the log normal distribution. Abramowitz and Stegun [44] show that the gamma function can be approximated by ln[Γ(k)] = (k − 0.5) ln(k) − k + 0.5 ln(2π) and based on this approximation, it follows that Г (k + 1) = kГ(k) and Г(1) = 1. As mentioned in the main text, Z is a standardised variable, and its mean (or expected value E) and its variance are given by where Var(Z) reads the variance of Z and the di (Ψ) and tri (Ψ / ) gamma functions are the first and second derivatives respectively of the log gamma function with respect to k and so The distribution of Z is degenerate in that as either k 1 or k 2 tends to infinity, Ψ / (k) tends to zero and so the variance of Z tends also to zero, i.e., the pdf for Z collapses to a single point. However, from the rules behind computing expected values, it follows that Var(W) = δ 2 Var(Z).
It follows from Equations (A7) and (A8) and the rules of expected values that: and Var(Y) = b 2 Var(Z) (A10) where Var(Y) is the variance of Y. Meeker and Escobar [45] have shown that, when k 2 > 1/β the mean and variance for T are given by There is no closed form expression for the survivor function, but the pth percentile for T is given by t p = exp(µ)(w k1,k2,p ) (b/δ) (A12) where w k1,k2,p is the pth quantile of an F distribution with (2k 1 , 2k 2 ) degrees of freedom. t p can thus be computed for various values of p and a plot of 1 − p against t p defines the survivor function for T.
The survivor function is therefore given by Substituting k 1 = 1 into Equation (A11) leads to the following expressions for the mean and variance of T Prentice [37] has shown that, when k 2 = ∞, Equation (9c) in Section 3.1.2 reduces to with δ = √ k 1 and thus W = Using the change of variable technique, it follows that f(y) = f(z) ∂z ∂y Again, using β = 1/b, λ = exp(−µ) and T = exp(Y), and the change of variable technique As k 2 tends to infinity it can be shown that Γ(k 2 − σ) Γ(k 2 ) (k 2 ) σ → 1 and Γ(k 2 − 2σ) Γ(k 2 ) (k 2 ) 2σ → 1 (A26) and substituting this expression and k 1 = 1 into Equation (A11) gives the formulas for the mean and variance of T when this variable is Weibull distributed Var[T] = e 2µ Γ(k 1 +2σ)Γ(k 2 −2σ) If σ also equals 1, then β = 1 and the Weibull distribution collapses to the exponential distribution σ = 1 implies β = 1/ √ k 1 and so