Skip to contents

Simulate cluster-mass statistics via bootstrapped permutations

Usage

permute_timewise_statistics(
  jlmer_spec,
  family = c("gaussian", "binomial"),
  statistic = c("t", "chisq"),
  nsim = 100L,
  predictors = NULL,
  ...
)

Arguments

jlmer_spec

Data prepped for jlmer from make_jlmer_spec()

family

A GLM family. Currently supports "gaussian" and "binomial".

statistic

Test statistic for calculating cluster mass. Can be one of "t" (default) from the regression model output or "chisq" from a likelihood ratio test (takes about twice as long to calculate).

nsim

Number of simulations description

predictors

(Optional) a subset of predictors to test. Defaults to NULL which tests all predictors.

...

Optional arguments passed to Julia for model fitting. Defaults to fast = TRUE (when family = "binomial") and progress = FALSE.

Value

A simulation-by-time-by-predictor 3D array of cluster statistics, of class timewise_statistics.

Examples

# \donttest{
# \dontshow{
options("jlmerclusterperm.nthreads" = 2)
jlmerclusterperm_setup(cache_dir = tempdir(), verbose = FALSE)
julia_progress(show = FALSE)
# }

library(dplyr, warn.conflicts = FALSE)

# Specification object
spec <- make_jlmer_spec(
  weight ~ 1 + Diet, filter(ChickWeight, Time <= 20),
  subject = "Chick", time = "Time"
)
spec
#> ── jlmer specification ───────────────────────────────────────── <jlmer_spec> ──
#> Formula: weight ~ 1 + Diet2 + Diet3 + Diet4
#> Predictors:
#>   Diet: Diet2, Diet3, Diet4
#> Groupings:
#>   Subject: Chick
#>   Trial:
#>   Time: Time
#> Data:
#> # A tibble: 533 × 6
#>   weight Diet2 Diet3 Diet4 Chick  Time
#>    <dbl> <dbl> <dbl> <dbl> <ord> <dbl>
#> 1     42     0     0     0 1         0
#> 2     51     0     0     0 1         2
#> 3     59     0     0     0 1         4
#> # ℹ 530 more rows
#> ────────────────────────────────────────────────────────────────────────────────

# Simulation x Time x Predictor array of t-statistics from regression output
reset_rng_state()
null_statistics <- permute_timewise_statistics(spec, nsim = 3)
round(null_statistics, 2)
#> , , Predictor = Diet2
#> 
#>    Time
#> Sim     0    2    4    6    8   10   12   14   16   18   20
#>   1 -0.90 0.41 0.39 1.12 1.58 2.13 2.34 2.58 2.47 2.67 2.61
#>   2  1.40 0.90 0.96 1.44 1.61 1.45 1.23 0.90 0.80 0.78 0.76
#>   3 -1.44 1.23 0.73 0.88 0.93 0.69 0.58 0.61 0.77 1.04 1.43
#> 
#> , , Predictor = Diet3
#> 
#>    Time
#> Sim    0     2     4     6     8    10    12    14    16    18    20
#>   1 0.00 -0.21  0.11 -0.68 -0.66 -0.38 -0.01 -0.02 -0.17 -0.20 -0.12
#>   2 0.70  0.56 -0.10  0.17  0.37  0.37  0.55  0.21  0.46  0.44  0.02
#>   3 0.72  1.90  0.85  0.20  0.19 -0.09 -0.34 -0.24  0.34  0.44  0.83
#> 
#> , , Predictor = Diet4
#> 
#>    Time
#> Sim     0     2     4     6     8    10    12   14    16    18    20
#>   1  0.45 -0.82 -0.33 -0.68 -0.36 -0.35 -0.41 0.19 -0.04 -0.16 -0.16
#>   2  2.09 -0.69 -0.10  0.34  0.56  0.63  0.35 0.05 -0.05 -0.26 -0.43
#>   3 -2.16 -1.38 -1.35 -0.85 -0.02 -0.12  0.01 0.59  0.71  0.67  0.53
#> 

# Collect as dataframe with `tidy()`
permuted_timewise_stats_df <- tidy(null_statistics)
permuted_timewise_stats_df
#> # A tibble: 99 × 4
#>    predictor  time statistic sim  
#>    <chr>     <dbl>     <dbl> <fct>
#>  1 Diet2         0 -8.98e- 1 1    
#>  2 Diet2         0  1.40e+ 0 2    
#>  3 Diet2         0 -1.44e+ 0 3    
#>  4 Diet3         0  2.34e-15 1    
#>  5 Diet3         0  6.98e- 1 2    
#>  6 Diet3         0  7.21e- 1 3    
#>  7 Diet4         0  4.49e- 1 1    
#>  8 Diet4         0  2.09e+ 0 2    
#>  9 Diet4         0 -2.16e+ 0 3    
#> 10 Diet2         2  4.12e- 1 1    
#> # ℹ 89 more rows

# Permutations ran under the same RNG state are identical
reset_rng_state()
null_statistics2 <- permute_timewise_statistics(spec, nsim = 3)
identical(null_statistics, null_statistics2)
#> [1] TRUE

get_rng_state()
#> [1] 111
null_statistics3 <- permute_timewise_statistics(spec, nsim = 3)
identical(null_statistics, null_statistics3)
#> [1] FALSE

# \dontshow{
JuliaConnectoR::stopJulia()
# }
# }