Likelihood Inference and Instability · Tutorial 133 of 499

Numerical Integration in Selection-Model Likelihoods

Set up likelihood integrals for missing outcomes and see how quadrature dimension and tail behavior can affect their numerical values and resulting estimates.

Advanced 12 min read

What You'll Learn

Write a selection-model likelihood contribution as an integral over the unobserved outcomes.
Identify the integration dimension for each missingness pattern.
Use Gaussian quadrature to approximate a one-dimensional normal expectation.
Explain why tensor-product quadrature can become costly as the number of missing outcomes increases.
Recognize how steep selection probabilities and tail behavior can challenge a quadrature approximation.
Distinguish numerical integration error from uncertainty in the statistical model.

Why Integration Enters the Likelihood

In a selection model, the probability that an outcome is observed may depend on the outcome itself. A subject whose outcome is missing therefore contributes information through the response mechanism, but the missing outcome is not available to plug into that mechanism. The likelihood must average over its possible values, weighting each value by the outcome model.

This averaging is an integral. In simple cases it is one-dimensional and can be approximated by a short weighted sum. With several missing outcomes, the integral becomes multidimensional; its dimension and the behavior of the integrand, especially in the tails, can make numerical likelihood evaluation demanding.

Definition: Let $Y_i$ be subject $i$'s full outcome vector, $R_i$ its response pattern, and $x_i$ observed covariates. Suppose the outcome model has density $f(y_i\mid x_i;\theta)$, and the selection model specifies $h(r_i\mid y_i,x_i;\gamma)=P(R_i=r_i\mid Y_i=y_i,x_i)$. If $y_{i,O}$ is observed and $y_{i,M}$ is missing, the likelihood contribution is $$ L_i(\theta,\gamma)=\int f(y_{i,O},y_{i,M}\mid x_i;\theta)\, h(r_i\mid y_{i,O},y_{i,M},x_i;\gamma)\,dy_{i,M}. $$ The integral is over the support of the missing components. Its dimension is the number of continuously valued missing components, unless model structure permits a reduction.

This expression assumes the specified outcome density and response probabilities describe the data-generating mechanism, and that subjects are independent conditional on their covariates, so their contributions multiply in the likelihood. The response probability must be evaluated for the actual observed pattern. If all outcomes for a subject are observed, there is no integration over missing outcomes. For a missing outcome, the likelihood combines its density with the probability of the response pattern that left it unobserved.

A useful rearrangement follows by factoring the outcome density into observed and conditional parts: $$ L_i=f(y_{i,O}\mid x_i;\theta) \int f(y_{i,M}\mid y_{i,O},x_i;\theta) h(r_i\mid y_{i,O},y_{i,M},x_i;\gamma)\,dy_{i,M}. $$ Thus the integral is an expectation under the conditional distribution of the missing outcomes given the observed ones. This form can make the integration dimension and the quantity being averaged especially clear.

From an Integral to a Weighted Sum

Numerical quadrature replaces an integral by a weighted sum of evaluations of its integrand. For example, if $Z\sim N(0,1)$, a three-point Gauss–Hermite rule for a standard normal expectation uses nodes $0,\pm\sqrt{3}$ and weights $2/3,1/6,1/6$. For a sufficiently well-behaved function $g$, it gives $$ E\{g(Z)\}=\int_{-\infty}^{\infty}g(z)\frac{e^{-z^2/2}}{\sqrt{2\pi}}\,dz \ \approx\ \frac{2}{3}g(0)+\frac{1}{6}g(-\sqrt{3})+\frac{1}{6}g(\sqrt{3}). $$

For $Y\sim N(\mu,\sigma^2)$, transform the nodes to $\mu$ and $\mu\pm\sqrt{3}\sigma$. If the missing outcomes have a multivariate normal conditional distribution, a linear transformation from independent standard normal variables can similarly turn the expectation into an integral over standard normal coordinates. A tensor-product rule applies a one-dimensional rule in each coordinate. With $K$ nodes per coordinate and $d$ missing continuous outcomes, it requires $K^d$ evaluations before any model-specific simplification.

Worked Example: One Missing Outcome With Logistic Selection

Suppose, hypothetically, that $Y\sim N(0,1)$, and that the probability an outcome is missing given its value is $q(y)=\operatorname{expit}(-0.5-y)=1/\{1+\exp(0.5+y)\}$. For a subject with no observed outcome, the likelihood contribution for the missingness pattern is $$ L_{\mathrm{miss}}=\int q(y)\frac{e^{-y^2/2}}{\sqrt{2\pi}}\,dy=E\{q(Y)\}. $$

The three-point rule has nodes $-\sqrt{3}\approx-1.732$, $0$, and $1.732$. The corresponding probabilities are $q(-1.732)\approx0.7742$, $q(0)\approx0.3775$, and $q(1.732)\approx0.09692$. Substitution gives $$ L_{\mathrm{miss}}\approx\frac{1}{6}(0.7742)+\frac{2}{3}(0.3775)+\frac{1}{6}(0.09692) =0.1290+0.2517+0.01615\approx0.3969. $$

The result is an approximation to the probability of this response pattern under the stated outcome and selection models, not a probability that the model is correct. The assumption $Y\sim N(0,1)$ and the logistic missingness probability define the example; the quadrature evaluates their implied integral.

Integration Dimension and Model Structure

Each additional missing continuous outcome generally adds an integration coordinate. A tensor-product rule with 5 nodes per coordinate requires 5 evaluations in one dimension, 25 in two dimensions, and 3,125 in five dimensions. This growth is exponential in dimension. The number of subjects is not the integration dimension: it is the number of missing components within a subject's contribution. Different response patterns can therefore have different-dimensional integrals.

The raw count is not the whole story. Conditional independence or other factorization may reduce the computation, while dependence among outcomes or a response probability involving their joint values may prevent such a reduction. It is important to retain the dependence specified by the model rather than factor an integral merely to make it easier.

Worked Example: Two Missing Outcomes and a Factorized Integral

Suppose, hypothetically, that $Y_1$ and $Y_2$ are independent standard normal outcomes. Given their values, their response indicators are conditionally independent, and the probability each outcome is missing is $q(y)=\operatorname{expit}(-0.5-y)$. For a subject with both outcomes missing, the likelihood contribution is $$ L_{00}=\int\!\!\int \phi(y_1)\phi(y_2)q(y_1)q(y_2)\,dy_1dy_2 =\{E[q(Y)]\}^2, $$ where $\phi$ is the standard normal density. Independence of the outcomes and conditional independence of the response indicators justify the factorization.

Using the one-dimensional approximation from the previous example, $E[q(Y)]\approx0.3969$. Hence $$ L_{00}\approx(0.3969)^2\approx0.1575. $$ Equivalently, a three-node tensor-product rule has nine pairs of nodes and weights equal to the products of the one-dimensional weights. Here the factorization shows why the nine-term sum equals the square of the three-term sum, apart from rounding.

If the outcomes were correlated, or the response probability for the joint pattern did not factor into separate functions of $y_1$ and $y_2$, this shortcut would not follow. The contribution would still be a two-dimensional integral, but its integrand would have to preserve the joint outcome density and response probability.

Tail Behavior Can Change the Approximation

Quadrature nodes and weights are chosen to approximate integrals, not to guarantee equal accuracy for every integrand. A selection probability can change sharply across the outcome range. In that case, the product of the outcome density and response probability may have important mass in a region that a small set of nodes represents poorly. An integral may be finite and the model well-defined while a particular numerical approximation is inadequate.

The following hypothetical example illustrates the issue using a probit selection probability. It also provides an exact benchmark for this particular integral, so the approximation error can be seen directly.

Worked Example: A Steep Probit Selection Curve

Let $Y\sim N(0,1)$, and suppose $P(R=1\mid Y=y)=\Phi(-2+3y)$, where $\Phi$ is the standard normal cumulative distribution function and $R=1$ means the outcome is observed. The contribution for a missing outcome is $$ L_{\mathrm{miss}}=E\{1-\Phi(-2+3Y)\}. $$

For this special case, introduce an independent $U\sim N(0,1)$. Conditional on $Y$, $\Phi(a+bY)=P(U\leq a+bY\mid Y)$. Averaging over $Y$ gives $P(U-bY\leq a)$. Since $U-bY\sim N(0,1+b^2)$, it follows that $E\{\Phi(a+bY)\}=\Phi\{a/\sqrt{1+b^2}\}$. With $a=-2$ and $b=3$, $$ E\{\Phi(-2+3Y)\}=\Phi\left(\frac{-2}{\sqrt{10}}\right)\approx\Phi(-0.6325)\approx0.2635, \qquad L_{\mathrm{miss}}\approx1-0.2635=0.7365. $$

Now apply the three-point rule to the observed-response probability. At the nodes $-1.732,0,1.732$, the selection probabilities are approximately $0,0.02275,0.9993$. Thus the rule gives $$ E\{\Phi(-2+3Y)\}\approx\frac{1}{6}(0)+\frac{2}{3}(0.02275)+\frac{1}{6}(0.9993) =0.01517+0.1666=0.1817. $$ It consequently approximates the missing-pattern contribution by $1-0.1817=0.8183$, rather than the benchmark $0.7365$. The sharp transition in the selection curve is poorly represented by these few nodes. This example does not imply that three-point quadrature is generally unreliable; it shows that node count alone does not establish accuracy for a particular integrand.

Computation, Diagnostics, and Alternatives

In likelihood estimation, the integral is evaluated repeatedly as $\theta$ and $\gamma$ vary. An approximation error that is nearly constant across parameter values may have little effect on the maximizing parameter, but parameter-dependent error can alter likelihood differences, scores, and curvature. It can therefore change the apparent location or shape of a maximum. This is numerical error, distinct from sampling uncertainty and from uncertainty about the model assumptions.

Warning signs include likelihood values that change noticeably when the integration setup is modified, unstable contributions for subjects with many missing outcomes, or a fitted likelihood that behaves irregularly as parameters change. These signs do not by themselves identify the cause: a difficult likelihood surface and a poor integral approximation can coexist. A useful first step is to inspect the contribution and response pattern responsible for the most demanding integral, including whether its response probability is concentrated in an extreme part of the outcome distribution.

The main alternatives depend on the integral's structure and dimension:

Analytic integration is preferable when a closed form is available, as in the probit-normal identity above. It avoids quadrature error for that integral, but such formulas are model-specific.
Gaussian quadrature is a deterministic weighted-sum approximation that can be effective for low-dimensional smooth integrals. Tensor products become costly as dimension grows.
Monte Carlo integration approximates an expectation by an average over simulated draws. It can be more practical in higher dimensions, but introduces simulation variability and may require many draws for stable likelihood comparisons.
Model-based factorization can reduce dimension when independence or conditional structure genuinely makes it valid. It is not a numerical trick: it follows from the stated probability model.

No method removes the need to specify a plausible outcome model and selection mechanism. Nor does a numerically precise integral establish that the assumed dependence structure or response probabilities are appropriate. Numerical integration addresses how to evaluate a stated likelihood, not whether that likelihood is scientifically adequate.

Common pitfall: Treating a missing outcome as if its response probability could be evaluated at a single “typical” value, such as the outcome-model mean. In general, $E\{h(r\mid Y)\}$ is not equal to $h(r\mid E[Y])$. The correct contribution averages the response probability over the outcome distribution, with the density and any observed outcomes included as specified by the likelihood.

Check Your Understanding

Answer each question, then compare with the brief solutions.

A subject has two missing continuous outcomes. What determines the dimension of that subject's likelihood integral?
With four quadrature nodes per coordinate, how many tensor-product evaluations are required for a three-dimensional integral?
For the two-outcome example above, why could the integral be squared as a one-dimensional expectation?
Why can a steep selection curve cause trouble for a low-node quadrature approximation even when the likelihood integral exists?

Solutions. (1) The number of missing continuous components being integrated over: here it is two, unless valid model structure reduces the integral. (2) $4^3=64$ evaluations. (3) The outcomes are independent and the conditional probability of both being missing factors into $q(y_1)q(y_2)$, so the joint expectation factors. (4) A small set of nodes may not represent the region where the density-weighted selection probability changes rapidly or contributes substantially; existence of the integral does not guarantee accuracy of that approximation.

Back to Tutorials