How do random effects adjust for confounding in a model? How do random effects adjust for confounding in a model? This post explains that they do just as fixed effects do. This makes sense intuitively. However, I don't understand how the inclusion of a random effect can change the estimated effect of a treatment on the outcome.
If a random intercept is assumed to have a mean of zero, wouldn't the other coefficients in the model be unaffected (only their standard errors would change)?
To make this more concrete: I want to estimate the effect of a medical treatment on death. I believe that the medical center is a confounder (the particular hospital affects both likelihood of the subject receiving the treatment and likelihood of the subject dying). Will including hospital as a random intercept adjust the estimated coefficient for treatment even though the expected value of the random effect for hospital is 0?
 A: I agree that this is can be a little confusing.

Will including hospital as a random intercept adjust the estimated coefficient for treatment even though the expected value of the random effect for hospital is 0?

Yes. Just because the random intercepts have a mean of zero does not mean that it doesn't control for confounding. Very often an analyst has difficulty deciding whether to model a factor as fixed or random. Very often there are different arguments in favour of both, but there is no argument as to whether one adjusts for confounding or not.
When we have clustered data with correlations within the clusters, we can control for this by eithe using a random intercept for the cluster ID, or fitting a fixed effect for the cluster ID [Generalised estimating equations are another option, but that's not relevant to this answer].
A simple simulation shows this:
set.seed(15)

n <- 50

X <- rbinom(n, 10, 0.5)

E <- (X/5) + rnorm(n)

Y <- E + X + rnorm(n)

Here we have an exposure E and an outcome Y, but the association is confounded by X. The "true" value for the estimate for E is 1:
X <- as.factor(X) 
lm1 <- lm(Y ~ E)
lm2 <- lm(Y ~ E + X)  
lmm <- lmer(Y ~ E + (1|X))


> summary(lm1)

E             1.5232

where I have omitted all but the the essential output. Evidently this is confounded.  But, if we include the confounder X as a fixed effect we get:
> summary(lm2)

E            1.0446

as expected. And we also we find
> summary(lmm)

E             1.0661

