R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Calculate reference distribution using parametric bootstrap

PBrefdist

R Documentation

Calculate reference distribution using parametric bootstrap

Description

Calculate reference distribution of likelihood ratio statistic in mixed effects models using parametric bootstrap

Usage

PBrefdist(largeModel, smallModel, nsim = 1000, seed=NULL, cl = NULL, details = 0)

Arguments

`largeModel`	A linear mixed effects model as fitted with the `lmer()` function in the lme4 package. This model muse be larger than `smallModel` (see below).
`smallModel`	A linear mixed effects model as fitted with the `lmer()` function in the lme4 package. This model muse be smaller than `largeModel` (see above).
`nsim`	The number of simulations to form the reference distribution.
`seed`	Seed for the random number generation.
`cl`	A vector identifying a cluster; used for calculating the reference distribution using several cores. See examples below.
`details`	The amount of output produced. Mainly relevant for debugging purposes.

Details

The model object must be fitted with maximum likelihood (i.e. with REML=FALSE). If the object is fitted with restricted maximum likelihood (i.e. with REML=TRUE) then the model is refitted with REML=FALSE before the p-values are calculated. Put differently, the user needs not worry about this issue.

Value

A numeric vector

Author(s)

Soren Hojsgaard sorenh@math.aau.dk

References

Ulrich Halekoh, S<c3><b8>ren H<c3><b8>jsgaard (2014)., A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models - The R Package pbkrtest., Journal of Statistical Software, 58(10), 1-30., http://www.jstatsoft.org/v59/i09/

Examples

data(beets)
head(beets)
beet0<-lmer(sugpct~block+sow+harvest+(1|block:harvest), data=beets, REML=FALSE)
beet_no.harv <- update(beet0, .~.-harvest)
rr <- PBrefdist(beet0, beet_no.harv, nsim=20)
rr

## Note clearly many more than 10 simulations must be made in practice.

## Computations can be made in parallel using several processors:
## Not run: 
cl <- makeSOCKcluster(rep("localhost", 4))
clusterEvalQ(cl, library(lme4))
clusterSetupSPRNG(cl)
rr <- PBrefdist(beet0, beet_no.harv, nsim=20)
stopCluster(cl)

## End(Not run)
## Above, 4 cpu's are used and 5 simulations are made on each cpu.