Package {mixpower}

Simulation-Based Power Analysis for Mixed-Effects Models

Description

mixpower is a simulation-based toolkit for power and sample-size analysis for linear and generalized linear mixed-effects models (LMMs and GLMMs). It is design-first (no pilot data required), supports Gaussian, binomial, Poisson, and negative binomial families via lme4; Wald and likelihood-ratio tests; multi-parameter sensitivity grids; power curves and minimum sample-size solvers; parallel evaluation with deterministic seeds; and full reproducibility (manifests, result bundling, export to CSV/JSON). Every run reports diagnostics (failure rate, singular-fit rate, effective N).

Details

Typical workflow: (1) define a design with mp_design() (cluster sizes, trials per cell) and effect-size assumptions with mp_assumptions(); (2) build a scenario with a backend constructor (e.g. mp_scenario_lme4() for Gaussian LMM); (3) run mp_power() for a single power estimate, mp_sensitivity() to vary parameters, or mp_power_curve() / mp_solve_sample_size() for curves and sample-size. Use a fixed seed for reproducibility. Failure and singular-fit rates are always reported and never suppressed.

Getting started

After loading the package, define a design and assumptions, build an lme4 scenario, then call mp_power(). See Examples below for a minimal runnable workflow.

Function overview

Design and assumptions: mp_design(), mp_assumptions().

Scenarios (LMM/GLMM): mp_scenario(), mp_scenario_lme4(), mp_scenario_lme4_binomial(), mp_scenario_lme4_poisson(), mp_scenario_lme4_nb().

Power and sensitivity: mp_power() (optional aggregate = "streaming"), mp_sensitivity(), mp_sensitivity_parallel(), mp_power_curve(), mp_power_curve_parallel(), mp_solve_sample_size(), mp_grid_sample_size(), mp_quick_power().

Backends and simulators: mp_backend(), validate_mp_backend(), mp_backend_lme4(), mp_backend_lme4_binomial(), mp_backend_lme4_poisson(), mp_backend_lme4_nb(), mp_backend_glmmtmb(), mp_scenario_glmmtmb_lmm(), simulate_glmm_binomial_data(), simulate_glmm_poisson_data(), simulate_glmm_nb_data().

Optional tidyverse (Suggests): tibble::as_tibble() and ggplot2::autoplot() methods.

Reproducibility and reporting: mp_manifest(), mp_bundle_results(), mp_report_table(), mp_write_results().

Plotting: plot() for mp_power_curve and mp_sensitivity objects.

Getting the most out of mixpower

• Use a fixed seed (e.g. seed = 123) so runs are reproducible.

• Check failure_rate and singular_rate in results; investigate if high.

• For nested model comparison use LRT with an explicit null_formula.

• Use mp_power_curve() or mp_solve_sample_size() to choose sample size.

• Vignettes give step-by-step guides: vignette("mixpower-intro", package = "mixpower"), vignette("mixpower-design", package = "mixpower"), vignette("mixpower-simulations", package = "mixpower"), vignette("mixpower-diagnostics", package = "mixpower"), vignette("mixpower-reproducibility", package = "mixpower"), vignette("mixpower-extending", package = "mixpower").

Author(s)

Maintainer: Alex Litovchenko al4877@columbia.edu

References

Bates D, Maechler M, Bolker B, Walker S (2015). "Fitting Linear Mixed-Effects Models Using lme4." Journal of Statistical Software, 67(1), 1–48. doi:10.18637/jss.v067.i01.

Green P, MacLeod CJ (2016). "SIMR: an R package for power analysis of generalized linear mixed models by simulation." Methods in Ecology and Evolution, 7(4), 493–498. doi:10.1111/2041-210X.12504.

See Also

Entry points: mp_design(), mp_power(), mp_quick_power(). Vignettes: vignette(package = "mixpower").

Examples

# Minimal workflow: design -> assumptions -> scenario -> power
d <- mp_design(clusters = list(subject = 30), trials_per_cell = 4)
a <- mp_assumptions(
  fixed_effects = list(`(Intercept)` = 0, condition = 0.3),
  residual_sd = 1
)
scn <- mp_scenario_lme4(y ~ condition + (1 | subject), design = d, assumptions = a)
res <- mp_power(scn, nsim = 50, seed = 123)
summary(res)

Coerce mixpower results to a tibble

Description

Requires the tibble package (installed with the tidyverse).

Usage

## S3 method for class 'mp_power'
as_tibble(x, ...)

Arguments

Value

A tibble::tibble().

ggplot2 diagnostic plot for sensitivity or power curve

Description

Requires ggplot2. Intended as an optional alternative to base plot.mp_sensitivity() / plot.mp_power_curve().

Usage

## S3 method for class 'mp_sensitivity'
autoplot(
  object,
  ...,
  y = c("estimate", "failure_rate", "singular_rate", "n_effective")
)

Arguments

Value

A ggplot2 object.

Effect-size converters for eliciting assumptions

Description

Helpers that translate standardized or published effect sizes into the raw coefficients and variance components that mp_assumptions() expects, so you do not have to hand-pick regression coefficients. They compose directly:

Usage

mp_d_to_beta(d, sd = 1)

mp_beta_to_d(beta, sd = 1)

mp_r2_to_beta(r2, sd_resid = 1, predictor_sd = 1)

mp_beta_to_r2(beta, sd_resid = 1, predictor_sd = 1)

mp_icc_to_sd(icc, sd_resid = 1)

mp_sd_to_icc(sd, sd_resid = 1)

mp_or_to_logodds(or)

mp_logodds_to_or(beta)

mp_t_to_beta(t, se)

mp_f_to_beta(f, se)

Arguments

Details

mp_assumptions(
  fixed_effects = list("(Intercept)" = 0, condition = mp_d_to_beta(0.5, sd = 1.12)),
  random_effects = list(subject = list(intercept_sd = mp_icc_to_sd(0.1, 1))),
  residual_sd = 1
)

• mp_d_to_beta() / mp_beta_to_d(): Cohen's d for a 0/1 predictor is the mean difference divided by sd, so the coefficient is d * sd.

• mp_r2_to_beta() / mp_beta_to_r2(): for a predictor with SD predictor_sd, a target (partial) variance-explained r2 against residual SD sd_resid implies beta = sqrt(r2/(1-r2)) * sd_resid / predictor_sd.

• mp_icc_to_sd() / mp_sd_to_icc(): an intraclass correlation icc with residual SD sd_resid implies a random-intercept SD sqrt(icc/(1-icc)) * sd_resid.

• mp_or_to_logodds() / mp_logodds_to_or(): a binomial GLMM coefficient is the log odds ratio.

• mp_t_to_beta() / mp_f_to_beta(): recover a coefficient from a published t (or one-numerator-df F) statistic and its standard error.

Value

A numeric scalar.

Examples

mp_d_to_beta(0.5, sd = 1.12)
mp_icc_to_sd(0.1, sd_resid = 1)
mp_r2_to_beta(0.02, sd_resid = 1)
mp_or_to_logodds(1.5)

Fit a model for a single simulated dataset

Description

Fit a model for a single simulated dataset

Usage

fit_model(data, formula)

Arguments

Value

A fitted model object.

Create modeling assumptions for simulation-based power

Description

Assumptions encode effect sizes and the variance components used to simulate data. Random-effect sizes are given as standard deviations on the linear predictor (the scale lme4 reports), via random_effects.

Usage

mp_assumptions(
  fixed_effects,
  random_effects = NULL,
  icc = NULL,
  residual_sd = NULL,
  notes = NULL
)

Arguments

Value

An object of class mp_assumptions.

Examples

a <- mp_assumptions(
  fixed_effects = list("(Intercept)" = 0, condition = 0.4),
  random_effects = list(subject = list(intercept_sd = 0.5)),
  residual_sd = 1
)
a

MixPower backend contract

Description

A backend is a list with three functions used by mp_power():

simulate_fun

First argument receives the mp_scenario (often named scenario); optional seed. Must return a data.frame.

fit_fun

Two arguments: simulated data.frame, then the mp_scenario (names may differ; mp_power() passes them positionally).

test_fun

Two arguments: fitted model, then the mp_scenario; returns ⁠list(p_value = <numeric scalar>)⁠.

Usage

mp_backend(
  simulate_fun,
  fit_fun,
  test_fun,
  name = "custom",
  version = NULL,
  notes = NULL,
  capabilities = NULL
)

Arguments

Details

Use mp_backend() to build a validated object of class mp_backend. Custom backends can also be plain lists with those three names; validate_mp_backend() checks the contract without requiring the class.

Value

An object of class c("mp_backend", "list") with the components above.

Build a glmmTMB backend for Gaussian LMM scenarios

Description

Fits with the glmmTMB function glmmTMB() (Gaussian). Useful for comparing simulation-based power to mp_backend_lme4() and for workflows that later extend to families supported by glmmTMB but not lme4.

Usage

mp_backend_glmmtmb(
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_method = c("wald", "lrt"),
  null_formula = NULL
)

Arguments

Value

An object of class mp_backend.

Build an lme4 backend for MixPower scenarios

Description

Build an lme4 backend for MixPower scenarios

Usage

mp_backend_lme4(
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_method = c("wald", "lrt", "satterthwaite", "kenward-roger", "pb"),
  null_formula = NULL,
  pb_nsim = 100L
)

Arguments

Value

A list containing simulate_fun, fit_fun, and test_fun.

Build an lme4 backend for binomial GLMM scenarios

Description

Build an lme4 backend for binomial GLMM scenarios

Usage

mp_backend_lme4_binomial(
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_method = c("wald", "lrt", "pb"),
  null_formula = NULL,
  pb_nsim = 100L
)

Arguments

Value

A list containing simulate_fun, fit_fun, and test_fun.

Build an lme4 backend for Negative Binomial GLMM scenarios

Description

Build an lme4 backend for Negative Binomial GLMM scenarios

Usage

mp_backend_lme4_nb(
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_method = c("wald", "lrt", "pb"),
  null_formula = NULL,
  pb_nsim = 100L
)

Arguments

Value

A list containing simulate_fun, fit_fun, and test_fun.

Build an lme4 backend for Poisson GLMM scenarios

Description

Build an lme4 backend for Poisson GLMM scenarios

Usage

mp_backend_lme4_poisson(
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_method = c("wald", "lrt", "pb"),
  null_formula = NULL,
  pb_nsim = 100L
)

Arguments

Value

A list containing simulate_fun, fit_fun, and test_fun.

Bundle results with manifest and optional labels

Description

Combines a single result object (mp_power, mp_sensitivity, or mp_power_curve), a reproducibility manifest, and optional user labels into one object. Diagnostics and result structure are retained unchanged.

Usage

mp_bundle_results(
  result,
  manifest,
  study_id = NULL,
  analyst = NULL,
  notes = NULL
)

Arguments

Value

An object of class mp_bundle with components result, manifest, and labels (list with study_id, analyst, notes).

Check the Type I error calibration of a scenario's test

Description

Runs the scenario under the null hypothesis — the focal fixed effect set to zero — and estimates the empirical Type I error rate, i.e. the proportion of replicates in which the (false) effect is declared significant. A trustworthy test rejects at approximately alpha. The estimate is compared to alpha using an exact (Clopper-Pearson) interval, giving a verdict:

Usage

mp_calibrate(
  scenario,
  term = NULL,
  nsim = 1000,
  alpha = 0.05,
  seed = NULL,
  conf_level = 0.95,
  failure_policy = c("count_as_nondetect", "exclude")
)

Arguments

Details

• "well-calibrated": the interval contains alpha.

• "anti-conservative": the interval lies entirely above alpha (the test rejects too often — e.g. a Wald test with few clusters, or a model that omits a random slope that is actually present in the data-generating process). Power computed with such a test is not trustworthy.

• "conservative": the interval lies entirely below alpha.

This is the recommended sanity check to run before trusting a power number: if a design/analysis combination does not control Type I error, its power is meaningless.

Value

An object of class mp_calibration: a list with term, alpha, type1 (empirical Type I rate), ci, conf_level, mcse, nsim, verdict, and the underlying power_result.

See Also

mp_recommend_method(), mp_power().

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  d <- mp_design(list(subject = 40), trials_per_cell = 6)
  a <- mp_assumptions(
    fixed_effects = list("(Intercept)" = 0, condition = 0.4),
    random_effects = list(subject = list(intercept_sd = 0.5)),
    residual_sd = 1
  )
  scn <- mp_scenario_lme4(y ~ condition + (1 | subject), design = d, assumptions = a)
  mp_calibrate(scn, nsim = 200, seed = 1)
}

Compare analysis models on the same simulated data

Description

Simulates data once per replicate from the first scenario, then fits and tests every supplied scenario on that same dataset. This is the analogue of powerlmm's sim_formula: because the models see identical data, differences in their rejection rates isolate the effect of the analysis choice. Use it to study power across competing specifications, or to expose Type I inflation from a misspecified model (e.g. dropping a random slope that is present in the data-generating process).

Usage

mp_compare_models(
  scenarios,
  nsim,
  alpha = 0.05,
  seed = NULL,
  conf_level = 0.95,
  failure_policy = c("count_as_nondetect", "exclude")
)

Arguments

Details

All scenarios must share the same data-generating process (design and assumptions); they should differ only in their analysis model (formula / random-effects structure / test). The first scenario's simulate_fun drives data generation.

Value

An object of class mp_model_comparison with a results data frame (one row per model: power, conf_low, conf_high, failure_rate, n_effective, nsim).

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  d <- mp_design(list(subject = 30), trials_per_cell = 8)
  a <- mp_assumptions(
    fixed_effects = list("(Intercept)" = 0, condition = 0),
    random_effects = list(subject = list(intercept_sd = 0.5,
                                          slopes = list(condition = 0.8))),
    residual_sd = 1
  )
  maximal <- mp_scenario_lme4(y ~ condition + (1 + condition | subject), d, a)
  reduced <- mp_scenario_lme4(y ~ condition + (1 | subject), d, a)
  mp_compare_models(list(maximal = maximal, reduced = reduced),
                    nsim = 50, seed = 1)
}

Create a study design specification

Description

mp_design() encodes how data will be collected: cluster sizes and repeated measurements. It does not encode effect sizes or analysis decisions.

Usage

mp_design(
  clusters,
  trials_per_cell = 1,
  predictors = NULL,
  nesting = NULL,
  notes = NULL
)

Arguments

Value

An object of class mp_design.

Examples

d <- mp_design(clusters = list(subject = 40), trials_per_cell = 10)
d

# Three-level: 8 sites, 5 subjects per site, 4 trials each.
d3 <- mp_design(
  clusters = list(site = 8, subject = 5),
  trials_per_cell = 4,
  nesting = c(subject = "site")
)

Scale a fitted-model scenario's sample size up or down

Description

Sets target level counts for one or more grouping factors of a scenario built with mp_from_fit(), so power can be evaluated at a sample size different from the pilot's. This is the analogue of simr::extend(): levels are cloned from the pilot's within-level structure with fresh ids, and fresh random effects are drawn from the fitted covariance, so the extended data represent a larger (or smaller) sample from the same population.

Usage

mp_extend(scenario, ...)

Arguments

Details

Pair with mp_power() for a single N, or with mp_power_curve() / mp_solve_sample_size() using the ⁠extend.<group>⁠ key for a curve over N.

Value

The scenario with its extend targets set.

See Also

mp_from_fit().

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  m <- lme4::lmer(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy)
  scn <- mp_from_fit(m, test_term = "Days")
  big <- mp_extend(scn, Subject = 40)
  mp_power(big, nsim = 20, seed = 1)
}

Build a power scenario from a fitted lme4 model

Description

Turns an existing lmer/glmer fit (e.g. from pilot or published data) into an mp_scenario, so its estimated effects and variance components inform a simulation-based power analysis. New responses are simulated from the fitted model with stats::simulate() (keeping the estimated random-effect structure and residual variance), the model is refit, and the focal term is tested.

Usage

mp_from_fit(
  fit,
  test_term = NULL,
  test_method = NULL,
  null_formula = NULL,
  pb_nsim = 100L,
  extend = NULL
)

Arguments

Details

The fixed effects used to simulate are read from the scenario's assumptions, which start at the fitted coefficients. This means mp_sensitivity() / mp_power_curve() can vary an effect size (e.g. fixed_effects.condition) for data-based vs smallest-effect-of-interest comparisons. Sample size can be scaled up or down from the pilot with mp_extend() (or the ⁠extend.<group>⁠ sensitivity key), which clones the pilot's structure with fresh levels and fresh random effects.

Value

An object of class mp_scenario.

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  m <- lme4::lmer(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy)
  scn <- mp_from_fit(m, test_term = "Days")
  mp_power(scn, nsim = 20, seed = 1)
}

Create a grid of values for sample-size search

Description

Returns a numeric vector suitable for mp_solve_sample_size()'s grid argument. Specify either the number of points (length.out) or the step size (by); bounds are always explicit (from, to).

Usage

mp_grid_sample_size(from, to, length.out = NULL, by = NULL)

Arguments

Value

A numeric vector. For integer cluster sizes, use round(mp_grid_sample_size(...)) or pass by as an integer (e.g. by = 10).

Examples

mp_grid_sample_size(20, 100, length.out = 9)
mp_grid_sample_size(20, 100, by = 10)

Reproducibility manifest for power analyses

Description

Captures scenario fingerprint, seed strategy, session info, timestamp, and optional git SHA so results can be reproduced or audited. Output is a plain list (and one-row data frame via as.data.frame()) suitable for saving alongside results.

Usage

mp_manifest(scenario, seed = NULL, session = TRUE)

Arguments

Value

A list with components: scenario_digest, seed, seed_strategy, timestamp, r_version, mixpower_version, session_info (if requested), git_sha (if in a git repo). Use as.data.frame() on the list for a single-row table (list components become columns where possible).

Generate a methods paragraph for a power analysis

Description

Produces a ready-to-edit prose description of a simulation-based power analysis from an mp_power() result: the design, model, effect tested, inference method, number of simulations, and the estimated power with its interval. Intended to seed the "Power analysis" paragraph of a methods section.

Usage

mp_methods_text(result, software = TRUE)

Arguments

Value

A length-1 character string with class mp_methods_text (its print() method word-wraps the paragraph).

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  d <- mp_design(list(subject = 30), trials_per_cell = 6)
  a <- mp_assumptions(list("(Intercept)" = 0, condition = 0.4),
                      random_effects = list(subject = list(intercept_sd = 0.5)),
                      residual_sd = 1)
  scn <- mp_scenario_lme4(y ~ condition + (1 | subject), design = d, assumptions = a)
  mp_methods_text(mp_power(scn, nsim = 50, seed = 1))
}

Add a missing-data / dropout mechanism to a scenario

Description

Wraps a scenario so that, on every replicate, observations are deleted from the simulated data before the model is fit. This lets a power analysis reflect realistic incomplete data (Gallop & Liu, 2017; Magnusson, 2018). Three mechanisms are supported:

Usage

mp_missing(
  scenario,
  mechanism = c("mcar", "mar", "dropout"),
  prob = NULL,
  on = NULL,
  slope = 0,
  time = NULL,
  dropout = NULL,
  subject = "subject"
)

Arguments

Details

• "mcar": each observation is deleted independently with probability prob (missing completely at random).

• "mar": each observation is deleted with a probability that depends on an observed column on through a logistic model plogis(qlogis(prob) + slope * on) (missing at random).

• "dropout": monotone longitudinal dropout along time within each subject — once a subject drops out it contributes no later observations. The dropout pattern is given either by dropout as a vector of cumulative dropout proportions (one per ordered timepoint) or as list(shape =, scale =) for a Weibull dropout time on the time scale.

Value

The scenario with its simulator wrapped to apply missingness.

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  d <- mp_design(list(subject = 30), trials_per_cell = 6,
                 predictors = list(time = "continuous"))
  a <- mp_assumptions(list("(Intercept)" = 0, time = 0.4),
                      random_effects = list(subject = list(intercept_sd = 0.5)),
                      residual_sd = 1)
  scn <- mp_scenario_lme4(y ~ time + (1 | subject), design = d,
                          assumptions = a, predictor = "time")
  scn_drop <- mp_missing(scn, "dropout", time = "time",
                         dropout = c(0, 0.1, 0.2, 0.35, 0.5, 0.6))
  mp_power(scn_drop, nsim = 20, seed = 1)
}

Create a PSOCK cluster with mixpower loaded on workers

Description

Caller is responsible for parallel::stopCluster() (typically via on.exit(..., add = TRUE)).

Usage

mp_parallel_cluster(workers)

Arguments

Value

A cluster object from parallel::makeCluster().

Simulation-based power estimation (engine-agnostic core)

Description

mp_power() runs repeated simulations under a scenario and estimates power for the scenario's test decision rule (typically p < alpha).

Usage

mp_power(
  scenario,
  nsim,
  alpha = 0.05,
  seed = NULL,
  failure_policy = c("count_as_nondetect", "exclude"),
  keep = c("minimal", "fits", "data"),
  conf_level = 0.95,
  ci_method = c("clopper-pearson", "wald"),
  aggregate = c("full", "streaming")
)

Arguments

Details

In Phase 4 core, the scenario must provide engine functions: simulate_fun, fit_fun, and test_fun. Later phases will supply defaults based on specific backends (e.g., lme4).

Diagnostics include Type S (wrong-sign) and Type M (exaggeration-ratio) errors among significant replicates (Gelman & Carlin, 2014), computed when the tested term's true effect is known and non-zero and the backend reports an estimate.

Value

An object of class mp_power.

Examples

# A tiny toy engine (not mixed models) just to demonstrate the workflow:
d <- mp_design(list(subject = 30), trials_per_cell = 1)
a <- mp_assumptions(list(condition = 0.3), residual_sd = 1)

sim_fun <- function(scn, seed) {
  n <- scn$design$clusters$subject
  x <- stats::rbinom(n, 1, 0.5)
  y <- scn$assumptions$fixed_effects$condition * x +
    stats::rnorm(n, sd = scn$assumptions$residual_sd)
  data.frame(y = y, condition = x)
}
fit_fun <- function(dat, scn) stats::lm(scn$formula, data = dat)
test_fun <- function(fit, scn) {
  sm <- summary(fit)
  p <- sm$coefficients["condition", "Pr(>|t|)"]
  list(p_value = as.numeric(p))
}

s <- mp_scenario(
  y ~ condition, d, a,
  simulate_fun = sim_fun,
  fit_fun = fit_fun,
  test_fun = test_fun
)
res <- mp_power(s, nsim = 50, seed = 1)
summary(res)

Resumable, checkpointed power simulation

Description

Runs mp_power() in batches, saving the accumulated per-replicate results to file (an .rds) after each batch. If file already exists for the same seed, the run resumes from where it left off and only the remaining replicates are simulated. This makes very large or long-running power analyses robust to interruption, and lets you grow nsim later without recomputing finished replicates.

Usage

mp_power_checkpoint(
  scenario,
  nsim,
  file,
  batch_size = 100,
  alpha = 0.05,
  seed = NULL,
  failure_policy = c("count_as_nondetect", "exclude"),
  conf_level = 0.95,
  ci_method = c("clopper-pearson", "wald"),
  progress = TRUE
)

Arguments

Details

Because replicate seeds are deterministic (seed + i - 1), the result is identical to a single mp_power(scenario, nsim, seed = seed) call: batching only changes when work happens, not what is computed.

Value

An object of class mp_power for all nsim replicates.

See Also

mp_power().

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  d <- mp_design(list(subject = 30), trials_per_cell = 6)
  a <- mp_assumptions(list("(Intercept)" = 0, condition = 0.4),
                      random_effects = list(subject = list(intercept_sd = 0.5)),
                      residual_sd = 1)
  scn <- mp_scenario_lme4(y ~ condition + (1 | subject), design = d, assumptions = a)
  f <- tempfile(fileext = ".rds")
  mp_power_checkpoint(scn, nsim = 60, file = f, batch_size = 20, seed = 1)
}

Power curve for a single design/assumption parameter

Description

Runs mp_power() across a one-dimensional grid of values for one parameter (e.g. cluster size) via mp_sensitivity(). Results include power estimates and per-grid-point diagnostics: failure rate, singular rate, and effective N.

Usage

mp_power_curve(
  scenario,
  vary,
  nsim = 100,
  alpha = 0.05,
  seed = NULL,
  failure_policy = c("count_as_nondetect", "exclude"),
  conf_level = 0.95
)

Arguments

Value

An object of class mp_power_curve with components vary, grid, results (estimate, mcse, conf_low, conf_high, failure_rate, singular_rate, n_effective, nsim, plus the varying parameter column), alpha, failure_policy, and conf_level.

Parallel power curve evaluation

Description

Evaluates power over a one-parameter grid by running mp_power() for each grid cell in parallel. Uses explicit per-cell seeds (seed + cell_index - 1L) so results are deterministic and match serial mp_power_curve() for the same seed. Does not modify mp_power(); parallelization is at the scenario-grid level only.

Usage

mp_power_curve_parallel(
  scenario,
  vary,
  workers = 2L,
  nsim = 100,
  alpha = 0.05,
  seed = NULL,
  failure_policy = c("count_as_nondetect", "exclude"),
  conf_level = 0.95,
  progress = FALSE,
  ...
)

Arguments

Value

An object of class mp_power_curve (same structure as mp_power_curve()).

Note

Parallel execution requires the parallel package (base R) and that mixpower is installed (e.g. install.packages() or devtools::install()) so that workers can load it.

Quick power run for a single LMM design

Description

One-call wrapper that builds design, assumptions, scenario (lme4 LMM), and runs mp_power(). Intended for the common case: one fixed effect, one random intercept. All arguments are explicit; pass further options to mp_power() via ... (e.g. failure_policy, conf_level, keep).

Usage

mp_quick_power(
  formula,
  clusters,
  trials_per_cell = 1,
  fixed_effects,
  residual_sd,
  nsim,
  alpha = 0.05,
  seed = NULL,
  random_effects = NULL,
  icc = NULL,
  test_method = c("wald", "lrt"),
  null_formula = NULL,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  ...
)

Arguments

Value

The result of mp_power() (object of class mp_power).

Examples

mp_quick_power(
  y ~ condition + (1 | subject),
  clusters = list(subject = 40),
  trials_per_cell = 8,
  fixed_effects = list(`(Intercept)` = 0, condition = 0.3),
  residual_sd = 1,
  nsim = 50,
  seed = 123
)

Recommend an inference method for a scenario

Description

Heuristic guidance on the test_method for a scenario, based on the number of levels of the random grouping factors. Wald (normal-approximation z/t) tests and, to a lesser extent, likelihood-ratio tests are known to be anti-conservative when the number of clusters is small (Luke, 2017): the degrees of freedom are overstated, so the test rejects too often. With few clusters, a degrees-of-freedom-corrected test (Satterthwaite or Kenward-Roger, for linear mixed models) or a parametric bootstrap (any family) controls Type I error far better.

Usage

mp_recommend_method(scenario, small_clusters = 30L)

Arguments

Details

This is a fast, design-based heuristic; to measure a specific design and method, use mp_calibrate().

Value

An object of class mp_recommendation: a list with method (the scenario's current method), n_groups (smallest grouping-factor size), is_lmm, caution (logical), recommended (character vector), and rationale.

See Also

mp_calibrate().

Examples

d <- mp_design(list(subject = 12), trials_per_cell = 8)
a <- mp_assumptions(
  fixed_effects = list("(Intercept)" = 0, condition = 0.4),
  random_effects = list(subject = list(intercept_sd = 0.5)),
  residual_sd = 1
)
scn <- mp_scenario_lme4(y ~ condition + (1 | subject), design = d, assumptions = a)
mp_recommend_method(scn)

Publication-ready summary table for power results

Description

Returns a flat data frame with power estimate, CI, failure/singularity rates, and effective simulation counts. Works with mp_power, mp_sensitivity, mp_power_curve, or the result of mp_bundle_results() (uses the bundled result).

Usage

mp_report_table(x, ...)

Arguments

Value

A data frame: for mp_power one row; for mp_calibration a one-row Type I summary; for sensitivity/curve one row per grid cell with parameter column(s), power_estimate, ci_low, ci_high, failure_rate, singular_rate, n_effective, nsim.

Safeguard (confidence-bound) effect size from a fitted model

Description

Computes a safeguard effect for power analysis: the bound of a confidence interval for a fitted effect that lies closest to zero (Perugini, Gallucci & Costantini, 2014). Planning power around this conservative, uncertainty-aware value protects against the optimism of using a noisy pilot point estimate.

Usage

mp_safeguard_effect(fit, term = NULL, conf_level = 0.9)

Arguments

Details

The interval is the Wald (normal-approximation) interval from the fitted coefficient and its standard error. Pair the result with mp_from_fit() and mp_sesoi() to run a safeguard-power simulation.

Value

An object of class mp_safeguard: a list with term, estimate, se, conf_level, the interval (lower, upper), and the safeguard bound (the interval limit nearest zero, in the direction of the estimate).

See Also

mp_sesoi(), mp_from_fit().

Examples

if (requireNamespace("lme4", quietly = TRUE)) {
  m <- lme4::lmer(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy)
  sg <- mp_safeguard_effect(m, term = "Days", conf_level = 0.90)
  sg
}

Create a power-analysis scenario

Description

A scenario combines: (1) a design, (2) assumptions, (3) a model specification, and (4) an analysis engine.

Usage

mp_scenario(
  formula,
  design,
  assumptions,
  test = c("wald", "lrt", "custom"),
  simulate_fun = NULL,
  fit_fun = NULL,
  test_fun = NULL,
  notes = NULL
)

Arguments

Details

In Phase 1, the engine is pluggable via three functions:

• simulate_fun(scenario, seed) returns a data.frame

• fit_fun(data, scenario) returns a fit object

• test_fun(fit, scenario) returns a list with at least p_value (numeric scalar)

This allows mp_power() to run before selecting a specific backend (e.g., lme4).

Value

An object of class mp_scenario.

Examples

d <- mp_design(list(subject = 20), trials_per_cell = 5)
a <- mp_assumptions(list(condition = 0.3), residual_sd = 1)
s <- mp_scenario(y ~ condition, d, a, test = "wald")
s

Gaussian LMM scenario using glmmTMB

Description

Same data-generating process as mp_scenario_lme4() but fits with glmmTMB (glmmTMB()).

Usage

mp_scenario_glmmtmb_lmm(
  formula,
  design,
  assumptions,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_term = predictor,
  test_method = c("wald", "lrt"),
  null_formula = NULL
)

Arguments

Create a fully specified MixPower scenario with the lme4 backend

Description

Create a fully specified MixPower scenario with the lme4 backend

Usage

mp_scenario_lme4(
  formula,
  design,
  assumptions,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_term = predictor,
  test_method = c("wald", "lrt", "satterthwaite", "kenward-roger", "pb"),
  null_formula = NULL,
  pb_nsim = 100L,
  contrast = NULL
)

Arguments

Value

An object of class mp_scenario.

Create a fully specified MixPower scenario with the binomial lme4 backend

Description

Create a fully specified MixPower scenario with the binomial lme4 backend

Usage

mp_scenario_lme4_binomial(
  formula,
  design,
  assumptions,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_term = predictor,
  test_method = c("wald", "lrt", "pb"),
  null_formula = NULL,
  pb_nsim = 100L
)

Arguments

Value

An object of class mp_scenario.

Create a fully specified MixPower scenario with the NB lme4 backend

Description

Create a fully specified MixPower scenario with the NB lme4 backend

Usage

mp_scenario_lme4_nb(
  formula,
  design,
  assumptions,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_term = predictor,
  test_method = c("wald", "lrt", "pb"),
  null_formula = NULL,
  pb_nsim = 100L
)

Arguments

Value

An object of class mp_scenario.

Create a fully specified MixPower scenario with the Poisson lme4 backend

Description

Create a fully specified MixPower scenario with the Poisson lme4 backend

Usage

mp_scenario_lme4_poisson(
  formula,
  design,
  assumptions,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  test_term = predictor,
  test_method = c("wald", "lrt", "pb"),
  null_formula = NULL,
  pb_nsim = 100L
)

Arguments

Value

An object of class mp_scenario.

Run power sensitivity analysis over a parameter grid

Description

Run power sensitivity analysis over a parameter grid

Usage

mp_sensitivity(
  scenario,
  vary,
  nsim = 100,
  alpha = 0.05,
  seed = NULL,
  failure_policy = c("count_as_nondetect", "exclude"),
  conf_level = 0.95
)

Arguments

Value

An object of class mp_sensitivity.

Run one sensitivity grid cell (for serial and parallel workers)

Description

Run one sensitivity grid cell (for serial and parallel workers)

Usage

mp_sensitivity_cell_run(
  i,
  scenario,
  vary,
  grid,
  nsim,
  alpha,
  seed,
  failure_policy,
  conf_level,
  checkpoint_dir,
  resume
)

Parallel sensitivity analysis over a parameter grid

Description

Like mp_sensitivity(), but evaluates each grid cell in parallel (or with a progress bar when progress = TRUE). Uses per-cell seeds seed + cell_index - 1L to match a serial ordering convention. Does not modify mp_power().

Usage

mp_sensitivity_parallel(
  scenario,
  vary,
  nsim = 100,
  alpha = 0.05,
  seed = NULL,
  failure_policy = c("count_as_nondetect", "exclude"),
  conf_level = 0.95,
  workers = 2L,
  progress = FALSE,
  checkpoint_dir = NULL,
  resume = TRUE,
  ...
)

Arguments

Value

An object of class mp_sensitivity (same structure as mp_sensitivity()).

Note

Parallel execution requires the parallel package and that mixpower can be loaded on workers (installed package).

Set a smallest effect size of interest (SESOI) on a scenario

Description

Convenience helper for planning power around a smallest effect size of interest rather than a single point estimate. It returns a copy of scenario with the focal fixed effect replaced, either by an explicit value or by scaling the current assumed effect (e.g. a conservative 15% reduction, multiplier = 0.85). This is the recommended way to avoid overoptimistic power based on a possibly inflated pilot/published effect (Anderson, Kelley & Maxwell, 2017; Kumle, Vo & Draschkow, 2021).

Usage

mp_sesoi(scenario, multiplier = 0.85, effect = NULL, term = NULL)

Arguments

Value

A modified mp_scenario object.

See Also

mp_safeguard_effect() for a data-driven conservative effect.

Examples

d <- mp_design(list(subject = 30), trials_per_cell = 8)
a <- mp_assumptions(
  fixed_effects = list("(Intercept)" = 0, condition = 0.5),
  random_effects = list(subject = list(intercept_sd = 0.5)),
  residual_sd = 1
)
scn <- mp_scenario_lme4(y ~ condition + (1 | subject), design = d, assumptions = a)
# Power for an effect 15% smaller than assumed:
scn_sesoi <- mp_sesoi(scn, multiplier = 0.85)
scn_sesoi$assumptions$fixed_effects$condition

Solve for minimum sample size achieving target power

Description

Evaluates power on a user-supplied grid of values for one parameter (e.g. cluster size) via mp_power_curve(), then returns the smallest grid value whose power estimate meets or exceeds the target. Diagnostics (failure rate, singular rate, n_effective) are exposed in the returned results table.

Usage

mp_solve_sample_size(
  scenario,
  parameter,
  grid,
  target_power = 0.8,
  nsim = 100,
  alpha = 0.05,
  seed = NULL,
  failure_policy = c("count_as_nondetect", "exclude"),
  conf_level = 0.95
)

Arguments

Value

A list with target_power, parameter, solution (numeric: minimum grid value achieving target power, or NA if none), and results (data frame with estimate, failure_rate, singular_rate, n_effective, etc., per grid point).

Write results or bundle to CSV or JSON

Description

Writes the report table (and for bundles, manifest/labels) to file. CSV writes the publication-ready table only; JSON writes report table plus manifest and labels when x is an mp_bundle.

Usage

mp_write_results(x, file, format = c("csv", "json"), ...)

Arguments

Value

Invisibly the path file.

Plot the p-value distribution of a power analysis

Description

A histogram of the per-replicate p-values from mp_power(), with the alpha threshold marked. The shaded area to the left of alpha is the estimated power. Requires a run with aggregate = "full" (the default), which retains per-replicate p-values.

Usage

## S3 method for class 'mp_power'
plot(x, ...)

Arguments

Value

Invisibly, the p-values plotted.

Plot a power curve

Description

Plot a power curve

Usage

## S3 method for class 'mp_power_curve'
plot(x, y = c("estimate", "failure_rate", "singular_rate", "n_effective"), ...)

Arguments

Value

Invisibly returns the plotted data.

Plot a sensitivity analysis

Description

For one varying parameter: line plot with optional CI segments when y = "estimate". For two varying parameters: heatmap. More than two parameters is not supported.

Usage

## S3 method for class 'mp_sensitivity'
plot(x, y = c("estimate", "failure_rate", "singular_rate", "n_effective"), ...)

Arguments

Value

Invisibly returns the plotted data (1D: ordered data frame; 2D: matrix).

Plot power results

Description

Plot power results

Usage

plot_power(results, ...)

Arguments

Value

Invisibly returns the plot data.

Placeholder for parallel execution

Description

Placeholder for parallel execution

Usage

run_parallel(fun, ...)

Arguments

Value

The result of fun.

Simulate binary outcome data for a GLMM with random effects

Description

Thin wrapper over the shared simulation engine (.mp_simulate_mixed) with a logit link and Bernoulli response. Random-effect sizes (intercept and optional slope on predictor) come from scenario$assumptions$random_effects.

Usage

simulate_glmm_binomial_data(
  scenario,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL
)

Arguments

Value

A data.frame with outcome and predictors.

Simulate count outcome data for a Negative Binomial GLMM with random effects

Description

Thin wrapper over the shared simulation engine (.mp_simulate_mixed) with a log link and negative-binomial response. Random-effect sizes (intercept and optional slope on predictor) come from scenario$assumptions$random_effects.

Usage

simulate_glmm_nb_data(
  scenario,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL,
  theta = NULL
)

Arguments

Value

A data.frame with outcome and predictors.

Simulate count outcome data for a Poisson GLMM with random effects

Description

Thin wrapper over the shared simulation engine (.mp_simulate_mixed) with a log link and Poisson response. Random-effect sizes (intercept and optional slope on predictor) come from scenario$assumptions$random_effects.

Usage

simulate_glmm_poisson_data(
  scenario,
  predictor = "condition",
  subject = "subject",
  outcome = "y",
  item = NULL
)

Arguments

Value

A data.frame with outcome and predictors.

Run a simple simulation-based power study

Description

Run a simple simulation-based power study

Usage

simulate_power(scenario, nsim = 100, seed = NULL)

Arguments

Value

A data.frame of simulated p-values.

Summarize simulation outputs

Description

Summarize simulation outputs

Usage

summarize_simulations(simulations)

Arguments

Value

A summary data.frame.

Extract a test statistic for a model term

Description

Extract a test statistic for a model term

Usage

test_effect(fit, term)

Arguments

Value

A data.frame with coefficient information.

Validate a MixPower backend

Description

Checks that simulate_fun, fit_fun, and test_fun are functions and that their formal arguments are compatible with what mp_power() and .run_one_rep() invoke. Does not run simulations.

Usage

validate_mp_backend(engine)

Arguments

Value

Invisibly TRUE if valid; otherwise throws an error.