# Experimental Designs for Accelerated Degradation Tests Based on Linear Mixed
Effects Models

Helmi Shat<EMAIL_ADDRESS>Rainer Schwabe<EMAIL_ADDRESS>Institute for
Mathematical Stochastics, Otto-von-Guericke University Magdeburg,
Universitätsplatz 2, 39106 Magdeburg, Germany

###### Abstract

Accelerated degradation tests are used to provide accurate estimation of
lifetime properties of highly reliable products within a relatively short
testing time. There data from particular tests at high levels of stress (e. g.
temperature, voltage, or vibration) are extrapolated, through a physically
meaningful model, to obtain estimates of lifetime quantiles under normal use
conditions. In this work, we consider repeated measures accelerated
degradation tests with multiple stress variables, where the degradation paths
are assumed to follow a linear mixed effects model which is quite common in
settings when repeated measures are made. We derive optimal experimental
designs for minimizing the asymptotic variance for estimating the median
failure time under normal use conditions when the time points for measurements
are either fixed in advance or are also to be optimized.

###### keywords:

Accelerated degradation test, linear mixed effects model, failure time
distribution, locally optimal design, destructive testing.

##  1 Introduction

Rapid advances in modern manufacturing technologies, along with market needs
for sustainable and reliable products, have motivated industrial companies to
design and manufacture products that can operate without failure for years or
even decades. As a consequence, manufacturers are demanded to provide their
customers with accurate information about the reliability of their products.
However, when the products get more reliable, it becomes more difficult or
even unfeasible to assess a sufficient amount of lifetime data on the basis of
traditional reliability testing in order to accurately estimate
characteristics of the lifetime distribution of the products because failure
or fatigue can hardly be observed under normal use conditions in a reasonable
time period for testing. As an alternative, for highly reliable and enduring
products, Accelerated Degradation Tests can be utilized to provide sufficient
information on the deterioration of the products to obtain a sufficiently
accurate estimate of lifetime properties within a relatively short testing
time period. In Accelerated Degradation Testing products are tested at various
elevated stress levels (for e. g. temperature, voltage, or vibration), and the
grade of degradation of the product under these stress levels is measured over
a tractable time period. The measured data are then extrapolated, through a
technically reasonable statistical model, to obtain estimates of lifetime
characteristics under normal use conditions. The precision of the estimates is
influenced by several factors, such as the number of units tested, the
duration of the testing period, the frequency of measurements, and, on
particular, the choice of the stress levels to which the units are exposed.

A vast amount of literature is devoted to the analysis of Accelerated
Degradation Tests, see, for example, (Meeker and Escobar, 2014) for a
comprehensive survey on various approaches in the literature used to assess
reliable information from degradation data. In more detail, (Bagdonavicius and
Nikulin, 2001) present models and methods of statistical analysis for
Accelerated Degradation Tests and further references can be found there. As
additional sources, (Nelson, 2005) provides an extensive list of references
related to accelerated test planning and (Limon et al., 2017) review prominent
methods for statistical inference and optimal design of accelerated testing
plans. There different types of test planning strategies are categorized
according to their merits and drawbacks and research trends are provided. (Li
and Kececioglu, 2006) presents an analytical method for the optimum planning
of Accelerated Degradation Tests with an application to the reliability of
Light-Emitting Diodes. There the author states that the variability of the
measured units have a substantial impact on the accuracy of estimation.
Therefore these random effects should be encountered in the choice of the
experimental settings for the Accelerated Degradation Tests. Based on the
observation that ignoring the variability in the normal use conditions may
lead to significant prediction errors, (Liao and Elsayed, 2006) extend
Accelerated Degradation Test models to predict field reliability by
considering variations in the stress levels by considering a degradation
process represented by a Brownian motion with linear drift via a stochastic
differential equation. (Wang et al., 2017) propose a $M$-optimality criterion
for designing constant stress Accelerated Degradation Tests when the
degradation path can be represented by an inverse Gaussian process with
covariates and random effects. This criterion focuses on a degradation
mechanism equivalence rather than on the evaluation precision or the
prediction accuracy which are usually employed in traditional optimization
criteria. Those authors prove that, with a slightly relaxed requirement of
prediction accuracy, the obtained optimum designs minimize the dispersion of
the estimated acceleration factor between the normal stress level and a higher
accelerated stress level. Wiener processes (Brownian motions) are intensively
used to represent degradation paths in Accelerated Degradation Testing, see
(Ye et al., 2015) and (Guan et al., 2016). For instance, (Lim and Yum, 2011)
develop optimal Accelerated Degradation Test plans assuming that the constant
stress loading method is employed and the degradation characteristics follows
a Wiener process. These authors determine the test stress levels and the
proportion of test units allocated to each stress level such that the
asymptotic variance of the maximum likelihood estimator of a particular
quantile of the lifetime distribution at the normal use condition is
minimized. In addition, compromise plans are also developed for checking the
validity of the relationship between the model parameters and the stress
variable. In a case study for random effects in degradation of semiconductors,
(Lu et al., 1997) propose a repeated measurements model with random regression
coefficients and a standard deviation function for analyzing linear
degradation data. The authors utilize several large sample interval estimation
procedures to estimate the failure time distribution and its quantiles.

On the other hand, the general theory of optimal design of experiments is well
developed in the mathematical context of approximate designs which allow for
analytical solutions (see e. g. (Silvey, 1980) or (Atkinson et al., 2007)). In
addition, (Schwabe, 1996b) deals with the theory of optimal designs for multi-
factor models which can be used here to treat more than one stress variable
and the choice of time plans simultaneously under various interaction
structures. In the presence of random effects, (Entholzner et al., 2005)
derive that for single samples the optimal designs for fixed effects models
retain their optimality for linear optimality criteria. (Debusho and Haines,
2008) show that this also holds for $D$-optimality in linear models when only
the intercept is random. However, in a multi-sample situation (Schmelter,
2007) and (Schwabe and Schmelter, 2008) exhibit that the variability of the
intercept has a non-negligible influence on the $D$-optimal design. In the
case of random slope effects this dependence already occurs in single samples
as outlines by (Schmelter et al., 2007). (Dette et al., 2010) consider the
problem of constructing $D$-optimal designs for linear and nonlinear random
effect models with applications in population pharmacokinetics. These authors
present a new approach to determine efficient designs for nonlinear least
squares estimation which addresses the problem of additional correlation
between observations within units. Based on geometrical arguments, (Graßhoff
et al., 2012) derive $D$-optimal designs for random coefficient regression
models when only one observation is available per unit, a situation which
occurs in destructive testing. (Mentre et al., 1997) present an approach to
optimal design of experiments for random effects regression models in the
presence of cost functions related to costs per unit and costs per measurement
with applications to toxicokinetics.

The present approach is based on the discussion paper by (Weaver and Meeker,
2014) in which two case studies are introduced for optimal planning of
repeated measures Accelerated Degradation Tests. There the authors consider
the influence of a single stress variable and use a criterion based on a
large-sample approximation of the precision for estimating a quantile of the
failure-time distribution under normal use conditions. We will adopt this
approach, generalize the results presented there to more general models, and
extend the design optimization also to generate optimal time plans.

The present paper is organized as follows. Section 2 starts with a motivation
example based on a case study in (Weaver and Meeker, 2014). In Sections 3, 4
and 5 we state the general model formulation, specify the maximum-likelihood
estimation and exhibit the corresponding information matrix. Basic concepts of
optimal design theory in the present context are collected in Section 6 while
Section 7 is devoted to the idea of soft failure due to degradation, where we
derive the design optimality criterion for estimating a quantile of the
failure time distribution under normal use conditions. In Section 8 optimal
designs are characterized when the time plan for repeated measurements at the
testing units is fixed in advance. In Section 9 also the measurement times are
optimized under the constraint that measurements are taken according to the
same time plan for all units, and in Section 10 optimal measurement times are
determined in the setting of destructive testing. The paper closes with a
short discussion in Section 11. Throughout the paper the theoretical concepts
are illustrated by two accompanying running examples. Some technical result is
deferred to an Appendix.

##  2 Introductory example

Before formulating our general degradation model in section 3, we start in
this section for motivation with the description of a simple introductory
example based on (Weaver and Meeker, 2014).

######  Example 1.

The model proposed in (Weaver and Meeker, 2014) is a linear mixed effect model
with a single stress variable $x$. In this model there are $n$ testing units
for which degradation $y_{ij}$ is observed at $k$ time points $t_{j}$,
$j=1,...,k$. The (standardized) stress variable $x$ can be chosen by the
experimenter from the design region $\mathcal{X}=[0,1]$.

On the unit level the response $y_{ij}$ for the degradation of testing unit
$i$ at time $t_{j}$ is represented by

$y_{ij}=\beta_{i,1}+\beta_{2}x_{i}+\beta_{i,3}t_{j}+\beta_{4}x_{i}t_{j}+\varepsilon_{ij},$
(2.1)

where the intercept $\beta_{i,1}$ is the mean degradation of unit $i$ at time
$t=0$ under the stress level $x=0$, $\beta_{2}$ is the common (not unit
specific) mean increase in degradation depending on the stress variable $x$,
$\beta_{i,3}$ is the mean increase in degradation of unit $i$ over time $t$
when the stress level is set to $x=0$, and $\beta_{4}$ is the interaction
effect between time and stress. The measurement errors $\varepsilon_{ij}$ are
assumed to be realizations of a normally distributed error variable with mean
zero and error variance $\sigma^{2}_{\varepsilon}$.

On the whole experiment level the unit specific parameters
$(\beta_{i,1},\beta_{i,3})^{T}$ of the units are assumed to be realizations of
a bivariate normal distribution with mean $(\beta_{1},\beta_{3})^{T}$ and a
variance covariance matrix
$\boldsymbol{\Sigma}=\left(\begin{array}[]{cc}\sigma_{1}^{2}&\rho\sigma_{1}\sigma_{2}\\\
\rho\sigma_{1}\sigma_{2}&\sigma_{2}^{2}\end{array}\right)$. All random effect
parameters and measurement errors are assumed to be independent both within as
well as between units. Under well controlled measuring testing conditions, the
variability of the response is completely described by both the unit to unit
variability $\boldsymbol{\Sigma}$ and the within unit variability of the
measurement errors.

To illustrate the situation some virtual degradation paths $y_{ij}$,
$j=1,...,k$, are depicted in Figure 1 (left panel) for three different values
of the stress variable $x$. There are three units shown at each value of the
stress level ($n=9$) and $k=11$ equally spaced measurement times $t_{j}$. The
roughness of the paths is due to the measurement errors $\varepsilon_{ij}$ The
corresponding underlying mean degradation paths
$\mu_{i}(x_{i},t)=\beta_{i,1}+\beta_{2}x_{i}+\beta_{i,3}t+\beta_{4}x_{i}t$,
corrected for the measurement errors, are shown in the right panel of Figure
1. These mean degradation paths are represented by straight lines over time,
where both the intercept and the slope may vary across units around an
aggregate value determined by the value $x_{i}$ of the stress variable.

Figure 1: Observed degradation paths (left panel) and corresponding mean
degradation paths (right panel)

For the analysis of degradation under normal use it is assumed that the model
equation 2.1 is also valid at the normal use condition $x_{u}$, where
typically $x_{u}<0$ on the standardized scale, i. e.
$\mu_{u}(x_{u},t)=\beta_{u,1}+\beta_{2}x_{u}+\beta_{u,3}t+\beta_{4}x_{u}t$
describes the mean degradation of a future unit $u$ at normal use condition
$x_{u}$ and time $t$, and the unit specific parameters
$(\beta_{u,1},\beta_{u,3})^{T}$ are bivariate normal with mean
$(\beta_{1},\beta_{3})^{T}$ and variance covariance matrix
$\boldsymbol{\Sigma}$.

The degradation is assumed to undergo soft failure, i. e. the unit $u$ will be
considered to fail due to degradation when its mean degradation exceeds a
predetermined threshold $y_{0}$. The corresponding time $t_{u}$, for which
$\mu_{u}(x_{u},t_{u})=y_{0}$, will be called the failure time of unit $u$
under normal use condition due to degradation. These failure times vary across
different unit because of the unit specific parameters $\beta_{u,1}$ and
$\beta_{u,3}$.

In both panels of Figure 1 the predetermined failure threshold $y_{0}=50$ is
indicated by a horizontal line. As typical for degradation studies failure
does not occur during the time of experiment even for the highest stress
level.

(Weaver and Meeker, 2014) derived locally $c$-optimal designs for certain
characteristics of the mean degradation curves (for future units) under normal
use condition $x_{u}$. The optimal design is always supported on the extremal
points of the design region $\mathcal{X}$ and the optimal weights depend on
the nominal values of the model parameters. In Table 1 we reproduce the
nominal values of Example 7.2 by (Weaver and Meeker, 2014) on scar width
growth after standardization for further use.

Table 1: Nominal values for Example 1 | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\sigma_{1}$ | $\sigma_{2}$ | $\rho$ | $\sigma_{\varepsilon}$ | ${x}_{u}$ | $y_{0}$
---|---|---|---|---|---|---|---|---|---|---
$f_{j}(x,t)$ | $1$ | $x$ | $t$ | $xt$ |  |  |  |  |  |
| $2.397$ | $1.629$ | $1.018$ | $0.0696$ | $0.114$ | $0.105$ | $-0.143$ | $0.048$ | $-0.056$ | $3.912$

In the following section we extend the work of (Weaver and Meeker, 2014) by
considering a general mixed effects model with multiple stress variables
$\mathbf{x}$ given that the general model can be applied to linear mixed
effect models, linear additive models or exponential models. Further, this
generalization allows for degradation models with and without interactions
between the multiple design variables $\mathbf{x}$.

##  3 Formulation of the model

In this section, we give a general formulation of a mixed effects regression
model incorporating a product-type structure with complete interactions
between the stress and the time variable. This general formulation can cover
the situation described in the introductory Example1 of the previous section,
but may readily be extended to more than one stress variable and to more
complex marginal models for both the stress variables and the time variable,
as indicated later. To become more specific we assume that there are $n$
testing units $i=1,...,n$, for which degradation $y_{ij}$ is to be measured at
$k$ subsequent time points $t_{j}$, $j=1,...,k$, $t_{1}<...<t_{k}$. Each unit
$i$ is observed under a value $\mathbf{x}_{i}$ of the stress variable(s),
which is kept fixed for each unit throughout the degradation process, but may
differ from unit to unit. The number $k$ of measurements and the time points
are the same for all units. The measurements $y_{ij}$ are regarded as
realizations of random variables $Y_{ij}$ which are described by a
hierarchical model.

For each unit $i$ the observation $Y_{ij}$ at time point $t_{j}$ is given by

$Y_{ij}=\mu_{i}(\mathbf{x}_{i},t)+\varepsilon_{ij},$ (3.1)

where $\mu_{i}(\mathbf{x},t)$ is the mean degradation of unit $i$ at time $t$,
when stress $\mathbf{x}$ is applied to unit $i$, and $\varepsilon_{ij}$ is the
associated measurement error at time point $t_{j}$. The mean degradation
$\mu_{i}(\mathbf{x},t)$ is assumed to be given by a linear model equation in
the stress variable $\mathbf{x}$ and time $t$,

$\mu_{i}(\mathbf{x},t)=\sum_{r=1}^{p}\beta_{i,r}f_{r}(\mathbf{x}_{i},t_{j})=\mathbf{f}(\mathbf{x}_{i},t_{j})^{T}\boldsymbol{\beta}_{i}$
(3.2)

where
$\mathbf{f}(\mathbf{x},t)=(f_{1}(\mathbf{x},t),...,f_{p}(\mathbf{x},t))^{T}$
is a $p$-dimensional vector of known regression functions
$f_{q}(\mathbf{x},t)$ in both the stress variable(s) $\mathbf{x}$ and the time
$t$, $\boldsymbol{\beta}_{i}=(\beta_{i,1},...,\beta_{i,p})^{T}$ is a
$p$-dimensional vector of unit specific parameters $\beta_{i,q}$. Hence, the
response is given by

$Y_{ij}=\mathbf{f}(\mathbf{x}_{i},t_{j})^{T}\boldsymbol{\beta}_{i}+\varepsilon_{ij}.$
(3.3)

The measurement error $\varepsilon_{ij}$ is assumed to be normally distributed
with zero mean and some potentially time dependent error variance
$\sigma_{\epsilon,j}^{2}$
($\varepsilon_{ij}\sim\mathrm{N}(0,\sigma_{\epsilon,j}^{2})$). Moreover, the
error terms may be correlated within a unit over time. So, in general the
vector
$\boldsymbol{\varepsilon}_{i}=(\varepsilon_{i1},...,\varepsilon_{ik})^{T}$ of
errors associated with the $k$ observations within one unit $i$ is
$k$-dimensional multivariate normally distributed with zero mean and positive
definite variance covariance matrix $\boldsymbol{\Sigma}_{\varepsilon}$
($\boldsymbol{\varepsilon}_{i}\sim\mathrm{N}(\mathbf{0},\boldsymbol{\Sigma}_{\varepsilon})$).

######  Example (Example 1 cont.).

In the introductory model

$Y_{ij}=\beta_{i,1}+\beta_{2}x_{i}+\beta_{i,3}t_{j}+\beta_{4}x_{i}t_{j}+\varepsilon_{ij}$

of Section 2 there is only one stress variable $\mathbf{x}=x$, the vector of
regression functions is given by $\mathbf{f}(x,t)=(1,x,t,xt)^{T}$, and the
unit specific vector of parameters is
$\boldsymbol{\beta}_{i}=(\beta_{i,1},\beta_{2},\beta_{i,3},\beta_{4})^{T}$
($p=4$). The error terms are assumed to be homoscedastic and independent, i.
e. $\boldsymbol{\Sigma}_{\varepsilon}=\sigma_{\varepsilon}^{2}\mathbf{I}_{k}$,
where $\mathbf{I}_{k}$ denotes the $k\times k$ identity matrix.

For the regression functions $\mathbf{f}(\mathbf{x},t)$ we suppose a product-
type structure with complete interactions between the stress variable
$\mathbf{x}$ and the time $t$, i. e. there are marginal regression functions
$\mathbf{f}_{1}(\mathbf{x})=(f_{11}(\mathbf{x}),...,f_{1p_{1}}(\mathbf{x}))^{T}$
and $\mathbf{f}_{2}(t)=(f_{21}(t),...,f_{2p_{2}}(t))^{T}$ of dimension $p_{1}$
and $p_{2}$ which only depend on the stress variable $\mathbf{x}$ and the time
$t$, respectively, and the vector
$\mathbf{f}(\mathbf{x},t)=\mathbf{f}_{1}(\mathbf{x})\otimes\mathbf{f}_{2}(t)$
of regression functions factorizes into its marginal counterparts
($p=p_{1}p_{2}$). Here “$\otimes$” denotes the Kronecker product of matrices
or vectors. Then the observation $Y_{ij}$ can be written as

$Y_{ij}=\sum_{r=1}^{p_{1}}\sum_{s=1}^{p_{2}}\beta_{i,rs}f_{1r}(\mathbf{x}_{i})f_{2s}(t_{j})+\varepsilon_{ij}=(\mathbf{f}_{1}(\mathbf{x})\otimes\mathbf{f}_{2}(t))^{T}\boldsymbol{\beta}_{i}+\varepsilon_{ij},$
(3.4)

where for notational convenience the entries of the vector
$\boldsymbol{\beta}_{i}=(\beta_{i,11},...,\beta_{i,1p_{2}},...,\beta_{i,p_{1}p_{2}})$
of parameters are relabeled lexicographically according to their associated
marginal regression functions ($q=(r-1)p_{2}+s$, $r=1,...,p_{1}$,
$s=1,...,p_{2}$).

######  Example (Example 1 cont.).

In the introductory model of Section 2 the rearranged vector of regression
functions $\mathbf{f}(x,t)=(1,t,x,xt)^{T}=((1,x)\otimes(1,t))^{T}$ factorizes
into the vector $\mathbf{f}_{1}(x)=(1,x)^{T}$ of a simple linear regression in
the stress variable $x$ ($p_{1}=2$) and the vector
$\mathbf{f}_{2}(t)=(1,t)^{T}$ of a simple linear regression in time $t$
($p_{2}=2$). The unit specific vector of parameters $\boldsymbol{\beta}_{i}$
is relabeled by $\beta_{i,11}=\beta_{i,1}$, $\beta_{i,12}=\beta_{i,3}$,
$\beta_{21}=\beta_{2}$, and $\beta_{22}=\beta_{4}$, where for $\beta_{21}$ and
$\beta_{22}$ the index $i$ of the unit is suppressed as these parameters do
not depend on the unit. Thus the model equation can be rewritten as

$Y_{ij}=\beta_{i,11}+\beta_{i,12}t_{j}+\beta_{21}x_{i}+\beta_{22}x_{i}t_{j}+\varepsilon_{ij}=((1,x_{i})\otimes(1,t_{j}))^{T}\boldsymbol{\beta}_{i}+\varepsilon_{ij}.$

Moreover, as in the introductory example, we will assume throughout that the
marginal regression function
$\mathbf{f}_{1}(\mathbf{x})=(f_{11}(\mathbf{x}),...,f_{1p_{1}}(\mathbf{x}))^{T}$
of the stress variable $\mathbf{x}$ contains a constant term,
$f_{11}(\mathbf{x})\equiv 1$ say, which is a common assumption in the majority
of situations, and that only the leading $p_{2}$ parameters
$\beta_{i,11},...,\beta_{i,1p_{2}}$ of $\boldsymbol{\beta}_{i}$ associated
with this constant term are unit specific. All other parameters in
$\boldsymbol{\beta}_{i}$ are assumed to take the same value $\beta_{rs}$,
$r=2,...,p_{1}$, $s=1,...,p_{2}$, for all individuals $i=1,...,n$. Hence, for
unit $i$ the model (3.4) can be rewritten as

$Y_{ij}=\left(\mathbf{f}_{1}(\mathbf{x}_{i})\otimes\mathbf{f}_{2}(t_{j})\right)^{T}\boldsymbol{\beta}+\mathbf{f}_{2}(t_{j})^{T}\boldsymbol{\gamma}_{i}+\varepsilon_{ij},$
(3.5)

where $\boldsymbol{\beta}=(\beta_{11},...,\beta_{p_{1}p_{2}})^{T}$ is the
vector of fixed effect (aggregate) parameters (averaged over the units)
associated with the constant term in the regression functions of the stress
variable $\mathbf{x}$ and
$\boldsymbol{\gamma}_{i}=(\gamma_{i1},...,\gamma_{ip_{2}})^{T}$ is the
$p_{2}$-dimensional vector of unit specific deviations
$\gamma_{is}=\beta_{i,1s}-\beta_{1s}$, $s=1,...,p_{2}$, from the corresponding
aggregate parameters.

On the aggregate level it is assumed that the units are representatives of a
larger entity. The deviations of the units from the aggregate value are then
modeled as random effects, i. e. they are $p_{2}$-dimensional multivariate
normal with zero mean and variance-covariance matrix
$\boldsymbol{\Sigma}_{\gamma}$
($\boldsymbol{\gamma}_{i}\sim\mathrm{N}(\mathbf{0},\boldsymbol{\Sigma}_{\gamma})$).
All vectors $\boldsymbol{\gamma}_{i}$ of random effects and all vectors
$\boldsymbol{\varepsilon}_{i}$ of random errors are assumed to be independent.

######  Example (Example 1 cont.).

In the introductory model of Section 2 the two-dimensional vector of
deviations associated with the regression functions $f_{1}(x,t)=1$ and
$f_{2}(x,t)=t$ is
$\boldsymbol{\gamma}_{i}=(\beta_{i,11}-\beta_{11},\beta_{i,12}-\beta_{12})^{T}$
with general $2\times 2$ variance-covariance matrix
$\boldsymbol{\Sigma}_{\gamma}=\boldsymbol{\Sigma}$. This leads to the mixed
effects model equation in standard notation

$Y_{ij}=\beta_{11}+\beta_{12}t_{j}+\beta_{21}x_{i}+\beta_{22}x_{i}t_{j}+\gamma_{i1}+\gamma_{i2}t_{j}+\varepsilon_{ij}=((1,x_{i})\otimes(1,t_{j}))^{T}\boldsymbol{\beta}+(1,t_{j})^{T}\boldsymbol{\gamma}_{i}+\varepsilon_{ij}.$

In vector notation the $k$-dimensional vector
$\mathbf{Y}_{i}=(Y_{i1},...,Y_{ik})^{T}$ of observations for unit $i$ can be
expressed as

$\mathbf{Y}_{i}=(\mathbf{f}_{1}(\mathbf{x}_{i})^{T}\otimes\mathbf{F}_{2})\boldsymbol{\beta}+\mathbf{F}_{2}\,\boldsymbol{\gamma}_{i}+\boldsymbol{\varepsilon}_{i},$

where $\mathbf{F}_{2}=(\mathbf{f}_{2}(t_{1}),...,\mathbf{f}_{2}(t_{k}))^{T}$
is the $k\times p_{2}$ marginal design matrix for the time variable. Then
$\mathbf{Y}_{i}$ is $k$-dimensional multivariate normally distributed with
mean
$(\mathbf{f}_{1}(\mathbf{x}_{i})^{T}\otimes\mathbf{F}_{2})\boldsymbol{\beta}$
and variance covariance matrix
$\mathbf{V}=\mathbf{F}_{2}{\boldsymbol{\Sigma}}_{\gamma}\mathbf{F}_{2}^{T}+{\boldsymbol{\Sigma}}_{\varepsilon}$.
The variance covariance matrix $\mathbf{V}$ is not affected by the choice of
the stress level $\mathbf{x}_{i}$ and, hence, equal for all units $i$.

######  Example (Example 1 cont.).

In the introductory model of Section 2 the marginal design matrix

$\mathbf{F}_{2}=\left(\begin{array}[]{ccc}1&...&1\\\
t_{1}&...&t_{k}\end{array}\right)^{T}$

is the common design matrix for simple linear regression over time. The
variance covariance matrix
$\mathbf{V}=\mathbf{F}_{2}{\boldsymbol{\Sigma}}\mathbf{F}_{2}^{T}+\sigma^{2}_{\varepsilon}\mathbf{I}_{k}$
has diagonal entries
$v_{jj}=\mathrm{Var}(Y_{ij})=\sigma_{1}^{2}+2\rho\sigma_{1}\sigma_{2}t_{j}+\sigma_{2}^{2}t_{j}^{2}+\sigma_{\varepsilon}^{2}$
and off-diagonal entries
$v_{jj^{\prime}}=\mathrm{cov}(Y_{ij},Y_{ij^{\prime}})=\sigma_{1}^{2}+\rho\sigma_{1}\sigma_{2}(t_{j}+t_{j^{\prime}})+\sigma_{2}^{2}t_{j}t_{j^{\prime}}$.

In total, for the observations of all $n$ units the stacked $nk$-dimensional
response vector $\mathbf{Y}=(\mathbf{Y}_{1}^{T},...,\mathbf{Y}_{n}^{T})^{T}$
can be represented in matrix notation as

$\mathbf{Y}=({\mathbf{F}_{1}\otimes\mathbf{F}_{2}})\boldsymbol{\beta}+(\mathbf{I}_{n}\otimes{\mathbf{F}_{2}})\,{\boldsymbol{\gamma}}+\boldsymbol{\varepsilon},$
(3.6)

where
$\mathbf{F}_{1}=(\mathbf{f}_{1}(\mathbf{x}_{1}),...,\mathbf{f}_{1}(\mathbf{x}_{n}))^{T}$
is the $n\times p_{1}$ marginal design matrix for the stress variables across
units,
${\boldsymbol{\gamma}}=({\boldsymbol{\gamma}}_{1}^{T},...,{\boldsymbol{\gamma}}_{n}^{T})^{T}$
is the $np_{2}$-dimensional stacked parameter vector of random effects and
${\boldsymbol{\varepsilon}}=({\boldsymbol{\varepsilon}}_{1}^{T},...,{\boldsymbol{\varepsilon}}_{n}^{T})^{T}$
is the $nk$-dimensional stacked vector of random errors. Such a model equation
is sometimes called the “marginal model” for the response $\mathbf{Y}$, but
should not be confused with models marginalized for the covariates
$\mathbf{x}$ and $t$, respectively (see the decomposition at the end of
Section 5). Note that the vectors
${\boldsymbol{\gamma}}\sim\mathrm{N}(\mathbf{0},\mathbf{I}_{n}\otimes{\boldsymbol{\Sigma}}_{\gamma})$
of all random effects and the vector
$\boldsymbol{\varepsilon}\sim\mathrm{N}(\mathbf{0},\mathbf{I}_{n}\otimes{\boldsymbol{\Sigma}}_{\varepsilon})$
are multivariate normal. Hence, the vector $\mathbf{Y}$ of all observations is
$nk$-dimensional multivariate normal,
$\mathbf{Y}\sim\mathrm{N}(\mathbf{0},\mathbf{I}_{n}\otimes\mathbf{V})$.

######  Example (Example 1 cont.).

In the introductory model of Section 2 the marginal design matrix

$\mathbf{F}_{1}=\left(\begin{array}[]{ccc}1&...&1\\\
x_{1}&...&x_{n}\end{array}\right)^{T}$

is the common design matrix for simple linear regression on the stress
variable $x$.

For the analysis of degradation under normal use we further assume that the
general model 3.5 is also valid at the normal use condition $\mathbf{x}_{u}$,
where typically $\mathbf{x}_{u}\not\in\mathcal{X}$, i. e.

$\mu(\mathbf{x}_{u},t)=(\mathbf{f}_{1}(\mathbf{x}_{u})\otimes\mathbf{f}_{2}(t))^{T}\boldsymbol{\beta}+\mathbf{f}_{2}(t)^{T}\boldsymbol{\gamma}_{u}$
(3.7)

describes the mean degradation of a future unit $u$ at normal use condition
$\mathbf{x}_{u}$ and time $t$, and the random effects
$\boldsymbol{\gamma}_{u}$ are $k$-dimensional multivariate normal with mean
zero and variance covariance matrix $\boldsymbol{\Sigma}_{\gamma}$.

In the following we will consider thoroughly a more complex example, where two
stress variables are involved.

######  Example 2.

In this example the degradation is influenced by two standardized accelerating
stress variables $x_{1}$ and $x_{2}$ which act linearly on the response with a
potential interaction effect associated with $x_{1}x_{2}$. The two stress
variables $x_{1}$ and $x_{2}$ can be chosen independently from marginal design
regions $\mathcal{X}_{1}=\mathcal{X}_{2}=[0,1]$, respectively. Also the time
is assumed to act linearly on the degradation and all interactions between
stress variables and time are present as in Example 1.

If, for testing unit $i$, the stress variables are set to $x_{i1}$ and
$x_{i2}$ the response $y_{ij}$ at time $t_{j}$ is given by

$y_{ij}=\beta_{i,1}+\beta_{2}x_{i1}+\beta_{3}x_{i2}+\beta_{4}x_{i1}x_{i2}+\beta_{i,5}t_{j}+\beta_{6}x_{i1}t_{j}+\beta_{7}x_{i2}t_{j}+\beta_{8}x_{i1}x_{i2}t_{j}+\varepsilon_{ij},$
(3.8)

where the intercept $\beta_{i,1}$ is the mean degradation of unit $i$ at time
$t=0$ under the stress levels $x_{1}=0$ and $x_{2}=0$, $\beta_{2}$ is the
common (not unit specific) mean increase in degradation depending on the
stress variable $x_{1}$ when $x_{2}=0$, $\beta_{3}$ is the common mean
increase in degradation depending on the stress variable $x_{2}$ when
$x_{1}=0$, and $\beta_{4}$ is the interaction effect between the two stress
variables. Accordingly $\beta_{i,5}$ is the mean increase in degradation of
unit $i$ over time $t$ when the stress levels are set to $x_{1}=0$ and
$x_{2}=0$, $\beta_{6}$ is the interaction effect between time and the stress
variable $x_{1}$ when $x_{2}=0$, $\beta_{7}$ is the interaction effect between
time and the stress variable $x_{2}$ when $x_{1}=0$, and $\beta_{8}$ is the
second-order interaction effect between time and the two stress variables.
Also here only the parameters $\beta_{i,1}$ and $\beta_{i,5}$ associated with
the constant term in the stress variables may vary across units. On the
aggregate level these two unit parameters are again assumed to be normally
distributed with means $\mathrm{E}(\beta_{i,1})=\beta_{1}$ and
$\mathrm{E}(\beta_{i,5})=\beta_{5}$ and variance covariance matrix
$\boldsymbol{\Sigma}=\left(\begin{array}[]{cc}\sigma^{2}_{1}&\rho\sigma_{1}\sigma_{2}\\\
\rho\sigma_{1}\sigma_{2}&\sigma_{2}^{2}\end{array}\right)$.

After rearranging terms and relabeling the parameters the model can be
rewritten as

$Y_{ij}=(\mathbf{f}_{11}(x_{i1})\otimes\mathbf{f}_{12}(x_{i2})\otimes\mathbf{f}_{2}(t_{j}))^{T}\boldsymbol{\beta}+\mathbf{f}_{2}(t_{j})^{T}\boldsymbol{\gamma}_{i}+\varepsilon_{ij},$
(3.9)

where $\mathbf{f}_{11}(x_{1})=(1,x_{1})^{T}$,
$\mathbf{f}_{12}(x_{2})=(1,x_{2})^{T}$ and $\mathbf{f}_{2}(t)=(1,t)^{T}$ are
the marginal regression functions for the stress variables $x_{1}$, $x_{2}$
and the time variable $t$, respectively,
$\boldsymbol{\beta}=(\beta_{111},\beta_{112},\beta_{121},\beta_{122},\beta_{211},\beta_{212},\\\
\beta_{221},\beta_{222})^{T}=(\beta_{1},\beta_{5},\beta_{3},\beta_{7},\beta_{2},\beta_{6},\beta_{4},\beta_{8})^{T}$
is the rearranged vector of aggregate parameters, and
$\boldsymbol{\gamma}_{i}=(\gamma_{i1},\gamma_{i2})^{T}=(\beta_{i,1}-\beta_{1},\beta_{i,5}-\beta_{5})^{T}$
is the vector of parameters for the deviations of unit $i$ from the aggregate
values. These deviations constitute again random effects with zero mean and
variance covariance matrix $\boldsymbol{\Sigma}$. With
$\mathbf{x}=(x_{1},x_{2})$,
$\mathcal{X}=\mathcal{X}_{1}\times\mathcal{X}_{2}=[0,1]^{2}$,
$\mathbf{f}_{1}(\mathbf{x})=(\mathbf{f}_{11}(x_{1})^{T},\mathbf{f}_{12}(x_{2})^{T})^{T}$,
$p_{1}=4$, $p_{2}=2$, and $p=8$ model 3.9 fits into the framework of the
general product-type model 3.4.

Later we will also consider shortly a model in two stress variables without
interaction terms between the stress variables.

##  4 Estimation of the model parameters

Under the distributional assumptions of normality for both the random effects
and the measurement errors the model parameters may be estimated by means of
the maximum likelihood method. Denote by
$\boldsymbol{\theta}=(\boldsymbol{\beta}^{T},\boldsymbol{\varsigma}^{T})^{T}$
the vector of all model parameters, where $\boldsymbol{\varsigma}$ collects
all variance covariance parameters from $\boldsymbol{\Sigma}_{\gamma}$ and
$\boldsymbol{\Sigma}_{\varepsilon}$ For the general model (3.6) the log-
likelihood is given by

$\ell(\boldsymbol{\theta};\mathbf{y})=-{\textstyle{\frac{nk}{2}}}\log(2\pi)-{\textstyle{\frac{n}{2}}}\log(\det(\mathbf{V}))-{\textstyle{\frac{1}{2}}}(\mathbf{y}-(\mathbf{F}_{1}\otimes\mathbf{F}_{2})\boldsymbol{\beta})^{T}(\mathbf{I}_{n}\otimes\mathbf{V})^{-1}(\mathbf{y}-(\mathbf{F}_{1}\otimes\mathbf{F}_{2})\boldsymbol{\beta}),$
(4.1)

where the variance covariance matrix
$\mathbf{V}=\mathbf{V}(\boldsymbol{\varsigma})$ of measurements per unit
depends only on $\boldsymbol{\varsigma}$. The maximum likelihood estimator of
$\boldsymbol{\beta}$ can be calculated as

$\displaystyle\widehat{\boldsymbol{\beta}}$ $\displaystyle=$
$\displaystyle((\mathbf{F}_{1}\otimes\mathbf{F}_{2})^{T}(\mathbf{I}_{n}\otimes\widehat{\mathbf{V}})^{-1}(\mathbf{F}_{1}\otimes\mathbf{F}_{2}))^{-1}(\mathbf{F}_{1}\otimes\mathbf{F}_{2})^{T}(\mathbf{I}_{n}\otimes\widehat{\mathbf{V}})^{-1}\mathbf{Y}$
(4.2) $\displaystyle=$
$\displaystyle((\mathbf{F}_{1}^{T}\mathbf{F}_{1})^{-1}\mathbf{F}_{1}^{T})\otimes((\mathbf{F}_{2}^{T}\widehat{\mathbf{V}}^{-1}\mathbf{F}_{2})^{T}\mathbf{F}_{2}^{T}\widehat{\mathbf{V}}^{-1})\mathbf{Y},$

if both $\mathbf{F}_{1}$ and $\mathbf{F}_{2}$ are of full column rank $p_{1}$
and $p_{2}$, respectively, and
$\widehat{\mathbf{V}}=\mathbf{V}(\widehat{\boldsymbol{\varsigma}})$, where
$\widehat{\boldsymbol{\varsigma}}$ is the maximum likelihood estimator of
$\boldsymbol{\varsigma}$. When $\mathbf{V}$ is known, at least up to a
multiplicative constant $\mathbf{V}=\sigma^{2}\mathbf{V}_{0}$, then
$\widehat{\boldsymbol{\beta}}$ is the best liner unbiased (general least
squares) estimator
$\widehat{\boldsymbol{\beta}}_{\mathrm{GLS}}=((\mathbf{F}_{1}^{T}\mathbf{F}_{1})^{-1}\mathbf{F}_{1}^{T})\otimes((\mathbf{F}_{2}^{T}\mathbf{V}_{0}^{-1}\mathbf{F}_{2})^{T}\mathbf{F}_{2}^{T}\mathbf{V}_{0}^{-1})\mathbf{Y}$
of $\boldsymbol{\beta}$. In particular, when the measurement errors are
uncorrelated and homoscedastic, i. e.
$\boldsymbol{\Sigma}_{\varepsilon}=\sigma_{\varepsilon}^{2}\mathbf{I}_{k}$,
then this estimator reduces to the ordinary least squares estimator
$\widehat{\boldsymbol{\beta}}_{\mathrm{OLS}}=((\mathbf{F}_{1}^{T}\mathbf{F}_{1})^{-1}\mathbf{F}_{1}^{T})\otimes((\mathbf{F}_{2}^{T}\mathbf{F}_{2})^{T}\mathbf{F}_{2}^{T})\mathbf{Y}$
because
$\mathbf{V}\mathbf{F}_{2}=\mathbf{F}_{2}(\boldsymbol{\Sigma}_{\gamma}\mathbf{F}_{2}^{T}\mathbf{F}_{2}+\sigma_{\varepsilon}^{2}\mathbf{I}_{p_{2}})$
by a result of (Zyskind, 1967). Hence, in the case of uncorrelated
homoscedastic measurement errors the maximum likelihood estimator of the
location parameters $\boldsymbol{\beta}$ does neither depend on the variance
covariance parameters nor on their estimates.

In general, the quality of the estimator $\widehat{\boldsymbol{\beta}}$ can be
measured in terms of its variance covariance matrix which is given by

$\mathrm{Cov}(\widehat{\boldsymbol{\beta}})=(\mathbf{F}_{1}^{T}\mathbf{F}_{1})^{-1}\otimes(\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2})^{-1}.$
(4.3)

By using the structure
$\mathbf{V}=\mathbf{F}_{2}\boldsymbol{\Sigma}_{\gamma}\mathbf{F}_{2}^{T}+\boldsymbol{\Sigma}_{\varepsilon}$
the last term can be calculated as

$(\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2})^{-1}=(\mathbf{F}_{2}^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F}_{2})^{-1}+\boldsymbol{\Sigma}_{\gamma}$
(4.4)

in terms of the variance covariance matrices $\boldsymbol{\Sigma}_{\gamma}$
and $\boldsymbol{\Sigma}_{\varepsilon}$ of the random effects and the
measurement errors, respectively.

##  5 Information

In general, the Fisher information matrix is defined as the variance
covariance matrix of the score function $\mathbf{U}$ which itself is defined
as the vector of first derivatives of the log likelihood with respect to the
components of the parameter vector $\boldsymbol{\theta}$. More precisely, let
$\mathbf{U}=(\frac{\partial}{\partial\theta_{1}}\ell(\boldsymbol{\theta};\mathbf{Y}),...,\frac{\partial}{\partial\theta_{q}}\ell(\boldsymbol{\theta};\mathbf{Y}))^{T}$,
where $q$ is the dimension of $\boldsymbol{\theta}$. Then for the full
parameter vector $\boldsymbol{\theta}$ the Fisher information matrix is
defined as $\mathbf{M}_{\boldsymbol{\theta}}=\mathrm{Cov}(\mathbf{U})$, where
the expectation is taken with respect to the distribution of $\mathbf{Y}$. The
Fisher information can also be computed as the expectations of the second
derivatives of the score function $U$, i. e.
$\mathbf{M}_{\boldsymbol{\theta}}=-\mathrm{E}\left(\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{T}}\ell(\boldsymbol{\theta};\mathbf{Y})\right)$.
Under common regularity conditions the maximum likelihood estimator
$\widehat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta}$ is consistent and
asymptotically normal with asymptotic variance covariance matrix equal to the
inverse $\mathbf{M}_{\boldsymbol{\theta}}^{-1}$ of the Fisher information
matrix $\mathbf{M}_{\boldsymbol{\theta}}$.

To specify the Fisher information matrix further, denote by
$\mathbf{{M}}_{\boldsymbol{\beta}}=-\mathrm{E}\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{T}}\ell(\boldsymbol{\theta};\mathbf{Y})\right)$,
$\mathbf{M}_{\boldsymbol{\varsigma}}=-\mathrm{E}\left(\frac{\partial^{2}}{\partial\boldsymbol{\varsigma}\partial\boldsymbol{\varsigma}^{T}}\ell(\boldsymbol{\theta};\mathbf{Y})\right)$,
$\mathbf{M}_{\boldsymbol{\beta}\boldsymbol{\varsigma}}=-\mathrm{E}\left(\frac{\partial^{2}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\varsigma}^{T}}\ell(\boldsymbol{\theta};\mathbf{Y})\right)$
and
$\mathbf{{M}}_{\boldsymbol{\varsigma}\boldsymbol{\beta}}=\mathbf{{M}}_{\boldsymbol{\beta}\boldsymbol{\varsigma}}^{T}$
the blocks of the Fisher information matrix corresponding to the second
derivatives with respect to $\boldsymbol{\beta}$ and $\boldsymbol{\varsigma}$
and the mixed derivatives, respectively. The mixed blocks can be seen to be
zero and the Fisher information matrix is block diagonal,

$\mathbf{M}_{\boldsymbol{\theta}}=\left(\begin{array}[]{cc}\mathbf{M}_{\boldsymbol{\beta}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{M}_{\boldsymbol{\varsigma}}\end{array}\right).$ (5.1)

Moreover, the block $\mathbf{M}_{\boldsymbol{\beta}}$ associated with the
aggregate location parameters $\boldsymbol{\beta}$ can be determined as

$\mathbf{M}_{\boldsymbol{\beta}}=(\mathbf{F}_{1}^{T}\mathbf{F}_{1})\otimes(\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2})$
(5.2)

which turns out to be the inverse of the variance covariance matrix for the
estimator $\widehat{\boldsymbol{\beta}}$ of $\boldsymbol{\beta}$, when
$\mathbf{V}$ is known. Actually, because the Fisher information matrix for
$\boldsymbol{\theta}$ is block diagonal, the inverse
$\mathbf{M}_{\boldsymbol{\beta}}^{-1}=(\mathbf{F}_{1}^{T}\mathbf{F}_{1})^{-1}\otimes(\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2})^{-1}$
of the block associated with $\boldsymbol{\beta}$ is the corresponding block
of the inverse of $\mathbf{M}_{\boldsymbol{\theta}}$ and is, hence, the
asymptotic variance covariance matrix of $\widehat{\boldsymbol{\beta}}$.

Accordingly the asymptotic variance covariance matrix for estimating the
variance parameters $\boldsymbol{\varsigma}$ is the inverse of the block
$\mathbf{M}_{\boldsymbol{\varsigma}}$. In the following we will call
$\mathbf{M}_{\boldsymbol{\beta}}$ and $\mathbf{M}_{\boldsymbol{\varsigma}}$
the information matrices for $\boldsymbol{\beta}$ and
$\boldsymbol{\varsigma}$, respectively, for short. The particular form of
$\mathbf{M}_{\boldsymbol{\varsigma}}$ will be not of interest here. However,
as the information matrix $\mathbf{M}_{\boldsymbol{\varsigma}}$ for the
variance parameters $\boldsymbol{\varsigma}$ is given by

$\textbf{M}_{\boldsymbol{\varsigma}}=\frac{n}{2}\left(\frac{\partial^{2}\log(\det(\mathbf{V}))}{\partial\boldsymbol{\varsigma}\partial\boldsymbol{\varsigma}^{T}}+\mathrm{tr}\left(\mathbf{V}\frac{\partial^{2}\mathbf{V}^{-1}}{\partial\boldsymbol{\varsigma}\partial\boldsymbol{\varsigma}^{T}}\right)\right).$

It is important to note that $\mathbf{M}_{\boldsymbol{\varsigma}}$ does not
depend on the settings $\mathbf{x}_{1},...,\mathbf{x}_{n}$ of the stress
variable in contrast to the information matrix
$\mathbf{M}_{\boldsymbol{\beta}}$ of the aggregate location parameters
$\boldsymbol{\beta}$. For the general product-type model (3.6) the information
matrix $\mathbf{M}_{\boldsymbol{\beta}}$ for the aggregate parameters
$\boldsymbol{\beta}$ factorizes according to

$\mathbf{M}_{\boldsymbol{\beta}}=\mathbf{M}_{1}\otimes\mathbf{M}_{2}$ (5.3)

into the information matrix $\mathbf{M}_{1}=\mathbf{F}_{1}^{T}\mathbf{F}_{1}$
in the marginal model

$Y_{i}^{(1)}=\mathbf{f}_{1}(\mathbf{x}_{i})^{T}\boldsymbol{\beta}^{(1)}+\varepsilon_{i}^{(1)},$
(5.4)

$i=1,...,n$, in the stress variable $\mathbf{x}$ with standardized
uncorrelated homoscedastic error terms, $\sigma^{2}_{\varepsilon^{(1)}}=1$,
and the information matrix
$\mathbf{M}_{2}=\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2}$ in the mixed
effects marginal model

$Y_{j}^{(2)}=\mathbf{f}_{2}(t_{j})^{T}\boldsymbol{\beta}^{(2)}+\mathbf{f}_{2}(t_{j})^{T}\boldsymbol{\gamma}^{(2)}+\varepsilon_{j}^{(2)},$
(5.5)

$j=1,...,k$, in the time variable $t$ with variance covariance matrices
$\boldsymbol{\Sigma}_{\gamma}$ and $\boldsymbol{\Sigma}_{\varepsilon}$ for the
random effects $\boldsymbol{\gamma}^{(2)}$ and measurement errors
$\boldsymbol{\varepsilon}^{(2)}=(\varepsilon_{1}^{(2)},...,\varepsilon_{k}^{(2)})^{T}$,
respectively. Then the information matrix $\mathbf{M}_{\boldsymbol{\theta}}$
in the full model depends on the settings $\mathbf{x}_{1},...,\mathbf{x}_{n}$
of the stress variable only through the information matrix $\mathbf{M}_{1}$ in
the first marginal model.

##  6 Design

The quality of the estimates will be measured in terms of the information
matrix and, hence, depends on both the settings of the stress variable and the
time points of measurements. When these variables are under the control of the
experimenter, then their choice will be called the design of the experiment.
Here we assume that the time plan $\mathbf{t}=(t_{1},...,t_{k})^{T}$ for the
time points of measurements within units is fixed in advance and is not under
disposition of the experimenter. Then only the settings
$\mathbf{x}_{1},...,\mathbf{x}_{n}$ of the stress variable $\mathbf{x}$ can be
adjusted to the units $i=1,...,n$. Their choice
$(\mathbf{x}_{1},...,\mathbf{x}_{n})$ is then called an “exact” design, and
their influence on the performance of the experiment is indicated by adding
them as an argument to the information matrices,
$\mathbf{M}_{\boldsymbol{\theta}}(\mathbf{x}_{1},...,\mathbf{x}_{n})$,
$\mathbf{M}_{\boldsymbol{\beta}}(\mathbf{x}_{1},...,\mathbf{x}_{n})$, and
$\mathbf{M}_{1}(\mathbf{x}_{1},...,\mathbf{x}_{n})$, where appropriate. Remind
that both $\mathbf{M}_{\boldsymbol{\varsigma}}$ and $\mathbf{M}_{2}$ do not
depend on the design for the stress variable.

As
$\mathbf{M}_{1}(\mathbf{x}_{1},...,\mathbf{x}_{n})=\sum_{i=1}^{n}\mathbf{f}(\mathbf{x}_{i})\mathbf{f}(\mathbf{x}_{i})^{T}$
it can easily be seen that the information matrices do not depend on the order
of the setting but only on their mutually distinct settings,
$\mathbf{x}_{1},...,\mathbf{x}_{m}$ say, and their corresponding frequencies
$n_{1},...,n_{m}$, such that $\sum_{i=1}^{m}n_{i}=n$, i. e.
$\mathbf{M}_{1}=\sum_{i=1}^{m}n_{i}\mathbf{f}(\mathbf{x}_{i})\mathbf{f}(\mathbf{x}_{i})^{T}$.
Finding optimal exact designs is, in general, a difficult task of discrete
optimization. To circumvent this problem we follow the approach of approximate
designs propagated by (Kiefer, 1959) in which the requirement of integer
numbers $n_{i}$ of testing units at a stress level $\mathbf{x}_{i}$ is
relaxed. Then continuous methods of convex optimization can be employed (see
e. g. (Silvey, 1980)) and efficient exact designs can be derived by rounding
the optimal numbers to nearest integers. This approach is, in particular, of
use when the number $n$ of units is sufficiently large. Moreover, the
frequencies $n_{i}$ will be replaced by proportions $w_{i}=n_{i}/n$, because
the total number $n$ of units does not play a role in the optimization. Thus
an approximate design $\xi$ is defined by a finite number of settings
$\mathbf{x}_{i}$, $i=1,...,m$, from an experimental region $\mathcal{X}$ with
corresponding weights $w_{i}\geq 0$ satisfying $\sum_{i=1}^{m}w_{i}=1$ and is
denoted by

$\xi=\left(\begin{array}[]{ccc}\mathbf{x}_{1}&...&\mathbf{x}_{m}\\\
w_{1}&...&w_{m}\end{array}\right),$ (6.1)

The corresponding standardized, per unit information matrices are accordingly
defined as

$\mathbf{M}_{1}(\xi)=\sum_{i=1}^{m}w_{i}\mathbf{f}_{1}(\mathbf{x}_{i})\mathbf{f}_{1}(\mathbf{x}_{i})^{T}$
(6.2)

for the marginal model on itself or by plugging (6.2) in into the
standardized, per unit information matrix

$\mathbf{M}_{\boldsymbol{\beta}}(\xi)=\mathbf{M}_{1}(\xi)\otimes\mathbf{M}_{2}$
(6.3)

for the aggregate parameters $\boldsymbol{\beta}$, where again
$\mathbf{M}_{2}=\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2}$, and

$\mathbf{M}_{\boldsymbol{\theta}}(\xi)=\left(\begin{array}[]{cc}\mathbf{M}_{1}(\xi)\otimes\mathbf{M}_{2}&\mathbf{0}\\\
\mathbf{0}&\widetilde{\mathbf{M}}_{\boldsymbol{\varsigma}}\end{array}\right)$
(6.4)

or the full parameter vector $\boldsymbol{\theta}$, where now
$\widetilde{\mathbf{M}}_{\boldsymbol{\varsigma}}=\frac{1}{n}\mathbf{M}_{\boldsymbol{\varsigma}}$
is the standardized, per unit information for the variance parameters
$\boldsymbol{\varsigma}$. If all $nw_{i}$ are integer, then these standardized
versions coincide with the information matrices of the corresponding exact
design up to the normalizing factor $1/n$ and are , hence, an adequate
generalization.

In order to optimize information matrices, some optimality criterion has to be
employed which is a real valued function of the information matrix and
reflects the main interest in the experiment.

##  7 Optimality criterion based on the failure time under normal use
condition

As in (Weaver and Meeker, 2014) we are interested in some characteristics of
the failure time distribution of soft failure due to degradation (see the
introductory example in Section 2). Therefore it is assumed that the model
equation (3.7)
$\mu_{u}(t)=\mu(\mathbf{x}_{u},t)=(\mathbf{f}_{1}(\mathbf{x}_{u})\otimes\mathbf{f}_{2}(t))^{T}\boldsymbol{\beta}+\mathbf{f}_{2}(t)^{T}\boldsymbol{\gamma}_{u}$
for the mean degradation paths is also valid under normal use condition
$\mathbf{x}_{u}$, where $\mu_{u}$ denotes the degradation path under normal
use condition for short. We further denote by
$\mu(t)=\mathrm{E}(\mu_{u}(t))=(\mathbf{f}_{1}(\mathbf{x}_{u})\otimes\mathbf{f}_{2}(t))^{T}\boldsymbol{\beta}$
the aggregate degradation path under normal use condition and by
$\boldsymbol{\delta}=\boldsymbol{\delta}(\boldsymbol{\beta})=(\delta_{1}(\boldsymbol{\beta}),...,\delta_{p_{2}}(\boldsymbol{\beta}))^{T}$
the vector of its coefficients
$\delta_{s}=\delta_{s}(\boldsymbol{\beta})=\sum_{r=1}^{p_{1}}f_{1r}(\mathbf{x}_{u})\beta_{rs}$,
$s=1,....,p_{2}$, in the regression functions $f_{2s}$ in $t$, i. e.
$\mu(t)=\boldsymbol{\delta}\mathbf{f}_{2}(t)=\sum_{s=1}^{p_{2}}\delta_{s}f_{2s}(t)$.

For the following it is assumed that the mean degradation paths are strictly
increasing over time. Then a soft failure due to degradation is defined as the
exceedance of the degradation over a failure threshold $y_{0}$. This
definition is based on the mean degradation path and not on a “real” path
subject to measurement errors. The failure time $T$ under normal use condition
is then defined as the first time $t$ the mean degradation path $\mu_{u}(t)$
reaches or exceeds the threshold $y_{0}$, i. e. $T=\min\\{t\geq
0;\,\mu_{u}(t)\geq y_{0}\\}$. As random effects $\boldsymbol{\gamma}_{u}$ are
involved in the mean degradation path, the failure time $T$ is random.
Actually, $T$ may become infinite, if the mean degradation path does not reach
the threshold, or may degenerate to $T=0$, if the degradation already exceeds
the threshold at time $t=0$, because of unfortunate values of the random
effects $\boldsymbol{\gamma}_{u}$, but this will happen only with low
probability and will not affect the further argumentation.

In order to describe certain characteristics of the distribution of the
failure time, we will determine the distribution function
$F_{T}(t)=\mathrm{P}(T\leq t)$. First note that $T\leq t$ if and only if
$\mu_{u}(t)\geq y_{0}$. Hence

$\displaystyle F_{T}(t)$ $\displaystyle=$
$\displaystyle\mathrm{P}(\mu_{u}(t)\geq y_{0})$ (7.1) $\displaystyle=$
$\displaystyle\mathrm{P}(\mu(t)+\mathbf{f}_{2}(t)^{T}\boldsymbol{\gamma}_{u}\geq
y_{0})$ $\displaystyle=$
$\displaystyle\mathrm{P}(-\mathbf{f}_{2}(t)^{T}\boldsymbol{\gamma}_{u}\leq\mu(t)-y_{0})$
$\displaystyle=$ $\displaystyle\Phi(h(t)),$

where

$h(t)=\frac{\mu(t)-y_{0}}{\sigma_{u}(t)},$ (7.2)

$\sigma_{u}^{2}(t)=\mathbf{f}_{2}(t)^{T}\boldsymbol{\Sigma}_{\gamma}\mathbf{f}_{2}(t)$
is the variance of the mean degradation path $\mu_{u}(t)$ at time $t$, and
$\Phi$ denotes the distribution function of the standard normal distribution.
Here it is tacitly assumed that the variance $\sigma_{\mu}^{2}(t)$ of the mean
degradation path is greater than zero for every $t\geq 0$. This condition is
satisfied, in particular, when the variance covariance matrix
$\boldsymbol{\Sigma}_{\gamma}$ of the random effects is positive definite.

We will be interested in quantiles $t_{\alpha}$ of the failure time
distribution, i. e. $\mathrm{P}(T\leq t_{\alpha})=\alpha$. For each $\alpha$
the quantile $t_{\alpha}$ gives the time up to which under normal use
conditions (at least) $\alpha\cdot 100$ percent of the units fail and (at
least) $(1-\alpha)\cdot 100$ percent of the units persist. The quantiles
$t_{\alpha}$ are increasing in $\alpha$. Note that this standard definition of
quantiles is in contrast to the“upper” quantiles ($t_{1-\alpha}$) used in
(Weaver and Meeker, 2014) where percentages of failures and persistence are
reversed. Of particular importance is the median $t_{0.5}$ up to which under
normal use conditions half of the units fails and half of the units persist
($\alpha=0.5$). Other characteristics of interest may be the five or ten
percent quantiles $t_{0.05}$ and $t_{0.1}$ which give the times up to which
$95$ or $90\,\%$ percent of the units persist, respectively. By (7.1) these
quantiles can be determined as the solutions of the equation

$h(t_{\alpha})=z_{\alpha},$ (7.3)

where $z_{\alpha}=\Phi^{-1}(\alpha)$ is the $\alpha$-quantile of the standard
normal distribution. For the median ($\alpha=1/2$) we have $z_{0.5}=0$ and,
hence the median failure time $t_{0.5}$ is the solution of $\mu(t)=y_{0}$, i.
e. the aggregate degradation path reaches the threshold at time $t_{0.5}$.
Note that the function $h$ represents the failure time distribution function
$F_{T}$ on a normal Q-Q-plot scale.

In the particular case of straight lines for the mean degradation paths, i. e.
$\mathbf{f}_{2}(t)=(1,t)^{T}$ as considered in Examples 1 and 2, the function
$h(t)$ specifies to

$h(t)=\frac{\delta_{2}t+\delta_{1}-y_{0}}{\sqrt{\sigma_{1}^{2}+2\rho\sigma_{1}\sigma_{2}t+\sigma_{2}^{2}t^{2}}},$
(7.4)

where $\delta_{1}=\sum_{r=1}^{p_{1}}f_{1r}(\mathbf{x}_{u})\beta_{r1}$ and
$\delta_{2}=\sum_{r=1}^{p_{1}}f_{1r}(\mathbf{x}_{u})\beta_{r2}$ are the
intercept and the slope of the aggregate degradation path
$\mu(t)=\delta_{1}+\delta_{2}t$ under normal use condition, respectively. The
median failure time is then given by $t_{0.5}=(y_{0}-\delta_{1})/\delta_{2}$
which provides a proper solution $t_{0.5}>0$ under the natural assumptions
that the aggregate degradation path is increasing, $\delta_{2}>0$, and that
the aggregate degradation at the beginning of the testing at time $t=0$ is
less than the threshold of soft failure, $\delta_{1}<y_{0}$.

######  Example (Example 1 cont.).

In the introductory example the aggregate degradation path
$\mu(t)=\delta_{1}+\delta_{2}t$ has intercept
$\delta_{1}=\beta_{1}+\beta_{2}x_{u}$ and slope
$\delta_{2}=\beta_{3}+\beta_{4}x_{u}$. Hence, the median failure time is given
by $t_{0.5}=(y_{0}-\beta_{1}-\beta_{2}x_{u})/(\beta_{3}+\beta_{4}x_{u})$. If
we use the standardized nominal values of Table 1 for Example 7.2 by (Weaver
and Meeker, 2014), the aggregate degradation path becomes
$\mu(t)=2.306+1.014t$ under normal use condition, and the median failure time
is $t_{0.5}=1.583$. Note that, as typical for degradation experiments, the
median failure time is larger than the maximal experimental time $t_{\max}=1$.

Under the additional assumption that the correlation of the random effects is
non-negative for the intercept and the slope of the mean degradation path,
$\rho\geq 0$, the function $h(t)$ can be seen to be strictly increasing,
$h^{\prime}(t)>0$, in $t>0$. This also remains true for small to moderate
negative correlations. However, the range of $h(t)$ is bounded and does not
cover the whole real line such that not all quantiles are non-degenerate. For
small $\alpha$ the $\alpha$-quantile to be positive requires
$z_{\alpha}>h(0)=-(y_{0}-\delta_{1})/\sigma_{1}$, i. e. the variance
$\sigma_{1}^{2}$ of the intercept of the mean degradation path has to be
sufficiently small compared to the distance from its mean $\delta_{1}$ to the
threshold $y_{0}$. In particular, in the case of the $5\,\%$-quantile
$\sigma_{1}<0.608(y_{0}-\delta_{1})$ is needed for $t_{0.05}>0$. For large
$\alpha$ the $\alpha$-quantile is finite if
$z_{\alpha}<\lim_{t\to\infty}h(t)=\delta_{2}/\sigma_{2}$, i. e. the variance
$\sigma_{2}^{2}$ of the slope of the mean degradation path has to be
sufficiently small compared to its mean $\delta_{2}$. Note that
$\Phi(h(0))=1-\Phi((y_{0}-\delta_{1})/\sigma_{1})$ is the probability that
under normal use condition the mean degradation path exceeds the threshold
$y_{0}$ already at the initial time $t=0$. Note also that formally
$1-\Phi(\lim_{t\to\infty}h(t))=\Phi(-\delta_{2}/\sigma_{2})$ is the
probability that the mean degradation path has a negative slope which may be
interpreted as the probability that soft failure due to degradation will not
occur at all under normal use condition. When the $\alpha$-quantile is non-
degenerate ($0<t_{\alpha}<\infty$), then $t_{\alpha}$ is a solution of the
quadratic equation

$(\delta_{2}t+\delta_{1}-y_{0})^{2}=z_{\alpha}^{2}(\sigma_{1}^{2}+2\rho\sigma_{1}\sigma_{2}t+\sigma_{2}^{2}t^{2}),$

as indicated by (Weaver and Meeker, 2014) for the situation of Example 1. In
the special case of only a random intercept in the random effects, i. e.
$\sigma_{2}^{2}=0$, all $\alpha$-quantiles $t_{\alpha}$ finitely exist for
$\alpha\geq\Phi(-(y_{0}-\delta_{1})/\sigma_{1})$ and can be determined as the
solution of a linear equation to
$t_{\alpha}=(y_{0}-\delta_{1}+z_{\alpha}\sigma_{1})/\delta_{2}$.

######  Example (Example 1 cont.).

For the introductory example the function $h(t)$ is plotted in the right panel
of Figure 2 under the standardized nominal values of (Weaver and Meeker, 2014)
given in Table 1. The defining function $h(t)$ is seen to be strictly
increasing although the correlation is moderately negative ($\rho=-0.143$).
Thus the distribution function $F_{T}(t)=\Phi(h(t))$ is well-defined, and it
is represented in the left panel of Figure 2. In both plots the median failure
time $t_{0.5}=1.583$ is indicated by a dashed vertical line. Moreover, as
$h(0)=-14.03$ and $\lim_{t\to\infty}h(t)=9.67$, the range of $h$ covers all
reasonable quantiles.

Figure 2: Failure time distribution $F_{T}$ (left panel) and defining function
$h$ (right panel) for Example 1

It has to be noted that in the case of the standardized nominal values of
Table 1 for high stress levels ($x_{h}=1$) the mean degradation path $\mu_{i}$
exceeds the threshold $y_{0}$ for soft failure due to degradation with high
probability ($\mathrm{P}(\mu_{i}(1,0)\geq y_{0})>1/2$) already at the initial
experimental time $t_{\min}=0$. Hence, care has to be taken that the model
equation for the mean degradation paths is also valid beyond the threshold, i.
e. in the case that soft failure has already occurred. To avoid this
complication we consider in Example 2 nominal values which guarantee that soft
failure occurs during the experiment only with negligible probability.

######  Example (Example 2 cont.).

For the model with two interacting stress variables $x_{1}$ and $x_{2}$ we
consider the virtual nominal values for the parameters, normal use conditions
and threshold given in Table 2.

Table 2: Nominal values for Example 2 | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | $\beta_{6}$ | $\beta_{7}$ | $\beta_{8}$ | $\sigma_{\gamma}$ | $\sigma_{\varepsilon}$ | ${x}_{u1}$ | ${x}_{u2}$ | $y_{0}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$f_{j}(\mathbf{x},t)$ | $1$ | $x_{1}$ | $x_{2}$ | $x_{1}x_{2}$ | $t$ | $x_{1}t$ | $x_{2}t$ | $x_{1}x_{2}t$ |  |  |  |  |
| $4.0$ | $1.5$ | $0.75$ | $1.8$ | $0.5$ | $0.25$ | $0.25$ | $4.03$ | $0.7$ | $0.85$ | $-0.5$ | $-0.4$ | $14.39$

The aggregate degradation path $\mu(t)$ is a straight line with intercept
$\delta_{1}=\beta_{1}+\beta_{2}x_{u1}+\beta_{3}x_{u2}+\beta_{4}x_{u1}x_{u2}=3.31$
and slope
$\delta_{2}=\beta_{5}+\beta_{6}x_{u1}+\beta_{7}x_{u2}+\beta_{8}x_{u1}x_{u2}=1.08$,
i. e. $\mu(t)=3.31+1.08t$ under normal use condition. With a threshold of
$y_{0}=14.39$ for soft failure the median failure time results in
$t_{0.5}=(y_{0}-\delta_{1})/\delta_{2}=10.25$ which is substantially larger
than the maximal experimental time $t_{\max}=1$. For the characterization of
other quantiles the function $h(t)$ is plotted in the right panel of Figure 3
together with the corresponding distribution function $F_{T}(t)=\Phi(h(t))$
(left panel). The median failure time $t_{0.5}=10.25$ is indicated in both
plots by a dashed vertical line. As $\rho=0$ the function $h(t)$ is strictly
increasing and ranges from $h(0)=-15.83$ to
$h_{\max}=\lim_{t\to\infty}h(t)=1.54$. Thus, quantiles $t_{\alpha}$ are non-
degenerate as long as $\alpha\leq\alpha_{\max}$, where
$\alpha_{\max}=\Phi(h_{\max})=0.939$, and $(1-\alpha_{\max})\cdot 100=6.1$
percent of the mean degradation paths do not lead to a soft failure. Both
$\alpha_{\max}$ and $h_{\max}$ are indicated in the respective plots by a
dashed horizontal line.

Figure 3: Failure time distribution $F_{T}$ (left panel) and defining function
$h$ (right panel) for Example 2

Note that for the nominal values of Table 2 the mean degradation
$\mu(\mathbf{x},t)$ under experimental conditions attains its maximum $13.08$
for the maximal stress levels ($x_{1}=x_{2}=1$) and maximal experimental time
($t=1$), and, hence, the mean degradation paths do not exceed the threshold
for all experimental settings.

In any case the quantile $t_{\alpha}=t_{\alpha}(\boldsymbol{\theta})$ is a
function of both the aggregate location parameters $\boldsymbol{\beta}$ and
the variance parameters $\boldsymbol{\varsigma}$, in general. Hence, the
maximum likelihood estimator of the quantile $t_{\alpha}$ is given by
$\widehat{t}_{\alpha}=t_{\alpha}(\widehat{\boldsymbol{\theta}})$ in terms of
the maximum likelihood estimator $\widehat{\boldsymbol{\theta}}$ of
$\boldsymbol{\theta}$. The task of designing the experiment will now be to
provide an as precise estimate of the $\alpha$-quantile as possible.

By the delta-method $\widehat{t}_{\alpha}$ is seen to be asymptotically normal
with asymptotic variance

$\mathrm{aVar}(\widehat{t}_{\alpha})=\mathbf{c}^{T}\mathbf{M}_{\boldsymbol{\theta}}^{-1}\mathbf{c},$
(7.5)

where $\mathbf{c}=\frac{\partial}{\partial\boldsymbol{\theta}}t_{\alpha}$ is
the gradient vector of partial derivatives of $t_{\alpha}$ with respect to the
components of the parameter vector ${\boldsymbol{\theta}}$. The asymptotic
variance depends on the design of the experiment through the information
matrix $\mathbf{M}_{\boldsymbol{\theta}}$ and will be chosen as the optimality
criterion for the design.

The gradient $\mathbf{c}$ can be seen to be equal to

$\mathbf{c}=-c_{0}\left({\textstyle{\frac{\partial}{\partial\boldsymbol{\theta}}}}\mu(t)|_{t=t_{\alpha}}-z_{\alpha}{\textstyle{\frac{\partial}{\partial\boldsymbol{\theta}}}}\sigma_{u}(t)|_{t=t_{\alpha}}\right),$
(7.6)

in view of (7.2) and (7.3) by the implicit function theorem (see e. g. (Krantz
and Parks, 2012)), where
$c_{0}=1/(\mu^{\prime}(t_{\alpha})-z_{\alpha}\sigma_{u}^{\prime}(t_{\alpha}))$
is the inverse of the derivative of the defining function $h$ with respect to
$t$.

As the aggregate mean degradation $\mu(t)$ only depends on the aggregate
location parameters $\boldsymbol{\beta}$ and the variance $\sigma_{u}^{2}(t)$
only depends on the variance parameters $\boldsymbol{\varsigma}$ the gradient
simplifies to
$\mathbf{c}=-c_{0}(\mathbf{c}_{\boldsymbol{\beta}}^{T},\mathbf{c}_{\boldsymbol{\varsigma}}^{T})^{T}$,
where

$\mathbf{c}_{\boldsymbol{\beta}}={\textstyle{\frac{\partial}{\partial\boldsymbol{\beta}}}}\mu(t)|_{t=t_{\alpha}}=\mathbf{f}(\mathbf{x}_{u},t_{\alpha})$

is the gradient of $\mu(t)$ with respect to $\boldsymbol{\beta}$ and

$\mathbf{c}_{\boldsymbol{\varsigma}}=-z_{\alpha}{\textstyle{\frac{\partial}{\partial\boldsymbol{\varsigma}}}}\sigma_{u}(t)|_{t=t_{\alpha}}$

is $-z_{\alpha}$ times the gradient of $\sigma_{u}(t)$ with respect to
$\boldsymbol{\varsigma}$. The particular shape of
$\mathbf{c}_{\boldsymbol{\varsigma}}$ does not play a role here, in general.
But note that $\mathbf{c}_{\boldsymbol{\varsigma}}=\mathbf{0}$ in the case of
the median ($\alpha=0.5$).

By the block diagonal form (5.1) of the information matrix the asymptotic
variance (7.5) of $\widehat{t}_{\alpha}$ becomes

$\mathrm{aVar}(\widehat{t}_{\alpha})=c_{0}^{2}\left(\mathbf{c}_{\boldsymbol{\beta}}^{T}\mathbf{M}_{\boldsymbol{\beta}}^{-1}\mathbf{c}_{\boldsymbol{\beta}}+\mathbf{c}_{\boldsymbol{\varsigma}}^{T}\mathbf{M}_{\boldsymbol{\varsigma}}^{-1}\mathbf{c}_{\boldsymbol{\varsigma}}\right)$
(7.7)

which simplifies to

$\mathrm{aVar}(\widehat{t}_{0.5})=c_{0}^{2}\mathbf{c}_{\boldsymbol{\beta}}^{T}\mathbf{M}_{\boldsymbol{\beta}}^{-1}\mathbf{c}_{\boldsymbol{\beta}}$
(7.8)

in the case of the median.

For the product-type model (3.6) the expression related to the aggregate
parameters $\boldsymbol{\beta}$ further decomposes,

$\mathbf{c}_{\boldsymbol{\beta}}^{T}\mathbf{M}_{\boldsymbol{\beta}}^{-1}\mathbf{c}_{\boldsymbol{\beta}}=\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})\cdot\mathbf{f}_{2}(t_{\alpha})^{T}\mathbf{M}_{2}^{-1}\mathbf{f}_{2}(t_{\alpha}),$
(7.9)

as
$\mathbf{c}_{\boldsymbol{\beta}}=\mathbf{f}_{1}(\mathbf{x}_{u})\otimes\mathbf{f}_{2}(t_{\alpha})$
factorizes.

##  8 Optimal designs in the case of predetermined measurement times

From (7.7) and (7.9) it can be seen that for obtaining a minimal asymptotic
variance for $\widehat{t}_{\alpha}$ only
$\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})$
has to be minimized, because all other terms do not depend on the experimental
settings $\mathbf{x}_{1},...,\mathbf{x}_{n}$ of the stress variable, when the
measurement times $t_{1},...,t_{k}$ are predetermined. The optimality
criterion of minimization of the asymptotic variance of $\widehat{t}_{\alpha}$
thus reduces to a $c$-criterion
$\mathbf{c}_{1}^{T}\mathbf{M}_{1}(\xi)^{-1}\mathbf{c}_{1}$ for extrapolation
of the marginal response at normal use condition $\mathbf{x}_{u}$ in the first
marginal model (5.4), $\mathbf{c}_{1}=\mathbf{f}_{1}(\mathbf{x}_{u})$, which
is a well-known problem from the literature (see (Kiefer and Wolfowitz,
1964)). It is remarkable that this criterion and, hence, the corresponding
optimal design is the same whatever the value of $\alpha$ is, as long as there
is a proper solution $0<t_{\alpha}<\infty$ for the $\alpha$-quantile of the
failure time.

######  Proposition 8.1.

If the design $\xi^{*}$ is $c$-optimal for extrapolation of the mean response
at the normal use condition in the marginal model (5.4) for the stress
variable, then $\xi^{*}$ minimizes the asymptotic variance for the estimator
$\widehat{t}_{\alpha}$ of the $\alpha$-quantile of the failure time for every
$\alpha$ when $0<t_{\alpha}<\infty$ (for predetermined measurement times
$t_{1},...,t_{k}$).

Although the normal use condition is typically outside the experimental
region, the above proposition also would hold for interpolation, i. e.
$\mathbf{x}_{u}\in\mathcal{X}$. The result of Proposition 8.1 is next used to
derive optimal designs for the situation in Examples 1 and 2.

######  Example (Example 1 cont.).

In the introductory model of Section 2 the marginal model for the stress
variable $x$ is given by a simple linear regression,
$\mathbf{f}_{1}(x)=(1,x)^{T}$. In this marginal model the $c$-optimal design
$\xi^{*}$ for extrapolation of the mean response
$\mu^{(1)}(x_{u})=\beta_{1}^{(1)}+\beta_{2}^{(1)}x_{u}$ under normal use
condition $x_{u}<0$ assigns weight $w^{*}=|x_{u}|/(1+2|x_{u}|)$ to the highest
stress level $x_{h}=1$ and weight $1-w^{*}=(1+|x_{u}|)/(1+2|x_{u}|)$ to the
lowest stress level $x_{l}=0$ on the standardized scale $\mathcal{X}=[0,1]$
(see (Kiefer and Wolfowitz, 1964)). Note that larger weight $1-w^{*}>w^{*}$ is
assigned to the lowest stress level $x_{l}=0$ which is closer to $x_{u}<0$
than $x_{h}=1$ and that the weight $1-w^{*}$ at $x_{l}$ decreases from $1$ to
$1/2$, when the distance between the normal use condition and the experimental
region gets larger, i. e. $x_{u}$ decreases. For the standardized value
$x_{u}=-0.056$ of the normal use condition from Table 1 the optimal weights
for extrapolation at $x_{u}=-0.056$ are $w^{*}=0.05$ at $x_{h}=1$ and
$1-w^{*}=0.95$ at $x_{l}=0$, and the optimal design is

$\xi^{*}=\left(\begin{array}[]{cc}0&1\\\ 0.95&0.05\end{array}\right).$

Further examples for extrapolation at $x_{u}=-0.4$, $-0.5$, and $-1$ give
optimal weights $w^{*}=0.22$, $0.25$, and $0.33$ at $x_{h}=1$, and
$1-w^{*}=0.78$, $0.75$, and $0.67$ at $x_{l}=0$, respectively.

By Proposition 8.1 the design $\xi^{*}$ is also optimal for minimization of
the asymptotic variance for estimating the $\alpha$-quantile $t_{\alpha}$ of
the failure time for soft failure due to degradation under normal use
condition $x_{u}$, when $0<t_{\alpha}<\infty$ and the measurement times
$t_{1},...,t_{k}$ are predetermined (see (Schwabe et al., 2014) for estimation
of the median, $\alpha=0.5$). In particular, for the standardized value
$x_{u}=-0.056$ of the normal use condition from Table 1 the optimal design for
estimating any $\alpha$-quantile $t_{\alpha}$ assigns weight $0.95$ to
$x_{l}=0$ and weight $0.05$ to $x_{h}=1$, as found numerically by (Weaver and
Meeker, 2014) in the case of the median $t_{0.5}$.

The so obtained designs for one stress variable can now be used to construct
$c$-optimal designs in the presence of two stress variables with interactions.

######  Example (Example 2 cont.).

In the model with two interacting stress variables $x_{1}$ and $x_{2}$ the
marginal model for the combined stress variable $\mathbf{x}=(x_{1},x_{2})$ is
given itself by a product-type structure,
$\mathbf{f}_{1}(\mathbf{x})=\mathbf{f}_{11}(x_{1})\otimes\mathbf{f}_{12}(x_{2})$,
where both components $x_{1}$ and $x_{2}$ are specified by
$\mathbf{f}_{1v}(x_{v})=(1,x_{v})^{T}$ as simple linear regressions in their
corresponding submarginal models
$Y_{i}^{(1v)}=\beta_{1}^{(1v)}+\beta_{2}^{(1v)}x_{v}+\varepsilon_{i}^{(1v)}$,
$v=1,2$, with standardized homoscedastic and uncorrelated error terms.
Moreover, the experimental region $\mathcal{X}=[0,1]^{2}$ for the combined
stress variable $\mathbf{x}$ is the Cartesian product of the marginal
experimental regions $\mathcal{X}_{1}=\mathcal{X}_{2}=[0,1]$ for the
components $x_{1}$ and $x_{2}$, respectively. The vector $\mathbf{c}$ for
extrapolation of the mean response
$\mathbf{c}^{T}\boldsymbol{\beta}^{(1)}=\mu^{(1)}(\mathbf{x}_{u})=\beta_{1}^{(1)}+\beta_{2}^{(1)}x_{u1}+\beta_{3}^{(1)}x_{u2}+\beta_{4}^{(1)}x_{u1}x_{u2}$
under normal use condition $\mathbf{x}_{u}=(x_{u1},x_{u2})$,
$x_{u1},x_{u2}<0$, is given by
$\mathbf{f}_{1}(\mathbf{x}_{u})=\mathbf{f}_{11}(x_{u1})\otimes\mathbf{f}_{12}(x_{u2})$
and, hence, also factorizes as
$\mathbf{c}=\mathbf{c}_{1}\otimes\mathbf{c}_{2}$, where
$\mathbf{c}_{v}=\mathbf{f}_{1v}(x_{uv})$. In this setting the $c$-optimal
design $\xi^{*}$ for extrapolation at $\mathbf{x}_{u}$ can be obtained as the
product $\xi^{*}=\xi_{1}^{*}\otimes\xi_{2}^{*}$ of the $c$-optimal designs
$\xi_{v}^{*}$ for extrapolation at $x_{uv}$ in the submarginal models (see
Theorem 4.4 in (Schwabe, 1996b)).

The submarginal $c$-optimal designs $\xi_{v}^{*}$ can be derived as in Example
1. They assign weight $w_{v}^{*}=|x_{uv}|/(1+2|x_{uv}|)$ to $x_{v}=1$ and
weight $1-w_{v}^{*}=(1+|x_{uv}|)/(1+2|x_{uv})$ to $x_{v}=0$. Hence, the
$c$-optimal design $\xi^{*}=\xi_{1}^{*}\otimes\xi_{2}^{*}$ for extrapolation
at $\mathbf{x}_{u}$ is given by

$\xi^{*}=\left(\begin{array}[]{cccc}(0,0)&(0,1)&(1,0)&(1,1)\\\
(1-w_{1}^{*})(1-w_{2}^{*})&(1-w_{1}^{*})w_{2}^{*}&w_{1}^{*}(1-w_{2}^{*})&w_{1}^{*}w_{2}^{*}\end{array}\right).$

Then, by Proposition 8.1, the design $\xi^{*}$ is also optimal for
minimization of the asymptotic variance for estimating the $\alpha$-quantile
$t_{\alpha}$ of the failure time for soft failure due to degradation, when
$0<t_{\alpha}<\infty$ and the measurement times $t_{1},...,t_{k}$ are
predetermined. For example, when the normal use conditions are ${x}_{u1}=-0.5$
for the first component and ${x}_{u2}=-0.4$ for the second component as
specified in Table 2, then by the results in Example 1 the optimal marginal
weights are $w_{1}^{*}=0.25$ and $w_{2}^{*}=0.22$, and the optimal design
$\xi^{*}=\xi_{1}^{*}\otimes\xi_{2}^{*}$ is given by

$\xi^{*}=\left(\begin{array}[]{cccc}(0,0)&(0,1)&(1,0)&(1,1)\\\
0.58&0.17&0.19&0.06\end{array}\right).$

The $c$-optimal extrapolation designs can be obtained in both Examples 1 and 2
by Elfving’s theorem ((Elfving, 1952)) which provides a geometrical
construction of a $c$-optimal design (see (Schwabe, 1996b), Theorem 2.13). To
give a rough idea of this construction one has to consider the Elfving set
which is the convex hull

$\mathcal{E}=\mathrm{conv}(\\{\mathbf{f}(\mathbf{x});\,\mathbf{x}\in\mathcal{X}\\}\cup\\{-\mathbf{f}(\mathbf{x});\,\mathbf{x}\in\mathcal{X}\\})$

of the union of the so-called induced design region
$\\{\mathbf{f}(\mathbf{x});\,\mathbf{x}\in\mathcal{X}\\}$ and its image
$\\{-\mathbf{f}(\mathbf{x});\,\mathbf{x}\in\mathcal{X}\\}$ under reflection at
the origin $\mathbf{0}$ in $\mathbb{R}^{p}$. Here $\mathbf{x}$ and
$\mathbf{f}$ denote variables and regression functions associated with a
generic model
$Y_{i}=\mathbf{f}(\mathbf{x}_{i})^{T}\boldsymbol{\beta}+\varepsilon_{i}$. In
Example 1 we have $\mathbf{x}=x$, $\mathbf{f}(\mathbf{x})=(1,x)^{T}$ and
$\mathcal{X}=[0,1]$, and the Elfving set $\mathcal{E}$ is given as a
parallelogram in $\mathbb{R}^{2}$ with one edge from $(1,0)^{T}$ to
$(1,1)^{T}$ representing the induced design region and the opposite edge from
$(-1,0)^{T}$ to $(-1,-1)^{T}$ representing its image under reflection. For
another illustration of an Elfving set see Figure 4 in Section 10, where an
optimal design is sought for the time variable $t$.

The $c$-optimal design for estimating $\mathbf{c}^{T}\boldsymbol{\beta}$ can
then be constructed as follows: Determine the intersection point of the ray
$\lambda\mathbf{c}$, $\lambda>0$, with the boundary of the Elfving set,
$\lambda_{c}\mathbf{c}$ say. This point can be represented as a convex
combination of (extremal) points $\pm\mathbf{f}(\mathbf{x}_{i})$ of the
induced design region and its reflection,

$\lambda_{c}\mathbf{c}=\sum_{i=1}^{m}w_{i}z_{i}\mathbf{f}(\mathbf{x}_{i}),$

where $z_{i}=1$, when the (extremal) point $\mathbf{f}(\mathbf{x}_{i})$ is
from the induced design region, and $z_{i}=-1$, when the point
$-\mathbf{f}(\mathbf{x}_{i})$ is from the reflection, and the weights $w_{i}$
of the convex combination satisfy $w_{i}>0$ and $\sum_{i=1}^{m}w_{i}=1$. Then
Elfving’s theorem states that the design $\xi^{*}$ which assigns weights
$w_{i}$ to the settings $\mathbf{x}_{i}$ is $c$-optimal (for $\mathbf{c}$).
Moreover, this construction also provides the value of the $c$-criterion,
$\mathbf{c}^{T}\mathbf{M}(\xi^{*})^{-1}\mathbf{c}=1/\lambda_{c}^{2}$.

In Example 1 the ray $\lambda(1,x_{u})^{T}$ intersects the boundary of the
Elfving set $\mathcal{E}$ at the connecting line from
$(1,0)^{T}=\mathbf{f}(0)$ to $(-1,-1)^{T}=-\mathbf{f}(1)$ at
$\lambda_{c}(1,x_{u})^{T}=w_{1}(-1,-1)^{T}+w_{2}(1,0)^{T}$ with
$\lambda_{c}=1/(1-2x_{u})$, $w_{1}=-x_{u}/(1-2x_{u})>0$ and
$w_{2}=1-w_{1}=(1-x_{u})/(1-2x_{u})>0$. Hence, the optimality of the given
design follows.

We will use Elfving’s theorem next to characterize optimal designs for the
situation with two non-interacting stress variables.

######  Example 3.

In the case of two non-interacting stress variables $x_{1}$ and $x_{2}$ we
consider the model equation

$y_{ij}=\beta_{i,1}+\beta_{2}x_{i1}+\beta_{3}x_{i2}+\beta_{i,4}t_{j}+\beta_{5}x_{i1}t_{j}+\beta_{6}x_{i2}t_{j}+\varepsilon_{ij}$
(8.1)

for the combined stress variable $\mathbf{x}=(x_{1},x_{2})$ and the time
variable $t$. This model contains all terms of the full interaction model 3.8)
with the exception of the terms $x_{1}x_{2}$ and $x_{1}x_{2}t$ related to
potential interactions between the stress variables. The interpretation of all
other terms in (8.1) is the same as in Example 2. Model (8.1) is constructed
from the marginal model

$Y_{i}^{(1)}=\beta_{1}^{(1)}+\beta_{2}^{(1)}x_{i1}+\beta_{3}^{(1)}x_{i2}+\varepsilon_{i}^{(1)}$

which is additive in the effects of the stress variables $x_{1}$ and $x_{2}$,
i. e. $\mathbf{f}_{1}(\mathbf{x})=(1,x_{1},x_{2})^{T}$.

For this marginal model of multiple regression, the Elfving set is an oblique
prism with quadratic base with vertices $(1,0,0)^{T}$, $(1,0,1)^{T}$,
$(1,1,0)^{T}$, $(1,1,1)^{T}$ and quadratic top with vertices $(-1,0,0)^{T}$,
$(-1,0,-1)^{T}$, $(-1,-1,0)^{T}$, $(-1,-1,-1)^{T}$. To find the $c$-optimal
extrapolation design at the normal use condition
$\mathbf{x}_{u}=(x_{u1},x_{u2})$ by Elfving’s theorem we have to determine the
intersection point of the ray $\lambda(1,x_{u1},x_{u2})^{T}$ with the surface
of the Elfving set. For $x_{u1}<x_{u2}<0$ the ray intersects the surface at
the quadrangular face of the prism spanned by $(1,0,0)^{T}$, $(1,0,1)^{T}$,
$(-1,-1,0)^{T}$, and $(-1,-1,-1)^{T}$ when $\lambda_{c}=1/(1+2|x_{u1}|)$. The
representation of the intersection point $\lambda_{c}(1,x_{u1},x_{u2})^{T}$ by
the vertices of the quadrangle is not unique. There are two $c$-optimal
designs

$\xi_{0}^{*}=\left(\begin{array}[]{ccc}(0,0)&(0,1)&(1,1)\\\
(1+|x_{u2}|)\lambda_{c}&(|x_{u1}|-|x_{u2}|)\lambda_{c}&|x_{u1}|\lambda_{c}\end{array}\right)$

and

$\xi_{1}^{*}=\left(\begin{array}[]{ccc}(0,0)&(1,0)&(1,1)\\\
(1+|x_{u1}|)\lambda_{c}&(|x_{u1}|-|x_{u2}|)\lambda_{c}&|x_{u2}|\lambda_{c}\end{array}\right)$

which are supported on three vertices. As a consequence, also for all
coefficients $a$, $0<a<1$, the convex combination

$\displaystyle\xi_{a}^{*}$ $\displaystyle=(1-a)\xi_{0}^{*}+a\xi_{1}^{*}$
$\displaystyle=\mbox{\scriptsize{$\left(\begin{array}[]{cccc}(0,0)&(0,1)&(1,0)&(1,1)\\\
(1+a|x_{u1}|+(1-a)|x_{u2}|)\lambda_{c}&(1-a)(|x_{u1}|-|x_{u2}|)\lambda_{c}&a(|x_{u1}|-|x_{u2}|)\lambda_{c}&((1-a)|x_{u1}|+a|x_{u2}|)\lambda_{c}\end{array}\right)$}}$

supported on all four vertices is $c$-optimal for extrapolation at
$\mathbf{x}_{u}$, $0<t_{\alpha}<\infty$. Then, by Proposition 8.1, the designs
$\xi_{a}^{*}$ are also optimal for minimization of the asymptotic variance for
estimating the $\alpha$-quantile $t_{\alpha}$ of the failure time for soft
failure due to degradation, when $0<t_{\alpha}<\infty$ and the measurement
times $t_{1},...,t_{k}$ are predetermined, $0\leq\alpha\leq 1$. For example,
when the normal use conditions are $x_{u1}=-0.5$ for the first component and
$x_{u2}=-0.4$ for the second component as in Example 2, then the optimal
design $\xi_{a}^{*}$ is given by

$\xi_{a}^{*}=\left(\begin{array}[]{cccc}(0,0)&(0,1)&(1,0)&(1,1)\\\
0.70+0.05a&0.05-0.05a&0.05a&0.25-0.05a\end{array}\right)$

with the special cases

$\xi_{0}^{*}=\left(\begin{array}[]{ccc}(0,0)&(0,1)&(1,1)\\\
0.70&0.05&0.25\end{array}\right)\qquad\mathrm{and}\qquad\xi_{1}^{*}=\left(\begin{array}[]{ccc}(0,0)&(1,0)&(1,1)\\\
0.75&0.05&0.20\end{array}\right)$

supported on three vertices.

Note that there are also other designs which are $c$-optimal for extrapolation
at $\mathbf{x}_{u}$, but which are not solely supported on the vertices. For
example, for $x_{u1}<x_{u2}<0$ the two-point design which assigns weight
$w=|x_{u1}|/(1+2|x_{u1}|)$ to $(0,0)$ and weight
$1-w=(1+|x_{u1}|)/(1+2|x_{u1}|)$ to $(1,x_{u2}/x_{u1})$ is $c$-optimal by
Elfving’s theorem. However, these designs can be used for estimating
$t_{\alpha}$ by means of maximum-likelihood only when the resulting
information matrix is non-singular, i. e. when the design has, at least, three
distinct support points.

For $x_{u2}<x_{u1}<0$ optimal designs can be obtained from the above case by
interchanging the roles of the two components $x_{1}$ and $x_{2}$.

In the case $x_{u1}=x_{u2}<0$ there is only one $c$-optimal design for
extrapolation. This design is supported on two vertices and assigns weight
$w=|x_{u1}|/(1+2|x_{u1}|)$ to $(0,0)$ and weight
$1-w=(1+|x_{u1}|)/(1+2|x_{u1}|)$ to $(1,1)$. As the resulting information
matrix is singular, this design cannot be used for estimating the
$\alpha$-quantile $t_{\alpha}$ of the failure time for soft failure due to
degradation. Hence, no suitable optimal design exists in this case, but the
$c$-optimal design may serve as a benchmark for judging the quality of a
competing design in terms of efficiency.

To quantify the quality of a standard design $\xi_{0}$ for estimating the
quantile $t_{\alpha}$ of the mean failure time under normal use condition we
make use of the efficiency

$\mathrm{eff}_{\mathrm{aVar}}(\xi_{0})=\frac{\mathrm{aVar}(\widehat{t}_{\alpha};\xi^{*})}{\mathrm{aVar}(\widehat{t}_{\alpha};\xi_{0})},$
(8.2)

where
$\mathrm{aVar}(\widehat{t}_{\alpha};\xi)=\mathbf{c}^{T}\mathbf{M}_{\boldsymbol{\theta}}(\xi)^{-1}\mathbf{c}$
denotes the standardized asymptotic variance for estimating $t_{\alpha}$ by
equation (6.4) when design $\xi$ is used, and $\xi^{*}$ is the corresponding
optimal design. The efficiency gives the proportion of units to be used under
the optimal design $\xi^{*}$ which provides (asymptotically) the same accuracy
(in terms of the asymptotic variance) compared to the standard design
$\xi_{0}$. For example, if the efficiency is $0.5$ twice the number of units
have to be used under $\xi_{0}$ than under the optimal design $\xi^{*}$ to get
the same accuracy. Note that both the asymptotic variance and the efficiency
may also depend on the parameter vector $\boldsymbol{\theta}$, at least,
through $t_{\alpha}$ and are, hence, local quantities (at
$\boldsymbol{\theta}$) without explicitly stated in the notation.

In the case of estimating the median $t_{0.5}$ the standardized asymptotic
variance factorizes as

$\mathrm{aVar}(\widehat{t}_{0.5};\xi)=\frac{1}{n}c_{0}^{2}\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}(\xi)^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})\cdot\mathbf{f}_{2}(t_{\alpha})^{T}\mathbf{M}_{2}^{-1}\mathbf{f}_{2}(t_{\alpha})$
(8.3)

by equations (7.8) and (7.9) for the general product-type model (3.6). Thus
the efficiency defined in (8.2) reduces to the $c$-efficiency

$\mathrm{eff}_{c}(\xi_{0})=\frac{\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}(\xi^{*})^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})}{\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}(\xi_{0})^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})}$

for extrapolation at the normal use condition $\mathbf{x}_{u}$ in the first
marginal model with uncorrelated homoscedastic errors and does not depend on
$\boldsymbol{\theta}$.

######  Example (Example 1 cont.).

In the introductory model of Section 2 the $c$-criterion for extrapolation at
$x_{u}$ is defined by
$\Phi_{c}(\xi)=\mathbf{f}_{1}(x_{u})^{T}\mathbf{M}_{1}(\xi)^{-1}\mathbf{f}_{1}(x_{u})$.
For $x_{u}<0$ it attains its minimal value value
$\Phi_{c}(\xi^{*})=(1+2|x_{u}|)^{2}$ for the optimal design $\xi^{*}$ which
assigns weight $w^{*}=|x_{u}|/(1+2|x_{u}|)$ to $x_{h}=1$ and weight
$1-w^{*}=(1+|x_{u}|)/(1+2|x_{u}|)$ to $x_{l}=0$. Common alternatives would be
uniform designs $\bar{\xi}_{m}$ which assign equal weights $w=1/m$ to $m$
experimental settings on an equidistant grid
$\\{x_{1},...,x_{m}\\}=\\{0,1/(m-1),...,1\\}$ of the experimental region
$\mathcal{X}=[0,1]$. For these designs the $c$-criterion for extrapolation at
$x_{u}<0$ can be calculated as
$\Phi_{c}(\bar{\xi}_{m})=1+a_{m}(1+2|x_{u}|)^{2}$, where $a_{m}=3(m-1)/(m+1)$.
Their $c$-efficiency for extrapolation at $x_{u}$ and, hence, their efficiency
for estimating the median failure time under normal use condition $x_{u}$ is
equal to
$\mathrm{eff}_{c}(\bar{\xi}_{m})=\Phi_{c}(\xi^{*})/\Phi_{c}(\bar{\xi}_{m})=(1-1/(1+a_{m}(1+2|x_{u}|)^{2}))/a_{m}$
which increases from $1/(1+a_{m})=(m+1)/(4m-2)$ for $x_{u}$ close to the
lowest stress level $x_{l}=0$ to $1/a_{m}=(m+1)/(3m-3)$ when $x_{u}$ tends to
minus infinity. Moreover, for fixed $x_{u}$, the efficiency decreases when $m$
increases, i. e. when the grid becomes more dense. For selected values of the
normal use condition $x_{u}$ and numbers $m$ of grid points numerical values
of the efficiency are reported in Table 3

Table 3: Efficiency of uniform designs $\bar{\xi}_{m}$ for various normal use conditions $x_{u}$ in Example 1 $m$ | $x_{u}$
---|---
| $\ \quad 0\quad\ $ | $-0.056$ | $-0.400$ | $-0.500$ | $-1.000$ | $\quad-\infty\quad$
$2$ | $0.50$ | $0.55$ | $0.76$ | $0.80$ | $0.90$ | $1.00$
$3$ | $0.40$ | $0.43$ | $0.55$ | $0.57$ | $0.62$ | $0.67$
$4$ | $0.36$ | $0.38$ | $0.47$ | $0.49$ | $0.52$ | $0.56$
$5$ | $0.33$ | $0.36$ | $0.43$ | $0.44$ | $0.47$ | $0.50$
$\infty$ | $0.25$ | $0.26$ | $0.30$ | $0.31$ | $0.32$ | $0.33$

Note that in Table 3 the row $m=\infty$ corresponds to a continuous uniform
design as an approximation to large numbers $m$ of grid points, while the
columns $x_{u}=0$ and $x_{u}=-\infty$ give approximations for normal use
conditions $x_{u}$ close to the lowest experimental stress level or far away,
respectively.

For the particular case $m=2$, where the design $\bar{\xi}_{2}$ assigns equal
weights $w=1-w=1/2$ to both the highest and the lowest stress level $x_{h}=1$
and $x_{l}=0$, we have $a_{2}=1$ and, hence,
$\Phi_{c}(\xi_{0})=1+(1+2|x_{u}|)^{2}$ for the $c$-criterion. The
$c$-efficiency of $\bar{\xi}_{2}$ for extrapolation at $x_{u}$ and, thus, its
efficiency for estimating the median failure time under normal use condition
$x_{u}$ is equal to $\mathrm{eff}_{c}(\bar{\xi}_{2})=1-1/(1+(1+2|x_{u}|)^{2})$
which ranges from $1/2$ for $x_{u}$ close to the lowest stress level $x_{l}=0$
to $1$ when $x_{u}$ tends to minus infinity.

For the nominal value $x_{u}=-0.056$ of the normal use condition in Table 1
the efficiency of the equidistant grid designs $\bar{\xi}_{m}$ is reported in
the third column of Table 3. In particular, for the uniform design
$\bar{\xi}_{2}$ on the endpoints of the experimental region this efficiency is
$0.55$ which means that
$\mathrm{eff}_{c}(\bar{\xi}_{2})^{-1}-1=1/(1+2|x_{u}|)^{2}=81\,\%$ more units
would have to be used for design $\bar{\xi}_{2}$ to obtain the same quality
for estimating the median failure time than for the optimal design $\xi^{*}$.

The efficiencies in the case of one stress variable can be used to compute the
efficiency in the presence of two stress variables with interactions.

######  Example (Example 2 cont.).

In the model with two interacting stress variables $x_{1}$ and $x_{2}$ the
$c$-criterion
$\Phi_{c}(\xi)=\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}(\xi)^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})$
for extrapolation at $\mathbf{x}_{u}=(x_{u1},x_{u2})$ factorizes into its
counterparts in the submarginal models,
$\Phi_{c}(\xi)=\mathbf{f}_{11}({x}_{u1})^{T}\mathbf{M}_{11}(\xi_{1})^{-1}\mathbf{f}_{11}({x}_{u1})\cdot\mathbf{f}_{12}({x}_{u2})^{T}\mathbf{M}_{12}(\xi_{2})^{-1}\mathbf{f}_{12}({x}_{u2})$.
Because also the $c$-optimal design $\xi^{*}=\xi_{1}^{*}\otimes\xi_{2}^{*}$
has product-type structure, the $c$-efficiency for extrapolation at
$\mathbf{x}_{u}$ and, hence, the efficiency for estimating the median failure
time factorizes,
$\mathrm{eff}_{c}(\xi_{1}\otimes\xi_{2})=\mathrm{eff}_{c1}(\xi_{1})\cdot\mathrm{eff}_{c2}(\xi_{2})$,
where $\mathrm{eff}_{cv}(\xi_{v})$ is the corresponding efficiency in the
$v$th submarginal model, $v=1,2$.

The design $\bar{\xi}$ which assigns equal weights $1/4$ to the four vertices
$(0,0)$, $(0,1)$, $(1,0)$, and $(1,1)$ of the experimental region serves as a
natural standard design. This design can be seen to be the product
$\bar{\xi}=\bar{\xi}_{2}\otimes\bar{\xi}_{2}$ of submarginal designs
$\bar{\xi}_{2}$ which assign equal weights $1/2$ to the lowest and highest
stress level $x_{vl}$and $x_{vh}$ in the submarginal models, $v=1,2$. Hence,
from Example 1 we get the efficiency of $\bar{\xi}$ as
$\mathrm{eff}_{c}(\bar{\xi})=((1+2|x_{u1}|)(1+2|x_{u2}|))^{2}/((1+(1+2|x_{u1}|)^{2})(1+(1+2|x_{u2}|)^{2}))$
which ranges from $1/4$ for $\mathbf{x}_{u}$ close to the combination
$(x_{1l},x_{2l})$ of lowest stress levels $x_{1l}=0$ and $x_{2l}=0$ to $1$
when both normal use conditions $x_{u1}$ and $x_{u2}$ tend to minus infinity.

For example, when the normal use conditions are ${x}_{u1}=-0.5$ for the first
component and ${x}_{u2}=-0.4$ for the second component as specified in Table
2, then according to Table 3 the efficiency $\mathrm{eff}_{cv}(\bar{\xi}_{2})$
of $\bar{\xi}_{2}$ is $0.80$ and $0.76$ in the respective submarginal models,
$v=1,2$. By the above considerations the efficiency of $\bar{\xi}$ is
$\mathrm{eff}_{c}(\bar{\xi})=0.80\cdot 0.76=0.61$. This means that
$\mathrm{eff}_{c}(\bar{\xi})^{-1}-1=0.39/0.61=64\,\%$ more units have to be
used for design $\bar{\xi}$ to obtain the same quality for estimating the
median failure time than for the optimal design $\xi^{*}$. Hence, the optimal
design $\xi^{*}$ performs much better than the standard design $\bar{\xi}$ in
this situation.

Even more prominent results can be obtained for the marginal model without
interactions between the stress variables.

######  Example (Example 3 cont.).

In the model with two non-interacting stress variables $x_{1}$ and $x_{2}$ the
value of the $c$-criterion for the locally $c$-optimal design $\xi^{*}$ for
extrapolation at $\mathbf{x}_{u}=(x_{u1},x_{u2})$, $x_{u1}<x_{u2}<0$, is given
by $\Phi_{c}(\xi^{*})=1/\lambda_{c}^{2}=(1+2|x_{u1}|)^{2}$ as seen before. The
uniform design $\bar{\xi}$ which assigns equal weights $1/4$ to the four
vertices $(0,0)$, $(0,1)$, $(1,0)$, and $(1,1)$ of the experimental region has
a value of $\Phi_{c}(\bar{\xi})=1+(1+2|x_{u1}|)^{2}+(1+2|x_{u2}|)^{2}$. Hence,
the uniform design $\bar{\xi}$ has efficiency
$\mathrm{eff}_{c}(\bar{\xi})=((1+2|x_{u1}|)^{2}/(1+(1+2|x_{u1}|)^{2}+(1+2|x_{u2}|)^{2})$
which ranges from $1/3$ for $\mathbf{x}_{u}$ close to the combination
$(x_{1l},x_{2l})$ of lowest stress levels $x_{1l}=0$ and $x_{2l}=0$ to $1$
when the lower normal use condition $x_{u1}$ tends to minus infinity while
$x_{u2}$ remains fixed. Moreover, the efficiency approaches $1/2$ when
$x_{u2}\approx x_{u1}$ and both normal use conditions tend to minus infinity
simultaneously.

For example, when the normal use conditions are ${x}_{u1}=-0.5$ for the first
component and ${x}_{u2}=-0.4$ for the second component as specified in Table
2, then the values of the $c$-criterion are $\Phi_{c}(\xi^{*})=4.00$ for the
optimal design $\xi^{*}$ and $\Phi_{c}(\bar{\xi})=8.24$ for the uniform design
$\bar{\xi}$, respectively. Hence, the efficiency of the uniform design
$\bar{\xi}$ is $\mathrm{eff}_{c}(\bar{\xi})=4.00/8.24=0.49$. This means that
more than twice as many units have to be used for design $\bar{\xi}$ to obtain
the same quality for estimating the median failure time than for the optimal
design $\xi^{*}$. This highlights that the optimal design $\xi^{*}$ performs
substantially better than the standard design $\bar{\xi}$ in the current model
of two non-interacting stress variables.

The above efficiency calculations are all related to estimating the median
failure time for soft failure due to degradation under normal use conditions
$\mathbf{x}_{u}$. For estimating any other quantile $t_{\alpha}$ of the
failure time distribution, the efficiency of a design $\xi$ can be written as

$\mathrm{eff}_{\mathrm{aVar}}(\xi)=\mathrm{eff}_{c}(\xi)+(1-\mathrm{eff}_{c}(\xi))\mathbf{c}_{\boldsymbol{\varsigma}}^{T}\widetilde{\mathbf{M}}_{\boldsymbol{\varsigma}}^{-1}\mathbf{c}_{\boldsymbol{\varsigma}}/\mathrm{aVar}(\widehat{t}_{\alpha};\xi)$
(8.4)

by equations (7.7) and (6.4). This efficiency depends on the variance
parameters, but it is bounded from below by the $c$-efficiency
$\mathrm{eff}_{c}(\xi)$ of $\xi$ for extrapolation at $\mathbf{x}_{u}$. Hence,
designs with a high efficiency for estimating the median failure time are also
suitable for estimating any other reasonable quantile $t_{\alpha}$,
$0<t_{\alpha}<\infty$.

##  9 Optimization of the time plan

In Section 8, we considered the situation when there is a fixed time plan
$\mathbf{t}=(t_{1},...,t_{k})^{T}$, and only the stress variables
$\mathbf{x}\in\mathcal{X}$ are to be optimized across units. In contrast to
that we consider now the situation when also the settings $t_{1},...,t_{k}$
for the standardized time variable $t\in[0,1]$ may be optimized within units.
As in Section 3 we still assume that the same time plan
$\mathbf{t}=(t_{1},...,t_{k})^{T}$ is used for all units.

Similar to the case of fixed time points $t_{1},...,t_{k}$, the standardized
asymptotic variance is given by

$\displaystyle\mathrm{aVar}(\widehat{t}_{\alpha};\xi,t_{1},...,t_{k})$ (9.1)
$\displaystyle=$ $\displaystyle
c_{0}^{2}\left(\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}(\xi)^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})\cdot\mathbf{f}_{2}(t_{\alpha})^{T}\mathbf{M}_{2}(t_{1},...,t_{k})^{-1}\mathbf{f}_{2}(t_{\alpha})+\mathbf{c}_{\boldsymbol{\varsigma}}^{T}\widetilde{\mathbf{M}}_{\boldsymbol{\varsigma}}(t_{1},...,t_{k})^{-1}\mathbf{c}_{\boldsymbol{\varsigma}}\right)$

by equations (7.7) and (7.9). But here the dependence of the asymptotic
variance on the settings $t_{1},...,t_{k}$ for the time variable is explicitly
stated which comes through the information matrix
$\mathbf{M}_{2}(t_{1},...,t_{k})$ in the second marginal model as well as the
standardized information matrix
$\widetilde{\mathbf{M}}_{\boldsymbol{\varsigma}}(t_{1},...,t_{k})$ for the
variance parameters.

As has been seen in Section 8 the optimization with respect to the stress
variables $\mathbf{x}_{u}$ can be done independently of the second marginal
model and the settings $t_{1},...,t_{k}$ of the time variable. The reverse
does not hold true for the time points $t_{1},...,t_{k}$, in general, by the
extra dependence of the asymptotic variance on the information matrix
$\widetilde{\mathbf{M}}_{\boldsymbol{\varsigma}}(t_{1},...,t_{k})$ for the
variance parameters.

To circumvent this problem we restrict to the case of estimating the median
failure time $t_{0.5}$. There the second term on the right hand side of
equation (9.1) vanishes and the standardized asymptotic variance simplifies to

$\mathrm{aVar}(\widehat{t}_{0.5};\xi,t_{1},...,t_{k})=c_{0}^{2}\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}(\xi)^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})\cdot\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}(t_{1},...,t_{k})^{-1}\mathbf{f}_{2}(t_{0.5})\,.$
(9.2)

Hence, for the median failure time also the optimization of the measurement
times $t_{1},..,t_{k}$ can be performed independently of the marginal model of
the stress variables and their settings. Then, only the marginal $c$-criterion
$\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}(t_{1},...,t_{k})^{-1}\mathbf{f}_{2}(t_{0.5})$
has to be minimized.

Remember that the marginal information matrix
$\mathbf{M}_{2}(t_{1},...,t_{k})$ is given by
$\mathbf{M}_{2}(t_{1},...,t_{k})=\mathbf{F}_{2}(t_{1},...,t_{k})^{T}\mathbf{V}(t_{1},...,t_{k})^{-1}\mathbf{F}_{2}(t_{1},...,t_{k})$
in the second marginal model for the time variable $t$, where the variance
covariance matrix $\mathbf{V}(t_{1},...,t_{k})$ of the vector $\mathbf{Y}_{i}$
of measurements for each unit $i$ is given by
$\mathbf{V}(t_{1},...,t_{k})=\mathbf{F}_{2}(t_{1},...,t_{k}){\boldsymbol{\Sigma}}_{\gamma}\mathbf{F}_{2}(t_{1},...,t_{k})^{T}+{\boldsymbol{\Sigma}}_{\varepsilon}$
and ${\boldsymbol{\Sigma}}_{\gamma}$ and ${\boldsymbol{\Sigma}}_{\varepsilon}$
are the variance covariance matrices for the random effects
$\boldsymbol{\gamma}$ and for the measurement errors
$\boldsymbol{\varepsilon}$, respectively.

(Schmelter, 2007) proposed a representation of the inverse of the information
matrix in random effects models in the case of uncorrelated homoscedastic
errors
($\boldsymbol{\Sigma}_{\varepsilon}=\sigma^{2}_{\varepsilon}\mathbf{I}_{k}$).
This can be readily extended to a general, non-singular variance covariance
structure $\boldsymbol{\Sigma}_{\varepsilon}$ for the measurement errors
$\boldsymbol{\varepsilon}$ (see Lemma A.1 in the Appendix). From this we
obtain for the marginal information matrix $\mathbf{M}_{2}(t_{1},...,t_{k})$
that its inverse can be decomposed to
$\mathbf{M}_{2}(t_{1},...,t_{k})^{-1}=\mathbf{F}_{2}(t_{1},...,t_{k})^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F}_{2}(t_{1},...,t_{k})+\boldsymbol{\Sigma}_{\gamma}$.

As a consequence, the marginal $c$-criterion
$\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}(t_{1},...,t_{k})^{-1}\mathbf{f}_{2}(t_{0.5})$
can be split up into

$\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}(t_{1},...,t_{k})^{-1}\mathbf{f}_{2}(t_{0.5})=\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}^{(0)}(t_{1},...,t_{k})^{-1}\mathbf{f}_{2}(t_{0.5})+\mathbf{f}_{2}(t_{0.5})^{T}\boldsymbol{\Sigma}_{\gamma}\mathbf{f}_{2}(t_{0.5}),$

where
$\mathbf{M}_{2}^{(0)}(t_{1},...,t_{k})=\mathbf{F}_{2}(t_{1},...,t_{k})^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F}_{2}(t_{1},...,t_{k})$
is the information matrix in the marginal fixed effect model

$Y_{j}^{(2,0)}=\mathbf{f}_{2}(t_{j})^{T}\boldsymbol{\beta}^{(2)}+\varepsilon_{j}^{(2)},$
(9.3)

$j=1,...,k$, in the time variable $t$ with variance covariance matrix
$\boldsymbol{\Sigma}_{\varepsilon}$ for the vector
$\boldsymbol{\varepsilon}^{(2)}=(\varepsilon_{1}^{(2)},...,\varepsilon_{k}^{(2)})^{T}$
of error terms. Hence, as for any linear criterion, the optimization of the
asymptotic variance (9.2) of the median failure time ${t}_{0.5}$ with respect
to the time plan $t_{1},...,t_{k}$ does not depend on the variance covariance
matrix $\boldsymbol{\Sigma}_{\gamma}$ of the random effects. For optimization
only the term
$\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}^{(0)}(t_{1},...,t_{k})^{-1}\mathbf{f}_{2}(t_{0.5})$
has to be minimized which is the $c$-criterion for extrapolation of the mean
response at $t_{0.5}$ in the marginal fixed effect model in $t$ with
$\mathrm{Cov}(\boldsymbol{\varepsilon}^{(2)})=\boldsymbol{\Sigma}_{\varepsilon}$.
This leads to the following result which is similar to Proposition 8.1 for
optimization with respect to the stress variable $\mathbf{x}$.

######  Proposition 9.1.

If the time plan $t_{1}^{*},...,t_{k}^{*}$ is $c$-optimal for extrapolation of
the mean response at $t_{0.5}$ in the marginal fixed effect model (9.3) with
covariance $\boldsymbol{\Sigma}_{\varepsilon}$ for the time variable $t$, then
$t_{1}^{*},...,t_{k}^{*}$ minimize the asymptotic variance for the estimator
$\widehat{t}_{0.5}$ of the median failure time $t_{0.5}$ under normal use
condition.

Note that in degradation experiments the median failure time $t_{0.5}$ under
normal use condition is typically much larger than the time horizon of the
experiment, i. e. $t_{0.5}>1$ on the standardized time scale. However, the
above proposition also holds for interpolation, $t_{0.5}\in[0,1]$.

The optimal time plan depends on the location parameters $\boldsymbol{\beta}$
through $t_{0.5}$, but also on the variance covariance structure
$\boldsymbol{\Sigma}_{\varepsilon}$ of the measurement errors within units,
and are, hence, local. In Examples 1 and 2 we assumed uncorrelated
homoscedastic measurement errors,
$\boldsymbol{\Sigma}_{\varepsilon}=\sigma_{\varepsilon}^{2}\mathbf{I}_{k}$. In
that case, the $c$-criterion for extrapolation can be reduced to
$\mathbf{f}_{2}(t_{0.5})^{T}(\mathbf{F}_{2}(t_{1},...,t_{k})^{T}\mathbf{F}_{2}(t_{1},...,t_{k}))^{-1}\mathbf{f}_{2}(t_{0.5})$,
and the optimization does not depend on the error variance
$\sigma_{\varepsilon}^{2}$.

If the number $k$ of measurement times $t_{1},...,t_{k}$ is large, one might
be tempted to use the concept of approximate designs also here. In that case,
the approximate design $\tau$ will be defined by mutually distinct time points
$t_{1},...,t_{\ell}$ from the standardized experimental time interval
$\mathcal{T}=[0,1]$ with corresponding proportions $\pi_{1},...,\pi_{\ell}>0$
satisfying $\sum_{j=1}^{\ell}\pi_{j}=1$ with the interpretation that
(approximately) $\pi_{j}k$ measurements are performed at time point $t_{j}$,
$j=1,...,\ell$, for each unit. In the situation of a linear time trend,
$\mathbf{f}_{2}(t)=(1,t)^{T}$, as in Examples 1 and 2, this leads to
essentially the same extrapolation problem as for the stress variable in
Example 1. The $c$-optimal design for extrapolation at $t_{0.5}>1$ is
concentrated on the $\ell=2$ endpoints $t_{1}=0$ and $t_{2}=1$ of the
standardized experimental time interval $\mathcal{T}=[0,1]$ with corresponding
proportions $\pi_{1}=\frac{t_{0.5}-1}{2t_{0.5}-1}$ and
$\pi_{2}=\frac{t_{0.5}}{2t_{0.5}-1}$ (see (Schwabe, 1996b), Example 2.1).
Similarly, as for extrapolation of the stress variable in Example 1, the
proportion $\pi_{2}$ at the endpoint $t_{2}=1$ in the direction of the
extrapolation time $t_{0.5}$ decreases from $1$, when $t_{0.5}$ is close to
$1$, to $0.5$, when the distance gets large ($t_{0.5}\to\infty$).

However, from a practical point of view, it does not seem meaningful to have
replications, i. e. to have more than one measurement at a time at the same
unit. Moreover, even if this would be possible, these measurements would be
expected to be strongly correlated with a correlation beyond that caused by
the random effects. To avoid these technical and modeling problem additional
constraints have to be imposed on the time plan to guarantee independence of
the measurement errors. For instance, the different measurement times
$t_{1},...,t_{k}$ have to be at least some sufficiently large $\Delta t>0$
apart, or more specifically, they are restricted to a grid on the standardized
experimental time interval with grid size $\Delta t$. Nevertheless approximate
design theory can be used to determine optimal time plans. To do so the
standardized experimental time interval is discretized to a sufficiently
coarse grid, $\mathcal{T}=\\{j\Delta t;\,j=0,1,...,J\\}$, where $\Delta
t=1/J$, i. e. $\mathcal{T}=\\{0,\Delta t,2\Delta t,...,1\\}$. Additionally
constraints are imposed on the proportions $\pi_{j}$ that none of these
proportions is larger than $1/k$, i. e. the number of measurements at a time
is bounded by one and, hence, there are at least $k$ different time points.
Optimal approximate time plans can then be obtained under these constraints by
using standard algorithms for design optimization. Actually, the so obtained
proportions may be smaller than one for some of the time points. But a
theoretical result based on an equivalence theorem under constraints (convex
optimization with Kuhn-Tucker conditions) guarantees that this only occurs for
a small number of time points: (Sahm and Schwabe, 2001) proved for the
$D$-optimality criterion that the experimental region splits up in subregions
where the weight of the optimal design attains the upper bound $1/k$ and where
it is $0$. Only at the boundary of these subregions weights different from $0$
and $1/k$ may occur. In particular, for straight line regression on $[0,1]$,
the subregions with maximal weight $1/k$ are adjacent to the endpoints $0$ and
$1$ of the interval while in the interior the optimal design has zero weight.
This result carries over directly to other criteria like the $c$-criterion
under appropriate conditions which are met in the present case of
extrapolation. From this approach efficient exact designs can be obtained by
adjusting the weights to the admissible values $0$ and $1/k$ under the
constraint of total weight $1$.

######  Example (Example 2 cont.).

In the setup of Example 2 we deal with straight line regression

$Y_{j}^{(2,0)}=\beta_{1}^{(2)}+\beta_{2}^{(2)}t_{j}+\varepsilon_{j}^{(2)}$

for the fixed effects marginal model in the time variable $t$. The measurement
errors $\varepsilon_{j}^{(2)}$ are supposed to be uncorrelated and
homoscedastic
($\boldsymbol{\Sigma}_{\varepsilon}=\sigma^{2}_{\varepsilon}\mathbf{I}_{k}$).
We are searching for a locally $c$-optimal design $\tau^{*}$ for extrapolation
at the median failure time $t_{0.5}$ which is equal to $10.25$ under the
nominal values of Table 2. As constraints for the design we assume that $k=6$
observation can be taken on a grid with increment $\Delta t=0.05$ of the
standardized experimental time interval $[0,1]$, i. e. $J=20$ and
$\mathcal{T}=\\{0,0.05,0.10,...,1\\}$.

By the multiplicative algorithm (see e.g. (Silvey et al., 1978) ) adapted to
the present constraint situation the following numerical solution is obtained
for the locally $c$-optimal design

$\tau^{*}=\left(\begin{array}[]{ccccccc}0.00&0.05&0.10&0.85&0.90&0.95&1.00\\\
0.166&0.166&0.130&0.040&0.166&0.166&0.166\end{array}\right)\,.$

The optimal approximate design $\tau^{*}$ is supported on seven time points
$t^{*}_{1},...,t_{7}^{*}$ which are concentrated to the ends of the
standardized experimental time interval $[0,1]$ with maximal admitted weight
$1/k$ for all but the two boundary points $t_{3}^{*}=0.10$ and
$t_{4}^{*}=0.85$ separating grid points $t_{1}^{*}=0.00$, $t_{2}^{*}=0.05$,
$t_{5}^{*}=0.90$, $t_{6}^{*}=0.95$, and $t_{7}^{*}=1.00$ with full weight
$1/k$ from those with zero weight. This shape of the optimal design is in
accordance with the findings of (Sahm and Schwabe, 2001) in the case of
$D$-optimality. In view of Proposition 9.1 the design $\tau^{*}$ is also
optimal for the estimation of the median failure time.

For practical use, the optimal approximate design $\tau^{*}$ may be adjusted
to

$\tau_{0}=\left(\begin{array}[]{cccccc}0.00&0.05&0.10&0.90&0.95&1.00\\\
0.166&0.166&0.166&0.166&0.166&0.166\end{array}\right)$

which is supported on exactly $k=6$ time points by transferring the weight
from the boundary point $t_{4}^{*}=0.85$ with the lower weight to the boundary
point $t_{3}^{*}=0.10$ with the higher weight (see (Dorfleitner and Klein,
1999)). As a consequence all weights in the adjusted design $\tau_{0}$ are
equal to $1/6$. Then the design $\tau_{0}$ can be realized as an exact design
by taking one measurement at each of the six time points $0.00$, $0.05$,
$0.10$, $0.90$, $0.95$, and $1.00$. To quantify what might have got lost, the
quality of the adjusted design $\tau_{0}$ may be measured in terms of the
local $c$-efficiency

$\mathrm{eff}_{c}(\tau_{0})=\frac{\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}(\tau^{*})^{-1}\mathbf{f}_{2}(t_{0.5})}{\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}(\tau_{0})^{-1}\mathbf{f}_{2}(t_{0.5})}=\frac{464.42}{469.91}=98.70\%$

for extrapolation at $t_{0.5}$ in the marginal mixed effects model and, hence,
for estimation of the median failure time. This indicates that the adjusted
design $\tau_{0}$ is highly efficient and can be recommended as the time plan
for conducting an accelerated degradation testing experiment.

Note that by Proposition 9.1 the optimal design $\tau^{*}$ does not depend on
the variance covariance structure $\boldsymbol{\Sigma}_{\gamma}$ of the random
effects, but the efficiency of the adjusted design $\tau_{0}$ may be affected
by the random effects. Nevertheless, similar to equation (8.4), also here the
$c$-efficiency for extrapolation at $t_{0.5}$ provides a lower bound for the
efficiency in estimating the median failure time. The assumption
$\boldsymbol{\Sigma}_{\varepsilon}=\sigma_{\varepsilon}^{2}\mathbf{I}_{k}$ of
uncorrelated measurement errors could be replaced by equal correlations
$\rho_{\varepsilon}$ between all measurements, i. e.
$\boldsymbol{\Sigma}_{\varepsilon}$ is compound symmetric with all diagonal
entries equal to $\sigma_{\varepsilon}^{2}$ and all off-diagonal entries equal
to $\rho_{\varepsilon}\sigma_{\varepsilon}^{2}$. For estimating the variance
covariance parameters an additional identifiability condition would be
required in this case to distinguish the correlation $\rho_{\varepsilon}$ of
the measurement errors from the variance $\sigma_{1}^{2}$ of the random
intercept. Then estimation of the location parameters $\boldsymbol{\beta}$ is
not affected, and the design $\tau^{*}$ retains its optimality.

For further approaches allowing for correlations depending on the distance of
the measurement times, e. g. when measurements follow a random process, we
refer to (Näther, 1985) and (Müller, 2007) for optimization of the measurement
times.

##  10 Experimental design with a cross-sectional time plan

In contrast to the previous sections we will allow here for time plans
$t_{i1},...,t_{ik}$ differing across units,

$Y_{ij}=(\mathbf{f}_{1}(\mathbf{x}_{i})\otimes\mathbf{f}_{2}(t_{ij}))^{T}\boldsymbol{\beta}+\mathbf{f}_{2}(t_{ij})^{T}\boldsymbol{\gamma}_{i}+\varepsilon_{ij}\,,$
(10.1)

$j=1,...,k$, $i=1,...,n$, where all other expressions have the same meaning as
in the general model equation 3.5. Differing time plans are, for example,
required when the number $k$ of time points is less than the number $p_{2}$ of
parameters in the marginal model for the degradation paths. In particular, in
the case of destructive testing only $k=1$ measurement will be available per
unit,

$Y_{i}=(\mathbf{f}_{1}(\mathbf{x}_{i})\otimes\mathbf{f}_{2}(t_{i}))^{T}\boldsymbol{\beta}+\mathbf{f}_{2}(t_{i})^{T}\boldsymbol{\gamma}_{i}+\varepsilon_{i}\,,$
(10.2)

$i=1,...,n$. We will restrict to the case of $k=1$ measurements per unit in
the remainder of this section. In that case all measurements $Y_{i}$ are
assumed to be independent. Denote by

$\sigma^{2}(t)=\mathbf{f}_{2}(t)^{T}\boldsymbol{\Sigma}_{\gamma}\mathbf{f}_{2}(t)+\sigma^{2}_{\varepsilon}>0$

the variance function for measurements at time $t$, i. e.
$\mathrm{Var}(Y_{i})=\sigma^{2}(t_{i})$. As discussed at the end of the
previous section also here an identifiability condition has to be imposed on
the variance parameters $\boldsymbol{\varsigma}$ to distinguish between the
variance $\sigma^{2}_{\varepsilon}$ of the measurement error and the variance
$\sigma_{1}^{2}$ of the random intercept, see (Graßhoff et al., 2012) who
derived $D$-optimal designs.

Under the assumption of normally distributed random effects and measurement
errors the information matrix $\mathbf{M}_{\boldsymbol{\beta}}$ for the
location parameter $\boldsymbol{\beta}$ can be written as

$\mathbf{M}_{\boldsymbol{\beta}}=\sum_{i=1}^{n}(\mathbf{f}_{1}(\mathbf{x}_{i})\mathbf{f}_{1}(\mathbf{x}_{i})^{T})\otimes\left({\textstyle{\frac{1}{\sigma^{2}(t_{i})}}}\mathbf{f}_{2}(t_{i})\mathbf{f}_{2}(t_{i})^{T}\right)$

according to the product-type structure of the model (10.2). For estimating
the median failure time $t_{0.5}$ the asymptotic variance is given by
$\mathrm{aVar}(\widehat{t}_{0.5})=c_{0}^{2}\mathbf{c}^{T}\mathbf{M}_{\boldsymbol{\beta}}^{-1}\mathbf{c}$
as in Section 7. Hence, we are searching for a $c$-optimal design, where the
vector $\mathbf{c}$ factorizes according to
$\mathbf{c}=\mathbf{c}_{1}\otimes\mathbf{c}_{2}$ into components
$\mathbf{c}_{1}=\mathbf{f}_{1}(\mathbf{x}_{u})$ and
$\mathbf{c}_{2}=\mathbf{f}_{2}(t_{0.5})$ associated with the marginal models
for the stress and the time variables, see (7.9). Moreover, we suppose, as
before, that the experimental settings $\mathbf{x}_{i}$ for the stress
variable and $t_{i}$ for the time of measurement may be chosen independently
from their standardized experimental regions $\mathcal{X}$ and $\mathcal{T}$,
respectively, i. e. $(\mathbf{x}_{i},t_{i})\in\mathcal{X}\times\mathcal{T}$.

Also here we make use of the concept of approximate designs denoted by $\zeta$
which are specified by $m$ mutually distinct combinations
$(\mathbf{x}_{i},t_{i})\in\mathcal{X}\times\mathcal{T}$ with corresponding
weights $\eta_{i}>0$, $i=1,...,\nu$, $\sum_{i=1}^{\nu}\eta_{i}=1$. Accordingly
the corresponding normalized information matrix is defined as

$\mathbf{M}(\zeta)=\sum_{i=1}^{\nu}\eta_{i}(\mathbf{f}_{1}(\mathbf{x}_{i})\mathbf{f}_{1}(\mathbf{x}_{i})^{T})\otimes(\tilde{\mathbf{f}}_{2}(t_{i})\tilde{\mathbf{f}}_{2}(t_{i})^{T})\,,$

where $\tilde{\mathbf{f}}_{2}(t)=\mathbf{f}_{2}(t)/\sigma(t)$ is the weighted
version of the marginal regression function for the time variable
(standardized by the standard deviation for measurement at time $t$). Then a
$c$-optimal design for $\mathbf{c}=\mathbf{c}_{1}\otimes\mathbf{c}_{2}$ can be
obtained as the product-type design generated from $c$-optimal designs for
$\mathbf{c}_{v}$ in the associated marginal models.

To be more specific, let

$\xi=\left(\begin{array}[]{ccc}\mathbf{x}_{1}&...&\mathbf{x}_{m}\\\
w_{1}&...&w_{m}\end{array}\right)\qquad\mathrm{and}\qquad\tau=\left(\begin{array}[]{ccc}t_{1}&...&t_{\ell}\\\
\pi_{1}&...&\pi_{\ell}\end{array}\right)$

be approximate designs in the marginal models for the stress variables and the
time, respectively. Then the product-type design $\zeta=\xi\otimes\tau$ is
supported on the $\nu=m\ell$ mutually distinct combinations
$(\mathbf{x}_{i},t_{j})$ with corresponding weights
$\eta_{ij}=w_{i}\pi_{\ell}$, $i=1,...,m$, $j=1,...,\ell$, and its standardized
information matrix

$\mathbf{M}(\xi\otimes\tau)=\mathbf{M}_{1}(\xi)\otimes\tilde{\mathbf{M}}_{2}(\tau)$

factorizes into the standardized information matrix
$\mathbf{M}_{1}(\xi)=\sum_{i=1}^{m}w_{i}\mathbf{f}_{1}(\mathbf{x}_{i})\mathbf{f}_{1}(\mathbf{x}_{i})^{T}$
of $\xi$ in the marginal model for the stress variables and the standardized
information matrix
$\tilde{\mathbf{M}}_{2}(\tau)=\sum_{j=1}^{k}\pi_{j}\tilde{\mathbf{f}}_{2}(t_{j})\tilde{\mathbf{f}}_{2}(t_{j})^{T}$
of $\tau$ in the weighted marginal model

$Y_{j}^{(2)}=\tilde{\mathbf{f}}_{2}(t_{j})^{T}\boldsymbol{\beta}^{(2)}+\varepsilon^{(2)}$
(10.3)

for the time variable. Then also the $c$-criterion

$\mathbf{c}^{T}\mathbf{M}(\xi\otimes\tau)^{-1}\mathbf{c}=\mathbf{c}_{1}^{T}\mathbf{M}_{1}(\xi)^{-1}\mathbf{c}_{1}\cdot\mathbf{c}_{2}^{T}\tilde{\mathbf{M}}_{2}(\tau)^{-1}\mathbf{c}_{2}$

for $\mathbf{c}=\mathbf{c}_{1}\otimes\mathbf{c}_{2}$ factorizes into the
$c$-criteria for $\mathbf{c}_{1}$ and $\mathbf{c}_{2}$ in the corresponding
marginal models. This expression is minimized if $\zeta^{*}$ is the product
$\zeta^{*}=\xi^{*}\otimes\tau^{*}$ of the $c$-optimal marginal designs
$\xi^{*}$ and $\tau^{*}$ for $\mathbf{c}_{1}$ and $\mathbf{c}_{2}$ in their
marginal models, respectively. The design $\zeta^{*}=\xi^{*}\otimes\tau^{*}$
is not only $c$-optimal for $\mathbf{c}=\mathbf{c}_{1}\otimes\mathbf{c}_{2}$
in the class of product-type designs $\xi\otimes\tau$, but also compared to
all designs $\zeta$ on $\mathcal{Z}$ (see (Schwabe, 1996a), Theorem 4.4) which
establishes the following result.

######  Theorem 10.1.

If $\xi^{*}$ is $c$-optimal for extrapolation at normal use condition
$\mathbf{x}_{u}$ in the marginal model for the stress variables $\mathbf{x}$
and $\tau^{*}$ is $c$-optimal for extrapolation at the median failure time
$t_{0.5}$ in the weighted marginal model (10.3) for the time variable $t$,
then the design $\zeta^{*}=\xi^{*}\otimes\tau^{*}$ is optimal for estimating
the median failure time $t_{0.5}$ in the combined model (10.2).

Note that for the second marginal model we utilize that
$\tilde{\mathbf{f}}_{2}(t_{0.5})$ and $\mathbf{c}_{2}=\mathbf{f}_{2}(t_{0.5})$
only differ by the factor $1/\sigma(t_{0.5})>0$ which does not affect
$c$-optimality.

The $c$-optimal designs $\xi^{*}$ for extrapolation at $\mathbf{x}_{u}$ are
the same as obtained in Section 8, and the $c$-optimal designs for
extrapolation at $t_{0.5}$ can be obtained similarly by Elfving’s theorem
applied to the weighted regression functions
$\tilde{\mathbf{f}}_{2}(t)=\mathbf{f}_{2}(t)/\sigma(t)$.

######  Example (Example 1 cont.).

In the introductory model of Section 2 we have straight line regression in
both the stress variable $x$ and the time variable $t$ with an interaction
term $xt$. For the stress variable, the $c$-optimal marginal design $\xi^{*}$
for extrapolation at $x_{u}<0$ assigns weight $w^{*}=|x_{u}|/(1+2|x_{u}|)$ to
$x_{h}=1$ and weight $1-w^{*}=(1+|x_{u}|)/(1+2|x_{u}|)$ to $x_{l}=0$, see
Section 8. For the time variable $t$, the $c$-optimal marginal design
$\tau^{*}$ for extrapolation at $t_{0.5}>1$ can be similarly obtained by
Elfving’s theorem applied to the vector of weighted regression functions
$\tilde{\mathbf{f}}_{2}(t)$. The shape of the Elfving set is exhibited in
Figure 4.

Figure 4: Elfving set for $t$ in Example 1: Induced design region (right solid
line), negative image (left solid line), boundary (dashed lines),
$\mathbf{c}_{2}=(1,t_{0.5})^{T}$ (arrow) and corresponding ray (dotted line)

Note that the induced design region
$\\{\tilde{\mathbf{f}}_{2}(t);\,t\in[0,1]\\}$ and its negative image
constitute non-overlapping arc segments of the ellipse defined by
$\mathbf{z}^{T}\boldsymbol{\Sigma}\mathbf{z}=1$,
$\mathbf{z}\in\mathbb{R}^{2}$, centered at $(0,0)^{T}$, where
$\boldsymbol{\Sigma}=\boldsymbol{\Sigma}_{\gamma}+\left(\begin{array}[]{cc}\sigma_{\varepsilon}^{2}&0\\\
0&0\end{array}\right)$ is the variance covariance matrix of the random effects
augmented by the variance of the measurement error. The ray
$\lambda\mathbf{f}_{2}(t_{0.5})$ intersects the boundary of the Elfving set at
the line segment connecting $\tilde{\mathbf{f}}_{2}(1)=(1,1)^{T}/\sigma(1)$
and $-\tilde{\mathbf{f}}_{2}(0)=(1,0)^{T}/\sigma(0)$ when $t_{0.5}>1$. Hence,
the $c$-optimal design $\tau^{*}$ is supported by the endpoints of the
standardized experimental time interval $\mathcal{T}=[0,1]$, and the optimal
weights can be calculated by Elfving’s theorem to
$\pi^{*}=t_{0.5}\sigma(1)/(t_{0.5}\sigma(1)+(t_{0.5}-1)\sigma(0))$ at
$t_{\max}=1$ and
$1-\pi^{*}=(t_{0.5}-1)\sigma(0)/(t_{0.5}\sigma(1)+(t_{0.5}-1)\sigma(0))$ at
$t_{\min}=0$.

By Theorem 10.1, the design

$\zeta^{*}=\xi^{*}\otimes\tau^{*}=\left(\begin{array}[]{cccc}(x_{l},t_{\min})&(x_{l},t_{\max})&(x_{h},t_{\min})&(x_{h},t_{\max})\\\
(1-w^{*})(1-\pi^{*})&(1-w^{*})\pi^{*}&w^{*}(1-\pi^{*})&w^{*}\pi^{*}\end{array}\right)$

is optimal for estimating the median failure time $t_{0.5}$. Under the nominal
values of Table 1 the optimal weights are $w^{*}=0.05$ for extrapolation at
$x_{u}=-0.056$ in the marginal model for the stress variable, see Section 8,
and $\pi^{*}=0.77$ for extrapolation at $t_{0.5}=1.583$ in the weighted
marginal model for the time variable. Hence, the optimal design for estimating
the median failure time is given by

$\zeta^{*}=\left(\begin{array}[]{cccc}(0,0)&(0,1)&(1,0)&(1,1)\\\
0.22&0.73&0.01&0.04\end{array}\right)\,.$

According to this design $73\,\%$ of the units should be exposed to the lowest
stress level $x_{l}=0$ and measured at the maximum allowed time $t_{\max}=1$
while only $1\,\%$ of the units should be exposed to the highest stress level
$x_{h}=1$ and measured at the initial time $t_{\min}=0$. Accordingly $22\,\%$
should be exposed to low stress and measured initially, and $4\,\%$ should be
measured at the end of the experimental time interval under high exposure.

A similar results can be obtained for the situation of Example 2.

######  Example (Example 2 cont.).

In the model with two interacting stress variables (Example 2) the $c$-optimal
marginal design $\xi^{*}$ for extrapolation at
$\mathbf{x}_{u}=(x_{u1},x_{u2})^{T}$ is the same as in Section 8,
$\xi^{*}=\xi_{1}^{*}\otimes\xi_{2}^{*}$. For the time variable $t$, the
$c$-optimal marginal design $\tau^{*}$ for extrapolation at $t_{0.5}>1$ is
obtained as in Example 1: $\tau^{*}$ assigns weights
$\pi^{*}=t_{0.5}\sigma(1)/(t_{0.5}\sigma(1)+(t_{0.5}-1)\sigma(0))$ to
$t_{\max}=1$ and $1-\pi^{*}$ to $t_{\min}=0$. By Theorem 10.1, the optimal
design for estimation of the median failure time is then
$\zeta^{*}=\xi_{1}^{*}\otimes\xi_{2}^{*}\otimes\tau^{*}$.

Under the nominal values of Table 2 the optimal weights are $w_{1}^{*}=0.25$
and $w_{2}^{*}=0.22$ for extrapolation at ${x}_{u1}=-0.5$ and ${x}_{u2}=-0.4$
in the submarginal models of the stress variables, respectively, see Section
8. For the time variable the optimal weight results in $\pi^{*}=0.57$ for
extrapolation at $t_{0.5}=10.25$ in the weighted marginal model. Hence, the
optimal design for estimating the median failure time is given by

$\xi_{1}^{*}\otimes\xi_{2}^{*}\otimes\tau^{*}=\left(\begin{array}[]{cccccccc}(0,0,0)&(0,0,1)&(0,1,0)&(0,1,1)&(1,0,0)&(1,0,1)&(1,1,0)&(1,1,1)\\\
0.25&0.33&0.07&0.10&0.09&0.11&0.02&0.03\end{array}\right).$

The locally optimal designs for estimating the median failure time are
influenced by both the location and the variance parameters
$\boldsymbol{\beta}$ and $\boldsymbol{\varsigma}$. Therefore we make a
sensitivity analysis, how the optimal designs change with the parameters and
how well they perform under parameter misspecification. For this we note first
that the optimal marginal design $\xi^{*}$ for extrapolation at
$\mathbf{x}_{u}$ is globally optimal and does not depend on the parameters.
Moreover, in the case of straight line degradation paths as in Examples 1 and
2 only the weight
$\pi^{*}=t_{0.5}\sigma(1)/(t_{0.5}\sigma(1)+(t_{0.5}-1)\sigma(0))$ of the
optimal marginal design $\tau^{*}$ depends on the location parameters
$\boldsymbol{\beta}$ through $t_{0.5}$ and on the variance parameters
$\boldsymbol{\varsigma}$ through the ratio $\sigma(1)/\sigma(0)$ of the
standard deviations at the endpoints of the experimental time region. Note
that the $c$-efficiency for extrapolation at $t_{0.5}$ in the weighted
marginal model for $t$ coincides with the efficiency for estimating the median
failure time when the optimal design $\xi^{*}$ is used for the stress
variables.

######  Example (Example 1 cont.).

For the introductory model with one stress variable (Example 1), the optimal
weight $\pi^{*}=t_{0.5}\sigma(1)/(t_{0.5}\sigma(1)+(t_{0.5}-1)\sigma(0))$ is
plotted in Figure 5 as a function of $t_{0.5}$ while the ratio
$\sigma(1)/\sigma(0)$ is held fixed to $1.22$ induced by the nominal values in
Table 1 (left panel) and as a function of the ratio $\sigma(1)/\sigma(0)$
while $t_{0.5}$ is held fixed to the nominal value $1.583$ (right panel).

Figure 5: Optimal weights $\pi^{*}$ in dependence on $t_{0.5}$ (left panel)
and $\sigma(1)/\sigma(0)$ (right panel) for Example 1

When $t_{0.5}$ increases, the optimal weight $\pi^{*}$ decreases from $1$ for
$t_{0.5}$ close to the maximal experimental time $t_{\max}=1$ to
$\sigma(1)/(\sigma(1)+\sigma(0))=0.55$ for $t_{0.5}\to\infty$. On the other
hand the optimal weight $\pi^{*}$ increases in the ratio $\sigma(1)/\sigma(0)$
of standard deviations from $0$ for the ratio close to $0$ to $1$ when the
ratio tends to infinity. The limits for $\pi^{*}$ are indicated in Figure 5 by
horizontal dashed lines in both panels while the nominal values for $t_{0.5}$
and $\sigma(1)/\sigma(0)$ are indicated by vertical dotted lines,
respectively.

To assess the influence of the variation of the optimal weights we consider
the efficiency of a locally optimal design when the underlying parameters are
misspecified. Define by

$\mathrm{eff}_{\mathrm{aVar}}(\zeta;\boldsymbol{\theta})=\frac{\mathrm{aVar}_{\boldsymbol{\theta}}(\widehat{t}_{0.5};\zeta)}{\mathrm{aVar}_{\boldsymbol{\theta}}(\widehat{t}_{0.5};\zeta^{*}_{\boldsymbol{\theta}})}$

the asymptotic efficiency of the design $\zeta$ for estimating $t_{0.5}$,
where $\mathrm{aVar}_{\boldsymbol{\theta}}(\widehat{t}_{0.5};\zeta)$ denotes
the asymptotic variance of $\widehat{t}_{0.5}$, when $\zeta$ is used and when
$\boldsymbol{\theta}$ is the true vector of underlying parameters, and by
$\zeta^{*}_{\boldsymbol{\theta}}$ the locally optimal design at
$\boldsymbol{\theta}$. Similar to the situation of submarginal models in
Example 2 (see Section 8), the asymptotic efficiency

$\mathrm{eff}_{\mathrm{aVar}}(\xi\otimes\tau;\boldsymbol{\theta})=\mathrm{eff}_{c1}(\xi)\cdot\mathrm{eff}_{c2}(\tau;\boldsymbol{\theta})\,,$

factorizes for product-type designs $\zeta=\xi\otimes\tau$ into the marginal
efficiency $\mathrm{eff}_{c1}(\xi)$ of $\xi$ for extrapolation at
$\mathbf{x}_{u}$ in the marginal model for the stress variable and the
marginal efficiency $\mathrm{eff}_{c2}(\tau;\boldsymbol{\theta})$ of $\tau$
for extrapolation at $t_{0.5}$ in the weighted marginal model for the time
variable. Remember that, for the stress variable, the marginal information
matrix $\mathbf{M}_{1}(\xi)$ and, hence, the marginal efficiency
$\mathrm{eff}_{c1}(\xi)$ does not depend on $\boldsymbol{\theta}$. Thus, also
the $c$-optimal marginal design $\xi^{*}$ does not rely on the parameters and
can be used throughout in the comparison. The asymptotic efficiency of a
product-type design $\zeta=\xi^{*}\otimes\tau$ with optimal marginal $\xi^{*}$
reduces to the $c$-efficiency of the second marginal $\tau$,

$\mathrm{eff}_{\mathrm{aVar}}(\xi^{*}\otimes\tau;\boldsymbol{\theta})=\mathrm{eff}_{c2}(\tau;\boldsymbol{\theta})=\frac{\mathbf{c}_{2}(\boldsymbol{\theta})^{T}\mathbf{M}_{2}(\tau;\boldsymbol{\theta})\mathbf{c}_{2}(\boldsymbol{\theta})}{\mathbf{c}_{2}(\boldsymbol{\theta})^{T}\mathbf{M}_{2}(\tau^{*}_{\boldsymbol{\theta}};\boldsymbol{\theta})\mathbf{c}_{2}(\boldsymbol{\theta})}\,,$

where $\mathbf{c}_{2}(\boldsymbol{\theta})=\mathbf{f}_{2}(t_{0.5})$ and
$\tau^{*}_{\boldsymbol{\theta}}$ is the locally $c$-optimal design at
$\boldsymbol{\theta}$ for extrapolation at $t_{0.5}$ in the weighted marginal
model for the time variable. The $c$-vector
$\mathbf{c}_{2}(\boldsymbol{\theta})=\mathbf{f}_{2}(t_{0.5})$ depends on the
location parameters $\boldsymbol{\beta}$ merely through $t_{0.5}$, but not on
$\boldsymbol{\varsigma}$. In contrast to that, the information matrices
$\mathbf{M}_{2}(\tau;\boldsymbol{\theta})$ are only affected by the variance
parameters $\boldsymbol{\varsigma}$, but not by $\boldsymbol{\beta}$. For
straight line regression in the time variable (cf. Examples 1 and 2), only the
ratio $\sigma(1)/\sigma(0)$ of the standard deviations for measurements at the
endpoints of the experimental time interval has an effect on the information
matrix $\mathbf{M}_{2}(\tau;\boldsymbol{\theta})$, when the design $\tau$ is
supported by these endpoints, as is the case for the locally $c$-optimal
design $\tau^{*}_{\boldsymbol{\theta}}$.

######  Example (Example 1 cont.).

For the setting of the introductory model with one stress variable (Example
1), we examine the efficiency of the design $\zeta*=\xi^{*}\otimes\tau^{*}$
which is locally optimal for estimation of the median failure time under the
nominal values of Table 1 when the nominal parameters are misspecified. In
Figure 6 we plot the efficiency of the locally optimal design
$xi^{*}\otimes\tau^{*}$ at the nominal values (solid line) together with the
efficiency of the designs $xi^{*}\otimes\bar{\tau}_{2}$ (dashed line) and
$xi^{*}\otimes\bar{\tau}_{6}$ (dashed and dotted line) for which the marginal
design $\bar{\tau}_{k}$ is uniform on $k$ equally spaced time points, i. e.
$\bar{\tau}_{2}$ assigns weight $1/2$ to the endpoints $0$ and $1$ of the time
interval, and $\tau_{6}$ assigns weight $1/6$ to the time points $0.0$, $0.2$,
$0.4$, $0.6$, $0.8$, and $1.0$ (cf. the definition of $\bar{\xi}_{m}$ in
Section 8).

Figure 6: Efficiency of $xi^{*}\otimes\tau^{*}$ (solid line),
$\xi^{*}\otimes\bar{\tau}_{2}$ (dashed line) and
$\xi^{*}\otimes\bar{\tau}_{6}$ (dashed and dotted line) in Example 1

In the left panel of Figure 6 the efficiency is displayed in dependence on the
true value of the median failure time $t_{0.5}$ while the ratio
$\sigma(1)/\sigma(0)$ of standard deviations is held fixed to the nominal
value $1.22$. In the right panel the efficiency is shown in dependence on the
true value of the ratio $\sigma(1)/\sigma(0)$ of standard deviations while the
median failure time $t_{0.5}$ is held fixed to $1.583$ derived from the
nominal values. Also here the nominal values for $t_{0.5}$ and
$\sigma(1)/\sigma(0)$ are indicated by vertical dotted lines in the
corresponding panel, respectively.

The efficiency seems to be more sensitive with respect to deviations in the
ratio than in the median failure time $t_{0.5}$. However, commonly neither
small values of $t_{0.5}$ nor small values of $\sigma(1)/\sigma(0)$ seem to be
reasonable. In particular, we may expect that the variance $\sigma(1)^{2}$ at
the end of the experimental time interval is larger than the variance
$\sigma(0)^{2}$ at the initial time, i. e. $\sigma(1)/\sigma(0)\geq 1$. This
is satisfied when $\rho\geq-\sigma_{1}/(2\sigma_{2})$. In total, the design
$\xi^{*}\otimes\tau^{*}$ which is locally optimal for the nominal values seems
to perform quite well over a wide range of parameters while the design
$\xi^{*}\otimes\bar{\tau}_{2}$ with the same number of units at the endpoints
of the intervals is only preferable for larger values of the median failure
time $t_{0.5}$. The design $\xi^{*}\otimes\bar{\tau}_{6}$ with six equally
spaced time points performs substantially worse throughout for reasonable
parameter values. It has to be noted in this context that the efficiency of
$\bar{\tau}_{6}$ depends on the variance parameters $\boldsymbol{\varsigma}$
not only through the ratio $\sigma(1)/\sigma(0)$ because measurements are also
to be taken in the interior of the interval. The current plot has been
generated by fixing $\sigma_{2}^{2}=\sigma_{1}^{2}+\sigma_{\varepsilon}^{2}$
and letting $\rho$ vary. However, other choices of the variance parameters do
not change much in the performance of the design.

Similar results can be obtained for the settings of Example 2.

##  11 Discussion and conclusion

During the design stage of highly reliable systems it is extremely important
to assess the reliability related properties of the product. One method to
handle this issue is to conduct accelerated degradation testing. Accelerated
degradation tests have the advantage to provide an estimation of lifetime and
reliability of the system under study in a relatively short period of time. To
account for variability between units in accelerated degradation tests, it is
assumed that the degradation function can be described by a mixed-effects
linear model. This also leads to a non-degenerate distribution of the failure
time, due to soft failure by exceedance of the expected (conditionally per
unit) degradation path over a threshold, under normal use conditions.