Optimal designs for nonlinear mixed-effects models using competitive swarm optimizer with mutated agents

Nature-inspired meta-heuristic algorithms are increasingly used in many disciplines to tackle challenging optimization problems. Our focus is to apply a newly proposed nature-inspired meta-heuristics algorithm called CSO-MA to solve challenging design problems in biosciences and demonstrate its flexibility to find various types of optimal approximate or exact designs for nonlinear mixed models with one or several interacting factors and with or without random effects. We show that CSO-MA is efficient and can frequently outperform other algorithms either in terms of speed or accuracy. The algorithm, like other meta-heuristic algorithms, is free of technical assumptions and flexible in that it can incorporate cost structure or multiple user-specified constraints, such as, a fixed number of measurements per subject in a longitudinal study. When possible, we confirm some of the CSO-MA generated designs are optimal with theory by developing theory-based innovative plots. Our applications include searching optimal designs to estimate (i) parameters in mixed nonlinear models with correlated random effects, (ii) a function of parameters for a count model in a dose combination study, and (iii) parameters in a HIV dynamic model. In each case, we show the advantages of using a meta-heuristic approach to solve the optimization problem, and the added benefits of the generated designs.

Drs. Wong and Zhang were partially supported by a grant from the National Institute of General Medical Sciences of the National Institutes of Health under Award Number R01GM107639. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. Dr. Wong was also partially supported by the Yushan Scholarship Award from the Ministry of Education of Taiwan.

Authors and Affiliations

1. Department of Biostatistics, UCLA, 650 Charles E Young Dr S, Los Angeles, CA, 90095, USA

   Elvis Han Cui, Zizhao Zhang & Weng Kee Wong

2. Alibaba Group, 969 West Wen Yi Road, Yuhang District, Hangzhou, 311121, Zhejiang, China

   Zizhao Zhang

Contributions

Cui and Zhang wrote the main manuscript text and Wong polished the writing. All authors reviewed the manuscript. Zhang prepared Figure 1 and Tables 1-3; wrote the CSOMA codes in Matlab. Cui prepared Figure 2 and Tables 4-5; wrote the CSOMA codes in Python.

Corresponding author

Correspondence to Elvis Han Cui.

Conflict of interest

The author declares no competing interest.

Publisher's Note

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

Appendices

This appendix has 5 parts. In Appendix A, we list notations used in the paper and give their explanation. Appendix B provides an overview of nature-inspired meta-heuristic algorithms, review the CSO algorithm and show how CSO-MA is derived as a variant of CSO. In Appendix C, we mention online tools available to use CSO-MA and generate optimal designs for the examples in this paper and beyond. In Appendix D, we provide an illustrative derivation of the information matrix for a negative binomial model with a random intercept and a single factor. Extension to more than one factor is quite straightforward. Appendix E provides additional information to fill in some technical details for Application 6. In addition, two algorithms are included in the appendix. Algorithm 1 describes how to use CSO-MA to find optimal designs and Algorithm 2 provides further details on finding Bayesian optimal designs for the nonlinear model used to study HIV dynamics.

Appendix A Table of notations

Appendix B Nature-inspired meta-heuristic algorithms

Nature-inspired meta-heuristic algorithms, such as swarm-based algorithms, now form a dominant component in the optimization literature. Whitacre (2011a, 2011b) documented their meteoric rise in popularity in solving real, high-dimensional and complex optimization problems, and their interest and growth has continued unabated. They are already widely used in engineering and computer science, and recently also in other disciplines for tackling high-dimensional and complex real optimization problems. There are many reviews of such algorithms and their many multi-faceted applications to solve very different types of challenging optimization problems; see, for example, Korani and Mouhoub (2021), and Kashif et al. (2019), among several others. They have several appealing features over current algorithms and the main ones are that their codes are widely and freely available, they are easy to implement and use; they are fast, essentially assumptions free and are general-purpose optimization algorithms. While they do not guarantee that they will find the optimal solution, they frequently do so or find a nearly optimal solution quickly. Some recent studies have shown that swarm-based algorithms can search for hard-to-find optimal designs that were previously deemed difficult. For instance, Chen et al. (2017) found optimal designs for standardized maximin optimal designs for nonlinear models and Xu et al. (2019) found optimal designs for logistic models with many interacting factors.

There are many monographs on nature-inspired meta-heuristic algorithms at various levels; for instance, Yang (2010) monograph is at an introductory level and Heaton (2019) is at a bit more advanced level. Some monographs target specific disciplines; for example, Mendes et al. (2019) focus on agricultural applications and Sharma and Kaur (2021) is interested in feature selection problems and has numerous applications in finance and reinforce learning. meta-heuristics has undoubtedly emerged as a set of dominant optimization tools in industry and academia (Whitacre 2011a) and in artificial intelligence research as well Ezugwu and Shukla (2021). Their applications are plentiful and have been shown capable of solving different types of complex and high-dimensional optimization problems (Khan et al. 2019; Hesami et al. 2020; Cui et al. 2022). Websites like R, Python and MATLAB codes are available at youtube.com, github.com, including at professional websites, such as Matlab.com and SAS.com. For example, the website https://pyswarms.readthedocs.io/en/latest/ has a set of PSO tools written in Python (Miranda 2018). See also https://github.com/ElvisCuiHan/CSOMA, where codes are available for using CSO-MA to solve different types of estimation and design problems in statistics, including imputing missing data optimally or in the optimal selection of variables in ecology. Haidar et al. (2021) provided a high-level Python package for selecting machine learning algorithms and their parameters using PSO.

1.1 B.1 Competitive swarm optimizer

Competitive swarm optimizer (CSO) is a relatively novel swarm-based algorithm that has been proven to be very effective in solving different types of optimization problems (Cheng and Jin 2015). CSO has had successful applications to solve large and hard optimization problems. For example, Gu et al. (2018) applied CSO to select variables for high-dimensional classification models, and Xiong and Shi (2018) used CSO to study a power system economic dispatch, which is typically a complex nonlinear multivariable strongly coupled optimization problem with equality and inequality constraints.

Competitive swarm optimizer algorithm, or CSO for short, minimizes a given function \(f(\textbf{x})\) over a user-specified compact space \(\varvec{\Omega }\) by first generating a set of candidate solutions. In our case, they take the form of a swarm of n particles at positions \(\textbf{x}_1, \;\ldots , \;\textbf{x}_n\), along with their corresponding random velocities \(\textbf{v}_1, \;\ldots , \;\textbf{v}_n\).

After the initial swarm is generated, at each iteration we randomly divide the swarm into \(\left\lfloor \frac{n}{2} \right\rfloor \) pairs and compare their objective function values. We identify \(\textbf{x}^t_i\) as the winner and \(\textbf{x}^t_j\) as the loser if these two are competed at the iteration t and \(f(\textbf{x}^t_i) < f(\textbf{x}^t_j)\). The winner retains the status quo and the loser learns from the winner. The two defining equations for CSO are

where \(\textbf{R}_1, \;\textbf{R}_2, \;\textbf{R}_3\) are all random vectors whose elements are drawn from U(0, 1); operation \(\otimes \) also represents element-wise multiplication; vector \(\bar{\textbf{x}}^t\) is simply the swarm center at iteration t; social factor \(\phi \) controls the influence of the neighboring particles to the loser and a large value is helpful for enhancing swarm diversity (but possibly impacts convergence rate). This process iterates until some stopping criteria are met.

Simulation results have repeatedly shown that CSO either outperforms or is competitive with several state-of-the-art evolutionary algorithms, including several enhanced versions of PSO. This conclusion was arrived at after comparing CSO performance with state-of-the-art EAs using a variety of benchmark functions with dimensions up to 5000 and found that CSO was frequently the fastest and winner in terms of the quality of the solution (Cheng and Jin 2015; Mohapatra et al. 2017; Sun et al. 2016; Zhang et al. 2016, 2017).

1.2 B.2 Competitive swarm optimizer with mutated agents

Our experience with CSO is that many of the generated design points are at the boundary of the design space \(\varvec{\Omega }\). This is helpful because many optimal designs, especially D-optimal designs, have their support points at the boundary \(\varvec{\Omega }\).

To this end, we propose an improvement on CSO and call the enhanced version, competitive swarm optimizer with Mutated Agents or, in short, CSO-MA. After pairing up the swarm in groups of two at each iteration, we randomly choose a loser particle p as an agent, randomly pick a variable indexed as q and then randomly change the value of \(\textbf{x}_{pq}\) to either \(\textbf{xmax}_{q}\) or \(\textbf{xmin}_q\), where \(\textbf{xmax}_q\) and \(\textbf{xmin}_q\) represent, respectively, the upper bound and lower bound of the q-th variable. If the current optimal value is already close to the global optimum, this change will not hurt since we implement this experiment on a loser particle, which is not leading the movement for the whole swarm; otherwise, this chosen agent restarts a journey from the boundary and has a chance to escape from a local optimum. In an algorithmic form, the proposed CSO-MA is given below.

CSO-MA Algorithm

Full size image

The computational complexity of CSO is \(\mathcal {O}(nD)\), where n is the swarm size and D is the dimension of the problem. Since our modification only adds one coordinate mutation operation to each particle, its computational complexity is the same as that of CSO. The improved performance of CSO-MA over CSO-MA to find the optimum for many complex multi-dimensional benchmark functions has been validated (Zhang et al. 2020).

Appendix C On the online CSOMA package

The CSOMA algorithm was originally written in Matlab and is available online at https://github.com/ElvisCuiHan/CSOMA/tree/main/CSOMA%20for%20Matlab. We extend it to Python, utilizing the pyswarms package (Miranda 2018) and it is available at https://github.com/ElvisCuiHan/CSOMA/tree/main/CSOMA%20for%20Python. On the same page, an example for finding D-optimal design in a trinomial dose-response model is given.

In addition, codes for application 3 and 7 are also given in https://github.com/ElvisCuiHan/CSOMA/ and users can reproduce the results in the paper by running the run_fractional.m (in the folder High_Dim_D_Optimal_Example) and run_glm_fisher.m (in the folder Fractional_Polynomial_Regression_Example) files in Matlab, respectively.

Appendix D Derivation of the information matrix

Suppose we have a negative binomial model with a random intercept and a single fixed factor. Suppose further that there are n subjects, each with T independent observations and \(\mu _{ij}\) is the mean of the count outcome \(y_{ij}\). Then

Lawless (1987) showed that the total likelihood function for such a mixed model is

where a is the dispersion factor. Apart from an unimportant additive constant, the log-likelihood is

Noting that

we calculate the following quantities for the Fisher information matrix:

The information matrix for a mixed negative binomial model with a random intercept and multiple fixed factors can be derived similarly.

Appendix E Details on the Bayesian optimal design for nonlinear mixed models applied to HIV dynamics

In this subsection, we give details on computing the utility function \(\textbf{E}[\textbf{Var}(\mu _{2}|\textbf{y})]\) for the HIV dynamics model. Details for computing \(\textbf{E}[\textbf{Var}(\mu _{3}|\textbf{y})]\) are similar and therefore omitted. Mathematically, the utility function can be written as

Recall that \(\varvec{\mu }\sim \mathcal {N}(\varvec{\eta },\varvec{\Lambda })\), so the first term does not involve \(t_j\) because the second moment of \(\mu _2\) is just \(\eta _2^2+\varvec{\Lambda }_{22}\) where \(\eta _2\) is the second element of \(\varvec{\eta }\) and \(\varvec{\Lambda }_{22}\) is the \((2,2)^{th}\) element of \(\varvec{\Lambda }\). Hence optimization with respect to \(t_j\) only involves the second term and according to theorem 1 in Han and Chaloner (2004), it is finite. It requires the computation of

Recall that for fixed i and j, we have \(f(y_{ij}|\varvec{\mu }, \varvec{\Sigma }^{-1}, \sigma ^{-2}) = \int f(y_{ij}|\varvec{\theta }_i,t_j, \sigma ^{-2})f(\varvec{\theta }_i|\varvec{\mu }, \varvec{\Sigma })d\varvec{\theta }_i,\) where \(\varvec{\theta }_i=(\log V_{0i},\log c_i,\log \delta _i)^T\) is defined in equation (6) in the paper. The density \(f(\textbf{y}_{i}|\varvec{\mu },\varvec{\Sigma }^{-1},\sigma ^{-2})\) can be approximated as follows (for a fixed design \(t_{ij},j=1,\ldots ,T\)):

where T is the number of time points (T=16 in our paper) and K is the number of quadratures. We found that \(K=5\) is sufficient to obtain stable results. We denote the above approximation by \(F_1\)(\(\textbf{t}\), \(\textbf{y}_{i}\), \(\varvec{\mu }\), \(\varvec{\Sigma }^{-1}\), \(\sigma ^{-2}\), K). Next, we use a sampling scheme to derive i.i.d. samples \((\varvec{\mu }_i,\varvec{\Sigma }_i^{-1},\sigma _i^{-2})\) for \(i=1,\cdots ,s\), and hence

where \(N=5\) gives stable results and

and

We denote the above approximation (E4) as \(F_2\)(\(\textbf{t}\), \(\textbf{y}_i\), K, \(\varvec{\Phi }\), \(\gamma \), \(\varvec{\eta }\), \(\varvec{\Lambda }\), \(\alpha \), \(\beta \), N). We then draw S different y samples to approximate the outer expectation:

To sum up, we use the following algorithm to generate a Bayesian design:

• We start with a sequence of random y, say, generated uniformly from \(-\,10\) to 10.

• Given the initial \(\textbf{y}\), we apply CSO-MA to find the optimal \(\textbf{t}\).

• Given the optimal \(\textbf{t}\) we obtain in Step 2, we resample \(\textbf{y}\).

• Iterate between Steps 2 and 3 until some user-specified stopping criterion is met. In our case, we run 30 times in total and we set \(T=16\), \(N=5, K=5\) and \(S=50\).

Bayesian Optimal Design for Nonlinear Mixed Models to Study HIV Dynamics

Full size image

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