Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
72,184 characters · 25 sections · 22 citation commands
A Bayes-Factor-Guided Approach to Post-Double Selection with Bootstrapped Multiple Imputation
\doublespacing
Estimating the effect of a variable of interest on an outcome in the presence of many potential control variables is a common problem in empirical research. When the number of candidate covariates is large, variable selection methods such as LASSO are routinely used to construct parsimonious models. However, data-driven selection can distort estimation if relevant variables are omitted, particularly when they are associated with both the variable of interest and the outcome Hayashi2000.
The post-double selection (PDS) approach of BelloniCH2014b,BelloniCH2014a addresses this issue by performing variable selection in both the outcome and the variable-of-interest equations and combining the resulting sets of selected variables. This construction increases the likelihood that relevant controls are retained and has become a standard tool for estimation after selection in high-dimensional settings.
In empirical applications, variable selection is often combined with additional procedures that introduce further stochastic variability. Missing covariate data are commonly handled via multiple imputation (MI) Rubin1987,LittleR2019, while sampling uncertainty is frequently addressed through bootstrap resampling Efron1979,EfronT1994. When selection methods are applied repeatedly across such perturbed datasets, the set of selected variables typically varies across iterations. Standard aggregation rules, such as the union of all selected variables, may therefore lead to overly large models, while frequency-based rules may discard variables with weaker but persistent signals.
Although perfect variable selection is not required for reliable estimation after selection BelloniCH2014b, models containing many irrelevant covariates are often undesirable in practice. Moreover, researchers are often interested not only in obtaining stable estimates but also in identifying a meaningful subset of relevant control variables. In the presence of missing data and resampling procedures, assessing variable relevance therefore becomes particularly challenging.
This paper proposes a sequential evidence-based aggregation procedure for variable-selection results generated by repeated bootstrap and multiple-imputation (BOOT-MI) replications that use PDS screening within each iteration. The procedure uses PDS to construct within-iteration candidate sets, but applies an additional detection and aggregation layer across BOOT-MI replications. We implement a sequential BOOT-MI procedure in which each iteration consists of bootstrapping the incomplete dataset, performing stochastic imputation, constructing a PDS candidate set via LASSO in the outcome and variable-of-interest equations, and selectively confirming outcome-selected variables in a final regression step.
Rather than aggregating selected variables using heuristic rules, we construct a likelihood-ratio–based evidence measure with a Bayes-factor interpretation from a working detection model with empirically calibrated null and alternative detection probabilities. We then accumulate evidence for variable relevance sequentially across iterations. The resulting evidence measure provides a decision criterion for variable inclusion and induces a natural stopping rule for the sequential BOOT-MI procedure.
Our contribution is threefold. First, we cast the aggregation of variable-selection outcomes across BOOT-MI datasets as a sequential decision problem. Specifically, we model binary detection indicators as outcomes of a stochastic detection process with hypothesis-dependent probabilities, leading to a likelihood-ratio-type evidence accumulation process. Second, we derive a stopping rule and a variable inclusion criterion based on this evidence process, eliminating the need to fix the number of perturbation iterations ex ante. Third, we propose a pilot-based calibration strategy for the detection probabilities and examine its performance across a range of simulation scenarios.
The proposed procedure is not intended as a new version of PDS itself, but as a calibrated sequential evidence-aggregation framework built around repeated BOOT-MI implementations that use PDS candidate screening. It combines an asymmetric within-iteration detection rule, pilot-based calibration of detection probabilities, and sequential log-evidence accumulation across perturbations. The resulting method is designed to yield interpretable and stable inclusion decisions in settings where repeated perturbations induce substantial variability in selected models.
The remainder of the paper is structured as follows. Section (ref) reviews the related literature. Section (ref) develops the proposed sequential evidence-aggregation framework. Section (ref) reports the Monte Carlo study, and Section (ref) provides an empirical illustration. Section (ref) provides concluding remarks and discusses directions for future research.
This paper relates to four strands of the literature: variable selection with multiply imputed data, bootstrap and resampling approaches to model selection, Bayesian evidence and Bayes factors, and sequential testing procedures. Each strand addresses a component of the problem considered here, but the combination of missing data, resampling-based variability, and principled aggregation of variable-selection outcomes remains only partially explored.
Variable selection in the presence of missing covariate data has received increasing attention in the statistical literature. MI provides a principled framework for handling incomplete datasets by replacing missing values with draws from their predictive distribution Rubin1987,WoodWR2008. By generating multiple completed datasets, MI allows standard complete-data methods to be applied while accounting for imputation uncertainty.
Combining variable selection with MI introduces additional complexity, since selection results may differ across imputed datasets. Several approaches have been proposed to address this issue. ChenW2013 introduce MI-LASSO, which applies a group LASSO penalty across imputations to encourage consistent selection. Alternative approaches rely on stacked or grouped penalized regression formulations DuB2020, while comparative studies evaluate the performance of LASSO-based and Bayesian procedures under missing data BainterM2023. Bayesian methods, such as stochastic search variable selection and reversible jump Markov chain Monte Carlo, have also been applied in this setting GeorgeM1997,OHaraS2009.
A related framework is PDS BelloniCH2014b,BelloniCH2014a, originally developed in a causal inference context. The method performs variable selection in both the outcome equation and an auxiliary equation for the variable of interest and combines the resulting sets of selected variables. While this construction increases the likelihood that relevant controls are retained, it is typically combined with simple aggregation rules, such as the union of selected variables, which may lead to overly large models when applied across multiple imputations or resampled datasets.
Overall, this strand of literature highlights the difficulty of achieving stable and interpretable variable selection in the presence of missing data and provides limited guidance on how to aggregate selection outcomes across multiple stochastic completions of the data.
Resampling methods are widely used to assess sampling variability in variable selection procedures. Bootstrap-based approaches generate perturbed versions of the data and apply selection methods repeatedly, thereby providing information on the stability of selected variables. Several studies combine bootstrap resampling with MI to jointly account for sampling and imputation uncertainty LongJ2015,MusoroTORB2014.
A closely related approach is stability selection MeinshausenB2010, which evaluates variable importance based on selection frequencies across subsamples. Stability selection provides theoretical guarantees on the expected number of false selections under suitable assumptions and can be applied as a wrapper around a wide range of selection methods. However, its aggregation mechanism is based on frequency thresholds, which treat all detections symmetrically and do not account for the sequential structure of the selection process.
Recent work has also explored machine learning approaches to variable selection under missing data and resampling. For example, GunnR2022 study procedures that combine cross-validation, sample splitting, and MI in predictive modeling settings. These approaches emphasize predictive performance and model stability, but they typically rely on heuristic aggregation rules.
Taken together, resampling-based approaches highlight the intrinsic variability of variable selection under data perturbations. However, they generally rely on aggregation rules such as frequency thresholds or union rules, which may either include too many irrelevant variables or exclude variables with weaker but persistent signals. In contrast, the approach proposed in this paper aggregates selection outcomes through a sequential evidence-accumulation mechanism that accounts for the asymmetric informational content of detections and non-detections.
Complementary to resampling-based approaches, the literature has developed evidence-based frameworks for model comparison and variable selection. Bayes factors provide a formal mechanism for quantifying evidence in favor of competing hypotheses or models KassR1995. Unlike classical hypothesis testing based on $p$-values, they allow direct comparison between alternative explanations and quantify evidence for both the null and the alternative.
The approach proposed in this paper differs from a fully Bayesian variable selection method. Rather than specifying a complete probabilistic model over parameters and model space, we construct a likelihood-ratio--based evidence measure with a Bayes-factor interpretation based on a working detection model with empirically calibrated parameters. This strategy is related to a broader literature on approximate Bayes factors, including the Schwarz criterion Schwarz1978 and calibrated Bayes factors that reinterpret frequentist test statistics in evidential terms Sellke2001.
Unlike these approaches, which operate on model likelihoods or test statistics, our method is based on binary detection outcomes generated by repeated perturbation and selection. This perspective yields a tractable and interpretable evidence measure tailored to the aggregation of variable-selection results across BOOT-MI datasets.
Sequential hypothesis testing, introduced by Wald1945, provides a framework for accumulating evidence through repeated observations and terminating data collection once sufficient evidence has been reached. The Sequential Probability Ratio Test (SPRT) is based on cumulative likelihood ratios and admits both frequentist and Bayesian interpretations. In particular, SPRT decision boundaries can be related to Bayes-factor thresholds under suitable prior specifications EdwardsLJ1963.
The approach proposed in this paper draws on this connection by interpreting the cumulative likelihood ratio of the detection model as a sequential evidence measure with a Bayes-factor interpretation. However, unlike classical sequential testing settings, the observations in our framework are not independent draws from a sampling distribution but rather detection outcomes generated by stochastic perturbations of a fixed dataset. The resulting evidence accumulation mechanism is therefore motivated by, but not formally equivalent to, the SPRT framework.
Despite these connections, previous works provide limited guidance on how to apply sequential evidence accumulation to the aggregation of variable-selection outcomes across BOOT-MI datasets. The present paper addresses this gap by developing a sequential evidence-based aggregation procedure tailored to this setting.
The methodology uses perturbation-based PDS screening as a within-iteration device and adds a detection model together with a sequential evidence-aggregation rule across perturbations.
We consider estimation in a sparse linear model with a large set of potential control variables, some of which may contain missing values. Let $Y$ denote the outcome variable, $D$ the variable of interest, and $X = (X_1,\dots,X_p)$ a vector of candidate covariates. Our parameter of interest is the coefficient $\alpha$ in the partially linear regression model
Following PDS BelloniCH2014b,BelloniCH2014a, variable selection is based on two sparse regression problems corresponding to the outcome equation and an auxiliary equation for $D$. In particular, we consider LASSO regressions of $Y$ on $(D, X)$ and of $D$ on $X$, that is,
The outcome equation identifies variables predictive of $Y$, while the equation for $D$ identifies variables predictive of the variable of interest. In the standard PDS approach, the selected control set within a given dataset is formed by the union rule
where $S_{Y,t}$ and $S_{D,t}$ denote the sets of variables selected in the outcome equation and the auxiliary equation for $D$, respectively, in perturbation iteration $t$.
In each perturbation iteration, we generate a single completed dataset via stochastic imputation. Within the perturbation stage, each bootstrap draw is completed by one stochastic imputation and then passed to the selection/detection step. Thus, the repeated perturbation loop does not pool over multiple imputations within a given iteration. Instead, imputation uncertainty enters the selection stage through repeated stochastic completions across perturbations. Classical multiple-imputation pooling is used only after the final variable set has been determined, when the parameter of interest is re-estimated across $M$ completed datasets and combined using Rubin's rules Rubin1987. Specifically, if $Q^{(m)}$ denotes the estimate of the parameter of interest and $U^{(m)}$ its complete-data variance estimate in imputed dataset $m=1,\dots,M$, the pooled point estimate is \[ \bar{Q} = \frac{1}{M}\sum_{m=1}^M Q^{(m)}. \] The average within-imputation variance is \[ \bar{U} = \frac{1}{M}\sum_{m=1}^M U^{(m)}, \] and the between-imputation variance is \[ B = \frac{1}{M-1}\sum_{m=1}^M \left(Q^{(m)}-\bar{Q}\right)^2. \] The total variance is then given by \[ T = \bar{U} + \left(1+\frac{1}{M}\right)B. \] Hence, Rubin's rules combine estimation uncertainty within each imputed dataset and additional uncertainty arising from variation across imputations.
Missing covariate data are handled using MI, while sampling uncertainty is addressed through bootstrap resampling. In this paper we focus on the BOOT-MI approach, in which bootstrap resampling is performed before imputation. Let $t = 1,\dots,T$ index perturbation iterations. In each iteration, we proceed as follows:
Thus, each perturbation iteration produces a selected set $S_t$ of controls. Figure (ref) illustrates the workflow of the perturbation-based variable selection procedure.
Because both bootstrap sampling and stochastic imputation introduce randomness, the selected sets $S_t$ vary across perturbation iterations. While each individual set is typically sparse, the union across iterations, \[ S_{\mathrm{union}} = \bigcup_{t=1}^T S_t, \] can become large when selection is unstable across perturbation iterations, thereby producing overly inclusive aggregated control sets.
To summarize variable-specific outcomes, we define for each variable $j$ and iteration $t$ a binary detection indicator $Z_{jt}$ (formalized in Section (ref)). The sequence $\{Z_{jt}\}_{t=1}^{T}$ records whether variable $j$ is detected in each perturbed dataset and forms the basis for aggregation.
Rather than applying fixed aggregation rules such as unions or frequency thresholds, we interpret $\{Z_{jt}\}$ as realizations of a stochastic detection mechanism and aggregate them through a sequential evidence process.
The perturbation-based procedure produces, for each variable $j$, a sequence of binary indicators $\{Z_{jt}\}_{t=1}^{T}$ summarizing variable-specific detection outcomes across perturbation iterations.
In contrast to pure selection indicators, we define a two-stage detection event that combines PDS screening with outcome-model confirmation. Specifically, let $S_{Y,t}$ denote the set of variables selected by LASSO in the outcome equation and let $S_{D,t}$ denote the set of variables selected in the auxiliary equation for $D$ at perturbation iteration $t$. We form the PDS candidate set \[ S_t = S_{Y,t} \cup S_{D,t}, \] which corresponds to the standard PDS union rule applied within each perturbation iteration.
To aggregate variable-selection outcomes across perturbation iterations, we require a variable-specific representation of the selection history. We therefore map the within-iteration PDS candidate set into a binary detection indicator for each variable and iteration. The next paragraph defines this detection rule, clarifies how it relates to the standard PDS inclusion rule, and then introduces the working probabilistic model that underlies sequential evidence accumulation.
\paragraph{Asymmetric detection rule.} For each variable $j$ and iteration $t$, we define the detection indicator
Thus, variables selected in the auxiliary equation ($j \in S_{D,t}$) are retained unconditionally, while variables selected only in the outcome equation ($j \in S_{Y,t}$) must additionally pass a significance test in the post-selection outcome regression.
This asymmetric rule reflects the structure of the PDS framework: variables predictive of $D$ are retained automatically, whereas variables selected only in the outcome equation are filtered more strictly to reduce overfitting and instability across perturbations. The resulting procedure should be interpreted as a detection rule built on the PDS candidate set rather than as the standard PDS inclusion rule. Accordingly, the formal inferential guarantees of standard PDS do not automatically carry over to the full procedure considered here. Because this detection rule is not identical to the standard PDS union rule, it is useful to state its relation to PDS explicitly.
\paragraph{Relation to PDS.} The proposed detection rule departs from the standard PDS inclusion rule by applying an additional confirmation step to variables selected only in the outcome equation. In particular, defining \[ S_t^{\mathrm{ASYM}} = \{ j : Z_{jt} = 1 \}, \] we have \[ S_t^{\mathrm{ASYM}} \subseteq S_t = S_{Y,t} \cup S_{D,t}. \] That is, the procedure can be interpreted as a data-dependent pruning of the PDS candidate set.
While the PDS framework provides formal guarantees under suitable conditions, these guarantees do not carry over to the modified rule. In particular, the procedure may omit weak but relevant controls, which can affect estimation of the coefficient on the variable of interest. The magnitude of this risk is therefore assessed empirically in Section (ref).
Having clarified how the detection indicator is constructed and how it departs from standard PDS, we now use the resulting sequence $\{Z_{jt}\}$ as the basic input for the sequential evidence framework.
\paragraph{Detection model.} For the purpose of sequential aggregation, we reduce the variable-specific decision problem to a binary distinction between relevance and irrelevance. We consider the hypotheses \[ H_0: \text{Variable } j \text{ is irrelevant}, \qquad H_1: \text{Variable } j \text{ is relevant}. \]
Under this formulation, the sequence $\{Z_{jt}\}$ represents repeated outcomes of a stochastic detection mechanism induced by the combined selection-and-testing procedure across perturbation iterations.
We model this detection mechanism using a Bernoulli working model with effective detection probabilities \[ Z_{jt} \mid H_0 \sim \text{Bernoulli}(\pi_0), \qquad Z_{jt} \mid H_1 \sim \text{Bernoulli}(\pi_1), \] where \[ \pi_0 = \operatorname*{Prob}(Z_{jt}=1 \mid H_0), \qquad \pi_1 = \operatorname*{Prob}(Z_{jt}=1 \mid H_1), \] and $\pi_1 > \pi_0$.
This model is a working approximation: the detection indicators are neither independent across iterations nor generated by a true Bernoulli process. Rather, the model provides a parsimonious representation of the marginal detection behavior under relevance and irrelevance, which enables a tractable likelihood-ratio-type evidence construction in the next subsection. The conditional independence assumption is adopted as a simplifying device to obtain a tractable likelihood-ratio representation; the resulting evidence measure is therefore an approximation whose practical validity is assessed empirically.
The parameters $\pi_0$ and $\pi_1$ summarize the effective detection behavior of the full two-stage procedure. In particular, $\pi_0$ reflects both the probability that an irrelevant variable enters the candidate set and the probability that it is subsequently confirmed. These quantities are therefore treated as working parameters characterizing the empirical detection behavior of the procedure.
For each variable $j$, the sequential aggregation problem is formulated as a comparison between the hypotheses \[ H_{0j}: \text{variable } j \text{ is irrelevant}, \qquad H_{1j}: \text{variable } j \text{ is relevant}. \] Under the working detection model from Section (ref), the sequence of binary detection indicators $\{Z_{js}\}_{s=1}^{t}$ has hypothesis-dependent working distributions. For a single perturbation iteration $s$, define the likelihood ratio \[ L_{js} = \frac{\operatorname*{Prob}(Z_{js}\mid H_{1j})}{\operatorname*{Prob}(Z_{js}\mid H_{0j})} = \frac{\pi_1^{Z_{js}}(1-\pi_1)^{1-Z_{js}}} {\pi_0^{Z_{js}}(1-\pi_0)^{1-Z_{js}}}. \]
Under the Bernoulli working model, $L_{js}$ compares the relevance hypothesis $H_{1j}$ to the irrelevance hypothesis $H_{0j}$ for the single observation $Z_{js}$. If the detection indicators were conditionally independent across perturbation iterations, the joint likelihood ratio for $\{Z_{js}\}_{s=1}^t$ would factorize as the product of the per-iteration likelihood ratios. This motivates defining the cumulative quantity \[ E_{jt} = \prod_{s=1}^{t} L_{js}, \] which aggregates evidence sequentially across perturbations. Under equal prior weights, each $L_{js}$ corresponds to a Bayes factor under the working model, so $E_{jt}$ admits a Bayes-factor–type interpretation.
For numerical stability, we work with the logarithmic form
initialized at $\log E_{j0}=0$. Each perturbation iteration contributes additively: a positive detection adds $\log(\pi_1/\pi_0)>0$, while a non-detection adds $\log((1-\pi_1)/(1-\pi_0))<0$. Hence the evidence process can be updated recursively as \[ \log E_{jt} = \log E_{j,t-1} + \log L_{jt}. \]
This representation shows that the evidence process is a weighted cumulative sum of detection outcomes, where detections and non-detections contribute asymmetric increments determined by $(\pi_0,\pi_1)$.
A useful summary of the working evidence model is the break-even detection frequency, that is, the detection probability at which the expected increment of the log-evidence process is exactly zero. If a variable has effective marginal detection probability $q$, then the expected log-evidence increment under the working model is \[ m(q;\pi_0,\pi_1) = q\log\!\left(\frac{\pi_1}{\pi_0}\right) + (1-q)\log\!\left(\frac{1-\pi_1}{1-\pi_0}\right). \] This quantity is positive if detections occur sufficiently often and negative otherwise. The critical value at which the drift changes sign is \[ q^\ast(\pi_0,\pi_1) = \frac{ \log\!\left(\frac{1-\pi_0}{1-\pi_1}\right) }{ \log\!\left(\frac{\pi_1(1-\pi_0)}{\pi_0(1-\pi_1)}\right) }. \] Thus, $q^\ast(\pi_0,\pi_1)$ is the detection frequency for which the expected log-evidence increment is zero. Variables with effective detection rates above $q^\ast$ have positive expected log-evidence drift, whereas variables with detection rates below $q^\ast$ have negative expected drift under the working model. In the implementation, this quantity is used as a diagnostic for whether a given calibration of $(\pi_0,\pi_1)$ yields an overly permissive evidence process.
Because $(\pi_0,\pi_1)$ are calibrated empirically and perturbation iterations are not independent, the resulting evidence measure should be interpreted as a working-model approximation rather than as an exact Bayesian posterior-odds update. The construction nevertheless provides a tractable and interpretable score that mimics likelihood-ratio accumulation.
The speed of evidence accumulation is governed by the separation between $\pi_1$ and $\pi_0$. Larger separation implies larger expected increments in absolute value and therefore faster threshold crossing. Accordingly, the proposed procedure is best viewed as a likelihood-ratio–motivated scoring and stopping rule for repeated detection outcomes, rather than as a formally calibrated Bayesian or frequentist testing procedure. Relative to simple frequency thresholding, the likelihood-ratio formulation provides asymmetric weighting of detections and non-detections and yields a natural stopping rule based on cumulative evidence rather than a fixed number of iterations.
The evidence process also admits a decision-theoretic interpretation under the working model. In particular, the next proposition shows that thresholding the cumulative evidence is equivalent to the Bayes-optimal classification rule when false inclusion and false exclusion are assigned explicit losses and prior relevance probabilities are specified.
A proof is provided in Appendix (ref). As with the other theoretical results in this section, the proposition is derived under the idealized working model and should therefore be interpreted as an approximation to the behavior of the empirical procedure rather than as a finite-sample guarantee.
Proposition (ref) shows that, under the working detection model, the cumulative evidence measure $E_{jT}$ fully determines the Bayes-optimal classification of variable relevance. Thus, the proposed threshold rule is not merely heuristic: it coincides with the optimal decision rule under asymmetric misclassification costs and prior relevance probabilities.
In particular, symmetric loss and prior specifications imply a threshold at zero on the log-evidence scale, while more conservative inclusion rules correspond to higher relative loss for false inclusion or higher prior odds in favor of irrelevance. The result therefore provides a direct mapping between threshold choice and the implied decision-theoretic trade-off between false inclusion and false exclusion.
To provide intuition for the proposed evidence process, consider an idealized setting in which the detection indicators $\{Z_{jt}\}_{t\geq 1}$ are independent and identically distributed with \[ Z_{jt} \sim \mathrm{Bernoulli}(q_j), \] where $q_j = \operatorname*{Prob}(Z_{jt}=1)$ denotes the marginal detection probability.
Under the working model, the expected increment of the log-evidence process is \[ m(q_j;\pi_0,\pi_1) = q_j\log\!\left(\frac{\pi_1}{\pi_0}\right) + (1-q_j)\log\!\left(\frac{1-\pi_1}{1-\pi_0}\right). \] This quantity determines the direction of evidence accumulation. In particular, the drift is positive if $q_j > q^\ast(\pi_0,\pi_1)$ and negative otherwise, where $q^\ast(\pi_0,\pi_1)$ is the break-even detection probability defined in Section (ref).
The following proposition makes this intuition precise by showing that the sign of the drift determines the asymptotic classification outcome under the idealized working model.
A proof is provided in Appendix (ref).
The result provides intuition for why the proposed score can separate persistent detections from persistent non-detections under the working model. Variables with stronger separation from the break-even detection rate are classified more rapidly, while variables near the boundary require more iterations.
Under the same idealized working model, the expected stopping time is approximately inversely proportional to the absolute drift of the log-evidence process; a formal statement is provided in Appendix (ref).
The formal proofs and additional results on boundary-crossing probabilities are provided in Appendix (ref). The practical behavior of the evidence process therefore depends on how the working detection probabilities are calibrated, which is discussed next.
The evidence process requires specification of the working detection probabilities $(\pi_0,\pi_1)$. We calibrate these parameters in a short pilot phase and hold them fixed during subsequent evidence accumulation.
Let $T_{\mathrm{pilot}}$ denote the number of pilot iterations and define the pilot detection frequencies \[ \hat{\pi}_{j,\mathrm{pilot}} = \frac{1}{T_{\mathrm{pilot}}} \sum_{t=1}^{T_{\mathrm{pilot}}} Z_{jt}. \]
We estimate $\pi_0$ from variables with low pilot detection frequencies and $\pi_1$ from variables with high pilot detection frequencies, using lower- and upper-quantile subsets of $\{\hat{\pi}_{j,\mathrm{pilot}}\}_{j=1}^p$.
To ensure stable behavior of the evidence process, we regularize $\hat{\pi}_0$ by shrinking it toward the nominal level and imposing a lower bound. This prevents excessively aggressive evidence accumulation when the raw pilot estimate is very small.
The resulting calibrated values $(\hat{\pi}_0,\hat{\pi}_1)$ define the likelihood-ratio increments used in the sequential evidence process. Additional calibration strategies and sensitivity analyses are discussed in Appendix (ref).
The sequential evidence process provides a variable-specific rule for classifying controls as relevant or irrelevant as perturbation iterations proceed. Let $E_{jt}$ denote the cumulative evidence in favor of $H_1$ for variable $j$ after $t$ evidence iterations. Because $E_{jt}$ is constructed as a likelihood-ratio--based evidence measure under the working detection model, variable classification can be formulated through evidence thresholds.
Specifically, for thresholds $0<\tau_0<1<\tau_1$, a variable is classified as relevant once \[ E_{jt}\geq \tau_1, \] and as irrelevant once \[ E_{jt}\leq \tau_0. \] Thus, $\tau_1$ is the upper evidence threshold required for classifying a variable as relevant, whereas $\tau_0$ is the lower evidence threshold for classifying a variable as irrelevant.
Equivalently, in logarithmic form, \[ \log E_{jt}\geq c_1 \qquad\text{or}\qquad \log E_{jt}\leq c_0, \] where $c_1=\log\tau_1$ and $c_0=\log\tau_0$.
In the baseline implementation, we use symmetric log-thresholds $c_1=c$ and $c_0=-c$ for some $c>0$, so that \[ \tau_1=e^c \qquad\text{and}\qquad \tau_0=e^{-c}. \] Thus, $c$ is the required magnitude of cumulative log-evidence for classification: a variable is classified as relevant once its log-evidence reaches $+c$ and as irrelevant once its log-evidence reaches $-c$. Equivalently, on the evidence scale, $e^c$ is the upper evidence threshold in favor of relevance and $e^{-c}$ is the lower threshold in favor of irrelevance.
Hence, larger values of $c$ require stronger accumulated evidence before a classification is made, leading to more conservative decisions and typically longer perturbation runs, whereas smaller values of $c$ allow earlier classifications but may be less stable. In the simulation study, we consider values such as $c=\log(3)$, $c=\log(10)$, and $c=\log(30)$, corresponding to increasingly stringent evidence requirements.
The proposed procedure relies on a working model in which detection indicators are treated as approximately independent across perturbation iterations. In practice, bootstrap resampling and stochastic imputation induce dependence, so the resulting evidence measure should be interpreted as an approximate likelihood-ratio score rather than as an exact Bayesian quantity. The calibration of $(\pi_0,\pi_1)$ is empirical and may affect finite-sample behavior; its performance is therefore assessed through simulation and sensitivity analysis. Additional implementation details are provided in Appendix (ref).
The preceding discussion developed the proposed procedure in modular form by introducing the detection rule, the working evidence model, the calibration strategy, and the stopping rule in turn. Algorithm 1 brings these elements together and summarizes the full practical implementation of the BOOT-MI sequential evidence procedure.
This section evaluates the proposed sequential evidence aggregation procedure in a Monte Carlo simulation study. The simulation is designed to address four questions: whether the method improves variable-selection accuracy relative to standard aggregation rules, whether the asymmetric detection rule preserves estimation performance, how sensitive the procedure is to the calibration of $(\pi_0,\pi_1)$ and the threshold $c$, and whether sequential stopping yields computational gains relative to fixed-budget alternatives.
The simulation builds on the data-generating processes of BelloniCH2014a, extended to incorporate missing data and perturbation-based estimation. We generate a synthetic population of size $N=10{,}000$ and draw samples of size $n \in \{100,500,1000\}$. Each scenario is evaluated using 500 Monte Carlo replications.
The main simulation design follows the partially linear model
where $(\varepsilon_i,\nu_i)$ are independent standard normal random variables and $X_i \sim \mathcal{N}(0,\Sigma)$ with $\Sigma_{jl}=0.5^{|j-l|}$. The covariate vector has dimension $p=50$, of which $k_0=5$ variables are relevant. The coefficient vectors satisfy \[ \beta_j = \gamma_j =
\] where $\kappa$ is chosen to achieve target coefficients of determination $R^2 \in \{0.2,0.6\}$. We consider both homoscedastic and heteroscedastic error specifications. Missingness is introduced under MCAR, MAR, and MNAR mechanisms at rates of $20\%$, $40\%$, and $60\%$. As a robustness check, we additionally consider an extended design with group-specific nuisance components in the outcome equation.
Table (ref) summarizes the resulting simulation scenarios. The design spans both favorable and difficult settings for perturbation-based variable selection by varying sample size, signal strength, error structure, missing-data mechanism, and missingness severity. This allows us to assess not only average selection performance, but also the robustness of the proposed calibration and stopping rule under increasingly challenging conditions. Because Design 2 is included as a robustness check rather than as part of the full factorial design, the simulation evidence is driven primarily by the 108 scenarios from Design 1.
We implement the proposed method with pilot calibration ($T_{\mathrm{pilot}}=20$), maximum evidence iterations $T_{\max}=200$, and minimum classification delay $t_{\min}=5$. Unless stated otherwise, the simulation results use the baseline threshold $c=\log(10)$. Sensitivity analyses additionally consider $c\in\{\log(3),\log(10),\log(30)\}$. The baseline implementation uses the stabilized pilot-based calibration described in Section (ref); nominal and permutation-based calibrations are considered only as robustness checks.
We compare the proposed method to three benchmark aggregation rules: the union rule and frequency thresholding at 50% and 75%. Benchmarks are evaluated both at a fixed budget ($T=200$) and at matched budgets corresponding to the stopping time of the proposed method.
We evaluate performance in terms of variable selection accuracy (TPR, FPR, precision, model size), estimation performance (bias, RMSE, coverage), and computational efficiency (number of perturbation iterations until stopping).
We organize the results around variable selection, treatment-effect estimation, and computational efficiency. Throughout, we compare the proposed sequential evidence method to the union rule and to frequency-threshold aggregation at the 50% and 75% levels. For the frequency-based benchmarks, we report both fixed-budget and matched-budget comparisons.
Table (ref) summarizes variable selection performance under the fixed-budget comparison. The union rule achieves perfect sensitivity by construction but selects all variables. The 75% threshold is highly selective but suffers from low sensitivity. The proposed method achieves the highest TPR among the selective methods, at the cost of a higher FPR. The 50% threshold provides a more conservative alternative with lower FPR but reduced sensitivity.
Results are nearly unchanged under the matched-budget comparison in Table (ref), indicating that the gains are not driven solely by longer runs of the benchmark methods.
Figure (ref) shows the aggregate TPR--FPR trade-off under the fixed-budget comparison, while Figure (ref) reports the matched-budget analogue. In both cases, the proposed method occupies a distinctly higher-TPR position than the 50% threshold, while remaining far less inclusive than the union rule.
Figure (ref) shows that the proposed method consistently attains a higher TPR than the 50% threshold across all scenario dimensions. The gain is largest in more difficult settings, especially small samples and MNAR missingness, while the higher FPR reflects the method's more permissive treatment of persistent but weaker signals.
Table (ref) reports treatment-effect performance under the fixed-budget comparison. The proposed method achieves bias and RMSE comparable to the 50% threshold while providing the highest coverage among the selective methods. By contrast, the 75% threshold exhibits the largest bias and lowest coverage, reflecting the inferential cost of omitting relevant controls. Results under the matched-budget comparison are closely similar (Table (ref)).
Table (ref) shows that the proposed stopping rule requires on average 50.7 iterations, compared with the fixed budget of $T=200$ iterations used by the benchmark methods. This corresponds to a reduction of approximately 75% in the number of perturbation iterations.
Figure (ref) illustrates the sequential evidence accumulation for a representative realization. Relevant variables exhibit upward drift and irrelevant variables downward drift, with classification occurring once the paths cross the decision boundaries.
Across all simulation scenarios, the pilot calibration is stable and the fallback rule is not triggered.
The simulation study yields three main findings. First, the proposed method achieves the highest sensitivity among the selective methods and translates this into favorable estimation performance, especially relative to stricter frequency thresholding. Second, the TPR advantage is robust across simulation settings and is most pronounced in difficult scenarios such as small samples, MNAR missingness, and weak signals. Third, the sequential stopping rule substantially reduces computational cost relative to fixed-budget procedures.
Additional disaggregated results by sample size, missingness rate, missing-data mechanism, signal strength, and error structure are reported in Appendix (ref).
We illustrate how the proposed BOOT-MI sequential evidence procedure behaves in a realistic applied setting using data from Round 10 of the European Social Survey (ESS), integrated file, edition 3.3. The empirical analysis is based on an analytic subsample of women constructed from the ESS Round 10 integrated dataset. We restrict the sample to observations with nonmissing values for the employment indicator and the household child indicator; the remaining covariates are allowed to contain missing values and are handled through multiple imputation. After applying the sample restrictions described in Appendix (ref), the resulting analysis sample contains $n=8{,}543$ observations.
Each perturbation iteration consists of bootstrap resampling of the incomplete dataset, stochastic imputation using random forests, and variable selection using the PDS framework. Within each iteration, LASSO is applied to the outcome equation and the auxiliary equation for the variable of interest to construct the PDS candidate set. The asymmetric detection rule from Section (ref) is then applied: variables selected in the auxiliary equation are retained automatically, while variables selected only in the outcome equation must additionally be significant at level $\alpha=0.05$ in the post-selection outcome regression.
The working detection probabilities $(\pi_0,\pi_1)$ are calibrated from a pilot phase with $T_{\mathrm{pilot}}=20$ iterations. In the ESS application, this yields $\hat{\pi}_{0,\mathrm{raw}}=0.075$, a stabilized value $\hat{\pi}_0=0.0687$, and $\hat{\pi}_1=0.9900$, implying a break-even detection frequency of $q^\ast=0.6296$. In the main empirical specification, we use the more conservative threshold $c=\log(1000)$, check threshold crossing only after $t_{\min}=5$ evidence iterations, and impose a maximum of $T_{\max}=200$ evidence iterations. We use a stricter threshold in the empirical application than in the simulation baseline because the goal there is not method comparison but a more conservative, high-evidence variable classification in a single realized dataset. Below, we report sensitivity to less conservative threshold choices.
As benchmarks, we consider three fixed-iteration aggregation rules based on the raw within-iteration PDS union history. Let $S_t^{\mathrm{PDS}} = S_{Y,t} \cup S_{D,t}$ denote the PDS candidate set in perturbation iteration $t$. The Union benchmark includes a variable if it appears in $S_t^{\mathrm{PDS}}$ in at least one iteration, while the Freq(50%) and Freq(75%) benchmarks include a variable if its empirical inclusion frequency across 200 perturbation iterations is at least $0.50$ or $0.75$, respectively.
\paragraph{Variable Selection} Table (ref) reports the variables selected under each aggregation rule. The benchmark methods aggregate repeated raw PDS candidate sets, whereas the proposed method aggregates asymmetric detection indicators through the calibrated sequential evidence process.
The union rule selects all 55 candidate variables and therefore produces the densest specification. The 50% and 75% frequency rules yield more parsimonious models with 41 and 31 variables, respectively. The proposed method selects 34 variables and thus lies between the two frequency-based rules in terms of model size.
Substantively, the proposed method retains a stable core set of controls that is also selected by the benchmark rules, including age, years of education, health, household income, subjective financial situation, unemployment experience, several education and marital-status indicators, and a subset of country indicators. Differences across methods arise primarily for weaker attitudinal variables, additional missing-value indicators, urbanicity indicators, and the breadth of included country effects. Overall, the proposed procedure delivers a data-driven compromise between the inclusiveness of the union rule and the stricter frequency thresholds.
Figure (ref) displays the corresponding log-evidence paths. Variables with sustained support across perturbations cross the upper threshold quickly, while variables with weak or unstable support remain near the continuation region or drift downward. In the baseline ESS specification, all 55 candidate variables are classified within 10 perturbation iterations.
\paragraph{Estimated Coefficient on the Variable of Interest}
Table (ref) reports the estimated coefficient on the variable of interest across aggregation methods. For each rule, the selected control set is fixed first, and the coefficient is then re-estimated across $M=10$ multiply imputed datasets using Rubin's rules.
The estimated coefficient is close to zero and statistically insignificant under all four aggregation rules. Point estimates range from $-0.001$ to $0.002$, standard errors range from $0.011$ to $0.012$, and the confidence intervals overlap almost completely. Thus, in this application, the estimated conditional association is highly robust to the choice of aggregation rule despite substantial differences in model size. A more detailed substantive discussion is provided in Appendix (ref).
\paragraph{Sensitivity Analysis}
Table (ref) summarizes the sensitivity of the empirical results to alternative evidence thresholds and calibration choices.
The empirical findings are qualitatively stable across reasonable changes in the evidence threshold and calibration method. Varying $c$ changes the selected model size only modestly, while larger thresholds mechanically require more iterations. Permutation calibration yields a noticeably more conservative specification because it implies a much larger null detection probability. Additional discussion is provided in Appendix (ref).
This section concludes by summarizing the main contributions of the paper, discussing the limitations of the proposed procedure, and outlining promising directions for future research.
This paper proposes a sequential evidence-aggregation procedure for repeated stochastic variable-selection outcomes generated by bootstrap resampling and stochastic imputation, using PDS as the within-iteration screening device. By modeling detection outcomes across perturbation iterations with a working Bernoulli detection model and accumulating likelihood-ratio–based evidence, the procedure provides a structured inclusion rule together with a data-driven stopping criterion.
The Monte Carlo study evaluates the procedure across 126 simulation scenarios varying sample size, signal strength, missing-data mechanism, and missingness level. The results indicate that the proposed method achieves a favorable trade-off between true positive and false positive rates relative to the union rule and frequency-based aggregation, while producing sparser models and competitive treatment effect estimates. The sequential stopping rule substantially reduces the number of perturbation iterations compared to fixed-iteration procedures.
Several limitations should be noted. First, the Bernoulli working model treats detection indicators as conditionally independent across perturbation iterations, whereas bootstrap resampling and multiple imputation from the same observed sample induce positive dependence. The simulation study provides evidence on the practical consequences of this approximation, but the procedure should not be interpreted as delivering exact Bayesian posterior probabilities. Second, the calibration of the detection probabilities $(\hat{\pi}_0,\hat{\pi}_1)$ relies on a short pilot phase and heuristic estimators. While the simulation results suggest that the procedure is robust across the calibration strategies considered, the sensitivity of the evidence process to pilot length and calibration method warrants further investigation. Third, the procedure is evaluated here in a linear model setting with LASSO-based selection; its performance in nonlinear or high-dimensional nonparametric settings remains to be studied.
The evidence framework is not restricted to LASSO-based variable selection. It applies to any procedure that produces binary detection indicators with different detection probabilities under relevance and irrelevance. For methods with standard errors, the detection event can be defined via a significance threshold. For machine learning methods without standard errors, detection events can be defined through selection stability, with calibration based on empirical detection frequencies or permutation of the outcome variable to estimate the null detection rate. This generality suggests several directions for future work, including the application of the framework to random forests and gradient boosting methods, the development of adaptive calibration strategies that update $(\hat{\pi}_0,\hat{\pi}_1)$ during the sequential run, and the extension to settings with multiple treatment variables or high-dimensional treatment effect heterogeneity.
An important direction for future work is fuller uncertainty propagation across the entire pipeline, for example by embedding the complete perturbation-selection-estimation procedure in an outer resampling scheme. Such an extension would be particularly relevant for causal applications in which uncertainty from model selection is itself substantively important.
The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 35/1597-1 FUGG. The authors gratefully acknowledge the support of the HERMES network (Higher Education and Research in Management of European Universities) for a travel grant and for facilitating this study through its network.