Multivariate Sample Selection Model

• The multivariate sample selection model is defined as a framework to correct for endogenous observability by modeling multiple selection patterns and correlated outcomes.
• Estimation methods, including plug-in and double machine learning, leverage efficient influence functions to reduce bias and achieve robust, root-n consistent inference.
• Extensions to unordered categorical outcomes and matrix-variate systems provide practical tools for addressing complex empirical applications with multiple sources of selection.

A multivariate sample selection model is a class of models in which the outcome, or a collection of outcomes, is observed only under a nontrivial selection mechanism that is itself multivariate in structure. In the recent literature, this multivariate character appears in several distinct but related forms: the outcome may be observed only when a multivariate selection pattern is realized; treatment assignment and outcome observability may jointly induce a “double selection” problem; the latent outcome may be categorical with selection-specific dependence that varies by category; or several correlated outcomes may each be subject to their own selection equation within a matrix-variate system (Dolgikh et al., 16 Nov 2025, Bia et al., 2020, Boussim, 7 Oct 2025, Lim et al., 3 May 2026). These formulations share a common objective: identification and estimation of substantively meaningful parameters when observability is endogenous rather than incidental.

1. Conceptual scope and basic structure

The classical reference point is the Heckman sample selection model, which specifies an outcome equation

$Y_{1i}=x_i^\top \beta+\epsilon_{1i},$

a latent selection equation

$Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$

and a selection indicator

$C_i=\mathbb{I}(Y_{2i}>0).$

The outcome is observed only when $C_i=1$, and under bivariate normal errors the conditional mean of the observed outcome contains the inverse Mills ratio correction,

$$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$

This is the standard sample selection bias correction and the scalar benchmark from which multivariate extensions proceed (Lim et al., 3 May 2026).

In one multivariate formulation, the outcome is only observed when a certain multivariate selection pattern is realized. The recent causal-inference literature describes this as a multivariate sample selection model with ordinal selection equations, and the conclusion of the 2025 paper on double machine learning indicates that selection can be encoded by a binary vector $z^{(t)}$, where each coordinate indicates whether the outcome is observed in a given period. The same source describes the framework as generalizing standard sample selection by allowing selection to depend on multiple components (Dolgikh et al., 16 Nov 2025).

A second formulation treats multivariate selection as arising from multiple sources of selectivity rather than from multiple observed outcomes. In double machine learning for sample selection models, the observed data are $(Y\cdot S, D, S, X)$ in the MAR case, and $(Y\cdot S, D, S, X, Z)$ or $(Y\cdot S, D, S, X, M)$ in IV or dynamic-confounding cases. Here treatment assignment and outcome observability are both selective, so the model is explicitly described as a “double selection” problem (Bia et al., 2020).

A third formulation extends selection correction to several correlated outcomes observed under possibly different selection processes. The matrix-variate Multiple Heckman Selection Model arranges, for each subject, the latent outcome row and the latent selection row into a $2\times R$ matrix and uses a matrix normal distribution to capture dependence both within the outcome/selection block and across outcomes (Lim et al., 3 May 2026).

These formulations define a common field but not a single canonical specification. This suggests that “multivariate” in sample selection models operates along at least three axes: multivariate selection patterns, multiple sources of selection, and multiple correlated outcomes.

Formulation	Observability structure	Key feature
MSSM with ordinal selection equations	Outcome observed only under a multivariate selection pattern	Causal estimands such as ATE, ATET, and LATE
DML sample selection model	Outcome observed only when $Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$0, with selective treatment assignment and observability	“Double selection” and orthogonal scores
Categorical-outcome selection model	Latent unordered category observed only if selected	Local Logistic Representation
Multiple Heckman Selection Model	Multiple correlated outcomes, each with its own selection equation	Matrix normal and SUN correction

2. Identification under multivariate selection

Identification in multivariate sample selection models depends on restrictions that separate the latent outcome mechanism from the selection mechanism. In the 2025 causal MSSM paper, the conclusion explicitly highlights the importance of an exclusion restriction assumption and notes robustness to violations of it as a topic for future work. The same conclusion also states that, for nested sample selection in the estimation of ATE and ATET, one can replace the paper’s main assumption 1 with a multivariate probit model with sample selection, indicating that the baseline identification result is not simply a direct restatement of multivariate probit (Dolgikh et al., 16 Nov 2025).

In the DML treatment-evaluation framework, identification is decomposed into assumptions for treatment assignment and assumptions for selection. For treatment selection, the model assumes selection on observables,

$Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$1

For sample selection, it considers either selection on observables / MAR,

$Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$2

or nonignorable selection with an instrument $Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$3, combined with a threshold-crossing model

$Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$4

The dynamic-confounding extension replaces static conditioning with the sequential conditional independence restriction

$Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$5

allowing covariates that jointly affect selection and the outcome to be influenced by treatment (Bia et al., 2020).

For unordered categorical outcomes, identification is developed through a Local Logistic Representation rather than through latent thresholding. The latent outcome is

$Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$6

selection is driven by latent selection variables $Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$7 and a binary instrument $Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$8, and the outcome is observed only when $Y_{2i}=w_i^\top \gamma+\epsilon_{2i},$9. The central decomposition writes

$C_i=\mathbb{I}(Y_{2i}>0).$0

where $C_i=\mathbb{I}(Y_{2i}>0).$1 is a marginal outcome parameter, $C_i=\mathbb{I}(Y_{2i}>0).$2 is a marginal selection parameter, and $C_i=\mathbb{I}(Y_{2i}>0).$3 is a category-specific local association parameter. Without restrictions, the model is underidentified. The key identifying restriction is

$C_i=\mathbb{I}(Y_{2i}>0).$4

which says that the way selection sorts people across categories does not depend on the instrument. Under the IV assumptions and the interiority condition $C_i=\mathbb{I}(Y_{2i}>0).$5, Theorem 1 states that $C_i=\mathbb{I}(Y_{2i}>0).$6 is point-identified as the unique interior solution to the nonlinear system (Boussim, 7 Oct 2025).

Across these literatures, the common logic is that selection probabilities must vary in a way that is informative about observability, but the dependence structure linking selection and the latent outcome must remain sufficiently stable or sufficiently structured for the latent object of interest to be recoverable.

3. Estimation: plug-in, influence functions, and double machine learning

Recent work on multivariate sample selection has moved from correction formulas toward semiparametric and machine-learning-based estimators. The 2025 MSSM paper proposes both plug-in (PI) and double machine learning (DML) estimators of ATE, ATET, and LATE, and states that the DML estimators are doubly-robust and based on the efficient influence functions. The same paper studies finite sample properties on simulated data and reports that, without addressing multivariate sample selection, estimates of the causal parameters may be highly biased, whereas the proposed estimators allow these biases to be avoided (Dolgikh et al., 16 Nov 2025).

In the PI approach, the estimand is constructed by first estimating the needed conditional means, probabilities, or nonparametric regression functions and then substituting those estimates into the identified formula. In the DML approach, estimation proceeds through efficient influence-function-based scores. This suggests a shift from direct substitution toward orthogonal score construction, with robustness to regularization bias provided by semiparametric orthogonality.

The 2020 DML sample-selection paper gives the score structure explicitly. Under MAR, the mean potential outcome is identified by

$C_i=\mathbb{I}(Y_{2i}>0).$7

with

$C_i=\mathbb{I}(Y_{2i}>0).$8

where

$C_i=\mathbb{I}(Y_{2i}>0).$9

The paper states that this score is equivalent to the regression and IPW representations and is doubly robust: it identifies $C_i=1$0 if either the outcome regression $C_i=1$1 or the propensity components $C_i=1$2 are correctly specified. Analogous orthogonal scores are derived for IV-based selection and for the dynamic-confounding case, where the score also contains the nested regression $C_i=1$3 (Bia et al., 2020).

The same DML framework uses K-fold sample splitting / cross-fitting: nuisance functions are estimated on auxiliary folds, predicted on held-out folds, and inserted into the orthogonal score before aggregation. The paper states that a rate of $C_i=1$4 for the nuisance estimators is sufficient for root-$C_i=1$5 inference, and proves root-$C_i=1$6 consistency and asymptotic normality for the main estimators under bounded moments, overlap, and product-rate conditions (Bia et al., 2020).

The methodological significance is that multivariate sample selection is treated not merely as a correction problem but as a semiparametric estimation problem in which nuisance learning, orthogonality, and efficiency are jointly relevant.

4. Outcome spaces beyond scalar continuous responses

A common misconception is that sample selection models are intrinsically tied to latent continuous outcomes and Gaussian corrections. The categorical-outcome literature shows that this is not the case. The 2025 paper on correcting sample selection bias with categorical outcomes explicitly departs from the usual latent continuous-outcome logic and extends sample selection modeling to unordered categorical outcomes by developing a local representation that applies to joint probabilities, thereby eliminating the need to impose an artificial ordering on categories (Boussim, 7 Oct 2025).

Its core object is the Local Logistic Representation

$C_i=1$7

Lemma 1 states that for any measurable sets $C_i=1$8 and $C_i=1$9 with positive probabilities, there exists a unique local association parameter $$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$0 such that

$$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$1

and the sign of $$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$2 matches the sign of the covariance of the corresponding indicators. In the categorical setting, each joint category-selection probability is decomposed into a marginal outcome term, a marginal selection term, and a category-specific association parameter:

$$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$3

The paper then turns this identification result into a semiparametric multinomial logit model with selection. For each non-baseline category,

$$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$4

the selection equation is

$$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$5

and the category-specific local association is parameterized as

$$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$6

Estimation is carried out by a computationally tractable two-step estimator: first estimate $$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$7 by ordinary logit regression of $$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$8 on $$\mathbb{E}[Y_{1i}\mid C_i=1,x_i,w_i] = x_i^\top \beta+\rho\sigma\,\lambda(w_i^\top\gamma).$$9, then estimate the outcome and sorting parameters through the second-stage objective $z^{(t)}$0. Under standard regularity conditions, the second-stage estimator is consistent and asymptotically normal as a two-step M-estimator with a generated regressor (Boussim, 7 Oct 2025).

This line of work broadens multivariate sample selection in two ways. First, it shows that sample selection correction need not be tied to ordered support. Second, it locates multivariate structure in the category-specific dependence between selection and the latent outcome, rather than only in multivariate continuous response vectors.

5. Multiple outcomes, matrix-normal dependence, and the SUN correction

A distinct branch of the literature extends sample selection from one outcome to many. The matrix-variate Multiple Heckman Selection Model introduces a model for multiple correlated outcomes observed under possibly different selection processes by arranging the latent variables for subject $z^{(t)}$1 as a $z^{(t)}$2 matrix,

$z^{(t)}$3

with outcome equations

$z^{(t)}$4

selection equations

$z^{(t)}$5

and indicators

$z^{(t)}$6

The $z^{(t)}$7th outcome is observed only when the corresponding latent selection variable is positive (Lim et al., 3 May 2026).

The error matrix is assumed matrix normal,

$z^{(t)}$8

where $z^{(t)}$9 captures the dependence between the outcome row and the selection row, and $(Y\cdot S, D, S, X)$0 captures dependence across the $(Y\cdot S, D, S, X)$1 outcomes. The paper fixes the selection variance through

$(Y\cdot S, D, S, X)$2

thereby resolving the matrix-normal scale nonidentifiability (Lim et al., 3 May 2026).

Selection bias is represented through truncation of the latent selection row. For each subject, the selected index set is

$(Y\cdot S, D, S, X)$3

The observed-data likelihood decomposes into a density part for the observed components and a rectangle probability for the censored or latent components. A key theoretical result is that the conditional distribution of the observed outcomes belongs to the multivariate unified skew-normal (SUN) family. The multivariate selection correction is then

$(Y\cdot S, D, S, X)$4

where $(Y\cdot S, D, S, X)$5 is the multivariate analogue of the inverse Mills ratio. When $(Y\cdot S, D, S, X)$6, this reduces to the scalar Heckman correction (Lim et al., 3 May 2026).

Estimation is based on an ECM algorithm with closed-form updates for all model parameters. The E-step computes first and second moments of a truncated multivariate normal distribution, and the CM-steps update the regression and covariance parameters. The paper reports that the approach is implemented in the R package mvHeckman (Lim et al., 3 May 2026).

This multiple-outcome line of research places the multivariate aspect in cross-outcome dependence itself. It is therefore different from the causal MSSM formulation, in which multivariate structure is attached to the selection pattern, and different again from the DML formulation, in which multivariate structure appears through multiple selection mechanisms and high-dimensional nuisance structure.

6. Empirical implications, adjacent frameworks, and open directions

The empirical implications of multivariate sample selection are consistent across the recent literature: selection matters, and ignoring its multivariate form can distort both point estimation and uncertainty quantification. In the 2025 MSSM paper, the main simulation takeaway is explicit: without addressing multivariate sample selection, the estimates of the causal parameters may be highly biased, whereas the proposed PI and DML estimators avoid these biases. The same conclusion also points to future extensions to ATEG, nested sample selection, and robustness to violations of the exclusion restriction (Dolgikh et al., 16 Nov 2025).

Applications span several data types. In the DML sample-selection paper, the Job Corps application evaluates the effect of academic training or vocational training on hourly wages that are only observed conditional on employment. The paper estimates short-run effects one year after assignment using MAR and IV approaches and longer-run effects four years after assignment using the dynamic-confounding approach with post-treatment covariates $(Y\cdot S, D, S, X)$7 (Bia et al., 2020). In the categorical-outcome paper, the proposed framework is motivated by occupational choice, marital status, political affiliation, educational field, health status, and survey response categories, and the empirical illustration studies health care utilization in Côte d’Ivoire, where provider choice is observed only for people who seek care (Boussim, 7 Oct 2025). In the matrix-variate Heckman paper, real-data illustrations include the Mroz labor supply data with joint outcomes log wage and log hours worked, and NHANES blood pressure data with systolic and diastolic blood pressure. The paper states that ignoring cross-outcome dependence can distort coefficient estimates and uncertainty quantification (Lim et al., 3 May 2026).

An important interpretive boundary is provided by the conformal-inference literature. “Multivariate Conformal Selection” is described explicitly as not a classical parametric statistical selection model and not a latent-variable sample selection model in the econometric sense. Instead, it is a conformal prediction / multiple-testing-based inferential selection framework built from multivariate nonconformity scores, conformal p-values, regional monotonicity, and Benjamini–Hochberg FDR control (Bai et al., 1 May 2025). This distinction matters because the phrase “multivariate selection” is used in more than one technical tradition. In econometric sample selection models, the object is the mechanism that governs observation or attrition; in conformal selection, the object is a selection rule with finite-sample FDR control.

Taken together, the current literature shows that multivariate sample selection is not a single model but a family of identification and estimation problems generated by endogenous observability under richer dependence structures than the scalar Heckman benchmark. The main directions already visible in recent work are flexible semiparametric identification, orthogonal-score-based causal estimation, categorical and unordered outcomes, matrix-variate dependence, and extensions to nested sample selection and endogenous groups.