The R package ALTopt has been developed with the aim of creating and evaluating optimal experimental designs of censored accelerated life tests (ALTs). This package takes the generalized linear model approach to ALT planning, because this approach can easily handle censoring plans and derive information matrices for evaluating designs. Three types of optimality criteria are considered: D-optimality for model parameter estimation, U-optimality for reliability prediction at a single use condition, and I-optimality for reliability prediction over a region of use conditions. The Weibull distribution is assumed for failure time data and more than one stress factor can be specified in the package. Several graphical evaluation tools are also provided for the comparison of different ALT test plans.
Accelerated life testing (ALT) is commonly used for obtaining a product’s failure time data by subjecting it to elevated stress conditions, such as temperature, humidity, and voltage. As a result, the product fails in a shorter time period than would be expected under normal stress conditions. The failure data obtained from ALTs can then be extrapolated to the normal use stress level to estimate the product’s lifetime distribution. Nelson (2005a,2005b) provides a comprehensive review of ALT papers up to 2005.
To avoid poor experimental results and to obtain more accurate inference on the acceleration model and on reliability prediction, it is necessary to have an effective ALT test plan. A well-designed ALT test plan often aims to achieve some statistical optimality. However, conventional experimental designs (e.g., factorial designs) are not effective as ALT test plans because of the following features of ALTs:
Extrapolation – Test stress levels are typically higher than the normal use stress levels. As failure time data will be collected at these higher stress levels, extrapolating them to the normal use stress level is needed for reliability prediction. Nonlinear relationships between failure time and stress levels are expected.
Non-normal distributions of failure times – The failure time distribution is typically positively skewed, e.g., the Weibull distribution.
Censoring of failure time data – Censoring occurs when the exact failure times of test units are not observed. There are several reasons for censoring. In some cases, test units do not fail by the end of the test period, in which case the data becomes right-censored. In other cases, test units are periodically inspected, so the only information available is the time interval of failure, while the exact failure time is unknown. The latter case is called interval-censoring.
In this article, we introduce an R package, ALTopt (Seo and Pan 2015), that constructs optimal test plans for ALTs with right- and interval-censored data. This package is based on the work done by Monroe et al. (2011) and Yang and Pan (2013), where generalized linear models (GLMs) were used to model censored ALT data.
The ALTopt package accommodates three optimality criteria: D-optimality, U-optimality and I-optimality. A D-optimal design minimizes the generalized variance of parameter estimates, while a U-optimal or I-optimal design minimizes, respectively, the (average) variance of response prediction at a single use condition or over a region of use conditions.
In the GLM context, the D-optimal design is defined by
Using the same notation, the U-optimal and I-optimal designs can be
defined, respectively, as
A function that links failure time and stress variables is needed in order to extrapolate the results obtained in the test region to the use region. The GLM formulation for ALT is built upon the Cox’s proportional hazard (PH) assumption. This section provides the derivation of these formulations for right-censored and interval-censored ALT data.
The PH model assumes that, given the vector of explanatory variables
where
Note that the baseline hazard function is a function of time only. From Eq. (1) we can derive that
where
where
The baseline hazard function of a Weibull distribution is given by
With the proportional hazard assumption, the failure density function is
given by
For a failure time data set that includes right-censored survival times,
each observation can be expressed as a pair
From Eqs. (1) and (3) the log-likelihood function can be
written as
Note that the last two terms of sum on the right-hand side of Eq.
(4) are the same as the kernel of the log-likelihood function of
The response variables,
The linear predictor is
The link function is given by
This offset term in the link function is
Since the log link function is the canonical link function for the
Poisson distribution, the asymptotic variance-covariance matrix of
Monroe et al. (2011) has shown that
where
For a Weibull distribution with the known shape parameter
where
For interval-censored data, the whole testing period is divided into
multiple time intervals such as
and the conditional probability of surviving at the beginning of the
It can be shown that
Define an indicator variable,
The likelihood function can be expressed as
which is equivalent to
where
The likelihood function of Eq. (5) has the same form as the
likelihood function of independent binomial random variables. We treat
Now, we examine the probability
By Eq. (3) it becomes
The second term of the right hand side of Eq. (7) does not
depend on the regression coefficient
The response variables,
The linear predictor is
The link function is given by
Since the log-log link is not a canonical link for the binomial
distribution, we need to introduce
Then, the asymptotic variance-covariance matrix of
where
In a binomial distribution,
Assuming a Weibull distribution for a product’s lifetime, we have
Assume all time intervals have the same length,
We also have
Substituting Eq. (9) and Eq. (10) into Eq. (8) yields the weight matrix for the interval-censored Weibull failure time data:
The main purpose of the ALTopt package is to construct D-, U-, and
I-optimal ALT test plans. Two main functions, altopt.rc
and
altopt.ic
, are developed respectively for the right-censoring and
interval-censoring cases. The following assumptions are required for
using this package:
Failure time data follows the Weibull distribution where the shape parameter is specified by the user.
Log-linear functions are used to model the relationship between a failure time distribution parameter and the stress factors.
For interval-censored data, all the intervals have the same length.
Lastly, the package can accommodate many stress factors, but 5 or fewer is recommended for computational efficiency.
This package also provides two functions for evaluating existing test
plans – alteval.rc
and alteval.ic
. These functions can be used for
comparing test plans generated by this ALTopt package or any other
methods. Graphical displays of prediction variance are made available.
These plotting features enhance the usefulness of this package for
comparing and selecting test plans.
The syntax of altopt.rc
and altopt.ic
functions are as follows:
altopt.rc(optType, N, tc, nf, alpha, formula, coef,
useCond, useLower, useUpper, nOpt = 1, nKM = 30, nCls = NULL)
altopt.ic(optType, N, t, k, nf, alpha, formula, coef,
useCond, useLower, useUpper, nOpt = 1, nKM = 30, nCls = NULL)
The arguments within these functions are:
optType
– Choice of "D"
, "U"
, and "I"
optimality.
N
– The number of test units.
tc
(altopt.rc
only) – The planned right censoring time.
t
(altopt.ic
only) – The planned total testing time (i.e., the
end point of the last interval).
k
(altopt.ic
only) – The number of time intervals.
nf
– The number of stress factors.
alpha
– The value of the shape parameter of the Weibull
distribution.
formula
– The object of class “formula
” expressing the linear
predictor model.
coef
– The numeric vector containing the coefficients of each
term in formula
.
useCond
– The numeric vector of use condition. It should be
provided when optType
is "U"
. The length of the vector should be
same as the number of stress factors.
useLower
– The numeric vector of lower bound of use region in
coded units. It should be provided when optType
is "I"
. The
length of the vector should be the same as the number of stress
factors.
useUpper
– The numeric vector of upper bound of use region in
coded units. It should be provided when optType
is "I"
. The
length of the vector should be the same as the number of stress
factors.
nOpt
– The number of repetitions of optimization processes. The
default value is 1. If nOpt
is larger than 1, each optimization
process starts from randomly chosen design points in the design
region and each solution may slightly differ. The output shows the
best solution overall.
nKM
– The number of repetitions of k-means clustering, which is
used to generate the optimal design clustered by kmeans
. The
default value is 30.
nCls
– The number of clusters used for k-means clustering. If not
specified, it is set as the number of parameters in the linear
predictor model.
We use the function stats::optim
with the "L-BFGS-B"
method to
perform optimization. This function allows box constraints on design
variables. In our case, we have a cuboidal design region where the
levels of each stress factor are coded to be between 0 and 1. More
details about the "L-BFGS-B"
method are available in Byrd et al. (1995).
The output of these functions are given as a list with the following components:
call
– The matched call.
opt.design.rounded
– The optimal design clustered by rounding.
opt.value.rounded
– The objective function value of
opt.design.rounded
.
opt.design.kmeans
– The optimal design clustered by k-means
algorithm.
opt.value.kmeans
– the objective function value of
opt.design.kmeans
.
The procedure begins by generating an initial test plan with nCls
. By carefully selecting
the number of clusters, it is possible to reduce the number of distinct
design points without significantly affecting the value of the objective
function. The final recommended test plans are provided by the elements,
opt.design.rounded
, or opt.design.kmeans
, presented by a table
containing each design point location and the number of test units
allocated at each design point. The corresponding values of the
objective function of these plans are also stored in opt.value.rounded
and opt.value.kmeans
.
The design.plot
function displays the recommended test plan as a
bubble plot on a two-dimensional design region of user-specified stress
factors. The size of each bubble represents the relative size of the
test unit allocations. The arguments of design.plot
are as follows:
design.plot(design, xAxis, yAxis)
design
– A data frame containing the coordinates and the test
unit allocation at each design point. The components,
opt.design.rounded
or opt.design.kmeans
, of an output created by
altopt.rc
or altopt.ic
can be given for this argument directly
and any other design with the same form of those can also be given.
xAxis
– The name of the factor to be displayed on the
yAxis
– The name of the factor to be displayed on the
This package provides several methods to evaluate an ALT test plan. The
first method is the numerical evaluation of a given test plan using
alteval.rc
or alteval.ic
. These functions return the value of the
objective function of the test plan. The arguments are as follows:
alteval.rc(designTable, optType, tc, nf, alpha, formula, coef,
useCond, useLower, useUpper)
alteval.ic(designTable, optType, t, k, nf, alpha, formula, coef,
useCond, useLower, useUpper)
The existing test plan is specified in the argument designTable
. The
other arguments of alteval.rc
and alteval.ic
are similar to the
arguments in altopt.rc
and altopt.ic
.
ALTopt also provides three different graphs for evaluating a test plan
– the prediction variance (PV) contour plot, the fraction of use space
(FUS) plot, and the variance dispersion of use space (VDUS) plot. These
graphical tools are useful when visualizing the prediction variance
throughout the entire use-space region (Myers et al. 2009 8). The
PV contour plot displays the contours of the estimated prediction
variance from the design region to the use region of a two-dimensional
user-specified stress factor space. Functions pv.contour.rc
and
pv.contour.ic
generate the PV contour plot of an ALT test plan with
right and interval censoring, respectively. The FUS plot is an extension
of the fraction of design space (FDS) proposed by Zahran et al. (2003).
The vertical axis of a FUS plot is the fraction of the use space region
that has prediction variance less than or equal to the given values in
the horizontal axis. Functions pv.fus.rc
and pv.fus.ic
create the
FUS plot of right and interval censoring ALT plans, respectively. In
addition, the FUS curves of multiple designs can be overlaid on one
graph by using compare.fus
, so these designs can be compared
graphically. The VDUS plot is an extension of the variance dispersion
graphs (VDGs) of Giovannitti-Jensen and Myers (1989) to the cuboidal use space
region. It shows plots of minimum, average and maximum prediction
variance from the center to the boundary of the use region. The
comparison of multiple VDUS is also available through compare.vdus
.
The arguments of these functions are omitted here, because they are
similar to previously described functions.
In this section, we demonstrate the use of ALTopt using the
right-censored ALT data set from Yang and Pan (2013). In this experiment, an
ALT of 100 test units is conducted with two stress factors –
temperature and humidity. The lowest and highest stress levels in the
test region are (60 , 60 %) and (110 , 90 %), respectively. The normal
use condition is (30 , 25 %), while the typical use region has the range
from (20 , 20 %) to (40 , 30 %). The natural stress variables of these
two factors are defined by
R> NuseCond <- c(11605 / (30 + 273.15), log(25))
R> NuseLow <- c(11605 / (20 + 273.15), log(20))
R> NuseHigh <- c(11605 / (40 + 273.15), log(30))
R> NdesLow <- c(11605 / (60 + 273.15), log(60))
R> NdesHigh <- c(11605 / (110 + 273.15), log(90))
Next, we apply a coding scheme on these natural variables so that the highest stress level becomes (0, 0) and the lowest stress level becomes (1, 1).
Here, convert.stress.level
, to convert the natural stress level to the coded
stress level, and vice versa. The use condition and the use stress
region are accordingly coded as follows:
R> library(ALTopt)
R> (useCond <- as.numeric(convert.stress.level(NdesLow, NdesHigh,
+ actual = NuseCond))) # Coded use condition
[1] 1.758337 3.159172
R> (useLower <- as.numeric(convert.stress.level(NdesLow, NdesHigh,
+ actual = NuseHigh))) # Coded use region's lower bound
[1] 1.489414 2.709511
R> (useUpper <- as.numeric(convert.stress.level(NdesLow, NdesHigh,
+ actual = NuseLow))) # Coded use region's upper bound
[1] 2.045608 3.709511
We assume that the failure times follow an exponential distribution,
i.e., alpha = 1
, and the pre-specified linear predictor is given by
Suppose the total testing time is 30 time units. The D-optimal test plan is generated by the following lines of code:
R> set.seed(10)
R> DR <- altopt.rc("D", N = 100, tc = 30, nf = 2, alpha = 1,
+ formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01))
R> DR
$call
altopt.rc(optType = "D", N = 100, tc = 30, nf = 2, alpha = 1,
formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01))
$opt.design.rounded
x1 x2 allocation
1 0.000 0 21
2 0.835 0 28
3 0.000 1 26
4 0.639 1 25
$opt.value.rounded
[1] 27153.91
$opt.design.kmeans
x1 x2 allocation
1 0.8353075 0 28
2 0.0000000 0 21
3 0.6390136 1 25
4 0.0000000 1 26
$opt.value.kmeans
[1] 27153.92
While the formula does not include the intercept term explicitly, the value of the intercept parameter still needs to be specified (in this case, it is 0). From the final design output, we noticed that the designs generated by rounding and clustering are almost the same.
We can also generate the U-optimal and I-optimal designs using the following lines of code:
R> set.seed(50)
R> UR <- altopt.rc("U", N = 100, tc = 30, nf = 2, alpha = 1,
+ formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01),
+ useCond = useCond)
R> set.seed(100)
R> IR <- altopt.rc("I", N = 100, tc = 30, nf = 2, alpha = 1,
+ formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01),
+ useLower = useLower, useUpper = useUpper)
Using design.plot
function, we can draw the bubble plots of these test
plans.
R> design.plot(DR$opt.design.rounded, xAxis = x1, yAxis = x2)
R> design.plot(UR$opt.design.rounded, xAxis = x1, yAxis = x2)
R> design.plot(IR$opt.design.rounded, xAxis = x1, yAxis = x2)
|
|
|
(a) | (b) | (c) |
design.plot
function.
From Figure 1, one can see that the U- and I-optimal test plans resemble each other, while the D-optimal test plan is very different from the other two. This is expected because the objective functions of U- and I-optimal designs involve the variance of reliability prediction, while the D-optimal design involves the variance of parameter estimation. From the U- and I-optimal test plans, it is noticeable that a large number of test units is allocated at the lowest stress level. This type of test unit allocation scheme is common in ALTs (e.g., Meeker and Nelson 1975).
To compare the D-optimal test plan and the U-optimal test plan, the
pv.contour.rc
function generates the contour plot of prediction
variance using the following lines of code:
R> pv.contour.rc(DR$opt.design.rounded, xAxis = x1, yAxis = x2,
+ tc = 30, nf = 2, alpha = 1, formula = ~ x1 + x2 + x1:x2,
+ coef = c(0, -4.086, -1.476, 0.01), useCond = useCond)
R> pv.contour.rc(UR$opt.design.rounded, xAxis = x1, yAxis = x2,
+ tc = 30, nf = 2, alpha = 1, formula = ~ x1 + x2 + x1:x2,
+ coef = c(0, -4.086, -1.476, 0.01), useCond = useCond)
|
|
(a) | (b) |
pv.contour.rc
function.
|
(a) |
|
(b) |
Figure 2 shows that the U-optimal test plan has lower prediction variance than the D-optimal test plan at the normal use condition. The FUS and VDUS plots can also be used for further comparison of these test plans. These plots are shown in Figure 3.
R> fusDR <- pv.fus.rc(DR$opt.design.rounded, tc = 30, nf = 2, alpha = 1,
+ formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01),
+ useLower = useLower, useUpper = useUpper)
R> fusUR <- pv.fus.rc(UR$opt.design.rounded, tc = 30, nf = 2, alpha = 1,
+ formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01),
+ useLower = useLower, useUpper = useUpper)
R> compare.fus(fusDR, fusUR)
R> vdusDR <- pv.vdus.rc(DR$opt.design.rounded, tc = 30, nf = 2, alpha = 1,
+ formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01),
+ useLower = useLower, useUpper = useUpper)
R> vdusUR <- pv.vdus.rc(UR$opt.design.rounded, tc = 30, nf = 2, alpha = 1,
+ formula = ~ x1 + x2 + x1:x2, coef = c(0, -4.086, -1.476, 0.01),
+ useLower = useLower, useUpper = useUpper)
R> compare.vdus(vdusDR, vdusUR)
Figure 3 shows that the U-optimal test plan performs better in majority of the use-space region with respect to the prediction variance.
Finally, the convert.stress.level
function is useful for converting
the
R> convert.stress.level(NdesLow, NdesHigh, stand = UR$opt.design.rounded)
x1 x2 allocation
1 30.28840 4.499810 17
2 34.53414 4.499810 26
3 30.28840 4.094345 16
4 33.81136 4.094345 41
This paper describes the ALTopt package in R for constructing optimal ALT test plans for right- and interval-censored data. The package accommodates three statistical optimality criteria – D-optimal, U-optimal and I-optimal. It applies the GLM approach to the modeling of failure/censoring times and the derivation of the asymptotic variance-covariance matrix of regression coefficients. Failure times are assumed to follow a Weibull distribution. To use the package effectively, users are required to specify the linear predictor of the GLM and the shape parameter of the Weibull distribution. An example demonstrated the construction of optimal test plans for an ALT with two stress factors and right-censored data. This package also provides graphical functions for evaluating and comparing various test plans.
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Seo & Pan, "ALTopt: An R Package for Optimal Experimental Design of Accelerated Life Testing", The R Journal, 2015
BibTeX citation
@article{RJ-2015-029, author = {Seo, Kangwon and Pan, Rong}, title = {ALTopt: An R Package for Optimal Experimental Design of Accelerated Life Testing}, journal = {The R Journal}, year = {2015}, note = {https://rjournal.github.io/}, volume = {7}, issue = {2}, issn = {2073-4859}, pages = {177-188} }