R code for constructing likelihood based confidence intervals for the cumulative probabilities and quantiles of an Accelerated Failure Time model

Failure times may be modeled as a function of explanatory variables. The Accelerated Failure Time model (AFT model) is often used for finding the relationship between failure times and explanatory variables.
In a reliability engineering context, for instance, an Accelerated Life Test is often used for determining the effect of variables (such as temperature or voltage) on the durability of some component. For relating the variables to the durability of the component, the reliability engineer usually employs an AFT model.
The following R code computes the likelihood based confidence intervals at specific values of the explanatory variables for 1) the cumulative probabilities, and 2) the quantiles. Note that a reliability engineer may refer to the cumulative probabilities as failure probabilities, and to the quantiles as life time quantiles.

AFT model quantiles and cumulative probabilities

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 normal-approximation confidence intervals for the cumulative probabilities and quantiles at specific values of the explanatory variables.

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

#Device-A data from Meeker & Escobar (1998), Examples 19.3-19.8, pp. 496-504
DevA <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/DeviceA.txt", header=T)
#failure=1, censored=0
DevA$Cens <- as.numeric(DevA$Status)-1
head(DevA)



##fit lognormal-model to data
#see Meeker & Escobar, pp. 498-499

#Arrhenius temperature scale
DevA$tempArr <- 11605/(DevA$Temp + 273.15) #note: temp is in Celsius

#Arrhenius-lognormal regression model (which is an AFT model)
#note: use x=TRUE
mod <- survreg(Surv(Hours, Cens) ~ tempArr, data=DevA, weights=Weight,
               x=TRUE, dist="lognormal")
summary(mod)

#MLE for all regression coefficients
coefsmod <- coef(mod)

#MLE for regression coefficients of predictor variables
betaPv <- coef(mod)[-1]

#MLE for location of fitted distribution
sigma <- mod$scale
names(sigma) <- "sigma"



##FAILURE PROBABILITIES

##likelihood based confidence interval for F(te)
#F(te) is the proportion of units that fails by time te

#specify time
te <- 30000

#function for predicting failure probabilities (=F[te])
predFT <- function (coefs, sigma, te, newdata, dist){
  x <- c(1, newdata)
  mu <- x%*%coefs
  if (dist=="lognormal") {Ft <- pnorm( (log(te)-mu ) / sigma)}
  else if (dist=="weibull") {Ft <- psev( (log(te)-mu ) / sigma)}
  as.vector(Ft)
}

#log-likelihood function
llikmodf <- function (theta, y, x, d, wt=rep(1,length(y)), newdata, dist, fte) {
  sigma <- exp(theta[1])
  betas <- theta[-1]
  
  dotprod <- betas%*%newdata
  
  if (dist=="lognormal") {beta0 <- (log(te) - qnorm(fte)*sigma) - dotprod}
  else if (dist=="weibull") {beta0 <- (log(te) - qsev(fte)*sigma) - dotprod}
  
  betaNew <- c(beta0, betas)
  
  mu <- x%*%betaNew
  
  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 F(te)
profileFt <- function(y, x, d, wt=rep(1,length(y)), parms, newdata, dist, ftes) {
  plke <- rep(0, length(ftes))
  
  for (i in 1:length(ftes)) {
    temp <- optim(parms, llikmodf,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, x=x, d=d, wt=wt, newdata=newdata, dist=dist, fte=ftes[i] )
    plke[i] <- temp$value
  }
  plke
}

#MLE for F(te) at specified temperature
newTemp <- 10 #temperature in Celsius
newArrTemp <- 11605/(newTemp + 273.15) #Arrhenius temperature
(Fte <- predFT(coefs=coefsmod, sigma=sigma, te=te,
               newdata=newArrTemp, dist="lognormal"))

#generate fixed values for the failure probabilities
Fts <- seq(.001, .15, length.out=100)

#likelihood profile
ll <- profileFt(ftes=Fts, y=DevA$Hours, x=mod$x, d=DevA$Cens, wt=DevA$Weight,
                parms=c(log(sigma), betaPv), newdata=newArrTemp, dist="lognormal")

#plot profile likelihood for F(te) including a 95% confidence interval
plot(Fts, ll, type='l', xlab=paste("F(",te,") at DegreesC=", newTemp, sep=""),
     ylab="log-likelihood")
#log-likelihood
loglikw <- mod$loglik[2]
limit <- loglikw - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for Ft at specified temperature
abline(v=Fte, col="red", lty=2)


#compute limits of confidence interval by hand
limitlkh <- function(y, xobs, d, parms, newdata, dist, ftes, wt=rep(1,length(y))) {
  profileFt(y, x=xobs, d, wt, parms, newdata, dist, ftes) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(.001, Fte), y=DevA$Hours, xobs=mod$x, d=DevA$Cens,
                 wt=DevA$Weight, parms=c(log(sigma), betaPv), newdata=newArrTemp,
                 dist="lognormal")
upperci <- fzero(limitlkh, c(Fte, .15), y=DevA$Hours, xobs=mod$x, d=DevA$Cens,
                 wt=DevA$Weight, parms=c(log(sigma), betaPv), newdata=newArrTemp,
                 dist="lognormal")
round(c(lowerci$x, upperci$x), 6)

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



#LIFE TIME QUANTILES

##likelihood based confidence interval for tp
#tp is the time (t) at which a specified proportion (p) of the units fails

#specify proportion
p <-.2

#function for predicting life time quantiles (=tp)
predQp <- function (coefs, sigma, p, newdata, dist){
  x <- c(1, newdata)
  mu <- x%*%coefs
  if (dist=="lognormal") {Qp <- exp(mu + qnorm(p)*sigma)}
  else if (dist=="weibull") {Qp <- exp(mu + qsev(p)*sigma)}
  as.vector(Qp)
}

#log-likelihood function
llikmodq <- function (theta, y, x, d, wt=rep(1,length(y)), newdata, dist, tp) {
  sigma <- exp(theta[1])
  betas <- theta[-1]
  
  dotprod <- betas%*%newdata
  
  if (dist=="lognormal") {beta0 <- (log(tp)-qnorm(p)*sigma) - dotprod}
  else if (dist=="weibull") {beta0 <- (log(tp)-qsev(p)*sigma) - dotprod}
  
  betaNew <- c(beta0, betas)
  
  mu <- x%*%betaNew
  
  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 tp
profileTp <- function(y, x, d, wt=rep(1,length(y)), parms, newdata, dist, tps) {
  plke <- rep(0, length(tps))
  
  for (i in 1:length(tps)) {
    temp <- optim(parms, llikmodq,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, x=x, d=d, wt=wt, newdata=newdata, dist=dist, tp=tps[i] )
    plke[i] <- temp$value
  }
  plke
}

#MLE for tp at specified temperature
newTemp <- 10 #temperature in Celsius
newArrTemp <- 11605/(newTemp + 273.15) #Arrhenius temperature
(Tp <- predQp(coefs=coefsmod, sigma=sigma, p=p,
              newdata=newArrTemp, dist="lognormal"))

#generate fixed values for the life time quantiles
Tps <- seq(35000, 300000, length.out=100)

#likelihood profile
ll <- profileTp(tps=Tps, y=DevA$Hours, x=mod$x, d=DevA$Cens, wt=DevA$Weight,
                parms=c(log(sigma), betaPv), newdata=newArrTemp, dist="lognormal")

#plot profile likelihood for tp including a 95% confidence interval
plot(Tps, ll, type='l', xlab=bquote(t[.(p)]~" at DegreesC="~.(newTemp)),
     ylab="log-likelihood")
#log-likelihood
loglikw <- mod$loglik[2]
limit <- loglikw - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for tp at specified temperature
abline(v=Tp, col="red", lty=2)


#compute limits of confidence interval by hand
limitlkh <- function(y, xobs, d, parms, newdata, dist, tps, wt=rep(1,length(y))) {
  profileTp(y, x=xobs, d, wt, parms, newdata, dist, tps) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(35000, Tp), y=DevA$Hours, xobs=mod$x, d=DevA$Cens,
                 wt=DevA$Weight, parms=c(log(sigma), betaPv), newdata=newArrTemp,
                 dist="lognormal")
upperci <- fzero(limitlkh, c(Tp, 300000), y=DevA$Hours, xobs=mod$x, d=DevA$Cens,
                 wt=DevA$Weight, parms=c(log(sigma), betaPv), newdata=newArrTemp,
                 dist="lognormal")
round(c(lowerci$x, upperci$x), 6)

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

In the following example the failure times are modeled as a function of temperature, voltage, and the interaction between temperature and voltage.

##data

#Capacitor data from Meeker & Escobar (1998), pp. 512-515
Cap <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/Tantalum.txt", header=T)
#failure=1, censored=0
Cap$Cens <- as.numeric(Cap$Status)-1
head(Cap)



##fit lognormal-model to data
#see Meeker & Escobar, pp. 498-499

#Arrhenius temperature scale
Cap$ArrTemp <- 11605/(Cap$DegreesC + 273.15) #note: temp is in Celsius

#log transform of Volts
Cap$logVolts <- log(Cap$Volts)

#fit model to data (which is an AFT model)
#note: use x=TRUE
mod <- survreg(Surv(Hours, Cens) ~ ArrTemp + logVolts + ArrTemp:logVolts,
               data=Cap, weights=Weight, x=TRUE, dist="weibull")
summary(mod)

#MLE for all regression coefficients
coefsmod <- coef(mod)

#MLE for regression coefficients of predictor variables
betaPv <- coef(mod)[-1]

#MLE for location of fitted distribution
sigma <- mod$scale
names(sigma) <- "sigma"



##FAILURE PROBABILITIES

##likelihood based confidence interval for F(te)
#F(te) is the proportion of units that fails by time te

#specify time
te <- 10000

#MLE for F(te) at specified temperature and voltage
newTemp <- 45 #temperature in Celsius
newArrTemp <- 11605/(newTemp + 273.15) #Arrhenius temperature
newVolts <- 57
newlogVolts <- log(newVolts) #log transform of voltage
(Fte <- predFT(coefs=coefsmod, sigma=sigma, te=te,
               newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
               dist="weibull"))

#generate fixed values for the failure probabilities
Fts <- seq(.03, .12, length.out=100)

#likelihood profile
ll <- profileFt(ftes=Fts, y=Cap$Hours, x=mod$x, d=Cap$Cens, wt=Cap$Weight,
                parms=c(log(sigma), betaPv),
                newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
                dist="weibull")

#plot profile likelihood for F(te) including a 95% confidence interval
plot(Fts, ll, type='l', xlab=paste("F(",te,") at DegreesC=", newTemp,
                                   " and Volts=", newVolts, sep=""),
     ylab="log-likelihood")
#log-likelihood
loglikw <- mod$loglik[2]
limit <- loglikw - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for Ft at specified temperature and voltage
abline(v=Fte, col="red", lty=2)


#compute limits of confidence interval by hand
limitlkh <- function(y, xobs, d, parms, newdata, dist, ftes, wt=rep(1,length(y))) {
  profileFt(y, x=xobs, d, wt, parms, newdata, dist, ftes) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(.03, Fte), y=Cap$Hours, xobs=mod$x, d=Cap$Cens,
                 wt=Cap$Weight, parms=c(log(sigma), betaPv),
                 newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
                 dist="weibull")
upperci <- fzero(limitlkh, c(Fte, .12), y=Cap$Hours, xobs=mod$x, d=Cap$Cens,
                 wt=Cap$Weight, parms=c(log(sigma), betaPv),
                 newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
                 dist="weibull")
round(c(lowerci$x, upperci$x), 6)

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



#LIFE TIME QUANTILES

##likelihood based confidence interval for tp
#tp is the time (t) at which a specified proportion (p) of the units fails

#specify proportion
p <-.2

#MLE for tp at specified temperature and voltage
newTemp <- 45 #temperature in Celsius
newArrTemp <- 11605/(newTemp + 273.15) #Arrhenius temperature
newVolts <- 57
newlogVolts <- log(newVolts) #log transform of voltage
(Tp <- predQp(coefs=coefsmod, sigma=sigma, p=p,
              newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
              dist="weibull"))

#generate fixed values for the life time quantiles
Tps <- seq(40000, 1100000, length.out=100)

#likelihood profile
ll <- profileTp(tps=Tps, y=Cap$Hours, x=mod$x, d=Cap$Cens, wt=Cap$Weight,
                parms=c(log(sigma), betaPv),
                newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
                dist="weibull")

#plot profile likelihood for tp including a 95% confidence interval
plot(Tps, ll, type='l', xlab=bquote(t[.(p)]~" at DegreesC="~.(newTemp)~
                                      "and Volts="~.(newVolts)),
     ylab="log-likelihood")
#log-likelihood
loglikw <- mod$loglik[2]
limit <- loglikw - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for tp at specified temperature and voltage
abline(v=Tp, col="red", lty=2)


#compute limits of confidence interval by hand
limitlkh <- function(y, xobs, d, parms, newdata, dist, tps, wt=rep(1,length(y))) {
  profileTp(y, x=xobs, d, wt, parms, newdata, dist, tps) - limit
}

#upper and lower limits of two-sided confidence interval
lowerci <- fzero(limitlkh, c(40000, Tp), y=Cap$Hours, xobs=mod$x, d=Cap$Cens,
                 wt=Cap$Weight, parms=c(log(sigma), betaPv),
                 newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
                 dist="weibull")
upperci <- fzero(limitlkh, c(Tp, 1100000), y=Cap$Hours, xobs=mod$x, d=Cap$Cens,
                 wt=Cap$Weight, parms=c(log(sigma), betaPv),
                 newdata=c(newArrTemp, newlogVolts, newArrTemp*newlogVolts),
                 dist="weibull")
round(c(lowerci$x, upperci$x), 6)

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