#  $D$-optimal Approximate Design for Binary Regression and Quantal Response
in Toxicology Studies

Elvis Han Cui
Department of Biostatistics, UCLA
<EMAIL_ADDRESS>

###### Abstract

We provide a systematic treatment of $D$-optimal design for binary regression
and quantal response models in toxicology studies. For the two-parameter case,
we provide an analytical equation (WC equation) for computing the $D$-optimal
design quickly and when analytical solution is not available, we apply
particle swarm optimization to solve for the $D$-optimal design. Examples with
various link functions are given as well as the sensitivity functions. We
extend the two-parameter case to three-parameter case by providing a neat
formula for the determinant of the information matrix. We also suggest
practitioners to work with the neat formula to derive optimal designs for
three-parameter binary regression models.

_K_ eywords $D$-optimality $\cdot$ Quantal response $\cdot$ Binary regression
$\cdot$ WC equation

##  1 Preliminaries

In this section, we first briefly discuss quantal response in toxicology
studies and then introduce the motivation and the basic concepts of optimal
approximate design and review binary regression for toxicology studies.

###  1.1 Quantal Response in Toxicology Studies

The outcome of interest in many toxicological studies is qualitative in
nature. In a vasr majority of these experiments, the outcome is quantal and
binary [Raz22], to name a few: beetle mortality and embryogenic anthers in
toxicology studies [DB18]; tumor progression status in cancer studies
[ASE+17]; low-density lipo-protein (LDL) cholesterol levels (desirable vs
undesirable) [CCHW22]. Further, in developmental toxicity studies, pregnant
female animals are being exposed to a teratogen during a specific time of
pregnancy [Raz22]. To estimate the potential risk from a defined source of
hazard (e.g., exposure to a teratogen) quantitatively and qualitatively, we
usually starts with fitting a dose-response curve to the data. For quantal and
binary outcome, the range of the curve is within 0 to 1. In Figure 1, we
demonstrate the dose-response curve using 7 different potential “link
functions” (for details, see 1.3), the $x$-axis refers to potential dose level
(potential means the range can be shifted and re-scaled) and the $y$-axis
refers to the binary outcome probability. As Figure 1 has shown, when the dose
increases, the probability of the quantal / binary outcome is monotonically
increasing.

Figure 1: Dose response curves for different link functions.

We are particularly interested estimating the curve as accurate as possible so
that the downstream risk assessment [Raz22] can proceed smoothly. Hence, it is
desirable for us to estimate the curve as accurate as possible and
simultaneously reduce the required total number of doses [WL96, ZAKW98, ZW00,
BZWW06]. To achieve this goal, we apply the theory of optimal design and we
give a brief introduction in section 1.2.

###  1.2 Motivation and Basic Concepts of Optimal Design

In a dose–response experiment, decisions regarding the dose range, the number
of doses, the dose levels, and the number of experimental units at each dose
are sometimes made predicated on nebulous criteria. These are design issues
that can potentially have a substantial impact on the quality of the
statistical inference at the end of the study, yet they are decided in some
cases on an ad-hoc basis. Frequently, an equal number of experimental units
are assigned at each dose. When the doses are equally spaced, these are called
uniform designs in the statistical literature and while they are appealing and
intuitive, it has been shown that they can be inefficient, depending on the
goal of the study and the underlying model assumed. For example, [WL96] showed
that performance of such designs can depend sensitively on the choice of the
number of doses in a uniform design, the model, and the optimality criteria.
Therefore, each aspect in the design of the study must be carefully considered
to realize maximum accuracy in the information. Such attention to detail will
enhance reproducibility, thus addressing a current issue in animal
experimentation [Gil06] and reducing the overall cost of experiments. More
specifically, if the current cost for producing a new drug is 10 dollars per
dose, then using optimal design theory, one is able to reduce the cost to 5
dollars per dose.

To optimally design an experiment, model assumptions are required to work out
the mathematical and statistical details. Invariably, the goal is formulated
as an objective function defined on the user-specified dose range (or design
interval) that depends on the statistical model and the design. The
optimization of the criterion can then be performed among a specific class of
designs, for example, among all designs with five doses, or among all designs
on a given dose interval. The resulting optimal design is therefore model-
based and, as a consequence, can be highly model-dependent, suggesting that
choice of a statistical model for the dose–response study is also important.

Optimal approximate designs in clinical trials can help investigators achieve
higher quality results for the given resource constraints [SSQ+06, JM13,
SRW20, ZWY21]. The creation of this field can be traced back to [Smi18]. From
1950s to 1980s, the field of approximate design has witnessed a booming
development [Fed72, Kie74, Páz86, ADT07, Sil13] and we give a brief
introduction below.

Consider a linear model $\mathbb{E}(y)=\mu=\beta^{T}f(x)$ where $\mu$ is the
expectation of $y$. A $k$-point design $\xi$ is a $2\times k$ matrix of the
following form

$\displaystyle\xi$
$\displaystyle=\left(\begin{matrix}x_{1}&x_{2}&\cdots&x_{k-1}&x_{k}\\\
p_{1}&p_{2}&\cdots&p_{k-1}&p_{k}\end{matrix}\right)$

where $x_{i}$’s are called design points and $p_{i}$’s are non-negative
weights that sum to 1. In practice, we usually have a total of $n$
observations and $np_{i}$ is not an integer in general. Hence, we choose the
closest integer of $np_{i}$ and assign the corresponding dose $x_{i}$ to these
individuals. For this reason, $\xi$ is referred as approximate design in
literature. The information matrix associated with design $\xi$ is

$\mathfrak{M}(\xi)=\int_{\mathcal{X}}f(x)f(x)^{T}\xi(dx)=\sum_{i=1}^{k}p_{i}f(x_{i})f(x_{i})^{T}$

where $\xi(dx)$ is the measure induced by $p_{i}$’s. The optimal design seeks
to find a $\xi^{*}$ that minimizes $\phi(\mathfrak{M}(\xi))$ where
$\phi(\cdot)$ is a real-valued function. Commons choices of $\phi(\cdot)$ and
their terminologies are given in the following table

Table 1: List of Common Choices of $\phi(\cdot)$ Optimalitya |  Choice of $\phi$ | Remarks
---|---|---
A | $\text{Tr}(\mathfrak{M}^{-1})\textsuperscript{b}$ | Sum of variances
c | $\text{Var}(g(\widehat{\beta}))\textsuperscript{c}$ | Variance of $g(\widehat{\beta})$
D | $\log\det\mathfrak{M}^{-1}$ | Log-volume of the ellipse
E | $\min\lambda_{i}(\mathfrak{M})=\lambda_{min}\textsuperscript{d}$ | Length of minor axis
G | $\max_{i}(\mathfrak{M}^{-1})_{ii}$ | Maximum of $\text{Var}(\widehat{\beta})$
I | $\int_{\mathcal{X}}f(x)^{T}\mathfrak{M}^{-1}f(x)\mu(dx)$ | Integrated variance

a For example, the first line reads $A$-optimality.

b Tr refers to the trace function of a matrix.

c $g$ is a function of $\beta$, $\widehat{\beta}$ is the MLE of $\beta$ and
the variance of $g(\widehat{\beta})$ can be derived using Delta method, i.e.,
$\widehat{\text{Var}}(g(\widehat{\beta}))=\nabla
g(\widehat{\beta})^{T}\mathfrak{M}^{-1}\nabla g(\widehat{\beta})$ and $\nabla$
refers to the gradient operator.

d $\lambda$’s refer to eigenvalues of $\mathfrak{M}$.

Further, if we extend the linear model framework to generalized linear models,
then the information matrix depends on parameters (see section 2). In this
case, one usually plug-in plausible parameter values and then calculate the
optimal design $\xi^{*}$. We call it $\phi$-optimal design, where $\phi$
refers to $D$-, $A$-, $c$-, $E$-, etc. In practice, researchers would like to
consider different criteria simultaneously, leading to the so-called compound
criteria. For a comprehensive review of different optimality criteria, see
Chapter 10 of [ADT07] or the review paper by [Fed10].

To verify that the resulting design is globally optimal, that is, optimal
among all possible designs, one needs to apply the equivalence theorem [Kie74]
and plot the sensitivity functions to check. Different optimality criteria
corresponds to different types of equivalence theorems and sensitivity
functions, hence for brevity, we only state the equivalence theorem for
$D$-optimal design in section 1.3. In addition, [CCW22] provides a
comprehensive review on the application of PSO in optimal approximate design.

###  1.3 Basic Concepts of Binary Regression

In the following, we first give a review on binary regression and then provide
the equivalence theorem for $D$-optimal design. Assume that for
$i=1,\cdots,n$, the response $y_{i}$ is a binary outcome with covariate
$x_{i}\in\mathbb{R}^{d}$. The $y_{i}$’s are independently distributed and the
density is

$p(y|x,\pi)=p(y|\eta(x,\beta))=\exp\left(y\ln\frac{\pi}{1-\pi}+\ln(1-\pi)\right)$

where $\beta$ is a $p$-dimensional parameter of interest and
$\pi=F(\eta)=\int_{-\infty}^{\eta}f(s)ds$. Let $\mathcal{X}$ be the design
space and the design $\xi$ be

$\displaystyle\xi$
$\displaystyle=\left(\begin{matrix}x_{1}&x_{2}&\cdots&x_{n-1}&x_{n}\\\
p_{1}&p_{2}&\cdots&p_{n-1}&p_{n}\end{matrix}\right)$

where for all $i$, $x_{i}\in\mathcal{X}$, $p_{i}\geq 0$ and
$\sum_{i=1}^{n}p_{i}=1$. Let
$g(x,\beta)=\mathbb{E}\left(-\frac{\partial^{2}\log
p(y|x,\pi)}{\partial\beta\partial\beta^{T}}\right)$ be the Fisher information
associated with a single point $x$, then

$g(x,\beta)=\omega xx^{T}$

where $\omega=F^{\prime}(\eta)^{2}/(\pi(1-\pi))=f(\eta)^{2}/(\pi(1-\pi))$. The
information matrix associated with the design $\xi$ is

$\displaystyle\mathfrak{M}(\xi)$
$\displaystyle=\int_{\mathcal{X}}g(x,\beta)\xi(dx)$
$\displaystyle=\sum_{i=1}^{n}p_{i}\omega_{i}x_{i}x_{i}^{T},$ (1.1)

where

$\omega_{i}=\frac{f(\eta_{i})^{2}}{(\pi_{i}(1-\pi_{i}))},\
\sum_{i=1}^{n}p_{i}=1,\ p_{i}\geq 0.$

Here are some examples of commonly used models in practice.

  * •

(Logit) The most famous model is the logistic regression with logit link.

$\displaystyle f(\eta_{i})$
$\displaystyle=\frac{\exp(\eta_{i})}{(1+\exp(\eta_{i}))^{2}},$
$\displaystyle\eta_{i}$ $\displaystyle=\ln\frac{\pi_{i}}{1-\pi_{i}},$
$\displaystyle\pi_{i}$
$\displaystyle=\frac{\exp(\eta_{i})}{1+\exp(\eta_{i})},$
$\displaystyle\mathfrak{M}(\xi)$
$\displaystyle=\sum_{i=1}^{n}\frac{p_{i}\exp(\eta_{i})}{(1+\exp(\eta_{i}))^{2}}x_{i}x_{i}^{T}.$

  * •

(Probit) Prior to the presense of logit link, one uses the probit link.

$\displaystyle f(\eta_{i})$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}\eta_{i}^{2}\right),$
$\displaystyle\eta_{i}$ $\displaystyle=\Phi^{-1}(\pi_{i}),$
$\displaystyle\pi_{i}$ $\displaystyle=\Phi(\eta_{i}),$
$\displaystyle\mathfrak{M}(\xi)$
$\displaystyle=\sum_{i=1}^{n}\frac{p_{i}\exp(-\eta_{i}^{2})}{2\pi\Phi(-\eta_{i})\Phi(\eta_{i})}x_{i}x_{i}^{T}.$

where $\Phi(\cdot)$ is the cumulative function of standard normal.

  * •

(Laplace) If we want the rate of decay is faster than student $t$ but slower
than probit, then Laplace density is an alternative:

$\displaystyle f(\eta_{i})$
$\displaystyle=\frac{1}{2}\exp\left(-|\eta_{i}|\right),$
$\displaystyle\eta_{i}$ $\displaystyle=F^{-1}(\pi_{i}),$
$\displaystyle\pi_{i}$
$\displaystyle=F(\eta_{i})=\frac{1}{2}+\frac{1}{2}\text{sgn}(\eta_{i})\left(1-\exp(-|\eta_{i}|)\right),$
$\displaystyle\mathfrak{M}(\xi)$
$\displaystyle=\sum_{i=1}^{n}\frac{p_{i}\exp(-2|\eta_{i}|)}{4F(\eta_{i})S(\eta_{i})}x_{i}x_{i}^{T}$

where $\text{sgn}(\cdot)$ is the sign function and $S(\cdot)=1-F(\cdot)$.

  * •

(Cox regression) The associated distribution is called the Gumbel extreme
value distribution, and the link function is called the complementary log-log:

$\displaystyle f(\eta_{i})$
$\displaystyle=\exp\left(\eta_{i}-\exp(\eta_{i})\right),$
$\displaystyle\eta_{i}$ $\displaystyle=\ln\left(-\ln(1-\pi_{i})\right),$
$\displaystyle\pi_{i}$ $\displaystyle=1-\exp(-\exp(\eta_{i})),$
$\displaystyle\mathfrak{M}(\xi)$
$\displaystyle=\sum_{i=1}^{n}\frac{p_{i}\exp(2\eta_{i}-2\exp(\eta_{i}))}{\exp(-\exp(\eta_{i})-\exp(-2\exp(\eta_{i})))}x_{i}x_{i}^{T}.$

  * •

(Student t) The Student t distribution with degrees of freedom $k$ is useful
in hypothesis testing and its density is symmetric with respect to 0:

$\displaystyle f(\eta_{i})$
$\displaystyle=\frac{\Gamma(\frac{k+1}{2})}{\sqrt{k\pi}\Gamma(\frac{k}{2})}\left(1+\frac{\eta_{i}^{2}}{k}\right)^{-\frac{k+1}{2}}$
$\displaystyle\eta_{i}$ $\displaystyle=F^{-1}(\pi_{i})$ $\displaystyle\pi_{i}$
$\displaystyle=F(\eta_{i})=1-\frac{1}{2}\mathbf{I}_{x(\eta_{i})}\left(\frac{k}{2},\frac{1}{2}\right)$

where $\mathbf{I}$ is the regularized incomplete Beta function and
$x(s)=\frac{k}{s^{2}+k}$.

A $D$-optimal design seeks to find a design $\xi^{*}$ such that
$\det\mathfrak{M}(\xi)$ is maximized. It is well known that if we know the
$D$-optimal design for a $k$-parameter model is supported at $k$-points, then
all points are equally weighted [ADT07, Won21].

######  Lemma 1 (Property of $D$-optimal design).

If we know in advance that a $D$-optimal design for binary regression (Formula
1.1) has $k$ support points, then all design points have design weight
$\frac{1}{k}$.

###### Proof.

Suppose $n=k$ then we can write the determinant of formula 1.1 as

$\displaystyle\det\mathfrak{M}(\xi)$
$\displaystyle=\det\left(\sum_{i=1}^{k}p_{i}\omega_{i}x_{i}x_{i}^{T}\right)$
$\displaystyle=\det\left(\left(\begin{matrix}\sqrt{\omega_{1}}x_{1}&\cdots&\sqrt{\omega_{k}}x_{k}\end{matrix}\right)\left(\begin{matrix}p_{1}&0&\cdots&0\\\
0&p_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&p_{k}\end{matrix}\right)\left(\begin{matrix}\sqrt{\omega_{1}}x_{1}\\\
\cdots\\\ \sqrt{\omega_{k}}x_{k}\end{matrix}\right)\right)$
$\displaystyle=\det\left(\sum_{i=1}^{k}\omega_{i}x_{i}x_{i}^{T}\right)\det\left(\begin{matrix}p_{1}&0&\cdots&0\\\
0&p_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&p_{k}\end{matrix}\right)$
$\displaystyle=\det\left(\sum_{i=1}^{k}\omega_{i}x_{i}x_{i}^{T}\right)\left(\prod_{i=1}^{k}p_{i}\right)$
$\displaystyle\leq\det\left(\sum_{i=1}^{k}\omega_{i}x_{i}x_{i}^{T}\right)\left(\frac{\sum_{i=1}^{k}p_{i}}{k}\right)^{k}$
$\displaystyle=\det\left(\sum_{i=1}^{k}\omega_{i}x_{i}x_{i}^{T}\right)\left(\frac{1}{k}\right)^{k}$

where the inequality is due to the AM–GM inequality [HLP52]. The “=” is
attained if and only if $p_{1}=p_{2}=\cdots=p_{k}=\frac{1}{k}$. Hence, we
conclude that all design points have equal weights. ∎

Finally, to check a design $\xi^{*}$ is whether globally optimal or not (i.e.,
optimal among all possible designs), we use the following theorem and plot the
sensitivity function $\psi$.

######  Theorem 1.1 (Equivalence theorem [ADT07, Won21]).

Let $\mathfrak{M}(\xi)$ be the information matrix associated with design
$\xi$, then the following are equivalent ($\dim(\mathcal{X})=d$),

  1. 1.

The design $\xi^{*}$ is $D$-optimal, i.e.,
$\xi^{*}=\arg\min_{\xi}\det\mathfrak{M}(\xi)$.

  2. 2.

The inequality
$\psi(\mathbf{x},\beta)=w(\mathbf{x})(\mathbf{x}^{T}\mathfrak{M}(\xi^{*})^{-1}\mathbf{x})-d\leq
0$ holds for all $\mathbf{x}\in\mathcal{X}\subset\mathbb{R}^{d}$ where
$w(\mathbf{x})$ is a weight depending on the link $\eta(\mathbf{x})$ and
$\psi$ is called the sensitivity function.

The equivalence theorem says that if a design $\xi^{*}$ is $D$-optimal, the
the sensitivity function $\psi(\mathbf{x},\beta)$ is less or equal to 0 within
the design space $\mathcal{X}$. Further, the sensitivity function attains 0 at
the design points. Some preliminary work on $D$-optimal design for binary
regression is given in [KW00, BZWW06, HKO07, ADT07, KH12]. However, there lack
a detailed and unified framework for binary regression with different types of
link functions under $D$-optimality. Hence, we provide a systematic treatment
in the next two sections.

##  2 Two-parameter Binary Regression

In this section, we always assume $k=2$ so that the resulting design always
has $2$ equally support points. Then by formula 1.1, we have

$\displaystyle\det\mathfrak{M}(\xi)$
$\displaystyle=\frac{1}{4}\det\left(\begin{matrix}\omega_{1}+\omega_{2}&\omega_{1}x_{1}+\omega_{2}x_{2}\\\
\omega_{2}x_{2}+\omega_{1}x_{1}&\omega_{1}x_{1}^{2}+\omega_{2}x_{2}^{2}\end{matrix}\right)$
$\displaystyle\propto\left(\omega_{1}+\omega_{2}\right)\left(\omega_{1}x_{1}^{2}+\omega_{2}x_{2}^{2}\right)-\left(\omega_{1}x_{1}+\omega_{2}x_{2}\right)^{2}$
$\displaystyle=\left[\omega_{1}^{2}x_{1}^{2}+\omega_{1}\omega_{2}(x_{1}^{2}+x_{2}^{2})+\omega_{2}^{2}x_{2}^{2}\right]-$
$\displaystyle\ \ \ \
\left[\omega_{1}^{2}x_{1}^{2}+2\omega_{1}\omega_{2}x_{1}x_{2}+\omega_{2}^{2}x_{2}^{2}\right]$
$\displaystyle=\omega_{1}\omega_{2}\left(x_{1}^{2}-2x_{1}x_{2}+x_{2}^{2}\right)$
$\displaystyle\propto\frac{f(\eta_{1})^{2}f(\eta_{2})^{2}}{F(\eta_{1})S(\eta_{1})F(\eta_{2})S(\eta_{2})}(x_{1}-x_{2})^{2}$
(2.1)

where for $i=1,2$,
$\eta_{i}=\beta_{0}+\beta_{1}x_{i},F(\eta_{i})=\int_{-\infty}^{\eta_{i}}f(s)ds,S(\eta_{i})=1-F(\eta_{i})$.
Plug-in $\eta=\beta_{0}+\beta_{1}x$, then we have

$\det\mathfrak{M}(\xi)=\frac{f(\eta_{1})^{2}f(\eta_{2})^{2}}{F(\eta_{1})S(\eta_{1})F(\eta_{2})S(\eta_{2})}\left(\frac{\eta_{1}-\eta_{2}}{\beta_{1}}\right)^{2}$

and taking the logarithm of $\det\mathfrak{M}(\xi)$ and setting the derivative
w.r.t. $\eta_{1}$ equal to zero gives

$\displaystyle\frac{2f^{\prime}(\eta_{1})}{f(\eta_{1})}-\frac{f(\eta_{1})}{F(\eta_{1})}+\frac{f(\eta_{1})}{S(\eta_{1})}+\frac{2}{\eta_{1}-\eta_{2}}=0$
(2.2)

We denote the above the key equation and call it the WC equation where $WC$
stands for Wong and Cui. Now it is natural to consider, if we are given a
maximizer $(\eta_{1}^{*},\eta_{2}^{*})$ of Formula 2.1, is it unique? The
short answer is no unless $\eta_{1}^{*}$ and $\eta_{2}^{*}$ is symmetric
around $2a$ where $a$ is a real number and we provide a lemma below.

######  Lemma 2.

If $f(s)$ is symmetric, i.e., $f(a+s)=f(a-s)$ for some $a$, then for any
design

$\xi^{*}=\left(\begin{matrix}\frac{\eta_{1}^{*}-\beta_{0}}{\beta_{1}}&\frac{\eta_{2}^{*}-\beta_{0}}{\beta_{1}}\\\
0.5&0.5\end{matrix}\right)$

we have

$\det\mathfrak{M}(\xi^{*})=\det\mathfrak{M}(\xi^{\prime})$

where

$\xi^{\prime}=\left(\begin{matrix}\frac{2a-\eta_{2}^{*}-\beta_{0}}{\beta_{1}}&\frac{2a-\eta_{1}^{*}-\beta_{0}}{\beta_{1}}\\\
0.5&0.5\end{matrix}\right)$

###### Proof.

For $i=1,2$, let $x_{i}^{*}=\frac{\eta_{1}^{*}-\beta_{0}}{\beta_{1}}$ and
$x_{i}^{\prime}=\frac{2a-\eta_{2}^{*}-\beta_{0}}{\beta_{1}}$, then

$(x_{1}^{*}-x_{2}^{*})^{2}=\left(\frac{\eta_{1}^{*}-\eta_{2}^{*}}{\beta_{1}}\right)^{2}=(x_{1}^{\prime}-x_{2}^{\prime})^{2}$

Next, let $\eta_{1}^{\prime}=2a-\eta_{2}^{*}$ and
$\eta_{2}^{\prime}=2a-\eta_{1}^{*}$, we have

$f(\eta_{1}^{*})=f(a+(\eta_{1}^{*}-a))=f(a-(\eta_{1}^{*}-a))=f(2a-\eta_{1}^{*})=f(\eta_{2}^{\prime})$

and similarly, $f(\eta_{2}^{*})=f(\eta_{1}^{\prime})$.

Finally, we have

$\displaystyle F(\eta_{1}^{*})$ $\displaystyle=F(a+(\eta_{1}^{*}-a))$
$\displaystyle=S(\eta_{2}^{\prime})$

and $S(\eta_{1}^{*})=F(\eta_{2}^{\prime})$,
$F(\eta_{2}^{*})=S(\eta_{1}^{\prime})$,
$S(\eta_{2}^{*})=F(\eta_{1}^{\prime})$. ∎

###  2.1 Symmetric Densities

WLOG, we may assume that $f(s)=f(-s)$. As an example, one such $f$ is logistic
density $f(s)=e^{-s}/(1+e^{-s})^{2}$. For a symmetric two-point design, we let
$\eta_{1}=-\eta_{2}$, then

$\det\mathfrak{M}(\xi)=\frac{f(\eta_{1})^{4}}{F(\eta_{1})^{2}S(\eta_{1})^{2}}\left(\frac{2\eta_{1}}{\beta_{1}}\right)^{2}$

and taking the logarithm of $\det\mathfrak{M}(\xi)$ and setting the derivative
w.r.t. $\eta_{1}$ equal to zero gives

$\displaystyle\frac{2f^{\prime}(\eta_{1})}{f(\eta_{1})}-\frac{f(\eta_{1})}{F(\eta_{1})}+\frac{f(\eta_{1})}{S(\eta_{1})}+\frac{1}{\eta_{1}}=0$
(2.3)

The resulting solutions provide the design points of the $D$-optimal and
$G$-optimal designs among all $2$-point designs. To verify that it is optimal
among all possible designs, we need to calculate the sensitivity function
based on the Theorem 1.1.

In the following, we apply the key equation to a few examples and verify the
results using PSO . In short, we write $\eta$ for $\beta_{0}+\beta_{1}x$ and
$x$ can be solved by $x=(\eta-\beta_{0})/\beta_{1}$.

######  Example 2.1 (Logit).

For this problem, $f(\eta)={\exp(\eta)}/{(1+\exp(\eta))^{2}}$ and
$w=1/((1+\exp\eta)(1-\exp\eta))$. Plug-in all necessary elements, the key
equation is

$2-\frac{4\exp\eta}{1+\exp\eta}+\frac{2}{\eta}=0$

Solving it numerically, we obtain $\eta_{1}=+1.5434$ and $\eta_{2}=-1.5434$.
Hence, the resulting design is

$\displaystyle\begin{cases}x_{1}=(+1.5434-\beta_{0})/\beta_{1}\\\
x_{2}=(-1.5434-\beta_{0})/\beta_{1}\end{cases}$ (2.4)

The Figure 2 demonstrates the sensitivity functions of two locally D-optimal
designs with logit link and specified parameter values.

Figure 2: Sensitivity functions for logit link.

######  Example 2.2 (Probit).

For this problem, $f(\eta)=\exp\left(-\frac{1}{2}\eta^{2}\right)/\sqrt{2\pi}$
and $\Phi(\eta)=\int_{\infty}^{\eta}f(s)ds$ and
$w=\exp(-\eta^{2})/(2\pi\Phi(\eta)(1-\Phi(\eta)))$. Plug-in all necessary
elements, the key equation is

$\frac{2}{\eta}-4\eta-2f(\eta)\left(\frac{1}{\Phi(\eta)}-\frac{1}{1-\Phi(\eta)}\right)=0$

Solving it numerically, we obtain $\eta_{1}=+1.1381$ and $\eta_{2}=-1.1381$.
Hence, the resulting design is

$\displaystyle\begin{cases}x_{1}=(+1.1381-\beta_{0})/\beta_{1}\\\
x_{2}=(-1.1381-\beta_{0})/\beta_{1}\end{cases}$ (2.5)

The two panels of Figure 3 demonstrates the sensitivity functions of two
locally D-optimal designs with probit link and specified parameter values.

Figure 3: Sensitivity functions for probit link.

######  Example 2.3 (Laplace).

For this problem, we have $f(\eta)=\exp(-|\eta|)/2$ and
$w=1/(2\exp(|\eta|)-1)$. Plug-in all necessary elements, the key equation is

$-4\
\text{sgn}(\eta)-\frac{2\exp(-|\eta|)}{1+\text{sgn}(\eta)(1-\exp(-|\eta|))}+\frac{2\exp(-|\eta|)}{1-\text{sgn}(\eta)(1-\exp(-|\eta|))}+\frac{2}{\eta}=0$

Solving it numerically, we obtain $\eta_{1}=+0.7680$ and $\eta_{2}=-0.7680$.
Hence, the resulting design is

$\displaystyle\begin{cases}x_{1}=(+0.7680-\beta_{0})/\beta_{1}\\\
x_{2}=(-0.7680-\beta_{0})/\beta_{1}\end{cases}$ (2.6)

However, the left panel of Figure 4 has shown that $(+0.7680,-0.7680)$ is NOT
a locally D-optimal design. According to Federov’s algorithm [ADT07], it
suggests that we need to add a design point at 0. This is empirically verified
by Particle Swarm Optimization (PSO) using the Python package “pyswarms”
[Mir18]. The right panel of Figure 4 has shown the sensitivity function of the
three point design generated by PSO.

Figure 4: Sensitivity functions for Laplace link.

Now is it natural to ask that when does the two point design is indeed
globally optimal, i.e., is it possible that we have a three-or-more point
$D$-optimal design? Due to Caratheodory’s theorem (for example, see appendix
in [Sil13]), the $D$-optimal design for a two parameter binary regression
model has AT MOST three support points. To illustrate this idea, we add
another example with an artificial regression function.

######  Example 2.4 (Logit with an artificial regression function).

Instead of $\eta=\beta_{0}+\beta_{1}x$, suppose now we have

$\eta=\beta_{0}+\beta_{1}x+\frac{1}{|x|+1}$

In this case, the WC equation boils down to
$1.5434=\beta_{0}+\beta_{1}x+\frac{1}{|x|+1}$ and
$-1.5434=\beta_{0}+\beta_{1}x+\frac{1}{|x|+1}$. Alternatively, we solve the
$D$-optimal design using PSO (50 particles with 500 iterations, hyper-
parameters are $c_{1}=0.5,c_{2}=0.3$ and $w=0.9$) and the results are shown in
left panel of Figure 5. Note that still we only have 2 parameters $\beta_{0}$
and $\beta_{1}$ but the $D$-optimal has $3$ support points
($x_{1}=0,x_{2}=-1.9918\text{ and }x_{3}=1.1242$) and $3$ is, by
Caratheodory’s theorem, the most number of points for a two parameter binary
regression model. In this case, the three points are not equally weighted and
they have weights $0.2648,0.4289,0.3063$ respectively.

In addition, it is not necessary that a $D$-optimal design for this artificial
regression has $3$ support points. For example, if we let $\beta_{0}=0.1$,
$\beta_{1}=1$ and restrict $\eta$ to $[-10,10]$, then the $D$-optimal design
has only $2$ support points (right panel of Figure 5). We have
$x_{1}=-0.01861$ and $x_{2}=-10$.

Figure 5: Sensitivity functions for logit link with a erratic regression
function.

######  Example 2.5 (Student t).

For this problem,
$f(\eta)=\frac{\Gamma(\frac{k+1}{2})}{\sqrt{k\pi}\Gamma(\frac{k}{2})}\left(1+\frac{\eta^{2}}{k}\right)^{-\frac{k+1}{2}}$
and $\omega=f(\eta)^{2}/(F(\eta)S(\eta))$. The Figure 6 demonstrates the
sensitivity functions of two locally D-optimal designs with Student t link
(degrees of freedom is 2) and specified parameter values. The left panel has
design space $[-10,10]$ and the $D$-optimal design is
$x_{1}=-10,x_{2}=-0.5247$. The right panel has design space $[0,1]$ and the
$D$-optimal design is $x_{1}=1.7121,x_{2}=0$.

Figure 6: Sensitivity functions for Student t link.

###  2.2 Asymmetric Densities

If $f(\cdot)$ is asymmetric, then the resulting optimal design is not
symmetric in general. In this case, we have a system of non-linear WC
equations:

$\displaystyle\frac{2f^{\prime}(\eta_{1})}{f(\eta_{1})}-\frac{f(\eta_{1})}{F(\eta_{1})}+\frac{f(\eta_{1})}{S(\eta_{1})}+\frac{2}{\eta_{1}-\eta_{2}}$
$\displaystyle=0$
$\displaystyle\frac{2f^{\prime}(\eta_{2})}{f(\eta_{2})}-\frac{f(\eta_{2})}{F(\eta_{2})}+\frac{f(\eta_{2})}{S(\eta_{2})}+\frac{2}{\eta_{2}-\eta_{1}}$
$\displaystyle=0$

Define

$W(\eta)=\frac{2f^{\prime}(\eta)}{f(\eta)}-\frac{f(\eta)}{F(\eta)}+\frac{f(\eta)}{S(\eta)}$

and substitute it into the above system of non-linear equations:

$\displaystyle W(\eta_{1})$
$\displaystyle=-W\left(\frac{2}{W(\eta_{1})}+\eta_{1}\right)$
$\displaystyle\eta_{2}$ $\displaystyle=\eta_{1}+\frac{2}{W(\eta_{1})}$ (2.7)

Hence, in practice, we solve the first equation for $\eta_{1}$ and then plug-
in it to the second one to derive $\eta_{2}$.

######  Example 2.6 (Cox Regression).

For this problem, $f(\eta)=\exp\left(\eta-\exp(\eta)\right)$. Plug-in all
necessary elements and solving numerically, we obtain $\eta_{1}=+0.9796$ and
$\eta_{2}=-1.3378$.

The Figure 7 demonstrates the sensitivity functions of two locally D-optimal
designs with complementary log-log link and specified parameter values.

Figure 7: Sensitivity functions for complementary log-log link.

######  Example 2.7 (One-Hit Model).

The one-hit model is also known as exponential regression in toxicology
studies [Raz22] and the density is $f(\eta)=\exp\left(-\eta\right)$ for all
$\eta\geq 0$. The other terms associated with the one-hit model are
$f^{\prime}(\eta)=-\exp(-\eta),\pi=F(\eta)=1-\exp(-\eta)\text{ and
}S(\eta)=\exp(-\eta)$. However, the non-linear WC equation is numerically
unstable in this case. That is, if we let
$h(\eta)=W(\eta)+W(\frac{2}{W(\eta)}+\eta)$, then

$\lim_{\eta\downarrow 0}h(\eta)=-\infty$

To show this, we first write $W(\eta)$ and $h(\eta)$ explicitly, that is,

$\displaystyle W(\eta)$ $\displaystyle=-\frac{1}{1-\exp(-\eta)}$
$\displaystyle h(\eta)$
$\displaystyle=-\frac{1}{1-\exp(-\eta)}-\frac{1}{1-\exp(2-2\exp(-\eta)-\eta)}$

If $\eta\downarrow 0$, then both terms go to negative infinity. Next, if we
plug-in $\eta=1$, then $h(\eta)\approx 5.2263$. By continuity of $h(\eta)$,
there is at least one zero point between 0 and 1, and we plot the behavior of
$h(\eta)$ near $0$ below.

Figure 8: Behavior of $h(\eta)$ near 0.

As we can see from Figure 8, there are multiple zero points of $h(\eta)$ when
$\eta$ is close to zero. Therefore, as a take home message for practitioners,
it is better for us to enforce $\eta$ no less than a strictly positive value
$\eta_{\text{low}}$, say $\eta_{\text{low}}=0.5$. In other words, it means
that if the dose level is 0 ($x=0$), then the baseline response probability is

$\pi=1-\exp(-\eta_{\text{low}})=1-\exp(-0.5)=0.3935$

and the choice of $\eta_{\text{low}}$ should come from toxicologists.

The Figure 9 demonstrates the sensitivity functions of two locally D-optimal
designs with exponential link and specified parameter values. The results are
generated by PSO with 50 particles and 1,000 iterations, and the hyper-
parameters are $c_{1}=0.5,c_{2}=0.3$ and $w=0.9$. The left panel has design
space $\eta\in[0.5,1]$ and the optimal design is $x_{1}=0.5$, $x_{2}=1$. The
right panel has design space $\eta\in[1,\infty]$ and the optimal design is
$x_{1}=0$, $x_{2}=1.7978$.

Figure 9: Sensitivity functions for exponential link.

######  Example 2.8 (Weibull).

The Weibull distribution includes the exponential distribution as a special
case. The density of Weibull distribution is
$f(x)=\exp(-\beta_{0}-\beta_{1}x^{\alpha})$ where $\alpha$ is a positive
parameter. If $\alpha=1$, then it is the exponential density. Let
$\eta=\beta_{0}+\beta_{1}x^{\alpha}$ then the whole framework for the one-hit
model (exponential regression) can be copied almost completely. Hence, in
practice, we suggest practitioners to start with a small positive $\beta_{0}$
to avoid the numerical issue.

The Figure 10 demonstrates the sensitivity functions of two locally D-optimal
designs with Weibull link ($\alpha=3$) and specified parameter values. The
results are generated by PSO with 50 particles and 500 iterations, and the
hyper-parameters are $c_{1}=0.5,c_{2}=0.3$ and $w=0.9$. The left panel has
design space $\eta\in[0,3]$ and the optimal design is $x_{1}=0$,
$x_{2}=0.7686$. The right panel has design space $\eta\in[0,1]$ and the
optimal design is $x_{1}=0$, $x_{2}=1$.

Figure 10: Sensitivity functions for Weibull link.

##  3 Extension to Three-parameter Binary Regression

Let $\pi=c+(1-c)F(\eta)$ where $c$ is the background response probability
within 0 to 1 [Raz22]. Such model is particularly useful when we have a
potential baseline response probability that is strictly positive. Then the
parameter now becomes $\theta^{T}=(c,\beta)$ where $\beta$ is a
$d$-dimensional vector and in this section we always assume $d=2$. Given a
single design point $x$, the Fisher information matrix of $\theta$ is
($\mathcal{L}=\log p(y|x,\pi)$)

$\displaystyle\mathbb{E}\left(\frac{\partial\mathcal{L}}{\partial\theta}\right)\left(\frac{\partial\mathcal{L}}{\partial\theta^{T}}\right)$
$\displaystyle=\mathbb{E}\left(\frac{\partial\mathcal{L}}{\partial\pi}\frac{\partial\pi}{\partial\theta}\right)\left(\frac{\partial\mathcal{L}}{\partial\pi}\frac{\partial\pi}{\partial\theta^{T}}\right)$

The term $\frac{\partial\pi}{\partial\theta}$ is a $(d+1)$-vector and we have

$\displaystyle\frac{\partial\mathcal{L}}{\partial\pi}$
$\displaystyle=\frac{y-\pi}{\pi(1-\pi)}=\frac{y-F(\eta)}{F(\eta)S(\eta)}$
$\displaystyle\frac{\partial\pi}{\partial c}$
$\displaystyle=1-F(\eta)=S(\eta)$
$\displaystyle\frac{\partial\pi}{\partial\beta}$ $\displaystyle=(1-c)f(\eta)x$

Plug-in them into the Fisher information matrix, we have

$\displaystyle\mathbb{E}\left(\frac{\partial\mathcal{L}}{\partial\theta}\right)\left(\frac{\partial\mathcal{L}}{\partial\theta^{T}}\right)$
$\displaystyle=\mathbb{E}\left(\frac{y-F(\eta)}{F(\eta)S(\eta)}\right)^{2}\left(\begin{matrix}S(\eta)\\\
(1-c)f(\eta)x\end{matrix}\right)\left(\begin{matrix}S(\eta)\\\
(1-c)f(\eta)x\end{matrix}\right)^{T}$
$\displaystyle=\frac{1}{F(\eta)S(\eta)}\left(\begin{matrix}S(\eta)^{2}&(1-c)S(\eta)f(\eta)x^{T}\\\
(1-c)f(\eta)S(\eta)x&(1-c)^{2}f(\eta)^{2}xx^{T}\end{matrix}\right)$

Hence, for a $k$-point design, the information matrix is

$\displaystyle\mathfrak{M}(\xi)=\sum_{i=1}^{k}p_{i}\left(\begin{matrix}\frac{S(\eta_{i})}{F(\eta_{i})}&\frac{f(\eta_{i})}{F(\eta_{i})}(1-c)x_{i}^{T}\\\
\\\
\frac{f(\eta_{i})}{F(\eta_{i})}(1-c)x_{i}&\frac{f(\eta_{i})^{2}}{F(\eta_{i})S(\eta_{i})}(1-c)^{2}x_{i}x_{i}^{T}\end{matrix}\right)$

where $p_{i}$’s are weights summing to $1$. For example, suppose we are
interested in $3$-point $D$-optimal design, i.e., $k=3$, then by theorem 1 we
have $p_{1}=p_{2}=p_{3}=1/3$. The determinant of the resulting information
matrix has a neat representation and we state and prove it below.

######  Lemma 3.

The determinant of $\mathfrak{M}(\xi)$ for $k=3$ and $p_{i}=1/3,i=1,2,3$ has
the following representation

$\displaystyle\det\mathfrak{M}(\xi)\propto\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\right)^{3}\det\left(\text{Var}(\widetilde{x})\right)$
(3.1)

where $\widetilde{x}_{i}=\frac{f(\eta_{i})}{S(\eta_{i})}(1-c)x_{i}$ (a
$2$-dimensional vector) and $\text{Var}(\widetilde{x})$ is the $2\times 2$
variance-covariance matrix of $\widetilde{x}$ with respect to the tilted
probability measure [Dab19, Wai19] supported at
$(\widetilde{x}_{1},\widetilde{x}_{2},\widetilde{x}_{3})$ with probability
proportional to
$\left(\frac{S(\eta_{1})}{F(\eta_{1})},\frac{S(\eta_{2})}{F(\eta_{2})},\frac{S(\eta_{3})}{F(\eta_{3})}\right)$.
Therefore, we expression the determinant of a $3\times 3$ matrix in terms of a
$2\times 2$ matrix.

###### Proof.

By the block matrix determinant formula [BV04], we have

$\displaystyle\det\mathfrak{M}(\xi)\propto\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\right)^{-1}\det\left(\underbrace{\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\right)\left(\sum_{i=1}^{3}\frac{f(\eta_{i})^{2}(1-c)^{2}}{F(\eta_{i})S(\eta_{i})}x_{i}x_{i}^{T}\right)-\left(\sum_{i=1}^{3}\frac{f(\eta_{i})}{F(\eta_{i})}(1-c)x_{i}\right)^{\otimes
2}}_{(\Delta)}\right)$

where $y^{\otimes 2}=yy^{T}$ and we write $\otimes 2$ for convenience. Next,
we re-write $(\Delta)$ as

$\displaystyle(\Delta)$
$\displaystyle=\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\right)\left(\sum_{i=1}^{3}\frac{S(\eta_{i})f(\eta_{i})^{2}}{F(\eta_{i})S(\eta_{i})^{2}}(1-c)^{2}x_{i}x_{i}^{T}\right)-\left(\sum_{i=1}^{3}\frac{S(\eta_{i})f(\eta_{i})}{F(\eta_{i})S(\eta_{i})}(1-c)x_{i}\right)^{\otimes
2}$
$\displaystyle=\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\right)\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\left(\frac{f(\eta_{i})}{S(\eta_{i})}(1-c)x_{i}\right)^{\otimes
2}\right)-\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\left(\frac{f(\eta_{i})}{S(\eta_{i})}(1-c)x_{i}\right)\right)^{\otimes
2}$
$\displaystyle=\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\right)^{2}\left[\left(\sum_{i=1}^{3}\frac{S(\eta_{i})/F(\eta_{i})}{\sum_{j=1}^{3}S(\eta_{j})/F(\eta_{j})}\widetilde{x}^{\otimes
2}\right)-\left(\sum_{i=1}^{3}\frac{S(\eta_{i})/F(\eta_{i})}{\sum_{j=1}^{3}S(\eta_{j})/F(\eta_{j})}\widetilde{x}\right)^{\otimes
2}\right]$
$\displaystyle=\left(\sum_{i=1}^{3}\frac{S(\eta_{i})}{F(\eta_{i})}\right)^{2}\text{Var}(\widetilde{x})$

where $\text{Var}(\widetilde{x})$ is the variance-covariance matrix with
respect to the tilted probability measure [Wai19] supported at
$(\widetilde{x}_{1},\widetilde{x}_{2},\widetilde{x}_{3})$ with probability
$(S(\eta_{1})/F(\eta_{1}),S(\eta_{2})/F(\eta_{2}),S(\eta_{3})/F(\eta_{3}))$
and $\widetilde{x}_{i}=\frac{f(\eta_{i})}{S(\eta_{i})}x_{i}$. Substituting the
new expression of $(\Delta)$ into $\det\mathfrak{M}(\xi)$, we have the neat
formula 3.1 for the determinant of $\mathfrak{M}(\xi)$. ∎

Unfortunately, an analytical maximizer (or a similar WC equation for formula
3.1) is difficult to obtain and not useful in practice. Hence, we suggest to
apply formula 3.1 and optimization tools such as PSO to solve for the optimal
design in practice. Further, if the interest is in estimating $c$ in
particular, then the corresponding $D_{s}$-optimality is the Schur complement
[BV04] of the first element of $\mathfrak{M}(\xi)$, that is, we seek a design
$\xi$ that maximizes

$\left(\sum_{i=1}^{3}p_{i}\frac{S(\eta_{i})}{F(\eta_{i})}\right)-\left(\sum_{i=1}^{3}p_{i}\frac{f(\eta_{i})}{F(\eta_{i})}x_{i}^{T}\right)\left(\sum_{i=1}^{3}\frac{f(\eta_{i})^{2}}{F(\eta_{i})S(\eta_{i})}x_{i}x_{i}^{T}\right)^{-1}\left(\sum_{i=1}^{3}p_{i}x_{i}\frac{f(\eta_{i})}{F(\eta_{i})}\right)$

It is worth noting that $D_{s}$-optimality may lead to singular information
matrix $\mathfrak{M}(\xi)$ and in this case, only certain linear combinations
of the parameters are estimable [Páz86, CW21, Sil13]. Further, the sensitivity
function for model $\pi=c+(1-c)F(\eta)$ is different from the one provided in
1.1 and by the linearity technique [Páz86, Sil13], we define

$\displaystyle\psi(\mathbf{x},\theta)$
$\displaystyle=\mathbf{u}^{T}\mathfrak{M}(\xi)^{-1}\mathbf{u}-d$ (3.2)
$\displaystyle\mathbf{u}^{T}$
$\displaystyle=\left(\begin{matrix}S(\eta),&(1-c)\sqrt{\omega(\mathbf{x})}\mathbf{x}\end{matrix}\right)$

where $\omega(\mathbf{x})$ is, again, a weight depending on the link
$\eta(\mathbf{x})$. The global optimality of $\xi$ can be verified using the
analogue of the equivalence theorem 1.1.

######  Example 3.1.

(Logit) We provide two examples regarding the logit link here. In Figure 11,
the left panel shows the sensitivity function when $\beta_{0}=1,\beta_{1}=0.5$
and $c=0.1$ and the constraint is $[0,1]$. The locally $D$-optimal design is
$\xi^{*}=[0,0.4643,1]$ with equal weights. The right panel shows the
sensitivity function when $\beta_{0}=0,\beta_{1}=1$ and $c=0.2$ and the
constraint is $[-10,10]$. The locally $D$-optimal design is
$\xi^{*}=[-10,-1.4555,1.6137]$ with equal weights. Interestingly, both designs
include the lower bound (0 and -10) as a design point, this is probably due to
the accurate estimation of the baseline parameter $c$.

Figure 11: Sensitivity functions for logit link.

##  4 Applications to Toxicology Studies

In this section, we apply the developed WC theory to a dataset which comes
from toxicology studies using sea urchins [CCHW22] (Figure 14). There are two
endpoints (failure types): EDA/D and Radial:Ab and we use the second endpoint
for illustration. The concentration level for the second endpoint is within 0
to 450 $\mu M$, and we re-scale it to $[0,0.45]$ by dividing 1000. We run two
binary regression models using logit and complementary log-log (Cox
regression) link functions respectively. The results are generated by
‘gtsummary‘ package in R [SWC+21] and given in Table 2.

Table 2: Binary regression with two different link functions using sea urchin
data Logit link
---
Characteristic | Estimation | 95% CI |  $p$-value
$\beta_{0}$ | -4.5 | (-4.7,-4.4) | < 0.001
$\beta_{1}$ | 20 | (19,21) | < 0.001
Cox regression
Characteristic | Estimation | 95% CI |  $p$-value
$\beta_{0}$ | -3.7 | (-3.8, -3.6) | < 0.001
$\beta_{1}$ | 14 | (13, 14) | < 0.001

The fitted dose-response curve (in this case, concentration-response curve) is
given in Figure 12: the orange and dodgerblue curves correspond to Cox
regression and logit link respectively. The dots represent the true
observations from Table 14 with concetration level greater than 450 removed.

Figure 12: The fitted concentration-response curves.

Hence, by the WC equation 2.2 for two-parameter binary regression, the
$D$-optimal designs are

$\displaystyle\xi_{\text{logit}}$
$\displaystyle=\left(\begin{matrix}0.1478&0.3022\\\
0.5&0.5\end{matrix}\right)$ $\displaystyle\xi_{\text{Cox}}$
$\displaystyle=\left(\begin{matrix}0.1687&0.3343\\\
0.5&0.5\end{matrix}\right)$

and the sensitivity functions are given in Figure 13 (left panel: logit link;
right panel: Cox regression).

Figure 13: Sensitivity functions for the sea urchin data.

Multiplying by 1000, the resulting $D$-optimal designs at the original scale
are

$\displaystyle\xi^{*}_{\text{logit}}$
$\displaystyle=\left(\begin{matrix}147.8&302.2\\\ 0.5&0.5\end{matrix}\right)$
(4.1) $\displaystyle\xi^{*}_{\text{Cox}}$
$\displaystyle=\left(\begin{matrix}168.7&334.3\\\ 0.5&0.5\end{matrix}\right)$
(4.2)

Comparing them with the original design given in [CCHW22]:

$\displaystyle\xi^{*}_{\text{original}}$
$\displaystyle=\left(\begin{matrix}0&100&125&150&175&180&200&225&300&450\\\
0.254&0.148&0.0129&0.169&0.0263&0.0338&0.128&0.0370&0.155&0.0360\end{matrix}\right)$
(4.3)

we find that the $D$-optimal design reduces the number of required
concentration levels significantly.

##  5 Discussion

In this paper, we have systematically discussed the $D$-optimal design in two-
parameter binary regression model with various link functions. As an
extension, we have provided a analytical formula of determinant for handling
three-parameter binary regression. PSO, a type of metaheuristics [CCW22], is
applied to derive $D$-optimal designs when analytical solution is not
available by the WC equation. For four- and more parameter designs, the
analytical formula is neither easy to derive nor useful in practice.
Therefore, we suggest practitioners to use PSO to find optimal designs instead
of working with the analytical solutions. We have also provided a real data
example in toxicology studies illustrating the use of the developed results.

To handle more complicated situations (i.e., multi-hills of dose-response
curve), we can further extend the function $\pi=F(\eta)$ to multistage
[Car81], multi-hit [RVR81] and dichotomous hill models [GRM+16]. We leave
these as future work and emphasis that metaheuristics can be applied to these
models conveniently compared with analytical solutions.

Figure 14: The dataset from toxicology studies using sea urchins.

## References

  * [ADT07] Anthony Atkinson, Alexander Donev, and Randall Tobias.  Optimum experimental designs, with SAS, volume 34.  Oxford University Press, 2007.
  * [ASE+17] Joseph P Antonios, Horacio Soto, Richard G Everson, Diana Moughon, Joey R Orpilla, Namjo P Shin, Shaina Sedighim, Janet Treger, Sylvia Odesa, Alexander Tucker, et al.  Immunosuppressive tumor-infiltrating myeloid cells mediate adaptive immune resistance via a pd-1/pd-l1 mechanism in glioblastoma.  Neuro-oncology, 19(6):796–807, 2017.
  * [BV04] Stephen P Boyd and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [BZWW06] InYoung Baek, Wei Zhu, Xiangfeng Wu, and Weng Kee Wong.  Bayesian optimal designs for a quantal dose-response study with potentially missing observations.  Journal of Biopharmaceutical Statistics, 16(5):679–693, 2006.
  * [Car81] FW Carlborg.  Multi-stage dose-response models in carcinogenesis.  Food and cosmetics toxicology, 19:361–365, 1981.
  * [CCHW22] Michael D Collins, Elvis Han Cui, Seung Won Hyun, and Weng Kee Wong.  A model-based approach to designing developmental toxicology experiments using sea urchin embryos.  Archives of toxicology, pages 1–14, 2022.
  * [CCW22] Ping-Yang Chen, Ray-Bing Chen, and Weng Kee Wong.  Particle swarm optimization for searching efficient experimental designs: A review.  Wiley Interdisciplinary Reviews: Computational Statistics, page e1578, 2022.
  * [CW21] Elvis Cui and Weng Kee Wong.  Lecture Notes for Biostat 250AB.  Unpublished Manuscript at UCLA, 2021.
  * [Dab19] Dorota M Dabrowska.  Elements of real analysis and advanced probability, volume 1.  Lecture Notes at UCLA, 2019.
  * [DB18] Annette J Dobson and Adrian G Barnett.  An introduction to generalized linear models.  Chapman and Hall/CRC, 2018.
  * [Fed72] VV Federov.  Theory of Optimal Experiments, translated and edited by WJ.  JSTOR, 1972.
  * [Fed10] Valerii Fedorov.  Optimal experimental design.  Wiley Interdisciplinary Reviews: Computational Statistics, 2(5):581–589, 2010.
  * [Gil06] Jim Giles.  Animal experiments under fire for poor design.  Nature, 444(7122):981–982, 2006.
  * [GRM+16] Bradford W Gutting, Andrey Rukhin, David Marchette, Ryan S Mackie, and Brandolyn Thran.  Dose-response modeling for inhalational anthrax in rabbits following single or multiple exposures.  Risk Analysis, 36(11):2031–2038, 2016.
  * [HKO07] Linda M Haines, Gaëtan Kabera, and Timothy E O’Brien.  D-optimal designs for logistic regression in two variables.  In mODa 8-Advances in Model-Oriented Design and Analysis, pages 91–98. Springer, 2007.
  * [HLP52] Godfrey Harold Hardy, John Edensor Littlewood, and George Pólya.  Inequalities.  Cambridge university press, 1952.
  * [JM13] Katarzyna Jóźwiak and Mirjam Moerbeek.  Podse: A computer program for optimal design of trials with discrete-time survival endpoints.  Computer methods and programs in biomedicine, 111(1):115–127, 2013\.
  * [KH12] M Gaëtan Kabera and Linda M Haines.  A note on the construction of locally d-and ds-optimal designs for the binary logistic model with several explanatory variables.  Statistics & Probability Letters, 82(5):865–870, 2012.
  * [Kie74] Jack Kiefer.  General equivalence theory for optimum designs (approximate theory).  The annals of Statistics, pages 849–879, 1974.
  * [KW00] Joy King and Weng-Kee Wong.  Minimax d-optimal designs for the logistic model.  Biometrics, 56(4):1263–1267, 2000.
  * [Mir18] Lester James Miranda.  Pyswarms: a research toolkit for particle swarm optimization in python.  Journal of Open Source Software, 3(21):433, 2018.
  * [Páz86] Andrej Pázman.  Foundations of optimum experimental design, volume 14.  Springer, 1986.
  * [Raz22] Mehdi Razzaghi.  Statistical Models in Toxicology.  CRC, 2022.
  * [RVR81] Kamta Rai and John Van Ryzin.  A generalized multihit dose-response model for low-dose extrapolation.  Biometrics, pages 341–352, 1981.
  * [Sil13] Samuel Silvey.  Optimal design: an introduction to the theory for parameter estimation, volume 1.  Springer Science & Business Media, 2013.
  * [Smi18] Kirstine Smith.  On the standard deviations of adjusted and interpolated values of an observed polynomial function and its constants and the guidance they give towards a proper choice of the distribution of observations.  Biometrika, 12(1/2):1–85, 1918.
  * [SRW20] Oleksandr Sverdlov, Yevgen Ryeznik, and Weng Kee Wong.  On optimal designs for clinical trials: an updated review.  Journal of Statistical Theory and Practice, 14(1):1–29, 2020.
  * [SSQ+06] Marcio Schwaab, Fabrício M Silva, Christian A Queipo, Amaro G Barreto Jr, Márcio Nele, and José Carlos Pinto.  A new approach for sequential experimental design for model discrimination.  Chemical engineering science, 61(17):5791–5806, 2006.
  * [SWC+21] Daniel D Sjoberg, Karissa Whiting, Michael Curry, Jessica A Lavery, and Joseph Larmarange.  Reproducible summary tables with the gtsummary package.  R Journal, 13(1), 2021.
  * [Wai19] Martin J Wainwright.  High-dimensional statistics: A non-asymptotic viewpoint, volume 48.  Cambridge University Press, 2019.
  * [WL96] Weng Kee Wong and Peter A Lachenbruch.  Designing studies for dose response.  Statistics in Medicine, 15(4):343–359, 1996.
  * [Won21] Weng Kee Wong.  Lecture Notes for Biostat 279.  Unpublished Manuscript at UCLA, 2021.
  * [ZAKW98] Wei Zhu, Hongshik Ahn, and Weng Kee Wong.  Multiple-objective optimal designs for the logit model.  Communications in Statistics-Theory and Methods, 27(6):1581–1592, 1998.
  * [ZW00] Wei Zhu and Weng Kee Wong*.  Multiple-objective designs in a dose-response experiment.  Journal of Biopharmaceutical Statistics, 10(1):1–14, 2000.
  * [ZWY21] Xiao-Dong Zhou, Yun-Juan Wang, and Rong-Xian Yue.  Optimal designs for discrete-time survival models with random effects.  Lifetime Data Analysis, 27(2):300–332, 2021.