Small sample performance of bias‐corrected sandwich estimators for cluster‐randomized trials with binary outcomes

The sandwich estimator in generalized estimating equations (GEE) approach underestimates the true variance in small samples and consequently results in inflated type I error rates in hypothesis testing. This fact limits the application of the GEE in cluster‐randomized trials (CRTs) with few clusters. Under various CRT scenarios with correlated binary outcomes, we evaluate the small sample properties of the GEE Wald tests using bias‐corrected sandwich estimators. Our results suggest that the GEE Wald z‐test should be avoided in the analyses of CRTs with few clusters even when bias‐corrected sandwich estimators are used. With t‐distribution approximation, the Kauermann and Carroll (KC)‐correction can keep the test size to nominal levels even when the number of clusters is as low as 10 and is robust to the moderate variation of the cluster sizes. However, in cases with large variations in cluster sizes, the Fay and Graubard (FG)‐correction should be used instead. Furthermore, we derive a formula to calculate the power and minimum total number of clusters one needs using the t‐test and KC‐correction for the CRTs with binary outcomes. The power levels as predicted by the proposed formula agree well with the empirical powers from the simulations. The proposed methods are illustrated using real CRT data. We conclude that with appropriate control of type I error rates under small sample sizes, we recommend the use of GEE approach in CRTs with binary outcomes because of fewer assumptions and robustness to the misspecification of the covariance structure. Copyright © 2014 John Wiley & Sons, Ltd.     

Covariance estimators for generalized estimating equations (GEE) in longitudinal analysis with small samples

Generalized estimating equations (GEE) is a general statistical method to fit marginal models for longitudinal data in biomedical studies. The variance–covariance matrix of the regression parameter coefficients is usually estimated by a robust “sandwich” variance estimator, which does not perform satisfactorily when the sample size is small. To reduce the downward bias and improve the efficiency, several modified variance estimators have been proposed for bias‐correction or efficiency improvement. In this paper, we provide a comprehensive review on recent developments of modified variance estimators and compare their small‐sample performance theoretically and numerically through simulation and real data examples. In particular, Wald tests and t‐tests based on different variance estimators are used for hypothesis testing, and the guideline on appropriate sample sizes for each estimator is provided for preserving type I error in general cases based on numerical results. Moreover, we develop a user‐friendly R package “geesmv” incorporating all of these variance estimators for public usage in practice. Copyright © 2015 John Wiley & Sons, Ltd.     

Sample size calculation in three-level cluster randomized trials using generalized estimating equation models

Three-level cluster randomized trials (CRTs) are increasingly used in implementation science, where 2fold-nested-correlated data arise. For example, interventions are randomly assigned to practices, and providers within the same practice who provide care to participants are trained with the assigned intervention. Teerenstra et al proposed a nested exchangeable correlation structure that accounts for two levels of clustering within the generalized estimating equations (GEE) approach. In this article, we utilize GEE models to test the treatment effect in a two-group comparison for continuous, binary, or count data in three-level CRTs. Given the nested exchangeable correlation structure, we derive the asymptotic variances of the estimator of the treatment effect for different types of outcomes. When the number of clusters is small, researchers have proposed bias-corrected sandwich estimators to improve performance in two-level CRTs. We extend the variances of two bias-corrected sandwich estimators to three-level CRTs. The equal provider and practice sizes were assumed to calculate number of practices for simplicity. However, they are not guaranteed in practice. Relative efficiency (RE) is defined as the ratio of variance of the estimator of the treatment effect for equal to unequal provider and practice sizes. The expressions of REs are obtained from both asymptotic variance estimation and bias-corrected sandwich estimators. Their performances are evaluated for different scenarios of provider and practice size distributions through simulation studies. Finally, a percentage increase in the number of practices is proposed due to efficiency loss from unequal provider and/or practice sizes.     

A determinant‐based criterion for working correlation structure selection in generalized estimating equations

In generalized estimating equations (GEE), the correlation between the repeated observations on a subject is specified with a working correlation matrix. Correct specification of the working correlation structure ensures efficient estimators of the regression coefficients. Among the criteria used, in practice, for selecting working correlation structure, Rotnitzky‐Jewell, Quasi Information Criterion (QIC) and Correlation Information Criterion (CIC) are based on the fact that if the assumed working correlation structure is correct then the model‐based (naive) and the sandwich (robust) covariance estimators of the regression coefficient estimators should be close to each other. The sandwich covariance estimator, used in defining the Rotnitzky‐Jewell, QIC and CIC criteria, is biased downward and has a larger variability than the corresponding model‐based covariance estimator. Motivated by this fact, a new criterion is proposed in this paper based on the bias‐corrected sandwich covariance estimator for selecting an appropriate working correlation structure in GEE. A comparison of the proposed and the competing criteria is shown using simulation studies with correlated binary responses. The results revealed that the proposed criterion generally performs better than the competing criteria. An example of selecting the appropriate working correlation structure has also been shown using the data from Madras Schizophrenia Study. Copyright © 2015 John Wiley & Sons, Ltd.     

A bias‐corrected covariance estimate for improved inference with quadratic inference functions

The method of quadratic inference functions (QIF) is an increasingly popular method for the analysis of correlated data because of its multiple advantages over generalized estimating equations (GEE). One advantage is that it is more efficient for parameter estimation when the working covariance structure for the data is misspecified. In the QIF literature, the asymptotic covariance formula is used to obtain standard errors. We show that in small to moderately sized samples, these standard error estimates can be severely biased downward, therefore inflating test size and decreasing coverage probability. We propose adjustments to the asymptotic covariance formula that eliminate finite‐sample biases and, as shown via simulation, lead to substantial improvements in standard error estimates, inference, and coverage. The proposed method is illustrated in application to a cluster randomized trial and a longitudinal study. Furthermore, QIF and GEE are contrasted via simulation and these applications. Copyright © 2012 John Wiley & Sons, Ltd.     

Small sample GEE estimation of regression parameters for longitudinal data

Longitudinal (clustered) response data arise in many bio‐statistical applications which, in general, cannot be assumed to be independent. Generalized estimating equation (GEE) is a widely used method to estimate marginal regression parameters for correlated responses. The advantage of the GEE is that the estimates of the regression parameters are asymptotically unbiased even if the correlation structure is misspecified, although their small sample properties are not known. In this paper, two bias adjusted GEE estimators of the regression parameters in longitudinal data are obtained when the number of subjects is small. One is based on a bias correction, and the other is based on a bias reduction. Simulations show that the performances of both the bias‐corrected methods are similar in terms of bias, efficiency, coverage probability, average coverage length, impact of misspecification of correlation structure, and impact of cluster size on bias correction. Both these methods show superior properties over the GEE estimates for small samples. Further, analysis of data involving a small number of subjects also shows improvement in bias, MSE, standard error, and length of the confidence interval of the estimates by the two bias adjusted methods over the GEE estimates. For small to moderate sample sizes (), either of the bias‐corrected methods GEEBc and GEEBr can be used. However, the method GEEBc should be preferred over GEEBr, as the former is computationally easier. For large sample sizes, the GEE method can be used. Copyright © 2014 John Wiley & Sons, Ltd.     

Meta‐analysis of gene‐environment interaction exploiting gene‐environment independence across multiple case‐control studies

Multiple papers have studied the use of gene‐environment (G‐E) independence to enhance power for testing gene‐environment interaction in case‐control studies. However, studies that evaluate the role of G‐E independence in a meta‐analysis framework are limited. In this paper, we extend the single‐study empirical Bayes type shrinkage estimators proposed by Mukherjee and Chatterjee (2008) to a meta‐analysis setting that adjusts for uncertainty regarding the assumption of G‐E independence across studies. We use the retrospective likelihood framework to derive an adaptive combination of estimators obtained under the constrained model (assuming G‐E independence) and unconstrained model (without assumptions of G‐E independence) with weights determined by measures of G‐E association derived from multiple studies. Our simulation studies indicate that this newly proposed estimator has improved average performance across different simulation scenarios than the standard alternative of using inverse variance (covariance) weighted estimators that combines study‐specific constrained, unconstrained, or empirical Bayes estimators. The results are illustrated by meta‐analyzing 6 different studies of type 2 diabetes investigating interactions between genetic markers on the obesity related FTO gene and environmental factors body mass index and age.     

Modified robust variance estimator for generalized estimating equations with improved small-sample performance

Generalized estimating equations (GEE (Biometrika 1986; 73(1):13-22) is a general statistical method to fit marginal models for correlated or clustered responses, and it uses a robust sandwich estimator to estimate the variance-covariance matrix of the regression coefficient estimates. While this sandwich estimator is robust to the misspecification of the correlation structure of the responses, its finite sample performance deteriorates as the number of clusters or observations per cluster decreases. To address this limitation, Pan (Biometrika 2001; 88(3):901-906) and Mancl and DeRouen (Biometrics 2001; 57(1):126-134) investigated two modifications to the original sandwich variance estimator. Motivated by the ideas underlying these two modifications, we propose a novel robust variance estimator that combines the strengths of these estimators. Our theoretical and numerical results show that the proposed estimator attains better efficiency and achieves better finite sample performance compared with existing estimators. In particular, when the sample size or cluster size is small, our proposed estimator exhibits lower bias and the resulting confidence intervals for GEE estimates achieve better coverage rates performance. We illustrate the proposed method using data from a dental study.     

Comparison of the risk difference, risk ratio and odds ratio scales for quantifying the unadjusted intervention effect in cluster randomized trials

This paper evaluates methods for unadjusted analyses of binary outcomes in cluster randomized trials (CRTs). Under the generalized estimating equations (GEE) method the identity, log and logit link functions may be specified to make inferences on the risk difference, risk ratio and odds ratio scales, respectively. An alternative, 'cluster-level', method applies the t-test to summary statistics calculated for each cluster, using proportions, log proportions and log odds, to make inferences on the respective scales. Simulation was used to estimate the bias of the unadjusted intervention effect estimates and confidence interval coverage, generating data sets with different combinations of number of clusters, number of participants per cluster, intra-cluster correlation coefficient rho and intervention effect. When the identity link was specified, GEE had little bias and good coverage, performing slightly better than the log and logit link functions. The cluster-level method provided unbiased point estimates when proportions were used to summarize the clusters. When the log proportion and log odds were used, however, the method often had markedly large bias for two reasons: (i) bias in the modified summary statistic used for cluster-level estimation when a cluster has zero cases with the outcome of interest (arising when the number of participants sampled per cluster is small and the outcome prevalence is low) and (ii) asymptotically, the method estimates the ratio of geometric means of the cluster proportions or odds, respectively, between the trial arms rather than the ratio of arithmetic means.     

Analysis of counts for cluster randomized trials: Negative controls and test-negative designs

In cluster randomized trials (CRTs), the outcome of interest is often a count at the cluster level. This occurs, for example, in evaluating an intervention with the outcome being the number of infections of a disease such as HIV or dengue or the number of hospitalizations in the cluster. Standard practice analyzes these counts through cluster outcome rates using an appropriate denominator (eg, population size). However, such denominators are sometimes unknown, particularly when the counts depend on a passive community surveillance system. We consider direct comparison of the counts without knowledge of denominators, relying on randomization to balance denominators. We also focus on permutation tests to allow for small numbers of randomized clusters. However, such approaches are subject to bias when there is differential ascertainment of counts across arms, a situation that may occur in CRTs that cannot implement blinded interventions. We suggest the use of negative control counts as a method to remove, or reduce, this bias, discussing the key properties necessary for an effective negative control. A current example of such a design is the recent extension of test-negative designs to CRTs testing community-level interventions. Via simulation, we compare the performance of new and standard estimators based on CRTs with negative controls to approaches that only use the original counts. When there is no differential ascertainment by intervention arm, the count-only approaches perform comparably to those using debiasing negative controls. However, under even modest differential ascertainment, the count-only estimators are no longer reliable.     

On the estimation of intracluster correlation for time‐to‐event outcomes in cluster randomized trials

Cluster randomized trials (CRTs) involve the random assignment of intact social units rather than independent subjects to intervention groups. Time‐to‐event outcomes often are endpoints in CRTs. Analyses of such data need to account for the correlation among cluster members. The intracluster correlation coefficient (ICC) is used to assess the similarity among binary and continuous outcomes that belong to the same cluster. However, estimating the ICC in CRTs with time‐to‐event outcomes is a challenge because of the presence of censored observations. The literature suggests that the ICC may be estimated using either censoring indicators or observed event times. A simulation study explores the effect of administrative censoring on estimating the ICC. Results show that ICC estimators derived from censoring indicators or observed event times are negatively biased. Analytic work further supports these results. Observed event times are preferred to estimate the ICC under minimum frequency of administrative censoring. To our knowledge, the existing literature provides no practical guidance on the estimation of ICC when substantial amount of administrative censoring is present. The results from this study corroborate the need for further methodological research on estimating the ICC for correlated time‐to‐event outcomes. Copyright © 2016 John Wiley & Sons, Ltd.     

Stepped wedge designs could reduce the required sample size in cluster randomized trials

ObjectiveThe stepped wedge design is increasingly being used in cluster randomized trials (CRTs). However, there is not much information available about the design and analysis strategies for these kinds of trials. Approaches to sample size and power calculations have been provided, but a simple sample size formula is lacking. Therefore, our aim is to provide a sample size formula for cluster randomized stepped wedge designs.Study Design and SettingWe derived a design effect (sample size correction factor) that can be used to estimate the required sample size for stepped wedge designs. Furthermore, we compared the required sample size for the stepped wedge design with a parallel group and analysis of covariance (ANCOVA) design.ResultsOur formula corrects for clustering as well as for the design. Apart from the cluster size and intracluster correlation, the design effect depends on choices of the number of steps, the number of baseline measurements, and the number of measurements between steps. The stepped wedge design requires a substantial smaller sample size than a parallel group and ANCOVA design.ConclusionFor CRTs, the stepped wedge design is far more efficient than the parallel group and ANCOVA design in terms of sample size.     

Measures of between‐cluster variability in cluster randomized trials with binary outcomes

Cluster randomized trials (CRTs) are increasingly used to evaluate the effectiveness of health‐care interventions. A key feature of CRTs is that the observations on individuals within clusters are correlated as a result of between‐cluster variability. Sample size formulae exist which account for such correlations, but they make different assumptions regarding the between‐cluster variability in the intervention arm of a trial, resulting in different sample size estimates. We explore the relationship for binary outcome data between two common measures of between‐cluster variability: k, the coefficient of variation and ρ, the intracluster correlation coefficient. We then assess how the assumptions of constant k or ρ across treatment arms correspond to different assumptions about intervention effects. We assess implications for sample size estimation and present a simple solution to the problems outlined. Copyright © 2009 John Wiley & Sons, Ltd.     

Optimal designs in three-level cluster randomized trials with a binary outcome

Cluster randomized trials (CRTs) were originally proposed for use when randomization at the subject level is practically infeasible or may lead to a severe estimation bias of the treatment effect. However, recruiting an additional cluster costs more than enrolling an additional subject in an individually randomized trial. Under budget constraints, researchers have proposed the optimal sample sizes in two-level CRTs. CRTs may have a three-level structure, in which two levels of clustering should be considered. In this paper, we propose optimal designs in three-level CRTs with a binary outcome, assuming a nested exchangeable correlation structure in generalized estimating equation models. We provide the variance of estimators of three commonly used measures: risk difference, risk ratio, and odds ratio. For a given sampling budget, we discuss how many clusters and how many subjects per cluster are necessary to minimize the variance of each measure estimator. For known association parameters, the locally optimal design is proposed. When association parameters are unknown but within predetermined ranges, the MaxiMin design is proposed to maximize the minimum of relative efficiency over the possible ranges, that is, to minimize the risk of the worst scenario.     

Exact and approximate unconditional confidence intervals for proportion difference in the presence of incomplete data

Confidence interval (CI) construction with respect to proportion/rate difference for paired binary data has become a standard procedure in many clinical trials and medical studies. When the sample size is small and incomplete data are present, asymptotic CIs may be dubious and exact CIs are not yet available. In this article, we propose exact and approximate unconditional test‐based methods for constructing CI for proportion/rate difference in the presence of incomplete paired binary data. Approaches based on one‐ and two‐sided Wald's tests will be considered. Unlike asymptotic CI estimators, exact unconditional CI estimators always guarantee their coverage probabilities at or above the pre‐specified confidence level. Our empirical studies further show that (i) approximate unconditional CI estimators usually yield shorter expected confidence width (ECW) with their coverage probabilities being well controlled around the pre‐specified confidence level; and (ii) the ECWs of the unconditional two‐sided‐test‐based CI estimators are generally narrower than those of the unconditional one‐sided‐test‐based CI estimators. Moreover, ECWs of asymptotic CIs may not necessarily be narrower than those of two‐sided‐based exact unconditional CIs. Two real examples will be used to illustrate our methodologies. Copyright © 2008 John Wiley & Sons, Ltd.     

A comparison of strategies for analyzing dichotomous outcomes in genome-wide association studies with general pedigrees

Genome-wide association studies (GWAS) have been frequently conducted on general or isolated populations with related individuals. However, there is a lack of consensus on which strategy is most appropriate for analyzing dichotomous phenotypes in general pedigrees. Using simulation studies, we compared several strategies including generalized estimating equations (GEE) strategies with various working correlation structures, generalized linear mixed model (GLMM), and a variance component strategy (denoted LMEBIN) that treats dichotomous outcomes as continuous with special attentions to their performance with rare variants, rare diseases, and small sample sizes. In our simulations, when the sample size is not small, for type I error, only GEE and LMEBIN maintain nominal type I error in most cases with exceptions for GEE with very rare disease and genetic variants. GEE and LMEBIN have similar statistical power and slightly outperform GLMM when the prevalence is low. In terms of computational efficiency, GEE with sandwich variance estimator outperforms GLMM and LMEBIN. We apply the strategies to GWAS of gout in the Framingham Heart Study. Based on our results, we would recommend using GEE ind-san in the GWAS for common variants and GEE ind-fij or LMEBIN for rare variants for GWAS of dichotomous outcomes with general pedigrees.     

Stepped‐wedge cluster randomised controlled trials: a generic framework including parallel and multiple‐level designs

Stepped‐wedge cluster randomised trials (SW‐CRTs) are being used with increasing frequency in health service evaluation. Conventionally, these studies are cross‐sectional in design with equally spaced steps, with an equal number of clusters randomised at each step and data collected at each and every step. Here we introduce several variations on this design and consider implications for power. One modification we consider is the incomplete cross‐sectional SW‐CRT, where the number of clusters varies at each step or where at some steps, for example, implementation or transition periods, data are not collected. We show that the parallel CRT with staggered but balanced randomisation can be considered a special case of the incomplete SW‐CRT. As too can the parallel CRT with baseline measures. And we extend these designs to allow for multiple layers of clustering, for example, wards within a hospital. Building on results for complete designs, power and detectable difference are derived using a Wald test and obtaining the variance–covariance matrix of the treatment effect assuming a generalised linear mixed model. These variations are illustrated by several real examples. We recommend that whilst the impact of transition periods on power is likely to be small, where they are a feature of the design they should be incorporated. We also show examples in which the power of a SW‐CRT increases as the intra‐cluster correlation (ICC) increases and demonstrate that the impact of the ICC is likely to be smaller in a SW‐CRT compared with a parallel CRT, especially where there are multiple levels of clustering. Finally, through this unified framework, the efficiency of the SW‐CRT and the parallel CRT can be compared. © 2014 The Authors. Statistics in Medicine Published by John Wiley & Sons Ltd.     

Properties of analysis methods that account for clustering in volume-outcome studies when the primary predictor is cluster size

In recent years health services researchers have conducted 'volume-outcome' studies to evaluate whether providers (hospitals or surgeons) who treat many patients for a specialized condition have better outcomes than those that treat few patients. These studies and the inherent clustering of events by provider present an unusual statistical problem. The volume-outcome setting is unique in that 'volume' reflects both the primary factor under study and also the cluster size. Consequently, the assumptions inherent in the use of available methods that correct for clustering might be violated in this setting. To address this issue, we investigate via simulation the properties of three estimation procedures for the analysis of cluster correlated data, specifically in the context of volume-outcome studies. We examine and compare the validity and efficiency of widely-available statistical techniques that have been used in the context of volume-outcome studies: generalized estimating equations (GEE) using both the independence and exchangeable correlation structures; random effects models; and the weighted GEE approach proposed by Williamson et al. (Biometrics 2003; 59:36-42) to account for informative clustering. Using data generated either from an underlying true random effects model or a cluster correlated model we show that both the random effects and the GEE with an exchangeable correlation structure have generally good properties, with relatively low bias for estimating the volume parameter and its variance. By contrast, the cluster weighted GEE method is inefficient.     

Interval estimation of the over‐dispersion parameter in the analysis of one‐way layout of count data

The over‐dispersion parameter is an important and versatile measure in the analysis of one‐way layout of count data in biological studies. For example, it is commonly used as an inverse measure of aggregation in biological count data. Its estimation from finite data sets is a recognized challenge. Many simulation studies have examined the bias and efficiency of different estimators of the over‐dispersion parameter for finite data sets (see, for example, Clark and Perry, Biometrics 1989; 45:309–316 and Piegorsch, Biometrics 1990; 46:863–867), but little attention has been paid to the accuracy of the confidence intervals (CIs) of it. In this paper, we first derive asymptotic procedures for the construction of confidence limits for the over‐dispersion parameter using four estimators that are specified by only the first two moments of the counts. We also obtain closed‐form asymptotic variance formulae for these four estimators. In addition, we consider the asymptotic CI based on the maximum likelihood (ML) estimator using the negative binomial model. It appears from the simulation results that the asymptotic CIs based on these five estimators have coverage below the nominal coverage probability. To remedy this, we also study the properties of the asymptotic CIs based on the restricted estimates of ML, extended quasi‐likelihood, and double extended quasi‐likelihood by eliminating the nuisance parameter effect using their adjusted profile likelihood and quasi‐likelihoods. It is shown that these CIs outperform the competitors by providing coverage levels close to nominal over a wide range of parameter combinations. Two examples to biological count data are presented. Copyright © 2010 John Wiley & Sons, Ltd.     

Sample size considerations for stratified cluster randomization design with binary outcomes and varying cluster size

Stratified cluster randomization trials (CRTs) have been frequently employed in clinical and healthcare research. Comparing with simple randomized CRTs, stratified CRTs reduce the imbalance of baseline prognostic factors among different intervention groups. Due to the popularity, there has been a growing interest in methodological development on sample size estimation and power analysis for stratified CRTs; however, existing work mostly assumes equal cluster size within each stratum and uses multilevel models. Clusters are often naturally formed with random sizes in CRTs. With varying cluster size, commonly used ad hoc approaches ignore the variability in cluster size, which may underestimate (overestimate) the required number of clusters for each group per stratum and lead to underpowered (overpowered) clinical trials. We propose closed-form sample size formulas for estimating the required total number of subjects and for estimating the number of clusters for each group per stratum, based on Cochran-Mantel-Haenszel statistic for stratified cluster randomization design with binary outcomes, accounting for both clustering and varying cluster size. We investigate the impact of various design parameters on the relative change in the required number of clusters for each group per stratum due to varying cluster size. Simulation studies are conducted to evaluate the finite-sample performance of the proposed sample size method. A real application example of a pragmatic stratified CRT of a triad of chronic kidney disease, diabetes, and hypertension is presented for illustration.