lme4 (version 0.999999-2)

simulate-mer: Simulate based on mer fits

Description

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

Usage

## S3 method for class 'mer':
simulate(object, nsim = 1 , seed = NULL, \dots)
## S3 method for class 'mer':
refit(object, newresp, \dots)

Arguments

object

An object of a suitable class - usually an "mer" object.

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, seed is used in a call to set.seed before simulating the response vectors. If set

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.

Value

• 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).

Examples

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)))
## 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))
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)

cbpp_PB <- replicate(500,pboot(gm0,gm1))

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(...)
       })