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

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

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

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