huge difference between estimates of binomial regresssion when including random effect vs when not

Question

I'm trying to estimate the average score for two groups of students. I use a binomial regression model. The total_ans is the total question they've have answered, which may be different for different students.

Model 1 directly estimates

model <- glm(cbind(total_correct, total_ans-total_correct) ~ student_type,family= binomial, data = df)

Call:  glm(formula = cbind(total_correct, total_ans - total_correct) ~ student_type, family = binomial, data = df)

Coefficients:
                  (Intercept)  student_group_2  
                      -1.9684           0.2139  

Degrees of Freedom: 1552 Total (i.e. Null);  1551 Residual Null
Deviance:       1480  Residual Deviance: 1477   AIC: 1764


lsmeans(model,~ student_type, type="response")

 student_type           prob         SE df asymp.LCL asymp.UCL
 student_group_1   0.1225627 0.00654160 NA 0.1103074 0.1359715
 student_group_2   0.1474774 0.01275231 NA 0.1241918 0.1742602

In model 2, I use a random effect to better account for individual variances.

model <- glmer(cbind(total_correct, total_ans-total_correct) ~ (1|student)  + student_type, family= binomial, data = sub_df, control=glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE))

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(total_correct, total_ans - total_correct) ~ (1 | student) +  
    student_type
   Data: sub_df
      AIC       BIC    logLik  deviance  df.resid 
1653.9049 1669.9488 -823.9525 1647.9049      1550 
Random effects:
 Groups Name        Std.Dev.
 student (Intercept) 1.881   
Number of obs: 1553, groups:  author, 1553
Fixed Effects:
                  (Intercept)  student_group_2  
                      -3.0571   0.3915  


lsmeans(model,~ student_type, type="response")

 student_type            prob          SE df  asymp.LCL asymp.UCL
 student_group_1   0.04491007 0.004626728 NA 0.03666574 0.0549025
 student_group_2   0.06503249 0.015117905 NA 0.04097546 0.1017156

I'm surprised that there is such a big difference between the results in the two groups. What might be the reason for this?

more info:group 1 has 1434 students, group 2 has 119 students. these are naturally occurring groups

Answer 1

There is enough information in the question to clear this up. lsmeans is simply using the coefficients to obtain the group predicted probabilities.

So for the GLM, OP's implied model is:

\begin{equation} \hat\pi=\frac{1}{1+e^{-(-1.9684+0.2139\times d)}} \end{equation}

where $d$ is an indicator for membership in group 2. So the predicted probabilities are: $(1+e^{-(-1.9684)})^{-1}$ and $(1+e^{-(-1.9684+0.2139)})^{-1}$ for groups 1 and 2 respectively. These result in predicted probabilities of about $12.25\%$ and $14.75\%$ respectively.

For the multilevel model (or GLMM), OP's implied model is:

\begin{equation} \hat\pi=\frac{1}{1+e^{-(-3.0571+0.3915\times d+1.881\times \hat{u})}} \end{equation}

where $\hat{u}$ is the random intercept assumed to be standard normal. The predicted probabilities from lsmeans assume a random intercept value of zero $(\hat{u}=0)$ resulting in: $(1+e^{-(-3.0571)})^{-1}$ and $(1+e^{-(-3.0571+0.3915)})^{-1}$ for groups 1 and 2 respectively. These result in predicted probabilities of about $4.49\%$ and $6.50\%$ respectively. These are the lsmeans GLMM results.

An issue is the intercept in GLMM is the expected log-odds of success for someone who is "average" relative to others. So using this as a basis for reporting the whole model is an issue. About the other coefficient, one suggestion for why the group difference coefficient increases is that the model quality is better so the coefficient increases, search for collapsible on this website or see Is meta-analysis of odds ratios essentially hopeless?.

To make the lsmeans GLMM results comparable to the lsmeans GLM results. We have to use the observed values of the random intercept, $\hat{u}$. One can simulate random intercepts to match OP's specific model:

set.seed(12345)
u_hat <- rnorm(1553, 0, 1.881)
d <- c(rep(0, 1434), rep(1, 119))
# predicted log odds:
pred_log_odds <- -3.0571 + 0.3915 * d + u_hat
pred_prob <- plogis(pred_log_odds)
coef(lm(pred_prob ~ 0 + factor(d)))

factor(d)0 factor(d)1 
 0.1189244  0.1490395

In this simulated example, these values are a lot closer to the lsmeans GLM results. If you run something like the syntax below on your data:

lm(fitted(model) ~ 0 + df$student_type)

where model is the GLMM, you should get values quite close to the lsmeans values. I'm assuming that when you call fitted() on glmer(), it also includes the random intercept and the values that are returned are probabilities.

---

In your situation where the groups are naturally occurring, an additional thing to explore in the data is the differing group variances on the random intercept, so a model like:

glmer(cbind(total_correct, total_ans-total_correct) ~
  (0 + student_type || student) + student_type,
  family = binomial, data = sub_df,
  control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE))

might be worth exploring since right now, you're assuming the random intercept variances do not differ by group. I used the || so lme4 does not attempt to correlate the two random intercepts.

---

I did not know that one could add a random intercept for each student when they essentially have only one row in the data. But I am rationalizing it away by assuming that the trial vs. failures per row amount to several rows in long form.

Answer 2

Model 1, Logistic regression without random effect:

$$\mathrm{Logit}(\Pr(Y=1)) = X\beta$$

We know MLE $\hat \beta$ is asymptotically unbiased. But $$\hat Pr(Y=1) = \mathrm{Logit}^{-1}(X\hat\beta)$$ is biased estimate of $\Pr(Y=1)$, because the non-linearity of the logit function. But asymptotically it is unbiased. So for model 1, $\hat Pr(Y=1) = \mathrm{Logit}^{-1}(X\hat\beta)$ is acceptable.

Model 2, Logistic regression with random intercept:

$$\mathrm{Logit}(\Pr(Y_{ij}=1)) = X_{ij}\beta + \gamma_i$$ $$\gamma \sim N(0,\sigma^2)$$

In this situation, the common mistake is $$\Pr(Y=1|X) = \mathrm{E}(Y|X) = \mathrm{E}(\mathrm{Logit}^{-1}(X\beta + \gamma_i)) =\mathrm{Logit}^{-1}( \mathrm{E}(X\beta + \gamma_i)) = \mathrm{Logit}^{-1}(X\beta) $$

In this process the non-linearity of the logit function is totally igored.

So the following output is very misleading. It gives a idear that for group 1 student, Probability of correction is 4.5%.

    student_type            prob          SE df  asymp.LCL asymp.UCL
    student_group_1   0.04491007 0.004626728 NA 0.03666574 0.0549025
    student_group_2   0.06503249 0.015117905 NA 0.04097546 0.1017156

The mistake above cannot be contributed to mixed effect logistic regression, should come from the misinterpretation of the results.

Let look how to correctly derived the marginal mean of $Y$. $$\Pr(Y=1|X) = \mathrm{E}(Y|X) = \mathrm{E}(\mathrm{Logit}^{-1}(X\beta + \gamma_i)) =\int_{-\infty}^{\infty} \mathrm{Logit}^{-1}(z)\phi(z)dz$$ where $z=X\beta + \gamma_i$ and $\phi$ is pdf of $N(X\beta, \sigma^2)$

For Group 1 student, the results indicate that $X\beta= -3.0571$, $\sigma = 1.881$. Using Gauss–Hermite quadrature method with 20 points, I got $\Pr(Y=1|X=0) = 0.1172$. similarly, $\Pr(Y=1|X=0) = 0.1492$.