In a probit or logit model the coefficient is not the quantity anyone
interprets. Conclusions are stated in probability points: the treatment
raises take-up by ten percentage points. That quantity is the
average partial effect (APE), and it depends on all
coefficients and on the covariate distribution, so power tools that work
in coefficient units cannot say how many observations are needed to
detect it. powerape states the effect in APE units and
simulates the whole study: it draws covariates and outcomes, fits the
model, computes the APE with a delta-method confidence interval, and
checks a claim on that interval (Riesthuis, 2024):
-
minimum: the lower bound of the interval exceeds the smallest effect size of interest (SESOI), so the effect is meaningfully large. This is the default claim. -
detect: the interval excludes zero. -
equivalence: the interval lies inside (−SESOI, +SESOI), so the effect is too small to matter.
The Stata command powerape is a port of the R package of
the same name and is documented side by side with it on this site. This
page takes the six steps of the R article Powering the minimum-effect claim through
Stata, in the same world, and then shows what is specific to Stata: who
fits the simulated studies, and how to start from pilot data.
Installation
net install copies the command, its help file and its
Mata library; mata: mata mlib index makes the new library
known to the running Stata session (needed once, right after
installing); powerape version checks that the command and
the library belong together. To update, add replace to
net install. net get powerape, from(...) with
the same address copies two files into the current folder:
showcase.do, a guided tour of every subcommand, and
build.do, which recompiles the Mata library.
Stata 16 or later is required. The control-function route for
endogenous treatments uses cfprobit (Stata 18 or later),
and computing interaction effects with ginteff itself, an
option, needs that command from the Stata Journal
(net install st0711, from(http://www.stata-journal.com/software/sj23-2)).
The command has been tested on StataNow 19.5. After installation,
help powerape opens the help file, which is also on this site.
Step 1: describe the world
One command describes the data-generating process (DGP) the study
will sample from and stores it under a name, the way
estimates store stores results. Here: a probit model, a
binary treatment given to half the sample, two covariates, an outcome
rate of 0.30 without treatment, and covariates that explain 15% of the
latent variation.
powerape dgp d, focal(treat, binary p(0.5)) ///
covariates((age, normal mean(45) sd(12)) (female, binary p(0.55))) ///
correlation(0.2) baseline(0.30) signal(0.15)powerape DGP -- probit
focal: treat (binary, p = 0.50)
covariates: age (normal, mean 45.00, sd 12.00), female (binary, p = 0.55)
baseline P(Y=1 | reference) = 0.300, nuisance signal (latent pseudo-R2) = 0.150
target effect not set yet - run powerape setape or powerape setaie.Each variable is written as (name, type options):
binary p() or normal mean() sd(). The nuisance
side is set on two friendly scales, the baseline outcome
rate (baseline()) and the covariates’
latent pseudo-R² (signal()), and
powerape solves for the coefficients that produce them.
correlation() sets the latent correlation among the focal
variable and the covariates. A pilot data set or a fitted pilot model
can take the place of these guesses; see Pilot
data below.
Step 2: pin the true APE
Suppose prior evidence suggests a true effect of ten percentage points, and effects below five points would not justify the program: the planning value is 0.10 and the SESOI is 0.05.
powerape DGP -- probit
focal: treat (binary, p = 0.50)
covariates: age (normal, mean 45.00, sd 12.00), female (binary, p = 0.55)
baseline P(Y=1 | reference) = 0.300, nuisance signal (latent pseudo-R2) = 0.150
true APE pinned at + 0.1000 (beta_focal = 0.2943, implied P(Y=1 | focal=1) = 0.400)powerape setape solves for the coefficient of the
treatment that makes the true APE exactly 0.10, given everything else,
and refuses with an explanation when the target cannot exist (an APE of
+0.40 at a baseline rate of 0.85, say). For interaction designs
powerape setaie pins the average interaction effect
instead; see Interaction
effects.
Step 3: power for the claim
powerape -- minimum-effect claim (CI lower bound > 0.050)
probit, n = 1500, assumed true APE + 0.1000, 90% CI (alpha = .05), nsim = 500
fit: probit via engine(stata), effects(closed)
power = 0.652 (MCSE 0.021)
outcomes: minimum 0.652 | detect-only 0.338 | inconclusive 0.010 | equivalence 0.000 | failed 0.000In each of the 500 replications powerape drew a sample
of 1,500 observations, fitted it, and asked whether the 90 percent
confidence interval of the APE lies above the SESOI (a one-sided test at
the 5 percent level). It did in 65.2 percent of them. The line beginning
fit: says who fitted the simulated studies: Stata’s own
probit, with the APE and its standard error computed from
e(b) and e(V) exactly as margins
computes them.
Read the outcome distribution, not only the power. About a third of the studies would show an effect different from zero but could not claim that it clears the SESOI (“detect-only”): exactly the inconclusive result that a power analysis for detection alone hides. The R article finds the same split in the same world, up to simulation error.
Stored results
Every subcommand stores its results in r().
scalars:
r(is_margin) = 0
r(mde) = .
r(nsim_confirm) = .
r(confirmed) = 0
r(goal) = .
r(n_failed) = 0
r(min_cell) = 206.746
r(separated) = 0
r(target) = .1
r(sesoi) = .05
r(conf) = .9
r(alpha) = .05
r(nsim) = 500
r(n) = 1500
r(mcse) = .0213023942316351
r(power) = .652
macros:
r(fit) : "probit via engine(stata), effects(closed)"
r(effects) : "closed"
r(engine) : "stata"
r(dgp) : "d"
r(separation) : "fail"
r(se) : "model"
r(model) : "probit"
r(estimand) : "ape"
r(claim) : "minimum"
r(kind) : "power"
r(name) : "_last"
matrices:
r(outcomes) : 1 x 5Results also stay available by name: powerape list shows
the stored DGPs and results, and the option name() keeps a
result under a name of your choice (see Justified alpha for a use).
Step 4: required sample size and the power curve
powerape n searches for the sample size that reaches the
target power, then confirms it with a high-precision run and moves it
upward, never down, if the confirmation falls short.
powerape required sample size -- minimum claim
n = 2103 for 80% target power (confirmed 0.808, MCSE 0.009)
assumed true APE + 0.1000, sesoi 0.050, 90% CI (alpha = .05), probit
fit: Mata engine (Fisher scoring)
search: 1 step(s); confirmed in 1 round(s) at nsim = 2000.A search runs thousands of simulated studies (here the confirmation
stage alone ran 2,000), so this call uses the Mata engine,
engine(mata), a self-contained port of the R package’s
fitting routine that is many times faster than calling
probit. Engines below
explains the choice. The search took 15 seconds.
A power curve evaluates a list of sample sizes on the same footing.
powerape power curve -- minimum claim, probit (APE), nsim = 300 per point
fit: Mata engine (Fisher scoring)
n alpha power mcse failed separated
500 0.0500 0.3267 0.0271 0.0000 0.0000
1000 0.0500 0.5267 0.0288 0.0000 0.0000
1500 0.0500 0.6667 0.0272 0.0000 0.0000
2000 0.0500 0.7567 0.0248 0.0000 0.0000
2500 0.0500 0.8700 0.0194 0.0000 0.0000
3000 0.0500 0.8800 0.0188 0.0000 0.0000r(results) holds the curve as a matrix, ready to
plot:
matrix C = r(results)
svmat C, names(col)
twoway connected power n, yline(0.8, lpattern(dash)) ///
ytitle("Power") xtitle("Sample size") ylabel(0(0.2)1)
When n is fixed: the minimum detectable effect
Often the sample is given: the archive has the firms it has.
powerape mde turns the question around and finds the
smallest true APE that the claim reaches with 80 percent power, with the
same confirmation stage as powerape n.
powerape minimum detectable APE -- minimum claim
MDE = 0.1094 at n = 1500 for 80% target power (confirmed 0.797, MCSE 0.009)
smallest effect demonstrably above sesoi 0.050
fit: Mata engine (Fisher scoring)
search: 1 step(s); confirmed in 1 round(s) at nsim = 2000.With 1,500 observations, the study can show an effect to be meaningfully above 0.05 if the true effect is about 0.109 or larger.
Step 5: how fragile is this to the contextual guesses?
The baseline rate and the covariate signal are guesses.
powerape robust repeats the analysis over a grid of them
and reports the worst case (Hancock & Feng, 2026). By default the
APE is re-pinned in every scenario (pin(ape)), so only the
precision of the estimate varies.
powerape robust d, n(2200) claim(minimum) sesoi(0.05) ///
vary(baseline(0.20 0.40)) gridpoints(3) nsim(300) seed(4) nonmax engine(mata)powerape robustness sweep -- minimum claim, APE, pin = ape
n = 2200, nsim = 300 per scenario, 3 scenario(s) over: baseline
power range [0.793, 0.870]; worst scenario: row 3
fit: Mata engine (Fisher scoring)
row baseline implied_eff power mcse separated
1 .2 0.1000 0.8700 0.0194 0.0000
2 .3 0.1000 0.8100 0.0226 0.0000
3 .4 0.1000 0.7933 0.0234 0.0000Dropping nonmax adds n_max, the largest
required sample size across the scenarios, at 90 percent power by
default (nmaxpower()): the sample size that insures the
design against the guesses. mode(mde) reports the minimum
detectable effect per scenario instead of power.
Step 6: the citable paragraph
powerape statement writes the last result as a methods
paragraph for a grant application or a preregistration.
powerape -- minimum-effect claim (CI lower bound > 0.050)
probit, n = 2200, assumed true APE + 0.1000, 90% CI (alpha = .05), nsim = 500
fit: probit via engine(stata), effects(closed)
power = 0.816 (MCSE 0.017)
outcomes: minimum 0.816 | detect-only 0.180 | inconclusive 0.004 | equivalence 0.000 | failed 0.000
We conducted a simulation-based power analysis for the average partial effect
(APE) of treat using the powerape package for Stata (version 0.1.3 (Stata port
of powerape 1.11.3)), following the confidence-interval approach of Riesthuis
(2024). The assumed data-generating process was a probit model with focal
variable treat (binary, prevalence 0.50); parametric covariates age
(continuous, mean 45, SD 12), female (binary, prevalence 0.55), with an
exchangeable latent (Gaussian-copula) correlation of 0.20 among the focal
variable and the covariates; baseline outcome rate 0.300 with the focal at
reference; nuisance covariates contribute a latent pseudo-R-squared of 0.15.
The assumed true APE was 0.100 (10.0 percentage points). At a sample size of n
= 2200, simulated power for the minimum-effect claim (the 90% confidence
interval's lower bound exceeding the smallest effect size of interest, 0.050;
a one-sided test at alpha = 5%) was 0.816 (Monte Carlo SE 0.017; 500
replications). Across replications, the probability of concluding a meaningful
effect was 0.816, of detection without meaningfulness 0.180, of an
inconclusive result 0.004, and of equivalence 0.000. Estimation used maximum
likelihood with delta-method Wald confidence intervals; replications without a
maximum-likelihood estimate, from a failed fit or from separation (a focal
cell with no events or no non-events), counted against the claim (0.0% of
replications).Engines: who fits the simulated studies
The simulation draws the data in Mata. Fitting each simulated study can be left to Stata or done in Mata:
-
engine(stata), the default, writes every simulated study to a temporary frame and fits it with Stata’s ownprobitorlogit, withvce(cluster)for panels andcfprobitfor endogenous treatments. The effects are computed frome(b)ande(V)by the closed-form delta method thatmarginsuses. -
engine(mata)fits each study with a Mata port of the R package’s routine (Fisher scoring, expected information).
On the same seed both engines draw the same studies, so they can be compared directly:
powerape -- minimum-effect claim (CI lower bound > 0.050)
probit, n = 1500, assumed true APE + 0.1000, 90% CI (alpha = .05), nsim = 500
fit: Mata engine (Fisher scoring)
power = 0.652 (MCSE 0.021)
outcomes: minimum 0.652 | detect-only 0.340 | inconclusive 0.008 | equivalence 0.000 | failed 0.000The power is the same, 0.652. The outcome split differs by one study
of 500, because probit computes standard errors from the
observed information and the Mata engine from the expected information,
which moves an interval that sits right at a boundary. The Mata engine
took 1.9 seconds, probit 23.8. Use the Mata engine for
searches, curves and sweeps, and confirm the final design with Stata’s
estimator.
effects(margins) goes one step further and computes
every APE with margins itself (or with ginteff
for interaction effects). On the same seed it reproduces the closed-form
results exactly:
powerape power d, n(1500) claim(minimum) sesoi(0.05) nsim(100) seed(7)
powerape power d, n(1500) claim(minimum) sesoi(0.05) nsim(100) seed(7) ///
effects(margins)powerape -- minimum-effect claim (CI lower bound > 0.050)
probit, n = 1500, assumed true APE + 0.1000, 90% CI (alpha = .05), nsim = 100
fit: probit via engine(stata), effects(closed)
power = 0.720 (MCSE 0.045)
outcomes: minimum 0.720 | detect-only 0.260 | inconclusive 0.020 | equivalence 0.000 | failed 0.000
powerape -- minimum-effect claim (CI lower bound > 0.050)
probit, n = 1500, assumed true APE + 0.1000, 90% CI (alpha = .05), nsim = 100
fit: probit via engine(stata), effects(margins)
power = 0.720 (MCSE 0.045)
outcomes: minimum 0.720 | detect-only 0.260 | inconclusive 0.020 | equivalence 0.000 | failed 0.000It is also slower (14.1 against 3.6 seconds for these 100
replications), so its use is as a check. engine() and
effects() given on one command apply to that command only;
powerape engine shows and sets the defaults for the
session:
engine(stata) effects(closed)
engine(mata) effects(closed) set as the session default
engine(stata) effects(closed) set as the session defaultPilot data
When pilot data are available, the covariates can come from them instead of from parametric guesses. Here is a small simulated pilot:
set seed 20261016
set obs 500
generate treat = runiform() < .5
generate age = rnormal(45, 12)
generate female = runiform() < .55
generate y = runiform() < normal(-0.9 + 0.3*treat + 0.01*age + 0.2*female)empirical resamples whole rows of the pilot data in
memory and calibrates the outcome model to the stated baseline rate and
signal. A listed focal variable is resampled jointly with the
covariates; one specified as (name, binary p()) is drawn
independently of them instead.
powerape dgp pilot_rows, empirical focal(treat) covariates(age female) ///
baseline(0.30) signal(0.15)
powerape setape pilot_rows, target(0.10)
powerape power pilot_rows, n(1500) claim(minimum) sesoi(0.05) nsim(500) seed(11)powerape DGP -- probit (empirical covariates)
focal: treat (binary, p = 0.49) [resampled jointly with covariates]
covariates: age, female
pilot rows: 500 (resampled with replacement)
baseline P(Y=1 | reference) = 0.300, nuisance signal (latent pseudo-R2) = 0.150
target effect not set yet - run powerape setape or powerape setaie.
powerape DGP -- probit (empirical covariates)
focal: treat (binary, p = 0.49) [resampled jointly with covariates]
covariates: age, female
pilot rows: 500 (resampled with replacement)
baseline P(Y=1 | reference) = 0.300, nuisance signal (latent pseudo-R2) = 0.150
true APE pinned at + 0.1000 (beta_focal = 0.2940, implied P(Y=1 | focal=1) = 0.400)
powerape -- minimum-effect claim (CI lower bound > 0.050)
probit, n = 1500, assumed true APE + 0.1000, 90% CI (alpha = .05), nsim = 500
fit: probit via engine(stata), effects(closed)
power = 0.656 (MCSE 0.021)
outcomes: minimum 0.656 | detect-only 0.330 | inconclusive 0.014 | equivalence 0.000 | failed 0.000fromfit takes the world from a fitted
probit or logit model: the nuisance
coefficients and the pilot rows come from the fit, and
setape re-solves the focal coefficient for the planning
value. baseline() can override the pilot’s outcome
rate.
probit y treat age female, nolog
powerape dgp pilot_fit, fromfit focal(treat)
powerape setape pilot_fit, target(0.10)
powerape power pilot_fit, n(1500) claim(minimum) sesoi(0.05) nsim(500) seed(12)Probit regression Number of obs = 500
LR chi2(3) = 14.74
Prob > chi2 = 0.0021
Log likelihood = -332.77504 Pseudo R2 = 0.0217
------------------------------------------------------------------------------
y | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
treat | .3269278 .1142048 2.86 0.004 .1030904 .5507652
age | .0132259 .0051028 2.59 0.010 .0032246 .0232272
female | .0701507 .1152546 0.61 0.543 -.1557442 .2960455
_cons | -1.00213 .2559281 -3.92 0.000 -1.50374 -.5005203
------------------------------------------------------------------------------
powerape DGP -- probit (empirical covariates)
focal: treat (binary, p = 0.49) [resampled jointly with covariates]
covariates: age, female
pilot rows: 500 (resampled with replacement)
pilot model: probit y treat age female (nuisance coefficients lifted; focal coefficient re-solved)
baseline P(Y=1 | reference) = 0.358, nuisance signal (latent pseudo-R2) = 0.024
target effect not set yet - run powerape setape or powerape setaie.
powerape DGP -- probit (empirical covariates)
focal: treat (binary, p = 0.49) [resampled jointly with covariates]
covariates: age, female
pilot rows: 500 (resampled with replacement)
pilot model: probit y treat age female (nuisance coefficients lifted; focal coefficient re-solved)
baseline P(Y=1 | reference) = 0.358, nuisance signal (latent pseudo-R2) = 0.024
true APE pinned at + 0.1000 (beta_focal = 0.2614, implied P(Y=1 | focal=1) = 0.458)
powerape -- minimum-effect claim (CI lower bound > 0.050)
probit, n = 1500, assumed true APE + 0.1000, 90% CI (alpha = .05), nsim = 500
fit: probit via engine(stata), effects(closed)
power = 0.616 (MCSE 0.022)
outcomes: minimum 0.616 | detect-only 0.362 | inconclusive 0.022 | equivalence 0.000 | failed 0.000The R functions for the two routes are
ape_dgp_empirical() and ape_dgp_from_fit().
Panel data and endogenous treatments have their own routes,
panel and iv; see Panel designs and IV designs.
Where next
- The workflows take one design question each, with the R and the Stata code side by side: the classic SESOI analysis, equivalence, interaction effects, panel designs, IV designs and justified alpha.
- The help file documents every subcommand, option and stored result.
- R and Stata side by side matches the R functions to the Stata subcommands.
-
net get powerapecopiesshowcase.do, a commented tour of every subcommand that runs in five to ten minutes.
References
Hancock, G. R., & Feng, Y. (2026). nmax and the quest to restore caution, integrity, and practicality to the sample size planning process. Psychological Methods, 31(1), 174–199. https://doi.org/10.1037/met0000776
Riesthuis, P. (2024). Simulation-based power analyses for the smallest effect size of interest: A confidence-interval approach for minimum-effect and equivalence testing. Advances in Methods and Practices in Psychological Science, 7(2). https://doi.org/10.1177/25152459241240722