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).
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.datall-analyse.nl/blog_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.datall-analyse.nl/blog_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.datall-analyse.nl/blog_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)