Using remaining useful life predictions to assess model fit with right censored data

Question

I am modelling time to failure of some units. All units are made my the same manufacturer and there are no recorded covariates to distinguish between the units.

This means that the central, $95\%$ say, remaining useful life (RUL) prediction intervals for unit $i$ and unit $j$ are the same if unit $i$ and unit $j$ are the same age.

I would like to assess the predictive capability of my model. I initially thought I could investigate the model's ability to predict remaining useful life.

Suppose all units are new (age 0) at the time of the first prediction and suppose I have 1000 units and suppose initially that I have observed all failure times (no right censoring).

If my model is well calibrated I should expect the coverage probability of the central (or any other) $95\%$ prediction interval to be approximately 0.95. I.e., I should expect the true observations to lie in the central $95\%$ prediction interval approximately $95\%$ of the time.

More generally, I should expect the probability integral transform (PIT) to look uniform. This would mean that approximately $p\%$ of my observations fall below the $p$th quantile of my model.

This is a common way to assess model fit.

Now suppose I have 1000 units (all age 0) but I only observe 50 failures (perhaps I end the test early) and hence 950 units are right censored.

I can only assess whether the true value lies within the prediction interval for the 50 units that fail, and hence I have discarded the 950 right censored units.

I find that my PIT is not uniform and consequently the coverage probability of the central $95\%$ prediction interval is not approximately 0.95.

Initially, I worried that my model fit poorly to the data. However, I think this is to be expected.

All prediction intervals are the same (since all units are new and identical). Suppose the $95\%$ central prediction interval is $[L,U]$. If I had observed all failures, I should expect approximately 25 units ($2.5\%$ of units) to fall below $L$ [i.e. below the lower $2.5\%$ quantile].

I have only observed the earliest 50 failures. So, when I discard the 950 units with no failure times, I find that the coverage probability is very low, which I now think is to be expected.

More specifically, I think I should expect the coverage probability to be 0.5. I.e., I should expect 25 units to fall below $L$ (the earliest 25 failures) and I should expect 25 failures to lie within the prediction interval. Then, since I discarded the 950 right censored units, the coverage probability is expected to be 25/50 = 0.5.

Does this make sense? If yes, are there adjustments I can make to incorporate the unobserved right censored failure times?

Answer 1

You might consider checking that the whole predicted survival distribution fits the proposed model by computing Cox-Snell residuals. This is related to the PIT. See the R rms package val.surv function.

Answer 2

For parametric survival models, comparing the distribution of residuals against what's expected from the chosen parametric family can extend your idea while incorporating censoring. Accelerated failure time models without covariates can be written as

$$\log T = \mu + \sigma W ,$$

where $T$ is survival time, $\mu$ and $\sigma$ are location and scale parameters, and $W$ has a distribution consonant with the chosen parametric family. For such a model, scaled residuals estimated from the model,

$$ \frac{\log T_i - \hat\mu}{\hat\sigma},$$

should be distributed as $W$. For a Weibull model $W$ should be (minimum) extreme value; for a lognormal model, as standard Gaussian.

The trick is to recognize that right censoring of event times also means right censoring of those residuals. You can superimpose a plot of the survival curve of $W$, $1-F(w)$, on the corresponding right-censored survival curve of scaled residuals to go beyond assessing selected quantiles and, in principle, examine the entire distribution.

Start with 1000 event times drawn from a standard exponential distribution, and achieve the 95% censoring you propose in two ways. The first is to censor all but the first 50 event times. The second is to define potential censoring times with an exponential distribution with a rate of 20, and choose the lower of the potential censoring or actual event times. Fit the default Weibull model to each.

library(survival)
set.seed(101)
failTimes <- rexp(1000)
low50 <- failTimes < failTimes[order(failTimes)][51]
low50times <- ifelse(low50, failTimes, failTimes[order(failTimes)][51])
censTimes <- rexp(1000,20)
timesWithCens <- pmin(failTimes,censTimes)
notCens <- failTimes < censTimes 
modW1 <- survreg(Surv(failTimes)~1)
modW2 <- survreg(Surv(low50times,low50)~1)
modW3 <- survreg(Surv(timesWithCens,notCens)~1)

The last model shows how this scaled residuals approach can incorporate right censoring. Plot the scaled residuals, with censored values indicated by crosses at censoring times. Superimpose the corresponding minimum extreme value distribution function from survreg.disributions.

plot(survfit(Surv((log(timesWithCens)-coef(modW3))/modW3$scale,notCens)~1),
    bty="n", main="Residuals distribution", mark.time=TRUE,
    xlab="Scaled residuals", ylab="1 - residuals CDF", ylim=c(0.5,1))
extremeDist <- function(x) survreg.distributions$extreme$density(x,1)[,2]
curve(extremeDist,from = -10, to = 1,lwd=3,col="blue",add=TRUE)
legend("bottomleft",legend="95% random censoring",bty="n")

The model based on censoring the last 950 event times, modW2, only shows censoring of residuals at that highly shared censoring time (not shown). You might want to play similarly with other model types like lognormal or loglogistic.

This evaluation of scaled residuals is equivalent to evaluation of Cox-Snell residuals for parametric models; see Klein and Moeschberger, page 415.

With very heavy right censoring after very early event times, as you posit in the question, you run a risk of fitting the true distribution poorly. Although all 3 models agree reasonably well at very early times, they diverge at later times.

plot(survfit(Surv(failTimes)~1),bty="n",
    xlab="Time", ylab="Fraction surviving",
    xlim=c(0,2),main="Longer survival")   
lines(predict(modW1,newdata=data.frame(failTimes=1),
    type="quantile",p=seq(0.001,0.99,by=0.002)),
    seq(0.999,0.01,by=-.002),type="l",lwd=2)     
lines(predict(modW2,newdata=data.frame(failTimes=1),
    type="quantile",p=seq(0.001,0.99,by=0.002)),
    seq(0.999,0.01,by=-.002),type="l",lwd=2,col="red")
lines(predict(modW3,newdata=data.frame(failTimes=1),
    type="quantile",p=seq(0.001,0.99,by=0.002)),
    seq(0.999,0.01,by=-.002),type="l",lwd=2,col="blue")
legend(-0.14,0.2,
    legend="Black, full fit\nRed, first 50 fit, rest censored\nBlue, random censoring",bty="n")

That's a danger you face with heavy censoring.