Worked Binary-Outcome Analysis With Bivariate Probit

Why Model Treatment and Outcome Jointly?

When both treatment and outcome are binary, an ordinary probit outcome model can describe outcome probabilities, but treating an instrumented treatment as if it were exogenous can give misleading estimates. Bivariate probit addresses this by specifying a latent treatment-selection equation and a latent outcome equation together. Their errors may be correlated, allowing unobserved factors that affect treatment receipt to be associated with unobserved factors that affect the outcome.

This is a parametric IV approach: its causal interpretation depends not only on the instrument but also on the specified latent-index and error-distribution assumptions. The model yields intervention risks under treatment and no treatment, from which a risk difference can be calculated for a defined covariate distribution.

Model and Identifying Assumptions

The treatment equation includes the instrument; the outcome equation excludes it. This exclusion is the model’s representation of the claim that, conditional on \(X_i\), the instrument affects the outcome only by changing treatment. The parameter \(\rho\) represents residual dependence between the latent treatment and outcome equations. When \(\rho\ne0\), a probit outcome model that simply includes \(D_i\) while assuming it is exogenous is generally not adequate.

For a causal interpretation, the model requires more than successful numerical estimation. The key conditions are:

• Instrument relevance: conditional on \(X_i\), variation in \(Z_i\) changes the treatment probability. In the specified equation, the instrument coefficients must provide meaningful variation in \(Z_i\gamma\).
• Conditional instrument exogeneity: given \(X_i\), the instrument is independent of the latent errors \((u_i,\epsilon_i)\). This includes the absence of unmeasured common causes of instrument assignment and outcome.
• Exclusion: \(Z_i\) has no direct effect on the outcome index. It appears in the treatment equation but not the outcome equation.
• Correct model form: the latent-index equations and the joint normal distribution, including the constant correlation \(\rho\), are adequate descriptions of the data.
• Well-defined treatment and outcome: consistency and no interference are assumed, so each person’s outcome under a treatment level is well defined and does not depend on other people’s treatment.
• Useful support: the observed covariates and instrument provide enough overlap and variation to estimate the model. A nominally nonzero first-stage effect can still be too weak for stable inference.

The bivariate normality assumption is especially consequential. It supplies the functional form that allows the model to separate treatment effects from correlated latent errors. An instrument’s relevance and exclusion do not, by themselves, validate that distributional assumption.

From the Latent Model to the Likelihood

Write \(\eta_i=Z_i\gamma+X_i\delta\) for the treatment index and \(q_{id}=X_i\beta+\tau d\) for the outcome index if treatment were set to \(d\). Let \(\Phi\) be the standard normal cumulative distribution function and let \(\Phi_2(a,b;r)\) be the probability that two standard normal variables with correlation \(r\) are at most \(a\) and \(b\), respectively.

For an observed pair \((D_i,Y_i)=(d,y)\), the data imply that \(u_i\) lies above or below \(-\eta_i\), and that \(\epsilon_i\) lies above or below \(-q_{id}\). Define \(s_d=2d-1\) and \(s_y=2y-1\), each equal to \(-1\) or \(1\). The observed event is \(s_d(u_i+\eta_i)>0\) and \(s_y(\epsilon_i+q_{id})>0\). Changing signs to express this as a lower-tail bivariate normal probability gives

The correlation changes sign when exactly one of the two latent variables is sign-reversed. This sign change is important: using the same correlation for all four treatment-outcome cells gives incorrect likelihood contributions when \(\rho\ne0\). Conditional independence across sampled individuals gives the log likelihood

where \(\theta\) collects the coefficients and \(\rho\). Maximum likelihood estimation chooses parameter values that maximize this log likelihood. In practice, software evaluates the bivariate normal probabilities numerically and uses numerical optimization. The result should be checked for convergence and a well-behaved likelihood, not accepted solely because software returned coefficient estimates.

Intervention Risks and the Average Treatment Contrast

Under the outcome equation, setting treatment to \(d\) gives the latent outcome index \(X_i\beta+\tau d+\epsilon_i\). Since the marginal distribution of \(\epsilon_i\) is standard normal, the intervention risk conditional on \(X_i\) is

The correlation \(\rho\) affects the joint distribution of observed treatment and outcome, and hence the estimation of the model parameters. It does not enter this marginal intervention-risk expression: under an intervention setting \(D_i=d\), the treatment-selection error is not used to determine treatment, while the marginal outcome error remains standard normal.

For a target population represented by \(n\) covariate profiles, the model-standardized average treatment contrast on the risk-difference scale is

This is the average predicted risk under treatment minus the average predicted risk under no treatment, using the same covariate profiles in both predictions. It is measured in probability points (or can be multiplied by 100 and reported in percentage points). It is not generally equal to \(\widehat{\tau}\): \(\widehat{\tau}\) is a shift in the latent probit index, whereas the risk difference depends on the covariate profile and the normal CDF.

Estimation Workflow and Uncertainty

A confidence interval for the average contrast should reflect uncertainty in all parameters used to predict the two risks, including their covariance. A coefficient’s standard error alone is not a standard error for the standardized contrast. For a regular, well-identified maximum-likelihood fit, the delta method uses the gradient of the contrast with respect to the fitted parameters and their estimated joint covariance matrix. A bootstrap can provide an alternative, but it must repeat model fitting rather than hold the fitted coefficients fixed.

Diagnostics, Alternatives, and Limitations

Diagnostics cannot establish exclusion or conditional instrument exogeneity from the observed data alone. They can, however, reveal problems that make the analysis untrustworthy or highly dependent on its assumptions.

• Weak or sparse first-stage variation: compare treatment frequencies across instrument values within relevant covariate patterns. Near-identical probabilities, sparse cells, or near-separation can produce imprecise estimates and unstable predictions.
• Numerical instability: nonconvergence, a nearly singular information matrix, estimates that change substantially with starting values, or \(\widehat{\rho}\) near \(-1\) or \(1\) warrant investigation. A boundary estimate may signal weak information or a poorly fitting distributional model.
• Observed-data fit: compare predicted and observed joint treatment-outcome frequencies across meaningful covariate and instrument groups. Large discrepancies can indicate misspecified indices, omitted interactions, or a poor joint-error assumption.
• Prediction support: check that standardized covariate profiles are represented in the data and that predicted risks are not driven by extrapolation or near-certain fitted treatment probabilities.
• Assumption sensitivity: compare reasonable covariate specifications and assess whether conclusions depend heavily on the bivariate-normal latent-error structure. Stability is reassuring but cannot prove the assumptions.

An exogenous probit model is simpler, but assumes treatment can be entered without modeling its endogeneity; it does not become a valid IV analysis merely because an instrument is available. Control-function approaches model endogeneity differently and may be preferable under other specifications, but they also require their own assumptions and correct implementation. Bivariate probit is attractive when the binary latent-index structure and joint normality are defensible. Its limitations are the strong distributional assumptions, sensitivity to weak instruments, the difficulty of validating exclusion, and the possibility that a constant latent-error correlation is inadequate. Report the model as a parametric IV analysis, not as a design-based guarantee of causality.

Common Mistake