Design optimisation and post-trial analysis in group sequential stepped-wedge cluster randomised trials
Michael Grayling David Robertson James Wason Adrian Mander





Design optimisation and post-trial analysis in group sequential stepped-wedge cluster randomised trials
    M. J. Grayling1, D. S. Robertson1, J. M. S. Wason1,2, A. P. Mander1

1. Hub for Trials Methodology Research, MRC Biostatistics Unit, Cambridge, UK, 
 2. Institute of Health and Society, Newcastle University, Newcastle, UK.
    
==================================================================================================================================================================================================================================



Corresponding Author: M. J. Grayling, MRC Biostatistics Unit, Forvie Site, Robinson Way, Cambridge CB2 0SR, UK; Tel: +44-(0)1223-330300; E-mail: mjg211@cam.ac.uk.


Abstract

Recently, methodology was presented to facilitate the incorporation of interim analyses in stepped-wedge (SW) cluster randomised trials (CRTs). Here, we extend this previous discussion. We detail how the stopping boundaries, allocation sequences, and per-cluster per-period sample size of a group sequential SW-CRT can be optimised. We then describe methods by which point estimates, p-values, and confidence intervals, which account for the sequential nature of the design, can be calculated. We demonstrate that optimal sequential designs can reduce the expected required number of measurements under the null hypothesis, compared to the classical design, by up to 30%, with no cost to the maximal possible required number of measurements. Furthermore, the adjusted analysis procedure almost universally reduces the average bias in the point estimate, and consistently provides a confidence interval with coverage close to the nominal level. In contrast, the coverage of a naive 95% confidence interval is observed to range between 92 and 98%. Methodology is now readily available for the efficient design and analysis of group sequential SW-CRTs. In scenarios in which there are substantial ethical or financial reasons to terminate a SW-CRT as soon as possible, trialists should strongly consider a group sequential approach.



Keywords: Cluster randomized trial; Group sequential; Interim analysis; Stepped-wedge.





§ INTRODUCTION



In a cluster randomised trial (CRT), groups of subjects rather than individuals are randomised to treatment conditions. Today, numerous reasons are recognised for why one may employ a CRT design. Primarily, CRTs aid in the control of contamination between subjects, where participants from the same group would likely influence each other if randomised to different treatments. Furthermore, they can bring increased administrative efficiency, helping to overcome the barriers of contacting and recruiting large numbers of participants, whilst they additionally allow the study of interventions that could not be directed toward a single individual. Unfortunately though, there are also several noted disadvantages to these designs. In particular, since participants from a cluster are likely to be more similar than those from different clusters, there is a potential for correlation between observations from within each cluster. This must be accounted for not only in the analysis, but also at the design stage, as standard sample size calculation formulae are no longer applicable. See, for example, Eldridge and Kerry (2012) <cit.>, Campbell and Walters (2014) <cit.>, or Hayes and Moulton (2017) <cit.> for further discussions of these issues. 

One sub-category of CRT design that has been the subject of substantial research in recent years is the stepped-wedge (SW) CRT. The key feature of a SW-CRT is that an experimental intervention is introduced across the clusters over a number of time periods using one-directional switching only. That is, once a cluster switches to the experimental intervention in a SW-CRT, it does not revert back to the control intervention for the remainder of the trial. Today, methodology is well established for sample size calculation <cit.>, and for data analysis <cit.>, in fixed sample SW-CRTs.

In addition, techniques for applying classical group sequential trial design methodology to SW-CRTs with interim analyses were recently described <cit.>. This work followed previous suggestions that early stopping for efficacy may be important in SW-CRTs given the design is regularly associated with scenarios in which the experimental intervention is expected to do more good than harm <cit.>. Whilst a review of the outcome of completed SW-CRTs led to the proposal that early stopping for futility may be helpful to guard against the over-enthusiastic use of the design <cit.>. Ultimately, in this work, it was demonstrated that the incorporation of interim analyses in to a SW-CRT could bring substantial efficiency gains in terms of the expected number of required measurements, compared to utilising a classical fixed sample design.

However, one factor left unexamined was how the design of a group sequential SW-CRT could be optimised. Optimal fixed sample SW-CRT design has been a topic of much recent interest <cit.>, with the same also true of the optimisation of stopping boundaries in conventional individually randomised group sequential trials <cit.>. For both fixed sample SW-CRTs and group sequential trials, design optimisation was demonstrated to provide substantial efficiency gains. We may therefore expect that it would be of great value to establish methodology for optimising the design of group sequential SW-CRTs.

In addition, methodology for post-trial analysis in group sequential SW-CRTs was not considered by Grayling et al. (2017a) <cit.>. It is well known that if a final analysis is performed on data acquired within a sequential trial, that does not take in to account its sequential nature, then biased point estimates, and p-values and confidence intervals with incorrect coverage will be routinely acquired <cit.>. Accordingly, the lack of available guidance on how to analyse data accrued in a group sequential SW-CRT currently limits the application of such designs in practice.

In this work, we address these problems. Specifically, we first detail how the stopping boundaries, allocation sequences, and per-cluster per-period sample size of a group sequential SW-CRT can be optimised. Following this, we describe an approach for acquiring a median unbiased estimator for the treatment effect, along with methodology for constructing adjusted p-values and confidence intervals. Finally, for two example trial design scenarios, we explore the efficiency gains optimal designs could bring, and examine the performance of the proposed procedures for data analysis.

We now proceed by describing the notation utilised in the rest of the paper. The aforementioned methodology is then presented along with the results of a simulation study. We conclude with a short discussion. 



§ METHODS
 



 §.§ Notation, hypotheses and analysis in group sequential SW-CRTs


Define ℕ_a={k∈ℕ : k≤ a} for any a∈ℕ, and ℕ_a^+=ℕ_a\{0}. Then, we consider a SW-CRT conducted in C∈ℕ\ℕ_1 clusters over T∈ℕ\ℕ_1 time periods, with m∈ℕ\ℕ_1 measurements taken per-cluster per-period, to compare a single experimental intervention to some control. We do not restrict ourselves to cross-sectional designs (in which the measurements are assumed to be taken on different patients in each time period) or cohort designs (in which the measurements are assumed to be taken on the same patients in each time period). However, for simplicity, we focus only on complete block SW-CRTs, where measurements are accrued on every cluster in every time period. Given this, we denote by X, dim(X)=C× T, the matrix of binary treatment indicators X_c,t, where X_c,t=1 if cluster c receives the intervention in time period t, and X_c,t=0 otherwise.

Denoting the vector of responses available after time period t∈ℕ_T^+ by Y_t|m,X, we assume that Y_t|m,X∼ N(D_t|m,Xβ,Σ_t|m,X), with Σ_T|m,X non-singular. Here the subscript t∈ℕ_T^+ indicates we consider the subsections of D_T|m,X and Σ_T|m,X pertaining to the first t time periods only, where D_T|m,X and Σ_T|m,X are designated design and covariance matrices for the vector of responses Y_T|m,X. Furthermore, β=(β_1,…,β_p)^⊤ is a specified vector of fixed effects. Note that we have explicitly stated the distribution conditional on the values of m and X, since these will later be treated as variables to be determined.

We assume that we would like to make inference about the final element of β, β_p, which we denote for simplicity by τ. Formally, we test the one-sided hypothesis H_0 : τ≤0, against H_1 : τ>0.

A group sequential SW-CRT can then be completely defined given

	
  * A set of integers 𝒯={t_1,…,t_|𝒯|}, with t_i∈ℕ_T^+ for i∈ℕ_|𝒯|^+, t_i>t_i-1 for i∈ℕ_|𝒯|\ℕ_1, and t_|𝒯| =T. 
	
  * Futility and efficacy boundaries f=(f_1,…,f_|𝒯|)^⊤ and e=(e_1,…,e_|𝒯|)^⊤, with (f,e)∈ℬ_|𝒯|={(a,b)∈ℝ^|𝒯|×ℝ^|𝒯| : ∀ i∈ℕ_|𝒯|-1^+ a_i<b_i, a_|𝒯|=b_|𝒯|}.


With this, |𝒯| analyses are carried out in total, with analysis i∈ℕ_|𝒯|^+ conducted after time period t_i∈ℕ_T^+. Explicitly, for each i∈ℕ_|𝒯|^+ we compute an estimate for β, β̂_i|m,X, using generalised least squares

    β̂_i|m,X   = (β̂_1,i|m,X,…,β̂_p,i|m,X)^⊤,
       = (D_t_i|m,X^⊤Σ_t_i|m,X^-1D_t_i|m,X)^-1D_t_i|m,X^⊤Σ_t_i|m,X^-1Y_t_i|m,X.

From this, we use the estimate of τ,  τ̂_i|m,X=β̂_p,i|m,X, to form the following standardised test statistic

    Z_i|m,X=τ̂_i|m,X/{Var(τ̂_i|m,X)}^1/2 = τ̂_i|m,X{(D_t_i|m,X^⊤Σ_t_i|m,X^-1D_t_i|m,X)_[p,p]^-1}^-1/2 = τ̂_i|m,X I_i|m,X^1/2.


Note that for this to be possible (i.e., for τ to be estimable at every interim analysis), it is important that at least one cluster receives the intervention in one of the time periods t∈ℕ_t_1^+. This can be verified using the matrix X by checking ∑_c∈ℕ_C^+∑_t∈ℕ_t_1^+X_c,t>0.

An algorithm for conducting a group sequential SW-CRT can then be defined as follows


	
  * Set t=1.
	
  * Conduct period t of the trial, allocating treatments according to the values X_c,t for c∈ℕ_C^+, and accruing m measurements in each cluster.
	
  * If t∉𝒯 set t=t+1 and return to 2. Otherwise, if t∈𝒯, t=t_i say, compute Z_i|m,X. Then
	
		
  * If Z_i|m,X≤ f_i stop the trial, and do not reject H_0.
		
  * If Z_i|m,X> e_i stop the trial, and reject H_0.
		
  * If f_i < Z_i|m,X≤ e_i, set t=t+1 and return to 2.
	


Here, the trial is guaranteed to terminate, with a decision also made on H_0, by the end of time period T because of the stipulation that f_|𝒯|=e_|𝒯|.

Importantly, (Z_1|m,X,…,Z_|𝒯||m,X)^⊤ is multivariate normal with

    𝔼(Z_i|m,X)    = τ I_i|m,X^1/2,      i∈ℕ_|𝒯|^+,
    Cov(Z_i|m,X,Z_j|m,X)    = (I_i|m,X/I_j|m,X)^1/2,      i,j∈ℕ_|𝒯|^+, i≤ j,

and so multivariate normal integration can be used to evaluate the statistical operating characteristics of the design <cit.>.

To see this clearly, define the random variable Γ∈ℕ_|𝒯|^+ as the analysis at which the trial is terminated, and Ψ∈ℕ_1 as the random variable with Ψ=1 if H_0 is rejected, and Ψ=0 otherwise. The possible outcomes to the trial are then categorised by the pair {Ψ,Γ}, and for any {ψ,γ}∈ℕ_1×ℕ_|𝒯|^+

	
    ℙ(Γ=γ,Ψ=ψ | τ,f,e,m,X) = ∫_l(1,γ,ψ,f,e)^u(1,γ,ψ,f,e)…∫_l(γ,γ,ψ,f,e)^u(γ,γ,ψ,f,e)ϕ(y,τI_γ|m,X^1/2,Λ_γ|m,X)dy_γ…dy_1,


where 

	
  * ϕ{y,μ,Λ} is the probability density function of a multivariate normal distribution with mean μ and covariance matrix Λ, evaluated at vector y,
	
  * I_γ|m,X^1/2 = (I_1|m,X^1/2,…,I_γ|m,X^1/2)^⊤.
	
  * l and u are functions that provide the lower and upper integration limits for each integral. Explicitly, for any (f,e)∈ℬ_|𝒯|
	
    l(i,γ,ψ,f,e) =  f_i    : i<γ, 
     e_i    : i=γ, ψ=1, 
     -∞   : i=γ, ψ=0,     u(i,γ,ψ,f,e) =  e_i    : i<γ, 
    ∞   : i=γ, ψ=1, 
     f_i     : i=γ, ψ=0.

	
  * Λ_γ|m,X is the covariance matrix of the standardized test statistics across the interim analyses 1,…,γ
	
    Λ_γ|m,X =
    	[ Cov(Z_1|m,X,Z_1|m,X)                    … Cov(Z_1|m,X,Z_γ|m,X);                    ⋮                    ⋱                    ⋮; Cov(Z_γ|m,X,Z_1|m,X)                    … Cov(Z_γ|m,X,Z_γ|m,X) ].


	
With this, the probability that H_0 is rejected, and the expected number of required measurements (ENM), can be computed as


    P(τ|f,e,m,X)    = ∑_γ∈ℕ_|𝒯|^+ℙ(Γ=γ,Ψ=1 | τ,f_γ,e_γ,m,X),
    
    ENM(τ|f,e,m,X)    = ∑_γ∈ℕ_|𝒯|^+∑_ψ∈ℕ_1mCt_γℙ(Γ=γ,Ψ=ψ | τ,f_γ,e_γ,m,X).




 §.§ Optimal group sequential designs


Section <ref> describes a method for conducting a group sequential SW-CRT given choices for several design parameters. However, at the design stage, we require a method for choosing values for these parameters to provide desired operating characteristics. One simple approach to this, after designating a value for X, is the error spending method, as discussed in Grayling et al. (2017a) <cit.>. Here, we consider instead how X, f, e, and either C or m, can be optimised. 

For this, much depends on whether the choice is to fix C or m. In the former case, the matrix X can be readily optimised, as will be described shortly. In the latter case, optimisation of X is complex since its row dimension is variable in the search for the required value of C. In this case, X should be specified implicitly through rules on what proportion of clusters should switch to the intervention in each time period, reducing the problem to optimising only f, e, and C. Then, having determined an appropriate value for C, one can employ the approach for fixed C to optimise X if desired.

Thus, from here, we focus on the problem of simultaneously optimising X, f, e, and m, when C has been fixed. To this end, define the vector S=(S_1,…,S_C)^⊤∈𝒮_𝒯 of switching times, where

    𝒮_𝒯={s∈{ℕ_T+1^+}^C : min_c∈ℕ_C^+ s_c≤ t_1, ∀ t∈ℕ_T^+ ∏_c∈ℕ_C^+𝕀(s_c=t)=0 }.

That is, X_c,t=𝕀(t ≥ S_c) for (c,t)∈ℕ_C^+×ℕ_T^+. Here, we have used S_c=T+1 to indicate that cluster c does not receive the intervention during the trial, whilst the additional constraints ensure τ is estimable at every interim analysis, and that at least two distinct sequences of allocations are used across the clusters. We will use the notation X(S) to indicate the matrix X implied by a particular choice for S.

We then consider an optimal design framework where we wish to identify the solution to a problem which can be cast in the following form


    _{(f,e),m,S}∈ℬ_|𝒮|×(ℕ\ℕ_1)×𝒮_𝒯 O(f,e,m,S),

subject to

    P{0|f,e,m,X(S)}≤α,    
    P{δ|f,e,m,X(S)}≥1 - β.

Here, the constraints are interpreted as controlling the type-I error-rate to some α∈(0,1) when τ=0, and to have power of at least 1-β∈(0,1) when τ=δ>0.

In practice, the objective function O(·) could be specified to be as complex as desired. For example, as discussed in Grayling et al. (2017c) <cit.>, it could penalise designs which roll-out the intervention beyond a certain speed, or impose a cost on designs where clusters do not receive the experimental intervention. Here though, adapting previous relevant works on optimal group sequential trial design <cit.>, we consider

    O(f,e,m,S)=w_1 ENM{0|f,e,m,X(S)}+w_2 ENM{δ|f,e,m,X(S)}+w_3 mCT.

With this, the w_i∈ℝ^+∪{0} are weights indicating the relative strengths of the desires to minimise the ENM under the null, the ENM under the alternative hypothesis, and the maximal required number of observations. Later, we will use the notation w=(w_1,w_2,w_3).

Unfortunately, the complexity of the constraints on the trial's operating characteristics means they cannot easily be incorporated in to a conventional numerical optimisation routine that we may seek to use to solve the above problem. Therefore, following Wason and Jaki (2012) <cit.> and Wason et al. (2012) <cit.>, we translate our problem to be

	
    _{(f,e),m,S}∈ℬ_|𝒮|×(ℕ\ℕ_1)×𝒮_𝒯    O(f,e,m,S)+M_SW( 𝕀[P{0|f,e,m,X(S)}>α][P{0|f,e,m,X(S)}-α/α] .
               +. 𝕀[1-P{δ|f,e,m,X(S)}>β][1-P{δ|f,e,m,X(S)}-β/β] ).


Here M_SW is the number of observations required by a corresponding fixed-sample near-balanced SW-CRT (i.e., a design with as near-equal as possible numbers of clusters switching to the intervention in each time period). The final factor is thus a penalty for designs that do not meet the desired operating characteristics, to assist in converging towards a solution meeting these requirements.

The final challenge therefore is to identify the solution of the problem (<ref>). For this, a numerical optimisation routine capable of handling the fact that m and the S_c are discrete, but the components of f and e are continuous, is required. This algorithm must also be capable of handling the many local optima the search space will likely contain. One such routine, available in , is the  function from the package of the same name <cit.>. We utilise this for all examples given here, and in the Appendix we elaborate on our choices for the algorithms initialisation. 

Our full procedure for optimising the design of a group sequential SW-CRT has then been prescribed. In Section <ref> we present numerous such optimal designs.



 §.§ Post trial analysis: Point estimates, p-values, and confidence intervals


In this section, we turn our attention to the final analysis of data accrued in a completed group sequential SW-CRT. For notational simplicity, we therefore suppose that the trial was conducted using design 𝒟={𝒯,f, e,m,X} (i.e., these factors are all fixed and known in what follows). Additionally, we assume that the trial terminated with {Γ,Z_Γ|m,X}={γ,z_γ}∈Ξ_𝒟 = {{i,z_i}∈𝒯×ℝ : z_i∉(f_i,e_i]}.

With this, a simple method for specifying a point estimate for the treatment effect τ, is to use the maximum likelihood estimate (MLE) τ̂_N = τ̂_γ|m,X=β̂_p,γ|m,X,
as defined by Equation (<ref>). Additionally, one could designate a p-value as p_N=1 - Φ(z_γ), where Φ(z) is the cumulative density function of the standard normal distribution. Finally, a one-sided confidence interval for τ, could be taken as
C_N=(τ̂_LCI,N,∞) for τ̂_LCI,N=τ̂_N-Φ^-1(1-α)I_γ|m,X^-1/2.

The problem with these methods however is that they do not take into account the sequential nature of the design. As was noted in Section 1, such naive analysis procedures, as the subscript N is meant to indicate, will routinely result in biased point estimates, and p-values and confidence intervals with incorrect coverage. For example, Tsiatis et al. (1984) <cit.> demonstrated that the coverage of a naive 90% confidence interval at the end of a four-stage group sequential trial could be as low as 0.846.

Accordingly, methodology for performing an adjusted analysis is required. Unfortunately, there is no method that is optimal for analysis following a sequential trial in all situations. Whitehead (1986) <cit.> proposed one of the first methods for analysing data accrued in a group sequential trial, suggesting that the point estimate be taken as the solution to a simple equation involving the computation of the expected bias in the estimated treatment effect. This estimator was demonstrated to typically have extremely small bias. However, its evaluation is complex and it provides no easy means of nominating an associated confidence interval with the desired coverage. Emerson and Fleming (1990) <cit.> suggested use of a uniform minimum variance unbiased estimator, formed through Rao-Blackwellisation of the unbiased estimate of τ from the first stage data. Unfortunately, this estimator was found to have large variance.
 
Consequently, we instead present an approach that is both easy to apply, and generally has good performance. We use the so-called stage-wise ordering in a group sequential design to designate p-values and confidence intervals with exact coverage probabilities, and to specify a median unbiased point estimate for τ. These computations will require evaluation only of integrals of the form in Equation (<ref>).

The stage-wise ordering was first proposed by Armitage (1957) <cit.>, and has subsequently been utilised in many papers on individually randomised group sequential trials (see, for example, Siegmund (1978) <cit.> and Tsiatis et al. (1984) <cit.>). For any {i,z_i},{j,z_j}∈Ξ_𝒟, the stage-wise ordering considers {i,z_i} to be more extreme than {j,z_j} if i<j or if i=j and z_i≥ z_j. This ordering can then be used to derive our adjusted p-values and confidence intervals, as well as the median unbiased point estimate. Specifically, define

    E(τ| i,z_i,f,e,m,X)    = ∑_j∈ℕ_i-1^+ℙ(Γ=j,Ψ=1 | τ,f_j,e_j,m,X),
            +ℙ{Γ=i,Ψ=1 | τ,f_i,(e_i-1^⊤,z_i)^⊤,m,X},

which is the probability, for the given stage-wise ordering, of seeing a result more extreme than (i,z_i) for a particular value of τ.

Using this, an overall p-value for H_0 can be taken as p_SO = E(0|γ,z_γ,f,e,m,X). Similarly, the associated adjusted confidence interval is C_SO = (τ̂_LCI,SO,∞), where τ̂_LCI,SO is the solution of  α = E(τ̂_LCI,SO|γ,z_γ,f,e,m,X). Finally, a median unbiased estimate can be designated as the solution τ̂_SO to 0.5 = E(τ̂_SO|γ,z_γ,f,e,m,X).

Both τ̂_LCI,SO and τ̂_SO can be determined using a simple one-dimensional root finding algorithm, and thus we have established a potential means of performing an improved final analysis on data in group sequential SW-CRTs.

To examine the statistical performance of these analysis procedures, we simulate R=10^5 trials for several optimised group sequential designs 𝒟 and values for τ. The details of the data degeneration procedure in the simulation study are described in the Appendix. From this, for each considered design we compute, for r∈ℕ_R^+, C_x,r(τ|𝒟) and τ̂_x,r(τ|𝒟) for x∈{N,SO}. Here, τ̂_x,r(τ|𝒟) is, in a trial using design 𝒟, the point estimate at the final analysis for simulation r using estimator x, when the true treatment effect is τ, with an analogous definition holding for C_x,r(τ).

The performance of the two point estimation procedures is then assessed using the average bias and average root mean square error (RMSE), computed as

    b(τ| x,𝒟) = R^-1∑_r∈ℕ_R^+{τ̂_x,r(τ|𝒟)-τ},    
    RMSE(τ| x,𝒟) = √(R^-1∑_r∈ℕ_R^+{τ̂_x,r(τ|𝒟)-τ}^2).

Similarly, the confidence intervals are compared using their respective coverages, given by

    COV(θ| x,𝒟) = R^-1∑_r∈ℕ_R^+𝕀{τ∈ C_x,r(τ|𝒟)}.




 §.§ Examples


We follow Grayling et al. (2017a) <cit.> and consider two example trial design scenarios (TDSs). In both, we assume a cross-sectional design, with analysis conducted utilising the following linear mixed model, as proposed by Hussey and Hughes (2007) <cit.>

    y_ijk=μ+π_j+τ X_ij+c_i+ϵ_ijk,

for i∈ℕ_C^+, j∈ℕ_T^+, and k∈ℕ_m^+. Here

	
  * y_ijk is the response of the kth individual in the ith cluster and jth time period.
	
  * μ is an intercept term.
	
  * π_j is a fixed effect for the jth period (with π_1=0 for identifiability).
	
  * c_i ∼ N(0,σ_c^2) is a random effect for the ith cluster.
	
  * e_ijk∼ N(0,σ_e^2) is the residual error.


Now, for TDS1, we consider the SW-CRT reported by Bashour et al. (2013) <cit.>, on the effect of training doctors in communication skills on women's satisfaction with doctor-woman relationship during labour and delivery. In this trial, a balanced four cluster five period design was used, and the final analysis estimated that σ_c^2=0.02 and σ_e^2=0.51. Based on this, a sample size of 70 patients per-cluster per-period would have been required for the trial's desired operating characteristics (α=0.05, β=0.1, and δ=0.2), implying a total sample requirement of 1400 subjects. Thus, in TDS1 we fix C=4, T=5, σ_c^2=0.02, σ_e^2=0.51, α=0.05, β=0.1, δ=0.2, and M_SW=1400.

For TDS2, we base our design parameters on the average characteristics of completed SW-CRTs as reported in Grayling et al. (2017a) <cit.>. We set C=20 and T =9, to correspond to the median values used in-practice. Next, we suppose α=0.05, β=0.2,
and choose σ_e^2=1, σ_c^2=1/9, to provide a more moderate value for the intra-cluster correlation compared
to TDS1. In a near-balanced design (i.e., one where three clusters switch to the intervention in the second through fifth periods, and two clusters in each of the remaining time periods), a design with m=7 and δ=0.24 provides the desired power and type-I error-rate. This will therefore be our designated value of δ, and M_SW=1260.

Note that in these scenarios, D_T|m,X and Σ_T|m,X can be specified easily for any m and X. Moreover, recalling that C and T are designated for any X by the relationship dim(X)=C× T, we can then easily determine I_t|m,X using the result of Hussey and Hughes (2007) <cit.>

    I_t|m,X = (σ^2+tσ_c^2)(CU-W) + σ_c^2(U^2-CV)/Cσ^2(σ^2+tσ_c^2),

for σ^2=σ_e^2/m, and where

    U=∑_c∈ℕ_C^+∑_j∈ℕ_t^+X_c,j,     V=∑_c∈ℕ_C^+(∑_j∈ℕ_t^+X_c,j)^2,     W=∑_j∈ℕ_t^+(∑_c∈ℕ_C^+X_c,j)^2.


Code to replicate our determination of optimal designs, and to explore the performance of the considered analysis procedures, is available from https://github.com/mjg211/article_code.



§ RESULTS




 §.§ Optimal designs


For the TDSs described in Section <ref>, a numerical search was performed to identify the optimal values of e, f, m, and X for w∈{(1/3,1/3,1/3), (1/2,1/2,0), (0,1/2,1/2)}, with 𝒯={3,5} in TDS1 and 𝒯={3,6,9} in TDS2. The determined optimal designs, and a summary of their statistical operating characteristics, are presented in Table <ref>. Their type-I error-rates and power are omitted as each was verified to be the desired level.


	


In all instances the ENM under the null and alternative hypotheses is substantially reduced compared to the corresponding fixed sample design. Specifically, in TDS1 it is by as much as 30.1% under the null hypothesis, with w=(1/2,0,1/2), and 24.6% under the alternative hypothesis, with w=(0,1/2,1/2).

For TDS1, we can observe that the ability to alter X has allowed the three designs to incorporate interim analyses without the need to increase the value of m. We also see that the optimal choice for X is the same for each value of w. However, as would be expected, the ENMs required under the null and alternative hypotheses is strongly dependent on the choice of w. For example, with w=(1/2,0,1/2) the value of ENM(0) is lower, and the value of ENM(δ) higher, than with w=(1/3,1/3,1/3). A similar statement holds for the design with w=(0,1/2,1/2). This is a consequence of the identified optimal boundaries. For example, the design with w=(0,1/2,1/2) has f_1=-5.05, almost ensuring the trial will not stop early for futility, as no concern is given to stopping early under the null hypothesis.

Similar statements hold for TDS2. However, here, the optimal S varies for each considered w. Additionally, in this case, setting w=(0,1/2,1/2) results in a design with approximately the same required sample size under the alternative hypothesis as for w=(1/3,1/3,1/3). It appears that this value is approximately the lowest that can be attained for this TDS. However, these two designs have noticably different features outside of ENM(δ). As for TDS1, w=(0,1/2,1/2) results in a design with lower values for the early efficacy stopping boundaries, and almost removes the possibility to stop early for futility. This is a feature that should be anticipated in general.



 §.§ Post-trial analysis


For the optimal designs given in Table <ref>, a simulation study was conducted as described in Section <ref> to evaluate the performance of the naive and adjusted analysis procedures. In each case, we set μ=π_1=…=π_T=0 for simplicity. Figure <ref> displays the average bias and RMSE of the point estimates, as well as the coverage of the confidence intervals, in TDS1, when τ∈[-0.3,0.5]. Figure <ref> displays the equivalent findings for TDS2.

For both TDSs, the naive point estimator has, for certain values of τ, large bias. The particular shape of the bias curve is largely dictated by the stopping boundaries. For example, in TDS2, the design with w=(1/2,0,1/2) (Design 2) will rarely stop early for efficacy given the optimal values of e_1 and e_2 in this instance. Accordingly, there is only ever a negative average bias in the point estimate. Similar statements hold for w=(0,1/2,1/2) (Design 3) in TDS2 also. These patterns are less clear in TDS1 however, as the optimal boundaries had a smaller magnitude on average. Importantly, the adjusted procedures almost universally reduce the average bias in the point estimate. However, this does come at a small cost to the RMSE in some instances. 

For the confidence intervals, the coverage of the naive procedures varies widely, reaching nearly 98% for the design with w=(1/2,0,1/2) (Design 2), and below 92% for w=(0,1/2,1/2) (Design 3), in TDS2 for certain values of τ. In contrast, the adjusted method provides approximately the desired coverage across all designs, and all values of τ.







§ DISCUSSION


In this paper, we have described methodology which allows the design of a group sequential SW-CRT to be optimised. In addition, we have detailed an approach through which post trial analysis can be conducted to account for the inclusion of interim analyses, reducing the bias in point estimates and improving the coverage over a naive final analysis. Ultimately, it is clear that our search procedure is able to identify highly efficient group sequential SW-CRT designs. In addition, the adjusted analysis procedures were able to substantially improve performance relative to the naive methods, reducing bias and providing coverage closer to the desired level.

One important point to discuss further is the influence of the choice of w. It is clear from Table <ref> that by altering w it is possible to improve performance, in terms of the ENM, under the null or alternative hypothesis. However, setting any of the w_i equal to zero can lead to undesirable results. For example, ENM(0) for the design with w=(0,1/2,1/2) in TDS1 was extremely large compared to the two other presented designs. This fact has been noted previously in other trial design settings <cit.>. It is advisable generally to set w_i≠0 for i=1,2,3, and arguably the balanced optimal design with w=(1/3,1/3,1/3) should usually be preferred.

Another important result arose as a consequence of permitting the optimisation of X. It is extremely common in group sequential trials that the maximal required sample size is larger than that of a corresponding fixed trial design. Here though, compared to a balanced fixed-sample SW-CRT, the optimal designs did not require any additional patients. This means that the optimal sequential designs can perform no worse than a balanced fixed-sample design in terms of the required number of measurements. Thus, as was discussed at length in Girling and Hemming (2016) <cit.> and Grayling et al. (2017c) <cit.>, if it is possible to utilise an optimised (typically non-balanced) X, this should be done, as large efficiency gains will likely be possible. 

In summary, we have presented methodology to allow efficient group sequential SW-CRTs to be determined, and also to permit a corresponding post trial analysis to be performed with desirable performance. In practice, these methods are also easy to implement. We have provided code to replicate our presented results, but all that is required is the ability to evaluate multivariate normal integrals. One drawback that should be noted though is that there is no guarantee that the search procedure will converge to the truly optimal design; more likely it will only be approximately optimal. Ideally, one should therefore check convergence by performing multiple searches from different randomised starting values.

An interesting avenue for future work would be to examine whether there is any pattern to the optimal value of X depending on the choices for w and 𝒯. In particular, whether these are similar to the findings obtained for fixed sample SW-CRTs would be of great interest, as this may then permit the run time of the search procedure for finding optimal sequential designs to be improved by utilising these results as the starting point.

Finally, it is important to recognise that the practical limitations of incorporating interim analyses in to a SW-CRT as discussed in Grayling et al. (2017a) <cit.> still stand. However, these notwithstanding, future trials should in our opinion consider incorporating at least one interim analysis in to a SW-CRT whenever possible. This will provide the majority of the efficiency gains from a sequential approach, whilst minimising the alterations required to the trials conduct.
	




§ APPENDIX







 §.§  initialisation


In Section <ref> we discussed our approach for simultaneously optimising X, f, e, and m, when C had been fixed. Ultimately, we were left to find the solution to the problem defined by Equation (<ref>), and proposed to use the function  in  for this. Here, we describe our initialisation of this function.

The problem (<ref>) imposes the following constraints on {(f,e),m,S}

	
  * m∈ℕ\ℕ_1.
	
  * S_c∈ℕ_T+1^+ for c∈ℕ_C^+.
	
  * min_c∈ℕ_C^+S_c ≤ t_1.
	
  * ∏_c∈ℕ_C^+𝕀(S_c=t)=0 for t∈ℕ_T^+.
	
  * {f,e}∈ℝ^|𝒯|×ℝ^|𝒯|
	
  * f_i<e_i for i∈ℕ_|𝒯|-1.
	
  * f_|𝒯|=e_|𝒯|.

Here, constraints 2-4 are to ensure that S∈𝒮_𝒯, whilst constraints 5-7 guarantee {f,e}∈ℬ_|𝒯|.

In , discrete parameters can be incorporated for by direct specification of a set of allowed values. Consequently, constraint 2 was included easily, whilst constraint 1 could be practically achieved by setting m∈ℕ_m_max for some suitably large m_max∈ℕ\ℕ_1. For our example, we set m_max=10N_SW/CT. Moreover, to ensure constraint 3 was adhered to, we designated without loss of generality that S_1∈ℕ_t_1^+.

To include constraints 5-7, rather than optimising the parameters f and e, we transformed our problem to require optimisation of f and  r=(e_1-f_1,…,e_|𝒯|-1-f_|𝒯|-1)^⊤. For, {f,r} completely defines a pair {f,e} conforming to constraints 5-7 if {f,r}∈ℝ^|𝒯|×{ℝ^+}^|𝒯-1|, which are linear constraints supported by .

Unfortunately, constraint 4 could not be incorporated using the input parameters that can be passed to . Nor would a convenient transformation of the problem work as above. Our solution therefore was, for each parameter set {(f,e),m,S} considered in the search, to manually check whether constraint 4 was met. In the event that it was not, O(f,e,m,S) was assigned to be infinite to ensure the design was rejected.

Note that in treating m as discrete, we deviated from the approach taken in, for example, Wason (2015) <cit.>. This was advantageous in that then only a single optimisation run was required to determine the optimal implementable design.

Now whilst the above ensured our determined optimal solution would confirm to all requirements placed on the design parameters, this is only one aspect of utilising . Like many other global continuous optimisation routines,  allows the specification of several control parameters of the underlying algorithm's execution. These include an initial probability mass function for the discrete parameters, initial mean values and standard deviations for the continuous parameters, a sample size for each iteration of the algorithm, and a rarity parameter used in `elite' design selection. For more information on these, we refer the reader to Benham et al. (2017) <cit.>. 

Fortunately, the cross entropy method is relatively robust to its initialisation. Therefore we placed uniform distributions on the discrete parameters, as well as means of 0 and a (large) standard deviation of 10 on the continuous parameters. We also heeded the advice of Benham et al. (2017) <cit.>, and set the rarity parameter ρ_CE=0.01. Moreover, it was advised that the sample size for each iteration, N_CE, be chosen so that N_CEρ_CE is of the order 10 times the dimension of the optimisation problem. The dimension of our problem is C+2|𝒯|. Accordingly, we set N_CE=10000(C+2|𝒯|). Of course, these choices may be sub-optimal, and attempts could be made to improve upon them. Practically however, the search procedure is robust and efficient enough that this is not necessary. Indeed, as was seen, for the considered examples our utilised approach led to efficient design determination. 

Finally, it is important to note the redundancy which exists in the specification of the vector S, which we were forced to accept within our search procedure. This is a result of the fact that, for example, O{f,e,m,(1,2,3,4)^⊤}=O{f,e,m,(2,1,3,4)^⊤}. We may anticipate that we could remove the presence of this redundancy by enforcing S_c≥ S_c-1 for c∈ℕ_C\ℕ_2, using either a manual check on candidate design proposal (as for constraint 4 above), or an approach similar to that taken for constraints 5-7 with the introduction of r.

Unfortunately, the latter approach fails because there is then no effective means of ensuring that S_c∈ℕ_T+1^+ for c∈ℕ_C^+. Moreover, a manual check is only a viable solution when C is small. To see this, consider the probability of randomly choosing an S with S_c≥ S_c-1 for c∈ℕ_C\ℕ_2, which is given by


    ∑_s_1=1^t_1∑_s_c=s_1^T+1…∑_s_c=s_c-1^T+1∏_c=1^Cℙ(S_c=s_c).


Suppose for example that S_1 ∼ U{1,t_1} and S_c ∼ U{1,T+1} for c∈ℕ_C\ℕ_2, as is the case at the beginning of our search given the specifications above. Then, Table <ref> presents the probability (<ref>) when t_1=T/2 for several C∈ℕ\ℕ_1 and T∈ℕ\ℕ_1. Clearly, without utilising an intractably large value for N_CE, this approach is in general not feasible as the chance of selecting a design that would pass our manual check is unacceptably small.





 §.§ Simulation study


Here, we describe how we generated trial data during the simulation study. The assumed correlation structure of the accrued responses in a SW-CRT means  more complex methods are required for this, compared to the simulation of individually randomised trials for example.

First, recall that Y_t|m,X∼ N(D_t|m,Xβ,Σ_t|m,X). For simplicity, from here drop explicit reference to m and X and write Y_t∼ N(D_tβ,Σ_t). Then, for any set 𝒫={t_1,…,t_|𝒫|} with t_i∈ℕ_T^+ for i∈ℕ_|𝒫|^+, define Y_𝒫=(Y_t_1^⊤,…,Y_t_|𝒫|^⊤)^⊤. Here for any j∈ℕ_T^+, Y_j, dim(Y_j)=mC×1, is the vector of responses from time period j∈ℕ_T^+.

Using this, write

    [ Y_ℕ_t-1^+;       Y_t ]∼ N{[ D_ℕ_t-1^+;       D_t ]β, [ Σ_ℕ_t-1^+,ℕ_t-1^+       Σ_ℕ_t-1^+,t;       Σ_t,ℕ_t-1^+             Σ_t,t ]},

with for any sets 𝒫 and 𝒬, D_𝒫 and Σ_𝒫,𝒬 defined similarly as for Y_𝒫 above.

Finally, we generate the Y_t as follows. Suppose that Y_ℕ_t-1^+ has been generated for some t∈ℕ_T\ℕ_1, Y_ℕ_t-1^+=y_ℕ_t-1^+ say. Then

    (Y_t|Y_ℕ_t-1^+=y_ℕ_t-1^+) ∼ N{D_tβ + Σ_t,ℕ_t-1^+Σ_ℕ_t-1^+,ℕ_t-1^+^-1(y_ℕ_t-1^+-D_ℕ_t-1^+β),Σ_t,t - Σ_t,ℕ_t-1^+Σ_ℕ_t-1^+,ℕ_t-1^+^-1Σ_ℕ_t-1^+,1},

meaning Y_t can be generated using standard methods for generating random vectors following specified multivariate normal distributions.



§ ACKNOWLEDGEMENTS


This work was supported by the Wellcome Trust [grant number 099770/Z/12/Z to M.J.G.]; and the Medical Research Council [grant number MC_UU_00002/3 to M.J.G. and A.P.M., and grant number MC_UU_00002/6 to D.S.R. and J.M.S.W.].

99
	biblabel#1
	@bibitem#1@bibitem#1-
	eldridge[1] Eldridge, S., and Kerry, S. (2012). A Practical Guide to Cluster Randomised Trials in Health Services Research, Chichester: Wiley.
	
	campbell[2] Campbell, M. J., and Walters, S. J. (2014). How to Design, Analyse and Report Cluster Randomised Trials in Medicine and Health Related Research, Chichester: Wiley.
	
	hayes[3] Hayes, R. J., and Moulton, L. H. (2017). Cluster Randomised Trials, Boca Raton: CRC Press, Second Edition.
	
	hussey[4] Hussey, M. A., and Hughes, J. P. (2007). Design and analysis of stepped wedge cluster randomized  trials, Contemp Clin Trials 28: 182-191. DOI: 10.1016/j.cct.2006.05.007.
	
	woertman[5] Woertman, W., de Hoop, E., Moerbeek, M., Zuidema, S. U., Gerritsen, D. L, and Teerenstra, S. (2013). Stepped wedge designs could reduce the required sample size in cluster randomized trials, J Clin Epidemiol 66: 752-758. DOI: 10.1016/j.jclinepi.2013.01.009.
	 
	baio[6] Baio, G., Copas, A., Ambler, G., Hargreaves, J., Beard, E., and Omar, R. Z. (2015). Sample size calculation for a stepped wedge trial, Trials 16: 354. DOI: 10.1186/s13063-015-0840-9.
	
	hooper[7] Hooper, R., Teerenstra, S., de Hoop, E., and Eldridge, S. (2016). Sample size calculation for stepped wedge and other longitudinal cluster randomised trials, Stat Med 35: 4718-4728. DOI: 10.1002/sim.7028.
	
	hemming[8] Hemming, K., Taljaard, M., and Forbes, A. (2017). Analysis of cluster randomised stepped wedge  trials with repeated cross-sectional samples, Trials 18: 101. DOI: 10.1186/s13063-017-1833-7.
	
	wang[9] Wang, R., and De Gruttola, V. (2017). The use of permutation tests for the analysis of parallel and stepped-wedge cluster-randomized trials, Stat Med 36: 2831-2843. DOI: 10.1002/sim.7329.
	
	kasza[10] Kasza, J., Hooper, R., Matthews, J. N. S., and Forbes, A. B. (2017). Impact of non-uniform correlation structure on sample size and power in multiple-period cluster randomised trials, Stat Meth Med Res DOI: 10.1177/0962280217734981.
	
	grayling1[11] Grayling, M. J., Wason, J. M. S., and Mander, A. P. (2017a). Group sequential designs for stepped-wedge cluster randomised trials, Clin Trials 14: 507-517. DOI: 10.1177/1740774517716937.
	 
	dehoop[12] de Hoop, E., van der Tweel, I., van der Graaf, R., Moons, K. G. M., van Delden, J. J. M., Reitsma, J. B., and Koffijberg, H. (2015). The need to balance merits and limitations from different disciplines when considering the stepped wedge cluster randomized trial design, BMC Med Res Methodol 15: 93. DOI: 10.1186/s12874-015-0090-2.
	 
	grayling2[13] Grayling, M. J., Wason, J. M. S., and Mander, A. P. (2017b). Stepped wedge cluster randomized controlled  trial designs: a review of reporting quality and design features, Trials 18: 33. DOI: 10.1186/s13063-017-1783-0.
	
	lawrie[14] Lawrie, J., Carlin, J. B., and Forbes, A. B. (2015). Optimal stepped wedge designs, Stat Probabil Lett 99: 210-214. DOI: 10.1016/j.spl.2015.01.024.
	
	girling[15] Girling, A. J., and Hemming, K. (2016). Statistical efficiency and optimal design for stepped cluster studies under linear mixed effects models, Stat Med 35: 2149-2166. DOI: 10.1002/sim.6850.
	
	thompson[16] Thompson, J. A., Fielding, K., Hargreaves, J., and Copas, A. (2017). The optimal design of stepped wedge trials with equal allocation to sequences and a comparison to other trial designs, Clin Trials 14: 639-647. DOI: 10.1177/1740774517723921.
	
	li[17] Li, F., Turner, E. L., and Preisser, J. S. (2018). Optimal allocation of clusters in cohort stepped wedge designs, Stat Probabil Lett 137: 257-263. DOI: 10.1016/j.spl.2018.02.002.
	
	wason1[18] Wason, J. M. S., and Jaki, T. (2012). Optimal Design of Multi-arm Multi-stage Trials, Stat Med 31: 4269-4279. DOI: 10.1002/sim.5513 .
	
	wason2[19] Wason, J. M. S., Mander, A. P., and Thompson, S. G. (2012). Optimal Multistage Designs for Randomised Clinical Trials with Continuous Outcomes, Stat Med 31: 301-312. DOI: 10.1002/sim.4421.
	
	an[20] An, Q. and Chang, M. (2016). Computation of optimal group sequential designs, Sequential Anal 35: 453-464. DOI: 10.1080/07474946.2016.1238256.
	
	jennison[21] Jennison, C., and Turnbull, B. W. (2000). Group Sequential Methods with Applications to Clinical Trials, Boca Raton: Chapman and Hall/CRC.
	
	grayling3[22] Grayling, M. J., Mander, A. P., and Wason, J. M. S. (2017c). Admissible multi-arm stepped-wedge cluster randomized trial designs, arXiv:1710.03496.
	
	benham[23] Benham, T., Duan, Q., Kroese, D., and Liquet, B. (2017). CEoptim: Cross-Entropy R Package for Optimization, J Stat Soft 76: 1-29. DOI: 10.18637/jss.v076.i08.
	
	tsiatis[24] Tsiatis, A. A., Rosner, G. L., and Mehta, C. R. (1984). Exact confidence intervals following a group sequential test, Biometrics 40: 797-803. DOI: 10.2307/2530924.
	
	whitehead[25] Whitehead, J. (1986). On the bias of maximum likelihood estimation following a sequential test, Biometrika 73: 573-581. DOI: 10.1093/biomet/73.3.573.
	
	emerson[26] Emerson, S. S., and Fleming, T. R. (1990). Parameter estimation following group sequential hypothesis testing, Biometrika 77: 875-892. DOI: 10.1093/biomet/77.4.875.
	
	armitage[27] Armitage, P. (1957). Restricted sequential procedures, Biometrika 44: 9-56. DOI: 10.1093/biomet/44.1-2.9.
	
	siegmund[28] Siegmund, D. (1978). Estimation following sequential tests, Biometrika 65: 341-349. DOI: 10.1093/biomet/65.2.341.
	
	bashour[29] Bashour, H. N., Kanaan, M., Kharouf, M. H., Abdulsalam, A. A., Tabbaa, M. A., and Cheikha, S. A. (2013). The effect of training doctors in communication skills on women's satisfaction with doctor-woman relationship during labour and delivery: a stepped wedge cluster randomised trial in Damascus, BMJ Open 3: e002674. DOI: 10.1136/bmjopen-2013-002674.
	
	mander[30] Mander, A. P., Wason, J. M. S., Sweeting, M. J., and Thompson, S. G. (2012). Admissable two-stage designs for phase II cancer clinical trials that incorporate the expected sample size under the alternative, Pharm Stat 11: 91-96.DOI: 10.1002/pst.501.
	
	wason3[31] Wason, J. (2015). OptGS: An R Package for Finding Near-Optimal Group-Sequential Designs, J Stat Soft 66: 1-13. DOI: 10.18637/jss.v066.i02.



