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)

                  (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) +  
   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

  • 3
    $\begingroup$ @PeterFlom that's just a way to specify number of successes, number of failures in the glm. the sum of the two is the total trials $\endgroup$ – user926321 Nov 1 '18 at 11:57
  • 3
    $\begingroup$ The main reason the results differ so much is that the second model demonstrates the first is grossly wrong. The first assumes a binomial response with constant probability within each group, whereas the large standard deviation of the random component (1.8) in the second model shows the data are inconsistent with that assumption of constant probability. Histograms of the proportions correct by group ought to reveal that clearly. A few well chosen plots of the data would help clear up this issue. $\endgroup$ – whuber Nov 1 '18 at 13:45
  • 2
    $\begingroup$ Interesting thread, though one key piece of information seems to be missing: Have students in each group answered the test multiple times, such that each time you recorded their number of correct answers relative to the number of questions answered? If not, I don't see why you would include a random student effect in your model. $\endgroup$ – Isabella Ghement Nov 1 '18 at 14:13
  • 1
    $\begingroup$ The other issue is that these are naturally occurring groups, so they are likely to differ in many ways (which may explain the differences seen between groups). If you can control for any confounding variables that you may have measured, that might help with the interpretation of your results. $\endgroup$ – Isabella Ghement Nov 1 '18 at 14:15
  • 2
    $\begingroup$ Short comment: lsmeans is deprecated. You should switch to emmeans. The syntax will most likely be nearly identical. $\endgroup$ – COOLSerdash Nov 1 '18 at 20:46

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:

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.


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$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.