update confint names/method=Wald also returns v-cov params

R/profile.R

@@ -1,5 +1,19 @@
 ## --> ../man/profile-methods.Rd
 
+profnames <- function(object,signames) {
+    useSc <- as.logical(object@devcomp$dims[["useSc"]])
+    ntp <- length(getME(object,"theta"))
+    nn <- if (signames) {
+        sprintf(".sig%02d",seq(ntp))
+    } else {
+        tnames(object,old=FALSE,prefix=c("sd","cor"))
+    }
+    if (useSc) {
+        nn <- c(nn,if (signames) ".sigma" else "sigma")
+    }
+    return(nn)
+}
+
 ##' @importFrom splines backSpline interpSpline periodicSpline
 ##' @importFrom stats profile
 ##' @method profile merMod
@@ -382,6 +396,17 @@ profile.merMod <- function(fitted,
 
 ##' Transform 'which' \in {parnames | integer | "beta_" | "theta_"}
 ##' into integer indices
+##' @param which numeric or character vector
+##' @param nvp number of variance-covariance parameters
+##' @param nptot total number of parameters
+##' @param parnames vector of parameter names
+##' @param verbose print messages?
+##' @examples
+##' fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
+##' tn <- names(getME(fm1,"theta"))
+##' nn <- c(tn,names(fixef(fm1)))
+##' get.which("theta_",length(tn),length(nn),nn, verbose=TRUE)
+##' 
 get.which <- function(which, nvp, nptot, parnames, verbose=FALSE) {
     if (is.null(which))
         seq_len(nptot)
@@ -450,19 +475,10 @@ devfun2 <- function(fm, useSc, signames)
     ## opt <- c(pp$theta*sig, sig)
     if (useSc) {
         opt <- Cv_to_Sv(pp$theta, n=vlist, s=sig)
-        names(opt) <- if (signames) {
-            c(sprintf(".sig%02d", seq(length(opt)-1)), ".sigma")
-        } else {
-            c(tnames(fm,old=FALSE,prefix=c("sd","cor")),"sigma")
-        }
     } else {
         opt <- Cv_to_Sv(pp$theta, n=vlist)
-        names(opt) <- if (signames) {
-            sprintf(".sig%02d", seq_along(opt))
-        } else {
-            tnames(fm,old=FALSE,prefix=c("sd","cor"))
-        }
     }
+    names(opt) <- profnames(fm, signames)
     opt <- c(opt, fixef(fm))
     resp <- fm@resp$copy()
     np <- length(pp$theta)
@@ -736,39 +752,47 @@ confint.merMod <- function(object, parm, level = 0.95,
 {
     method <- match.arg(method)
     boot.type <- match.arg(boot.type)
-    if (!missing(parm) && !is.numeric(parm) && method %in% c("profile","boot"))
-        stop("for method='",method,"', 'parm' must be specified as an integer")
-    if (!quiet && method %in% c("profile","boot")) {
-        mtype <- switch(method, profile="profile", boot="bootstrap")
-        message("Computing ",mtype," confidence intervals ...")
-        flush.console()
+    ## 'parm' corresponds to 'which' in other contexts
+    if (method=="boot" && !is.null(FUN)) {
+        ## custom boot function, don't expand parameter names
+    } else {
+        vn <- profnames(object,oldNames)
+        an <- c(vn,names(fixef(object)))
+        if (missing(parm)) {
+            parm <- seq(length(an))
+        } else {
+            parm <- get.which(parm, nvp=length(vn), nptot=length(an),
+                              parnames=an)
+        }
+        if (!quiet && method %in% c("profile","boot")) {
+            mtype <- switch(method, profile="profile", boot="bootstrap")
+            message("Computing ",mtype," confidence intervals ...")
+            flush.console()
+        }
     }
     switch(method,
 	   "profile" =
            {
-               pp <- if(missing(parm)) {
-                   profile(object, signames=oldNames, ...)
-               } else {
-                   profile(object, which=parm, signames=oldNames, ...)
-               }
+               pp <- profile(object, which=parm, signames=oldNames, ...)
                confint(pp, level=level, zeta=zeta)
            },
 	   "Wald" =
            {
+               a <- (1 - level)/2
+               a <- c(a, 1 - a)
+               ci.vcov <- array(NA,dim = c(length(vn), 2L),
+                                dimnames = list(vn, format.perc(a,3)))
                ## copied with small changes from confint.default
                cf <- fixef(object)   ## coef() -> fixef()
                pnames <- names(cf)
-               if (missing(parm))
-                   parm <- pnames
-               else if (is.numeric(parm))
-                   parm <- pnames[parm]
-               a <- (1 - level)/2
-               a <- c(a, 1 - a)
                ## N.B. can't use sqrt(...)[parm] (diag() loses names)
-               ses <- sqrt(diag(vcov(object)[parm,parm, drop=FALSE]))
-               array(cf[parm] + ses %o% qnorm(a), dim = c(length(parm), 2L),
-                     dimnames = list(parm, format.perc(a, 3)))
-               ## only gives confidence intervals on fixed effects ...
+               ses <- sqrt(diag(vcov(object)))
+               ci.fixed <- array(cf + ses %o% qnorm(a),
+                                 dim = c(length(pnames), 2L),
+                                 dimnames = list(pnames, format.perc(a, 3)))
+               vnames <- tnames(object)
+               ci.all <- rbind(ci.vcov,ci.fixed)
+               ci.all[parm,,drop=FALSE]
            },
 	   "boot" =
            {
@@ -778,13 +802,12 @@ confint.merMod <- function(object, parm, level = 0.95,
                    scaleTh <- (isLMM(x) || isNLMM(x))
                    useSc <- as.logical(x@devcomp$dims[["useSc"]])
 		   ## FIXME: still ugly.  Best cleanup via Cv_to_Sv ...
-		   tn <- tnames(x, old=FALSE, prefix=c("sd","cor"))
 		   ss <- if (scaleTh) {	 ## scale variances by sigma and include it
-		       setNames(Cv_to_Sv(th,n=nvec,s=sigma(x)), c(tn,"sigma"))
+		       setNames(Cv_to_Sv(th,n=nvec,s=sigma(x)), vn)
 		   } else if (useSc) { ## don't scale variances but do include sigma
-		       setNames(c(Cv_to_Sv(th,n=nvec),sigma(x)), c(tn,"sigma"))
+		       setNames(c(Cv_to_Sv(th,n=nvec),sigma(x)), vn)
 		   } else {  ## no scaling, no sigma
-		       setNames(Cv_to_Sv(th,n=nvec), tn)
+		       setNames(Cv_to_Sv(th,n=nvec), vn)
 		   }
                    c(ss, fixef(x))
                }
@@ -800,11 +823,10 @@ confint.merMod <- function(object, parm, level = 0.95,
                a <- (1 - level)/2
                a <- c(a, 1 - a)
                dimnames(citab) <- list(names(bb[["t0"]]), format.perc(a, 3))
-               pnames <- rownames(citab)
-               if (missing(parm))
-                   parm <- pnames
-               else if (is.numeric(parm))
-                   parm <- pnames[parm]
+               if (missing(parm)) {
+                   ## only happens if we have custom boot method
+                   parm <- rownames(citab)
+               }
                citab[parm, , drop=FALSE]
            },
            ## should never get here ...

inst/NEWS.Rd

@@ -45,6 +45,13 @@
       \item \code{lmList} now returns objects of class \code{lmList4},
       to avoid overwriting \code{lmList} methods from the recommended
       \code{nlme} package
+      \item names of random effects parameters in \code{confint}
+      changed (modified for consistency across methods);
+      \code{oldNames=TRUE} (default) gives \code{".sig01"}-style names,
+      \code{oldNames=FALSE} gives \code{"sd_(Intercept)|Subject"}-style names
+      \item \code{confint(.,method="Wald")} result now contains rows
+      for random effects parameters (values set to \code{NA})
+      as well as for fixed-effect parameters
     }
   }
   \subsection{BUG FIXES}{

inst/tests/test-methods.R

@@ -8,6 +8,7 @@ if (getRversion() > "3.0.0") {
 
 fm0 <- fit_sleepstudy_0
 fm1 <- fit_sleepstudy_1
+fm2 <- fit_sleepstudy_2
 gm1 <- fit_cbpp_1
 gm2 <- fit_cbpp_2
 ## More objects to use in all contexts :
@@ -17,6 +18,11 @@ dNA <- data.frame(
     xcov= runif(100),
     resp= rnorm(100))
 dNA[sample(1:100, 10), "xcov"] <- NA
+
+CI.boot <- function(fm, nsim=10, seed=101, ...)
+    suppressWarnings(confint(fm, method="boot", nsim=nsim,
+                             quiet=TRUE, seed=seed, ...))
+
 ##
 rSimple <- function(rep = 2, m.u = 2, kinds = c('fun', 'boring', 'meh')) {
     stopifnot(is.numeric(rep), rep >= 1,
@@ -118,29 +124,30 @@ test_that("lmer", {
 context("bootMer confint()")
 set.seed(47)
 test_that("bootMer", {
-    CI.boot <- function(fm, nsim=10, ...)
-        suppressWarnings(confint(fm, method="boot", nsim=nsim, quiet=TRUE, ...))
     ## testing bug-fix for ordering of sd/cor components in sd/cor matrix with >2 rows
     m1 <- lmer(strength~1+(cask|batch),Pastes)
     ci <- CI.boot(m1)
     corvals <- ci[grep("^cor_",rownames(ci)),]
     expect_true(all(abs(corvals) <= 1))
     ## test bootMer with GLMM, multiple RE
-    ci1 <- CI.boot(gm2, nsim=3)
-    ci2 <- CI.boot(gm2, nsim=5, parm=4:6)
-    expect_true(nrow(ci2) == 3)
-    expect_equal(ci1[4:6,], ci2, tolerance = 0.4)# 0.361
+    ci1 <- CI.boot(gm2, nsim=5)
+    ci2 <- CI.boot(gm2, nsim=5, parm=3:6)
+    ci2B <- CI.boot(gm2, nsim=5, parm="beta_")
+    ## previously tested with nsim=5 vs nsim=3
+    expect_true(nrow(ci2) == 4)
+    expect_equal(ci2,ci2B)
+    expect_equal(ci1[3:6,], ci2) ## , tolerance = 0.4)# 0.361
     ## bootMer with NA values
     PastesNA <- Pastes
-    PastesNA$Sample[1:3] <- NA
+    PastesNA$cask[1:3] <- NA
+    ## previously set 'Sample' (sic) -- no effect!
     m2 <- update(m1, data=PastesNA)
-    ci3 <- CI.boot(m2, seed=101)
-    expect_equal(ci, ci3, tol = 0.06)# 0.0425
-    m3 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+    ci3 <- CI.boot(m2)
+    expect_equal(ci, ci3, tol=0.2)
     sleepstudyNA <- sleepstudy
     sleepstudyNA$Days[1:3] <- NA
-    m4 <- update(m3, data = sleepstudyNA)
-    expect_true(nrow(ci4 <- CI.boot(m4, seed = 101)) == 6) # could check more
+    m4 <- update(fm2, data = sleepstudyNA)
+    expect_true(nrow(ci4 <- CI.boot(m4)) == 6) # could check more
     ##
     ## semipar bootstrapping
     m5 <- lmer(Yield ~ 1|Batch, Dyestuff)
@@ -156,17 +163,42 @@ test_that("bootMer", {
                                  quiet=TRUE)))),
                  c(243.7551,256.9104),tol=1e-3)
 })
-
-context("confint")
+context("confint_other")
 test_that("confint", {
     load(system.file("testdata", "gotway_hessianfly.rda", package = "lme4"))
+    ## generated via:
     ## gotway_hessianfly_fit <- glmer(cbind(y, n-y) ~ gen + (1|block),
     ##              data=gotway.hessianfly, family=binomial,
     ##              control=glmerControl(check.nlev.gtreq.5="ignore"))
     ## gotway_hessianfly_prof <- profile(gotway_hessianfly_fit,which=1)
     ## save(list=ls(pattern="gotway"),file="gotway_hessianfly.rda")
     expect_equal(confint(gotway_hessianfly_prof)[1,1],0)
     ## FIXME: should add tests for {-1,1} bounds on correlations as well
+    expect_equal(c(confint(fm1,method="Wald",parm="beta_")),
+                 c(232.301892,8.891041,270.508318,12.043531),
+                 tol=1e-5)
+    ## Wald gives NA for theta values
+    expect_true(all(is.na(confint(fm1,method="Wald",parm="theta_"))))
+
+    ## check names
+    ci1.p <- suppressWarnings(confint(fm1,quiet=TRUE))
+    ci1.w <- confint(fm1,method="Wald")
+    ci1.b <- CI.boot(fm1, nsim=2)
+    expect_equal(dimnames(ci1.p),
+                 list(c(".sig01", ".sigma", "(Intercept)", "Days"),
+                      c("2.5 %", "97.5 %")))
+    expect_equal(dimnames(ci1.p),dimnames(ci1.w))
+    expect_equal(dimnames(ci1.p),dimnames(ci1.b))
+    ci1.p.n <- suppressWarnings(confint(fm1,quiet=TRUE,oldNames=FALSE))
+    ci1.w.n <- confint(fm1,method="Wald", oldNames=FALSE)
+    ci1.b.n <- CI.boot(fm1, nsim=2, oldNames=FALSE)
+    expect_equal(dimnames(ci1.p.n),
+                         list(c("sd_(Intercept)|Subject", "sigma", "(Intercept)", "Days"),
+                              c("2.5 %", "97.5 %")))
+    expect_equal(dimnames(ci1.p.n),dimnames(ci1.w.n))
+    expect_equal(dimnames(ci1.p.n),dimnames(ci1.b.n))
+
+
 })
 
 context("refit")

man/confint.merMod.Rd

@@ -15,10 +15,12 @@
 }
 \arguments{
   \item{object}{a fitted [ng]lmer model}
-  \item{parm}{parameters of which intervals are sought. Specified by
-    integer position, or for \code{method = "profile"}, e.g., via a
-    \code{\link{character}} vector, see Note below, and
-    \code{\link[=profile.merMod]{profile}}.}
+  \item{parm}{parameters for which intervals are sought. Specified by an
+    integer vector of positions, \code{\link{character}} vector of
+    parameter names, or (unless doing parametric bootstrapping with a
+    user-specified bootstrap function) \code{"theta_"} or \code{"beta_"}
+    to specify variance-covariance or fixed effects parameters only: see the
+    \code{which} parameter of \code{\link[=profile.merMod]{profile}}.}
   \item{level}{confidence level \eqn{< 1}, typically above 0.90.}
   \item{method}{a \code{\link{character}} string determining the method
     for computing the confidence intervals.}
@@ -34,9 +36,9 @@
     are unavailable because they require additional components to be
     calculated.)}
   \item{quiet}{(logical) suppress messages about computationally intensive profiling?}
-  \item{oldNames}{(logical) use old-style names for
-    \code{method="profile"}? (See \code{signames} argument to
-    \code{\link{profile}}).}
+  \item{oldNames}{(logical) use old-style names for variance-covariance
+    parameters, e.g. \code{".sig01"}, rather than newer (more informative) names such as
+    \code{"sd_(Intercept)|Subject"}? (See \code{signames} argument to \code{\link{profile}}).}
   \item{\dots}{additional parameters to be passed to
     \code{\link{profile.merMod}} or \code{\link{bootMer}}, respectively.}
 }
@@ -64,7 +66,9 @@
       based on the likelihood ratio test;}
     \item{\code{"Wald"}:}{approximating
       the confidence intervals (of fixed-effect parameters
-      only) based on the estimated local curvature of the
+      only; all variance-covariance parameters
+      CIs will be returned as \code{NA})
+      based on the estimated local curvature of the
       likelihood surface;}
     \item{\code{"boot"}:}{performing parametric
       bootstrapping with confidence intervals computed from the