When we need population-level insights without the distributional baggage of random effects

Picture this: We’ve just finished implementing that elegant mixed effects model from last week’s deep dive, and the hospital executives are thrilled with our patient length-of-stay predictions. But then the epidemiologist on our team asks a different question: “What’s the average effect of this new treatment protocol across all our hospitals?”

Our mixed effects model gives us subject-specific predictions – how much longer will this particular patient stay at this specific hospital? The policy question requires population-average inference – if we implement this protocol system-wide, what’s the expected change in average length of stay across the entire health system?

Now, our mixed effects model could answer this population-average question. We could average predictions across the distribution of hospital random effects, use marginal standardization techniques, or employ other post-hoc methods. But doing so requires additional computational steps, strong distributional assumptions about those random effects, and careful handling of the integration over the random effects distribution.

This moment crystallizes a fundamental choice in statistical modeling: Do we want population-average effects as a byproduct of a more complex individual-level model, or do we want them directly as our primary target? Welcome to Generalized Estimating Equations (GEE) – a different philosophical approach that puts population-average inference front and center.

Disclaimer: Python codes provided in this article are all to understand algorithmic steps, rather than fully functioning implementations

The Question That Changes Everything

Last week, we explored mixed effects models that decompose patient outcomes into universal patterns plus hospital-specific deviations. The mathematical beauty lay in modeling the distribution of those random effects, creating a hierarchical structure that enables both population and individual predictions.

GEE takes a radically different approach with a deceptively simple question: What if we don’t care about modeling individual group effects at all?

Instead of asking “How does Hospital A differ from the population?” GEE asks “What’s the average relationship between treatment and outcome across all hospitals, accounting for the fact that patients within hospitals are correlated?”

This philosophical shift – from conditional to marginal inference – has profound implications that ripple through methodology, computation, and application.

The Marginal vs. Conditional Divide

Consider this concrete example from our medication adherence analysis:

Mixed Effects Question: For a patient with diabetes managed at Penn Presbyterian, what’s their predicted medication adherence given their demographic profile and hospital-specific care patterns? (Primary output: conditional predictions; population averages require additional computation)

GEE Question: Across the entire Penn Medicine health system, what’s the population-average improvement in medication adherence for patients who receive automated text reminders? (Primary output: marginal effects directly)

The mixed effects model provides conditional inference as its natural output – predictions conditional on being at a specific hospital with specific random effects. Population-average effects are obtainable but require integrating over the random effects distribution. GEE provides marginal inference directly – population-averaged effects without requiring distributional assumptions about group-level variation.

This distinction affects both computational efficiency and robustness. Mixed effects models give us more information (both individual and population effects) but at the cost of stronger assumptions and more complex computation. GEE gives us exactly what we need for population-level questions, more robustly and directly.

Why We Need Another Approach: The Directness and Robustness Revolution

The skeptic in us might ask: Why not just use mixed effects models and compute population averages from them? After all, mixed effects models can provide population-level inferences through marginal standardization or averaging over the random effects distribution.

The answer reveals why GEE has become indispensable: computational directness and distributional robustness.

Mixed effects models are elegant and informative, but they require additional steps to obtain population averages. More critically, they demand strong assumptions about the world:

  • Random effects follow a specific distribution (usually normal)
  • Correct specification of the random effects structure
  • Proper modeling of the variance components
  • The relationship between individual and population effects

When these assumptions hold, mixed effects models provide rich information about both individual predictions and population averages. But real healthcare data has an uncomfortable habit of violating theoretical assumptions. When our random effects aren’t actually normally distributed (spoiler: they often aren’t), when our variance structures are misspecified, or when we’re unsure about the hierarchical relationships, the population-average inferences derived from mixed effects models can be biased.

GEE sidesteps this complexity entirely by targeting population averages directly. It focuses only on two things:

  • The mean structure (how covariates relate to outcomes)
  • The correlation structure (how observations within groups are related)

If our correlation assumptions are wrong in GEE, our population-average inferences remain consistent. If our random effects assumptions are wrong in a mixed model, both our individual predictions and our derived population averages are in trouble.

This robustness isn’t theoretical luxury – it’s practical necessity when we need reliable population-level estimates without the distributional baggage.

The Mathematical Foundation: Quasi-Likelihood and Working Correlations

Let’s build the mathematical machinery that makes this robustness possible, piece by piece.

For clustered data, we have observations $y$ for subject $i$ in cluster $j$. GEE models the marginal expectation:

$$E[y_{ij}] = \mu_{ij} = g^{-1}(X_{ij}\beta)$$

where $g$ is a link function and $\beta$ represents population-average effects.

The breakthrough insight is the working correlation matrix that captures within-cluster correlation without requiring distributional assumptions:

$$\text{Var}(Y_j) = V_j = A_j^{1/2} R_j(\alpha) A_j^{1/2}$$

where $A$ denote diagonal matrices of variance functions.

The “working” designation is crucial – we don’t need this correlation structure to be correct, just reasonable. This is where GEE’s elegance becomes apparent.

Working Correlation Structures: The Art of Being Approximately Right

Common working correlation structures include:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# Working independence: assume uncorrelated observations
def independence(cluster_size):
    return identity_matrix(cluster_size)

# Exchangeable: constant correlation within cluster  
def exchangeable(cluster_size, alpha):
    R = alpha * ones_matrix(cluster_size, cluster_size)
    diagonal(R) = 1.0
    return R

# AR(1): correlation decays with time/distance
def ar1(cluster_size, alpha):
    R = zeros_matrix(cluster_size, cluster_size)
    for i in range(cluster_size):
        for j in range(cluster_size):
            R[i, j] = alpha ** abs(i - j)
    return R

The beauty is that even if we choose the wrong structure, our population-average estimates remain consistent thanks to the sandwich estimator – the mathematical safety net that makes GEE so robust.

Building GEE from Scratch: The Quasi-Likelihood Revolution

Imagine we implement GEE step by step, revealing the algorithmic elegance that makes it both robust and computationally reliable.

The Foundation

We first need a different kind of likelihood.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

class GEEFromScratch:
    def __init__(self, family='gaussian', link='identity', correlation='exchangeable'):
        self.family = family
        self.link = link  
        self.correlation = correlation
        
    def link_function(self, eta):
        if self.link == 'identity': return eta
        elif self.link == 'logit': return 1 / (1 + exp(-eta))
        elif self.link == 'log': return exp(eta)
    
    def variance_function(self, mu):
        if self.family == 'gaussian': return ones_like(mu)
        elif self.family == 'binomial': return mu * (1 - mu)
        elif self.family == 'poisson': return mu

The Iterative Heart: Where Robustness Emerges

The core of GEE is an iterative algorithm that alternates between estimating regression parameters and correlation parameters. Unlike mixed effects models that wrestle with complex likelihood surfaces, GEE’s iterations are remarkably well-behaved:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

def fit(self, y, X, clusters, max_iter=50, tol=1e-6):
    # Initialize with OLS estimates
    beta = least_squares_solution(X, y)
    
    # Group observations by cluster
    cluster_indices = group_observations_by_cluster(clusters)
    
    for iteration in range(max_iter):
        beta_old = beta.copy()
        
        # E-step: Update working correlation parameters
        alpha = estimate_correlation_parameters(y, X, beta, cluster_indices)
        
        # M-step: Update regression parameters using weighted least squares
        beta = update_regression_parameters(y, X, clusters, alpha, cluster_indices)
        
        # Check convergence
        if max(abs(beta - beta_old)) < tol:
            break
    
    # Compute robust variance estimate (sandwich estimator)
    self.beta = beta
    self.robust_cov = compute_sandwich_variance(y, X, clusters, cluster_indices)
    return self

def estimate_correlation_parameters(self, y, X, beta, cluster_indices):
    # Compute residuals from current beta estimates
    eta = X @ beta
    mu = link_function(eta)
    residuals = y - mu
    
    if self.correlation == 'exchangeable':
        # Method of moments for exchangeable correlation
        sum_products = 0
        total_pairs = 0
        
        for cluster_idx in cluster_indices.values():
            if len(cluster_idx) > 1:
                cluster_residuals = residuals[cluster_idx]
                # Compute all pairwise products
                for i in range(len(cluster_residuals)):
                    for j in range(i + 1, len(cluster_residuals)):
                        sum_products += cluster_residuals[i] * cluster_residuals[j]
                        total_pairs += 1
        
        alpha = sum_products / total_pairs if total_pairs > 0 else 0.0
        return clip(alpha, -0.99, 0.99)  # Keep correlation in valid range

The Sandwich Estimator: GEE’s Secret Weapon

The sandwich estimator is the mathematical innovation that makes GEE robust to correlation misspecification. Its name comes from its mathematical form, but its importance comes from what it accomplishes:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

def compute_sandwich_variance(self, y, X, clusters, cluster_indices):
    # Current fitted values and derivatives
    eta = X @ self.beta
    mu = link_function(eta)
    d_mu_d_eta = derivative_of_link(eta)
    var_mu = variance_function(mu)
    
    # Build matrices for sandwich formula: (X'WX)^-1 (X'ΩX) (X'WX)^-1
    
    # Working weight matrix W  
    W = diagonal_matrix(d_mu_d_eta**2 / var_mu)
    
    # Information matrix (bread)
    information_matrix = X.T @ W @ X
    information_inv = inverse(information_matrix)
    
    # Empirical variance matrix Ω (meat) - cluster-robust
    empirical_variance = zeros_matrix(X.shape[1], X.shape[1])
    
    for cluster_idx in cluster_indices.values():
        X_cluster = X[cluster_idx]
        residuals_cluster = (y - mu)[cluster_idx]
        
        # Each cluster's contribution to empirical variance
        cluster_score = X_cluster.T @ residuals_cluster
        empirical_variance += outer_product(cluster_score, cluster_score)
    
    # Sandwich formula: bread × meat × bread
    sandwich_cov = information_inv @ empirical_variance @ information_inv
    return sandwich_cov

The sandwich estimator gets its name from its mathematical form:

$(X'WX)^{-1}(X'\Omega X)(X'WX)^{-1}$

The middle term $\Omega$ captures the true covariance structure empirically, while the outer terms use the working correlation structure.

This is why GEE is robust – even if our working correlation is wrong, the empirical middle term corrects for it asymptotically. It’s mathematical insurance against our modeling assumptions.

We Observe Different Coefficients!

Suppose we run mixed effects and GEE with the same data and identical regressors and grouping. We are very likely to end up with different – if slightly – results.

Here’s where our understanding about the mixed effects vs. GEE distinction becomes critical for practical applications. The difference isn’t just methodological – it’s interpretational and has real implications for how we apply our results. In certrain circumstances, this difference may have real implications for policy decisions. If we’re designing system-wide clinical protocols, we want the GEE interpretation. If we’re making predictions for specific patients at specific facilities, we want the mixed effects interpretation.

Now you may ask: Then what do the coefficients in our regression table actually mean?

What Regression Tables Actually Tell Us

In academic research, we’re often less interested in prediction and more interested in understanding relationships – the marginal effect of a key variable while controlling for others. This is where the ME vs. GEE distinction becomes absolutely critical for interpretation.

Consider the following hypothetical example: We want to understand the effect of a new diabetes medication on HbA1c levels across our health system.

Mixed Effects Regression Table Interpretation

1
2
3
4
5

Fixed Effects:
                    Estimate    Std.Error    t-value    p-value
(Intercept)         7.2         0.15         48.0       <0.001
New_Medication     -0.8         0.12         -6.7       <0.001
Age                 0.02        0.005         4.0       <0.001

What -0.8 means in ME: “For the average group (zero random effect), switching from old to new medication is associated with a 0.8 point reduction in HbA1c, controlling for age.”

This is a conditional effect – the average marginal effect conditioned on average group characteristics. We decouple individual-level effects from group-level variation, then average them separately.

GEE Regression Table Interpretation

1
2
3
4
5

Population-Average Effects:
                    Estimate    Robust.SE    Z-value    p-value
(Intercept)         7.1         0.18         39.4       <0.001
New_Medication     -0.7         0.14         -5.0       <0.001
Age                 0.02        0.006         3.3       0.001

What -0.7 means in GEE: “Across our entire health system, patients receiving the new medication have, on average, 0.7 points lower HbA1c than those receiving the old medication, controlling for age.”

This is a marginal effect – averaged across the actual distribution of hospitals and their varying characteristics, without decoupling individual and group effects.

Why The Numbers (May) Differ

The coefficients differ (-0.8 vs -0.7) because of how we handle the averaging:

  • ME: Average the marginal effect conditioning on average group characteristics (u=0)
  • GEE: Average across the actual distribution of group characteristics

In ME, we separate individual-level relationships from group-level variation, then examine the individual-level relationship at the group average. In GEE, we average everything together without this decoupling.

Clinical and Policy Implications

These different interpretations have real consequences:

For Clinical Decision-Making (ME interpretation): “If patients are treated at an average hospital, I expect their HbA1c to improve by 0.8 points on average.”

For Health System Policy (GEE interpretation): “If we switch all eligible patients system-wide to the new medication, we expect average HbA1c to improve by 0.7 points across all facilities.”

The policy decision depends on the population-average effect. The clinical decision might depend on the hospital-specific conditional effect.

When The Differences Matter Most

The divergence between conditional and marginal effects is largest when:

  • Random effects have large variance (high between-cluster heterogeneity)
  • Non-linear link functions are used (logistic, Poisson, etc.)
  • Effect modification exists (treatment effects vary by cluster characteristics)

In linear models with random intercepts only, conditional and marginal effects are often similar. But in generalized linear mixed models, they can be substantially different.

This is why many epidemiological studies prefer GEE for exposure effect estimation – the coefficients directly answer the population-level questions that drive public health recommendations.

The Convergence Advantage: Why GEE Almost Always Works

Before we dive deeper into applications, there’s a critical practical advantage of GEE that every applied researcher needs to understand: GEE almost always converges and produces results, while mixed effects models frequently struggle with real-world data.

If we’ve ever spent hours debugging why our lme4 model won’t converge, or why SAS PROC MIXED throws cryptic warnings about non-positive definite G matrices, we’ve experienced this firsthand. Meanwhile, GEE implementations typically converge smoothly and quickly.

Why GEE Converges When Mixed Effects Fail

The theoretical reasons reveal fundamental differences in computational complexity:

1. Quasi-Likelihood vs. Full Likelihood
Mixed effects models maximize complex likelihood functions that integrate over random effects distributions. These high-dimensional optimization problems often have flat regions, multiple local optima, or numerical instabilities. GEE uses quasi-likelihood methods that avoid integration entirely – the optimization landscape is much more well-behaved.

2. Parameter Space Complexity
Mixed effects models simultaneously estimate fixed effects, random effects, and variance components – with constraints ensuring variance components remain positive. GEE only estimates fixed effects and simple correlation parameters, with fewer constraints and a more tractable parameter space.

3. Distributional Robustness Translates to Numerical Robustness
When our random effects aren’t actually normally distributed (spoiler: they often aren’t in real healthcare data), mixed effects likelihood surfaces can become pathological. GEE’s distribution-free approach maintains numerical stability even when underlying assumptions fail.

4. Identifiability Guarantees
Complex random effects structures can create identifiability issues – multiple parameter combinations yielding identical likelihoods. GEE’s focus on marginal means and simple correlation structures rarely encounters these problems.

The Speed Advantage

Beyond convergence reliability, GEE is typically 3-10x faster than mixed effects models on identical datasets. The quasi-likelihood updates are computationally simpler than the Newton-Raphson or EM algorithms required for mixed effects likelihood maximization.

This speed difference compounds when we need to:

  • Fit models across multiple bootstrap samples
  • Conduct cross-validation with clustered data
  • Explore different model specifications
  • Scale to datasets with hundreds of clusters

We’ll explore the deep mathematical reasons behind these convergence properties – including the theory of estimating equations, asymptotic behavior of quasi-likelihood methods, and numerical optimization landscapes – in upcoming posts on the computational foundations.

Beyond Continuous Outcomes

GEE really shines when we extend to generalized linear models. The framework extends naturally to logistic regression for binary outcomes. Let me sketch how:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

class LogisticGEE(GEEFromScratch):
    def __init__(self, correlation='exchangeable'):
        super().__init__(family='binomial', link='logit', correlation=correlation)
    
    def update_regression_parameters(self, y, X, clusters, alpha, cluster_indices):
        # Newton-Raphson update for logistic regression
        eta = X @ self.beta
        mu = link_function(eta)
        
        score = zeros_vector(X.shape[1])
        hessian = zeros_matrix(X.shape[1], X.shape[1])
        
        for cluster_idx in cluster_indices.values():
            cluster_size = len(cluster_idx)
            X_cluster = X[cluster_idx]
            y_cluster = y[cluster_idx]
            mu_cluster = mu[cluster_idx]
            
            # Working correlation matrix
            R = exchangeable_correlation(cluster_size, alpha)
            
            # Variance matrix  
            V = diagonal_matrix(mu_cluster * (1 - mu_cluster))
            
            # Working covariance
            W = inverse(V @ R @ V)
            
            # Score and Hessian contributions
            residuals = y_cluster - mu_cluster
            score += X_cluster.T @ W @ residuals
            hessian += X_cluster.T @ W @ V @ X_cluster
        
        # Newton-Raphson update
        beta_update = solve_linear_system(hessian, score)
        return self.beta + beta_update

A Strategic Decision Framework: When to Choose What

The choice between mixed effects and GEE isn’t just methodological – it’s strategic. Our research objectives should drive our modeling decisions.

Choose Mixed Effects When:

1. Individual Prediction Matters We want subject-specific predictions that account for both population-level relationships and group-specific deviations.

2. Understanding Group Heterogeneity is the Goal We want to identify which hospitals have systematically different outcomes, understand the distribution of facility effects, or make predictions conditional on facility characteristics.

3. Hierarchical Structure is Well-Understood We have strong theoretical reasons to believe in specific random effects structures, or we’re willing to invest in model selection across different hierarchical specifications.

Choose GEE When:

1. Population-Average Policy Questions We need marginal effects that directly inform population-level interventions and policy decisions.

2. Robustness is Critical We’re unsure about distributional assumptions, have limited theoretical guidance about correlation structures, or need results that are robust to model misspecification.

3. Computational Efficiency and Reliability Matter GEE often converges faster and scales better to large datasets than mixed effects models, especially with complex random effects structures.

The Subtle Art of Working Correlations

One of GEE’s most elegant features is that we can be “wrong” about correlation and still get correct population-average estimates. But being “more right” improves efficiency:

Key Insight: All correlation structures give consistent estimates of the treatment effect, but the correctly specified structure provides the most efficient estimates. We trade some efficiency for robustness, but we maintain consistency.

Implementation Nuances That Matter in Practice

Handling Unbalanced Clusters

Real healthcare data rarely has perfectly balanced clusters. Our implementation needs to handle this gracefully. Let me show one example of some possible approaches.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

def handle_unbalanced_clusters(self, y, X, clusters):
    # Adjust correlation estimation for varying cluster sizes
    weighted_alpha = 0
    total_weight = 0
    
    for cluster_id in unique_clusters:
        cluster_size = count_observations_in_cluster(cluster_id)
        
        if cluster_size > 1:
            # Weight by number of pairs in cluster
            n_pairs = cluster_size * (cluster_size - 1) / 2
            cluster_alpha = estimate_cluster_correlation(cluster_id)
            
            weighted_alpha += n_pairs * cluster_alpha
            total_weight += n_pairs
    
    return weighted_alpha / total_weight if total_weight > 0 else 0

Ensuring Numerical Stability

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

def ensure_numerical_stability(self, correlation_matrix):
    # Ensure correlation matrix is positive definite
    eigenvals, eigenvecs = eigendecomposition(correlation_matrix)
    
    # Clip negative eigenvalues
    eigenvals_clipped = maximum(eigenvals, 1e-8)
    
    # Reconstruct matrix
    stable_correlation = eigenvecs @ diagonal_matrix(eigenvals_clipped) @ eigenvecs.T
    
    # Ensure diagonal is 1
    diagonal(stable_correlation) = 1.0
    return stable_correlation

Beyond Standard Applications: Modern GEE Extensions

Like we discussed in the previous article, understanding GEE mechanics also enables powerful extensions that go beyond traditional statistics packages:

Time-Varying Correlation Structures

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

class AdaptiveGEE(GEEFromScratch):
    def estimate_time_varying_correlation(self, y, X, time_var, clusters):
        # Allow correlation to change over time
        correlation_params = {}
        
        for time_period in unique_time_periods:
            period_mask = time_var == time_period
            period_alpha = estimate_correlation_for_period(period_mask, clusters)
            correlation_params[time_period] = period_alpha
        
        return correlation_params

Regularized GEE for High-Dimensional Settings

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

class RegularizedGEE(GEEFromScratch):
    def penalized_score_function(self, beta, y, X, clusters, lambda_reg):
        # Modified score function with regularization
        standard_score = compute_standard_score(beta, y, X, clusters)
        
        # Add L1 penalty (LASSO)
        l1_penalty = lambda_reg * sign(beta)
        l1_penalty[0] = 0  # Don't penalize intercept
        
        return standard_score - l1_penalty

Machine Learning Integration

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

def gee_inspired_feature_engineering(df, outcome_col, cluster_col, feature_cols):
    # Create population-average features using GEE principles
    features = {}
    
    for feature in feature_cols:
        # Fit simple GEE for this feature
        simple_gee = GEEFromScratch()
        X = create_design_matrix([intercept, feature])
        
        simple_gee.fit(y=outcomes, X=X, clusters=clusters)
        
        # Use GEE coefficient as population-average feature importance
        features[f'{feature}_gee_importance'] = simple_gee.beta[1]
        
        # Create residuals adjusted for clustering
        predicted = X @ simple_gee.beta
        features[f'{feature}_cluster_adjusted_residual'] = outcomes - predicted
    
    return features

Embracing Uncertainty

At its core, GEE represents a fundamentally different methodological philosophy: embrace uncertainty about correlation structure while maintaining focus on population-average effects.

This philosophy has served us well across diverse applications:

  • Clinical Trial Analysis: Estimating average treatment effects across study sites without getting bogged down in site-specific modeling
  • Health Policy Research: Understanding population-level intervention impacts that translate directly to policy recommendations
  • Quality Improvement: Measuring system-wide changes while accounting for facility clustering
  • Epidemiological Studies: Robust estimation of exposure effects across populations with varying correlation structures

The implementation details we’ve explored – working correlations, sandwich estimators, iterative algorithms – aren’t just technical requirements. They’re the machinery that enables a more robust, policy-relevant approach to clustered data analysis.

The Bigger Picture

When our mixed effects model assumptions feel shaky, when we care more about population effects than individual predictions, or when computational robustness trumps distributional elegance, GEE provides a principled alternative.

But more than that, understanding both approaches – and their implementation details – gives us flexibility to adapt when standard packages fall short. We can:

  • Modify convergence criteria for domain-specific requirements
  • Integrate population-average concepts into modern ML architectures
  • Develop hybrid approaches that combine the strengths of both methods
  • Create computational optimizations for specific data characteristics

The choice between mixed effects and GEE isn’t about better or worse. It’s about matching our methodology to our research questions. Now we have both tools – and the implementation knowledge to transcend their standard limitations.

Next time we face clustered data, we won’t just reflexively reach for mixed effects. We’ll ask the crucial question: “Do we need subject-specific predictions or population-average effects?” That distinction will guide us toward the right methodological path – and the right answers for our stakeholders.

In the end, GEE teaches us that sometimes the most robust approach is the one that makes the fewest assumptions while still capturing what matters most. It’s a lesson that extends far beyond statistical modeling into the heart of how we approach uncertainty in science and policy.

The mathematical elegance is appealing, but the practical wisdom is what makes GEE indispensable in our modern toolkit for analyzing complex, clustered data.