Skip to content

Power and design analysis

pilotr turns a specification into evidence for study planning: simulate from the ground truth, fit the analysis model and summarise across replicates. Alongside power it reports the Type S (sign) and Type M (magnitude, or exaggeration) errors of Gelman and Carlin (2014), computed over the significant replicates.

Two-group Gaussian power

power handles the two-group Gaussian design with a two-sample t-test. It needs scipy (install the power extra):

from pilotr import power

spec = {
    "name": "d", "seed": 1,
    "units": {"subject": {"n": 64}},
    "factors": [{"name": "group", "levels": ["a", "b"],
                 "contrasts": {"effect": [-0.5, 0.5]}, "between": "subject"}],
    "fixed": {"intercept": 100, "coefficients": {"effect": 5}},
    "response": {"family": "gaussian", "name": "score", "sigma": 10},
}

res = power(spec, n_sims=500)
print(table([{k: res[k] for k in (
    "power", "type_s", "type_m", "true_effect", "mean_estimate"
)}]))
power type_s type_m true_effect mean_estimate
0.504 0 1.4 5 4.99

At roughly 50% power the Type M ratio is well above 1. Conditional on significance the estimated effect is exaggerated, even though the average estimate over all replicates is unbiased. This is the statistical-significance filter that design analysis is meant to expose.

Power over sample size

power_curve sweeps the number of subjects and reports power at each:

import matplotlib.pyplot as plt
from math import sqrt
from pilotr import power_curve, target_n

n_sims = 200
curve = power_curve(spec, subject_ns=[16, 32, 48, 64, 96, 128, 160, 192], n_sims=n_sims)
ns = [p["n_subject"] for p in curve]
pw = [p["power"] for p in curve]
# Each power estimate is a proportion over n_sims replicates, so it carries a
# binomial Monte Carlo standard error. The shaded band is the 95% interval.
se = [sqrt(p * (1 - p) / n_sims) for p in pw]
lo = [max(0, p - 1.96 * s) for p, s in zip(pw, se)]
hi = [min(1, p + 1.96 * s) for p, s in zip(pw, se)]

# The sample size the curve solves to, drawn as a vertical band. The next
# section explains the call.
solved = target_n(curve, target=0.8)

fig, ax = plt.subplots(figsize=(6, 3.4))
ax.axhline(0.8, ls="--", color="grey")
ax.axvspan(solved["lo"], solved["hi"], color="grey", alpha=0.12)
ax.axvline(solved["value"], ls="--", color="grey")
ax.fill_between(ns, lo, hi, color=BLUE, alpha=0.15)
ax.plot(ns, pw, "-o", color=BLUE)
ax.set_ylim(0, 1)
ax.set_xlabel("$N$ subjects")
ax.set_ylabel("Power")
print(show(fig))
2026-08-23T12:11:18.750806 image/svg+xml Matplotlib v3.11.1, https://matplotlib.org/

The shaded horizontal band along the curve is the Monte Carlo interval, the binomial standard error of each power estimate over the n_sims replicates widened to a 95% envelope. The vertical band is the solved sample size, with its own interval.

The sample size the curve implies

The point of a power curve is the sample size at which power reaches the target, and that is the number a preregistration quotes. Judging the crossing by eye gives a bare figure over points whose Monte Carlo intervals overlap. target_n fits the curve and inverts the fit, so the crossing arrives with the uncertainty that a simulated curve carries.

solved = target_n(curve, target=0.8)
print(table([{"target": solved["target"], "n": solved["n"],
              "n_lo": solved["n_lo"], "n_hi": solved["n_hi"]}]))
target n n_lo n_hi
0.8 133 123 142

This design has a closed-form answer to check against. With an effect of 5 and a residual standard deviation of 10, power.t.test in R puts 0.80 power at 127.5 subjects in total, which falls inside the interval above. The interval spans some twenty subjects, the honest report at 200 replicates a point, and raising n_sims narrows it.

Nothing is extrapolated. A curve that never reaches the target within the sizes swept is refused, and the refusal reports the range the sweep did cover. Had the sweep above stopped at 128 subjects, where power was still short of 0.80, that is what would have happened.

solve_curve is the general form. It takes any curve with a decision rate against a swept value, so an effect-size sweep or a sweep over items is solved the same way, with transform="identity" for an effect-size axis.

Crossed mixed-effects power

power_mixed fits a crossed mixed model with statsmodels (needs statsmodels and pandas). The design has one within-subject factor and crossed by-subject and by-item random intercepts and slopes:

from pilotr import power_mixed

spec_mixed = {
    "name": "priming", "seed": 1,
    "units": {"subject": {"n": 12}, "item": {"n": 8}},
    "factors": [{"name": "condition", "levels": ["related", "unrelated"],
                 "contrasts": {"cond": [-0.5, 0.5]}, "vary_within": "subject"}],
    "fixed": {"intercept": 6, "coefficients": {"cond": 0.1}},
    "random": {
        "subject": {"intercept_sd": 0.12, "slopes": {"cond": 0.04},
                    "correlations": {"intercept~cond": 0.2}},
        "item": {"intercept_sd": 0.08, "slopes": {"cond": 0.02},
                 "correlations": {"intercept~cond": -0.1}},
    },
    "response": {"family": "shifted_lognormal", "name": "RT",
                 "sigma": 0.3, "shift": 200},
}

# A tiny replicate count keeps the docs build fast. Use 200 or more for real planning.
res = power_mixed(spec_mixed, n_sims=12)
print(table([{k: res[k] for k in (
    "power", "n_converged", "true_effect", "mean_estimate", "type_s", "type_m"
)}]))
power n_converged true_effect mean_estimate type_s type_m
0.5 12 0.1 0.1 0 1.33

n_converged reports how many replicates the mixed model actually fit. Convergence problems are common in small crossed designs, so this is a useful diagnostic in its own right, and it is the denominator of power, the significant proportion among the converged replicates.

Even at this tiny n_sims, the fixed effect is recovered (mean_estimate is close to true_effect). The statsmodels variance-component fit overstates random-slope variance, so the power estimate is conservative for random-slope designs, and the R package's lme4-based power_mixed is the reference for them. Data generation is identical across the two languages, and the discrepancy is in the estimator alone.

power_mixed carries its own simulation loop over the portable specification, with no other power package underneath it. It covers territory pioneered by simr (Green and MacLeod, 2016) and mixedpower (Kumle, Vo and Draschkow, 2021). pilotr differs in being driven by the portable cross-language specification, in reporting Type S and Type M errors alongside power and in built-in parallelisation.

Parallel execution

Every analysis on this page takes a workers argument that spreads the Monte Carlo replicates across local processes with concurrent.futures. Each replicate takes its own seed from replicate_seeds(), so the results are identical to a serial run whatever the worker count, and parallelisation costs nothing in reproducibility. The model fits dominate the running time, which makes the speed-up close to linear in the number of processes. power_curve starts one process pool and reuses it across the whole sweep.

from pilotr import power, power_curve, power_mixed

power(spec, n_sims=2000, workers=8)          # same numbers as workers=1, sooner
power_curve(spec, subject_ns=[16, 32, 48, 64, 96, 128], n_sims=2000, workers=8)
power_mixed(spec_mixed, n_sims=500, workers=8)

When calling a parallel analysis from a script on Windows or macOS, put the call inside an if __name__ == "__main__": block, the standard requirement for Python's spawn-based multiprocessing. Interactive sessions and notebooks need no guard.

This design answers a serial bottleneck familiar from simr::powerCurve() in R, which pilotr's maintainer previously worked around by splitting the sample-size grid across separate jobs by hand and recombining the results afterwards (Bernabeu, 2021). In pilotr the same gain takes one argument.