# Puzzling predicted values in generalized multilevel model

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)[]))
(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[], 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.

• Try using predict( mod_1, random.only = TRUE) and look at estimates using only the random effects, don't the estimate look suspiciously good? Then try using defining mod_0 = glmer(y ~ (1 | id), dat, binomial), don't the estimate look even more suspiciously good? The random effect has sucked the life out of your model. Especially for person who answered a single question all their variability is transferred to the random component and treated as noise. – usεr11852 Jul 22 '18 at 20:33
• +1 As an answer, could you explain more of the reasoning behind: "The random effect has sucked the life out of your model. Especially for person who answered a single question all their variability is transferred to the random component and treated as noise." What makes this the case? How does only one observation for a person transfer to the random component? And why do you think MCMCglmm works better? – Mark White Jul 22 '18 at 21:31
• @usεr11852 I make that previous comment (forgot to tag you in it) because I tried simulating some data with the same distribution of 1, 2, and 3 responses and still got OK inferences: gist.github.com/markhwhiteii/4987ed886017eadc1785a5cdef298b3d. Granted, these simulated data were generated by exactly how the model is specified, which is undoubtedly supremely optimistic. You can see the distribution of predictions, broken out by whether or not there are people with 1, 2, or 3, 2 or 3, and all 3 observations: s22.postimg.cc/qv80p008h/… – Mark White Jul 22 '18 at 21:45
• Ultimately, the random effects part of an MEM is there because we recognise there is some (possibly id-related) structure in our errors/residuals. If we allow as single error parameter of each residual point then the remaining variability is minor. As such our fixed effects are left with next to nothing to model. This results to the very small coefficient you see. (Note that in logistic regression there is no error term so in comparison with the simple LMER the effect will be even more pronounced) – usεr11852 Jul 22 '18 at 22:00

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

• Good one (+1 from me). Before posting my comment, I tried using MCMCglmm, it too gives reasonable results when predicting. – usεr11852 Jul 22 '18 at 20:36
• +1 I'm afraid I'm still not following why 1 observation for 60% of the sample causes numeric conflicts (mainly, people use the terms random and fixed effects to mean lots of things, so I'm not sure what people are referring to often); I get in retrospect that it makes sense, but I'm not sure why. Thanks for the GEE suggestion—I'm reading up on them because I'm curious why the numerical conflicts don't arise there but do in a mixed model. – Mark White Jul 24 '18 at 1:12