Using random effects to adjust for cluster-level confounding?

Question

There is a usage of random intercepts to adjust for unobserved cluster-level confounding, as for example argued here:

Are random effects confounding variables?

How do random effects adjust for confounding in a model?

Based on this advice and examples from the literature in a similar spirit, one would imagine that random effects can be used for adjustment in a DAG like this, where there is an unobserved confounder on the cluster level:

For example, imagine a clinical study where hospitals differ in their propensity to enroll high-risk patients (more likely to experience the adverse outcome) and also in their propensity to give the treatment under study, due to an unobserved structural characteristic.

On the other hand, a core assumption of random effects models is that the predictor (here: Treatment) is uncorrelated with the random intercepts, see for example Verbeek (2008):

"...it may be the case that $𝛼_i$ [random effects] and $x_{it}$ [predictor] are correlated, in which case the random effects approach, ignoring this correlation, leads to inconsistent estimators. We saw an example of this previously, where $𝛼_i$ included management quality and was argued to be correlated with the other inputs included in the production function. The problem of correlation between the individual effects $𝛼_i$ and the explanatory variables in $x_{it}$ can be handled by using the fixed effects approach, which essentially eliminates the $𝛼_i$ from the model, and thus eliminates any problems that they may cause."

or Setodji and Shwartz (2013):

"...base their choice of model type on whether unobserved time-invariant omitted variables, which are captured in $\phi_j$ [random effects], are uncorrelated with the main predictor of interest. If uncorrelated (an assumption that can be assessed using the Hausman test), random-effect models are appropriate; otherwise, fixed-effect models are used."

If, by definition, a confounder is correlated with exposure, and random effects models assume uncorrelatedness of random effects and exposure, how can random effects then be used to adjust for confounding?

References

• Verbeek, M. (2008). A guide to modern econometrics. John Wiley & Sons.
• Setodji, C. M., & Shwartz, M. (2013). Fixed-effect or random-effect models: what are the key inference issues?. Medical care, 51(1), 25-27.

Answer 1

The thing about assumptions, is that they are there to be violated. It is rare, if not impossible in observational studies for 2 variables to have a correlation of zero. Correlation is expected, even if it is just due to random sampling and not confounding or some other causal mechanism. The interesting questions are: to what extent is an assumption volitated, and how robust is a particular model to such violations. The first point is subjective and the latter can be quite difficult to establish in all but simple models. As usual simulation can be your friend, so let's take a look using your example:

Here we will simulate data so that the confounder X is highly correlated with the exposure E, with correlations ranging from 0.5 to 0.95

set.seed(15)
N <- 100
n.sim <- 100
simvec.E <- numeric(n.sim)
rhos <- seq(0.5, 0.95, by = 0.05)
simvec.rho <- numeric(length(rhos))

for (j in 1:length(rhos)) {

  Sigma = matrix(c(1, rhos[j], rhos[j], 1), byrow = TRUE, nrow = 2)

  for(i in 1:n.sim) {
    dt <- data.frame(mvrnorm(N, mu = c(0,0), Sigma = Sigma, empirical = TRUE))  

    # put them on a bigger scale, so it's easy to create the group factor
    dt1 <- dt + 5
    dt1 <- dt1 * 10
  
    X <- as.integer(dt1$X1)

    E <- dt1$X2

    Y <- E + X + rnorm(N)  # so we expect estimate for E that we want to recover is 1
  
    X <- as.factor(X) 
    lmm <- lmer(Y ~ E + (1|X))
    simvec.E[i] <- summary(lmm)$coef[2]
  }
  simvec.rho[j] <- mean(simvec.E)
}

ggplot(data.frame(rho = rhos, E = simvec.rho), aes(x = rho, y = E)) + geom_line()

This produces:

So, yes, the is some bias introduced when the correlation becomes large, but at correlations below 0.85 or so, this is fairly negligable. In other words the mixed model seems quite robust. Note that the way I simulated the grouping factor here leads to quite small cluster sizes. Increasing N will lead to larger clusters, though this takes longer to run of course. With N <- 1000 I get:

which is a considerable improvement. Of course we could also look at standard errors, and other sample size / designs, random slopes etc. but I'll leave that for another day.

With real data where this problem arose I would always want to compare a fixed effects model as well as random effects.

Answer 2

A random-effects model does not control for unobserved invariant unit-level heterogeneity ($\alpha_i$ in your excerpt from Verbeek). If your intention is to make causal claims from the model and you have reasons to believe that $\alpha_i$ is correlated with the causal variable of interest, your model will be rejected by the scientific community because it is not the best possible evidence on the issue. Why? Because if you can run a random-effects model, it implies you have multiple observations for the same unit. In such a situation, you can easily adjust for $\alpha_i$ and thus you did not produce the best possible evidence for the question in hand.

To fix ideas, assume your models is: $y_{it} = \beta_0 + B_1 X_{it} + \beta_2 D_{it} + \alpha_i + \epsilon_{it}$

Assume that $i$ represents the unit and $t$ represents the time period, $y_{it}$ is the observed outcome for unit $i$ at time $t$, $X_{it}$ is a vector of covariates, $D_{it}$ is the causal variable, which varies over time for some units, and $\alpha_i$ is the time-invariant unobserved heterogeneity. The quantity we are interested in estimating is $\beta_2$, which is the treatment effect. Further, assume that $\alpha_i$ is correlated with $D_{it}$. One easy fix for $\alpha_i$ is to take the difference between two observations for each unit and use it to estimate the model (this time without $\alpha_i$, which gets differenced out).

$\Delta y_{it} = B_1 \Delta X_{it} + \beta_2 \Delta D_{it} + \Delta \epsilon_{it}$

Now, we can consistently estimate $\beta_2$ assuming that we have no unmeasured confounding conditional on $X$. The cost for first differencing is the loss of observations but we get the gain far outweighs the cost.