How do you get the random intercept effects estimators from a lme4 result object?
set.seed(247)
# Create Data
n=1000
x = runif(n)
id = rep(NA,n)
for (i in 1:10) {
id_s = (i-1)*100+1
id_e = i*100
id[id_s:id_e] = i
}
effects = rnorm(10)
lp = -0.5+0.5*x + effects[id]
probs = exp(lp)/(1+exp(lp))
Y2 = rbinom(n, 1, probs)
library(lme4)
fit_glmm2 = glmer(Y2 ~ x + (1|id), family = "binomial",control = glmerControl(calc.derivs = FALSE))
I thought maybe they are the u's but there's a slight difference between them:
yy = coef(fit_glmm2) # looking only at the intercept
fit_glmm2@u + fit_glmm2@beta[1]
If you want the random effects, ranef() is the best way to get them:
r <- ranef(fit_glmm2)
str(r)
## List of 1
## $ id:'data.frame': 10 obs. of 1 variable:
## ..$ (Intercept): num [1:10] -0.693 0.297 0.54 -0.467 0.755 ...
## ..- attr(*, "postVar")= num [1, 1, 1:10] 0.0463 0.0385 0.0392 0.0436 0.0409 ...
## - attr(*, "class")= chr "ranef.mer"
raw <- unname(unlist(ranef(fit_glmm2)$id))
identical(raw, fit_glmm2@u*fit_glmm2@theta) ## TRUE
As described in vignette("lmer", package = "lme4"), the @u values are the spherical random effects, i.e. they're iid N(0,1) and need to be transformed to get to the random effects b used in the formula X %*% beta + Z %*% b. In this case (an intercept-only RE), theta corresponds to the standard deviation of the random effect. u*theta won't work for more complicated cases ... in this case you need getME(fit_glmm2, "Lambda") %*% getME(fit_glmm2, "u").
getME(., "b") will also work, but again for more complex models you'll have to work out how the b-vector is split into random intercepts, slopes, different RE terms, etc..