lme4中的GLMM随机截尾估计

Question

如何从lme4结果对象获得随机拦截效果估计器？

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

我以为它们可能是u的，但它们之间有一些细微的区别：

yy = coef(fit_glmm2) # looking only at the intercept
fit_glmm2@u + fit_glmm2@beta[1]

Answer 1

如果您想要随机效果，ranef()是获得它们的最佳方法：

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

正如vignette("lmer", package = "lme4")所描述的，@u值是球形随机效应，即它们是iid N(0,1)，需要变换才能得到公式X %*% beta + Z %*% b中使用的随机效应b。在这种情况下(纯截距RE)，theta对应于随机效应的标准差.u*theta不能处理更复杂的案子.在这种情况下，您需要getME(fit_glmm2, "Lambda") %*% getME(fit_glmm2, "u")。

getME(., "b")也可以工作，但是对于更复杂的模型，您必须知道如何将b-vector分解为随机拦截、斜率、不同的RE项等等。

Answer 2

结果表明，您可以通过将u参数与theta参数相乘，或者通过调用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]