[R-sig-ME] confint.merMod() problem.

Zaslavsky, Alan M. zaslavsk at hcp.med.harvard.edu
Tue Apr 29 00:09:37 CEST 2014


I am running into what might be a bug in confint.merMod()  with boot.type="norm".  The error occurs on the line near the end of confint.merMod() that labels the value:          dimnames(citab) <- list(names(bb[["t0"]]), pct)
but appears to arise a few lines up, where the expression that creates citab is hardwired to extract from the structure created when boot.type="perc":
        citab <- t(sapply(bci, function(x) x[["percent"]][4:5]))
More details appear below.  I don't understand the structures enough to propose  a fix other than suggesting that the line above should extract in a way dependent on the value of boot.type.

Package:            lme4
Version:            1.1-6

platform       x86_64-w64-mingw32          
version.string R version 3.1.0 (2014-04-10)
=================================================================================================

Here are the last few lines of confint.merMod():   
     bb <- bootMer(object, bootFun, nsim = nsim, ...)
        bci <- lapply(seq_along(bb$t0), boot.out = bb, boot::boot.ci, 
            type = boot.type, conf = level)
        citab <- t(sapply(bci, function(x) x[["percent"]][4:5]))
        a <- (1 - level)/2
        a <- c(a, 1 - a)
        pct <- format.perc(a, 3)
        dimnames(citab) <- list(names(bb[["t0"]]), pct)
        pnames <- rownames(citab)
        if (missing(parm)) parm <- pnames else if (is.numeric(parm)) parm <- pnames[parm]
        citab[parm, ]
    }, stop("unknown confidence interval method"))

=================================================================================================
This is the model and the call to confint() followed by verbose output and error message:
> xmodsp
Linear mixed model fit by REML ['lmerMod']
Formula: im.pneum ~ 0 + group + (0 + group | aco.id)
   Data: mydata
REML criterion at convergence: 10837.66
Random effects:
 Groups   Name    Std.Dev. Corr
 aco.id   groupFH 0.03173      
          groupFL 0.05405  1.00
 Residual         0.41305      
Number of obs: 10022, groups: aco.id, 152
Fixed Effects:
groupFH  groupFL  
  1.135    1.239  
> confint(xmodsp,parm=1:5,method="boot",boot.type="norm",nsim=5,verbose=T,.progress="win")
Computing bootstrap confidence intervals ...
    1 : Named num [1:6] 0.0208 1 0.0603 0.4118 1.1406 ...
 - attr(*, "names")= chr [1:6] "sd_groupFH|aco.id" "cor_groupFL.groupFH|aco.id" "sd_groupFL|aco.id" "sigma" ...
    2 : Named num [1:6] 0.0477 0.9133 0.0657 0.4121 1.138 ...
 - attr(*, "names")= chr [1:6] "sd_groupFH|aco.id" "cor_groupFL.groupFH|aco.id" "sd_groupFL|aco.id" "sigma" ...
    3 : Named num [1:6] 0.0302 1 0.0598 0.4136 1.1177 ...
 - attr(*, "names")= chr [1:6] "sd_groupFH|aco.id" "cor_groupFL.groupFH|aco.id" "sd_groupFL|aco.id" "sigma" ...
    4 : Named num [1:6] 0.0745 1 0.0512 0.4159 1.1364 ...
 - attr(*, "names")= chr [1:6] "sd_groupFH|aco.id" "cor_groupFL.groupFH|aco.id" "sd_groupFL|aco.id" "sigma" ...
    5 : Named num [1:6] 0.055 1 0.0533 0.4148 1.1247 ...
 - attr(*, "names")= chr [1:6] "sd_groupFH|aco.id" "cor_groupFL.groupFH|aco.id" "sd_groupFL|aco.id" "sigma" ...
Error in dimnames(citab) <- list(names(bb[["t0"]]), pct) : 
  length of 'dimnames' [1] not equal to array extent
=================================================================================================
Here are values of the variables involved just before the offending line:
Browse[1]> bb
Call:
bootMer(x = object, FUN = bootFun, nsim = nsim, verbose = ..1, 
    .progress = "win")
Bootstrap Statistics :
      original       bias    std. error
t1* 0.03172667 -0.003491017 0.020826059
t2* 1.00000000 -0.135380451 0.302719857
t3* 0.05404599 -0.003365170 0.008439177
t4* 0.41304973 -0.002945966 0.003552465
t5* 1.13522212  0.008405088 0.014985916
t6* 1.23883999  0.002633327 0.004492324
Browse[1]> str(bci)
List of 6
 $ :List of 4
  ..$ R     : int 5
  ..$ t0    : Named num 0.0317
  .. ..- attr(*, "names")= chr "sd_groupFH|aco.id"
  ..$ call  : language FUN(boot.out = ..1, conf = ..3, type = ..2, index = 1:6[[1L]])
  ..$ normal: num [1, 1:3] 0.95 -0.0056 0.076
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr "sd_groupFH|aco.id"
  .. .. ..$ : chr [1:3] "conf" "" ""
  ..- attr(*, "class")= chr "bootci"
		[5 more components, omitted for legibility]
Browse[1]> citab
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] NULL NULL NULL NULL NULL NULL
=================================================================================================



More information about the R-sig-mixed-models mailing list