R code for constructing likelihood based confidence intervals for parametric survival models

The following R code may be used for constructing likelihood based intervals for parametric survival models (such as the Weibull model). These likelihood based intervals are also known as likelihood ratio bounds, or profile likelihood intervals.

The code constructs confidence intervals for the two distribution parameters of the parametric surival model (location μ and scale σ), life time quantiles tp, and failure probabilities F(te).

The R code focuses on the Weibull distribution, but can easily be adapted for modeling with other distributions (e.g., lognormal distribution).

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)
library(pracma)

##data
#Shock absorber failure data (see: Meeker & Escobar, Examples 8.4-8.7, pp. 180-183)
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 data
mod <- survreg(Surv(Distance,Cens)~1, data=ShAbs, dist="weibull")
summary(mod)

#extract MLE for mu and sigma
mu <- coef(mod)[1]
names(mu) <- "mu"
sigma <- mod$scale
names(sigma) <- "sigma"

#extract log-likelihood
loglikmod <- mod$loglik[2]

Likelihood based confidence intervals for distribution parameter σ

##profile likelihood sigma
#sigma is the scale of the distribution

#log-likelihood function
#note: the wt argument weights observations
lliksigma <- function (mu, y, d, wt=rep(1, length(y)), dist, sigma) {
  z <- (log(y)-mu)/sigma
  
  #log-likelihood for exact (not censored) and right-censored observations
  if (dist=="lognormal") {
    sum(wt*log((( 1/(sigma*y) * dnorm(z)) ^d)*( (1-pnorm(z)) ^(1-d) )))}
  else if (dist=="weibull") {
    sum(wt*log((( 1/(sigma*y) * dsev(z)) ^d)*( (1-psev(z)) ^(1-d) )))}
}

#profile function for sigma
#note: the wt argument weights observations
profile.sigma <- function(y, d, wt=rep(1, length(y)), dist, sigmas) {
  plkmu <- rep(0, length(sigmas))
  
  for (i in 1:length(sigmas)) {
    temp <- optim(mu, lliksigma,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, d=d, wt=wt, dist=dist, sigma=sigmas[i] )
    plkmu[i] <- temp$value
  }
  plkmu
}

#MLE for sigma
sigma
#fixed values of sigma
sigmas <- seq(0.2, 0.6, length.out=50)

#likelihood profile
ll <- profile.sigma(sigmas=sigmas, y=ShAbs$Distance, d=ShAbs$Cens, dist="weibull")

#plot profile likelihood for sigma including a 95% confidence interval
plot(sigmas,ll, type='l', xlab=expression(sigma), ylab="log-likelihood")
#include lower limit for log-likelihood values
limit <- loglikmod - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for sigma
abline(v=sigma, col="red", lty=2)
#gray lines represent the boundaries of the 95% likelihood based interval
#these boundaries were reported by Meeker & Escobar, Table 8.1, p. 178
abline(v=c(.21,.527), col="gray", lty=3)


#compute limits of confidence interval by hand
limitlkh <- function(y, d, dist, sigmas, wt=rep(1, length(y))) {
  profile.sigma(y, d, wt, dist, sigmas) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(.2, sigma), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
upperci <- fzero(limitlkh, c(sigma, 0.6), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
round(c(lowerci$x, upperci$x), 6)

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

Likelihood based confidence interval for distribution parameter μ

##profile likelihood mu
#mu is the location of the distribution

#log-likelihood function
#note: the wt argument weights observations
llikmu <- function (sigma, y, d, wt=rep(1, length(y)), dist, mu) {
  sigma <- exp(sigma)
  z <- (log(y)-mu)/sigma
  
  #log-likelihood for exact (not censored) and right-censored observations
  if (dist=="lognormal") {
    sum(wt*log((( 1/(sigma*y) * dnorm(z)) ^d)*( (1-pnorm(z)) ^(1-d) )))}
  else if (dist=="weibull") {
    sum(wt*log((( 1/(sigma*y) * dsev(z)) ^d)*( (1-psev(z)) ^(1-d) )))}
}

#profile function for sigma
#note: the wt argument weights observations
profile.mu <- function(y, d, wt=rep(1, length(y)), dist, mus) {
  plkmu <- rep(0, length(mus))
  
  for (i in 1:length(mus)) {
    temp <- optim(log(sigma), llikmu,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, d=d, wt=wt, dist=dist, mu=mus[i] )
    plkmu[i] <- temp$value
  }
  plkmu
}

#MLE for mu
mu
#fixed values of mu
mus <- seq(10, 10.6, length.out=50)

#likelihood profile
ll <- profile.mu(mus=mus, y=ShAbs$Distance, d=ShAbs$Cens, dist="weibull")

#plot profile likelihood for mu including a 95% confidence interval
plot(mus,ll, type='l', xlab=expression(mu), ylab="log-likelihood")
#include lower limit for log-likelihood values
limit <- loglikmod - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for mu
abline(v=mu, col="red", lty=2)
#boundaries of the 95% interval as reported by Meeker & Escobar, Table 8.1, p. 178
abline(v=c(10.06,10.54), col="gray", lty=3)


#compute limits of confidence interval by hand
limitlkh <- function(y, d, dist, mus, wt=rep(1, length(y))) {
  profile.mu(y, d, wt, dist, mus) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(10, mu), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
upperci <- fzero(limitlkh, c(mu, 10.6), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
round(c(lowerci$x, upperci$x), 6)

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

Likelihood based confidence interval for life time quantile tp

##profile likelihood tp
#tp is the time (t) at which a specified proportion (p) of the units fails

#specify proportion
p <-.1

#log-likelihood function
#note: the wt argument weights observations
lliktp <- function (sigma, y, d, wt=rep(1, length(y)), dist, tp) {
  sigma <- exp(sigma)
  
  #log-likelihood for exact (not censored) and right-censored observations
  if (dist=="lognormal") {
    mu <- log(tp)-qnorm(p)*sigma
    z <- (log(y)-mu)/sigma
    sum(wt*log((( 1/(sigma*y) * dnorm(z)) ^d)*( (1-pnorm(z)) ^(1-d) )))}
  else if (dist=="weibull") {
    mu <- log(tp)-qsev(p)*sigma
    z <- (log(y)-mu)/sigma
    sum(wt*log((( 1/(sigma*y) * dsev(z)) ^d)*( (1-psev(z)) ^(1-d) )))}
}

#profile function for tp
#note: the wt argument weights observations
profile.tp <- function(y, d, wt=rep(1, length(y)), dist, tps) {
  plkp <- rep(0, length(tps))
  
  for (i in 1:length(tps)) {
    temp <- optim(log(sigma), lliktp,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, d=d, wt=wt, dist=dist, tp=tps[i] )
    plkp[i] <- temp$value
  }
  plkp
}

#MLE for tp
#note: for lognormal distribution use: (tpMLE <- exp(mu+qnorm(p)*sigma))
(tpMLE <- exp(mu+qsev(p)*sigma))
#fixed values of tp
tps <- seq(8000, 18000, length.out=50)

#likelihood profile
ll <- profile.tp(tps=tps, y=ShAbs$Distance, d=ShAbs$Cens, dist="weibull")

#plot profile likelihood for tp including a 95% confidence interval
plot(tps,ll, type='l', xlab=bquote(t[.(p)]), ylab="log-likelihood")
#include lower limit for log-likelihood values
limit <- loglikmod - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for tp
abline(v=tpMLE, col="red", lty=2)
#boundaries of the 95% interval as reported by Meeker & Escobar, Table 8.1, p. 178
abline(v=c(9400,17300), col="gray", lty=3)


#compute limits of confidence interval by hand
limitlkh <- function(y, d, dist, tps, wt=rep(1, length(y))) {
  profile.tp(y, d, wt, dist, tps) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(8000, tpMLE), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
upperci <- fzero(limitlkh, c(tpMLE, 18000), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
round(c(lowerci$x, upperci$x), 6)

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



##it is also possible to use the Brent optimization algorithm

#specify proportion
p <-.001

#profile function for tp (using Brent)
#note: the wt argument weights observations
sigmaSE <- sqrt(diag(vcov(mod)))[2] #extract standard error for log(sigma)

profile.tp <- function(y, d, wt=rep(1, length(y)), dist, tps) {
  plkp <- rep(0, length(tps))
  
  for (i in 1:length(tps)) {
    temp <- optim(log(sigma), lliktp,
                  method="Brent", control=list(fnscale=-1),
                  lower=log(sigma)-20*sigmaSE, upper=log(sigma)+20*sigmaSE,
                  y=y, d=d, wt=wt, dist=dist, tp=tps[i] )
    plkp[i] <- temp$value
  }
  plkp
}

#MLE for tp
#note: for lognormal distribution use: (tpMLE <- exp(mu+qnorm(p)*sigma))
(tpMLE <- exp(mu+qsev(p)*sigma))
#fixed values of tp
tps <- seq(700, 7000, length.out=50)

#likelihood profile
ll <- profile.tp(tps=tps, y=ShAbs$Distance, d=ShAbs$Cens, dist="weibull")

#plot profile likelihood for tp including a 95% confidence interval
plot(tps,ll, type='l', xlab=bquote(t[.(p)]), ylab="log-likelihood")
#include lower limit for log-likelihood values
limit <- loglikmod - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for tp
abline(v=tpMLE, col="red", lty=2)



#compute limits of confidence interval by hand
limitlkh <- function(y, d, dist, tps, wt=rep(1, length(y))) {
  profile.tp(y, d, wt, dist, tps) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(700, tpMLE), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
upperci <- fzero(limitlkh, c(tpMLE, 7000), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
round(c(lowerci$x, upperci$x), 6)

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

Likelihood based confidence interval for failure probability F(te)

##profile likelihood F(te)
#F(te) is the proportion of units that fails by time te

#specify time
te <- 10000

#log-likelihood function
#note: the wt argument weights observations
llikfte <- function (sigma, y, d, wt=rep(1, length(y)), dist, fte) {
  sigma <- exp(sigma)
  
  #log-likelihood for exact (not censored) and right-censored observations
  if (dist=="lognormal") {
    mu <- log(te) - qnorm(fte)*sigma
    z <- (log(y)-mu)/sigma
    sum(wt*log((( 1/(sigma*y) * dnorm(z)) ^d)*( (1-pnorm(z)) ^(1-d) )))}
  else if (dist=="weibull") {
    mu <- log(te) - qsev(fte)*sigma
    z <- (log(y)-mu)/sigma
    sum(wt*log((( 1/(sigma*y) * dsev(z)) ^d)*( (1-psev(z)) ^(1-d) )))}
}

#profile function for F(te)
#note: the wt argument weights observations
profile.fte <- function(y, d, wt=rep(1, length(y)), dist, ftes) {
  plke <- rep(0, length(ftes))
  
  for (i in 1:length(ftes)) {
    temp <- optim(log(sigma), llikfte,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, d=d, wt=wt, dist=dist, fte=ftes[i] )
    plke[i] <- temp$value
  }
  plke
}

#MLE for F(te)
#note: for lognormal distribution use: (fteMLE <- pnorm( (log(te)-mu) / sigma))
(fteMLE <- psev( (log(te)-mu) / sigma))
#fixed values of F(te)
ftes <- seq(.001, .14, length.out=50)

#likelihood profile
ll <- profile.fte(ftes=ftes, y=ShAbs$Distance, d=ShAbs$Cens, dist="weibull")

#plot profile likelihood for F(te) including a 95% confidence interval
plot(ftes,ll, type='l', xlab=paste("F(",te,")", sep=""), ylab="log-likelihood")
#include lower limit for log-likelihood values
limit <- loglikmod - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for F(te)
abline(v=fteMLE, col="red", lty=2)
#boundaries of the 95% interval as reported by Meeker & Escobar, Table 8.1, p. 178
abline(v=c(.0092,.1136), col="gray", lty=3)


#compute limits of confidence interval by hand
limitlkh <- function(y, d, dist, ftes, wt=rep(1, length(y))) {
  profile.fte(y, d, wt, dist, ftes) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(.001, fteMLE), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
upperci <- fzero(limitlkh, c(fteMLE, .14), y=ShAbs$Distance, d=ShAbs$Cens,
                 dist="weibull")
round(c(lowerci$x, upperci$x), 6)

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