Lack of within-cluster variability

Question

I am working on patients' data. I want to do multilevel logistic regression. The cluster is hospital, exposure variable is treatment (A, B, C), and independent variables include sex, age and others. I am concerned that some hospitals (i.e. clusters) may use one treatment for all their patients. What would happen in this case to the results? How can I test it? and how can I overcome this problem?!

Answer 1

Assuming you are using a mixed model, i.e. model clusters with random effects, even if every hospital uses only one treatment it will be fine, as long as you have enough hospitals. (See here: https://stats.stackexchange.com/a/479806/341520). If you don't have enough you might end up with convergence problems or uselessly wide confidence intervals.

What is enough depends on a lot. A general recommendation is more than 20 (https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#should-i-treat-factor-xxx-as-fixed-or-random), but I recommend simulation. Even with 30 hospitals 10 patients per hospital and admittedly large random effect variance the model does not always converge:

Setup with linear response followed by binomial response. Model converges

library(lme4)
set.seed("000")
treatment <- rep(c("A", "B", "C"), each = 100)
hospital  <- rep(1:30, each = 10)
plot(hospital, as.factor(treatment))

random_effect <- rnorm(30)
y <- ifelse(treatment == "C", 2, 0) + random_effect[hospital] + rnorm(300)

my_mod <- lmer(y ~ treatment + (1|hospital))
summary(my_mod)
confint(my_mod)


y2 <- rbinom(300, 1, prob = plogis(ifelse(treatment == "C", 2, 0) + random_effect[hospital]))
my_mod <- glmer(y2 ~ treatment + (1|hospital), family=binomial(link = "logit"))
summary(my_mod)
confint(my_mod)

Result on y2

> summary(my_mod)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: y2 ~ treatment + (1 | hospital)

     AIC      BIC   logLik deviance df.resid 
   350.7    365.5   -171.3    342.7      296 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9749 -0.7596  0.2678  0.6439  1.5497 

Random effects:
 Groups   Name        Variance Std.Dev.
 hospital (Intercept) 0.6944   0.8333  
Number of obs: 300, groups:  hospital, 30

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.09238    0.34366  -0.269    0.788    
treatmentB  -0.13102    0.48225  -0.272    0.786    
treatmentC   2.26221    0.55826   4.052 5.07e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) trtmnB
treatmentB -0.712       
treatmentC -0.619  0.437
> confint(my_mod)
Computing profile confidence intervals ...
                 2.5 %    97.5 %
.sig01       0.4195519 1.3892689
(Intercept) -0.8129415 0.6209828
treatmentB  -1.1374003 0.8717255
treatmentC   1.2119191 3.5084389

Different seed and it doesn't converge

set.seed("004")
random_effect <- rnorm(30)
y2 <- rbinom(300, 1, prob = plogis(ifelse(treatment == "C", 2, 0) + random_effect[hospital]))
my_mod <- glmer(y2 ~ treatment + (1|hospital), family=binomial(link = "logit"))

warning:

In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.0233217 (tol = 0.002, component 1)