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

Question

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?

Answer 1

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.