A readily available improvement over method of moments for intra-cluster correlation estimation in the context of cluster randomized trials and fitting a GEE

Abstract

Background/aims

Cluster randomized trials are popular in health-related research due to the need or desire to randomize clusters of subjects to different trial arms as opposed to randomizing each subject individually. As outcomes from subjects within the same cluster tend to be more alike than outcomes from subjects within other clusters, an exchangeable correlation arises that is measured via the intra-cluster correlation coefficient. Intra-cluster correlation coefficient estimation is especially important due to the increasing awareness of the need to publish such values from studies in order to help guide the design of future cluster randomized trials. Therefore, numerous methods have been proposed to accurately estimate the intra-cluster correlation coefficient, with much attention given to binary outcomes. As marginal models are often of interest, we focus on intra-cluster correlation coefficient estimation in the context of fitting such a model with binary outcomes using generalized estimating equations. Traditionally, intra-cluster correlation coefficient estimation with generalized estimating equations has been based on the method of moments, although such estimators can be negatively biased. Furthermore, alternative estimators that work well, such as the analysis of variance estimator, are not as readily applicable in the context of practical data analyses with generalized estimating equations. Therefore, in this article we assess, in terms of bias, the readily available residual pseudo-likelihood approach to intra-cluster correlation coefficient estimation with the GLIMMIX procedure of SAS (SAS Institute, Cary, NC). Furthermore, we study a possible corresponding approach to confidence interval construction for the intra-cluster correlation coefficient.

Methods

We utilize a simulation study and application example to assess bias in intra-cluster correlation coefficient estimates obtained from GLIMMIX using residual pseudo-likelihood. This estimator is contrasted with method of moments and analysis of variance estimators which are standards of comparison. The approach to confidence interval construction is assessed by examining coverage probabilities.

Results

Overall, the residual pseudo-likelihood estimator performs very well. It has considerably less bias than moment estimators, which are its competitor for general generalized estimating equation–based analyses, and therefore, it is a major improvement in practice. Furthermore, it works almost as well as analysis of variance estimators when they are applicable. Confidence intervals have near-nominal coverage when the intra-cluster correlation coefficient estimate has negligible bias.

Conclusion

Our results show that the residual pseudo-likelihood estimator is a good option for intra-cluster correlation coefficient estimation when conducting a generalized estimating equation–based analysis of binary outcome data arising from cluster randomized trials. The estimator is practical in that it is simply a result from fitting a marginal model with GLIMMIX, and a confidence interval can be easily obtained. An additional advantage is that, unlike most other options for performing generalized estimating equation–based analyses, GLIMMIX provides analysts the option to utilize small-sample adjustments that ensure valid inference.

Keywords

Empirical covariance matrix generalized estimating equations group randomized trial method of moments residual pseudo-likelihood small-sample corrections

Introduction

Increasingly in health-related research, cluster randomized trials, or group randomized trials, are used to randomize entire clusters, or groups, of subjects to trial arms as opposed to individually randomizing subjects.^1–6 Such trials are conducted, for instance, because of the need or desire to apply an intervention to the entire cluster due to natural grouping or to avoid contamination. Popular examples of possible clusters include, but are not limited to, communities, schools, and medical care practices.

The number of clusters and their sizes depend on the setting. For instance, Van Breukelen and Candel⁷ summarize results from reviews on primary care cluster randomized trials reported by Adams et al.⁸ and Eldridge et al.,⁹ with the median number of clusters ranging from 31 to 41, and a median of the mean number of subjects per cluster equaling 23⁸ or 32.⁹ Alternatively, in many cluster randomized trials, the number of clusters is small, whereas the size of each cluster can be large. For example, Roetzheim et al.^10,11 and Lee et al.¹² report results from a repeated cross-sectional cluster randomized trial with focus on three cancer screening outcomes, which were subject-level binary outcomes of whether or not a given subject had received a mammography, a papanicolaou (pap) smear, or a fecal occult blood test. The intervention and control conditions were only randomized to four clinics each, and data were obtained from 51 to 127 subjects per clinic, depending on the time period and outcome type.¹²

In this article, we focus on situations in which subjects contribute binary outcomes; for example, an indicator for having been screened for a particular type of cancer. Because outcomes from subjects within the same cluster tend to be more alike than outcomes from subjects within other clusters, for example, due to differences in clinical practices, an exchangeable correlation arises that is measured via the intra-cluster, or intra-class, correlation coefficient (ICC). For more detail, a review on different definitions for the ICC can be found in Eldridge et al.¹³ In this article, we focus on ICCs corresponding to the proportions scale.

The ICC is important because it must be appropriately accounted for in the design of cluster randomized trials and during the statistical analysis of data resulting from such a trial. According to the extension of the Consolidated Standards of Reporting Trials (CONSORT) statement to cluster randomized trials,¹⁴ ICC estimates should be reported, which are very helpful for planning future cluster randomized trials.⁶ For example, Hade et al.¹⁵ reported 138 unadjusted and adjusted ICC estimates, from 12 different studies, corresponding to different types of cancer screening, and Taljaard et al.¹⁶ reported 86 estimated ICC values for variables corresponding to maternal and perinatal health. In addition, there has been an increase in reporting of ICCs in primary care settings.¹⁷ Furthermore, databases of ICC values have even been created for certain types of studies; for example, a surgical-based ICC database.^6,18

Due to the importance of reporting reliably estimated ICC values, numerous ICC estimation methods have been developed and studied against each other. For summaries, see, for instance, the references cited here.^6,19,20 Although multiple estimators have been shown to work well, such as the well-known ANOVA estimator^19–23 which was used by Hade et al.,¹⁵ to our knowledge many of these estimators cannot be used in the context of practical data analyses in which marginal models, fit using generalized estimating equations (GEE),²⁴ do or do not adjust for covariates.

With GEE, the exchangeable correlation, or the ICC in our context, has traditionally been estimated using the method of moments. This is the case, for instance, in the popular R package geepack^25–28 as well as SAS GENMOD.²⁹ However, this approach is well known to produce negatively biased ICC estimates, as demonstrated in, for instance, Wu et al.²⁰ and later in our simulation study. The resulting problem is that the use of such negatively biased estimates to plan future cluster randomized trials could result in the calculation of a sample size (number of clusters or the sizes of the clusters) that is too small or alternatively an over-calculation of power for a fixed number of clusters and their sizes.

Due to the limitations of well-known and routinely used estimators, in this article we demonstrate the ability of SAS GLIMMIX²⁹ to estimate the ICC, on the proportions scale, in a GEE-based analysis. In short, residual pseudo-likelihood is used via SAS GLIMMIX to estimate the ICC as opposed to the method of moments. Although our main focus in this article is demonstrating the ability of this procedure to estimate the ICC, our desire to demonstrate how well GLIMMIX works in this regard is also due to this procedure’s ability to implement small-sample bias corrections to empirical standard error estimates and the allowance of user-defined degrees of freedom in order to ensure valid inference. We discuss this issue in detail later in the “Discussion” section and in an online supplementary appendix.

We note that the construction of confidence intervals for the ICC has received a considerable amount of attention in the cluster randomized trial literature. For example, see Chakrabory and Hossain³⁰ and Braschel et al.³¹ Methods for confidence interval construction have been established for the ANOVA and methods of moments estimators,^30,32 but not for the residual pseudo-likelihood estimator. Therefore, we study a general confidence interval construction approach of Demetrashvili et al.³³ that is naturally applicable for use when estimating the ICC with residual pseudo-likelihood.

Next, we briefly give details on GEE-based analyses and discuss ICC estimation with the method of moments and residual pseudo-likelihood. Furthermore, as a standard of comparison, we also describe ANOVA estimators due to their popularity, their proven ability to estimate the ICC well, as well as the capability to use this type of estimator even when adjusting for cluster-level covariates in the marginal model. Via a simulation study, we then demonstrate the utility of the residual pseudo-likelihood method of GLIMMIX in terms of estimating the ICC relative to ANOVA-based and method of moment estimators. Furthermore, we study the confidence interval construction approach of Demetrashvili et al.³³ to determine how well it works when the ICC is estimated with residual pseudo-likelihood. The ICC estimation approaches are then contrasted via application to analyses of data from the cancer screening study of Roetzheim et al.^10,11 Finally, we provide concluding remarks, including a discussion on the additional practical advantages of using GLIMMIX.

Methods

Notation

For simplicity, we assume two trial arms, each consisting of $N / 2$ clusters. Outcomes from $n_{i}$ subjects are obtained from cluster i, $i = 1, \dots, N$ . The corresponding outcome vector is given by $Y_{i} = [Y_{i 1}, \dots, Y_{i n_{i}}]^{T}$ , with marginal mean denoted by $E (Y_{i}) = μ_{i}$ . At the subject level, we have the marginal model $f (μ_{ij}) = x_{ij}^{T} β$ , $j = 1, \dots, n_{i}$ ; $i = 1, \dots, N$ , in which $x_{ij} = [1, x_{1 ij}, \dots, x_{(p - 1) ij}]^{T}$ is the vector of covariate values for the $j th$ subject in cluster i and the p marginal regression parameters are $β = [β_{0}, β_{1}, \dots, β_{p - 1}]^{T}$ . In our context of binary outcomes, we assume f is the logit function and $μ_{ij} = π_{ij}$ is a marginal probability.

The working covariance structure of $Y_{i}$ is given by $V_{i} = A_{i}^{1 / 2} R_{i} (ρ) A_{i}^{1 / 2}$ , $i = 1, \dots, N$ . Here, $A_{i} = diag [π_{i 1} (1 - π_{i 1}), \dots, π_{i n_{i}} (1 - π_{i n_{i}})]$ is a diagonal matrix of marginal variances, $ρ$ is the ICC which is typically non-negative for cluster randomized trials,¹³ and $R_{i}$ is an $n_{i} \times n_{i}$ working correlation matrix with 1 on the diagonal and $ρ$ elsewhere. Although the ICC value may vary across clusters, for simplicity, we currently assume a constant value as is often done in practice.³⁴

Marginal modeling via GEE

To fit a marginal model utilizing GEE, the estimate for $β$ , $\hat{β}$ , is found by iteratively solving

\sum_{i 1 =}^{N} D_{i}^{T} A_{i}^{- 1 / 2} R_{i}^{- 1} (\hat{ρ}) A_{i}^{- 1 / 2} (Y_{i} - μ_{i}) = 0

in which $D_{i} = \partial μ_{i} / \partial β^{T}$ and $\hat{ρ}$ is the estimated ICC. We note that when there are only cluster-level covariates, that is, $x_{i 1} = \dots = x_{i n_{i}} = x_{i}$ and thus $π_{i 1} = \dots = π_{i n_{i}} = π_{i}$ , these equations can be simplified to the quasi-likelihood³⁵ form given by

\sum_{i = 1}^{N} x_{i} \frac{Y_{i} - n_{i} π_{i}}{1 + (n_{i} - 1) \hat{ρ}} = 0

in which $Y_{i} = \sum_{j = 1}^{n_{i}} Y_{ij}$ is the number of events out of $n_{i}$ subjects in the ith cluster. In order to obtain consistent regression parameter estimates, the marginal mean must be correctly modeled such that the estimating equations are unbiased. Although we assume a logit link function as is often done in practice, in general the true and working link functions are not required to be the logit function.

ICC estimation

With our focus on GEE, the ICC corresponds to a common exchangeable correlation parameter which, in general, is given by $Corr (Y_{ij}, Y_{ik}) = Cov (Y_{ij}, Y_{ik}) / \sqrt{π_{ij} (1 - π_{ij}) π_{ik} (1 - π_{ik})}$ . Due to the use of marginal variances, focus is on the proportions scale. For instance, assume there are only cluster-level covariates, and let $p_{i}$ denote the true probability of an event for any subject in cluster i, $i = 1, \dots, N$ , that has a marginal mean of $E [p_{i}] = π_{i}$ . The constant ICC is then given as $ρ = var (p_{i}) / [π_{i} (1 - π_{i})]$ ,¹³ $i = 1, \dots, N$ , which is a formula we later utilize to generate data in our simulation study.

The traditional approach to obtaining $\hat{ρ}$ with GEE is the use of method of moments. Although slight modifications can exist, a common estimator used by SAS [29] for binary outcomes is given by

\hat{ρ} = {[\hat{ϕ} (0.5 \sum_{i = 1}^{N} n_{i} (n_{i} - 1) - p)]}^{- 1} \sum_{i = 1}^{N} \sum_{j = 1}^{n_{i} - 1} \sum_{k = j + 1}^{n_{i}} e_{ij} e_{ik}

(1)

where $e_{ij} = (Y_{ij} - π_{ij}) / \sqrt{π_{ij} (1 - π_{ij})}$ and $\hat{ϕ} = \sum_{i = 1}^{N} \sum_{j = 1}^{n_{i}} e_{ij}^{2} / (\sum_{i = 1}^{N} n_{i} - p)$ . Alternatively, with an applied emphasis given to cluster randomized trials, Lu et al.³⁶ proposed a general method of moments approach that incorporates bias corrections within estimating equations in order to reduce bias. Due to the complexity of this approach, we refer the reader to this article or Preisser et al.³² for more detail. We also note that the subtraction of p within equation (1) is similar in nature to this bias correction.

Due to the popularity and utility of the ANOVA estimator,^19–23 as well as the possibility to use it with GEE when no subject-level covariates are incorporated within the marginal model, we later utilize it as a standard of comparison in our simulation study. When a trial arm indicator is the only covariate in the marginal model, the estimate for $ρ$ using the ANOVA method is given by

\hat{ρ} = \frac{MSB - MSE}{MSB + (K - 1) MSE}

(2)

Let $m_{1} = \sum_{i = 1}^{N / 2} n_{i}$ and $m_{2} = \sum_{i = N / 2 + 1}^{N} n_{i}$ be the total number of subjects in the first and second trial arms, respectively. In the above ANOVA formula³⁷

\begin{matrix} MSB = \frac{1}{N - 2} \sum_{i = 1}^{N} n_{i} {(\frac{Y_{i}}{n_{i}} - {\hat{π}}_{i})}^{2} \\ MSE = \frac{1}{\sum_{i = 1}^{N} n_{i} - N} \sum_{i = 1}^{N} Y_{i} (1 - \frac{Y_{i}}{n_{i}}) \\ K = \frac{1}{N - 2} [\sum_{i = 1}^{N} n_{i} - \sum_{i = 1}^{N / 2} \frac{n_{i}^{2}}{m_{1}} - \sum_{i = N / 2 + 1}^{N} \frac{n_{i}^{2}}{m_{2}}] \end{matrix}

When other cluster-level covariates are included in the marginal model, an alternative formula can be derived based on the ANOVA method³⁸ that still has the form given in equation (2), although formulas differ as follows

\begin{matrix} MSB = \frac{1}{N} \sum_{i = 1}^{N} \frac{1}{n_{i}} \frac{{(Y_{i} - n_{i} {\hat{π}}_{i})}^{2}}{1 - h_{i}} \\ MSE = \frac{1}{\sum_{i = 1}^{N} n_{i} - N} \sum_{i = 1}^{N} Y_{i} (1 - \frac{Y_{i}}{n_{i}}) \\ K = \frac{1}{ρ} \frac{\sum_{i = 1}^{N} n_{i} - N}{B - A} [\frac{1}{N} [A + ρ (B - A)]] + 1 - \frac{1}{ρ} \\ A = \sum_{i = 1}^{N} {\hat{π}}_{i} (1 - {\hat{π}}_{i}), B = \sum_{i = 1}^{N} n_{i} {\hat{π}}_{i} (1 - {\hat{π}}_{i}) \end{matrix}

and $h_{i}$ is the ith diagonal element of $H = W^{1 / 2} X (X' WX)^{- 1} X' W^{1 / 2}$ , $X' = [x_{1}, \dots, x_{N}]$ , $W = diag (n_{i} π_{i} (1 - π_{i}) / [1 + (n_{i} - 1) ρ])$ . K is both a function of $ρ$ and is used to estimate $ρ$ , and therefore an iterative procedure should be used.³⁸ Furthermore, as opposed to the use of $(N - 2)^{- 1}$ in the prior version of the MSB, the latter version incorporates leverage values, $h_{i}$ , $i = 1, \dots, N$ , because $E [(Y_{i} - n_{i} {\hat{π}}_{i})^{2}] \approx n_{i} π_{i} (1 - π_{i}) [1 + (n_{i} - 1) ρ] (1 - h_{i}) = (1 - h_{i}) Var (Y_{i})$ .^39–41 We note that the two different ANOVA formulas provide equivalent estimates if trial arm is the only covariate in the model and clusters are all of the same size, although distinct estimates can occur when clusters vary in size.

In this article, we propose the use of residual pseudo-likelihood via the GLIMMIX procedure in SAS.²⁹ In short, Taylor series expansions are implemented such that pseudo, or transformed, response data are created to be used in a weighted linear mixed model, and therefore SAS documentation refers to the corresponding marginal model as a “GEE-type” model.²⁹ Due to the use of this weighting, the estimated compound symmetry and residual covariance parameters, ${\hat{σ}}_{CS}^{2}$ and ${\hat{σ}}_{Residual}^{2}$ , respectively, that SAS produces are common. Therefore, the estimated common ICC is given by

\hat{ρ} = \frac{{\hat{σ}}_{CS}^{2}}{{\hat{σ}}_{CS}^{2} + {\hat{σ}}_{Residual}^{2}}

(3)

As most cluster randomized trials do not involve a large number of clusters, we suggest the use of the default restricted maximum likelihood⁴² for covariance parameter estimation, that is, estimation of the ICC. Details on residual pseudo-likelihood are provided in the online supplementary appendix, and in-depth information is provided by Fitzmaurice et al.⁴² and SAS’s User’s Guide.²⁹

We note that methods for constructing confidence intervals corresponding to ICC estimates from the residual pseudo-likelihood approach are not established. However, Demetrashvili et al.³³ proposed an approach to obtain confidence intervals with approximately correct coverage in the context of variance component models. Therefore, due to the use of restricted maximum likelihood, and based on the form for the ICC estimate given in equation (3), the approach of Demetrashvili et al.³³ naturally can be applied when utilizing residual pseudo-likelihood to estimate the ICC. For details, we refer the reader to Demetrashvili et al.³³ and the SAS code that is available in Supplemental Material.

Simulation study

We conduct a simulation study to demonstrate how well the residual pseudo-likelihood approach works relative to the method of moments and ANOVA methods in terms of estimating the ICC. Specifically, we compare empirical biases. We note that for each of these methods, if a negative ICC estimate is produced, we use an ICC value of 0 as is often done in practice.⁶ Furthermore, we assess the approach of Demetrashvili et al.³³ to constructing a confidence interval for the ICC based on the residual pseudo-likelihood estimation method. Specifically, we compare empirical coverage probabilities of these 95% confidence intervals with the nominal 0.950 value.

For each setting, 10,000 datasets are generated from either an unadjusted or adjusted marginal logistic regression model. The unadjusted model is given by logit $(π_{i}) = β_{0} + β_{1} I_{i} (Intervention)$ with corresponding results given in Tables 1 and 2, and the adjusted model, logit $(π_{i}) = β_{0} + β_{1} I_{i} (Intervention) + β_{2} X_{i}$ , includes an additional cluster-level covariate, $X_{i}$ , with corresponding results given in Table 3. Here, $I_{i} (Intervention)$ is an indicator for the ith cluster having been randomized to intervention (1 – Intervention, 0 – Control).

Table 1.

Empirical biases of estimated intra-cluster correlation coefficients (ICCs) from 10,000 datasets generated for each setting (row) from the unadjusted marginal logistic regression model with given coefficient of variation (CV) of cluster sizes (mean size is 100), true ICC values, number of clusters per trial arm $(N / 2)$ , and marginal probabilities $(π)$ for each arm.

CV	ICC	$N / 2$	$π$	ANOVA 1	MM ANOVA 2	MM Common	Corrected	RPL	RPL CI CP
			0.1	0.000	0.000	−0.003	−0.003	0.000	0.962
		5	0.2	0.000	0.000	−0.003	−0.002	0.000	0.962
	0.01		0.5	0.000	0.000	−0.003	−0.002	0.000	0.962
			0.1	0.000	0.000	−0.002	−0.002	0.000	0.936
		10	0.2	0.000	0.000	−0.002	−0.001	0.000	0.917
			0.5	0.000	0.000	−0.002	−0.001	0.000	0.927
0
			0.1	−0.0007	−0.0007	−0.030	−0.026	−0.013	0.819
		5	0.2	−0.002	−0.002	−0.024	−0.022	−0.005	0.887
	0.10		0.5	0.000	0.000	−0.020	−0.018	0.000	0.905
			0.1	−0.003	−0.003	−0.016	−0.013	−0.006	0.847
		10	0.2	0.000	0.000	−0.011	−0.010	−0.002	0.914
			0.5	0.001	0.001	−0.009	−0.008	0.001	0.938
			0.1	0.002	0.000	−0.004	−0.005	0.001	0.923
		5	0.2	0.001	0.000	−0.004	0.005	0.000	0.933
	0.01		0.5	0.001	0.000	−0.004	−0.004	0.000	0.942
			0.1	0.000	0.000	−0.002	−0.003	0.000	0.914
		10	0.2	0.000	0.000	−0.002	−0.003	0.000	0.922
			0.5	0.000	0.000	−0.002	−0.002	0.000	0.924
1
			0.1	0.002	−0.011	−0.037	−0.046	−0.010	0.812
		5	0.2	0.008	−0.005	−0.029	−0.041	−0.003	0.869
	0.10		0.5	0.008	−0.004	−0.025	−0.048	−0.002	0.893
			0.1	−0.001	−0.006	−0.022	−0.028	−0.006	0.826
		10	0.2	0.003	−0.002	−0.015	−0.022	−0.002	0.900
			0.5	0.004	−0.001	−0.012	−0.020	0.000	0.930

Empirical coverage probabilities (CPs) of 95% confidence intervals (CIs) corresponding to the residual pseudo-likelihood estimator and construction approach of Demetrashvili et al.³³ are also presented in the last column.

ANOVA 1: The traditional ANOVA estimator that cannot account for covariate adjustment.

ANOVA 2: The ANOVA estimator that can account for covariate adjustment and utilizes leverage values.

MM Common: The common method of moments estimator given in equation (1).

MM Corrected: The corrected method of moments estimator of Lu et al.³⁶

RPL: residual pseudo-likelihood from GLIMMIX.

Table 2.

Empirical biases of estimated intra-cluster correlation coefficients (ICCs) from 10,000 datasets generated for each setting (row) from the unadjusted marginal logistic regression model with given coefficient of variation (CV) of cluster sizes (mean size is 30), true ICC values, number of clusters per trial arm $(N / 2)$ , and marginal probabilities $(π)$ for each arm.

CV	ICC	$N / 2$	$π$	ANOVA 1	ANOVA 2	MM Common	MM Corrected	RPL	RPL CI CP
			0.1	0.001	0.001	−0.002	−0.001	0.001	0.945
		15	0.2	0.001	0.001	−0.001	−0.001	0.001	0.934
			0.5	0.001	0.001	−0.001	0.000	0.001	0.951
	0.01
			0.1	0.001	0.001	−0.001	−0.001	0.000	0.937
		20	0.2	0.001	0.001	−0.001	0.000	0.001	0.933
			0.5	0.001	0.001	−0.001	0.000	0.001	0.939
0
			0.1	−0.003	−0.003	−0.013	−0.010	−0.006	0.858
		15	0.2	−0.001	−0.001	−0.010	−0.008	−0.002	0.916
			0.5	0.000	0.000	−0.008	−0.005	0.000	0.951
	0.10
			0.1	−0.002	−0.002	−0.010	−0.008	−0.004	0.854
		20	0.2	−0.001	−0.001	−0.007	−0.006	−0.001	0.924
			0.5	0.000	0.000	−0.006	−0.004	0.000	0.952
			0.1	0.001	0.001	−0.002	−0.002	0.001	0.939
		15	0.2	0.001	0.001	−0.002	−0.002	0.001	0.945
			0.5	0.001	0.001	−0.002	−0.001	0.001	0.950
	0.01
			0.1	0.001	0.001	−0.002	−0.002	0.000	0.924
		20	0.2	0.001	0.001	−0.002	−0.001	0.000	0.928
			0.5	0.001	0.000	−0.002	−0.001	0.000	0.937
1
			0.1	−0.001	−0.004	−0.017	−0.020	−0.005	0.849
		15	0.2	0.000	−0.002	−0.012	−0.016	−0.002	0.913
			0.5	0.002	−0.001	−0.01	−0.014	−0.001	0.938
	0.10
			0.1	−0.001	−0.003	−0.013	−0.015	−0.003	0.851
		20	0.2	0.001	−0.001	−0.009	−0.013	−0.001	0.914
			0.5	0.001	−0.001	−0.007	−0.010	0.000	0.943

ANOVA 1: The traditional ANOVA estimator that cannot account for covariate adjustment.

ANOVA 2: The ANOVA estimator that can account for covariate adjustment and utilizes leverage values.

MM Common: The common method of moments estimator given in equation (1).

MM Corrected: The corrected method of moments estimator of Lu et al.³⁶

RPL: Residual pseudo-likelihood from GLIMMIX.

Table 3.

Empirical biases of estimated intra-cluster correlation coefficients (ICCs) from 10,000 datasets generated for each setting (row) from the adjusted marginal logistic regression model with given mean cluster size $(\bar{n})$ , coefficient of variation (CV) of cluster sizes, true ICC values, and number of clusters per trial arm $(N / 2)$ .

$\bar{n}$	CV	ICC	$N / 2$	ANOVA 2	MM Common	MM Corrected	RPL	RPL CI CP
			5	0.000	−0.005	−0.004	0.000	0.961
		0.01
			10	0.000	−0.003	−0.002	0.000	0.929
	0
			5	−0.003	−0.035	−0.032	−0.007	0.873
		0.10
			10	−0.001	−0.017	−0.016	−0.003	0.904
100
			5	0.001	−0.006	−0.006	0.001	0.923
		0.01
			10	0.000	−0.003	−0.004	0.000	0.913
	1
			5	−0.009	−0.043	−0.058	−0.006	0.865
		0.10
			10	−0.003	−0.021	−0.031	−0.002	0.903
			15	0.001	−0.003	−0.001	0.001	0.941
		0.01
			20	0.001	−0.002	−0.001	0.001	0.938
	0
			15	−0.001	−0.014	−0.011	−0.003	0.921
		0.10
			20	−0.001	−0.011	−0.009	−0.002	0.924
30
			15	0.001	−0.003	−0.003	0.001	0.940
		0.01
			20	0.001	−0.003	−0.002	0.000	0.929
	1
			15	−0.002	−0.017	−0.023	−0.002	0.918
		0.10
			20	−0.001	−0.012	−0.017	−0.002	0.920

ANOVA 2: The ANOVA estimator that can account for covariate adjustment and utilizes leverage values.

MM Common: The common method of moments estimator given in equation (1).

MM Corrected: The corrected method of moments estimator of Lu et al.³⁶

RPL: Residual pseudo-likelihood from GLIMMIX.

We only present results from the unadjusted model settings in which $β_{1} = 0$ . Therefore, both arms have the same marginal probability, $π$ , for the unadjusted model. We use true marginal probabilities of 0.10, 0.20, and 0.50 such that we study ICC estimation relatively near the boundary of the parameter space as well as at its center. Not allowing an intervention effect allows us to clearly study results at these different places in the parameter space. For the adjusted model, not only do we allow an intervention effect, we more generally allow marginal probabilities to greatly vary across clusters. Specifically, we utilize the model in Westgate⁴³ such that $β = [- 1.5, 1, 0.5]^{T}$ and $X_{i}$ is independently generated from $Uniform (- 1.4, 1)$ , allowing marginal probabilities to vary anywhere from 0.10 to 0.50. We note that generating $X_{i}$ from this distribution is meant to represent an analysis of covariance model as detailed later in our application example.

Similar to our application example and the numbers of clusters and their sizes reported in Adams et al.⁸ and Eldridge et al.,⁹ we focus on general cluster randomized trial settings with either a small number of large clusters or a medium number of smaller clusters. Specifically, the former (latter) settings comprise either 5 or 10 (15 or 20) clusters in each trial arm, with a mean cluster size of 100 (30). Due to computational time, we do not study larger clusters or a larger number of clusters. However, we do note that as the number of clusters increases, the degree of bias in the ICC estimation approaches decreases. We utilize true ICC values of 0.01 and 0.10, as most estimated ICCs reported from such cluster randomized trials tend to not be large.^6,8,9,17 Furthermore, the coefficients of variation of cluster sizes were 0 and 1, with sizes generated from a normal distribution and then rounded to the nearest integer.^44,45 More details are provided in the online supplementary appendix.

Data are generated from the given marginal models via the use of a beta-binomial distribution. R version 3.4.3²⁵ and SAS²⁹ are used to conduct the study. Details on data generation and software use are provided in the online supplementary appendix.

Application

To illustrate the different ICC estimation methods assessed in our simulation study, we analyze data on females from the Roetzheim et al.^10,11 and Lee et al.¹² cluster randomized trial. As a reminder, there are eight total clinics, or clusters, such that four are in each of the control and intervention arms, and three different outcomes are subject-level binary indicators of whether or not a given subject had received a mammography, a papanicolaou smear, or a fecal occult blood test. We utilize data from baseline and 2 years post baseline. Although this study utilizes a repeated cross-sectional design, some subjects are observed at both points in time simply by chance. The number of subjects in each clinic contributing data with respect to a papanicolaou smear, mammography, and fecal occult blood test at 2 years post baseline ranges from 66 to 85, 111 to 122, and 112 to 126, respectively. The proportion of women receiving a papanicolaou smear, a mammography, and a fecal occult blood test ranges from 0.38 to 0.55, 0.50 to 0.71, and 0.04 to 0.22, respectively, in the control arm and 0.15 to 0.71, 0.45 to 0.77, and 0.14 to 0.37, respectively, in the intervention arm.

For each outcome, we fit two marginal logistic regression models:

Unadjusted: logit $(π_{2 i}) = β_{0} + β_{1} I_{i} (Intervention)$ ;

Adjusted: logit $(π_{2 i}) = β_{0} + β_{1} I_{i} (Intervention) + β_{2}$ logit $({\hat{p}}_{0 i})$ .

Within these models, $π_{2 i}$ is the marginal probability for cluster i at 2 years post baseline, and logit $({\hat{p}}_{0 i})$ is the natural log of the baseline sample odds of having the given outcome in cluster i. We note that the cluster-level adjustment covariate, logit $({\hat{p}}_{0 i})$ , is utilized because its use results in an analysis of covariance model that has previously been shown to be preferable.⁴⁶

Results

Simulation study

As expected, the ANOVA estimators perform very well overall, as they often provide unbiased ICC estimates or estimates with little bias. Furthermore, both types of moment estimators typically provide negatively biased ICC estimates. Overall, the residual pseudo-likelihood estimator performs very well. It notably outperforms the moment estimators that are traditionally used with GEE analyses in practice. Furthermore, it tends to work almost as well as the ANOVA estimators. When the true ICC value is 0.01, the residual pseudo-likelihood estimator yields unbiased results with the exception in some settings of the slight simulated bias of only 0.001. Although there are a few settings in which small negative bias exists with the residual pseudo-likelihood estimator when the true ICC value is 0.10, this bias is not nearly as great as the bias in the moment estimators, and the bias in the ANOVA estimators is only slightly smaller. Furthermore, there are settings in which the residual pseudo-likelihood estimator performs best.

For all estimators, bias tends to decrease as the number of clusters increases. Also, the closer marginal probabilities are to the edge of the parameter space, that is, the further probabilities are from 0.5, the greater the potential bias in any of these estimators. We note that this bias is more notable with the larger true ICC value of 0.1. We also note that allowing marginal probabilities to vary across clusters, as opposed to having equal marginal probabilities across clusters, does not impact our findings.

The confidence interval construction approach of Demetrashvili et al.³³ tends to work well, but not perfectly, when applied with the residual pseudo-likelihood estimator. When bias is negligible, for example, $\leq 0.001$ , 95% confidence interval coverage probabilities are often relatively near-nominal. Specifically, for such negligible bias, the majority of empirical coverage probabilities are within 0.02 of the desired 0.95 value, but the lowest observed coverage probability is as low as 0.905. As expected with confidence intervals in general, when the residual pseudo-likelihood estimator is biased, coverage probabilities are subnominal. In short, the greater the bias, the lower the coverage probability.

Application

ICC estimates are provided in Table 4, corresponding to each outcome type, model, and estimation method. In general, results are similar to findings in the simulation study. The method of moments estimators yield the smallest ICC values. Alternatively, the residual pseudo-likelihood estimates are similar to the estimates produced by the ANOVA methods. Based on the results from our simulation study, this is the result we are both expecting and hoping for.

Table 4.

Estimates of the intra-cluster correlation coefficient from analyses of the cancer screening data.

Covariate adjustment?	Outcome	ANOVA 1	ANOVA 2	MM Common	MM Corrected	RPL
	Pap Smear	0.126	0.127	0.096	0.099	0.123
No	Mammography	0.063	0.063	0.047	0.049	0.063
	FOBT	0.046	0.046	0.034	0.034	0.048
	Pap Smear	–	0.089	0.047	0.051	0.084
Yes	Mammography	–	0.010	0.003	0.009	0.009
	FOBT	–	0.056	0.034	0.034	0.058

Unadjusted marginal logistic regression model: logit $(π_{2 i}) = β_{0} + β_{1} I_{i} (Intervention)$ .

Adjusted marginal logistic regression model: logit $(π_{2 i}) = β_{0} + β_{1} I_{i} (Intervention) + β_{2} logit ({\hat{p}}_{0 i})$ .

$π_{2 i}$ : The marginal probability for cluster i of having the given outcome at 2 years post baseline.

Logit $({\hat{p}}_{0 i})$ : The natural log of the baseline sample odds of having the given outcome in cluster i.

ANOVA 1: The traditional ANOVA estimator that cannot account for covariate adjustment.

ANOVA 2: The ANOVA estimator that can account for covariate adjustment and utilizes leverage values.

MM Common: The common method of moments estimator given in equation (1).

MM Corrected: The corrected method of moments estimator of Lu et al.³⁶

RPL: Residual pseudo-likelihood from GLIMMIX.

FOBT: Fecal occult blood test.

Discussion

In this article, we demonstrate that ICC estimation via residual pseudo-likelihood using the GLIMMIX procedure in SAS works very well. Specifically, it results in either unbiased estimates or estimates with small bias. Furthermore, any bias is notably smaller in magnitude than the negative bias resulting from the method of moments approaches. Furthermore, residual pseudo-likelihood estimates are often comparable to ANOVA estimates.

The construction of confidence intervals for the ICC has also received a considerable amount of attention in the cluster randomized trial literature. Therefore, we studied the applicability of the construction approach of Demetrashvili et al.³³ in our context of estimating the ICC with residual pseudo-likelihood. In our simulation study, this approach is shown to often, but not always, produce near-nominal coverage probabilities as long as the bias in the residual pseudo-likelihood estimator is negligible. Therefore, this approach is recommended for use in practice, although it should be kept in mind that coverage may not be perfect. This is especially the case when notable bias exists with the estimator, although in practice it will not be known if any bias is present. We also note that this is a limitation with any ICC estimator and confidence interval construction approach.

In addition to providing better ICC estimates than other approaches that are commonly utilized with GEE in practice, GLIMMIX has additional practical advantages. First, user-defined degrees of freedom and small-sample corrections to empirical standard errors can be employed. Second, the ICC structure can be flexibly modeled.⁴² A discussion on the importance of these advantages, with details, is provided in an online supplementary appendix and is based on the references cited here.^{20,32,34,36,38,40–45,47–58}

In Supplemental Material, we provide SAS code to demonstrate how to conduct a GEE-type analysis with GLIMMIX. We focus on ICC estimation via residual pseudo-likelihood, and we provide code to calculate a confidence interval for the ICC based on the method of Demetrashvili et al.³³ Furthermore, we show how to utilize corrections to the empirical sandwich estimator and user-defined degrees of freedom. Furthermore, we demonstrate how to obtain a separate ICC estimate with corresponding confidence interval for each trial arm.

A limitation of estimating a separate ICC for each trial arm is that fewer clusters will be utilized to estimate each ICC, and therefore estimation will be less precise. Therefore, estimates for a common ICC as well as separate ICCs for each trial arm could be reported. Furthermore, both crude and adjusted estimates can be reported, as in Hade et al.^6,15 We note that ICC estimates will often, but not always, be smaller when adjusting for covariates in the regression model.^8,15,17,59 For example, this is seen with both pap smears and mammographies in our application example. Studies should also report the method for obtaining ICC estimates in order for researchers to judge if the ICC was estimated using a reliable approach; for example, an estimate from the method of moments will likely be negatively biased, whereas an estimate from the use of residual pseudo-likelihood will often be approximately unbiased. Finally, due to the potentially high variability in ICC estimates, approaches to utilize information from reported ICC estimates in order to design future cluster randomized trials have been developed. For instance, see Turner et al.⁶⁰

Although the use of generalized linear mixed models⁴² is also popular, in this article, we focus on the use of GEE to fit a marginal model, as population-average interpretations from statistical models are more natural in cluster randomized trial settings in which clusters are only observed within one trial arm during the study.^4,32,61 Furthermore, with GEE a correct specification of the distribution for the outcomes is not required. We point out that the ICCs for binary data can be on the proportions scale or the log odds scale if using a logit link function, and ICCs for different scales are not equivalent.^13,17 Specifically, the ICC calculated from a GEE-type marginal model is on the proportions scale and is not equivalent to the ICC calculated from a mixed-effects logistic model. More detail is provided in the online supplementary appendix.

Supplemental Material

803635_supp_mat – Supplemental material for A readily available improvement over method of moments for intra-cluster correlation estimation in the context of cluster randomized trials and fitting a GEE–type marginal model for binary outcomes

Supplemental material, 803635_supp_mat for A readily available improvement over method of moments for intra-cluster correlation estimation in the context of cluster randomized trials and fitting a GEE–type marginal model for binary outcomes by Philip M Westgate in Clinical Trials

Footnotes

Acknowledgements

We thank an anonymous associate editor and two reviewers for their constructive comments that helped improve this paper.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Murray

DM.

Design and analysis of group-randomized trials. New York: Oxford University Press, 1998.

Donner

Klar

Design and analysis of cluster randomization trials in health research. London: Edward Arnold, 2000.

Hayes

Moulton

LH.

Cluster randomised trials. Boca Raton, FL: Chapman and Hall/CRC Press, 2009.

Eldridge

Kerry

A practical guide to cluster randomised trials in health services research. Hoboken, NJ: John Wiley & Sons, 2012.

Campbell

Walters

SJ.

How to design, analyse and report cluster randomised trials in medicine and health related research. Hoboken, NJ: John Wiley & Sons, 2014.

Chakraborty

Lyons

Cluster randomized trials: considerations for design and analysis. J Stat Theory Pract2015; 9: 685–698.

Van Breukelen

Candel

. Efficient design of cluster randomized and multicentre trials with unknown intraclass correlation. Stat Methods Med Res2015; 24: 540–556.

Adams

Gulliford

Ukoumunne

et al . Patterns of intra-cluster correlation from primary care research to inform study design and analysis. J Clin Epidemiol2004; 57: 785–794.

Eldridge

Ashby

Feder

et al . Lessons for cluster randomized trials in the twenty-first century: a systematic review of trials in primary care. Clin Trials2004; 1: 80–90.

10.

Roetzheim

Christman

Jacobsen

et al . A randomized controlled trial to increase cancer screening among attendees of community health centers. Ann Fam Med2004; 2: 294–300.

11.

Roetzheim

Christman

Jacobsen

et al . Long-term results from a randomized controlled trial to increase cancer screening among attendees of community health centers. Ann Fam Med2005; 3: 109–114.

12.

Lee

Schell

Roetzheim

Analysis of group randomized trials with multiple binary endpoints and small number of groups. PLoS ONE2009; 4: e7265.

13.

Eldridge

Ukoumunne

Carlin

The intra-cluster correlation coefficient in cluster randomized trials: a review of definitions. Int Stat Rev2009; 77: 378–394.

14.

Campbell

Elbourne

Altman

et al . CONSORT statement: extension to cluster randomised trials. Brit Med J2004; 328: 702–708.

15.

Hade

Murray

Pennell

et al . Intraclass corrrelation estimates for cancer screening outcomes: estimates and applications in the design of group-randomized cancer screening studies. J Natl Cancer Inst Monograph2010; 40: 97–103.

16.

Taljaard

Donner

Villar

et al . Intracluster correlation coefficients from the 2005 WHO global survey on maternal and perinatal health: implications for implementation research. Paediatr Perinat Ep2008; 22: 117–125.

17.

Yelland

Salter

Ryan

et al . Adjusted intraclass correlation coefficients for binary data: methods and estimates from a cluster-randomized trial in primary care. Clin Trials2011; 8: 48–58.

18.

Cook

Bruckner

MacLennan

et al . Clustering in surgical trials—database of intracluster correlations. Trials2012; 13: 2.

19.

Ridout

Demetrio

CGB

Firth

Estimating intraclass correlation for binary data. Biometrics1999; 55: 137–148.

20.

Crespi

Wong

WK.

Comparison of methods for estimating the intraclass correlation coefficient for binary responses in cancer prevention cluster randomized trials. Contemp Clin Trials2012; 33: 869–880.

21.

Donner

Donald

The statistical analysis of multiple binary measurements. J Clin Epidemiol1988; 41: 899–905.

22.

Reed

JF.

Eliminating bias in randomized cluster trials with correlated binomial outcomes. Comput Meth Prog Bio2000; 61: 119–123.

23.

Jung

Kang

Ahn

Sample size calculations for clustered binary data. Stat Med2001; 20: 1971–1982.

24.

Liang

Zeger

SL.

Longitudinal data analysis using generalized linear models. Biometrika1986; 73: 13–22.

25.

R Development Core Team. R: a language and environment for statistical computing, Vienna, Austria, 2011, http://www.R-project.org/

26.

Højsgaard

Halekoh

Yan

The R package geepack for generalized estimating equations. J Stat Softw2006; 15: 1–11.

27.

Yan

Fine

Estimating equations for association structures. Stat Med2004; 23: 859–874.

28.

Yan

Fine

geepack: yet another package for generalized estimating equations. R-News2002; 2: 12–14.

29.

SAS and Institute Inc,. SAS/STAT 12.3 user’s guide. Cary, NC: SAS Institute Inc., 2013.

30.

Chakraborty

Hossain

R package to estimate intracluster correlation coefficient with confidence interval for binary data. Comput Meth Prog Bio2018; 155: 85–92.

31.

Braschel

Svec

Darlington

et al . A comparison of confidence interval methods for the intraclass correlation coefficient in community-based cluster randomization trials with a binary outcome. Clin Trials2016; 13: 180–187.

32.

Preisser

Qaqish

BF.

Finite sample adjustments in estimating equations and covariance estimators for intracluster correlations. Stat Med2008; 27: 5764–5785.

33.

Demetrashvili

Wit

Van Den Heuvel

ER.

Confidence intervals for intraclass correlation coefficients in variance components models. Stat Methods Med Res2016; 25: 2359–2376.

34.

Crespi

Wong

Mishra

SI.

Using second-order generalized estimating equations to model heterogeneous intraclass correlation in cluster-randomized trials. Stat Med2009; 28: 814–827.

35.

Wedderburn

RWM

. Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method. Biometrika1974; 61: 439–447.

36.

Preisser

Qaqish

et al . A comparison of two bias-corrected covariance estimators for generalized estimating equations. Biometrics2007; 63: 935–941.

37.

Westgate

Braun

TM.

Improving small-sample inference in group randomized trials with binary outcomes. Stat Med2011; 30: 201–210.

38.

Westgate

PM.

On small-sample inference in group randomized trials with binary outcomes and cluster-level covariates. Biometrical J2013; 55: 789–806.

39.

Williams

DA.

Extra-binomial variation in logistic-linear models. J Roy Stat Soc C: App1982; 31: 144–148.

40.

Kauermann

Carroll

RJ.

A note on the efficiency of sandwich covariance matrix estimation. J Am Stat Assoc2001; 96: 1387–1396.

41.

Mancl

DeRouen

TA.

A covariance estimator for GEE with improved small-sample properties. Biometrics2001; 57: 126–134.

42.

Fitzmaurice

Laird

Ware

JH.

Applied longitudinal analysis. 2nd ed.Hoboken, NJ: John Wiley & Sons, 2011.

43.

Westgate

PM.

Intra-cluster correlation selection for cluster randomized trials. Stat Med2016; 35: 3272–3284.

44.

Redden

DT.

Small sample performance of bias-corrected sandwich estimators for cluster-randomized trials with binary outcomes. Stat Med2015; 34: 281–296.

45.

Ford

Westgate

PM.

Improved standard error estimator for maintaining the validity of inference in cluster randomized trials with a small number of clusters. Biometrical J2017; 59: 478–495.

46.

Ukoumunne

Thompson

SG.

Analysis of cluster randomized trials with repeated cross-sectional binary measurements. Stat Med2001; 20: 417–433.

47.

Huang

Fiero

Bell

ML.

Generalized estimating equations in cluster randomized trials with a small number of clusters: review of practice and simulation study. Clin Trials2016; 13: 445–449.

48.

Murray

Varnell

Blitstein

JL.

Design and analysis of group-randomized trials: a review of recent methodological developments. Am J Public Health2004; 94: 423–432.

49.

Turner

Heagerty

et al . An evaluation of constrained randomization for the design and analysis of group-randomized trials with binary outcomes. Stat Med2017; 36: 3791–3806.

50.

Crespi

Maxwell

Cluster randomized trials of cancer screening interventions: are appropriate statistical methods being used?

Contemp Clin Trials 2011; 32: 477–484.

51.

Murray

Pals

Blitstein

et al . Design and analysis of group-randomized trials in cancer: a review of current practices. J Natl Cancer I2008; 100: 483–491.

52.

Ivers

Taljaard

Dixon

et al . Impact of CONSORT extension for cluster randomised trials on quality of reporting and study methodology: review of random sample of 300 trials, 2000–8. Brit Med J2011; 343: d5886.

53.

Fay

Graubard

BI.

Small-sample adjustments for Wald-type tests using sandwich estimators. Biometrics2001; 57: 1198–1206.

54.

Morel

Bokossa

Neerchal

NK.

Small sample correction for the variance of GEE estimators. Biometrical J2003; 45: 395–409.

55.

Omar

Thompson

SG.

Analysis of a cluster randomized trial with binary outcome data using a multi-level model. Stat Med2000; 19: 2675–2688.

56.

Turner

Omar

Thompson

SG.

Bayesian methods of analysis for cluster randomized trials with binary outcome data. Stat Med2001; 20: 453–472.

57.

Gulliford

Adams

Ukoumunne

et al . Intraclass correlation coefficient and outcome prevalence are associated in clustered binary data. J Clin Epidemiol2005; 58: 246–251.

58.

Crespi

Wong

A new dependence parameter approach to improve the design of cluster randomized trials with binary outcomes. Clin Trials2011; 8: 687–698.

59.

Feng

Diehr

Yasui

et al . Explaining community-level variance in group randomized trials. Stat Med1999; 18: 539–556.

60.

Turner

Thompson

Spiegelhalter

DJ.

Prior distributions for the intracluster correlation coefficient, based on multiple previous estimates, and their application in cluster randomized trials. Clin Trials 2005; 2: 108–118.

61.

Leyrat

Caille

Donner

et al . Propensity score methods for estimating relative risks in cluster randomized trials with low-incidence binary outcomes and selection bias. Stat Med 2014; 33: 3556–3575.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.24 MB