Why Bootstrap a Selection Model?
In a selection model, the probability that an outcome is observed may depend on the outcome itself. This makes uncertainty assessment challenging: the likelihood is built from observed outcomes and response indicators, but the outcomes missing from the data still affect the likelihood through integration. In finite samples, likelihoods can also be asymmetric, biased, or weakly curved. A parametric bootstrap offers a direct way to investigate the sampling behavior of an estimator under a fitted model: simulate complete outcomes and response indicators, hide outcomes where the simulated response indicator is zero, and refit the model to each resulting incomplete data set.
The procedure does not create new information about unobserved outcomes. It asks a conditional, model-based question: if the fitted selection model generated data like these, how much would the estimator vary across repeated samples? Its answer depends on the assumed outcome and response distributions, and on whether the fitted likelihood can be reliably refit.
Model and Observed-Data Likelihood
Let \(Y_i\) be a continuous outcome, \(X_i\) observed covariates, and \(R_i\) an indicator equal to 1 when \(Y_i\) is observed and 0 otherwise. A useful general specification has an outcome density \(f(y\mid x;\theta)\) and a response probability \(\pi(y,x;\gamma)=P(R=1\mid Y=y,X=x)\). The parameters \(\theta\) describe the outcome distribution and \(\gamma\) describes selection. For example, a normal outcome model and logistic response model are
where \(\operatorname{logit}(p)=\log\{p/(1-p)\}\). If \(\gamma_1\ne 0\), response probability depends on the outcome even after conditioning on the covariates, which is a nonignorable selection mechanism. The model must specify the joint data-generating assumptions, not just a regression for observed outcomes.
Assume subjects are independent conditional on their observed covariates, the outcome density and response probability are correctly specified, and the fitted parameter is a suitable plug-in approximation to the data-generating parameter. Covariates are commonly held fixed at their observed values. If the covariates themselves are sampled and their sampling variability is part of the target, they must instead be simulated from an appropriate model. Each replicate must be analyzed with the same observed-data likelihood used for the original data.
For one subject, the likelihood contribution conditional on \(x_i\) is
The missing subject contributes a probability integrated over possible outcomes; the refitting step must not treat a simulated missing outcome as observed. The integration and optimization should be accurate enough that numerical error does not dominate the bootstrap variation.
Simulation, Masking, and Refitting
With fitted parameters \((\hat\theta,\hat\gamma)\), one bootstrap replicate is generated as follows. For each subject, draw a latent outcome from \(f(y\mid x_i;\hat\theta)\). Given that outcome and the covariates, draw \(R_i^*\) from a Bernoulli distribution with success probability \(\pi(Y_i^*,x_i;\hat\gamma)\). Retain \(Y_i^*\) only if \(R_i^*=1\). Then maximize the observed-data likelihood for the replicate, obtaining \((\hat\theta^*,\hat\gamma^*)\). Repeat this process independently \(B\) times.
Maximize the selection-model observed-data likelihood and record the estimate, convergence information, and numerical settings.
At each fixed covariate value, simulate an outcome from the fitted outcome model and then a response indicator from the fitted selection model.
Keep the simulated outcome for \(R_i^*=1\); treat it as unobserved for \(R_i^*=0\).
Fit the same likelihood and parameterization, recording convergence, boundary estimates, and integration diagnostics for each replicate.
Use the empirical distribution of the successful, adequately checked refits to describe bias and uncertainty, while investigating any failures.
Worked Example: Generating One Replicate
Suppose a hypothetical fitted model has \(Y\sim N(10,4)\), so the outcome standard deviation is 2, and \(P(R=1\mid Y=y)=\operatorname{expit}(-1+0.2y)\), where \(\operatorname{expit}(a)=1/(1+e^{-a})\). For illustration, take the following simulated standard-normal draws \(z_i\) and uniform draws \(u_i\). Generate \(Y_i^*=10+2z_i\), calculate \(p_i^*=\operatorname{expit}(-1+0.2Y_i^*)\), and set \(R_i^*=1\) exactly when \(u_i<p_i^*\).
| Subject | \(z_i\) | \(u_i\) | \(Y_i^*\) | \(p_i^*\) | \(R_i^*\) |
|---|---|---|---|---|---|
| 1 | -1.0 | 0.30 | 8 | 0.6457 | 1 |
| 2 | 0.5 | 0.75 | 11 | 0.7685 | 1 |
| 3 | 1.0 | 0.90 | 12 | 0.8022 | 0 |
| 4 | -0.5 | 0.20 | 9 | 0.6900 | 1 |
For subject 3, the response linear predictor is \(-1+0.2(12)=1.4\), giving \(p_3^*=1/(1+e^{-1.4})=0.8022\). Since \(0.90\) is greater than \(0.8022\), this subject has \(R_3^*=0\). The replicate therefore contains observed outcomes 8, 11, and 9, plus a missing outcome for subject 3. Although 12 was generated to determine the response indicator, it is hidden from the refit. The missing subject contributes the integral in the observed-data likelihood, not a likelihood term evaluated at 12.
What the Bootstrap Distribution Estimates
Write \(\hat\psi\) for an estimate of a parameter or other prespecified scalar target \(\psi\), and \(\hat\psi_b^*\) for its value in bootstrap replicate \(b\). The empirical distribution of the \(B\) refitted estimates approximates the sampling distribution of the estimator when data are generated under the fitted model. For example, the bootstrap estimate of bias is
A percentile interval uses quantiles of the bootstrap estimates directly. This is a convenient summary, not a guarantee of good coverage: skewness, bias, boundary estimates, and weak identification can all affect its performance. A basic interval reflects bootstrap errors around \(\hat\psi\); Wald intervals use an estimated standard error and a symmetric reference distribution. These alternatives can behave differently near likelihood ridges, where a symmetric approximation may be especially unpersuasive. A profile-likelihood interval is another likelihood-based comparison. No interval method removes dependence on the selection model’s assumptions.
Worked Example: A Percentile Interval and Bootstrap Bias
Suppose a hypothetical analysis estimates the response coefficient as \(\hat\gamma_1=0.18\). In \(B=2{,}000\) adequately checked bootstrap refits, the mean estimate is 0.21, the estimate at sorted rank 50 is \(-0.07\), and the estimate at sorted rank 1,950 is 0.46. Define the empirical quantile as the order statistic at rank \(\lceil Bq\rceil\), for quantile level \(q\).
The estimated bootstrap bias is \(0.21-0.18=0.03\). For a 95% percentile interval, the lower rank is \(\lceil 2{,}000(0.025)\rceil=50\), and the upper rank is \(\lceil 2{,}000(0.975)\rceil=1{,}950\). The resulting interval is \([-0.07,0.46]\) on the response-coefficient scale. It includes zero and is not symmetric around 0.18, which is compatible with an asymmetric bootstrap distribution. This interval is a model-based uncertainty summary; it does not mean there is a 95% probability that the fixed coefficient lies in this particular interval.
Justification and Practical Diagnostics
The justification is a plug-in approximation. If the fitted parameter were the true parameter, repeatedly generating data under it and refitting would reproduce the estimator’s sampling variation under that model. The true parameter is unknown, so the fitted value is substituted. Formally, if \(\mathcal{L}_{\eta}\) denotes the sampling distribution of the estimator when the model parameter is \(\eta=(\theta,\gamma)\), the bootstrap approximates \(\mathcal{L}_{\eta}\) by \(\mathcal{L}_{\hat\eta}\). Thus, the approximation requires both a reasonable fitted model and an estimate \(\hat\eta\) that represents the relevant data-generating region. A bootstrap cannot repair a poorly identified or misspecified likelihood.
The repeated refits are themselves a diagnostic. Examine the fraction that converge, whether estimates accumulate at boundaries, whether different starting values lead to distinct solutions, and whether the bootstrap distribution has long tails or separated modes. For each concerning replicate, check the optimizer’s status, likelihood value, gradient or Hessian information, and integration accuracy. A successful optimizer flag alone does not establish that the global maximum has been found.
Worked Example: Interpreting Failed Bootstrap Refits
In a hypothetical run of \(B=1{,}000\) replicates, 972 refits converge under the initial fitting procedure and 28 do not. The apparent success fraction is \(972/1{,}000=0.972\), or 97.2%. That figure describes the computation, not a statistical confidence level. If the 28 failures occur mainly in replicates with very few observed outcomes, simply discarding them could remove precisely the difficult cases that reveal instability.
Suppose further investigation shows that 20 of the 28 fits converge after trying dispersed starting values and tighter integration tolerances, while 8 still fail. The revised success fraction is \(992/1{,}000=99.2\%\). The remaining 8 should still be investigated; their estimates cannot be silently replaced by convenient values or omitted without comment. If failures persist or successful refits split into distinct likelihood modes, the bootstrap summary is not a dependable description until the numerical and inferential problem is addressed. Report the number and pattern of failures alongside any interval.
A nonparametric bootstrap that resamples observed records is not an automatic substitute: it resamples observed outcomes only among responders and does not generate the missing outcomes needed to reproduce a selection mechanism. A parametric bootstrap does generate them, but only by relying on the specified outcome and response models. Compared with Wald intervals, it can represent skewness and finite-sample behavior more directly; compared with profile likelihood, it requires repeated simulation and refitting. The earlier material on finite-sample bias motivates examining the bootstrap mean, but the bootstrap bias estimate itself can be unstable when the likelihood is weakly identified.
A further limitation is that the bootstrap is centered on one fitted solution. If the likelihood has a ridge or multiple local maxima, a bootstrap around one solution may understate uncertainty across other plausible parameter regions. Likewise, a fitted response model that predicts near-zero response probabilities for some outcomes can generate replicates with little information about those regions. The resulting bootstrap distribution is conditional on the fitted model and fitting procedure; it is not a general correction for nonidentification or sensitivity to assumptions.
Common Mistakes and Interpretation
- Using simulated missing outcomes in the refit. The latent outcome is used to generate its response indicator, then hidden when \(R_i^*=0\). The likelihood contribution for that case must integrate over the unobserved outcome.
- Simulating response independently of the outcome. This discards the selection mechanism whenever \(\pi(y,x;\gamma)\) depends on \(y\). Simulate \(Y_i^*\) first and use it to calculate the response probability.
- Resampling only observed outcomes as if they were a complete sample. This does not reproduce the fitted selection process or the distribution of missing outcomes.
- Treating convergence as proof of validity. Check the likelihood solution, numerical integration, boundary behavior, and sensitivity to starting values.
- Interpreting a percentile interval as assumption-free. Its coverage is only as credible as the model and the bootstrap approximation that generated it.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In a selection model, why must the simulated response indicator be generated after the simulated outcome when response depends on outcome?
- For a subject with \(R_i^*=0\), should the refit use the simulated latent outcome as an observed value? What likelihood contribution is appropriate?
- A bootstrap run has 500 replicates and the mean of the refitted estimates is 1.34, while the original estimate is 1.25. Calculate the estimated bootstrap bias.
- Why might a nonparametric bootstrap of observed records fail to reproduce uncertainty under a selection model?
- What should be investigated if many bootstrap refits fail or converge to distinct solutions?
Solutions. (1) The response probability is conditional on the outcome, so simulate \(Y_i^*\) first and then draw \(R_i^*\) using \(\pi(Y_i^*,x_i;\hat\gamma)\). (2) No. Hide that outcome and use \(\int f(y\mid x_i;\theta)\{1-\pi(y,x_i;\gamma)\}\,dy\) in the observed-data likelihood. (3) The estimated bias is \(1.34-1.25=0.09\), in the units of the estimated quantity. (4) Resampling observed records does not generate outcomes for nonresponders under the fitted outcome and response models, and observed outcomes are selected by response. (5) Check optimizer behavior, starting values, boundary estimates, likelihood modes, gradients or curvature, integration accuracy, and whether failures concentrate in especially uninformative replicates; do not silently discard the failures.