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

The following R-code may be used for computing likelihood based confidence intervals for the median residual life.

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

\frac{1}{2}*S(t) = S(t + mdrl)

where mdrl is the median residual life, and S(ยท) the survival function of the failure times.

Computing likelihood based intervals for MDRL
The method for computing the likelihood based confidence interval for the median residual life is similar to the one I previously used for the mean residual life (this blog post explains the method for computing the likelihood based confidence intervals for the mean residual life).

median residual life

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

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

library(SPREDA)

##Median Residual Life (MDRL) for Weibull distribution

##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
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 Median Residual Life (MDRL)

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

#function for computing abs(S(t+mdrl)/S(t)-.5)
XmuMin1 <- function (Xmu, Xsigma, t, mdrl){
  shape <- exp(Xmu)
  scale <- 1/Xsigma
  c1 <- -((t+mdrl)/shape)^scale
  c2 <- -((t)/shape)^scale
  abs((exp(c1))/(exp(c2))-.5)
}

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

llik.mdrl <- function (modpar, x, xd, t, mdrlhat) {
  xsigma <- exp(modpar[1])
  
  xmu <- optim(Xmu, XmuMin1, method='Brent',
               lower=Xmu-10*XmuSE, upper=Xmu+10*XmuSE,
               Xsigma=xsigma, t=t, mdrl=mdrlhat)$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.mdrl <- function(x, xd, t, mdrlhats) {
  plkmdrlhat <- rep(0, length(mdrlhats))
  
  for (i in 1:length(mdrlhats)) {
    temp <- optim(c(log(Xsigma)), llik.mdrl,
                  method='Brent', control=list(fnscale=-1),
                  lower=log(Xsigma)-5*XsigmaSE, upper=log(Xsigma)+5*XsigmaSE,
                  x=x, xd=xd, t=t, mdrlhat=mdrlhats[i] )
    plkmdrlhat[i] <- temp$value
  }
  plkmdrlhat
}


#MLE for MDRL at specified t
mdrl <- function (mdrl, Xmu, Xsigma, t){
  shape <- exp(Xmu)
  scale <- 1/Xsigma
  c1 <- -((t+mdrl)/shape)^scale
  c2 <- -((t)/shape)^scale
  (exp(c1))/(exp(c2))-.5
}

(mdrlhat <- uniroot(mdrl, c(0, 1e+9), t=t, Xmu=Xmu, Xsigma=Xsigma)$root)

#generate fixed values for MDRLhat
MDRLhats <- 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.mdrl function
#that is, change the bounds to, for instance, Xmu+/-5*XmuSE
#should this not help, then decrease the width of the interval
#for possible values of Xsigma in the profile.mdrl function
#that is, change the bounds to, for instance, log(Xsigma)+/-4*XsigmaSE
ll <- profile.mdrl(mdrlhats=MDRLhats, x=Xdata, xd=Xd, t=t)

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

#include lower limit for log-likelihood values
loglikw <- llik.mdrl(modpar=log(Xsigma), x=Xdata, xd=Xd, mdrlhat=mdrlhat, 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 MDRLhat
abline(v=mdrlhat, col="red", lty=2)


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

lowerci <- uniroot(limitlkh, c(11000, mdrlhat), x=Xdata, xd=Xd, t=t)
upperci <- uniroot(limitlkh, c(mdrlhat, 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)





##Median Residual Life (MDRL) for lognormal distribution

library(survival)

#fit lognormal model to 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 Median Residual Life (MDRL)

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

#function for computing abs(S(t+mdrl)/S(t)-.5)
XmuMin1 <- function (Xmu, Xsigma, t, mdrl){
  zz <- (log(t)-Xmu)/Xsigma
  zm <- (log(t+mdrl)-Xmu)/Xsigma
  abs((1-pnorm(zm))/(1-pnorm(zz))-.5)
}

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

llik.mdrl <- function (modpar, x, xd, t, mdrlhat) {
  xsigma <- exp(modpar[1])
  
  xmu <- optim(Xmu, XmuMin1, method='Brent',
               lower=Xmu-10*XmuSE, upper=Xmu+10*XmuSE,
               Xsigma=xsigma, t=t, mdrl=mdrlhat)$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 Xsigma

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


#MLE for MDRL at specified t
mdrl <- function (mdrl, Xmu, Xsigma, t){
  zz <- (log(t)-Xmu)/Xsigma
  zm <- (log(t+mdrl)-Xmu)/Xsigma
  (1-pnorm(zm))/(1-pnorm(zz))-.5
}

(mdrlhat <- uniroot(mdrl, c(0, 1e+9), t=t, Xmu=Xmu, Xsigma=Xsigma)$root)


#generate fixed values for MDRLhat
MDRLhats <- seq(11000, 30000, 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.mdrl function
#that is, change the bounds to, for instance, Xmu+/-5*XmuSE
#should this not help, then decrease the width of the interval
#for possible values of Xsigma in the profile.mdrl function
#that is, change the bounds to, for instance, log(Xsigma)+/-4*XsigmaSE
ll <- profile.mdrl(mdrlhats=MDRLhats, x=Xdata, xd=Xd, t=t)

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

#include lower limit for log-likelihood values
loglikw <- llik.mdrl(modpar=log(Xsigma), x=Xdata, xd=Xd, mdrlhat=mdrlhat, 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 MDRLhat
abline(v=mdrlhat, col="red", lty=2)


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

lowerci <- uniroot(limitlkh, c(11000, mdrlhat), x=Xdata, xd=Xd, t=t)
upperci <- uniroot(limitlkh, c(mdrlhat, 30000), 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)