# Bayesian $D$-Optimal Design of Experiments with Quantitative and Qualitative
Responses

Lulu Kang1 1Illinois Institute of Technology, Chicago, IL, 2Virginia Tech,
Blacksburg, VA.  Xinwei Deng2 1Illinois Institute of Technology, Chicago, IL,
2Virginia Tech, Blacksburg, VA.  Ran Jin2 1Illinois Institute of Technology,
Chicago, IL, 2Virginia Tech, Blacksburg, VA.

###### Abstract

Systems with both quantitative and qualitative responses are widely
encountered in many applications. Design of experiment methods are needed when
experiments are conducted to study such systems. Classic experimental design
methods are unsuitable here because they often focus on one type of response.
In this paper, we develop a Bayesian $D$-optimal design method for experiments
with one continuous and one binary response. Both noninformative and conjugate
informative prior distributions on the unknown parameters are considered. The
proposed design criterion has meaningful interpretations regarding the
$D$-optimality for the models for both types of responses. An efficient point-
exchange search algorithm is developed to construct the local $D$-optimal
designs for given parameter values. Global $D$-optimal designs are obtained by
accumulating the frequencies of the design points in local $D$-optimal
designs, where the parameters are sampled from the prior distributions. The
performances of the proposed methods are evaluated through two examples.
KEYWORDS: Bayesian $D$-optimal design, Conjugate prior, Generalized linear
model, Multivariate responses, Noninformative prior, Point-exchange.

##  1 Introduction

In many applications, both quantitative and qualitative responses are often
collected for evaluating the quality of the system. Often, the two types of
responses are mutually dependent. We call such a system with both types of
quality responses quantitative-qualitative system. Such systems are widely
encountered in practice (Kang et al., 2022, 2021, 2018). In Kang et al.
(2018), the authors studied an experiment of the lapping stage of the wafer
manufacturing process. The qualitative response is the conformity of the site
total indicator reading (STIR) of the wafer, which has two possible outcomes:
whether or not the STIR of a wafer is within the tolerance. The quantitative
response is the total thickness variation (TTV) of the wafer. Kang et al.
(2021) focused on the birth records and examined the mutual dependency of
birth weight and preterm birth. The birth weight of an infant is a
quantitative outcome and the preterm birth is a binary indicator of whether an
infant is born before 36 gestational weeks. The two types of outcomes are
correlated as an infant is usually underweight if the infant is born preterm.
In Kang et al. (2022), two case studies of quantitative-qualitative systems
from material sciences and gene expressions are illustrated. In the gene
expression study, the qualitative response has three possible outcomes:
healthy individuals, patients with Crohn’s disease, and patients with
Ulcerative colitis.

This work is motivated by a study of an etching process in a wafer
manufacturing process. In the production of silicon wafers, the silicon ingot
is sliced into wafers in fine geometry parameters. Inevitably, this step
leaves scratches on the wafers’ surface. An etching process is used to improve
the surface finish, during which the wafers are submerged in the container of
etchant for chemical reaction. The quality of the wafers after etching is
measured by two response variables: the total thickness variation of the wafer
(TTV) and the binary judgment that whether the wafer has cloudy stains in its
appearance. The two responses measure the quality from different but connected
aspects. There is a hidden dependency between the continuous TTV and binary
judgment of stains. To improve the etching quality, expensive experiments are
to be carried out to reveal this hidden dependency and to model the quality-
process relationship. Therefore, an ideal experimental design for continuous
and binary responses should be able to extract such useful information with
economic run size.

The classic design of experiments methods mainly focus on experiments with a
single continuous response. There have been various methods developed for a
single discrete response too, including Sitter and Torsney (1995); Woods et
al. (2006); Russell et al. (2009); Woods (2010); Woods and Van de Ven (2012).
For multiple responses, Draper and Hunter (1966) proposed the seminal work for
continuous responses, Denman et al. (2011) developed a design method for
bivariate binary responses modeled by Copula functions. In the case of mixed
types of responses, the literature is very scarce. A naive design method is to
combine the two designs that are separately constructed for each type of
response. However, such a naive strategy could be reasonable for one type of
response but problematic for the other by ignoring the dependency between the
types of responses, as shown in Example 1.

######  Example 1.

Denote by $Y$ and $Z$ a continuous response and a binary response,
respectively. Assume that the true model of the binary response $Z$ is
$\mathbb{E}(Z|x)=\pi(x)=\exp(1+x)/(1+\exp(1+x))$. The true model of $Y$ is
related to $Z$ in the form $Y|Z=z\sim N(1-(1-z)x^{2},0.3^{2})$. Thus,
$\mathbb{E}(Y|Z=1)=1$ and $\mathbb{E}(Y|Z=0)=1-x^{2}$. Using the naive design
method, a 14-point design is constructed, which consists of an 8-point local
$D$-optimal design for the model of $Z$ with
$\log(\pi(x)/(1-\pi(x))=\eta_{0}+\eta_{1}x$, and a 6-point $D$-optimal design
for linear regression with a quadratic model of $x$. Given the design, we
generate the responses from the true models of $Y$ and $Z$. Figure 1 (a)-(c)
show the data $(x_{i},y_{i},z_{i})$ from the 8-point, 6-point and their
combined 14-point design, respectively.

Figure 1: (a) Observations from design for $Z$; (b) Observations from design
for $Y$; (c) Observations from the combined design. Dashed line “- - -”
denotes $\mathbb{E}(Y|Z=1)=1$; solid line “—” denotes
$\mathbb{E}(Y|Z=0)=1-x^{2}$; point “+” denotes $(x_{i},y_{i})$ with $z_{i}=1$;
point “o” denotes $(x_{i},y_{i})$ with $z_{i}=0$.

In this example, there is a strong dependency between the two responses since
the true underlying models of $\mathbb{E}(Y|Z)$ are different when $Z=1$ and
$Z=0$. In both designs for a single response shown in Figure 1 (a) and (b),
the design points are balanced and reasonably distributed for the targeted
response. However, since there are no $Y$ observations for $Z=0$ at $x=1.0$
shown in Figure 1 (c), the quadratic model for $Y|Z=0$ is not estimable.
Clearly, the combined design is not suitable here. Note that this problem is
not caused by outliers, since all the points for $Z=1$ (with “+”) are varied
around $Y=1$ and the points for $Z=0$ (with “o”) are around $Y=1-x^{2}$. In
fact $P(Z=0|x=1)=0.12$, which is relatively small. Thus it is less likely to
observe $Y$ with $Z=0$ at $x=1.0$. A simple solution is to add more
replications at $x=1.0$, but it is not clear how many replications should be
sufficient. It becomes more difficult to spot a direct solution when the
experiments get more complicated.

Such experiments call for new experimental design methods to account for both
continuous and binary responses. Note that under the experimental design
framework, the linear model is often considered for modeling the continuous
response, and the generalized linear model (GLM) is often considered for
modeling the qualitative response. A joint model must be developed to
incorporate both types of responses. Compared to the classic design methods
for linear models or GLMs, design for the joint model is more challenging due
to the following aspects. First, the design criterion for the joint model is
more complicated, as the joint model is more complicated than the separate
ones. Second, experimental design for the GLM itself is more difficult than
that for the linear model, which is naturally inherited by the design for the
joint model. Third, efficient design construction algorithms are needed to
handle the complexity of the design criterion based on the joint model. Kang
and Huang (2019) proposed an $A$-optimal design for the experiments with both
quantitative and qualitative responses. The $A$-optimality was derived under a
Bayesian framework proposed in Kang et al. (2018). Although Kang and Huang
(2019) addressed the three challenges to a degree, the $A$-optimality is not a
commonly used criterion. More importantly, only informative prior is
considered, which circumvented some difficulties brought by noninformative
prior of the parameters.

In this paper, we choose the most commonly used $D$-optimal design criterion
and propose a novel Bayesian design method for the continuous and binary
responses. The proposed method considers both cases of noninformative priors
and informative priors. With the noninformative priors, the Bayesian framework
is equivalent to the frequentist approach. In this case, we also establish
some regularity conditions on the experimental run sizes. With the informative
priors, we develop the $D$-optimal design using conjugate priors. The derived
design criterion has meaningful interpretations in terms of the $D$-optimality
criteria for the models of both continuous and binary responses. Moreover, we
develop an efficient point-exchange algorithm to construct the proposed
designs. The construction algorithm can be applied to more general settings
other than factorial designs.

The rest of the paper is organized as follows. Section 2 reviews the general
Bayesian quantitative-qualitative (QQ) model and the optimal design criterion.
The Bayesian $D$-optimal design criterion is derived using noninformative
prior distributions in Section 3. In Section 4, the design criterion is
derived with conjugate informative priors. Efficient algorithms for
constructing optimal designs are elaborated in Section 5. One artificial
example and the etching experimental design are shown in Section 6. Section 7
concludes this paper with some discussions.

##  2 General Bayesian QQ Model and Design

We first review the general Bayesian QQ model introduced in Kang et al. (2018)
and focus on the scenario that $Y$ is a continuous response and $Z$ is a
binary response. The input variable
$\bm{x}=(x_{1},\ldots,x_{p})^{\prime}\in\mathbb{R}^{p}$ contains $p$
dimensions. Denote the data as $(\bm{x}_{i},y_{i},z_{i}),i=1,\ldots,n$, where
$y_{i}\in\mathbb{R}$ and $z_{i}\in\\{0,1\\}$. The vectors
$\bm{y}=(y_{i},\ldots,y_{n})^{\prime}$ and
$\bm{z}=(z_{1},\ldots,z_{n})^{\prime}$ are the vectors of response
observations. To jointly model the continuous response $Y$ and the binary
response $Z$ given $\bm{x}$, consider the joint probability of $Y|Z$ and $Z$.
The conditional model on $Y|Z$ is assumed to be a linear regression model,
while the model of $Z$ is a logistic regression model. Specifically, we
consider joint modeling of $Y$ and $Z$ as follows,

$\displaystyle Z=\left\\{\begin{array}[]{cl}1,&\textrm{ with probability
}\pi(\bm{x})\\\ 0,&\textrm{ with probability
}1-\pi(\bm{x})\end{array}\right.\textrm{ with
}\pi(\bm{x},\bm{\eta})=\frac{\exp(\bm{f}(\bm{x})^{\prime}\bm{\eta})}{1+\exp(\bm{f}(\bm{x})^{\prime}\bm{\eta})},$
(3)

where $\bm{f}(\bm{x})=\left(f_{1}(\bm{x}),\ldots,f_{q}(\bm{x})\right)$
contains $q$ modeling effects including the intercept, the main, interaction
and quadratic effects, etc., and
$\bm{\eta}=(\eta_{1},\ldots,\eta_{q})^{\prime}$ is a vector of parameter
coefficients. Conditioning on $Z=z$, the quantitative variable $Y$ has the
distribution

$Y|Z=z\sim
N(\mu_{0}+z\bm{f}(\bm{x})^{\prime}\bm{\beta}^{(1)}+(1-z)\bm{f}(\bm{x})^{\prime}\bm{\beta}^{(2)},\sigma^{2}),$
(4)

where
$\bm{\beta}^{(i)}=(\beta^{(i)}_{1},\ldots,\beta^{(i)}_{q})^{\prime},i=1,2$ are
the corresponding coefficients of the model effects. The parameter $\mu_{0}$
is the mean and $\sigma^{2}$ is the noise variance. The above conditional
model (4) indicates that $Y|Z=1\sim
N(\mu_{0}+\bm{f}(x)^{\prime}\bm{\beta}^{(1)},\sigma^{2})$ and $Y|Z=0\sim
N(\mu_{0}+\bm{f}(\bm{x})^{\prime}\bm{\beta}^{(2)},\sigma^{2})$. We assume the
same variance $\sigma^{2}$ for the two conditional distributions of $Y|Z=1$
and $0$. The design method developed in the paper can be easily adapted to the
case with different variances.

The association between the two responses $Y$ and $Z$ is represented using the
conditional model $Y|Z$. When the two linear models for $Y|Z=0$ and $Y|Z=1$
are different, i.e., $\bm{\beta}^{(1)}\neq\bm{\beta}^{(2)}$, then it is
important to take account of the influence of the qualitative response $Z$
when modeling the quantitative response $Y$. Let
$\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{n})^{\prime}$ be the $n\times p$ design
matrix with $\bm{x}_{i}$ as the $i$th design point. Based on the CB model, we
can express the sampling distributions as

$\displaystyle\bm{y}|\bm{z},\bm{\beta}^{(1)},\bm{\beta}^{(2)},\mu_{0},\sigma^{2},\bm{X}\sim
N(\mu_{0}{\bf
1}+\bm{V}_{1}\bm{F}\bm{\beta}^{(1)}+\bm{V}_{2}\bm{F}\bm{\beta}^{(2)},\sigma^{2}\bm{I}_{n}),$
(5)
$\displaystyle\bm{z}|\bm{\eta},\bm{X}\sim\textrm{Bernoulli}(\pi(\bm{x}_{i},\bm{\eta}))\textrm{
for }i=1,\ldots,n,\textrm{ and }$ (6) $\displaystyle
p(\bm{z}|\bm{\eta},\bm{X})\propto\exp\left\\{\sum_{i=1}^{n}\left(z_{i}\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}-\log(1+e^{\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}})\right)\right\\},$

where $p(\cdot)$ denotes a general density function. Here
$\bm{V}_{1}=\mbox{diag}\\{z_{1},\ldots,z_{n}\\}$ is a diagonal matrix,
$\bm{I}_{n}$ is the $n\times n$ identity matrix and
$\bm{V}_{2}=\bm{I}_{n}-\bm{V}_{1}$, $\bm{F}$ is the model matrix with the
$i$th row as $\bm{f}(\bm{x}_{i})^{\prime}$, and $\bf 1$ is a vector of ones.

Denote $p(\bm{\beta}^{(1)})$, $p(\bm{\beta}^{(2)})$, and $p(\bm{\eta})$ as the
prior distributions of the parameters $\bm{\beta}^{(1)}$, $\bm{\beta}^{(2)}$,
and $\bm{\eta}$. Note that we focus on the estimation accuracy of the three
groups of parameters. The mean $\mu_{0}$ and variance $\sigma^{2}$ are
considered nuisance parameters and thus excluded from the optimal design
criterion. In this work, we assume that the priors of $\bm{\beta}^{(1)}$,
$\bm{\beta}^{(2)}$, and $\bm{\eta}$ are independent. Under this assumption,
the conditional posterior distribution of $\bm{\eta}$, $\bm{\beta}^{(1)}$, and
$\bm{\beta}^{(2)}$ are also independent as explained in Sections 3 and 4.
Under the Bayesian framework, the conditional posterior distribution of the
parameters $(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta})$ can be derived as

$\displaystyle
p(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X})$
$\displaystyle\propto
p(\bm{y}|\bm{z},\bm{\beta}^{(1)},\bm{\beta}^{(2)},\mu_{0},\sigma^{2},\bm{X})p(\bm{\beta}^{(1)})p(\bm{\beta}^{(2)})p(\bm{z}|\bm{\eta},\bm{X})p(\bm{\eta})$
(7)

Using (7) we develop the general Bayesian optimal design criterion. Let $\psi$
be a criterion function on the conditional posterior distribution of the
parameters. For example, it can be the Shannon information (or equivalently,
the Kullback-Leibler distance), $A$/$I$-optimality (Fedorov, 1972), or other
design criteria. However, $\psi(\cdot)$ cannot be directly used as the final
optimal design criterion because its value depends on the random parameters
$(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta})$ and the experimental outputs
$(\bm{y},\bm{z})$ that are not yet observed. The randomness of
$(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta})$ can be removed by calculating
the mean of $\psi$ with respect to these parameters. The uncertainty of
$(\bm{y},\bm{z})$ can be removed by calculating the mean
$\mathbb{E}(\mathbb{E}(\psi|\bm{y},\bm{z}))$. Therefore, the general Bayesian
optimal design criterion on the design matrix $\bm{X}$ is

$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})=\int
p(\bm{y},\bm{z}|\mu_{0},\sigma^{2},\bm{X})$
$\displaystyle\times\left(\int\psi(p(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X}))\times\right.$
(8)
$\displaystyle\left.p(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X})d\bm{\beta}^{(1)}d\bm{\beta}^{(2)}d\bm{\eta}\right)d\bm{y}d\bm{z}.$

It is well-known that the Bayesian $D$-optimal design is equivalent to the
Shannon information criterion (Chaloner and Verdinelli, 1995), omitting the
constant terms to $\bm{X}$. The criterion function $\psi(\cdot)$ of Shannon
information is $\log(\cdot)$. Next, we develop the Bayesian $D$-optimal design
criteria (8) under different prior distributions.

##  3 Optimal Design under Noninformative Priors

When lacking domain knowledge or proper historical data, experimenters often
favor the frequentist approach as no priors need to be specified. The
frequentist approach can be seen as the Bayesian approach using noninformative
priors. In this section, we derive the optimal design criterion and the
regularity conditions for noninformative priors.

###  3.1 Design Criterion

Assume the non-informative priors $p(\bm{\beta}^{(i)})\propto 1$ for $i=1,2$
and $p(\bm{\eta})\propto 1$. The conditional posterior distribution in (7) is
the same as the joint distribution of the data. It can be further factorized
into

$p(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X})\propto
p(\bm{\eta}|\bm{z},\bm{X})\prod_{i=1}^{2}p(\bm{\beta}^{(i)}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X}),$

with the posterior distributions

$\displaystyle\bm{\beta}^{(i)}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X}\sim
N\left((\bm{F}^{\prime}\bm{V}_{i}\bm{F})^{-1}\bm{F}^{\prime}\bm{V}_{i}(\bm{y}-\mu_{0}{\bf
1}),\sigma^{2}(\bm{F}^{\prime}\bm{V}_{i}\bm{F})^{-1}\right)\textrm{ for
}i=1,2,$ (9) $\displaystyle
p(\bm{\eta}|\bm{z},\bm{X})\propto\exp\left\\{\sum_{i=1}^{n}\left(z_{i}\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}-\log(1+e^{\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}})\right)\right\\}.$
(10)

Conditioning on $\bm{z}$, the posterior distributions of
$(\bm{\beta}^{1},\bm{\beta}^{2})$ and $\bm{\eta}$ are independent. Note that
the noninformative prior $p(\bm{\eta})\propto 1$ is proper because it leads to
proper posterior $p(\bm{\eta}|\bm{z},\bm{X})$. Under the noninformative
priors, the Bayesian estimation is identical to the frequentist estimation.

Using the posterior distributions (9)–(10) and the criterion function
$\psi(\cdot)=\log(\cdot)$ in the general Bayesian optimal design criterion
(8), we obtain the Bayesian $D$-optimal design criterion (11).

$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})$
$\displaystyle=\mathbb{E}_{\bm{z},\bm{\eta}}\left\\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\\}$
(11)
$\displaystyle+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\\{\log\det\\{(\bm{F}^{\prime}\bm{V}_{i}\bm{F})\\}\right\\}+\textrm{constants}.$

The derivation is in Appendix A1. The first additive term in (11) is exactly
the Bayesian $D$-optimal design criterion for GLMs. Unfortunately, its exact
integration is not tractable. The common approach in experimental design for
GLMs is to use a normal approximation for the posterior distribution
$p(\bm{\eta}|\bm{z},\bm{X})$ (Chaloner and Verdinelli, 1995; Khuri et al.,
2006). Such an approximation leads to

$\mathbb{E}_{\bm{z},\bm{\eta}}\left\\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\\}\approx\mathbb{E}_{\bm{\eta}}\\{\log\det\bm{I}(\bm{\eta}|\bm{X})\\}+\textrm{constant},$
(12)

where $\bm{I}(\bm{\eta}|\bm{X})$ is the Fisher information matrix. We can
easily show that

$\bm{I}(\bm{\eta}|\bm{X})=-\mathbb{E}_{\bm{z}}\left(\frac{\partial^{2}l(\bm{z},\bm{\eta}|\bm{X})}{\partial\bm{\eta}\partial\bm{\eta}^{T}}\right)=\sum_{i=1}^{n}\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))=\bm{F}^{\prime}\bm{W}_{0}\bm{F},$

where $\bm{W}_{0}$ is a diagonal weight matrix

$\bm{W}_{0}=\mbox{diag}\\{\pi(\bm{x}_{1},\bm{\eta})(1-\pi(\bm{x}_{1},\bm{\eta})),\ldots,\pi(\bm{x}_{n},\bm{\eta})(1-\pi(\bm{x}_{n},\bm{\eta}))\\}.$

Omitting the irrelevant constant, we approximate the exact criterion
$\Psi(\bm{X}|\mu_{0},\sigma^{2})$ in (11) as follows.

$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})\approx\mathbb{E}_{\bm{\eta}}\\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\\}+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\\{\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F})\right\\}.$
(13)

To construct the optimal design, we consider maximizing the approximated
$\Psi(\bm{X}|\mu_{0},\sigma^{2})$ in (13). But this is not trivial, because it
involves the expectation on $Z_{i}$’s in the second additive term. To overcome
this challenge, Hainy et al. (2013) constructed optimal designs by simulating
samples from the joint distribution of responses and the unknown parameters.
But this method can be computationally expensive for even slightly larger
dimensions of experimental factors. Instead of simulating $Z_{i}$’s, we derive
the following Theorem 1 that gives a tractable upper bound $Q(\bm{X})$. Thus
we propose using the upper bound $Q(\bm{X})$ as the optimal criterion.

######  Theorem 1.

Assume that the matrices $\bm{F}^{\prime}\bm{W}_{0}\bm{F}$,
$\bm{F}^{\prime}\bm{V}_{1}\bm{F}$, and $\bm{F}^{\prime}\bm{V}_{2}\bm{F}$ are
all nonsingular. Omitting the irrelevant constant, an upper bound of the
approximated $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ is

$Q(\bm{X})=\mathbb{E}_{\bm{\eta}}\left\\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})+\frac{1}{2}\sum_{i=1}^{2}\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F})\right\\},$
(14)

where
$\bm{W}_{1}=\mbox{diag}\\{\pi(\bm{x}_{1},\bm{\eta}),\ldots,\pi(\bm{x}_{n},\bm{\eta})\\}$
and $\bm{W}_{2}=\bm{I}_{n}-\bm{W}_{1}$.

The proof of Theorem 1 is in Appendix A1. Note that Theorem 1 requires that
$\bm{F}^{\prime}\bm{W}_{i}\bm{F}$ for $i=0,1,2$ are all nonsingular. Obviously
$\bm{W}_{0}=\bm{W}_{1}\bm{W}_{2}$. It is easy to see that
$\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ is nonsingular if and only if both
$\bm{F}^{\prime}\bm{W}_{1}\bm{F}$ and $\bm{F}^{\prime}\bm{W}_{2}\bm{F}$ are
nonsingular.

The matrices $\bm{F}^{\prime}\bm{V}_{1}\bm{F}$ and
$\bm{F}^{\prime}\bm{V}_{2}\bm{F}$ involve the responses $Z_{i}$’s that are not
yet observed at the experimental design stage. We can only choose the
experimental run size and the design points to avoid the singularity problem
with a larger probability for given values of $\bm{\eta}$. Once the run size
is chosen, the design points can be optimally arranged by maximizing
$Q(\bm{X})$. The weight matrix $\bm{W}_{1}$ (or $\bm{W}_{2}$) gives more
weight to the feasible design points that are more likely to lead to $Z=1$ (or
$Z=0$) observations so that the parameters $\bm{\beta}^{(1)}$ (or
$\bm{\beta}^{(2)}$) of the linear model $Y|Z=1$ (or $Y|Z=0$) are more likely
to be estimable. Next, we introduce some regularity conditions on the run size
and number of replications to alleviate the singularity problem.

###  3.2 Regularity Conditions

Let $m$ be the number of distinct design points in the design matrix $\bm{X}$,
$n_{i}$ be the number of repeated point $\bm{x}_{i}$ in $\bm{X}$ for
$i=1,\ldots,m$. Thus $n=\sum_{i=1}^{m}n_{i}$ and $n_{i}\geq 1$ for
$i=1,\ldots,m$. First, it is necessary that $m\geq q$ for the linear
regression model to be estimable under the noninformative priors. The if-and-
only-if condition for $\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ to be nonsingular is
that $\textrm{rank}(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\geq q$. If $m\geq q$ and
$\pi(\bm{x}_{i},\bm{\eta})\in(0,1)$ for $i=1,\ldots,m$, then
$\bm{F}^{\prime}\bm{W}_{i}\bm{F}$ for $i=0,1,2$ are all nonsingular and thus
positive definite. To make sure $\pi(\bm{x}_{i},\bm{\eta})\in(0,1)$ for
$i=1,\ldots,m$, it is sufficient to assume that $\bm{\eta}$ is finitely
bounded. This condition is typically used for the frequentist $D$-optimal
design for GLMs. For instance, Woods et al. (2006) suggested using the
centroids of the finite bounded space of $\bm{\eta}$ to develop the local
$D$-optimal design for GLMs. Dror and Steinberg (2006) clustered different
local $D$-optimal designs with $\bm{\eta}$ randomly sampled from its bounded
space.

The if-and-only-if condition for $\bm{F}^{\prime}\bm{V}_{1}\bm{F}$ and
$\bm{F}^{\prime}\bm{V}_{2}\bm{F}$ being nonsingular is

$\sum_{i=1}^{m}I(\sum_{j=1}^{n_{i}}Z_{ij}>0)\geq
q\quad\textrm{and}\quad\sum_{i=1}^{m}I(\sum_{j=1}^{n_{i}}Z_{ij}<n_{i})\geq q,$
(15)

where $Z_{ij}$ is the $j$th random binary response at the unique design point
$\bm{x}_{i}$ and $I(\cdot)$ is the indicator function. In the following, we
discuss how to choose sample sizes $n_{i}$ and $n$ under two scenarios (i)
$m=q$ and (ii) $m>q$.

######  Proposition 1.

Assume $m=q$. Both $\bm{F}^{\prime}\bm{V}_{1}\bm{F}$ and
$\bm{F}^{\prime}\bm{V}_{2}\bm{F}$ are nonsingular if and only if
$I(0<\sum_{j=1}^{n_{i}}Z_{ij}<n_{i})=1$ for $i=1,2,\ldots,m$. For any given
$\kappa\in(0,1)$, a sufficient condition on $n_{i}$ for
$\Pr(0<\sum_{j=1}^{n_{i}}Z_{ij}<n_{i})\geq\kappa$ is

$n_{i}\geq
1+\left\lceil\frac{\log(1-\kappa)}{\log\left(\max\left\\{\pi(\bm{x}_{i},\bm{\eta}),1-\pi(\bm{x}_{i},\bm{\eta})\right\\}\right)}\right\rceil\quad\textrm{for
}i=1,2,\ldots,m,$ (16)

and a necessary condition is

$n_{i}\geq\left\lceil\frac{2\log\left(\frac{1-\kappa}{2}\right)}{\log\pi(\bm{x}_{i},\bm{\eta})+\log(1-\pi(\bm{x}_{i},\bm{\eta}))}\right\rceil\quad\textrm{for
}i=1,2,\ldots,m.$ (17)

######  Proposition 2.

Assume $m>q$. To make both $\bm{F}^{\prime}\bm{V}_{1}\bm{F}$ and
$\bm{F}^{\prime}\bm{V}_{2}\bm{F}$ nonsingular with large probability, or
equivalently,

$\mathbb{E}\left\\{\sum_{i=1}^{m}I(\sum_{j=1}^{n_{i}}Z_{ij}>0)\right\\}\geq
q\quad\textrm{and}\quad\mathbb{E}\left\\{\sum_{i=1}^{m}I(\sum_{j=1}^{n_{i}}Z_{ij}<n_{i})\right\\}\geq
q,$

  * (i)

a sufficient condition is

$n_{0}\geq\max\left\lceil\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\min})},\frac{\log(1-q/m)}{\log\pi_{\max}}\right\\}\right\rceil,$
(18)

which is the same as

$n\geq\left\lceil
m\cdot\max\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\min})},\frac{\log(1-q/m)}{\log\pi_{\max}}\right\\}\right\rceil,$
(19)

where $n_{0}=\min\\{n_{1},\ldots,n_{m}\\}$,
$\pi_{\min}=\min_{i=1}^{m}\pi(\bm{x}_{i},\bm{\eta})$,
$\pi_{\max}=\max_{i=1}^{m}\pi(\bm{x}_{i},\bm{\eta})$ and $\bm{x}_{i}$’s are
the unique design points;

  * (ii)

a necessary condition is

$n_{0}\geq\left\lceil\max\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\max})},\frac{\log(1-q/m)}{\log\pi_{\min}}\right\\}\right\rceil$
(20)

which is the same as

$n\geq\left\lceil
m\cdot\max\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\max})},\frac{\log(1-q/m)}{\log\pi_{\min}}\right\\}\right\rceil.$
(21)

Proposition 1 gives a sufficient condition on the lower bound of $n_{i}$ when
saturated design ($m=q$) is used. Under the sufficient condition, the
nonsingularity of $\bm{F}\bm{V}_{1}\bm{F}$ and $\bm{F}\bm{V}_{2}\bm{F}$ holds
with a probability larger than $\kappa^{m}$. For Example 1 in Section 1,
suppose that if the possible values of $x$ can only be $-1,0,1$, then $m=q=3$.
Let $\bm{\eta}=(1,1)^{\prime}$. If $\kappa=0.5$, then the numbers of
replications for $x=-1,0,1$ need to satisfy $n_{1}\geq 2$, $n_{2}\geq 4$, and
$n_{3}\geq 7$, respectively. If $\kappa=0.9$, then $n_{1}\geq 5$, $n_{2}\geq
9$, and $n_{3}\geq 20$. Proposition 1 is useful in Step 1 to construct the
initial design in Algorithm 2 in Section 5.

Proposition 2 provides one sufficient condition and one necessary condition
when $m>q$ on the smallest number of replications and the overall run size.
But these conditions only examine the nonsingularity of the two matrices with
large probability, which is weaker than Proposition 1. For given $\bm{\eta}$
value, Algorithm 2 in Section 5.1 can return the local $D$-optimal design.
Proposition 2 can be useful to check the local $D$-optimal design, as
$\pi_{\min}$ and $\pi_{\max}$ depend on $\bm{\eta}$. Take the artificial
example in Section 6.1 for instance. The local $D$-optimal design for $\rho=0$
(Table A1 in Appendix A2) has $m=50$ unique design points and there are $q=22$
effects. According to Proposition 2, the sufficient condition requires
$n_{0}\geq 7$ and the necessary condition requires $n_{0}\geq 1$. The local
$D$-optimal design in Table A1 only satisfies the necessary condition. To meet
the sufficient condition $n$ has to be much larger. For the global optimal
design considering all possible $\bm{\eta}$ values, Proposition 2 can provide
some guidelines for the design construction when $\bm{\eta}$ is varied in a
relatively small range.

##  4 Optimal Design under Conjugate Priors

When prior information for parameters $\bm{\beta}^{(1)}$, $\bm{\beta}^{(2)}$,
and $\bm{\eta}$ is available, it would be desirable to consider the optimal
design under the informative priors. In this section, we detail the proposed
Bayesian $D$-optimal design using the conjugate priors.

###  4.1 Design Criterion

For the parameters $\bm{\beta}^{(1)}$ and $\bm{\beta}^{(2)}$, the conjugate
priors are normal distribution since $Y|Z$ follows normal distribution. Thus
we consider their priors as

$\bm{\beta}^{(1)}\sim N(\bm{0},\tau^{2}\bm{R}_{1}),\quad\bm{\beta}^{(2)}\sim
N(\bm{0},\tau^{2}\bm{R}_{2}).$

where $\tau^{2}$ is the prior variance and $\bm{R}_{i}$ is the prior
correlation matrix of $\bm{\beta}^{(i)}$ for $i=1,2$. Here we use the same
prior variance $\tau^{2}$ only for simplicity. The matrix $\bm{R}_{i}$ can be
specified flexibly such as using $(\bm{F}^{\prime}\bm{F})^{-1}$, or those in
Joseph (2006) for factorial designs.

For the parameters $\bm{\eta}$, we choose the conjugate prior derived in Chen
and Ibrahim (2003). It takes the form

$\bm{\eta}\sim
D(s,\bm{b})\propto\exp\left\\{\sum_{i=1}^{n}s\left(b_{i}\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}-\log(1+e^{\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}})\right)\right\\},$
(22)

where $D(s,\bm{b})$ is the distribution with parameters $(s,\bm{b})$. Here $s$
is a scalar factor and $\bm{b}\in(0,1)^{n}$ is the marginal mean of $\bm{z}$
as shown in Diaconis and Ylvisaker (1979). The value of $\bm{b}$ can be
interpreted as a prior prediction (or guess) for $\mathbb{E}(\bm{Z})$. Based
on the priors for $(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta})$ we can
derive the posteriors as follows.

######  Proposition 3.

For priors $\bm{\beta}^{(i)}\sim N(\bm{0},\tau^{2}\bm{R}_{i})$, $i=1,2$ , and
$\bm{\eta}\sim D(s,\bm{b})$, the posterior distributions of
$\bm{\beta}^{(1)}$, $\bm{\beta}^{(2)}$ and $\bm{\eta}$ are independent of each
other with the following forms.

$\displaystyle\bm{\beta}^{(i)}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X}\sim N$
$\displaystyle\left(\bm{H}_{i}^{-1}\bm{F}^{\prime}\bm{V}_{i}(\bm{y}-\mu_{0}{\bf
1}),\sigma^{2}\bm{H}_{i}^{-1}\right),\qquad\textrm{for }i=1,2.$ (23)
$\displaystyle\bm{\eta}|\bm{z},\bm{X}\sim D$
$\displaystyle\left(1+s,\frac{\bm{z}+s\bm{b}}{1+s}\right),$ (24)

where $\bm{H}_{i}=\bm{F}^{\prime}\bm{V}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}$ with
$\rho=\frac{\sigma^{2}}{\tau^{2}}$.

The proof of Proposition 3 can be derived following the standard Bayesian
framework, thus is omitted. To derive the Bayesian $D$-optimal design
criterion, we take the posterior distributions in Proposition 3 to (8) and set
$\psi(\cdot)=\log(\cdot)$. The derivation is very similar to that in (11), and
thus we obtain the exact design criterion as

$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})$
$\displaystyle=\mathbb{E}_{\bm{z},\bm{\eta}}\left\\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\\}$
(25)
$\displaystyle+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\\{\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F}+\rho\bm{R}_{i}^{-1})\right\\}+\textrm{constant}.$

As the integration of
$\mathbb{E}_{\bm{z},\bm{\eta}}\left\\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\\}$
is not tractable, we adopt the same normal approximation of the posterior
distribution $p(\bm{\eta}|\bm{z},\bm{X})$ as in (12). A straightforward
calculation leads to getting Fisher information matrix
$\bm{I}(\bm{\eta}|\bm{X})=(1+s)\bm{F}^{\prime}\bm{W}_{0}\bm{F}$. Thus we have

$\mathbb{E}_{\bm{z},\bm{\eta}}\left\\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\\}\approx\mathbb{E}_{\bm{\eta}}\\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\\}+\textrm{constant}.$
(26)

Disregarding the constant, we can approximate the exact
$\Psi(\bm{X}|\mu_{0},\sigma^{2})$ by

$\Psi(\bm{X}|\mu_{0},\sigma^{2})\approx\mathbb{E}_{\bm{\eta}}\\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\\}+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\\{\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F}+\rho\bm{R}_{i}^{-1})\right\\},$

The following Theorem 2 gives an upper bound of the approximated criterion
$\Psi(\bm{X}|\mu_{0},\sigma^{2})$ to avoid the integration with respect to
$\bm{z}$, which plays the same role as Theorem 1.

######  Theorem 2.

Assume that the prior distributions of $\bm{\beta}^{(i)}$ are
$\bm{\beta}^{(i)}\sim N(\bm{0},\tau^{2}\bm{R}_{i})$ for $i=1,2$ and
$\bm{\eta}$ has either the conjugate prior $\bm{\eta}\sim D(s,\bm{b})$ or the
noninformative prior $p(\bm{\eta})\propto 1$. If
$\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ is nonsingular, an upper bound of the
approximated $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ is

$Q(\bm{X})=\mathbb{E}_{\bm{\eta}}\left(\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})+\frac{1}{2}\sum_{i=1}^{2}\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1})\right).$
(27)

For the same argument as in Section 3, we use the upper bound in (27) as the
optimal design criterion. Note that since $\rho\bm{R}_{i}^{-1}$ is added to
$\bm{F}^{\prime}\bm{V}_{i}\bm{F}$ and $\bm{F}^{\prime}\bm{W}_{i}\bm{F}$,
$\bm{F}^{\prime}\bm{V}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}$ and
$\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}$ are nonsingular. The
derivation of $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ and $Q(\bm{X})$ only needs
$\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ to be nonsingular, which requires $m\geq q$
and $\bm{\eta}$ to be finitely bounded as in Section 3.

###  4.2 Interpretation

Note that the criterion in (27) has a similar formulation with $Q(\bm{X})$ in
(14). The only difference is that (14) does not involve $\rho\bm{R}_{i}^{-1}$.
For consistency, we use the formula (27) as the design criterion $Q(\bm{X})$
for both cases. When noninformative priors for $\bm{\beta}^{(1)}$ and
$\bm{\beta}^{(2)}$ are used, we set $\rho=0$. From another point of view, as
$\tau^{2}\rightarrow\infty$, $\rho\rightarrow 0$, the variances in the priors
$p(\bm{\beta}_{1})$ and $p(\bm{\beta}_{2})$ diffuse and result in a
noninformative priors.

The criterion $Q(\bm{X})$, consisting of three additive terms, can be
interpreted intuitively. The first additive term
$\mathbb{E}_{\bm{\eta}}\\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\\}$ is
known as the Bayesian $D$-optimal criterion for logistic regression and
$\mathbb{E}_{\bm{\eta}}\\{\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1})\\}$
is the Bayesian $D$-optimal criterion for the linear regression model of $Y$.
To explain the weights, we rewrite $Q(\bm{X})$ as follows.

$Q(\bm{X})=1\cdot\mathbb{E}_{\bm{\eta}}\left(\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\right)+1\cdot\left(\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}(\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}))\right)$

Since there are equal numbers of binary and continuous response observations,
the design criterion should put the same weight (equal to 1) on both design
criteria for $Z$ and $Y$. For the two criteria for the linear regression
models, the same weight 1/2 is used. This is also reasonable because we assume
$\pi_{i}\in(0,1)$. Then none of the diagonal entries of $\bm{W}_{1}$ and
$\bm{W}_{2}$ are zero, so the two terms should split the total weight 1
assigned for the entire linear regression part. Therefore, even though
$Q(\bm{X})$ are derived analytically, all the additive terms and their weights
make sense intuitively.

###  4.3 Prior Parameters

Note that the conjugate prior $p(\bm{\eta})$ requires prior parameters
$(s,\bm{b})$ to be specified. Moreover, the prior distribution (22) contains
$f(\bm{x}_{i})$, which depends on the design points. When sampling $\bm{\eta}$
from the prior (22), it does not matter whether $f(\bm{x}_{i})$’s are actually
from the design points. If relevant historical data is available, we can
simply sample $\bm{\eta}$ from the likelihood of the data. Alternatively, one
can adopt the method in Chen and Ibrahim (2003) to estimate the parameters
$(s,\bm{b})$. Without the relevant data, we would use the noninformative prior
for $\bm{\eta}$, i.e., $p(\bm{\eta})\propto 1$ in the bounded region for
$\bm{\eta}$.

The design criterion $Q(\bm{X})$ contains some unknown parameters, including
the noise-to-signal ratio $\rho=\sigma^{2}/\tau^{2}$ and the correlation
matrices $\bm{R}_{i}$’s for $i=1,2$. The value of $\rho$ has to be specified
either from the historical data or from the domain knowledge. Typically we
would assume $\rho<1$ such that the measurement error has a smaller variance
than the signal variance.

The setting of $\bm{R}_{i}$ can also be specified flexibly. If historical data
are available, $\bm{R}_{i}$ can be set as the estimated correlation matrix of
$\bm{\beta}^{(i)}$. Otherwise, we can use the correlation matrix in Joseph
(2006) and Kang and Joseph (2009), which is targeted for factorial design.
Specifically, let $\bm{\beta}$ be the unknown coefficients of the linear
regression model and the prior distribution is $\bm{\beta}\sim
N(\bm{0},\tau^{2}\bm{R})$. For 2-level factor coded in $-1$ and 1, Joseph
(2006) suggests that $\bm{R}$ is a diagonal matrix and the priors for
individual $\beta_{j}$ is

$\displaystyle\beta_{0}\sim N(0,\tau^{2}),$ (28) $\displaystyle\beta_{j}\sim
N(0,\tau^{2}r),\qquad i=1,\ldots,p,$ $\displaystyle\beta_{j}\sim
N(0,\tau^{2}r^{2}),\qquad i=p+1,\ldots,p+\binom{p}{2},$ $\displaystyle\vdots$
$\displaystyle\beta_{2^{p}-1}\sim N(0,\tau^{2}r^{p}),$

where $\beta_{j}$ $i=1,\ldots,p$ are main effects, $\beta_{j}$
$j=p+1,\ldots,p+\binom{p}{2}$ are 2-factor-interactions and up to the
$p$-factor-interaction $\beta_{2^{p}-1}$. The variance of $\beta_{j}$
decreases exponentially with the order of their corresponding effects by
$r\in(0,1)$, thus it incorporates the _effects hierarchy principle_ (Wu and
Hamada, 2011). Joseph (2006) showed that if $\bm{f}(\bm{x})$ contains all the
$2^{p}$ effects of all $p$ orders, $\tau^{2}\bm{R}$ can be represented
alternatively by Kronecker product as
$\tau^{2}\bm{R}=\varsigma^{2}\bigotimes_{j=1}^{p}\bm{F}_{j}(x_{j})^{-1}\bm{\Psi}_{j}(x_{j})(\bm{F}_{j}(x_{j}))^{-1}$.
The model matrix for the 2-level factor and the correlation matrix are

${\bm{F}_{j}(x_{j})}=\left(\begin{array}[]{cc}1&-1\\\
1&1\end{array}\right)\;\;\textrm{and}\;\;\bm{\Psi}_{j}(x_{j})=\left(\begin{array}[]{cc}1&\zeta\\\
\zeta&1\end{array}\right).$ (29)

To keep the two different presentations equivalent, let
$\zeta=\frac{1-r}{1+r}$ and $\tau^{2}=(\frac{1+\zeta}{2})^{p}\varsigma^{2}$.
For the mixed-level of 2- and 3-level experiments, Kang and Joseph (2009) have
extended the 2-level case to the format

$\tau^{2}\bm{R}=\varsigma^{2}\bigotimes_{j=1}^{p_{2}+p_{3,c}+p_{3,q}}\bm{F}_{j}(x_{j})^{-1}\bm{\Psi}_{j}(x_{j})(\bm{F}_{j}(x_{j})^{-1})^{\prime},$
(30)

where $p_{2}$ is the number of 2-level factors, $p_{3,c}$ is the number of
3-level qualitative (categorical) factors, and $p_{3,q}$ is the number of
3-level quantitative factors. For all the 3-level factors, the model matrix is

${\bm{F}_{j}(x_{j})}=\left(\begin{array}[]{ccc}1&-\sqrt{\frac{3}{2}}&\sqrt{\frac{1}{2}}\\\
1&0&-\sqrt{2}\\\ 1&\sqrt{\frac{3}{2}}&\sqrt{\frac{1}{2}}\end{array}\right).$
(31)

An isotropic correlation function is recommended for the 3-level qualitative
factors and a Gaussian correlation function for quantitative factors. Thus,
the correlation matrices for the 3-level qualitative and quantitative factors
are

${\bm{\Psi}_{j}(x_{j})}=\left(\begin{array}[]{ccc}1&\zeta&\zeta\\\
\zeta&1&\zeta\\\
\zeta&\zeta&1\end{array}\right)\;\;\textrm{and}\;\;\bm{\Psi}_{j}(x_{j})=\left(\begin{array}[]{ccc}1&\zeta&\zeta^{4}\\\
\zeta&1&\zeta\\\ \zeta^{4}&\zeta&1\end{array}\right),$ (32)

respectively. To keep the covariance (30) consistent with the 2-level case we
still set $\zeta=\frac{1-r}{1+r}$. To keep the variance of the intercept equal
to $\tau^{2}$ (Kang and Joseph, 2009), we set

$\tau^{2}=\varsigma^{2}\left(\frac{1+\zeta}{2}\right)^{p_{2}}\left(\frac{1+2\zeta}{3}\right)^{p_{3,c}}\left(\frac{3+4\zeta+2\zeta^{4}}{9}\right)^{p_{3,q}},$

and thus

$\bm{R}=\left(\left(\frac{1+\zeta}{2}\right)^{p_{2}}\left(\frac{1+2\zeta}{3}\right)^{p_{3,c}}\left(\frac{3+4\zeta+2\zeta^{4}}{9}\right)^{p_{3,q}}\right)^{-1}\times\bigotimes_{j=1}^{p_{2}+p_{3,c}+p_{3,q}}\bm{F}_{j}(x_{j})^{-1}\bm{\Psi}_{j}(x_{j})(\bm{F}_{j}(x_{j})^{-1})^{\prime}.$

It is straightforward to prove that $\bm{R}$ is a diagonal matrix if only
2-level and 3-level qualitative factors are involved, but not so if any
3-level quantitative factors are involved, and the first diagonal entry of
$\bm{R}$ is always 1.

To specify different prior distributions for $\bm{\beta}^{(1)}$ and
$\bm{\beta}^{(2)}$, we only need to use different values $r_{1}$ (or
$\zeta_{1}$) and $r_{2}$ (or $\zeta_{2}$) to construct the prior correlation
matrix. If the prior knowledge assumes that the two responses $Z$ and $Y$ are
independent, one can set $r_{1}=r_{2}=r$ so that the two correlation matrices
$\bm{R}_{i}$’s are the same, denoted as $\bm{R}$. Kang and Joseph (2009) has
used $r=1/3$ (equivalently $\zeta=1/2$) according to a meta-analysis of 113
data sets from published experiments (Li et al., 2006). Thus we also use
$r=1/3$ in all the examples. The readers can specify different values for
$r_{1}$ and $r_{2}$ if needed.

In computation, we construct $\bm{R}$ using the Kronecker product in (30). But
such $\bm{R}$ is for $\bm{f}(\bm{x})$ containing effects of all possible
orders. Usually, we would assume the model just contains lower-order effects.
So we just pick the rows and columns that correspond to the lower-order
effects assumed in the model as the correlation matrix.

##  5 Design Search Algorithm

In this work, we focus on the construction of optimal design based on
factorial design, which is suited for the prior distribution introduced in
Section 4.3. For optimizing the design criterion $Q(\bm{X})$ we consider two
cases. First, for fixed $\bm{\eta}$ value, we develop a point-exchange
algorithm to construct a _local optimal design_ that maximizes the criterion
$Q(\bm{X}|\bm{\eta})$. Second, we construct a _global optimal design_ based on
the prior distribution of $\bm{\eta}$. Specifically, we construct the local
optimal designs for different $\bm{\eta}$’s sampled from its prior
distribution. Then the global optimal continuous design is obtained by
accumulating the frequencies of design points selected into those local
optimal designs.

###  5.1 Local Optimal Design for Fixed $\eta$

For a fixed $\bm{\eta}$, we adapt the point-wise exchange algorithm to
maximize the criterion

$Q(\bm{X}|\bm{\eta})=\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})+\frac{1}{2}\sum_{i=1}^{n}\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}).$

The point-wise exchange algorithm is commonly used to construct $D$-optimal
designs. It was first introduced by Fedorov (1972) and then widely used in
many works (Cook and Nachtrheim, 1980; Nguyen and Miller, 1992).

The point-wise exchange algorithm finds the optimal design from a candidate
set. Here the candidate set is chosen to be the full factorial design without
replicates. For now, we develop the method for 2- and 3-level factors, but it
can be generalized to factors of more levels. Use previous notation that
$p_{2}$, $p_{3,c}$, $p_{3,q}$ as the number of 2-level, 3-level categorical,
and 3-level quantitative factors. The total number of full factorial design
points is $N=2^{p_{2}}3^{p_{3,c}+p_{3,q}}$, which can be large if the
experiment involves many factors. To make the algorithm efficient, we filter
out the candidate points that are unlikely to be the optimal design points.
Following the suggestion from Dror and Steinberg (2005), we exclude the
candidate design points whose corresponding probabilities
$\pi(\bm{x},\bm{\eta})$ is outside of $[0.15,0.85]$. This range is used
because the approximate variance of
$\log\left(\frac{\pi_{i}}{1-\pi_{i}}\right)$ is nearly constant for
$\pi_{i}\in(0.2,0.8)$ but increases rapidly if $\pi_{i}$ is outside that range
(Wu, 1985). Denote the reduced candidate set as $\bm{X}_{c}$ with size
$N^{\prime}$.

Next we construct the initial design of size $n$, such that
$\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ is nonsingular, and so should be
$\bm{F}^{\prime}\bm{W}_{i}\bm{F}$ if $\rho=0$ for $i=1,2$. If $N^{\prime}\geq
q$, we construct the initial design by reduction. Starting the initial design
as $\bm{X}_{c}$, we remove the design points one by one until there are $q$
points left. The remaining $n-q$ design points are then sampled from these $q$
initial design points with probabilities proportional to the lower bounds in
the sufficient condition in Proposition 1. For removing one design point, we
select the one having the smallest deletion function $d(\bm{x})$ defined in
(33). Shortcut formulas are developed in Appendix for updating the inverse of
the matrices $\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ and
$\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}$ for $i=1,2$ after one design
point is removed. If $N^{\prime}\leq q$, we have to restore the candidate set
back to the full factorial design and construct the initial design in the same
reduction fashion.

To simplify the notation for $d(\bm{x})$, we define
$v_{i}(\bm{x}_{1},\bm{x}_{2})=\bm{f}(\bm{x}_{1})^{\prime}\bm{M}_{i}\bm{f}(\bm{x}_{2})$
and $v_{i}(\bm{x})=\bm{f}(\bm{x})^{\prime}\bm{M}_{i}\bm{f}(\bm{x})$ for
$i=0,1,2$, where
$\bm{M}_{0}=\left(\bm{F}^{\prime}\bm{W}_{0}\bm{F}\right)^{-1}$ and
$\bm{M}_{i}=\left(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}\right)^{-1}$
for $i=1,2$. Denote $\bm{X}$ as the current design and $\bm{X}_{-i}$ the
design of $\bm{X}$ with the $i$th row removed. Then the deletion function can
be derived as

$\displaystyle d(\bm{x}_{i})$
$\displaystyle=Q(\bm{X}|\bm{\eta})-Q(\bm{X}_{-i}|\bm{\eta})$ (33)
$\displaystyle=-\left\\{\log\left[1-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))v_{0}(\bm{x}_{i})\right]+\frac{1}{2}\log\left[1-\pi(\bm{x}_{i},\bm{\eta})v_{1}(\bm{x}_{i})\right]\right.$
$\displaystyle\left.+\frac{1}{2}\log\left[1-(1-\pi(\bm{x}_{i},\bm{\eta}))v_{2}(\bm{x}_{i})\right]\right\\}.$

The smaller $d(\bm{x}_{i})$ is, the less contribution the corresponding point
makes for the overall objective $Q(\bm{X}|\bm{\eta})$.

One key of the point-wise exchange algorithm is to compute
$\Delta(\bm{x},\bm{x}_{i})=Q(\bm{X}^{*}|\bm{\eta})-Q(\bm{X}|\bm{\eta})$, the
change in the criterion after the candidate design point $\bm{x}$ replaces
$\bm{x}_{i}$ in the current design $\bm{X}$. Here $\bm{X}^{*}$ is the new
design matrix after the exchange. To compute $\Delta(\bm{x},\bm{x}_{i})$
efficiently, we can obtain the following formula.

$\displaystyle\Delta(\bm{x},\bm{x}_{i})$
$\displaystyle=Q(\bm{X}^{*}|\bm{\eta})-Q(\bm{X}|\bm{\eta})$ (34)
$\displaystyle=\log\Delta_{0}(\bm{x},\bm{x}_{i})+\frac{1}{2}\sum_{i=1}^{2}\log\Delta_{i}(\bm{x},\bm{x}_{i}),$

where $\Delta_{i}(\bm{x},\bm{x}_{i})$ for $i=0,1,2$ are derived in Appendix.
The matrices $\bm{M}_{i}$ for $i=0,1,2$ need to be updated after the exchange
of design points. Denote the updated matrices as $\bm{M}_{i}^{*}$ for the
updated design $\bm{X}^{*}$. We derive the shortcut formulas to easily compute
$\bm{M}_{i}^{*}$ for $i=0,1,2$ as shown in Appendix.

Given the initial design, we can iteratively exchange the current design
points with candidate design points to improve the objective
$Q(\bm{X}|\bm{\eta})$. The details are listed in the following Algorithm 1.

Algorithm 1 Exchange-Point Algorithm for Local $D$-Optimal Design

Step 0: Generate the candidate design set from full factorial design. Filter
out the points with probabilities $\pi(\bm{x},\bm{\eta})$ outside of
$[0.15,0.85]$.

Step 1: Generate the initial design. Based on the initial design $\bm{X}$,
update the matrices $\bm{F}$, $\bm{W}_{i}$, and $\bm{M}_{i}$ for $i=0,1,2$.
Compute the current objective value $Q(\bm{X}|\bm{\eta})$.

Step 2: Compute the deletion function $d(\bm{x}_{i})$ for each $\bm{x}_{i}$ in
$\bm{X}$. Randomly sample one design point with probability inversely
proportional to $d(\bm{x}_{i})$’s. Denote it as $\bm{x}_{i_{0}}$.

Step 3: Find $\bm{x}^{*}$ as the candidate point having the largest
$\Delta(\bm{x},\bm{x}_{i_{0}})$. If $\Delta(\bm{x}^{*},\bm{x}_{i_{0}})>0$,
exchange $\bm{x}^{*}$ with $\bm{x}_{i_{0}}$ in $\bm{X}$ and update the
objective function value to
$Q(\bm{X}|\bm{\eta})+\Delta(\bm{x}^{*},\bm{x}_{i_{0}})$.

Step 4: Repeat _Step 2_ and _3_ until the objective function has been
stabilized or the maximum number of iterations is reached.

The Algorithm 1 can return different optimal designs due to different initial
designs and the random sampling in Step 2. Thus, we run Algorithm 1 a few
times and return the design with the best optimal value. We have several
remarks regarding the algorithm. (1) The initial design generated via
reduction does not have singularity issues. (2) The updated design from point-
exchange does not have the singularity problem either, based on the way
$\bm{x}^{*}$ is selected and $\bm{M}_{i}^{*}$ for $i=0,1,2$ are computed. (3)
To avoid being trapped in a local maximum, in Step 2 we randomly sample the
design point for an exchange instead of deterministically picking the “worst”
point. (4) Different from some other point-exchange algorithms, the candidate
set here remains the same through Steps 1-4 since no points are deleted if
they are selected in the design. It enables the resultant optimal design
having replicated design points.

###  5.2 Global Optimal Design

Based on Algorithm 1 for local $D$-optimal design, we can use the following
Algorithm 2 to construct global optimal design.

Algorithm 2 Algorithm for Global $D$-Optimal Design

Step 0: If $p(\bm{\eta})$ is informative, simulate $\bm{\eta}_{j}\sim
p(\bm{\eta})$ for $j=1,\ldots,B$. Otherwise, $\bm{\eta}$ is uniformly
distributed in a rectangular high-dimensional space.

Step 1: For each $\bm{\eta}_{j}$, call _Algorithm 1_ to construct the local
optimal design $\bm{X}_{j}$.

Step 2: For each point in the candidate set, count its frequency of being
selected in the local optimal designs. The continuous optimal design is formed
by the normalized frequency as a discrete distribution.

Step 3: To obtain a discrete optimal design, sample $n$ design points from the
continuous optimal design.

In Step 1 of generating $\bm{\eta}$ uniformly, we can use uniform design (Fang
et al., 2000), maximin Latin hypercube design (Morris and Mitchell, 1995), or
other space filling design methods (Joseph and Hung, 2008; Lin et al., 2010;
Qian, 2012) to select samples $\bm{\eta}_{j}$ for $j=1,\ldots,B$. From
Algorithm 2, it is likely that the discrete design obtained in Step 3 has some
design points with $n_{i}=1$. When experimenters prefer to have replications
at every design point, they can choose a saturated design by sampling $m=q$
unique design points in Step 3. Then sample some $\bm{\eta}$ values as in Step
0. Compute the lower bounds for $n_{i}$ for every $\bm{\eta}$ sample according
to Proposition 1 and use the averaged lower bounds to set $n_{i}$. If
$\sum_{i=1}^{m}n_{i}$ exceeds $n$, the experimenters have to either increase
the experiment budget or reduce the $\kappa$ value.

##  6 Examples

In this section, we use two examples to demonstrate the proposed Bayesian
$D$-optimal design and the construction method. For both examples, we set
$r=1/3$ (equivalently $\zeta=1/2$) as explained in Section 4.3. Since there
are few existing works on experimental design for continuous and binary
responses, we compare the proposed method with three alternative designs: the
optimal designs for the quantitative-only response, the optimal design for the
binary-only response, and the naively combined design method as mentioned in
Example 1.

###  6.1 Artificial Example

In this artificial experiment, there are three 2-level factors $x_{1}\sim
x_{3}$, one 3-level categorical factor $x_{4}$, and one 3-level quantitative
factor $x_{5}$. The underlying model assumed is the complete quadratic model
and $\bm{f}(\bm{x})$ contains $q=22$ model effects including the intercept and
the following model effects.

First order effects: $\displaystyle
x_{1},x_{2},x_{3},x_{4,1},x_{4,2},x_{5,l},$ Second order effects:
$\displaystyle
x_{1}x_{2},x_{1}x_{3},x_{1}x_{4,1},x_{1}x_{4,2},x_{1}x_{5,l},x_{2}x_{3},x_{2}x_{4,1},x_{2}x_{4,2},$
$\displaystyle
x_{2}x_{5,l},x_{3}x_{4,1},x_{3}x_{4,2},x_{3}x_{5,l},x_{4,1}x_{5,l},x_{4,2}x_{5,l},{\color[rgb]{0,0,0}x_{5,\text{quad}}}.$

Here for the 3-level factors $x_{4}$ and $x_{5}$, the effects $x_{4,1}$ (1st
comparison) and $x_{5,l}$ (linear effect) have values
$\left\\{-\sqrt{\frac{3}{2}},0,\sqrt{\frac{3}{2}}\right\\}$, and $x_{4,2}$
(2nd comparison) and $x_{5,\text{quad}}$ (quadratic effect) have values
$\left\\{-\sqrt{\frac{1}{2}},\sqrt{2},\sqrt{\frac{1}{2}}\right\\}$. For the
2-level factors, the effects $x_{i}$ $i=1,2,3$ have the same values as the
design settings $\\{-1,1\\}$. We consider independent uniform distribution for
each $\eta_{i}$. Specifically, $\eta_{i}\sim\textrm{Uniform}[-1,1]$ for the
intercept and the first order effects and Uniform$[-0.5,0.5]$ for the second
order effects. The ranges of $\eta_{i}$’s satisfy the effect hierarchy
principle.

Table 1: An example of $\bm{\eta}$ value for the local $D$-optimal design. Effect | $\eta$ | Effect | $\eta$ | Effect | $\eta$ | Effect | $\eta$
---|---|---|---|---|---|---|---
Intercept | -0.0153 | $x_{1}$ | -0.6067 | $x_{2}$ | 0.7212 | $x_{1}x_{2}$ | 0.0080
$x_{3}$ | -0.1682 | $x_{1}x_{3}$ | 0.0010 | $x_{2}x_{3}$ | 0.1349 | $x_{4,1}$ | 0.0283
$x_{1}x_{4,1}$ | 0.0594 | $x_{2}x_{4,1}$ | -0.1719 | $x_{3}x_{4,1}$ | 0.1492 | $x_{4,2}$ | -0.1468
$x_{1}x_{4,2}$ | 0.0553 | $x_{2}x_{4,2}$ | -0.0634 | $x_{3}x_{4,2}$ | -0.2629 | $x_{5,l}$ | -0.0660
$x_{1}x_{5,l}$ | -0.1054 | $x_{2}x_{5,l}$ | -0.0857 | $x_{3}x_{5,l}$ | -0.0807 | $x_{4,1}x_{5,l}$ | -0.1198
$x_{4,2}x_{5,l}$ | -0.0292 | $x_{5,\text{quad}}$ | -0.1336 |

We set the experimental run size to be $n=66$. Table 1 illustrates the values
of a randomly chosen $\bm{\eta}$. Using Algorithm 2 with this $\bm{\eta}$, we
construct the proposed local $D$-optimal designs for QQ model $D_{QQ}$ for
$\rho=0$ and $\rho=0.3$, respectively. For comparison, we also generate three
alternative designs via R package _AlgDesign_ developed by Wheeler (2014).
Specifically, they are (i) the 66-run classic $D$-optimal design for linear
regression model, denoted as $D_{L}$, (ii) the 66-run local $D$-optimal design
for logistic regression model given the $\bm{\eta}$, denoted as $D_{G}$, (iii)
and the naively combined design of $44$-run local $D$-optimal design for
logistic regression model and 22-run $D$-optimal design for the linear
regression model, denoted as $D_{C}$. The details of these designs can be
found in Table A1 in Appendix.

To evaluate the performance of the proposed design in comparison with
alternative designs, we consider the efficiency between two designs (Woods et
al., 2006) as

$\displaystyle\textrm{eff}(D_{1},D_{2}|\bm{\eta})=\exp\left\\{\frac{1}{q}\left(Q(D_{1}|\bm{\eta})-Q(D_{2}|\bm{\eta})\right)\right\\}.$
(35)

Table 2 reports the efficiency of $D_{QQ}$ compared with $D_{L}$, $D_{G}$, and
$D_{C}$, respectively. The proposed QQ optimal design $D_{QQ}$ gains the best
efficiency over the three alternative designs. It appears that the combined
design $D_{C}$ has the second-best design efficiency.

Table 2: Design efficiency between the proposed local design $D_{QQ}$ and three alternative designs.  $\rho$ | $\textrm{eff}(D_{QQ},D_{L}|\bm{\eta})$ | $\textrm{eff}(D_{QQ},D_{G}|\bm{\eta})$ | $\textrm{eff}(D_{QQ},D_{C}|\bm{\eta})$
---|---|---|---
0 | 1.08 | 1.11 | 1.05
0.3 | 1.10 | 1.14 | 1.07

Next, we focus on the comparison of $D_{QQ}$ with $D_{C}$ under different
$\bm{\eta}$ values. We generate a maximin Latin hypercube design of $B=500$
runs (R package _lhs_ by Carnell (2012)) for $\bm{\eta}$ with the lower and
upper bounds specified earlier. For each of the $\bm{\eta}$ values, we
construct a local QQ optimal design $D_{QQ}$ and the combined design $D_{C}$.
Figure 2 shows the histogram of the $\textrm{eff}(D_{QQ},D_{C}|\bm{\eta})$ for
different $\bm{\eta}$ value. All $\textrm{eff}(D_{QQ},D_{C}|\bm{\eta})$ values
are larger than 1, indicating that the local QQ optimal design outperforms the
combined design.

Figure 2: Artificial example: the efficiency between each local $D_{QQ}$ and
$D_{C}$ under different $\bm{\eta}$ values: (a) $\rho=0$ (b) $\rho=0.3$.

Based on Algorithm 2, we accumulate frequencies of locally optimal designs and
obtain the global $D$-optimal designs shown in Figure 3. Denote $d_{QQ}$ and
$d_{C}$ are the proposed global optimal design for the QQ model and the global
optimal combined design, respectively. The bar plots show the normalized
frequencies for all the candidate points with the largest 22 frequencies
colored blue. From Figure 3, for $d_{QQ}$, the points in the middle have much
smaller frequencies than the other points. It is known that these points in
the middle correspond to the points with $x_{5}=0$ in Table A1. Note that
these points are only necessary for estimating the coefficient for $x_{5,q}$,
whose variances are the smallest in the prior for $\bm{\beta}^{(i)}$’s based
on the effects hierarchy principle. In contrast, such a pattern is not
observed for points with $x_{4}=0$. The reason is that $x_{4}$ is a
categorical variable and $x_{4}=-1,0,1$ are equally necessary to estimate the
effects $x_{4,1}$ and $x_{4,2}$. For $d_{C}$, the points with the largest 22
frequencies correspond to the 22-run $D$-optimal design for the linear
regression model, which is independent of $\bm{\eta}$ and remains the same
every time. The points with $x_{5}=1$ and $-1$ only have slightly higher
frequencies than the ones with $x_{5}=0$, due to the way we specify the prior
distribution of $\bm{\eta}$.

To compare the performances of the global designs, the design efficiencies in
(35) is used with $Q(d|\bm{\eta})$ adapted as

$\displaystyle Q(d|\bm{\eta})$
$\displaystyle=\log\det\left(n\sum_{i=1}^{N}d(\bm{x}_{i})\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))\bm{f}(\bm{x}_{i})f(\bm{x}_{i})^{\prime}\right)$
$\displaystyle+\frac{1}{2}\log\det\left(n\sum_{i=1}^{N}d(\bm{x}_{i})\pi(\bm{x}_{i},\bm{\eta})\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}+\rho\bm{R}\right)$
$\displaystyle+\frac{1}{2}\log\det\left(n\sum_{i=1}^{N}d(\bm{x}_{i})(1-\pi(\bm{x}_{i},\bm{\eta}))\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}+\rho\bm{R}\right)$

for a global optimal design $d$ given $\bm{\eta}$ value. Here $d(\bm{x}_{i})$
is the probability frequency for candidate design point $\bm{x}_{i}$ specified
by the design $d$ and $\sum_{i=1}^{N}d(\bm{x}_{i})=1$. For the global optimal
designs $d_{QQ}$ and $d_{C}$ obtained previously, Figure 4 shows the
histograms of the $\textrm{eff}(d_{QQ},d_{C}|\bm{\eta})$ values, where the
$\bm{\eta}$ values are generated from another 100-run maximin Latin hypercube
design. It is clear that $d_{QQ}$ is universally better than $d_{C}$, thus the
proposed design is more robust to different values of $\bm{\eta}$.

Figure 3: Artificial example: global QQ optimal designs for (a) $\rho=0$ and
(b) $\rho=0.3$ and (c) global combined design. Figure 4: Artificial example:
the efficiency between each global QQ optimal design and combined design under
different $\bm{\eta}$ values: (a) $\rho=0$ (b) $\rho=0.3$.

###  6.2 Etching Experiment

In the etching process described in Section 1, the etchant is circulating at a
certain flow rate. The wafers are rotated and swung horizontally and
vertically. Meanwhile, the air is blown in the etchant with certain pressure.
There are five factors involved in the etching process, the wafer rotation
speed ($x_{1}$), the pressure for blowing the bubbles ($x_{2}$), the
horizontal and vertical frequencies for swinging wafers ($x_{3}$, $x_{4}$),
and the flow rate of circulating the etchant ($x_{5}$). The engineers intend
to experiment to study the relationship between these factors and the two QQ
responses.

Because of the newly developed process, the historical data on similar
processes are not directly applicable to this experiment. Based on some
exploratory analysis, we set $\rho=0.5$. Both domain knowledge and data have
shown that the wafer appearance is the worst when both the rotating speed
($x_{1}$) and bubble pressure ($x_{2}$) are low. Accordingly, we set the prior
of $\bm{\eta}$ as follows. For intercept $\eta_{0}\sim\textrm{Uniform}[0,6]$.
The linear effects of rotating speed and bubble pressure follow
$\textrm{Uniform}[1,5]$. The other linear effects follow
$\textrm{Uniform}[-1,1]$ and the second order interactions and quadratic
effects $\textrm{Uniform}[-0.3,0.3]$. The experimental run size is set to be
$n=21\times 6=126$, 6 times the number of effects $q=21$.

We generate a maximin Latin hypercube design of $B=500$ runs for $\bm{\eta}$
values. For each $\bm{\eta}$ value, we obtain the local optimal designs
$D_{QQ}$ and $D_{C}$. Here the local combined design $D_{C}$ has $2/3$ of the
runs generated from the local $D$-optimal design for logistic regression and
$1/3$ of the runs from the $D$-optimal design for linear regression. The
efficiency between each pair of local designs $D_{QQ}$ and $D_{C}$ is reported
in Figure 6(a). We can see that almost every local design $D_{QQ}$ is better
than $D_{C}$. Moreover, we obtain the global optimal designs $d_{QQ}$ and
$d_{C}$ by accumulating the frequencies of the local designs. To compare
$d_{QQ}$ and $d_{C}$, we generate another 100-run maximin Latin hypercube
design for $\bm{\eta}$ values and compute the efficiencies between $d_{QQ}$
and $d_{C}$ under different $\bm{\eta}$ values, which are shown in Figure 6
(b). Clearly, $d_{QQ}$ is universally better and more robust to $\bm{\eta}$
than $d_{C}$.

Fractional factorial design (Wu and Hamada, 2011) is another commonly used
design in practice. We compare the proposed design with a $3^{5-2}$ minimum
aberration (MA) fractional factorial design by defining the contrast subgroup
as $I=ABD^{2}=AB^{2}CE^{2}=AC^{2}DE=BCDE^{2}$. Each design point is replicated
5 times and the overall run is $3^{5-2}\times 5=135$. Figure 6(c) shows the
histogram of the efficiencies between $d_{QQ}$ and the MA design, and the
proposed global optimal design is still superior.

Figure 5: Etching experiment: (a) the global Bayesian QQ $D$-optimal design
for $\rho=0.5$ and (b) the global combined design. Figure 6: Etching
experiment: histograms of the efficiencies (a) efficiencies between local
designs $D_{QQ}$ and $D_{C}$ for different $\bm{\eta}$’s; (b) efficiencies
between global designs $d_{QQ}$ and $d_{C}$; (c) efficiencies between global
design $d_{QQ}$ and the $3^{5-2}$ MA fractional factorial design.

##  7 Discussion

In this paper, we propose the Bayesian $D$-optimal design criterion for QQ
responses. The adoption of the Bayesian approach allows us to consider both
the non-informative priors as the frequentist approach and informative priors
when domain knowledge or historical data are available. A new point-exchange
algorithm is developed for efficiently constructing the proposed designs. This
algorithm can also be used to construct non-factorial designs when the
candidate set is not a full factorial design. Moreover, the proposed method
can be directly generalized for the sequential design with the QQ response. In
the following, we discuss some other scenarios for the proposed method that
are not considered in detail previously.

Non-conjugate Prior for $\bm{\eta}$

Other than the conjugate prior $p(\bm{\eta})$, we can also use a non-conjugate
prior distribution $\bm{\eta}\sim N(\bm{0},\tau_{0}^{2}\bm{R}_{0})$. In this
situation, one can consider the normal approximation in the posterior
distribution $p(\bm{\eta}|\bm{z})$. Then the design criterion for the binary
response becomes (Chaloner and Verdinelli, 1995)
$\mathbb{E}_{\bm{z},\bm{\eta}}\\{\log(p(\bm{\eta}|\bm{z})\\}\approx\mathbb{E}_{\bm{\eta}}\left\\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F}+\rho_{0}\bm{R}_{0}^{-1})\right\\}$.
The overall design criterion $Q(\bm{X})$ can be updated as

$Q(\bm{X})=\sum_{i=0}^{2}\mathbb{E}_{\bm{\eta}}\left\\{\log\det\left(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\tau_{i}\bm{R}_{i}^{-1}\right)\right\\},$

where $\rho_{i}=\sigma^{2}/\tau_{i}$, $\tau_{i}$ and $\bm{R}_{i}$ are the
prior variance and correlation of $\bm{\eta}$, $\bm{\beta}^{(1)}$, and
$\bm{\beta}^{(2)}$ respectively. The proposed design construction algorithm
can still be applied with minor modifications.

Multiple QQ Responses

In this paper, we focus on optimal designs for one quantitative response and
one qualitative response. The proposed method can also be generalized to
accommodate the QQ models with multiple quantitative responses
$Y_{1},\ldots,Y_{l}$ and binary responses $Z_{1},\ldots,Z_{k}$. For example, a
multi-level qualitative response can be transformed into a set of dummy binary
responses. One idea is to generalize the QQ models in (3) and (4) for both
$l\geq 1$ and $k\geq 1$. For $l\geq 1$ and $k=1$, we generalize the QQ model
by introducing the correlation matrix between $Y_{1},\ldots,Y_{l}$ as in the
standard multi-response regression (Breiman and Friedman, 1997). Then the
corresponding optimal design can be established by studying its likelihood
function. For $l\geq 1$ and $k>1$ with multiple binary responses, considering
all $2^{k}$ conditional models for
$(Y_{1},\ldots,Y_{l}|Z_{1}=1,\ldots,Z_{k}=1)$,…,
$(Y_{1},\ldots,Y_{l}|Z_{1}=0,\ldots,Z_{k}=0)$ only works for a small $k$.
Moreover, the construction algorithm can be more complicated as it needs to
involve the multi-logit model (McCullagh, 1980) for modeling the multiple
binary responses. When $k$ is relatively large, we are going to pursue an
alternative QQ model and develop its corresponding optimal design method as a
future research topic.

Continuous Design

The point-exchange algorithm is to construct the exact discrete optimal
designs, which are different from the theoretical continuous optimal designs.
As described in Sections 5 and 6, the way of generating the frequency as the
local optimal design is heuristic. The rigorous definition of the local
continuous $D$-optimal design criterion is the probability measure $\xi$ on
the design space $\Omega$ that maximizes

$\displaystyle
Q(\bm{X}|\bm{\eta})=\log\det\left(\int\pi(\bm{x})(1-\pi(\bm{x}))f(\bm{x})f(\bm{x})^{\prime}d\xi(\bm{x})\right)+$
$\displaystyle\log\det\left(\int(\pi\left(\bm{x})f(\bm{x})f(\bm{x})^{\prime}+\rho\bm{R}_{1}^{-1}\right)d\xi(\bm{x})\right)+\log\det\left(\int\left((1-\pi(\bm{x}))f(\bm{x})f(\bm{x})^{\prime}+\rho\bm{R}_{2}^{-1}\right)d\xi(\bm{x})\right).$

Yang et al. (2013) developed a method to obtain the optimal $\xi$ for the
nonlinear models. It will be interesting to extend their framework and develop
the method to obtain the optimal $\xi$ for QQ models.

Different QQ Models

The proposed design is not restricted to the logit model for the binary
response. For example, if the probit model is used, the Bayesian $D$-optimal
design criterion can be directly obtained by replacing the logit
transformation with the probit transformation in both $p(\bm{z}|\bm{\eta})$
and $p(\bm{\eta})$. The design criterion can be derived similarly with minor
modifications. The criterion formula remains the same with the following
different diagonal matrices,

$\displaystyle\bm{W}_{0}$
$\displaystyle=\mbox{diag}\left\\{\frac{\phi^{2}(\bm{f}(\bm{x}_{1})^{\prime}\bm{\eta})}{\Phi(\bm{f}(\bm{x}_{1})^{\prime}\bm{\eta})\left(1-\Phi(\bm{f}(\bm{x}_{1})^{\prime}\bm{\eta})\right)},\ldots,\frac{\phi^{2}(\bm{f}(\bm{x}_{n})^{\prime}\bm{\eta})}{\Phi(\bm{f}(\bm{x}_{n})^{\prime}\bm{\eta})\left(1-\Phi(\bm{f}(\bm{x}_{n})^{\prime}\bm{\eta})\right)}\right\\},$
$\displaystyle\bm{W}_{1}$
$\displaystyle=\mbox{diag}\left\\{\Phi\left(\bm{f}(\bm{x}_{1})^{\prime}\bm{\eta}\right),\ldots,\Phi\left(\bm{f}(\bm{x}_{n})^{\prime}\bm{\eta}\right)\right\\},\quad\bm{W}_{2}=\bm{I}-\bm{W}_{1},$

where $\Phi$ and $\phi$ are CDF and PDF of the standard normal distribution.

It is worth pointing out that the design criterion in the work is based on the
QQ model constructed by the joint model of $Y|Z$ in (3) and $Z$ in (4). Kang
et al. (2021) created a new QQ model based on $Z|U$ where $U$ is a latent
continuous variable that is assumed to be correlated with the observed
continuous response variable $Y$. Besides the conditional model structures,
other model structures such as mixed graphical models (Yang et al., 2014) can
also be used as long as the $D$-optimality can be derived.

Acknowledgement

The authors were partly supported by U.S. National Science Foundation for this
research project. Dr. Lulu Kang was supported by grants CMMI-1435902,
DMS-1916467, and DMS-2153029, Dr. Xinwei Deng by CMMI-1233571 and
CMMI-1435996, and Dr. Ran Jin by CMMI-1435996.

## References

  * (1)
  * Breiman and Friedman (1997) Breiman, L., and Friedman, J. H. (1997), “Predicting multivariate responses in multiple linear regression,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59(1), 3–54.
  * Carnell (2012) Carnell, R. (2012), lhs: Latin Hypercube Samples.  R package version 0.10.
http://CRAN.R-project.org/package=lhs

  * Chaloner and Verdinelli (1995) Chaloner, K., and Verdinelli, I. (1995), “Bayesian experimental design: A review,” Statistical Science, 10(3), 273–304.
  * Chen and Ibrahim (2003) Chen, M. H., and Ibrahim, J. G. (2003), “Conjugate priors for generalized linear models,” Statistica Sinica, 13(2), 461–476.
  * Cook and Nachtrheim (1980) Cook, R. D., and Nachtrheim, C. J. (1980), “A comparison of algorithms for constructing exact D-optimal designs,” Technometrics, 22(3), 315–324.
  * Denman et al. (2011) Denman, N., McGree, J., Eccleston, J., and Duffull, S. B. (2011), “Design of experiments for bivariate binary responses modelled by Copula functions,” Computational Statistics & Data Analysis, 55(4), 1509–1520.
  * Diaconis and Ylvisaker (1979) Diaconis, P., and Ylvisaker, D. (1979), “Conjugate priors for exponential families,” The Annals of Statistics, 7(2), 269–281.
  * Draper and Hunter (1966) Draper, N. R., and Hunter, W. G. (1966), “Design of experiments for parameter estimation in multiresponse situations,” Biometrika, 53(3-4), 525–533.
  * Dror and Steinberg (2005) Dror, H. A., and Steinberg, D. M. (2005), Approximate Local D-optimal Experimental Design for Binary Response, Technical Report RP SOR-0501, Tel Aviv University.
  * Dror and Steinberg (2006) Dror, H. A., and Steinberg, D. M. (2006), “Robust experimental design for multivariate generalized linear models,” Technometrics, 48(4), 520–529.
  * Fang et al. (2000) Fang, K.-T., Lin, D. K., Winker, P., and Zhang, Y. (2000), “Uniform design: theory and application,” Technometrics, 42(3), 237–248.
  * Fedorov (1972) Fedorov, V. V. (1972), Theory of Optimal Experiments, Vol. 12, New York: Academic Press.
  * Hainy et al. (2013) Hainy, M., Müller, W. G., and Wynn, H. P. (2013), “Approximate Bayesian computation design (ABCD), an introduction,” in mODa 10–Advances in Model-Oriented Design and Analysis, Switzerland: Springer, pp. 135–143.
  * Joseph (2006) Joseph, V. R. (2006), “A Bayesian approach to the design and analysis of fractionated experiments,” Technometrics, 48(2), 219–229.
  * Joseph and Hung (2008) Joseph, V. R., and Hung, Y. (2008), “Orthogonal-maximin Latin hypercube designs,” Statistica Sinica, 18(1), 171–186.
  * Kang and Huang (2019) Kang, L., and Huang, X. (2019), “Bayesian A-Optimal Design of Experiment with Quantitative and Qualitative Responses,” Journal of Statistical Theory and Practice, 13(4), 64.
https://doi.org/10.1007/s42519-019-0063-6

  * Kang and Joseph (2009) Kang, L., and Joseph, V. R. (2009), “Bayesian optimal single arrays for robust parameter design,” Technometrics, 51(3), 250–261.
  * Kang et al. (2018) Kang, L., Kang, X., Deng, X., and Jin, R. (2018), “A Bayesian hierarchical model for quantitative and qualitative responses,” Journal of Quality Technology, 50(3), 290–308.
https://doi.org/10.1080/00224065.2018.1489042

  * Kang et al. (2022) Kang, X., Kang, L., Chen, W., and Deng, X. (2022), “A generative approach to modeling data with quantitative and qualitative responses,” Journal of Multivariate Analysis, 190, 104952.
https://www.sciencedirect.com/science/article/pii/S0047259X22000045

  * Kang et al. (2021) Kang, X., Ranganathan, S., Kang, L., Gohlke, J., and Deng, X. (2021), “Bayesian auxiliary variable model for birth records data with qualitative and quantitative responses,” Journal of Statistical Computation and Simulation, 91(16), 3283–3303.
https://doi.org/10.1080/00949655.2021.1926459

  * Khuri et al. (2006) Khuri, A. I., Mukherjee, B., Sinha, B. K., and Ghosh, M. (2006), “Design issues for generalized linear models: A review,” Statistical Science, 21(3), 376–399.
  * Li et al. (2006) Li, X., Sudarsanam, N., and Frey, D. D. (2006), “Regularities in data from factorial experiments,” Complexity, 11(5), 32–45.
  * Lin et al. (2010) Lin, C. D., Bingham, D., Sitter, R. R., and Tang, B. (2010), “A new and flexible method for constructing designs for computer experiments,” The Annals of Statistics, 38, 1460–1477.
  * McCullagh (1980) McCullagh, P. (1980), “Regression models for ordinal data,” Journal of the Royal Statistical Society. Series B (Methodological), 42, 109–142.
  * Morris and Mitchell (1995) Morris, M. D., and Mitchell, T. J. (1995), “Exploratory designs for computational experiments,” Journal of Statistical Planning and Inference, 43(3), 381–402.
  * Nguyen and Miller (1992) Nguyen, N.-K., and Miller, A. J. (1992), “A review of some exchange algorithms for constructing discrete D-optimal designs,” Computational Statistics and Data Analysis, 14(4), 489–498.
  * Qian (2012) Qian, P. Z. (2012), “Sliced Latin hypercube designs,” Journal of the American Statistical Association, 107(497), 393–399.
  * Russell et al. (2009) Russell, K. G., Woods, D. C., Lewis, S., and Eccleston, J. (2009), “D-optimal designs for Poisson regression models,” Statistica Sinica, 19, 721–730.
  * Sitter and Torsney (1995) Sitter, R. R., and Torsney, B. (1995), “Optimal designs for binary response experiments with two design variables,” Statistica Sinica, 5, 405–419.
  * Wheeler (2014) Wheeler, B. (2014), AlgDesign: Algorithmic Experimental Design.  R package version 1.1-7.2.
http://CRAN.R-project.org/package=AlgDesign

  * Woods (2010) Woods, D. (2010), “Robust designs for binary data: applications of simulated annealing,” Journal of Statistical Computation and Simulation, 80(1), 29–41.
  * Woods and Van de Ven (2012) Woods, D. C., and Van de Ven, P. (2012), “Blocked designs for experiments with correlated non-normal response,” Technometrics, .
  * Woods et al. (2006) Woods, D., Lewis, S., Eccleston, J., and Russell, K. (2006), “Designs for generalized linear models with several variables and model uncertainty,” Technometrics, 48(2), 284–292.
  * Wu (1985) Wu, C. F. J. (1985), “Efficient Sequential Designs With Binary Data,” Journal of the American Statistical Association, 80(392), 974–984.
  * Wu and Hamada (2011) Wu, C. J., and Hamada, M. S. (2011), Experiments: planning, analysis, and optimization, Vol. 552, New Jersey: John Wiley & Sons.
  * Yang et al. (2014) Yang, E., Baker, Y., Ravikumar, P., Allen, G., and Liu, Z. (2014), Mixed Graphical Models via Exponential Families,, in Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, pp. 1042–1050.
  * Yang et al. (2013) Yang, M., Biedermann, S., and Tang, E. (2013), “On optimal designs for nonlinear models: a general and efficient algorithm,” Journal of the American Statistical Association, 108(504), 1411–1420.

Appendix: Proofs, Derivations, and Extra Table

##  A1 Proofs and Derivations

Proof of (11)

###### Proof.

$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})=\int
p(\bm{y},\bm{z}|\mu_{0},\sigma^{2})$
$\displaystyle\int\log\left(p(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta}|\bm{y},\bm{z},\mu_{0},\sigma^{2})\right)$
$\displaystyle
p(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta}|\bm{y},\bm{z},\mu_{0},\sigma^{2})d\bm{\beta}^{(1)}d\bm{\beta}^{(2)}d\bm{\eta}d\bm{y}d\bm{z}$
$\displaystyle=\int p(\bm{y},\bm{z}|\mu_{0},\sigma^{2})$
$\displaystyle\left\\{\int\log\left(p(\bm{\eta}|\bm{z})\right)p(\bm{\eta}|\bm{z})d\bm{\eta}\right.$
$\displaystyle+\left.\sum_{i=1}^{2}\int\log\left(p(\bm{\beta}^{(i)}|\bm{y},\bm{z},\mu_{0},\sigma^{2})\right)p(\bm{\beta}^{(i)}|\bm{y},\bm{z},\mu_{0},\sigma^{2})d\bm{\beta}^{(i)}\right\\}d\bm{y}d\bm{z}$
$\displaystyle=\int p(\bm{y},\bm{z}|\mu_{0},\sigma^{2})$
$\displaystyle\left\\{\int\log\left(p(\bm{\eta}|\bm{z})\right)p(\bm{\eta}|\bm{z})d\bm{\eta}\right.$
$\displaystyle-\left.\frac{1}{2}\sum_{i=1}^{2}\log\det\\{\sigma^{-2}(\bm{F}^{\prime}\bm{V}_{i}\bm{F})^{-1}\\}-n\log(2\pi)-n\right\\}d\bm{y}d\bm{z}$
$\displaystyle=\int\log(p(\bm{\eta}|\bm{z}))p(\bm{\eta}|\bm{z})$
$\displaystyle p(\bm{z})d\bm{\eta}d\bm{z}$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{2}\int\log\det\\{\sigma^{-2}(\bm{F}^{\prime}\bm{V}_{i}\bm{F})^{-1}\\}p(\bm{z})d\bm{z}+\textrm{constants}$
$\displaystyle=\int\log(p(\bm{\eta}|\bm{z}))p(\bm{z},\bm{\eta})$
$\displaystyle d\bm{\eta}d\bm{z}$
$\displaystyle+\frac{1}{2}\sum_{i=1}^{2}\int\log\det\\{(\bm{F}^{\prime}\bm{V}_{i}\bm{F})\\}p(\bm{z}|\bm{\eta})p(\bm{\eta})d\bm{z}d\bm{\eta}+\textrm{constants}.$

Write the integration into the form of expectation,

$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})$
$\displaystyle=\mathbb{E}_{\bm{z},\bm{\eta}}\left\\{\log(p(\bm{\eta}|\bm{z}))\right\\}$
$\displaystyle+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\\{\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F})\right\\}+\textrm{constants}.$

∎

Proof of Theorem 1

###### Proof.

To show (14), we just need to show that
$\mathbb{E}_{\bm{z}|\bm{\eta}}\left(\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F})\right)\leq\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F}^{\prime})$
for $i=1,2$.

First we need to show that
$\log\det(\bm{F}^{\prime}\bm{A}\bm{F})=\log\det(\sum_{i=1}^{n}a_{i}\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime})$
with $\bm{A}=\mbox{diag}\\{a_{1},\ldots,a_{n}\\}$ is a concave function of
$\bm{a}=(a_{1},\ldots,a_{n})^{\prime}$ for $\bm{a}\in[0,1]^{n}$. Denote
$\bm{B}=\mbox{diag}\\{b_{1},\ldots,b_{n}\\}$ for
$\bm{b}=(b_{1},\ldots,b_{n})^{\prime}\in[0,1]^{n}$. We assume that
$\bm{F}^{\prime}\bm{A}\bm{F}$ is nonsingular, thus
$\bm{F}^{\prime}\bm{A}\bm{F}$ is positive definite. For any scalar value of
$t$, define function $g(t)$

$\displaystyle g(t)$
$\displaystyle=\log\det\left(\bm{F}^{\prime}\bm{A}\bm{F}+\bm{F}^{\prime}(t\bm{B})\bm{F}\right)$
$\displaystyle=\log\det(\bm{F}^{\prime}\bm{A}\bm{F})+\log\det(\bm{I}_{q}+t(\bm{F}^{\prime}\bm{A}\bm{F})^{-1/2}(\bm{F}^{\prime}\bm{B}\bm{F})(\bm{F}^{\prime}\bm{A}\bm{F})^{-1/2})$
$\displaystyle=\log\det(\bm{F}^{\prime}\bm{A}\bm{F})+\sum_{i=1}^{q}\log(1+t\lambda_{i}),$

where $\lambda_{i}$’s are the eigenvalues of the positive definite matrix
$(\bm{F}^{\prime}\bm{A}\bm{F})^{-1/2}(\bm{F}^{\prime}\bm{B}\bm{F})(\bm{F}^{\prime}\bm{A}\bm{F})^{-1/2}$.
Thus $g(t)$ is a concave function in $t$ for any choice of $\bm{a}$, which is
the sufficient and necessary condition that
$\log\det(\bm{F}^{\prime}\bm{A}\bm{F})$ is a concave function of $\bm{a}$.
According to Jensen’s inequality, if $\pi(\bm{x}_{i},\bm{\eta})\in(0,1)$ for
$i=1,\ldots,n$, then

$\displaystyle\mathbb{E}_{\bm{z}|\bm{\eta}}\left(\log\det(\bm{F}^{\prime}\bm{V}_{1}\bm{F})\right)$
$\displaystyle\leq\log\det\left(\sum_{i=1}^{n}\mathbb{E}(Z_{i}|\bm{\eta})\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\right)$
$\displaystyle=\log\det\left(\sum_{i=1}^{n}\pi(\bm{x}_{i},\bm{\eta})\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\right),$
$\displaystyle\mathbb{E}_{\bm{z}|\bm{\eta}}\left(\log\det(\bm{F}^{\prime}\bm{V}_{2}\bm{F})\right)$
$\displaystyle\leq\log\det\left(\sum_{i=1}^{n}\mathbb{E}((1-Z_{i})|\bm{\eta})\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\right)$
$\displaystyle=\log\det\left(\sum_{i=1}^{n}(1-\pi(\bm{x}_{i},\bm{\eta}))\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\right).$

So we have proved (14). ∎

Proof of Proposition 1

###### Proof.

$m=\sum_{i=1}^{m}I\left(0<\sum_{j=1}^{n_{i}}Z_{ij}<n_{i}\right)+\sum_{i=1}^{m}I\left(\sum_{j=1}^{n_{i}}Z_{ij}=0\right)+\sum_{i=1}^{m}I\left(\sum_{j=1}^{n_{i}}Z_{ij}=n_{i}\right)\geq
q$

If $m=q$, (15) is equivalent to
$I\left(0<\sum_{j=1}^{n_{i}}Z_{ij}<n_{i}\right)=1$ for $i=1,2,\ldots,m$. A
sufficient condition for any
$\Pr\left(0<\sum_{j=1}^{n_{i}}Z_{ij}<n_{i}\right)\geq\kappa$ can be derived as
follows.

$\displaystyle\Pr\left(0<\sum_{j=1}^{n_{i}}Z_{ij}<n_{i}\right)$
$\displaystyle=1-\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}-(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}\geq\kappa$
$\displaystyle\Longleftrightarrow\quad$
$\displaystyle\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}+(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}\leq
1-\kappa.$

Next develop an upper bound for
$\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}+(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}$. If
$\pi(\bm{x}_{i},\bm{\eta})\geq 1/2$, denote
$a=(1-\pi(\bm{x}_{i},\bm{\eta}))/\pi(\bm{x}_{i},\bm{\eta})$ and $a\leq 1$.
Then

$\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}+(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}=\frac{1+a^{n_{i}}}{(1+a)^{n_{i}}}\leq\frac{1+a}{(1+a)^{n_{i}}}.$

If

$\frac{1+a}{(1+a)^{n_{i}}}\leq 1-\kappa\quad\Longleftrightarrow\quad n_{i}\geq
1+\frac{\log(1-\kappa)}{\log\pi(\bm{x}_{i},\bm{\eta})},$

then
$\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}+(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}\leq
1-\kappa$. If $\pi(\bm{x}_{i},\bm{\eta})<1/2$, denote
$a=\pi(\bm{x}_{i},\bm{\eta})/(1-\pi(\bm{x}_{i},\bm{\eta}))$ and $a<1$. Then

$\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}+(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}=\frac{1+a^{n_{i}}}{(1+a)^{n_{i}}}<\frac{1+a}{(1+a)^{n_{i}}}.$

The sufficient condition becomes

$\frac{1+a}{(1+a)^{n_{i}}}\leq 1-\kappa\quad\Longleftrightarrow\quad n_{i}\geq
1+\frac{\log(1-\kappa)}{\log(1-\pi(\bm{x}_{i},\bm{\eta}))}.$

Combining the two cases, the sufficient condition on $n_{i}$ for
$i=1,\ldots,m$ is (16). It is known that

$2\left(\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta})\right)^{n_{i}/2}\leq\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}+(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}\leq
1-\kappa.$

The necessary condition is

$2\left(\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta})\right)^{n_{i}/2}\leq
1-\kappa\quad\Longleftrightarrow\quad
n_{i}\geq\frac{2\log\left(\frac{1-\kappa}{2}\right)}{\log\pi(\bm{x}_{i},\bm{\eta})+\log(1-\pi(\bm{x}_{i},\bm{\eta}))}.$

∎

Proof of Proposition 2

###### Proof.

If $m>q$, (15) is equivalent to

$\sum_{i=1}^{m}I\left(\sum_{j=1}^{n_{i}}Z_{ij}=0\right)\leq m-q,\textrm{ and
}\sum_{i=1}^{m}I\left(\sum_{j=1}^{n_{i}}Z_{ij}=n_{i}\right)\leq m-q.$

For the two inequalities to hold with large probability,

$\displaystyle\sum_{i=1}^{m}\mathbb{E}\left\\{I\left(\sum_{j=1}^{n_{i}}Z_{ij}=0\right)\right\\}=\sum_{i=1}^{m}(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}\leq
m-q,$ (36) and
$\displaystyle\sum_{i=1}^{m}\mathbb{E}\left\\{I\left(\sum_{j=1}^{n_{i}}Z_{ij}=n_{i}\right)\right\\}=\sum_{i=1}^{m}\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}\leq
m-q.$ (37)

It is known that

$m\left(\prod_{i=1}^{m}(1-\pi_{\max})^{n_{i}}\right)^{1/m}\leq\sum_{i=1}^{m}(1-\pi_{\max})^{n_{i}}\leq\sum_{i=1}^{m}(1-\pi(\bm{x}_{i},\bm{\eta}))^{n_{i}}\leq\sum_{i=1}^{m}(1-\pi_{\min})^{n_{i}}\leq
m(1-\pi_{\min})^{n_{0}}.$

Thus one sufficient condition for (36) is

$m(1-\pi_{\min})^{n_{0}}\leq m-q\quad\Longleftrightarrow\quad
n_{0}\geq\frac{\log(1-q/m)}{\log(1-\pi_{\min})}.$

One necessary condition for (36) is

$m\left(\prod_{i=1}^{m}(1-\pi_{\max})^{n_{i}}\right)^{1/m}\leq
m-q\quad\Longleftrightarrow\sum_{i=1}^{m}n_{i}\geq
m\frac{\log(1-q/m)}{\log(1-\pi_{\max})}.$

Similarly,

$m\left(\prod_{i=1}^{m}\pi_{\min}^{n_{i}}\right)^{1/m}\leq\sum_{i=1}^{m}\pi_{\min}^{n_{i}}\leq\sum_{i=1}^{m}\pi(\bm{x}_{i},\bm{\eta})^{n_{i}}\leq\sum_{i=1}^{m}\pi_{\max}^{n_{i}}\leq
m\cdot\pi_{\max}^{n_{0}}.$

Thus one sufficient condition for (37) is

$m\cdot\pi_{\max}^{n_{0}}\leq m-q\quad\Longleftrightarrow\quad
n_{0}\geq\frac{\log(1-q/m)}{\log\pi_{\max}}.$

One necessary condition for (37) is

$m\left(\prod_{i=1}^{m}\pi_{\min}^{n_{i}}\right)^{1/m}\leq
m-q\quad\Longleftrightarrow\sum_{i=1}^{m}n_{i}\geq
m\frac{\log(1-q/m)}{\log\pi_{\min}}.$

Thus the sufficient condition for (36) and (37) derived here is

$n_{0}\geq\max\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\min})},\frac{\log(1-q/m)}{\log\pi_{\max}}\right\\},$

or equivalently,

$\sum_{i=1}^{m}n_{i}\geq m\cdot n_{0}\geq
m\cdot\max\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\min})},\frac{\log(1-q/m)}{\log\pi_{\max}}\right\\}.$

The necessary condition for (36) and (37) derived here is

$\sum_{i=1}^{n}n_{i}\geq
m\cdot\max\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\max})},\frac{\log(1-q/m)}{\log\pi_{\min}}\right\\}.$

It is easy to see that this lower bound in the necessary condition is smaller
than the one in the sufficient condition,

$\sum_{i=1}^{m}n_{i}\geq m\cdot n_{0}\geq
m\cdot\max\left\\{1,\frac{\log(1-q/m)}{\log(1-\pi_{\min})},\frac{\log(1-q/m)}{\log\pi_{\max}}\right\\}.$

∎

Proof of Theorem 2

###### Proof.

We only need to show that
$\log\det(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1})=\log\det(\sum_{i=1}^{n}a_{i}\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}+\rho\bm{R}_{i}^{-1})$
with $\bm{A}=\mbox{diag}\\{a_{1},\ldots,a_{n}\\}$ is a concave function of
$\bm{a}=(a_{1},\ldots,a_{n})^{\prime}$ for $\bm{a}\in[0,1]^{n}$. Denote
$\bm{B}=\mbox{diag}\\{b_{1},\ldots,b_{n}\\}$ for
$\bm{b}=(b_{1},\ldots,b_{n})^{\prime}\in[0,1]^{n}$ and
$\bm{K}=\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1}$. Define $g(t)$ as
follows.

$\displaystyle g(t)$
$\displaystyle=\log\det\left(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1}+\bm{F}^{\prime}(t\bm{B})\bm{F}+\rho\bm{R}_{i}^{-1}\right)$
$\displaystyle=\log\det(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1})$
$\displaystyle+\log\det(\bm{I}_{q}+(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1})^{-1/2}(\bm{F}^{\prime}t\bm{B}\bm{F}+\rho\bm{R}_{i}^{-1})(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1})^{-1/2})$
$\displaystyle=\log\det(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1})$
$\displaystyle+\log\det(\bm{I}_{q}+t\bm{K}^{-1/2}\bm{F}^{\prime}\bm{B}\bm{F}\bm{K}^{-1/2}+\rho\bm{K}^{-1/2}\bm{R}_{i}^{-1}\bm{K}^{-1/2})$
$\displaystyle=\log\det(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1})+\log\det(\bm{I}_{q}+\rho\bm{K}^{-1/2}\bm{R}_{i}^{-1}\bm{K}^{-1/2})+\sum_{i=1}^{q}\log(1+t\lambda_{i}),$

where $\lambda_{i}$’s are the eigenvalues of

$(\bm{I}_{q}+\rho\bm{K}^{-1/2}\bm{R}_{i}^{-1}\bm{K}^{-1/2})^{-1/2}\bm{K}^{-1/2}\bm{F}^{\prime}\bm{B}\bm{F}\bm{K}^{-1/2}(\bm{I}_{q}+\rho\bm{K}^{-1/2}\bm{R}_{i}^{-1}\bm{K}^{-1/2})^{-1/2}.$

Therefore, $\lambda_{i}\geq 0$ for $i=1,\ldots,q$. Thus $g(t)$ is a concave
function in $t$ for any choice of $\bm{a}\in[0,1]^{n}$, which is the
sufficient and necessary condition for
$\log\det(\bm{F}^{\prime}\bm{A}\bm{F}+\rho\bm{R}_{i}^{-1})$ to be a concave
function of $\bm{a}$. ∎

Proof of the deletion function (33)

###### Proof.

Denote $\bm{F}_{-i}$ as the model matrix for $\bm{X}_{-i}$, which is the
design matrix without the $i$th design point, and $\bm{W}_{0,-i}$,
$\bm{W}_{1,-i}$, and $\bm{W}_{2,-i}$ the weight matrices accordingly.
According to the properties of matrix determinants, we can show the following.

$\displaystyle\det\left(\bm{F}^{\prime}_{-i}\bm{W}_{0,-i}\bm{F}_{-i}\right)$
$\displaystyle=\det\left(\sum_{j\neq
i}\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))\bm{f}(\bm{x}_{j})\bm{f}(\bm{x}_{j})^{\prime}\right)$
$\displaystyle=\det\left(\bm{F}^{\prime}\bm{W}_{0}\bm{F}-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\right)$
$\displaystyle=\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\left[1-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))\bm{f}(\bm{x}_{i})^{\prime}(\bm{F}^{\prime}\bm{W}_{0}\bm{F})^{-1}\bm{f}(\bm{x}_{i})\right]$
$\displaystyle=\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\left[1-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))v_{0}(\bm{x}_{i})\right].$
$\displaystyle\det\left(\bm{F}^{\prime}_{-i}\bm{W}_{1,-i}\bm{F}_{-i}+\rho\bm{R}^{-1}\right)$
$\displaystyle=$ $\displaystyle\det\left(\sum_{j\neq
i}\pi(\bm{x}_{i},\bm{\eta})\bm{f}(\bm{x}_{j})\bm{f}(\bm{x}_{j})^{\prime}+\rho\bm{R}^{-1}\right)$
$\displaystyle=$
$\displaystyle\det\left(\bm{F}^{\prime}\bm{W}_{1}\bm{F}+\rho\bm{R}^{-1}-\pi(\bm{x}_{i},\bm{\eta})\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\right)$
$\displaystyle=$
$\displaystyle\det\left(\bm{F}^{\prime}\bm{W}_{1}\bm{F}+\rho\bm{R}^{-1}\right)\left[1-\pi(\bm{x}_{i},\bm{\eta})\bm{f}(\bm{x}_{i})^{\prime}(\bm{F}^{\prime}\bm{W}_{1}\bm{F}+\rho\bm{R}^{-1})^{-1}\bm{f}(\bm{x}_{i})\right]$
$\displaystyle=$
$\displaystyle\det\left(\bm{F}^{\prime}\bm{W}_{1}\bm{F}+\rho\bm{R}^{-1}\right)\left[1-\pi(\bm{x}_{i},\bm{\eta})v_{1}(\bm{x}_{i})\right].$

Similarly,

$\det\left(\bm{F}^{\prime}_{-i}\bm{W}_{2,-i}\bm{F}_{-i}+\rho\bm{R}_{2}^{-1}\right)=\det\left(\bm{F}^{\prime}\bm{W}_{2}\bm{F}+\rho\bm{R}^{-1}\right)\left[1-(1-\pi(\bm{x}_{i},\bm{\eta}))v_{2}(\bm{x}_{i})\right].$

The deletion function (33) is then obtained. ∎

Shortcut formulas for $\bm{M}_{i,-j}$.

###### Proof.

Denote $\bm{M}_{i,-j}$ for $i=0,1,2$ as the updated $\bm{M}_{i}$ when the
$j$th design point is removed. The following shortcut formulas are used in
constructing the initial design.

$\bm{M}_{i,-j}=\bm{M}_{i}+\left[\frac{\pi(\bm{x}_{j},\bm{\eta})(1-\pi(\bm{x}_{j},\bm{\eta}))}{1-\pi(\bm{x}_{j},\bm{\eta})(1-\pi(\bm{x}_{j},\bm{\eta}))v_{0}(\bm{x}_{j})}\right]\bm{M}_{i}\bm{f}(\bm{x}_{j})\bm{f}(\bm{x}_{j})^{\prime}\bm{M}_{i},\quad\textrm{
for }i=0,1,2.$

∎

Proof of $\Delta(\bm{x},\bm{x}_{i})$ in (34).

###### Proof.

Denote $\bm{F}^{*}$ as the updated model matrix for updated design
$\bm{X}^{*}$, and $\bm{W}_{i}^{*}$ for $i=0,1,2$ as the updated weight
matrices accordingly. Let

$\bm{G}_{0}=[\sqrt{\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta})}\bm{f}(\bm{x}),i\sqrt{\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta})}\bm{f}(\bm{x}_{i})],$

where $i=\sqrt{-1}$.

$\displaystyle\det(\bm{F}^{*^{\prime}}\bm{W}_{0}^{*}\bm{F}^{*})=\det\left(\sum_{j\neq
i}\pi(\bm{x}_{j},\bm{\eta})(1-\pi(\bm{x}_{j},\bm{\eta}))\bm{f}(\bm{x}_{j})\bm{f}(\bm{x}_{j})^{\prime}\right.$
$\displaystyle\left.-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}+\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta}))\bm{f}(\bm{x})\bm{f}(\bm{x})^{\prime}\right)$
$\displaystyle=\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F}+\bm{G}_{0}\bm{G}_{0}^{\prime})$
$\displaystyle=\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\det\left(\bm{I}_{2}+\bm{G}_{0}^{\prime}\bm{M}_{0}\bm{G}_{0}\right).$

It can be computed that

$\det\left(\bm{I}_{2}+\bm{G}_{0}^{\prime}\bm{M}_{0}\bm{G}_{0}\right)=\Delta_{0}(\bm{x},\bm{x}_{i}),$

where $\Delta_{0}(\bm{x},\bm{x}_{i})$ is

$\displaystyle\Delta_{0}(\bm{x},\bm{x}_{i})$
$\displaystyle=\left[1+\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta}))v_{0}(\bm{x})\right]\left[1-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))v_{0}(\bm{x}_{i})\right]$
$\displaystyle+\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta}))\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))v_{0}(\bm{x},\bm{x}_{i})^{2}.$

Similarly, define

$\displaystyle\bm{G}_{1}$
$\displaystyle=[\sqrt{\pi(\bm{x},\eta)}\bm{f}(\bm{x}),i\sqrt{\pi(\bm{x}_{i},\eta)}\bm{f}(\bm{x}_{i})]\quad\textrm{and}$
$\displaystyle\bm{G}_{2}$
$\displaystyle=[\sqrt{(1-\pi(\bm{x},\bm{\eta})}\bm{f}(\bm{x}),i\sqrt{(1-\pi(\bm{x}_{i},\bm{\eta})}\bm{f}(\bm{x}_{i})].$

Following a similar derivation,

$\displaystyle\det(\bm{F}^{*^{\prime}}\bm{W}_{1}^{*}\bm{F}^{*}+\rho\bm{R})$
$\displaystyle=\det(\bm{F}^{\prime}\bm{W}_{1}\bm{F}+\rho\bm{R})\Delta_{1}(\bm{x},\bm{x}_{i}),$
$\displaystyle\det(\bm{F}^{*^{\prime}}\bm{W}_{2}^{*}\bm{F}^{*}+\rho\bm{R})$
$\displaystyle=\det(\bm{F}^{\prime}\bm{W}_{2}\bm{F}+\rho\bm{R})\Delta_{2}(\bm{x},\bm{x}_{i}),$

where $\Delta_{1}(\bm{x},\bm{x}_{i})$ and $\Delta_{2}(\bm{x},\bm{x}_{i})$ are

$\displaystyle\Delta_{1}(\bm{x},\bm{x}_{i})$
$\displaystyle=\left[1+\pi(\bm{x},\bm{\eta})v_{1}(\bm{x})\right]\left[1-\pi(\bm{x}_{i},\bm{\eta})v_{1}(\bm{x}_{i})\right]+\pi(\bm{x},\bm{\eta})\pi(\bm{x}_{i},\bm{\eta})v_{1}(\bm{x},\bm{x}_{i})^{2},$
$\displaystyle\Delta_{2}(\bm{x},\bm{x}_{i})$
$\displaystyle=\left[1+(1-\pi(\bm{x},\bm{\eta}))v_{2}(\bm{x})\right]\left[1-(1-\pi(\bm{x}_{i},\bm{\eta}))v_{2}(\bm{x}_{i})\right]$
$\displaystyle+(1-\pi(\bm{x},\bm{\eta}))(1-\pi(\bm{x}_{i},\bm{\eta}))v_{2}(\bm{x},\bm{x}_{i})^{2}.$

Thus $\Delta(\bm{x},\bm{x}_{i})$ is computed as in (34). ∎

Proof of update formulas $\bm{M}_{i}^{*}$ for $i=0,1,2$

###### Proof.

Use the same notation of $\bm{G}_{i}$ for $i=0,1,2$ as in the previous proof.
Define the functions

$\displaystyle\bm{S}_{i}(\bm{x})=\bm{M}_{i}\bm{f}(\bm{x})\bm{f}(\bm{x})^{\prime}\bm{M}_{i},\quad\textrm{for
}i=0,1,2,$
$\displaystyle\bm{S}_{i}(\bm{x}_{1},\bm{x}_{2})=\bm{M}_{i}\bm{f}(\bm{x}_{1})\bm{f}(\bm{x}_{2})^{\prime}\bm{M}_{i},\quad\textrm{for
}i=0,1,2.$

It is straightforward to derive

$\displaystyle\bm{M}_{0}^{*}$
$\displaystyle=\left(\bm{F}^{\prime}\bm{W}_{0}\bm{F}+\bm{G}_{0}\bm{G}_{0}^{\prime}\right)^{-1}$
$\displaystyle=\bm{M}_{0}-\bm{M}_{0}\bm{G}_{0}\left(\bm{I}_{2}+\bm{G}_{0}^{\prime}\bm{M}_{0}\bm{G}_{0}\right)^{-1}\bm{G}^{\prime}\bm{M}_{0}.$

For simpler notation, denote
$a(\bm{x})=\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta}))$ and
$a(\bm{x}_{i})=\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))$.

$\displaystyle\left(\bm{I}_{2}+\bm{G}_{0}^{\prime}\bm{M}_{0}\bm{G}_{0}\right)^{-1}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{ll}1+a(\bm{x})v_{0}(\bm{x}),&i\sqrt{a(\bm{x})a(\bm{x}_{i})}v_{0}(\bm{x},\bm{x}_{i})\\\
i\sqrt{a(\bm{x})a(\bm{x}_{i})}v_{0}(\bm{x},\bm{x}_{i}),&1-a(\bm{x}_{i})v_{0}(\bm{x}_{i})\end{array}\right)^{-1}$
$\displaystyle=$
$\displaystyle\frac{1}{\Delta_{0}(\bm{x},\bm{x}_{i})}\left(\begin{array}[]{ll}1-a(\bm{x}_{i})v_{0}(\bm{x}_{i}),&-i\sqrt{a(\bm{x})a(\bm{x}_{i})}v_{0}(\bm{x},\bm{x}_{i})\\\
-i\sqrt{a(\bm{x})a(\bm{x}_{i})}v_{0}(\bm{x},\bm{x}_{i}),&1+a(\bm{x})v_{0}(\bm{x})\end{array}\right).$

$\displaystyle\bm{M}_{0}^{*}$
$\displaystyle=\bm{M}_{0}-\frac{1}{\Delta_{0}(\bm{x},\bm{x}_{i})}\bm{M}_{0}\left[\sqrt{a(\bm{x})}v_{0}(\bm{x}),i\sqrt{a(\bm{x}_{i})}v_{0}(\bm{x}_{i})\right]$
$\displaystyle\times$
$\displaystyle\left(\begin{array}[]{ll}1-a(\bm{x}_{i})v_{0}(\bm{x}_{i}),&-i\sqrt{a(\bm{x})a(\bm{x}_{i})}v_{0}(\bm{x},\bm{x}_{i})\\\
-i\sqrt{a(\bm{x})a(\bm{x}_{i})}v_{0}(\bm{x},\bm{x}_{i}),&1+a(\bm{x})v_{0}(\bm{x})\end{array}\right)\times\left[\begin{array}[]{c}\sqrt{a(\bm{x})}\bm{f}(\bm{x})^{\prime}\\\
i\sqrt{a_{i}(\bm{x}_{i})}\bm{f}(\bm{x}_{i})^{\prime}\end{array}\right]\bm{M}_{0}.$
$\displaystyle=\bm{M}_{0}-\frac{1}{\Delta_{0}(\bm{x},\bm{x}_{i})}\left\\{\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta}))\left[1-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))v_{0}(\bm{x}_{i})\right]\bm{S}_{0}(\bm{x})\right.$
$\displaystyle\left.-\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))\left[1+\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta}))v_{0}(\bm{x})\right]\bm{S}_{0}(\bm{x}_{i})\right.$
$\displaystyle+\left.\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))\pi(\bm{x},\bm{\eta})(1-\pi(\bm{x},\bm{\eta}))v_{0}(\bm{x},\bm{x}_{i})\left[\bm{S}_{0}(\bm{x},\bm{x}_{i})+\bm{S}_{0}(\bm{x}_{i},\bm{x})\right]\right\\}$

Update formulas for $\bm{M}_{1}^{*}$ and $\bm{M}_{2}^{*}$ can be derived
similarly as follows.

$\displaystyle\bm{M}_{1}^{*}=\bm{M}_{1}$
$\displaystyle-\frac{1}{\Delta_{1}(\bm{x},\bm{x}_{i})}\left\\{\pi(\bm{x},\bm{\eta})\left[1-\pi(\bm{x}_{i},\bm{\eta})v_{1}(\bm{x}_{i})\right]\bm{S}_{1}(\bm{x})-\pi(\bm{x}_{i},\bm{\eta})\left[1+\pi(\bm{x},\bm{\eta})v_{1}(\bm{x})\right]\bm{S}_{1}(\bm{x}_{i})\right.$
$\displaystyle+\left.\pi(\bm{x}_{i},\bm{\eta})\pi(\bm{x},\bm{\eta})v_{1}(\bm{x},\bm{x}_{i})\left[\bm{S}_{1}(\bm{x},\bm{x}_{i})+\bm{S}_{1}(\bm{x}_{i},\bm{x})\right]\right\\},$
$\displaystyle\bm{M}_{2}^{*}=\bm{M}_{2}$
$\displaystyle-\frac{1}{\Delta_{2}(\bm{x},\bm{x}_{i})}\left\\{(1-\pi(\bm{x},\bm{\eta}))\left[1-(1-\pi(\bm{x}_{i},\bm{\eta}))v_{2}(\bm{x}_{i})\right]\bm{S}_{2}(\bm{x})\right.$
$\displaystyle\left.-(1-\pi(\bm{x}_{i},\bm{\eta}))\left[1+(1-\pi(\bm{x},\bm{\eta}))v_{2}(\bm{x})\right]\bm{S}_{2}(\bm{x}_{i})\right.$
$\displaystyle+\left.(1-\pi(\bm{x}_{i},\bm{\eta}))(1-\pi(\bm{x},\bm{\eta}))v_{2}(\bm{x},\bm{x}_{i})\left[\bm{S}_{2}(\bm{x},\bm{x}_{i})+\bm{S}_{2}(\bm{x}_{i},\bm{x})\right]\right\\}.$

∎

##  A2 Appendix A2. Table

Table A1: Local Bayesian $D$-optimal Designs for $\rho=0$ and 0.3 and other three alternative designs. The values in columns 7-11 are frequencies.  Point | $x_{1}$ | $x_{2}$ | $x_{3}$ | $x_{4}$ | $x_{5}$ |  $D_{QQ}$ $\rho=0$ |  $D_{QQ}$ $\rho=0.3$ | $D_{G}$ | $D_{L}$ | $D_{n}$
---|---|---|---|---|---|---|---|---|---|---
1 | -1 | -1 | -1 | -1 | -1 | 1 | 2 | 1 | 1 | 2
2 | 1 | -1 | -1 | -1 | -1 | 1 | 1 | 2 | 1 | 1
3 | -1 | 1 | -1 | -1 | -1 | 1 | 1 | 2 | 1 | 1
4 | 1 | 1 | -1 | -1 | -1 | 2 | 1 | 2 | 1 | 2
5 | -1 | -1 | 1 | -1 | -1 | 1 | 0 | 2 | 1 | 1
6 | 1 | -1 | 1 | -1 | -1 | 1 | 1 | 0 | 1 | 1
7 | -1 | 1 | 1 | -1 | -1 | 2 | 2 | 1 | 1 | 2
8 | 1 | 1 | 1 | -1 | -1 | 1 | 2 | 1 | 1 | 1
9 | -1 | -1 | -1 | 0 | -1 | 2 | 2 | 1 | 1 | 1
10 | 1 | -1 | -1 | 0 | -1 | 0 | 0 | 1 | 1 | 1
11 | -1 | 1 | -1 | 0 | -1 | 1 | 1 | 2 | 1 | 2
12 | 1 | 1 | -1 | 0 | -1 | 1 | 2 | 1 | 1 | 1
13 | -1 | -1 | 1 | 0 | -1 | 1 | 1 | 1 | 1 | 2
14 | 1 | -1 | 1 | 0 | -1 | 2 | 2 | 1 | 1 | 1
15 | -1 | 1 | 1 | 0 | -1 | 1 | 1 | 0 | 1 | 0
16 | 1 | 1 | 1 | 0 | -1 | 1 | 1 | 1 | 1 | 2
17 | -1 | -1 | -1 | 1 | -1 | 1 | 0 | 1 | 1 | 2
18 | 1 | -1 | -1 | 1 | -1 | 2 | 2 | 1 | 1 | 1
19 | -1 | 1 | -1 | 1 | -1 | 2 | 2 | 1 | 1 | 1
20 | 1 | 1 | -1 | 1 | -1 | 1 | 1 | 1 | 1 | 1
21 | -1 | -1 | 1 | 1 | -1 | 2 | 2 | 1 | 1 | 1
22 | 1 | -1 | 1 | 1 | -1 | 0 | 1 | 1 | 1 | 1
23 | -1 | 1 | 1 | 1 | -1 | 0 | 1 | 1 | 1 | 2
24 | 1 | 1 | 1 | 1 | -1 | 2 | 1 | 1 | 1 | 1
25 | -1 | -1 | -1 | -1 | 0 | 1 | 0 | 1 | 1 | 0
26 | 1 | -1 | -1 | -1 | 0 | 0 | 0 | 1 | 1 | 1
27 | -1 | 1 | -1 | -1 | 0 | 0 | 0 | 0 | 1 | 0
28 | 1 | 1 | -1 | -1 | 0 | 0 | 1 | 1 | 1 | 0
29 | -1 | -1 | 1 | -1 | 0 | 0 | 1 | 1 | 1 | 1
30 | 1 | -1 | 1 | -1 | 0 | 0 | 0 | 0 | 1 | 0
31 | -1 | 1 | 1 | -1 | 0 | 0 | 0 | 0 | 1 | 0
32 | 1 | 1 | 1 | -1 | 0 | 1 | 0 | 1 | 1 | 0
33 | -1 | -1 | -1 | 0 | 0 | 0 | 0 | 1 | 1 | 0
34 | 1 | -1 | -1 | 0 | 0 | 1 | 1 | 1 | 1 | 1
35 | -1 | 1 | -1 | 0 | 0 | 0 | 0 | 0 | 1 | 0
36 | 1 | 1 | -1 | 0 | 0 | 1 | 0 | 1 | 1 | 0
37 | -1 | -1 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0
38 | 1 | -1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0
39 | -1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0
40 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0
41 | -1 | -1 | -1 | 1 | 0 | 0 | 2 | 1 | 0 | 0
42 | 1 | -1 | -1 | 1 | 0 | 0 | 0 | 1 | 1 | 0
43 | -1 | 1 | -1 | 1 | 0 | 0 | 0 | 0 | 1 | 0
44 | 1 | 1 | -1 | 1 | 0 | 1 | 0 | 1 | 0 | 0
45 | -1 | -1 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0
46 | 1 | -1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 | 0
47 | -1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 0
48 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 1
49 | -1 | -1 | -1 | -1 | 1 | 1 | 1 | 1 | 1 | 1
50 | 1 | -1 | -1 | -1 | 1 | 2 | 2 | 1 | 1 | 2
51 | -1 | 1 | -1 | -1 | 1 | 1 | 1 | 0 | 1 | 1
52 | 1 | 1 | -1 | -1 | 1 | 1 | 1 | 1 | 1 | 1
53 | -1 | -1 | 1 | -1 | 1 | 2 | 1 | 1 | 1 | 2
54 | 1 | -1 | 1 | -1 | 1 | 0 | 0 | 0 | 1 | 0
55 | -1 | 1 | 1 | -1 | 1 | 1 | 1 | 1 | 1 | 1
56 | 1 | 1 | 1 | -1 | 1 | 1 | 2 | 1 | 1 | 2
57 | -1 | -1 | -1 | 0 | 1 | 1 | 1 | 1 | 1 | 2
58 | 1 | -1 | -1 | 0 | 1 | 1 | 1 | 1 | 1 | 0
59 | -1 | 1 | -1 | 0 | 1 | 2 | 2 | 1 | 1 | 1
60 | 1 | 1 | -1 | 0 | 1 | 1 | 1 | 1 | 1 | 2
61 | -1 | -1 | 1 | 0 | 1 | 2 | 2 | 1 | 0 | 1
62 | 1 | -1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 2
63 | -1 | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 1
64 | 1 | 1 | 1 | 0 | 1 | 2 | 2 | 1 | 0 | 1
65 | -1 | -1 | -1 | 1 | 1 | 2 | 1 | 1 | 1 | 2
66 | 1 | -1 | -1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
67 | -1 | 1 | -1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
68 | 1 | 1 | -1 | 1 | 1 | 1 | 2 | 1 | 1 | 2
69 | -1 | -1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
70 | 1 | -1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1
71 | -1 | 1 | 1 | 1 | 1 | 2 | 2 | 1 | 1 | 2
72 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1