Bootstrapping emmeans 来自多级回归，但是“t.star[r, ] 中的错误 <- res[[r]] : number of items to replace is not..."

Question

我有一个数据集，其中参与者经历了 5 种不同的情况，并在每种情况下测量了他们的行为。这是一个具有类似结构的示例数据集：

id1<-1:10
id<-rep(id1, each=5)
socbeh<-as.integer(rnorm(50, 3, 1))
Index1<-1:5
situation<-rep(Index1, times=10)
df<-data.frame(id, situation, socbeh)

我正在使用 R lme4 运行多级回归，参与者 ID 作为随机效应，情况作为 5 级因子预测变量，如下所示：

library(lme4)
model<-lmer(socbeh ~ (1|id)+factor(situation), data=df)
#and then derive the estimated marginal means via emmeans
library(emmeans)
em<-emmeans(model, ~situation)

我想对估计的边际均值及其置信区间进行非参数引导。但是，我所有的尝试都会产生错误：

Error in t.star[r, ] <- res[[r]] : 
  number of items to replace is not a multiple of replacement length

我首先尝试了以下内容：

fun<-function(data, idx){
     model<-lmer(socbeh ~ (1|id) + factor(situation), 
       data=data[idx,])
    rg<-ref_grid(model, mult.levs = rm_levels)
    em_<-emmeans(rg, ~situation)
}

B<-boot(df, fun, R=1000)

然而，这会产生错误："Error in t.star[r, ] <- res[[r]] : number of items to replace is not a multiple of replacement length"

我尝试从数据中删除所有 NA -> 同样的错误。

我试着仔细阅读这个回复中的建议：：

df2 <- model.matrix(~socbeh + situation + id - 1, data=df) 

然后使用 df2 作为数据再次运行启动，但我仍然得到同样的错误。

我也尝试了一个简单版本的函数：

fun<-function(data, idx){
     model<-lmer(socbeh ~ (1|id) + factor(situation), 
       data=data[idx,])
       em<-emmeans(model, ~situation)
}

但是又出现了同样的错误。我还尝试通过

仅引导回归系数

fun<-function(data, idx) {
coef(lmer(socbeh ~ (1|id)+factor(situation), data=data[idx,]))
}
B<-boot(df2, fun, R=1000)

以前对我有用（使用连续预测变量），但现在我收到此错误“model.frame.default(data = data[idx, ], drop.unused.levels = TRUE, : 可变长度不同（针对“因素（情况）”找到）“

不言而喻，我在编程方面非常缺乏经验，根本不了解发生了什么。谁能帮忙？非常感谢！

Answer 1

emmeans()

所做的部分工作是从模型拟合调用中重建数据。因此，您需要数据作为对象存在于模型安装的环境中，并且在安装模型和运行

emmeans

之间不会更改。也许在函数的开头添加一行就足够了，比如

dat <- data[idx, ]

，然后在模型拟合中使用它。也可以将

data = dat

添加到

ref_grid

调用中。

我想知道的另一件事是函数应该返回什么。现在它正在返回一个

emmGrid

对象，这是一件相当复杂的事情。我猜

boot()

期待一个估计向量。如果是这样，我认为你应该在末尾添加一行

predict(em_)

——这将导致它只返回 EMM。

附录

我试过一个类似的例子。毫无疑问，返回值是一个问题，您应该返回一个估计向量。我发现的另一个问题是，并非所有 bootstrap 样本都具有相关因素的所有水平，在这些情况下，您会得到不同数量的估计值。所以你需要非常小心才能正确排列它们。

这是一个例子，我必须使用

fiber

数据，emmeans

function(data, idx) {
    dat = data[idx, ]
    fit = lm(strength ~ machine + diameter, data = dat)
    em = emmeans(fit, "machine", data = dat)
    rtn = c(A = NA, B = NA, C = NA)
    rtn[em@levels$machine] = predict(em)
    rtn
}
> boot(fiber, fun, R = 100)

ORDINARY NONPARAMETRIC BOOTSTRAP

Call:
boot(data = fiber, statistic = fun, R = 100)

Bootstrap Statistics :
    original      bias    std. error
t1* 40.38241 -0.18833737    1.154697
t2* 41.41922 -0.07372249    1.250100
t3* 38.79836 -0.11139516    1.305454

注意

fun

的最后3行。首先我们将所有的返回值设置为

NA

，对应

machine

的三个层次。 然后我们将我们实际估计的那些元素（如在

em@levels$machine

中找到的）设置为我们得到的估计。然后我们返回

rtn

。这些代码行确保每次返回 3 个值。