R code for constructing likelihood based confidence intervals for MTTF

The following R code may be used for constructing likelihood based intervals for the Mean Time To Failure (MTTF). These likelihood based intervals are also known as likelihood ratio bounds, or profile likelihood intervals.

The MTTF is defined as the mean of a failure-time distribution. The R code calculates the MTTF and its likelihood based confidence interval for the lognormal and Weibull failure-time distribution.

MTTF confidence interval

See this blog post for constructing likelihood based intervals for the Median Time To Failure.

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



##MTTF lognormal distribution

##obtain start values

#fit lognormal model
mod <- survreg(Surv(Distance,Cens)~1, data=ShAbs, dist="lognormal")
summary(mod)

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



##profile likelihood MTTF
#MTTF is equal to the mean of the lognormal distribution,
#which is defined as exp(mu + sigma^2/2)

#log-likelihood function
lliknorm.mttf <- function (sigma, y, d, MTTF) {
  mu <- log(MTTF) - sigma^2/2
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log((( 1/(sigma*y) * dnorm( ((log(y)-mu)/sigma) )) ^d)*
            ( (1-pnorm( ((log(y)-mu)/sigma) )) ^(1-d) )))
}

#profile function for MTTF
profile.mttf <- function(y, d, MTTFs) {
  plkmttf <- rep(0, length(MTTFs))
  
  for (i in 1:length(MTTFs)) {
    temp <- optim(sigma, lliknorm.mttf,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, d=d, MTTF=MTTFs[i] )
    plkmttf[i] <- temp$value
  }
  plkmttf
}

#MLE for MTTF
(MLEmttf <- exp(mu + sigma^2/2))
#fixed values of MTTF
MTTFs <- seq(20000, 60000, 500)

#likelihood profile
ll <- profile.mttf(MTTFs, y=ShAbs$Distance, d=ShAbs$Cens)

#plot profile likelihood for MTTF
plot(MTTFs,ll, type='l', xlab="MTTF", ylab="log-likelihood")
#include lower limit for log-likelihood values
alpha <- .05 #95% confidence interval
limit <- loglikw - .5*qchisq(1-alpha, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for MTTF
abline(v=MLEmttf, col="red", lty=2)

#calculate the 95% confidence interval by hand
limitlkh <- function(y, d, mttfs) {
  profile.mttf(y, d, mttfs) - limit
}

lowerci <- uniroot(limitlkh,c(20000,MLEmttf), y=ShAbs$Distance, d=ShAbs$Cens)
upperci <- uniroot(limitlkh,c(MLEmttf,60000), y=ShAbs$Distance, d=ShAbs$Cens)
c(lowerci$root, upperci$root)



##MTTF Weibull distribution

##obtain start values

#fit Weibull model
mod <- survreg(Surv(Distance,Cens)~1, data=ShAbs, dist="weibull")
summary(mod)

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



##profile likelihood MTTF
#MTTF is equal to the mean of the Weibull distribution,
#which is defined as exp(mu)*gamma(1+sigma)

#log-likelihood function
llikweibull.mttf <- function (sigma, y, d, MTTF) {
  mu <- log(MTTF) - log(gamma(1+sigma))
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log((( 1/(sigma*y) * dsev( ((log(y)-mu)/sigma) )) ^d)*
            ( (1-psev( ((log(y)-mu)/sigma) )) ^(1-d) )))
}

#profile function for MTTF
profile.mttf <- function(y, d, MTTFs) {
  plkmttf <- rep(0, length(MTTFs))
  
  for (i in 1:length(MTTFs)) {
    temp <- optim(sigma, llikweibull.mttf,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, d=d, MTTF=MTTFs[i] )
    plkmttf[i] <- temp$value
  }
  plkmttf
}

#MLE for MTTF
(MLEmttf <- exp(mu)*gamma(1+sigma))
#fixed values of MTTF
MTTFs <- seq(20000, 60000, 500)

#likelihood profile
ll <- profile.mttf(MTTFs, y=ShAbs$Distance, d=ShAbs$Cens)

#plot profile likelihood
plot(MTTFs,ll, type='l', xlab="MTTF", ylab="log-likelihood")
#include lower limit for log-likelihood values
alpha <- .05 #95% confidence interval
limit <- loglikw - .5*qchisq(1-alpha, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for MTTF
abline(v=MLEmttf, col="red", lty=2)

#calculate the 95% confidence interval by hand
limitlkh <- function(y, d, mttfs) {
  profile.mttf(y, d, mttfs) - limit
}

lowerci <- uniroot(limitlkh,c(20000,MLEmttf), y=ShAbs$Distance, d=ShAbs$Cens)
upperci <- uniroot(limitlkh,c(MLEmttf,40000), y=ShAbs$Distance, d=ShAbs$Cens)
c(lowerci$root, upperci$root)