R code for fitting an Accelerated Life Test Model

The following code may be used for fitting an accelerated life test model in R.

Currently the R code implements the lognormal and Weibull distribution for failure times. However, the code can be easily adapted to implement other distributions as well (such as the Gumbel distribution).

accelerated life test model

As an example, temperature (on the Arrhenius scale) is used as acceleration factor in the code. But it also possible to use other acceleration factors (for instance, voltage stress). Moreover, the code can be extended to include multiple accelerators in the accelerated life test model, and even their interaction.

An excellent introduction to accelerated life test models is given by Meeker and Escobar in their 1998 book Statistical methods for reliability data.

Note that the calculation of confidence intervals for failure probabilities is based on Nelson’s method. This method is described in this paper by Hong, Meeker, and Escobar.

See this blog post for computing likelihood based confidence intervals for the failure probabilities and life time quantiles of the Accelerated Life Test model.

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)

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



##examine scatter plot of failure times versus accelarating variable
plot(DevA$Temp, DevA$Hours, pch=as.numeric(DevA$Cens),
     col=as.numeric(DevA$Cens)+3, log="y",
     xlab="Degrees Celsius", ylab="Hours")
legend("bottomleft", pch=c(1,0), col=c(4,3), legend=c("failure","censored"))



##multiple probability plot
#see Meeker & Escobar, Figure 19.3, p. 497

#function for calculating probability plotting positions
ProbPos <- function (times, event, weight=rep(1, length(times)), dist){
  timesSorted <- times[order(times)]
  eventSorted <- event[order(times)]
  timesi <- unique(timesSorted[which(eventSorted==1)])
  
  cumSurvRaw <- survfit(Surv(times, event)~ 1, weights=weight)
  cumSurv <- unique(rep(cumSurvRaw$surv, cumSurvRaw$n.event))
  cumFail <- 1-cumSurv
  lagFail <- c(0, cumFail[-length(cumFail)])
  Prob <- .5*(cumFail+lagFail)
  if (dist=="lognormal") {yi <- qnorm(Prob)}
  else if (dist=="weibull") {yi <- qsev(Prob)}
  cbind(timesi, yi, Prob)
}

#calculate probability plotting positions
#assume lognormal distribution for failure times
deg40 <- DevA[which(DevA$Temp==40),]
deg40pp <- ProbPos(deg40$Hours, deg40$Cens, weight=deg40$Weight, dist="lognormal")
deg60 <- DevA[which(DevA$Temp==60),]
deg60pp <- ProbPos(deg60$Hours, deg60$Cens, weight=deg60$Weight, dist="lognormal")
deg80 <- DevA[which(DevA$Temp==80),]
deg80pp <- ProbPos(deg80$Hours, deg80$Cens, weight=deg80$Weight, dist="lognormal")

#multiple probability plot

#determine plot ranges
minProb <- min(rbind(deg40pp, deg60pp, deg80pp)[,2])
maxProb <- max(rbind(deg40pp, deg60pp, deg80pp)[,2])
minTime <- min(rbind(deg40pp, deg60pp, deg80pp)[,1])
maxTime <- max(rbind(deg40pp, deg60pp, deg80pp)[,1])
#construct plot
plot(0, type="n", log="x", xlim=c(minTime, maxTime), ylim=c(minProb, maxProb),
     xlab="Time", ylab="Linearized CDF")
points(deg40pp[,1], deg40pp[,2], col="red")
points(deg60pp[,1], deg60pp[,2], col="blue")
points(deg80pp[,1], deg80pp[,2], col="darkgreen")
legend("topleft", pch=1, col=c("red", "blue", "darkgreen"),
       legend=c("40 Deg C", "60 Deg C", "80 Deg C"))



##adding ML estimates to multiple probability plot
#see Meeker & Escobar, Figure 19.3, p. 497

#MLEs for mu and sigma (lognormal distribution) at each temperature
mleDeg40 <- survreg(Surv(Hours,Cens)~1, data=deg40, weights=Weight, dist="lognormal")
mleDeg60 <- survreg(Surv(Hours,Cens)~1, data=deg60, weights=Weight, dist="lognormal")
mleDeg80 <- survreg(Surv(Hours,Cens)~1, data=deg80, weights=Weight, dist="lognormal")

#function for predicting quantiles
predictTime <- function (p, mu, sigma, dist) {
  if (dist=="lognormal") {exp(qnorm(p)*sigma + mu)}
  else if (dist=="weibull") {exp(qsev(p)*sigma + mu)}
}

#generate predictions (lognormal distribution)
ps <- seq(0,1,length.out=200)
plotposps <- qnorm(ps) #caution: in case of Weibull distribution use qsev(ps)
predTime40deg <- sapply(ps, predictTime, mu=coef(mleDeg40),
                        sigma=mleDeg40$scale, dist="lognormal")
predTime60deg <- sapply(ps, predictTime, mu=coef(mleDeg60),
                        sigma=mleDeg60$scale, dist="lognormal")
predTime80deg <- sapply(ps, predictTime, mu=coef(mleDeg80),
                        sigma=mleDeg80$scale, dist="lognormal")

#include ML estimates in probability plot, and assess
#whether the slopes of the lines are similar (=assumption of
#the accelerated life test model)
plot(0, type="n", log="x", xlim=c(minTime, maxTime), ylim=c(minProb, maxProb),
     xlab="Time", ylab="Linearized CDF")
points(deg40pp[,1], deg40pp[,2], col="red")
lines(predTime40deg, plotposps, col="red")
points(deg60pp[,1], deg60pp[,2], col="blue")
lines(predTime60deg, plotposps, col="blue")
points(deg80pp[,1], deg80pp[,2], col="darkgreen")
lines(predTime80deg, plotposps, col="darkgreen")
legend("topleft", pch=1, col=c("red", "blue", "darkgreen"),
       legend=c("40 Deg C", "60 Deg C", "80 Deg C"))



##MLEs and their confidence intervals at each temperature
#see Meeker & Escobar, Table 19.1, p. 498

#40 degrees Celsius
#MLEs for mu and sigma
summary(mleDeg40)
#95% confidence interval for mu
confint(mleDeg40, parm=1, level=.95)
#95% confidence interval for sigma
SEparms40deg <- sqrt(diag(vcov(mleDeg40))) #extract standard errors
exp(log(mleDeg40$scale) + c(1,-1)*qnorm(.05/2)*SEparms40deg[2])

#60 degrees Celsius
#MLEs for mu and sigma
summary(mleDeg60)
#95% confidence interval for mu
confint(mleDeg60, parm=1, level=.95)
#95% confidence interval for sigma
SEparms60deg <- sqrt(diag(vcov(mleDeg60))) #extract standard errors
exp(log(mleDeg60$scale) + c(1,-1)*qnorm(.05/2)*SEparms60deg[2])


#80 degrees Celsius
#MLEs for mu and sigma
summary(mleDeg80)
#95% confidence interval for mu
confint(mleDeg80, parm=1, level=.95)
#95% confidence interval for sigma
SEparms80deg <- sqrt(diag(vcov(mleDeg80))) #extract standard errors
exp(log(mleDeg80$scale) + c(1,-1)*qnorm(.05/2)*SEparms80deg[2])



##fit acceleration relationship
#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
accmod <- survreg(Surv(Hours, Cens) ~ tempArr, data=DevA, weights=Weight,
                  dist="lognormal")
summary(accmod)



##confidence intervals for parameters of Arrhenius model
#see Meeker & Escobar, Table 19.2, p. 499

#95% confidence interval for beta0 and beta1
confint(accmod, level=.95)
#95% confidence interval for sigma
SEparmsAcc <- sqrt(diag(vcov(accmod))) #extract standard errors
exp(log(accmod$scale) + c(1,-1)*qnorm(.05/2)*SEparmsAcc[3])



##construct variance covariance matrix of model parameters
#note: apply the delta method for "undoing" the log transformation of sigma
vcovAcc <- vcov(accmod)
vcovAcc[1,3] <- vcovAcc[3,1] <- t(c(0,0,exp(log(accmod$scale))))%*%vcovAcc%*%c(1,0,0)
vcovAcc[2,3] <- vcovAcc[3,2] <- t(c(0,0,exp(log(accmod$scale))))%*%vcovAcc%*%c(0,1,0)
vcovAcc[3,3] <- t(c(0,0,exp(log(accmod$scale))))%*%vcovAcc%*%c(0,0,exp(log(accmod$scale)))
rownames(vcovAcc) <- colnames(vcovAcc) <- c("beta0", "beta1", "sigma")
#variance covariance matrix (estimates are identical to those of Meeker & Escobar, p. 499)
vcovAcc



##checking model assumptions (diagnostics)

##standardized residuals
linFit <- predict(accmod, type="lp")
sderr <- (log(DevA$Hours)-linFit)/accmod$scale
#lognormal distribution
sderrpp <- ProbPos(sderr, DevA$Cens, weight=DevA$Weight, dist="lognormal")
#probability plot of standardized residuals
plot(sderrpp[,1], sderrpp[,2],
     xlab="Standardized Residuals", ylab="Linearized CDF")
#visual aid for detecting linearity
abline(a=0, b=1, lty=2, col="blue")



##Martingale residuals
#for more information on Martingale residuals
#see: http://blogs2.datall-analyse.nl/2016/02/19/rcode_martingale_residuals_aft_model/

#null model for lognormal distribution
sr <- survreg(Surv(Hours, Cens) ~ 1, data=DevA, weights=Weight,
              dist="lognormal")

#standardized residuals
linFit <- predict(sr, type="lp")
sderr <- (log(DevA$Hours)-linFit)/sr$scale

#Cox-Snell residuals
CoxSnellResidual <- function (standRes, weight=1, dist){
  standRes <- standRes[rep(seq_len(length(standRes)), weight)]
  if (dist=="lognormal") {csr <- -log(1-pnorm(standRes))}
  else if (dist=="weibull") {csr <- -log(exp(-exp(standRes)))}
}

cxsn <- CoxSnellResidual(standRes=sderr, weight=DevA$Weight, dist="lognormal")

#Martingale residuals
MartingaleResidual <- function (CoxSnellResid, event, weight=1) {
  delta <- event[rep(seq_len(length(event)), weight)]
  martingale <- delta-CoxSnellResid
  data.frame(Martingale=martingale, Event=delta)
}

mgale <- MartingaleResidual(CoxSnellResid=cxsn, event=DevA$Cens, weight=DevA$Weight)

#Martingale residuals plot detects an approximately linear functional form
#for Arrhenius temperature, but has (as usual) the wrong (opposite) sign
#that is, the slope of the plot is negative, whereas the
#coefficient for Arrhenius temperature in the lognormal model is positive
tempArrExpanded <- DevA$tempArr[rep(seq_len(length(DevA$tempArr)), DevA$Weight)]
plot(tempArrExpanded, mgale$Martingale, pch=mgale$Event,
     xlab="Arrhenius temperature", ylab="Martingale Residual")
abline(h=0, lty=2, col="gray40")
lines(lowess(tempArrExpanded, mgale$Martingale, iter=0), col="blue")
legend("topright", pch=c(1,0), legend=c("failure","censored"))



##predictions of lognormal model
#comparable to Meeker & Escobar, Figure 19.4, p. 499

#failure plot

#generate sequence of proportions
ps <- seq(0, 1, length.out=200)
predMu <- unique(predict(accmod, type="lp"))
#predictions of lognormal model
predTime <- sapply(ps, predictTime, mu=predMu, sigma=accmod$scale,
                   dist="lognormal")

#plot predictions
matplot(x=t(predTime), y=ps,
        type="l",
        xlab="Hours", ylab="Proportion Failing",
        lty=c(1,1,1,1), lwd=c(1,1,1,1),
        col=c("cyan", "red", "blue", "darkgreen"),
        xlim=c(0,50000), ylim=c(0,1))
#include observed failure times
points(deg40pp[,1], deg40pp[,3], col="red")
points(deg60pp[,1], deg60pp[,3], col="blue")
points(deg80pp[,1], deg80pp[,3], col="darkgreen")
legend("topright", pch=1, col=c("cyan", "red", "blue", "darkgreen"),
       legend=c("10 Deg C", "40 Deg C", "60 Deg C", "80 Deg C"))
#draw raster
abline(h = seq( 0, 1, .05 ), lty=3, col="gray")
abline(v = seq( 0, 50000, 5000), lty=3, col="gray")



#for predicting failure probabilities and life time quantiles
#at a specified temperature, it is necessary to compute
#the variance covariance for mu and sigma at this temperature

#specify temperature
Temp <- 10
#transform to Arrhenius scale
ArrTemp <- 11605 / (Temp + 273.15)

#function for computing variance covariance matrix for mu and sigma
#at the specified temperature
vcovMuSigma <- function (vcov, newdata){
  #variance(mu):
  #partial derivatives of the function g=beta0+beta1*x with respect to
  #beta0, beta1, and sigma, respectively
  dg.dtheta <- c(1, newdata, 0)
  #apply delta method
  varMu <- t(dg.dtheta)%*%vcov%*%dg.dtheta
  
  #covariance(mu,sigma):
  #partial derivatives of the function g1=beta0+beta1*x with respect to
  #beta0, beta1, and sigma, respectively
  dg1.dtheta <- c(1, newdata, 0)
  #partial derivatives of the function g2=sigma with respect to
  #beta0, beta1, and sigma, respectively
  dg2.dtheta <- c(0, 0, 1)
  #apply delta method
  covMuSigma <- t(dg1.dtheta)%*%vcov%*%dg2.dtheta
  
  #construct variance covariance matrix
  vcovMS <- matrix(nrow=2, ncol=2)
  vcovMS[1,1] <- varMu
  vcovMS[1,2] <- vcovMS[2,1] <- covMuSigma
  vcovMS[2,2] <- vcov[3,3]
  rownames(vcovMS) <- colnames(vcovMS) <- c("mu", "sigma")
  vcovMS
}
#variance covariance matrix (estimates are identical to those of Meeker & Escobar, p. 503)
(vcMS <- vcovMuSigma(vcov=vcovAcc, newdata=ArrTemp))



##predict failure probabilities at specified temperature
##and calculate the confidence intervals

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

#Nelson's method for computing confidence intervals
#function for computing standard error for failure probability
seZhat <- function (vcov, mu, sigma, te, newdata, dist){
  if (dist=="lognormal"){
    Ze <- qnorm(predFT(mu, sigma, te, newdata, dist))}
  else if (dist=="weibull"){
    Ze <- qsev(predFT(mu, sigma, te, newdata, dist))}
  #partial derivative of the function z=(log(t)-mu)/sigma
  #with respect to mu
  dz.dmu <- -1/sigma
  #partial derivative of the function z=(log(t)-mu)/sigma
  #with respect to sigma
  dz.dsigma <- (-1/sigma)*Ze
  #compute variance covariance matrix at specified temperature
  vcMS <- vcovMuSigma(vcov, newdata)
  #delta method
  seZ <- t(c(dz.dmu, dz.dsigma))%*%vcMS%*%c(dz.dmu, dz.dsigma)
  as.vector(sqrt(seZ))
}

#function for computing confidence intervals
confintFtz <- function (vcov, mu, sigma, te, newdata, alpha, dist){
  Fe <- predFT(mu, sigma, te, newdata, dist)
  seZ <- seZhat(vcov, mu, sigma, te, newdata, dist)
  w <- qnorm(1-alpha/2)*seZ
  if (dist=="lognormal"){
    Ze <- qnorm(Fe); lcl <- (pnorm(Ze-w)); ucl <- (pnorm(Ze+w))}
  else if (dist=="weibull"){
    Ze <- qsev(Fe); lcl <- (psev(Ze-w)); ucl <- (psev(Ze+w))}
  c(Te=te, Fhat=Fe, std.err=seZ, lcl=lcl, ucl=ucl)
}

#confidence interval for F(10000) and F(30000) at 10 deg Celsius
tes <- c(10000, 30000) #predict failure probability at 10000 and 30000 hours
t(sapply(tes, confintFtz, vcov=vcovAcc, newdata=ArrTemp, alpha=.05,
         dist="lognormal", mu=coef(accmod), sigma=accmod$scale))

#failure plot at 10 deg Celsius (including confidence intervals)
#comparable to Meeker & Escobar, Figure 19.4, p. 499
tes <- seq(1000, 2000000, length.out=200)
Fes <- data.frame(t(sapply(tes, confintFtz, vcov=vcovAcc,
                           newdata=ArrTemp, alpha=.05, dist="lognormal",
                           mu=coef(accmod), sigma=accmod$scale)))
plot(Fes$Te, Fes$Fhat, type="l", col="blue", xlab="Time",
     ylab="Proportion Failing", ylim=c(0,1), xlim=c(0,2000000))
lines(Fes$Te, Fes$lcl, lty=2, col="red")
lines(Fes$Te, Fes$ucl, lty=2, col="red")



##predict life time quantiles at specified temperature
##and calculate the confidence intervals

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

#function for computing standard error for life time quantile
seQphat <- function (vcov, mu, sigma, p, newdata, dist) {
  Qp <- predQp(mu, sigma, p, newdata, dist)
  dg.dmu <- 1
  if (dist=="lognormal") {dg.dsigma <- qnorm(p)}
  else if (dist=="weibull") {dg.dsigma <- qsev(p)}
  vcMS <- vcovMuSigma(vcov, newdata)
  seQp <- t(c(dg.dmu, dg.dsigma))%*%vcMS%*%c(dg.dmu, dg.dsigma)
  as.vector(Qp*sqrt(seQp))
}

#function for computing confidence intervals
confintQp <- function (vcov, mu, sigma, p, newdata, alpha, dist){
  Qp <- predQp(mu, sigma, p, newdata, dist)
  seQp <- seQphat(vcov, mu, sigma, p, newdata, dist)
  w <- exp((qnorm(1-alpha/2)*seQp) / Qp)
  c(p=p, Quantile=Qp, std.err=seQp, lcl=Qp/w, ucl=Qp*w )
}

#confidence interval for p=.2 and p=.7 at 10 deg Celsius
pes <- c(.2, .7) #predict life time quantile at p=.2 and p=.7
t(sapply(pes, confintQp, vcov=vcovAcc, newdata=ArrTemp, alpha=.05,
         mu=coef(accmod), sigma=accmod$scale, dist="lognormal"))

#failure plot at 10 deg Celsius (including confidence intervals)
#comparable to Meeker & Escobar, Figure 19.4, p. 499
pes <- seq(.001, 1, length.out=200)
Qes <- data.frame(t(sapply(pes, confintQp, vcov=vcovAcc,
                           newdata=ArrTemp, alpha=.05, dist="lognormal",
                           mu=coef(accmod), sigma=accmod$scale)))
plot(Qes$Quantile, Qes$p, type="l", col="blue", xlab="Time",
     ylab="Proportion Failing", ylim=c(0,1), xlim=c(0,2000000))
lines(Qes$lcl, Qes$p, lty=2, col="red")
lines(Qes$ucl, Qes$p, lty=2, col="red")