A Bayes-Factor-Guided Approach to Post-Double Selection with Bootstrapped Multiple Imputation

1. Problem Statement

In empirical research involving high-dimensional data, researchers often face two simultaneous challenges:

1. Variable Selection: Identifying a parsimonious set of control variables from a large pool of candidates to estimate the effect of a variable of interest ($D$) on an outcome ($Y$). The Post-Double Selection (PDS) method (Belloni et al., 2014) is a standard approach that selects controls via LASSO in both the outcome equation and an auxiliary equation for $D$.
2. Data Uncertainty: Real-world data often contain missing values (requiring Multiple Imputation, MI) and require sampling uncertainty assessment (requiring Bootstrap resampling).

The Core Issue: When PDS is applied repeatedly across bootstrapped and multiply imputed datasets (BOOT-MI), the set of selected variables varies across iterations.

• Union Rule: Aggregating results by taking the union of all selected variables across iterations leads to overly dense models containing many irrelevant covariates.
• Frequency Thresholds: Aggregating based on fixed frequency thresholds (e.g., selecting variables appearing in >50% of iterations) can be arbitrary, potentially discarding variables with persistent but weaker signals or retaining noise.
• Lack of Stopping Rules: Existing methods typically require fixing the number of bootstrap/imputation iterations ex ante, which is computationally inefficient and lacks a principled stopping criterion.

The authors propose a solution that moves beyond heuristic aggregation rules to a sequential evidence-based framework that determines variable relevance dynamically.

2. Methodology

The proposed method integrates PDS screening within a sequential evidence accumulation framework. It treats the variable selection process across perturbation iterations as a series of Bernoulli trials.

A. The Within-Iteration Procedure (BOOT-MI)

For each iteration $t$:

1. Bootstrap: Draw a sample from the incomplete dataset.
2. Imputation: Perform stochastic imputation to create a completed dataset.
3. PDS Screening:
   • Run LASSO on the outcome equation ($Y \sim D + X$) to get set $S_{Y,t}$.
   • Run LASSO on the auxiliary equation ($D \sim X$) to get set $S_{D,t}$.
   • Form the candidate set $S_t = S_{Y,t} \cup S_{D,t}$.
4. Asymmetric Detection Rule: Define a binary indicator $Z_{jt}$ for variable $j$:
   • If $j \in S_{D,t}$ (selected in auxiliary equation), $Z_{jt} = 1$ (unconditionally retained).
   • If $j \in S_{Y,t}$ only, $Z_{jt} = 1$ only if the coefficient of $X_j$ is statistically significant in the post-selection outcome regression. Otherwise, $Z_{jt} = 0$.
   • This creates a stricter filter for outcome-only variables to reduce overfitting.

B. Sequential Evidence Aggregation

Instead of aggregating $Z_{jt}$ via union or frequency, the authors model the sequence $\{Z_{jt}\}_{t=1}^T$ as a stochastic detection process.

• Hypotheses:
  • $H_0$: Variable $j$ is irrelevant.
  • $H_1$: Variable $j$ is relevant.
• Working Model: Assume $Z_{jt} \sim \text{Bernoulli}(\pi_h)$, where $\pi_1 > \pi_0$.
  • $\pi_1$: Probability of detection given relevance.
  • $\pi_0$: Probability of detection given irrelevance (false positive rate).
• Likelihood Ratio (Bayes Factor): The evidence for $H_1$ vs. $H_0$ is accumulated via the log-likelihood ratio:
  $$\log E_{jt} = \log E_{j,t-1} + \left[ Z_{jt} \log\left(\frac{\pi_1}{\pi_0}\right) + (1-Z_{jt}) \log\left(\frac{1-\pi_1}{1-\pi_0}\right) \right]$$
  • A detection ($Z_{jt}=1$) adds positive evidence.
  • A non-detection ($Z_{jt}=0$) adds negative evidence.

C. Calibration and Stopping

• Pilot Calibration: A short pilot phase (e.g., $T_{pilot}=20$) is used to estimate $\hat{\pi}_0$ and $\hat{\pi}_1$ empirically from the detection frequencies of variables. $\hat{\pi}_0$ is regularized to prevent overly aggressive accumulation.
• Stopping Rule: The procedure runs sequentially until:
  • $\log E_{jt} \geq c$ (Classify as Relevant).
  • $\log E_{jt} \leq -c$ (Classify as Irrelevant).
  • $t = T_{max}$ (Max iterations reached).
• Decision-Theoretic Link: The threshold $c$ corresponds to a Bayes-optimal decision rule under specific loss functions for false inclusion/exclusion.

3. Key Contributions

1. Sequential Decision Framework: Casts variable selection aggregation across BOOT-MI replications as a sequential hypothesis testing problem, eliminating the need to fix the number of iterations ex ante.
2. Bayes-Factor Interpretation: Derives a likelihood-ratio-based evidence measure that admits a Bayes-factor interpretation, providing a principled way to weigh detections vs. non-detections asymmetrically.
3. Asymmetric Detection Rule: Introduces a two-stage detection rule that filters outcome-selected variables more strictly than auxiliary-selected variables, improving stability.
4. Empirical Calibration Strategy: Proposes a pilot-based method to estimate detection probabilities ($\pi_0, \pi_1$) directly from the data, making the method adaptive to specific dataset characteristics.

4. Results

The authors evaluated the method via a Monte Carlo study (126 scenarios) and an empirical application using European Social Survey (ESS) data.

Monte Carlo Simulation

• Variable Selection Accuracy:
  • The Union Rule achieved perfect sensitivity (TPR=1.0) but selected all 50 variables (FPR=1.0), resulting in overly dense models.
  • Frequency Thresholds (50%/75%) were too conservative, missing relevant variables (low TPR) in difficult settings (small samples, MNAR missingness).
  • Proposed Method: Achieved the highest TPR among selective methods (0.774) while maintaining a reasonable model size (21.1 variables vs. 50 for Union). It outperformed frequency thresholds, particularly in small samples and high-missingness scenarios.
• Estimation Performance: The proposed method yielded bias and RMSE comparable to the 50% frequency threshold but with superior coverage (85.1% vs. 84.5%), indicating more reliable inference.
• Computational Efficiency: The sequential stopping rule reduced the average number of iterations from a fixed 200 to 50.7, a 75% reduction in computational cost.

Empirical Illustration (ESS Data)

• Applied to a study of women's employment ($N=8,543$).
• Model Size: The proposed method selected 34 variables, lying between the Union rule (55 variables) and the 75% frequency rule (31 variables).
• Robustness: The estimated coefficient for the variable of interest was statistically insignificant and stable across all aggregation methods, but the proposed method provided a more interpretable, data-driven subset of controls without the arbitrariness of fixed thresholds.

5. Significance and Implications

• Practical Utility: The method solves the "aggregation dilemma" in high-dimensional inference with missing data, offering a middle ground between the inclusiveness of the Union rule and the strictness of frequency thresholds.
• Efficiency: By providing a natural stopping rule, it saves significant computational resources, which is crucial for bootstrapped multiple imputation workflows.
• Theoretical Rigor: It bridges frequentist resampling methods with Bayesian evidence accumulation, offering a decision-theoretic justification for variable inclusion that is robust to the specific calibration of parameters.
• Generalizability: While demonstrated with LASSO and PDS, the framework is applicable to any selection method that produces binary detection indicators, including machine learning algorithms.

Limitations: The method relies on a "working model" assumption that detection indicators are conditionally independent across iterations (which is an approximation given bootstrap dependence). However, simulations suggest this approximation performs well in practice. Future work could explore adaptive calibration during the sequential run and extensions to non-linear models.