Confidence Interval of lmer model producing NA

Question

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

Ben Bolker Ben Bolker · Accepted Answer · 2015-08-09T18:52:44

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

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

Confidence Interval of lmer model producing NA

1 Answers