D-optimal population designs in linear mixed effects models for multiple longitudinal data

Abstract

The main purpose of this paper is to investigate D-optimal population designs in multi-response linear mixed models for longitudinal data. Observations of each response variable within subjects are assumed to have a first-order autoregressive structure, possibly with observation error. The equivalence theorems are provided to characterise the D-optimal population designs for the estimation of fixed effects in the model. The semi-Bayesian D-optimal design which is robust against the serial correlation coefficient is also considered. Simulation studies show that the correlation between multi-response variables has tiny effects on the optimal design, while the experimental costs are important factors in the optimal designs.

Keywords:

• D-optimal designs
• longitudinal data
• multi-response linear mixed model
• equivalence theorem

MSC 2010:

• 62K05

1. Introduction

In many statistical studies, it is quite common that repeated measured data are collected over time and are often referred to as longitudinal data. Mixed effects models are widely used for the analysis of such correlated data. However, the responses are usually time-structured when they are measured at different time points. Modeling of autocorrelations through use of autoregressive error structure, in addition with random effects terms, may be more appropriate (Chi & Reinsel, 1989). For a list of correlation structures, we also refer to the monographs by Verbeke and Molenberghs (2000) and Molenberghs and Verbeke (2005).

The aim of optimal experimental design is to arrange a reasonable design schedule, so that the experimenter can obtain the maximal information of the experiment at minimal cost. In the case of single-response with correlated errors under the AR(1) structure assumption, Tan and Berger (1999) studied optimal designs for linear mixed effects models with longitudinal data, while Berger et al. (2003) and Azurduy (2009) studied the robust designs in a maximin approach. Nyberg et al. (2012) considered the optimal designs for nonlinear mixed effects models with serial correlation. Liu et al. (2019) discussed the D-optimal design for the multi-response linear mixed effects models by modelling the serial correlation with totally mixed effects. For the case of longitudinal data generated from the multi-response models with repeated measurements, correlated errors also play an important role in model description and estimation. Roy (2006) discussed how the correlation coefficient between two response variables was effected by the autocorrelation structure of repeated measures. Wang and Fan (2010) applied the multi-response linear mixed effects model (MLMM) with correlated errors to model repeated continuous data and presented ECM procedure to obtain the maximum likelihood estimates of model parameters. A general description of the statistical analysis of the multi-response longitudinal models can be found in Verbeke et al. (2014).

In this paper, we are interested in deducing the approximate designs of the fixed effects parameters in the multi-response linear mixed effects models with AR(1) autoregressive error structure. The generalised estimating equation (GEE) approach (Liang & Zeger,1986) is used to estimate the model. The effects of the experimental cost on the optimal population design is also discussed.

The rest of the paper is organised as follows. Section 2 introduces the notation and provides the Fisher information matrices of individual design and population design. Section 3 gives the equivalence theorem for the D-optimal population design with respect to the estimation of the fixed effects and defines the relative D-efficiency by considering a linear experimental cost. The methodology is illustrated with a real example concerning the AIDS Clinical Trials Group Study in Section 4. Concluding remarks are given in Section 5.

2. Model description and information matrix

2.1. Assumptions and preliminaries

Suppose there are N independent individuals in a longitudinal study and for each individual there are r response variables observed over time. Although each response variable may be observed under different design schedules, index matrices can be used (Ogungbenro et al., 2007) to describe it. Hence, a balanced design (Schmelter, 2007b), in which all r responses are measured under the same design schedule, is considered.

Let ${\mathit{y}}_{ij}=(y_{ij1},y_{ij2},\dots ,y_{ijr}{)}^T$ be an $r\times 1$ response vector observed at the jth measurement of individual i, and let ${\mathit{e}}_{ij}=(e_{ij1},e_{ij2},\dots ,e_{ijr}{)}^T$ be the within individual observation error vector with zero mean, $j=1,2,\dots ,m_i$. Let $\mathit{G}$ be an $r\times r$ covariance matrix describing the correlation of the r different responses. The multi-response linear regression model with mixed effects is defined as (1) ${\mathit{y}}_{ij}={\mathit{X}}_{ij}\mathit{\beta }+{\mathit{Z}}_{ij}{\mathit{b}}_i+{\mathit{e}}_{ij},$(1) where ${\mathit{X}}_{ij}$ is an $r\times p$ fixed-effects design matrix, ${\mathit{Z}}_{ij}$ is an $r\times l$ random-effects design matrix, $\mathit{\beta }$ is a p-dimensional fixed-effects vector and ${\mathit{b}}_i$ is an l-dimensional ($l\le p$) random-effects vector with mean zero and $l\times l$ variance–covariance matrix $\mathit{D}$ (Molenberghs & Verbeke, 2005). It is assumed that $\mathit{D}$ is a positive definite matrix.

In general, in a longitudinal study, at least three repeated observations of the responses are required over time. We assume ${\mathit{t}}_i=(t_{i1},t_{i2},\dots ,t_{i{m}_i})$ ($m_i\ge 3$) is an $m_i$-point design (without repetition) for individual i with equal weight $1/m_i$ on each point, $i=1,2,\dots ,N$. The observations of different subjects are independent, but those within subject i conditioned on ${\mathit{b}}_i$ are dependent. Let ${\mathit{R}}_i$ be a working autocorrelation matrix among all $m_i$ possible observations $y_{i1k},y_{i2k},\dots ,y_{i{m}_{i}k}$ of the kth response variable. The common working correlation structure could be as simple as independent, equicorrelated, first-order autoregressive AR(1), or could be unspecified. In this paper, the autoregressive process of order one, AR(1) is used (Diggle, 1988; Wang & Fan, 2010), and then the autocorrelation matrix of the individual i is ${\mathit{R}}_i=({\rho }^{|t_{ij}-t_{i{j}^{\prime }}|})$, where ρ ($0\le \rho \le 1$) is the autocorrelation coefficient that gives the correlation between two time points that are one unit of time apart. This correlation structure implies that repeated measurements are more highly correlated in closer time and become less correlated as the difference in time points increases.

For each individual i, denoting ${\mathit{y}}_i=({\mathit{y}}_{i1}^T,{\mathit{y}}_{i2}^T,\dots ,{\mathit{y}}_{i{m}_i}^T{)}^T$ and ${\mathit{e}}_i=({\mathit{e}}_{i1}^T,{\mathit{e}}_{i2}^T$, …, ${\mathit{e}}_{i{m}_i}^T{)}^T$, the multi-response linear mixed model in (1(1) ${\mathit{y}}_{ij}={\mathit{X}}_{ij}\mathit{\beta }+{\mathit{Z}}_{ij}{\mathit{b}}_i+{\mathit{e}}_{ij},$(1) ) is expressed as (2) ${\mathit{y}}_i={\mathit{X}}_i\mathit{\beta }+{\mathit{Z}}_i{\mathit{b}}_i+{\mathit{e}}_i,$(2) where ${\mathit{X}}_i=({\mathit{X}}_{i1}^T,{\mathit{X}}_{i2}^T,\dots ,{\mathit{X}}_{i{m}_i}^T{)}^T$ and ${\mathit{Z}}_i=({\mathit{Z}}_{i1}^T,{\mathit{Z}}_{i2}^T,\dots ,{\mathit{Z}}_{i{m}_i}^T{)}^T$. It is assumed that the random-effects vector ${\mathit{b}}_i$ is independent with ${\mathit{e}}_i$. By letting ${\mathit{A}}_i$ be the diagonal marginal variance matrix of ${\mathit{y}}_i$, the covariance matrix ${\mathit{H}}_i$ of ${\mathit{y}}_i$ in model (2(2) ${\mathit{y}}_i={\mathit{X}}_i\mathit{\beta }+{\mathit{Z}}_i{\mathit{b}}_i+{\mathit{e}}_i,$(2) ) can be expressed as (Wang & Fan, 2010) (3) ${\mathit{H}}_i=\mathrm{Cov}({\mathit{y}}_i)={\mathit{Z}}_i\mathit{D}{\mathit{Z}}_i^T+{\mathit{A}}_i^{1/2}({\mathit{R}}_i\otimes \mathit{G}){\mathit{A}}_i^{1/2}.$(3) According to the GEE method, the Fisher information matrix of a single individual design ${\mathit{t}}_i=(t_{i1},t_{i2},\dots ,t_{i{m}_i})$ for individual i is (4) ${\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_i)={\mathit{X}}_i^T{\mathit{H}}_i^{-1}{\mathit{X}}_i.$(4) The Fisher information matrix in (4(4) ${\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_i)={\mathit{X}}_i^T{\mathit{H}}_i^{-1}{\mathit{X}}_i.$(4) ) depends on the value of ρ, which is often not explicitly known in an experiment. Assuming that ρ has a range of possible values, we investigate a robust D-optimal design criterion (see Section 3).

2.2. Information matrix of the population design

Consider now a population design comprising Q distinct individual designs with $n_q$ individuals allocated to the design with $m_q$ time points ${\mathit{t}}_q=(t_{q1},t_{q2},\dots ,t_{q{m}_q})$ for $q=1,2,\dots ,Q$ (Schmelter, 2007a, 2007b). Denote this population design as follows: (5) $\zeta =\left\{\begin{array}{cccc}{\mathit{t}}_1& {\mathit{t}}_2& \cdots & {\mathit{t}}_Q\\ w_1& w_2& \cdots & w_Q\end{array}\right\},$(5) where $w_q=n_q/N$. To turn the discrete optimisation problem into a continuous one, we allow approximate designs, where the weights $w_q$'s can be arbitrary numbers in the interval $(0,1)$ summing up to 1. Then, the Fisher information matrix for the fixed-effects vector $\mathit{\beta }$ at the population design ζ is defined by (6) $\mathit{M}(\zeta )=\sum _{q=1}^Q{w}_q{\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_q),$(6) where ${\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_q)$ is determined as in (4(4) ${\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_i)={\mathit{X}}_i^T{\mathit{H}}_i^{-1}{\mathit{X}}_i.$(4) ). By defining $\mathit{X}=({\mathit{X}}_1^T,{\mathit{X}}_2^T,\dots ,{\mathit{X}}_Q^T{)}^T$, $\mathit{H}=\mathrm{diag}({\mathit{H}}_1,{\mathit{H}}_2,\dots ,{\mathit{H}}_Q)$ and $\mathit{W}=\mathrm{diag}(w_1,w_2,\dots ,w_Q)$, the matrix $\mathit{M}(\zeta )$ in (6(6) $\mathit{M}(\zeta )=\sum _{q=1}^Q{w}_q{\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_q),$(6) ) can be expressed as (7) $\mathit{M}(\zeta )={\mathit{X}}^T{\mathit{H}}^{-1/2}\mathit{W}{\mathit{H}}^{-1/2}\mathit{X}.$(7)

3. The D-optimal design and robust design

3.1. The general equivalence theorem

A D-optimal population design, ${\zeta }^{\ast }$, for estimating the fixed effects in the MLMM (1(1) ${\mathit{y}}_{ij}={\mathit{X}}_{ij}\mathit{\beta }+{\mathit{Z}}_{ij}{\mathit{b}}_i+{\mathit{e}}_{ij},$(1) ) maximises the logarithm of the determinant of the Fisher information matrix $\mathit{M}(\zeta )$ given in (7(7) $\mathit{M}(\zeta )={\mathit{X}}^T{\mathit{H}}^{-1/2}\mathit{W}{\mathit{H}}^{-1/2}\mathit{X}.$(7) ), i.e., (8) ${\zeta }^{\ast }=\arg \underset{\zeta }{max}\log |\mathit{M}(\zeta )|.$(8) The General Equivalence Theorem relating to the approximate D-optimal population design is given in the following, which is based on very mild assumptions that the set of all population designs ζ is compact on the design space and the D-optimality criterion $\mathrm{\Psi }(\zeta )=-\log |\mathit{M}(\zeta )|$ is convex and differentiable. The Frechet derivative of $\mathrm{\Psi }(\zeta )$ in the direction of any individual design $\mathit{t}$ is $\varphi (\mathit{t},\zeta )=p-\mathrm{tr}({\mathit{M}}^{-1}(\zeta ){\mathit{M}}_{\mathrm{ind}}(\mathit{t}))$.

Theorem 3.1

A population design ${\zeta }^{\ast }$ is D-optimal for the estimation of the fixed-effects vector $\mathit{\beta }$ if and only if

1. The design ${\zeta }^{\ast }$ maximises $\log |\mathit{M}(\zeta )|.$

2. The design ${\zeta }^{\ast }$ minimise $\underset{\mathit{t}}{max}\mathrm{tr}({\mathit{M}}^{-1}(\zeta ){\mathit{M}}_{\mathrm{ind}}(\mathit{t})).$

3. The maximum of $\mathrm{tr}({\mathit{M}}^{-1}({\zeta }^{\ast }){\mathit{M}}_{\mathrm{ind}}(\mathit{t}))$ is attained at the support individual designs of ${\zeta }^{\ast },$ and the maximum value is p.

This theorem can be used to construct or check the optimality of a population design (5(5) $\zeta =\left\{\begin{array}{cccc}{\mathit{t}}_1& {\mathit{t}}_2& \cdots & {\mathit{t}}_Q\\ w_1& w_2& \cdots & w_Q\end{array}\right\},$(5) ), which provides a basis for the following discussion of optimal designs considering the experimental spending.

3.2. Semi-Bayesian D-optimal design

Both the individual Fisher information matrix in (4(4) ${\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_i)={\mathit{X}}_i^T{\mathit{H}}_i^{-1}{\mathit{X}}_i.$(4) ) and the population Fisher information matrix in (7(7) $\mathit{M}(\zeta )={\mathit{X}}^T{\mathit{H}}^{-1/2}\mathit{W}{\mathit{H}}^{-1/2}\mathit{X}.$(7) ) depend on ρ, which is always unknown before starting an experiment. Assume that ρ obeys a distribution with probability density function $f(\rho )$ defined on an interval $[a,b]$ ($0\le a<b\le 1$). A design is said to be semi-Bayesian D-optimal, denoted by ${\zeta }_B^{\ast }$, if it maximises (9) ${\mathrm{\Phi }}_B(\zeta )={\int }_a^b\log |\mathit{M}(\zeta )|f(\rho )\mathrm{d}\rho .$(9) It is known that the semi-Bayesian D-optimal design is robust against a variety of ρ.

In order to calculate the integral defined by (9(9) ${\mathrm{\Phi }}_B(\zeta )={\int }_a^b\log |\mathit{M}(\zeta )|f(\rho )\mathrm{d}\rho .$(9) ), a set of 901 points is drawn from the prior distribution of ρ. Numerical approximations are applied by sampling parameter values from the prior distribution of $f(\rho )$ and replacing the integral in (9(9) ${\mathrm{\Phi }}_B(\zeta )={\int }_a^b\log |\mathit{M}(\zeta )|f(\rho )\mathrm{d}\rho .$(9) ) with the sample mean. All computations are carried out in MATLAB R2010a (64bt) using a computer with 2.40GHz Intel processor and 4GB RAM. The Matlab function ‘fmincon’ is used to find the maximum of (9(9) ${\mathrm{\Phi }}_B(\zeta )={\int }_a^b\log |\mathit{M}(\zeta )|f(\rho )\mathrm{d}\rho .$(9) ) and the optimal design. The aim of Bayesian design criterion is to find a design which maximises ${\mathrm{\Phi }}_B(\zeta )$, which is equivalent to find a design which minimises the opposite of ${\mathrm{\Phi }}_B(\zeta )$. ‘Fmincon’ just attempts to find a constrained minimum of a scalar function of several variables starting at an initial estimate. The equidistant design is used as the starting points in our program, and the optimal design is not sensitive to the starting points. The computing time to generate a Bayesian optimal design is about 3 hours, and the running time increases with the increasing number of time points.

For longitudinal data, the more measurements per individual are taken, the more additional information provides. However, the number of measurements per individual is usually limited because of the experimental sampling budget or logistical restrictions. The linear cost function defined by Tekle et al. (2008) is expanded to the population design ζ for the model with r responses (10) $C=\sum _{q=1}^Q{n}_q(c_1+c_{2}r{m}_q),$(10) where $c_1$ is the cost of recruiting a new individual, $c_2$ represents the average cost per measurement, and $\sum _{q=1}^Q{n}_q=N$.

Thus, as proposed in Fedorov et al. (2002), the Fisher information matrix in (7(7) $\mathit{M}(\zeta )={\mathit{X}}^T{\mathit{H}}^{-1/2}\mathit{W}{\mathit{H}}^{-1/2}\mathit{X}.$(7) ) under the population design ζ of the form (5(5) $\zeta =\left\{\begin{array}{cccc}{\mathit{t}}_1& {\mathit{t}}_2& \cdots & {\mathit{t}}_Q\\ w_1& w_2& \cdots & w_Q\end{array}\right\},$(5) ) should be normalised for the cost constraint in (10(10) $C=\sum _{q=1}^Q{n}_q(c_1+c_{2}r{m}_q),$(10) ), that is (11) $\tilde{\mathit{M}}(\zeta )={\displaystyle \frac{1}{C}\sum _{q=1}^Q{w}_q{\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_q)\triangleq \sum _{q=1}^Q{\tilde{w}}_q{\tilde{\mathit{M}}}_{\mathrm{ind}}({\mathit{t}}_q),}$(11) where ${\tilde{\mathit{M}}}_{\mathrm{ind}}({\mathit{t}}_q)={\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_q)/(c_1+c_{2}r{m}_q)$ is the standard information matrix at the individual design ${\mathit{t}}_q$ and ${\tilde{w}}_q={\displaystyle \frac{n_q(c_1+c_{2}r{m}_q)}{CN}={\displaystyle \frac{n_q(c_1+c_{2}r{m}_q)}{N\sum _{q=1}^Q{n}_q(c_1+c_{2}r{m}_q)}.}}$ The optimal population design problem is to identify the optimal weight, and the locations of the time points of each individual design. The value of frequency $n_q$ of the individual design ${\mathit{t}}_q$ can be obtained by rounding ${\tilde{w}}_{q}CN/(c_1+c_{2}r{m}_q)$ to the nearest integer subject to $\sum _{q=1}^Q{n}_q(c_1+c_{2}r{m}_q)\le C$. For comparing the efficiencies of two population designs ${\zeta }_1$ and ${\zeta }_2$ under the cost constraint, the D-efficiency of ${\zeta }_1$ relative to ${\zeta }_2$ is defined by (12) $RE({\zeta }_1;{\zeta }_2)={\left({\displaystyle \frac{|\tilde{\mathit{M}}({\zeta }_1)|}{|\tilde{\mathit{M}}({\zeta }_2)|}}\right)}^{1/p}.$(12) For the normalised Fisher information matrix $\tilde{\mathit{M}}(\zeta )$, the general equivalence theorem still holds, and the verification is given in Section 4 by a simulation study.

Theorem 3.2

For a given linear cost function, a population design ${\zeta }^{\ast }$ is D-optimal for the estimation of the fixed-effects vector $\mathit{\beta }$ if and only if

1. The design ${\zeta }^{\ast }$ maximises $\log |\tilde{\mathit{M}}(\zeta )|.$

2. The design ${\zeta }^{\ast }$ minimise $\underset{\mathit{t}}{max}\mathrm{tr}({\tilde{\mathit{M}}}^{-1}(\zeta ){\tilde{\mathit{M}}}_{\mathrm{ind}}(\mathit{t})).$

3. The maximum of $\mathrm{tr}({\tilde{\mathit{M}}}^{-1}(\zeta ){\tilde{\mathit{M}}}_{\mathrm{ind}}(\mathit{t}))$ is attained at the support individual designs of ${\zeta }^{\ast },$ and the maximum value is p.

4. Numerical studies

In this section, a linear and quadratic mixed effects model is considered to show the effects of autocorrelation coefficient ρ on optimal designs. The optimal designs for the bivariate AIDS model studied in Wang and Fan (2010) is also discussed.

4.1. Simulation example

Consider a linear and quadratic mixed effects model (13) $\begin{array}{r}\{\begin{array}{l}{y}_1={\beta }_{10}+{\beta }_{11}t+b_{10}+e_1,\\ y_2={\beta }_{20}+{\beta }_{21}t+{\beta }_{22}{t}^2+b_{20}+b_{21}t+e_2.\end{array}\end{array}$(13) Assuming that each individual is observed under time schedule ${\mathit{t}}_1=(t_1,t_2,\dots ,t_m)$. The unknown parameters are given previously to find the optimal population design. Let $\mathit{\beta }=(1.5,-3,-1,2,3{)}^T$, and the covariance matrix of the two responses $y_1$ and $y_2$ is $\mathit{G}=\left(\begin{array}{cc}5& {\gamma }_{12}\sqrt{10}\\ {\gamma }_{12}\sqrt{10}& 2\end{array}\right).$ The covariance matrix of the three mixed effects variables $b_{10},b_{20},b_{21}$ is $\mathit{D}=\left(\begin{array}{ccc}2& 0.1& 0.5\\ 0.1& 4& 0.9\\ 0.5& 0.9& 3\end{array}\right).$ The design region is $[-1,1]$.

We first show that the correlation between the two characteristics has slightly effect on the optimal designs. Fixing $\rho =0.2$, the individual D-optimal design for different values of m = 3, 4, 5, 6 and ${\gamma }_{12}=0$, $\pm 0.02$, $\pm 0.1$, $\pm 0.3$, $\pm 0.6$, $\pm 0.9$ are shown in Figure 1. All the optimal designs have design points located in the left and right endpoints of the design area, for value of ${\gamma }_{12}$ in the neighbourhood of 0 that the optimal designs are almost the same, with the increase of absolute value of ${\gamma }_{12}$, there are some subtle differences between the allocation of the optimal designs. In addition, the relative D-efficiency of the individual optimal designs at ${\gamma }_{12}\ne 0$ compared to ${\gamma }_{12}=0$ are all almost 1, whatever the number of design points is. Hence, the correlation between characteristics has slight effects on the D-optimal designs, which is consistent with the results in Atkinson et al. (2007). Hence, we take ${\gamma }_{12}=0.02$ in the model (13(13) $\begin{array}{r}\{\begin{array}{l}{y}_1={\beta }_{10}+{\beta }_{11}t+b_{10}+e_1,\\ y_2={\beta }_{20}+{\beta }_{21}t+{\beta }_{22}{t}^2+b_{20}+b_{21}t+e_2.\end{array}\end{array}$(13) ) and consider the effects of ρ on D-optimal designs in what follows.

Figure 1. The individual D-optimal design allocations of time points for the linear and quadratic mixed effects model (13(13) $\begin{array}{r}\{\begin{array}{l}{y}_1={\beta }_{10}+{\beta }_{11}t+b_{10}+e_1,\\ y_2={\beta }_{20}+{\beta }_{21}t+{\beta }_{22}{t}^2+b_{20}+b_{21}t+e_2.\end{array}\end{array}$(13) ) for different values of m and ${\gamma }_{12}$.

The individual D-optimal allocations of time points for different autocorrelation coefficient ρ $(\rho =0,0.001,0.01,0.1,0.2,\dots ,0.9)$ are presented in Figure 2, the number of design points are $m=3,4,5,$ and 6, respectively. By assuming that ρ obeys an uniform distribution on $[0.05,0.95]$, the Bayesian optimal designs (9(9) ${\mathrm{\Phi }}_B(\zeta )={\int }_a^b\log |\mathit{M}(\zeta )|f(\rho )\mathrm{d}\rho .$(9) ) are also calculated, and the results are also shown in Figure 2. The asterisks show the D-optimal allocations of time points when $\rho =0$, the filled dots are the design position for $\rho \ne 0$ and the plus signs are the design position of Bayesian optimal designs. It can be observed that when the number of design points is three, the allocation of the optimal designs are almost the same. While for $m\ge 4$, when ρ equals zero, the optimal design points clustered at the left, middle and right endpoints of the design area, which maybe infeasible and less informative in practice. With the increase of ρ, the locations of the D-optimal design points are moderately affected by the autocorrelation coefficient ρ, especially for large values of ρ. And there isn't much difference between the Bayesian optimal design and locally D-optimal design when $\rho >0.1$.

Figure 2. The individual optimal design allocations of time points for the linear and quadratic mixed effects model (13(13) $\begin{array}{r}\{\begin{array}{l}{y}_1={\beta }_{10}+{\beta }_{11}t+b_{10}+e_1,\\ y_2={\beta }_{20}+{\beta }_{21}t+{\beta }_{22}{t}^2+b_{20}+b_{21}t+e_2.\end{array}\end{array}$(13) ) for different values of m and ρ.

4.2. Optimal designs for the AIDS model

The bivariate longitudinal AIDS model is based on a set of data in the AIDS Clinical Trials Group Study 175, and only the post-week-12 data were considered, in which the CD4 and CD8 cell counts (per cubic millimetre) of a total of 20 patients were recorded per protocol about every 12 weeks thereafter on each patience visit. Let ${\mathit{y}}_{ij}=({\log }_{10}(\text{CD4}{)}_{ij},{\log }_{10}(\text{CD8}{)}_{ij}{)}^T$ denote the ${\log }_{10}$-transformed CD4 and CD8 cell counts for patient i at the scaled time $t_{ij}=({\text{week}}_{ij}-8)/12\in [1,12]$, $i=1,2,\dots ,20$, $j=1,2,\dots ,n_i$. The model considered here is assumed to be the following two-response linear mixed model with AR(1) errors (14) $\{\begin{array}{l}{y}_{ij1}={\beta }_{10}+{\beta }_{11}{t}_{ij}+b_{i10}+b_{i11}{t}_{ij}+e_{i1j},\\ y_{ij2}={\beta }_{20}+{\beta }_{21}{t}_{ij}+b_{i20}+b_{i21}{t}_{ij}+e_{i2j}.\end{array}$(14) According to Table 1 in Wang and Fan (2010), the maximum likelihood estimates of the fixed-effects vector $\mathit{\beta }$ via the ECM algorithm based on the data is $\hat{\mathit{\beta }}=(2.5708,-0.0236,2.8627,-0.0110{)}^T$, the estimate of the autocorrelation coefficient ρ is 0.3821, and the estimate of the covariance matrix of two response variables is $\mathit{G}=\left(\begin{array}{cc}cc0.0083& 0.0061\\ 0.0061& 0.0111\end{array}\right)$.

We now consider the D-optimal single-group design for model (14(14) $\{\begin{array}{l}{y}_{ij1}={\beta }_{10}+{\beta }_{11}{t}_{ij}+b_{i10}+b_{i11}{t}_{ij}+e_{i1j},\\ y_{ij2}={\beta }_{20}+{\beta }_{21}{t}_{ij}+b_{i20}+b_{i21}{t}_{ij}+e_{i2j}.\end{array}$(14) ), where all individuals are observed under the same settings, that is, ${\zeta }_m^{\ast }=\left\{\begin{array}{c}{\mathit{t}}_m^{\ast }\\ 1\end{array}\right\},$ where ${\mathit{t}}_m^{\ast }$ is a single sequence of m time points in the design region $[1,12]$. Suppose the experimental cost to be considered in an experiment. The relative efficiency defined in (12(12) $RE({\zeta }_1;{\zeta }_2)={\left({\displaystyle \frac{|\tilde{\mathit{M}}({\zeta }_1)|}{|\tilde{\mathit{M}}({\zeta }_2)|}}\right)}^{1/p}.$(12) ) is calculated for the cost of recruiting a new individual, $c_1=100$, and the average cost per measurement, $c_2$, ranging from 0.01 to 20. The D-optimal designs ${\zeta }_m^{\ast }$ are determined by Theorem 3.2. Figure 3 shows the relative efficiencies of ${\zeta }_m^{\ast }$ ($3\le m\le 12$) compared to ${\zeta }_{12}^{\ast }$ for $c_1=100$ and $c_2\in \{0.01,0.1,1,2,5,20\}$. In each curve corresponding to a value of $c_2$, the optimal number of time points is the number m at which the relative efficiency is the highest. It is observed that the number of optimal design points is m = 3 for $c_2=1$, 2, 5, 20. And the performances of the designs, ${\zeta }_m^{\ast }$ ($3\le m\le 12$), are almost the same for $c_2=0.01$ and 0.1.

Figure 3. Relative efficiencies of ${\zeta }_m^{\ast }$ relative to ${\zeta }_{12}^{\ast }$ for the two-response linear mixed model (14(14) $\{\begin{array}{l}{y}_{ij1}={\beta }_{10}+{\beta }_{11}{t}_{ij}+b_{i10}+b_{i11}{t}_{ij}+e_{i1j},\\ y_{ij2}={\beta }_{20}+{\beta }_{21}{t}_{ij}+b_{i20}+b_{i21}{t}_{ij}+e_{i2j}.\end{array}$(14) ) with costs $c_1=100$ and $c_2\in \{0.01,0.1,1,2,5,20\}$.

We now consider more-group designs. The set of possible sampling times considered is a predefined set T (Castañeda & López-Rios, 2016). For example, we let T be the 12-point single group balanced D-optimal design for model (14(14) $\{\begin{array}{l}{y}_{ij1}={\beta }_{10}+{\beta }_{11}{t}_{ij}+b_{i10}+b_{i11}{t}_{ij}+e_{i1j},\\ y_{ij2}={\beta }_{20}+{\beta }_{21}{t}_{ij}+b_{i20}+b_{i21}{t}_{ij}+e_{i2j}.\end{array}$(14) ), that is $\begin{array}{rl}{\mathit{t}}_{12}^{\ast }& =\{1,1.563,2.162,2.809,3.519,4.319,5.268,6.531,\\ & 9.007,10.249,11.196,12\}.\end{array}$ In practice, experimenter can choose a set of possible sampling times by their experience or some prior information (Nyberg et al., 2012). We assume that three or four measurements are available for each individual, and there is no replication at a certain time point. It is noted from Figure 3 that a large value of $c_2$ results in a small sample size, and we then consider the cost function (10(10) $C=\sum _{q=1}^Q{n}_q(c_1+c_{2}r{m}_q),$(10) ) where $c_1=100$, $c_2=1.7$, for example. We find the D-optimal population design as follows: ${\zeta }^{\ast }=\left\{\begin{array}{cc}(1,2.809,12)& (1,2.162,4.319,12)\\ 0.2352& 0.7648\end{array}\right\},$ which can be verified by Theorem 3.2. Therefore, according to Gagnon and Leonov (2005) and (11(11) $\tilde{\mathit{M}}(\zeta )={\displaystyle \frac{1}{C}\sum _{q=1}^Q{w}_q{\mathit{M}}_{\mathrm{ind}}({\mathit{t}}_q)\triangleq \sum _{q=1}^Q{\tilde{w}}_q{\tilde{\mathit{M}}}_{\mathrm{ind}}({\mathit{t}}_q),}$(11) ), the ratio of the frequencies is $\frac{N_1}{N_2}=\frac{w_1}{w_2}\cdot \frac{c_1+2c_2{m}_2}{c_1+2c_2{m}_1}=0.3170,$ which means that about $24.07\mathrm{\%}$ of the population are observed at $\{1,2.8090,12\}$, and about $75.93\mathrm{\%}$ should be measured at $\{1,2.1623,4.3190,12\}$. Obviously, this population optimal design is more efficient than the single design ${\xi }_{12}^{\ast }$ when considering the cost.

On the other hand, the result is a little different from Schmelter (2007a, 2007b), and each individual design involved in ${\zeta }^{\ast }$ is no longer D-optimal. When the value of $c_2$ varies, the individual designs remain the same, while the weight of each individual design changes, and the population design tends to single three or four points single-group design.

5. Conclusion

The problem of choosing the optimal population design is investigated for the multi-response linear mixed effects models with correlated errors. AR(1) structure is used to model the time series autocorrelation of the errors. A higher order AR may be worth consideration, however, it brings more unknown parameters into the model. Simulation study is used to show the effect of autocorrelation on designs and parameter estimations. The equivalence theorems are provided to characterise the optimal population designs for the estimation of the fixed effects.

It should be noticed that the autocorrelation coefficient and the variance parameters in the random effect distribution are pre-specified in finding the D-optimal designs. In practice, the experimenter may get some prior information of the parameters, such as the range of parameter values. Then, a semi-Bayesian criterion is used in combination with the D-optimal criterion to find the optimal design robust to the autocorrelation coefficient. A topic for further research may be to construct optimal designs with respect to other optimality criteria, such as T-optimality criterion, or compound criterion. The result may be extended to polynomial mixed models with measurement error, and also to generalised linear mixed effects models.

Acknowledgments

This work was partly supported by the National Natural Science Foundation of China (Nos. 11971318, 11871143) and Shanghai Rising-Star Program (No. 20QA1407500).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was partly supported by the National Natural Science Foundation of China (Nos. 11971318, 11871143) and Shanghai Rising-Star Program (No. 20QA1407500).

Notes on contributors

Hongyan Jiang

Hongyan Jiang is a Ph.D. in Statistics, Shanghai Normal University, 2019, Associate Professor in Huaiyin Institute of Technology, Jiangsu, 2019–Present.

Rongxian Yue

Rongxian Yue is a Ph.D. in Statistics, Hong Kong Baptist University, 1997. Postdoctoral Research Fellow in Statistics, East China Normal University, 1997–1999. Associate Professor (1999–2001) and Professor (2001–Present), Shanghai Normal University.