The Bivariate Probit Model For Endogenous Binary Treatment

Why Model the Outcome and Treatment Jointly?

When both an outcome and a treatment are binary, a probit model for the outcome alone does not account for the possibility that unobserved factors influencing treatment also influence the outcome. The bivariate probit model represents the two binary variables as thresholded latent indices and allows their latent errors to be correlated. An instrument can affect the treatment index while being excluded from the outcome index.

This is a joint parametric model: its likelihood uses the probabilities of observing each combination of outcome and treatment. The model’s assumptions are substantial, especially the joint normality of latent errors. Its coefficient on treatment is a parameter of the latent outcome index; it is not, by itself, a probability-scale treatment effect.

Latent-Index Specification and Assumptions

Definition: Let $Y\in\{0,1\}$ be a binary outcome, $D\in\{0,1\}$ a binary treatment, $Z$ a vector of excluded instruments, and $X$ a vector of included exogenous covariates. The bivariate probit model is $$ Y^* = X'\beta + \tau D + \varepsilon,\qquad Y = \mathbf{1}(Y^*>0), $$ $$ D^* = X'\gamma + Z'\pi + u,\qquad D = \mathbf{1}(D^*>0). $$ Here $Y^*$ and $D^*$ are latent continuous indices, and $\mathbf{1}(\cdot)$ equals 1 when its condition is true and 0 otherwise. Conditional on $X,Z$, the errors satisfy $$ \begin{pmatrix}\varepsilon\\u\end{pmatrix}\,\middle|\,X,Z \sim N\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\\rho&1\end{pmatrix} \right),\qquad -1<\rho<1. $$

The unit variances set the scales of the latent indices: binary observations reveal whether an index is above or below zero, not its unrestricted scale. The parameter $\rho$ is the correlation between the latent outcome error and the latent treatment error, conditional on $X,Z$. A nonzero $\rho$ allows unobserved determinants of treatment to be associated with unobserved determinants of the outcome. It is not the treatment coefficient $\tau$, nor is it itself an effect on the outcome probability.

The specification relies on several conditions. The threshold equations and their linear indices must be correctly specified; observations are treated as independent; and the conditional error distribution must be bivariate normal with the stated zero means, unit variances, and constant correlation. The included covariates $X$ enter both equations as written. The instrument enters the treatment index but not the outcome index. For an instrumental-variable interpretation, the instrument must also be conditionally independent of the structural errors given $X$, and its exclusion from the outcome equation must be credible. A statistical fit cannot establish that exclusion restriction.

How the Exclusion Restriction Enters

The exclusion restriction is visible in the equations: $Z$ appears in $D^*$, through $Z'\pi$, but does not appear directly in $Y^*$. Thus the model does not allow the instrument to shift the latent outcome index while holding treatment and $X$ fixed. The instrument may still be associated with the observed outcome through its effect on treatment and the joint distribution of the latent errors.

This restriction concerns the model’s structural outcome equation, not whether the instrument is statistically associated with the outcome in the data. An association can arise through treatment. Conversely, leaving a variable out of the outcome equation does not prove that it has no direct or otherwise prohibited connection to the outcome. Substantive knowledge is needed to assess that claim.

Worked Example: Writing the Model and Checking Exclusion

Suppose a hypothetical study records whether a patient receives a binary intervention, $D$, and whether a binary health outcome occurs, $Y$. Let $Z$ be a binary instrument and $X$ include age and baseline severity. A specification could be $$ Y^*=\beta_0+\beta_1\text{Age}+\beta_2\text{Severity}+\tau D+\varepsilon, $$ $$ D^*=\gamma_0+\gamma_1\text{Age}+\gamma_2\text{Severity}+\pi Z+u. $$

The exclusion restriction is represented by the absence of a $Z$ term in the first equation. If $\pi\ne 0$, the model allows $Z$ to shift the latent treatment index; it does not thereby establish that the instrument is valid. For the joint-error assumption, the investigator must also posit that, conditional on age, severity, and $Z$, $(\varepsilon,u)$ has a bivariate normal distribution with unit marginal variances and correlation $\rho$. Those are model assumptions, not facts verified by writing the equations.

The coefficient $\tau$ shifts the latent outcome index by $\tau$ units when $D$ changes from 0 to 1, holding $X$ and the latent error fixed. Because $Y$ is determined by whether $Y^*$ crosses zero, this statement does not mean that the outcome probability changes by $\tau$ percentage points.

Deriving the Likelihood Contribution

For an observation with covariates $x,z$, define the two index values $a=x'\beta+\tau d$ and $b=x'\gamma+z'\pi$, where $d$ is the observed treatment. Given $D=d$, the outcome equation has index $a+\varepsilon$; the treatment equation has index $b+u$. The observed $y,d$ specify whether each error is above or below its threshold:

\begin{aligned} Y=1 &\Longleftrightarrow \varepsilon>-a, &Y=0 &\Longleftrightarrow \varepsilon\le -a,\\ D=1 &\Longleftrightarrow u>-b, &D=0 &\Longleftrightarrow u\le -b. \end{aligned}

Write $s_y=2y-1$ and $s_d=2d-1$, so each is $+1$ for a value of 1 and $-1$ for a value of 0. The observed pair requires $s_y\varepsilon>-s_ya$ and $s_du>-s_db$. The transformed errors $(-s_y\varepsilon,-s_du)$ are standard bivariate normal with correlation $s_ys_d\rho$. Therefore the probability of the observed pair is a bivariate normal cumulative probability.

Likelihood contribution: Let $\Phi_2(r,q;c)$ denote the cumulative distribution function of a standard bivariate normal vector with correlation $c$, evaluated at $r,q$. For one observation, $$ P(Y=y,D=d\mid x,z) = \Phi_2\!\left(s_y a,\ s_d b;\ s_ys_d\rho\right). $$ For $n$ independent observations, the conditional likelihood is the product of these contributions, and the log-likelihood is the sum of their logarithms.

The sign on the correlation matters. When the observed outcome and treatment agree (both 1 or both 0), $s_ys_d=1$, so the relevant correlation is $\rho$. When they disagree, the transformed correlation is $-\rho$. The formula accommodates all four possible observed pairs without fitting four unrelated probabilities.

Worked Example: A Likelihood Contribution at Zero Indices

Consider a hypothetical observation for which $a=0$ and $b=0$, with a model value $\rho=0.5$. Suppose the observed values are $Y=1,D=1$. Then $s_y=s_d=1$, and the contribution is $$ P(Y=1,D=1\mid x,z)=\Phi_2(0,0;0.5). $$

For a standard bivariate normal distribution, at two zero thresholds this probability is $\frac14+\frac{\arcsin(c)}{2\pi}$, where $c$ is the correlation. Substituting $c=0.5$, $$ \Phi_2(0,0;0.5)=\frac14+\frac{\arcsin(0.5)}{2\pi} =\frac14+\frac{\pi/6}{2\pi} =\frac13\approx 0.3333. $$

For the same zero indices and correlation, the discordant pair $Y=1,D=0$ has transformed correlation $-0.5$, so its probability is $\Phi_2(0,0;-0.5)=\frac14-\frac{\pi/6}{2\pi}=\frac16\approx0.1667$. The four cell probabilities are $1/3,1/6,1/6,1/3$, and sum to 1. This verifies that the sign change in the likelihood formula distinguishes concordant from discordant observations.

Worked Example: Computing a Row Probability Under Independence

For a hypothetical observation, let the fitted indices be $a=0.5$ and $b=-0.2$, let $\rho=0$, and suppose $Y=1,D=0$. Then $s_y=1$, $s_d=-1$, and the formula gives $$ P(Y=1,D=0\mid x,z)=\Phi_2(0.5,0.2;0). $$

At correlation zero, the two standard normal variables are independent, so the bivariate CDF factors: $$ \Phi_2(0.5,0.2;0)=\Phi(0.5)\Phi(0.2) \approx (0.6915)(0.5793)=0.4006. $$

Thus, under these specified indices and zero error correlation, the model assigns about 0.4006 probability to this particular outcome-treatment pair. The calculation is a joint probability, not the outcome probability alone. If $\rho\ne0$, the product of univariate probabilities is generally not the correct contribution.

What the Correlation Does—and Does Not—Say

If $\rho=0$, the model’s conditional bivariate normal errors are independent, and the joint probability factors into the treatment and outcome probabilities for a given treatment value. If $\rho\ne0$, the errors are dependent, so the likelihood must use their joint distribution. This is the model’s way of allowing latent unobservables affecting treatment to be related to latent unobservables affecting the outcome.

A nonzero fitted $\rho$ is not proof that treatment is causally endogenous, and a fitted value near zero does not prove that the instrument is valid or that the entire model is correctly specified. Its meaning depends on the assumed latent-index form and error distribution. Likewise, $\tau$ is a coefficient on the latent outcome index, not a risk difference or a probability ratio. Probability-scale effects require translating the model’s probabilities and are a separate interpretive task.

Diagnostics, Alternatives, and Limitations

The likelihood is only as credible as the assumptions that produce it. Review whether the outcome and treatment indices have defensible covariates and functional forms, whether the instrument’s exclusion can be justified, and whether the joint normality and constant-correlation assumptions are plausible. Sparse combinations of $Y$ and $D$, especially within covariate or instrument groups, can make the fitted model sensitive to a small number of observations. Convergence problems, extreme estimates, or large changes after modest specification changes are warnings to investigate rather than evidence that the model is correct.

A bivariate probit uses a single joint likelihood and directly specifies the dependence between latent errors. Separate probit fits do not incorporate that shared error correlation. Residual-inclusion approaches instead add a first-stage-derived control to an outcome model; they are a different estimation strategy and require their own assumptions. Neither an alternative estimator nor a good fit automatically resolves instrument validity. The bivariate probit’s parametric distributional assumptions can be consequential, so results should not be presented as though the joint-normal model followed from the observed binary data alone.

A common mistake is to treat a nonzero $\rho$ as the estimated treatment effect, or to claim that omitting $Z$ from the outcome equation proves the exclusion restriction. The correction is to keep three ideas distinct: $\tau$ is a latent-index treatment coefficient; $\rho$ is latent-error correlation; and exclusion is a substantive restriction on the outcome equation that must be defended, not merely imposed in software.

Write both threshold equations.
Specify the outcome and treatment indices and identify the observed binary variables as indicators that the latent indices exceed zero.

State the joint-error assumptions.
Give the conditional means, variances, correlation, and distribution, along with the independence-across-observations assumption.

Check the instrument placement.
Confirm that the instrument appears in the treatment index and is excluded from the outcome index; separately justify its conditional exogeneity and exclusion.

Form the joint likelihood.
Use the observed signs of $Y,D$ to select the appropriate bivariate normal CDF probability and multiply contributions across observations.

Check Your Understanding

Answer each question, then compare with the solutions below.

In the bivariate probit model, what does $\rho$ represent, and what does it not represent?
For an observation with $Y=0,D=1$, what is $s_ys_d$, and what correlation enters its likelihood contribution?
Suppose an instrument $Z$ appears directly in the outcome index as well as the treatment index. Which model restriction has changed?
At zero indices and $\rho=0$, what is the probability of observing any one of the four $Y,D$ pairs?
Why is $\tau$ not automatically a change in outcome probability measured in percentage points?

Solutions. (1) $\rho$ is the conditional correlation of the latent outcome and treatment errors; it is neither the treatment coefficient nor a probability-scale treatment effect. (2) $s_y=-1$, $s_d=1$, so $s_ys_d=-1$, and the likelihood uses correlation $-\rho$. (3) The instrument exclusion restriction is no longer imposed: the model now permits a direct instrument term in the latent outcome index. (4) The zero-correlation bivariate normal CDF factors, so each pair has probability $\Phi(0)\Phi(0)=0.5(0.5)=0.25$. (5) $\tau$ shifts the latent index; the observed outcome is a thresholded indicator, so its probability change depends on the index values and the normal distribution, not just on $\tau$.