Inference With Estimated Control Functions

Why Estimated Residuals Affect Inference

In a linear control-function analysis, the first stage produces residuals that are then included in an outcome regression. Those residuals are not directly observed: they depend on estimated first-stage coefficients. The outcome regression therefore uses generated regressors. Its coefficient estimates can be consistent under suitable assumptions, but uncertainty from estimating the residuals also belongs in the variance calculation.

This tutorial considers a single endogenous regressor in a linear outcome equation. The central inference question is not whether to include a control function, but how to calculate uncertainty after estimating it. In particular, a second-stage regression calculation that holds the generated residuals fixed need not describe the sampling variation of the complete two-step estimator.

Here \(v_i\) is the population first-stage disturbance and \(\hat v_i=d_i-w_i'\hat\pi\) is its estimated counterpart. The second-stage regression uses \(\hat v_i\), not \(v_i\). The coefficient \(\rho\) captures the part of the outcome disturbance represented by the control function under the stated model. The coefficient of primary interest is often \(\beta_d\), interpreted as a conditional linear effect only under the model’s substantive assumptions; a control-function regression does not by itself establish causality.

Where the Additional Variance Enters

Let \(q_i(\pi)=(x_i',d_i,d_i-w_i'\pi)'\) be the second-stage regressor vector, and let \(\gamma=(\beta_x',\beta_d,\rho)'\). At the true first-stage parameter, the second-stage error is \(u_i=y_i-q_i(\pi)'\gamma\). The two stages can be written together as estimating equations: the outcome equation sets the sample average of \(q_i(\pi)u_i\) to zero, while the first stage sets the sample average of \(w_i(d_i-w_i'\pi)\) to zero.

The first equation depends on \(\pi\) through the estimated residual. To see the consequence, let \(e_r\) be a vector with a one in the control-function coordinate of \(q_i\) and zeros elsewhere. Then \(\partial q_i/\partial\pi'=-e_r w_i'\). Differentiating the outcome moment \(q_i u_i\) with respect to \(\pi\) gives

The first term reflects the change in the residual regressor itself; the second reflects the resulting change in the fitted outcome error. Thus, uncertainty in \(\hat\pi\) affects the second-stage moments even when the outcome regression is otherwise correctly specified.

A First-Order Variance Derivation

The derivation makes the correction explicit. Define \(Q=E[q_iq_i']\), \(R=E[w_iw_i']\), and \(D=E[D_i]\), evaluated at the true parameters. Assume \(Q\) and \(R\) are invertible. A first-order expansion of both sample moment equations around \((\gamma,\pi)\) gives the following influence function for the second-stage coefficient vector:

To verify the expression, the first-stage equation yields \(\operatorname{IF}_{\pi,i}=R^{-1}w_i v_i\). The outcome equation has derivative \(-Q\) with respect to \(\gamma\), and derivative \(D\) with respect to \(\pi\). Its first-order expansion is therefore \(0=q_i u_i-Q\,\operatorname{IF}_{\gamma,i}+D\,\operatorname{IF}_{\pi,i}\). Solving for \(\operatorname{IF}_{\gamma,i}\) gives the displayed formula. The variance of the estimator is approximately the variance of this influence function divided by \(n\). This derivation identifies precisely the first-stage contribution: it is not an arbitrary penalty, but the effect of first-stage estimation on the outcome moments.

This is a heteroskedasticity-robust form when the moment covariance is estimated from observation-level outer products. Under independent clusters, form cluster-level sums of moments before calculating \(\hat B\). If a homoskedastic model is imposed, a specialized covariance expression may be available, but that assumption must be justified rather than inferred from the use of ordinary least squares. In all cases, the first-stage and second-stage moments must be treated as parts of the same estimator.

Sandwich Inference in Practice

A practical implementation should retain the dependence between stages rather than treating \(\hat v_i\) as measured data. One route is to evaluate the stacked moments and their derivative matrix directly. Another is to implement the influence-function calculation using estimated quantities. The resulting covariance estimate is then used for standard errors, confidence intervals, and tests.

For example, if \(\hat\beta_d=2.40\) outcome units per unit of \(d_i\) and its full-sandwich standard error is \(0.30\), an approximate 95% large-sample confidence interval is \(2.40\pm1.96(0.30)\), or \([1.812,2.988]\) outcome units per unit of \(d_i\). Under the model and regularity conditions, this procedure has approximately 95% coverage across repeated samples. It does not mean there is a 95% probability that the fixed parameter lies in this particular interval.

Bootstrap Inference: Repeat the Whole Procedure

A bootstrap can estimate the sampling variation without relying solely on a plug-in sandwich formula, provided each resample reproduces the estimator that would be applied to new data. For independent observations, draw observations with replacement; in each bootstrap sample, estimate the first stage, construct new residuals, and fit the outcome equation again. For clustered data, resample independent clusters instead. The bootstrap standard error is the sample standard deviation of the bootstrap coefficient estimates.

The essential rule is to re-estimate the first stage in every resample. Resampling only the outcome regression while keeping the original \(\hat v_i\) fixed conditions on the estimated first stage and fails to reproduce the full two-step sampling variation. Bootstrap inference still requires suitable sampling and regularity conditions; it does not repair weak identification, a misspecified control-function model, or an invalid instrument.

When the Fixed-Regressor Standard Error Misleads

An ordinary second-stage regression calculation conditions on the residual values supplied to that regression. It can therefore omit the change in the outcome moments caused by a different first-stage estimate. When that omitted component is material, the reported standard error can be too small and confidence intervals too narrow. As the variance decomposition showed, the full correction also includes covariance terms, so it is not a theorem that every fixed-residual standard error must be smaller in every sample or model.

In the linear setting, there is an important connection to instrumental-variables estimation: with the same included exogenous regressors and first-stage projection, the coefficient on \(d_i\) in the residual-inclusion regression equals the corresponding just-identified projected-regressor coefficient. That algebraic equivalence does not justify using an arbitrary second-stage OLS covariance. Inference must still use a covariance estimator appropriate to the full estimator and its assumptions.

Common Pitfall

• Wrong reasoning: “The residuals are now columns in my outcome data set, so the usual second-stage standard errors account for them.”
• Correction: The residual values depend on \(\hat\pi\), which was estimated from the same sample. Use a joint sandwich or a bootstrap that re-estimates the first stage and reconstructs the residuals. State the sampling unit and the assumptions behind the covariance calculation.