Understanding how to tell whether random effects assumption is sufficiently violated to pose a problem in practice

Question

Consider a situation where I want to predict a binary health outcome for patients with various medical conditions, who are treated at different hospitals. I want to use patients' medical conditions as predictors, and it is the coefficients & confidence intervals for these conditions that I care most about (i.e. I don't care so much about differences between these specific hospitals). So it seems like a perfect situation to treat hospital as a random effect (random-intercepts model), e.g.

glmer(outcome ~ conditionA + conditionB + ... + (1 | hospital), family="binomial")

However, patients' medical conditions may very well be correlated with hospital, because patients in the most serious condition are more likely to be sent to some hospitals than others. The amount of multicollinearity here is not super-strong -- the VIF of 'hospital' in a model where hospital is treated as a fixed effect is 3.25 if all potentially relevant conditions are included as independent variables, and under 1.5 if LASSO or stepwise regression is used to exclude nuisance variables -- but it's not nothing.

With that background, I am trying to determine whether it makes more sense to treat hospital as a random or fixed effect in this instance. As noted in this question,

The random effects assumption is that the individual unobserved heterogeneity is uncorrelated with the independent variables. The fixed effect assumption is that the individual specific effect is correlated with the independent variables.

An answer to how to test for this advised extracting the random effects in R via ranef and "plotting them against the predictors". To be clear: is it true that in my case, a recommended approach would be to get the random effect for each hospital from a model where hospital is treated as a random effect; then to get the coefficients for each hospital from a model where hospital is treated as a fixed effect; and then checking if the hospital (fixed-effect) coefficients are significantly correlated with the random effects? Is this roughly equivalent to conducting a Hausman test to decide whether to treat a variable as a fixed or random effect, as described on slide 16 here?

Finally, irrespective of correlations between observed variables, is this a situation for a fixed-effects rather than random-effects model merely on the theoretical grounds that there are unobserved variables underlying the fact that people with more serious medical conditions are more likely to wind up at some hospitals than others (e.g., differential availability of acute services at different hospitals)?

Answer 1

In my experience, the issue of the correlation of predictors / exposures with the random effects only becomes a problem when

1. the correlation is very high - typically in the region of 0.8 or higher.

2. when the cluster sizes are small.

3. when the goal of the analysis is inference rather than prediction.

Regarding 1, in healthcare settings, this is fairly implausible.

Regarding 2, even with small cluster sizes, mixed models are quite robust as we will see from the simulations below

Regarding 3, you specifically mention prediction as the goal of your analyis so again, we will see below that predictions from mixed models with correlated fixed and random effects are not greatly affected by the degree of corelation.

It is also worth noting here, that in this kind of applied setting, we are not talking about a problem of confounding - it is mediation. The exposure causes the outcome, and also the group (hospital) assignment, and the hospital has a causal effect on the outcome. So, in a causal framework if we were interested in the total effect of the exposure on the outcome we would not adjust for the hospital effect, either as fixed effects or random effects, but we would do so if we were only interested in the direct effect. Again, if we are interested in prediction instead, rather than inference, then this problem wanes.

So here is a simple simulation were we look at varying levels of correlation between an exposure E and grouping variable X from 0.5 to 0.95 and we look at the impact of this on the estimate for E and the mean squared error of predictions:

library(MASS)
set.seed(15)
N <- 100
n.sim <- 100
simvec.E <- numeric(n.sim)          # a vector to hold the estimates for E
simvec.mse <- numeric(n.sim)        # a vector to hold the mse for the predictions
rhos <- seq(0.5, 0.95, by = 0.05)
simvec.rho <- numeric(length(rhos))    # vector for the mean estimates at each rho
simvec.rho.mse <- numeric(length(rhos))  # vector for mse at each rho

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 the 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.mse[i] <- sum((Y - predict(lmm))^2)
  }
  simvec.rho[j] <- mean(simvec.E)
  simvec.rho.mse[j] <- mean(simvec.mse)
}

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

So here we see that the estimates for E (simulate with a value of 1) are largely unbiased up to correlations of around 0.8. Even at 0.95 the bias is only 6%

Here we see no marked effect on mean squared error of prediction.

As mentioned above, small cluster sizes exacerbate the bias. In these simulations each dataset has only 100 observations with 35-40 groups, so the cluster sizes are small.

We can easily create more clusters by increasing N to 1000 which results in around 50-60 groups

Here we see that the bias is smaller.

And here again we see no discernable impact of correlation on mean squared error of prediction.

I would encourage you to play around with these or similar simulations, there are many parameters that can be changed, as well as changing the way the data are simulated to better reflect your actual use case.