lmer模型产生NA的置信区间

问题描述

NA对于lmer模型的置信区间正在发生?我该如何摆脱呢?

NA is occurring for confidence interval of lmer model ? How can I get rid of it ?

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,sqrt(0.5))
     sim_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,control=lmerControl(optCtrl=list(maxfun=20000)))

   ci=confint.merMod(fit,oldName=FALSE,c("sd_(Intercept)|group","sd_X|group","sigma"))

    ci.u0 = as.numeric(ci[1,])
    nc.u0 = ifelse((ci.u0[1]<sqrt(sig2_0) & ci.u0[2]>sqrt(sig2_0)),0,1)

    ci.u1 = as.numeric(ci[2,])
    nc.u1 = ifelse((ci.u1[1]<sqrt(sig2_1) & ci.u1[2]>sqrt(sig2_1)),0,1)

    ci.e = as.numeric(ci[3,])
    nc.e = ifelse((ci.e[1]<sqrt(0.5) & ci.e[2]>sqrt(0.5)),0,1)

    nc = data.frame(nc.u0=nc.u0,nc.u1=nc.u1,nc.e=nc.e)

  }

 fit <- replicate(10,noncoverage(10,5,1,.3,.3,.3,(1/18),0,(1/18)))

 fit
, , 1

  nc.u0 nc.u1 nc.e
 1 0     0     0

 , , 2

  nc.u0 nc.u1 nc.e
 1 0     0     0

 , , 3

 nc.u0 nc.u1 nc.e
1 1     0     0

, , 4

 nc.u0 nc.u1 nc.e
1 NA    0     0

 , , 5

  nc.u0 nc.u1 nc.e
1 0     NA    0

, , 6

 nc.u0 nc.u1 nc.e
1 1     0     0

 , , 7

 nc.u0 nc.u1 nc.e
1 0     0     1

 , , 8

  nc.u0 nc.u1 nc.e
1 0     0     0

 , , 9

 nc.u0 nc.u1 nc.e
1 0     0     0

, , 10

 nc.u0 nc.u1 nc.e
1 0     NA    0

推荐答案

此处的问题(在lme4的开发版本中已解决……)问题是，似然曲线是使用样条拟合来构造的.如果轮廓太平坦，则花键拟合将失败.现在，开发版本尝试在这种情况下替代线性插值.

The problem here (which has been fixed in the development version of lme4 ...) is that likelihood profiles are constructed using spline fits. If the profile is too flat, the spline fit will fail. The development version now tries to substitute linear interpolation in this case.

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   <- MASS::mvrnorm(J,mu=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,sqrt(0.5))
    sim_data <- data.frame(Y=y,X=x,Z=rep(z,each=n_j),
                           group=rep(1:J,each=n_j))
}
set.seed(102)
dat <- simfun(10,5,1,.3,.3,.3,(1/18),0,(1/18))
library("lme4") ## version 1.1-9
fit <- lmer(Y~X+Z+X:Z+(X||group),data=dat)
pp <- profile(fit,"theta_",quiet=TRUE)  ## warnings
cc <- confint(pp)  ## warning
##            2.5 %    97.5 %
## .sig01 0.0000000 0.2880457
## .sig02 0.0000000 0.5427609
## .sigma 0.5937762 0.8802176

您还应该注意，这些置信区间是基于ML而不是REML拟合的 ...

You should also note that these confidence intervals are based on ML, not REML, fits ...

这篇关于lmer模型产生NA的置信区间的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！