When are cluster robust standard errors a valid "alternative" to mixed models? 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?

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.
