Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The powercurve function from the simr package in R (Green & MacLeod, 2016) can incur very long running times when the method used for the calculation of p values is Kenward-Roger or Satterthwaite (see Luke, 2017). Here I suggest three ways for cutting down this time.

1. Where possible, use a high-performance (or high-end) computing cluster. This removes the need to use personal computers for these long jobs.

2. In case you’re using the fixed() parameter of the powercurve function, and calculating the power for different effects, run these at the same time (‘in parallel’) on different machines, rather than one after another.

3. Parallelize the breaks argument. The breaks argument of the powercurve function allows the calculation of power for different levels of the grouping factor passed to along. Some grouping factors are participant, trial and item. The breaks argument sets the different sample sizes for which power will be calculated. Parallelizing breaks is done by running each number of levels in a separate function. When each has been run and saved, they are combined to allow the plotting. This procedure is demonstrated below.

## Parallelizing breaks

Let’s do a minimal example using a toy lmer model. A power curve will be craeted for the fixed effect of x along different sample sizes of the grouping factor g.

Notice that the six sections of the power curve below are serially arranged, one after another. In contrast, to enable parallel processing, each power curve would be placed in a single script, and they would all be run at the same time.

Although the power curves below run in a few minutes, the settings that are often used (e.g., a larger model; fixed('x', 'sa') instead of fixed('x'); nsim = 500 instead of nsim = 50) take far longer. That is where parallel processing becomes useful.1

library(lme4)
library(simr)

# Toy model
fm = lmer(y ~ x + (x | g), data = simdata)

# Extend sample size of g
fm_extended_g = extend(fm, along = 'g', n = 12)

# Parallelize breaks by running each number of levels in a separate function.

# 4 levels of g
pwcurve_4g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 4,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 6 levels of g
pwcurve_6g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 6,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 8 levels of g
pwcurve_8g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 8,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 10 levels of g
pwcurve_10g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 10,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 12 levels of g
pwcurve_12g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 12,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

Having saved each section of the power curve, we must now combine them to be able to plot them together.

# Create a destination object using any of the power curves above.
all_pwcurve = pwcurve_4g

# Combine results
all_pwcurve$ps = c(pwcurve_4g$ps, pwcurve_6g$ps, pwcurve_8g$ps,
pwcurve_10g$ps, pwcurve_12g$ps)

# Combine the different numbers of levels.
all_pwcurve$xval = c(pwcurve_4g$nlevels, pwcurve_6g$nlevels, pwcurve_8g$nlevels,
pwcurve_10g$nlevels, pwcurve_12g$nlevels)

print(all_pwcurve)
## Power for predictor 'x', (95% confidence interval),
## by number of levels in g:
##       4: 46.00% (31.81, 60.68) - 40 rows
##       6: 74.00% (59.66, 85.37) - 60 rows
##       8: 92.00% (80.77, 97.78) - 80 rows
##      10: 98.00% (89.35, 99.95) - 100 rows
##      12: 100.0% (92.89, 100.0) - 120 rows
##
## Time elapsed: 0 h 0 m 7 s
plot(all_pwcurve, xlab = 'Levels of g') # For reproducibility purposes
sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
##  en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
##  stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
##  simr_1.0.5    lme4_1.1-27.1 Matrix_1.3-4
##
## loaded via a namespace (and not attached):
##   Rcpp_1.0.7        lattice_0.20-44   tidyr_1.1.3       assertthat_0.2.1
##   digest_0.6.27     utf8_1.2.2        plyr_1.8.6        R6_2.5.1
##   cellranger_1.1.0  backports_1.2.1   evaluate_0.14     blogdown_1.5
##  pillar_1.6.2      rlang_0.4.11      curl_4.3.2        readxl_1.3.1
##  minqa_1.2.4       data.table_1.14.0 car_3.0-11        nloptr_1.2.2.2
##  jquerylib_0.1.4   rmarkdown_2.11    splines_4.1.1     stringr_1.4.0
##  foreign_0.8-81    broom_0.7.9       compiler_4.1.1    xfun_0.26
##  pkgconfig_2.0.3   mgcv_1.8-36       htmltools_0.5.2   tidyselect_1.1.1
##  tibble_3.1.4      binom_1.1-1       bookdown_0.24     rio_0.5.27
##  fansi_0.5.0       crayon_1.4.1      dplyr_1.0.7       MASS_7.3-54
##  grid_4.1.1        nlme_3.1-153      jsonlite_1.7.2    lifecycle_1.0.0
##  DBI_1.1.1         magrittr_2.0.1    zip_2.2.0         stringi_1.7.4
##  carData_3.0-4     bslib_0.3.0       ellipsis_0.3.2    vctrs_0.3.8
##  generics_0.1.0    boot_1.3-28       openxlsx_4.2.4    iterators_1.0.13
##  tools_4.1.1       forcats_0.5.1     glue_1.4.2        purrr_0.3.4
##  hms_1.1.0         plotrix_3.8-2     parallel_4.1.1    abind_1.4-5
##  pbkrtest_0.5.1    fastmap_1.1.0     yaml_2.2.1        RLRsim_3.1-6
##  knitr_1.34        haven_2.4.3       sass_0.4.0

### Just the code

library(lme4)
library(simr)

# Toy model
fm = lmer(y ~ x + (x | g), data = simdata)

# Extend sample size of g
fm_extended_g = extend(fm, along = 'g', n = 12)

# Parallelize breaks by running each number of levels in a separate function.

# 4 levels of g
pwcurve_4g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 4,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 6 levels of g
pwcurve_6g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 6,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 8 levels of g
pwcurve_8g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 8,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 10 levels of g
pwcurve_10g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 10,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# 12 levels of g
pwcurve_12g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 12,
nsim = 50, seed = 123,
# No progress bar
progress = FALSE)

# Create a destination object using any of the power curves above.
all_pwcurve = pwcurve_4g

# Combine results
all_pwcurve$ps = c(pwcurve_4g$ps, pwcurve_6g$ps, pwcurve_8g$ps,
pwcurve_10g$ps, pwcurve_12g$ps)

# Combine the different numbers of levels.
all_pwcurve$xval = c(pwcurve_4g$nlevels, pwcurve_6g$nlevels, pwcurve_8g$nlevels,
pwcurve_10g$nlevels, pwcurve_12g$nlevels)

print(all_pwcurve)

plot(all_pwcurve, xlab = 'Levels of g')


1. The number of simulations set by nsim should be larger (Brysbaert & Stevens, 2018; Green & MacLeod, 2016). In addition, the effect size for x should be adjusted to the value that best fits with the planned study (Kumle et al., 2021).↩︎