Case Study: Mechanical Switches
Data
Table 1 shows failure times8 in Millions of Operations of 40 randomly selected mechanical switches tested in a facility. 3 switches had not failed by the end of the test, leading to right censored observations. This dataset is from an article by Nair (1984)9 Nair, V. N. (1984). Confidence bands for survival functions with censored data: A compar- ative study. Technometrics 26, 265–275.. It is also publicly available on Datashare, Iowa State University’s open data repository and can be accessed through this link.
Table 1. Mechanical Switches Life Test Data
| Millions of Operations | Failure Mode |
|---|---|
| 1.151 | Spring B |
| 1.170 | Spring B |
| 1.248 | Spring B |
| 1.331 | Spring B |
| 1.381 | Spring B |
| 1.499 | Spring A |
| 1.508 | Spring B |
| 1.534 | Spring B |
| 1.577 | Spring B |
| 1.584 | Spring B |
| 1.667 | Spring A |
| 1.695 | Spring A |
| 1.710 | Spring A |
| 1.955 | Spring B |
| 1.965 | Spring A |
| 2.012 | Spring B |
| 2.051 | Spring B |
| 2.076 | Spring B |
| 2.109 | Spring A |
| 2.883 | Censored |
| 2.883 | Censored |
| 3.793 | Censored |
| 2.116 | Spring B |
| 2.119 | Spring B |
| 2.135 | Spring A |
| 2.197 | Spring A |
| 2.199 | Spring B |
| 2.227 | Spring A |
| 2.250 | Spring B |
| 2.254 | Spring A |
| 2.261 | Spring B |
| 2.349 | Spring B |
| 2.369 | Spring A |
| 2.547 | Spring A |
| 2.548 | Spring A |
| 2.738 | Spring B |
| 2.794 | Spring A |
| 2.910 | Spring A |
| 3.015 | Spring A |
| 3.017 | Spring A |
There are 2 modes of failure - Spring A and Spring B10 The 2 springs are identical in construction, but are subject to different stresses, thereby leading to different time-to-failure probability distributions.. We are interested in constructing a prediction interval for the time to failure of a mechanical switch, irrespective of mode of failure. Note that a proper analysis of this dataset should treat the two modes of failures as competing risks and apply appropriate methods that are capable of handling competing risks. We are not going to do that in this article. The purpose of this blog article is to discuss prediction intervals and their calibration. A simple parametric model ignoring competing risks will be assumed in the interest of a succinct article focusing on one topic i.e prediction intervals.
Model
The log location scale distribution model that we consider for this dataset is shown in Eq. A below.
\[ \hspace{-2cm} \begin{align*} F(t;\mu,\sigma) & = Pr(T\le t)\\ & = \Phi \left[\frac{log(t)-\mu}{\sigma}\right] \tag{Eq. A} \end{align*} \]
The distribution of \(\Phi\) determines the time-to-failure distribution. For example, if \(t\) is assumed to have a Weibull distribution, \(\Phi\) is the cumulative distribution function (CDF) of the standard smallest extreme value (SEV) distribution. Similarly, if \(t\) is assumed to have a log-normal distribution, then \(\Phi\) is the CDF of the standard normal distribution.
Weibull and log-normal distributions are widely used to describe time-to-failure distributions in reliability applications. Other distributions11 Eg: exponential, log-logistic, generalized gamma, normal etc can also be considered, depending on the application and motivated by the failure process. Only log-normal and Weibull distributions are considered as candidate models in this article.
Figure 1 shows a Weibull probability plot and Figure 2 shows a log-normal probability plot of the mechanical switches dataset.
Figure 1: Weibull probability plot
Figure 2: log-normal probability plot
The log-normal model looks like a better fit to the data as compared to the Weibull model, especially at the lower tail. The AIC of the Weibull model is 83.01 and the AIC of the log-normal model is 77.38, further reinforcing that the log-normal is a better choice for this data. The maximum likelihood (ML) estimates of the parameters of the log-normal distribution and their standard errors are shown in the table 2 below.
Table 2. Maximum Likelihood estimates of the log-normal model
| Estimate | |
|---|---|
| \(\hat{\mu}\) | 0.723 |
| \(\mathrm SE_{\hat{\mu}}\) | 0.048 |
| \(\hat{\sigma}\) | 0.300 |
| \(\mathrm SE_{\hat{\sigma}}\) | 0.036 |
Plug-in Prediction Interval
Now that we have the ML estimates of the log-normal distribution parameters, we can easily obtain the 90% approximate12 As opposed to an “exact” prediction interval where the confidence level is exactly met. prediction interval. A one sided lower 95% approximate prediction bound and a one sided upper 95% approximate prediction bound are combined to form an equal tailed two sided 90% approximate prediction interval. Although it might be possible to find a narrower interval with unequal upper and lower tail probabilities that add to \(\alpha\) = 10%, two-sided intervals are often reported in reliability applications even when the primary interest is on one side or the other. With the appropriate adjustment to the confidence level, a two-sided interval can be correctly interpreted as a one-sided prediction bound. For the mechanical switches, a two sided approximate 90% prediction interval calculated using the ML estimates of the parameters is [1.26,3.38].
Calibrating Prediction Intervals
The method of finding the approximate prediction interval described in the previous section is a “naive” method because it simply uses the ML estimates of the parameters. It ignores the uncertainty in the estimated parameter values, which then leads to the coverage probability13 Coverage probability is the probability that a prediction interval procedure gives an interval containing the future prediction. being generally smaller than the nominal confidence level. When there is a lot of information in the data set i.e large number of failures and little censoring, the naive plug-in method may suffice. If this is not the case, it is necessary to refine/calibrate the interval to ensure better coverage probability.
Escobar and Meeker (1999)14 Escobar, L. A. and W. Q. Meeker (1999). Statistical prediction based on censored life data. Technometrics 41, 113–124., discuss different prediction interval calibration methods in their paper. Here, we use the method where calibration is done by simulating B15 B is typically chosen between 10,000 and 50,000. realizations of the data and averaging the conditional coverage probabilities. A comprehensive discussion on this method is outside the scope of this blog article and I encourage you to read the paper or the textbook that I have mentioned previously. A summary of the steps involved in this method is shown below:
Step 1: Choose a particular value of the one sided confidence level \(1 - \alpha_{c}\). In our case study, the one sided confidence level is 5%. So, we can start with \(\alpha_{c} = 0.05\).
Step 2: Simulate a realization of the available data using fractional random weight bootstrap16 Note on fractional random weight bootstrap method is included in the next section. sampling.
Step 3: Find the maximum likelihood estimates of the distribution parameter for the data simulated in Step 2.
Step 4: Find the one sided upper 95% prediction interval from the parameter estimates in step 3.
Step 5: Find the coverage probability of the interval from Step 4 using the original ML estimates of the parameters in Table 2.
Step 6: Repeat Steps 2,3,4 and 5 \(B\) times. \(B = 10,000\) has been used in this case.
Step 7: Take the average of the 10,000 coverage probabilities from Step 5. If this value is less than our target confidence level of 95%, decrease the value of \(\alpha_{c}\) in Step 1 and repeat the whole procedure until the average in Step 7 is at or above our target confidence level.
The above procedure will give us a calibrated one sided approximate 95% upper prediction bound. Repeat the same procedure to find the calibrated one sided approximate 95% lower prediction bound by modifying Step 4 to calculate the lower bound instead of the upper bound.
Combining the two one-sided 95% bounds gives a two-sided 90% interval. Table 3 below shows the calibrated interval along with the naive plug in interval.
Table 3: Approximate 90% Prediction Intervals
| Naive 90% Prediction Interval | Calibrated 90% Prediction Interval |
|---|---|
| [1.26, 3.38] | [1.22, 3.49] |
Notice that the width of the calibrated interval is greater than the width of the naive interval. The calibrated approximate prediction interval is expected to have better coverage properties than the naive approximate prediction interval.
Note on Fractional Random Weight Bootstrap
With a traditional bootstrap approach, we randomly sample data (from Table 1), with replacement. In a given bootstrap sample, some rows from table 1 may appear more than once while other may not appear at all. Think of the number of times that a particular row appears in a bootstrap sample as the “weight” for that observation. These weights are integer values and can take on values between 0 and \(n\), where \(n\) is the number of data points in a bootstrap sample. Note that these weights are positive “integer” values. When there is heavy censoring in the data, a bootstrap sample might consist of only censored observations, or very few failures. This leads to estimation problems when we try to estimate parameters using maximum likelihood. An alternative to these “integer weights” is “fractional weights”, where each row in Table 1 gets assigned random non-integer weights. These non-integer weights can be generated from a continuous distribution of a random variable that has the same mean and standard deviation as the integer resampling weights. Xu et al. (2020)17 Xu, L., C. Gotwalt, Y. Hong, C. B. King, and W. Q. Meeker (2020). Applications of the fractional‐random‐weight bootstrap. The American Statistician 74, 345–358. discusses an application of this method and shows how to generate the random weights. Although we don’t have heavy censoring in our case and would probably not have encountered any estimation problems, we have still used fractional random weight bootstrapping to be on the safer side.