Optimal experimental designs for treatment contrasts in heteroscedastic models with covariates

TL;DR

This paper develops optimal experimental design strategies for treatment contrast estimation in heteroscedastic models with covariates, extending known optimality results and proposing a linear programming approach for practical design implementation.

Contribution

It introduces new optimal design methods for heteroscedastic models with covariates, including eigenvalue-based criteria and a sparsification technique for practical application.

Findings

Derived A- and E-optimal product designs for heteroscedastic models.

Proposed a linear programming method for sparse, practical designs.

Validated methods through illustrative examples.

Abstract

In clinical trials, the response of a given subject often depends on the selected treatment as well as on some covariates. We study optimal approximate designs of experiments in the models with treatment and covariate effects. We allow for the variances of the responses to depend on the chosen treatments, which introduces heteroscedasticity into the models. For estimating systems of treatment contrasts and linear functions of the covariates, we extend known results on D-optimality of product designs by providing product designs that are optimal with respect to general eigenvalue-based criteria. In particular, A- and E-optimal product designs are obtained. We then formulate a method based on linear programming for constructing optimal designs with smaller supports from the optimal product designs. The sparser designs can be more easily converted to practically applicable exact designs. The provided results and the proposed sparsification method are demonstrated on some examples.

Tables

Table 1. Table 1 : Optimal product design for Example 1 . The indices k = 1 , … , 8 𝑘 1 … 8 k=1,\ldots,8 denote the active covariate values 𝐳 = ( − 1 , − 1 , − 1 ) T , ( − 1 , − 1 , 1 ) T , … , ( 1 , 1 , 1 ) T 𝐳 superscript 1 1 1 𝑇 superscript 1 1 1 𝑇 … superscript 1 1 1 𝑇 \mathbf{z}=(-1,-1,-1)^{T},(-1,-1,1)^{T},\ldots,(1,1,1)^{T} arranged in the lexicographical order; the treatments are denoted by i = 1 , 2 , 3 𝑖 1 2 3 i=1,2,3 .

$i\backslash k$	1	2	3	4	5	6	7	8
1	0.0295	0.0295	0.0295	0.0295	0.0295	0.0295	0.0295	0.0295
2	0.0477	0.0477	0.0477	0.0477	0.0477	0.0477	0.0477	0.0477
3	0.0477	0.0477	0.0477	0.0477	0.0477	0.0477	0.0477	0.0477

Table 2. Table 2 : Optimal non-product design for Example 1 ; the notation is the same as in Table 1 .

$i\backslash k$	1	2	3	4	5	6	7	8
1	0.0378	0	0.0212	0.0591	0.0212	0.0591	0.0378	0
2	0	0.1909	0	0	0	0	0.1909	0
3	0.1909	0	0	0	0	0	0	0.1909

Table 3. Table 3 : Exact design for n = 48 𝑛 48 n=48 trials constructed from the optimal non-product design given in Table 2 .

$i\backslash k$	1	2	3	4	5	6	7	8
1	2	0	1	3	1	3	2	0
2	0	9	0	0	0	0	9	0
3	9	0	0	0	0	0	0	9

Table 4. Table 4 : Optimal non-product design for Example 2 . The table contains the values ξ ​ ( i , j , k ) 𝜉 𝑖 𝑗 𝑘 \xi(i,j,k) , where i 𝑖 i is the treatment, j 𝑗 j is the first covariate level (expressed in rows) and k 𝑘 k is the second covariate level (expressed in columns). For example, the value 0.0121 0.0121 0.0121 in row 2 and column 3 of the first “block” of the table means that ξ ​ ( 1 , 2 , 3 ) = 0.0121 𝜉 1 2 3 0.0121 \xi(1,2,3)=0.0121 .

$\xi (1,j,k)$	1	2	3	4	5
1	0.0303	0.0121	0.0303	0.0121	0.0061
2	0.0061	0.0303	0.0121	0.0242	0.0182
3	0.0182	0.0121	0.0121	0.0182	0.0303
$\xi (2,j,k)$	1	2	3	4	5
1	0	0	0	0.0485	0.0727
2	0.0242	0	0.0727	0.0242	0
3	0.0485	0.0727	0	0	0
$\xi (3,j,k)$	1	2	3	4	5
1	0	0.0727	0	0.0242	0.0242
2	0.0727	0	0	0	0.0485
3	0	0	0.0727	0.0485	0

Table 5. Table 5 : Exact design for n = 40 𝑛 40 n=40 trials constructed from the optimal non-product design given in Table 4 .

$\xi (1,j,k)$	1	2	3	4	5
1	1	1	1	1	1
2	1	1	1	1	1
3	1	1	1	1	1
$\xi (2,j,k)$	1	2	3	4	5
1	0	0	0	2	2
2	1	0	2	1	0
3	2	2	0	0	0
$\xi (3,j,k)$	1	2	3	4	5
1	0	2	0	1	1
2	3	0	0	0	2
3	0	0	2	2	0

Table 6. Table 6 : Optimal non-product design for Example 3 . i 𝑖 i – treatment, k 𝑘 k – time moment.

$i\backslash k$	1	2	3	4	5	6
1	0.0880	0.1667	0.0552	0.0166	0	0.1048
2	0.0787	0	0	0	0.1667	0.0036
3	0	0	0	0.1501	0	0.0260
4	0	0	0.1115	0	0	0.0323

Table 7. Table 7 : Exact design constructed from the optimal non-product design given in Table 6 .

$i\backslash k$	1	2	3	4	5	6
1	1	1	0	0	0	1
2	0	0	0	0	1	0
3	0	0	0	1	0	0
4	0	0	1	0	0	0

Equations

$$y(i,k)={\tau }_i+\mu +{\mathbf{g}}^T(k)\beta +\epsilon (i,k),$$

$$\mathbf{M}(\xi )=\left[\begin{array}{cc}{\mathbf{M}}_{11}(\xi )& {\mathbf{M}}_{12}(\xi )\\ {\mathbf{M}}_{12}^T(\xi )& {\mathbf{M}}_{22}(\xi )\end{array}\right],$$

$$\mathbf{A}=\left[\begin{array}{cc}{\mathbf{Q}}_1& 0_{v_1\times s_2}\\ 0_{(v_2+1)\times s_1}& {\mathbf{Q}}_2\end{array}\right]=\mathrm{diag}({\mathbf{Q}}_1,{\mathbf{Q}}_2)$$

$${\Phi }_p(\mathbf{N})=\{\begin{array}{ll}(\det (\mathbf{N}){)}^{1/s},& p=0;\\ (\frac{1}{s}\mathrm{tr}({\mathbf{N}}^p){)}^{1/p},& p\in (-\infty ,0);\\ {\gamma }_{\min }(\mathbf{N}),& p=-\infty ,\end{array}$$

$$y(i)={\tau }_i+\epsilon (i),$$

$$y(k)=\mu +{\mathbf{g}}^T(k)\beta +{\epsilon }_k.$$

$${\mathbf{M}}_2(\alpha )=\sum _{k=1}^d{\alpha }_k\mathbf{h}(k){\mathbf{h}}^T(k)=\left[\begin{array}{cc}1& {\sum }_k{\alpha }_k{\mathbf{g}}^T(k)\\ {\sum }_k{\alpha }_k\mathbf{g}(k)& {\sum }_k{\alpha }_k\mathbf{g}(k){\mathbf{g}}^T(k)\end{array}\right].$$

$$\mathbf{B}=\left[\begin{array}{cc}{\mathbf{B}}_{11}& {\mathbf{B}}_{12}\\ {\mathbf{B}}_{12}^T& {\mathbf{B}}_{22}\end{array}\right]$$

$$\mathbf{G}=\left[\begin{array}{cc}{\mathbf{B}}_{11}^{-1}+{\mathbf{B}}_{11}^{-1}{\mathbf{B}}_{12}{\mathbf{B}}_{\mathbf{S}}^{-}{\mathbf{B}}_{12}^T{\mathbf{B}}_{11}^{-1}& -{\mathbf{B}}_{11}^{-1}{\mathbf{B}}_{12}{\mathbf{B}}_{\mathbf{S}}^{-}\\ -{\mathbf{B}}_{\mathbf{S}}^{-}{\mathbf{B}}_{12}^T{\mathbf{B}}_{11}^{-1}& {\mathbf{B}}_{\mathbf{S}}^{-}\end{array}\right],$$

$${\mathbf{N}}_{\mathbf{A}}(w\otimes \alpha )=\left[\begin{array}{cc}{\mathbf{N}}_{{\mathbf{Q}}_1}(w)& 0\\ 0& ({\sum }_i{\lambda }_i{w}_i){\mathbf{N}}_{\mathbf{K}}(\alpha )\end{array}\right].$$

$${\mathbf{S}}_1(\xi )$$

$$\displaystyle=\Big{(}\sum_{i=1}^{v_{1}}\lambda_{i}w_{i}\Big{)}\Big{(}\sum_{k=1}^{d}\alpha_{k}\mathbf{h}(k)\mathbf{h}^{T}(k)-(\sum_{k=1}^{d}\alpha_{k}\mathbf{h}(k))(\sum_{k=1}^{d}\alpha_{k}\mathbf{h}^{T}(k))\Big{)}$$

$$\displaystyle=\Big{(}\sum_{i=1}^{v_{1}}\lambda_{i}w_{i}\Big{)}\begin{bmatrix}0&\mathbf{0}^{T}\\ \mathbf{0}&\mathbf{S}(\alpha)\end{bmatrix}.$$

$${\mathbf{M}}_{11}^{-1}(\xi ){\mathbf{M}}_{12}(\xi )$$

$$=1_{v_1}\sum _{k=1}^d{\alpha }_k{\mathbf{h}}^T(k).$$

$${\mathbf{N}}_{\mathbf{A}}^{-1}(\xi )=\left[\begin{array}{cc}{\mathbf{Q}}_1^T{\mathbf{M}}_1^{-1}(w){\mathbf{Q}}_1& 0\\ 0& {\mathbf{Q}}_2^T{\mathbf{S}}_1^{-}(\xi ){\mathbf{Q}}_2\end{array}\right],$$

$$\tilde{\mathbf{M}}(\xi )=\left[\begin{array}{cc}{\mathbf{M}}_{11}(\xi )& -{\mathbf{M}}_{12}(\xi )\\ -{\mathbf{M}}_{12}^T(\xi )& {\mathbf{M}}_{22}(\xi )\end{array}\right].$$

$${\tilde{\mathbf{N}}}_{\mathbf{A}}(\xi )=\left[\begin{array}{cc}-\mathbf{I}& 0\\ 0& \mathbf{I}\end{array}\right]{\mathbf{N}}_{\mathbf{A}}(\xi )\left[\begin{array}{cc}-\mathbf{I}& 0\\ 0& \mathbf{I}\end{array}\right].$$

$$\Phi (\frac{1}{2}({\mathbf{N}}_{\mathbf{A}}(\xi )+{\tilde{\mathbf{N}}}_{\mathbf{A}}(\xi )))\ge \frac{1}{2}\Phi ({\mathbf{N}}_{\mathbf{A}}(\xi ))+\frac{1}{2}\Phi ({\tilde{\mathbf{N}}}_{\mathbf{A}}(\xi ))=\Phi ({\mathbf{N}}_{\mathbf{A}}(\xi )).$$

$$\det ({\mathbf{N}}_{\mathbf{A}}(w\otimes \alpha ))=(\sum _{i=1}^{v_1}{\lambda }_i{w}_i{)}^{s_2}\det ({\mathbf{N}}_{{\mathbf{Q}}_1}(w))\det ({\mathbf{N}}_{\mathbf{K}}(\alpha )).$$

$$\Phi_{p}(\mathbf{N}_{\mathbf{A}}(w\otimes\alpha))=\Big{(}s^{-1}\big{(}\mathrm{tr}(\mathbf{N}_{\mathbf{Q}_{1}}^{p}(w))+(\sum_{i}\lambda_{i}w_{i})^{p}\mathrm{tr}(\mathbf{N}_{\mathbf{K}}^{p}(\alpha))\big{)}\Big{)}^{1/p}.$$

$$\mathbf{G}=\mathrm{diag}({\mathbf{M}}_1^{-1}(w^{\ast }),(\sum _i{\lambda }_i{w}_i^{\ast }{)}^{-1}{\mathbf{M}}_2^{-}({\alpha }^{\ast }))$$

$$w_i=w_i^{\ast },i=1,\dots ,v_1$$

$$\left[\frac{1}{w_1^{\ast }}\sum _{k=1}^d\xi (1,k)\mathbf{g}(k),\dots ,\frac{1}{w_v^{\ast }}\sum _{k=1}^d\xi (v_1,k)\mathbf{g}(k)\right]{\mathbf{Q}}_1=0_{v_2\times s_1},$$

$$\left[\begin{array}{c}{\sum }_k(\xi (1,k)-w_1^{\ast }{\alpha }_k^{\ast }){\mathbf{g}}^T(k)\\ ⋮\\ {\sum }_k(\xi (v_1,k)-w_{v_1}^{\ast }{\alpha }_k^{\ast }){\mathbf{g}}^T(k)\end{array}\right]{\mathbf{S}}^{-}({\alpha }^{\ast })\mathbf{K}=0_{v_1\times s_2},$$

$$\Big{(}\sum_{i=1}^{v_{1}}\lambda_{i}w_{i}^{*}\Big{)}^{-1}\Big{(}(\sum_{i,k}\lambda_{i}\xi(i,k)\mathbf{g}(k)\mathbf{g}^{T}(k))-(\sum_{i,k}\lambda_{i}\xi(i,k)\mathbf{g}(k))\sum_{k=1}^{v_{2}}\alpha^{*}_{k}\mathbf{g}^{T}(k)\Big{)}\mathbf{S}^{-}(\alpha^{*})\mathbf{K}=\mathbf{K},$$

$$\mathbf{G}=\mathrm{diag}({\mathbf{M}}_1^{-1}(w^{\ast }),(\sum _i{\lambda }_i{w}_i^{\ast }{)}^{-1}{\mathbf{M}}_2^{-}({\alpha }^{\ast })).$$

$${\mathbf{M}}_2^{-}({\alpha }^{\ast })=\left[\begin{array}{cc}1+({\sum }_k{\alpha }_k^{\ast }{\mathbf{g}}^T(k)){\mathbf{S}}^{-}({\alpha }^{\ast })({\sum }_k{\alpha }_k^{\ast }{\mathbf{g}}^T(k))& -{\sum }_k{\alpha }_k^{\ast }{\mathbf{g}}^T(k){\mathbf{S}}^{-}({\alpha }^{\ast })\\ -{\mathbf{S}}^{-}({\alpha }^{\ast }){\sum }_k{\alpha }_k^{\ast }\mathbf{g}(k)& {\mathbf{S}}^{-}({\alpha }^{\ast })\end{array}\right].$$

$${\mathbf{M}}_{11}(\xi ){\mathbf{M}}_1^{-1}(w^{\ast }){\mathbf{Q}}_1={\mathbf{Q}}_1$$

$${\mathbf{M}}_{12}(\xi ){\mathbf{M}}_2^{-}({\alpha }^{\ast }){\mathbf{Q}}_2=0$$

$${\mathbf{M}}_{12}^T(\xi ){\mathbf{M}}_1^{-1}(w^{\ast }){\mathbf{Q}}_1=0$$

$$(\sum _{i=1}^{v_1}{\lambda }_i{w}_i^{\ast }{)}^{-1}{\mathbf{M}}_{22}(\xi ){\mathbf{M}}_2^{-}({\alpha }^{\ast }){\mathbf{Q}}_2={\mathbf{Q}}_2$$

Full text

1 Introduction

In the experiments performed to estimate the effects of a selected set of treatments, it is common that the responses are also affected by some other experimental conditions (the covariates – e.g., time trend, block effects, the ages or the gender of the subjects in a clinical trial). The effects of the covariates are traditionally considered to be nuisance effects in the experimental design literature (e.g., Cox (1951), Majumdar and Notz (1983), Atkinson and Donev (1996), Jacroux et al. (1997), Rosa and Harman (2016)). Recently, Atkinson (2015) suggested that particularly with the growing importance of personalized medicine, the covariate effects may also be of prominent interest.

Moreover, the interest often lies in some functions of the model parameters rather than in the parameters themselves. The comparisons of test treatments with a control or with the placebo is a common example, especially in clinical trials. The design and analysis are sometimes further complicated by the presence of heteroscedasticity – the responses under different treatments may have different variances. In the present paper, we study optimal approximate designs in heteroscedastic models with treatment effects and covariate effects where there is interest in a set of treatment contrasts and in a set of linear combinations of the covariate effects.

For the abovementioned settings, Atkinson (2015) studied $D$-optimality in a model with only two treatments. Wang and Ai (2016) extended his results by providing $D$-optimal product designs for arbitrary numbers of treatments. The $D$-optimality is the most popular criterion in optimal design literature and is particularly nice to work with analytically. However, if the interest lies in a set of treatment contrasts (e.g., in the test treatment-control comparisons, which are a natural choice when placebo is included in the clinical trial), $D$-optimality tends to ignore the special interest in the chosen contrasts; i.e., this criterion tends to select designs that do not provide more information on the contrasts of interest compared to other treatment contrasts. As such, $D$-optimality is generally not recommended for such experimental interests (cf. Hedayat et al. (1988), Morgan and Wang (2010)). This also corresponds to the observation by Wang and Ai (2016) that if the experimental objective is to estimate the test treatment-control comparisons and all covariate effects, then the $D$-optimal designs are the same as if there was interest in all covariate effects and a uniform interest in all the treatments. That is, in such a case, the $D$-optimality does not place any special emphasis on the test treatment-control comparisons.

In contrast, the $A$-optimality criterion possesses a natural statistical interpretation for the considered settings – it minimizes the average variance for the linear functions of interest. Hence, $A$-optimality is very popular for the comparisons with the control (e.g., see Hedayat et al. (1988)). Recently, $E$-optimality was also argued to be meaningful for estimating treatment contrasts (e.g., Morgan and Wang (2011), Rosa (2018)). In this paper, we therefore extend the results by Wang and Ai (2016) to other eigenvalue-based optimality criteria (see Section 2). In particular, we obtain $\Phi_{p}$-optimal product designs for the other Kiefer’s $\Phi_{p}$-optimality criteria besides $D$-optimality, including $A$- and $E$-optimality.

The observation that the product designs are generally optimal for multi-factor models (like the treatment-covariate one) is extensively used in the literature; e.g., see Schwabe and Wierich (1995), Schwabe (1996), Rodriguez and Ortiz (2005), Graßhoff et al. (2007). One drawback of the optimal product designs is that they have large supports, and therefore it is sometimes difficult to construct designs for the actual experiments from the product designs by rounding procedures111The product designs, like other approximate designs, specify only proportions of trials to be performed for the particular experimental conditions, as will be formalized later. To obtain actual integer numbers of trials, some rounding is generally performed. (e.g., the well-known efficient rounding procedure by Pukelsheim and Rieder (1992)). However, in Section 3, we provide an entire class of $\Phi$-optimal designs characterized by linear constraints. This allows us to formulate a linear programming method for constructing optimal designs with smaller supports from the optimal product designs. Similar approach was employed by Rosa and Harman (2016) in a homoscedastic model, where the covariate effects were considered to be nuisance parameters.

The application of the theoretical results, and in particular the construction of the optimal designs with sparser supports and their usefulness in obtaining efficient exact designs of experiments is demonstrated on some examples in Section 4.

1.1 Notation

By $\mathbf{1}_{n}$ and $\mathbf{0}_{n}$, we denote the vector of ones and the vector of zeros in $\mathbb{R}^{n}$, respectively. The vector $\mathbf{e}_{i}$ has 1 on the $i$th position and zeros elsewhere. The identity matrix is denoted by $\mathbf{I}_{n}$ and the $m\times n$ matrix of zeros is denoted by $\mathbf{0}_{m\times n}$. For brevity, we sometimes omit the subscripts expressing the dimensions. The expression $\mathrm{diag}(a_{1},\ldots,a_{n})$, where $a_{i}\in\mathbb{R}$, denotes the diagonal matrix with $a_{1},\ldots,a_{n}$ on diagonal. If $\mathbf{A}$ and $\mathbf{B}$ are matrices, $\mathrm{diag}(\mathbf{A},\mathbf{B})$ denotes the corresponding block diagonal matrix. The smallest eigenvalue of a nonnegative definite matrix $\mathbf{A}$ is denoted by $\gamma_{\min}(\mathbf{A})$. Given a matrix $\mathbf{A}$, the symbol $\mathbf{A}^{-}$ denotes a generalized inverse of $\mathbf{A}$, and $\mathbf{A}^{+}$ is the Moore-Penrose pseudoinverse of $\mathbf{A}$.

1.2 The model

Consider the model

[TABLE]

where $i\in\{1,\ldots,v_{1}\}$ is the chosen treatment and $k\in\mathfrak{T}$, $|\mathfrak{T}|=d$, represents the chosen covariates. The vector $\tau=(\tau_{1},\ldots,\tau_{v_{1}})$ represents the treatment effects, $\mu$ is the constant term, $\beta=(\beta_{1},\ldots,\beta_{v_{2}})^{T}$ are the covariate effects, and $\mathbf{g}(k)\in\mathbb{R}^{v_{2}}$ is the regression function for the covariates $k$. The errors $\varepsilon$ are uncorrelated with zero mean, and their variance depends on the treatment chosen for the $i$th trial: $\mathrm{Var}(\varepsilon(i,k))=\sigma^{2}/\lambda_{i}$. The positive function $\lambda:\{1,\ldots,v_{1}\}\to\mathbb{R}_{++}$ is called the efficiency function and is assumed to be known.

We assume that the set $\mathfrak{T}$ of all potential covariates is finite, which is often the case in practice: the covariates are either naturally discrete (e.g., blocks) or the continuous covariates are discretized. For ease of notation, we number the covariates: $\mathfrak{T}=\{1,\ldots,d\}$. Model (1) can be expressed in a compact form as $y(i,k)=\mathbf{f}^{T}(i,k)\theta+\varepsilon(i,k)$, where $\theta=(\tau^{T},\mu,\beta^{T})^{T}$, and $\mathbf{f}(i,k)=(\mathbf{e}_{i}^{T},1,\mathbf{g}^{T}(k))^{T}$ is the regression function for the given $i$ and $k$.

The approximate design $\xi$ (or the design $\xi$, in short) is a probability measure on the design space $\mathfrak{X}=\{1,\ldots,v_{1}\}\otimes\{1,\ldots,d\}$; i.e., $\xi$ is a nonnegative function from $\mathfrak{X}$ to $\mathbb{R}$ that satisfies $\sum_{i,k}\xi(i,k)=1$. The value $\xi(i,k)$ represents the proportion of all trials that are performed with treatment $i$ and covariates $k$. The exact design $\xi_{e}$, which describes an actual experiment consisting of $n$ trials, specifies the number of trials $\xi_{e}(i,k)$ that are performed with treatment $i$ and covariates $k$ for each $i$ and $k$. The limit on the number of trials means that $\sum_{i,k}\xi_{e}(i,k)=n$. In this paper, we consider approximate designs, unless specified otherwise.

The moment matrix of a design $\xi$ is $\mathbf{M}(\xi)=\sum_{i,k}\xi(i,k)\lambda_{i}\mathbf{f}(i,k)\mathbf{f}^{T}(i,k)$. Let $\mathbf{h}(k):=(1,\mathbf{g}^{T}(k))^{T}$; then, $\mathbf{M}(\xi)$ can be expressed in the block form

[TABLE]

where $\mathbf{M}_{11}(\xi)=\mathrm{diag}(\lambda_{1}w_{1},\ldots,\lambda_{v_{1}}w_{v_{1}})$, \mathbf{M}_{12}(\xi)=\big{(}\sum_{k}\lambda_{1}\xi(1,k)\mathbf{h}(k),\ldots,\sum_{k}\lambda_{v_{1}}\xi(v_{1},k)\mathbf{h}(k)\big{)}^{T} and $\mathbf{M}_{22}(\xi)=\sum_{i,k}\lambda_{i}\xi(i,k)\mathbf{h}(k)\mathbf{h}^{T}(k)$.

The experimental interest in a set of $s_{1}$ treatment contrasts $\mathbf{Q}_{1}^{T}\tau$ and in a set of $s_{2}$ linear functions of the covariate effects $\mathbf{K}^{T}\beta$, where $\mathbf{Q}_{1}\in\mathbb{R}^{v_{1}\times s_{1}}$ and $\mathbf{K}\in\mathbb{R}^{v_{2}\times s_{2}}$, can be expressed as $\mathbf{A}^{T}\theta$, where

[TABLE]

and $\mathbf{Q}_{2}^{T}=(\mathbf{0}_{s_{2}},\mathbf{K}^{T})$. Because $\mathbf{Q}_{1}^{T}\tau$ is a system of contrasts, the matrix $\mathbf{Q}_{1}$ satisfies $\mathbf{Q}^{T}\mathbf{1}_{v_{1}}=\mathbf{0}_{s_{1}}$. The subsystem $\mathbf{A}^{T}\theta$ is estimable under $\xi$ if $\mathcal{C}(\mathbf{A})\subseteq\mathcal{C}(\mathbf{M}(\xi))$. In such a case, we say that $\xi$ is feasible for $\mathbf{A}^{T}\theta$. The constant term $\mu$ is not estimable because of the presence of the treatment effects; therefore, the first row of matrix $\mathbf{Q}_{2}$ is a row of zeros, as formulated above. Let both $\mathbf{Q}_{1}$ and $\mathbf{K}$ be of full column rank. We also suppose that there is interest in all treatments; i.e., no row of $\mathbf{Q}_{1}$ is a row of zeros. The information matrix for $\mathbf{A}^{T}\theta$ of a feasible $\xi$ is $\mathbf{N}_{\mathbf{A}}(\xi)=(\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A})^{-1}$.

A design $\xi^{*}$ is $\Phi$-optimal if it maximizes $\Phi(\mathbf{N}_{\mathbf{A}}(\xi))$ for a given functional $\Phi$ of the information matrix. For example, $\xi^{*}$ is $D$-optimal if it maximizes $\det(\mathbf{N}_{\mathbf{A}}(\xi))$, $A$-optimal if it maximizes $1/\mathrm{tr}(\mathbf{N}_{\mathbf{A}}^{-1}(\xi))$ and $E$-optimal if it maximizes $\gamma_{\min}(\mathbf{N}_{\mathbf{A}}(\xi))$. The mentioned optimality criteria can be extended to an entire class of the so-called Kiefer’s $\Phi_{p}$-optimality criteria, $p\in[-\infty,0]$ (see Pukelsheim (1993), Chapter 6):

[TABLE]

where $\mathbf{N}$ is an $s\times s$ positive definite matrix. The criteria of $D$-, $A$- and $E$-optimality are obtained by setting $p=0$, $-1$ and $-\infty$, respectively.

We require that the optimality criteria possess some basic properties so that they correctly measure the amount of the obtained information. In particular, the criteria must be positively homogeneous, concave, nonnegative, nonconstant and upper semicontinuous; such criteria are called information functions (Pukelsheim (1993), Chapter 5). The optimality criteria are usually eigenvalue based, i.e., they depend only on the eigenvalues of the information matrix. In the present paper, we restrict ourselves to eigenvalue-based information functions. One common class of such functions is the class of the $\Phi_{p}$-criteria.

If the system of interest is rank deficient (i.e., if $\mathbf{A}$ is not of full column rank), then $\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A}$ is never non-singular, and therefore $(\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A})^{-1}$ does not exist. In such cases, especially for eigenvalue-based optimality criteria, the information matrix for $\mathbf{A}^{T}\theta$ is set to be $\mathbf{N}_{\mathbf{A}}(\xi)=(\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A})^{+}$; see Section 8.18 by Pukelsheim (1993). The $\Phi_{p}$-optimality criteria are then defined on the positive eigenvalues of $\mathbf{N}_{\mathbf{A}}(\xi)$. The system $\tau_{i}-\bar{\tau}$ ($i=1,\ldots,v_{1}$) for estimating centered treatment effects, represented by $\mathbf{Q}_{1}=\mathbf{I}_{v_{1}}-\mathbf{1}_{v_{1}}\mathbf{1}_{v_{1}}^{T}/v_{1}$, is a simple example of a rank-deficient system of treatment contrasts.

1.3 Marginal models

For the first marginal model of (1), we consider the heteroscedastic model for treatment effects

[TABLE]

where $\mathrm{Var}(\varepsilon(i))=\sigma^{2}/\lambda_{i}$. We denote the marginal treatment design $w$ as the approximate design in (2). That is, $w$ specifies $v_{1}$ nonnegative treatment weights that sum to one; we denote these weights as $w_{1},\ldots,w_{v_{1}}$. The moment matrix of $w$ is $\mathbf{M}_{1}(w)=\mathrm{diag}(\lambda_{1}w_{1},\ldots,\lambda_{v_{1}}w_{v_{1}})$, and $w$ is feasible for $\mathbf{Q}_{1}^{T}\tau$ if $w_{i}>0$ for all $i=1,\ldots,v_{1}$, denoted as $w>0$. The information matrix for $\mathbf{Q}_{1}^{T}\tau$ of $w>0$ is $\mathbf{N}_{\mathbf{Q}_{1}}(w)=(\mathbf{Q}_{1}^{T}\mathbf{M}_{1}^{-1}(w)\mathbf{Q}_{1})^{-1}$.

For the second marginal model, we consider the homoscedastic model for covariates with constant term:

[TABLE]

Then, the marginal covariate design $\alpha$ is an approximate design in (3). The marginal covariate design is analogously characterized by the $d$ covariate weights $\alpha_{1},\ldots,\alpha_{d}$, and its moment matrix is

[TABLE]

The subsystem $\mathbf{K}^{T}\beta$ is estimable if $\mathcal{C}(\mathbf{Q}_{2})\subseteq\mathcal{C}(\mathbf{M}_{2}(\alpha))$, which is equivalent to $\mathcal{C}(\mathbf{K})\subseteq\mathbf{S}(\alpha)$, where $\mathbf{S}(\alpha)=\sum_{k}\alpha_{k}\mathbf{g}(k)\mathbf{g}^{T}(k)-(\sum_{k}\alpha_{k}\mathbf{g}(k))(\sum_{k}\alpha_{k}\mathbf{g}^{T}(k))$ is the Schur complement of $\mathbf{M}_{2}(\alpha)$. The Schur complement of a nonnegative definite matrix in the block form

[TABLE]

is $\mathbf{B}_{\mathbf{S}}=\mathbf{B}_{22}-\mathbf{B}_{12}^{T}\mathbf{B}_{11}^{-1}\mathbf{B}_{12}$. More precisely, $\mathbf{B}_{\mathbf{S}}$ is the Schur complement of $\mathbf{B}_{11}$ in $\mathbf{B}$.

The information matrix for $\mathbf{K}^{T}\beta$ of a feasible $\alpha$ is $\mathbf{N}_{\mathbf{K}}(\alpha)=(\mathbf{Q}_{2}^{T}\mathbf{M}_{2}^{-}(\alpha)\mathbf{Q}_{2})^{-1}$. The information matrix can be expressed using the Schur complement, because the matrix

[TABLE]

is a generalized matrix of an arbitrary nonnegative definite matrix $\mathbf{B}$ given by (4); see Theorem 9.6.1 by Harville (1997). Formula (5) implies for $\mathbf{M}_{2}(\alpha)$ that $\mathbf{N}_{\mathbf{K}}(\alpha)=(\mathbf{K}^{T}\mathbf{S}^{-}(\alpha)\mathbf{K})^{-1}$.

For a design $\xi$ in model (1), the marginal treatment design $w$ of $\xi$ is given by $w_{i}=\sum_{k}\xi(i,k)$ ($i=1,\ldots,v_{1}$) and the marginal covariate design $\alpha$ of $\xi$ is given by $\alpha_{k}=\sum_{i}\xi(i,k)$ ($k=1,\ldots,d$).

2 Optimal product designs

The product design $\xi=w\otimes\alpha$ of the marginal designs $w$ and $\alpha$ satisfies $\xi(i,k)=w_{i}\alpha_{k}$ ($i=1,\ldots,v_{1}$, $k=1,\ldots,d$).

Proposition 1**.**

If $w>0$ and if $\alpha$ is feasible for $\mathbf{K}^{T}\beta$, then $\xi=w\otimes\alpha$ is feasible for $\mathbf{A}^{T}\theta$ and its information matrix is

[TABLE]

Proof.

For $\xi=w\otimes\alpha$, a generalized inverse of $\mathbf{M}(\xi)$ can be expressed as in (5) using the Schur complement

[TABLE]

Observe that $\mathbf{M}_{11}(\xi)=\mathbf{M}_{1}(w)$ and that

[TABLE]

Because $\mathbf{Q}_{1}^{T}\mathbf{1}=\mathbf{0}$, we obtain $\mathbf{Q}_{1}^{T}\mathbf{M}_{11}^{-1}(\xi)\mathbf{M}_{12}(\xi)=\mathbf{0}$. Let us calculate $\mathbf{M}(\xi)\mathbf{M}^{-}(\xi)\mathbf{A}$, where $\mathbf{M}^{-}(\xi)$ is given by (5). By employing the abovementioned observations, it is straightforward to show that $\mathbf{M}(\xi)\mathbf{M}^{-}(\xi)\mathbf{A}=\mathrm{diag}(\mathbf{Q}_{1},\mathbf{S}_{1}(\xi)\mathbf{S}_{1}^{-}(\xi)\mathbf{Q}_{2})$. Because $\alpha$ is feasible for $\mathbf{K}^{T}\beta$, we have $\mathcal{C}(\mathbf{K})\subseteq\mathcal{C}(\mathbf{S}(\alpha))$, which implies $\mathcal{C}(\mathbf{Q}_{2})\subseteq\mathcal{C}(\mathbf{S}_{1}(\xi))$. Then, $\mathbf{S}_{1}(\xi)\mathbf{S}_{1}^{-}(\xi)\mathbf{Q}_{2}=\mathbf{Q}_{2}$, which yields $\mathbf{M}(\xi)\mathbf{M}^{-}(\xi)\mathbf{A}=\mathbf{A}$; i.e., $\xi$ is feasible for $\mathbf{A}^{T}\theta$.

The information matrix of $\xi$ can be expressed using the same generalized inverse $\mathbf{M}^{-}(\xi)$ and the fact that $\mathbf{M}_{11}^{-1}(\xi)\mathbf{M}_{12}(\xi)\mathbf{Q}_{1}=\mathbf{0}$:

[TABLE]

which is equivalent to (6). ∎

The following preliminary lemma shows that no design is better with respect to any eigenvalue-based criterion than the product of its marginals.

Lemma 1**.**

Let $\Phi$ be an eigenvalue-based information function, let $\xi$ be a design in model (1) and let $w$ and $\alpha$ be its marginal designs. Then, $\Phi(\mathbf{N}_{\mathbf{A}}(\xi))\leq\Phi(\mathbf{N}_{\mathbf{A}}(w\otimes\alpha))$.

Proof.

Similarly to the proof of Lemma 3.1 by Schwabe (1996), let us consider the modification of model (1), where $\tau_{i}$ is changed to $-\tau_{i}$. This corresponds to the change in regressors from $\mathbf{f}(i,k)=(\mathbf{e}_{i}^{T},1,\mathbf{g}^{T}(k))$ to $\tilde{\mathbf{f}}(i,k)=(-\mathbf{e}_{i}^{T},1,\mathbf{g}^{T}(k))$. The moment matrix $\tilde{\mathbf{M}}(\xi)$ in the modified model is

[TABLE]

Then, the information matrix changes from $\mathbf{N}_{\mathbf{A}}(\xi)$ to

[TABLE]

Because the matrix $\mathrm{diag}(-\mathbf{I},\mathbf{I})$ is orthogonal, the information matrices $\mathbf{N}_{\mathbf{A}}(\xi)$ and $\tilde{\mathbf{N}}_{\mathbf{A}}(\xi)$ have the same eigenvalues; hence $\Phi(\mathbf{N}_{\mathbf{A}}(\xi))=\Phi(\tilde{\mathbf{N}}_{\mathbf{A}}(\xi))$. From concavity of $\Phi$ it follows that

[TABLE]

Moreover, $(\mathbf{N}_{\mathbf{A}}(\xi)+\tilde{\mathbf{N}}_{\mathbf{A}}(\xi))/2=\mathrm{diag}(\mathbf{N}_{11}(\xi),\mathbf{N}_{22}(\xi))=:\mathbf{N}_{*}(\xi)$. By expressing the generalized inverse $\mathbf{M}^{-}(\xi)$ through the Schur complement $\mathbf{S}_{1}(\xi)$ (see (5)), we obtain that $\mathbf{N}_{*}(\xi)=(\mathbf{A}^{T}\mathrm{diag}(\mathbf{M}_{11}^{-1}+\mathbf{M}_{11}^{-1}\mathbf{M}_{12}\mathbf{S}_{1}^{-}\mathbf{M}_{12}^{T}\mathbf{M}_{11}^{-1},\mathbf{S}_{1}^{-})\mathbf{A})^{-1}$. Moreover, $\mathrm{diag}(\mathbf{M}_{11}^{-1}+\mathbf{M}_{11}^{-1}\mathbf{M}_{12}\mathbf{S}_{1}^{-}\mathbf{M}_{12}^{T}\mathbf{M}_{11}^{-1},\mathbf{S}_{1}^{-})\succeq\mathrm{diag}(\mathbf{M}_{11}^{-1},\mathbf{S}_{1}^{-})$, which implies that $\mathbf{N}_{*}(\xi)\preceq(\mathbf{A}^{T}\mathrm{diag}(\mathbf{M}_{11}^{-1},\mathbf{S}_{1}^{-})\mathbf{A})^{-1}=\mathbf{N}_{\mathbf{A}}(w\otimes\alpha)$. Hence, $\Phi(\mathbf{N}_{\mathbf{A}}(w\otimes\alpha))\geq\Phi(\mathbf{N}_{*}(\xi))\geq\Phi(\mathbf{N}_{\mathbf{A}}(\xi))$. ∎

For any eigenvalue-based $\Phi$, Lemma 1 implies that if there exists a $\Phi$-optimal design, then there exists a product design that is $\Phi$-optimal among all designs. In particular, if $\xi^{*}$ is $\Phi$-optimal with marginal treatment design $w^{*}$ and marginal covariate design $\alpha^{*}$, then $w^{*}\otimes\alpha^{*}$ is $\Phi$-optimal.

Theorem 1**.**

Let $\Phi$ be an eigenvalue-based information function. Then, if there exists a $\Phi$-optimal design for $\mathbf{A}^{T}\theta$ in (1), then there exists a product design in (1) that is $\Phi$-optimal for $\mathbf{A}^{T}\theta$.

Theorem 1 shows that to obtain a $\Phi$-optimal design it suffices to find a $\Phi$-optimal product design. For $D$-optimality, the problem of finding a $\Phi$-optimal product design simplifies because of its multiplicative form:

[TABLE]

This means that the optimal $w$ and $\alpha$ can be computed separately, which was extensively used by Wang and Ai (2016).

For the other $\Phi_{p}$-optimality criteria, we reduce the complexity of the problem as follows: For $p\in(-\infty,0)$, we have

[TABLE]

The additive form allows one to first calculate the covariate design $\alpha^{*}$, which is $\Phi_{p}$-optimal in the marginal model (3), by maximizing $\Phi_{p}(\mathbf{N}_{\mathbf{K}}(\alpha))$ and then calculate the corresponding optimal $w^{*}$ by maximizing $\Phi_{p}(\mathbf{N}_{\mathbf{A}}(w\otimes\alpha^{*}))$ over all $w$. Then, from (7) it follows that the product design $w^{*}\otimes\alpha^{*}$ is $\Phi_{p}$-optimal. Analogously, to calculate the $E$-optimal ($p=-\infty$) product design, one can first find the $E$-optimal $\alpha^{*}$ in (3), which maximizes $\gamma_{\min}(\mathbf{N}_{\mathbf{K}}(\alpha))$, and then the optimal $w^{*}$ that maximizes $\gamma_{\min}(\mathbf{N}_{\mathbf{A}}(w\otimes\alpha^{*}))$. Therefore, for any $\Phi_{p}$-optimality criterion, the optimization problem of size $v_{1}\cdot d$ can be split into two maximization problems of sizes $d$ and $v_{1}$, respectively, even for $p$ other than [math].

Theorem 2**.**

(i) Let $p\in(-\infty,0)$, let $\alpha^{*}$ be $\Phi_{p}$-optimal for $\mathbf{K}^{T}\beta$ in (3) and denote $\varphi^{*}=\mathrm{tr}(\mathbf{N}^{p}_{\mathbf{K}}(\alpha^{*}))$. Let $w^{*}$ maximize $\Phi_{p}(\mathbf{N}_{\mathbf{A}}(w\otimes\alpha^{*}))=(s^{-1}(\mathrm{tr}(\mathbf{N}^{p}_{\mathbf{Q}_{1}}(w))+(\sum_{i}\lambda_{i}w_{i})^{p}\varphi^{*}))^{1/p}$. Then, $\xi^{*}=w^{*}\otimes\alpha^{*}$ is $\Phi$-optimal for $\mathbf{A}^{T}\theta$ in (1).

(ii) Let $\alpha^{*}$ be $E$-optimal for $\mathbf{K}^{T}\beta$ in (3) and let $w^{*}$ maximize $\gamma_{\min}(\mathbf{N}_{\mathbf{A}}(w\otimes\alpha^{*}))$. Then, $\xi^{*}=w^{*}\otimes\alpha^{*}$ is $E$-optimal for $\mathbf{A}^{T}\theta$ in (1).

Theorem 2 shows that $\alpha^{*}$ in the $\Phi_{p}$-optimal product design $w^{*}\otimes\alpha^{*}$ is the $\Phi_{p}$-optimal marginal covariate design, but $w^{*}$ is generally not $\Phi_{p}$-optimal in (2). Due to heteroscedasticity, a slight correction (represented by $\varphi^{*}$) needs to be present to calculate the “optimal” marginal treatment design $w^{*}$.

3 Optimal non-product designs

In the previous section, we showed that to calculate an optimal design, it suffices to restrict oneself to the product designs. However, in model (1), besides an optimal product design, there usually exists a rich class of optimal designs that are not of the product form. Fortunately, once a $\Phi$-optimal product design $\xi^{*}$ is found, other $\Phi$-optimal designs can generally be constructed from $\xi^{*}$. The following theorem provides such a construction.

Theorem 3**.**

Let $\Phi$ be an eigenvalue-based information function and let $\xi^{*}=w^{*}\otimes\alpha^{*}$ be the $\Phi$-optimal product design for $\mathbf{A}^{T}\theta$. Let $\xi$ satisfy $\mathbf{M}(\xi)\mathbf{G}\mathbf{A}=\mathbf{A}$, where

[TABLE]

and where $\mathbf{M}_{2}^{-}(\alpha^{*})$ is any generalized inverse of $\mathbf{M}_{2}(\alpha^{*})$. Then, $\xi$ is $\Phi$-optimal for $\mathbf{A}^{T}\theta$.

In particular, let $\xi$ satisfy

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $w$ is the marginal treatment design of $\xi$. Then, $\xi$ is $\Phi$-optimal for $\mathbf{A}^{T}\theta$.

Proof.

Lemma 1 by Rosa and Harman (2016) shows that if $\tilde{\mathbf{M}}$ is a non-negative definite matrix and if a design $\xi$ satisfies $\mathbf{M}(\xi)\tilde{\mathbf{M}}^{-}\mathbf{A}=\mathbf{A}$, then $\xi$ is feasible for $\mathbf{A}^{T}\theta$ and $\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A}=\mathbf{A}^{T}\tilde{\mathbf{M}}^{-}\mathbf{A}$. We choose $\tilde{\mathbf{M}}=\mathrm{diag}(\mathbf{M}_{1}(w^{*}),(\sum_{i}\lambda_{i}w_{i}^{*})\mathbf{M}_{2}(\alpha^{*}))$ and

[TABLE]

Then $\mathbf{G}$ is a generalized inverse of $\tilde{\mathbf{M}}$, and any $\xi$ satisfying $\mathbf{M}(\xi)\mathbf{G}\mathbf{A}=\mathbf{A}$ satisfies $\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A}=\mathbf{A}^{T}\mathbf{G}\mathbf{A}=\mathbf{N}_{\mathbf{A}}^{-1}(\xi^{*})$; i.e., such $\xi$ is $\Phi$-optimal.

For the second part, for the generalized inverse of $\mathbf{M}_{2}(\alpha^{*})$ we choose

[TABLE]

The condition $\mathbf{M}(\xi)\mathbf{G}\mathbf{A}=\mathbf{A}$ consists of the following equalities:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Because no row of $\mathbf{Q}_{1}$ is a row of zeros, condition (13) implies that $w_{i}=w_{i}^{*}$. Equality (14) can be expressed as

[TABLE]

which can be simplified to (10). Equality (15) can be simplified to

[TABLE]

Finally, condition (16) can be expressed by the following equalities:

[TABLE]

[TABLE]

However, (17) can be simplified to

[TABLE]

which follows from (10). ∎

Because the conditions in Theorem 3 are linear, optimal designs satisfying these conditions can be calculated via linear programming (LP) once $w^{*}$ and $\alpha^{*}$ are computed:

[TABLE]

where $\mathbf{x}$ represents the $v_{1}d$ design values, $\mathbf{C}\mathbf{x}=\mathbf{b}$ consists of the conditions given by Theorem 3 and of the design constraint $\sum_{j}\mathbf{x}_{j}=1$. The vector $\mathbf{c}$ can be chosen arbitrarily, as the objective is to find any $\mathbf{x}\geq\mathbf{0}$ that satisfies $\mathbf{C}\mathbf{x}=\mathbf{b}$.

The possibility of employing linear programming has two crucial advantages. The more obvious one is that it is relatively simple to solve the LP problems, and most mathematical software packages (e.g., MATLAB and R) contain reliable and fast LP solvers. The other advantage is that the vertices of the set of the feasible solutions of (18) have high numbers of zeroes among all feasible solutions; for technical details, see, e.g., Theorem 2.4 by Bertsimas and Tsitsiklis (1997). In fact, any vertex solution of (18) is guaranteed to have at least $v_{1}d-r$ zeros, where $r=\mathrm{rank}(\mathbf{C})$. This is beneficial, because such designs with small supports can then be obtained by solving the LP problem via the simplex method, as this method provides vertex solutions. Note that the MATLAB implementation of the interior point method also seems to provide vertex solutions. Theorem 3 therefore allows us to formulate a linear programming “sparsification” method of the product designs based on solving (18). In Section 4, we demonstrate the applicability of this method.

Remarks

• •

For simpler settings (e.g., in the homoscedastic case or if the interest lies in $\mathbf{Q}_{1}^{T}\tau$ only), the obtained results simplify correspondingly. For instance, if the interest lies only in a set of treatment contrasts (i.e., $\mathbf{A}=(\mathbf{Q}_{1}^{T},\mathbf{0})^{T}$), then the product design $w^{*}\otimes\alpha$ is $\Phi$-optimal for any marginal covariate design $\alpha$, where $w^{*}$ is a $\Phi$-optimal marginal treatment design for $\mathbf{Q}_{1}^{T}\tau$. Moreover, for such $\mathbf{A}^{T}\theta$, any design $\xi$ that satisfies (9) and whose marginal treatment design is $w^{*}$ is $\Phi$-optimal. The designs satisfying (9) were denoted as resistant to nuisance effects in a slightly different context by Rosa and Harman (2016). In the present settings, designs satisfying (8) and (9) may be called covariate resistant: if the interest lies only in $\mathbf{Q}_{1}^{T}\tau$, then no relevant information is lost under such designs due to the presence of covariates.

• •

In the case of a rank-deficient system, both the results of Section 2 and of Section 3 hold, as the proofs in Section 2 can be easily adapted to $\mathbf{N}_{\mathbf{A}}(\xi)=(\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A})^{+}$ and Theorem 3 ensures that $\mathbf{A}^{T}\mathbf{M}^{-}(\xi)\mathbf{A}=\mathbf{A}^{T}\tilde{\mathbf{M}}^{-}\mathbf{A}$, which also implies the equality of the rank-deficient information matrices.

• •

If a $\Phi$-optimal design $\xi^{*}$ is known, other $\Phi$-optimal designs can trivially be found by solving the linear equality $\sum_{i,k}\xi(i,k)\mathbf{f}(i,k)\mathbf{f}^{T}(i,k)=\mathbf{M}(\xi^{*})$ that guarantees that the design $\xi$ has the same moment matrix as $\xi^{*}$. The conditions in Theorem 3 are more general, as it can be shown that any design satisfying $\mathbf{M}(\xi)=\mathbf{M}(\xi^{*})$ also satisfies $\mathbf{M}(\xi)\mathbf{G}\mathbf{A}=\mathbf{A}$ with $\mathbf{G}$ given by Theorem 3. Conditions (8)-(9) can generally be satisfied also by designs $\xi$ that do not have the same moment matrix as $\xi^{*}$; only the information matrices of $\xi$ and $\xi^{*}$ are guaranteed to coincide.

4 Examples

Example 1**.**

Consider a model with effects of $v_{2}$ continuous covariates

[TABLE]

where $\mathbf{z}\in[-1,1]^{v_{2}}$, as in Example 2 by Wang and Ai (2016). By a discretization of the continuous covariates, say $z_{i}\in\{\pm j/10\ |\ j=0,\ldots,10\}$, the model becomes a special case of (1). Suppose that we wish to obtain $A$-optimal designs for estimating the comparisons with the control $\tau_{i}-\tau_{1}$, $i=2,\ldots,v_{1}$ and all the covariate effects; i.e., $\mathbf{Q}=(-\mathbf{1}_{v_{1}-1},\mathbf{I}_{v_{1}-1})^{T}$ and $\mathbf{K}=\mathbf{I}_{v_{2}}$. Then, the marginal covariate design $\alpha^{*}$ that is uniform on $\{-1,1\}^{v_{2}}$ is $A$-optimal for $\mathbf{K}^{T}\beta$. The corresponding information matrix is $\mathbf{N}_{\mathbf{K}}(\alpha^{*})=\mathbf{I}_{v_{2}}$, and hence $\varphi^{*}=\mathrm{tr}(\mathbf{I}_{v_{2}}^{-1})=v_{2}$. Then, the optimal marginal treatment design $w^{*}$ can be obtained by minimizing $\mathrm{tr}(\mathbf{Q}^{T}\mathbf{M}_{1}^{-1}(w)\mathbf{Q})+v_{2}(\sum_{i}\lambda_{i}w_{i})^{-1}$.

Let $v_{1}=v_{2}=3$ and let $\lambda_{1}=9$, $\lambda_{2}=\lambda_{3}=1$; i.e., the observations under the control treatment have smaller variances. Then, the vector of the optimal treatment weights given by $w^{*}$ is $(0.236,0.382,0.382)^{T}$, and hence the product design $\xi^{*}=w^{*}\otimes\alpha^{*}$ (see Table 1) is $A$-optimal for $\mathbf{A}^{T}\theta$. As observed in Theorem 2, the optimal treatment design $w^{*}$ depends on $\varphi^{*}$, which is equal to the number of covariates $v_{2}$ in the current example. Figure 1 depicts the dependence of $w_{1}^{*}$ on $v_{2}$ for the abovementioned settings.

An $A$-optimal design, say $\xi_{s}^{*}$, that is supported on a smaller number of design points was computed by solving the linear program (18). Although $\xi_{s}^{*}$ was allowed to attain non-zero values even outside of the support of $\xi^{*}$, it turns out that the support of $\xi_{s}^{*}$ is a subset of the support of $\xi^{*}$. Therefore, $\xi_{s}^{*}$ can also be expressed by considering only $\mathbf{z}\in\{-1,1\}^{3}$, see Table 2. Whereas the support of the optimal product design is of size 24, the support of the optimal non-product design is only of size 10. One can also observe that the marginal treatment design of $\xi_{s}^{*}$ is $w^{*}$, but $\alpha^{*}$ is not its marginal covariate design.

Let us demonstrate the usefulness of the sparsely supported designs for constructing the exact designs by the efficient rounding procedure by Pukelsheim and Rieder (1992). For $n=48$ trials the rounded product design $\xi_{ep}$ satisfies $\xi_{ep}(i,\mathbf{z})=2$ for each $i$ and for each $\mathbf{z}\in\{-1,1\}^{3}$; the rounded non-product design $\xi_{es}$ is given in Table 3. To compare the quality of the obtained designs, we calculate their efficiencies: $\mathrm{eff}(\xi_{e})=\Phi(\mathbf{N}_{\mathbf{A}}(\xi_{e}/n))/\Phi(\mathbf{N}_{\mathbf{A}}(\xi^{*}))$, where $\xi_{e}$ is an exact design for $n$ trials and $\xi^{*}$ is the $\Phi$-optimal approximate design. We have $\mathrm{eff}(\xi_{ep})=0.9641$ for the rounded product design and $\mathrm{eff}(\xi_{es})=0.9991$ in the non-product case. The very high efficiency of the exact design $\xi_{es}$ relative to the optimal approximate design means that $\xi_{es}$ is either an optimal exact design, or very nearly optimal.

The efficient rounding method requires that the number of trials $n$ be equal or greater than the support size of the approximate design. Therefore, the efficient rounding of $\xi^{*}$ cannot be performed for $n\leq 23$, as the support size of $\xi^{*}$ is $24$. However, the rounding of the non-product $\xi_{s}^{*}$ can be done even for $10\leq n\leq 23$. This is demonstrated in Figure 2, which shows efficiencies of the respective rounded designs for varying $n$. Note that the rounded non-product design is not always more efficient than the rounded product design.

Example 2**.**

Suppose that the responses are affected, besides the treatment effects, by two qualitative covariates:

[TABLE]

where $j\in\{1,\ldots,r\}$ is the level of the first covariate, $k\in\{1,\ldots,c\}$ is the level of the second covariate, and $\eta_{j}$ and $\rho_{k}$ are the respective covariate effects. In the optimal design literature, such a model is usually known as two-way elimination of heterogeneity (e.g., see Hedayat et al. (1988)) or a row-column model (e.g., see Jacroux (1982)), where the trials are split into $r$ rows and $c$ columns. If the effects of one of the covariates is disregarded, we obtain the well-known model with block effects. Suppose that all the treatments are of the same interest, which can be expressed by estimating the treatment effects corrected for the mean $\tau_{i}-\bar{\tau}$; i.e., $\mathbf{Q}_{1}=\mathbf{I}_{v_{1}}-\mathbf{1}_{v_{1}}\mathbf{1}_{v_{1}}^{T}/{v_{1}}$. Similarly, let $\mathbf{K}=\mathrm{diag}(\mathbf{I}_{r}-\mathbf{1}_{r}\mathbf{1}_{r}^{T}/r,\mathbf{I}_{c}-\mathbf{1}_{c}\mathbf{1}_{c}^{T}/c)$. Let $v_{1}=3$, $r=3$, $c=5$ and $\lambda=(4,1,1)^{T}$. We shall provide $E$-optimal designs for estimating $\mathbf{A}^{T}\theta$.

The $E$-optimal marginal covariate design for $\mathbf{K}^{T}\beta$ is uniform: $\alpha^{*}(j,k)=1/15$ for all $j,k$. The optimal criterion value $\Phi_{E}(\mathbf{N}_{\mathbf{Q}_{2}}(\alpha^{*}))$ is $0.2$, and then the corresponding optimal $w^{*}$ calculated using Theorem 2 is given by the vector of treatment weights $(0.273,0.364,0.364)^{T}$. Then, the $E$-optimal product design $w^{*}\otimes\alpha^{*}$ is supported on all the $v_{1}\times rc=45$ design points, with values $\xi(1,j,k)=0.0182$ and $\xi(2,j,k)=\xi(3,j,k)=0.0242$ for each $j$ and $k$. By solving the linear program (18), an $E$-optimal design $\xi_{s}^{*}$ with a smaller number of support points can be obtained. Such a design is given in Table 4. The non-product design $\xi_{s}^{*}$ has a smaller support of only size 28, as expressed by the large number of zeroes in Table 4. Table 5 gives the rounded non-product design $\xi_{s}^{*}$ for $n=40$ trials; for this number of trials, the efficient rounding of the optimal product design cannot be performed because of its large support size of $48>n$. The efficiency of the exact design obtained by the rounding of $\xi_{s}^{*}$ is $0.8493$.

Example 3**.**

Let us demonstrate the obtained theoretical results in settings, where there is no interest in the covariates. Consider a model with exponential trend effect, which is considered to be a nuisance:

[TABLE]

where $g(k)=e^{k}/\sum_{j}e^{j}$, $k=1,\ldots,d$. In each time $k=1,\ldots,d$, the same number of trials must be performed, which results in the design constraint $\sum_{i}\xi(i,k)=1/d$ for each $k$. Suppose that the interest lies in the treatment-control comparisons only; i.e., $\mathbf{Q}^{T}=(-\mathbf{1}_{v_{1}-1},\mathbf{I}_{v_{1}-1})$ and $\mathbf{A}^{T}=(\mathbf{Q}^{T},\mathbf{0})$. Let $v_{1}=4$, $d=6$, $\lambda=(1,1,2,3)^{T}$ and consider the $A$-optimality criterion.

The $A$-optimal marginal treatment design $w^{*}$ can be found by maximizing $\Phi_{A}(\mathbf{N}_{\mathbf{Q}}(w))$, and its values are given by the vector $(0.431,0.249,0.176,0.144)^{T}$. The $A$-optimal product design is then $\xi^{*}=w^{*}\otimes\alpha$, where the marginal covariate design $\alpha_{k}=1/d$ for $k=1,\ldots,d$ is implied by the design constraints $\sum_{i}\xi(i,k)=1/d$. An optimal non-product design can be obtained as in Theorem 3 by solving the linear program with constraints $\xi(i,k)\geq 0$ for all $i,k$, $\sum_{i}\xi(i,k)=1/d$ for all $k$, and $\mathbf{M}(\xi)\mathbf{G}\mathbf{A}=\mathbf{A}$, where $\mathbf{G}=\mathrm{diag}(\mathbf{M}_{1}^{-1}(w^{*}),\mathbf{0})$. The resulting design $\xi_{s}^{*}$, which is supported on 12 design points, is given in Table 6; the optimal product design is supported on all the 24 points.

Because of the requirement that exactly one trial should be performed in each time moment, the usual rounding methods cannot be used to obtain efficient exact designs. A natural rounding method for the sparsely supported non-product designs is to select in each time the treatment with the highest design value; such a rounded design based on $\xi_{s}^{*}$ has efficiency 0.8871 and is given in Table 7. Note that it is unclear how the optimal product designs can be rounded in the current example – these designs do not provide any information on the suggested time sequence of the treatments, only the optimal treatment proportions are obtained. For instance, the aforementioned method would result for $\xi^{*}$ in a singular design that selects the first treatment in each time moment.

5 Discussion

Although $D$-optimality is very beneficial analytically and computationally in the considered model, we showed that other eigenvalue-based optimality criteria like $A$- and $E$-optimality are not much more difficult to work with (Theorems 1 and 2). Because the latter place actual emphasis on the selected system of interest, they seem to be more appropriate in the current settings.

The use of optimal product designs is also computationally very beneficial, because they allow one to reduce the $v_{1}d$-dimensional optimal design problem in the multi-factor model (1) to two simpler problems of sizes $v_{1}$ and $d$ for the two marginal models. However, the product designs suffer from large support sizes, as noted earlier. The support sizes of the product designs can in some settings be reduced by combinatorial approaches; e.g., by replacing full factorials by orthogonal arrays as in Graßhoff et al. (2004). Such reductions of support sizes are however very dependent on the particular well-studied models with a very regular structure that can be made use of in the combinatorial arguments.

To reduce the support sizes in more general settings, we provide an entire class of optimal designs (see Theorem 3) characterized by linear constraints. Such a characterization allows one to use linear programming for constructing optimal non-product designs. That the obtained constraints for optimal designs are linear is particularly useful, because the linear programming routines tend to provide vertex solutions, which have high numbers of zeros. As a result, the solvers “automatically” provide designs with small supports. As we demonstrated on examples in Section 4, it seems that the method tends to allow for significant reductions of the support sizes. The use of linear programming also means that the algorithm is reliable and fast. Unlike the combinatorial approaches, the proposed method is not tailored for specific settings – it can be used for any model of the form (1).

We would like to emphasize that the proposed support reduction method is not particularly tied to the current model (1). As noted earlier, the method was applied for a similar model by Rosa and Harman (2016). Moreover, the proposed approach can be used for any linear regression model where an optimal design $\xi^{*}$ with a large support is found, either analytically or by design algorithms. Then one can find designs that satisfy either $\mathbf{M}(\xi)=\mathbf{M}(\xi^{*})$ or $\mathbf{M}(\xi)\mathbf{G}\mathbf{A}=\mathbf{A}$ with a suitably chosen $\mathbf{G}$ as in Theorem 3 via linear programming to obtain vertex solutions. As shown in the present paper, the vertex solutions correspond to designs that tend to have smaller supports and are therefore generally more useful for practical purposes.