simulate-mer {lme4}R Documentation

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

## S4 method for signature 'mer'
simulate(object, nsim = 1 , seed = NULL, ...)
## S4 method for signature 'mer'
refit(object, newresp, ...)

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, the value is saved as the 'seed' attribute of the returned value. The default, 'NULL', will not change the random generator state, and return '.Random.seed' as the 'seed' attribute

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

Methods

refit

signature(object = "mer", newresp = "numeric"): Update the response vector only and refit the model. See refit.

simulate

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.

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

[Package lme4 version 0.999375-42 Index]