Bayesian power and design analysis by simulation

Simulation-based design analysis for multilevel SCD models. Data are generated under an assumed treatment effect, the model is refit by update() each replication (reusing one compiled Stan model), and the proportion of fits meeting the chosen posterior decision rule is the power. The same multilevel structure covers a single study (cases) and a meta-analysis of single-case data (cases nested in studies, set n_studies > 1 with between-study effect heterogeneity study_sd).

Usage

run_power_simulation(
  n_sims = 200,
  n_cases = 5,
  n_studies = 1,
  n_phase_a = 5,
  n_phase_b = 5,
  effect_sizes = c(0, 0.5, 1, 1.5),
  autocorr = 0.2,
  icc = 0.2,
  study_sd = 0,
  slope_sd = 0,
  fit_autocor = TRUE,
  iter = 1000,
  chains = 2,
  design = "ab",
  family = "gaussian",
  coefficient = "level",
  decision_rule = "prob_direction",
  prob_threshold = 0.975,
  rope = 0.1,
  standardized = TRUE,
  prior = NULL
)

Arguments

n_sims

Replications per effect magnitude. 200+ is reasonable; 500 gives a Monte Carlo SE below 0.02.

n_cases

Cases per study.

n_studies

Studies. 1 for a single-study design; > 1 adds a study-level random effect (meta-analysis of single-case data).

n_phase_a, n_phase_b

Time points per phase.

effect_sizes

Treatment-effect magnitudes on the model link scale (log-rate ratio for count families). For a Gaussian family with standardized = TRUE the value is multiplied by the baseline (within-case) SD (a Cohen's-d-style planning grid). The element 0 yields the Type I rate.

autocorr

Lag-1 autocorrelation of within-case residuals.

icc

Intraclass correlation; sets the between-case SD via tau = sigma * sqrt(icc / (1 - icc)).

study_sd

Between-study SD of the treatment effect (link scale); used when n_studies > 1.

slope_sd

Between-case SD of the treatment effect, adding case-level random-slope heterogeneity to the data-generating process and fitting a matching uncorrelated case random slope on the tested coefficient. On the link scale for count families; in SD units for a standardized Gaussian family. Default 0 (homogeneous effect): the fitted model then uses only a random intercept, so the effect_size = 0 row reflects zero between-case heterogeneity. Set > 0 to plan for a realistic spread of case-level effects.

fit_autocor

Logical; add a level-1 AR(1) term (ar(time, gr = case)) to the fitted model so it matches the autocorrelated errors the simulator generates. Gaussian family only (default TRUE).

iter, chains

MCMC iterations and chains.

design

One of "ab", "reversal", "multiple_baseline".

family

One of "gaussian", "poisson", "negbinomial".

coefficient

Coefficient carrying the effect and being tested: "level" (b_phase_num), "trend" (b_time_center), or "slope_change" (b_time_since_change).

decision_rule

One of "prob_direction" (probability of direction pd >= prob_threshold), "prob_meaningful" (posterior probability of an effect beyond the ROPE, P(|effect| >= rope) >= prob_threshold), or "rope" (posterior probability inside the ROPE <= 1 - prob_threshold). These mirror the pd, prob_meaningful_change, and rope_prob quantities reported by the analysis decision summaries.

prob_threshold

Posterior-probability threshold the chosen rule must meet to flag an effect. The default 0.975 corresponds, for the prob_direction rule, to a 5\ (pd >= 0.975 mirrors a two-sided alpha of 0.05); 0.95 mirrors a 10\ two-sided rate.

rope

Half-width of the region of practical equivalence, on the same scale as effect_sizes.

standardized

For a Gaussian family, multiply the magnitude by the baseline (within-case) SD (planning in SD units). Ignored for non-Gaussian families.

prior

A brms prior. Defaults to a weakly-informative prior on fixed effects scaled to the outcome (normal(0, residual SD) for Gaussian, normal(0, 1) on the log scale for count families); diffuse priors inflate the Type I rate at small samples.

Value

A data frame with one row per magnitude: effect_size, true_beta (magnitude on the coefficient scale), power (the Type I rate when the magnitude is 0), mc_se (Agresti-Coull Monte Carlo SE of power, which stays positive when power is 0 or 1), type_s (sign-error rate, conditional on the posterior-mean sign) and type_m (posterior-mean exaggeration ratio among fired replications, in the spirit of Gelman and Carlin (2014); note the posterior mean is shrunk by the prior, so this is not their MLE-based quantity), n_valid, n_nonconverged (replications dropped for Rhat > 1.01, which can bias the estimate if non-convergence correlates with effect size), and n_sims.

References

Gelman, A., & Carlin, J. (2014). Beyond power calculations: Assessing Type S (sign) and Type M (magnitude) errors. Perspectives on Psychological Science, 9(6), 641-651.

Kruschke, J. K. (2018). Rejecting or accepting parameter values in Bayesian estimation. Advances in Methods and Practices in Psychological Science, 1(2), 270-280.

Examples

if (interactive()) {
run_power_simulation(n_sims = 50, n_cases = 4,
                     effect_sizes = c(0, 0.5, 1.0),
                     decision_rule = "prob_direction")
}