R code for constructing likelihood based confidence intervals for the mean residual life

The following R-code may be used for computing likelihood based confidence intervals for the mean (or expected) residual life.

Definition of the mean residual life
The mean residual life (or expected remaining lifetime) of a unit at time t is given by:

MRL(t) = \frac{\int^{\infty}_{t} (u-t)f(u)\ du}{S(t)} = \frac{\int^{\infty}_{t} S(u) \ du}{S(t)}

where MRL(t) is the mean residual life at time t, f(·) the probability density function of the failure times, and S(·) the survival function.

Computing likelihood based intervals for MRL
For computing the likelihood, a location-scale based parametric distribution will be employed. The location of this distribution will be denoted by μx, and its scale by σx.
A consequence of using a location-scale based distribution is that both f(·) and S(·) in the above equation for MRL will be some function of μx and σx. Furthermore, the likelihood function (denoted by L) also will be some function of μx and σx.

The likelihood function L will be used in computing the likelihood based intervals for MRL. However, for computing these MRL confidence intervals we need to reparameterize L. For this reparameterization to work, notice that MRL(t) is defined as a function of the parameters μx and σx , denoted by MRL=f(μx, σx). As a consequence, it is also possible to express μx as some function of MRL and σx , denoted by μx=g(MRL, σx).

The reparameterization of the likelihood function L requires that μx is substituted by g(MRL, σx) in the original likelihood function for L. In that way, L becomes a function of MRL and σx. This latter parameterization of L makes it possible to construct likelihood based confidence intervals for MRL.

For our specific reparameterization of L, it is thus needed to express μx as a function of the parameters MRL and σx. In other words, we need to compute μx given any values for MRL and σx.
The definition of MRL(t) above showed that MRL=f(μx, σx), and therefore that MRL – f(μx, σx) = 0. This latter equation may yield our sought for value for μx, given any values for the remaining parameters in the equation (i.e., MRL and σx).
One way of obtaining this value for μx is by relying on numerical optimization. More specifically, μx can be found for any given values MRL and σx by minimizing abs(MRL – f(μx, σx)) over all possible values of μx.

Mean residual life ci

A description of likelihood based confidence intervals can be found in Meeker and Escobar’s 1998 book Statistical methods for reliability data (Chapter 8, pp. 173-186).

See this blog post for computing likelihood based confidence intervals for the median residual life.

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).
See this page for an overview of all of Stefan’s R code blog posts.

R code

##Mean Residual Life (MRL) for Weibull distribution

library(SPREDA)

##data

#Shock absorber failure data (see: Meeker & Escobar, Example 3.8, p. 59)
ShAbs <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/ShockAbsorber.txt", header=T)

#failure=1, censored=0
ShAbs$Cens <- as.numeric(ShAbs$Censor)-1
head(ShAbs)

#fit Weibull model to shock absorber failure data
modX <- survreg(Surv(Distance,Cens)~1, data=ShAbs, dist="weibull")
summary(modX)

#obtain MLE for mu and sigma
Xmu <- coef(modX)[1]
names(Xmu) <- "Xmu"
Xsigma <- modX$scale
names(Xsigma) <- "Xsigma"
#log-likelihood
(loglikw <- modX$loglik[1])

#construct variables for analysis
Xdata <- ShAbs$Distance
Xd <- ShAbs$Cens



##likelihood based interval for MRL

#specify time at which to compute the MRL and its confidence interval
t <- 10000

#specify integrand for numerical integration
integrand <- function(u, Xmu, Xsigma) {
  shape <- exp(Xmu)
  scale <- 1/Xsigma
  c1 <- -((u)/shape)^scale
  exp(c1)
}

#function for computing abs(MRLhat-f(Xmu,Xsigma)/(S[t]))
XmuMin1 <- function (XmuPar, Xsigma, MRLhat, t){
  #for stability of the optimization algorithm, replace the upper limit (=Inf) by:
  uplimit <- exp(qsev(.99999999)*3*Xsigma+XmuPar)
  #compute abs(MRLhat-f(Xmu,Xsigma)/(S[t]))
  shape <- exp(XmuPar)
  scale <- 1/Xsigma
  c1 <- -((t)/shape)^scale
  abs(MRLhat-integrate(integrand, t, uplimit, Xmu=XmuPar, Xsigma=Xsigma)$value/
        (exp(c1)))}

#log-likelihood function
XmuSE <- sqrt(diag(vcov(modX)))[1] #extract standard error for Xmu

llik.mrl <- function (modpar, x, xd, MRLhat, t) {
  xsigma <- exp(modpar[1])
  
  xmu <- optim(par=Xmu, XmuMin1, method="Brent",
               lower=Xmu-20*XmuSE, upper=Xmu+20*XmuSE,
               MRLhat=MRLhat, Xsigma=xsigma, t=t)$par
  
  z <- (log(x)-xmu)/xsigma
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log((( 1/(xsigma*x) * dsev(z)) ^xd)*( (1-psev(z)) ^(1-xd) )))
}

#profile function
XsigmaSE <- sqrt(diag(vcov(modX)))[2] #extract standard error for Xsigma

profile.mrl <- function(x, xd, t, MRLhats) {
  plkmrlhat <- rep(0, length(MRLhats))
  
  for (i in 1:length(MRLhats)) {
    temp <- optim(c(log(Xsigma)), llik.mrl,
                  method='Brent', control=list(fnscale=-1),
                  lower=log(Xsigma)-5*XsigmaSE, upper=log(Xsigma)+5*XsigmaSE,
                  x=x, xd=xd, t=t, MRLhat=MRLhats[i] )
    plkmrlhat[i] <- temp$value
  }
  plkmrlhat
}

#MLE for MRL at specified t
uplimit <- exp(qsev(.99999999)*3*Xsigma+Xmu)
shapeT <- exp(Xmu)
scaleT <- 1/Xsigma
c1T <- -((t)/shapeT)^scaleT
(MRLhat <- as.vector(integrate(integrand, t, uplimit, Xmu=Xmu, Xsigma=Xsigma)$value/
                       (exp(c1T))))
#compare this estimate of MRL with that where the MRL is computed by
#using upper=Inf as limit for the integral
as.vector(integrate(integrand, t, Inf, Xmu=Xmu, Xsigma=Xsigma)$value/(exp(c1T)))

#generate fixed values for MRLhat
MRLhats <- seq(11000, 25000, length.out=100)

#likelihood profile
#note: in case of warning messages, decrease the width of the interval
#for possible values of Xmu in the llik.mrl function
#that is, change the bounds to, for instance, Xmu+/-10*XmuSE, or Xmu+/-5*XmuSE
#should this not help, then decrease the width of the interval
#for possible values of Xsigma in the profile.mrl function
#that is, change the bounds to, for instance, log(Xsigma)+/-4*XsigmaSE
ll <- profile.mrl(MRLhats=MRLhats, x=Xdata, xd=Xd, t=t)

#plot profile likelihood for MRLhat
plot(MRLhats,ll, type='l', xlab=paste("Mean Residual Life (at t=",t,")", sep=""), ylab="log-likelihood")

#include lower limit for log-likelihood values
loglikw <- llik.mrl(modpar=log(Xsigma),
                    x=Xdata, xd=Xd, MRLhat=MRLhat, t=t)
alpha <- .10 #90% confidence interval
limit <- loglikw - .5*qchisq(1-alpha, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for MRLhat
abline(v=MRLhat, col="red", lty=2)


#compute limits of confidence interval by hand
limitlkh <- function(x, xd, t, MRLhats) {
  profile.mrl(x, xd, t, MRLhats) - limit
}

lowerci <- uniroot(limitlkh, c(11000, MRLhat), x=Xdata, xd=Xd, t=t)
upperci <- uniroot(limitlkh, c(MRLhat, 25000), x=Xdata, xd=Xd, t=t)
c(lowerci$root, upperci$root)

#include confidence bounds in plot
abline(v=lowerci$root, col="darkgreen", lty=2)
abline(v=upperci$root, col="darkgreen", lty=2)





##Mean Residual Life (MRL) for lognormal distribution

library(survival)

#fit lognormal model to shock absorber failure data
modX <- survreg(Surv(Distance,Cens)~1, data=ShAbs, dist="lognormal")
summary(modX)

#obtain MLE for mu and sigma
Xmu <- coef(modX)[1]
names(Xmu) <- "Xmu"
Xsigma <- modX$scale
names(Xsigma) <- "Xsigma"
#log-likelihood
(loglikw <- modX$loglik[1])

#construct variables for analysis
Xdata <- ShAbs$Distance
Xd <- ShAbs$Cens



##likelihood based interval for MRL

#specify time at which to compute the MRL and its confidence interval
t <- 10000

#specify integrand for numerical integration
integrand <- function(u, Xmu, Xsigma) {
  1-pnorm((log(u)-Xmu)/Xsigma)
}

#function for computing abs(MRLhat-f(Xmu,Xsigma)/(S[t]))
XmuMin1 <- function (XmuPar, Xsigma, MRLhat, t){
  uplimit <- exp(qnorm(.99999999)*Xsigma+XmuPar)
  abs(MRLhat-integrate(integrand, t, uplimit, Xmu=XmuPar, Xsigma=Xsigma)$value/
        (1-pnorm((log(t)-XmuPar)/Xsigma)))}

#log-likelihood function
XmuSE <- sqrt(diag(vcov(modX)))[1] #extract standard error for Xmu

llik.mrl <- function (modpar, x, xd, MRLhat, t) {
  xsigma <- exp(modpar[1])
  
  xmu <- optim(par=Xmu, XmuMin1, method="Brent",
               lower=Xmu-20*XmuSE, upper=Xmu+20*XmuSE,
               MRLhat=MRLhat, Xsigma=xsigma, t=t)$par
  
  z <- (log(x)-xmu)/xsigma
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log((( 1/(xsigma*x) * dnorm(z)) ^xd)*( (1-pnorm(z)) ^(1-xd) )))
}

#profile function
XsigmaSE <- sqrt(diag(vcov(modX)))[2] #extract standard error for Xmu

profile.mrl <- function(x, xd, t, MRLhats) {
  plkmrlhat <- rep(0, length(MRLhats))
  
  for (i in 1:length(MRLhats)) {
    temp <- optim(c(log(Xsigma)), llik.mrl,
                  method='Brent', control=list(fnscale=-1),
                  lower=log(Xsigma)-5*XsigmaSE, upper=log(Xsigma)+5*XsigmaSE,
                  x=x, xd=xd, t=t, MRLhat=MRLhats[i] )
    plkmrlhat[i] <- temp$value
  }
  plkmrlhat
}

#MLE for MRL at specified t
uplimit <- exp(qnorm(.99999999)*Xsigma+Xmu)
(MRLhat <- as.vector(integrate(integrand, t, uplimit, Xmu=Xmu, Xsigma=Xsigma)$value/
                       (1-pnorm((log(t)-Xmu)/Xsigma))))
#compare this estimate of MRL with that where the MRL is computed by
#using upper=Inf as limit for the integral
as.vector(integrate(integrand, t, Inf, Xmu=Xmu, Xsigma=Xsigma)$value/
            (1-pnorm((log(t)-Xmu)/Xsigma)))

#generate fixed values for MRLhat
MRLhats <- seq(11000, 40000, length.out=100)

#likelihood profile
#note: in case of warning messages, decrease the width of the interval
#for possible values of Xmu in the llik.mrl function
#that is, change the bounds to, for instance, Xmu+/-10*XmuSE, or Xmu+/-5*XmuSE
#should this not help, then decrease the width of the interval
#for possible values of Xsigma in the profile.mrl function
#that is, change the bounds to, for instance, log(Xsigma)+/-4*XsigmaSE
ll <- profile.mrl(MRLhats=MRLhats, x=Xdata, xd=Xd, t=t)

#plot profile likelihood for MRLhat
plot(MRLhats,ll, type='l', xlab=paste("Mean Residual Life (at t=",t,")", sep=""), ylab="log-likelihood")

#include lower limit for log-likelihood values
loglikw <- llik.mrl(modpar=log(Xsigma),
                    x=Xdata, xd=Xd, MRLhat=MRLhat, t=t)
alpha <- .10 #90% confidence interval
limit <- loglikw - .5*qchisq(1-alpha, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for MRLhat
abline(v=MRLhat, col="red", lty=2)


#compute limits of confidence interval by hand
limitlkh <- function(x, xd, t, MRLhats) {
  profile.mrl(x, xd, t, MRLhats) - limit
}

lowerci <- uniroot(limitlkh, c(11000, MRLhat), x=Xdata, xd=Xd, t=t)
upperci <- uniroot(limitlkh, c(MRLhat, 40000), x=Xdata, xd=Xd, t=t)
c(lowerci$root, upperci$root)

#include confidence bounds in plot
abline(v=lowerci$root, col="darkgreen", lty=2)
abline(v=upperci$root, col="darkgreen", lty=2)