2025-10-26
Source:vignettes/articles/99-model-validation.qmd
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))
library(rctbayespower)
# Set number of cores and simulations:
n_cores <- parallel::detectCores() - 2
n_sims <- 500\[ y_{ij} = \beta_0 + d \cdot \text{arm}_{ij} + \rho \cdot \text{covariate}_{ij} + \epsilon_{ij},\\ \epsilon_{ij} \sim \mathcal{N}(0,1)\]
where \(i\) is the subject and \(j\) is the arm (1 or 2).
True effect size and covariate effect, and sample size range.
rho <- .3
d <- .5
sample_size <- seq(40, 260, 20)Set up model, design, and conditions for Monte Carlo power analysis.
# Create ANCOVA model with uninformative priors
model <- build_model_ancova_cont_2arms(
b_arm_treat = d,
b_covariate = rho,
prior_treatment = brms::set_prior("normal(0, 1e3)", class = "b"),
prior_covariate = brms::set_prior("normal(0, 1e3)", class = "b", coef = "covariate"),
prior_intercept = brms::set_prior("constant(0)", class = "Intercept"),
prior_sigma = brms::set_prior("normal(0, 100)", class = "sigma", lb = 0),
link_sigma = "identity"
)Compiling the brms model ...
Model compilation done!
# Create design
design <- build_design(
model = model,
target_params = "b_arm2",
n_interim_analyses = 0,
thresholds_success = 0,
thresholds_futility = 0,
p_sig_success = 0.95,
p_sig_futility = 0.5
)
# Create conditions
conditions <- build_conditions(
design = design,
condition_values = list(
n_total = sample_size
),
static_values = list())Run the Monte Carlo power analysis.
power <- rctbayespower::power_analysis(
conditions = conditions,
n_sims = n_sims,
n_cores = n_cores,
brms_args = list(
chains = 4,
iter = 1000,
warmup = 250,
init = .5
),
verbose = TRUE
)
=== Power Analysis ===
Conditions:
# A tibble: 12 × 2
id_cond n_total
<int> <dbl>
1 1 40
2 2 60
3 3 80
4 4 100
5 5 120
6 6 140
7 7 160
8 8 180
9 9 200
10 10 220
11 11 240
12 12 260
Conditions to test: 12
Simulations per condition: 500
Total simulations: 6000 Functions exported from package namespace
Combined results dimensions: 6000 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: 5.29 minutes
power_analytic <- analytical_power_ancova_cont_2arms(
n = sample_size,
d = d,
beta_cov = rho,
sigma = 1,
alpha = .05,method = "theory",
covariate_method = "expected",
alternative = "greater"
)
dplyr::tibble(
"N total" = sample_size,
"MC estimate" = power_mc,
"95% CI" = power_ci_95 ,
"Analytical estimate" = power_analytic
) |>
dplyr::mutate(
`MC estimate` = round(`MC estimate`, 3),
`Analytical estimate` = round(`Analytical estimate`, 3),
Difference = round(`MC estimate` - `Analytical estimate`, 3)
) |>
dplyr::arrange(`N total`) |>
kableExtra::kable("html",align = c("c","c", "c", "c", "c"))| N total | MC estimate | 95% CI | Analytical estimate | Difference |
|---|---|---|---|---|
| 40 | 0.422 | [0.379, 0.465] | 0.455 | -0.033 |
| 60 | 0.576 | [0.533, 0.619] | 0.600 | -0.024 |
| 80 | 0.676 | [0.635, 0.717] | 0.711 | -0.035 |
| 100 | 0.788 | [0.752, 0.824] | 0.795 | -0.007 |
| 120 | 0.814 | [0.78, 0.848] | 0.857 | -0.043 |
| 140 | 0.876 | [0.847, 0.905] | 0.901 | -0.025 |
| 160 | 0.920 | [0.896, 0.944] | 0.932 | -0.012 |
| 180 | 0.952 | [0.933, 0.971] | 0.954 | -0.002 |
| 200 | 0.964 | [0.948, 0.98] | 0.969 | -0.005 |
| 220 | 0.976 | [0.963, 0.989] | 0.979 | -0.003 |
| 240 | 0.980 | [0.968, 0.992] | 0.986 | -0.006 |
| 260 | 0.996 | [0.99, 1.002] | 0.991 | 0.005 |
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), processx(v.3.8.6), inline(v.0.3.21), lattice(v.0.22-7), callr(v.3.7.6), ps(v.1.9.1), 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)