Optimal Designs in Multiple Group Random Coefficient Regression Models
Maryna Prus




Optimal Designs in Multiple Group Random Coefficient Regression Models 

Maryna Prus[Maryna Prus: mailto:maryna.prus@ovgu.demaryna.prus@ovgu.de]


Otto-von-Guericke University Magdeburg, Institute for Mathematical Stochastics, 


PF 4120, D-39016 Magdeburg, Germany





Abstract: The subject of this work is multiple group random coefficients regression models with several treatments and one control group. Such models are often used for studies with cluster randomized trials. We investigate A-, D- and E-optimal designs for estimation and prediction of fixed and random treatment effects, respectively, and illustrate the obtained results by numerical examples.

Keywords: Optimal design, treatment and control, random effects, cluster randomization, mixed models, estimation and prediction 





§ INTRODUCTION


Random coefficient regression (RCR) models were originally introduced in plant and animal breeding and used for selection purposes (see e. g. <cit.>). 
The subject of this paper is the multiple group random coefficient models, in which observational units (individuals) are allocated in several treatment groups and one control group. In each treatment group some group-special kind of treatment is available, in the control group there is no treatment. Such models are typically used for cluster randomization or cluster randomized trials (see e. g. <cit.> or <cit.>). 

RCR models with known population parameters were considered in detail by <cit.>. 
<cit.> investigated optimal designs for estimation of unknown population mean parameters in RCR models, where all individuals are observed under the same regime. Analytical results for designs, which are optimal for the prediction of individual random parameters in hierarchical models with the same treatment for all individuals, have been presented in <cit.>. A practical approach for computation of optimal approximate and exact designs was proposed by <cit.>.

The fixed effects version of the multiple group models considered in this paper may be recognized as the well known one way layout model. Classical one-way layout models have been well discussed in the literature. Results for the optimal designs can be found e. g. in <cit.>, <cit.>, <cit.>, <cit.> or <cit.>.

Optimal designs for estimation of fixed parameters in multiple group models with random coefficients were considered e. g. in <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. <cit.> worked on optimal designs, which minimize a loss function in multicentre trials models. In <cit.> models with the same fixed number of observations in all groups were investigated. <cit.> proposed a design optimization method based on the generalized least squares estimation. <cit.> considered models with carryover effects. In <cit.> optimal designs were computed for the maximum likelihood estimation.

Optimal designs for prediction of random effects in models, where the population mean parameters differ from group to group, were briefly discussed in <cit.>. 
 
In this paper we investigate multiple group models with the same unknown population parameters across all groups. We present analytical results for A-, D- and E-optimal designs (-optimal group sizes) based on the best linear unbiased estimation or prediction of fixed or random treatment effects, respectively.

The paper has the following structure: In Section 2 the model will be specified. In Section 3 the best linear unbiased estimation for the population parameters (mean treatment effects) and a best linear unbiased prediction for the random treatment effects of the observational units will be discussed. Section 4 provides analytical results for designs, which are optimal for estimation or prediction. The results will be illustrated by a numerical example, in which we compare the optimal group sizes in the model under investigation with optimal group sizes in the fixed effects model (one-way-layout). The paper will be concluded by a discussion of the obtained results and possible directions for the future research in Section 5.



§ MODEL SPECIFICATION

We consider here a multiple group model with J groups and N individuals. In the first J-1 (treatment) groups individuals get group-special kinds of treatment: 1, …, J-1. Each treatment group includes n>0 individuals. The last group (group J) is a control group (no treatment) with m=N-(J-1)n, m>0 individuals. 

The k-th observation at individual i in the j-th treatment group is given by the following formula:

    Y_jik=μ_i+α_ji + ε_jik,   k=1, …, K,   i=(j-1)n+1, …, jn,   j=1, …, J-1,

while in the control group the k-th observation at the i-th individual is given by

    Y_Jik=μ_i + ε_Jik,   k=1, …, K,   i=(J-1)n+1, …, N,

where K is the number of observations per individual, which is assumed to be the same for all individuals across all groups, ε_jik are observational errors with zero expected value and common variance σ^2>0.

In this work we optimize the numbers of individuals in treatment and control groups. Therefore the group allocation of individuals is not completely clear. For this reason we define the individual treatment effects α_ji for all individuals and all treatments, i. e. i=1, …, N and j=1, …, J-1. For individual i in the control group the treatment effects α_ji would appear if the individual would be treated with treatment j. For individual i in treatment group j the treatment effect α_j'i would appear if the individual would get treatment j' instead of j for j, j'=1, …, J-1, j≠ j'.

Individual intercepts μ_i are defined for all individual, i. e. i=1, …, N. μ_i and α_ji have unknown expected values E(μ_i)=μ, E(α_ji)=α_j and variances var(μ_i)=u σ^2, var(α_ji)=v σ^2 for some given positive values of u and v. All individual parameters μ_i and α_ji' and all observational errors ε_j'i”k, for k=1, …, K, i, i'=1, …, N and j, j', j”=1, …, J, are uncorrelated.























For further considerations we define the regression functions

    𝐟(j):={[   (1,𝐞_j^⊤)^⊤ ,     j=1, …, J-1; (1,0_J-1^⊤)^⊤ ,             j=J ].,

where 𝐞_j and 0_J-1 denote the j-th unit vector and the zero vector of length J-1, respectively, and the total number of individuals

    N_j:={[         nj, j=1, …, J-1;          N,         j=J ].

in the first j groups, and we fix N_0 by 0. Then for the vector _i:= (μ_i, α_1 i, …, α_J-1 i)^⊤ of individual random parameters the multiple group model (defined by (<ref>) and (<ref>)) may be rewritten in the following form:

    Y_jik=𝐟(j)^⊤_i + ε_jik,   k=1, …, K,   i=N_j-1+1, …, N_j,   j=1, …, J.

The parameter vectors _i have the expected value E(_i)=_0:=(μ, α_1, …, α_J-1)^⊤ and the covariance matrix Cov(_i)=σ^2block-diag(u, v 𝕀_J-1), where 𝕀_p denotes the p× p identity matrix and block-diag(𝐀_1, …, 𝐀_p) is the block diagonal matrix with l× l blocks 𝐀_1, …, 𝐀_p.

On the individual level the vectors 𝐘_j,i=(Y_ji1, …, Y_jiK)^⊤ for individuals in the j-th group can be specified as

    𝐘_j,i=1_K 𝐟(j)^⊤_i + _j,i,   i=N_j-1+1, …, N_j,   j=1, …, J,

where 1_K denotes the vector of length K with all entries equal to 1 and _j,i=(ε_ji1, …, ε_jiK)^⊤.


On the group level for the group parameters _j:= (_N_j-1+1^⊤, …, _N_j^⊤)^⊤ the vector 𝐘_j=(𝐘_j, (j-1)n+1^⊤, …, 𝐘_j, N_j^⊤)^⊤ of all observations at all individuals in the j-th group is given by




    𝐘_j=(𝕀_r_j⊗(1_K 𝐟(j)^⊤))_j + _j,   j=1, …, J,

where 

    r_j:={[         n , j=1, …, J-1;         m ,         j=J ].,

is the number of individuals in the j-th group, “⊗” denotes the Kronecker product and _j=(_j, (j-1)n+1^⊤, …, _j, N_j^⊤)^⊤.


Finally, we introduce the vector of all individual random parameters = (_1^⊤, …,_J^⊤)^⊤ (or, equivalently, := (_1^⊤, …, _N^⊤)^⊤), for which 



the full vector 𝐘=(𝐘_1^⊤, …, 𝐘_J^⊤)^⊤ of observations at all individuals in all groups is given by the following formula:

    𝐘=Diag_j=1^J(𝕀_r_j⊗(1_K 𝐟(j)^⊤)) + ,

where Diag_s=1^p(𝐀_s) is the block diagonal matrix with l_s× t_s blocks 𝐀_s, s=1, …, p,






 and =(_1^⊤, …, _J^⊤)^⊤.



Alternatively, for the random vector :=-(1_N⊗𝕀_J) _0 the model (<ref>) can be represented in the form

    𝐘=Vec_j=1^J(1_r_j⊗(1_K 𝐟(j)^⊤))_0+Diag_j=1^J(𝕀_r_j⊗(1_K 𝐟(j)^⊤)) + ,

where Vec_s=1^p(𝐀_s)=(𝐀_1^⊤, …, 𝐀_p^⊤)^⊤ for some l_s× t matrices 𝐀_s, s=1, …, p. 

Note that = (_1^⊤, …, _N^⊤)^⊤ for _i=_i-_0, i=1, …, N. 
The expected value of  is zero and Cov()=σ^2 𝕀_N⊗block-diag(u,  v 𝕀_J-1).



§ ESTIMATION AND PREDICTION


In this chapter we determine the best linear unbiased estimator (BLUE) for the population mean parameters α_j and μ and the best linear unbiased predictor (BLUP) for the individual random parameters α_ji and μ_i for j=1, …, J-1 and i=1, …, N. 

We denote the mean observation in group j by Y̅_j=1/r_j∑_i=N_j-1+1^N_jY̅_j,i and the mean observation at individual i in group j by Y̅_j,i=1/K∑_k=1^KY_jik for all j=1, …, J, i=N_j-1+1, …, N_j. Then we obtain the following results for the BLUEs of the mean parameters μ and α_j and the BLUPs of the treatment effects μ_i and α_ji.

The BLUE for the population intercept parameter μ is given by 

    μ̂=Y̅_J

and the BLUE for the population treatment effect α_j is given by

    α̂_̂ĵ= Y̅_j-Y̅_J,   j=1, …, J-1.



If the i-th individual is in the control group, the BLUP for the individual intercepts μ_i is given by

    μ̂_i=Ku/Ku+1Y̅_J,i+1/Ku+1Y̅_J,   i=N_J-1+1, …, N,

otherwise, if the i-th individual is in the j-th treatment group, the BLUP is given by

    μ̂_i=Ku/K(v+u)+1 (Y̅_j,i-Y̅_j)+Y̅_J,   j=1, …, J-1,    i=N_j-1+1, …, N_j.

If the i-th individual is in the j-th treatment group, the BLUP for the individual treatment effect α_ji is given by

    α̂_ji=Kv/K(v+u)+1 (Y̅_j,i-Y̅_J)+Ku+1/K(v+u)+1 (Y̅_j-Y̅_J),   j=1, …, J-1,   i=N_j-1+1, …, N_j

otherwise the BLUP is given by

    α̂_ji=Y̅_j-Y̅_J,   j=1, …, J-1,   i=1, …, N,    i≠ N_j-1+1, …, N_j.




For the vector _0:=(α_1, …, α_J-1)^⊤ of the mean treatment effects the BLUE is given by _0=(α̂_1, …, α̂_J-1)^⊤. 








We denote the vector of all individual treatment effects by :=(_1^⊤, …, _N^⊤)^⊤, where _i:=(α_1 i, …, α_J-1 i)^⊤ is the vector of individual treatment effects for the i-th individual, i=1, …, N. Then the BLUP of  is given by =(_1^⊤, …, _N^⊤)^⊤, where _i=(α̂_1 i, …, α̂_J-1 i)^⊤ is the BLUP for _i, i=1, …, N.  


The next theorems provide the covariance matrix of _0 and the MSE matrix of .













The covariance matrix of the BLUE _0 is given by

    Cov(_0)=σ^2(Ku+1/Km 1_J-1 1_J-1^⊤+K(v+u)+1/Kn 𝕀_J-1).



The mean squared error matrix of the BLUP  is given by

    Cov(-)= 𝐁_1+𝐁_2+𝐁_2^⊤+𝐁_3,

where

    𝐁_1=σ^21_N1_N^⊤⊗(Ku+1/Km 1_J-1 1_J-1^⊤+K(v+u)+1/Kn 𝕀_J-1),


    𝐁_2=-σ^2v 1_N⊗(tVec_j=1^J-1(1/n1_n^⊤⊗(𝐞_j𝐞_j^⊤))  ⋮  0_(J-1)× m(J-1)),

where 0_p× q denotes the p× q zero matrix and tVec_s=1^p(𝐀_s)=(Vec_s=1^p(𝐀_s^⊤))^⊤, and

    𝐁_3=σ^2v block-diag(𝕀_n(J-1)^2-Kv/K(v+u)+1Diag_j=1^J-1((𝕀_n-1/n1_n1_n^⊤)⊗(𝐞_j𝐞_j^⊤) ), 𝕀_m(J-1)).






For the proofs of Theorems <ref>-<ref> see Appendix <ref>.







§ OPTIMAL DESIGN


In this chapter we optimize the numbers n and m of individuals in the treatment and control groups, respectively. We define the exact experimental design as


    ξ = ([   1   … J-1   J;   n   …   n   m ]),

where the indexes 1, …, J-1 denote the treatment groups and the index J is used for the control group. 

For analytical purposes, we also define the approximate design:

    ξ = ([        1        …      J-1        J;        w        …        w 1-(J-1)w ]),   w∈(0, 1/J-1),

where w=n/N is the weight of a treatment group and 1-(J-1)w=m/N is the weight of the control group. Then only the optimal weight w^* of a treatment group has to be determined. 




 §.§ A-criterion


For the estimation of the population treatment effects _0 the A-criterion for an exact design is defined as the trace of the covariance matrix of the BLUE _0:

    A_Ψ_0(ξ) = tr( Cov(_0) ).

We determine the trace of the covariance matrix (<ref>),

replace n by Nw and m by N(1-w) and obtain the next results for an approximate design.

The A-criterion for the estimation of the population treatment effects _0 is given for an approximate design by

    A_Ψ_0(w) = σ^2(J-1)/K N(K(v+u)+1/w+Ku+1/1-(J-1)w).



The A-optimal weight w^*_A,Ψ_0 for the estimation of the population treatment effects is given by

    w^*_A,Ψ_0 = 1/J-1+√(J-1)√(Ku+1/K(v+u)+1).


Note that for large values of the intercepts variance (u→∞) the optimal weight (<ref>) tends to the value w^*_A,Ψ_0 = 1/J-1+√(J-1), which coincides with the optimal weight in the fixed effects model (see e.g. <cit.>, ch. 3 or <cit.>, ch. 3) . If the treatment effects variance takes a very large value (v→∞), the limiting optimal design assigns all observations to be taken in the treatment groups: w^*_A,Ψ_0 = 1/J-1. If both variances are large and the variance ratio b=v/u is fixed, the limiting optimal design is given by 

    w^*_A,Ψ_0 = 1/J-1+√(J-1)√(1/1+b).


The A-criterion for the prediction of the individual treatment effects _i is defined for an exact design as the trace of the mean squared error matrix of the BLUP :

    A_Ψ(ξ) = tr( Cov(-) ).

The next theorem presents the A-criterion for an approximate design.

The A-criterion for the prediction of the individual treatment effects _i is given for an approximate design by

    A_Ψ(w) = σ^2(J-1)(K(v+u)+1/Kw+Ku+1/K(1-(J-1)w)+v(N-2-Kv (Nw-1)/K(v+u)+1)).



The result of Theorem <ref> follows from (<ref>) and

    tr(𝐁_1) = σ^2(J-1)/K(K(v+u)+1/w+Ku+1/(1-(J-1)w)),


    tr(𝐁_2) = -σ^2(J-1) v

and

    tr(𝐁_3) = σ^2(J-1) v(N-Kv (Nw-1)/K(v+u)+1).


Note that there is no explicit formula for the optimal weight in this case. However, it is easy to see that there is a unique solution w^*, which may be determined numerically for given values of N, J, K, u and v. To illustrate the behavior of optimal designs, we consider a numerical example.

Example 1. 
Let the total number of individuals be N=100 and the number of observations per individual be K=10 and the variance ratio b=v/u be fixed by the values 2, 0.6 and 0.001. The next graphics (Figure 1 and Figure 2) illustrate the behavior of the A-optimal weight in dependence on the treatment effects variance for the special cases of one (left panel) and two (right panel) treatment groups (J=2 and J=3, respectively). The parameter ρ=v/(1+v) is used instead of the variance parameter v to cover all values of the treatment effects variance by a finite interval. On the graphics the solid, dashed and dotted lines present the optimal weight for the values 2, 0.6 and 0.001 of the ratio b. Note that the optimal weight w^* takes all its values in the intervals (0,1) and (0,0.5) in the models with one and two treatment groups, respectively. 
 
			
For the models with one and two treatment groups the optimal the optimal weights start (for ρ→ 0) at points w^*=0.5 and w^*≈0.29, respectively. This may be explained by the fact that the optimal designs for the fixed effects models are equal to w^*_A,fix = 1/J-1+√(J-1) and result in w^*_A,fix = 0.5 for J=2 and w^*_A,fix = 1/2+√(2)≈0.29 for J=3. The optimal weights increase with increasing variance of the individual treatment effects with limiting values (for ρ→ 1) w^*=0.91, w^*=0.80 and w^*≈0.50 for J=2 and w^*=0.44, w^*=0.38 and w^*≈0.29 for J=3 for b=2, b=0.6 and b=0.001, respectively.

Figure 3 and Figure 4 present the efficiency of the optimal weight w^*_A,fix from the fixed effects model for the present model for one (left panel) and two (right panel) treatment groups for the values 2, 0.6 and 0.001 of the ratio b.

 
		
For both particular models the efficiencies start at point 1 and decrease with limits eff=0.65, eff=0.90 and eff≈1 for J=2 and eff=0.89, eff=0.98 and eff≈1 for J=3 for b=2, b=0.6 and b=0.001.



 §.§ D- and E-criterion

In this section we consider D- and E-optimality criteria for the estimation and the prediction in multiple group models. We consider the general case of model (<ref>) for the estimation of population parameters and we restrict ourselves to the special case J=2 for the prediction of individual treatment effects.

For further considerations we will use the following result.


The eigenvalues of the covariance matrix of the BLUE _0 are

    λ_1=σ^2/K((Ku+1)(J-1)/(N-(J-1) n)+K(v+u)+1/n)

with algebraic multiplicity 1 and

    λ_2=σ^2(K(v+u)+1)/K n

with algebraic multiplicity J-2.



To determine the eigenvalues of the covariance matrix Cov(_0) we solve the equation

    det( Cov(_0)-λ 𝕀_J-1)=0,

where λ denotes an eigenvalue of Cov(_0).

    det( Cov(_0)-λ 𝕀_J-1)
            =   det(σ^2(Ku+1)/K(N-(J-1) n)1_J-1 1_J-1^⊤+(σ^2(K(v+u)+1)/K n-λ) 𝕀_J-1)
            =   (σ^2/K((Ku+1)(J-1)/N-(J-1) n+K(v+u)+1/n)-λ)(σ^2(K(v+u)+1)/K n-λ)^J-2.

Then

    λ_1=σ^2/K((Ku+1)(J-1)/(N-(J-1) n)+K(v+u)+1/n)

and

    λ_2=σ^2(K(v+u)+1)/K n

are the solutions of (<ref>).


For the estimation of the population treatment effects the D-criterion is defined as the logarithm of the determinant of the covariance matrix of the BLUE _0:

    D_Ψ_0(ξ) = ln det( Cov(_0) ).

for an exact design. We compute the determinant as the product of the eigenvalues, which are given in Lemma <ref>, and receive using n=Nw the following result for approximate designs.


The D-criterion for the estimation of the population treatment effects _0 is given for an approximate design by

    D_Ψ_0(w) = c+ln(K(v+u)+1/w+(J-1)(Ku+1)/1-(J-1)w)+(J-2) ln(K(v+u)+1/w),

where c=(J-1) ln (σ^2/K N).











The D-optimal weight w^*_D,Ψ_0 for the estimation of the population treatment effects is given by

    w^*_D,Ψ_0 = J-2+z/(J-1)(J+z),

where z=√(4(J-1)Kv/Ku+1+J^2).





The E-criterion for the estimation of the population treatment effects is defined for an exact design as the largest eigenvalue of the covariance matrix of the BLUE _0

    E_Ψ_0(ξ) = λ_max( Cov(_0) ),

where λ_max(A) denotes the largest eigenvalue of the matrix A.

Using Lemma <ref> we receive the following form of the E-criterion for approximate designs.

The E-criterion for the estimation of the population treatment effects _0 is given for an approximate design by

    E_Ψ_0(w) =σ^2/K N(K(v+u)+1/w+(J-1)(Ku+1)/1-(J-1)w).




The E-optimal weight w^*_E,Ψ_0 for the estimation of the population treatment effects is given by

    w^*_E,Ψ_0 = 1/(J-1)(1+√(Ku+1/K(v+u)+1)).



Note that also for the D- and E-criteria the optimal weights for the estimation of the population parameters ((<ref>) and (<ref>)) tend to the optimal weights in the fixed effects model: w^*_D,Ψ_0→1/J and w^*_E,Ψ_0→1/2(J-1), for u→∞. For large values of the treatment effects variance (v→∞) all observations should to be taken in the treatment groups: w^*_D,Ψ_0 =w^*_E,Ψ_0 = 1/J-1. If both variances are large and b=v/u, the limiting values for the optimal weights are

    w^*_D,Ψ_0→J-2+√(4(J-1)b+J^2)/(J-1)(J+√(4(J-1)b+J^2))
 
and 

    w^*_E,Ψ_0→1/(J-1)(1+√(1/1+b)).


For the prediction of the individual treatment effects we consider the particular multiple group model with one treatment group and one control group (J=2). In this case the mean squared error matrix (<ref>) of the prediction simplifies to

    Cov(-)= ([   𝐇_11   𝐇_12; 𝐇_12^⊤   𝐇_22 ]),

where

    𝐇_11=σ^2(Ku+1)(N/Kn m1_n1_n^⊤+v/K(u+v)+1(𝕀_n-1/n1_n1_n^⊤)),


    𝐇_12=σ^2(Ku+1)N/Kn m1_n1_m^⊤,


    𝐇_22=σ^2((K(u+v)+1/Kn+Ku+1/Km)1_m1_m^⊤+v 𝕀_m).

The approximate design (<ref>) simplifies to

    ξ = ([   1   2;   w 1-w ]).

The next Lemma provides the eigenvalues of the mean squared error matrix (<ref>).

The eigenvalues of the mean squared error matrix of the BLUP  are

    λ_1=σ^2 v (Ku+1)/K(v+u)+1,


    λ_2=σ^2v,


    λ_3=σ^2/2(v/w+Ku+1/Kw(1-w)+√(s_w)/Kw(1-w)),

where s_w=K^2(1-w)^2v^2+2K(1-w)(1-2w)(Ku+1)v+(Ku+1)^2, and

    λ_4=σ^2/2(v/w+Ku+1/Kw(1-w)-√(s_w)/Kw(1-w))

with algebraic multiplicities n-1, m-1, 1 and 1, respectively.

For the proof see Appendix <ref>.

We define the D-criterion for the prediction as the logarithm of the determinant of the mean squared error matrix of :

    D_Ψ(ξ) = ln det( Cov(-) ).

We compute the determinant using the results of Lemma <ref> and obtain the following criterion for approximate designs.

The D-criterion for the prediction of the individual treatment effects _i is given for an approximate design by

    D_Ψ(w) = d+N w ln(Ku+1/K(u+v)+1)-ln(w (1-w)),

where d=ln((σ^2)^Nv^N-1 (K(u+v)+1)1/K).


The D-optimal weight w^*_D,Ψ for the prediction of the individual treatment effects is given by

    w^*_D,Ψ = 1/N t+1/2+√(1/N^2 t^2+1/4) ,

where t=ln(Ku+1/K(u+v)+1).

We define the E-criterion for the prediction as the largest eigenvalue of the mean squared error matrix:

    E_Ψ(ξ) = λ_max( Cov(-) ).

The E-criterion for approximate designs follows directly from Lemma <ref>.

The E-criterion for the prediction of the individual treatment effects _i is given for an approximate design by

    E_Ψ(w) = σ^2/2(v/w+Ku+1/Kw(1-w)+√(s_w)/Kw(1-w)).





The E-optimal weight w^*_E,Ψ for the prediction of the individual treatment effects is given by

    w^*_E,Ψ =K(2u+v)+2/K(4u+v)+4.



Note that for both D- and E-criteria for the prediction the optimal weights tend to w^*_D,Ψ=w^*_E,Ψ=0.5, which is optimal for the fixed effects model (see e.g. <cit.>, p. 44), for u→∞ and to w^*_D,Ψ=w^*_E,Ψ=1 (all observations are being taken in the treatment group) for v→∞. For large u, large v and fixed ratio b=v/u we observe 

    w^*_D,Ψ→ -1/N ln(1+b)+1/2+√(1/N^2(ln(1+b))^2+1/4)

and

    w^*_E,Ψ→2/b+1/4/b+1.


Example 2. We consider again the model with N=100 individuals and K=10 observations per individual. The variance ratio b takes the values 2, 0.6 and 0.001 for both D-criterion and E-criteria. The number of treatment groups is fixed by one (J=2). 
The next picture (Figure 5 and Figure 6) illustrates the behavior of the D- (left panel) and E- (right panel) optimal weights in dependence on the variance parameter ρ. 

 

As we can observe on the picture the optimal weights start with value w^*=0.5, which is optimal for the fixed effects model for both criteria: w^*_D,fix=w^*_E,fix=0.5, and increase with the variance parameter ρ with limiting values w^*=0.99, w^*=0.98 and w^*=0.51 for the D-criterion and w^*=0.67, w^*=0.57 and w^*≈0.5 for the E-criterion for b=2, b=0.6 and b=0.001, respectively. 

On the graphics the values of the optimal designs for the D-criterion are larger than those for the E-criterion for all ρ and b. The A-optimal designs illustrated by Figure 1 turn out to take all their values between the corresponding D- and E-optimal weights. 

Figure 7 and Figure 8 present the efficiency of the optimal designs from the fixed effects model for the present model for the D- and E- criteria.

 

The efficiencies on the picture start at point eff=1 and decrease with increasing variance of the individual treatment effects with limiting values eff=0.60, eff=0.82 and eff≈1 for the D- and eff=0.88, eff=0.98 and eff≈1 for the E-criterion for b=2, b=0.6 and b=0.001, respectively.

According to these graphics and Figure 2 the efficiencies for the D-criterion are smaller than those for the E-criterion and the A-criterion is in between for all values of ρ and b.	

Note also that for small values of the variance ratio b the optimal designs are very close to those in the fixed effects model for all design criteria. Hence, the corresponding efficiencies are close to 1.
		



§ DISCUSSION AND CONCLUSIONS


In the present work multiple group RCR models with several treatment groups and a control group have been considered. We have obtained A-, D- and E-optimality criteria for the estimation of population parameters and for the prediction of individual treatment effects using the covariance matrix of the BLUE and the mean squared error matrix of the BLUP, respectively. The optimal designs (optimal group sizes) turned out to be different for the estimation and the prediction and do not coincide with those in the corresponding fixed-effects-model (one-way-layout). For large values of the treatment effects variance the optimal designs assign almost all observations to be taken in the treatment groups. If the variance of the individual intercepts is large, the optimal groups sizes tend to those in the fixed effects model. 

The optimal group sizes are locally optimal, i. e. they depend on the variance parameters. To avoid this, minimax-optimal designs, which minimize the worst case for the criterion function over some reasonable region of values of the covariance matrix, or some other design criteria, which are robust with respect to the variance parameter, may be considered in the next step of this research. 

For the investigated models we assumed a diagonal covariance matrix of random effects. Models with more complicated covariance structure may also be considered in the future.







§ ACKNOWLEDGMENT
 This research has been supported by grant SCHW 531/16-1 of the German Research Foundation (DFG).



§ APPENDIX



 §.§ Proofs of Theorems <ref>-<ref>

To make use of the available results for estimation and prediction we recognize the model (<ref>) as a special case of the linear mixed model (see e. g. <cit.>)

    𝐘=𝐗 + 𝐙 + ,

where 𝐗 and 𝐙 are known fixed and random effects design matrices,  and  are vectors of fixed and random effects, respectively. The random effects  and the observational errors  are assumed to be uncorrelated and to have zero means and non-singular covariance matrices 𝐆= () and 𝐑= (). 

According to <cit.> for full column rank design matrix 𝐗 the BLUE  for  and the BLUP  for  are provided by the mixed model equations

    ( [ ;  ])
    = ( [      𝐗^⊤𝐑^-1𝐗      𝐗^⊤𝐑^-1𝐙;      𝐙^⊤𝐑^-1𝐗 𝐙^⊤𝐑^-1𝐙+𝐆^-1 ])^-1( [ 𝐗^⊤𝐑^-1𝐘; 𝐙^⊤𝐑^-1𝐘 ]),

which can be rewritten in the alternative form

    =(𝐗^⊤ (𝐙𝐆𝐙^⊤+𝐑)^-1𝐗)^-1𝐗^⊤ (𝐙𝐆𝐙^⊤+𝐑)^-1𝐘,


    =𝐆𝐙^⊤ (𝐙𝐆𝐙^⊤ + 𝐑)^-1(𝐘-𝐗).


The mean squared error matrix of the estimator and predictor (^⊤, ^⊤)^⊤ is given by (see <cit.>)

    ( [  ; - ])
    =( [      𝐗^⊤𝐑^-1𝐗      𝐗^⊤𝐑^-1𝐙;      𝐙^⊤𝐑^-1𝐗 𝐙^⊤𝐑^-1𝐙+𝐆^-1 ])^-1
 
and can be represented as the partitioned matrix

    Cov ( [  ; - ])=(
    [   𝐂_11   𝐂_12; 𝐂_12^⊤   𝐂_22 ]),

where 𝐂_11=Cov(), 𝐂_22=Cov(-),

    𝐂_11=(𝐗^⊤(𝐙𝐆𝐙^⊤+𝐑)^-1𝐗)^-1,


    𝐂_22=(𝐙^⊤𝐑^-1𝐙 +𝐆^-1-𝐙^⊤𝐑^-1𝐗 (𝐗^⊤𝐑^-1𝐗)^-1𝐗^⊤𝐑^-1𝐙)^-1,


    𝐂_12= -𝐂_11 𝐗^⊤𝐑^-1𝐙(𝐙^⊤𝐑^-1𝐙 +𝐆^-1)^-1.


For =_0, =, 𝐗=Vec_j=1^J(1_r_j⊗(1_K 𝐟(j)^⊤)), 𝐙=Diag_j=1^J(𝕀_r_j⊗(1_K 𝐟(j)^⊤)),𝐆=σ^2 𝕀_N⊗block-diag(u, 𝕀_J-1) and 𝐑=Cov()=σ^2 𝕀_NK our model (<ref>) is of form (<ref>). 









Using formulas (<ref>) and (<ref>) and after employing some linear algebra we obtain the following 
BLUE and BLUP for the fixed and random effects _0 and :

    _0=([                   𝐘̅_J; Vec_j=1^J-1(𝐘̅_j-𝐘̅_J) ]),

and

    =([ K/K(v+u)+1 Vec_j=1^J-1Vec_i=N_J-1+1^N_j(([     u; v 𝐞_j ])(𝐘̅_j,i-𝐘̅_j));                K/Ku+1 Vec_i=N_J-1+1^N_j(([     u; 0_J-1 ])(𝐘̅_J,i-𝐘̅_J)) ]).

Now the results (<ref>)-(<ref>) of Theorems <ref> and <ref> are straightforward to verify.


















To proof Theorems <ref> and <ref> we firstly compute the blocks 𝐂_11, 𝐂_12 and 𝐂_22 of the mean squared error matrix (<ref>):





    𝐂_11=σ^2(Ku+1)/Km( [                                     1                              -1_J-1^⊤;                                -1_J-1 (K(u+v)+1)m/(Ku+1)n𝕀_J-1+1_J-11_J-1^⊤ ])


    𝐂_22=σ^2([ 𝐂_221     0;     0 𝐂_222 ]),

where

    𝐂_221   =   𝕀_n(J-1)⊗block-diag(u, 𝕀_J-1)
          - K/K(v+u)+1Diag_j=1^J-1((𝕀_n-1/n1_n1_n^⊤)⊗(([     u; v 𝐞_j ])([     u; v 𝐞_j ])^⊤) ),


    𝐂_222=𝕀_m⊗block-diag(u, 𝕀_J-1)-k/Ku+1(𝕀_m-1/m1_m1_m^⊤)⊗(([     u; 0_J-1 ])([     u; 0_J-1 ])^⊤),

and

    𝐂_12=-σ^2(𝐂_121  ⋮  𝐂_122),

where

    𝐂_121= tVec_j=1^J-1(1/n1_n^⊤⊗(([   0; 𝐞_j ])([    u; v𝐞_j ])^⊤)),


    𝐂_122=1/m1_m^⊤⊗(([      1; -1_J-1 ])([     u; 0_J-1 ])^⊤).


We can observe that _0=(0_J-1 ⋮ 𝕀_J-1)_0 and _0=(0_J-1 ⋮ 𝕀_J-1)_0 is the BLUE of _0. Then the covariance matrix of _0 can be determined using the formula

    Cov(_0)=(0_J-1 ⋮ 𝕀_J-1)𝐂_11(0_J-1 ⋮ 𝕀_J-1)^⊤,

which implies result (<ref>).





For the vector  of all individual treatment effects it cam be verified that





    =(1_N⊗(0_J-1 ⋮ 𝕀_J-1))_0+(𝕀_N⊗(0_J-1 ⋮ 𝕀_J-1))

and the BLUP of  is given by 

    =(1_N⊗(0_J-1 ⋮ 𝕀_J-1))_0+(𝕀_N⊗(0_J-1 ⋮ 𝕀_J-1)).

Then the mean squared error matrix of  is of general form (<ref>) with

    𝐁_1=(1_N⊗(0_J-1 ⋮ 𝕀_J-1))𝐂_11(1_N⊗(0_J-1 ⋮ 𝕀_J-1))^⊤,


    𝐁_2=(1_N⊗(0_J-1 ⋮ 𝕀_J-1))𝐂_12(𝕀_N⊗(0_J-1 ⋮ 𝕀_J-1))^⊤

and

    𝐁_3=(𝕀_N⊗(0_J-1 ⋮ 𝕀_J-1))𝐂_22(𝕀_N⊗(0_J-1 ⋮ 𝕀_J-1))^⊤.

After applying (<ref>)-(<ref>) we obtain the result of Theorem <ref>. 































 §.§ Proof of Lemma <ref>

To determine the eigenvalues of Cov(-), we have to solve the equation

    det( Cov(-)-λ 𝕀_N)=0.

From (<ref>) it follows that

    Cov(-)-λ 𝕀_N =:([  𝐇̃_11   𝐇_12; 𝐇_12^⊤  𝐇̃_22 ]),

where

    𝐇̃_11=(a_1-λ) 1/n1_n1_n^⊤+(a_2-λ)(𝕀_n-1/n1_n1_n^⊤)

for a_1=σ^2 N(Ku+1)/K m and a_2=σ^2 v (Ku+1)/K(v+u)+1,

    𝐇̃_22=a_3 1/m1_m1_m^⊤+(σ^2v-λ)𝕀_m

for a_3=σ^2(K(v+u)+1)m/K n+a_1m/N,
and 𝐇_12 is the same as in (<ref>).

Then we compute the determinant of Cov(-)-λ 𝕀_N as

    det( Cov(-)-λ 𝕀_N)=det(𝐇̃_11)det(𝐇̃_22-𝐇_12^⊤ 𝐇̃_11^-1𝐇_12),

where

    det(𝐇̃_11)=(a_1-λ)(a_2-λ)^n-1,


    𝐇̃_11^-1=1/a_1-λ 1/n1_n1_n^⊤+1/a_2-λ(𝕀_n-1/n1_n1_n^⊤),


    𝐇̃_22-𝐇_12^⊤ 𝐇̃_11^-1𝐇_12=(a_3-a_1^2 m/(a_1-λ ) n)1/m1_m1_m^⊤+(σ^2v-λ)𝕀_m,


    det(𝐇̃_22-𝐇_12^⊤ 𝐇̃_11^-1𝐇_12)=(σ^2v-λ)^m-1(a_3-a_1^2 m/(a_1-λ ) n+σ^2v-λ).

Then we obtain

    det( Cov(-)-λ 𝕀_N)=(a_2-λ)^n-1(σ^2v-λ)^m-1((a_3+σ^2v-λ)(a_1-λ)-a_1^2 m/n),

which results in the following solutions of equation (<ref>):

    λ_1=σ^2 v (Ku+1)/K(v+u)+1,


    λ_2=σ^2v,


    λ_3=σ^2N/2 Kn m(Km v+N(Ku+1)+√(s_n,m)),

where s_n,m=K^2m^2v^2+2Km(m-n)(Ku+1)v+N^2(Ku+1)^2,
and

    λ_4=σ^2N/2 Kn m(Km v+N(Ku+1)-√(s_n,m)).





After applying n=N w and m=N(1-w) we can see that s_n,m=N^2s_w and obtain the results of the lemma.

natbib













