pbkrtest (version 0.4-7)

pb-refdist: 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)
# S3 method for lm
PBrefdist(largeModel, smallModel, nsim = 1000, seed = NULL,
  cl = NULL, details = 0)
# S3 method for merMod
PBrefdist(largeModel, smallModel, nsim = 1000, seed = NULL,
  cl = NULL, details = 0)
# S3 method for mer
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.

Value

A numeric vector

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.

References

Ulrich Halekoh, S<U+00F8>ren H<U+00F8>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/

See Also

PBmodcomp, KRmodcomp

Examples

# NOT RUN {
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: Many more 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)
# }
# NOT RUN {
## Above, 4 cpu's are used and 5 simulations are made on each cpu.

# }