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

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

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
#failure=1, censored=0
CapG\$Cens <- as.numeric(CapG\$Status)-1

##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
#failure=1, censored=0
Cap\$Cens <- as.numeric(Cap\$Status)-1

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

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