Estimating And Interpreting Bivariate Probit Effects

From a Fitted Model to Probability-Scale Effects

A bivariate probit model describes a binary treatment and a binary outcome jointly, allowing their latent errors to be correlated. Fitting the model produces coefficients on latent indices. Those coefficients are not, by themselves, changes in outcome probabilities. To report effects on a clinically or substantively meaningful scale, calculate predicted probabilities under each treatment condition and contrast them.

The focus here is the translation from estimated model to counterfactual risks and average treatment effects (ATEs). The earlier material on identification in bivariate probit models establishes why excluded instruments and model structure can identify the parameters. The interpretation below still depends on the model and causal assumptions being adequate.

Model: Let $D_i$ be a binary treatment, $Y_i$ a binary outcome, $X_i$ a vector of included covariates, and $Z_i$ a vector containing the included covariates and at least one excluded instrument. The bivariate probit model is $$ D_i=1[Z_i^\top\gamma+u_i>0], \qquad Y_i=1[X_i^\top\theta+\beta D_i+\epsilon_i>0]. $$ The errors have a standard bivariate normal distribution: $$ \begin{pmatrix}u_i\\ \epsilon_i\end{pmatrix} \sim N\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\ \rho&1\end{pmatrix} \right). $$ The latent-error variances are fixed at one to set the probit scales. The parameter $\rho$ represents dependence between unobserved determinants of treatment and outcome.

For a causal interpretation, the treatment must be well-defined, with consistency and no interference between units. The instrument must be relevant, excluded from the outcome equation conditional on $X_i$, and independent of the structural errors conditional on included covariates. The specified latent-index forms and joint normal error distribution must also be adequate. Identification additionally requires sufficient variation and the rank or support conditions described in the earlier material. These are substantive and statistical assumptions, not consequences of estimating the model.

If $\rho\ne 0$, unobserved factors affecting treatment are associated with unobserved factors affecting outcome. This is why a probit outcome model that simply treats $D_i$ as exogenous may not estimate the structural outcome parameters consistently. The joint model instead estimates treatment and outcome equations together.

How the Joint Likelihood Uses the Data

Write the treatment index as $\eta_i=Z_i^\top\gamma$, and the outcome index under treatment value $d$ as $q_{id}=X_i^\top\theta+\beta d$. Let $\Phi$ be the standard normal cumulative distribution function and let $\Phi_2(a,b;\rho)$ be the bivariate standard normal probability that both variables are at most their respective thresholds $a,b$, with correlation $\rho$.

The model assigns a probability to each observed pair $(D_i,Y_i)$. For the treated, the probability of the outcome being one is a joint probability involving the two latent errors. For the untreated, the outcome equation uses $q_{i0}$, not $q_{i1}$. The four cell probabilities are:

\begin{aligned} P(D_i=1,Y_i=1\mid X_i,Z_i) &=\Phi_2(\eta_i,q_{i1};\rho),\\ P(D_i=1,Y_i=0\mid X_i,Z_i) &=\Phi(\eta_i)-\Phi_2(\eta_i,q_{i1};\rho),\\ P(D_i=0,Y_i=1\mid X_i,Z_i) &=\Phi(q_{i0})-\Phi_2(\eta_i,q_{i0};\rho),\\ P(D_i=0,Y_i=0\mid X_i,Z_i) &=1-\Phi(\eta_i)-\Phi(q_{i0})+\Phi_2(\eta_i,q_{i0};\rho). \end{aligned}

These expressions can be seen by changing signs on both latent errors: the events $D_i=1$ and $Y_i=1$ become lower-tail threshold events with thresholds $\eta_i$ and $q_{i1}$. The other cell probabilities follow by subtracting from the corresponding marginal probabilities or taking the remaining probability. They are nonnegative under the model and sum to one.

Maximum likelihood chooses $\gamma,\theta,\beta,\rho$ to maximize the product of the probability assigned to each unit’s observed treatment-outcome pair. Equivalently, it maximizes the sum of their log probabilities. This is joint estimation: it is not ordinary probit after replacing treatment with a first-stage fitted value.

Worked Example: A Joint Likelihood Contribution

Consider a hypothetical unit with fitted treatment index $\eta_i=0$, outcome indices $q_{i0}=q_{i1}=0$, and fitted error correlation $\rho=0.5$. This illustrates a likelihood contribution; it is not a claim that these values were estimated from a real data set.

At zero thresholds, the bivariate standard normal probability is $\Phi_2(0,0;0.5)=\frac14+\frac{\arcsin(0.5)}{2\pi}=\frac14+\frac1{12}=\frac13$. Also, $\Phi(0)=0.5$. Thus $P(D_i=1,Y_i=1)=1/3$ and $P(D_i=1,Y_i=0)=0.5-1/3=1/6$.

For the untreated cells, $P(D_i=0,Y_i=1)=0.5-1/3=1/6$, and $P(D_i=0,Y_i=0)=1-0.5-0.5+1/3=1/3$. The four probabilities sum to $1/3+1/6+1/6+1/3=1$. If the observed unit has $D_i=0,Y_i=1$, its likelihood contribution is $1/6$, and its log-likelihood contribution is $\log(1/6)$. A different observed cell uses that cell’s probability.

The positive correlation places relatively more probability in the concordant cells $(1,1)$ and $(0,0)$ than in the discordant cells here. This dependence is part of what the joint likelihood models.

Counterfactual Risks and the Average Treatment Effect

The structural outcome equation defines potential outcomes under each treatment value as $Y_i(d)=1[X_i^\top\theta+\beta d+\epsilon_i>0]$, for $d\in\{0,1\}$. Because the marginal distribution of $\epsilon_i$ is standard normal, the model-implied counterfactual probability conditional on covariates is $P(Y_i(d)=1\mid X_i)=\Phi(X_i^\top\theta+\beta d)$. The instrument and treatment-assignment equation are essential to estimating the structural parameters when treatment is endogenous, but the marginal risk formula follows from the outcome equation.

For a particular covariate profile, the probability-scale treatment effect is the risk difference $\Phi(X_i^\top\theta+\beta)-\Phi(X_i^\top\theta)$. It is generally not equal to $\beta$. For a target population represented by $n$ units with covariates $X_i$, the model-based ATE is obtained by averaging their predicted risk differences:

\widehat{\mathrm{ATE}} =\frac{1}{n}\sum_{i=1}^{n} \left[ \Phi(X_i^\top\hat{\theta}+\hat{\beta}) -\Phi(X_i^\top\hat{\theta}) \right].

This is a risk difference: for a binary outcome it is measured in probability units, often reported as a percentage-point difference. The average depends on the covariate distribution used. Averaging over the estimation sample targets that sample’s covariate distribution; another target population requires averaging over its covariates instead.

Worked Example: A Profile-Specific Risk Difference

Suppose a hypothetical fitted model gives an untreated outcome index of $-0.2$ for a specified covariate profile and a treatment coefficient $\hat{\beta}=0.7$. The treated index is therefore $-0.2+0.7=0.5$.

The untreated counterfactual risk is $\Phi(-0.2)=0.42074029$, and the treated counterfactual risk is $\Phi(0.5)=0.69146246$. Their difference, calculated before rounding, is $0.69146246-0.42074029=0.27072217$, or approximately $0.2707$. The model-implied risk is about 27.07 percentage points higher under treatment for this profile.

The coefficient $0.7$ is a shift in the latent probit index, not a 70-percentage-point increase in outcome probability. The probability change depends on the starting index: the standard normal curve has different slopes at different points.

Why the Coefficient Is Not the ATE

For a fixed covariate profile, the risk difference uses a nonlinear transformation of both indices. If the treatment coefficient is small, a local approximation is $\Phi(q+\beta)-\Phi(q)\approx\phi(q)\beta$, where $q=X_i^\top\theta$ and $\phi$ is the standard normal density. Since $\phi(q)$ varies with $q$, even this approximation varies across profiles. The exact contrast is preferable when calculating effects.

A marginal effect obtained by differentiating a probit probability with respect to a continuous regressor is also not automatically the treatment effect for a binary treatment. For binary $D_i$, calculate the discrete contrast between $D_i=1$ and $D_i=0$, then average those contrasts if an ATE is wanted.

Worked Example: Averaging Counterfactual Risks

Suppose a hypothetical fitted model has $\hat{\beta}=0.5$. For three target-population profiles, the untreated indices $X_i^\top\hat{\theta}$ are $-1,0,1$. The treated indices are consequently $-0.5,0.5,1.5$.

For the first profile, the risks are $\Phi(-1)=0.1586553$ and $\Phi(-0.5)=0.3085375$, giving a risk difference of $0.1498822$. For the second, they are $\Phi(0)=0.5000000$ and $\Phi(0.5)=0.6914625$, giving $0.1914625$. For the third, they are $\Phi(1)=0.8413447$ and $\Phi(1.5)=0.9331928$, giving $0.0918481$.

The average untreated risk is $(0.1586553+0.5000000+0.8413447)/3=0.5000000$. The average treated risk is $(0.3085375+0.6914625+0.9331928)/3=0.6443976$. Thus the ATE is $0.6443976-0.5000000=0.1443976$, approximately $0.1444$, or 14.44 percentage points. It is also the average of the three profile-specific contrasts. It is not the latent-index coefficient $0.5$.

Observed-Group Risks Are Different Quantities

A counterfactual risk under intervention, $P(Y_i(d)=1\mid X_i)$, should not be confused with an observed-group probability such as $P(Y_i=1\mid D_i=1,X_i,Z_i)$. The latter conditions on treatment assignment, which selects units partly through $u_i$. When $\rho\ne0$, that selection is associated with $\epsilon_i$, so observed treated and untreated groups can have different outcome-error distributions.

For example, the model gives $P(Y_i=1\mid D_i=1,X_i,Z_i)=\Phi_2(\eta_i,q_{i1};\rho)/\Phi(\eta_i)$. For untreated units it gives $P(Y_i=1\mid D_i=0,X_i,Z_i)=1-\{1-\Phi(\eta_i)-\Phi(q_{i0})+\Phi_2(\eta_i,q_{i0};\rho)\}/\{1-\Phi(\eta_i)\}$. These are selection-conditioned predictions, not the intervention risks used in the ATE formula. The distinction matters because an observed treated-versus-untreated risk contrast need not equal a causal effect.

Uncertainty, Diagnostics, and Limitations

An estimated ATE is itself a function of estimated parameters and, for a sample average, of the covariate values used for standardization. Its uncertainty should therefore reflect uncertainty in the fitted model. A delta-method calculation can propagate the estimated parameter covariance through the ATE formula; a bootstrap that refits the full joint model and recalculates the ATE is another option when its sampling design is respected. Report an interval and the target population, outcome, treatment contrast, and units.

Useful checks include confirming model convergence, examining whether $\hat{\rho}$ is near a boundary, checking that fitted cell probabilities are valid, and comparing predicted and observed treatment-outcome patterns across relevant covariate or instrument groups. Investigate sparse treatment-outcome cells, weak instrument variation, sensitivity to covariate specification, and whether estimates change substantially under plausible alternative distributional specifications. These checks can reveal trouble; they cannot prove the exclusion restriction or joint normality.

A naïve probit outcome model is simpler but does not account for endogenous treatment. A linear IV model works on a different functional and effect scale and does not impose this joint bivariate normal structure. Bivariate probit can deliver model-based counterfactual probabilities, but its conclusions may rely strongly on distributional and functional-form assumptions, especially when exclusion-based variation is limited. A precisely reported estimate does not remove those limitations.

Common pitfall: “The treatment coefficient is $0.7$, so treatment raises the outcome probability by 70 percentage points.” This is incorrect: $0.7$ shifts the latent probit index. Calculate $\Phi(X_i^\top\hat{\theta}+\hat{\beta})-\Phi(X_i^\top\hat{\theta})$ for each profile, then average the contrasts for an ATE.

Check Your Understanding

Use the definitions and calculations above to answer these questions. Solutions follow in the same box.

If the untreated index is $0$ and $\hat{\beta}=0.4$, what are the two counterfactual risks and their risk difference?
Why can the model-implied ATE differ from $\hat{\beta}$, even when the coefficient is estimated precisely?
In the joint-likelihood example, what is the likelihood contribution for an observed unit with $D_i=1,Y_i=0$?
Why is $P(Y_i=1\mid D_i=1,X_i,Z_i)$ not generally the treated counterfactual risk?
What does the ATE average depend on besides the fitted parameters?

Solutions. (1) The untreated risk is $\Phi(0)=0.5000$; the treated risk is $\Phi(0.4)=0.6554$. The risk difference is $0.6554-0.5000=0.1554$, or 15.54 percentage points. (2) The coefficient shifts a latent index, while the probability contrast depends on the starting index through the nonlinear $\Phi$ function; the ATE also averages profile-specific contrasts. Precision does not change the estimand’s scale. (3) It is $1/6$, as calculated from $\Phi(0)-\Phi_2(0,0;0.5)=0.5-1/3$. (4) Observed treatment conditions on the treatment-selection process; correlated latent errors can make the outcome error distribution in the treated group differ from its marginal distribution under intervention. (5) It depends on the covariate distribution over which profile-specific counterfactual risks are averaged, and therefore on the chosen target population.