# 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))
library(rctbayespower)
# Set number of cores and simulations:
n_cores <- parallel::detectCores() - 2
n_sims <- 500Setup:
1 Introduction
The rctbayespower package provides tools for conducting Bayesian power analysis for randomized controlled trials (RCTs) using the brms package and Stan.
2 Basic Power Analysis Using ANCOVA Design
2.1 Model
Show predefined models available in the package:
[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.2 Design
Find possible target parameters in the model.
model_ancova@parameter_names_brms[1] "b_Intercept" "b_covariate" "b_arm2" Specify the study design with target parameter, thresholds for success and futility, and significance levels for frequentist-like power.
design <- build_design(
model = model_ancova,
target_params = "b_arm2",
thresholds_success = 0.1,
thresholds_futility = 0,
p_sig_success = 0.9,
p_sig_futility = 0.5
)
print(design)
S7 Object of class: 'rctbp_design'
--------------------------------------------------
=== Model Specifications ===
Number of endpoints: 1
Endpoint types: continuous
Number of arms: 2
Number of repeated measures: 0
Parameter names - simulation function: n_total, n_arms, contrasts, p_alloc, intercept, b_arm_treat, b_covariate, sigma
Parameter names - brms model: b_Intercept, b_covariate, b_arm2
=== Design Specifications ===
Design name: NULL
Target parameters: b_arm2
Number of interim analyses: 0
Thresholds for success: 0.1
Thresholds for futility: 0
Probability of success significance: 0.9
Probability of futility significance: 0.5
Parameter names - interim function: NULL
=== 'brms' Model ===
Family: gaussian
Links: mu = identity
Formula: outcome ~ 1 + covariate + arm
Data: mock_data_ancova (Number of observations: 20)
Draws: 1 chains, each with iter = 500; warmup = 250; thin = 1;
total post-warmup draws = 250
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.09 0.29 -0.47 0.64 1.03 317 151
covariate 0.02 0.33 -0.57 0.69 1.01 177 147
arm2 -0.24 0.36 -1.00 0.49 1.00 418 276
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.08 0.19 0.83 1.50 1.00 148 173
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).2.2.1 Conditions
Find parameters that need user specification for the design.
required_fn_args(design)
Arguments that need user specification.
Simulation function:
n_total, b_arm_treat, b_covariate
Interim function:
Specify conditions.
conditions <- build_conditions(
design = design,
condition_values = list(
# two sample sizes
n_total = seq(80, 160, by = 40),
# two effect sizes
b_arm_treat = c(0, 0.3)
),
static_values = list(
# baseline effect
b_covariate = 0
)
)
print(conditions, n = 100)
S7 Object of class: 'rctbp_conditions'
--------------------------------------------------
Number of conditions: 6
Number of varying parameters: 2
Number of static parameters: 1
Condition Grid:
# A tibble: 6 × 3
id_cond n_total b_arm_treat
<int> <dbl> <dbl>
1 1 80 0
2 2 80 0.3
3 3 120 0
4 4 120 0.3
5 5 160 0
6 6 160 0.32.3 Run Power Analysis
Run the analysis using the power_analysis function. This will run the model for each condition and return the results.
power <- power_analysis(
conditions = conditions,
n_cores = n_cores,
n_sims = n_sims,
verbose = TRUE,
brms_args = list(
chains = 4,
iter = 300,
warmup = 100
)
)
=== Power Analysis ===
Conditions:
# A tibble: 6 × 3
id_cond n_total b_arm_treat
<int> <dbl> <dbl>
1 1 80 0
2 2 80 0.3
3 3 120 0
4 4 120 0.3
5 5 160 0
6 6 160 0.3
Conditions to test: 6
Simulations per condition: 500
Total simulations: 3000 Functions exported from package namespace
Combined results dimensions: 3000 18
Column names: parameter, threshold_success, threshold_futility, success_prob, futility_prob, power_success, power_futility, median, mad, mean, sd, rhat, ess_bulk, ess_tail, id_sim, id_cond, converged, error
Total analysis time: 2.81 minutes2.4 Power for Effect of 0.3
power@summarized_results |>
dplyr::filter(b_arm_treat != 0) |>
select(
n_total,
b_arm_treat,
prob_success,
prob_success_se,
power_success,
power_success_se,
prob_futility,
prob_futility_se,
power_futility,
power_futility_se
) |>
dplyr::arrange(desc(b_arm_treat)) |>
kableExtra::kable(digits = 3, format = "html")| n_total | b_arm_treat | prob_success | prob_success_se | power_success | power_success_se | prob_futility | prob_futility_se | power_futility | power_futility_se |
|---|---|---|---|---|---|---|---|---|---|
| 80 | 0.3 | 0.741 | 0.011 | 0.356 | 0.021 | 0.168 | 0.009 | 0.088 | 0.013 |
| 120 | 0.3 | 0.781 | 0.010 | 0.418 | 0.022 | 0.118 | 0.007 | 0.036 | 0.008 |
| 160 | 0.3 | 0.826 | 0.008 | 0.478 | 0.022 | 0.079 | 0.005 | 0.020 | 0.006 |
2.5 Alpha Error Rate for Null Effect
The power for success indicates the alpha error rate for the null effect (b_arm_treat = 0). This is the probability of rejecting the null hypothesis when it is true, which should be close to the significance level (0.05) if the design is well specified.
power@summarized_results |>
dplyr::filter(b_arm_treat == 0) |>
select(
n_total,
b_arm_treat,
prob_success,
prob_success_se,
power_success,
power_success_se,
power_futility,
power_futility_se,
prob_futility,
prob_futility_se
) |>
dplyr::arrange(desc(b_arm_treat)) |>
kableExtra::kable(digits = 3,
format = "html")| n_total | b_arm_treat | prob_success | prob_success_se | power_success | power_success_se | power_futility | power_futility_se | prob_futility | prob_futility_se |
|---|---|---|---|---|---|---|---|---|---|
| 80 | 0 | 0.383 | 0.012 | 0.034 | 0.008 | 0.492 | 0.022 | 0.491 | 0.013 |
| 120 | 0 | 0.352 | 0.012 | 0.028 | 0.007 | 0.492 | 0.022 | 0.496 | 0.013 |
| 160 | 0 | 0.315 | 0.012 | 0.018 | 0.006 | 0.516 | 0.022 | 0.512 | 0.013 |
2.6 Run Time
Show the code
Total run time: 2.8 minutes for 3000 total simulation repetitions using 14 cores.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), reshape2(v.1.4.4), 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), plyr(v.1.8.9), R6(v.2.6.1) and systemfonts(v.1.3.1)