1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
|
\name{bootMer}
\alias{bootMer}
\title{Model-based (Semi-)Parametric Bootstrap for Mixed Models}
\usage{
bootMer(x, FUN, nsim = 1, seed = NULL, use.u = FALSE, re.form=NA,
type = c("parametric", "semiparametric"),
verbose = FALSE, .progress = "none", PBargs = list(),
parallel = c("no", "multicore", "snow"),
ncpus = getOption("boot.ncpus", 1L), cl = NULL)
}
\arguments{
\item{x}{a fitted \code{merMod} object: see
\code{\link{lmer}}, \code{\link{glmer}}, etc.}
\item{FUN}{a function taking a fitted
\code{merMod} object as input and returning the
\emph{statistic} of interest, which must be a (possibly named)
numeric vector.}
\item{nsim}{number of simulations, positive integer; the
bootstrap \eqn{B} (or \eqn{R}).}
\item{seed}{optional argument to \code{\link{set.seed}}.}
\item{use.u}{logical, indicating whether the spherical
random effects should be simulated / bootstrapped as
well. If \code{TRUE}, they are not changed, and all
inference is conditional on these values. If
\code{FALSE}, new normal deviates are drawn (see
Details).}
\item{re.form}{formula, \code{NA} (equivalent to \code{use.u=FALSE}),
or \code{NULL} (equivalent to \code{use.u=TRUE}):
alternative to \code{use.u} for
specifying which random effects to incorporate.
See \code{\link{simulate.merMod}} for details.}
\item{type}{character string specifying the type of
bootstrap, \code{"parametric"} or
\code{"semiparametric"}; partial matching is allowed.}
\item{verbose}{logical indicating if progress should
print output}
\item{.progress}{character string - type of progress bar
to display. Default is \code{"none"}; the function will
look for a relevant \code{*ProgressBar} function, so
\code{"txt"} will work in general; \code{"tk"} is
available if the \pkg{tcltk} package is loaded; or
\code{"win"} on Windows systems. Progress bars are
disabled (with a message) for parallel operation.}
\item{PBargs}{a list of additional arguments to the
progress bar function (the package authors like
\code{list(style=3)}).}
\item{parallel}{The type of parallel operation to be used (if any).
If missing, the
default is taken from the option \code{"boot.parallel"} (and if that
is not set, \code{"no"}).}
\item{ncpus}{integer: number of processes to be used in parallel operation:
typically one would choose this to be the number of available CPUs.}
\item{cl}{An optional \pkg{parallel} or \pkg{snow} cluster for use if
\code{parallel = "snow"}. If not supplied, a cluster on the
local machine is created for the duration of the \code{boot} call.}
}
\value{
an object of S3 \code{\link{class}} \code{"boot"},
compatible with \CRANpkg{boot} package's
\code{\link[boot]{boot}()} result. (See Details for information on how
to retrieve information about errors during bootstrapping.)
}
\description{
Perform model-based (Semi-)parametric bootstrap for mixed
models.
}
\note{
If you are using \code{parallel="snow"}, you will need to run
\code{clusterEvalQ(cl,library("lme4"))} before calling
\code{bootMer} to make sure that the
\code{lme4} package is loaded on all of the workers; you may
additionally need to use \code{\link[parallel]{clusterExport}}
if you are using a summary function that calls any objects
from the environment.
}
\details{
The semi-parametric variant is only partially implemented, and
we only provide a method for \code{\link{lmer}} and
\code{\link{glmer}} results.
Information about warning and error messages incurred during the
bootstrap returns can be retrieved via the attributes
\describe{
\item{bootFail}{number of failures (errors)}
\item{boot.fail.msgs}{error messages}
\item{boot.all.msgs}{messages, warnings, and error messages}
}
e.g. \code{attr("boot.fail.msgs")} to retrieve error messages
The working name for bootMer() was
\dQuote{simulestimate()}, as it is an extension of \code{simulate}
(see \link{simulate.merMod}), but we want to emphasize its potential
for valid inference.
\itemize{
\item If \code{use.u} is \code{FALSE} and \code{type} is
\code{"parametric"}, each simulation generates new values of both
the \dQuote{\emph{spherical}} random effects \eqn{u} and the
i.i.d. errors \eqn{\epsilon}, using \code{\link{rnorm}()}
with parameters corresponding to the fitted model \code{x}.
\item If \code{use.u} is \code{TRUE} and \code{type=="parametric"},
only the i.i.d. errors (or, for GLMMs, response values drawn from
the appropriate distributions) are resampled, with the values of
\eqn{u} staying fixed at their estimated values.
\item If \code{use.u} is \code{TRUE} and \code{type=="semiparametric"},
the i.i.d. errors are sampled from the distribution of (response)
residuals. (For GLMMs, the resulting
sample will no longer have the same properties as the original
sample, and the method may not make sense; a warning is generated.)
The semiparametric bootstrap is currently an experimental feature,
and therefore may not be stable.
\item The case where \code{use.u} is \code{FALSE} and
\code{type=="semiparametric"} is not implemented; Morris (2002)
suggests that resampling from the estimated values of \eqn{u} is not
good practice.
} %% itemize
} %% details
\references{
\itemize{
\item \insertRef{davison1997bootstrap}{lme4}
\item \insertRef{morris2002blups}{lme4}
}
}
\seealso{
\itemize{
\item \code{\link{confint.merMod}},
for a more specific approach to bootstrap confidence
intervals on parameters.
\item \code{\link{refit}()}, or \code{\link[pbkrtest]{PBmodcomp}()}
from the \CRANpkg{pbkrtest} package, for parametric bootstrap comparison
of models.
\item \code{\link[boot]{boot}()}, and then
\code{\link[boot]{boot.ci}}, from the \pkg{boot} package.
\item \code{\link{profile-methods}}, for likelihood-based inference,
including confidence intervals.
\item \code{\link{pvalues}},
for more general approaches to inference and p-value computation
in mixed models.
}
}
\examples{
if (interactive()) {
fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
## see ?"profile-methods"
mySumm <- function(.) { s <- sigma(.)
c(beta =getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) }
(t0 <- mySumm(fm01ML)) # just three parameters
## alternatively:
mySumm2 <- function(.) {
c(beta=fixef(.),sigma=sigma(.), sig01=sqrt(unlist(VarCorr(.))))
}
set.seed(101)
## 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 (requireNamespace("boot")) {
boo01
## note large estimated bias for sig01
## (~30\% low, decreases _slightly_ for nsim = 1000)
## extract the bootstrapped values as a data frame ...
head(as.data.frame(boo01))
## ------ Bootstrap-based confidence intervals ------------
## warnings about "Some ... intervals may be unstable" go away
## for larger bootstrap samples, e.g. nsim=500
## intercept
(bCI.1 <- boot::boot.ci(boo01, index=1, type=c("norm", "basic", "perc")))# beta
## Residual standard deviation - original scale:
(bCI.2 <- boot::boot.ci(boo01, index=2, type=c("norm", "basic", "perc")))
## Residual SD - transform to log scale:
(bCI.2L <- boot::boot.ci(boo01, index=2, type=c("norm", "basic", "perc"),
h = log, hdot = function(.) 1/., hinv = exp))
## Among-batch variance:
(bCI.3 <- boot::boot.ci(boo01, index=3, type=c("norm", "basic", "perc"))) # sig01
confint(boo01)
confint(boo01,type="norm")
confint(boo01,type="basic")
## Graphical examination:
plot(boo01,index=3)
## Check stored values from a longer (1000-replicate) run:
(load(system.file("testdata","boo01L.RData", package="lme4")))# "boo01L"
plot(boo01L, index=3)
mean(boo01L$t[,"sig01"]==0) ## note point mass at zero!
} %% if boot package available
} %% interactive
}
\keyword{htest}
\keyword{models}
|
