ForLion: a new algorithm for D-optimal designs under general parametric statistical models with mixed factors

In this paper, we address the problem of designing an experimental plan with both discrete and continuous factors under fairly general parametric statistical models. We propose a new algorithm, named ForLion, to search for locally optimal approximate designs under the D-criterion. The algorithm performs an exhaustive search in a design space with mixed factors while keeping high efficiency and reducing the number of distinct experimental settings. Its optimality is guaranteed by the general equivalence theorem. We present the relevant theoretical results for multinomial logit models (MLM) and generalized linear models (GLM), and demonstrate the superiority of our algorithm over state-of-the-art design algorithms using real-life experiments under MLM and GLM. Our simulation studies show that the ForLion algorithm could reduce the number of experimental settings by 25% or improve the relative efficiency of the designs by 17.5% on average. Our algorithm can help the experimenters reduce the time cost, the usage of experimental devices, and thus the total cost of their experiments while preserving high efficiencies of the designs.

No datasets were generated or analysed during the current study.

The authors gratefully acknowledge the authors of Ai et al. (2023), Lukemire et al. (2019) and Lukemire et al. (2022) for kindly sharing their source codes, which we used to implement and compare their methods with ours. The authors gratefully acknowledge the support from the U.S. NSF grants DMS-1924859 and DMS-2311186.

Authors and Affiliations

1. Department of Mathematics, Statistics, and Computer Science, University of Illinois at Chicago, Chicago, IL, 60607, USA

   Yifei Huang & Jie Yang

2. Department of Mathematics, University of Alabama at Birmingham, Birmingham, AL, 35294, USA

   Keren Li

3. Department of Statistics, University of Georgia, Athens, GA, 30602, USA

   Abhyuday Mandal

Contributions

Conceptualization and methodology, all authors; software, Y.H., K.L., and J.Y.; validation, Y.H., A.M., and J.Y.; formal analysis, Y.H. and K.L.; investigation, Y.H. and J.Y.; resources, all authors; data curation, Y.H.; writing-original draft preparation, all authors; writing-review and editing, all authors.; supervision, J.Y. and A.M.; project administration, J.Y. and A.M.; funding acquisition, J.Y. and A.M.. All authors have read and agreed to the published version of the manuscript.

Corresponding author

Correspondence to Jie Yang.

Conflict of interest

The authors declare no conflict of interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Computing \(u_{st}^{{\textbf{x}}}\) in Fisher information \({\textbf{F}}_{{\textbf{x}}}\)

In this section, we provide more technical details for Sect. 3.1 and Theorem 2.

For MLM (1), Corollary 3.1 in Bu et al. (2020) provided an alternative form \({{\textbf{F}}}_{{{\textbf{x}}}_i} = {\textbf{X}}_i^T {{\textbf{U}}}_i {{\textbf{X}}}_i\), which we use for computing the Fisher information \({{\textbf{F}}}_{{\textbf{x}}}\) at an arbitrary \({{\textbf{x}}} \in {{\mathcal {X}}}\). More specifically, first of all, the corresponding model matrix at \({{\textbf{x}}}\) is

where \({{\textbf{h}}}^T_j(\cdot ) = (h_{j1}(\cdot ), \ldots , h_{jp_j}(\cdot ))\) and \({{\textbf{h}}}^T_c(\cdot )\) \( = \) \((h_1(\cdot ), \ldots , h_{p_c}(\cdot ))\) are known predictor functions. We let \(\varvec{\beta }_j\) and \(\varvec{\zeta }\) denote the model parameters associated with \({{\textbf{h}}}^T_j({{\textbf{x}}})\) and \({{\textbf{h}}}^T_c({{\textbf{x}}})\), respectively, then the model parameter vector \(\varvec{\theta }=(\varvec{\beta }_{1},\varvec{\beta }_{2},\cdots ,\varvec{\beta }_{J-1},\varvec{\zeta })^T \in {\mathbb {R}}^p\), and the linear predictor \(\varvec{\eta }_{{\textbf{x}}} = {{\textbf{X}}}_{{\textbf{x}}} \varvec{\theta }= (\eta _1^{{\textbf{x}}}, \ldots , \eta _{J-1}^{{\textbf{x}}}, 0)^T \in {\mathbb {R}}^J\), where \(\eta _j^{{\textbf{x}}} = {{\textbf{h}}}_j^T({{\textbf{x}}}) \varvec{\beta }_j + {{\textbf{h}}}_c^T({{\textbf{x}}}) \varvec{\zeta }\), \(j=1, \ldots , J-1\).

According to Lemmas S.10, S.12 and S.13 in the Supplementary Material of Bu et al. (2020), the categorical probabilities \(\varvec{\pi }_{{\textbf{x}}} = (\pi _1^{{\textbf{x}}}, \ldots , \pi _J^{{\textbf{x}}})^T \in {\mathbb {R}}^J\) at \({{\textbf{x}}}\) for baseline-category, adjacent-categories and continuation-ratio logit models can be expressed as follows:

for \(j=1, \ldots , J-1\), where \(D_j = \exp \{\eta _{J-1}^{{\textbf{x}}} + \cdots + \eta _1^{{\textbf{x}}}\} + \exp \{\eta _{J-1}^{{\textbf{x}}} + \cdots + \eta _2^{{\textbf{x}}}\} + \cdots +\exp \{\eta _{J-1}^{{\textbf{x}}}\} + 1\), and

where \(D_J = \exp \{\eta _{J-1}^{{\textbf{x}}} + \cdots + \eta _1^{{\textbf{x}}}\} + \exp \{\eta _{J-1}^{{\textbf{x}}} + \cdots + \eta _2^{{\textbf{x}}}\} + \cdots +\exp \{\eta _{J-1}^{{\textbf{x}}}\} + 1\). Note that we provide the expression of \(\pi _J^{{\textbf{x}}}\) for completeness while \(\pi _J^{{\textbf{x}}} = 1 - \pi _1^{{\textbf{x}}} - \cdots - \pi _{J-1}^{{\textbf{x}}}\) is an easier way for numerical calculations.

As for cumulative logit models, the candidate \({{\textbf{x}}}\) must satisfy \(-\infty< \eta _1^{{\textbf{x}}}< \eta _2^{{\textbf{x}}}< \cdots< \eta _{J-1}^{{\textbf{x}}} < \infty \). Otherwise, \(0< \pi _j^{\textbf{x}} < 1\) might be violated for some \(j=1, \ldots , J\). In other words, the feasible design region should be

which depends on the regression parameter \(\varvec{\theta }\) (see Section S.14 in the Supplementary Material of Bu et al. (2020) for such an example). For cumulative logit models, if \({{\textbf{x}}} \in {{\mathcal {X}}}_{\varvec{\theta }}\), then

according to Lemma S.11 of Bu et al. (2020).

Once \(\varvec{\pi }_{{\textbf{x}}} \in {\mathbb {R}}^J\) is obtained, we can calculate \(u_{st}^{{\textbf{x}}} = u_{st}({\varvec{\pi }}_{\textbf{x}})\) based on Theorem A.2 in Bu et al. (2020) as follows:

1. (i)

   \(u_{st}^{{\textbf{x}}} = u_{ts}^{{\textbf{x}}}\), \(s,t=1, \ldots , J\);

2. (ii)

   \(u_{sJ}^{{\textbf{x}}} = 0\) for \(s=1, \ldots , J-1\) and \(u_{JJ}^{{\textbf{x}}} = 1\);

3. (iii)

   For \(s=1, \ldots , J-1\), \(u_{ss}^{{\textbf{x}}}\) is

   $$\begin{aligned} \left\{ \begin{array}{ll} \pi _s^{{\textbf{x}}} (1-\pi _s^{{\textbf{x}}}) &{} \text{ for } \text{ baseline-category },\\ (\gamma _s^{{\textbf{x}}})^2(1-\gamma _s^{{\textbf{x}}})^2((\pi _s^{{\textbf{x}}})^{-1} + (\pi _{s+1}^{{\textbf{x}}})^{-1}) &{} \text{ for } \text{ cumulative },\\ \gamma _s^{{\textbf{x}}}(1-\gamma _s^{{\textbf{x}}}) &{} \text{ for } \text{ adjacent-categories },\\ \pi _s^{{\textbf{x}}}(1-\gamma _s^{{\textbf{x}}})(1-\gamma _{s-1}^{\textbf{x}})^{-1} &{} \text{ for } \text{ continuation-ratio }; \end{array}\right. \end{aligned}$$
4. (iv)

   For \(1\le s < t \le J-1\), \(u_{st}^{{\textbf{x}}}\) is

   $$\begin{aligned} \left\{ \begin{array}{ll} -\pi _s^{{\textbf{x}}} \pi _t^{{\textbf{x}}} &{} \text{ for } \text{ baseline-category },\\ -\gamma _s^{{\textbf{x}}}\gamma _t^{{\textbf{x}}}(1-\gamma _s^{{\textbf{x}}})(1-\gamma _t^{{\textbf{x}}})(\pi _t^{{\textbf{x}}})^{-1} &{} \text{ for } \text{ cumulative }, t-s=1,\\ 0 &{} \text{ for } \text{ cumulative }, t-s>1,\\ \gamma _s^{{\textbf{x}}}(1-\gamma _t^{{\textbf{x}}}) &{} \text{ for } \text{ adjacent-categories },\\ 0 &{} \text{ for } \text{ continuation-ratio }; \end{array}\right. \end{aligned}$$

where \(\gamma _j^{{\textbf{x}}} = \pi _1^{{\textbf{x}}} + \cdots + \pi _j^{{\textbf{x}}}\), \(j=1, \ldots , J-1\); \(\gamma _0^{{\textbf{x}}}\equiv 0\) and \(\gamma _J^{{\textbf{x}}}\equiv 1\).

Example that \({{\textbf{F}}}_{{{\textbf{x}}}} = {{\textbf{F}}}_{{{\textbf{x}}}'}\) with \({{\textbf{x}}} \ne {{\textbf{x}}}'\)

Consider a special MLM (1) with proportional odds (po) (see Section S.7 in the Supplementary Material of Bu et al. (2020) for more technical details). Suppose \(d=2\) and a feasible design point \({{\textbf{x}}} = (x_1, x_2)^T \in [a, b]\times [-c, c] = \mathcal{X}\), \(c > 0\), \(J\ge 2\), \({{\textbf{h}}}_c({{\textbf{x}}}) = (x_1, x_2^2)^T\). Then the model matrix at \({{\textbf{x}}} = (x_1, x_2)^T\) is

Then \(p=J+1\). Let \(\varvec{\theta }= (\beta _1, \ldots , \beta _{J-1}, \zeta _1, \zeta _2)^T \in {\mathbb {R}}^{J+1}\) be the model parameters (since \(\varvec{\theta }\) is fixed, we may assume that \(\mathcal{X} = \mathcal{X}_{\varvec{\theta }}\) if the model is a cumulative logit model). Let \({{\textbf{x}}}' = (x_1, -x_2)^T\). Then \({{\textbf{X}}}_{{\textbf{x}}} = {{\textbf{X}}}_{{{\textbf{x}}}'}\) and thus \(\varvec{\eta }_{{\textbf{x}}} = \varvec{\eta }_{{{\textbf{x}}}'}\). According to (A2) (or (A4)), we obtain \(\varvec{\pi }_{{\textbf{x}}} = \varvec{\pi }_{{{\textbf{x}}}'}\) and then \({{\textbf{U}}}_{{\textbf{x}}} = {{\textbf{U}}}_{{{\textbf{x}}}'}\) . The Fisher information matrix at \({{\textbf{x}}}\) is \({\textbf{F}}_{{\textbf{x}}} = {{\textbf{X}}}_{{\textbf{x}}}^T {{\textbf{U}}}_{{\textbf{x}}} {{\textbf{X}}}_{{\textbf{x}}} = {{\textbf{X}}}_{{{\textbf{x}}}'}^T {\textbf{U}}_{{{\textbf{x}}}'} {{\textbf{X}}}_{{{\textbf{x}}}'} = {{\textbf{F}}}_{{\textbf{x}}'}\). Note that \({{\textbf{x}}} \ne {{\textbf{x}}}'\) if \(x_2 \ne 0\).

First-order derivative of sensitivity function

As mentioned in Sect. 3.1, to apply Algorithm 1 for MLM, we need to calculate the first-order derivative of the sensitivity function \(d({\textbf{x}}, \varvec{\xi })\).

Recall that the first k (\(1\le k\le d\)) factors are continuous. Given \({{\textbf{x}}} = (x_1, \ldots , x_d)^T \in \mathcal{X}\), for each \(i=1, \ldots , k\), according to Formulae 17.1(a), 17.2(a) and 17.7 in Seber (2008),

where

\(\frac{\partial {{\textbf{U}}}_{{\textbf{x}}}}{\partial x_i} = \left( \frac{\partial u^{{\textbf{x}}}_{st}}{\partial x_i}\right) _{s,t=1, \ldots , J}\) with

\({{\textbf{C}}}\) and \({{\textbf{L}}}\) defined as in (1), and \({{\textbf{D}}}_{{\textbf{x}}} = \textrm{diag}({{\textbf{L}}} \varvec{\pi }_{{\textbf{x}}})\). Explicit formula of \(({{\textbf{C}}}^T {{\textbf{D}}}_{{\textbf{x}}}^{-1} {{\textbf{L}}})^{-1}\) can be found in Section S.3 in the Supplementary Material of Bu et al. (2020) with \({{\textbf{x}}}_i\) replaced by \( {{\textbf{x}}}\). As for \(\frac{\partial u^{{\textbf{x}}}_{st}}{\partial \varvec{\pi }_{{\textbf{x}}}^T}\), we have the following explicit formulae

1. (i)

   \(\frac{\partial u_{st}^{{\textbf{x}}}}{\partial \varvec{\pi }_{{\textbf{x}}}} = \frac{\partial u_{ts}^{{\textbf{x}}}}{\partial \varvec{\pi }_{{\textbf{x}}}}\), \(s,t=1, \ldots , J\);

2. (ii)

   \(\frac{\partial u_{sJ}^{{\textbf{x}}}}{\partial \varvec{\pi }_{{\textbf{x}}}} = {{\textbf{0}}} \in {\mathbb {R}}^J\) for \(s=1, \ldots , J\);

3. (iii)

   For \(s=1, \ldots , J-1\), \(\frac{\partial u_{ss}^{{\textbf{x}}}}{\partial \varvec{\pi }_{{\textbf{x}}}}\) is

   $$\begin{aligned} \left\{ \begin{array}{cl} \left( \pi _s^{{\textbf{x}}} {{\textbf{1}}}_{s-1}^T, 1-\pi _s^{{\textbf{x}}}, \pi _s^{{\textbf{x}}} {{\textbf{1}}}_{J-s}^T\right) ^T &{} \text{ for } \text{ baseline-category }\\ u_{ss}^{{\textbf{x}}} \left[ \left( \frac{2}{\gamma _s^{{\textbf{x}}}} {{\textbf{1}}}_s^T, \frac{2}{1-\gamma _s^{{\textbf{x}}}} {{\textbf{1}}}_{J-s}^T\right) ^T\right. &{} \\ \left. - \frac{\pi _{s+1}^{{\textbf{x}}} {{\textbf{e}}}_s}{\pi _s^{{\textbf{x}}} (\pi _s^{{\textbf{x}}} + \pi _{s+1}^{{\textbf{x}}})} - \frac{\pi _s^{{\textbf{x}}} {{\textbf{e}}}_{s+1}}{\pi _{s+1}^{{\textbf{x}}} (\pi _s^{{\textbf{x}}} + \pi _{s+1}^{{\textbf{x}}})}\right] &{} \text{ for } \text{ cumulative }\\ \left( (1-\gamma _s^{{\textbf{x}}}) {{\textbf{1}}}_s^T, \gamma _s^{{\textbf{x}}} {{\textbf{1}}}_{J-s}^T\right) ^T &{} \text{ for } \text{ adjacent-categories }\\ \left( {{\textbf{0}}}_{s-1}^T, \frac{(1-\gamma _s^{\textbf{x}})^2}{(1-\gamma _{s-1}^{{\textbf{x}}})^2}, \frac{(\pi _s^{{\textbf{x}}})^2 {{\textbf{1}}}_{J-s}^T}{(1-\gamma _{s-1}^{{\textbf{x}}})^2}\right) ^T &{} \text{ for } \text{ continuation-ratio } \end{array}\right. \end{aligned}$$

   where \({{\textbf{e}}}_s\) is the \(J\times 1\) vector with the sth coordinate 1 and all others 0, \({{\textbf{1}}}_s\) is the \(s\times 1\) vector of all 1, and \({{\textbf{0}}}_s\) is the \(s\times 1\) vector of all 0.

4. (iv)

   For \(1\le s < t \le J-1\), \(\frac{\partial u_{st}^{{\textbf{x}}}}{\partial \varvec{\pi }_{{\textbf{x}}}}\) is

   $$\begin{aligned} \left\{ \begin{array}{cl} \left( {{\textbf{0}}}_{s-1}^T, -\pi _t^{{\textbf{x}}}, {{\textbf{0}}}_{t-s-1}^T, -\pi _s^{{\textbf{x}}}, {{\textbf{0}}}_{J-t}^T\right) ^T &{} \text{ for } \text{ baseline-category }\\ \left( -(1-\gamma _s^{{\textbf{x}}})(1-\gamma _t^{{\textbf{x}}})\left( 1 + \frac{2\gamma _s^{{\textbf{x}}}}{\pi _t^{{\textbf{x}}}}\right) {{\textbf{1}}}_s^T, \right. &{} \\ -\gamma _s^{{\textbf{x}}} (1-\gamma _t^{{\textbf{x}}}) \left[ 1 - \frac{\gamma _s^{{\textbf{x}}} (1-\gamma _t^{{\textbf{x}}})}{(\pi _t^{{\textbf{x}}})^2}\right] , &{} \\ \left. -\gamma _s^{{\textbf{x}}} \gamma _t^{{\textbf{x}}} \left[ 1 + \frac{2(1-\gamma _t^{{\textbf{x}}})}{\pi _t}\right] {{\textbf{1}}}_{J-s-1}^T\right) ^T &{} \text{ for } \text{ cumulative }, t-s=1\\ {{\textbf{0}}}_J &{} \text{ for } \text{ cumulative }, t-s>1\\ \left( (1-\gamma _t^{{\textbf{x}}}) {{\textbf{1}}}_s^T, {{\textbf{0}}}_{t-s}^T, \gamma _s^{{\textbf{x}}} {{\textbf{1}}}_{J-t}^T\right) ^T &{} \text{ for } \text{ adjacent-categories }\\ {{\textbf{0}}}_J &{} \text{ for } \text{ continuation-ratio } \end{array}\right. \end{aligned}$$

where \(\gamma _j^{{\textbf{x}}} = \pi _1^{{\textbf{x}}} + \cdots + \pi _j^{{\textbf{x}}}\), \(j=1, \ldots , J-1\); \(\gamma _0^{{\textbf{x}}}\equiv 0\) and \(\gamma _J^{{\textbf{x}}}\equiv 1\).

Thus the explicit formulae for \(\frac{\partial d({{\textbf{x}}}, \varvec{\xi })}{\partial x_i}\), \(i=1, \ldots , k\) can be obtained via (C5). Only \(\frac{\partial {\textbf{X}}_{{\textbf{x}}}}{\partial x_i}\) is related to i, which may speed up the computations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions