如何使用 lmeresampler 从嵌套随机效应模型中获取系数的 p 值

Question

我使用lme4 包的lmer 命令估计了一个具有嵌套随机效应结构（参与者在不同组中）的混合效应模型。
mixed.model <- lmer(ln.v ~ treatment*level+age+income+(1 | group/participant),data=data)
然后，由于嵌套结构，我从 lmeresampler 包中引导了 bootstrap 命令。我使用了半参数引导程序。
boot.mixed.model <- bootstrap(model = mixed.model, type = "cgr", fn = extractor, B = 10000, resample=c(data$group,data$participant))
我可以通过boot.ci（包boot）获得自举置信区间，但此外我想报告系数的p值。自举模型boot.mixed.model 的输出仅提供偏差和标准误差：

Bootstrap Statistics :
         original        bias    std. error
t1*   0.658442415 -7.060056e-02  2.34685668
t2*  -0.452128438 -2.755208e-03  0.17041300
…

根据这些值计算 p 值的最佳方法是什么？

Answer 1

我不知道名为 lmeresampler 的软件包，由于兼容性问题（Cran 检查失败），它似乎已从 cran 中删除。

此外，问题不包含数据且未定义 extractor，因此示例不是 reproducible. 但是输出与使用 lme4 中的 bootMer 函数得到的输出相同，因此产生以及使用内置函数的示例。

这基本上遵循help(bootMer) 页面中的示例，但针对特定问题进行了扩展。如果lmeresampler包返回的对象相似，就会包含使用的对象。

可重现的例子

library(lme4)
data(Dyestuff, package = "lme4")
fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)

现在bootMer 函数只需要一个函数，它会输出一个有趣的参数向量。

StatFun <- function(merMod){
    pars <- getME(merMod, c("fixef", "theta", "sigma"))
    c(beta = pars$fixef, theta = unname(pars$theta * pars$sigma), sigma = pars$sigma) ### <<== Error corrected
}

我们可以使用bootMer 执行我们的引导，它还包含type 中的参数选项（我建议阅读help(bootMer) 页面中的详细信息以获取更多信息）

boo01 <- bootMer(fm01ML, StatFun, nsim = 100, seed = 101)

现在，对于更精确的 p 值，我建议 p 值更接近 1000，但由于时间原因，这可能并非在所有情况下都可行。

不管输出是否存储在矩阵 t 中，我们都可以使用它来执行简单的 Kolmogorov-supremum 测试：

H0 <- c(0, 0, 0)
Test <- sweep(abs(boo01$t), 2, H0, "-") <= H0 ###<<=== Error corrected
pVals <- colSums(Test)/nrow(Test)
print(pVals)
#output#
beta.(Intercept)            theta            sigma 
            0.00             0.12             0.00