Finite-Sample Bias in Selection-Model Estimates

Why assess finite-sample performance?

A selection-model likelihood can contain a long, shallow ridge: quite different combinations of an outcome parameter and a response parameter may fit the observed data nearly equally well. An optimizer may return a numerical estimate, but that does not guarantee that the estimate is centered near the truth in repeated samples or that its confidence interval has its advertised coverage. Simulation makes those finite-sample properties measurable under a specified data-generating model.

The target is not a universal property of an estimator. Bias and coverage depend on the model, the true parameter values, the sample size, the response pattern, and the estimation and interval procedures. A simulation therefore asks a conditional question: if this particular model and design were repeated many times, how would this particular procedure perform?

Model and simulation targets

For a concrete setup, let \(Y_i\) be a continuous outcome and let \(R_i=1\) indicate that \(Y_i\) is observed. Suppose the outcome model and response model are

where \(\operatorname{expit}(u)=1/(1+\exp(-u))\). Here \(\mu\) is the outcome mean, and \(\gamma\) describes how response probability changes with the outcome; \(\alpha\) is the response intercept. The response mechanism depends on \(Y_i\), so this is a missing-not-at-random selection model when \(\gamma\ne 0\). In each simulated data set, \(R_i\) is generated for every subject, but \(Y_i\) is available to the analyst only when \(R_i=1\).

Under this model, an observed subject contributes \(f(y_i;\mu,\sigma^2)\pi(y_i)\) to the likelihood, where \(\pi(y)=\operatorname{expit}(\alpha+\gamma y)\). A subject whose outcome is unobserved contributes \(\int [1-\pi(y)]f(y;\mu,\sigma^2)\,dy\). Thus the simulation must generate and fit the same observed-data likelihood that the analyst would use. The model assumptions include the stated outcome distribution, the stated response model, independent subjects, and correct implementation of the likelihood and interval procedure. Simulated performance does not validate these assumptions for a real data set.

For \(B\) simulation replicates, let \(\hat\theta_b\) be the estimate and \(C_b\) the confidence interval from replicate \(b\). When all \(B\) generated replicates are included and the procedure defines an estimate and interval for every replicate, the Monte Carlo estimates are

The first quantity has the units of \(\theta\); the second is a proportion. The empirical standard deviation of the estimates describes their simulated sampling variability. It is not the Monte Carlo standard error of the estimated bias. If the replicate estimates have standard deviation \(s_{\hat\theta}\), the Monte Carlo standard error of their mean is approximately \(s_{\hat\theta}/\sqrt{B}\). For estimated coverage, the Monte Carlo standard error is approximately \(\sqrt{\hat C(1-\hat C)/B}\).

How likelihood shape can produce bias

A short Taylor expansion shows why a flat likelihood can make finite-sample bias more consequential. Consider a scalar log likelihood \(\ell(\theta)\), with score \(U(\theta)=\ell'(\theta)\), and suppose the maximum-likelihood estimate solves \(U(\hat\theta)=0\). Write \(J=-U'(\theta)\) and \(\delta=\hat\theta-\theta\). A second-order expansion around the true \(\theta\) gives

If \(J>0\) and the quadratic term is a small correction, first substitute \(\delta\approx U(\theta)/J\) into that term. Solving to this order yields

The score has mean zero at the true parameter under the usual regularity conditions, so the leading linear term is centered when \(J\) is treated as fixed. But the curvature term need not have mean zero. In an actual likelihood, \(J\) is random and correlated with the score and higher derivatives, so this expression is a local explanation, not a complete general bias formula. It does show the key vulnerability: when curvature is weak, the inverse-curvature factors can amplify random variation and nonlinear effects.

In a selection model, weak identification often involves a ridge coupling \(\mu\) and \(\gamma\). Movement along the ridge changes both parameters while changing the likelihood little. The estimate can then move substantially in response to small data fluctuations, and the likelihood may be asymmetric or have a boundary maximum. A simulation can reveal the resulting bias and undercoverage, but it should also retain likelihood diagnostics: a poor optimizer, inaccurate integration, or an unintended boundary constraint can create apparent “statistical” bias that is actually a computational defect.

Comparing designs and procedures

A useful simulation varies a feature expected to affect identification, such as sample size, while holding the true parameters and the estimation procedure fixed. The following hypothetical results use 2,000 replicates per setting and 95% intervals. All fits and intervals are available for all replicates. The sample size is \(n=200\) in the weak-information setting and \(n=1{,}000\) in the more informative setting.

For example, the weak-information coverage for \(\gamma\) is \(1{,}750/2{,}000=0.875\); the corresponding Monte Carlo standard error is \(\sqrt{0.875(0.125)/2{,}000}=0.00740\). The more informative estimate, \(0.951\), is close to the nominal 0.95. These results are consistent with improved finite-sample behavior as information increases, but they do not prove that sample size alone caused the difference unless the simulation changed only sample size. The likelihood profiles and convergence records should be reviewed alongside the summary table.

Simulation is not the only way to assess an interval procedure. Wald intervals use an estimate and an estimated standard error with an approximately normal reference distribution; near a ridge, that approximation can be poor. Profile-likelihood intervals use likelihood changes and can better reflect asymmetry, although they too can have inaccurate coverage in finite samples. Simulation evaluates the actual repeated-sampling behavior of either procedure under the chosen model and design. It does not make a flawed likelihood or interval method valid by itself.

Failures, diagnostics, and limits

A failed fit is part of the performance of a computational procedure, not a replicate that may automatically be discarded. If only successful fits are averaged, the resulting bias and coverage estimates are conditional on successful fitting. They do not estimate unconditional performance over all generated data sets. Report the failure rate and label successful-only summaries explicitly. Alternatively, prespecify an operational rule that supplies an estimate and interval for every replicate, and evaluate that complete rule; the rule itself then defines the procedure being assessed.

For example, check whether the estimates cluster at parameter boundaries, whether profile likelihoods are unusually flat or asymmetric, whether different starting values lead to materially different solutions, and whether numerical integration is accurate enough to stabilize the likelihood. Also inspect the distribution of estimates, not only its mean: a long tail or a small group of extreme estimates can produce important finite-sample effects. If failures occur disproportionately in data sets with certain response patterns, a successful-only average may describe a selected subset of the simulation rather than typical repeated use.

The main limitation is dependence on the simulation assumptions. A model can show low bias and nominal coverage when it generates the data and still be wrong for the application. Simulation also cannot establish identification from the observed data alone: it describes performance over the parameter and design settings chosen by the analyst. Results should therefore include the true parameter values, sample size, response mechanism, number of replicates, estimation and interval methods, failure handling, and Monte Carlo uncertainty. The earlier material on likelihood ridges and profile curvature helps interpret why performance changes across settings.