Comparative Algorithm Performance in ‘rctbayespower’ • rctbayespower

Setup:

# Set CRAN mirror
options(repos = c(CRAN = "https://cloud.r-project.org/"))

# Check and install required packages
packages <- c("dplyr", "kableExtra", "pander")
install.packages(setdiff(packages, rownames(installed.packages())))
invisible(lapply(packages, require, character.only = TRUE))

# Load the package
library(rctbayespower)

# Set number of cores and simulations:
n_cores <- parallel::detectCores() - 2
n_sims <- 500

1 Algorithm Performance and Computational Considerations

This vignette provides guidance on algorithm selection, performance optimization, and computational considerations when using the rctbayespower package for Bayesian power analysis.

1.1 Algorithm Overview

The rctbayespower package supports multiple algorithms via brms / Stan:

meanfield (variational inference)
fullrank (variational inference)
sampling (MCMC)

2 Performance Comparison

2.1 Setup

2.1.1 Model

Show predefined models available in the package:

list_predefined_models()

[1] "ancova_cont_2arms" "ancova_cont_3arms"

Get the predefined ANCOVA model for continuous outcomes with two arms:

model_ancova <- build_model(predefined_model = "ancova_cont_2arms")
model_ancova@parameter_names_brms

[1] "b_Intercept" "b_covariate" "b_arm2"

2.1.2 Design

design <- build_design(
  model = model_ancova,
  target_params = "b_arm2",
  thresholds_success = 0.1,
  thresholds_futility = 0,
  p_sig_success = 0.975,
  p_sig_futility = 0.5
)
# check the required parameters for the design
required_fn_args(design)


Arguments that need user specification.

Simulation function:
n_total, b_arm_treat, b_covariate

Interim function:

2.1.3 Conditions

conditions <- build_conditions(
  design = design,
  condition_values = list(
    # two sample sizes
    n_total = 400,
      # two effect sizes
    b_arm_treat = c(0,0.3)
  ),
  static_values = list(
    # equal allocation
    p_alloc = c(0.5, 0.5),
          # baseline effect
    b_covariate = 0
  )
)

2.2 Run Simulations

# Meanfield
  power_result_meanfield <-
    power_analysis(
      conditions = conditions,
      n_cores = n_cores,
      n_sims = n_sims,
      verbose = FALSE,
      brms_args = list(
        algorithm = "meanfield",
        importance_resampling = TRUE,
        iter = 1e4,
        output_samples = 2e3
      )
    )

Chain 1: ------------------------------------------------------------
Chain 1: EXPERIMENTAL ALGORITHM:
Chain 1:   This procedure has not been thoroughly tested and may be unstable
Chain 1:   or buggy. The interface is subject to change.
Chain 1: ------------------------------------------------------------
Chain 1:
Chain 1:
Chain 1:
Chain 1: Gradient evaluation took 2.6e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Begin eta adaptation.
Chain 1: Iteration:   1 / 250 [  0%]  (Adaptation)
Chain 1: Iteration:  50 / 250 [ 20%]  (Adaptation)
Chain 1: Iteration: 100 / 250 [ 40%]  (Adaptation)
Chain 1: Iteration: 150 / 250 [ 60%]  (Adaptation)
Chain 1: Iteration: 200 / 250 [ 80%]  (Adaptation)
Chain 1: Success! Found best value [eta = 1] earlier than expected.
Chain 1:
Chain 1: Begin stochastic gradient ascent.
Chain 1:   iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes
Chain 1:    100          -37.151             1.000            1.000
Chain 1:    200          -36.361             0.511            1.000
Chain 1:    300          -36.550             0.342            0.022
Chain 1:    400          -36.492             0.257            0.022
Chain 1:    500          -36.403             0.206            0.005   MEDIAN ELBO CONVERGED
Chain 1:
Chain 1: Drawing a sample of size 2000 from the approximate posterior...
Chain 1: COMPLETED.

# Fullrank
  power_result_fullrank <-
    power_analysis(
      conditions = conditions,
      n_cores = n_cores,
      n_sims = n_sims,
      verbose = FALSE,
      brms_args = list(
        algorithm = "fullrank",
        importance_resampling = TRUE,
        iter = 1e4,
        output_samples = 2e3
      )
    )

Chain 1: ------------------------------------------------------------
Chain 1: EXPERIMENTAL ALGORITHM:
Chain 1:   This procedure has not been thoroughly tested and may be unstable
Chain 1:   or buggy. The interface is subject to change.
Chain 1: ------------------------------------------------------------
Chain 1:
Chain 1:
Chain 1:
Chain 1: Gradient evaluation took 8e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Begin eta adaptation.
Chain 1: Iteration:   1 / 250 [  0%]  (Adaptation)
Chain 1: Iteration:  50 / 250 [ 20%]  (Adaptation)
Chain 1: Iteration: 100 / 250 [ 40%]  (Adaptation)
Chain 1: Iteration: 150 / 250 [ 60%]  (Adaptation)
Chain 1: Iteration: 200 / 250 [ 80%]  (Adaptation)
Chain 1: Success! Found best value [eta = 1] earlier than expected.
Chain 1:
Chain 1: Begin stochastic gradient ascent.
Chain 1:   iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes
Chain 1:    100          -61.435             1.000            1.000
Chain 1:    200          -38.327             0.801            1.000
Chain 1:    300          -38.504             0.536            0.603
Chain 1:    400          -38.890             0.404            0.603
Chain 1:    500          -38.210             0.327            0.018
Chain 1:    600          -37.965             0.274            0.018
Chain 1:    700          -37.698             0.236            0.010   MEDIAN ELBO CONVERGED
Chain 1:
Chain 1: Drawing a sample of size 2000 from the approximate posterior...
Chain 1: COMPLETED.

# Sampling
power_result_sampling <- power_analysis(
  conditions = conditions,
  n_cores = n_cores,
  n_sims = n_sims,
  verbose = FALSE,
  brms_args = list(
    chains = 4,
    iter = 700,
    warmup = 200
  )
)

Extract the summarized results from each power analysis:

results_summary<- rbind(
  power_result_meanfield@summarized_results |> 
    mutate(Algorithm = "Meanfield"),
  power_result_fullrank@summarized_results |> 
    mutate(Algorithm = "Fullrank"),
  power_result_sampling@summarized_results |> 
    mutate(Algorithm = "Sampling")
)

2.3 Comparison of Results

2.3.1 Probability of Success and Futility and Power

results_summary |> 
  select(
    Algorithm,
    n_total,
    b_arm_treat,
    prob_success,
    prob_futility,
    power_success,
    power_futility
  ) |>
  dplyr::arrange(desc(b_arm_treat)) |> 
  kableExtra::kable(digits = 2)

Algorithm	n_total	b_arm_treat	prob_success	prob_futility	power_success	power_futility
Meanfield	400	0.3	0.89	0.03	0.41	0.01
Fullrank	400	0.3	0.90	0.03	0.40	0.01
Sampling	400	0.3	0.92	0.02	0.55	0.00
Meanfield	400	0.0	0.27	0.48	0.00	0.47
Fullrank	400	0.0	0.27	0.49	0.00	0.48
Sampling	400	0.0	0.24	0.49	0.00	0.49

2.3.2 Convergence

results_summary |> 
  select(
    Algorithm,
    n_total,
    b_arm_treat,
    rhat,
    ess_bulk,
    ess_tail,
    conv_rate
  ) |>
  kableExtra::kable(digits = 3)

Algorithm	n_total	b_arm_treat	rhat	ess_bulk	ess_tail	conv_rate
Meanfield	400	0.0	1.001	959.844	918.042	1
Meanfield	400	0.3	1.001	966.890	917.698	1
Fullrank	400	0.0	1.000	967.312	914.260	1
Fullrank	400	0.3	1.001	960.745	921.348	1
Sampling	400	0.0	1.002	2026.675	1457.781	1
Sampling	400	0.3	1.002	2053.590	1475.123	1

2.3.3 Run Time

# extract elapsed time
data.frame(
  Alorithm = c(
    "Meanfield",
    "Fullrank",
    "Sampling"
  ),
    "Elapsed Time" = c(
      power_result_meanfield@elapsed_time,
      power_result_fullrank@elapsed_time,
      power_result_sampling@elapsed_time
    )
  ) |> 
  kableExtra::kable(digits = 2)

Alorithm	Elapsed.Time
Meanfield	1.94
Fullrank	1.88
Sampling	2.26

3 Session Information

pander::pander(sessionInfo())

R version 4.5.1 (2025-06-13 ucrt)

Platform: x86_64-w64-mingw32/x64

locale: LC_COLLATE=German_Germany.utf8, LC_CTYPE=German_Germany.utf8, LC_MONETARY=German_Germany.utf8, LC_NUMERIC=C and LC_TIME=German_Germany.utf8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: rctbayespower(v.0.0.0.9000), pander(v.0.6.6), kableExtra(v.1.4.0) and dplyr(v.1.1.4)

loaded via a namespace (and not attached): gtable(v.0.3.6), tensorA(v.0.36.2.1), QuickJSR(v.1.8.1), xfun(v.0.53), ggplot2(v.4.0.0), htmlwidgets(v.1.6.4), inline(v.0.3.21), lattice(v.0.22-7), vctrs(v.0.6.5), tools(v.4.5.1), generics(v.0.1.4), stats4(v.4.5.1), parallel(v.4.5.1), tibble(v.3.3.0), pkgconfig(v.2.0.3), brms(v.2.23.0), Matrix(v.1.7-4), data.table(v.1.17.8), checkmate(v.2.3.3), RColorBrewer(v.1.1-3), S7(v.0.2.0), distributional(v.0.5.0), RcppParallel(v.5.1.11-1), lifecycle(v.1.0.4), compiler(v.4.5.1), farver(v.2.1.2), stringr(v.1.5.2), textshaping(v.1.0.4), Brobdingnag(v.1.2-9), codetools(v.0.2-20), htmltools(v.0.5.8.1), bayesplot(v.1.14.0), yaml(v.2.3.10), lazyeval(v.0.2.2), plotly(v.4.11.0), pillar(v.1.11.1), tidyr(v.1.3.1), StanHeaders(v.2.32.10), bridgesampling(v.1.1-2), abind(v.1.4-8), nlme(v.3.1-168), posterior(v.1.6.1), rstan(v.2.32.7), tidyselect(v.1.2.1), digest(v.0.6.37), mvtnorm(v.1.3-3), stringi(v.1.8.7), purrr(v.1.0.4), fastmap(v.1.2.0), grid(v.4.5.1), cli(v.3.6.5), magrittr(v.2.0.3), loo(v.2.8.0), pkgbuild(v.1.4.8), withr(v.3.0.2), scales(v.1.4.0), backports(v.1.5.0), rmarkdown(v.2.30), httr(v.1.4.7), matrixStats(v.1.5.0), gridExtra(v.2.3), pbapply(v.1.7-4), coda(v.0.19-4.1), evaluate(v.1.0.5), knitr(v.1.50), pwr(v.1.3-0), viridisLite(v.0.4.2), rstantools(v.2.5.0), rlang(v.1.1.6), Rcpp(v.1.1.0), glue(v.1.8.0), xml2(v.1.4.0), svglite(v.2.2.2), rstudioapi(v.0.17.1), jsonlite(v.2.0.0), R6(v.2.6.1) and systemfonts(v.1.3.1)