G-optimal designs for hierarchical linear models: an equivalence theorem and a nature-inspired meta-heuristic algorithm.

Hierarchical linear models are widely used in many research disciplines and estimation issues for such models are generally well addressed. Design issues are relatively much less discussed for hierarchical linear models but there is an increasing interest as these models grow in popularity. This paper discusses the G-optimality for predicting individual parameters in such models and establishes an equivalence theorem for confirming the G-optimality of an approximate design. Because the criterion is non-differentiable and requires solving multiple nested optimization problems, it is much harder to find and study G-optimal designs analytically. We propose a nature-inspired meta-heuristic algorithm called competitive swarm optimizer (CSO) to generate G-optimal designs for linear mixed models with different means and covariance structures. We further demonstrate that CSO is flexible and generally effective for finding the widely used locally D-optimal designs for nonlinear models with multiple interacting factors and some of the random effects are correlated. Our numerical results for a few examples suggest that G and D-optimal designs may be equivalent and we establish that D and G-optimal designs for hierarchical linear models are equivalent when the models have only a random intercept only. The challenging mathematical question of whether their equivalence applies more generally to other hierarchical models remains elusive.

Introduction

Hierarchical models are widely used to analyze data in various disciplines, such as in psychology, medicine, manufacturing industry and education. Such models are especially appealing for analyzing longitudinal analysis because they allow for the presence of missing data, time-varying or invariant covariates, and subjects measured at different time points. In educational research, hierarchical models are commonly used to evaluate the effectiveness of teaching methods using data from students nested within classrooms and classrooms are nested within schools that use different teaching methods. The distinguishing feature of these models is that they account for both individual-level and population-level effects. In the literature across disciplines, hierarchical models are variously referred to as multilevel models, nested data models, mixed models, random coefficient, random-effects models or random parameter models. In our work, we make no distinction among them and refer to them as statistical models with parameters that vary at one or more levels.

When the experimental settings are under the control of the investigator, design issues arise and they must be carefully addressed for maximal precision in the statistical inference at a minimal cost. The basic questions to answer are given a statistical model defined on a compact design space , an optimality criterion and a fixed amount of resources to take N observations for the study, what are the optimal number of points, where these points are and how many replicates to take at each of these points. We denote these quantities by k, and , respectively where . Optimal exact designs are challenging to find and study because there are no theoretical tools for finding them in general and in addition to depending on the design criterion and the model, the optimal exact design also depends on the value of N.

An alternative option is to formulate the design problem to find an optimal approximate design, where we determine k, as before, and now the optimal weights defined by with . In practice, there are implemented by taking rounding each to the nearest integer and taking observations at subject to . Approximate designs were proposed in Kiefer (1959) and they are appealing because when the criterion is a concave functional, there is a theoretical tool called an equivalence theorem for confirming the optimality of an approximate design among all designs on the given design space . There are algorithms with proof of convergence for searching some types of optimal designs and there are also tools to assess the proximity of an approximate design from the optimum without knowing the latter. For these reasons, we focus on approximate designs in the rest of the paper.

A primary interest of our paper is to find G-optimal designs for hierarchical linear models, with G standing for global. These designs are best for estimating the overall response surface with minimal variance across the design space and so protect against the worst-case scenario. This is unlike the much simpler situation where there is interest in finding an optimal design for estimating the predicted response at a single point. G-optimality is a minimax-type of design criterion, it is non-differentiable and requires at least solving two nested layers of optimization problems over different search spaces. Consequently, even for fixed-effects models, G-optimal designs are among the most difficult to study mathematically. Compounding the optimization problem is that current algorithms for searching optimal designs cannot find G-optimal designs effectively, let alone one with proof of convergence.

The main aims of this paper are to find G-optimal designs for hierarchical linear models, propose a nature-inspired meta-heuristic algorithm to find them and develop an equivalence theorem to confirm the G-optimality of an approximate design for hierarchical linear models. We demonstrate the flexibility and effectiveness of the algorithm for finding G-optimal designs when the mean response is modeled by fractional polynomials and D-optimal designs for hierarchical nonlinear models with multiple interacting factors with possibly correlated random effects. We also establish that D and G-optimal designs for hierarchical linear models are equivalent when the models have only a random intercept only and pose the challenging mathematical question of whether their equivalence applies more generally to other hierarchical models.

The rest of the paper is organized as follows. Section 2 provides the background of hierarchical linear models and recent literature on constructing D-optimal designs for such models. Section 3 establishes an equivalence theorem for G-optimality and Sect. 4 presents G-optimal designs for various types of hierarchical models for relatively simple models. To find G-optimal designs for more complicated models, Sect. 5 introduces a nature-inspired meta-heuristic algorithm called competitive swarm optimizer (CSO) to find G-optimal designs for more complicated hierarchical linear models, such as, when the mean function is a fractional polynomial and the random effects may be correlated. We further show CSO can find locally D-optimal designs for estimating model parameters in a Poisson hierarchical model with multiple factors defined on a user-specified design space and some random effects are correlated. Section 6 concludes with a discussion of possible equivalence between D and G-optimal designs for linear mixed models.

Preliminaries and model specification

This section gives the background, a brief literature review on constructing optimal designs for hierarchical linear models before we describe our statistical models.

Equivalence theorem

Constructing optimal designs for a given model under a specified criterion is challenging if the model is complex, and more so if the criterion is complicated. An analytical approach for finding optimal designs is limiting and is possible only for simple models with a couple of parameters, see examples in design monographs, such as, Fedorov (1972) and Berger and Wong (2009). A common practice is to find an optimal design under a set of restrictive conditions and hope that the design remains optimal under a broader setup. For instance, one may find the best two-point optimal design for the simple logistic model. Is the same design still optimal among all three-point designs? Equivalence theorem provides an answer, and more generally, is able to confirm whether the design is optimal among all designs of interest.

Specifically, suppose the given design space is compact and the optimality criterion is formulated as a concave (or convex) function of approximate designs or, equivalently, as a function of the normalized Fisher Information matrix . The elements of this matrix are the expectation of the negative of the second derivatives of the log-likelihood function with respect to the parameters. As an example, suppose we have a standard linear regression model with the mean responsewhere is a univariate function and is the vector of unknown parameters. The expectation is over the model errors assumed to have independent normal variates with zero means and constant variance. A direct calculation shows the Fisher Information matrix is proportional toA typical design goal is find a design so that the estimated has a minimal covariance matrix. The most common design criterion is D-optimality and a D-optimal design satisfiesover all approximate designs on . Here is the determinant of and is the set of all approximate designs on . After the D-optimal design is found, the worth of another design in terms of D-optimality is measured by its D-efficiency defined byClearly, the above ratio is between 0 and 1 and if it is equal to one half, the interpretation is that the design has to be replicated twice to do as well as the D-optimal design . If in the denominator is replaced by another design , the ratio is the D-efficiency of relative to . The efficiency or relative efficiency of a design under another criterion is similarly defined.

The function is concave on and using a standard convex analysis argument (Fedorov 1972; Pazman 1986), it can be shown that the equivalence theorem for D-optimality is as follows: is D-optimal if and only if for all ,The function on the left-hand side of the inequality is frequently called the sensitivity function of the design. As an example, if , the design space is and we want to estimate all the three parameters as accurately as possible. A direct calculation shows the D-optimal design takes equal proportion of observations at and and one can directly verify its sensitivity function satisfies (1).

Each concave functional has a directional derivative which is used to derive its own unique equivalence theorem; see details in design monographs, such as Fedorov (1972) and Berger and Wong (2005). When is an interval or a two-dimensional space, the optimality of a design can be readily checked by plotting the sensitivity function on the left-hand side of the above inequality across the design space and visually ascertain whether the conditions in the equivalence theorem are satisfied. Equivalence theorems are frequently used to confirm the optimality of an approximate design, and more importantly, to ascertain the efficiency of any approximate design. An example is Wong and Cook (1993), where they proposed an algorithm to find G-optimal designs for fixed-effects linear models and assess the quality of the generated design using an equivalence theorem.

Brief review of optimal designs

There is a lot of work on finding optimal designs for various models in different fields and most concerned finding D-optimal designs for linear models; see references in Berger and Wong (2005) and those cited below. For nonlinear models, the information matrix depends on the unknown parameters and nominal values are required to replace them in the information matrix before the criterion is optimized. Because the resulting designs depend on the nominal values, they are called locally D-optimal designs and are the simplest to construct for nonlinear models and commonly used in practice when a single best guess for the parameters is available. However, where there are conflicting opinions from experts or prior information from previous studies, locally optimal designs become problematic to implement.

Two common design approaches for such a scenario are to adopt a Bayesian paradigm or a minimax approach. The latter design strategy was used by Berger et al. (2000) where they searched for the best D-optimal design over a plausible set of values for the nominal values. The minimax D-optimal design from Berger et al. (2000) minimizes the maximum inefficiencies that arise from misspecified values in the plausible set; variations of the theme are possible; see Chen et al. (2018). An appeal of this approach is that practitioners are likely able to provide a set of plausible nominal values for the model parameters than having to provide a prior distribution to implement the Bayesian optimal design. A disadvantage of the minimax approach is that minimax optimal designs are both theoretically and numerically more challenging to find than Bayesian optimal designs.

Random-effects models have more complicated equivalence theorems and a limited number of them are available. For instance, Fedorov and Hackl (1972) derived an equivalence theorem to confirm the D-optimality of an approximate design for a random coefficient regression model and Entholzner et al. (2005) obtained optimal or efficient designs for mixed models. Schmelter (2007b) showed that the search for an optimal design for hierarchical linear models could be restricted to the class of group-wise identical individual designs and Schmelter (2007a) noted that optimal designs found in the class of single-group designs remain optimal in a larger class having more group designs. Entholzner et al. (2005), Schmelter (2007a) and Schmelter (2007b) investigated optimal designs under random intercept models, random slope models and random coefficient cubic regression models, respectively. Debusho and Haines (2008) and Debusho and Haines (2011) constructed V- and D-optimal population designs for linear and quadratic regression models with a random intercept term. More recently, Prus and Schwabe (2016b) constructed optimal designs for predicting individual parameters in hierarchical models and called them D-optimal designs. They also extended their work to find interpolation and extrapolation optimal designs for random coefficient regression models in Prus and Schwabe (2016a). Further, Prus (2019) proposed a design criterion called G-optimality for predicting the response surface in hierarchical models and noted technical difficulties in finding such an optimal design, analytically or numerically, for more complicated models.

Statistical models

Throughout, we denote the jth response from the ith subject by . The hierarchical model is defined on a user-defined compact design space and given bywhere the jth observation from individual i is taken at the experimental setting , n is the number of individuals, is the number of observations from individual i, is the vector of known regression functions, and is the individual parameter vector specifying the individual response. The observational errors are assumed to be centred with zero mean and homoscedastic and uncorrelated with common variance var. We assume that E()= and Cov()= (), where is known, and all ’s are uncorrelated with all ’s. We note that the matrix can be singular, which happens when some of the individual parameters are non-random.

In practice, experiments are usually conducted with identical regimes for all individuals, i.e., all individuals i have the same number of observations at the same values of the experimental settings. Such designs are popular because they are simple to implement and analyze (Prus and Schwabe 2016b). In what is to follow, we assume such a setting.

Let be the vector of observations from individual i, let be the average response across all individuals and let be the individual common design matrix for all individuals. Prus and Schwabe (2016b) showed that the best linear unbiased estimator of the population parameter isand the best linear unbiased predictor of the individual parameter is a weighted average of the individualized estimate based on the observations of subject i only and the estimator for the population parameter. Specifically, we haveThe quality of the prediction in (4) for may be measured by the mean squared error matrix of :Here is the identity matrix, is the matrix with all entries equal to 1 and “” denotes the Kronecker product of matrices.

An equivalence theorem for G-optimality

G-optimality for predicting individual parameters

The G-optimality criterion was proposed for prediction of the individual parameters in Prus (2019) and a G-optimal design minimizes the maximal prediction mean squared error over the experimental region. The standardized information matrix of an approximate design for the above model without individual effects isThe matrix stands for the information obtained per observation and corresponds to the information contributed by the observations at the experimental settings per individual. In this paper, we are only concerned with approximate designs on with non-singular information matrices and denote this set by . If , its MSE-matrix in (5) iswhere . When is non-singular, this expression simplifies toFor predicting purposes, Prus (2019) used model (2) and proposed the G-criterion as the maximal sum of the expected squared differences of the predicted and real response across all individuals with respect to all possible observational settings:This criterion can be rewritten aswhereIn the appendix, we show that the criterion is a convex function on and a design that minimizes the G-criterion over all designs is G-optimal. The next two analytical results give a lower bound for the function for any design and an equivalence theorem to confirm the optimality of a design. The technical arguments are deferred to the appendix.

Theorem 1

For any design , we have

An implication of the above theorem is that a design is G-optimal for model (2) if it satisfies (10) with equality at the support points. In what is to follow, we establish below an equivalence theorem for G-optimality to characterize G-optimal designs.

Theorem 2

For model (2), let , let , let , and letwhere is a probability measure on defined bywith

A design on is G-optimal if and only if there exists a probability measure on such that for all

Similar to (1), the function on the left hand side of the inequality is the sensitivity function of the design . The equivalence theorem is more complex than the one for D-optimality because the G-optimality criterion is not differentiable. The theorem has a form similar to those for confirming G-optimal designs for heteroscedastic linear models in Wong and Cook (1993).

In Theorem 2, the probability measure exists only if the design under investigation is G-optimal, whereupon the inequality becomes an equality at all design points of the optimal design . This follows because if there is a support point x of such that the inequality in Theorem 2 is strict, then integrating both sides with respect to , we haveimplying , which is impossible. It follows that all support points of are roots of the sensitivity function. In the next section, we show how this information is used to generate the G-optimal design.

G-optimal designs

We now provide examples of G-optimal designs for several types of linear mixed models and use Theorem 2 to confirm their optimality. Here and in the rest of the paper, to fix ideas, we assume that but other values for n and m can be similarly used. The first example is relatively simple and a formula for the G-optimal is available. When the model is slightly generalized to two additive factors with uncorrelated random factors and no intercept term, an analytical description for the G-optimal design for two-parameter models is no longer possible and numerical methods must be used. The implication is that G-optimal designs are difficult to find and they have to be found numerically. This motivates us to use competitive swarm optimizer (CSO) to find G-optimal designs and note that some of the results in this section are found by CSO.

Examples

G-optimal design for the simple linear model with a random slope

Suppose we wish to find the G-optimal design for model (2) with and on the experimental region . The model isand we assume that the slope parameter is random with mean and the dispersion matrix is diagonal and equal to .

Let and let be the two-point design supported at 0 and 1 with weight at 1 equal to . A direct application of Theorem 2 shows that the G-optimal design for model (12) hasand for the design , (9) becomesThis maximum of this function is attained at and and so Let be the two-point probability measure defined on supported at 0 and 1 and its weight at 1 isThe sensitivity function of the design in terms of and iswith a maximum value of at and . This value can be shown to equal toand so by Theorem 2, the design is G-optimal.

G-optimal designs for a two additive factor model without an intercept

Regression models with no intercept are quite common and they either arise naturally or from constraints imposed on the variables, see for example, Huang et al. (1995). For such a model with two additive random factors, we havewhere . The dispersion matrix of is diagonal and equal to . To find G-optimal designs, first consider designs of the formand let be the associated probability measure on defined byThe table below shows the optimal weights and the weights of the measure for model (16) for selected values of and These designs have been verified to be numerically G-optimal by Theorem 2.

G-optimal designs for quadratic and fractional polynomial mixed models

We show here that even for a relatively simple mixed model, an analytic description of the G-optimal design can be problematic. We consider two linear mixed models with independent errors having zero means and common variance. These two models areandModel (17) is quadratic and model (18) is an example of a fractional polynomial, which is increasingly used in the biomedical sciences to model a univariate continuous response. Fractional polynomial models were proposed by Royston and Altman (1994) who showed that they are more effective for modeling a continuous outcome than using polynomials (Royston et al. 1999; Royston and Sauerbrei 2008). We recall a fractional polynomial (FP) is given by , where are the real-valued coefficients and are defined sequentially, ,The powers are given by the Box-Tidwell transformation with if , otherwise and for practical applications, ‘powers’ in a FP are selected from the set (Royston and Altman 1994). Many software statistical packages now provide an option for fitting FP models, suggesting that FP models are gaining recognition as a modelling tool in statistics. Interestingly, optimal designs for FP models have never been reported in the literature.

There are no closed-form descriptions for the optimal designs for these two relatively simple models with random effects. A practical way is to find them using an algorithm. There are many traditional algorithms for generating many types of optimal designs and our experience is that many of them do not work well in such a setting where the criterion is not differentiable and the optimization problem has two levels. To this end, we use a meta-heuristic algorithm called competitive swarm optimizer to find the optimal designs.

Competitive swarm optimizer

An analytical approach is generally unable to determine a G-optimal design and we need an effective algorithm to generate a G-optimal design, or search for a design with sufficiently high G-efficiency for practical applications. The last few G-optimal designs in Sect. 4 were found by an algorithm that we now describe.

Nature-inspired meta-heuristic algorithms are commonly used in engineering and computer science to tackle hard-to-solve optimization problems (Yang 2010). Examples are particle swarm optimization (PSO), differential evolutionary (DE), cuckoo search (CS) and imperialist competitive algorithm (ICA). Their main appeal is that they are general purpose optimization tools, they tend to be assumptions-free, easy to implement and use and frequently able to find high-quality solutions quickly. Their meteoric rise in popularity is well documented in Whitacre (2011a, 2011b) with reasons. Recently, Chen et al. (2018), Phoa et al. (2016), Kim and Wong (2018), Masoudi et al. (2019), and Storn and Price (1997) used meta-heuristic algorithms to tackle different types of optimal design problems. For example, Chen et al. (2018) applied a version of PSO to find standardized maximin optimal designs for several enzyme-kinetic inhibition models by solving multilevel nested optimization problems over different types of search spaces, and Kim and Wong (2018) likewise applied PSO and solved an adaptive clinical trial design problem by solving a complex discrete optimization problem by determining the optimal choice of ten integers with multiple constraints.

Evolutionary algorithms are continuously evolving and nature-inspired meta-heuristic algorithms are a major component of evolutionary algorithms. Competitive swarm optimizer (CSO) is popular because many simulation results in the literature show that it either outperforms or is competitive with several state-of-the-art evolutionary algorithms. This conclusion was arrived at after comparing CSO’s performance with several state-of-the-art evolutionary algorithms using a variety of benchmark functions with dimensions up to 5000 (Cheng and Jin 2014; Zhou et al. 2016; Sun et al. 2018; Mohapatra et al. 2017; Zhang et al. 2016). They showed that CSO was frequently not only the winner but also required significantly less run time. CSO has also been successfully applied to solve many different types of complex optimization problems; see, for example, Gu et al. (2018), Kumarappan and Arulraj (2016) and Xiong and Shi (2018).

CSO was initially proposed by Cheng and Jin (2014) to tackle the premature convergence issue met by many evolutionary algorithms. CSO first generates a swarm of n particles at positions with random velocities in . In each iteration, we randomly divide them into pairs and compare their objective function values. If we have a minimization problem, at iteration t, we compare each pair of particles and and identify as the winner and as the loser if the objective function has a smaller value at than at . Winner retains status quo and the loser learns from the winner. The two defining equations for CSO arewhere are all random vectors whose elements are drawn from U(0, 1); operation also represents element-wise multiplication; vector is simply the swarm center at iteration t; social factor controls the influence of the neighboring particles to the loser and a large value is helpful for enhancing the swarm diversity (but possibly impacts convergence rate). This process iterates until some stopping criteria are met. Algorithm 1 displays a pseudo code of CSO.

The tuning parameters we used in the CSO algorithm for finding the optimal designs are similar to those suggested by Cheng and Jin (2014). For example, because the optimization problems we have dimensions fewer than 100, we set and used 128 particles in the search. This is because the examples to follow have about 100 variables or fewer to optimize. The maximum number of iterations for each search was set to be 350. For our examples, the algorithm typically converges in 200 or fewer iterations or less than 1 second of CPU time. The hardware we used is a Windows PC with 3.20GHz Intel i7-8700 CPU, 32GB DDR4 2666MHz memory and 512G SSD storage. Here convergence means that successive values of the objective function do not differ by less than in absolute value. On average, we observe that each search usually converged in 200 or fewer iterations.

Application of CSO to find G-optimal designs

To search for a G-optimal design, we set up a two-layer optimization structure for min-max the objective function. The inner optimization step is a low-dimensional maximization problem represented in formula (9) and the outer loop is to minimize formula (8), which, in our cases, is a multi-dimensional function with less than 100, or more specifically, around 10–20 variables to optimize. Accordingly, we set , 32 particles, 128 iterations for the inner optimization and , 128 particles and 200 iterations for the outer optimization task. During the search, the target matrix may sometimes become singular or close to being singular and affect the floating-point calculation accuracy. We found adding a diagonal matrix with very small positive diagonal elements to the original matrix can make inverting an ill-conditioned matrix more stable. In our work, we set

After finding a design , we determine its answering set and the probability measure that meet the conditions in the equivalence theorem. To efficiently find all x’s that maximize , we suggest to split the design space into two or more subspaces and then run CSO on each subspace to search for all x’s that maximize . For instance, if the design space , we may split the search space into [0, 0.5] and [0.5, 1]. There are no firm rules for the number of subspaces and our suggestion is to first try with two subspaces, then aggregate all such points to obtain . We then sequentially increase the number of subspaces and stop the process when splitting the design space into more subspace does not enlarge the size of .

To find the probability measure that meets the conditions in the equivalence theorem, we proceed iteratively as follows. For each generated design from the algorithm, we first determine elements in and find a candidate for among probability measures supported at the k points. There are 2k variables in the pairs , where is the weight of at . Since the support points of must be roots of the sensitivity function when is optimal, one may apply CSO to minimize the following function with respect to the variable weights : We note that in (19), for each design , we have to first find x that maximizes . This can be done by calling CSO again to tackle this optimization problem. However, after a lot of experiments, we find that only minimizing the second term in (19) frequently suffices to determine , i.e. we suggest given a design, first find by simply minimizingbefore solving the more complicated problem (19) to determine . We ran CSO with , 64 particles for 1200 iterations. We stop the algorithm if and when the objective value attains a small user pre-specified value of, say, .

The CSO-generated designs for the next few examples were found by optimizing (20) and we confirm their G-optimality via the equivalence theorem. On average, the CPU time required to solve (20) was about 5 seconds or less for our examples and we expect a longer time is required to solve (19). This suggests that if we want to find G-optimal designs for more complicated problems, an efficient strategy to find G-optimal designs is to first solve the simpler optimization problem in (20) before solving (19).

The next few tables list CSO-generated designs for various linear mixed models with different assumptions on the covariance structure of the random effects and their G-optimality criterion values. We display the sensitivity plot of each CSO-generated design across the design space and it confirms the G-optimality of the CSO-generated design. All parameters for the design problems, such as the design space, the elements in the covariance matrix , n and m are chosen randomly for illustrative purposes. The error variance is a nuisance parameter and does not affect the optimization process, so we set in all examples.

Sensitivity function of the CSO-generated for model (17) with uncorrelated random effects in Table 1

Table 1

CSO-generated design and its features for model (17) for specific choices of the covariance matrix of the uncorrelated random effects, design space and (n, m)

Model	y_ij=β_i0+β_i1x_j+β_i2xj2+ϵ_ij
D	diag(0.2, 0.2, 0.3)
Design space	[0, 2]
(n, m)	(10, 5)
G-optimal design	0.0000.9662.0000.1460.1400.714
G-criterion value	13.480
μ^∗ on A(ξ)	0.0000.9482.0000.2310.1750.595
Plot of ϕ_G(x,ξ)	Figure 1

Sensitivity function of the CSO-generated for model (17) with correlated random effects in Table 2

Table 2

CSO-generated design and its features for model (17) for specific choices of the covariance matrix of the correlated random effects, design space and (n, m)

Model	y_ij=β_i0+β_i1x_j+β_i2xj2+ϵ_ij
D	0.800.300.100.300.500.080.100.080.40
Design space	[0, 3]
(n, m)	(11, 4)
G-optimal design	0.0001.2743.0000.1650.2700.565
G-criterion value	19.000
μ^∗ on A(ξ)	0.0001.2433.0000.2070.3380.455
Plot of ϕ_G(x,ξ)	Figure 2

Sensitivity function of the CSO-generated design for model (18) with uncorrelated random effects in Table 3

Table 3

CSO-generated design and its features for model (18) for specific choices of the covariance matrix of the uncorrelated random effects, design space and (n, m)

Model	y_ij=β_i0+β_i1xj1/2+β_i2x_j+β_i3xj2+ϵ_ij
D	diag(0.3, 0.5, 0.8, 0.2)
Design space	[1, 3]
(n, m)	(10, 5)
G-optimal design	1.0001.4402.3423.0000.1860.1900.1200.504
G-criterion value	17.939
μ^∗ on A(ξ)	1.0001.4312.3383.0000.2330.2240.1000.443
Plot of ϕ_G(x,ξ)	Figure 3

Sensitivity function of the CSO-generated design for model (18) with correlated random effects in Table 4

Table 4

CSO-generated design and its features for model (18) for specific choices of the covariance matrix of the correlated random effects, design space and (n, m)

Model	y_ij=β_i0+β_i1xj1/2+β_i2x_j+β_i3xj2+ϵ_ij
D	0.800.300.100.050.300.500.080.040.100.080.400.020.050.040.020.30
Design space	[1, 3]
(n, m)	(8, 4)
G-optimal design	1.0001.4192.3723.0000.2300.1860.1220.462
G-criterion value	15.546
μ^∗ on A(ξ)	1.0001.4102.3603.0000.2740.1980.1060.422
Plot of ϕ_G(x,ξ)	Figure 4

CSO-generated locally D-optimal designs for Poisson models with mixed factors

This subsection demonstrates the utility and flexibility of CSO to find other types of optimal designs, such as locally D-optimal designs for estimating parameters in a Poisson regression model with possibly interacting factors and some factors are random.

Poisson models are commonly used to study count data in a regression setting even though they have restrictive assumptions, such as, requiring the mean and variance of the response to be equal. Negative binomial regression models extend Poisson models when the variance is larger or smaller than the mean response. These models are frequently used in clinical trials when the primary outcome is a count variable, such as the number of falls by an elderly patient, or number of CT scans received in the last three months or modeling the numbers of insurance claims and results. Walters (2007a, 2007b) used a negative binomial regression model to account for the number of aggressive incident reports in the following 12 months after subjects were put in an institutional correction center. Like all meta-heuristic algorithms, they can be modified to search for an optimum more effectively for specific problems. An enhanced version of CSO algorithm is available to find various types of optimal designs for the negative binomial regression models (Zhang et al. 2020).

Locally D-optimal designs for two-factor mixed Poisson models with an intercept term have been reported for two models, one with and the other without an interaction term. By using the quasi-information matrix defined in Niaparast and Schwabe (2013) and assuming that the random effects are uncorrelated, i.e., the covariance matrix is diagonal, Naderi et al. (2018) provided theoretical details for finding minimally supported D-optimal designs, including an equivalence theorem to confirm the D-optimality of an approximate design. We recall that minimally supported designs have the number of support points equal to the number of parameters in the mean function and so cannot be used to perform a lack of fit test to check model adequacy. Thus their approach can be restrictive in practice and it is also not clear whether their numerical procedure works for finding D-optimal designs when the model has more interacting factors with random effects or under another design criterion.

We applied CSO and found the same locally D-optimal designs in Tables 1, 2 and 3 of Naderi et al. (2018) for different models with various numbers of uncorrelated random coefficients. The next two sets of results in Tables 5 and 6 show that CSO-generated designs for a model with two interacting factors when the random effects are uncorrelated or correlated. The accompanying plots in Figs. 5 and 6 display their sensitivity plots and they confirm the CSO-generated designs are locally D-optimal. Interestingly, the results show different design spaces can produce D-optimal designs that are not minimally supported.The usefulness of the CSO algorithm can also be seen when we applied it find D-optimal designs for mixed models when some random coefficients are correlated and able to find D-optimal designs that may not be minimally supported. We demonstrate using two Poisson models,one with four and the other with five mixed factors and both have some interaction terms. This is helpful because the bulk of theoretical optimal designs in the literature for nonlinear models have only a couple of additive factors and numerical results for models with interaction terms are also very limited.

Table 5

A locally D-optimal design for a mixed Poisson model with an interaction term and uncorrelated random effects

Model	Y_ij∼P(λ_ij),λ_ij=exp(β_0i+β_1ix_ij,1+β_2ix_ij,2+β_03x_ij,1x_ij,2)
Covariance matrix D	diag(0.3, 0.2, 0.5, 0.6)
Coefficient mean β	[-0.5,0.2,-0.3,0.4]
Design space	[-0.5,1.7]×[-1.0,0.6]
D-optimal design	-0.500-0.5001.2181.7001.700-1.0000.600-1.000-0.5780.6000.2630.3820.0520.1440.159
log|M(ξ)|	− 2.660

Table 6

A locally D-optimal design for a mixed Poisson model with an interaction term and correlated random effects

Model	Y_ij∼P(λ_ij),λ_ij=exp(β_0i+β_1ix_ij,1+β_2ix_ij,2+β_03x_ij,1x_ij,2)
Covariance matrix D	0.300.020.100.000.021.100.600.000.100.601.200.000.000.000.000.00
Coefficient mean β	[-0.9,1.0,1.2,-1.5]
Design space	[-0.8,1.4]×[-1.3,0.5]
D-optimal design	-0.8000.5280.7011.1791.4001.4000.500-1.3000.500-0.280-1.3000.5000.3650.3200.1350.0950.0750.010
log|M(ξ)|	− 4.211

Fig. 5

The sensitivity function of the CSO-generated design under the D-optimality criterion for the two-factor Poisson model with an interaction term and uncorrelated random effects shown in Table 5

Fig. 6

The sensitivity function of the CSO-generated design under the D-optimality criterion for a two-factor Poisson model with an interaction term and correlated random effects shown in Table 6

For the model with four factors, we assume thatthe covariance matrix is and the coefficients in the linear predictor function is . The 12-point CSO-generated design is shown below and criterion value is . Here and elsewhere, the last row in the design for a multi-factor model shows the mass of the design point above it.

The second example concerns a Poisson regression model with five factors () and three interaction terms () defined on . As an illustrative example, suppose the covariance matrix isand the nominal coefficients in the linear predictor function . The criterion value of the 12-point CSO-generated is and the design is shown below. This is a more complicated model and as expected, CSO, like other meta-heuristic algorithms, can also encounter problems when the optimization problem becomes more complex. Convergence and numerical stability issues can arise and repeated reruns of the algorithm with different tuning parameters and swarm size did not produce an optimal design, which is the case here.

A common strategy to try to overcome the above problem is to hybridize the algorithm with another algorithm, such that the hybridized version performs better than either of the algorithms. Another option is to use the equivalence theorem to derive an efficiency lower bound for the generated approximate design without knowing the optimum (Pazman 1986). For D-efficiency, the lower bound depends on the maximum value of the sensitivity function of the design. Since CSO is a stochastic algorithm and depends on the initial design and choice of tuning parameters, we ran the algorithm multiple times and reported the minimum of the D-efficiency of the generated design. For this example, the CSO generated design is and its computed minimum D-efficiency is , which may suffice in practice.

For models with three or more factors, like in the last two examples, the CSO-generated designs appear to be D-optimal or highly D-efficient. Unlike previous examples, the sensitivity functions are now high dimensional and so it is difficult to plot them and visually appreciate their features. One can discretize the design space using a fine grid evaluate the values of their sensitivity function at every grid point but this can be time-consuming. Another way is to apply CSO to optimize the multi-dimensional sensitivity function and confirm all peaks occur at the design points of the CSO-generated design. We chose the latter option and have verified that the CSO-generated design is numerically optimal and is not.

Are D and G-optimal designs equivalent for hierarchical linear models?

It is well known that for regression models with homoscedastic errors, the D and G-optimal approximate designs are equivalent (Kiefer and Wolfowitz 1960). The two optimality criteria look different and have different purposes and so it is an intriguing result. Does the equivalence of the two types of designs apply when we have linear mixed regression models?

Let be the n-dimensional square matrix whose elements are all ones, let , and let denote the Kronecker product. D-optimal design for model (2) minimizes the log whereTables 7 and 8 display CSO-generated designs for the D and G-optimality criteria on different settings for the mixed quadratic model with random components having various covariance matrices. Tables 9 and 10 show corresponding results for a fractional polynomial model with random components having various covariance matrices. The numerical results suggest that D and G-optimal designs for these models are almost equivalent since their D and G-efficiencies relative to the other are all very close to 1. Additional numerical results not shown here for space consideration, for other models we have investigated with more polynomial or fractional polynomial terms and different types of covariance structures also show that the two types of optimal designs have relative efficiencies very close to 1.

Table 7

D and G-optimal designs for a quadratic mixed model with uncorrelated random effects

Model	y_ij=β_i0+β_i1x_j+β_i2xj2+ϵ_ij
D	diag(0.2, 0.2, 0.3)
Design space	[0, 2]
(n, m)	(10, 5)
G-optimal design	0.0000.9662.0000.1460.1400.714
G-criterion value	13.480
D-efficiency	99%
D-optimal design	0.0000.9462.0000.1570.1400.703
D-criterion value	− 4.666
G-efficiency	99%

Table 8

D and G-optimal designs for a quadratic mixed model with correlated random effects

Model	y_ij=β_i0+β_i1x_j+β_i2xj2+ϵ_ij
D	0.800.300.100.300.500.080.100.080.40
Design space	[0, 3]
(n, m)	(11, 4)
G-optimal design	0.0001.2743.0000.1650.2700.565
G-criterion value	19.000
D-efficiency	99%
D-optimal design	0.0001.2393.0000.1630.2720.565
D-criterion value	− 8.798
G-efficiency	99%

Table 9

D and G-optimal designs for a fractional polynomial mixed model with uncorrelated random effects

Model	y_ij=β_i0+β_i1xj1/2+β_i2x_j+β_i3xj2+ϵ_ij
D	diag(0.3, 0.5, 0.8, 0.2)
Design space	[1, 3]
(n, m)	(10, 5)
G-optimal design	1.0001.4402.3423.0000.1860.1900.1200.504
G-criterion value	17.939
D-efficiency	99%
D-optimal design	1.0001.4122.3613.0000.2300.1860.1220.462
D-criterion value	− 2.116
G-efficiency	99%

Table 10

D and G-optimal designs for a fractional polynomial mixed model with correlated random effects

Model	y_ij=β_i0+β_i1xj1/2+β_i2x_j+β_i3xj2+ϵ_ij
D	0.8000.3000.1000.0500.3000.5000.0800.0400.1000.0800.4000.0200.0500.0400.0200.300
Design space	[1, 4]
(n, m)	(10, 5)
G-optimal design	1.0001.6022.9144.0000.2060.2180.1090.467
G-criterion value	20.271
D-efficiency	99%
D-optimal design	1.0001.5992.8844.0000.2050.2170.1100.468
D-criterion value	− 89.985
G-efficiency	99%

Our numerical results in the tables suggest that the celebrated theorem of Kiefer and Wolfowitz may also apply to hierarchical linear models, but a general proof is elusive at this time. However, we are able to show the equivalence of the two types of designs holds for hierarchical linear models when they contain only a random intercept. To see this, we recall D-optimality discussed in Prus and Schwabe (2016b) for the hierarchical linear model (2). Assuming that the dispersion matrix has rank q, the D-criterion for prediction is the logarithm of the product of the largest eigenvalues of the MSE-matrix:where are the q largest eigenvalues of . Corollary 6 of Prus and Schwabe (2016b) used this definition and gave an equivalence theorem for D-optimality, and showed that the D-optimal design in the fixed-effects model is also D-optimal for prediction in the random intercept model.

We next consider using the G-optimal design for prediction in a hierarchical linear model with a random intercept. Assume in model (2). The dispersion matrix can be written as , where denotes the first unit vector in . For an approximate design , the MSE-matrix for prediction in this random-intercepts model iswhere as defined in (13).

It follows that the lower bound in (10) iswhich does not depend on . If is D-optimal for the fixed-effects model, we haveand it follows that the design is G-optimal for prediction of individual parameters. Therefore, the D-optimal and G-optimal designs are equivalent for prediction in the random-intercept model.

Conclusions

G-optimal designs are challenging to determine and study because the criterion is not differentiable and they require solving two or more layers of nested optimization problems over different spaces. However, the criterion is compelling and should appeal to researchers interested to design an experiment to estimate the overall response surface. To facilitate greater use of such designs for linear models with random effects, we proposed an effective meta-heuristic algorithm call CSO to find G-optimal designs and developed an equivalence theorem to confirm whether a design is G-optimal.

Additionally, we showed CSO is flexible and can also search for locally D-optimal designs for Poisson mixed regression models with several interacting factors, where some random effects may be correlated. We also provide R-codes freely and the interested reader can request them from the third author.

d_1	d_2	w1∗	w2∗	w3∗	wμ1∗	wμ2∗	wμ3∗
1	1	0.2788	0.2788	0.4423	0.3272	0.3272	0.3456
	5	0.2312	0.3674	0.4014	0.3268	0.3355	0.3377
	10	0.2251	0.3789	0.3960	0.3269	0.3361	0.3371
5	1	0.3674	0.2312	0.4014	0.3355	0.3268	0.3377
	5	0.3215	0.3215	0.3569	0.3331	0.3331	0.3337
	10	0.3156	0.3333	0.3511	0.3331	0.3333	0.3336
10	1	0.3789	0.2251	0.3960	0.3361	0.3269	0.3371
	5	0.3333	0.3156	0.3511	0.3333	0.3331	0.3336
	10	0.3274	0.3274	0.3452	0.3333	0.3333	0.3334