Score-Based Estimation and Inference for the Heckman Model

From the Likelihood to Its Score

The full-information Heckman model uses both the observed outcomes of selected units and the selection indicators for all units. Its likelihood therefore has distinct contributions from selected and unselected observations. The score is the vector of derivatives of the log-likelihood with respect to its parameters. It is central to numerical maximum-likelihood estimation and provides the observation-level quantities needed for sandwich inference.

Let $S_i$ indicate whether unit $i$ is selected, let $x_i$ be its outcome regressors, and let $z_i$ be its selection regressors. The outcome $y_i$ is observed only when $S_i=1$. Write the outcome equation as $y_i=x_i'\beta+\varepsilon_i$, and the selection rule as $S_i=1$ when $z_i'\gamma+u_i>0$. The standard Gaussian specification assumes that $(\varepsilon_i,u_i)$ is jointly normal, with variances $\sigma^2$ and 1 and correlation $\rho$, where $\sigma>0$ and $-1<\rho<1$. Observations are independent across $i$, and the regressors are treated as given for the likelihood.

Define $e_i=y_i-x_i'\beta$, $q_i=e_i/\sigma$, and $a_i=z_i'\gamma$. Let $r=\sqrt{1-\rho^2}$. For selected observations, the likelihood uses the density of the observed outcome together with the probability of selection conditional on that outcome. For unselected observations, it uses only the probability of nonselection.

Definition: Define $\phi(\cdot)$ and $\Phi(\cdot)$ as the standard normal density and distribution functions. The observation-level log-likelihood contributions are $$ \ell_i = \begin{cases} -\log(\sigma)+\log\phi(q_i)+\log\Phi(t_i), & S_i=1,\\ \log\Phi(-a_i), & S_i=0, \end{cases} \qquad t_i=\frac{a_i+\rho q_i}{r}. $$ The sample log-likelihood is $\ell(\theta)=\sum_{i=1}^n\ell_i$, where $\theta=(\beta',\gamma',\sigma,\rho)'$.

The $\log\Phi(t_i)$ term is what links selection to the observed outcome: conditional on the outcome residual, selection depends on the correlated selection error. By contrast, an unselected unit has no observed $y_i$, so it contributes no outcome residual or outcome-density term. It still contributes information about $\gamma$ through its observed selection status.

Score Contributions

For a compact expression, define the inverse Mills ratio $\lambda(v)=\phi(v)/\Phi(v)$, and for an unselected observation define $m_i=\phi(a_i)/\Phi(-a_i)$. Differentiating each log-likelihood contribution gives the following scores in the original parameterization. A score for an individual parameter is the derivative of that observation's log-likelihood with respect to that parameter.

Result: For a selected observation, the score components are $$ \begin{aligned} s_{\beta i} &=\frac{x_i}{\sigma}\left(q_i-\frac{\rho}{r}\lambda(t_i)\right),\\ s_{\gamma i} &=\frac{z_i}{r}\lambda(t_i),\\ s_{\sigma i} &=\frac{q_i^2-1}{\sigma} -\frac{\rho q_i}{\sigma r}\lambda(t_i),\\ s_{\rho i} &=\lambda(t_i)\frac{q_i+\rho a_i}{r^3}. \end{aligned} $$ For an unselected observation, $$ s_{\beta i}=0,\qquad s_{\gamma i}=-z_i m_i,\qquad s_{\sigma i}=0,\qquad s_{\rho i}=0. $$

These expressions follow by the chain rule. For example, $\partial q_i/\partial\beta=-x_i/\sigma$, and $\partial\log\Phi(t_i)/\partial t_i=\lambda(t_i)$. Thus the derivative with respect to $\beta$ combines the normal outcome-density score $x_iq_i/\sigma$ and the selection-probability derivative $-x_i\rho\lambda(t_i)/(\sigma r)$. For the selection index, $\partial t_i/\partial\gamma=z_i/r$, giving $s_{\gamma i}=z_i\lambda(t_i)/r$ when selected. When unselected, differentiating $\log\Phi(-a_i)$ gives $-z_i\phi(a_i)/\Phi(-a_i)$. The remaining components follow by differentiating $q_i$, $t_i$, and the scale term $-\log\sigma$.

Worked Example: Score for a Selected Observation

Consider one hypothetical selected observation with scalar $x_i=z_i=1$, $y_i=2.0$, and parameter values $\beta=1.5$, $\gamma=0.2$, $\sigma=1.0$, and $\rho=0.4$. These values are for illustration, not data from a real study.

First, $e_i=2.0-1.5=0.5$, $q_i=0.5/1.0=0.5$, and $a_i=0.2$. Also, $r=\sqrt{1-0.4^2}=\sqrt{0.84}=0.9165$, and $t_i=(0.2+0.4(0.5))/0.9165=0.4364$. Using $\lambda(0.4364)\approx0.5424$, the outcome-coefficient score is $$ s_{\beta i}=\frac{1}{1.0}\left(0.5-\frac{0.4}{0.9165}(0.5424)\right)\approx0.263. $$ The selection-coefficient score is $s_{\gamma i}=0.5424/0.9165\approx0.592$. The scale score is $$ s_{\sigma i}=\frac{0.5^2-1}{1.0} -\frac{0.4(0.5)}{1.0(0.9165)}(0.5424)\approx-0.868. $$ Finally, $r^3\approx0.7699$, so $s_{\rho i}=0.5424(0.5+0.4(0.2))/0.7699\approx0.409$. These are derivatives at the stated parameter values, not estimates of the parameters themselves.

The log-likelihood contribution is also directly evaluable: $\ell_i=-\log(1.0)+\log\phi(0.5)+\log\Phi(0.4364)\approx-1.446$.

Worked Example: Score for an Unselected Observation

For a hypothetical unselected unit with $z_i=1$ and $\gamma=0.2$, its selection index is $a_i=0.2$. Its log-likelihood contribution is $\log\Phi(-0.2)\approx\log(0.4207)=-0.8657$. The selection score is $$ s_{\gamma i}=-\frac{\phi(0.2)}{\Phi(-0.2)} \approx-\frac{0.3910}{0.4207}=-0.9294. $$

This unit's outcome is unobserved, so its scores for $\beta$, $\sigma$, and $\rho$ are zero in its individual likelihood contribution. Its selection status nevertheless informs the selection equation. It would be incorrect to omit unselected units from the likelihood just because their outcomes are missing.

Numerical Maximum-Likelihood Estimation

The maximum-likelihood estimate is the parameter value that maximizes the summed log-likelihood. Except in special cases, the nonlinear $\log\Phi(t_i)$ terms mean the score equations do not have a closed-form solution. A numerical optimizer searches for a parameter vector where the likelihood is high and the total score is close to zero.

In practice, it is useful to optimize over unconstrained parameters $\tau=\log\sigma$ and $\kappa=\operatorname{atanh}(\rho)$, recovering $\sigma=\exp(\tau)>0$ and $\rho=\tanh(\kappa)\in(-1,1)$. The chain rule converts derivatives in the original parameterization: $\partial\ell/\partial\tau=\sigma\,\partial\ell/\partial\sigma$ and $\partial\ell/\partial\kappa=(1-\rho^2)\,\partial\ell/\partial\rho$. This avoids proposing inadmissible scale or correlation values. Stable numerical routines for $\log\Phi(v)$ and its associated ratio are important when $v$ is far into a tail, where direct evaluation can underflow or lose precision.

Construct the contributions.
For every unit, use its selection indicator; include the outcome and its density contribution only when $S_i=1$.

Choose starting values and parameterization.
Initialize $\beta,\gamma$, and unconstrained scale and correlation parameters. Multiple reasonable starting values can help reveal sensitivity to local solutions.

Maximize the summed log-likelihood.
Use a gradient-based or other suitable optimizer, supplying analytic scores when available or checking numerical derivatives carefully.

Check convergence and curvature.
Confirm that the log-likelihood has stabilized, the total score is small, parameter values are admissible, and the estimated curvature is suitable for inference.

A convergence flag alone is not enough. A nearly singular Hessian can indicate weak identification, poor scaling, or a flat likelihood; a solution at an extreme correlation can also merit scrutiny. Comparing starting values and inspecting the score and curvature help distinguish a reliable numerical optimum from an optimizer's stopping point.

Model-Based and Sandwich Covariance

Let $s_i(\hat\theta)$ denote the score vector for observation $i$ at the maximum-likelihood estimate. Define the observed information $A$ as the negative Hessian of the sample log-likelihood, evaluated at $\hat\theta$. Under correct specification and regularity conditions, the model-based covariance estimate is $A^{-1}$. This uses the likelihood's curvature to estimate uncertainty.

Definition: The empirical score outer-product matrix is $B=\sum_{i=1}^n s_i(\hat\theta)s_i(\hat\theta)'$. The sandwich covariance estimate is $$ \widehat{\operatorname{Var}}_{\mathrm{sand}}(\hat\theta) =A^{-1}BA^{-1}. $$ For independent observations, this allows the score variance to differ from the curvature implied by the assumed likelihood. If observations are clustered or otherwise dependent, the meat $B$ must instead be constructed from appropriately aggregated scores.

The justification is the first-order expansion of the estimating equation. Near the probability limit $\theta_0$, the total score is approximately its score at $\theta_0$ minus the curvature matrix times $\hat\theta-\theta_0$. Thus the estimator's variation is approximately the inverse curvature applied to the variation in the summed scores. If the scores have covariance $B$, this yields inverse-curvature times score-covariance times inverse-curvature. Under correct likelihood specification, the information identity makes the population score covariance equal to the expected information, so model-based and sandwich estimates agree asymptotically. In finite samples they need not be equal.

The sandwich estimate protects standard errors against some likelihood misspecification when the estimator still targets a well-defined pseudo-true parameter and regularity conditions hold. It does not make the Gaussian model correct, eliminate selection bias, or resolve weak identification. Both covariance estimates require a well-behaved estimate with adequate information in the data.

Worked Example: Comparing Standard Errors and Intervals

Suppose a hypothetical converged fit reports two parameter estimates, an outcome slope $\hat\beta=1.20$ and a selection coefficient $\hat\gamma=-0.40$. For a standalone two-parameter illustration, suppose the following are the full observed information and empirical score outer-product matrices: $$ A= \begin{pmatrix} 25&0\\ 0&16 \end{pmatrix}, \qquad B= \begin{pmatrix} 227.256&0\\ 0&69.809 \end{pmatrix}. $$

For these blocks, the model-based standard errors are the square roots of the diagonal entries of $A^{-1}$: $SE(\hat\beta)=\sqrt{1/25}=0.2000$ and $SE(\hat\gamma)=\sqrt{1/16}=0.2500$. The sandwich variances are $227.256/25^2=0.3636$ and $69.809/16^2=0.2727$, giving standard errors $0.6030$ and $0.5222$, respectively.

Using a normal approximation, the 95% interval for the outcome slope is $1.20\pm1.96(0.6030)=[0.0181,\,2.382]$. The corresponding interval for the selection coefficient is $-0.40\pm1.96(0.5222)=[-1.424,\,0.6235]$.

The first interval is expressed in outcome units per one-unit increase in the outcome regressor, and the second in selection-index units per one-unit increase in its regressor. Under the stated large-sample approximation, the interval procedure has 95% coverage across repeated samples under its assumptions; it is not a statement that the fixed parameter has a 95% probability of lying in this particular computed interval.

Common Pitfall and Practical Interpretation

A frequent error is to calculate standard errors from the selected outcome observations alone, treating the selection equation as if it had been estimated independently. Full-information estimation uses the likelihood contributions from both selected and unselected units, and the selected contribution itself couples the outcome and selection parameters. The relevant Hessian and score vectors must therefore correspond to the full likelihood actually maximized.

A second error is to call a sandwich covariance estimate a correction for all problems in a selection model. It adjusts the estimated covariance for a difference between score variation and likelihood curvature under suitable conditions. It does not repair an invalid observation rule, guarantee adequate support, or make a point estimate meaningful when the model is poorly identified. Numerical convergence, identification, and the substantive assumptions remain separate questions.

Key takeaway: The Heckman score has different contributions for selected and unselected observations: selected units inform all model components, while unselected units directly inform the selection equation. Use the full likelihood for estimation, inspect convergence and curvature, and choose covariance estimates with their assumptions in view.

Check Your Understanding

Try these questions before moving on.

Why does an unselected observation contribute to the likelihood even though its outcome is not observed?
In the selected-observation score, what two parts contribute to the derivative with respect to $\beta$?
Why can optimizing over $\log\sigma$ and $\operatorname{atanh}(\rho)$ be numerically convenient?
What matrices form the curvature and score-variation parts of the sandwich covariance estimate?
Give one limitation of sandwich standard errors in a sample-selection analysis.