MCPower LME Validation Process

1. What is MCPower?

MCPower is a Monte Carlo power analysis library for complex statistical models. It estimates the probability (power) that a study will detect a true effect at a given sample size, using simulation rather than closed-form formulas.

Its LME (Linear Mixed-Effects) module handles clustered/hierarchical data — students nested in schools, patients nested in hospitals, repeated measures within subjects. Instead of relying on R or Python's statsmodels, MCPower uses a custom C++ solver based on the profiled-deviance approach from Bates et al. (2015), achieving roughly 200× speedup over pure Python.

Statistical inference uses Likelihood Ratio (LR) tests for overall significance and Wald z-tests (with REML) for individual fixed effects — the same testing framework that lmerTest provides in R.

2. Why validate against lme4?

R's lme4 package (Bates, Mächler, Bolker & Walker, 2015) is the gold standard for fitting linear mixed-effects models. It is used in tens of thousands of published studies across psychology, education, ecology, and medicine. By validating MCPower's solver against lme4, we establish that our custom implementation produces statistically equivalent results.

3. Validation framework overview

The framework uses four independent, complementary strategies that together test the solver, the data-generating process (DGP), and end-to-end power estimation:

4. Strategy 1 — External Data Agreement

This strategy isolates the solver from everything else. Data is generated by a third party (NumPy) so neither MCPower nor R influences the data-generating process. Both solvers receive identical datasets and we compare their conclusions.

# Fit REML model (lmerTest provides Satterthwaite p-values)
fit <- lmer(as.formula(formula_str), data = d, REML = TRUE,
            control = lmerControl(optimizer = "bobyqa",
                                    optCtrl = list(maxfun = 50000)))

# Extract fixed effects (skip intercept)
coefs <- summary(fit)$coefficients
beta    <- coefs[-1, "Estimate"]
se      <- coefs[-1, "Std. Error"]
pvalues <- coefs[-1, "Pr(>|t|)"]

# Variance components
vc    <- as.data.frame(VarCorr(fit))
tau2  <- vc$vcov[vc$grp == "g"]
sigma2 <- sigma(fit)^2

# Likelihood Ratio Test: full vs null (intercept-only + RE)
null_fit <- lmer(y ~ 1 + (1 | g), data = d, REML = FALSE,
                 control = lmerControl(optimizer = "bobyqa",
                                         optCtrl = list(maxfun = 50000)))
full_ml <- update(fit, REML = FALSE)
lr_test <- anova(null_fit, full_ml)
lr_pvalue <- lr_test$`Pr(>Chisq)`[2]
significant <- !is.na(lr_pvalue) && lr_pvalue < 0.05

5. Strategy 2 — MCPower DGP Validation

Strategy 1 tests the solver in isolation using third-party data. Strategy 2 adds MCPower's data-generating process to the test. Data is generated through MCPower's full pipeline (the same code path used during power analysis), then both MCPower and R fit models on these datasets.

6. Strategy 3 — Parallel Power Simulation

The ultimate test: do MCPower and R produce the same power estimates when each generates data and fits models completely independently?

Both systems use the same statistical model specification (ICC, effect sizes, sample size, number of clusters) but use independent RNG streams and independent DGP implementations.

# Compute variance components from ICC (same formula as MCPower)
sigma2_fixed <- sum(effect_values^2)
sigma2_within <- 1.0 + sigma2_fixed
tau2 <- icc / (1 - icc) * sigma2_within

for (sim in 1:n_sims) {
  # Generate cluster assignments
  g <- rep(1:n_clusters, each = cluster_size)

  # Generate predictors ~ N(0, 1)
  x1 <- rnorm(n_obs)

  # Random intercepts ~ N(0, tau2)
  b <- rnorm(n_clusters, 0, sqrt(tau2))

  # Response: y = X*beta + b[group] + error
  y <- beta * x1 + b[g] + rnorm(n_obs)

  # Fit & test via LRT
  fit  <- lmer(y ~ x1 + (1 | g), REML = FALSE)
  null <- lmer(y ~  1 + (1 | g), REML = FALSE)
  lr_p <- anova(null, fit)$`Pr(>Chisq)`[2]

  if (!is.na(lr_p) && lr_p < 0.05) n_significant += 1
}
power <- n_significant / n_converged

7. Strategy 4 — Statistical Power Comparison (Z-test + FDR)

Strategy 3 compares MCPower and R power estimates using a fixed threshold (|diff| ≤ 5 percentage points). With 500 Monte Carlo simulations, sampling noise alone can cause ~2–3pp deviations, leading to occasional false failures in sensitivity scenarios near the threshold boundary. Strategy 4 replaces this fixed threshold with a proper statistical test.

Two-Proportion Z-test

For each scenario, we test H₀: Power_MCPower = Power_R using the standard two-proportion z-test:

where the pooled standard error under H₀ is:

and p̂_pool = (p̂_MC · n_MC + p̂_R · n_R) / (n_MC + n_R) is the pooled proportion. The two-sided p-value is computed via the complementary error function: p = erfc(|z| / √2).

Edge case: When both tools agree on extreme power (0% or 100%), the pooled SE is zero. In this case z = 0, p = 1.0, and the scenario trivially passes.

Benjamini-Hochberg FDR Correction

With ~50 scenarios tested simultaneously, multiple testing inflation would produce false positives. The BH procedure controls the False Discovery Rate at level α = 0.05:

Interpretation

8. Scenario Matrix

The validation uses two scenario categories: ~27 core scenarios (tested with all four strategies S1–S4) and ~50 sensitivity scenarios (tested with S4 only). Core scenarios verify solver correctness with well-powered designs. Sensitivity scenarios span the full 5%–90% power range to validate power curve accuracy.

Random Intercept Scenarios (~20)

Random Slopes Scenarios (~4)

Nested Random Effects Scenarios (~3)

Sensitivity Scenarios (~50, S4 only)

9. Pass/Fail Thresholds

Each scenario gets a PASS or FAIL verdict based on the following thresholds:

10. R Implementation Details

For R users who want to understand exactly what lme4 configuration is being compared against.

R Packages

Model Fitting Configuration

• Estimation: REML for parameter extraction, ML refit for Likelihood Ratio Tests
• Optimizer: bobyqa (bound optimization by quadratic approximation)
• Max iterations: 50,000 function evaluations via optCtrl = list(maxfun = 50000)
• Convergence failures: caught with tryCatch(), marked as unconverged, excluded from comparisons

Statistical Tests

# Null model: intercept + random effects only
null_fit <- lmer(y ~ 1 + (1 | g), data = d, REML = FALSE)
full_ml  <- update(fit, REML = FALSE)
lr_test  <- anova(null_fit, full_ml)
# Chi-squared test with df = number of fixed effects

R↔Python Communication

R scripts are invoked via subprocess.run(["Rscript", "--vanilla", script_path, ...]) with a 600-second timeout. Data passes as CSV files; results return as JSON on stdout. Stderr is reserved for R warnings/messages and discarded unless an error occurs.

A smart caching layer hashes the scenario parameters + R script contents (SHA-256) so unchanged scenarios skip re-running R, which is critical since a full validation run with all strategies can take 30+ minutes.

11. Running the Validation

Prerequisites

• Python 3.10+ with MCPower installed (pip install -e MCPower/)
• R with packages: lme4, lmerTest, jsonlite
• Rscript available on $PATH

Commands

# Run all strategies on all scenarios (full validation)
python LME4-validation/run_validation.py

# Run only Strategy 1 (solver comparison)
python LME4-validation/run_validation.py --strategy 1

# Run specific scenarios by name prefix
python LME4-validation/run_validation.py --scenarios intercept_icc0.2

# Run with more datasets (default: 100) or simulations (default: 500)
python LME4-validation/run_validation.py --n-datasets 200 --n-sims 1000

# Clear cached R results and re-run from scratch
python LME4-validation/run_validation.py --clear-cache

# Disable caching entirely
python LME4-validation/run_validation.py --no-cache

# Set custom random seed (default: 42)
python LME4-validation/run_validation.py --seed 123

Example output

Running validation: 27 scenarios, strategies [1, 2, 3, 4]
  n_datasets=100, n_sims_s3=500, seed=42

============================================================
Strategy 1: External Data Agreement
============================================================
  [S1] intercept_icc0.2_cl20_n1000_small ... PASS (agree=100.0%, 12.4s)
  [S1] intercept_icc0.2_cl20_n1000_medium ... PASS (agree=100.0%, 11.8s)
  ...

============================================================
Strategy 3: Parallel Power Simulation
============================================================
  [S3] intercept_icc0.2_cl20_n1000_medium ... PASS (MC=0.847, R=0.832, 45.2s)
  ...

============================================================
Strategy 4: Statistical Power Comparison (Z-test + FDR)
============================================================
  [S4] intercept_icc0.2_cl20_n1000_medium: PASS (z=0.672, adj_p=0.8921)
  ...

============================================================
DONE: 27/27 scenarios passed (892.3s total)
Report: LME4-validation/reports/validation_report_20260222_143201.html
============================================================

12. Reading the Reports

After each run, a self-contained HTML report is generated in LME4-validation/reports/. The report includes:

Reports are fully self-contained (inline CSS, no external dependencies) and can be shared directly. Each cell is color-coded green (pass) or red (fail) for quick scanning.