Cluster Bootstrap Inference for GEE

Why Resample Clusters?

In GEE, observations within a cluster may be correlated, while clusters are treated as independent sampling units. A bootstrap that samples individual observations separately breaks the observed within-cluster patterns: it may combine outcomes and covariates that never occurred together, and it treats correlated observations as if they were independent. The pairs cluster bootstrap instead resamples entire clusters, keeping every sampled cluster’s rows, covariates, outcomes, and internal ordering together.

The procedure approximates repeated sampling of independent clusters from a population of clusters. It is most directly justified when clusters are independent and identically distributed sampling units and the number of clusters is sufficiently large. It is not a universal repair for few-cluster inference: if only a small number of independent clusters were observed, resampling those same clusters cannot create new independent information.

Definition: Suppose the data consist of $K$ independent clusters. Cluster $k$ contains $n_k$ observations, with its full data record denoted $D_k$. A pairs cluster bootstrap replicate draws $K$ indices independently, with replacement, from $\{1,\ldots,K\}$, and fits the GEE to the $K$ selected records $D_k$, including all rows in each selected record. A cluster selected more than once appears more than once in the replicate.

The sampling unit is the cluster, not necessarily a person. For example, if patients are grouped within clinics and the independent sampling units are clinics, resample clinics and retain all sampled patients’ records. If cluster sizes differ, keep each selected cluster’s original size; do not truncate clusters or sample the same number of rows from each one. Thus a bootstrap data set can contain a different total number of rows from the original data, while still containing $K$ sampled clusters.

Implementing the Pairs Cluster Bootstrap

Let $\hat{\beta}$ be the fitted GEE coefficient vector and let $\widehat{V}$ be the covariance estimator chosen for inference. For each bootstrap replicate $b$, draw cluster indices $I_{b1},\ldots,I_{bK}$, independently and uniformly from $\{1,\ldots,K\}$. Refit the GEE to $D_{I_{b1}},\ldots,D_{I_{bK}}$, obtaining $\hat{\beta}^{*(b)}$. If studentized inference is planned, also obtain the matching replicate covariance estimate $\widehat{V}^{*(b)}$.

Identify the independent unit.
Specify which records form one cluster and verify that the analysis data preserve the full cluster membership and within-cluster structure.

Draw cluster indices.
Sample $K$ cluster labels with replacement. Keep all observations in each selected cluster; repeated labels create repeated copies of that cluster.

Refit the full GEE procedure.
Re-estimate the coefficients and any working-correlation parameters as prescribed by the original analysis. Compute the intended covariance estimate if the interval requires it.

Record and diagnose.
Save the coefficient estimates and, for studentization, their standard errors. Record failures rather than silently discarding them.

Form intervals and assess stability.
Use the bootstrap coefficient distribution or studentized statistics, then examine convergence, endpoint stability, and the number of distinct clusters supporting the result.

The bootstrap replicate must repeat the coefficient-fitting procedure, not merely resample a previously calculated list of residuals or score contributions. If the analysis estimates a working correlation, that estimation is generally repeated in each replicate too. The covariance estimator used for a studentized interval should be specified in advance and applied consistently; for example, it may be the ordinary robust sandwich estimator or a chosen small-sample adjustment.

Worked Example: One Cluster-Resampled GEE Fit

Consider an intercept-only logistic GEE with independence working correlation and six clusters, each containing two binary outcomes. The cluster outcome pairs are $D_1=(0,0)$, $D_2=(0,1)$, $D_3=(0,1)$, $D_4=(1,1)$, $D_5=(0,0)$, and $D_6=(1,0)$. The original data contain six successes among twelve observations, so the fitted marginal probability is $6/12=0.5$, and the intercept estimate is $\hat{\beta}=\log(0.5/0.5)=0$.

Suppose one bootstrap draw selects cluster labels $4,2,2,1,5,4$. The replicate contains six clusters and twelve rows: the complete records $D_4,D_2,D_2,D_1,D_5,D_4$. The duplicated clusters are retained as duplicated pairs, not broken apart into individual outcomes.

The selected records contain $2+1+1+0+0+2=6$ successes out of 12 observations. Therefore the replicate’s fitted probability is $6/12=0.5$, and $\hat{\beta}^{*}=\log(0.5/0.5)=0$. Its sandwich standard error, if needed, is computed from these six resampled cluster records. It is not obtained by treating the twelve rows as independent clusters.

Percentile and Studentized Intervals

For a scalar coefficient $\beta_j$, let $\hat{\beta}_j^{*(1)},\ldots,\hat{\beta}_j^{*(B)}$ be the successful bootstrap estimates from $B$ replicates. Write $\hat{q}_p$ for the empirical $p$-quantile of these estimates, using a stated quantile convention. The two-sided percentile interval with nominal level $1-\alpha$ is

[\hat{q}_{\alpha/2},\ \hat{q}_{1-\alpha/2}].

This uses the location and shape of the bootstrap coefficient distribution directly. It is not centered automatically at $\hat{\beta}_j$. Quantile conventions can differ slightly in a finite bootstrap sample, so software settings should be reported when exact reproducibility matters.

A studentized interval scales each replicate’s coefficient difference by that replicate’s estimated standard error. If $\widehat{\mathrm{se}}^{*(b)}_j$ is the standard error for replicate $b$, define

T^{*(b)}_j= \frac{\hat{\beta}^{*(b)}_j-\hat{\beta}_j} {\widehat{\mathrm{se}}^{*(b)}_j}.

Let $\hat{t}_{p}$ denote the empirical $p$-quantile of these $B$ statistics, and let $\widehat{\mathrm{se}}_j$ be the standard error from the original fit. The studentized interval is

\left[ \hat{\beta}_j-\hat{t}_{1-\alpha/2}\widehat{\mathrm{se}}_j,\ \hat{\beta}_j-\hat{t}_{\alpha/2}\widehat{\mathrm{se}}_j \right].

The reversed quantile order in the endpoint formula follows from solving the bootstrap pivot inequalities for $\beta_j$. Studentization can account for changing standard errors across replicates, but it requires those standard errors to be finite and meaningful. It can be especially unstable when a replicate has very few event-bearing clusters.

Worked Example: A Percentile Interval

Suppose nine successful cluster bootstrap replicates give coefficient estimates, in sorted order, $-0.30,-0.20,-0.10,0,0.05,0.10,0.20,0.30,0.50$. The original estimate is $\hat{\beta}=0.08$. Construct a 95% percentile interval using the linear-interpolation quantile convention with index $h=1+(B-1)p$.

Here $B=9$. At $p=0.025$, $h=1+8(0.025)=1.2$, so the quantile is 20% of the way from the first value to the second: $ -0.30+0.2[-0.20-(-0.30)]=-0.28$.

At $p=0.975$, $h=1+8(0.975)=8.8$, so the quantile is 80% of the way from the eighth to the ninth value: $0.30+0.8(0.50-0.30)=0.46$. The percentile interval is therefore $[-0.28,0.46]$. This calculation illustrates the definition, not a recommendation to use only nine replicates: such a small number makes tail quantiles highly coarse and unstable.

Worked Example: A Studentized Interval

Suppose the original coefficient estimate is $\hat{\beta}=0.40$, with $\widehat{\mathrm{se}}=0.20$. Three bootstrap replicates have estimates and standard errors $(0.10,0.15)$, $(0.30,0.20)$, and $(0.70,0.25)$. Their studentized statistics are

\frac{0.10-0.40}{0.15}=-2.00,\qquad \frac{0.30-0.40}{0.20}=-0.50,\qquad \frac{0.70-0.40}{0.25}=1.20.

For illustration, suppose a larger set of replicates yields the 2.5th and 97.5th percentile studentized quantiles $\hat{t}_{0.025}=-1.88$ and $\hat{t}_{0.975}=1.38$. The lower endpoint is $0.40-1.38(0.20)=0.124$; the upper endpoint is $0.40-(-1.88)(0.20)=0.776$. Thus the studentized interval is $[0.124,0.776]$. The three displayed replicates show how individual pivots are calculated; they are not enough to estimate 2.5th and 97.5th percentile quantiles.

What the Bootstrap Approximation Requires

The cluster bootstrap treats the observed clusters as an empirical approximation to the population distribution of clusters. Its basic rationale is clearest when the clusters are independent and identically distributed and the number $K$ grows. For a regular estimator with an asymptotically linear representation, the coefficient error is approximately an average of cluster-level influence contributions. Resampling clusters then resamples the units contributing to that average. This rationale does not guarantee good finite-sample coverage, especially when $K$ is small, the estimator is near a boundary, or some resamples cannot be fitted.

Result (conditional variance of resampled cluster contributions): For fixed scalar cluster contributions $a_1,\ldots,a_K$, let $I_1,\ldots,I_K$ be independent uniform draws from $\{1,\ldots,K\}$, and set $\bar{a}^{*}=K^{-1}\sum_{r=1}^K a_{I_r}$ and $\bar{a}=K^{-1}\sum_{k=1}^K a_k$. Conditional on the observed contributions, $$ E[\bar{a}^{*}\mid a_1,\ldots,a_K]=\bar{a},\qquad \operatorname{Var}(\bar{a}^{*}\mid a_1,\ldots,a_K) =\frac{1}{K}\left\{\frac{1}{K}\sum_{k=1}^K(a_k-\bar{a})^2\right\}. $$

Proof. Each draw $a_{I_r}$, conditional on the observed values, has mean $K^{-1}\sum_k a_k=\bar{a}$ and variance $K^{-1}\sum_k(a_k-\bar{a})^2$. The $K$ draws are conditionally independent. Therefore the mean of their average is $\bar{a}$, and the variance of their average is $K^{-2}$ times the sum of their $K$ identical variances. This gives the displayed expression. $\square$

This identity explains why whole-cluster resampling reproduces cluster-level sampling variation in an average. For a GEE coefficient, the relevant contributions are generally vectors and depend on the fitted estimating procedure; the identity alone is not a proof of bootstrap validity for every GEE estimator. Validity also requires suitable regularity, including a stable fit and enough independent clusters for the large-sample approximation.

Result (percentile quantile consistency): Let $F$ be the limiting distribution function for a suitably centered and scaled estimator, and let $q_p$ be its $p$-quantile. Suppose $F$ is continuous and strictly increasing in a neighborhood of $q_p$. If the conditional bootstrap distribution function converges in probability to $F$ at every point in that neighborhood, then its conditional $p$-quantile converges in probability to $q_p$.

Proof. Fix any $\varepsilon>0$ small enough that both $q_p-\varepsilon$ and $q_p+\varepsilon$ lie in the neighborhood. Strict increase and continuity at the quantile give $F(q_p-\varepsilon)<p<F(q_p+\varepsilon)$. By the assumed conditional convergence, with probability tending to one the bootstrap distribution function at $q_p-\varepsilon$ is below $p$, and its value at $q_p+\varepsilon$ is above $p$. On that event, the conditional bootstrap $p$-quantile lies between $q_p-\varepsilon$ and $q_p+\varepsilon$. Since this holds for every such $\varepsilon$, the conditional quantile converges in probability to $q_p$. $\square$

The result clarifies the role of approximation quality: the quantile construction is reliable only when the bootstrap distribution approximates the relevant sampling distribution. In a GEE analysis, that approximation can be poor with few clusters or unstable estimating equations, even if many bootstrap replicates are generated.

Failure Modes: Few Clusters and Sparse Outcomes

A large number of bootstrap replicates reduces Monte Carlo error in the estimated bootstrap distribution; it does not increase the number of independent clusters in the original study. With few clusters, the empirical cluster distribution has limited diversity. Many resamples repeat the same clusters, and intervals can be sensitive to one influential cluster. This is why the bootstrap should be considered alongside the few-cluster limitations and small-sample methods already relevant to GEE inference, not as a substitute for them.

Sparse binary outcomes create a more direct computational problem. A logistic GEE replicate can contain no events, or no non-events, for a covariate pattern. In such a replicate, a finite coefficient estimate may not exist. Even when the fit converges, extremely large coefficients and standard errors can make percentile endpoints or studentized statistics erratic. A report should state how many replicates failed, why they failed when known, and what rule was used for the analysis. Silently dropping failures changes the bootstrap distribution to one conditional on successful fitting.

Worked Example: Chance of a No-Event Replicate

Suppose $K=6$ clusters are observed, exactly one of which contains any events; the other five have all-zero outcomes. In a cluster bootstrap replicate of six draws, the probability of selecting no event-bearing cluster is the probability that all six draws come from the five event-free clusters:

P(\text{no event-bearing cluster}) =\left(\frac{5}{6}\right)^6 =\frac{15625}{46656}\approx 0.335.

Thus about one third of replicates are expected to contain only zero outcomes. An intercept-only logistic fit to such a replicate has no finite maximum-likelihood intercept: the fitted event probability approaches zero as the intercept tends to $-\infty$. Increasing the number of bootstrap replicates will reveal this failure more precisely; it will not remove the underlying lack of event-bearing clusters.

Practical Reporting and a Common Mistake

A useful report states the number of independent clusters, the resampling unit, the number of attempted and successful replicates, the interval type and quantile convention, and the covariance estimator used for studentization. It should also describe convergence problems and assess whether interval endpoints change materially with a larger number of replicates. If failed replicates are common, that is substantive evidence of instability, not merely a software inconvenience.

Common mistake: “There are hundreds of observations, so I can bootstrap individual rows and obtain a precise interval.” This reasoning confuses the number of rows with the number of independent sampling units. If observations are correlated within clusters, individual-row resampling destroys that dependence and exaggerates the independent information. The correction is to resample the independent clusters as intact records and to judge inferential reliability using the number and diversity of those clusters. A very large row count does not compensate for having only a few clusters.

Key takeaway: A pairs cluster bootstrap resamples independent clusters with replacement and refits the GEE to every sampled set of complete cluster records. Percentile and studentized intervals summarize that refitting distribution, but neither method can manufacture independent clusters or guarantee stable inference under sparse outcomes.

Check Your Understanding

Answer each question by identifying both the resampling operation and the inferential issue it addresses.

A study has 30 clinics with unequal numbers of patients. What is sampled in one pairs cluster bootstrap replicate, and what happens to a clinic selected twice?
State the percentile interval endpoints for a scalar coefficient in terms of bootstrap quantiles.
Why does a studentized interval use a replicate-specific standard error, and why are its endpoint quantiles reversed in the interval formula?
A logistic GEE replicate contains no events. Why may its coefficient estimate fail to be finite, and why should the replicate not be silently ignored?
Why does increasing the number of bootstrap replicates not solve the problem of having only a few independent clusters?