R code for fitting an Accelerated Degradation Model

This R code may be used for fitting the Accelerated Degradation Model described by Meeker and Escobar in their 1998 book Statistical Methods for Reliability Data (on pages 564-567).
The same degradation model is also described in this paper by Meeker, Escobar, and Lu.

The degradation model was fitted in R using the nlme package (version 3.1-118).

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

##Accelerated Degradation Model
#Meeker & Escobar, pp. 564-567


#accelerated degradation test data of Device-B
devB <- read.table("http://www.datall-analyse.nl/blog_data/DeviceB.txt", header=T)

#grouped data format
devBG <- groupedData(PowerDrop ~ Hours | Device, data=devB)

#inspect data (Meeker & Escobar, p. 564 , figure 21.1)
plot(devBG, outer= ~1)
plot(devBG, outer= ~ DegreesC, aspect=2)

#temperature differential factor (Meeker & Escobar, p. 472)
devBG$TDF <- (11605/(195+273.15))-(11605/(devBG$DegreesC+273.15))

##nonlinear mixed effects model

#obtain start values with nlsList (ignore error)
fmdevB <- nlsList(SSasympOrig, data=devBG)

#obtain start values for fixed effects in nlme model
nlmeDevB <- nlme(fmdevB)

#check relationship between "rate of reaction" (lrc) and
#the covariate "TDF" (temperature differential factor)
fmRE <- ranef(nlmeDevB, aug=T)
plot(fmRE, form = lrc ~ TDF)

#rate of reaction (lrc) is related to TDF
#thus, include TDF as covariate in the model
nlmeDevBCov <- update(nlmeDevB, fixed = list(Asym ~ 1, lrc ~ TDF), start = c(nlmeDevB$coefficients$fixed[1], nlmeDevB$coefficients$fixed[2], 0))

#final model
#note: the model's estimates are similar to those reported by Meeker & Escobar on p. 567
#also note that their estimate of the asymptote is -1*exp(.3510)=-1.420, which is close
#to the model's estimate for "Asym"
#applying the delta method, Meeker & Escobar's estimate of the "Asym" variance is
#Meeker & Escobar's estimate of the covariance between "Asym" and "lrc" is
#-exp(.3510)*1*-.02918=0.04145 (applying the delta method), and the
#model's estimate for this covariance is:

#confidence intervals
intervals(nlmeDevBCov, which="all")

##diagnostics final model

#covariance structure of random effects
#outliers are marked by numbers in the plot
pairs(nlmeDevBCov, ~ranef(.), id=.10)

#normality of random effects
qqnorm(nlmeDevBCov, ~ ranef(., standard=T), abline=c(0,1))

#testing for the fitted bivariate normal distribution of the random effects
#(comparable to Meeker & Escobar, figure 21.3, p. 568)
#extract variance-covariance matrix of fitted random effects
VarCovNLME <- matrix(nrow=2, ncol=2)
VarCovNLME[1,1] <- as.numeric(VarCorr(nlmeDevBCov)[1,1])
VarCovNLME[2,2] <- as.numeric(VarCorr(nlmeDevBCov)[2,1])
VarCovNLME[1,2] <- VarCovNLME[2,1] <- as.numeric(VarCorr(nlmeDevBCov)[2,3])*
#contour representing the estimated bivariate normal distribution
plot(data.frame(ranef(nlmeDevBCov)), xlim=c(-.5,.5), ylim=c(-1,1))
ellipse(mu=c(0,0), sigma=VarCovNLME, alpha=.05, npoints=250, col="red")

#standardized residuals
plot(nlmeDevBCov, Device~resid(., type="p"), abline=0)
plot(nlmeDevBCov, resid(., type="p")~fitted(.), abline=0)
qqnorm(nlmeDevBCov, ~resid(., type="p"), abline=c(0,1))

#within-group autocorrelation of standardized residuals
plot(ACF(nlmeDevBCov, maxLag=10, resType="p"), alpha=.01)

#predictions model
plot(augPred(nlmeDevBCov, level=0:1))