# Why is the random intercept variance so much larger in R than in SPSS in my model and how do I interpret the results?

I am new to Cross Validated so please forgive me if this question has been asked before. However, I did not find any post that answered my question, so here it is:

I am running a 3 level multilevel binary logistic regression (one binary outcome variable and one binary predictor variable) with 839 observations nested in 171 study participants nested in 29 groups. I am using the glmer() function of the lme4 package in R. When I am specifying the empty model and testing it against a “normal” logistic regression without a random intercept for groups and participants, the results clearly tell me that my data is clustered at the participant level and that I do need to use multilevel modeling.

Models:
M0_simple: Outcome ~ 1
M0:        Outcome ~ (1 | Group/Person)

Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
M0_simple        1 975.15 979.87 -486.57   973.15
M0               3 831.07 845.25 -412.54   825.07 148.07      2  < 2.2e-16
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Moreover, when I look at the estimates of the empty multilevel model, the random intercept variance at the participant level (level 2) is very high. And when I calculated the VPC for the participant level, the result of .975 is also extremely high.

Random effects:
Groups   Name               Variance Std.Dev.
Person:Group (Intercept)     127.4    11.29
Group        (Intercept)       0.0     0.00
Number of obs: 839, groups:  Person:Group, 171; Group, 29

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   8.8693     0.8254   10.74 <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


These are the results of the random intercept model once I put in my predictor variable:

M1 <- glmer(Outcome ~ Predictor + (1 | Group/Person), family =
binomial("logit"), data = data_M)
Random effects:
Groups   Name                Variance Std.Dev.
Person:Group   (Intercept)   129.1    11.36
Group          (Intercept)     0.0     0.00
Number of obs: 839, groups:  Person:Group, 171; Group, 29

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                         9.5390     0.9430  10.116   <2e-16 ***
Predictor                          -0.8154     0.4899  -1.664   0.0961 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


How should I interpret these huge random intercept variances at the participant level? I realize that the data is obviously strongly clustered at the participant level. Is this variance so high because the level 1 variance is fixed to 3.29 for multilevel logistic regressions? And does this large variance also affect the fixed effects?

I tried to calculate the predicted probabilities for the random intercept level and ended up with 99.9 %. Moreover, odds-ratios for the intercept of 13891.05 do seem weird. Did I misspecify the model somehow or what might be the issue here? When I run the same model with SPSS 23 it gives out much more reasonable results:

GENLINMIXED
/DATA_STRUCTURE SUBJECTS=Group*Person
/FIELDS TARGET=Outcome TRIALS=NONE OFFSET=NONE
/TARGET_OPTIONS DISTRIBUTION=BINOMIAL LINK=LOGIT
/FIXED  EFFECTS=Predictor USE_INTERCEPT=TRUE
/RANDOM USE_INTERCEPT=TRUE SUBJECTS=Group
COVARIANCE_TYPE=VARIANCE_COMPONENTS
/RANDOM USE_INTERCEPT=TRUE SUBJECTS=Group*Person
COVARIANCE_TYPE=VARIANCE_COMPONENTS
/BUILD_OPTIONS TARGET_CATEGORY_ORDER=DESCENDING
INPUTS_CATEGORY_ORDER=DESCENDING
MAX_ITERATIONS=100 CONFIDENCE_LEVEL=95 DF_METHOD=RESIDUAL COVB=ROBUST
PCONVERGE=0.000001(ABSOLUTE)
SCORING=0 SINGULAR=0.000000000001
/EMMEANS_OPTIONS SCALE=ORIGINAL PADJUST=LSD.

Random effects:
Groups   Name                Variance Std.Error
Person:Group   (Intercept)    5.052     .831
Group          (Intercept)    0.0       …
Number of obs: 839, groups:  Person:Group, 171; Group, 29

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                         2.571     0.2413  10.652    .000
Predictor                          -0.445     0.2400  -1.853    .064
---


While the p-values do seem in the same ball park, the results from SPSS make much more sense to interpret. When I calculate the predicted probabilities from the fixed effects of the intercept and the predictor I get 89.3 % and 92.9 % and the odds-ratio for the intercept of 13.076 seem much more likely than 13891.05.

So what it comes down to, I guess, are the following questions:

I have read that if I want to use likelihood ratio tests to determine the significance of a predictor, I have to use a statistical program that uses Maximum Likelihood (ML) and not Restricted Maximum Likelihood (REML). This is why I use glmer() in R.

However, once I have established that a certain model (with the predictor) is correct or not, how do I interpret the results? Can I simply look at the estimates and interpret the odds-ratio and calculate VPCs and predicted probabilities? Or is this susceptible to mistakes, since the variance at level 1 for multilevel logistic regressions is fixed at 3.29 and the “higher” random variances are scaled accordingly?

Am I even allowed to calculate predicted probabilities from the random intercept model and why are the results of SPSS and R so different?

• I suspect you have mostly zeroes or mostly ones and also that many participants scored all zero or all one but that is only a guess. Sep 13, 2018 at 17:09
• Thanks for the comment @mdewey! Yes, I already checked the data and out of the 171 persons in the study, 141 had only either zeros or ones. So I guess there just is very little variation at the observation level, which rescales the person-level variance into this obscure high number, right? For another variable I investigated, only 93 out of the 141 persons had always either zeros or ones and here the variances seemed much more reasonable. Sep 13, 2018 at 19:13

## 1 Answer

Differences may potentially stem from the fact that SPSS is using the Penalized Quasi Likelihood (PQL) method to fit generalized linear mixed effects models. However, it is known in the literature that the PQL is not a good algorithm, especially for Bernoulli data and Poisson data with small counts. The glmer() function uses by default the Laplace approximation, which is better but still can be inferior to the "gold standard" which is the adaptive Gaussian quadrature (AGQ). In case you have random intercepts and only a single grouping factor, you can fit the model using AGQ with glmer() by setting a higher value into the nAGQ argument. If you want to include more than random intercepts, e.g., random slopes, you can have a look at the GLMMadaptive package.

• Thanks @Dimitris! This answer makes a lot of sense and I have found similar information elsewhere, too. However, in my case I can't use higher numerical integration with the nAGQ argument and will have to stay with the Laplace approximation because I have to estimate random intercepts for persons AND groups, right? Also, just to clarify I am not doing something wrong: I have read that using likelihood ratio tests to test for fixed effects is superior to just simply looking at the p-level for a predictor in the model. Is this correct? And can I still interpret the odds-ratio from the model? Sep 13, 2018 at 19:21
• @Sebastian yes, if you want to include random effects for both persons & groups, you can only do Laplace in R. With regard to the estimated coefficients, you have to be aware of the fact that in GLMMs because of the nonlinear link function used in the specification of the model, they have an interpretation conditional on the random effects. For more on this, check the discussion in this question: stats.stackexchange.com/questions/365907/… Sep 13, 2018 at 19:36
• Hi @Dimitris, thank you so much. You're helping me a lot! :) To make sure I unterstand correctly: The coefficients have to be interpreted by essentially saying "holding co-variables (which I don't have) and the random intercepts for persons and groups constant, the odds ratio for the coefficient is X and therefore the predictor increases the odds by...", right? So would the odds ratio give me the change in odds for the average person in the average group? I guess I'm still a bit unsure how to interpret it exactly. Sep 14, 2018 at 9:53
• Moreover, I tried using your GLMMadaptive package and it gives similar results than the glmer package I've used so far. So I guess the odds ratio from the coefficients from these models would be the change in odds for the same person in the same group, right? And the marginal_coefs() function doing in GLMMadaptive is giving me the change in odds across persons and groups , right? Sep 14, 2018 at 9:56
• @Sebastian yes, the interpretation will be conditional on the person. Most often you're interested in marginal interpretation. That is, what is the odds ratio between the group of persons with predictor value $x$ and the group of persons with predictor value $x + 1$. For example, what is the odds ratio between males and females (i.e., groups of people) not the odds ratio if you changed the sex of a specific person. For a summary of these points, check slide 332 of my course notes: drizopoulos.com/courses/EMC/CE08.pdf Sep 14, 2018 at 10:16