CODEWITHSUNDEEP
rlme4glmm
Maverick Meerkat
Maverick Meerkat

Reputation: 6404

GLMM Random Intercept estimators in lme4

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]

Upvotes: 0

Views: 104

Answers (2)

Ben Bolker
Ben Bolker

Reputation: 226097

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

Upvotes: 2

Maverick Meerkat
Maverick Meerkat

Reputation: 6404

Turns out you can get them by multiplying the u parameter with the theta parameter, or by calling getME(.,"b"):

yy = coef(fit_glmm2) # looking only at the intercept
fit_glmm2@u*fit_glmm2@theta + fit_glmm2@beta[1] # or 
# getME(fit_glmm2,"b") + fit_glmm2@beta[1]

Upvotes: 0

Related Questions