Longitudinal studies are common in clinical studies and involve multiple measurements of the subject over time. For example, subjects in a clinical trial are recruited and randomized into treatment groups and responses from each patient are observed multiple times over a user-selected period of time. There is a huge literature on analyzing longitudinal models using different methods but design issues for such studies lagged. What seems to be generally known is that having excessive time points does not necessarily improve the quality of the statistical inference and having too few observations per patient to facilitate calculation may not meet the scientific requirements in the study.

Given a statistical model and a design criterion, design issues for a longitudinal study concern the optimal number of time points to observe the outcomes, where the time points are over the study period and the number of replicates at each time point. For making statistical inference as accurately as possible, these decisions have to be made judiciously and at minimum cost. Analytical approaches are difficult because they invariably involve number-theory to solve the optimization problem. Even when it is possible to derive the results mathematically, they can be limiting because they are valid to a single statistical model and does not apply to even a slightly changed model. Software and efficient computational methods are therefore useful for answering such design questions.

We propose using nature-inspired metaheuristic algorithms to tackle these challenging questions. This class of algorithms is commonly used in engineering and computer science to tackle hard to solve optimization problems but surprisingly under-used in mainstream statistical research. We focus on particle swarm optimization (PSO), which is a popular member of this class of algorithms and apply it to design longitudinal studies with different mean functions and various correlation structures. Such algorithms typically have variants, which are enhancements motivated by different desires to improve certain aspects of its performance. We also compare results from PSO with its various enhancements and determine whether the latter provide markedly superior results than the results from the original PSO.

To fix the idea, we discuss the proposed methodology for the 2-parameter nonlinear Michaelis-Menten model which is widely used to study enzyme-substrate dose relationship in kinetic biological systems. Its simplicity and usefulness make it one of the most widely used models across many disciplines, such as, in agriculture [37], biochemistry [19], biology [4], environmental study [36]. [35] notes that the Michaelis-Menten model is a special case of the 3-parameter logistic model, which is more widely used in practice.

In its simplest form, the Michaelis-Menten model for a chemical reaction is E ( υ ) = η ( θ , t ) = a t b + t , t ∈ T , θ = ( a , b ) ⊤ , where υ is the observed velocity of the reaction when the substrate concentration applied is t ∈ T, where S is the user-selected range from which concentration is selected. The mean response is a nonlinear function η ( θ , t ) with 2 parameters a and b. The parameter b is called the Michaelis-Menten constant, which controls the rate of the reaction and so between the two parameters, it is the more biologically meaningful parameter. The maximum velocity attainable is a, which is reached when the substrate concentration is increased without bound. Errors in this model are assumed to be independent with mean 0 and constant variance σ 2 .

This model can also be used to study the growth curves of animals or subjects by taking repeated measurements on the subjects [17]. Suppose that there are n subjects, and for each subject, m repeated observations are allowed over a period of time. [7] addressed design issues when the responses are correlated using the following model (1.1) υ = η ( θ , t i , j ) + ϵ i , j = a t i , j b + t i , j + ϵ i , j , j = 1 , … , m , i = 1 , … , n , where the error terms ϵ i , j are normally distributed, each with mean 0 and constant variance and ϵ i , j and ϵ i ′ , j ′ are independent if i ≠ i ′ , and each pair of ϵ i , j and ϵ i , j ′ has the correlation coefficient λ ∣ t i , j − t i , j ′ ∣ , j , j ′ = 1 , … , m. Here λ ∈ ( 0 , 1 ) and t i , j ′ ∈ S. Given n and an optimality criterion, design issues for the model concern how to optimally select m, the number of time points and the values of the s to observe outcomes so that the criterion is optimized. A typical criterion is to estimate the parameters a and b or some function thereof as accurately as possible.

Throughout, we assume the following setup. We are given a total of N observations for the study and a regression model defined on a given design space S. Given an optimality criterion, finding an optimal exact design means finding the optimal value of the criterion by choosing the optimal number of points k to observe the response, where the time points and the number of replicates n i at each time point t i so that n 1 + n 2 + ⋯ + n k = N. A complicating feature of such an optimization problem is that there is neither a general theory for finding an optimal exact design nor a theoretical tool for confirming the optimality of an exact design. Consequently, optimal exact designs are rarely studied in the literature, especially for nonlinear models. As an early example, [3] found D-optimal exact designs numerically for homoscedastic polynomial models of low degrees.

Our goal is to investigate the effectiveness of using a nature-inspired metaheuristic algorithms, such as particle swarm optimization (PSO), for finding D-optimal exact designs for the Michaelis-Menten model and related models with correlated errors, c-optimal designs for estimating a function of the parameters in a growth curve model, and a maximin optimal exact design for optimizing efficiencies in a HIV dynamic model. We show our results agree with the theoretical designs in [7], which are only available for small values of m and n and able to find the optimum in all problems quite quickly. However, there is no general theory for finding optimal exact designs or confirming optimality of an exact design, which may, in part, explain why exact designs are much less studied than approximate designs, where there is a general framework and general methods for finding and conforming optimality of an approximate design when the criterion is convex functional.

The rest of the article is organized as follows. Section 2 briefly reviews particle swarm optimization. Section 3 applies PSO to search for locally D-optimal exact designs for the Michaelis-Menten model with various correlation structures. The Michaelis-Menten model is a nonlinear model, and, as we explain below, its information matrix and hence the design criterion depends on unknown model parameters that we wish to estimate. The simplest design strategy is to use results from a pilot study or similar studies to have a good guess of the true parameters that we wish to estimate. These nominal values, which may be also available from experts, are then substituted into the design criterion and the problem then involves optimizing the design parameters only. We next discuss a few PSO enhancements in Section 3.2 and compare their performance for finding locally D-optimal designs for the Michaelis-Menten model with various error structures. Additionally, we apply PSO to find locally c-optimal designs for estimating one or more functions for a generalized version of the Michaelis-Menten model for studying growth rates in Section 4. Section 5 considers an alternative to locally optimal design and constructs maximin locally exact designs for a longitudinal study to estimate parameters for a HIV dynamic model. Section 6 concludes with a short discussion.

Recently a class of algorithms called nature-inspired metaheuristic algorithms has proved very popular in the optimization literature. [31, 32] provided reasons for the rapid rise and interest in these algorithms. Early users of such algorithms to find optimal exact designs for linear models include [12] who used an annealing algorithm to search for optimal designs for linear regression models, and [20] who used a genetic algorithm to construct exact D-optimal designs. Of particular note is the particle swarm optimization (PSO) introduced by [8] for tackling optimization problems. PSO is increasingly used across disciplines to solve hard optimization problems. PSO is essentially assumptions free and it searches in a simple and effective way. For example, unlike many algorithms, PSO does not require the objective function to be differentiable or convex and can solve non-convex high-dimensional optimization problems.

PSO is a metaheuristic optimization algorithm inspired from the way animals, such as birds and fishes, search for food. The birds fly continuously in the sky to look for food on the ground. Each has its own perception where the food is (local optimum) but it communicates with the rest and collectively the flock decides where the food is (global optimum). Accordingly, each bird flies toward the global optimum in the direction of the local optimum (not giving up completely where it thinks the food is). Birds are referred as particles and each bird represents a candidate solution to the optimization problem. Velocities and locations of each bird are adjusted at each iteration and if and when the flock converges, the perceived global optimum is found. In order to efficiently identify the optimal points, we initiate a flock of birds in the pre-defined search space. Let X ( k ) be the locations of particles at the k-th iteration. Define t L ( k − 1 ) to be the locations with the best objective function values discovered by each particle before the k-th iteration, and t G ( k − 1 ) to be the locations of the best value found by the whole swarm before the k-th iteration. At the k-th iteration, the particles are updated by (2.1) X ( k ) = X ( k − 1 ) + V ( k ) , and (2.2) V ( k ) = w V ( k − 1 ) + c 1 R 1 ⊗ [ t L ( k − 1 ) − X ( k − 1 ) ] + c 2 R 2 ⊗ [ t G ( k − 1 ) − X ( k − 1 ) ] , where V ( k ) is the velocity of the particle. There are several parameters in (2.2). The inertia weight represents how active the birds are and is denoted by w. This parameter may be chosen to be a positive constant but more typically its value changes over iteration and eventually decrease to 0. The parameters c 1 and c 2 are two positive constants which are recommended to be 2, and R 1 and R 2 are two random vectors whose components are independently drawn from the uniform variate on [ 0 , 1 ]. In practice, the number of iterations and the swarm size are the most influential parameters in PSO. Large swarm size gives a better vision on the search area such that PSO could achieve the global optimum with a higher chance, and the more iterations allows the particles having search experience due to random perturbation. More details on PSO and the related metaheuristic optimization algorithms are available in [34].

For an individual patient with n = 1 and N = m, we determine the optimal N time points to observe responses from the patient. We allow replications, i.e. the N time points need not be distinct. To search for a N-point optimal exact design, we search for an optimal N × 1 vector, ( t 1 , t 2 , … , t N ) and allow some of the t i ’s may be the same. Here N is user-specified and the objective function is the design criterion Φ ( · ), which may not be convex function. A key challenge is the dimension of the problem when N is large. When there are multiple optimal solutions, and we only report one of them. Hence at the beginning of the PSO, we randomly generate initial particles (designs) on the design space. Then at each iteration, we update the particles (designs) based on the (2.1) and (2.2). We also watch out for those that flew outside the search boundaries and when this happens, we need to adjust those particles to make sure they are in the design space properly. The hope is that after a number of iterations, the particles will converge to a point and this point is supposedly the global best solution or the optimal design we are after. Algorithm 1 below is a pseudo code of the PSO algorithm for finding an optimal design.

In the next section and beyond, we apply PSO-based algorithms to find various types of optimal exact designs for models with uncorrelated errors and errors with various correlation structures and compare the designs. Throughout, the hardware we used is a PC with 3.50 GHz Intel(R) Core(TM) i7-4770K CPU.

Suppose in a longitudinal study, we have n subjects and each is observed at m time points t i j , j = 1 , … , m from a time interval T. We assume T, m and n are given. The normalized information matrix is proportional to M ( ξ , θ 0 ) = 1 N ∑ i = 1 n f i ⊤ Σ i − 1 f i , where N = n m, θ 0 ⊤ = ( a , b , σ 2 , λ ) and f i is the m × 2 matrix with the j-th row containing the derivative of f i , which is t i , j b + t i , j − a t i , j b + t i , j 2 , j = 1 , … , m .

The m × m correlation matrix of the m responses from the ith subject is Σ i and its ( j , k )-th element is σ 2 λ ∣ t i , j − t i , k ∣ , j , k = 1 , … , m . Thus, given a fixed parameters θ 0 = ( a , b , σ 2 , λ ) ⊤ , the locally D-optimal design for estimating the model parameters as accurately as possible is the one that maximizes Φ ( ξ , θ 0 ) = log | M ( ξ , θ 0 ) | over all designs ξ on T. Following the argument in [7], we may, without loss of generality, set T = [ 0 , 1 ]. For the Michaelis-Menten model, locally D-optimal designs do not depend on a and σ 2 and so, we choose a = σ 2 = 1. We also assume that the experimental conditions for each subject are the same and each subject is observed at the same set of time points, t 1 j = t 2 j = ⋯ = t n j , for j = 1 , … , m. Thus, to search for the locally D-optimal exact design via PSO, we can set n = 1, i.e. N = 1 × m = m, and the target design is represented as an m × 1 vector.

We first consider the case when we have 2 observations from a single subject and the design questions are which two time points to take observations and how the spread out the observations between the two points. [7] theoretically identified that the locally D-optimal exact design is supported at two points, u and 1, when b ≥ 1 3 , and the design point u solves the equation (3.1) b − ( 2 b + 1 ) u u ( 1 − u ) ( b + u ) = log ( λ ) λ 2 ( 1 − u ) 1 − λ 2 ( 1 − u ) . The solution has no closed form when b < 1 / 3. When there are 3 or 4 observations to be taken from each subject, locally D-optimal exact designs cannot be described analytically and only numerical results were provided. The numerical approach they used is to check all possible designs using a pre-specified grid.

We next use Algorithm 1 to systematically find the locally D-optimal exact designs. We used 256 particles and set the maximum number of iterations equals to 500. The PSO parameters c 1 and c 2 were set to their default values with both equal to 2. For the inertia weight, w, we started from 0.95 and linearly decreased it to 0.4 in the first 350 iterations and then fixed w as 0.4 for the remaining iterations. We implemented the algorithm and found PSO-generated D-optimal designs for several values of b and λ, i.e. b = 0.2 , 0.7 , 1.2 , 1.7 , 2.2 , 2.7 and λ = 0.1 , 0.5 , 0.9. We first consider the cases when the numbers of observations per subject were m = 2 , 3 , 4 and a direct application of PSO shows the generated designs coincided with the analytical designs in [7] when m = 2 and b ≥ 1 / 3. The PSO-generated designs when we have 3 or 4 observations per subject also agree with those numerically found optimal designs from [7]. For space consideration, we do not display them. We next use Algorithm 1 to search for D-optimal designs when there are more observations per subject. Table 1 displays the PSO-generated optimal exact designs for 5 and 6 observations per subject. Unlike the global numerical search approach in [7], we are not constrained by the grid size imposed on the experimental region, S, and PSO can still identify the best designs efficiently. The CPU time required by our Algorithm 1 to find the optimal design is short. On average, our MATLAB codes take around 4.11, 5.54, 6.72, 7.83 and 8.71 seconds to find the optimal designs when there are 2, 3, 4, 5 and 6 observations, respectively.

In summary, Algorithm 1 seems efficient for finding locally D-optimal exact designs for the Micahelis-Menten model. Tables 2 to 4 show PSO-generated designs have more design points as the value of m is increased. To conserve space, we consider particular cases of ( λ , b ) ⊤ when λ is fixed at 0.1, 0.5 and 0.9, and b = 0.5 , 1.0 , 1.5 , 2.0 and 2.5.

A growth curve describes the change of an outcome over time. To estimate a growth curve model, we take multiple observations at different time points from the same subject and so the observed responses are correlated.

A generalized version of the Michaelis-Menten model to study the biomass W of an animal over time was introduced by [17]. At time t, the model is as follows: (4.1) W ( t ; θ ) = W 0 K h + W f t h K h + t h , t > 0 , θ = ( W 0 , W f , K , h ) ⊤ where W 0 and W f are the zero- and infinite-time values of the biomass, respectively. If K > 0, it is the time when half-maximal growth is achieved. When h = 1 and W 0 = 0, model (4.1) reduces to the usual Michaelis-Menten model.

We consider the data set from [17] consisting of 17 weight (kg) records of one castrated male Percheron horse from birth to 220 weeks of age. Its growth pattern has been described by the generalized Michaelis-Menten model with nominal values θ 0 = ( W 0 , W f , K , h ) ⊤ = ( 85.50 , 731.00 , 56.70 , 1.39 ) ⊤ in [17]. Using θ 0 , they found a locally D-optimal exact design, ξ D ∗ for estimating parameters in (4.1) by maximizing the determinant of the Fisher information matrix. For model (4.1), the elements in the covariance matrix are Σ i j = σ 2 exp { ( − λ | t i − t j | } and f i = k h k h + t i h , t i h k h + t i h , h k h − 1 W 0 k h + t i h − h k h − 1 ( k h W 0 + t i h W f ) ( k h + t i h ) 2 , ( k h W 0 log k + t i h W f log t i ) k h + t i h − ( k h log k + t i h log t i ) ( k h W 0 + t i h W f ) ( k h + t i h ) 2 ⊤ .

Human immunodeficiency viruses (HIV) are a subgroup of retrovirus that can cause infection, resulting in progressive failure of the immune system. The presence of viruses is usually detected when the CD4 cell counts are low. Monitoring the cell counts continuously is key to understanding disease progression. Various statistical models have been proposed to model CD4 cell count longitudinally. A common model is to assume a fixed effect model defined on a time line in hours over a pre-specified time interval, say, t ∈ T = [ 0 , 6 11 12 ] = [ 0 , 6.917 ], (5.1) E ( Y j ∣ V 0 , c , δ , t j ) = η ( V 0 , c , δ , t j ) = log V 0 + log c 2 ( c − δ ) 2 e − δ t j − c 2 − ( c − δ ) 2 ( c − δ ) 2 e − c t j − c δ c − δ t e − c t j .

The above time interval is taken from Table 1 in [13] that describes HIV RNA copies per milliliter plasma and time is measured since treatment initiation and the predictor, t, is time after pharmacologic delay. Frequently interest is in estimating one or more of the parameters in the model parameters since not every parameter has a meaningful biological interpretation. When all parameters are interesting, the D-criterion discussed previously is an appropriate design objective; otherwise, if only a subset of the parameters is interesting, the design criterion is c-optimality. In the above model, we are interested to find an optimal design to estimate θ T = ( log V 0 , log c , log δ ) and optimal designs to best estimate each of its components. Additionally, we also find a maximin optimal design goal that maximizes the lowest efficiency among the three efficiencies for estimating each of the three parameters.

For such studies, there is usually a standard protocol recommended by physicians for implementation. The protocol specifies how many times and when to sample blood to determine cd4 counts and for other laboratory tests. However, the protocol is usually not based on statistical considerations. We now apply PSO to find various optimal designs and ascertain how efficient is the protocol design under various scenarios, including situations where there are multiple objectives and there are different emphases in each of the study objectives.

We compute various optimal exact designs and evaluate the efficiencies of the recommended or protocol design or baseline design for such a study. For a given number of points, uniform designs have them evenly distributed on the time interval T. For example, such a design with 8 observations is ξ 8 = { 0.000 , 0.917 , 1.917 , 2.917 , 3.917 , 4.917 , 5.917 , 6.917 }. A direct calculation shows its Fisher information matrix is proportional to (5.2) M ( ξ , θ ) = 1 N F ⊤ F where the ( i j ) t h element of the matrix F is F i j ( t ) = ∂ η ( θ i , t ) / ∂ θ j and η ( θ , t ) is the mean response at time t ∈ T. For estimating all parameters in the model, D-optimality is commonly used and is defined by the determinant of M ( ξ , θ ). Because the information matrix contains unknown parameters, we replace θ by its nominal values before we maximize the determinant using PSO. Using nominal values θ ⊤ = ( log V 0 , log c , log δ ) ⊤ = ( 11.0 , 1.1 , − 1.0 ), the resulting PSO-generated locally D-optimal design with 8 observations is ξ D = 0.000 , 0.000 , 0.000 , 2.083 , 2.083 , 6.917 , 6.917 , 6.917 . This implies that the implemented design requires 3 replicates at 0, 2 replicates at 2.08 and 3 replicates at 6.92.

To estimate log c or log δ, we find an exact design that minimizes the asymptotic variance of the estimated function. Setting c ⊤ = ( 0 , 1 , 0 ) to estimate log c, c ⊤ = ( 0 , 0 , 1 ) to estimate log δ in (4.2) and using the given nominal values for the model parameters, the PSO-generated designs for estimating log c- and log δ are ξ log c = { 0.000 , 0.000 , 0.000 , 2.113 , 2.113 , 2.113 , 2.113 , 6.917 } and ξ log δ = { 0.000 , 1.923 , 1.923 , 1.923 , 1.923 , 6.917 , 6.917 , 6.917 } .