R:从引导混合效应模型结果中获取系数&CI

R: obtain coefficients&CI from bootstrapping mixed-effect model results

工作数据如下:

set.seed(1234)
df <- data.frame(y = rnorm(1:30), 
                 fac1 = as.factor(sample(c("A","B","C","D","E"),30, replace = T)),
                 fac2 = as.factor(sample(c("NY","NC","CA"),30,replace = T)),
                 x = rnorm(1:30))

lme模型拟合为:

library(lme4)
mixed <- lmer(y ~ x + (1|fac1) + (1|fac2), data = df)

我使用 bootMer 到 运行 参数自举,我可以成功获得固定和随机效应的系数(截距)和 SE:

mixed_boot_sum <- function(data){s <- sigma(data)
c(beta = getME(data, "fixef"), theta = getME(data, "theta"), sigma = s)}

mixed_boot <- bootMer(mixed, FUN = mixed_boot_sum, nsim = 100, type = "parametric", use.u = FALSE)

我的第一个问题是如何从引导结果中获得两个随机效应的每个单独水平的系数(斜率)mixed_boot

我使用 broom 包中的 augment 函数从 mixed 模型中提取系数(斜率)没问题,见下文:

library(broom)
mixed.coef <- augment(mixed, df)

不过,好像broom 不能处理boot class 对象。我不能直接在 mixed_boot.

上使用上述功能

我还尝试通过添加 mmList 来修改 mixed_boot_sum(我认为这就是我要找的),但 R 抱怨为:

Error in bootMer(mixed, FUN = mixed_boot_sum, nsim = 100, type = "parametric",  : 
  bootMer currently only handles functions that return numeric vectors

此外,是否可以通过指定FUN同时获得固定和随机效果的CI?

现在,为了满足我的需求,我对 FUN 的正确规格感到非常困惑。任何有关我的问题的帮助将不胜感激!

My first question is how to obtain the coefficients(slope) of each individual levels of the two random effects from the bootstrapping results mixed_boot ?

我不确定 "coefficients(slope) of each individual level" 是什么意思。 broom::augment(mixed, df) 给出了 每个观测值 的预测(残差等)。如果你想要每个级别的预测系数,我会尝试

mixed_boot_coefs <- function(fit){
   unlist(coef(fit))
}

原始模型给出

mixed_boot_coefs(mixed)
## fac1.(Intercept)1 fac1.(Intercept)2 fac1.(Intercept)3 fac1.(Intercept)4 
##        -0.4973925        -0.1210432        -0.3260958         0.2645979 
## fac1.(Intercept)5           fac1.x1           fac1.x2           fac1.x3 
##        -0.6288728         0.2187408         0.2187408         0.2187408 
##           fac1.x4           fac1.x5 fac2.(Intercept)1 fac2.(Intercept)2 
##         0.2187408         0.2187408        -0.2617613        -0.2617613 
##  ...

如果你想让生成的对象更清楚地命名,你可以使用:

flatten <- function(cc) setNames(unlist(cc),
                                outer(rownames(cc),colnames(cc),
                                      function(x,y) paste0(y,x)))

mixed_boot_coefs <- function(fit){
   unlist(lapply(coef(fit),flatten))
}

当 运行 到 bootMer/confint/boot::boot.ci 时,这些函数将为每个值提供置信区间(请注意,所有斜率 fac<strong>W</strong>.x<strong>Z</strong> 在各组中是相同的,因为模型仅假设截距随机变化)。换句话说,无论您知道如何从拟合模型中提取信息(条件 modes/BLUPs [ranef]、分组变量每个级别的预测截距和斜率 [coef]、参数估计 [fixefgetME]、随机效应方差 [VarCorr]、特定条件下的预测 [predict] ...) 可用于 bootMerFUN参数,只要你能把它的结构扁平化成一个简单的数值向量.