Puzzling predicted values in generalized multilevel model

Question

I ran a multilevel model with a dichotomous variable as the outcome using the logistic link function. I did this with the lme4 package in R. An approximation of my data can be found at this GitHub Gist: https://gist.github.com/markhwhiteii/3a954c4da82efdad5d8abee601b8a6a8

The design of the data collection process is pretty simple: Participants are randomly assigned to a treatment (t) or control (c) condition. They are then asked three questions; they may answer one, two, or three of them. For each of the questions (q1, q2, q3), they can respond positively (1) or negatively (0). I do not necessarily care about the effect of variable—I am considering each question as indicators of more or less the same construct.

The question is simple: Does the treatment increase the probability of responding positively (1)? Importantly, the effect size I want to use is a model-implied percentage point difference between the two groups. This requires transforming predicted values from logits to probability.

The model I specify is an intercept-only model, with the response at Level 1 (denoted by $i$) and the person at Level 2 (donated by $j$):

$\text{logit}(\pi_{ij}) = \beta_{0j}$ and $\beta_{0j} = \gamma_{00} + \gamma_{01}X_j + u_{0j}$

In the syntax of lme4, this is:

mod_1 <- glmer(y ~ x + (1 | id), dat, binomial)

(If you read.csv() in the .csv file in the Gist link above, this code will estimate the model.)

Here is where the confusion rises. I can transform the predicted values for responses in the control condition and those in the treatment condition into probability space. They give me:

inv_logit <- function(x) exp(x) / (1 + exp(x))
(control <- inv_logit(fixef(mod_1)[[1]]))
(treatment <- inv_logit(sum(fixef(mod_1))))
treatment - control

This will show about a 0.0003 probability of responding positively in the control condition as well as the treatment condition, or a .03% chance of being positive. The difference between conditions is quite small.

These conditional probabilities in each condition seem too small, especially when we look at the naive percentages of positive responses in both conditions:

with(dat, tapply(y, x, mean))

This will show that about 22% of the responses in the control condition are positive, while 23% in the treatment condition are positive. My fundamental question is: How is there such a big difference between the 22% or so positive outcomes empirically and only .03% implied from the model?

We get the same thing if we first average by person (to not allow those responding more to have more weight) and then by condition:

library(dplyr)
dat %>% 
  group_by(id) %>% 
  mutate(p = mean(y)) %>% 
  slice(1) %>% 
  group_by(x) %>% 
  summarise(p = mean(p))

Both conditions show about 24% of the responses being positive. Again, this is drastically different than the .03% implied by mod_1.

Some other information that may be helpful:

• 604 people had 1 response, 301 had 2 responses, and 95 had 3 responses. I think part of the issue is that the majority of the people only had one observation? As far as I can tell, though, we don't need to estimate a different normal distribution of responses for each cluster (in this case, a person), so I don't get why only 1 response would be giving me hang-ups.

• The distribution of random intercepts was decidedly not normal as can be seen by hist(ranef(mod_1)$id[[1]], breaks = 50):

The model clearly seems to be misspecified from that horrific distribution of random intercepts, but I would like to know:

1. Why does my naive look at that about 22-24% of the responses being positive differ so much from the model-implied .03%? I'm looking for an answer beyond "the model is misspecified," because I am curious as to why this is occurring, plus that doesn't get me any closer to specifying a proper model. Which leads me to...

2. How can I specify a proper model that answers my question regarding the experimental condition? I suppose I could also try to use some type of sandwich covariance estimation that handles dependent responses, such as sandwich::vcovCL()? That approach doesn't change the SEs too much in the current circumstance. This also doesn't satisfy my curiosity as to why the mixed model is giving me strange predictions.

Answer 1

The answer might at least partly be hidden in your answer:

The predicted random effects are all non-negative. Adding them to the linear predictor of the fixed effects will push up the average prediction to the desired level.

So what could be the reason for the non-centered distribution of the eBLUPs? I think it is hidden in the extremely unbalanced data situation: 604 of the 1000 ids provide just one single measurement, which can easily lead to numeric conflicts between random and fixed effects.

In such situations, I often run as a sensitivity analysis a normal mixed-model. If we do this with your data, we get:

mod_2 <- lmer(y ~ x + (1 | id), dat)
summary(mod_2)

# Output
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.06961  0.2638  
 Residual             0.10825  0.3290  
Number of obs: 1491, groups:  id, 1000

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.228315   0.022716  10.051
xt          0.007224   0.026922   0.268

Correlation of Fixed Effects:
   (Intr)
xt -0.844

The fixed effects now look as expected by the descriptive analysis and very different from what the GLMER found. Now, the eBLUPs are even centered (while of course still far away from a normal):

summary(unlist(ranef(mod_2, drop = TRUE)))

# Output
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.15513 -0.12844 -0.09218  0.00000  0.14878  0.50347 

Depending on your exact research questions and goals of your analysis, you could try generalized estimation equations GEE, in R:

library(gee)
fit_gee <- gee(y ~ x, 
               id = id,         
               data = dat %>% arrange(id), 
               family = binomial, 
               corstr = "exchangeable")
summary(fit_gee)

# Output
Coefficients:
               Estimate Naive S.E.    Naive z Robust S.E.   Robust z
(Intercept) -1.22561668  0.1257257 -9.7483411   0.1235148 -9.9228358
xt           0.04325299  0.1484306  0.2914022   0.1479156  0.2924167

Estimated Scale Parameter:  0.9859774
Number of Iterations:  3

Working Correlation
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.3133096 0.3133096
[2,] 0.3133096 1.0000000 0.3133096
[3,] 0.3133096 0.3133096 1.0000000

# Distribution of predictions
summary(fitted(fit_gee))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2269  0.2269  0.2346  0.2324  0.2346  0.2346