# Using random effects to adjust for cluster-level confounding?

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.
• Standard ways to adjust for an unobserved confounder would be either front-door adjustment if you can insert a variable between the treatment and the outcome; or, if you have continuous variables, you can try an instrumental variable approach. – Adrian Keister Aug 18 '20 at 15:28
• Yes, thanks, I don't yet see how this relates to the question though? – stefgehrig Aug 18 '20 at 17:16
• Doesn't the sinulstion in your second link describe your situation very closely ? – Robert Long Aug 18 '20 at 20:16
• Yes, it does, and it seems to me that the random effects independence assumption is violated in that case. (Of course, the model can still be fit and give reasonable results) Am I right or making a mistake here? – stefgehrig Aug 19 '20 at 5:44
• I also get inconsistency in that model (by @RobertLong) depending on how strong the correlation is between X and E is, it surprises me that a mixed model does so well. – Jesper for President Aug 19 '20 at 12:18

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.

• Yes, I am still very puzzled by this because it clearly does violate the assumption of a mixed model (to my understanding of mixed models which i admit is rudimentary). I would be inclined to say that theoretically random effects as defined in the question cannot make sense of confounding. – Jesper for President Aug 19 '20 at 16:41
• @JesperforPresident just because an assumption is violated doesn't mean the model isn't useful. All models are wrong, remember. And in many cases, some assumptions only apply in small samples. Context is key. – Robert Long Aug 19 '20 at 18:39
• Thanks for the simulation, it clarifies things. In more complex cases, as you say, it probably makes sense then to compare RE and FE to have a better idea about the amount of bias. As @JesperforPresident, I am also a bit puzzled about the theoretical side of things, because the estimate will be more biased the more important (i.e., correlated with exposure) the group-level confounder is - so there are two different sources of bias weighed against each other. It would be great to see some theoretical discussion, or at least appreciation, of this topic! Literature suggestions welcome – stefgehrig Aug 20 '20 at 7:27
• I agree. Perhaps you can ask a question on the theory behind the assumption. I had a brief look at various Bates papers and did not see an explicit assumption. Implicitly it is there, perhaps. – Robert Long Aug 20 '20 at 7:48
• These authors for example are really explicit about the assumption, even calling it omitted variable bias when cluster level confounders are not included/observed but "only" the REs are specified. jyx.jyu.fi/bitstream/handle/123456789/66704/2/Antonakisym.pdf – stefgehrig Aug 21 '20 at 14:34

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.