Optimal Design in Hierarchical Models with application in Multi-center Trials
Maryna Prus Norbert Benda Rainer Schwabe





Optimal Design in Hierarchical Models with application in Multi-center Trials
    Maryna Prus, Norbert Benda, Rainer Schwabe
    Received 22 May 2018 / Accepted 17 July 2018
=============================================================================





Hierarchical random effect models are used for different purposes in clinical research and other areas. In general, the main focus is on population parameters related to the expected treatment effects or group differences among all units of an upper level (e.g. subjects in many settings). Optimal design for estimation of population parameters are well established for many models.  
However, optimal designs for the prediction for the individual units may be different. Several settings are identified in which individual prediction may be of interest. In this paper we determine optimal designs for the individual predictions, e.g. in multi-center trials, and compare them to a conventional balanced design with respect to treatment allocation.
Our investigations show, that balanced designs are far from optimal if the treatment effects vary strongly as compared to the residual error and more subjects should be recruited to the active (new) treatment in multi-center trials. 
Nevertheless, efficiency loss may be limited resulting in a moderate sample size increase when individual predictions are foreseen with a balanced allocation.





§ INTRODUCTION

Hierarchical random effect models are used for different purposes. Common applications, e.g. in clinical research, are random effect meta analyses of several clinical trials using individual patient data (IPD), multi-center trials assuming random center effects or mixed or random effect models for repeated measurements (MMRM) for longitudinal data within subjects. Clinical trials or meta-analyses of clinical trials may involve several populations or subgroups of patients (e.g. defined by genetic biomarkers), where in some settings between-population variability may be modeled using a hierarchical approach. Apart from clinical trials or observational studies data on quality parameters of drugs arise from hierarchical settings with multiple layers, i.e. with multiple sources of variability according to the underlying manufacturing process, e.g. given by manufacturing sites, production batches and samples within batches. In general, the main focus is on population parameters related to the expected treatment effects or group differences among all units of an upper level (e.g. trials in IPD meta-analyses, centers in multi-center trials, patients in longitudinal trials, batches in quality control, etc.). Several authors considered optimal design for estimation of population parameters (expected values of random effects) in similar models, see, e.g. <cit.>, <cit.>, <cit.> and <cit.>.

However, prediction of the outcome in the individual units may also be of interest in several settings, as for treatment effects in single centers to assess qualification of individual clinics, manufacturing sites in manufacturing control, treatment effects in different subpopulations of patients. In these cases, the question arises, whether optimal designs for populations parameters can also be used for individual predictions of random effects, whether optimal designs for individual predictions differ from those for populations parameters and to which extent and which efficiency loss (or sample size increase) can be anticipated if another conventional design is chosen.

Therefore, we investigated optimal designs for random effects to be applied, e.g. in multi-center trials and compared them to a conventional balanced design with respect to treatment allocation.

The structure of the paper is as follows: 
In the second section a model of a multi-center trial is specified and the best linear unbiased predictions of the individual center parameters (intercepts and treatment effects) are derived. 
In Section 3 analytical results are presented to characterize optimal designs for prediction in multi-center trials. 
The results are illustrated by some numerical examples. 
The paper is concluded by a short discussion.



§ MODEL SPECIFICATION

We consider a multicenter trial with K different centers. 
In all of these centers individuals are allocated to two treatment groups. 
In the first group (denoted by x=1) the individuals receive an active treatment while in the second group (denoted by x=0) a placebo or a control treatment is applied. 

Denote by μ_i and α_i the intercept (mean response at placebo or control) and the effect of the active treatment (compared to placebo or control), respectively, in center i=1,...,K, which both may vary across the centers.
The response Y_ij of an individual j=1,...,N_i in center i can be described as

    Y_ij=μ_i+α_ix_ij + ε_ij ,

where x_ij is equal to 1, if the individual belongs to the treatment group, and x_ij is equal to 0 for the control (or placebo) group and ε_ij denotes the random variation in the response of the individuals.
The individual variations ε_ij are assumed ot have zero mean and to be homoscedastic with common variance σ^2.

The centers are assumed to be similar and, hence, representatives of a larger entity of centers. 
Thus the center specific intercepts μ_i and treatment effects α_i can be assumed as random with (unknown) expected values E(μ_i)=μ and E(α_i)=α characterizing the mean intercept and mean treatment effect across the centers and covariance structure Cov((μ_i, α_i)^⊤)=σ^2𝐃 for some 2× 2 positive definite dispersion matrix 𝐃.
All random effects and all individual variations are assumed to be uncorrelated.

For the sake of simplicity we further assume that the total number N of individuals is the same for all centers (N_i=N) and that the allocation rate is constant across the centers, i. e.  the number n of individuals in the treatment group is the same for all centers. 
The design problem can then be formulated in terms of finding the optimal allocation rate w=n/N for the treatment group.
 
Because of exchangeability of the individuals within each center we may sort them in the analysis regardless of randomization in such a way that the first n individuals j=1,...,n to be analyzed receive the active treatment and the remaining N-n individuals j=n+1,...,N are in the control group.
Then the experimental settings x_ij in (<ref>) can be specified by

    x_ij=x_j={[           1,   j=1,..., n;           0, j=n+1,..., N ].

and are independent of the center i.

Hence, the multi-center model (<ref>) can be identified as a particular case of the random coefficient regression model 

    Y_ij=𝐟(x_j)^⊤_i+ε_ij,   i=1, …, K,    j=1, …, N

investigated by <cit.>
when the regression functions and the center parameters are specified by 𝐟(x)=(1, x)^⊤ and
_i=(μ_i,α_i)^⊤, respectively.
In general this model can be written in vector notation as

    𝐘_i=𝐅_i+ε_i,

where 𝐘_i=(Y_i1,...,Y_iN)^⊤ and  ε_i=(ε_i1,...,ε_iN)^⊤ are the N-dimensional vectors of observations and individual variations at center i, respectively, and 𝐅=(𝐟(x_1), ..., 𝐟(x_N))^⊤ is the within center design matrix which is equal across all centers. 

For the present multicenter model the design matrix 𝐅 simplifies to

    𝐅=([    1_n    1_n;  1_N-n 0_ N-n ]),

where 1_ℓ and 0_ ℓ denote the ℓ-dimensional vectors with all entries equal to 1 and 0, respectively. 

According to <cit.> (see also <cit.>) in model (<ref>) the best linear unbiased predictors (BLUP) 

    _i=(𝐅^⊤𝐅+𝐃^-1)^-1(𝐅^⊤𝐅 _i; ind+𝐃^-1_0)

of the random parameters _i are weighted combinations of the individual estimates  _i;ind=(𝐅^⊤𝐅)^-1𝐅^⊤𝐘_i based only on the observations in  center i and the the best linear unbiased estimator (BLUE)  _0=(𝐅^⊤𝐅)^-1𝐅^⊤𝐘̅ of the population parameter _0=E(_i), where 𝐘̅=1/K∑_i=1^K𝐘_i is the mean observational vector averaged across the centers..

We additionally assume that the center intercepts μ_i and the center treatment effects α_i are uncorrelated for all centers, i. e. 𝐃=diag(u,v), where u=σ_μ^2/σ^2>0 and v=σ_α^2/σ^2>0 are the variance ratios of the intercepts and the treatment effects in relation to the observational variance of the individuals. 

With the standard notations Y_i ·^(T)=1/n∑_j=1^nY_ij and Y_i ·^(C)=1/N-n∑_j=n+1^NY_ij for the mean response in the treatment (“T”) and the control (“C”) groups in center i, Y_· ·^(T)=1/K∑_i=1^KY_i ·^(T) and Y_· ·^(C)=1/K∑_i=1^KY_i ·^(C) for the overall mean of the treatment and the control groups, respectively, the BLUPs for the center parameters _i=(μ̂_i,α̂_i)^⊤ of the random intercepts and the random treatment effects _i=(μ_i,α_i)^⊤ in model (<ref>) can be written as weighted averages

    μ̂_i=c_0(Y_i ·^(T)-Y_· ·^(T))+cY_i ·^(C)+(1-c)Y_· ·^(C),

where c_0=nu/(Nu+1)(nv+1)-n^2uv and c=u(N-n)(nv+1)/(Nu+1)(nv+1)-n^2uv,
and

    α̂_i=c_1Y_i ·^(T)+(1-c_1)Y_· ·^(T)-c_2Y_i ·^(C)-(1-c_2)Y_· ·^(C)

with weights c_1=(Nu+1)nv-n^2uv/(Nu+1)(nv+1)-n^2uv and c_2=Nu+nv+1/(Nu+1)(nv+1)-n^2uv.

To measure the quality of a design we will use the mean squared error (MSE) matrix of the BLUP =(_1^⊤,..., _K^⊤)^⊤ of the complete vector =(_1^⊤,..., _K^⊤)^⊤ of all random parameters. 
The mean squared error (MSE) matrix for the BLUP  can be computed by means of the following formula :

    Cov(-)
       =   σ^2(1/K(1_K1_K^⊤ )⊗  (𝐅^⊤𝐅)^-1.
           + .(I_K-1/K1_K1_K^⊤ )⊗ (𝐅^⊤𝐅+𝐃^-1)^-1)

 (see <cit.>), where I_ℓ denotes the ℓ×ℓ identity matrix and ⊗ is the symbol for the Kronecker product of matrices or vectors.

Further denote by _α=(_1,..., _K)^⊤ the vector of treatment effects for all centers. 
Then _α= (𝕀_K⊗(0,  1) ) and, hence, 

    Cov(_α-_α)=(𝕀_K⊗(0,  1) )Cov(-)(𝕀_K⊗(0,  1) )^⊤
 
for the MSE matrix of the BLUP _α=(α̂_1,..., α̂_K)^⊤ of _α.
Using this and formula (<ref>) we obtain the MSE matrix 

    Cov(_α-_α)   =   σ^2(N/n(N-n)1/K1_K1_K^⊤.
          +.v(Nu+1)/(Nu+1)(nv+1)-n^2uv(𝕀_K-1/K1_K1_K^⊤)).






§ OPTIMAL DESIGN


As individuals are interchangeable within treatment groups we may define an exact within center design

    ([   T   C;   n N-n ])

by the allocation numbers n and N-n to the treatment and control group T and C, respectively. 

For analytical purposes, we generalize this to the definition of an approximate design:

    ξ = ([   T   C;   w 1-w ]),

where w=n/N is the allocation rate to the treatment group and 1-w=N-n/N is the allocation rate to the control group. 
For finding an optimal design only the optimal allocation rate w^* to the treatment group has to be determined.

For an approximate design the definition of the MSE matrix (<ref>) of the BLUP _α is extended in a straightforward manner an can be rewritten (neglecting σ^2) as :

    MSE(w)   =   1/Nw(1-w)1/K1_K1_K^⊤
          +v(Nu+1)/(Nu+1)(Nwv+1)-N^2w^2uv(𝕀_K-1/K1_K1_K^⊤)

 in terms of the allocation rate w.
 The approach of approximate designs seems in so far to be appropriate as the total number N of individuals in each center should be sufficiently large.
 Otherwise optimal exact designs have to be obtained by adequate rounding of w^* to a multiple of 1/N (see below).

For the assessment of the MSE matrix we focus on the A-optimality criterion which averages the mean squared errors of the center treatment effects.
More specifically, the A-criterion Φ_α is the trace of the MSE matrix of the prediction _α of the center treatment effects. 
For an approximate design we get 

    Φ_α(w)=1/Nw(1-w)+(K-1)v(Nu+1)/(Nu+1)(Nwv+1)-N^2w^2uv

for the criterion function Φ_α in terms of the allocation rate w.
Because, in general, there is no explicit solution for the optimal allocation rate, which minimizes (<ref>), we will give an insight in the qualitative behavior by some numerical example below.

It is worth-while mentioning that the criterion (<ref>) is convex and, therefore, an optimal exact design may be obtained by choosing the best of two exact designs adjacent to an optimal approximate design.



















 §.§ Example.


For illustrative purposes we consider a numerical example with K=50 centers and N=10 individuals in each center. 
Figure 1 exhibits the behavior of the optimal allocation rate w^* to the treatment group in dependence of the variance ratio v of the treatment effects for some fixed values 0.01, 0.1, 0.25, 0.5 and 1,5 of the variance ratio u of the intercept. 
For reasons of presentation we plot the optimal allocation rate w^* against the rescaled variance ratio r_v=v/(1+v) in the spirit of intra-class correlation in order to cover all possible values of the treatment effects variance by a finite interval ((0,1)). 
Each value of the variance ratio u of the intercepts is represented by one solid line.
What can be seen from the picture is that for fixed values of u  the optimal allocation rate w^* is equal to 0.5 for v→ 0 and increases with increasing values of the variance ratio v of the treatment effects. 
The different lines associated to the different values of u appear in descending order which means that the optimal allocation rates decrease when the variance ratio of the intercepts gets larger. 

The next figure (Figure 2) shows the behavior of the optimal allocation rate in dependence of the variance ratio u of the intercepts for fixed values 0.01, 0.1, 0.2, 0.5 and 2 of the variance ratio v of the treatment effects, where again the variance ratio is rescaled (r_u=u/(u+1)). 
Also here it can be seen that the optimal allocation rate decreases with increasing values of u and increases with increasing values of v.
       
		
Finally, Figures 3 and 4 present the efficiency of the equal allocation rate w_0=0.5 which is optimal in the fixed effects model (u=v=0). The efficiency for the A-criterion  (A-efficiency) has been computed using the standard formula

    eff(w_0)=Φ_α(w^*)/Φ_α(w_0).

		    
		The efficiencies decrease with increasing values of the variance of the treatment effects and increase with increasing values of the variance of the intercepts.
		
		Example cont.

For illustrative purposes we consider numerical examples with K=50 or K=100 centers and N=10 or N=5 individuals in each center. 
Figures 5 and 6 exhibit the behavior of the optimal allocation rate w^* to the treatment group in dependence of the variance ratio q=v/u for some fixed values 0.01 (solid line), 0.1 (dashed line), 0.25 (dotted line), 0.5 (dashed-dotted line) and 1,5 (dashed line, long dashes) of the variance ratio u of the intercept. 
For reasons of presentation we plot the optimal allocation rate w^* against the rescaled variance ratio r_q=q/(1+q) in the spirit of intraclass correlation in order to cover all possible values of variance ratio by a finite interval ((0,1)). As we can observe on the graphics, the dependence on the variance ratio is more essential for K=100 and N=5 than for K=50 and N=10.
       
		
Figures 7 and 8 presents the efficiency of the equal allocation rate w_0=0.5 which is optimal in the fixed effects model. 
The efficiency is more sensible with respect to the variance ratio in the case K=100 and N=5 than for K=50 and N=10.
		    
		


§ DISCUSSION

As illustrated in the example, the larger the between-unit (between-trial) variability of the treatment effects the more the optimal weight deviates from equal allocation, especially if the variance of the units' intercepts is small. An increasing heterogeneity in the treatment effect leads to a decreased precision of the design that is optimal for population parameters: A balanced design is far from optimal if the treatment effects vary strongly as compared to the residual error and more subjects should be recruited to the active (new) treatment in multi-center trials. Nevertheless, it appears reassuring to the clinical trial practitioner that the efficiency loss may be limited as in the example resulting in a total sample size increase of about 10 - 20 % in the considered scenarios if individual predictions are foreseen with a balanced allocation. Usually, between-unit variability of treatment effects are considered to be rather small, indicating that equal allocation may suffice. However, using the results given in the paper, specific settings with different expectations can be assessed properly, in order to make optimally use of a limited number of patients or sample units to predict random effects of individual units.  
		

natbib





