1

我正在使用lmerfrom lme4package 来计算方差分量的置信区间。

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

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

4 回答 4

3

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

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

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

有许多其他工具可以在merTools. 但是，如果您想要实际结果的模拟，最好使用arm::sim.

于 2015-08-13T14:24:41.450 回答
1

您可以使用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 包于 2020-02-07 创建(v0.3.0)

于 2020-02-07T14:06:08.983 回答
1

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

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

因此，更改MyModel模型名称并将数字四舍五入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

于 2017-11-07T14:10:25.017 回答
0

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

于 2020-02-07T11:30:31.033 回答