lmer随机效应的置信区间答案

【问题标题】：Confidence interval of random effects with lmerlmer随机效应的置信区间
【发布时间】：2015-07-27 14:34:40
【问题描述】：

我正在使用 lme4 包中的 lmer 来计算方差分量的置信区间。

当我拟合模型时，会出现警告消息：

fit <- lmer(Y~X+Z+X:Z+(X|group),data=sim_data)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

我搜索了很多以了解为什么会发生错误，最终做出there is difference between error and warning in the R world的决定。

我想计算模型参数的置信区间并运行显示错误的代码：

 confint.merMod(fit,oldNames=FALSE)
 Computing profile confidence intervals ...
 Error in if (all(x[iu <- upper.tri(x)] == 0)) t(x[!iu]) else t(x)[!iu] : 
 missing value where TRUE/FALSE needed

还有其他方法可以用 lmer 获得随机效应的 CI 吗？

编辑：

simfun <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
    N <- sum(rep(n_j,J))  
    x <- rnorm(N)         
    z <- rnorm(J)        

    mu <- c(0,0)
    sig <- matrix(c(sig2_0,sig01,sig01,sig2_1),ncol=2)
    u   <- rmvnorm(J,mean=mu,sigma=sig)

    b_0j <- g00 + g01*z + u[,1]
    b_1j <- g10 + g11*z + u[,2]

    y <- rep(b_0j,each=n_j)+rep(b_1j,each=n_j)*x + rnorm(N,0,0.5)
    data <- data.frame(Y=y,X=x,Z=rep(z,each=n_j),group=rep(1:J,each=n_j))
  } 

noncoverage <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
    dat <- simfun(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1)
    fit <- lmer(Y~X+Z+X:Z+(X|group),data=dat)
 }

comb1 = replicate(1000,noncoverage(10,5,1,.3,.3,.3,(1/18),0,(1/18)))
comb26 = replicate(1000,noncoverage(100,50,1,.3,.3,.3,(1/8),0,(1/8)))

【问题讨论】：

如果您收到收敛警告，您可能会在后续分析中遇到麻烦。 confint 是在这种情况下通常不起作用的功能之一。你应该保持警惕，不应该相信这个模型。调查为什么它没有收敛。
我们可能需要找出警告的根本原因（尽管有助于您了解错误/警告的区别）。您是否阅读过?convergence 手册页，尤其是“如果您确实看到收敛警告”下面的内容？）您能否发布有关您的问题的更多信息（即可重现的示例）？

标签： r statistics lme4 mixed-models lmer

【解决方案1】：

这完全取决于您从置信区间中寻找什么，但 arm 包中的函数 sim 提供了一种从 lmer 或 glmer 对象的后部获取重复样本的好方法了解固定项和随机项的系数的可变性。

在merTools 包中，我们编写了一个包装器来简化提取这些值并与它们交互的过程：

library(merTools)
randomSims <- REsim(fit, n.sims = 500)
# and to plot it
plotREsim(REsim(fit, n.sims = 500))

merTools 中还有许多其他工具可用于探索这些内容。如果你想要实际的模拟结果，你最好使用arm::sim。

【讨论】：

刚刚遇到您的 merTools 包 - plotREsim(REsim(fit)) 选项非常棒 - 如果组不同，您是否可以选择将比例从 z 分布转换回概率保持较深的颜色符号。从平均数？
您应该在 SO 上打开一个新问题或在 GitHub 上发布一个问题（听起来我们可以对文档进行改进），然后有人会帮助您获得所需的内容。

【解决方案2】：

lmer 模型以矩阵形式显示结果。您可以访问模型的估计值和标准误来计算置信区间。

由于第一行是模型截距的估计值，因此第二行是显示操纵变量估计值的行。第一列是效果的估计值，第二列是标准误。

因此，将 MyModel 更改为模型名称并将数字四舍五入为 2，round(coef(summary(MyModel))[2,1],2)-round(coef(summary(MyModel))[2,2],2)*2 将为您提供置信区间的下限，而只需更改前面公式中加法的减法即可估计的上限： round(coef(summary(MyModel))[2,1],2)+round(coef(summary(MyModel))[2,2],2)*2

【讨论】：

不确定你们是否还在，你们的答案是否比@Daniel 的答案更通用？如果我想从固定或随机效应中找出CI，那我应该去他的答案吗？

【解决方案3】：

您可以使用parameters-package，它为您提供了计算固定效应 CI 的方法，包括小样本量的 Kenward-Roger 近似，或固定和随机效应的标准误差，或完整的模型摘要。请参阅下面的示例。

library(parameters)
library(lme4)
#> Loading required package: Matrix
model <- lmer(mpg ~ wt + (1 | gear), data = mtcars)

ci(model)
#>     Parameter CI   CI_low   CI_high
#> 1 (Intercept) 95 31.89749 40.482616
#> 2          wt 95 -6.29877 -3.791242

ci_kenward(model)
#>     Parameter CI    CI_low   CI_high
#> 1 (Intercept) 95 31.208589 41.171513
#> 2          wt 95 -6.573142 -3.516869

standard_error(model)
#>     Parameter        SE
#> 1 (Intercept) 2.1901246
#> 2          wt 0.6396873

standard_error(model, effects = "random")
#> $gear
#>   (Intercept)
#> 3   0.6457169
#> 4   0.6994964
#> 5   0.9067223

model_parameters(model, df_method = "satterthwaite")
#> Parameter   | Coefficient |   SE |         95% CI |     t |    df |      p
#> --------------------------------------------------------------------------
#> (Intercept) |       36.19 | 2.19 | [31.90, 40.48] | 16.52 | 13.85 | < .001
#> wt          |       -5.05 | 0.64 | [-6.30, -3.79] | -7.89 | 21.92 | < .001

^{由reprex package (v0.3.0) 于 2020 年 2 月 7 日创建}

【讨论】：

【解决方案4】：

这里的问题是您的随机效应是 (X|group)，我假设它应该是 (1|group)，随机截距模型或将两者都包含为 (1+X|group)。只有 (X|group) 会导致问题，因为您允许改变斜率，但不允许截距改变。

【讨论】：