Template - Power Simulation
The following is an R code that can be used to generate power calculations practically all statistical calculation. Just run the script, find the analysis you want to run and use it with the values you expect.
# =============================================================================
# power_sim_engine.R
# General-purpose simulation-based power analysis engine
#
# Defaults: alpha = 0.05, target power = 0.80
#
#
# Usage:
# source("power_sim_engine.R")
# list_scenarios()
# estimate_power("a01_t_test", n = 64, par = list(d = 0.5), nsims = 2000)
# find_n("f24_cluster_rct_lmm", par = list(m = 15, icc = 0.05, d = 0.3))
# power_curve("a08_ancova", n_values = seq(20, 120, 20))
#
# Dependencies (base/recommended R covers most): MASS, nnet, survival ship
# with R; lme4 + lmerTest needed for F-block; pscl needed for c18 only.
# =============================================================================
`%||%` <- function(a, b) if (is.null(a)) b else a
# ----------------------------------------------------------------------------
# Small helpers
# ----------------------------------------------------------------------------
# Multivariate normal via Cholesky (no external dependency)
.rmvnorm <- function(n, mu, Sigma) {
R <- chol(Sigma)
Z <- matrix(rnorm(n * length(mu)), nrow = n)
sweep(Z %*% R, 2, mu, `+`)
}
.require_pkgs <- function(scn) {
for (pkg in scn$requires) {
if (!requireNamespace(pkg, quietly = TRUE)) {
stop(sprintf("Scenario '%s' requires package '%s' (install.packages(\"%s\")).",
scn$name, pkg, pkg), call. = FALSE)
}
}
}
# ----------------------------------------------------------------------------
# Scenario constructor (public: use this to add your own scenarios)
# ----------------------------------------------------------------------------
make_scenario <- function(name, simulate, fit, extract,
n_label = "n",
description = "",
default_par = list(),
demo_n = 50,
default_n_range = c(4, 4000),
requires = character(0)) {
structure(list(name = name, simulate = simulate, fit = fit, extract = extract,
n_label = n_label, description = description,
default_par = default_par, demo_n = demo_n,
default_n_range = default_n_range, requires = requires),
class = "power_scenario")
}
print.power_scenario <- function(x, ...) {
cat(sprintf("Scenario: %s\n %s\n n = %s\n", x$name, x$description, x$n_label))
if (length(x$requires)) cat(" requires:", paste(x$requires, collapse = ", "), "\n")
cat(" default_par:\n")
for (nm in names(x$default_par)) {
v <- x$default_par[[nm]]
cat(sprintf(" %s = %s\n", nm,
paste(format(v, digits = 4), collapse = ", ")))
}
invisible(x)
}
# ----------------------------------------------------------------------------
# Core engine
# ----------------------------------------------------------------------------
.run_once <- function(scn, n, par, seed = NULL) {
if (!is.null(seed)) set.seed(seed)
had_warning <- FALSE
out <- withCallingHandlers(
tryCatch({
dat <- scn$simulate(n, par)
fit <- scn$fit(dat, par)
p <- scn$extract(fit, dat, par)
p <- as.numeric(p)[1]
if (!is.finite(p)) stop("non-finite p-value returned")
list(p = p, error = NA_character_)
}, error = function(e) list(p = NA_real_, error = conditionMessage(e))),
warning = function(w) { had_warning <<- TRUE; invokeRestart("muffleWarning") },
message = function(m) invokeRestart("muffleMessage")
)
out$warning <- had_warning
out
}
.resolve_scenario <- function(scn) {
if (is.character(scn)) {
if (!scn %in% names(SCENARIOS))
stop(sprintf("Unknown scenario '%s'. See list_scenarios().", scn), call. = FALSE)
scn <- SCENARIOS[[scn]]
}
stopifnot(inherits(scn, "power_scenario"))
scn
}
estimate_power <- function(scn, n, par = list(), nsims = 1000, alpha = 0.05,
seed = NULL, cores = 1) {
scn <- .resolve_scenario(scn)
.require_pkgs(scn)
par <- modifyList(scn$default_par, par)
if (!is.null(seed)) set.seed(seed)
rep_seeds <- sample.int(.Machine$integer.max, nsims)
runner <- function(i) .run_once(scn, n, par, seed = rep_seeds[i])
results <- if (cores > 1 && .Platform$OS.type == "unix") {
parallel::mclapply(seq_len(nsims), runner, mc.cores = cores)
} else {
lapply(seq_len(nsims), runner)
}
p <- vapply(results, function(r) r$p, numeric(1))
warn <- vapply(results, function(r) isTRUE(r$warning), logical(1))
errs <- vapply(results, function(r) r$error, character(1))
ok <- !is.na(p)
sig <- p[ok] < alpha
pw <- if (any(ok)) mean(sig) else NA_real_
structure(list(
scenario = scn$name, n = n, n_label = scn$n_label, par = par,
alpha = alpha, nsims = nsims, n_ok = sum(ok),
failure_rate = mean(!ok), warning_rate = mean(warn),
power = pw,
power_conservative = sum(sig) / nsims,
mc_se = if (any(ok)) sqrt(pw * (1 - pw) / sum(ok)) else NA_real_,
error_messages = if (any(!ok)) sort(table(errs[!ok]), decreasing = TRUE) else NULL,
pvalues = p
), class = "power_estimate")
}
print.power_estimate <- function(x, ...) {
cat(sprintf("Scenario : %s\n", x$scenario))
cat(sprintf("n : %s (%s)\n", format(x$n), x$n_label))
cat(sprintf("alpha : %.3f | nsims = %d\n", x$alpha, x$nsims))
cat(sprintf("Power : %.3f (MC SE %.3f)\n", x$power, x$mc_se))
if (x$failure_rate > 0) {
cat(sprintf("Fit failures : %.1f%% -> conservative power (failures = n.s.): %.3f\n",
100 * x$failure_rate, x$power_conservative))
cat("Top error messages:\n")
em <- head(x$error_messages, 3)
for (i in seq_along(em)) cat(sprintf(" [%dx] %s\n", em[[i]], names(em)[i]))
}
if (x$warning_rate > 0)
cat(sprintf("Replications w/ warnings (e.g. convergence): %.1f%%\n",
100 * x$warning_rate))
invisible(x)
}
power_curve <- function(scn, n_values, par = list(), nsims = 500, alpha = 0.05,
seed = NULL, cores = 1) {
scn <- .resolve_scenario(scn)
if (!is.null(seed)) set.seed(seed)
seeds <- sample.int(.Machine$integer.max, length(n_values))
rows <- lapply(seq_along(n_values), function(i) {
e <- estimate_power(scn, n_values[i], par, nsims, alpha,
seed = seeds[i], cores = cores)
data.frame(n = n_values[i], power = e$power, mc_se = e$mc_se,
failure_rate = e$failure_rate, warning_rate = e$warning_rate)
})
out <- do.call(rbind, rows)
attr(out, "scenario") <- scn$name
attr(out, "n_label") <- scn$n_label
out
}
# Bisection search for the smallest n reaching target power.
# Assumes power is (stochastically) monotone in n; nsims_search controls the
# Monte Carlo noise during the search, nsims_final the precision of the
# confirmatory estimate at the returned n.
find_n <- function(scn, par = list(), target_power = 0.80, alpha = 0.05,
n_range = NULL, nsims_search = 500, nsims_final = 2000,
seed = NULL, cores = 1, verbose = TRUE) {
scn <- .resolve_scenario(scn)
n_range <- n_range %||% scn$default_n_range
lo <- as.integer(ceiling(n_range[1])); hi <- as.integer(floor(n_range[2]))
if (!is.null(seed)) set.seed(seed)
say <- function(...) if (verbose) cat(sprintf(...))
e_hi <- estimate_power(scn, hi, par, nsims_search, alpha, cores = cores)
say("n = %5d : power = %.3f\n", hi, e_hi$power)
if (is.na(e_hi$power) || e_hi$power < target_power) {
warning(sprintf("Power at upper bound n = %d is %.3f < target %.2f; widen n_range.",
hi, e_hi$power, target_power))
return(invisible(list(n_star = NA_integer_, at_bound = TRUE, last = e_hi)))
}
e_lo <- estimate_power(scn, lo, par, nsims_search, alpha, cores = cores)
say("n = %5d : power = %.3f\n", lo, e_lo$power)
if (!is.na(e_lo$power) && e_lo$power >= target_power) {
say("Target already met at lower bound.\n")
hi <- lo
}
while (hi - lo > 1) {
mid <- (lo + hi) %/% 2
e <- estimate_power(scn, mid, par, nsims_search, alpha, cores = cores)
say("n = %5d : power = %.3f\n", mid, e$power)
if (!is.na(e$power) && e$power >= target_power) hi <- mid else lo <- mid
}
final <- estimate_power(scn, hi, par, nsims_final, alpha, cores = cores)
say("--> n* = %d (%s), confirmatory power = %.3f (MC SE %.3f, nsims = %d)\n",
hi, scn$n_label, final$power, final$mc_se, nsims_final)
invisible(list(n_star = hi, target_power = target_power, alpha = alpha,
confirmatory = final))
}
list_scenarios <- function() {
data.frame(
id = names(SCENARIOS),
n_label = vapply(SCENARIOS, function(s) s$n_label, character(1)),
requires = vapply(SCENARIOS, function(s)
paste(s$requires, collapse = ","), character(1)),
description = vapply(SCENARIOS, function(s) s$description, character(1)),
row.names = NULL
)
}
get_scenario <- function(id) .resolve_scenario(id)
# ----------------------------------------------------------------------------
# Reusable response simulators (fork 2: this is where generality lives)
# ----------------------------------------------------------------------------
# Latent-variable simulator for proportional-odds ordinal outcomes.
# base_probs: category probabilities when eta = 0. Returns integer 1..K.
.sim_ordinal <- function(eta, base_probs) {
K <- length(base_probs)
tau <- qlogis(cumsum(base_probs))[-K] # thresholds
ystar <- eta + rlogis(length(eta))
1L + rowSums(outer(ystar, tau, `>`)) # P(Y<=k|eta)=plogis(tau_k - eta)
}
# Survival times with administrative censoring (+ optional exponential dropout)
.censor <- function(T, follow_up, dropout_rate = 0) {
C <- rep(follow_up, length(T))
if (dropout_rate > 0) C <- pmin(C, rexp(length(T), rate = dropout_rate))
data.frame(time = pmin(T, C), event = as.integer(T <= C))
}
# ----------------------------------------------------------------------------
# THE 30 SCENARIOS
# ----------------------------------------------------------------------------
SCENARIOS <- list()
## ===================== A. Continuous / Gaussian ============================
SCENARIOS$a01_t_test <- make_scenario(
"a01_t_test",
description = "Independent two-sample t-test, equal variances (Cohen's d)",
n_label = "participants per group", demo_n = 64,
default_par = list(d = 0.5, sd = 1),
simulate = function(n, par) {
g <- rep(0:1, each = n)
data.frame(g = g, y = rnorm(2 * n, mean = par$d * g, sd = par$sd))
},
fit = function(dat, par) t.test(y ~ g, data = dat, var.equal = TRUE),
extract = function(fit, dat, par) fit$p.value
)
SCENARIOS$a02_welch <- make_scenario(
"a02_welch",
description = "Welch's t-test, unequal per-group SDs",
n_label = "participants per group", demo_n = 80,
default_par = list(m1 = 0, m2 = 0.5, sd1 = 1, sd2 = 2),
simulate = function(n, par) {
data.frame(g = rep(0:1, each = n),
y = c(rnorm(n, par$m1, par$sd1), rnorm(n, par$m2, par$sd2)))
},
fit = function(dat, par) t.test(y ~ g, data = dat, var.equal = FALSE),
extract = function(fit, dat, par) fit$p.value
)
SCENARIOS$a03_paired_t <- make_scenario(
"a03_paired_t",
description = "Paired t-test; bivariate-normal pre/post with correlation rho",
n_label = "pairs", demo_n = 50,
default_par = list(d = 0.4, sd1 = 1, sd2 = 1, rho = 0.5),
simulate = function(n, par) {
S <- matrix(c(par$sd1^2, par$rho * par$sd1 * par$sd2,
par$rho * par$sd1 * par$sd2, par$sd2^2), 2)
X <- .rmvnorm(n, c(0, par$d), S)
data.frame(y1 = X[, 1], y2 = X[, 2])
},
fit = function(dat, par) t.test(dat$y2, dat$y1, paired = TRUE),
extract = function(fit, dat, par) fit$p.value
)
SCENARIOS$a04_anova_omnibus <- make_scenario(
"a04_anova_omnibus",
description = "One-way ANOVA, omnibus F over k group means",
n_label = "participants per group", demo_n = 50,
default_par = list(means = c(0, 0.2, 0.5), sd = 1),
simulate = function(n, par) {
k <- length(par$means)
g <- factor(rep(seq_len(k), each = n))
data.frame(g = g, y = rnorm(n * k, mean = rep(par$means, each = n), sd = par$sd))
},
fit = function(dat, par) lm(y ~ g, data = dat),
extract = function(fit, dat, par) anova(fit)["g", "Pr(>F)"]
)
SCENARIOS$a05_contrast <- make_scenario(
"a05_contrast",
description = "Single planned contrast among k groups (t on contrast, pooled MSE)",
n_label = "participants per group", demo_n = 50,
default_par = list(means = c(0, 0.2, 0.5), sd = 1, contrast = c(-1, 0.5, 0.5)),
simulate = function(n, par) {
stopifnot(length(par$contrast) == length(par$means),
abs(sum(par$contrast)) < 1e-8)
k <- length(par$means)
g <- factor(rep(seq_len(k), each = n))
data.frame(g = g, y = rnorm(n * k, mean = rep(par$means, each = n), sd = par$sd))
},
fit = function(dat, par) lm(y ~ g, data = dat),
extract = function(fit, dat, par) {
gm <- tapply(dat$y, dat$g, mean)
ni <- tapply(dat$y, dat$g, length)
mse <- summary(fit)$sigma^2
cv <- par$contrast
tstat <- sum(cv * gm) / sqrt(mse * sum(cv^2 / ni))
2 * pt(-abs(tstat), df = sum(ni) - length(gm))
}
)
SCENARIOS$a06_slr <- make_scenario(
"a06_slr",
description = "Simple linear regression, slope of a continuous predictor",
n_label = "participants", demo_n = 90,
default_par = list(beta = 0.3, sd_x = 1, sd_e = 1),
simulate = function(n, par) {
x <- rnorm(n, sd = par$sd_x)
data.frame(x = x, y = par$beta * x + rnorm(n, sd = par$sd_e))
},
fit = function(dat, par) lm(y ~ x, data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["x", "Pr(>|t|)"]
)
SCENARIOS$a07_multiple_regression <- make_scenario(
"a07_multiple_regression",
description = "Multiple regression: target beta adjusted for correlated covariates",
n_label = "participants", demo_n = 120,
default_par = list(beta_target = 0.2, beta_covs = c(0.2, 0.2, 0.2),
r_predictors = 0.3, sd_e = 1),
simulate = function(n, par) {
p <- 1 + length(par$beta_covs)
R <- matrix(par$r_predictors, p, p); diag(R) <- 1
X <- .rmvnorm(n, rep(0, p), R)
colnames(X) <- paste0("x", seq_len(p))
y <- as.vector(X %*% c(par$beta_target, par$beta_covs)) + rnorm(n, sd = par$sd_e)
data.frame(y = y, X)
},
fit = function(dat, par) lm(y ~ ., data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["x1", "Pr(>|t|)"]
)
SCENARIOS$a08_ancova <- make_scenario(
"a08_ancova",
description = "ANCOVA: RCT post-test adjusted for baseline (pre-post corr rho)",
n_label = "participants per group", demo_n = 60,
default_par = list(d = 0.3, rho = 0.5),
simulate = function(n, par) {
trt <- rep(0:1, each = n)
pre <- rnorm(2 * n)
post <- par$d * trt + par$rho * pre + rnorm(2 * n, sd = sqrt(1 - par$rho^2))
data.frame(trt = trt, pre = pre, post = post)
},
fit = function(dat, par) lm(post ~ trt + pre, data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["trt", "Pr(>|t|)"]
)
SCENARIOS$a09_factorial_interaction <- make_scenario(
"a09_factorial_interaction",
description = "2x2 factorial: power for the interaction (moderation) term",
n_label = "participants per cell", demo_n = 60,
default_par = list(cell_means = rbind(c(0, 0), c(0, 0.4)), sd = 1),
simulate = function(n, par) {
d <- expand.grid(a = 0:1, b = 0:1)
d <- d[rep(seq_len(4), each = n), ]
mu <- par$cell_means[cbind(d$a + 1, d$b + 1)]
data.frame(A = factor(d$a), B = factor(d$b), y = rnorm(4 * n, mu, par$sd))
},
fit = function(dat, par) lm(y ~ A * B, data = dat),
extract = function(fit, dat, par) anova(fit)["A:B", "Pr(>F)"]
)
SCENARIOS$a10_correlation <- make_scenario(
"a10_correlation",
description = "Pearson correlation on bivariate-normal data",
n_label = "pairs", demo_n = 85,
default_par = list(rho = 0.3),
simulate = function(n, par) {
X <- .rmvnorm(n, c(0, 0), matrix(c(1, par$rho, par$rho, 1), 2))
data.frame(x = X[, 1], y = X[, 2])
},
fit = function(dat, par) cor.test(dat$x, dat$y),
extract = function(fit, dat, par) fit$p.value
)
## ===================== B. Binary outcomes ==================================
SCENARIOS$b11_two_proportions <- make_scenario(
"b11_two_proportions",
description = "Two-proportion chi-square (prop.test, no continuity correction)",
n_label = "participants per group", demo_n = 120,
default_par = list(p1 = 0.15, p2 = 0.30, correct = FALSE),
simulate = function(n, par) {
data.frame(x = c(rbinom(1, n, par$p1), rbinom(1, n, par$p2)), n = c(n, n))
},
fit = function(dat, par) prop.test(dat$x, dat$n, correct = par$correct),
extract = function(fit, dat, par) fit$p.value
)
SCENARIOS$b12_logistic_binary <- make_scenario(
"b12_logistic_binary",
description = "Logistic regression, binary predictor (odds ratio), Wald p",
n_label = "participants per group", demo_n = 150,
default_par = list(p0 = 0.15, OR = 2),
simulate = function(n, par) {
trt <- rep(0:1, each = n)
p <- plogis(qlogis(par$p0) + log(par$OR) * trt)
data.frame(trt = trt, y = rbinom(2 * n, 1, p))
},
fit = function(dat, par) glm(y ~ trt, family = binomial, data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["trt", "Pr(>|z|)"]
)
SCENARIOS$b13_logistic_continuous <- make_scenario(
"b13_logistic_continuous",
description = "Logistic regression, continuous predictor (OR per SD of x)",
n_label = "participants", demo_n = 200,
default_par = list(p_at_mean = 0.20, OR_per_sd = 1.5),
simulate = function(n, par) {
x <- rnorm(n)
p <- plogis(qlogis(par$p_at_mean) + log(par$OR_per_sd) * x)
data.frame(x = x, y = rbinom(n, 1, p))
},
fit = function(dat, par) glm(y ~ x, family = binomial, data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["x", "Pr(>|z|)"]
)
SCENARIOS$b14_logistic_interaction <- make_scenario(
"b14_logistic_interaction",
description = "Logistic 2x2: interaction OR (effect modification), Wald p",
n_label = "participants per cell", demo_n = 200,
default_par = list(p00 = 0.15, OR_a = 1.5, OR_b = 1.5, OR_ab = 2),
simulate = function(n, par) {
d <- expand.grid(a = 0:1, b = 0:1)
d <- d[rep(seq_len(4), each = n), ]
eta <- qlogis(par$p00) + log(par$OR_a) * d$a + log(par$OR_b) * d$b +
log(par$OR_ab) * d$a * d$b
data.frame(a = d$a, b = d$b, y = rbinom(4 * n, 1, plogis(eta)))
},
fit = function(dat, par) glm(y ~ a * b, family = binomial, data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["a:b", "Pr(>|z|)"]
)
SCENARIOS$b15_mcnemar <- make_scenario(
"b15_mcnemar",
description = "McNemar's test on paired binary data (discordant-pair DGP)",
n_label = "pairs", demo_n = 200,
default_par = list(p11 = 0.20, p10 = 0.15, p01 = 0.05, correct = FALSE),
simulate = function(n, par) {
p00 <- 1 - par$p11 - par$p10 - par$p01
stopifnot(p00 >= 0)
counts <- rmultinom(1, n, c(par$p11, par$p10, par$p01, p00))
data.frame(n11 = counts[1], n10 = counts[2], n01 = counts[3], n00 = counts[4])
},
fit = function(dat, par) {
tab <- matrix(c(dat$n11, dat$n10, dat$n01, dat$n00), 2, byrow = TRUE)
if (dat$n10 + dat$n01 == 0) stop("no discordant pairs")
mcnemar.test(tab, correct = par$correct)
},
extract = function(fit, dat, par) fit$p.value
)
## ===================== C. Count outcomes ===================================
SCENARIOS$c16_poisson_rate <- make_scenario(
"c16_poisson_rate",
description = "Poisson regression with offset: rate ratio, varying exposure",
n_label = "participants per group", demo_n = 100,
default_par = list(rate0 = 0.5, RR = 1.5, exposure_range = c(0.5, 1.5)),
simulate = function(n, par) {
trt <- rep(0:1, each = n)
expo <- runif(2 * n, par$exposure_range[1], par$exposure_range[2])
mu <- par$rate0 * par$RR^trt * expo
data.frame(trt = trt, exposure = expo, y = rpois(2 * n, mu))
},
fit = function(dat, par)
glm(y ~ trt + offset(log(exposure)), family = poisson, data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["trt", "Pr(>|z|)"]
)
SCENARIOS$c17_negbin <- make_scenario(
"c17_negbin",
description = "Negative binomial regression; dispersion theta is an explicit input",
n_label = "participants per group", demo_n = 100,
default_par = list(mu0 = 2, RR = 1.4, theta = 1.2),
requires = "MASS",
simulate = function(n, par) {
trt <- rep(0:1, each = n)
data.frame(trt = trt,
y = rnbinom(2 * n, mu = par$mu0 * par$RR^trt, size = par$theta))
},
fit = function(dat, par) MASS::glm.nb(y ~ trt, data = dat),
extract = function(fit, dat, par) summary(fit)$coefficients["trt", "Pr(>|z|)"]
)
SCENARIOS$c18_zip <- make_scenario(
"c18_zip",
description = "Zero-inflated Poisson: effect on the count component (pscl::zeroinfl)",
n_label = "participants per group", demo_n = 150,
default_par = list(pi_zero = 0.3, mu0 = 1.5, RR = 1.6),
requires = "pscl",
simulate = function(n, par) {
trt <- rep(0:1, each = n)
z <- rbinom(2 * n, 1, par$pi_zero)
y <- ifelse(z == 1, 0L, rpois(2 * n, par$mu0 * par$RR^trt))
data.frame(trt = trt, y = y)
},
fit = function(dat, par) pscl::zeroinfl(y ~ trt | 1, data = dat),
extract = function(fit, dat, par)
summary(fit)$coefficients$count["trt", "Pr(>|z|)"]
)
## ===================== D. Ordinal / categorical ============================
SCENARIOS$d19_prop_odds <- make_scenario(
"d19_prop_odds",
description = "Proportional-odds ordinal regression (latent-variable DGP), LRT",
n_label = "participants per group", demo_n = 100,
default_par = list(base_probs = c(0.2, 0.3, 0.3, 0.2), OR = 1.8),
requires = "MASS",
simulate = function(n, par) {
trt <- rep(0:1, each = n)
y <- .sim_ordinal(log(par$OR) * trt, par$base_probs)
data.frame(trt = trt,
y = factor(y, levels = seq_along(par$base_probs), ordered = TRUE))
},
fit = function(dat, par) {
list(f1 = MASS::polr(y ~ trt, data = dat, Hess = FALSE),
f0 = MASS::polr(y ~ 1, data = dat, Hess = FALSE))
},
extract = function(fit, dat, par)
pchisq(fit$f0$deviance - fit$f1$deviance, df = 1, lower.tail = FALSE)
)
SCENARIOS$d20_multinomial <- make_scenario(
"d20_multinomial",
description = "Multinomial logistic: LRT for predictor over all K-1 log-odds",
n_label = "participants per group", demo_n = 150,
default_par = list(probs0 = c(0.5, 0.3, 0.2), beta = c(0.4, 0.8)),
requires = "nnet",
simulate = function(n, par) {
K <- length(par$probs0)
stopifnot(length(par$beta) == K - 1)
eta0 <- log(par$probs0 / par$probs0[1]) # baseline-cat logits
eta1 <- eta0 + c(0, par$beta)
p1g <- exp(eta1) / sum(exp(eta1))
y <- c(sample.int(K, n, replace = TRUE, prob = par$probs0),
sample.int(K, n, replace = TRUE, prob = p1g))
data.frame(trt = rep(0:1, each = n), y = factor(y, levels = seq_len(K)))
},
fit = function(dat, par) {
list(f1 = nnet::multinom(y ~ trt, data = dat, trace = FALSE),
f0 = nnet::multinom(y ~ 1, data = dat, trace = FALSE))
},
extract = function(fit, dat, par) {
df <- length(levels(dat$y)) - 1
pchisq(fit$f0$deviance - fit$f1$deviance, df = df, lower.tail = FALSE)
}
)
SCENARIOS$d21_chisq_rxc <- make_scenario(
"d21_chisq_rxc",
description = "r x c chi-square test of independence on a contingency table",
n_label = "total N", demo_n = 300,
default_par = list(probs = rbind(c(0.25, 0.15, 0.10), c(0.15, 0.15, 0.20))),
simulate = function(n, par) {
counts <- rmultinom(1, n, as.vector(par$probs))
as.data.frame(matrix(counts, nrow = nrow(par$probs)))
},
fit = function(dat, par) chisq.test(as.matrix(dat), correct = FALSE),
extract = function(fit, dat, par) fit$p.value
)
## ===================== E. Time-to-event ====================================
SCENARIOS$e22_cox <- make_scenario(
"e22_cox",
description = "Cox PH: hazard ratio with administrative censoring + dropout",
n_label = "participants per group", demo_n = 150,
default_par = list(hr = 0.65, base_rate = 0.10, follow_up = 24,
dropout_rate = 0.01),
requires = "survival",
simulate = function(n, par) {
trt <- rep(0:1, each = n)
T <- rexp(2 * n, rate = par$base_rate * par$hr^trt)
cbind(trt = trt, .censor(T, par$follow_up, par$dropout_rate))
},
fit = function(dat, par)
survival::coxph(survival::Surv(time, event) ~ trt, data = dat),
extract = function(fit, dat, par)
summary(fit)$coefficients["trt", "Pr(>|z|)"]
)
SCENARIOS$e23_weibull_aft <- make_scenario(
"e23_weibull_aft",
description = "Parametric Weibull AFT with censoring (follow-up as design input)",
n_label = "participants per group", demo_n = 150,
default_par = list(shape = 1.3, scale0 = 20, aft_ratio = 1.5, follow_up = 24,
dropout_rate = 0),
requires = "survival",
simulate = function(n, par) {
trt <- rep(0:1, each = n)
T <- rweibull(2 * n, shape = par$shape,
scale = par$scale0 * par$aft_ratio^trt)
cbind(trt = trt, .censor(T, par$follow_up, par$dropout_rate))
},
fit = function(dat, par)
survival::survreg(survival::Surv(time, event) ~ trt, data = dat,
dist = "weibull"),
extract = function(fit, dat, par) summary(fit)$table["trt", "p"]
)
## ===================== F. Multilevel / clustered ===========================
SCENARIOS$f24_cluster_rct_lmm <- make_scenario(
"f24_cluster_rct_lmm",
description = "Cluster-randomized RCT, random-intercept LMM (Satterthwaite p)",
n_label = "clusters per arm", demo_n = 20,
default_n_range = c(3, 400),
default_par = list(m = 15, icc = 0.05, d = 0.30),
requires = c("lme4", "lmerTest"),
simulate = function(n, par) {
J <- 2 * n
cluster <- rep(seq_len(J), each = par$m)
trt <- rep(rep(0:1, each = n), each = par$m)
u <- rnorm(J, sd = sqrt(par$icc))
y <- par$d * trt + u[cluster] + rnorm(J * par$m, sd = sqrt(1 - par$icc))
data.frame(cluster = factor(cluster), trt = trt, y = y)
},
fit = function(dat, par) lmerTest::lmer(y ~ trt + (1 | cluster), data = dat),
extract = function(fit, dat, par)
summary(fit)$coefficients["trt", "Pr(>|t|)"]
)
SCENARIOS$f25_longitudinal_lmm <- make_scenario(
"f25_longitudinal_lmm",
description = "Longitudinal LMM: time x treatment slope difference, random slopes",
n_label = "subjects per arm", demo_n = 60,
default_n_range = c(6, 1000),
default_par = list(timepoints = 5, slope_time = 0.2, slope_diff = 0.5,
sd_int = 0.7, sd_slope = 0.4, cor_is = 0.2, sd_e = 0.6),
requires = c("lme4", "lmerTest"),
simulate = function(n, par) {
N <- 2 * n; Tn <- par$timepoints
time <- rep(seq(0, 1, length.out = Tn), times = N)
id <- rep(seq_len(N), each = Tn)
trt <- rep(rep(0:1, each = n), each = Tn)
S <- matrix(c(par$sd_int^2, par$cor_is * par$sd_int * par$sd_slope,
par$cor_is * par$sd_int * par$sd_slope, par$sd_slope^2), 2)
U <- .rmvnorm(N, c(0, 0), S)
y <- (par$slope_time + par$slope_diff * trt + U[id, 2]) * time +
U[id, 1] + rnorm(N * Tn, sd = par$sd_e)
data.frame(id = factor(id), trt = trt, time = time, y = y)
},
fit = function(dat, par)
lmerTest::lmer(y ~ time * trt + (time | id), data = dat),
extract = function(fit, dat, par)
summary(fit)$coefficients["time:trt", "Pr(>|t|)"]
)
SCENARIOS$f26_glmm_binary <- make_scenario(
"f26_glmm_binary",
description = "Clustered binary GLMM (glmer, Wald p); convergence rate reported",
n_label = "clusters per arm", demo_n = 25,
default_n_range = c(3, 400),
default_par = list(m = 20, p0 = 0.25, OR = 1.8, sd_u = 0.5),
requires = "lme4",
simulate = function(n, par) {
J <- 2 * n
cluster <- rep(seq_len(J), each = par$m)
trt <- rep(rep(0:1, each = n), each = par$m)
u <- rnorm(J, sd = par$sd_u)
p <- plogis(qlogis(par$p0) + log(par$OR) * trt + u[cluster])
data.frame(cluster = factor(cluster), trt = trt,
y = rbinom(J * par$m, 1, p))
},
fit = function(dat, par)
lme4::glmer(y ~ trt + (1 | cluster), data = dat, family = binomial),
extract = function(fit, dat, par)
summary(fit)$coefficients["trt", "Pr(>|z|)"]
)
SCENARIOS$f27_stepped_wedge <- make_scenario(
"f27_stepped_wedge",
description = "Stepped-wedge cluster design: trt adjusted for period effects",
n_label = "subjects per cluster-period", demo_n = 10,
default_n_range = c(2, 400),
default_par = list(n_clusters = 12, n_periods = 5, d = 0.30, icc = 0.05,
period_drift = 0.05),
requires = c("lme4", "lmerTest"),
simulate = function(n, par) {
J <- par$n_clusters; P <- par$n_periods; S <- P - 1
seq_id <- ((seq_len(J) - 1) %% S) + 1 # crossover sequence per cluster
d <- expand.grid(unit = seq_len(n), period = seq_len(P), cluster = seq_len(J))
d$trt <- as.integer(d$period > seq_id[d$cluster])
u <- rnorm(J, sd = sqrt(par$icc))
d$y <- par$d * d$trt + par$period_drift * (d$period - 1) +
u[d$cluster] + rnorm(nrow(d), sd = sqrt(1 - par$icc))
d$cluster <- factor(d$cluster); d$period <- factor(d$period)
d
},
fit = function(dat, par)
lmerTest::lmer(y ~ trt + period + (1 | cluster), data = dat),
extract = function(fit, dat, par)
summary(fit)$coefficients["trt", "Pr(>|t|)"]
)
SCENARIOS$f28_three_level <- make_scenario(
"f28_three_level",
description = "Three-level nested (pupils < classes < schools), school-level trt",
n_label = "schools per arm", demo_n = 12,
default_n_range = c(3, 300),
default_par = list(classes_per_school = 3, pupils_per_class = 12,
icc_school = 0.04, icc_class = 0.08, d = 0.30),
requires = c("lme4", "lmerTest"),
simulate = function(n, par) {
Js <- 2 * n; Jc <- par$classes_per_school; m <- par$pupils_per_class
school <- rep(seq_len(Js), each = Jc * m)
class <- rep(rep(seq_len(Jc), each = m), times = Js)
trt <- rep(rep(0:1, each = n), each = Jc * m)
us <- rnorm(Js, sd = sqrt(par$icc_school))
uc <- rnorm(Js * Jc, sd = sqrt(par$icc_class))
cid <- (school - 1) * Jc + class
e_sd <- sqrt(1 - par$icc_school - par$icc_class)
y <- par$d * trt + us[school] + uc[cid] + rnorm(Js * Jc * m, sd = e_sd)
data.frame(school = factor(school), class = factor(class), trt = trt, y = y)
},
fit = function(dat, par)
lmerTest::lmer(y ~ trt + (1 | school / class), data = dat),
extract = function(fit, dat, par)
summary(fit)$coefficients["trt", "Pr(>|t|)"]
)
## ===================== G. Robustness / non-standard ========================
SCENARIOS$g29_wilcoxon <- make_scenario(
"g29_wilcoxon",
description = "Wilcoxon rank-sum on a lognormal DGP (test != DGP by design)",
n_label = "participants per group", demo_n = 60,
default_par = list(sdlog = 0.8, delta_log = 0.4),
simulate = function(n, par) {
g <- rep(0:1, each = n)
data.frame(g = g, y = rlnorm(2 * n, meanlog = par$delta_log * g,
sdlog = par$sdlog))
},
fit = function(dat, par) wilcox.test(y ~ g, data = dat, exact = FALSE),
extract = function(fit, dat, par) fit$p.value
)
SCENARIOS$g30_mediation <- make_scenario(
"g30_mediation",
description = "Mediation a*b indirect effect; percentile-bootstrap p-value",
n_label = "participants", demo_n = 120,
default_par = list(a = 0.3, b = 0.3, cprime = 0.1, nboot = 500),
simulate = function(n, par) {
x <- rnorm(n)
m <- par$a * x + rnorm(n, sd = sqrt(max(1e-8, 1 - par$a^2)))
y <- par$b * m + par$cprime * x + rnorm(n)
data.frame(x = x, m = m, y = y)
},
fit = function(dat, par) {
ab <- function(idx) {
x <- dat$x[idx]; m <- dat$m[idx]; y <- dat$y[idx]
a <- cov(x, m) / var(x)
X <- cbind(1, x, m)
b <- solve(crossprod(X), crossprod(X, y))[3]
a * b
}
n <- nrow(dat)
boot <- vapply(seq_len(par$nboot),
function(i) ab(sample.int(n, n, replace = TRUE)),
numeric(1))
list(ab_hat = ab(seq_len(n)), boot = boot)
},
extract = function(fit, dat, par) {
B <- length(fit$boot)
# two-sided percentile-bootstrap p (add-one smoothed, never exactly 0)
min(1, 2 * min((sum(fit$boot <= 0) + 1) / (B + 1),
(sum(fit$boot >= 0) + 1) / (B + 1)))
}
)