# 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
/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

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

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.
• @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