Estimating intra-cluster correlation coefficients for planning longitudinal cluster randomized trials: a tutorial

Background

The cluster randomized trial (CRT) is an essential design for evaluating interventions delivered at the group level or when there are scientific or logistical reasons precluding individual randomization.¹ Unlike individually randomized trials, CRTs randomize intact groups of participants to either intervention or control conditions. A key characteristic of this design is that the randomization is at the level of the cluster, and individuals within the same cluster are more similar than individuals from different clusters; this should be taken into consideration in the design and analysis.² The similarity of outcomes within clusters is typically measured by the intra-cluster correlation coefficient (ICC).³ Substantial methodological development for longitudinal CRTs (also referred to as multiple-period CRTs) has taken place in recent years. In a longitudinal design, the outcomes are measured multiple times in each cluster, typically before and after intervention delivery. These outcomes can be measured on different participants over time (i.e. a repeated cross-sectional design) or the same participants (i.e. a cohort design).⁴ Three major types of longitudinal CRTs are multiple-period before and after parallel-arm (Figure 1a), multiple cross-over⁵ (Figure 1b) and stepped-wedge cluster randomized trials (SW-CRT) (Figure 1c).⁶

Two recent articles have summarized existing sample size calculation methods and statistical models for longitudinal CRTs.^7,8 A key feature differentiating existing methods is the type of correlation structure, i.e. how ICCs are assumed to change over time. Available sample size calculators require specification of an assumed correlation structure for the planned design and corresponding estimated correlation parameters. Parameters that must be specified are a within-period ICC, a cluster autocorrelation coefficient (CAC) and, in the case of a cohort design, a within-individual ICC.

A practical challenge for researchers planning a new trial is how to obtain estimates for these coefficients. The choice of correlation structure can affect the sample size requirements and the statistical inferences.^8,9 One possibility is to use estimates reported in previously published trials. However, despite the various CONSORT extensions for CRTs recommending that authors report the obtained ICC values for their trials,^10,11 adherence to this recommendation has been low.¹² Furthermore, even when reported, estimates may not be suitable for longitudinal correlation structures. A database reporting suitable ICC values for a selection of longitudinal CRTs with continuous outcomes has recently become available,¹³ but the range of outcomes remains limited. With the increasing availability of routinely collected databases that can be used for outcome assessment in a CRT, as well as data availability from published trials, researchers may be able to obtain access to real data for more reliable estimation of ICC values in advance of a new trial. However, researchers may lack knowledge of how to estimate complex correlation parameters from longitudinal clustered data.

In this tutorial, our goal is to provide readers with the tools to estimate correlation parameters required for longitudinal CRTs with continuous and binary outcomes. We first review the main types of mixed-effect statistical models and correlation structures available for longitudinal CRTs and define the required input parameters for sample size calculation under these models. We then describe the required data structure for estimating the relevant correlation values and highlight some practical considerations for implementation. We provide code in R, SAS and Stata to empirically estimate the correlation parameters under repeated cross-sectional and cohort designs and we make available an RShiny app that does not require access to statistical software. We conclude our tutorial by identifying some gaps in the literature.

Available correlation structures for longitudinal CRTs

A common analytical approach for longitudinal CRT is mixed-effects regression with fixed (discrete) effects for periods and treatment.¹⁴ In the methodological literature, substantial attention has been paid to the random-effects specification, with three main types of correlation structures defined by what has become known as the Hussey and Hughes (exchangeable) model,¹⁵ the Hooper and Girling (nested/block exchangeable) model,^16,17 and the Kasza (exponential decay) model.¹⁸ A visualization of these three correlation structures is presented in Supplementary File Section 1 (available as Supplementary data at IJE online).

The first statistical model we introduce for longitudinal CRTs was developed by Hussey and Hughes in 2007.¹⁵ Assuming a continuous outcome, this model can be expressed as:

where $y_{\mathit{ijk}}$ is the outcome of the $k$-th individual from the $i$-th cluster and $j$-th period. The baseline average response (average outcome in the first period under the control condition) is denoted as $\mu$; ${\alpha }_i$ is the random intercept for the cluster, which is assumed to follow the $N\left(0,{\sigma }_{\alpha }^2\right)$ distribution; and ${\beta }_j$ is the categorical time effect in the $j$-th period (with ${\beta }_1=0$ for identifiability). The binary treatment indicator is denoted $X_{ij}$ for the $i$-th cluster in the $j$-th period, where 0 is the control and 1 is the intervention condition. $\delta$ is the treatment effect. The error term is denoted by ${\epsilon }_{\mathit{ijk}}$ and assumed to be independently and identically distributed as $N\left(0,{\sigma }_{\epsilon }^2\right)$.

This model includes a random intercept for the cluster only, and thus, the ICC is assumed to be constant across different periods (called the exchangeable correlation structure). This may be a rather strong assumption, considering that factors that may influence outcomes can vary over time. Moreover, since the exchangeable model allows the variance of the treatment effect estimate to approach 0 with increasing cluster size,¹⁹ use of the exchangeable correlation for designing longitudinal CRTs is not recommended. An improvement is to allow the outcomes of two individuals within the same cluster but in different periods to weaken with increasing time separation; i.e. to allow for a within-period ICC (the correlation between two individuals in the same cluster and the same period) and a different between-period ICC (correlation between two individuals in the same cluster but different periods).

To allow for a different within- and between-period ICC, Hooper et al.¹⁶ and Girling and Hemming¹⁷ added an additional random cluster-by-period interaction term (${\gamma }_{ij} \sim N\left(0,{\sigma }_{\gamma }^2\right)),$ independent to the random intercept for clusters, which is written as:

This model defines a correlation structure referred to as the nested exchangeable correlation.²⁰ The additional random effect allows the cluster means to vary randomly across periods. The within-period and between-period ICCs between any pair of observations $k, m$ in any two periods $j, l$ can be defined as follows:

Within-period ICC (where $j=l$):

Between-period ICC:

Note that the between-period ICC is constant in this model regardless of the distance between $i$ and $j$. The ratio of the between-period and within-period ICCs is called the CAC. The CAC is typically less than 1; a CAC value of 1 implies the exchangeable correlation.

The Hooper and Girling model was later extended by Kasza et al.,¹⁸ where they further allowed the between-period ICC to decay over time. This model can be written as:

where the vector of random cluster-by-period effects in each cluster ${\mathit{\gamma }}_{\mathit{i}}={\left({\gamma }_{i1},\dots ,{\gamma }_{iJ}\right)}^{\prime }$ is assumed to follow the $N\left(0,{\sigma }_{\gamma }^2\tilde{Z}\right)$distribution. Kasza et al.¹⁸ focused on $\tilde{Z}$ as an autoregressive (AR(1)) structure where the between-period ICC decays exponentially at the rate of $r$ per period. Therefore, this model is referred to as the exponential decay model. For simplicity, the rate of decay per period, r, is called the CAC although it defines a different parameter compared with the Hooper and Girling model.

In the case of binary outcomes, data are often analysed using mixed-effects logistic regression. The mathematical expression of the model under each correlation structure is similar to the above except $y_{\mathit{ijk}} \mathrm{is} \mathrm{replaced} \mathrm{by} \mathit{log}\left(\frac{p_{\mathit{ijk}}}{1-p_{\mathit{ijk}}}\right), \mathrm{where} p_{\mathit{ijk}}$ where the probability of the binary outcome for the $k$-th individual from the $i$-th cluster and $j$-th period and the error term ${\epsilon }_{\mathit{ijk}}$ is omitted. The definition of the ICC for binary outcomes, however, is less straightforward. For mixed-effects logistic regression, one simple approach is to work on the latent response scale²¹ so that the ICC can still be defined as the ratio of the between-cluster variance and the total variance with the residual variance defined as ${\pi }^2/3$. Under this approach, the within-period ICC and CAC can be defined analogously to the linear mixed model case.

In the case of a cohort design, an additional random effect for the repeated measures on the same individual $(k$) in cluster $i$ (denoted as ${\varphi }_{ik}$) may be added to equations (1), (2) and (3) to model a constant correlation in repeated measures on the same individual over time. This random effect is assumed to follow the $N\left(0,{\sigma }_{\varphi }^2\right)$ distribution and defines the within-individual ICC. In this case, the three relevant ICCs can be derived under the nested exchangeable model as:

Within-period ICC (where $j=l$):

Between-period ICC:

Within-individual ICC:

In some sample size calculation methods,^16,22 an alternative correlation coefficient called the individual autocorrelation coefficient (IAC), defined as $\frac{{\sigma }_{\varphi }^2}{{{\sigma }_{\varphi }^2+\sigma }_{\epsilon }^2}$, is used to account for the correlation in repeated measures on the same individual in cohort designs.¹⁶ The fundamental difference between IAC and within-individual ICC is whether we consider the individual as being sampled from one given cluster or from any cluster. In practice, if using data from a single cluster, IAC may be a more suitable measure, and if using data from several different clusters, then within-individual ICC would be a better measure. In the literature, the cohort design version of a nested exchangeable correlation structure has been referred to as block exchangeable correlation.²⁰

Practical recommendations for estimating correlation parameters from available data

Sample size calculation methods are now available for all three models reviewed in the previous section,⁸ and have been implemented in major statistical software packages and web-based applications.^22–28 The input correlation parameters required under all three models are summarized in Table 1. In most cases, assuming an exchangeable correlation is inappropriate and anti-conservative at the design stage. Therefore, it is preferable to assume an alternative correlation structure and estimate both the within-period ICC and the CAC (in addition to the within-individual ICC or IAC in case of a cohort design). In this section, we provide practical recommendations for estimating values for these correlation parameters in advance of a planned longitudinal CRT. After obtaining estimates, readers can refer to a previously published tutorial for how to use these input parameters to determine the required sample size using any of the available sample size calculators including code and examples.⁸

Examples for illustration

In this section, we demonstrate how to estimate the correlation parameters in R, SAS and Stata using data from three real longitudinal CRTs. The code used to do the estimation is available via online code repository [https://rpubs.com/derek6561/estimateicc]. The estimation procedure we used in all three packages is restricted maximum likelihood (REML) estimation. As different software packages use different default algorithms (e.g. in our examples, SAS (PROC MIXED) uses Newton-Raphson, R (lmer) uses NLopt nonlinear optimization and Stata (mixed) uses expectation–maximization optimization) to optimize REML, it is possible that slightly different estimates may be obtained.

Discussion

In this tutorial, we considered three main correlation structures for longitudinal CRTs and provided practical recommendations and tools to allow investigators to obtain relevant estimates for a planned longitudinal CRT. Previously published trials may not have reported appropriate correlations, but data collected from such trials may be available under a data sharing agreement allowing estimates to be calculated. Furthermore, when the planned trial will use routinely collected databases for outcome assessment, investigators could potentially access these data (at individual level or in aggregate) to inform the design of their trial. We present three examples in this tutorial to demonstrate how available data can be used to estimate correlation parameters for future designs with continuous or binary outcomes, and with repeated cross-sectional or cohort designs.

For cohort designs, we considered models assuming a constant correlation in repeated measures on the same individual over time, although more complex correlation structures are available. In a proportional decay model for example, both the between-period ICC and the individual-level correlation are allowed to decay exponentially at the same or different rates.²⁰ Unfortunately, software packages that can fit proportional decay models are limited. Furthermore, the models considered in this tutorial assumed measurements are taken in discrete intervals of time; in reality, observations may be taken on a continuous basis.⁴⁹ Continuous time decay structures for trials with continuous recruitment are available,⁴⁹ but need further development. In practice, due to privacy issues, granular data allowing users to estimate such continuous time decay structures based on actual times of recruitment may not be available. Approximations may be obtained using a discrete time decay model by assuming participants arrive in equally spaced time intervals as described in Grantham et al.⁴⁹

In this manuscript, we did not consider covariates when estimating the ICCs. Currently available sample size formulae for longitudinal CRTs are based on models without baseline covariates; see, for example, the review of Li et al.⁷ In a real application, covariates may be available at the design stage. If the planned analysis includes adjustment for covariates, an adjusted ICC may be a more accurate input design parameter for sample size calculation, but if covariate information is not available at the design stage, using unadjusted ICC estimates may be appropriate and will typically yield a more conservative sample size estimate. When the planned analyses include accounting for treatment effect heterogeneity, multiple ICCs for example, the ICC for the outcome adjusting for covariates and the ICC for the covariates may be needed.^50,51 Readers may still be able to apply the methods described here to estimate each ICC and use available sample size formulae for powering the study to detect treatment effect heterogeneity.^50,52

In this tutorial, we focused on the situation where a suitable historical dataset has been identified. In the absence of suitable empirical data, researchers may need to use ‘rules of thumb’. For example, a within-period ICC of 0.05 and CAC of 0.8 may be reasonable values for continuous outcomes, although previous studies have suggested that higher values for process measures may be reasonable.³⁰ Several studies have examined patterns in databases of ICC values.^13,30–33 For binary outcomes, Campbell et al.³⁰ discussed the relationship between the ICC and prevalence, which can be helpful in choosing an appropriate ICC. Regardless of the choice of ICCs, a sensitivity analysis is strongly recommended to obtain a range of sample sizes under alternative assumptions about the ICC parameters.

There are several gaps in the literature requiring further research. First, estimating correlation parameters for longitudinal trials with binary outcomes is challenging and an under-explored area. Although obtaining estimates on the logistic scale may be straightforward, most sample size calculators based on a mixed-effects model for longitudinal CRTs require estimates on the proportions scale.⁸ Therefore, a practical strategy is to fit linear mixed-effect models to the data and treat the binary outcome as continuous; this way, one may be able to plug the resulting ICC estimates into the sample size formulae developed for continuous outcome for an approximate calculation (when the effect size of interest is the risk difference). Second, there is no consensus on how to choose the best correlation structure for longitudinal CRTs. In practice, researchers may calculate the sample sizes under both exponential decay and block/nested exchangeable models and choose the structure yielding the most conservative estimate. Third, we have focused on the mixed-effects model framework, which is most commonly used in practice and more accessible in existing software.^7,53,54 An alternative approach is generalized estimating equations (GEE).^55,56 The GEE approach can be more attractive for binary or count outcomes because the ICCs can be directly parameterized on the natural scale of the outcomes, even when a non-identity link function is used for the outcome model (that is, when the effect size of interest is a relative risk or odds ratio).^20,57,58 However, related software for implementing GEE models with longitudinal correlation structures is still being developed, and the associated computational challenges with large cluster sizes have only recently been addressed for cross-sectional designs.³⁵ Fourth, as discussed earlier, future work is required to determine guidelines for minimum required numbers of clusters and cluster period sizes to yield reliable estimates for ICC parameters in longitudinal CRT designs. Additionally, methods for obtaining confidence intervals around estimated within-period ICCs and CACs under nested/block exchangeable and exponential decay correlation structures, under a mixed-effects model framework, require further investigation. Finally, attrition may be present in a cohort design and may be informative, i.e. the missing value in a particular period may depend on the intervention and the outcome in a previous period, as well as other unmeasured factors. Even under non-informative missingness, failure to account for the complex nature of the correlations in the imputation model may lead to incorrect estimation of the correlation parameters based on imputation-completed data. Further work is required to develop appropriate imputation methods for longitudinal CRTs.

Ethics approval

Not applicable.

Supplementary Material

Data availability

There are no new data associated with this article. Some code for data manipulation and model fitting can be found at [https://rpubs.com/derek6561/estimateicc]. A RShiny app used to estimate ICC is available at [https://douyang.shinyapps.io/estimateICC].

Supplementary data

Supplementary data are available at IJE online.

Author contributions

Y.O. led the writing of the manuscript and developed the idea with M.T. M.T. conceived the ideas and co-led the writing. K.H. and F.L. participated in the development of the idea and the method. All authors contributed to writing, drafting and editing the manuscript.

Funding

M.T. and F.L. are supported by the National Institute of Aging (NIA) of the National Institutes of Health (NIH) under Award Number U54AG063546, which funds NIA Imbedded Pragmatic Alzheimer's Disease and AD-Related Dementias Clinical Trials Collaboratory (NIA IMPACT Collaboratory). Y.O. is funded by the Health System Impact Fellowship, which is supported by the Canadian Institutes of Health Research (CIHR). F.L. is supported through a Patient-Centered Outcomes Research Institute^® (PCORI^® Award ME-2020C3-21072). The statements presented in this article are solely the responsibility of the authors and do not necessarily represent the views of NIH, nor PCORI^®, its Board of Governors or Methodology Committee. K.H. is funded by an NIHR Senior Research Fellowship SRF-2017–10-002.

Conflict of interest

None declared.