# Optimal Time Plan in Accelerated Degradation Testing

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

###### Abstract

Many highly reliable products are designed to function for years without
failure. For such systems accelerated degradation testing may provide
significance information about the reliability properties of the system. In
this paper, we propose the $c$-optimality criterion for obtaining optimal
designs of constant-stress accelerated degradation tests where the degradation
path follow a linear mixed effects model. The present work is mainly concerned
with developing optimal desings in terms of the time variable rather than
considering an optimal design of the stress variable which is usually
considered in the majority of the literature. Finally, numerical examples are
presented and sensitivity analysis procedures are conducted for evaluating the
robustness of the optimal as well as standard designs against
misspecifications of the nominal values.

###### keywords:

Accelerated degradation test, linear mixed effects model, $c$-optimal design,
.

##  1 Introduction

The growing research focus in recent years on the analysis of degradation data
has played a major role in the field of reliability inference and management.
In particular, reliability data and characteristics for highly reliable
products can be collected through analyzing the obtained degradation data.
Consequently, statistical modeling and inference approaches have been
introduced on the basis of different degradation measuring techniques.

[8], for instance, provided a comprehensive reviews on several approaches in
the literature that are utilized to evaluate reliable information from
acclerated degradation data. [19] introduced an analytical optimal ADT design
prcedure for more efficient reliability demonstration by minimizing the
asymptotic variance of decision variable in reliability demonstration under
the constraints of sample size, test duration, test cost, and predetermined
decision risks. Considering linear mixed effects model (LMEM), [18] considered
the minimum asymptotic variance criterion to characterize optimal design as
well as compromise design plans for accelerated degradation tests. Further,
[17] developed an algorithm-based optimal ADT approach by minimizing the
asymptotic variance of the MLE of the mean failure time of a system, where the
sample size and termination time of each run of the ADT at a constant
measurement frequency were specified. [14] charachterized optimal designs for
repeated measures accelerated degradation tests with multiple stress
variables, where the degradation paths are assumed to follow a linear mixed
effects. [5] introduced a modeling precedure to simultaneously analyze linear
degradation data and traumatic failures with competing risks in step stress
ADT experiment.

The general theory of optimal design of experiments is efficiently
charachterized in terms of the mathematical context of approximate designs
(see e. g. [15]). For example, [13] suggested optimal designs for multi-factor
models which can be used to treat more than one stress variable and the choice
of time plans simultaneously under various interaction structures.

[1] considered the problem of constructing $D$-optimal designs for linear and
nonlinear random effect models with applications in population
pharmacokinetics. [4] presented $D$-optimal designs for random coefficient
regression models in the basis of geometrical arguments when only one
observation is available per unit, i.e. a situation which occurs in
destructive testing.

It can be noted that the vast majority of the existing literature focuses on
desinging the degradation stress variable(s). As opposed to this approach the
time variable will be considered in the present work in order to obtain
optimal experimental designs under various measurement plans.

The present paper is organized as follows. In Sections 2 we state the general
model formulation and in Section 3 we exhibit the corresponding information
matrix. Section 4 is devoted to derive the design optimality criterion for
estimating a quantile of the failure time distribution under normal use
conditions. In Section 5 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 6 optimal measurement times are determined in the
setting of destructive testing. The paper closes with a short discussion in
Section 7.

##  2 Formulation of the model

In this section, we introduce a mixed effects regression model incorporating a
product-type structure as in [14]. In accordance with [14] we assume that
there are $n$ testing units $i=1,...,n$, for which degradation $y_{ij}$ will
be measured according to a time plan $\mathbf{t}=(t_{1},...,t_{k})^{T}$ with
$k$ time points $t_{j}$, $j=1,...,k$, $0\leq t_{1}<...<t_{k}\leq t_{\max}$.
Each unit $i$ is measutred under a setting $\mathbf{x}_{i}$ of the stress
variable(s), which is kept fixed for each unit throughout the testing, but
varies between units. The number $k$ of measurements and the measurent times
$t_{1},...,t_{k}$ are the same for all units. The measurements $y_{ij}$ are
considered as realizations of random variables $Y_{ij}$ which can be
charachterized by a hierarchical model. For each unit $i$ the observation
$Y_{ij}$ at $t_{j}$ is expressed as

$Y_{ij}=\mu_{i}(\mathbf{x}_{i},t_{j})+\varepsilon_{ij},$ (2.1)

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

$\mu_{i}(\mathbf{x},t)=\sum_{q=1}^{p}\beta_{iq}f_{q}(\mathbf{x},t_{j})=\mathbf{f}(\mathbf{x},t_{j})^{T}\boldsymbol{\beta}_{i}$
(2.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$, and $\boldsymbol{\beta}_{i}=(\beta_{i,1},...,\beta_{i,p})^{T}$ is a
$p$-dimensional vector of unit specific parameters $\beta_{i,q}$. Thus, the
response is given by

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

Further, 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}$. The vector
$\boldsymbol{\varepsilon}_{i}=(\varepsilon_{i1},...,\varepsilon_{ik})^{T}$ of
errors within one unit $i$ is $k$-dimensional multivariate normal with zero
mean and positive definite variance covariance matrix
$\boldsymbol{\Sigma}_{\varepsilon}$. 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, such that 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. Then the
observation $Y_{ij}$ can be written as

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

where the entries of the vector
$\boldsymbol{\beta}_{i}=(\beta_{i11},...,\beta_{i1p_{2}},...,\beta_{ip_{1}p_{2}})$
of parameters are relabeled lexicographically according to the marginal
regression functions. Moreover, 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, and that only the leading $p_{2}$ parameters
$\beta_{i11},...,\beta_{i1p_{2}}$ of $\boldsymbol{\beta}_{i}$ associated with
this constant term are unit specific. The remaining parameters in
$\boldsymbol{\beta}_{i}$ are assumed to have the same value
$\beta_{irs}=\beta_{rs}$, $r=2,...,p_{1}$, $s=1,...,p_{2}$, for all units
$i=1,...,n$. Hence, for unit $i$ the model (2.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},$
(2.5)

where $\boldsymbol{\beta}=(\beta_{11},...,\beta_{p_{1}p_{2}})^{T}$ is the
vector of (aggregate) fixed effect parameters (averaged over the units) 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 where $\boldsymbol{\gamma}_{i}$ are independent and
identically distributed multivariate normal with mean 0 and covariance matrix
$\boldsymbol{\Sigma}_{\gamma}$. For the analysis of degradation under normal
use, we further assume that the general model 2.5 is also valid at the normal
use condition $\mathbf{x}_{u}$, i. e.

$\mu_{u}(\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}$
(2.6)

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 multivariate normal with mean zero and variance
covariance matrix $\boldsymbol{\Sigma}_{\gamma}$.

##  3 Information and design

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}$. Let
$\mathbf{F}_{1}=(\mathbf{f}_{1}(\mathbf{x}_{1}),...,\mathbf{f}_{1}(\mathbf{x}_{n}))^{T}$
be the $n\times p_{1}$ marginal design matrix for the stress variables across
units, and denote by
$\mathbf{F}_{2}=(\mathbf{f}_{2}(t_{1}),...,\mathbf{f}_{2}(t_{k}))^{T}$ the
$k\times p_{2}$ marginal design matrix for the time variable. The
corresponding per unit variance covariance matrix is
$\mathbf{V}=\mathbf{F}_{2}{\boldsymbol{\Sigma}}_{\gamma}\mathbf{F}_{2}^{T}+{\boldsymbol{\Sigma}}_{\varepsilon}$
where $\mathbf{V}$ depends only on $\boldsymbol{\varsigma}$, see [14] for
further details about estimation in the model parameters. 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)$ (3.1)

where $\mathbf{M}_{\boldsymbol{\beta}}$ and
$\mathbf{M}_{\boldsymbol{\varsigma}}$ are 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. By (2.6), the information matrix $\mathbf{M}_{\boldsymbol{\beta}}$
for the parameters $\boldsymbol{\beta}$ factorizes according to

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

into the information matrix $\mathbf{M}_{1}=\mathbf{F}_{1}^{T}\mathbf{F}_{1}$
in the marginal model in the stress variable and the information matrix
$\mathbf{M}_{2}=\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2}$ in the
marginal mixed effects model in the time variable $t$ where the time plan
$\boldsymbol{t}=(t_{1},...,t_{k})^{T}$ plays the role of a (exact) design
here.

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

In accordance with [18] we are focusing on some characteristics of the failure
time distribution of soft failure due to degradation. As already mentioned, we
assume that the model equation (2.6)
$\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}$. A soft failure due to degradation is expressed as the
exceedance of the mean degradation $\mu_{u}$ over a failure threshold $y_{0}$.
The failure time $T$ under normal use condition is 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. As in [14] we 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}=(\delta_{1},...,\delta_{p_{2}})^{T}$ the vector of its
coefficients $\delta_{s}=\sum_{r=1}^{p_{1}}f_{1r}(\mathbf{x}_{u})\beta_{rs}$,
in the regression functions $f_{2s}$ in $t$, i. e.
$\mu(t)=\mathbf{f}_{2}(t)^{T}\boldsymbol{\delta}=\sum_{s=1}^{p_{2}}\delta_{s}f_{2s}(t)$.

For the purpose of describing 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})$ (4.1) $\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)},$ (4.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 standard normal distribution function.

We will be interested in quantiles $t_{\alpha}$ of the failure time
distribution, $\mathrm{P}(T\leq t_{\alpha})=\alpha$, up to which under normal
use conditions $\alpha\cdot 100$ percent of the units fail and
$(1-\alpha)\cdot 100$ percent of the units persist. Of particular importance
is the median $t_{0.5}$. By (4.1) these quantiles can be determined as the
solutions of the equation $h(t_{\alpha})=z_{\alpha}$, where
$z_{\alpha}=\Phi^{-1}(\alpha)$ is the $\alpha$-quantile of the standard normal
distribution. For the median ($\alpha=0.5$) we have $z_{0.5}=0$ and, hence the
median failure time $t_{0.5}$ is the solution of $\mu(t)=y_{0}$. In the
particular case of straight lines for the mean degradation paths, i. e.
$\mathbf{f}_{2}(t)=(1,t)^{T}$ and
$\boldsymbol{\Sigma}_{\gamma}=\left(\begin{array}[]{cc}\sigma_{1}^{2}&\rho\sigma_{1}\sigma_{2}\\\
\rho\sigma_{1}\sigma_{2}&\sigma_{2}^{2}\end{array}\right)$ where
$\sigma_{1}^{2},\sigma_{2}^{2}$ and $\rho$ are defined, 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}}},$
(4.3)

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 gives a 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 time $t=0$ is less than the failure threshold,
$\delta_{1}<y_{0}$.

In general, under the 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$. Note that
$\Phi(h(0))=\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$ where
$\lim_{t\to\infty}h(t)=\delta_{2}/\sigma_{2}$ and $t_{\alpha}$ is
nondegenerate for
$\Phi(-(y_{0}-\delta_{1})/\sigma_{1})<\alpha<\Phi(\delta_{2}/\sigma_{2})$. 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
$t_{\alpha}=(y_{0}-\delta_{1}+z_{\alpha}\sigma_{1})/\delta_{2}$.

The quantile $t_{\alpha}=t_{\alpha}(\boldsymbol{\theta})$ is a function of
both the fixed effect parameters $\boldsymbol{\beta}$ and the variance
parameters $\boldsymbol{\varsigma}$, in general.

By the delta-method the maximum likelihood estimator $\widehat{t}_{\alpha}$ of
${t}_{\alpha}$ is asymptotically normal with asymptotic variance

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

where $\mathbf{c}=\frac{\partial}{\partial\boldsymbol{\theta}}t_{\alpha}$ is
the gradient vector of $t_{\alpha}$ with respect to 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}$ is given by

$\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),$
(4.5)

in view of (4.2) by the implicit function theorem (see e. g. [7]), where
$c_{0}\neq 0$ is a constant not dprending on the design. 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, up to a multiplicative constant, 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 it is important that $\mathbf{c}_{\boldsymbol{\varsigma}}=\mathbf{0}$ in
the case of the median ($\alpha=0.5$). By the block diagonal form (3.1) of the
information matrix the asymptotic variance (4.4) of $\widehat{t}_{\alpha}$
becomes

$\mathrm{aVar}(\widehat{t}_{\alpha})=c_{0}^{2}\left(\mathbf{f}(\mathbf{x}_{u},t_{\alpha})^{T}\mathbf{M}_{\boldsymbol{\beta}}^{-1}\mathbf{f}(\mathbf{x}_{u},t_{\alpha})+z_{\alpha}^{2}\frac{\partial}{\partial\boldsymbol{\varsigma}}\sigma_{u}(t_{\alpha})^{T}\mathbf{M}_{\boldsymbol{\varsigma}}^{-1}\frac{\partial}{\partial\boldsymbol{\varsigma}}\sigma_{u}(t_{\alpha})\right)$
(4.6)

which simplifies to
$\mathrm{aVar}(\widehat{t}_{0.5})=c_{0}^{2}\mathbf{f}(\mathbf{x}_{u},t_{\alpha})^{T}\mathbf{M}_{\boldsymbol{\beta}}^{-1}\mathbf{f}(\mathbf{x}_{u},t_{\alpha})$
in the case of the median. For the case of a product-type model the expression
related to the fixed effect parameters $\boldsymbol{\beta}$ decomposes into

$\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}).$
(4.7)

##  5 Optimization of the time plan

In [14] we treated 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. There we noticed
that the optimal design for the stress variable does not depend on the time
plan used. 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 2 we still assume that the same time
plan $\mathbf{t}=(t_{1},...,t_{k})^{T}$ is used for all units. By equations
(4.6) and (4.7), the asymptotic variance is given by

$\mathrm{aVar}(\widehat{t}_{\alpha})=c_{0}^{2}\left(\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})+\mathbf{c}_{\boldsymbol{\varsigma}}^{T}{\mathbf{M}}_{\boldsymbol{\varsigma}}^{-1}\mathbf{c}_{\boldsymbol{\varsigma}}\right)$
(5.1)

Here the dependence of the asymptotic variance on the settings at the time
plan $\mathbf{t}$ for the time variable comes in through the information
matrix $\mathbf{M}_{2}$ in the second marginal model as well as through the
information matrix ${\mathbf{M}}_{\boldsymbol{\varsigma}}$ for the variance
parameters. In general, the asymptotic variance is a compound criterion in
$\mathbf{f}$ which aims at minimizing a linear combination of a variance term
$\mathbf{f}_{2}(t_{\alpha})^{T}\mathbf{M}_{2}^{-1}\mathbf{f}_{2}(t_{\alpha})$
in the location parameters $\boldsymbol{\beta}$ and a variance term
$\mathbf{c}_{\boldsymbol{\varsigma}}^{T}\mathbf{M}_{\boldsymbol{\varsigma}}^{-1}\mathbf{c}_{\boldsymbol{\varsigma}}$
in the variance parameters. In this compound criterion the first coefficient
involves the marginal information matrix $\mathbf{M}_{1}$ for the stress
variable. As a consequence the criterion and, hence, the optimal time plan
depend on the design for the stress variable. Thus first the marginal desiogn
for the stress variable has to be optimized which can be done simultaneously
in all time plans. Then the corresponding value has to be inserted into the
compound criterion for optimizing the time plan. To facilitate further
calculations we restrict to the case of estimating the median failure time
$t_{0.5}$. Consequently, the asymptotic variance simplifies to

$\mathrm{aVar}(\widehat{t}_{0.5})=c_{0}^{2}\mathbf{f}_{1}(\mathbf{x}_{u})^{T}\mathbf{M}_{1}^{-1}\mathbf{f}_{1}(\mathbf{x}_{u})\cdot\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}^{-1}\mathbf{f}_{2}(t_{0.5})\,.$
(5.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. Only the marginal $c$-criterion
$\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}^{-1}\mathbf{f}_{2}(t_{0.5})$ has to
be minimized. Remember that the marginal information matrix $\mathbf{M}_{2}$
is given by $\mathbf{M}_{2}=\mathbf{F}_{2}^{T}\mathbf{V}^{-1}\mathbf{F}_{2}$
where the variance covariance matrix $\mathbf{V}$ of the vector
$\mathbf{Y}_{i}$ of measurements for each unit $i$ is given by
$\mathbf{V}=\mathbf{F}_{2}{\boldsymbol{\Sigma}}_{\gamma}\mathbf{F}_{2}^{T}+{\boldsymbol{\Sigma}}_{\varepsilon}$.
In the situation of uncorrelated and homoscedastic measurment erros
($\boldsymbol{\Sigma}_{\varepsilon}=\sigma^{2}_{\varepsilon}\mathbf{I}_{k}$),
[12] provided a representation of the invere information matrix in random
effects models. This can be readily extended through the following result to a
general, non-singular variance covariance structure
$\boldsymbol{\Sigma}_{\varepsilon}$ for the measurement errors
$\boldsymbol{\varepsilon}$.

######  Lemma 5.1.

Let $\mathbf{F}$ be a $k\times p$ matrix of rank $p$,
$\boldsymbol{\Sigma}_{\gamma}$ a non-negative definite $p\times p$ matrix,
$\boldsymbol{\Sigma}_{\varepsilon}$ a positive definite $k\times k$ matrix,
$\mathbf{V}=\mathbf{F}\boldsymbol{\Sigma}_{\gamma}\mathbf{F}^{T}+\boldsymbol{\Sigma}_{\varepsilon}$,
and $\mathbf{M}=\mathbf{F}^{T}\mathbf{V}^{-1}\mathbf{F}$. Then

$\mathbf{M}^{-1}=(\mathbf{F}^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F})^{-1}+\boldsymbol{\Sigma}_{\gamma}.$
(5.3)

where the proof of Lemma 5.1 is obtained in the Appendix. From this we obtain
for the marginal information matrix $\mathbf{M}_{2}$ that its inverse can be
decomposed to
$\mathbf{M}_{2}^{-1}=\mathbf{F}_{2}^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F}_{2}+\boldsymbol{\Sigma}_{\gamma}$.
As a consequence, the marginal $c$-criterion
$\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}^{-1}\mathbf{f}_{2}(t_{0.5})$ can be
split up into

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

where
$\mathbf{M}_{2}^{(0)}=\mathbf{F}_{2}^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F}_{2}$
is the information matrix in the marginal fixed effect model in the time
variable $t$. Hence, as for any linear criterion, the optimization of the
asymptotic variance (5.2) of the median failure time ${t}_{0.5}$ with respect
to the time plan $\mathbf{t}$ 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)}}^{-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$. This leads to
the following result for optimization with respect to the time variable.

######  Proposition 5.1.

If the time plan $\mathbf{t}$ is $c$-optimal for extrapolation of the mean
response at $t_{0.5}$ in the marginal fixed effect model with covariance
$\boldsymbol{\Sigma}_{\varepsilon}$ for the time variable $t$, then
$\mathbf{t}$ 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 is, hence, local. 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$, as introduced by[6], 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}$ this leads to
essentially the same extrapolation problem as for the stress variable in [14].
Under the additional assumptions of uncorrelated and homoscedastic measurement
errors and $t_{0.5}>=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 [13], Example 2.1). Similarly, as
for extrapolation of the stress variable in [14], 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 often 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. For instance, the different measurement times
$t_{1},...,t_{k}$ may be assumed 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. For the $D$-criterion [11]
showed 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 strictly between $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
consideration. 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$. For brevity we consider in the following
Example the model of Example 2 in [14] for numerical results.

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$

######  Example 1.

In this model we have straight line regression in both a single stress
variable $x$ and the time variable $t$ with an interaction term $xt$. The
response $y_{ij}$ of testing unit $i$ at time $t_{j}$ is given by

$y_{ij}=\beta_{i1}+\beta_{2}x_{i}+\beta_{i3}t_{j}+\beta_{4}x_{i}t_{j}+\varepsilon_{ij}$
(5.5)

under the nominal values of Table 1 which are reproduced from Example 7.2 in
[18] after standardization. The measurement errors are supposed to be
uncorrelated and homoscedastic
($\boldsymbol{\Sigma}_{\varepsilon}=\sigma^{2}_{\varepsilon}\mathbf{I}_{k}$),
and the measurement errors $\varepsilon_{ij}$ are realizations of a normally
distributed error variable with mean zero and error variance
$\sigma^{2}_{\varepsilon}$. Further, the unit specific parameters
$(\beta_{i1},\beta_{i3})^{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}_{\gamma}=\left(\begin{array}[]{cc}\sigma_{1}^{2}&\rho\sigma_{1}\sigma_{2}\\\
\rho\sigma_{1}\sigma_{2}&\sigma_{2}^{2}\end{array}\right)$. We are searching
for a locally $c$-optimal design $\tau^{*}$ for extrapolation at the median
failure time $t_{0.5}=1.583$. 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. [16] ) 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.055&0.166&0.166&0.166\end{array}\right)\,.$

The optimal approximate design $\tau^{*}$ is supported on seven time points
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=0.10$ and $t=0.85$ separating the grid points
$t=0.00,0.05,0.90,0.95$, and $1.00$ with full weight $1/k$ from those with
zero weight $(t=0.15,...,0.80).$ This shape of the optimal design is in
accordance with the findings of [11] for $D$-optimality. In view of
Proposition 5.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 with equal weights $1/k=1/6$
by transferring the weight from the boundary point $t=0.85$ to the boundary
point $t=0.10$ (see [2]). 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}$ can be measured in terms
of the local $c$-efficiency

$\mathrm{eff}_{c}(\tau_{0})=\frac{\mathrm{aVar}(\hat{t}_{0.5};\tau^{*})}{\mathrm{aVar}(\hat{t}_{0.5};\tau_{0})}=\frac{\mathbf{f}_{2}(\hat{t}_{0.5})^{T}\mathbf{M}_{2}(\tau^{*})^{-1}\mathbf{f}_{2}(\hat{t}_{0.5})}{\mathbf{f}_{2}(\hat{t}_{0.5})^{T}\mathbf{M}_{2}(\tau_{0})^{-1}\mathbf{f}_{2}(\hat{t}_{0.5})}=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 favorably used as the time
plan for conducting an accelerated degradation testing experiment.

Note that by Proposition 5.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, see equation (5.4). Nevertheless, the $c$-efficiency

$\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}^{(0)}(\tau^{*})^{-1}\mathbf{f}_{2}(t_{0.5})/\mathbf{f}_{2}(t_{0.5})^{T}\mathbf{M}_{2}^{(0)}(\tau_{0})^{-1}\mathbf{f}_{2}(t_{0.5})$

for extrapolation at $t_{0.5}$ in the fixed effect model provides a lower
bound for the efficiency in estimating the median failure time. 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 [10] and [9].

##  6 Cross-sectional time plans for destructive testing

In the case of destructive testing only $k=1$ measurement will be available
per unit. Hence, measurements have to be taken at different time points at
different units in order to guarantee estimability of the parameters,

$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}\,,$
(6.1)

$i=1,...,n$, where all other expressions have the same meaning as in the model
equation 2.5. 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$,
$\mathrm{Var}(Y_{i})=\sigma^{2}(t_{i})$. 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 [4]. 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 (6.1). 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 4, 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 (4.7). The experimental settings
$\mathbf{x}_{i}$ for the stress variable and $t_{i}$ for the time of
measurement may be chosen independently from the standardized experimental
regions $\mathcal{X}$ and $\mathcal{T}$, respectively. Here the approximate
designs $\zeta$ 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$. 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, in
accordance with [14], a $c$-optimal design $\zeta^{*}$ for
$\mathbf{c}=\mathbf{c}_{1}\otimes\mathbf{c}_{2}$ can be obtained as the
product-type design generated from the $c$-optimal designs $\xi^{*}$ for
$\mathbf{c}_{1}$ and $\tau^{*}$ for $\mathbf{c}_{2}$ in the associated
marginal models, respectively.

######  Theorem 6.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 with regression functions
$\tilde{\mathbf{f}}_{2}(t)=\mathbf{f}_{2}(t)/\sigma(t)$ 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 (6.1).

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 [14], 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)$.

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 the locally optimal designs 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. In the particular case of straight line degradation paths, as
in Example 1, the optimal marginal design $\tau^{*}$ is concentrated on the
endpoints 0 and 1 of the standardized time interval, as long as the median
$t_{0.5}$ is larger than the upper bound 1 of the interval. Then only the
optimal weight
$\pi^{*}=t_{0.5}\sigma(1)/(t_{0.5}\sigma(1)+(t_{0.5}-1)\sigma(0))$ at $t=1$
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 support points.

As in Example 1, the asymptotic eficiency

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

of the design $\zeta$ for estimating $t_{0.5}$, factorizes,

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

in the case of product-type designs $\zeta=\xi\otimes\tau$ into the marginal
efficiency $\mathrm{eff}_{c1}(\xi)$ of $\xi$ for extrapolation at
$\mathbf{x}_{u}$ and the marginal efficiency
$\mathrm{eff}_{c2}(\tau;\boldsymbol{\theta})$ of $\tau$ for extrapolation at
$t_{0.5}$. For a product-type design $\xi^{*}\otimes\tau$ with optimal first
marginal design $\xi^{*}$, the asymptotic efficiency reduces to the
$c$-efficiency of $\tau$ in the second marginal model,

$\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}};\boldsymbol{\theta})^{-1}\mathbf{c}_{2}(\boldsymbol{\theta})}{\mathbf{c}_{2}(\boldsymbol{\theta})^{T}\mathbf{M}_{2}(\tau;\boldsymbol{\theta})^{-1}\mathbf{c}_{2}(\boldsymbol{\theta})}\,.$

The vector $\mathbf{c}_{2}(\boldsymbol{\theta})=\mathbf{f}_{2}(t_{0.5})$
depends on the location parameters $\boldsymbol{\beta}$ through $t_{0.5}$, but
not on $\boldsymbol{\varsigma}$, while 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. Example 1), 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}}$. In the following we revisit Example 1
for obtaining an optimal design in regards to the current model of cross-
sectional time plans.

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

Under the current case of destructive testing, we provide in this example a
$c$-optimal design of the time variable for the degradation model presented in
Example 1. For the stress variable, the $c$-optimal marginal design $\xi^{*}$
is derived in [14]. In particular, the optimal marginal design $\xi^{*}$ is
concentrated on the endpoints $0$ and $1$ of the standardized design region
with weight $w^{*}=0.05$ at $x=1$. 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 ([3]) applied to the vector of weighted
regression functions $\tilde{\mathbf{f}}_{2}(t)=(1,t)^{T}/\sigma(t)$, where
the variance function $\sigma^{2}(t)$ is defined here by
$\sigma^{2}(t)=\sigma_{1}^{2}+\rho\sigma_{1}\sigma_{2}t+\sigma_{2}^{2}t^{2}+\sigma_{\varepsilon}^{2}$.
The shape of the Elfving set is exhibited in Figure 1.

Figure 1: 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)

The ray $\lambda\mathbf{f}_{2}(t_{0.5})=\lambda(1,t_{0.5})^{T}$ intersects the
boundary of the Elfving set at the upper 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 $0$ and $1$ of the
standardized experimental time interval, 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=1$
and $1-\pi^{*}$ at $t=0$.

By Theorem 6.1, the design

$\zeta^{*}=\xi^{*}\otimes\tau^{*}=\left(\begin{array}[]{cccc}(0,0)&(0,1)&(1,0)&(1,1)\\\
(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$ and $\pi^{*}=0.77$.
Then, 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)\,.$

In Figure 3, the optimal weight
$\pi^{*}=t_{0.5}\sigma(1)/(t_{0.5}\sigma(1)+(t_{0.5}-1)\sigma(0))$ is plotted
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 of Table 1.

Figure 2: Optimal weights $\pi^{*}$ in dependence on $t_{0.5}$ for Example 1

Figure 3: Optimal weights $\pi^{*}$ in dependence on $\sigma(1)/\sigma(0)$ 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$. In Figure 3 the
optimal weight $\pi^{*}$ is plotted as a function of the ratio
$\sigma(1)/\sigma(0)$ while $t_{0.5}$ is held fixed to $1.583$. On the other
hand the optimal weight $\pi^{*}$ increases in the ratio $\sigma(1)/\sigma(0)$
from $0$ for the ratio close to $0$ to $1$ when the ratio tends to infinity.
Those limiting values are indicated in Figures 3 and 3 by horizontal dashed
lines, and the nominal values for $t_{0.5}=1.583$ and
$\sigma(1)/\sigma(0)=1.22$ are indicated by vertical dotted lines,
respectively.

To assess the sensitivity of the optimal design with respect to parameter
misspecification, we plot the efficiency of the locally optimal design
$\zeta^{*}=\xi^{*}\otimes\tau^{*}$ in Figures 5 and 5. For a comparison with
standard marginal designs on the time scale, we also plot there the efficiency
of the designs $\xi^{*}\otimes\bar{\tau}_{2}$ (dashed line) and
$\xi^{*}\otimes\bar{\tau}_{6}$ (dashed and dotted line), 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 $t=0.0$,
$0.2$, $0.4$, $0.6$, $0.8$, and $1.0$.

Figure 4: Efficiency of $\zeta^{*}=\xi^{*}\otimes\tau^{*}$ (solid line),
$\xi^{*}\otimes\bar{\tau}_{2}$ (dashed line) and
$\xi^{*}\otimes\bar{\tau}_{6}$ (dashed and dotted line) in dependence on
$t_{0.5}$ for Example 1

Figure 5: Efficiency of $\zeta^{*}=\xi^{*}\otimes\tau^{*}$ (solid line),
$\xi^{*}\otimes\bar{\tau}_{2}$ (dashed line) and
$\xi^{*}\otimes\bar{\tau}_{6}$ (dashed and dotted line) in dependence on
$\sigma(1)/\sigma(0)$ for Example 1

In Figure 5 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 Figure 5 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 of the standard deviations 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^{2}(1)$ at the end of the experimental time interval is
larger than the variance $\sigma^{2}(0)$ at the initial time,
$\sigma(1)/\sigma(0)\geq 1$. This is satisfied in the case of positive or no
correlation, and it also holds for slight to moderate negative corrleations,
$\rho\geq-\sigma_{1}/(2\sigma_{2})$. In total, the locally optimal design
$\zeta^{*}=\xi^{*}\otimes\tau^{*}$ under the nominal values performs quite
well over a wide range of parameters while the design
$\xi^{*}\otimes\bar{\tau}_{2}$ which is concentrated on the same time points
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. It has to be noted for completeness
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. To
accomplish for this, the plot in Figure 5 has been generated by fixing
$\sigma_{1}^{2}=\sigma_{2}^{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.

##  7 Discussion and conclusion

An elaborate assessment of the reliability related characteristics is of
importance during the design stage of highly reliable systems. The variability
between units in accelerated degradation testing leads to a degradation
function that can be described by a mixed-effects linear model. This also
results in a nondegenerate distribution of the failure time under normal use
conditions such that it is desirable to estimate quantiles of this failure
time distribution as a characteristic of the reliability of the product. The
purpose of optimal experimental design of the time plan is then to obtain the
best settings for the time variable to obtain most accurate estimates for
these quantities.

In the present model for accelerated degradation testing, it is assumed that
in the experiment a cross-sectional design between units has to be chosen for
the stress variable while for repeated measurements the time variable varies
according to a longitudinal design within units.

Here we assumed a model with complete interactions between the time and the
stress variables and random effects only associated with time but not with
stress. Then the cross-sectional design for the stress variables and the
longitudinal design for the time variable can be optimized independently, and
the resulting common optimal design can be generated as the cross-product of
the optimal marginal designs for stress and time, respectively. In particular,
the same time plan for measurements can be used for all units in the test.

In contrast to that, in a model of destructive testing also the time variable
has to be chosen cross-sectionally. There the optimal choice of measurement
times may also be affected by the variance covariance parameters of the random
effects. In both cases (longitudinal and cross-sectional time settings) the
efficiency of the designs considered factorizes which facilitates to assess
their performance when the nominal values for these parameters are
misspecified at the design stage.

##  A Appendix

###### Proof.

We prove the statement of the Lemma 5.1 by showing that multiplication of
$\mathbf{M}$ with the right hand side
$\mathbf{C}=(\mathbf{F}^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F})^{-1}+\boldsymbol{\Sigma}_{\gamma}$
results in the $k\times k$ identity matrix $\mathbf{I}_{k}$. For this note
first that after premultiplication with $\mathbf{F}$ the right hand side can
be expanded to

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

Hence, by straightforward multiplication

$\mathbf{M}\mathbf{C}=\mathbf{F}^{T}\mathbf{V}^{-1}\mathbf{F}\mathbf{C}=\mathbf{F}^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F}(\mathbf{F}^{T}\boldsymbol{\Sigma}_{\varepsilon}^{-1}\mathbf{F})^{-1}=\mathbf{I}_{k}$

which proofs the lemma. ∎

## Acknowledgement

The work of the first author has been supported by the German Academic
Exchange Service (DAAD) under grant no. 2017-18/ID-57299294.

## References

## References

  * Dette et al. [2010] Dette, H., Pepelyshev, A., and Holland-Letz, T. (2010).  Optimal designs for random effect models with correlated errors with applications in population pharmacokinetics.  The Annals of Applied Statistics, 4(3):1430–1450.
  * Dorfleitner and Klein [1999] Dorfleitner, G. and Klein, T. (1999).  Rounding with multiplier methods: An efficient algorithm and applications in statistics.  Statistical Papers, 40(2):143–157.
  * Elfving [1952] Elfving, G. (1952).  Optimum allocation in linear regression theory.  The Annals of Mathematical Statistics, 23(2):255–262.
  * Graßhoff et al. [2012] Graßhoff, U., Doebler, A., Holling, H., and Schwabe, R. (2012).  Optimal design for linear regression models in the presence of heteroscedasticity caused by random coefficients.  Journal of Statistical Planning and Inference, 142(5):1108–1113.
  * Haghighi and Bae [2015] Haghighi, F. and Bae, S. J. (2015).  Reliability estimation from linear degradation and failure time data with competing risks under a step-stress accelerated degradation test.  IEEE Transactions on Reliability, 64(3):960–971.
  * Kiefer [1959] Kiefer, J. (1959).  Optimum experimental designs.  Journal of the Royal Statistical Society, Series B, 21(2):272–304.
  * Krantz and Parks [2012] Krantz, S. G. and Parks, H. R. (2012).  The implicit function theorem: history, theory, and applications.  Springer.
  * Limon et al. [2017] Limon, S., Yadav, O. P., and Liao, H. (2017).  A literature review on planning and analysis of accelerated testing for reliability assessment.  Quality and Reliability Engineering International, 33(8):2361–2383.
  * Müller [2007] Müller, W. G. (2007).  Collecting Spatial Data : Optimum Design of Experiments for Random Fields. 3rd rev. and extended ed. Springer.
  * Näther [1985] Näther, W. (1985).  Effective Observation of Random Fields.  Teubner.
  * Sahm and Schwabe [2001] Sahm, M. and Schwabe, R. (2001).  A note on optimal bounded designs.  In Optimum Design 2000, pages 131–140. Kluwer.
  * Schmelter [2007] Schmelter, T. (2007).  The optimality of single-group designs for certain mixed models.  Metrika, 65(2):183–193.
  * Schwabe [1996] Schwabe, R. (1996).  Optimum designs for multi-factor models.  Springer.
  * Shat and Schwabe [2021] Shat, H. and Schwabe, R. (2021).  Experimental designs for accelerated degradation tests based on linear mixed effects models.  arXiv preprint arXiv:2102.09446.
  * Silvey [1980] Silvey, S. D. (1980).  Optimal design, volume 1.  Chapman and Hall.
  * Silvey et al. [1978] Silvey, S. D., Titterington, D. H., and Torsney, B. (1978).  An algorithm for optimal designs on a design space.  Communications in Statistics – Theory and Methods, 7(14):1379–1389.
  * Tsai et al. [2016] Tsai, T.-R., Sung, W.-Y., Lio, Y. L., Chang, S. I., and Lu, J.-C. (2016).  Optimal two-variable accelerated degradation test plan for gamma degradation processes.  IEEE Transactions on Reliability, 65(1):459–468.
  * Weaver and Meeker [2014] Weaver, B. P. and Meeker, W. Q. (2014).  Methods for planning repeated measures accelerated degradation tests.  Applied Stochastic Models in Business and Industry, 30(6):658–671.
  * Zhang et al. [2015] Zhang, C., Lu, X., Tan, Y., and Wang, Y. (2015).  Reliability demonstration methodology for products with Gamma Process by optimal accelerated degradation testing.  Reliability Engineering and System Safety, 142(C):369–377.