Modeling clustering and treatment effect heterogeneity in parallel and stepped-wedge cluster randomized trials.

Cluster randomized trials are frequently used in health service evaluation. It is common practice to use an analysis model with a random effect to allow for clustering at the analysis stage. In designs where clusters are exposed to both control and treatment conditions, it may be of interest to examine treatment effect heterogeneity across clusters. In designs where clusters are not exposed to both control and treatment conditions, it can also be of interest to allow heterogeneity in the degree of clustering between arms. These two types of heterogeneity are related. It has been proposed in both parallel cluster trials, stepped-wedge, and other cross-over designs that this heterogeneity can be allowed for by incorporating additional random effect(s) into the model. Here, we show that the choice of model parameterization needs careful consideration as some parameterizations for additional heterogeneity induce unnecessary or implausible assumptions. We suggest more appropriate parameterizations, discuss their relative advantages, and demonstrate the implications of these model choices using a real example of a parallel cluster trial and a simulated stepped-wedge trial.

INTRODUCTION

Cluster randomized trials (CRTs) randomize entire clusters of individuals to treatment or control conditions.1 Another related design, which again randomizes entire clusters, is one in which all clusters are randomized to a treatment sequence that is either treatment followed by control, or control followed by treatment.2 Stepped‐wedge cluster randomized trials (SW‐CRTs) also randomize entire clusters but randomize clusters to a sequence of time periods spent in the control condition followed by time periods spent in the treatment condition.3 All cluster trials need to allow for the non‐independence of observations within the same cluster. This is typically done by using a linear or generalized linear mixed model with a random effect for cluster. 1

The effect of the treatment might vary across clusters, and it can be of interest to explore any treatment effect heterogeneity at the analysis stage. In stepped‐wedge designs or other cluster randomized designs in which clusters are exposed to both treatment and control, this treatment effect heterogeneity can be identified.4, 5, 6, 7 In the conventional parallel CRT, there may also be treatment effect heterogeneity. However, in trials such as the parallel CRT in which each cluster is exclusively exposed to either the treatment or control condition, this heterogeneity cannot be separated out from the cluster effect. But a related issue is where the variability between clusters may differ according to treatment arm in parallel cluster trials.8 While these two sources of heterogeneity are inherently different, they have similar implications for design and analysis when this heterogeneity is modeled using random effects. There is a confusing array of different parameterizations that have been proposed in the literature to model these different but related sources of heterogeneity, which all make different assumptions. This paper aims to review these models and make some recommendations for which models make fewer assumptions.

In this paper, we demonstrate that in cross‐sectional CRTs, the choice of model parameterization to incorporate treatment effect or cluster heterogeneity has important underlying assumptions. We demonstrate the implications of these model choices using practical examples of a parallel cluster trial and a simulated stepped‐wedge trial.

THEORETICAL MODELS

Differential clustering in parallel arm CRTs

In parallel arm CRTs, clusters are either fully exposed or unexposed to the treatment. Therefore, any treatment effect heterogeneity cannot be disentangled from the cluster effect. However, the correlation between observations within clusters might vary across treatment arms. The potential for such differential correlation is intuitive because different treatments might be expected to induce homogeneity—or even heterogeneity. This is particularly relevant in trials with differential clustering between arms—perhaps where a group therapy is compared with an individual therapy. Here, we review the parameterizations that have been proposed for these models in the context of parallel cluster trials and outline their respective assumptions. In Appendix A, we show derivations of the correlations for model 2b (one of the more complex models). In the next section, we show how these models for differential clustering are related to treatment effect heterogeneity models.

Basic model: single random effect

Let us consider a 2‐arm parallel CRT. We assume that there is a continuous outcome y , where i=1,…,m represents the individual and j=1,…,k represents the clusters. We assume that there are k clusters each of equal size m. One analysis model for this simple setup is1 where x represents the treatment indicator. We code this treatment indicator as 1 for the treatment and 0 for the control (coding of contrasts is important, a point to which we return in the discussion). Then, θ is the treatment effect, and α is a random effect for cluster j. Under this model, the correlation between two observations in the same cluster, the intra‐cluster correlation (ICC), will be

The model thus assumes that for any cluster, observations within that cluster share a common correlation ρ and this is the same for both treatment and control clusters.

Model extension 1: two separate random effects

To allow for different within‐ and between‐cluster variability in treatment and control clusters, two separate random effects, one for treatment and one for control, are incorporated: where α(T) and α(C) represent random cluster effects for the treatment and control clusters. We use the sub‐scripts capital T and capital C to denote the treatment and control clusters. This induces differing ICCs in control and treatment clusters. The ICC in the control clusters will be and the ICC in the treatment clusters will be Using this parameterization, there is no restriction on whether ρ is bigger or smaller than ρ . The model allows for random variation between the control clusters (variance ), and a different random variation between treatment clusters, .

Model extension 2: a random interaction

Seen perhaps as a more intuitive way to model heterogeneity in cluster variability across arms, an alternative to the above model is to use a parameterization that includes a random (simple multiplicative) interaction between the treatment covariate and cluster: where α(M) and α(I) are independent. We use the notation M to represent the main effect term, and I the interaction term. We have used different notation to that assumed in model extension 1 to make explicit that this is a different model. Under this assumption, the ICC in the control clusters will be and the ICC in the treatment clusters will be

Using this parameterization, then ρ <=ρ . That is, there is an implicit assumption that the total variance in the treatment clusters is greater than or equal to the total variance in the control clusters. This of course is not a tenable assumption in all situations. This means that if the treatment comparison was switched around, so that the control clusters were coded as 1 and the treatment clusters as 0 model results would differ, an aspect of a model fit which is clearly undesirable.

Treatment effect heterogeneity in cluster randomized designs

The previous section focused on differential clustering by arm (or heteroscedasticity), as appropriate in parallel designs as treatment is nested within clusters. We now move on to consider treatment effect heterogeneity in designs in which the treatment is crossed with cluster (cross‐over and stepped‐wedge trials).

We consider a conventional stepped‐wedge study, with S sequences to which the k clusters are randomly allocated, and where at each of the S+1 time points, a cross‐section of observations is taken.3, 9 Because each cluster in a SW‐CRT is exposed both to the treatment and control condition, treatment effects can be estimated within each cluster. While for simplicity, we focus on the SW‐CRT, models that follow will be generalizable to other cluster trials in which treatment is crossed with cluster. We outline various different models along with the corresponding correlations and give more details on relationships between the various parameterizations in Appendix C.

We extend the basic model (1) for parallel CRTs to the SW‐CRT, by incorporating fixed effects for each period9: where s=1,…,S denotes the period and π is a fixed effect for each period, and m now represents the size of each cluster at each period, y the outcome for individual i in cluster j at time s, and x the corresponding treatment indicator.

Model extension 1: two separate random effects

The first model we consider allows for different between‐cluster variability in treatment and control clusters and follows from model extension 1 above for CRTs. Here, we incorporate two separate random effects, one for treatment and one for control. In contrast to when the clustering is nested within treatment arm, when the clustering is crossed with treatment, the model with a non‐zero covariance between these two random effects becomes identifiable: again α(T) and α(C) represent random cluster effects for the treatment and control conditions, where and where σ is the possibly non‐zero covariance between these two random effects.

Model extension 1a: two separate independent random effects

In model extension 1a, we make the assumption that the correlation between these two separate random effects is zero (ie, σ =0). This induces differing ICCs in control and treated conditions. The correlation between two observations in the same cluster, both exposed to the control condition, will be and the correlation between two observations in the same cluster, both exposed to the treatment condition, will be

Using this parameterization, there is no restriction on whether ρ is bigger or smaller than ρ . However, an additional correlation also exists: the correlation between two observations within the same cluster, but one in the control condition and one in the treatment condition, which we call ρ . Under this model parameterization, this correlation, ρ , is assumed to be zero (clearly undesirable).

Model extension 1b: two separate nonindependent random effects

In model extension 1b, the correlation between these two separate random effects is allowed to be non‐zero. This induces differing ICCs in control and treatment clusters, and importantly a non‐zero correlation between observations within the same cluster but different treatment exposure. These correlations are

As in model 1a above, there are no restrictions on the relative magnitude of ρ and ρ , but in contrast to model 1a, this model no longer includes the restriction that ρ =0.

We now consider extensions that clearly include an interaction between treatment and cluster. Again, because the clustering is crossed with treatment, the model with a non‐zero covariance between these two random effects becomes identifiable: where the two random effects α(M) and α(I) (where M represents the main effect term, and I the interaction term) have the following distribution:

Model extension 2a: an independent random interaction

When the covariance term σ between the two random effects is assumed to be zero, the respective correlations are

Using this parameterization, we see that there is again the restriction that ρ is greater than or equal to ρ , as was the case when this model was used to allow for differential clustering between arms; furthermore, there is also the restriction that ρ is greater than or equal to ρ .

Model extension 2b: a nonindependent random interaction

Allowing the covariance term to be non‐zero, the correlations become

Inclusion of the non‐zero covariance term means that the model does not induce assumptions on the variance in the treated periods being greater than the variance in the control periods, nor that ρ is greater than or equal to ρ .

Model extension 3: an alternative random interaction

An alternative parameterization to allow for an interaction with treatment and cluster is one in which the random variation for the cluster effects is partitioned into a common part and an exposure specific part: and that Note that we have used a change of notation from α to γ for the random effects to avoid confusion with models 1 and 2.

Model extension 3a: an alternative random interaction with constrained total variance

When it is assumed that , the correlations between two observations in the same cluster and both exposed to the same treatment (whether that be treatment or control) will be and the correlations between two observations in the same cluster and both exposed to different treatments will be

This model assumes that the total variance under exposure to treatment is equivalent to the total variance under exposure to control. This parameterization therefore imposes a restriction on the correlations within clusters being identical for two observations that are both treated and two observations that are both not treated. Furthermore, this model also makes the assumption that ρ is less than or equal to ρ (=ρ ).

Model extension 3b: an alternative random interaction without constrained total variance

An alternative parameterization, to model 3a, which relaxes the assumption that the total variance in the treatment arm is equivalent to the total variance in the treatment arm, is one in which the parameters and are not constrained to be the same. The correlations then become

While this model does not make the assumption that one of the ICCs is larger than the other, and neither makes the assumption that the correlation between two observations, one treated and one control, is zero, this model does make assumptions that the three random effects are uncorrelated. This does induce additional assumptions (Appendix C) that there is a positive correlation between the treatment and control random effects (in terms of model 1 notation as increases so will see Appendix 3).

MODEL FITTING

We have fitted these models in Stata 14 using the xtmixed function, REML methods, and default settings. Careful model specification and understanding of the language syntax is required to avoid model misspecification. We have included sample Stata code (Appendix B).

EXAMPLE OF DIFFERENTIAL CLUSTERING IN A PARALLEL CRT

In this example, we illustrate results from fitting models 1 and 2 to a set of data from a parallel CRT. The trial we choose as an example is a parallel cluster trial conducted in 53 schools (clusters) of a behavioral treatment to prevent obesity in school‐aged children.10 The outcome we consider here is the child's body mass index measured at the end of the trial. There are a total of 689 observations in the treatment arm and 778 observations in the control arm, with an average cluster size of 24.

This dataset usefully highlights that one of the parameterizations (model 2) is not a good one. In this CRT, there is little variability between the treatment clusters, the estimated ICC from a simple analysis of variance on the treatment clusters only is very close to zero, whereas in the control clusters, it is higher at 0.057 (Table 1). Model 1 is able to capture this pattern in the data (estimating the ICCs in control and treatment arm as 0.050 and 0.000, respectively), whereas model 2 is unable to, estimating the ICCs in the control and treatment arm as 0.022 and 0.022 (Table 1). The estimated treatment effects differ slightly between the two approaches, but likelihoods do not indicate any preference for either model.

Table 1

Example of differential clustering in a parallel cluster randomized trial

			First random	Second random	Residual variance,
	Log‐likelihood	Treat effect, SE	effect, SE	effect, SE	SE	ICC, 95% CI
No differential
clustering	−2100.47	0.111 (0.092)	0.195 (0.052)	1.288 (0 .026)	0.022 (0 .008, 0.062)
					Control	Intervention
Stratified model–control arm			0.304 (0.072)		1.236 ( 0.034)	0.057 ( 0.023, 0.134)
Stratified model–treatment arm			0.000 (0.000)		1.337( 0.039)		0.000 (0.000, 0.000)
Differential clustering models
Model 1—two separate random effects	−2097.09	0.098 (0.094)	0.000 (0.000)	0.296 (0.073)	1.285 (0.026)	0.050	0.000
Model 2—random interaction	−2100.47	0.111 (0.092)	0.000 (0.000)	0.195 (0.052)	1.288 (0.026)	0.022	0.022

Abbreviations: CI, confidence interval; ICC, intra‐cluster correlation; SE, standard error.

To investigate the actual consequences on estimation of model parameters when modeling treatment effect heterogeneity in cluster trials in which treatment is crossed with cluster, we undertook a small simulation study. In this simulation study, we have compared model performance statistics (ie, bias and coverage) across the 6 models (1a to 3b above) for the estimates of treatments effects (and standard errors) and the estimates of the variance of the random effects parameters (ie, the ICCs in both treatment and control arms).

This simulation study is not exhaustive and does not cover the range of possible scenarios but serves to demonstrate only some of the possible implications of model choice. This simulation study is limited to a simple stepped‐wedge study. We assume a large number of clusters (100) and a large cluster size (1000) to avoid any issue of small sample sizes and assume a stepped‐wedge study with 4 steps (ie, 25 clusters each randomized to one of the 4 sequences). We assume a treatment effect of 0 and residual variance of 1 (fixed throughout), and no secular trend (ie, each π =0).

Data were generated from a linear mixed model as in model extension 1, adding cluster specific variation derived from three specified ICCs, one for the control condition, one for the treatment condition, and one for both. From this, we then derived variance components and a non‐zero covariance term. These ICCs (the scenarios) were chosen to represent a broad range of scenarios, some of which are concordant with the models presented here and some of which were discordant.

The first scenario (scenario 1) examines what happens when we fit a model that is more complex than warranted by the data. In the first scenario, there is no differential heterogeneity, and the ICC under both treatment conditions, and between treatment conditions, is set at a typical value of 0.01. In this scenario, the data are concordant with all model extensions presented in this paper, except model 1a (which assumes ρ =0).

The second scenario (scenario 2) illustrates what happens when we fit these models to data in which there is differential clustering—but differential clustering that is not too discordant with most models. In the second scenario, there is differential heterogeneity: the ICC in the control condition being 0.01, the ICC in the treatment condition being 0.05, and the ICC between control and treatment conditions 0.005. In this scenario, the data are discordant with model 2a(which assumes ρ >ρ ); model extension 3a (which assumes ρ =ρ ); and model 1a (which assumes ρ =0).

In the third scenario, there is again differential heterogeneity across both arms, but the direction is reversed: the ICC is 0.05 in the control condition and lower in the treatment condition (0.01), and between conditions 0.001. So in this situation, there is differential heterogeneity, and this is very discordant with some models. In this scenario, the data are discordant with model 2a (because ρ <ρ and ρ >ρ ); model extension 3a (because ρ ≠ρ ); and model 1a (because ρ ≠0).

We ran 1000 simulations for each scenario, which meant nominal coverage would be estimated to between 93.6% and 96.4%. We observed most models converged (8/1000 did not converge under model 1a); treatment effects showed little evidence of bias (Table 2). For the treatment effects, we observed some small departures from nominal coverage in some models, and some small differences in coverage between different models. However, these differences were small and likely to be due to sampling variation (values of coverage outside of the range 93.6 and 96.4 are possibly due to sampling variation). Model 1a exhibited higher than nominal coverage (100%) under scenario 1, explained by the zero covariance assumption of this model (Appendix C).

Table 2

Simulation study of impact of model choice for treatment effect heterogeneity in cluster randomized designs

	Model 1a	Model 1b	Model 2a	Model 2b	Model 3a	Model 3b
Scenario 1 (ICC = 0.01 in control clusters; ICC = 0.01 in intervention clusters; ICC = 0.01 in clusters crossed with treatment)
Treatment effects
Absolute bias	0.00	0.00	0.00	0.00	0.00	0.00
Coverage	100.0%	93.7%	95.4%	95.6%	95.6%	95.8%
ICC estimates (percentage bias)
ICC in treatment arm	−0.03	0.46	1.51	0.48	0.42	0.28
ICC in control arm	0.08	0.56	−0.77	0.51	0.42	0.44
Scenario 2 (ICC = 0.01 in control clusters; ICC = 0.05 in intervention clusters; ICC = 0.005 in clusters crossed with treatment)
Treatment effects
Absolute bias	0.00	0.00	0.00	0.00	0.00	0.00
Coverage	97.2%	96.1%	95.6%	96.1%	96.2%	96.1%
ICC estimates (percentage bias)
ICC in treatment arm	0.31	0.32	14.31	0.32	−40.4	0.24
ICC in control arm	−0.95	−0.95	−4.93	−0.96	210.8	−0.93
Scenario 3 (ICC = 0.05 in control clusters; ICC = 0.01 in intervention clusters; ICC = 0.001 in clusters crossed with treatment)
Treatment effects
Absolute bias	0.00	0.00	0.00	0.00	0.00	0.00
Coverage	96.0%	94.9%	93.0%	94.9%	94.6%	94.8%
ICC estimates (percentage bias)
ICC in treatment arm	−0.29	−0.29	907.7	−0.29	208.8	−0.29
ICC in control arm	−0.63	−0.63	−8.81	−0.63	−40.7	−0.63

Abbreviation: ICC, intra‐cluster correlation.

Model 1: two separate random effects one for cluster and one for treatment condition (a [b]: with a zero [non‐zero] covariance term).

Model 2: random interaction between treatment and cluster (a [b]: with a zero [non‐zero] covariance term).

Model 3: two separate random effects one for cluster and one for treatment condition, with a partition into a common part (a [b]: with same [different] variance in treatment and control arms).

Comparing the implications on the estimates of the ICCs, we observed that models 2a and 3a perform particularly poorly. Under scenario 2, model 3a shows a large percentage bias for the ICC in both the control and treatment conditions (as much as 200%); model 2a shows moderate degree of bias (up to 14%). For scenario 3, models 2a and 3a again perform particularly poorly (up to 900% bias for model 2a).

DISCUSSION

We have demonstrated that when heterogeneity in between‐cluster variability or treatment effect heterogeneity is included in models for the analysis of a CRT, the parameterization of the model and the resulting implicit assumptions are important for estimation of ICCs and in limited scenarios on the standard errors of treatment effects. We have identified and illustrated that some parameterizations make assumptions that are not only not obvious from the model itself but also will not always be tenable.

What is already known

Individually randomized trials

Where individually randomized studies have included random effects to model treatment effect heterogeneity, random interaction models allowing for a correlation between center main effects and treatment by center interactions have been used.11, 12, 13

Of note, Lee and Thompson discuss treatment effect heterogeneity where there is clustering in individually randomized trials and allow for non‐zero covariance terms (model 2b), and while they clearly express a preference for a model that makes fewer assumptions, the assumptions implicit in the parameterizations are not clear.13 The model proposed by Brown makes an assumption of equal variances in both arms (ie, is our model 3a).11 Others have not clearly outlined the model parameterization choice at all.12 Thus, even these scholarly articles we identified in the individually randomized literature did not make it apparent what model assumptions were being made.

Parallel CRTs

One situation where between‐cluster heterogeneity is likely to be important is in designs with differential clustering between arms.14, 15 However, to date, these models have concentrated on situations in which it is natural to expect that the correlation between the observations within the treatment clusters is greater than that in the control clusters. Situations like this arise, for example, when the treatment induces clustering, for example, group therapy.8 Model choice for such data very naturally fits in line with model 1a, and this has been the model of choice.8

CRTs with treatment crossed with cluster

In the cross‐over trial literature, Turner et al have used a parameterization very similar to our model 3b, and so while do not make assumptions about variances being larger in the treatment arm,6 the model does make other “hidden” assumptions (Appendix C). Yet others have followed more restrictive forms of parameterizations in the stepped‐wedge literature and have chosen to use model 2a without appreciating the assumptions implicit in this model.4, 5

Individual patient data meta‐analysis

In individual patient data meta‐analysis models, which also use mixed models and also allow for treatment effect heterogeneity (across studies rather than clusters), it has become common practice to include study effects as fixed rather than random. This means the issue of how to parameterize the additional random study by treatment heterogeneity is non‐problematic. Despite this, where random effects models have been used, model assumptions have been clearly articulated,16 which means clear guidance does exist for those wishing to model treatment effect heterogeneity using random effects. Of some note, it has been observed that the fully flexible model (model 2b here) can sometimes not converge, and so an alternative contrast parameterization (+1/2, −1/2) has been advocated for the model that has zero covariance terms that has a different constraint, namely, that the variances are equal between arms.17

Limitations

There are alternative ways to circumvent this issues raised in this paper. For example, an appropriate choice of contrast coding or selection of the appropriate arm to be coded as the arm under which heterogeneity is assumed to be greater. However, we do not endorse the choice of models that depend on choice of contrast coding, as this approach offers more potential for error and model misspecification on behalf of the user.

Our simulation study was small insofar, as it considered only a very limited number of scenarios. However, it illustrates that even in the case of large studies, where estimates of the ICCs are needed, model choice can have large implications on bias of these parameters and could induce unnecessarily large standard errors on treatment effects. It is likely that in other scenarios, where, for example, the cluster sizes are small or the number of clusters are small, the implications on the ICCs would be greater. We also limited our consideration to mixed models, whereas population average models offer an alternative approach.

While considering an interaction with cluster and treatment, we did not consider interactions between time and treatment or time and cluster.7, 18 Such interactions are likely to be of importance in trials that are longitudinal in nature as a decay in the strength of the correlations within a cluster will be expected.

Recommendations

In parallel CRTs, clusters are nested within treatment arms, whereas in the SW‐CRT and cross‐over designs, clusters are crossed with treatment arms. So in parallel CRTs, it does not make sense to consider clusters crossed with treatment, but rather clusters nested within treatment arms. This leads directly to model 1 as the natural approach for parallel cluster trials to model differential clustering between arms. It is also apparent from this that there are only two parameters here, and only two parameters estimable because each cluster gets only one treatment. Therefore, immediately any model with three parameters for the random effects structure is over‐parameterized.

In a SW‐CRT and cross‐over trial, however, clusters are now crossed with treatments, so interaction terms as well as potential for covariance between the two random effects are now estimable. Model 1a no longer becomes a good choice, as it is unlikely to be a reasonable assumption that the correlation between two observations in the same cluster but of different treatment status are independent. Model 2a is unlikely to be a good choice, as it assumes the correlation between two observations, one treated and the other not, is greater than two not treated. Model 3a is equally unlikely to be a good choice, as it assumes the total variance is equivalent across treatment conditions. Finally, model 3b, while not performing poorly in our simulation study, does induce additional assumptions. So, in the analysis of stepped‐wedge and cross‐over trials, it may be more natural to follow approaches that make minimal assumptions, that is, model 1b or 2b.12 Testing for simpler models (1a or 2a) will involve a different set of assumptions, which are clearly evident when placing the models into a common parameterization (Appendix C).