Difference between glmer and stan_glmer


#1

Hi noticed some clear differences between glmer and stan_glmer. The difference is in the estimated realizations of the “random effects” - in glmer I get values that appear to be a weighted average between grand mean and group-specific mean (as suggested by Gelman and Hill, 2007: 253). When I use stan_glmer I get very small (factor of 5 or 10 times smaller) estimated realizations.

Originally, I thought it might be because of the small sample size and that the prior is pushing the estimates towards 0. But I now set up a reproducible example with a large sample and the effect remains.

I would appreciate any suggestions or comments that would help me understand why there is a difference and potentially how to change the stan_glmer command to get estimates that would provide more sizable estimates for the realizations of the “random effects”.

Below is some code that allows to reproduce the issue (takes about an hour to run on my small MacBook and provides no warning messages at the end. All Rhat values are 1.0.

Best,
Lucas

rm(list=ls())
library(lme4)
library(blme)
library(rstanarm)
library(glmer2stan)



options(mc.cores = parallel::detectCores())


# Set up a MC analysis 
# 8008 indidivuals split across 25 level 2 units

n1 <- 8000
n2 <- 25
n1.in.2 <- n1/n2                # Equal number of individuals in level 2 units
err.ind <- rnorm(n1)            # normally distributed error across individuals
err.cont <- rnorm(n2)           # normally distributed error across level 2 units
# equal-sized clusters!
unit <- kronecker(c(1:n2),rep(1,n1.in.2))   # level 2 unit indicator
err.cont <- kronecker(err.cont,rep(1,n1.in.2))  # creating level 2 error vector with length n1
# generating two categorial features (one with four groups and one with 15 groups)
feature1.nr <- sample(c(1:4),n1,replace = TRUE)
feature2.nr <- sample(c(1:15),n1,replace = TRUE)
# Now, drawing their group-specific "effects"
feature1.vec <- rnorm(4,0,.3)
feature2.vec <- rnorm(15,0,.5)
# Assigning each observations its feature value 
feature1.x <- feature1.vec[feature1.nr]
feature2.x <- feature2.vec[feature2.nr]

y.star <- 0.3 + feature1.x + feature2.x + err.cont + err.ind
y <- rep(0,length(y.star))
y[y.star>0] <- 1

data.test <- data.frame(y,feature1.nr,feature2.nr,unit)

# Relying on glmer to estimate mixed model ("random effects" modeled for each feature)
m.glmer <- glmer(y ~ 1 + (1|feature1.nr) + (1|feature2.nr) + (1|unit), family=binomial("probit"), data=data.test)
summary(m.glmer)
ranef(m.glmer)


startL <- Sys.time()
m.stan <- stan_glmer(y ~ 1 + (1|feature1.nr) + (1|feature2.nr) + (1|unit), 
                     data=data.test,
                     family = binomial("probit"),
                     chains = 4,
                     cores = 2,
                     iter = 2000,
                     warmup = 800,
                     adapt_delta = 0.99)#,
#prior_covariance = decov(.1, .1, shape=1, scale=1))
endL <- Sys.time()
endL-startL

summary(m.stan)


# looking at feature 1 "random effect" realizations
ranef(m.glmer)$feature1.nr
ranef(m.stan)$feature1.nr
feature1.vec
# median absolute factor of glmer vs stan_glmer
median(abs(c(ranef(m.glmer)$feature1.nr/ranef(m.stan)$feature1.nr)$`(Intercept)`))

# looking at feature 2 "random effect" realizations
ranef(m.glmer)$feature2.nr
ranef(m.stan)$feature2.nr
feature2.vec
# median absolute factor of glmer vs stan_glmer
median(abs(c(ranef(m.glmer)$feature2.nr/ranef(m.stan)$feature2.nr)$`(Intercept)`))

#2

Which prior are you using for stan glmer? The code has a decov line commented out. Did your model include that line?

Could you include the values that are printed out? Just so we don’t have run some models for a long time to inspect the results…


#3

The code I posted is what I was using. The commented out version is from an older attempt. It does not affect the difference.

Here is the output:

looking at feature 2 “random effect” realizations

ranef(m.glmer)$feature2.nr
(Intercept)
1 -0.13601181
2 0.50738305
3 0.39911842
4 -0.64575670
5 0.23289099
6 -0.34473891
7 -0.85845488
8 0.43339556
9 0.35712270
10 -0.06913988
11 0.33540768
12 -0.05049033
13 -0.04830857
14 0.24739588
15 -0.36226043
ranef(m.stan)$feature2.nr
(Intercept)
1 1.911405e-03
2 2.114783e-04
3 2.057607e-03
4 2.074405e-03
5 8.985433e-04
6 4.627525e-03
7 4.923238e-03
8 -6.434273e-03
9 -9.464151e-04
10 7.018660e-03
11 -4.561079e-04
12 1.161595e-04
13 5.913821e-04
14 7.885671e-06
15 5.113488e-03
feature2.vec
[1] -0.2099732 0.2933298 0.3002041 -0.8878242 0.0293822
[6] -0.5403211 -0.9436911 0.3001267 0.2130320 -0.1806606
[11] 0.1751510 -0.2283298 -0.2106052 0.1150211 -0.5447367

median absolute factor of glmer vs stan_glmer

median(abs(c(ranef(m.glmer)$feature2.nr/ranef(m.stan)$feature2.nr)$(Intercept)))
[1] 193.9721


#4

What is packageVersion("rstanam")? I am getting very similar point estimates.


#5

packageVersion(“rstanarm”)
[1] ‘2.15.2’


#6

Windows?


#7

Mac OS


#8

clang4?


#9

Dear all,

Thank you very much for the help. I re-ran the examples with the updated package and the results are now essentially identical.

Thank you very much - Best,
Lucas