# pb-refdist: Calculate reference distribution using parametric bootstrap In hojsgaard/pbkrtest: Parametric Bootstrap and Kenward Roger Based Methods for Mixed Model Comparison

## Description

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

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10``` ```PBrefdist(largeModel, smallModel, nsim = 1000, seed = NULL, cl = NULL, details = 0) ## S3 method for class 'lm' PBrefdist(largeModel, smallModel, nsim = 1000, seed = NULL, cl = NULL, details = 0) ## S3 method for class 'merMod' 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` Used for controlling parallel computations. See sections 'details' and '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.

The argument 'cl' (originally short for 'cluster') is used for controlling parallel computations. 'cl' can be NULL (default), positive integer or a list of clusters.

Special care must be taken on Windows platforms (described below) but the general picture is this:

The recommended way of controlling cl is to specify the component `pbcl` in options() with e.g. `options("pbcl"=4)`.

If cl is NULL, the function will look at if the pbcl has been set in the options list with `getOption("pbcl")`

If cl=N then N cores will be used in the computations. If cl is NULL then the function will look for

A numeric vector

## Author(s)

Søren Højsgaard [email protected]

## References

Ulrich Halekoh, Søren Hø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/

`PBmodcomp`, `KRmodcomp`
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47``` ```data(beets) head(beets) beet0 <- lmer(sugpct ~ block + sow + harvest + (1|block:harvest), data=beets, REML=FALSE) beet_no.harv <- update(beet0, . ~ . -harvest) rd <- PBrefdist(beet0, beet_no.harv, nsim=20) rd ## Note: Many more simulations must be made in practice. # Computations can be made in parallel using several processors: # 1: On OSs that fork processes (that is, not on windows): # -------------------------------------------------------- if (Sys.info()["sysname"] != "Windows"){ N <- 2 ## Or N <- parallel::detectCores() # N cores used in all calls to function in a session options("mc.cores"=N) rd <- PBrefdist(beet0, beet_no.harv, nsim=20) # N cores used just in one specific call (when cl is set, # options("mc.cores") is ignored): rd <- PBrefdist(beet0, beet_no.harv, nsim=20, cl=N) } # In fact, on Windows, the approach above also work but only when setting the # number of cores to 1 (so there is to parallel computing) # In all calls: # options("mc.cores"=1) # rd <- PBrefdist(beet0, beet_no.harv, nsim=20) # Just once # rd <- PBrefdist(beet0, beet_no.harv, nsim=20, cl=1) # 2. On all platforms (also on Windows) one can do # ------------------------------------------------ library(parallel) N <- 2 ## Or N <- detectCores() clus <- makeCluster(rep("localhost", N)) # In all calls in a session options("pb.cl"=clus) rd <- PBrefdist(beet0, beet_no.harv, nsim=20) # Just once: rd <- PBrefdist(beet0, beet_no.harv, nsim=20, cl=clus) stopCluster(clus) ```