Why Include a Control Function?
A continuous exposure may be endogenous when unmeasured factors that affect exposure also affect the outcome. For example, people with greater unmeasured health need may receive more treatment and also have worse outcomes. In a nonlinear outcome model, simply replacing the exposure with its first-stage fitted value generally does not recover the desired structural relationship. A control function instead includes a first-stage residual, or another suitable control, in the outcome equation to account for the dependence between exposure and unobserved outcome determinants.
This tutorial develops the approach for a continuous exposure \(D\), an outcome \(Y\), an instrument \(Z\), and observed covariates \(X\). The main illustration uses a binary outcome and a logit response. The same logic applies with a different nonlinear response model only when its corresponding distributional and control-function assumptions are appropriate.
Model and Assumptions
Consider the first-stage model
Here \(V\) is the unobserved component of exposure after accounting for \(Z\) and \(X\). The first-stage conditional-mean restriction is \(E(V\mid Z,X)=0\). For a nonlinear first stage, an ordinary linear residual need not be an adequate control; the appropriate generalized residual must be defined from that first-stage model.
For a binary outcome, one structural specification that directly justifies a logit control-function equation is
The unit-scale logistic assumption fixes the scale of the latent outcome error. The model also requires a correctly specified first-stage conditional mean, instrument relevance (\(\pi_Z\ne 0\), or a nonzero coefficient vector for multiple instruments), and instrument validity: conditional on \(X\), \(Z\) must be independent of the outcome disturbances and must affect \(Y\) only through \(D\) and the modeled control \(V\). In addition, the control-function restriction must be substantively credible: once \(V\) and \(X\) are held fixed, no remaining unobserved factor may jointly drive \(D\) and \(Y\).
These assumptions are stronger than merely observing that the first-stage residual predicts the outcome. They assert that \(V\) captures the relevant dependence and that the latent error \(\eta\) has a common logistic distribution independent of \(V,Z,X\). Other distributions and links can be used, but they require their own assumptions; a logistic conditional response does not follow from arbitrary joint normality or from an unspecified residual decomposition.
Deriving the Nonlinear Outcome Equation
Given \(D,X,V\), the latent index excluding \(\eta\) is \(q=\beta_0+\beta_DD+X^\mathsf{T}\beta_X+\rho V\). Because \(\eta\) is independent and has the stated logistic distribution,
The third line uses the symmetry of the unit-scale logistic distribution. Thus the estimable outcome equation is
This is residual inclusion: estimate the first stage, calculate \(\hat V=D-\hat\pi_0-\hat\pi_ZZ-X^\mathsf{T}\hat\pi_X\), and include \(\hat V\) along with \(D\) and \(X\) in the nonlinear outcome model. The derivation explains why the residual enters the outcome index under the stated structural assumptions. It does not make those assumptions true, and it does not turn a misspecified first stage into a valid control.
Worked Example: Constructing the First-Stage Residual
Suppose hypothetical exposure data follow the specified first-stage conditional mean \(D=1+0.8Z+0.5X+V\), where \(Z\) and \(X\) are binary. For four observations, the exposure values are as follows. Treat the coefficients as known for this arithmetic illustration; in an analysis they would normally be estimated.
| Observation | \(Z\) | \(X\) | Observed \(D\) | Fitted exposure | Residual \(\hat V\) |
|---|---|---|---|---|---|
| A | 0 | 0 | 1.2 | 1.0 | 0.2 |
| B | 1 | 0 | 2.0 | 1.8 | 0.2 |
| C | 0 | 1 | 1.3 | 1.5 | -0.2 |
| D | 1 | 1 | 2.2 | 2.3 | -0.1 |
For observation C, for example, the fitted exposure is \(1+0.8(0)+0.5(1)=1.5\), so \(\hat V=1.3-1.5=-0.2\). For observation D, \(\hat V=2.2-[1+0.8(1)+0.5(1)]=-0.1\). The residual records how far exposure is above or below its fitted conditional mean; it is not itself an instrument and should not be interpreted as a directly measured confounder.
Before using these residuals, check that the first-stage specification is suitable and that \(Z\) predicts \(D\) conditional on \(X\). If \(Z\) has little conditional predictive power, the control-function analysis can be weakly informed even if a nonlinear outcome model converges.
Estimating and Interpreting Treatment Contrasts
The coefficient \(\beta_D\) is a change in the latent log-odds index per unit of exposure, holding \(X\) and \(V\) fixed. It is not a probability difference. A probability-scale contrast compares the fitted response at two exposure values using the same covariates and control-function value. Under the structural model, this is a contrast for those \(X,V\) values; an average contrast for a target population requires averaging such contrasts over that population’s distribution of \(X,V\).
Worked Example: A Conditional Risk Contrast
Suppose a hypothetical fitted model is \(P(Y=1\mid D,X,V)=\operatorname{logit}^{-1}(-1.2+0.4D+0.3X+0.7V)\). Calculate the probability contrast between \(D=2\) and \(D=4\) for \(X=1,V=0.5\), holding the latter values fixed.
At \(D=2\), the index is \(-1.2+0.4(2)+0.3(1)+0.7(0.5)=0.25\). Therefore, \(p_2=1/[1+\exp(-0.25)]=0.5622\).
At \(D=4\), the index is \(-1.2+0.4(4)+0.3(1)+0.7(0.5)=1.05\). Therefore, \(p_4=1/[1+\exp(-1.05)]=0.7408\).
Using the unrounded probabilities, the contrast is \(p_4-p_2=0.7407749-0.5621765=0.1785984\), or about 0.1786. Under the model and its causal assumptions, this is an estimated increase of 17.86 percentage points in outcome probability for this \(X,V\) profile when exposure changes from 2 to 4 units. It is not the effect for every individual or necessarily the population-average effect.
Uncertainty calculations must account for estimating the first stage and then using its residuals. A standard error that treats \(\hat V\) as known may be too small. Joint estimation or a bootstrap that repeats both stages can propagate this uncertainty; the bootstrap must resample the independent sampling units and refit the first and second stages in each replicate.
Support, Diagnostics, and Failure Modes
The model can produce a contrast at any numeric \(D\), but that does not mean the contrast is identified by the observed data. For a requested contrast, examine whether exposure values are supported in relevant covariate and control-function regions. In particular, overlap of the overall exposure ranges can conceal a lack of support at particular \(X,V\) values. A result that depends on extrapolating the fitted logit curve far beyond observed exposure is driven heavily by functional-form assumptions.
Worked Example: Detecting an Unsupported Contrast
Suppose a hypothetical first stage is \(D=1+0.8Z+V\), with binary \(Z\), and the observed residual range in each instrument group is approximately \(-0.5\) to \(0.5\). For \(Z=0\), the corresponding exposure range is \(1+[-0.5,0.5]=[0.5,1.5]\); for \(Z=1\), it is \(1.8+[-0.5,0.5]=[1.3,2.3]\). Their overall overlap is only \([1.3,1.5]\).
At a fixed \(V=0\), the first stage predicts \(D=1.0\) when \(Z=0\) and \(D=1.8\) when \(Z=1\). It supplies no instrument-induced observations at both \(D=1\) and \(D=2\) for that same residual value. A model-based contrast between \(D=1\) and \(D=2\) at \(V=0\) may still be calculated, but it relies on extrapolating the structural outcome form rather than on local support for that contrast. The small overall overlap is not evidence of broad conditional support.
Useful checks include plotting exposure against \(Z\) within covariate strata, examining the distribution of first-stage residuals across instrument values, and plotting fitted outcome probabilities against exposure with the control function accounted for. Inspect whether the residual has a systematic nonlinear relationship with the outcome index or whether its association changes across exposure levels. Such patterns can signal misspecified functional form, though no single diagnostic proves the assumptions.
A linear residual term \(\rho V\) is a substantive restriction. If the relationship between \(V\) and the outcome index is nonlinear, or if its role changes with \(D\), the simple model may not adequately control endogeneity. Flexible functions or interactions may be considered when justified and supported by data; those choices add assumptions and should be reported transparently. A weak instrument, a poorly specified first stage, instrument invalidity, limited support, or an incorrect outcome link can each undermine interpretation.
How the Approach Differs From Alternatives
| Approach | What enters the outcome model | Main distinction |
|---|---|---|
| Control-function residual inclusion | Observed \(D\), covariates \(X\), and first-stage control \(\hat V\) | Models the nonlinear response while controlling for specified dependence through the first stage. |
| Predictor substitution | First-stage fitted exposure \(\hat D\) in place of \(D\) | Generally does not reproduce the structural nonlinear response; fitted-value substitution is not interchangeable with residual inclusion. |
| Linear two-stage least squares | A linear outcome equation with fitted exposure in the second stage | Targets a linear structural relationship under its assumptions, not the nonlinear probability model developed here. |
| Joint likelihood | A joint model for exposure and outcome | Can estimate their dependence together, but requires a correctly specified joint distribution and may be sensitive to those assumptions. |
These methods do not become equivalent just because they use the same instrument. The control-function approach is useful when a nonlinear conditional outcome model is substantively important and a credible control can be constructed. It is not a general repair for invalid instruments, unmeasured dependence left outside \(V\), or unsupported treatment contrasts.
Common Mistakes and Practical Interpretation
- Reading \(\beta_D\) as a probability change. It is a log-odds index coefficient in the logit model. Compute probabilities at specified \(D,X,V\) values to obtain a risk contrast.
- Assuming any residual works. The first-stage residual is valid as a control only under the first-stage and control-function assumptions. A poor first-stage model can produce a misleading control.
- Using fitted exposure as if it were residual inclusion. Substitution and residual inclusion are different procedures in nonlinear models and generally have different implications.
- Claiming causal effects from prediction alone. A causal reading depends on instrument validity, the structural outcome and control-function restrictions, and adequate support for the requested contrast.
- Ignoring generated-residual uncertainty. The residual is estimated; inference should account for both stages.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- In a linear first stage, what does the residual \(V\) represent, and what condition should it satisfy given \(Z,X\)?
- In the stated logit model, is \(\beta_D\) a probability difference? What does it measure?
- Why does observing a residual-outcome association not by itself prove that the control-function assumptions hold?
- Why can overall overlap in exposure ranges fail to support a contrast at fixed \(X,V\)?
- What uncertainty is missed if a second-stage standard error treats the estimated first-stage residuals as known?
Solutions. (1) \(V\) is the exposure component remaining after the first-stage conditional mean is removed; the model requires \(E(V\mid Z,X)=0\). (2) No. It is the change in the latent log-odds index per exposure unit, holding \(X,V\) fixed. (3) The association is compatible with misspecification; validity requires the stated control-function restriction and instrument assumptions. (4) Exposure may overlap overall but have little or no overlap within particular covariate and residual strata, forcing extrapolation. (5) It misses the uncertainty from estimating first-stage coefficients and constructing the residuals, which can be propagated by joint estimation or a two-stage bootstrap.