R code for constructing likelihood based confidence intervals for the regression coefficients of an Accelerated Failure Time model

The following R code computes likelihood based confidence intervals for the regression coefficients of an Accelerated Failure Time model. The AFT model is a parametric survival model.

AFT model coefficients

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

#log-likelihood function
llikbeta <- function (theta, y, d, wt, censoring, xs, dist, betaFixed) {
  sigma <- exp(theta[1])
  coefs <- c(betaFixed, theta[-1])
  
  mu <- xs%*%coefs
  
  if (censoring=="right"){
    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) )))}
  }
  
  else if (censoring=="interval"){
    z1 <- (log(y[,1])-mu)/sigma
    z2 <- (log(y[,2])-mu)/sigma
    
    if (dist=="lognormal"){
      l <- ifelse(d==0, 1-pnorm(z1), #right censoring
                  ifelse(d==1, 1/(sigma*y[,1]) * dnorm(z1), #exact observations
                         ifelse(d==2, pnorm(z1), #left censoring
                                pnorm(z2)-pnorm(z1) )))} #interval censoring
    else if (dist=="weibull"){
      l <- ifelse(d==0, 1-psev(z1), #right censoring
                  ifelse(d==1, 1/(sigma*y[,1]) * dsev(z1), #exact observations
                         ifelse(d==2, psev(z1), #left censoring
                                psev(z2)-psev(z1) )))} #interval censoring
    -wt%*%log(l)}
}

#profiling using Nelder-Mead method
profileLogLik <- function(fixedValues, parms, model, xs, rsimplex){
  
  #weights
  if (is.null(model$weights)) wt <- rep(1,length(model$y[,1])) else
    wt <- model$weights
  
  #censoring type
  censoring <- ifelse(ncol(model$y)==2, "right", "interval")
  
  #regsimp=FALSE: initial simplex is not a regular simplex
  if (rsimplex==TRUE) controlNM <- list(tol=1e-10) else
    controlNM <- list(tol=1e-10, regsimp=FALSE)
  
  plke <- rep(0, length(fixedValues))
  
  if (censoring=="right") {
    for (i in 1:length(fixedValues)) {
      temp <- nmk(parms, llikbeta,
                  control=controlNM,
                  y=exp(model$y[,1]), d=model$y[,2], wt=wt, censoring=censoring,
                  xs=xs, dist=model$dist, betaFixed=fixedValues[i] )
      plke[i] <- temp$value}
  }
  else if (censoring=="interval") {
    for (i in 1:length(fixedValues)) {
      temp <- nmk(parms, llikbeta,
                  control=controlNM,
                  y=exp(model$y[,1:2]), d=model$y[,3], wt=wt, censoring=censoring,
                  xs=xs, dist=model$dist, betaFixed=fixedValues[i] )
      plke[i] <- temp$value}
  }
  
  -plke
}

#profile function for coefficients of an AFT model
profileBeta <- function(model, whichPar, range,
                        nvalues=50, alpha=.95, rsimplex=TRUE) {
  fixedValues <- seq(range[1], range[2], length.out=nvalues)
  
  beta <- coef(model)[whichPar]
  parms <- c(log(model$scale), coef(model)[-whichPar])
  
  xBeta <- model$x[, whichPar]
  xParms <- model$x[, -whichPar]
  
  xs <- cbind(xBeta, xParms)
  
  #profile regression coefficient (=beta)
  ll <- profileLogLik(fixedValues, parms, model, xs, rsimplex)
  
  #compute confidence bounds
  loglikw <- model$loglik[2]
  limit <- loglikw - .5*qchisq(alpha, df=1)
  
  limitlkh <- function(fixedValues, parms, model, xs, rsimplex) {
    profileLogLik(fixedValues, parms, model, xs, rsimplex) - limit}
  
  lowerci <- NA
  upperci <- NA
  
  try(lowerci <- round(fzero(limitlkh, c(range[1], beta), parms=parms,
                             model=model, xs=xs, rsimplex=rsimplex)$x, 6), TRUE)
  try(upperci <- round(fzero(limitlkh, c(beta, range[2]), parms=parms,
                             model=model, xs=xs, rsimplex=rsimplex)$x, 6), TRUE)
  
  #plot profile
  plot(fixedValues, ll, type='l', ylab="log-likelihood",
       xlab=paste("coefficient = ", names(coef(model)[whichPar]), sep=""))
  abline(h=limit, col="blue",lty=2) #log-likelihood limit
  abline(v=coef(model)[whichPar], col="red", lty=2) #include MLE of coefficient
  #include interval bounds
  if (!is.na(lowerci)) abline(v=lowerci, col="darkgreen", lty=2)
  if (!is.na(upperci)) abline(v=upperci, col="darkgreen", lty=2)
  #return bounds
  c(lcl=lowerci, ucl=upperci)
}



##data
#Glass capacitor data from Meeker & Escobar (1998), pp. 449-450
CapG <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/ZelenCap.txt", header=T)
#failure=1, censored=0
CapG$Cens <- as.numeric(CapG$Status)-1
head(CapG)



##fit an AFT model to data
#note: employ x=TRUE
mod <- survreg(Surv(Hours, Cens) ~ DegreesC + Volts,
               data=CapG, x=TRUE, weights=Weight, dist="weibull")
summary(mod)
confint(mod) #normal-approximation confidence intervals



##construct likelihood based confidence intervals for the regression coefficients
#compute the lower (lcl) and upper (ucl) confidence limits

#intercept
#note: inspect whether the curve of the log-likelihood function
#intersects the blue line
profileBeta(model=mod, whichPar=1, alpha=.95, range=c(9, 20))
#on the left, the curve does not intersect with the blue line and lcl=NA
#therefore, increase the width of the range
profileBeta(model=mod, whichPar=1, alpha=.95, range=c(8, 20))

#DegreesC
profileBeta(model=mod, whichPar=2, alpha=.95, range=c(-.06, 0))

#Volts
profileBeta(model=mod, whichPar=3, alpha=.95, range=c(-.01, 0))


##compute the likelihood based confidence intervals using profileModel
profLogLik <- function(restrFit) {-2*restrFit$loglik[2]}

prof.mod <- profileModel(mod, quantile=qchisq(0.95,1), objective=profLogLik,
                         stdErrors=summary(mod)$table[,2], gridsize=30)

(confints <- profConfint(prof.mod, method="zoom", endpoint.tolerance=1e-8))
plot(prof.mod, signed=FALSE)





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

#transform predictor variables
Cap$ArrTemp <- 11605/(Cap$DegreesC + 273.15) #Arrhenius temperature
Cap$logVolts <- log(Cap$Volts) #log transform of Volts



##fit an AFT model to data
#note: employ x=TRUE
mod <- survreg(Surv(Hours, Cens) ~ ArrTemp + logVolts + ArrTemp:logVolts,
               data=Cap, weights=Weight, x=TRUE, dist="weibull")
summary(mod)
confint(mod) #normal-approximation confidence intervals



##construct likelihood based confidence intervals for the regression coefficients
#compute the lower (lcl) and upper (ucl) confidence limits

#intercept
profileBeta(model=mod, whichPar=1, alpha=.95, range=c(-330, 150))

#ArrTemp
profileBeta(model=mod, whichPar=2, alpha=.95, range=c(-1.5, 15))

#logVolts
profileBeta(model=mod, whichPar=3, alpha=.95, range=c(-35, 80))

#interaction of ArrTemp and logVolts
profileBeta(model=mod, whichPar=4, alpha=.95, range=c(-3.5, .5))



##using profileModel
#note: works only for intercept, ArrTemp, and logVolts coefficients
profLogLik <- function(restrFit) {-2*restrFit$loglik[2]}

prof.mod <- profileModel(mod, quantile=qchisq(0.95,1), objective=profLogLik,
                         stdErrors=summary(mod)$table[,2], gridsize=30,
                         which=1:3)

(confints <- profConfint(prof.mod, method="zoom", endpoint.tolerance=1e-8))
plot(prof.mod, signed=FALSE)

#try the coefficient for the interaction
prof.mod <- profileModel(mod, quantile=qchisq(0.95,1), objective=profLogLik,
                         stdErrors=summary(mod)$table[,2], gridsize=30,
                         which=4)
#this generates an error
(confints <- profConfint(prof.mod, method="zoom", endpoint.tolerance=1e-8))





##data
#fit again an AFT model to the Glass capacitor data,
#but this time include the interaction between DegreesC and Volts

#AFT model (note: employ x=TRUE)
mod <- survreg(Surv(Hours, Cens) ~ DegreesC + Volts + DegreesC:Volts,
               data=CapG, x=TRUE, weights=Weight, dist="weibull")
summary(mod)
round(confint(mod),6) #normal-approximation confidence intervals



##construct likelihood based confidence intervals for the regression coefficients
#compute the lower (lcl) and upper (ucl) confidence limits

#intercept
#this generates an error
profileBeta(model=mod, whichPar=1, alpha=.95, range=c(-30, 50))
#but, the required likelihood based confidence intervals are computed when the
#Nelder-Mead algorithm does not use a regular simplex as initial simplex
profileBeta(model=mod, whichPar=1, alpha=.95, range=c(-30, 50), rsimplex=FALSE)

#DegreesC
profileBeta(model=mod, whichPar=2, alpha=.95, range=c(-.3, .3))

#Volts
profileBeta(model=mod, whichPar=3, alpha=.95, range=c(-.15, .15))

#interaction of DegreesC and Volts
profileBeta(model=mod, whichPar=4, alpha=.95, range=c(-.001, .001))

##note that profileModel does not work for this AFT model (try it yourself!)





##data
#IC device data from Meeker & Escobar (1998), pp. 508-511
ICd <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/ICDevice02.txt", header=T)
#using "interval2" coding approach for right censored observations
ICd$HoursU <- ifelse(ICd$Status=="Right", NA, ICd$HoursU<-ICd$HoursU)

#transform predictor variable
ICd$ArrTemp <- 11605/(ICd$DegreesC + 273.15) #Arrhenius temperature
head(ICd)



##fit an AFT model to the data
#notes:
#1) employ x=TRUE
#2) use "interval2" coding approach
mod <- survreg(Surv(HoursL, HoursU, type="interval2") ~ ArrTemp, data=ICd,
               weights=Count, x=TRUE, dist="lognormal")
mod
confint(mod) #normal-approximation confidence intervals



##construct likelihood based confidence interval for the regression coefficients
#compute the lower (lcl) and upper (ucl) confidence limits

#intercept
profileBeta(model=mod, whichPar=1, alpha=.95, range=c(-15, -6))

#ArrTemp
profileBeta(model=mod, whichPar=2, alpha=.95, range=c(.6, 1))



##using profileModel
profLogLik <- function(restrFit) {-2*restrFit$loglik[2]}

prof.mod <- profileModel(mod, quantile=qchisq(0.95,1), objective=profLogLik,
                         stdErrors=summary(mod)$table[,2], gridsize=30)

(confints <- profConfint(prof.mod, method="zoom", endpoint.tolerance=1e-8))
plot(prof.mod, signed=FALSE)