Package {PanelSelect}

PanelSelect: panel sample-selection models

Description

PanelSelect provides a unified and flexible framework for sample selection models with continuous, binary, and count outcomes in panel and cross-sectional data. In particular, it supports a broad family of panel selection models, which extend the classic Heckman framework to account for persistent unobserved heterogeneity. The package supports selection at two levels: the individual level, through correlated random effects, and the observation level, through correlated occasion-specific shocks. This allows researchers to distinguish non-random selection driven by persistent individual characteristics from selection driven by temporary shocks. Alongside model estimation, PanelSelect supports prediction of population-level outcomes.

Main functions

• panel_select() fits a model from an outcome formula, a selection formula (or NULL for a standalone model), the data, the panel ID column (or NULL for cross-sectional data) and the outcome family; its arguments also set restrictions, starting values, missing-data handling and the reported covariance. panel_control() sets the integration rules and the optimizer.

• Standard methods (panelselect-methods) summarize, compare and update fits; predict.panelselect(), residuals.panelselect() and random_effects() give predictions and posterior summaries.

• check_quadrature() compares likelihoods and scores at retained estimates under higher orders or an independent GL rule; agreement is evidence of stability, not an integration certificate.

• check_stationarity() reads stored scores and information, reporting Q_OPG and, with valid positive-definite stored observed information and its inverse, Q_Hessian. The default warning is strictly Q_OPG > 0.1. It does not compute a Hessian or alter backend convergence.

• add_information() explicitly returns a fit with observed information stored for reuse. Requested covariance can remain unavailable; optimizer convergence, stationarity, integration and inference require distinct checks.

Use panel_select() in place of the PanelSelect 1.x and PanelCount wrappers; vignette("migration", package = "PanelSelect") maps their calls and results to panel_select().

Quick start

fit = panel_select(y ~ x, observed ~ .outcome + z, data, id = "person",
                   family = "linear")
summary(fit)
check_stationarity(fit)[c("Q_OPG", "Q_Hessian")]
predict(fit, newdata, type = "selected_mean")
panel_select(count ~ x, NULL, data, "person", "poisson")  # no selection

Start with vignette("vignette", package = "PanelSelect"); users of earlier versions should read vignette("migration", package = "PanelSelect"). GL supports population predictions; history-conditioned predictions and random-effect summaries require a GH/AGH fit, with no prediction-time integration override. History-conditioned prediction requires explicit history; random_effects() defaults to the retained fitting rows.

Author(s)

Maintainer: Jing Peng jing.peng@uconn.edu

Authors:

• Jing Peng jing.peng@uconn.edu

Compute and store observed information for a fitted model

Description

Augment an OPG fit with the observed Hessian and its covariance estimators, without refitting. Assign the returned object to retain the calculation for later inference and stationarity diagnostics.

Usage

add_information(object)

Arguments

Details

The default vcov = "opg" fit does not compute a post-fit observed Hessian. add_information() differentiates the full score at the fitted free optimizer coordinates to obtain the negative log-likelihood Hessian. This is the observed information of the fitted numerical objective; it does not certify quadrature accuracy, a global maximum or calibrated intervals.

The returned fit stores the information, any integration-accuracy record, its computed/unavailable status and error reason, and the observed, OPG and sandwich covariances on both optimizer and natural parameter scales. Its coefficients, likelihood and default covariance choice are unchanged. check_stationarity() can then report Q_Hessian and Newton diagnostics when the stored information and its inverse are finite and positive definite. Missing or invalid curvature remains unavailable; no OPG or sandwich matrix substitutes for the observed inverse.

Explicit observed/sandwich requests in vcov(), confint(), summary() or predictions can calculate missing information for that request. They do not store it in the caller's original object. Use fit = add_information(fit) once to reuse the stored calculation in subsequent requests.

A fit whose information is already marked as computed is returned unchanged, including when its stored information is unavailable. This also applies to information stored by older versions, including before the Gaussian outcome-unit correction. Use update(fit, vcov = "observed") to refit and recompute information with current numerical steps.

A staged GL fit whose retained mesh cannot serve the information stencil warns and stores an unavailable result. It does not silently refine the mesh or replace the estimator. An explicit refit with vcov = "observed" prepares a mesh for the information calculation, subject to the requested accuracy and resource limits.

Value

A panelselect fit with stored information and covariance results, or the unchanged fit when information was already computed. Assign this return value to preserve the added information.

See Also

panelselect-methods

Other diagnostic functions: check_quadrature(), check_stationarity()

Examples

set.seed(1)
dat = data.frame(person = rep(1:60, each = 4), x = rnorm(240), z = rnorm(240))
u = rnorm(60); v = 0.5 * u + rnorm(60, sd = 0.8)  # correlated person effects
dat$observed = as.integer(0.3 + 0.5 * dat$z + 0.7 * u[dat$person] + rnorm(240) > 0)
dat$y = ifelse(dat$observed == 1, 1 + 0.5 * dat$x + v[dat$person] + rnorm(240), NA)
fit = panel_select(y ~ x, observed ~ x + z, dat, "person", "linear")
fit = add_information(fit)
sqrt(diag(vcov(fit, "observed")))
check_stationarity(fit)$Q_Hessian

Check integration sensitivity at the fitted estimates

Description

Re-evaluate a fitted model using higher quadrature orders or an alternative GL rule. The comparisons hold parameters fixed, so they isolate changes in the numerical likelihood and score. A higher order is not guaranteed to be more accurate, and agreement is evidence of stability rather than a certificate.

Usage

check_quadrature(
  object,
  nodes = NULL,
  observation_nodes = NULL,
  tolerance = c(loglik = 1e-05, score = 1e-04),
  refit = FALSE,
  compare = c("higher_order", "gl"),
  quadrature_options = list(),
  ...
)

## S3 method for class 'panelselect_quadrature_check'
print(x, digits = 3L, ...)

Arguments

Value

A panelselect_quadrature_check list with compare, fitted (a fitted-rule description), the comparison data frame, tolerance, optimizer_converged, covariance_valid, integration, failures, refit (NULL by default), and qualification. print() displays the main comparison columns and returns the object invisibly. Covariance validity describes the requested fitted estimator, not statistical calibration.

Choose a comparison

For PLN fits, person_effect changes only the outer integral, observation only the occasion-level integral, and both changes both where both are active. Fixed restrictions and exact reductions can remove an integral, so only effective changes produce rows. A fully closed-form model has no numerical integral to refine and the function stops with that explanation.

Read the results and decide what to do

• If discrepancies exceed tolerance, investigate larger orders or another method and assess changes in estimates/uncertainty after an explicit refit.

• If integration fails, inspect message, ⁠$failures⁠ and ⁠$integration⁠; change only the named controls whose limits or accuracy require attention.

• Even small changes do not prove both rules accurate. Parameter-dependent integration error can move an optimum despite correct finite-rule scores. Compare competing endpoints under a common checked rule when needed.

• Keep optimizer status, check_stationarity(), covariance availability, and substantive sensitivity separate from this integration comparison.

See Also

panel_control(), check_stationarity(), update.panelselect()

Other diagnostic functions: add_information(), check_stationarity()

Examples

set.seed(1)
dat = data.frame(person = rep(1:60, each = 4), x = rnorm(240), z = rnorm(240))
u = rnorm(60); v = 0.5 * u + rnorm(60, sd = 0.8)  # correlated person effects
dat$observed = as.integer(0.3 + 0.5 * dat$z + 0.7 * u[dat$person] + rnorm(240) > 0)
dat$y = ifelse(dat$observed == 1, 1 + 0.5 * dat$x + v[dat$person] + rnorm(240), NA)
fit = panel_select(y ~ x, observed ~ x + z, dat, "person", "linear")
check_quadrature(fit)
check_quadrature(fit, nodes = 20, tolerance = c(loglik = 1e-6, score = 1e-5))

Inspect score and curvature diagnostics without refitting

Description

Calculate stationarity diagnostics from the fitted object's stored scores and covariances. The primary screen is Q_OPG; when valid observed information is stored, Q_Hessian and local Newton diagnostics are also available. No likelihood, score or Hessian is evaluated by this function.

Usage

check_stationarity(object)

Arguments

Value

A list with top-level numeric Q_OPG and Q_Hessian (NA when unavailable), backend convergence/code/status, named optimizer score, the curvature/Newton quantities above, information_error, and retained rejected-trial summaries. opg contains Q_OPG, status, tolerance, warning, reason, and basis. Top-level status and reason concern the observed-information diagnostics; they can be unavailable while opg$status is available. The supplied fit is unchanged.

Diagnostic quantities

Let g be the total free-parameter score in optimizer coordinates, B the sum of independent-unit score outer products, and A=-\nabla^2\ell the observed information for the fitted numerical objective.

Optimizer coordinates are coefficients, log-SDs and Fisher-transformed correlations. The Hessian and Newton interpretation is local to these coordinates; it is not generally invariant to nonlinear reparameterization away from a stationary point. Neither Q is a p-value or a global-optimality test.

Warning policy and unavailable results

• The default warning cutoff is Q_OPG > 0.1, strictly. Change panel_control(stationarity_warning = ...) to use another Q cutoff; NULL disables the warning while preserving the statistic. This is a diagnostic policy, not an optimizer stopping criterion.

• There is no additional automatic cutoff for Q_Hessian.

• Missing, nonfinite, singular or indefinite observed information/inverse leaves Q_Hessian and the Newton quantities unavailable with a reason. OPG and sandwich covariance are never substituted for observed information.

• An all-fixed model has status "all_fixed" and no free-parameter stationarity quantities. An unavailable OPG screen is not an automatic pass.

• Native saved fits use current Q semantics and their stored scores/covariances. Previously stored Gaussian Hessians remain historical calculations; an explicit update/refit obtains corrected outcome-scale information.

A small Q, or a backend convergence code, does not certify integration accuracy, identification, calibrated intervals or a global maximum. Near parameter boundaries the screen can overstate remaining displacement. Use check_quadrature() for integration sensitivity, and explicit tighter optimizer settings or different starts when investigating optimization.

See Also

panel_control(), add_information(), check_quadrature()

Other diagnostic functions: add_information(), check_quadrature()

Examples

set.seed(1)
dat = data.frame(person = rep(1:60, each = 4), x = rnorm(240), z = rnorm(240))
u = rnorm(60); v = 0.5 * u + rnorm(60, sd = 0.8)  # correlated person effects
dat$observed = as.integer(0.3 + 0.5 * dat$z + 0.7 * u[dat$person] + rnorm(240) > 0)
dat$y = ifelse(dat$observed == 1, 1 + 0.5 * dat$x + v[dat$person] + rnorm(240), NA)
fit = panel_select(y ~ x, observed ~ x + z, dat, "person", "linear")
check_stationarity(fit)$opg
check_stationarity(add_information(fit))$max_standardized_step

Extract coefficients, covariance matrices and confidence intervals

Description

Inspect the natural-scale estimates of a panel_select() fit, select a covariance estimator, or obtain parameter confidence intervals. These methods do not refit the model. Observed and sandwich covariance can require an additional information calculation at the fitted estimates.

Usage

## S3 method for class 'panelselect'
coef(object, ...)

## S3 method for class 'panelselect'
vcov(object, type, vcov, ...)

## S3 method for class 'panelselect'
confint(object, parm, level = 0.95, vcov = object$control$vcov, ...)

Arguments

Value

coef() returns a named numeric vector. vcov() returns a square named natural-scale covariance matrix. confint() returns a two-column matrix of lower and upper limits for the requested parameters.

Choosing a covariance estimator

Let each independent unit contribute a score vector, let B be the sum of their outer products, and let A be the observed information. Independent units are people in a panel model and rows in a cross-sectional model.

Assign fit = add_information(fit) to retain information for repeated observed/sandwich requests. Information that is singular, indefinite, or numerically unavailable is not replaced by a different estimator. Staged GL uses the retained mesh; if information needs a refined mesh, explicitly update/refit with vcov = "observed" or "sandwich".

Parameter scales and intervals

• Regression coefficients are reported on the model's coefficient scale.

• SDs are positive natural SDs; intervals are constructed on the log scale and transformed back. Testing an SD against zero is a boundary problem, so summaries do not supply a usual zero-test p-value for SDs.

• Correlations are reported between -1 and 1; intervals and summary tests use the Fisher (atanh) scale and transform interval limits back.

• Fixed parameters remain in coef() and have zero covariance rows/columns. Their inferential entries, including confidence limits, are NA.

These are Wald intervals. Stored optimizer-coordinate covariance is used when available; supplied matrices and older fits without it use the natural-scale delta method. Numerical convergence alone does not establish interval coverage.

See Also

summary.panelselect(), add_information(), check_stationarity()

Examples

dat = data.frame(person = rep(1:20, each = 3), y = rep(c(0, 1, 3, 1, 2, 4), 10))
fit = panel_select(y ~ 1, NULL, dat, "person", "poisson",
                   fixed = c(sd_re_outcome = 0))
coef(fit)
vcov(fit)
confint(fit, parm = "outcome.(Intercept)")
fit = add_information(fit)
vcov(fit, type = "observed")

Recover fitted formulas, designs and omitted rows

Description

Inspect either equation of a fitted model using its retained data, factor levels and contrasts. The methods describe the fitted design rather than constructing a design for new prediction data.

Usage

## S3 method for class 'panelselect'
formula(x, equation = c("outcome", "selection"), ...)

## S3 method for class 'panelselect'
model.frame(formula, equation = c("outcome", "selection"), ...)

## S3 method for class 'panelselect'
terms(x, equation = c("outcome", "selection"), ...)

## S3 method for class 'panelselect'
na.action(object, ...)

## S3 method for class 'panelselect'
model.matrix(object, equation = c("outcome", "selection"), ...)

Arguments

Value

A formula, terms object, data frame, numeric design matrix, or omit/exclude row-index object, according to the generic.

Returned designs

• formula() returns the fitted equation formula, with .outcome selection shorthand expanded to the regressors used by the fit.

• terms() returns its terms object, including offset information.

• model.frame() rebuilds the retained rows in input order, using the fitted factor levels and allowing structural missing outcome entries.

• model.matrix() constructs the corresponding coefficient design with the fitted contrasts. Offsets belong to the formula/model frame, not additional estimated-coefficient columns.

• na.action() identifies rows omitted or excluded for required missing data; it returns NULL when none were dropped. Rows outside subset are separate.

Outcome values/regressors unused on unselected rows can remain missing in the returned frame. The likelihood's internal zero placeholders are not imputed outcomes. Categories represented only among unselected rows can produce an unidentified outcome coefficient and a rank-deficiency error; levels are not silently dropped as a model-selection step. Prepare an estimable outcome design and check new-data factor support explicitly.

See Also

panel_select(), predict.panelselect(), update.panelselect()

Examples

dat = data.frame(person = rep(1:20, each = 3), y = rep(c(0, 1, 3, 1, 2, 4), 10))
fit = panel_select(y ~ 1, NULL, dat, "person", "poisson",
                   fixed = c(sd_re_outcome = 0))
formula(fit)
head(model.frame(fit))
head(model.matrix(fit))
na.action(fit)

Inspect likelihood, sample size and information criteria

Description

Extract the fitted joint log likelihood or compare models fitted to the same sample and response definition. Within-person dependence is included in each panel's integrated likelihood contribution; observations within a person are not treated as independent likelihood units.

Usage

## S3 method for class 'panelselect'
logLik(object, ...)

## S3 method for class 'panelselect'
nobs(object, ...)

## S3 method for class 'panelselect'
AIC(object, ..., k = 2)

## S3 method for class 'panelselect'
BIC(object, ..., unit = c("panels", "observations"))

Arguments

Value

logLik() returns a logLik scalar with df, nobs and n_rows attributes. nobs() returns the retained row count. AIC/BIC return one number for one model or a data frame with df and the criterion for several.

Calculations and sample counts

If ell is the fitted log likelihood and p is the number of free parameters, AIC = -2 * ell + k * p and BIC = -2 * ell + log(n) * p.

Comparing models responsibly

• Use the same retained sample and response definition. Equal sample counts do not establish sample identity, and criteria from different likelihood definitions need not be comparable.

• The existing guard warns for different row counts, and for default BIC also different independent-unit counts. It does not compare row identities.

• Check integration sensitivity and optimization before interpreting small criterion differences. A higher-order fit need not be a more accurate fit.

• Use likelihood-ratio tests only with the appropriate nesting, identified parameters and boundary conditions; see panel_select().

See Also

panel_select(), update.panelselect(), check_quadrature()

Examples

dat = data.frame(person = rep(1:20, each = 3), y = rep(c(0, 1, 3, 1, 2, 4), 10))
fit = panel_select(y ~ 1, NULL, dat, "person", "poisson",
                   fixed = c(sd_re_outcome = 0))
logLik(fit)
c(rows = nobs(fit), people = attr(logLik(fit), "nobs"))
AIC(fit)
BIC(fit)
BIC(fit, unit = "observations")

Control numerical integration and optimization

Description

Collects numerical settings for a panel selection fit with panel_select(), also used by its cross-sectional selection and standalone models. Person effects use fixed Gauss-Hermite (GH), 10 nodes per effect, by default; the separate PLN observation integral uses fixed GH with 20 nodes. The default optimizer is optim BFGS with 500 iterations. These are starting settings to assess, not guarantees of accuracy. Missing-data handling, initialization and the reported covariance belong to panel_select() (na_action, init, vcov).

Usage

panel_control(
  quadrature = c("gh", "agh", "gl"),
  nodes = NULL,
  quadrature_options = list(),
  observation_quadrature = c("gh", "agh", "gl"),
  observation_nodes = NULL,
  observation_options = list(),
  mode_maxit = 50L,
  mode_tol = 1e-09,
  exact_reductions = TRUE,
  optimizer = c("optim", "maxLik"),
  method = "BFGS",
  maxit = 500L,
  reltol = sqrt(.Machine$double.eps),
  factr = 1e+07,
  parscale = NULL,
  trace = 0L,
  stationarity_warning = 0.1,
  ...
)

## S3 method for class 'panelselect_control'
print(x, ...)

Arguments

Value

A panelselect_control list of validated settings, printed by print().

Integration rules

Analytical scores differentiate the finite rule actually optimized. Correct differentiation and backend convergence do not establish accuracy of the continuous integral. Higher order is not uniformly more accurate; inspect likelihood and score sensitivity with check_quadrature(), especially for large scales, near-boundary correlations and rare response histories. Equal GH/GL orders are neither equal work nor matched accuracy, because GL uses that many nodes in each segment.

PLN has two integration levels. Choosing outer AGH or GL does not change the fixed-GH observation default. Coarse fixed observation nodes can make the person-effect integrand corrugated; observation_quadrature = "agh" differentiates its own moving observation geometry and can help in some cases, with additional computation. Select both rules explicitly for GL at both levels. Refine observation and person-effect orders separately and then jointly before interpreting changed estimates.

Exact reductions skip integrals made unnecessary by fixed zero scales or correlations. Cross-sectional selection without PLN is closed form. A selected rule can therefore do no integration work in a reduced model.

GL fitting and limits

Standalone panel probit and selected binary/Linear panels with cor_error fixed at zero use GL with analytically differentiated moving cuts. Other GL routes keep a positive mesh fixed within each optimization and information stage. A failed staged gate at any trial point ends that stage; the next stage starts from its best valid point on a refined mesh. At most three stages are allowed, including stages needed for requested covariance. Unresolved preparation or exhausted limits fail explicitly and retain diagnostic records. No rule, tolerance or covariance estimator is silently substituted.

Staged GL can be substantially more expensive than GH/AGH because it prepares, checks and refines meshes. Cost depends on the model, active dimensions, restrictions and accuracy target; the moving-cut routes and staged PLN routes do not share a universal timing ratio. Physical moving-cut GL can also be difficult when the person-effect correlation approaches -1 or 1. Refinement agreement remains empirical evidence, not a certified integration-error bound.

GL observed information must pass step/order checks and, when positive definite, inverse/SE stability checks; otherwise it remains unavailable. GL supports population prediction, serialization and explicit update/refit. History-conditioned prediction and random_effects() are unavailable when GL is selected for person effects or the PLN observation rule. There is no prediction-time integration override; use an explicit GH/AGH refit if those targets are required.

GL options

Both option lists accept the following positive finite scalar settings. Orders and capacity controls must be integers. These options are checked against the selected rule; unknown or misspelled names are errors.

Score budgets account for coordinate units: coefficient budgets use each regressor's root mean square. A small canonical/reference discrepancy is an empirical interior-error estimate; analytic omitted-tail envelopes are separate. More interior nodes cannot resolve an omitted-tail bound outside the chosen domain. Increase the named resource or domain limit explicitly when justified, then repeat the comparison without relaxing its tolerances. Observation time limits accumulate work across a person's observed rows, including observation-mesh refinement. They are not a fresh allowance for each row. Time guards run between bounded numerical blocks, so an indivisible block can finish after its deadline before the failure is raised. With outer GH/AGH and observation GL, the overall guard is the larger of 300 seconds and the supplied observation total limit. See the staged-quadrature mathematical chapter for the full budget accounting.

Optimizer methods

Optimizer coordinates are unrestricted coefficients, log SDs and Fisher correlations; supplied starts and reported estimates remain natural-scale. optim minimizes the negative log likelihood. The backend status reports that method's termination condition. SANN has no convergence test, and a zero-iteration run is reported as not run; neither is labelled converged. Newton methods may require curvature during optimization even with OPG covariance. No method guarantees a global maximum.

Without explicit scales, BFGS uses the reciprocal initial panel-score column norms, 1 / sqrt(colSums(score^2)), for regression coefficients with non-negligible initial scores. Other coordinates use 1 / pmax(1, sqrt(colSums(score^2))); all scales are bounded below by sqrt(.Machine$double.eps). This respects coefficient units under rescaled regressors and Linear outcomes. Fixed-GH panel binary selection with free cor_error and otherwise default optim BFGS settings instead uses parscale = 1/sqrt(n_panels) with reltol = 1e-12. Explicit parscale or reltol prevents that automatic pair. The random-effects auxiliaries used by init = "default" have fresh BFGS optimizer settings; the joint fit's optimizer options do not control them. Numerical auxiliaries inherit the equation's person-effect integration settings, and PLN also inherits the observation settings. The internal Linear auxiliary uses an exact Gaussian likelihood without quadrature. See panel_select() for supplied-start precedence and initialization fallback.

The fit's control keeps requested settings and optimizer$control records effective backend options. update() retains the fitted rules and settings unless explicitly replaced; a fresh panel_control(...) is a complete control replacement, not a partial patch to the old control.

See Also

check_quadrature()

Other fitting functions: panel_select()

Examples

panel_control()
panel_control(nodes = c(selection = 8, outcome = 12))
panel_control(quadrature = "agh", nodes = 15)
panel_control(quadrature = "gl", quadrature_options = list(bound = 16))
panel_control(optimizer = "optim", method = "L-BFGS-B", factr = 1e6)
panel_control(stationarity_warning = NULL)

Fit panel and cross-sectional sample-selection models

Description

panel_select() fits panel and cross-sectional sample-selection models with probit selection and linear, binary, or count outcomes. It also supports standalone random-effects panel models for binary and count outcomes.

Panel selection models allow dependence between selection and outcomes through correlated individual-level random effects, correlated observation-level errors, or both. These models are described in Bailey and Peng (2025) and Peng and Van den Bulte (2024).

Estimation uses maximum likelihood with analytical gradients and numerical integration of random effects. The mathematical engine is deeply optimized using C++. The table below summarizes the available models.

Usage

panel_select(
  outcome,
  selection,
  data,
  id,
  family,
  start = NULL,
  fixed = NULL,
  control = panel_control(),
  subset = NULL,
  na_action = c("fail", "omit", "exclude"),
  init = c("default", "re", "pooled", "cross_sectional", "zero", "unif", "norm"),
  vcov = c("opg", "observed", "sandwich")
)

Arguments

Value

An object of class panelselect. Use summary() to view the estimation results, including parameter estimates and standard errors. The object contains these main components:

• coefficients: named regression coefficients, SDs and correlations on natural scales, including fixed values. Extract them with coef().

• covariance: stored parameter covariance matrices. Use vcov() and summary() for uncertainty estimates with the chosen estimator.

• loglik: total log likelihood at the returned estimates; extract with logLik().

• optimizer: convergence and stationarity information; inspect with check_stationarity().

• integration: available numerical integration diagnostics. Assess sensitivity to the integration rule and node counts with check_quadrature().

• control and fixed: fitting settings and parameter restrictions. Modify them through update().

See panelselect-methods for the full set of fitted-model methods.

Models

In the selection and id columns, Yes means that the argument is required, while No means it must be NULL.

Regression coefficients use ⁠selection.<column>⁠ and ⁠outcome.<column>⁠. sd_re_selection and sd_re_outcome are individual-effect SDs, with correlation cor_re. sd_error is the residual SD for Linear outcomes or the observation-level random effect SD for PLN outcomes. cor_error measures selection-outcome error correlation within an observation. summary() describes the structural parameters for the fitted model.

Formulas and missing data

In selection models, indicators and selection regressors are required on every retained row. Outcomes and outcome regressors are required only on selected rows; unused entries may be NA. With selection = NULL, outcomes are required on every retained row, and omitting missing outcomes does not correct for selection.

Use .outcome to copy outcome regressors into the selection formula; observed ~ .outcome + z adds an exclusion restriction that can aid identification. Ordinary . shorthand is unsupported. Formulas support factors, interactions, transformations and offset().

subset is applied before missing-data checks. na_action handles missing required values: "fail" stops, "omit" drops rows, and "exclude" also restores NA positions in fitting-row predictions and residuals. See panelselect-design for design and factor-support details.

Restrictions with fixed

fixed holds named parameters, including regression coefficients, at natural-scale values. Use the parameter names listed under Models.

• fixed = c(cor_re = 0) removes dependence between individual effects.

• fixed = c(cor_error = 0) removes same-observation error dependence and can substantially speed up estimation. Use fixed = c(cor_re = 0, cor_error = 0) to combine both restrictions when both correlations are present.

• Setting an SD to zero removes that random effect and automatically fixes any unidentified associated correlation to zero. In PLN models, fixed = c(sd_error = 0) removes observation-level heterogeneity, reducing the outcome to Poisson.

Correlations must remain strictly between -1 and 1; Linear sd_error must be positive. Fixed parameters remain in the estimates, have no inferential entries, and do not count toward AIC/BIC.

Starting values

init supplies values not specified in start. The default uses separate random-effects fits for applicable panel selection equations; standalone and cross-sectional models use pooled starts. Failed, nonconvergent, or invalid auxiliary fits warn and retain pooled starts for the affected equation. Use init = "pooled" to bypass auxiliary fitting; see init for other choices.

Named partial starts override automatic values; a complete free-parameter vector bypasses initialization. Starts and fixed values use natural scales and must agree where both are supplied.

update() normally uses fitted coefficients as starts. To initialize afresh, use update(fit, start = NULL, init = "pooled"); changing init alone normally preserves fitted starts. See update.panelselect().

Inference

Use summary(), vcov() and confint() for standard errors and intervals. The default vcov = "opg" uses unit scores and avoids a Hessian calculation. "observed" uses observed information; "sandwich" combines it with unit scores. Both require an additional information calculation.

OPG and observed-information inference assume a correctly specified likelihood. Sandwich inference is asymptotic and requires independent units and consistent coefficient estimates; it does not correct coefficient bias or integration error. Independent units are panels, or rows for cross-sectional models.

fit = add_information(fit) stores observed information without refitting or changing the default covariance. Unavailable requested covariance is never replaced automatically. See panelselect-inference for formulas and interval construction, and add_information() for repeated requests.

Integration and numerical checks

panel_control() describes GH, AGH and GL settings. GL orders count nodes per integration segment, so equal GH/AGH and GL orders are not directly comparable. PLN's observation-level rule is set separately by observation_quadrature; changing the person-effect rule leaves it unchanged.

Use check_stationarity() to inspect optimizer convergence and Q_OPG; the package warns when Q_OPG > 0.1 by default. Stationarity does not establish integration accuracy. check_quadrature() compares likelihoods and scores with parameters fixed. Refit explicitly with update() or refit = TRUE to assess changes in estimates. Agreement supports numerical stability but does not guarantee accuracy or a global maximum.

GL supports population predictions, but history-conditioned predictions and random_effects() are unavailable when GL is used at either integration level. See predict.panelselect() for supported predictions.

References

Bailey, M., & Peng, J. (2025). A Random Effects Model of Non-Ignorable Nonresponse in Panel Survey Data. Available at SSRN doi:10.2139/ssrn.5475626

Peng, J., & Van den Bulte, C. (2024). Participation vs. Effectiveness in Sponsored Tweet Campaigns: A Quality-Quantity Conundrum. Management Science, 70(11):7961-7983. Available at doi:10.1287/mnsc.2019.01897

See Also

panelselect-methods, predict.panelselect(), vignette("vignette", package = "PanelSelect")

Other fitting functions: panel_control()

Examples

# 1. Fit a panel selection model
set.seed(1)
dat = data.frame(person = rep(1:60, each = 4), x = rnorm(240), z = rnorm(240))
u = rnorm(60); v = 0.5 * u + rnorm(60, sd = 0.8)  # correlated person effects
dat$observed = as.integer(0.3 + 0.5 * dat$z + 0.7 * u[dat$person] + rnorm(240) > 0)
dat$y = ifelse(dat$observed == 1, 1 + 0.5 * dat$x + v[dat$person] + rnorm(240), NA)

# z enters the selection equation but is excluded from the outcome.
fit = panel_select(y ~ x, observed ~ .outcome + z, dat, id = "person",
                   family = "linear")
summary(fit)

# 2. Check and refine the fit
check_stationarity(fit)[c("optimizer_converged", "Q_OPG")]
# Compare calculations with more nodes, holding the estimates fixed.
check_quadrature(fit, nodes = 20)

# The initial fit uses default GH; refit with 20 nodes.
fit20 = update(fit, control = panel_control(nodes = 20))
cbind(GH10 = coef(fit), GH20 = coef(fit20))


# 3. Predict population outcome means
targets = data.frame(x = c(-1, 0, 1), z = 0)
# Mean potential outcomes before selection, averaging over person effects.
predict(fit, targets, type = "response",
        conditioning = "population", interval = "confidence")


# 4. Set both selection-outcome correlations to zero
independent = update(fit, fixed = c(cor_re = 0, cor_error = 0))
coef(independent)

# 5. Fit a cross-sectional selection model
# Use one observation per independent person.
cross_dat = dat[!duplicated(dat$person), ]
cross = panel_select(y ~ x, observed ~ .outcome + z, data = cross_dat,
                     id = NULL, family = "linear")
coef(cross)

# 6. Fit a standalone Poisson model with exposure
# Counts are observed for everyone, so there is no selection equation.
dat$exposure = runif(240, 0.5, 2)
dat$count = rpois(240, dat$exposure * exp(0.2 + 0.3 * dat$x + 0.4 * u[dat$person]))
counts = panel_select(count ~ x + offset(log(exposure)), selection = NULL,
                      data = dat, id = "person", family = "poisson")
coef(counts)

Work with fitted PanelSelect models

Description

Use the standard R generics below to inspect or reuse a panel_select() fit. Reported coefficients and covariance matrices use natural parameter scales.

Details

See vignette("vignette", package = "PanelSelect") for a complete analysis workflow and vignette("migration", package = "PanelSelect") for earlier APIs.

See Also

panel_select(), panel_control()

Predict outcome means and selection probabilities

Description

Compute a clearly defined mean, probability or linear predictor at the fitted parameters. Population predictions average over person effects. History predictions instead update each person's effect distribution using explicitly supplied responses and selection indicators; they do not refit coefficients.

Usage

## S3 method for class 'panelselect'
predict(
  object,
  newdata = NULL,
  type = c("response", "selection_probability", "joint", "selected_mean", "link"),
  equation = c("outcome", "selection"),
  conditioning = c("population", "history"),
  history = NULL,
  se.fit = FALSE,
  vcov = NULL,
  aggregate = FALSE,
  weights = NULL,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

## S3 method for class 'panelselect'
fitted(object, ...)

Arguments

Value

Normally a named numeric vector; with aggregate = TRUE, a single estimate. interval = "confidence" returns a matrix with fit, lwr and upr. With se.fit = TRUE, the result is a list containing that fit vector/matrix, a numeric se.fit, and a gradient matrix with natural parameters in columns. Fixed-parameter gradient columns do not add variance.

Choose the target

Write D for selection and ⁠Y*⁠ for the potential outcome before selection. Means below condition on the supplied regressors and any explicit history.

Binary-outcome response means are probabilities. Linear-outcome means use outcome units, while Poisson/PLN means are expected counts, including exposure offsets when specified. Each target represents a fresh within-person shock.

Choose the conditioning information

• Population: integrate over the model's person-effect distribution. This is appropriate for a new person or a population-average target. It is available after GH, AGH and GL fitting and uses analytic target formulas.

• History: integrate over each person's posterior given exactly the rows supplied in history. IDs absent from usable history receive population predictions. The current GH/AGH history calculation uses the fit's stored integration settings; there is no prediction-specific method override.

The function does not infer chronology, past-only windows or leave-one-out samples. For forecasting, supply only information available before the target. Including a target's own observed response informs its person effects and produces a fresh-wave replicated prediction, not a held-out forecast or a prediction of that same realized shock. An absent history row is not an observed D = 0 event. Missing required history values follow the fitted na_action: "fail" stops, otherwise incomplete history rows are omitted.

For selected means the target selection event also reweights the history posterior. The result is a ratio of integrated moments, not an unweighted average of node-specific selected means.

Uncertainty and aggregation

• SEs use the delta method with uncertainty in the fitted parameters. For history targets, derivatives include posterior weights and moving nodes.

• Confidence bounds use a probit scale for probabilities, a log scale for count means, and the identity scale for Linear means and linear predictors, then transform back. Aggregation occurs before this interval calculation.

• A supplied matrix must name all natural parameters, be finite and symmetric, and be positive semidefinite (singular matrices are allowed). Symmetry and eigenvalue checks allow relative roundoff of 1e-10.

• Unavailable covariance produces unavailable SEs; a different estimator is never substituted. These intervals do not include realized-outcome variation.

Fitted values versus residuals

fitted() defaults to the potential-outcome mean, whereas residuals() uses the selected mean by default for selection models. Therefore adding their default results does not generally reconstruct observed selected outcomes.

See Also

random_effects(), residuals.panelselect(), panelselect-inference

Other prediction functions: random_effects(), residuals.panelselect()

Examples

set.seed(1)
dat = data.frame(person = rep(1:60, each = 4), x = rnorm(240), z = rnorm(240))
u = rnorm(60); v = 0.5 * u + rnorm(60, sd = 0.8)  # correlated person effects
dat$observed = as.integer(0.3 + 0.5 * dat$z + 0.7 * u[dat$person] + rnorm(240) > 0)
dat$y = ifelse(dat$observed == 1, 1 + 0.5 * dat$x + v[dat$person] + rnorm(240), NA)
fit = panel_select(y ~ x, observed ~ x + z, dat, "person", "linear")

targets = data.frame(person = 1:2, x = 0, z = 0)
predict(fit, targets, type = "selection_probability", interval = "confidence")
# Condition on each person's observed rows.
predict(fit, targets, type = "selected_mean", conditioning = "history",
        history = dat[dat$person %in% 1:2, ])
predict(fit, aggregate = TRUE, se.fit = TRUE)$fit

Summarize person effects conditional on observed history

Description

Calculate posterior means and SDs of person-specific random-intercept contributions at the fitted parameters. These summaries describe conditional heterogeneity in each equation, rather than re-estimating its coefficients.

Usage

random_effects(object, history = NULL)

Arguments

Value

A data frame with one row per history person and the columns above. A structurally fixed-zero effect has zero posterior mean and SD.

Interpretation and available summaries

The intercepts include their fitted person-effect SD multipliers. For count models they contribute to the log rate, not directly to the expected count. The posterior SDs measure remaining uncertainty about a person's effect at fixed fitted parameters; they are not parameter-estimation SEs or confidence intervals for regression coefficients.

Exactly the supplied history is used, without automatic date filtering. Default summaries use the person's own fitting responses and are descriptive. Only IDs represented in usable history are returned; an empty history returns an empty data frame with the appropriate columns. For predictions of people without observed history, use predict() and its population fallback.

See Also

predict.panelselect(), panel_select()

Other prediction functions: predict.panelselect(), residuals.panelselect()

Examples

set.seed(1)
dat = data.frame(person = rep(1:60, each = 4), x = rnorm(240), z = rnorm(240))
u = rnorm(60); v = 0.5 * u + rnorm(60, sd = 0.8)  # correlated person effects
dat$observed = as.integer(0.3 + 0.5 * dat$z + 0.7 * u[dat$person] + rnorm(240) > 0)
dat$y = ifelse(dat$observed == 1, 1 + 0.5 * dat$x + v[dat$person] + rnorm(240), NA)
fit = panel_select(y ~ x, observed ~ x + z, dat, "person", "linear")
head(random_effects(fit))

Calculate response residuals for observed outcomes

Description

Subtract a predicted mean from each observed outcome on the fitting sample. The default target is the selected mean in a selection model and the response mean in a standalone model. Unselected outcomes have no residual.

Usage

## S3 method for class 'panelselect'
residuals(object, ...)

Arguments

Value

A numeric vector in fitting-row order, with NA for unselected rows. Under na_action = "exclude", excluded fitting rows are also restored as NA; "omit" does not restore dropped rows.

Common choices

• residuals(fit) uses the population mean appropriate to the observed outcomes: selected mean for selection models, response mean otherwise.

• residuals(fit, type = "response") uses the potential-outcome mean; selection can make these residuals systematically nonzero among participants.

• residuals(fit, conditioning = "history", history = dat) uses person effects informed by dat. If it includes each row's own outcome, these are descriptive residuals and should not be interpreted as forecast errors.

fitted() has a different default target in selection models. To reconstruct selected outcomes, use predict(fit, type = "selected_mean") + residuals(fit).

See Also

predict.panelselect(), random_effects()

Other prediction functions: predict.panelselect(), random_effects()

Examples

set.seed(1)
dat = data.frame(person = rep(1:60, each = 4), x = rnorm(240), z = rnorm(240))
u = rnorm(60); v = 0.5 * u + rnorm(60, sd = 0.8)  # correlated person effects
dat$observed = as.integer(0.3 + 0.5 * dat$z + 0.7 * u[dat$person] + rnorm(240) > 0)
dat$y = ifelse(dat$observed == 1, 1 + 0.5 * dat$x + v[dat$person] + rnorm(240), NA)
fit = panel_select(y ~ x, observed ~ x + z, dat, "person", "linear")
summary(residuals(fit))
# Residuals from each person's posterior mean given their own rows.
head(residuals(fit, conditioning = "history", history = dat))

Summarize and print a fitted model

Description

Display the model specification, retained sample, estimates, uncertainty and numerical diagnostics. print(fit) gives a compact coefficient view; summary(fit) constructs the full coefficient table and its printed summary.

Usage

## S3 method for class 'panelselect'
summary(object, vcov = object$control$vcov, level = 0.95, ...)

## S3 method for class 'panelselect'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'summary.panelselect'
print(
  x,
  digits = max(3L, getOption("digits") - 3L),
  signif.stars = getOption("show.signif.stars"),
  ...
)

Arguments

Value

summary() returns a summary.panelselect list. coefficients contains estimate, se, z, p, lci, and uci; stationarity contains the OPG status/threshold/flag and the two named Q values. Other components describe the model, likelihood criteria, covariance and numerical settings. exclusion_restriction is a character vector (character(0) when none, NULL for standalone models). Print methods return their input invisibly.

What to inspect

Q_OPG > stationarity_warning is flagged as substantial (0.1 by default). There is no automatic warning cutoff for Q_Hessian. Neither a small Q nor optimizer convergence proves integration accuracy, identification or a global maximum. Saved native fits are displayed under the current Q convention using stored scores/covariances; inspecting them does not refit or overwrite them.

See Also

panelselect-methods, check_stationarity(), panelselect-inference

Examples

dat = data.frame(person = rep(1:20, each = 3), y = rep(c(0, 1, 3, 1, 2, 4), 10))
fit = panel_select(y ~ 1, NULL, dat, "person", "poisson",
                   fixed = c(sd_re_outcome = 0))
print(fit)
answer = summary(fit, vcov = "observed", level = 0.9)
print(answer, digits = 4, signif.stars = FALSE)
answer$stationarity

Refit a model with changed formulas, data or settings

Description

Refit a panel_select() model using retained data and fitted coefficients as warm starts. Settings and restrictions persist unless explicitly overridden. Use start = NULL to run an initializer afresh instead of warm-starting.

Usage

## S3 method for class 'panelselect'
update(object, formula., ..., evaluate = TRUE)

Arguments

Value

A newly fitted panelselect model, or a fitting call when evaluate = FALSE. The supplied fit is unchanged.

Common updates

Starting values and changed models

Explicit starts must be finite named numeric vectors with unique, valid parameter names; matrices, unnamed vectors and duplicate names are rejected before they can overwrite warm starts. Named partial starts remain supported. Parameters absent from a changed design/family leave inherited starts and restrictions. Newly freed SDs previously fixed at zero get automatic starts. An empty start vector adds no overrides; start = NULL requests automatic starts for all free parameters. Explicit start/fixed conflicts and nonpositive free SD starts are rejected when the fitting call is evaluated.

The stored na_action, init, vcov, integration rule and node orders persist unless overridden. Staged meshes are rebuilt for each update. No method is silently changed to improve convergence or integration accuracy.

Formula shorthand and saved fits

In a dotted selection update, .outcome is one parenthesized term, so . ~ . - x cannot remove just one regressor inside that shorthand; write the selection formula in full for that change. Earlier native 2.0 fits are read using current names, but removed integration/optimizer settings stop with guidance for an explicit replacement. See the migration vignette. The former initializer string "ProbitRE" can remain in a displayed saved fit, but it is not a supported update/refit setting. Supply init = "default" (or another current initializer) explicitly, even with evaluate = FALSE. update(fit, init = "default") keeps fitted warm starts; update(fit, init = "default", start = NULL) reruns the automatic initialization described in panel_select(). Earlier person-effect parameter names in a saved fit are read as sd_re_selection, sd_re_outcome and cor_re; explicit new start and fixed inputs must use these current names. The id data-column argument is unchanged. Previously stored Gaussian information is not repaired merely by printing or calling add_information() again; update/refit explicitly for the corrected outcome-scale differentiation steps.

See Also

panel_select(), panel_control(), panelselect-methods

Examples

dat = data.frame(person = rep(1:20, each = 3), y = rep(c(0, 1, 3, 1, 2, 4), 10))
fit = panel_select(y ~ 1, NULL, dat, "person", "poisson",
                   fixed = c(sd_re_outcome = 0))
observed = update(fit, vcov = "observed")
fresh = update(fit, start = NULL, init = "pooled")
call = update(fit, control = modifyList(fit$control, list(maxit = 1000)),
              evaluate = FALSE)
coef(fresh)