IV designs: endogenous focal variables via control functions
Source:vignettes/iv-designs.Rmd
iv-designs.RmdWhy an IV module?
Many focal variables of interest are not assigned: training take-up,
board structure, R&D intensity. When the focal variable is
endogenous – correlated with the structural error – a plain probit APE
is biased, and the modern flexible answer is the
control-function (CF) approach (Wooldridge, 2015):
model the endogenous variable in a first stage on instruments, and carry
its (generalized) residual into the outcome probit as an extra
regressor. Stata 18.5 ships this as cfprobit, with a clever
stacked method-of-moments variance that avoids bootstrapping
entirely.
powerape replicates that estimator exactly – point
estimates, margins-style ASF effects, and the stacked
robust and cluster-robust standard errors, verified against Stata to
about 1e-7 – and then wraps the usual power machinery around it. The
estimand is unchanged: the ASF-based APE (or AIE) of the focal variable,
which is what margins reports after
cfprobit.
Step 1: describe the world
Two knobs carry the IV structure, and both are quantities you can defend:
-
endogeneity– the correlation between the structural and first-stage errors. This is why you need instruments; 0 means plain probit would have been fine. -
iv_strength– the share of the focal’s variance (latent variance for a binary focal) the instruments explain. This is the first-stage strength you would defend with an F statistic in the paper.
library(powerape)
d <- ape_dgp_iv(
focal = pa_var("training", "binary", p = 0.4), # probit first stage
covariates = list(pa_var("size", "normal")),
instruments = pa_var("subsidy", "binary", p = 0.5),
endogeneity = 0.4,
iv_strength = 0.2,
baseline = 0.30,
signal = 0.10
)
d <- set_ape(d, target = 0.10)
d
#> powerape DGP -- probit
#> focal: training (binary, p = 0.40)
#> covariates: size (normal, mean 0.00, sd 1.00)
#> iv: endogenous focal (probit first stage; instruments: subsidy); endogeneity rho = 0.40, first-stage strength = 0.20; control-function probit, stacked robust SEs
#> baseline P(Y=1 | reference) = 0.300, nuisance signal (latent pseudo-R2) = 0.100
#> true APE pinned at +0.1000 (beta_focal = 0.2857, implied P(Y=1 | focal=1) = 0.400)powerape dgp d, iv focal(training, binary p(0.4)) covariates((size, normal)) ///
instruments((subsidy, binary p(0.5))) endogeneity(0.4) ivstrength(0.2) ///
baseline(0.30) signal(0.10)
powerape setape d, target(0.10)powerape DGP -- probit
focal: training (binary, p = 0.40)
covariates: size (normal, mean 0.00, sd 1.00)
iv: endogenous focal (probit first stage; instruments: subsidy); endogeneity rho = 0.40,
first-stage strength = 0.20; control-function probit, stacked robust SEs
baseline P(Y=1 | reference) = 0.300, nuisance signal (latent pseudo-R2) = 0.100
true APE pinned at + 0.1000 (beta_focal = 0.2857, implied P(Y=1 | focal=1) = 0.400)A continuous focal (pa_var("rd", "normal", ...)) gets a
linear first stage instead; the pinning step is identical.
Step 2: power – and the IV price
ape_power(d, n = 1200, claim = "detect", nsim = 400, seed = 1)
#> powerape -- detection claim (CI excludes 0, directional)
#> probit, n = 1200, assumed true APE +0.1000, 95% CI (alpha = 0.05), nsim = 400, stacked robust SEs
#> power = 0.210 (MCSE 0.020)
#> outcomes: detect 0.210 | inconclusive 0.790 | failed 0.000powerape -- detection claim (CI excludes 0, directional)
probit, n = 1200, assumed true APE + 0.1000, 95% CI (alpha = .05), nsim = 400, stacked robust SEs
fit: cfprobit via engine(stata), effects(closed)
power = 0.250 (MCSE 0.022)
outcomes: detect 0.250 | inconclusive 0.750 | failed 0.000Compare that with the same world without endogeneity, estimated by plain probit:
ds <- ape_dgp(
focal = pa_var("training", "binary", p = 0.4),
covariates = list(pa_var("size", "normal")),
baseline = 0.30,
signal = 0.10
)
ds <- set_ape(ds, target = 0.10)
ape_power(ds, n = 1200, claim = "detect", nsim = 400, seed = 2)
#> powerape -- detection claim (CI excludes 0, directional)
#> probit, n = 1200, assumed true APE +0.1000, 95% CI (alpha = 0.05), nsim = 400
#> power = 0.955 (MCSE 0.010)
#> outcomes: detect 0.955 | inconclusive 0.045 | failed 0.000powerape dgp ds, focal(training, binary p(0.4)) covariates((size, normal)) ///
baseline(0.30) signal(0.10)
powerape setape ds, target(0.10)
powerape power ds, n(1200) claim(detect) nsim(400) seed(2)powerape -- detection claim (CI excludes 0, directional)
probit, n = 1200, assumed true APE + 0.1000, 95% CI (alpha = .05), nsim = 400
fit: probit via engine(stata), effects(closed)
power = 0.970 (MCSE 0.009)
outcomes: detect 0.970 | inconclusive 0.030 | failed 0.000The gap is the IV price: the CF estimator identifies
the effect from the instrumented variation only, so its standard error
is larger by roughly sqrt(1 / iv_strength) – the familiar
2SLS variance logic. Studies that budget their sample size for an
exogenous analysis and then switch to IV at review time are,
mechanically, badly underpowered. This pair of calls prices that in
advance.
Since endogeneity and iv_strength are
guesses, sweep them:
ape_robust(d, n = 1200, claim = "detect",
vary = list(iv_strength = c(0.1, 0.3)), grid_points = 3,
nsim = 300, seed = 3, nmax = FALSE)
#> powerape robustness sweep -- detect claim, APE, pin = ape
#> n = 1200, nsim = 300 per scenario, 3 scenario(s) over: iv_strength
#> power range [0.133, 0.353]; worst scenario:
#> iv_strength implied_effect power mcse separated
#> 0.1 0.1 0.1333333 0.01962614 0
#> scenarios:
#> iv_strength implied_effect power mcse separated
#> 0.1 0.1 0.1333333 0.01962614 0
#> 0.2 0.1 0.2200000 0.02391652 0
#> 0.3 0.1 0.3533333 0.02759764 0powerape robust d, n(1200) claim(detect) vary(ivstrength(0.1 0.3)) gridpoints(3) ///
nsim(300) seed(3) nonmax engine(mata)powerape robustness sweep -- detect claim, APE, pin = ape
n = 1200, nsim = 300 per scenario, 3 scenario(s) over: ivstrength
power range [0.133, 0.340]; worst scenario: row 1
fit: Mata engine (Fisher scoring; control-function probit)
row ivstrength implied_eff power mcse separated
1 .1 0.1000 0.1333 0.0196 0.0000
2 .2 0.1000 0.2033 0.0232 0.0000
3 .3 0.1000 0.3400 0.0273 0.0000Endogenous interactions
Moderation with an endogenous focal is the fully general case: does a
(self-selected) training program work differently in public firms? The
estimating model gains the focal-by-moderator term, the control
function, and the control-function-by-moderator interaction (Stata:
interact() + mainonly()), and the first stage
includes the moderator and the instrument-by-moderator
interactions – as instruments must in interaction models.
di <- ape_dgp_iv(
focal = pa_var("training", "binary", p = 0.4),
moderator = pa_var("public", "binary", p = 0.5),
covariates = list(pa_var("size", "normal")),
instruments = pa_var("subsidy", "binary", p = 0.5),
endogeneity = 0.4,
iv_strength = 0.25,
baseline = 0.30,
signal = 0.10
)
di <- set_aie(di, target = 0.06, main_focal = 0.10, main_moderator = 0.05)
ape_power(di, n = 2000, claim = "detect", nsim = 250, seed = 4)
#> powerape -- detection claim (CI excludes 0, directional)
#> probit, n = 2000, assumed true AIE +0.0600, 95% CI (alpha = 0.05), nsim = 250, stacked robust SEs
#> power = 0.064 (MCSE 0.015)
#> outcomes: detect 0.064 | inconclusive 0.936 | failed 0.000powerape dgp di, iv focal(training, binary p(0.4)) ///
moderator(public, binary p(0.5)) covariates((size, normal)) ///
instruments((subsidy, binary p(0.5))) endogeneity(0.4) ///
ivstrength(0.25) baseline(0.30) signal(0.10)
powerape setaie di, target(0.06) mainfocal(0.10) mainmoderator(0.05)
powerape power di, n(2000) claim(detect) nsim(250) seed(4)powerape -- detection claim (CI excludes 0, directional)
probit, n = 2000, assumed true AIE + 0.0600, 95% CI (alpha = .05), nsim = 250, stacked robust SEs
fit: cfprobit via engine(stata), effects(closed)
power = 0.084 (MCSE 0.018)
outcomes: detect 0.084 | inconclusive 0.916 | failed 0.000The interaction tax and the IV price compound. Interaction claims from instrumented designs need very large samples, and it is far better to learn that at the design stage than from a null result.
What stands behind the numbers
On fixed datasets, powerape’s CF estimator reproduces Stata’s
cfprobit across first stages (linear, probit, fractional
probit, Poisson): coefficients and stacked SEs to ~1e-7 (robust and
clustered), and the margins ASF effects exactly, including
the interaction design run with Stata’s own syntax. One caveat we
discovered in the process: before the October 14, 2026 update of
StataNow 19.5, Stata’s margins rebuilt the control function
for fractional first stages with a binary-outcome formula that
ignores the fractional value – an internal inconsistency with its own
estimator, reproduced to 7 decimals in our validation battery and fixed
by StataCorp in that update. powerape uses the estimation-consistent
generalized residuals, as Stata’s margins now does.
Current scope: one endogenous focal variable; linear and probit first stages in the DGP route (fractional probit and Poisson are supported at the estimator level); binary exogenous moderators for AIE designs; robust (stacked sandwich) inference.