[Lme4-commits] r1631 - in pkg/lme4Eigen: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 29 14:51:51 CET 2012
Author: mmaechler
Date: 2012-02-29 14:51:51 +0100 (Wed, 29 Feb 2012)
New Revision: 1631
Modified:
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/NEWS
pkg/lme4Eigen/R/bootMer.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/man/bootMer.Rd
Log:
drop old bootMer() and rename new 'bootMer2' to 'bootMer'
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2012-02-29 11:57:05 UTC (rev 1630)
+++ pkg/lme4Eigen/NAMESPACE 2012-02-29 13:51:51 UTC (rev 1631)
@@ -1,5 +1,4 @@
export(bootMer)
-export(bootMer2)
exportClasses(glmerMod)
exportClasses(lmerMod)
exportClasses(lmList)
Modified: pkg/lme4Eigen/NEWS
===================================================================
--- pkg/lme4Eigen/NEWS 2012-02-29 11:57:05 UTC (rev 1630)
+++ pkg/lme4Eigen/NEWS 2012-02-29 13:51:51 UTC (rev 1631)
@@ -6,7 +6,7 @@
**************************************************
- CHANGES IN lme4 VERSION 0.999375-16
+ CHANGES IN lme4 VERSION 0.999375-16 [Jun 23, 2008]
This midsummer release has many, many changes, relative to earlier
versions. Be careful.
Modified: pkg/lme4Eigen/R/bootMer.R
===================================================================
--- pkg/lme4Eigen/R/bootMer.R 2012-02-29 11:57:05 UTC (rev 1630)
+++ pkg/lme4Eigen/R/bootMer.R 2012-02-29 13:51:51 UTC (rev 1631)
@@ -14,22 +14,22 @@
## instead of using rnorm()
##' Model-based (Semi-)Parametric Bootstrap for Mixed Models
-##'
+##'
##' Perform model-based (Semi-)parametric bootstrap for mixed models.
-##'
+##'
##' The semi-parametric variant is not yet implemented, and we only
##' provide a method for \code{\link{lmer}} and \code{\link{glmer}} results.
-##'
+##'
##' The working name for bootMer() was \dQuote{simulestimate()}, as it is an
##' extension of \code{\link{simulate}}, but we rather want to emphasize its
##' potential for valid inference.
-##'
+##'
##' In each of the \code{nsim} simulations --- that is what the
##' \emph{parametric} bootstrap does, both \dQuote{\emph{spherized}} random
##' effects \eqn{u} and the i.i.d errors \eqn{\epsilon} are generated, using
##' \code{\link{rnorm}()} with parameters corresponding to the fitted model
##' \code{x}.
-##'
+##'
##' @param x fitted \code{*lmer()} model, see \code{\link{lmer}},
##' \code{\link{glmer}}, etc.
##' @param FUN a \code{\link{function}(x)}, computating the \emph{statistic} of
@@ -53,7 +53,7 @@
##' \pkg{boot} package's \code{boot()} result.
##' @seealso For inference, including confidence intervals,
##' \code{\link{profile-methods}}.
-##'
+##'
##' \code{\link[boot]{boot}()}, and then \code{\link[boot]{boot.ci}} from
##' package \pkg{boot}.
##' @references Davison, A.C. and Hinkley, D.V. (1997) \emph{Bootstrap Methods
@@ -65,39 +65,38 @@
##' mySumm <- function(.) { s <- sigma(.)
##' c(beta =getME(., "beta"), sigma = s, sig01 = s * getME(., "theta")) }
##' (t0 <- mySumm(fm01ML)) # just three parameters
-##'
-##' \dontrun{%%--- fails for now --- FIXME
-##'
+##'
##' ## 3.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
##' system.time( boo01 <- bootMer(fm01ML, mySumm, nsim = 100) )
-##'
+##'
##' ## to "look" at it
##' if(need.boot <- is.na(match("package:boot", search())))
##' require("boot")## is a recommended package, i.e. *must* be there
##' boo01 # look at estimated bias for sig01 (-9.1, also when nsim = 1000)
-##'
+##'
##' ## ------ Bootstrap-based confidence intervals ------------
-##'
+##'
##' (bCI.1 <- boot.ci(boo01, index=1, type=c("norm", "basic", "perc")))# beta
-##'
+##'
##' ## Sigma - original scale, and interval calculation on log-scale:
##' (bCI.2 <- boot.ci(boo01, index=2, type=c("norm", "basic", "perc")))
##' (bCI.2l <- boot.ci(boo01, index=2, type=c("norm", "basic", "perc"),
##' h = log, hdot = function(.) 1/., hinv = exp))
-##'
+##'
##' (bCI.3 <- boot.ci(boo01, index=3, type=c("norm", "basic", "perc")))# sig01
-##' }
+##'
##' @export
-bootMer2 <- function(x, FUN, nsim = 1, seed = NULL, use.u = FALSE,
+bootMer <- function(x, FUN, nsim = 1, seed = NULL, use.u = FALSE,
type=c("parametric","semiparametric"),
verbose = FALSE, control = list(),
- .progress="none", PBargs=list()) {
- stopifnot((nsim <- as.integer(nsim[1])) > 0)
- if (.progress!="none") { ## progress bar
- pbfun <- get(paste(.progress,"ProgressBar",sep=""))
- setpbfun <- get(paste("set",.simpleCap(.progress),"ProgressBar",sep=""))
- pb <- do.call(pbfun,PBargs)
- }
+ .progress="none", PBargs=list())
+{
+ stopifnot((nsim <- as.integer(nsim[1])) > 0)
+ if (.progress!="none") { ## progress bar
+ pbfun <- get(paste(.progress,"ProgressBar",sep=""))
+ setpbfun <- get(paste("set",.simpleCap(.progress),"ProgressBar",sep=""))
+ pb <- do.call(pbfun,PBargs)
+ }
FUN <- match.fun(FUN)
type <- match.arg(type)
if(!is.null(seed)) set.seed(seed)
@@ -115,24 +114,24 @@
## FIXME: remove prefix when incorporated in package
if (type=="parametric") {
- ss <- simulate(x, nsim=nsim, use.u=use.u)
+ ss <- simulate(x, nsim=nsim, use.u=use.u)
} else {
- if (use.u) {
- ## FIXME: does this work for GLMMs???
- ss <- replicate(nsim,fitted(x)+sample(residuals(x,"response")),
- simplify=FALSE)
- } else {
- stop("not implemented")
- }
+ if (use.u) {
+ ## FIXME: does this work for GLMMs???
+ ss <- replicate(nsim,fitted(x)+sample(residuals(x,"response")),
+ simplify=FALSE)
+ } else {
+ stop("not implemented")
+ }
}
t.star <- matrix(t0, nrow = length(t0), ncol = nsim)
for(i in 1:nsim) {
- if (.progress!="none") { setpbfun(pb,i/nsim) }
- foo <- try(FUN(refit(x,ss[[i]])))
- if(verbose) { cat(sprintf("%5d :",i)); str(foo) }
- t.star[,i] <- if (inherits(foo, "error")) NA else foo
+ if (.progress!="none") { setpbfun(pb,i/nsim) }
+ foo <- try(FUN(refit(x,ss[[i]])))
+ if(verbose) { cat(sprintf("%5d :",i)); str(foo) }
+ t.star[,i] <- if (inherits(foo, "error")) NA else foo
}
- if (.progress!="none") { close(pb) }
+ if (.progress!="none") { close(pb) }
rownames(t.star) <- names(t0)
## boot() ends with the equivalent of
@@ -146,4 +145,4 @@
## these two are dummies
ran.gen = "simulate(<lmerMod>, 1, *)", mle = mle),
class = "boot")
-}## {bootMer}
+} ## {bootMer}
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-02-29 11:57:05 UTC (rev 1630)
+++ pkg/lme4Eigen/R/lmer.R 2012-02-29 13:51:51 UTC (rev 1631)
@@ -611,171 +611,9 @@
## }
## }
-### bootMer() --- <==> (TODO: semi-*)parametric bootstrap
-### -------------------------------------------------------
-## Doc: show how this is equivalent - but faster than
-## boot(*, R = nsim, sim = "parametric", ran.gen = simulate(.,1,.), mle = x)
-## --> return a "boot" object -- so we could use boot.ci() etc
-## TODO: also allow "semi-parametric" model-based bootstrap:
-## resampling the (centered!) residuals (for use.u=TRUE) or for use.u=FALSE,
-## *both* the centered u's + centered residuals
-## instead of using rnorm()
+## bootMer() ---> now in ./bootMer.R
-##' Model-based (Semi-)Parametric Bootstrap for Mixed Models
-##'
-##' Perform model-based (Semi-)parametric bootstrap for mixed models.
-##'
-##' Currently, the semi-parametric variant is not yet implemented, and we only
-##' provide a method for \code{\linkS4class{merMod}} classes, i.e.,
-##' \code{\link{lmer}} results.
-##'
-##' The working name for bootMer() was \dQuote{simulestimate()}, as it is an
-##' extension of \code{\link{simulate}}, but we rather want to emphasize its
-##' potential for valid inference.
-##'
-##' In each of the \code{nsim} simulations --- that is what the
-##' \emph{parametric} bootstrap does, both \dQuote{\emph{spherized}} random
-##' effects \eqn{u} and the i.i.d errors \eqn{\epsilon} are generated, using
-##' \code{\link{rnorm}()} with parameters corresponding to the fitted model
-##' \code{x}.
-##'
-##' @param x fitted \code{*lmer()} model, see \code{\link{lmer}},
-##' \code{\link{glmer}}, etc.
-##' @param FUN a \code{\link{function}(x)}, computating the \emph{statistic} of
-##' interest, which must be a numeric vector, possibly named.
-##' @param nsim number of simulations, positive integer; the bootstrap \eqn{B}
-##' (or \eqn{R}).
-##' @param seed optional argument to \code{\link{set.seed}}.
-##' @param use.u logical, indicating, if the spherized random effects should be
-##' simulated / bootstrapped as well. If \code{FALSE}, they are not changed,
-##' and all inference is conditional on these.
-##' @param verbose logical indicating if progress should print output
-##' @param control an optional \code{\link{list}}, to be passed to the minimizer
-##' (of the log-likelihood, or RE likelihood), which is currently set to
-##' \code{\link[minqa]{bobyqa}} in package \pkg{minqa}.
-##' @return an object of S3 \code{\link{class}} \code{"boot"}, compatible with
-##' \pkg{boot} package's \code{boot()} result.
-##' @seealso For inference, including confidence intervals,
-##' \code{\link{profile-methods}}.
-##'
-##' \code{\link[boot]{boot}()}, and then \code{\link[boot]{boot.ci}} from
-##' package \pkg{boot}.
-##' @references Davison, A.C. and Hinkley, D.V. (1997) \emph{Bootstrap Methods
-##' and Their Application}. Cambridge University Press.
-##' @keywords models htest
-##' @examples
-##' fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
-##' ## see ?"profile-methods"
-##' mySumm <- function(.) { s <- sigma(.)
-##' c(beta =getME(., "beta"), sigma = s, sig01 = s * getME(., "theta")) }
-##' (t0 <- mySumm(fm01ML)) # just three parameters
-##'
-##' \dontrun{%%--- fails for now --- FIXME
-##'
-##' ## 3.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
-##' system.time( boo01 <- bootMer(fm01ML, mySumm, nsim = 100) )
-##'
-##' ## to "look" at it
-##' if(need.boot <- is.na(match("package:boot", search())))
-##' require("boot")## is a recommended package, i.e. *must* be there
-##' boo01 # look at estimated bias for sig01 (-9.1, also when nsim = 1000)
-##'
-##' ## ------ Bootstrap-based confidence intervals ------------
-##'
-##' (bCI.1 <- boot.ci(boo01, index=1, type=c("norm", "basic", "perc")))# beta
-##'
-##' ## Sigma - original scale, and interval calculation on log-scale:
-##' (bCI.2 <- boot.ci(boo01, index=2, type=c("norm", "basic", "perc")))
-##' (bCI.2l <- boot.ci(boo01, index=2, type=c("norm", "basic", "perc"),
-##' h = log, hdot = function(.) 1/., hinv = exp))
-##'
-##' (bCI.3 <- boot.ci(boo01, index=3, type=c("norm", "basic", "perc")))# sig01
-##' }
-##' @export
-bootMer <- function(x, FUN, nsim = 1, seed = NULL, use.u = FALSE,
- verbose = FALSE, control = list()) {
- stopifnot((nsim <- as.integer(nsim[1])) > 0,
- is(x, "merMod"), is(x at resp, "lmerResp"))
- FUN <- match.fun(FUN)
- if(!is.null(seed)) set.seed(seed)
- else if(!exists(".Random.seed", envir = .GlobalEnv))
- runif(1) # initialize the RNG if necessary
- mc <- match.call()
- t0 <- FUN(x)
- if (!is.numeric(t0))
- stop("bootMer currently only handles functions that return numeric vectors")
-
- ## simplistic approach {original for old-lme4 by DB was much smarter}
-
- n <- nrow(X <- getME(x, "X"))
- if(use.u) {
- u <- getME(x,"u")
- } else {
- U <- getME(x,"Z") %*% getME(x, "Lambda")
- q <- ncol(U)
- }
- ## Zt <- x at Z
-
- beta <- getME(x, "beta")
- X.beta <- as.vector(X %*% beta) # fixed-effect contribution
- sigm.x <- sigma(x)
-
- ## Here, and below ("optimize"/"bobyqa") using the "logic" of lmer() itself:
- ## lmer..Update <- if(is(x, "lmerSp")) lmerSpUpdate else lmerDeUpdate
- ## devfun <- mkdevfun(x)
- ## oneD <- length(x at re@theta) < 2
- theta0 <- getME(x,"theta")
- ## just for the "boot" result -- TODOmaybe drop
- mle <- list(beta = beta, theta = theta0, sigma = sigm.x)
-
- t.star <- matrix(t0, nrow = length(t0), ncol = nsim)
- ## resp <- x at resp
- for(i in 1:nsim) {
- y <- {
- X.beta + sigm.x *
- ((if(use.u) u else as.vector(U %*% rnorm(q))) + rnorm(n))
- ## random effects contribution + Error
- }
- x @ resp <- new(Class=class(resp), REML=resp$REML, y=y, offset=resp$offset, weights=resp$weights)
-
- rho <- new.env(parent=parent.env(environment()))
- rho$pp <- x at pp
- rho$resp <- x @ resp ## FIXME: ???
-
- ## if (oneD) { # use optimize
- ## d0 <- devfun(0)
- ## opt <- optimize(devfun, c(0, 10))
- ## ## -------- <<< arbitrary
- ## ## FIXME ?! if optimal theta > 0, optimize will *not* warn!
- ## if (d0 <= opt$objective) ## prefer theta == 0 when close
- ## devfun(0) # -> theta := 0 and update the rest
- ## } else {
- devfun <- mkdevfun(rho, 0L)
- opt <- bobyqa(theta0, devfun, x at lower, control = control)
- xx <- mkMerMod(environment(devfun), opt, vals$reTrms, x at frame, mc)
- ## xx <- updateMod(x, opt$par, opt$fval)
- ## FIXME: also here, prefer \hat\sigma^2 == 0 (exactly)
- ## }
- foo <- tryCatch(FUN(xx), error = function(e)e)
- if(verbose) { cat(sprintf("%5d :",i)); str(foo) }
- t.star[,i] <- if (inherits(foo, "error")) NA else foo
- }
- rownames(t.star) <- names(t0)
-
- ## boot() ends with the equivalent of
- ## structure(list(t0 = t0, t = t.star, R = R, data = data, seed = seed,
- ## statistic = statistic, sim = sim, call = call,
- ## ran.gen = ran.gen, mle = mle),
- ## class = "boot")
- structure(list(t0 = t0, t = t(t.star), R = nsim, data = x at frame,
- seed = .Random.seed,
- statistic = FUN, sim = "parametric", call = mc,
- ## these two are dummies
- ran.gen = "simulate(<lmerMod>, 1, *)", mle = mle),
- class = "boot")
-}## {bootMer}
-
## Methods for the merMod class
## Anova for merMod objects
Modified: pkg/lme4Eigen/man/bootMer.Rd
===================================================================
--- pkg/lme4Eigen/man/bootMer.Rd 2012-02-29 11:57:05 UTC (rev 1630)
+++ pkg/lme4Eigen/man/bootMer.Rd 2012-02-29 13:51:51 UTC (rev 1631)
@@ -3,7 +3,9 @@
\title{Model-based (Semi-)Parametric Bootstrap for Mixed Models}
\usage{
bootMer(x, FUN, nsim = 1, seed = NULL, use.u = FALSE,
- verbose = FALSE, control = list())
+ type = c("parametric", "semiparametric"),
+ verbose = FALSE, control = list(), .progress = "none",
+ PBargs = list())
}
\arguments{
\item{x}{fitted \code{*lmer()} model, see
@@ -23,13 +25,23 @@
\code{FALSE}, they are not changed, and all inference is
conditional on these.}
+ \item{type}{character string specifying the type of
+ bootstrap, \code{"parametric"} or
+ \code{"semiparametric"}. Default is \code{"parametric"}.
+ Partial matching is allowed.}
+
\item{verbose}{logical indicating if progress should
print output}
\item{control}{an optional \code{\link{list}}, to be
passed to the minimizer (of the log-likelihood, or RE
- likelihood), which is currently set to
- \code{\link[minqa]{bobyqa}} in package \pkg{minqa}.}
+ likelihood).}
+
+ \item{.progress}{character string - type of progress bar
+ to display. Default is \code{"none"}.}
+
+ \item{PBargs}{a list of additional arguments to the
+ progress bar function.}
}
\value{
an object of S3 \code{\link{class}} \code{"boot"},
@@ -41,10 +53,9 @@
models.
}
\details{
- Currently, the semi-parametric variant is not yet
- implemented, and we only provide a method for
- \code{\linkS4class{merMod}} classes, i.e.,
- \code{\link{lmer}} results.
+ The semi-parametric variant is not yet implemented, and
+ we only provide a method for \code{\link{lmer}} and
+ \code{\link{glmer}} results.
The working name for bootMer() was
\dQuote{simulestimate()}, as it is an extension of
@@ -65,8 +76,6 @@
c(beta =getME(., "beta"), sigma = s, sig01 = s * getME(., "theta")) }
(t0 <- mySumm(fm01ML)) # just three parameters
-\dontrun{\%\%--- fails for now --- FIXME
-
## 3.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
system.time( boo01 <- bootMer(fm01ML, mySumm, nsim = 100) )
@@ -85,8 +94,8 @@
h = log, hdot = function(.) 1/., hinv = exp))
(bCI.3 <- boot.ci(boo01, index=3, type=c("norm", "basic", "perc")))# sig01
+
}
-}
\references{
Davison, A.C. and Hinkley, D.V. (1997) \emph{Bootstrap
Methods and Their Application}. Cambridge University
@@ -101,4 +110,3 @@
}
\keyword{htest}
\keyword{models}
-
More information about the Lme4-commits
mailing list