pbkrtest

https://hojsgaard.github.io/pbkrtest/
Admin
pbkrtest is an R package that implements parametric bootstrap, Kenward-Roger, and Satterthwaite...
Tokens:14,921
Snippets:107
Trust Score:9.7
Update:2 months ago
Show doc for...
Tokens:
Raw
### Parallel execution examples

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-refdist.html

Examples demonstrating how to run simulations in parallel on different operating systems.

```R
if (FALSE) {
## 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)
}
```

--------------------------------

### Basic usage example

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-refdist.html

A basic example of calculating the reference distribution using the beets dataset.

```R
data(beets)
head(beets)
#>   harvest  block  sow yield sugpct
#> 1   harv1 block1 sow3 128.0   17.1
#> 2   harv1 block1 sow4 118.0   16.9
#> 3   harv1 block1 sow5  95.0   16.6
#> 4   harv1 block1 sow2 131.0   17.0
#> 5   harv1 block1 sow1 136.5   17.0
#> 6   harv2 block2 sow3 136.5   17.0
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, cl=1)
rd
#>  [1] 3.120726e-01 1.530013e+00 7.464388e+00 2.615508e-02 2.557797e-01
#>  [6] 5.717452e-01 1.258622e+00 1.429834e-03 8.870589e-01 2.728823e-02
#> [11] 7.467763e-02 1.621034e-01 2.606911e+00 5.525129e+00 1.010020e-01
#> [16] 7.143766e-01 3.343534e-02 3.257146e+00 4.674691e+00 2.406896e-07
#> attr(,"cl")
#> [1] 1
#> attr(,"ctime")
#> elapsed 
#>   0.624 
#> attr(,"stat")
#>         tobs           df      p.value 
#> 1.291424e+01 1.000000e+00 3.260911e-04 
#> attr(,"samples")
#>        nsim        npos   n.extreme         pPB 
#> 20.00000000 20.00000000  0.00000000  0.04761905 
#> attr(,"class")
#> [1] "refdist"
```

--------------------------------

### Install pbkrtest

Source: https://hojsgaard.github.io/pbkrtest/index.html

Commands to install the package from CRAN or the development version from GitHub.

```R
## Install from CRAN:
install.packages('pbkrtest')
## Install from Github: Use the remotes package:
remotes::install_github("hojsgaard/pbkrtest", build_vignettes = TRUE)
```

--------------------------------

### Setup and Data Loading

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Load the beets dataset and define the number of simulations for the parametric bootstrap.

```R
data(beets, package="pbkrtest")
head(beets)

NSIM <- 50 ## Simulations in parametric bootstrap
```

--------------------------------

### Example Usage of compare_column_space

Source: https://hojsgaard.github.io/pbkrtest/reference/compare_column_space.html

Demonstrates how to compare column spaces between matrices of different dimensions.

```R
A1 <- matrix(c(1,1,1,1,2,3), nrow=3)
A2 <- A1[, 1, drop=FALSE]

compare_column_space(A1, A2)
#> [1] 1
compare_column_space(A2, A1)
#> [1] 0
compare_column_space(A1, A1)
#> [1] -1
```

--------------------------------

### Perform Satterthwaite F-tests

Source: https://hojsgaard.github.io/pbkrtest/reference/sat-modcomp.html

Examples demonstrating the use of SATmodcomp with restriction matrices, model objects, and formulas.

```R
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
L1 <- cbind(0,1)
SATmodcomp(fm1, L1)
```

```R
(fm2 <- lmer(Reaction ~ Days + I(Days^2) + (Days|Subject), sleepstudy))

## Test for no effect of Days. There are three ways of using the function:

## 1) Define 2-df contrast - since L has 2 (linearly independent) rows
## the F-test is on 2 (numerator) df:
L2 <- rbind(c(0, 1, 0), c(0, 0, 1))
SATmodcomp(fm2, L2)

## 2) Use two model objects 
fm3 <- update(fm2, ~. - Days - I(Days^2))
SATmodcomp(fm2, fm3)

## 3) Specify restriction as formula
SATmodcomp(fm2, ~. - Days - I(Days^2))
```

--------------------------------

### Advanced PBmodcomp Usage with Clusters

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Examples of using PBmodcomp with pre-computed reference distributions and parallel processing clusters.

```R
if (FALSE) {
## Vanilla
PBmodcomp(beet0, beet_no.harv, nsim=NSIM, cl=1)

## Simulate reference distribution separately:
refdist <- PBrefdist(beet0, beet_no.harv, nsim=1000)
PBmodcomp(beet0, beet_no.harv, ref=refdist, cl=1)

## Do computations with multiple processors:
## Number of cores:
(nc <- detectCores())
## Create clusters
cl <- makeCluster(rep("localhost", nc))

## Then do:
PBmodcomp(beet0, beet_no.harv, cl=cl)

## Or in two steps:
refdist <- PBrefdist(beet0, beet_no.harv, nsim=NSIM, cl=cl)
PBmodcomp(beet0, beet_no.harv, ref=refdist)

## It is recommended to stop the clusters before quitting R:
stopCluster(cl)
}
```

--------------------------------

### Calculating Denominator Degrees of Freedom Example

Source: https://hojsgaard.github.io/pbkrtest/reference/get_ddf_Lb.html

Demonstrates fitting linear mixed models, comparing them, and calculating the denominator degrees of freedom for a specific contrast.

```R
(fmLarge <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))

# removing Days
(fmSmall <- lmer(Reaction ~ 1 + (Days|Subject), sleepstudy))

anova(fmLarge, fmSmall)

KRmodcomp(fmLarge, fmSmall)  ## 17 denominator df's

get_Lb_ddf(fmLarge, c(0, 1)) ## 17 denominator df's

# Notice: The restriction matrix L corresponding to the test above
# can be found with
L <- model2restriction_matrix(fmLarge, fmSmall)
L
```

--------------------------------

### Visualize Sugar Beet Data

Source: https://hojsgaard.github.io/pbkrtest/index.html

Example using ggplot2 to visualize the relationship between sugar content and harvest time in the beets dataset.

```R
library(pbkrtest)
library(ggplot2)

## Sugar beets: Does suger content depend on harvest time?

beets |> ggplot(aes(x=sow, y=sugpct, group=harvest)) +
    geom_jitter(aes(color=harvest), width=0)
```

--------------------------------

### Perform Kenward-Roger F-test

Source: https://hojsgaard.github.io/pbkrtest/reference/kr-modcomp.html

Example demonstrating the comparison of two lmer models and the use of a restriction matrix for hypothesis testing.

```R
(fmLarge <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + (Days | Subject)
#>    Data: sleepstudy
#> REML criterion at convergence: 1743.628
#> Random effects:
#>  Groups   Name        Std.Dev. Corr
#>  Subject  (Intercept) 24.741       
#>           Days         5.922   0.07
#>  Residual             25.592       
#> Number of obs: 180, groups:  Subject, 18
#> Fixed Effects:
#> (Intercept)         Days  
#>      251.41        10.47  
## removing Days
(fmSmall <- lmer(Reaction ~ 1 + (Days|Subject), sleepstudy))
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ 1 + (Days | Subject)
#>    Data: sleepstudy
#> REML criterion at convergence: 1769.845
#> Random effects:
#>  Groups   Name        Std.Dev. Corr 
#>  Subject  (Intercept) 25.53         
#>           Days        11.93    -0.18
#>  Residual             25.59         
#> Number of obs: 180, groups:  Subject, 18
#> Fixed Effects:
#> (Intercept)  
#>       257.8  
anova(fmLarge, fmSmall)
#> refitting model(s) with ML (instead of REML)
#> Data: sleepstudy
#> Models:
#> fmSmall: Reaction ~ 1 + (Days | Subject)
#> fmLarge: Reaction ~ Days + (Days | Subject)
#>         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
#> fmSmall    5 1785.5 1801.4 -887.74   1775.5                         
#> fmLarge    6 1763.9 1783.1 -875.97   1751.9 23.537  1  1.226e-06 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
KRmodcomp(fmLarge, fmSmall)
#> large : Reaction ~ Days + (Days | Subject)
#> small : Reaction ~ 1 + (Days | Subject)
#>         stat    ndf    ddf F.scaling   p.value    
#> Ftest 45.853  1.000 17.000         1 3.264e-06 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## The same test using a restriction matrix
L <- cbind(0, 1)
KRmodcomp(fmLarge, L)
#> large : Reaction ~ Days + (Days | Subject)
#> L = 
#>      [,1] [,2]
#> [1,]    0    1
#>         stat    ndf    ddf F.scaling   p.value    
#> Ftest 45.853  1.000 17.000         1 3.264e-06 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```

--------------------------------

### Calculate adjusted covariance matrix

Source: https://hojsgaard.github.io/pbkrtest/reference/kr-vcovAdj.html

Example usage of vcovAdj compared against the standard vcov output for an lmer model.

```R
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
class(fm1)
#> [1] "lmerMod"
#> attr(,"package")
#> [1] "lme4"

## Here the adjusted and unadjusted covariance matrices are identical,
## but that is not generally the case:

v1 <- vcov(fm1)
v2 <- vcovAdj(fm1, details=0)
v2 / v1
#> 2 x 2 Matrix of class "dgeMatrix"
#>             (Intercept) Days
#> (Intercept)           1    1
#> Days                  1    1
```

--------------------------------

### Initialize data and linear models

Source: https://hojsgaard.github.io/pbkrtest/articles/a02-coercion.html

Create a sample dataset and define large and small linear models for comparison.

```R
N <- 4
dat <- data.frame(int=rep(1, N), x=1:N, y=rnorm(N))
```

```R
lg <- lm(y ~ x + I(x^2), data=dat)
sm <- lm(y ~ x, data=dat)
lg
```

```R
sm
```

--------------------------------

### Basic Model Comparison with PBmodcomp

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Initial comparison of two models using anova and PBmodcomp.

```R
anova(glm.D93, glm.D93.t, test="Chisq")
PBmodcomp(glm.D93, glm.D93.t, nsim=NSIM, cl=1)
```

--------------------------------

### Load and inspect the shoes dataset

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Load the shoes dataset from the MASS package and display its contents.

```R
data(shoes, package="MASS")
shoes
```

--------------------------------

### Compare with parametric bootstrap variance estimate

Source: https://hojsgaard.github.io/pbkrtest/reference/kr-vcovAdj.html

Demonstrates an alternative variance-covariance estimation approach using parametric bootstrap.

```R
if (FALSE) {
nsim <- 100
sim <- simulate(fm.ml, nsim)
B <- lapply(sim, function(newy) try(fixef(refit(fm.ml, newresp=newy))))
B <- do.call(rbind, B)
v3 <- cov.wt(B)$cov
v2/v1
v3/v1
}
```

--------------------------------

### Convert models and restriction matrices in R

Source: https://hojsgaard.github.io/pbkrtest/reference/model-coerce.html

Demonstrates constructing restriction matrices from linear models and subsequently reconstructing submodels from those matrices using the pbkrtest package.

```R
library(pbkrtest)
data("beets", package = "pbkrtest")
sug <- lm(sugpct ~ block + sow + harvest, data=beets)
sug.h <- update(sug, .~. - harvest)
sug.s <- update(sug, .~. - sow)

## Construct restriction matrices from models
L.h <- model2restriction_matrix(sug, sug.h); L.h
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,]    0    0    0    0    0    0    0   -1
L.s <- model2restriction_matrix(sug, sug.s); L.s
#>      [,1] [,2] [,3]       [,4]       [,5]       [,6]       [,7] [,8]
#> [1,]    0    0    0 -0.9176629  0.2294157  0.2294157  0.2294157    0
#> [2,]    0    0    0  0.0636285  0.8907987 -0.3181424 -0.3181424    0
#> [3,]    0    0    0  0.1048285  0.1048285  0.8386279 -0.5241424    0
#> [4,]    0    0    0 -0.3779645 -0.3779645 -0.3779645 -0.7559289    0

## Construct submodels from restriction matrices
mod.h <- restriction_matrix2model(sug, L.h); mod.h
#> 
#> Call:
#> lm(formula = sugpct ~ .X1 + .X2 + .X3 + .X4 + .X5 + .X6 + .X7 - 
#>     1, data = structure(list(.X1 = c(0, 0, 0, 0, 0, 1, 1, 1, 
#> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 
#> 0), .X2 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 
#> 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1), .X3 = c(0, 0, 0, 1, 0, 
#> 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
#> 0, 0, 1, 0), .X4 = c(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 
#> 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0), .X5 = c(0, 1, 
#> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 
#> 0, 0, 0, 1, 0, 0, 0), .X6 = c(0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 
#> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1), .X7 = c(1, 
#> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
#> 1, 1, 1, 1, 1, 1, 1, 1), sugpct = c(17.1, 16.9, 16.6, 17, 17, 
#> 17, 17, 16.7, 16.4, 16.8, 16.6, 16.9, 17.1, 16.8, 16.9, 17, 16.8, 
#> 16.5, 16.7, 17, 16.8, 16.9, 17, 17, 16.5, 16.7, 16.6, 16.9, 16.9, 
#> 16.4), block = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
#> 2L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
#> 3L, 3L, 3L, 3L, 3L), levels = c("block1", "block2", "block3"), class = "factor"), 
#>     sow = structure(c(3L, 4L, 5L, 2L, 1L, 3L, 2L, 4L, 5L, 1L, 
#>     5L, 2L, 3L, 4L, 1L, 2L, 1L, 5L, 4L, 3L, 4L, 1L, 3L, 2L, 5L, 
#>     1L, 4L, 3L, 2L, 5L), levels = c("sow1", "sow2", "sow3", "sow4", 
#>     "sow5"), class = "factor"), harvest = structure(c(1L, 1L, 
#>     1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
#>     2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), levels = c("harv1", 
#>     "harv2"), class = "factor")), class = "data.frame", row.names = c(NA, 
#> 30L)))
#> 
#> Coefficients:
#>     .X1      .X2      .X3      .X4      .X5      .X6      .X7  
#> -0.0500  -0.0800   0.1167   0.1667  -0.1000  -0.3500  16.8933  
#> 
mod.s <- restriction_matrix2model(sug, L.s); mod.s
#> 
#> Call:
#> lm(formula = sugpct ~ .X1 + .X2 + .X3 + .X4 - 1, data = structure(list(
#>     .X1 = c(-0.481123223875107, -0.481123223875107, -0.481123223875107, 
#>     -0.481123223875107, -0.481123223875107, -1.35446911968096, 
#>     -1.35446911968096, -1.35446911968096, -1.35446911968096, 
#>     -1.35446911968096, -0.557198119386703, -0.557198119386703, 
```

--------------------------------

### SATmodcomp(largeModel, smallModel, details = 0, eps = sqrt(.Machine$double.eps))

Source: https://hojsgaard.github.io/pbkrtest/reference/sat-modcomp.html

Performs an approximate F-test based on the Satterthwaite approach for lmerMod objects.

```APIDOC
## SATmodcomp

### Description
Performs an approximate F-test based on the Satterthwaite approach for linear mixed models.

### Arguments
- **largeModel** (lmerMod) - Required - An lmerMod model object.
- **smallModel** (lmerMod | matrix | formula) - Required - A smaller model, a restriction matrix, or a model formula to test against.
- **details** (numeric) - Optional - If greater than 0, timing details are printed.
- **eps** (numeric) - Optional - A small number used for numerical stability.

### Usage
SATmodcomp(largeModel, smallModel, details = 0, eps = sqrt(.Machine$double.eps))
```

--------------------------------

### Access the beets dataset

Source: https://hojsgaard.github.io/pbkrtest/reference/data-beets.html

Load the beets dataset into the R environment.

```R
beets
```

--------------------------------

### compare_column_space(X1, X2)

Source: https://hojsgaard.github.io/pbkrtest/reference/compare_column_space.html

Compares the column spaces of two matrices X1 and X2.

```APIDOC
## compare_column_space(X1, X2)

### Description
Compares the column spaces of two matrices with the same number of rows.

### Arguments
- **X1** (matrix) - Required - First matrix.
- **X2** (matrix) - Required - Second matrix.

### Value
- **-1**: C(X1)=C(X2), or the spaces are not nested.
- **0**: C(X1) is contained in C(X2).
- **1**: C(X2) is contained in C(X1).

### Example
```r
A1 <- matrix(c(1,1,1,1,2,3), nrow=3)
A2 <- A1[, 1, drop=FALSE]

compare_column_space(A1, A2)
# [1] 1
compare_column_space(A2, A1)
# [1] 0
compare_column_space(A1, A1)
# [1] -1
```
```

--------------------------------

### Extracting denominator degrees of freedom

Source: https://hojsgaard.github.io/pbkrtest/reference/get_modcomp.html

Demonstrates how to compute a model comparison and extract the denominator degrees of freedom using getKR.

```R
data(beets, package='pbkrtest')
lg <- lmer(sugpct ~ block + sow + harvest + (1|block:harvest), 
              data=beets, REML=FALSE)
sm <- update(lg, .~. - harvest)
modcomp <- KRmodcomp(lg, sm)
getKR(modcomp, "ddf") # get denominator degrees of freedom.
#> [1] 2
```

--------------------------------

### PBmodcomp(largeModel, smallModel, nsim = 1000, ref = NULL, seed = NULL, cl = NULL, details = 0)

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Performs model comparison of nested models using parametric bootstrap methods. It generates a reference distribution of the likelihood ratio test statistic under the null hypothesis.

```APIDOC
## PBmodcomp

### Description
Performs model comparison of nested models using parametric bootstrap methods. The function calculates p-values based on the fraction of simulated LRT-values that are larger or equal to the observed LRT value.

### Arguments
- **largeModel** (model object) - Required - A linear mixed effects model, generalized linear mixed effects model, linear normal model, or generalized linear model. Must be larger than smallModel.
- **smallModel** (model object or matrix) - Required - A model of the same type as largeModel or a restriction matrix.
- **nsim** (integer) - Optional - The number of simulations to form the reference distribution. Default is 1000.
- **ref** (vector) - Optional - Vector containing samples from the reference distribution. If NULL, this is generated using PBrefdist().
- **seed** (integer) - Optional - A seed for the simulation of new datasets.
- **cl** (vector) - Optional - A vector identifying a cluster for parallel computing.
- **details** (integer) - Optional - The amount of output produced, mainly for debugging.
```

--------------------------------

### Perform Mixed Model Hypothesis Tests

Source: https://hojsgaard.github.io/pbkrtest/index.html

Demonstrates fitting linear mixed models with lme4 and comparing them using ANOVA, parametric bootstrap, Kenward-Roger, and Satterthwaite methods.

```R
fm0 <- lmer(sugpct ~ block + sow + harvest + (1|block:harvest), data=beets)
fm1 <- update(fm0, .~. -harvest)

## Is there an effect of harvest time?
an <- anova(fm0, fm1)
pb <- PBmodcomp(fm0, fm1)
kr <- KRmodcomp(fm0, fm1)
sa <- SATmodcomp(fm0, fm1)

tidy(an)
#> # A tibble: 2 × 9
#>   term   npar   AIC   BIC logLik deviance statistic    df   p.value
#>   <chr> <dbl> <dbl> <dbl>  <dbl>    <dbl>     <dbl> <dbl>     <dbl>
#> 1 fm1       9 -69.1 -56.5   43.5    -87.1      NA      NA NA       
#> 2 fm0      10 -80.0 -66.0   50.0   -100.       12.9     1  0.000326
tidy(pb)
#> # A tibble: 2 × 4
#>   type    stat    df  p.value
#>   <chr>  <dbl> <dbl>    <dbl>
#> 1 LRT     12.9     1 0.000326
#> 2 PBtest  12.9    NA 0.0300
tidy(kr)
#> # A tibble: 1 × 6
#>   type   stat   ndf   ddf F.scaling p.value
#>   <chr> <dbl> <int> <dbl>     <dbl>   <dbl>
#> 1 Ftest  15.2     1  2.00         1  0.0599
tidy(sa)
#> # A tibble: 1 × 5
#>   type  statistic   ndf   ddf p.value
#>   <chr>     <dbl> <int> <dbl>   <dbl>
#> 1 Ftest      15.2     1  2.00  0.0599
```

--------------------------------

### Linear Model Comparison

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Comparing linear models using anova and PBmodcomp with different input formats.

```R
m11 <- lm(dist ~ speed + I(speed^2), data=cars)
m10 <- update(m11, ~.-I(speed^2))
anova(m11, m10)
PBmodcomp(m11, m10, cl=1, nsim=NSIM)
PBmodcomp(m11, ~.-I(speed^2), cl=1, nsim=NSIM)
PBmodcomp(m11, c(0, 0, 1), cl=1, nsim=NSIM)
```

--------------------------------

### Perform Asymptotic Likelihood Ratio Test

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Compare models using the ANOVA function with a Chi-squared test, tidying the output for readability.

```R
anova(lmm1.bal, lmm0.bal, test="Chisq")  |> tidy()
```

```R
anova(lmm1.imbal, lmm0.imbal, test="Chisq")  |> tidy()
```

--------------------------------

### seqPBmodcomp(largeModel, smallModel, h = 20, nsim = 1000, cl = 1)

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Performs sequential computing for bootstrap p-values, stopping the sampling process once a specified number of extreme cases is reached.

```APIDOC
## seqPBmodcomp

### Description
Performs sequential computing for bootstrap p-values to determine the number of extreme cases needed before stopping the sampling process.

### Arguments
- **largeModel** (model object) - Required - The larger model object.
- **smallModel** (model object) - Required - The smaller model object.
- **h** (integer) - Optional - The number of extreme cases needed to generate before the sampling process stops. Default is 20.
- **nsim** (integer) - Optional - The number of simulations. Default is 1000.
- **cl** (integer) - Optional - Cluster identifier for parallel computing. Default is 1.
```

--------------------------------

### Generalized Linear Model Comparison

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Comparing generalized linear models using anova and PBmodcomp.

```R
m21 <- glm(dist ~ speed + I(speed^2), family=Gamma("identity"), data=cars)
m20 <- update(m21, ~.-I(speed^2))
anova(m21, m20, test="Chisq")
PBmodcomp(m21, m20, cl=1, nsim=NSIM)
PBmodcomp(m21, ~.-I(speed^2), cl=1, nsim=NSIM)
PBmodcomp(m21, c(0, 0, 1), cl=1, nsim=NSIM)
```

--------------------------------

### Function Signature for compare_column_space

Source: https://hojsgaard.github.io/pbkrtest/reference/compare_column_space.html

The basic syntax for comparing the column spaces of two matrices.

```R
compare_column_space(X1, X2)
```

--------------------------------

### Compare column spaces

Source: https://hojsgaard.github.io/pbkrtest/articles/a02-coercion.html

Evaluate the relationship between the column spaces of two model matrices.

```R
## The first column space contains the second
compare_column_space(Xlg, Xsm)
```

```R
## The second column space contains the first
compare_column_space(Xsm, Xlg)
```

```R
## The two column spaces are identical
compare_column_space(Xlg, Xlg)
```

--------------------------------

### PBrefdist(largeModel, smallModel, nsim, seed, cl, details)

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-refdist.html

Calculates the reference distribution of the likelihood ratio statistic for comparing two nested linear mixed effects models.

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

### Description
Calculates the reference distribution of the likelihood ratio statistic in mixed effects models using parametric bootstrap. The function supports models fitted with lme4.

### Parameters
- **largeModel** (lmerMod) - Required - A linear mixed effects model (larger model).
- **smallModel** (lmerMod) - Required - A linear mixed effects model (smaller model).
- **nsim** (numeric) - Optional - The number of simulations to form the reference distribution (default: 1000).
- **seed** (numeric) - Optional - Seed for the random number generation.
- **cl** (numeric/list) - Optional - Used for controlling parallel computations (integer or cluster list).
- **details** (numeric) - Optional - The amount of output produced for debugging purposes.

### Value
A numeric vector representing the reference distribution.

### Example
```r
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, cl=1)
```
```

--------------------------------

### Summarize and tidy Kenward-Roger model comparison

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Use summary and tidy to extract F-test statistics from a Kenward-Roger model comparison object.

```R
summary(kr.bal) |> tidy()
```

--------------------------------

### Perform Parametric Bootstrap

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Conducts model comparison using parametric bootstrap with parallel processing support.

```R
pb.bal <- PBmodcomp(lmm1.bal, lmm0.bal, nsim=500, cl=2)
pb.bal |> tidy()
```

```R
summary(pb.bal) |> tidy()
```

--------------------------------

### Parallel computation configuration

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-refdist.html

Configuration options for parallel processing using the cl argument.

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

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

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

--------------------------------

### Generalized Linear Model Comparison

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Compare generalized linear models using PBmodcomp.

```R
counts    <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome   <- gl(3, 1, 9)
treatment <- gl(3, 3)
d.AD      <- data.frame(treatment, outcome, counts)
head(d.AD)
glm.D93   <- glm(counts ~ outcome + treatment, family = poisson())
glm.D93.o <- update(glm.D93, .~. -outcome)
glm.D93.t <- update(glm.D93, .~. -treatment)

anova(glm.D93, glm.D93.o, test="Chisq")
PBmodcomp(glm.D93, glm.D93.o, nsim=NSIM, cl=1)
```

--------------------------------

### Perform Parametric Bootstrap Model Comparison

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Executes a parametric bootstrap test on linear mixed models and tidies the output.

```R
pb.imbal <- PBmodcomp(lmm1.imbal, lmm0.imbal, nsim=500, cl=2)
pb.imbal |> tidy()
```

```R
summary(pb.imbal)  |> tidy()
```

--------------------------------

### Linear Normal Model Comparison

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Compare linear normal models using PBmodcomp.

```R
sug <- lm(sugpct ~ block + sow + harvest, data=beets)
sug.h <- update(sug, .~. -harvest)
sug.s <- update(sug, .~. -sow)

anova(sug, sug.h)
PBmodcomp(sug, sug.h, nsim=NSIM, cl=1)

anova(sug, sug.s)
PBmodcomp(sug, sug.s, nsim=NSIM, cl=1)
```

--------------------------------

### Apply Satterthwaite Approximation

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Uses SATmodcomp to perform model comparisons using the Satterthwaite approximation for balanced and unbalanced designs.

```R
sat.bal <- SATmodcomp(lmm1.bal, lmm0.bal)
sat.bal |> tidy()
```

```R
sat.imbal <- SATmodcomp(lmm1.imbal, lmm0.imbal)
sat.imbal |> tidy()
```

--------------------------------

### Function signatures for getKR and getSAT

Source: https://hojsgaard.github.io/pbkrtest/reference/get_modcomp.html

Defines the parameters for extracting components from a KRmodcomp object.

```R
getKR(
  object,
  name = c("ndf", "ddf", "Fstat", "p.value", "F.scaling", "FstatU", "p.valueU", "aux")
)

getSAT(object, name = c("ndf", "ddf", "Fstat", "p.value"))
```

--------------------------------

### getSAT(object, name)

Source: https://hojsgaard.github.io/pbkrtest/reference/get_modcomp.html

Extracts components from a KRmodcomp object using Satterthwaite approximation results.

```APIDOC
## getSAT(object, name)

### Description
Extracts specific components from a KRmodcomp object related to Satterthwaite approximation.

### Arguments
- **object** (KRmodcomp) - Required - The result of the KRmodcomp function.
- **name** (character) - Optional - The specific slot to return. Available options: "ndf", "ddf", "Fstat", "p.value".
```

--------------------------------

### Linear Mixed Effects Model Comparison

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-modcomp.html

Compare linear mixed effects models using PBmodcomp.

```R
sug   <- lmer(sugpct ~ block + sow + harvest + (1|block:harvest),
              data=beets, REML=FALSE)
sug.h <- update(sug, .~. -harvest)
sug.s <- update(sug, .~. -sow)

anova(sug, sug.h)
PBmodcomp(sug, sug.h, nsim=NSIM, cl=1)

anova(sug, sug.s)
PBmodcomp(sug, sug.s, nsim=NSIM, cl=1)
```

--------------------------------

### Apply Kenward-Roger Approximation

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Use the KRmodcomp function to perform the Kenward-Roger approximation for model comparison.

```R
kr.bal <- KRmodcomp(lmm1.bal, lmm0.bal)
kr.bal |> tidy()
```

--------------------------------

### Analyze sugar beets data with ANOVA

Source: https://hojsgaard.github.io/pbkrtest/reference/data-beets.html

Perform an analysis of variance on yield and sugar percentage using the interaction of block and harvest as the error term.

```R
data(beets)

beets$bh <- with(beets, interaction(block, harvest))
summary(aov(yield ~ block + sow + harvest + Error(bh), beets))
#> 
#> Error: bh
#>           Df Sum Sq Mean Sq  F value   Pr(>F)    
#> block      2   11.3     5.7    4.018 0.199292    
#> harvest    1 1725.2  1725.2 1225.000 0.000815 ***
#> Residuals  2    2.8     1.4                      
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Error: Within
#>           Df Sum Sq Mean Sq F value   Pr(>F)    
#> sow        4   8369  2092.3   119.7 1.13e-13 ***
#> Residuals 20    349    17.5                     
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary(aov(sugpct ~ block + sow + harvest + Error(bh), beets))
#> 
#> Error: bh
#>           Df  Sum Sq Mean Sq F value Pr(>F)  
#> block      2 0.03267 0.01633   2.579 0.2794  
#> harvest    1 0.09633 0.09633  15.211 0.0599 .
#> Residuals  2 0.01267 0.00633                 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Error: Within
#>           Df Sum Sq Mean Sq F value   Pr(>F)    
#> sow        4   1.01  0.2525     101 5.74e-13 ***
#> Residuals 20   0.05  0.0025                     
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
```

--------------------------------

### PBrefdist function signature

Source: https://hojsgaard.github.io/pbkrtest/reference/pb-refdist.html

The function signature for PBrefdist and its S3 methods for lm and merMod objects.

```R
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
)
```

--------------------------------

### Perform and tidy imbalanced data model comparison

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Execute KRmodcomp on imbalanced linear mixed models and tidy the resulting output.

```R
kr.imbal <- KRmodcomp(lmm1.imbal, lmm0.imbal)
kr.imbal |> tidy()
```

--------------------------------

### make_model_matrix

Source: https://hojsgaard.github.io/pbkrtest/reference/model-coerce.html

Creates a model matrix based on a provided matrix and restriction matrix.

```APIDOC
## make_model_matrix(X, L)

### Description
Creates a model matrix using a base matrix X and a restriction matrix L.

### Parameters
- **X** (matrix) - Required - The model matrix.
- **L** (matrix) - Required - The restriction matrix.
```

--------------------------------

### Prepare balanced dataframe for mixed models

Source: https://hojsgaard.github.io/pbkrtest/articles/a01-pbkrtest.html

Convert the list-based shoes data into a balanced dataframe suitable for mixed model analysis.

```R
boy <- rep(1:10, 2)
boyf<- factor(letters[boy])
material <- factor(c(rep("A", 10), rep("B", 10)))
## Balanced data:
shoe.bal <- data.frame(wear=unlist(shoes), boy=boy, boyf=boyf, material=material)
head(shoe.bal)
```
Tokens:
Raw