When are cluster robust standard errors a valid "alternative" to mixed models?

Question

Background

If we have data with multiple measurements per subject, I understand we must model or correct for the non-independence. I believe there are several approaches to doing this. One is a mixed effect model using subject ID as a random effect. Another models subject ID as a fixed effect (i.e. the "fixed effects model" in panel data terminology). Another involves fitting a regression model ignoring the clustering (i.e. a working independance model with no variable for ID), but correcting the standard errors after fitting the model, i.e., computing cluster robust standard errors. Although the latter approach corrects the standard errors, it does not correct the point estimates, i.e., estimated parameters. Therefore I cannot understand how this is a valid approach for estimating paramaters if, for example, some subjects appear more times than others in a dataset.

Example scenario

Supposing one subject with particularly high values of the response variable (Y) for a given predictor, appeared substantially more times in the dataset. Would this not cause the corresponding parameter estimate to be biased high? In the fixed or random effects model (explicitly modelling subject ID), I don't think this occurs since each subject has its intercept based on its average of Y (I think, though I might be missing some technical details). I'm aware that clustered standard errors can also be used after a fixed or random effects model.

Main question

If we have multiple measurments per subject, when is ONLY using the cluster standard errors approach valid i.e. NO fixed or random effect variable for subject ID, but correcting the SE's for clustering after model fitting? I am particularly concerned about bias in parameter estimates and any possible differences in approach w.r.t the type of regression model e.g. ols, Poisson, logistic, Cox etc.

Thoughts...

Perhaps we assume that, on average, the parameter estimate for Y will be unbiased since the subjects sampled more than once will (in the long run) balance out in terms of "low" and "high" Y values. I'm not sure if there is a way to "check" this after model fitting with the usual diagnostic plots?

Answer 1

The answer from @Eoin gets to the heart of the matter. One models marginal, population-level values while the other models conditional, subject-level values.

Here are a few additional suggestions, drawing from the McNeish et al. paper recommended by @Noah in a comment, "On the unnecessary ubiquity of hierarchical linear modeling," Psychological Methods 22: 114–140 (2017), and Chapter 7 of Frank Harrell's course notes and book, on modeling longitudinal responses.

First, your concern about a highly influential point leading to bias is an issue in all regression. If some unmodeled characteristic leads to an individual with "particularly high values of the response variable (Y) for a given predictor" in standard linear regression with equal numbers of observations per individual, you have a similar problem.

The fundamental assumption is that the regression model is capturing the associations of predictors with outcome, plus an error term to be estimated from the model in least squares or set by the family in a generalized linear model. If that assumption isn't met, then an individual with an unmodeled extra association with outcome could lead to undue influence on the coefficient estimates. Whether that ends up being more or less of an issue with multiple measurements per subject has to do with the relative number of observations associated with the subject. Checking influence is important, regardless.

Second, a hierarchical model doesn't necessarily avoid that problem. Remember that a strength of multi-level modeling is its partial pooling of observations among cases to obtain its estimates--cases with more observations contribute more to the model. If the high response value Y of an individual is due to an individually high slope with respect to predictor X but you only model random intercepts, you have much the same problem.

Third, a working independence model followed by robust standard errors isn't the only way to deal with repeated measures outside of a mixed linear regression model. Generalized least squares, unlike a working independence model, can provide best linear unbiased estimates with correlated observations, provided that you specify the correct correlation structure. The coefficient estimates will in general differ from those under a working independence assumption. Harrell discusses this in his chapters on longitudinal modeling and includes a helpful summary table of strengths and weaknesses of different approaches.

Although generalized least squares doesn't work other than with linear regression, the McNeish et al paper nicely outlines the approach of generalized estimating equation (GEE) models, which can start with several choices of initial correlation structure and iteratively account for correlations in the broader context of generalized linear models.

So there are several valid approaches, beyond mixed models or working independence followed by robust standard errors, for dealing with multiple measurements per subject. They are appropriate choices when the primary focus is on marginal, population-average effects and there isn't a multi-level hierarchy.

Answer 2

Good question!

I think the short answer is that the working independence linear model and the random effects model by design estimate different quantities, and in your example these quantities are not the same. In the simplest example...

m = lm(y ~ 1, data = data)
mixed_m = lmer(y ~ 1 + (1|id), data = data)

$$ \begin{align} \text{m: }&y \sim Normal(\mu_y, \sigma_y); \\ \text{mixed_m: }&y_s \sim Normal(\mu_s, \sigma) \\ &\mu_s \sim Normal(\mu_p, \sigma_p); \end{align} $$

(s for subject, p for population)

...the simple linear model m estimates the average across all observations ($\mu_y$), while the mixed model mixed_m estimates the average of each of the subject-level averages ($\mu_p$). As you say, these aren't necessarily the same thing, but that's doesn't mean one model is "biased", it's just being used to answer the wrong question.