Sample Selection Models for Panel and Cross-Sectional Data

1.1. Supported models

Table 1 summarizes the available models. Panels may be unbalanced, and participation may change across occasions.

Note: RE means random effects; — indicates an absent correlation channel. Standalone probit, Poisson and PLN random-effects models are also supported and provide default outcome-coefficient starting values for the corresponding panel-selection models.

1.2. What changed since 1.0.1?

Version 2.0.0 is a full refactor of the package, developed with assistance from Codex and Claude Code. We recommend 2.0.0 for new analyses and for refitting analyses conducted with 1.0.1. The four original panel-selection families are retained, with improvements that matter for everyday use:

• Corrected estimation. Corrections to likelihood and score calculations address some errors that can slightly distort parameter estimates and standard errors.
• A unified, user-friendly interface. panel_select() replaces separate model-specific fitting functions, with common controls and standard methods such as summary(), predict() and update().
• Broader model support. Cross-sectional selection models and standalone binary and count panel models complement the original panel-selection families.
• More control over numerical accuracy. Adaptive Gauss–Hermite (AGH) and Gauss–Legendre (GL) supplement fixed Gauss–Hermite (GH). Together with explicit numerical checks, they help detect and address inaccurate integration in difficult cases.
• Faster fitting in general. The numerical engine has been optimized for better performance.

The migration vignette gives the main changes needed to update existing code.

2.1. Continuous Outcome (Probit RE–Linear RE)

The selection and outcome equations are

\[ d_{it}=1\bigl(\mathbf{w}_{it}\boldsymbol{\alpha}+\delta u_i+\varepsilon_{it}>0\bigr), \]

\[ y^*_{it}=\mathbf{x}_{it}\boldsymbol{\beta}+\lambda v_i+\sigma\epsilon_{it}. \]

The correlation structure is

\[ \begin{pmatrix}u_i\\v_i\end{pmatrix} \sim\mathcal{N}_2\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\right). \]

\[ \begin{pmatrix}\varepsilon_{it}\\\epsilon_{it}\end{pmatrix} \sim\mathcal{N}_2\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\tau\\\tau&1\end{pmatrix}\right). \]

Use this model for continuous outcomes with conditionally normal errors. The scales \(\delta\), \(\lambda\) and \(\sigma\) are the SDs of selection effects, outcome effects and outcome errors. The correlations \(\rho\) and \(\tau\) capture persistent and same-occasion dependence. The selection shock has variance one.

2.2. Binary Outcome (Probit RE–Probit RE)

For a binary potential outcome,

\[ d_{it}=1\bigl(\mathbf{w}_{it}\boldsymbol{\alpha}+\delta u_i+\varepsilon_{it}>0\bigr), \]

\[ y^*_{it}=1\bigl(\mathbf{x}_{it}\boldsymbol{\beta}+\lambda v_i+\epsilon_{it}>0\bigr), \]

with

\[ \begin{pmatrix}u_i\\v_i\end{pmatrix} \sim\mathcal{N}_2\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\right). \]

\[ \begin{pmatrix}\varepsilon_{it}\\\epsilon_{it}\end{pmatrix} \sim\mathcal{N}_2\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\tau\\\tau&1\end{pmatrix}\right). \]

Use this model for binary outcomes. Both occasion-shock variances are normalized to one.

2.3. Poisson Count Outcome (Probit RE–Poisson RE)

For count outcomes with selection only through individual effects,

\[ d_{it}=1\bigl(\mathbf{w}_{it}\boldsymbol{\alpha}+\delta u_i+\varepsilon_{it}>0\bigr), \]

\[ y^*_{it}\mid\mathbf{x}_{it},v_i \sim\operatorname{Poisson}\!\left\{ \exp(\mathbf{x}_{it}\boldsymbol{\beta}+\lambda v_i)\right\}, \]

where

\[ \begin{pmatrix}u_i\\v_i\end{pmatrix} \sim\mathcal{N}_2\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\right). \]

\[ \varepsilon_{it}\sim\mathcal{N}(0,1). \]

The outcome random effect allows for overdispersion arising from unobserved individual heterogeneity. Conditional on the individual random effects, selection and the Poisson outcome are independent, so the model has no separate occasion-specific correlation parameter \(\tau\).

2.4. Poisson-Lognormal Count Outcome (Probit RE–PLN RE)

The Poisson-lognormal model extends the Poisson specification by adding an observation-level random effect to the log rate. This effect can be correlated with the selection shock \(\varepsilon_{it}\), allowing for same-occasion selection dependence. The equations are:

\[ d_{it}=1\bigl(\mathbf{w}_{it}\boldsymbol{\alpha}+\delta u_i+\varepsilon_{it}>0\bigr), \]

\[ y^*_{it}\mid\mathbf{x}_{it},v_i,\epsilon_{it} \sim\operatorname{Poisson}\!\left\{ \exp(\mathbf{x}_{it}\boldsymbol{\beta}+\lambda v_i+\gamma\epsilon_{it})\right\}, \]

with

\[ \begin{pmatrix}u_i\\v_i\end{pmatrix} \sim\mathcal{N}_2\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\right). \]

\[ \begin{pmatrix}\varepsilon_{it}\\\epsilon_{it}\end{pmatrix} \sim\mathcal{N}_2\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\tau\\\tau&1\end{pmatrix}\right). \]

Here \(\gamma\) measures the standard deviation of the additional random effect at the observation level. PLN accommodates overdispersion beyond that captured by individual effects and allows selection and outcome shocks to be correlated within occasions.

3.1. Simulate a panel with known parameters

The example follows the Linear specification. Both equations contain the same two covariates, with coefficients one, and all three SDs equal one. The correlations are \(\rho=-0.5\) and \(\tau=0.5\):

\[ d_{it}=1\bigl(x_{it}+w_{it}+u_i+\varepsilon_{it}>0\bigr), \]

\[ y^*_{it}=x_{it}+w_{it}+v_i+\epsilon_{it}. \]

We simulate 500 people over five occasions. Both covariates are independent standard normal variables, independent of the latent disturbances. The code retains unselected rows and records their outcomes as NA.

set.seed(100)
N = 500L
period = 5L
n = N*period
id = rep(1:N, each = period)
u = rnorm(N)
v = -0.5*u + sqrt(0.75)*rnorm(N)
eps = rnorm(n)
e = 0.5*eps + sqrt(0.75)*rnorm(n)
w = rnorm(n)
x = rnorm(n)
d = as.integer(x + w + u[id] + eps > 0)
y = x + w + v[id] + e
y[d == 0] = NA
dat = data.frame(id = id, time = rep(seq_len(period), N), x = x, w = w,
                 d = d, y = y)

3.2. Fit the model and compare estimates

Fit the equations jointly using 10-node GH integration.

fit = panel_select(
  outcome = y ~ x + w,
  selection = d ~ x + w,
  data = dat,
  id = "id",
  family = "linear",
  control = panel_control(quadrature = "gh", nodes = 10L))

Table 2 reports precomputed estimates from GH10, AGH10 and GL10 alongside the known parameters. The AGH and GL fits are not rerun when this vignette is built. For AGH or GL, keep nodes = 10L and change quadrature to "agh" or "gl", respectively. Standard errors use the default OPG covariance.

The three integration methods produce the same parameter estimates and standard errors to the reported precision. The estimated slopes and standard deviations are close to their true values, while the correlation estimates recover the opposing selection channels: negative dependence between individual effects and positive dependence between occasion-specific shocks.

3.3. Check and refine the fit

Step 1. Check convergence and stationarity. The optimizer should report convergence, and Q_OPG should be close to zero. The package warns when Q_OPG > 0.1 by default.

check_stationarity(fit)[c("optimizer_converged", "Q_OPG")]
#> $optimizer_converged
#> [1] TRUE
#> 
#> $Q_OPG
#> [1] 5.087207e-07

A small Q_OPG indicates that the fitted numerical objective is close to stationary. Integration accuracy requires a separate check.

Step 2. Compare calculations with more nodes. Re-evaluate the 10-node fit using 20 nodes while holding the parameter estimates fixed:

check_quadrature(fit, nodes = 20L)
#> Quadrature check (higher_order) of: Fixed Gauss-Hermite (gh), 10 x 10 nodes (selection x outcome)
#>           rule nodes loglik_change max_panel_change max_score_change
#>  person_effect    20       0.00131          0.00756           0.0329
#>  max_panel_score_change within_tolerance
#>                    0.06            FALSE
#> Tolerances: loglik 1e-05, score 1e-04; optimizer converged; covariance valid

Although the differences between the 10-node and 20-node GH calculations exceed the default comparison tolerances, refitting with 20 nodes changes the parameter estimates only slightly. To refine the fit, refit with 20 nodes and repeat the check using 30 nodes. If notable sensitivity persists, consider AGH or GL. Because GL applies the specified number of nodes within each integration segment, GL10 generally requires substantially more computation than GH10 or AGH10. GL10 or GL12 may be sufficient for many applications.

Step 3. Refit with more nodes if needed. Use update() to fit at the higher order:

fit_more = update(fit, control = modifyList(fit$control, list(nodes = 20L)))
coef(fit_more) - coef(fit)

update() uses the fitted coefficients as starting values and preserves the other controls. Compare the estimates, then repeat the checks on fit_more with a higher comparison order.