simulate-mer {lme4} | R Documentation |
These generic functions (1) generate simulations based on the estimated fitted models (conditional on the estimated values of both the random and fixed effects) and (2) efficiently refit a specified model based on a new vector, data frame, or (for binomial models) matrix of responses
## S4 method for signature 'mer' simulate(object, nsim = 1 , seed = NULL, ...) ## S4 method for signature 'mer' refit(object, newresp, ...)
object |
An object of a suitable class - usually an
|
nsim |
(integer) number of simulations to generate |
seed |
(integer or NULL) an object specifying if and how the
random number generator
should be initialized ('seeded').
If an integer, |
newresp |
a numeric vector or single-column data frame (for most models) or a matrix or two-column data frame (for binomial models where the response was specified as a matrix) containing new values for the response variable |
... |
Some methods may take additional, optional arguments. |
Data frame representing one or more simulations from the fitted model, conditional on both the fixed and random effect estimates. In cases of univariate responses (i.e. all cases except binomial models where the response is specified as a two-column matrix), each column of the data frame represents an independent simulation.
For binomial models where the response is specified as
a two-column matrix, the value is a non-standard data frame
where each column of the data frame (as accessed via x[,i]
or x[[i]]
) is actually a two-column matrix of
responses (successes and failures).
signature(object = "mer", newresp = "numeric")
:
Update the response vector only and refit the model. See
refit
.
signature(object = "mer")
: simulate
nsim
(defaults to 1) responses from the theoretical
distribution corresponding to the fitted model. The refit
method is particularly useful in combination with this method.
See also simulate
.
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy) ## generic parametric bootstrapping function; return a single simulated deviance ## difference between full (m1) and reduced (m0) models under the ## null hypothesis that the reduced model is the true model pboot <- function(m0,m1) { s <- simulate(m0) L0 <- logLik(refit(m0,s)) L1 <- logLik(refit(m1,s)) 2*(L1-L0) } obsdev <- c(2*(logLik(fm1)-logLik(fm2))) ## Not run: ## parametric bootstrap test of significance of correlation between ## random effects of `(Intercept)` and Days ## Timing: approx. 70 secs on a 2.66 GHz Intel Core Duo laptop set.seed(1001) sleepstudy_PB <- replicate(500,pboot(fm2,fm1)) ## End(Not run) library(lattice) ## null value for correlation parameter is *not* on the boundary ## of its feasible space, so we would expect chisq(1) qqmath(sleepstudy_PB,distribution=function(p) qchisq(p,df=1), type="l", prepanel = prepanel.qqmathline, panel = function(x, ...) { panel.qqmathline(x, ...) panel.qqmath(x, ...) }) ## classical test pchisq(obsdev,df=1,lower.tail=FALSE) ## parametric bootstrap-based test mean(sleepstudy_PB>obsdev) ## pretty close in this case ## cbpp data ## PB test of significance of main effect of period gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), family = binomial, data = cbpp) gm0 <- update(gm1, . ~. -period) ## Not run: cbpp_PB <- replicate(500,pboot(gm0,gm1)) ## End(Not run) obsdev <- c(2*(logLik(gm1)-logLik(gm0))) qqmath(cbpp_PB,distribution=function(p) qchisq(p,df=3), type="l", prepanel = prepanel.qqmathline, panel = function(x, ...) { panel.qqmathline(x, ...) panel.qqmath(x, ...) }) ## classical test pchisq(obsdev,df=3,lower.tail=FALSE) ## parametric bootstrap-based test mean(cbpp_PB>obsdev) ## alternative plot nsim <- length(cbpp_PB) pdat <- data.frame(PB=(1:nsim)/(nsim+1), nominal=pchisq(sort(cbpp_PB,decreasing=TRUE),3,lower.tail=FALSE)) xyplot(nominal~PB,data=pdat,type="l",grid=TRUE, scales=list(x=list(log=TRUE),y=list(log=TRUE)), panel=function(...) { panel.abline(a=0,b=1,col="darkgray") panel.xyplot(...) })