1
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ You don't only observe the failures, but also the fact that 950 units didn't fail., so obviously you should take it into account. Can you clarify what exactly you are asking besides that ? $\endgroup$
    – J. Delaney
    Commented May 30, 2022 at 17:15
  • $\begingroup$ The main question would be: Suppose I have a prediction interval $[L,U]$ for a right-censored observation, i.e., I know the failure time for this unit is $T>t$ (where $t$ is the age of the unit at the end of the study), how can I determine if the actual failure time is within the prediction interval? It is easy to calculate the coverage probability when I observe the actual failure time, $t$. I can just count the number of times they fall inside the interval, but how do I count the right-censored times? $\endgroup$
    – JLee
    Commented May 30, 2022 at 17:38
  • $\begingroup$ Well obviously you have no way of knowing what is the actual failure time if you don't observe it. You can only calculate upper/lower bounds by assuming extreme cases, but if some units have a much longer remaining useful time than your model predicts, you have no way of knowing it. $\endgroup$
    – J. Delaney
    Commented May 30, 2022 at 19:39

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
2
+100
$\begingroup$

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")

scaled residuals distribution

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")

compare 3 survival models on exponential data

That's a danger you face with heavy censoring.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.