Reproducible Aggregation of Sample-Split Statistics

titlepage\begin{spacing}{1} \begin{abstract} \smalltonormalsize{Statistical inference is often simplified by sample-splitting. This simplification comes at the cost of the introduction of randomness not native to the data. We propose a simple procedure for sequentially aggregating statistics constructed with multiple splits of the same sample. The user specifies a bound and a nominal error rate. If the procedure is implemented twice on the same data, the nominal error rate approximates the chance that the results differ by more than the bound. We illustrate the application of the procedure to several widely applied econometric methods.} \\ \\ Keywords: Sample-splitting, Cross-Fitting, Cross-Validation, Reproducibility \\ JEL: C01, C13, C52 \end{abstract} \end{spacing}

spacing{1.3} \section{Introduction} Sample-splitting is ubiquitous in modern econometric theory. Routine statistical tasks---model selection, dimension reduction, nuisance parameter estimation---can be implemented on a randomly selected subsample of a data set, without contaminating the validity of a statistical inference produced on its complement. This principle underlies the widely applied practices of cross-validation for predictive risk estimation stone1974cross,arlot2010survey and cross-fitting for adaptive estimation of semiparametric models bickel1982adaptive,schick1986asymptotically,chernozhukov2018double, among a growing set of additional applications. For a fixed data set, statistics constructed with sample-splitting are not deterministic. Two researchers can compute the same statistic on the same data and obtain different values. Researchers are incentivized to report significant results. If there is scope to materially alter the statistics that they report through the choice of the split of their sample, should this choice be left to chance? This paper makes two contributions. First, we show that many widely applied sample-split econometric methods exhibit significant residual randomness. We give examples from the applied economics literature where the randomness induced by sample-splitting determines the statistical significance of a treatment effect estimated with cross-fitting, the interpretation of a model selected with cross-validation, and the qualitative features of treatment targeting rules learned from cross-fit nuisance parameter estimates. Second, and accordingly, we propose an efficient method for removing the residual randomness from sample-split statistics. The procedure takes as input a bound and an error rate. The statistic of interest is sequentially aggregated over randomly drawn splits of the sample. The procedure is stopped after an estimate of the residual variation of the aggregate statistic falls below a pre-determined threshold. If the procedure were run twice, we show that the chance that the outputs differ by more than the bound is well-approximated by the error rate. That is, by setting the bound and error rate to be sufficiently small, sample-split statistics aggregated with the procedure are reproducible. We begin, in (ref), by discussing several widely applied sample-split econometric methods and demonstrating that, for each, the randomness induced by sample-splitting can substantively effect results. In (ref), we propose an efficient method for sequentially aggregating sample-split statistics that ensures that the residual randomness is small. We illustrate that, in each of our examples, the proposed method stabilizes results---ensuring reproducibility at a minimal computational expense. We establish that the procedure is valid, in a particular asymptotic sense, without imposing any restrictions on the data generating process or statistic of interest. Similarly, we show that, for a large class of applications, sample-split statistics aggregated with the procedure maintain (or, improve upon) unconditional statistical guarantees. That is, reproducibly aggregated sample-split statistics are still consistent and associated inferences are still valid. To implement the procedure, a user must make several choices that may affect performance, including the specification of a suitable bound and error rate. Additionally, in most applications, the procedure is applied to stabilize a statistic that is itself constructed with cross-splitting. In these cases, users must also specify how many folds to use for cross-splitting. To shed light on these choices, in (ref), we give an analysis of the performance of the procedure under a set of simplifying conditions. We give two main results. First, the computation needed to achieve a given bound on residual randomness is very sensitive to the desired error tolerance, but is insensitive to the number of folds that are used for cross-splitting. Second, and on the other hand, the accuracy of the nominal error rate of the procedure deteriorates as the number of folds used for cross-splitting increases. We conclude, in (ref), by synthesizing these results into a set of concrete recommendations for practice. We emphasize simple rules-of-thumb for choosing suitable error tolerances and approaches to cross-splitting. The main theoretical challenge posed by this analysis is the accommodation of the dependence between statistics computed on cross-splits of a sample. We address this through an application of the method of exchangeable pairs stein1986approximate, ross2011fundamentals,chen2011normal. To construct an appropriate exchangeable pair for our problem, we develop a novel application of a coupling argument due to chatterjee2005concentration, that may be of independent interest. \subsection{Related Literature} The procedure studied in this paper is applicable to a large variety of sample-split statistical methods. In (ref), we give a selective review of various sample-split methods that are frequently used in applied economics. There are many, additional, sample-split methods, proposed in the statistics literature, that have the potential to be useful in econometric applications. Generic sample-split methods for testing statistical hypotheses are studied in guo2017analysis, diciccio2020exact, and wasserman2020universal. Additional applications include sample-split procedures for selective inference rinaldo2019bootstrapping, inference on high-dimensional linear models meinshausen2010stability, conformal and predictive inference lei2018distribution, and knockoff tests of conditional independence barber2015controlling. chernozhukov2018double and chernozhukov2018generic advocate for the aggregation of estimators and $p$-values computed with sample splitting, over a pre-determined number of splits, in the context of applications related to semiparametric estimation and characterization of treatment effect heterogeneity, respectively. We second these recommendations and contribute a general purpose method that provides a statistical guarantee that residual randomness has been controlled up to a specified level of error, at a minimal computational cost. Our setting is related to a large literature that studies methods for constructing confidence intervals for cross-validated estimates of generalization error. Several examples include dietterich1998approximate, nadeau1999inference, lei2020cross, bayle2020cross, austern2020asymptotics, and bates2023cross. By contrast, we are interested in the randomness conditional on the data. Formally, our non-asymptotic results are most similar to the Berry-Esseen bounds given in austern2020asymptotics, who study the unconditional normal approximation of statistics similar to the those considered in (ref). Our bounds apply under weaker conditions and are substantially simpler. Some of our results build on a literature studying the role of algorithmic stability in the accuracy of cross-validation kale2011cross,kumar2013near. These papers are related to a broader literature that derives generalization bounds for stable algorithms, originating with bousquet2002stability. Some of the concentration inequalities that we derive can be compared to the results of cornec2010concentration and abou2019exponential. Again, the setting we study is different and our conditions, and resultant bounds, are substantially simpler. Although our emphasis is on statistics constructed with sample-splitting, the algorithmic and formal methods studied in this paper are potentially applicable to randomized algorithms more generally motwani1995randomized. See beran1987stochastic for a classical analysis of the asymptotics of randomized tests and estimators. \section{The Residual Randomness of Sample-Split Statistics} Sample-splitting has proliferated as a useful subroutine for simplifying various tasks associated with modern statistical inference. The high-level situation is as follows. Consider a researcher who observes the independent data $D = (D_i)_{i=1}^n$ and wishes to report the statistic \begin{equation} \Psi(D, \eta) , \end{equation} where $\eta$ is some unknown nuisance parameter. For example, each observation $D_i$ could contain a measurement of an outcome $Y_i$, a treatment $W_i$, and a vector collecting a large number of covariates $X_i$. The researcher could be interested in measuring the effect of $W_i$ on $Y_i$. To ensure that their estimate is not unnecessarily imprecise, they might like to include controls for only the subset of covariates that are correlated with the outcome. Here, formally, the nuisance parameter $\eta$ collects the indices of the subset of “relevant” controls and the statistic $\Psi(D, \eta)$ denotes a treatment effect estimate associated with, say, a regression of the outcome on the treatment with controls for covariates with indices in $\eta$. In practice, the researcher cannot compute the statistic (ref), as the nuisance parameter $\eta$ is unknown. Often, however, an estimator $\hat{\eta}(D)$ is available and so it may be tempting to report the feasible statistic \begin{equation} \Psi(D, \hat{\eta}(D)) . \end{equation} In the example, the estimator $\hat{\eta}(D)$ could collect the indices of the covariates whose absolute sample correlation with the outcome is greater than some pre-determined threshold. Statistical inferences that, counterfactually, would be appropriate if the infeasible statistic (ref) were available will not necessarily be valid if they were instead based on the feasible statistic (ref). In particular, the process of computing the nuisance parameter estimate $\hat{\eta}(D)$ might produce a confounding effect. In the example, screening control variables according to their correlation with the outcome produces a bias in the resultant treatment effect estimate. Sample-splitting solves this problem.\footnote{Of course, under certain conditions involving the sparsity of the covariance between the treatment, outcome, and controls, “double post selection” of control variables with a Lasso regression tibshirani1996regression can produce consistent treatment effect estimates belloni2014inference. However, consistent estimates can be obtained under weaker conditions through a closely related approach based on sample-splitting chernozhukov2018double,belloni2012sparse.} Let $\mathsf{s}$ denote a randomly drawn subset of the numbers $[n]=\{1,\ldots,n\}$ of size $b$ and let $\tilde{\mathsf{s}}$ denote its complement. Sample-split statistics take the form \begin{equation} T(\mathsf{s}, D) = \Psi(D_{\mathsf{s}}, \hat{\eta}(D_{\tilde{\mathsf{s}}})) , \end{equation} where the quantities $D_{\mathsf{s}} = (D_i)_{i\in\mathsf{s}}$ and $D_{\tilde{\mathsf{s}}} = (D_i)_{i\in\tilde{\mathsf{s}}}$ collect the data with indices in $\mathsf{s}$ and $\tilde{\mathsf{s}}$, respectively. Splitting the data $D$ into two independent subsamples $(D_{\mathsf{s}}, D_{\tilde{\mathsf{s}}})$ prevents the construction of nuisance parameter estimates from contaminating statistical inferences that the researcher may wish to make with the feasible statistic $\eqref{eq: T def intro}$. There are two, well-known, practical issues with this approach. First, by splitting the data, statistical precision may be meaningfully reduced. Second, the statistic (ref) is random through both the data $D$ and the choice of the random subset $\mathsf{s}$. That is, if the same sample-split statistic were computed on the same data by two different researchers, and the statistic was sensitive to the choice of the random subset $\mathsf{s}$, then the researchers could report meaningfully different results. To address these concerns, researchers often aggregate several replications of sample-split statistics through cross-splitting stone1974cross,schick1986asymptotically. In particular, researchers typically report aggregate statistics of the form \begin{equation} a(\mathsf{r}, D) = \frac{1}{k} \sum_{j=1}^k T(\mathsf{s}_j, D) , \end{equation} where the quantity $\mathsf{r} = (\mathsf{s}_j)^k_{j=1}$ denotes a random $k$-fold partition of $[n]$, i.e., a collection of $k$ mutually exclusive sets whose union is equal to $[n]$. By reusing subsamples of the data for computation of both the statistic of interest and the estimation of nuisance parameters, cross-split statistics mitigate potential losses in statistical precision. In this section, we demonstrate, in several real applications from applied economics, that the second concern---the residual randomness generated by sample-splitting---can substantively affect results. This residual variability is, often, not resolved by cross-splitting. Consequently, in (ref), we propose an approach to aggregating sample-split statistics that controls the scope of residual randomness at a minimal computational cost. Revisiting the applications, we demonstrate that sample-split statistics aggregated through this procedure are reproducible. \subsection{Cross-Validation} The most prevalent instance of sample-splitting in applied economics is the use of cross-validation for model selection. Often, in this setting, the data $D_i$ consist of a measurement of an outcome $Y_i$ and a vector $X_i$ collecting measurements of $p$ covariates. Interest is in choosing a parsimonious subset of the covariates that, together, best predicts the outcome. Lasso regression is a standard approach to this problem tibshirani1996regression,hastie2015statistical. The Lasso coefficient is the solution to the regularized least-squares regression \begin{equation} \hat{\eta}_{\lambda}(D) = \underset{\eta \in \mathbb{R}^p}{\arg\min} \left\{ \frac{1}{n} \sum_{i=1}^n (Y_i - \eta^\top X_i)^2 + \lambda \sum_{j=1}^p | \eta_j \vert \right\} , \end{equation} where $\lambda$ is some tuning parameter chosen by the user. Penalization of the $\ell_1$-norm of the coefficient $\eta$ encourages sparsity. That is, often, many elements of $\hat{\eta}_{\lambda}(D)$ are exactly equal to zero. The set of covariates associated with non-zero coefficients are referred to as the model “selected” by the Lasso. The selected model is sensitive to the choice of the tuning parameter $\lambda$. If $\lambda$ is sufficiently large, no covariates are selected. If $\lambda$ is sufficiently small, all covariates are selected.\footnote{For sufficiently small $\lambda$, a covariate is selected so long as it is not collinear with other covariates and its sample correlation with the outcome is not exactly equal to zero.} Cross-validation is a widely applied and, perhaps, uncontroversial approach for choosing tuning parameters. Here, cross-validation entails assigning $\lambda$ the value that minimizes the cross-split estimate of the out-of-sample mean-squared error \begin{equation} a_{\lambda}(\mathsf{r}, D) = \frac{1}{k} \sum_{j=1}^k T_{\lambda}(\mathsf{s}_j, D) , \quad where \quad T_{\lambda}(\mathsf{s},D)=\frac{1}{\vert \mathsf{s} \vert} \sum_{i\in\mathsf{s}} (Y_i - \hat{\eta}_{\lambda}(D_{\tilde{\mathsf{s}}})^\top X_i)^2 \end{equation} and, as before, $\mathsf{r}=(\mathsf{s}_j)^k_{j=1}$ is a random $k$-fold partition of $[n]$. By measuring error out-of-sample, i.e., in the independent subsample $D_{\mathsf{s}}$, the cross-validated risk estimate (ref) avoids any “over-fitting” bias induced by the estimation of the coefficient $\hat{\eta}_{\lambda}(D_{\tilde{\mathsf{s}}})$. This approach is widely applied throughout various subfields of applied economics.\footnote{See, for example, applications to Crime mastrobuoni2020crime, arnold2020measuring, Development casey2021experiment,blattman2024gang,sadka2024information, Economic Theory fudenberg2019predicting, Education ellison2021efficiency, Environment deryugina2019mortality,cicala2022imperfect, Finance koijen2024investors, Health abaluck2016determinants,cooper2020surprise, Industrial Organization kelly2023mechanics,dube2023personalized, Innovation chen2021notching,myers2022estimating, Labor muendler2010margins,card2020referees,adermon2021dynastic,derenoncourt2022can, Market Design agarwal2019market, Macroeconomics hansen2018transparency, and Political Economy gentzkow2019measuring,cantoni2022does.} The risk estimate (ref) is random both through the data $D$ and through the choice of the $k$-fold partition $\mathsf{r}$. That is, the choice of $\lambda$, and thereby, the selected model, have the potential to change for different choices of the random collection $\mathsf{r}$. To evaluate the scope of this sensitivity, we consider data from casey2021experiment. casey2021experiment study a large-scale experiment, implemented in Sierra Leone, in which a randomly selected subset of parliamentary elections were preceded by direct vote, party-specific primaries. They select covariates to include in various, downstream, econometric analyses with the cross-validated Lasso. In their setting, the outcome $Y_i$ is the vote share, in a poll of party officials, for each of 390 candidates. The vector $X_i$ collects measurements of 48 characteristics for each candidate. \begin{figure}[t] \begin{centering} \caption{Cross-Validation} \begin{tabular}{c} Panel A: Mean-Squared Error Quantiles\tabularnewline \tabularnewline Panel B: Model Variability\tabularnewline \tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) measures the residual randomness of the cross-validated Lasso, implemented in data from casey2021experiment. Panel A displays quantiles of the $10$-fold cross-validated estimate of the mean-squared error (ref) of the lasso regression (ref) over a grid of values of $\lambda$. In Panel B, the probabilities that each value of $\lambda$ minimize the cross-validated risk estimate are displayed with light green dots, relative to the left $y$-axis. The number of covariates selected at each value of $\lambda$ are displayed with dark blue dots, relative to the right $y$-axis. See (ref) for further details.} \hrulefill \end{figure} Panel A of (ref) displays quantiles, across random draws of the $10$-fold partition $\mathsf{r}$, of the cross-validated estimate of the mean-squared error (ref) over a grid of values of $\lambda$. The curve associated with the 5th quantile has a minimum around $\lambda = 0.02$. By contrast, the curve associated with the 95th quantile is monotonically decreasing. This induces instability in the value of $\lambda$ chosen by cross-validation. Panel B of (ref) displays, in light green, the probability---again, across random $10$-fold partitions $\mathsf{r}$---that each value of $\lambda$ is selected, i.e., minimizes the cross-validated estimate of the mean-squared error. The number of covariates in the model associated with each value of $\lambda$ is displayed in dark blue. The distribution of the selection probabilities has two maxima---both are associated with probabilities greater than 0.2. One entails selecting a model with 14 covariates. The other entails selecting a model with zero covariates. The random choice of the partition $\mathsf{r}$ substantively affects the size of the selected model.\footnote{casey2021experiment take measures to stabilize their results. In particular, they select covariates that are selected in more than 200 of 400 randomly drawn 10-fold partitions $\mathsf{r}$. We use this setting as an example, in part, because a serious attempt was made to address residual randomness. This is not commonplace in the literature surveyed in (ref) (a similar strategy is used in, e.g., hansen2018transparency, however).} \subsection{Cross-Fitting}Meaningful residual randomness is not particular to cross-validated risk estimation. A second class of sample-split methods, seeing increasing use in applied economics, are characterized by the application of “cross-fitting” to accommodate---or characterize---treatment effect heterogeneity. Often, in this case, the data $D_i$ consist of an outcome $Y_i$, a binary treatment $W_i$, and a vector of covariates $X_i$. Let $Y_i(1)$ and $Y_i(0)$ denote the potential outcomes induced by the treatment $W_i$. As before, let $\mathsf{s}$ denote a random a subset of $[n]$, with complement $\tilde{\mathsf{s}}$. The basic idea underlying this class of methods is to use the data from the units $i$ in $\tilde{\mathsf{s}}$ to predict the treatment effects $Y_i(1)-Y_i(0)$ for the units $i$ in $\mathsf{s}$. For example, a machine learning algorithm can be applied to the data $D_{\tilde{\mathsf{s}}}$ to construct an estimate of the conditional expectation \begin{equation} \mu_w(x) = \mathbb{E}[Y_i \mid X_i = x, W_i = w] \end{equation} for each $w$ in $\{0,1\}$. Collect these estimates into $\hat{\eta}(D_{\tilde{\mathsf{s}}}) = (\hat{\mu}_1, \hat{\mu}_0)$. Predictions of the treatment effects for the units $i$ in $\mathsf{s}$ can be constructed through the sample-split statistic \begin{equation} \psi(D_i, \hat{\eta}(D_{\tilde{\mathsf{s}}})) = \hat{\mu}_1(X_i) - \hat{\mu}_0(X_i) . \end{equation} These predictions can then be used to estimate average treatment effects and characterize treatment effect heterogeneity, among other related problems. The point is, by splitting the data, these “second-stage” analyses are not confounded by the construction of the “first-stage” estimators. “Double Machine Learning” (DML) estimates of average treatment effects are a leading example of a method that takes this structure chernozhukov2018double. Here, the sample-split statistics (ref) are aggregated with cross-splitting, through \begin{equation} a(\mathsf{r},D)= \frac{1}{k}\sum_{j=1}^k T(\mathsf{s}_j, D) , \quadwhere\quad T(\mathsf{s}, D) = \frac{1}{\vert \mathsf{s}\vert}\sum_{i\in\mathsf{s}} \psi(D_i, \hat{\eta}(D_{\tilde{\mathsf{s}}})) \end{equation} and, again, $\mathsf{r}$ is a random $k$-fold partition of $[n]$. An appropriate standard error for (ref) is itself given by the cross-split statistic \begin{equation} \mathsf{se}(\mathsf{r},D) = \frac{1}{n} \sqrt{ \sum_{j=1}^k \sum_{i \in \mathsf{s}_j} \left(\psi(D_i, \hat{\eta}(D_{\tilde{\mathsf{s}}_j})) - a(\mathsf{r},D) \right)^2 } . \end{equation} chernozhukov2018double show that an asymptotically efficient test of the null hypothesis that an average effect treatment effect is less than zero can be constructed by comparing (ref) to the critical value $\mathsf{cv}_{\alpha}(\mathsf{r},D) = z_{1-\alpha} \cdot \mathsf{se}(\mathsf{r},D)$, where $z_{1-\alpha}$ is the $1-\alpha$ quantile of the standard normal distribution.\footnote{This result does not apply to estimates of the form (ref), but rather to estimators constructed with a “Neyman Orthogonal” moment, which, in this case, additionally require a non-parametric estimate of the propensity score.} DML estimators have been increasingly used in applied economics, as they place weaker restrictions on treatment effect heterogeneity than, say, methods based on linear regression.\footnote{See, for example, applications in okunogbe2022technology, beraja2023ai, covert2023relinquishing, delfino2024breaking, farronato2024consumer.} The DML estimator (ref) and standard error (ref) are, again, random through both the data $D$ and the $k$-fold partition $\mathsf{r}$. To evaluate the extent of resultant variability, we consider data from chakravorty2024can. chakravorty2024can use DML to study the effect of a program, implemented in two Indian states, involving the provision of information concerning prospective jobs to vocational trainees, on employment outcomes. Here, the binary outcome $Y_i$ indicates employment five months after completing training for each of 890 trainees placed into jobs, $W_i$ denotes assignment to the program, and $X_i$ collects measurements of 77 pre-treatment covariates. Panel A of (ref) displays a heat map of the joint distribution of the $5$-fold cross-fit estimator (ref) and the associated critical value $\mathsf{cv}_{\alpha}(\mathsf{r},D)$ over random draws of the partition $\mathsf{r}$.\footnote{We follow the replication package associated with chakravorty2024can. Nuisance parameters are estimated with random forests using the “Ranger” R package wright2017ranger. Estimates and standard errors are constructed using the “DoubleML” R package bach2024double. See (ref) for further details.} A black line is placed at the threshold where the estimate is equal to the critical value. The residual variability in the estimate is large relative to estimates of the sampling variability, and is sufficient to determine the purported statistical significance. Concretely, the difference between the 5th and 95th quantiles of the distribution of the estimator (0.094 and 0.131, respectively) is equal to 68% of the median of the distribution of the standard error (0.054).\footnote{Some papers take measures to address analogous residual randomness (see e.g., covert2023relinquishing). In fact, chernozhukov2018double suggest taking the median over several partitions $\mathsf{r}$. Most implementations of DML in standard statistical software allow the user to aggregate estimates, although this aggregation does not occur by default (e.g., bach2024double in R or ahrens2024ddml in STATA).} \begin{figure}[p] \begin{centering} \caption{Cross-Fitting} \begin{tabular}{cc} \multicolumn{2}{c}{Panel A: Treatment Effect Estimation}\tabularnewline \multicolumn{2}{c}\tabularnewline \multicolumn{2}{c}{Panel B: Policy Evaluation}\tabularnewline \multicolumn{2}{c}\tabularnewline \multicolumn{2}{c}{\textit{Panel C: Testing for Treatment Effect Heterogeneity}}\tabularnewline \multicolumn{2}{c}\tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) displays discretized heat maps quantifying the residual randomness of three estimators constructed with cross-fitting. All three panels give the joint distribution of an estimator and an associated critical value. In each case, a black line has been placed at the threshold where the estimator is equal to the critical value. Panel A displays the distribution of the DML estimate (ref) cross $5$-fold cross-splits using data from chakravorty2024can. Panel B displays the difference between treatment effects estimates for the “most impacted” and “most deprived” individuals, following haushofer2022targeting, using data from egger2022general, again across $5$-fold cross-splits. Panel C displays the distribution of the estimate associated with a test of treatment effect heterogeneity using data from beaman2023selection. Here, each estimate and critical value is computed by averaging over 250 independently drawn half-samples $\mathsf{s}$. See (ref) for further details on the construction of each panel.} \end{figure} The same behavior is exhibited in related applications that use cross-fitting in more complicated ways. In these cases, measures to stabilize results are more common, although there is little guidance on how to best operationalize this stabilization. To illustrate this, we consider examples from haushofer2022targeting and beaman2023selection.\footnote{Additional, related, methods that uses cross-fitting to estimate and evaluate treatment targeting rules are proposed by athey2021policy and yadlowsky2024evaluating.} haushofer2022targeting use data from a randomized cash transfer implemented in Kenya. These data were originally considered in egger2022general. beaman2023selection use data from an experiment concerning agricultural lending in Mali. Both papers use sample-splitting to construct estimates, of the form (ref), of the treatment effects $Y_i(1) - Y_i(0)$ for each of the units $i$ in $\mathsf{s}$. haushofer2022targeting additionally construct sample-split estimates of the untreated outcome $Y_i(0)$ for each of the units in $\mathsf{s}$. haushofer2022targeting identify the 50% of units in $\mathsf{s}$ that have the largest predicted treatment effect as well as the 50% of units that have the smallest predicted untreated outcome. They refer to these groups as the “most impacted” and “most deprived,” respectively. They estimate the difference in the average treatment effect for the two groups, and construct a standard error, and associated critical value, for this difference with the bootstrap. See (ref) for details. Panel B of (ref) displays a heat map of the joint distribution of the $5$-fold cross-fit estimate of the difference between the treatment effect estimates for the two groups, and the associated critical value, over random draws of the partition $\mathsf{r}$.\footnote{We were not able to access a replication package associated with haushofer2022targeting, and so implement a simplified version of the exercise considered in that paper using data from the replication package associated with egger2022general. Treatment effects, and untreated outcomes, are estimated with random forests using the “GRF” R package athey2019generalized. Analogous estimates reported in haushofer2022targeting are statistically significant. We emphasize that, as we implement only a simplified version of the exercise conducted in haushofer2022targeting, our estimates should only be interpreted an an illustration of the scope of residual randomness in analyses of this type, rather than as a substantive characterization of underlying treatment effect heterogeneity. In particular, we make no attempt to reweigh observations according to their sampling probabilities and appear to be using a different measure of the time between the administration of the experiment and the measurement of post-treatment outcomes.} The residual randomness is considerable. Here, the difference between the 5th and 95th quantiles of the distribution of the estimator (37.82 and 109.69, respectively) is equal to 190% of the median of the distribution of the standard error (54.44). To stabilize these results, haushofer2022targeting average over 400 draws of the $5$-fold partition $\mathsf{r}$. Due to this increased computation, they only report confidence intervals for their main results. The methods developed in this paper ensure that stabilization of sample-split statistics, in this way, occurs at a minimal computational expense. Likewise, beaman2023selection implement a test of treatment effect heterogeneity proposed by chernozhukov2018generic. In particular, in data for units in $\mathsf{s}$, the outcome is regressed on the treatment, the treatment effect estimate, and an interaction between the treatment and the treatment effect estimate. The idea is that, if treatment effect estimates are well-calibrated, then the coefficient on the interaction should be statistically greater than zero. chernozhukov2018generic recommend aggregating $p$-values associated with the coefficients on the interaction over 250 replications of this sample-split test. Panel C of (ref) displays a heat map measuring the joint distribution of the average coefficient and critical value associated with this test.\footnote{We follow the details of the replication package associated with beaman2023selection. Nuisance parameters are estimated with random forests using the “GRF” R package athey2019generalized. Further details are given in (ref).} That is, each coefficient and critical value is computed by taking the average over 250 sample-splits. Despite this aggregation, the residual randomness remains meaningful.\footnote{chernozhukov2018generic recommend aggregating $p$-values the median to ensure robustness to outliers. In practice, if aggregation of a sample-split statistic with a median, rather than a mean, makes a material difference, we strongly encourage researcher to investigate why some sample-splits generate extreme (or highly skewed) estimates (e.g., outliers in the underlying data or poor overlap in an intervention).} The difference between the 5th and 95th quantiles of the distribution of the estimator (0.86 and 1.08, respectively) is equal to 28% of the median of the distribution of the standard error (0.80). Perhaps motivated by this instability, beaman2023selection aggregate over 1000 sample-splits. Two themes emerge from these examples. First, sample-splitting is a versatile tool for simplifying various tasks associated with modern statistical inference. Second, the residual randomness induced by sample splitting is often large and can affect the substantive interpretation of results. In the next section, we provide a general-purpose method for aggregating sample-split statistics that ensures that residual randomness is controlled at a minimal computational cost. \section{Reproducible Aggregation} We propose a sequential method for aggregating sample-split statistics. Our objective is to ensure that the auxiliary randomness induced by sample-splitting is small. To introduce the method, we require some additional notation. The set $\mathcal{S}_{n,b}$ consists of all subsets of $\left[n\right]=\{1,\ldots,n\}$ of size $b$. In turn, the set $\mathcal{R}_{n,k,b}$ contains all collections of $k$ mutually exclusive elements of $\mathcal{S}_{n,b}$. That is, if $n=k\cdot b$, then the set $\mathcal{R}_{n,k,b}$ collects all partitions of $[n]$ into $k$ mutually exclusive sets of size $b$. We refer to the elements of the set $\mathcal{R}_{n,k,b}$ as cross-splits. Throughout, the quantity $\mathsf{R}_{g,k}=(\mathsf{r}_i)_{i=1}^g$ collects $g$ cross-splits, each given by $\mathsf{r}_i = (\mathsf{s}_{i,j})_{j=1}^k$. Unless otherwise specified, the elements of $\mathsf{R}_{g,k}$ are random, sampled independently and uniformly from the collection $\mathcal{R}_{n,k,b}$. Suppose that we are interested in some sample-split statistic $T(\mathsf{s},D)$. We study the construction of aggregate statistics of the form \begin{equation} a(\mathsf{R}_{g,k}, D)= \frac{1}{g} \sum^g_{i=1} a(\mathsf{r_i}, D) , \quad\text{where}\quad a(\mathsf{r_i}, D) = \mathcal{A}(\{T(\mathsf{s}_{i,j},D)\}_{j=1}^k) \end{equation} and the function $\mathcal{A}(\cdot)$ aggregates the statistics $T(\mathsf{s}_{i,j},D)$ across the cross-split $\mathsf{r}_i$. To simplify exposition, for the time being, we will restrict attention to the case that the statistic $a(\mathsf{R}_{g,k}, D)$ is real-valued. This is natural, if, for example, the statistic $a(\mathsf{r}, D)$ is a treatment effect estimate or $p$-value constructed with cross-fitting. Later, in our application to cross-validated risk estimation, we treat vector-valued sample-split statistics, e.g., cross-validated risk estimates queried at a vector of values of a penalization parameter. Our task is to formulate a method for choosing the number of cross-splits $g$ to ensure that the residual variability of the aggregate statistic (ref) is small. We formalize this objective as follows. \begin{defn}[Reproducible Aggregation] Let the sequences $\{\mathsf{r}_i\}_{i=1}^\infty$ and $\{\mathsf{r}^\prime_{i}\}_{i=1}^\infty$ be drawn independently and uniformly, conditional on the data $D$, from the collection of cross-splits $\mathcal{R}_{n,k,b}$. Define $\mathsf{R}_{g,k} = \{\mathsf{r}_g\}_{i=1}^g $ and $\mathsf{R}^{\prime}_{g,k} = \{\mathsf{r}^{\prime}_g\}_{i=1}^g$ for each integer $g$. Suppose that the integers $\hat{g}$ and $\hat{g}^\prime$ are independent and identically distributed, again conditional on the data $D$. We say that the aggregate statistic $a(\mathsf{R}_{\hat{g},k}, D)$ is $(\xi,\beta)$-reproducible if \begin{equation} P\bigg\{ \big\vert a(\mathsf{R}_{\hat{g},k}, D) - a(\mathsf{R}^\prime_{\hat{g}^\prime,k}, D) \big\vert \geq \xi \mid D\bigg\} \leq \beta \end{equation} almost surely. \end{defn} (ref) is motivated by the following thought experiment. Suppose that the data $D$ are given to two researchers. Each researcher is tasked with producing an estimate of the form (ref). They generate the collections of splits, $\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^\prime_{\hat{g}^\prime,k}$, independently using the same, potentially data-dependent, procedure. That is, the two collections of splits are independent and identically distributed, conditional on the data $D$. If the estimate $a(\mathsf{R}_{\hat{g},k}, D)$ is $(\xi,\beta)$-reproducible, then the probability that the two researchers' estimates differ by more than $\xi$ is less than $\beta$. Of course, ensuring that estimates are reproducible, in the sense of (ref), does not preclude deceptive behavior. Researchers might compute a reproducibility aggregated statistic many times, until a desirable estimate is obtained. However, reporting reproducible statistics greatly increases the cost of searches of this form, as the scope of residual randomness has been reduced. We propose a sequential method for constructing reproducible sample-split statistics. The proposal is based on the fixed-length sequential confidence intervals of anscombe1952large and chow1965asymptotic. The procedure works by repeatedly drawing a cross-split $\mathsf{r}_g$ uniformly from $\mathcal{R}_{n,k,b}$, appending the cross-split to the collection $\mathsf{R}_{g,k} = (\mathsf{R}_{g-1,k}, \mathsf{r}_g)$, and estimating the conditional variance \begin{flalign} v_{g,k}\left(D\right) = \Var( a(\mathsf{R}_{g,k}, D) \mid D) \end{flalign} with the plug-in estimator \begin{flalign} \hat{v}\left(\mathsf{R}_{g,k}, D\right) & =\frac{1}{g}\frac{1}{g-1}\sum_{i=1}^{g} (a(\mathsf{r}_i,D) - a(\mathsf{R}_{g,k},D))^2 \end{flalign} until the condition \begin{equation} \hat{v}\left(\mathsf{R}_{g,k}, D\right) \leq \mathsf{cv}(\xi,\beta) = \frac{1}{2}\left(\frac{\xi}{z_{1-\beta/2}}\right)^2 \end{equation} is satisfied, where $z_{\alpha}$ denotes the $\alpha$ quantile of the standard normal distribution. In particular, let $g_{\mathsf{init}}\geq2$ denote some “burn-in” period chosen by the user. The number of cross-splits $\hat{g}$ chosen by the procedure is the smallest value of $g$, greater than $g_{\mathsf{init}}$, such that (ref) is satisfied. The procedure is summarized in (ref). \begin{algorithm}[t] \caption{Anscombe-Chow-Robbins Aggregation} \KwIn{Data $D$, tolerance $\xi$, error rate $\beta$, collection size $k$, split size $b$, initialization $g_{\mathsf{init}}$} Set $g \leftarrow g_{\mathsf{init}}$ Draw $\mathsf{r}_1$, ..., $\mathsf{r}_{g_{\mathsf{init}}}$ independently and uniformly from $\mathcal{R}_{n,k,b}$. Collect $\mathsf{R}_{g_{\mathsf{init}},k} = ( \mathsf{r}_j)_{j=1}^{g_{\mathsf{init}}}$. \While{$\hat{v}(\mathsf{R}_{g,k}, D) > \mathsf{cv}(\xi,\beta) $}{ Set $g \leftarrow g + 1$ Draw $\mathsf{r}_g$ uniformly from $\mathcal{R}_{n,k,b}$. Collect $\mathsf{R}_{g,k} = (\mathsf{R}_{g-1,k}, \mathsf{r}_g)$. } Set $\hat{g} \leftarrow g$ \Return $a(\mathsf{R}_{\hat{g},k}, D)$ \nonl \hrulefill \nonl {Notes: (ref) gives a method for sequentially aggregating sample-split statistics. The critical value $\mathsf{cv}(\xi,\beta)$ is defined in display (ref).} \end{algorithm} In (ref), we show that (ref) is applicable, off-the-shelf, to a large class of problems. In particular, without imposing any conditions on the data generating process or the statistic of interest, we show that sample-split statistics aggregated with (ref) are reproducible, in a particular asymptotic sense. Moreover, we show that, in many applications, statistics aggregated with (ref) maintain, or improve on, various unconditional statistical guarantees. We then revisit, in (ref), the examples considered in (ref). We show that, in each case, an application of (ref) produces a reproducible, stabilized, estimate. \subsection{Asymptotic Validity} The following theorem establishes that statistics aggregated with (ref) are reproducible, in a particular asymptotic sense. The proof is closely related to the arguments of anscombe1952large and chow1965asymptotic and is given in (ref). See e.g., Theorem 3.1 of gut2009stopped for a textbook treatment. \begin{theorem} Suppose that the conditional variance $\Var(a(\mathsf{r}, D)\mid D)$ is strictly positive, almost surely, where $\mathsf{r}$ denotes a random collection drawn uniformly from $\mathcal{R}_{n,k,b}$. If the collections $\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^\prime_{\hat{g}^\prime,k}$ are independently obtained using (ref), then \begin{equation} P\bigg\{\big\vert a(\mathsf{R}_{\hat{g},k}, D) - a(\mathsf{R}^\prime_{\hat{g}^\prime,k}, D) \big\vert \geq \xi \mid D\bigg\} \overset{\text{a.s}}{\to} \beta \quad \text{as}\quad\xi\to0 . \end{equation} \end{theorem} Two aspects of (ref) are worthy of emphasis. First, no restrictions on the data generating process or the statistic of interest are imposed. For instance, the data $D$ do not necessarily need to be i.i.d.\ and the statistic $a(\mathsf{r},D)$ can be random though quantities other than just $\mathsf{r}$ and $D$. The result, thus, provides a strong assurance that (ref) can be applied widely. Second, in the statement (ref), asymptotics are taken as $\xi \to 0$, with all other quantities, including the sample size $n$, fixed. This asymptotic framework is somewhat opaque, or at least nonstandard, as the parameter $\xi$ is a choice variable. The operational interpretation is that, if $\xi$ is chosen to be sufficiently small---at a point negligible relative to, say, the conditional variance $v_{1,k}(D)$---then the nominal reproducibility error $\beta$ is accurate. In (ref), by imposing some simplifying restrictions, we give a set of non-asymptotic results that make the dependence of the performance of (ref) on the choices of $\xi$, $\beta$, and the statistic $a(\mathsf{r}, D) $ more transparent. The intuition underlying (ref) is straightforward. The result follows directly from the observation that the cross-splits $\mathsf{r}_i$ are independent and identically distributed conditional on the data $D$. Thus, for large values of $g$, the variance estimator $\hat{v}\left(\mathsf{R}_{g,k}, D\right)$ will be close to the conditional variance $v_{g,k}(D)$. As a consequence, if $\xi$ is sufficiently small, then the number of cross-splits $\hat{g}$ chosen by the procedure will be close to the “oracle” stopping time \begin{align} g^{\star} & = \underset{g \geq g_{\mathsf{init}} }{\arg\min} \left\{ \mathsf{cv}(\xi,\beta) = \frac{1}{2}\left(\frac{\xi}{z_{1-\beta/2}}\right)^{2} \geq v_{g,k}(D) \right\} \\ & = \underset{g \geq g_{\mathsf{init}} }{\arg\min} \left\{ \xi \geq z_{1-\beta/2} \sqrt{\frac{2 v_{1,k}(D)}{g} } \right\} , \end{align} where we have used the fact that $v_{g,k}(D) = g^{-1}v_{1,k}(D)$ in writing (ref). This, then, ensures that the aggregate statistic is approximately $(\xi,\beta)$-reproducible, as \begin{flalign*} & P\left\{\big\vert a(\mathsf{R}_{g^\star,k},D)-a(\mathsf{R}_{g^\star,k}^{\prime},D)\big\vert \geq \xi\mid D\right\} \\ & =P\left\{ \bigg\vert\sqrt{\frac{g^{\star}}{2 v_{1,k}(D)}}\left(\frac{1}{g^{\star}}\sum_{i=1}^{g^{\star}}a(\mathsf{r}_{i},D)-a(\mathsf{r}_{i}^{\prime},D)\right)\bigg\vert \geq z_{1-\beta/2}\mid D\right\} \overset{\text{a.s}}{\to} \beta \end{flalign*} as $\xi\to0$, where the limit follows from the central limit theorem (as $g^{\star} \to \infty$ as $\xi\to0$).\footnote{The discrepancy between $a(\mathsf{R}_{\hat{g},k}, D)$ and $a(\mathsf{R}_{g^\star,k}, D)$ can be shown to be negligible by applying an appropriate maximal inequality---an idea due to renyi1957asymptotic.} In most cases, sample-split statistics are used because they have some desirable property. For example, a sample-split estimator $a(\mathsf{r}, D)$ might be asymptotically normal or might perform well in terms of some loss function. The following result---roughly, a sequential version of Jensen's inequality---can be applied to show that, in many situations, statistics aggregated with (ref) inherit these properties. The proof is given in (ref) and follows from an application of the martingale stopping theorem. \begin{theorem} If the collection $\mathsf{R}_{\hat{g},k}$ is obtained with (ref), the cross-split $\mathsf{r}$ is a random element of the set $\mathcal{R}_{n,k,b}$, and $f(\cdot)$ is any convex function, then \begin{equation} \mathbb{E}\left[ f( a(\mathsf{R}_{\hat{g},k}, D) ) \right] \leq \mathbb{E}\left[ f( a(\mathsf{r}, D) )\right] . \end{equation} \end{theorem} Suppose that the estimator $a(\mathsf{r}, D)$ has good performance in terms of some convex loss. (ref) demonstrates that the aggregate estimator $a(\mathsf{R}_{\hat{g},k}, D)$ will also perform well (and, in fact, may perform better). For example, chetverikov2021cross demonstrate that the cross-validated lasso has nearly optimal rates of convergence in sparse regression problems. (ref) shows that all higher-order moments of risk estimates obtained with our procedure are smaller than those of risk estimates obtained from a single cross-split. Thus, we expect the same result to hold for tuning parameters chosen with risk estimates aggregated with (ref). Alternatively, in many cases, sample-split statistics are shown to be asymptotically normal by demonstrating that the variance of the discrepancy \begin{equation} \sqrt{n}\left(a(\mathsf{r}, D) - \frac{1}{n} \sum_{i=1}^b\psi(D_i, \eta)\right) \end{equation} is small, where $\psi(D_i, \eta)$ is some score, or influence, function and $\eta$ is some unknown nuisance parameter.\footnote{That is, if second term in (ref) is asymptotically normal and the variance of (ref) converges to zero, then the statistic $a(\mathsf{r}, D)$ is also asymptotically normal by Chebychev's inequality and the continuous mapping theorem.} This is how chernozhukov2018double demonstrate that DML estimates of average treatment effects are asymptotically normal, for example. (ref) demonstrates that the same argument will apply to statistics aggregated with (ref). In some cases, substantial residual randomness in a sample-split statistic might cause concern for the validity of associated inferences. For example, for DML estimators, the variance of the term (ref) is bounded from below by the expectation of its variance conditional on the data. So, if the conditional variance is large, then the unconditional variance of (ref) is likely not small. On the other hand, the conditional variance can be reduced through the application of our procedure. In other cases, substantial residual randomness may not be a reason for concern, because inferences are based on estimating a nuisance parameter in one split of a data set, and then conducting inference in the other split by conditioning on this estimate. This is the case for the applications considered in beaman2023selection and haushofer2022targeting. There is a large statistical literature on the use and interpretation of average $p$-values. In general, a level $\alpha$ test can be constructed by comparing an average $p$-value to $2\cdot\alpha$ ruger1978maximale,vovk2020combining,diciccio2020exact.\footnote{A variety of methodological papers use this observation to construct sample-split hypothesis tests that control the Type I error rate under very general conditions. In (ref), we show that the sample-split hypothesis tests considered in e.g., diciccio2020exact, meinshausen2009p, and wasserman2020universal continue to be valid when they are constructed sequentially with (ref).} Using a result analogous to (ref), for a pre-determined number of sample-splits $g$, chernozhukov2018generic show that, under some high-level conditions, valid $p$-values aggregated by averaging over sample-splits are also valid $p$-values. Under the same conditions, (ref) can be applied to give an analogous result for $p$-values aggregated with (ref). \subsection{Reproducible Aggregation in Practice} Equipped with (ref), we return to the examples considered in (ref). We show that, in each case, an appropriate application of the procedure produces a stabilized, reproducible estimate. \subsubsection{Cross-Fitting} We begin by treating the three examples that use cross-fitting. The most important choice to make when implementing (ref) is the specification of the statistic $a(\mathsf{r}, D)$. There are many reasonable choices that one might make here. For example, one could choose to ensure that estimates or standard errors are stable up to a desired level of precision. For the sake of comparability across settings, we opt to stabilize the $p$-value \begin{equation} a(\mathsf{r_i}, D) = 1 - \Phi\left(\frac{\mathsf{est}(\mathsf{r}_i, D) }{\mathsf{se}(\mathsf{r}_i, D) }\right) , \end{equation} where the function $\Phi(\cdot)$ is the standard normal c.d.f.\ and the quantities $\mathsf{est}(\mathsf{r}_i, D)$ and $\mathsf{se}(\mathsf{r}_i, D)$ denote the cross-split estimate and standard error that corresponding to the axes of (ref). In effect, controlling the residual randomness of the average $p$-value ensures that the residual randomness in the aggregate estimate is small relative to the sampling error. (ref) summarizes the results. We choose $\xi$ equal to either 0.001 and 0.01, as these are the levels of precision relevant for the determination of statistical significance at levels $\alpha = 0.025$ and $\alpha = 0.10$, respectively. For now, we follow the approaches to cross-splitting taken by all three papers. In the applications to chakravorty2024can and haushofer2022targeting, we aggregate over $5$-fold cross splits. In the application to beaman2023selection, the split $\mathsf{r}_i$ contains a single half-sample, i.e., subset of $[n]$ of size $n/2$.\footnote{If it is desirable to aggregate the median $p$-value, (ref) will continue to apply with a small modification. In particular, the variance estimator (ref) can be replaced by an estimator constructed with the bootstrap. In this case, a result analogous to (ref) will continue to hold.} The estimates displayed in (ref) indicate that all three applications require aggregation over hundreds of cross-splits to ensure reproducibility, in the sense of (ref), at these levels of error tolerance. Moreover, we find that nominal error rate $\beta$ associated with (ref), here set to $0.05$, is very accurate. \begin{table} \begin{centering} \caption{Reproducible Aggregation in Practice: Cross-Fitting} \begin{tabular}{lcccc} \toprule Application & $\xi$ & $p$-value & Average $\hat{g}$ & Reproducibility\tabularnewline \midrule \midrule chakravorty2024can & 0.001 & 0.979 & 860.7 & 0.945 \tabularnewline haushofer2022targeting & 0.01 & 0.80 & 1360.6 & 0.947 \tabularnewline beaman2023selection & 0.01 & 0.83 & 2794.8 & 0.950 \tabularnewline \bottomrule \end{tabular} \end{centering} \justifying {Notes: (ref) summarizes the application of (ref) to three of the examples considered in (ref). The first column indicates the application. The second column specifies choices for the error tolerance $\xi$. In each case, we set the reproducibility error rate to $\beta = 0.05$ and the burn-in sample size $g_{\mathsf{init}}$ to 10. The third column gives the $p$-value produced by one-application of the procedure. The fourth column gives the average number of cross-splits $\hat{g}$ drawn in each implementation, taken across 2,000 replications of the aggregation procedure. The fifth column displays an estimate of the true reproducibility probability, i.e., the probability that two independent implementations of (ref) produce estimates that differ by less than $\xi$, computed with these replicates.} \hrulefill \end{table} \subsubsection{Cross-Validation} We now turn to our application to cross-validated risk estimation. In this setting, there is more ambiguity in how to best apply (ref). We find that the following approach works well in the data from casey2021experiment. Recall that, in this application, the sample-split statistic $a_{\lambda}(\mathsf{r}, D)$ denotes the cross-validated estimate of the mean-squared error, queried at a specified value of the regularization parameter $\lambda$. We are interested in stabilizing these estimates for each value of $\lambda$ in an increasing sequence $(\lambda_l)_{l=1}^p$, where $\lambda_p$ is the smallest value such that no covariates have non-zero coefficients when estimated using the full data.\footnote{Throughout, we use the grid $(\lambda_l)_{l=1}^p$ of values of the regularization parameter $\lambda$ chosen by default by the “glmnet” R package friedman2021package.} To do this, we apply (ref) to each component of the vector \begin{align} (b_{\lambda_1}(\mathsf{r}, D), \ldots, b_{\lambda_{p-1}}(\mathsf{r}, D)) , \quad\text{where}\quad b_{\lambda_i}(\mathsf{r},D) = ( a_{\lambda_i}(\mathsf{r}, D) - a_{\lambda_p}(\mathsf{r}, D) ) . \end{align} The rationale for this is that the relative, rather than absolute, values of the risk estimates are what are relevant for determining where the estimates takes their minimum value. We find that, in practice, consideration of the differences (ref) can substantially reduce the amount of aggregation needed for stabilization. We implement (ref), independently, for each component of the vector (ref). In other words, we ensure that each component of the risk estimate is marginally reproducible.\footnote{In principle, it is straightforward to ensure that the components of the vector (ref) are simultaneously reproducible. For example, simultaneous reproducibility at level $\beta$ can be obtained by ensuring that each component is reproducible at level $\beta / (p-1)$, mimicking the standard Bonferonni adjustment. More sophisticated schemes, based on the bootstrap, say, are also applicable. In practice, ensuring simultaneous reproducibility can substantially increase the required computation and tends not to lead to materially different risk estimates.} We let $\xi_i$ denote the error tolerance used for the $i$th component of the vector (ref). We use a simple approach for determining suitable values for these tolerances. Roughly speaking, we compute the statistic (ref) for each element of a small, initial sample of cross-splits (e.g., $g=20$). Using these estimates, we find that the level of precision needed to distinguish between the close-to-optimal values of $\lambda$ is approximately $\xi_i = 10^{-4}$. The error tolerances specified at values of $\lambda$ that can be determined to be sub-optimal are set to larger values. Further details are given in (ref). (ref) gives measurements of the performance of (ref) for aggregation of the statistic (ref), implemented in the data from casey2021experiment. Panel A displays quantiles of the reproducibly aggregated statistic (ref) across replications of the procedure. The $y$-axis has been truncated to focus attention on the close-to-optimal values of the regularization parameter. There is little residual randomness. Panel B displays quantiles of the number of $10$-fold cross-splits $\hat{g}$ chosen by the procedure at each value of the regularization parameter. Precise estimates at the close-to-optimal values require aggregation over thousands of cross-splits. Observe that less aggregation is used at small, sub-optimal values of $\lambda$, as we have used larger error tolerances for these values. We find that, if cross-validated model selection is implemented twice, and, in each case, aggregated at this level of error tolerance, then the same value of $\lambda$ is selected with probability 0.84 and the same 14 covariates are selected with probability 0.999. In (ref), we show that the nominal reproducability error, again $\beta = 0.05$, is very accurate over the full range of the regularization parameter. \begin{figure}[t] \begin{centering} \caption{Reproducible Aggregation in Practice: Cross-Validation} \begin{tabular}{c} \textit{Panel A: Mean-Squared Error Differences}\tabularnewline \tabularnewline \textit{Panel B: Computation}\tabularnewline \tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) displays the performance of (ref) in the application to cross-validated Lasso, implemented in data from casey2021experiment. (ref) is applied independently to each component of the vector (ref). The error tolerances $\xi_i$ are specified in (ref). We set the nominal reproducibility error to $\beta = 0.05$. Panel A displays quantiles of the statistic (ref) at each value of a grid of values of the regularization parameter $\lambda$. The $y$-axes is truncated to focus attention on large values of $\lambda$. We give an un-truncated version in (ref). Panel B displays quantiles of the number of $10$-fold cross-splits $\hat{g}$ chosen by the procedure at each value of the regularization parameter.} \hrulefill \end{figure} \section{Theoretical Analysis} The asymptotic results given in (ref) are quite general. In particular, (ref) holds in the absence of any restrictions on the data generating process or statistic under consideration. It is worth asking, however, whether these results confer a clear statistical understanding. At least two issues arise. First, (ref) relies entirely on the fact that, in (ref), successive cross-splits are sampled independently. That is, we have said nothing, yet, about the role of cross-splitting. Second, the interpretation of the asymptotic approximation with $\xi\to0$ is somewhat opaque. What would be a reasonable value of $\xi$ to choose to ensure that (ref) is accurate? And how do these choices impact the amount of computation required to implement the procedure? In this section, we give a non-asymptotic description of the performance of (ref). This is accomplished by placing some simplifying restrictions on the statistic of interest. The more specialized analysis that follows is aimed at providing qualitative and quantitative intuition for how the computational cost and reproducibility error of statistics aggregated with (ref) depend on the choices of $k$ and $\xi$. Proofs for results stated in this section are given in (ref). \subsection{Symmetry, Linearity, and Stability} We impose a set of simplifying restrictions. Recall that, in general, we are considering the aggregation of cross-split statistics of the form \begin{equation} a(\mathsf{r}, D) = \mathcal{A}(\{T(\mathsf{s}_{j},D)\}_{j=1}^k) , \quad\text{where}\quad T\left(\mathsf{s}, D \right)= \Psi \left(D_\mathsf{s}, \hat{\eta} \left(D_{\tilde{\mathsf{s}}}\right)\right) \end{equation} is a sample-split statistic and the function $\hat{\eta}(\cdot)$ is an estimator of an unknown nuisance parameter $\eta$. In the main text, we restrict attention to the case that $n = k\cdot b$, i.e., where each element $\mathsf{r}$ in $\mathcal{R}_{n,k,b}$ is a complete partition of $[n]$ into $k$ sets of size $b$. Each of the results given here will follow directly from more general results stated in (ref), where this restriction is not imposed. First, we assume that the sample-split statistic under consideration is symmetric and deterministic in each part of a split sample. \begin{assumption}[Symmetry and Determinism] For all sets $\mathsf{s}$ in $\mathcal{S}_{n,b}$ and data $D$, the statistic $T\left(\mathsf{s},D\right)$ is deterministic and invariant to permutations of the data with indices in $\mathsf{s}$ and $\tilde{\mathsf{s}}$, respectively. \end{assumption} The intention of (ref) is to restrict the residual randomness under consideration to the randomness introduced by sample-splitting. This holds in cases where $\hat{\eta}(\cdot)$ is deterministic, e.g., when $\hat{\eta}(\cdot)$ is a coefficient vector determined by a regularized regression or in the applications to hypothesis testing considered by diciccio2020exact or wasserman2020universal. (ref) rules out procedures where the estimator $\hat{\eta}(\cdot)$ is random conditional on the data. This excludes settings where, e.g., $\hat{\eta}(\cdot)$ is estimated with stochastic gradient descent, bagging or subsampling, or is itself constructed with data splitting. Second, we assume that the aggregation function $\mathcal{A}(\cdot)$ is an average and that the statistic $T\left(\mathsf{s}, D \right)$ is linearly separable in the first part of the split sample. \begin{assumption}[Linearity] For all cross-splits $\mathsf{r} = (\mathsf{s}_j)_{j=1}^k$ in $\mathsf{R}_{n,k,b}$, the cross-split statistic $a(\mathsf{r}, D)$ can be represented by \begin{equation} a(\mathsf{r}, D) = \frac{1}{k} \sum_{i\in\mathsf{s}} T\left(\mathsf{s}_j, D \right) . \end{equation} Moreover, for all sets $\mathsf{s}$ in $\mathcal{S}_{n,b}$, the sample-split statistic $T\left(\mathsf{s}, D \right)$ can be represented by \begin{equation} T\left(\mathsf{s},D\right) = \frac{1}{b} \sum_{i\in\mathsf{s}} \psi(D_i, \hat{\eta} \left(D_{\tilde{\mathsf{s}}}\right)) \end{equation} for some function $\psi(\cdot,\cdot)$. \end{assumption} We make these restrictions to ease exposition. (ref) is satisfied if, for example, the statistic under consideration is a cross-fit treatment effect or cross-validated mean-squared error estimate. In principle, (ref) rules out some applications of interest. In practice, so long as the linear representations (ref) and (ref) hold up to a suitable degree of approximation, the qualitative and quantitative predictions of the results that follow will continue to hold.\footnote{Extensions of our results to cases where the function $\mathcal{A}(\cdot)$ or the statistic $T\left(\mathsf{s}, D \right)$ satisfy component-wise Lipschitz or bounded differences conditions, say, are feasible, and will exhibit the same qualitative behavior.} In particular, we show below that the predictions of our results play out in the data from chakravorty2024can, where the statistic $T\left(\mathsf{s},D\right) $ is a $p$-value of the form (ref). The remaining assumptions, and our ensuing results, are expressed in terms of two objects that measure the sensitivity of the statistic under consideration to perturbations of the data and of the splits, respectively. We refer to these objects as stabilities. They are defined as follows. \begin{defn}[Sample Stability] Fix a set $\mathsf{s}\subseteq\mathcal{S}_{n,b}$ and let $i$ be an arbitrary element of $\mathsf{s}$. Let $D^\prime$ denote an independent and identical copy of the data $D$. For each $\mathsf{q}\subseteq [n]$, let $\tilde{D}^{(\mathsf{q})}$ be constructed by replacing $D_j$ with $D^\prime_j$ in $D$ for each $j$ in $\mathsf{q}$. Let $\mathsf{q}$ be a randomly selected subset of $\tilde{\mathsf{s}}$ of cardinality $q$. We refer to the quantity \begin{align} \sigma^{(r,q)} &=\mathbb{E}\left[ \big\vert \psi(D_i, \hat{\eta} \left(D_{\tilde{\mathsf{s}}}\right)) - \psi(D_i, \hat{\eta} (\tilde{D}^{(\mathsf{q})}_{\tilde{\mathsf{s}}})) \big\vert^{r}\right] \end{align} as the $(r,q)$-order sample stability. \end{defn} \begin{defn}[Split Stability] We refer to the quantity \begin{equation} \zeta^{(r)}=\mathbb{E}\left[\max_{\mathsf{s},\mathsf{s}^{\prime}\in\mathsf{S}_{n,b}}\left(T\left(\mathsf{s},D\right)-T\left(\mathsf{s}^{\prime},D\right)\right)^{r}\right]. \end{equation} as the $r$th-order split stability. \end{defn} We restrict attention to statistics whose sample stabilities decay in a suitable way with the sample size $n$. Throughout, we say $x\lesssim y$ if there exists a universal constant $C$ such that $x\leq Cy$. \begin{assumption}[Sample Stability Decay] The $(r,q)$-order sample stability satisfies the bound \begin{equation} \sigma^{(r,q)} \lesssim \left(\frac{\sqrt{q}}{n-b}\right)^r \end{equation} uniformly for each $q$ in $[b]$ and $r$ in $\{2,4\}$. \end{assumption} We call a statistic sample stable if it satisfies (ref).\footnote{The $(2,1)$-order sample stability $\sigma^{(2,1)}$ is a widely studied object in the statistical learning literature, where it is referred to as mean-square stability bousquet2002stability,kale2011cross,kumar2013near.} Many sample-split statistics of interest are sample stable. In (ref), we show that statistics satisfying (ref) are sample stable if the nuisance parameter estimator $\hat{\eta}(\cdot)$ is an empirical risk minimizer of a, potentially regularized, strictly convex loss. There is a large literature that gives analogous bounds for other standard machine learning estimators, including bagged or subsampled estimators, like random forests chen2022debiased,ritzwoller2024uniform, ensemble estimators elisseeff2005stability, and estimators computed with stochastic gradient descent hardt2016train.\footnote{chen2022debiased show that, under some regularity conditions, sample-splitting is unnecessary for the consistency and asymptotic normality of DML estimators, if a condition related to, but partially stronger than, (ref) is satisfied. We comment on the relationship between (ref) and the conditions considered in chen2022debiased in (ref).} Nevertheless, sample-stability should be viewed as a strong assumption, that is only applicable to highly regular estimators. Below, we show that the predictions that follow from the imposition of sample-stability, concerning the qualitative behavior of the residual randomness, play out in the applications to casey2021experiment and chakravorty2024can. These predications may not have the same quality in settings that use less well-behaved nuisance parameter estimators. The split stability $\zeta^{(r)}$ is a less frequently studied object. We will only require that it is finite for $r$ equal to 4 or 8, depending on the setting. This is a weak restriction that will hold, for example, if the statistic $T\left(\mathsf{s}, D \right)$ is bounded. \subsection{Computation and Concentration} (ref) entails sequentially computing the statistic $a(\mathsf{R}_{g,k}, D)$, as $g$ increases, until a stopping criteria is satisfied. How much computation should we expect to do? And how does this quantity depend on the parameters $k$ and $\xi$? Recall from (ref) that we should expect the total number of splits $\hat{m}=\hat{g} \cdot k$ used by (ref) to be close to the “oracle” quantity \begin{equation} m^\star = g^\star \cdot k \approx 2 k \cdot v_{1,k}(D)\left(\frac{z_{1-\beta/2}}{\xi}\right)^2 , \end{equation} where $v_{1,k}(D) = \Var(a(\mathsf{r}, D)\mid D)$ denotes the conditional variance of the statistic computed using a single cross-split. Two aspects of the expression (ref) are worth highlighting. First, the total number of splits depends on the error tolerance $\xi$ through the factor $\xi^{-2}$. In other words, in order to reduce the reproducibility error by a factor of $10$, e.g., to move from $\xi = 0.1$ to $\xi = 0.01$, the total number of splits must be increased by a factor of $100$. Second, the total number of splits depends on the parameter $k$ through the factor $k \cdot v_{1,k}(D)$. How should we expect this quantity to scale with $k$? We answer this question with the following result, which characterizes the rate of convergence of the statistic $a(\mathsf{R}_{g,k}, D)$ around its conditional mean, given by \begin{equation} \bar{a}\left(D\right) = \mathbb{E}\left[a(\mathsf{R}_{g,k}, D) \mid D\right] = \mathbb{E}\left[T(\mathsf{s}_{i,j},D) \mid D\right] \end{equation} under (ref). The nonstandard aspect of this result is that we account for the dependence in the summands in $a(\mathsf{R}_{g,k}, D)$ across cross-splits, i.e., the dependence induced by cross-fitting. This is accomplished by applying a coupling argument due to chatterjee2005concentration,chatterjee2007stein.\footnote{(ref) is closely related to the unconditional variance bounds given in kale2011cross and kumar2013near, who give bounds with the same dependence on $k$ for cross-validated risk estimation. Our result follows from a different method of argument. In particular, (ref) is a corollary of a more general result, presented in (ref), that gives an analogous large deviations bound. In particular, we show that, for each $\varepsilon>1$, the bound \begin{equation*} P\left\{ \vert a(\mathsf{R}_{g,k}, D)- \bar{a}(D) \vert \leq \sqrt{\frac{b-1}{n^2} \frac{1}{g} \frac{\log(\varepsilon^{-1})}{\delta}} \mid D \right\} \geq 1 - \varepsilon \end{equation*} holds with probability greater than $1-\delta$ as $D$ varies. That is, the rate of convergence suggested by (ref) holds for all higher-order moments as well. This large deviations bound is applied repeatedly to establish the result given in the following subsection.} \begin{theorem} Suppose that (ref) hold, the data $D$ are independently and identically distributed, and $n = k\cdot b$. If the 4th-order split stability $\zeta^{(4)}$ is finite, then for each $\delta > 0$, the inequality \begin{equation} v_{g,k}(D) = \mathbb{E}\left[\left(a\left(\mathsf{R}_{g,k}, D\right)-\bar{a}\left(D\right)\right)^{2} \mid D \right] \lesssim \frac{1}{\delta} \frac{b-1}{n^2} \frac{1}{g} . \end{equation} holds with probability greater than $1-\delta$ as $D$ varies. \end{theorem} The left-hand side of the inequality (ref) is random through the data $D$. (ref) says that, if $\mathcal{F}$ is the event that the inequality (ref) holds, then $P\{\mathcal{F}\} > 1-\delta$ unconditionally. This bound results from an application of Markov's inequality, at one point in the proof, to bound a complicated, data-dependent term with a term that depends on the sample stability. This strategy---bounding the conditional quantities of interest with unconditional quantities---is helpful because the resultant unconditional object, the sample stability, is tractable and has been characterized in many settings of interest. (ref) demonstrates that the variance $v_{g,k}(D)$ converges to zero at the rate \begin{equation} \frac{1}{n} \frac{b-1}{n} \frac{1}{g} \leq \frac{1}{n} \frac{1}{k} \frac{1}{g} . \end{equation} Consequently, for a fixed total number of splits $m=g\cdot k$, the rate of convergence is proportional to $(m n)^{-1}$. That is, the conditional randomness of aggregate statistics constructed with $m$ sample-splits concentrates like averages of $m$ i.i.d.\ random variables, despite the dependence across cross-splits. Observe, also, that if $b = 1$, corresponding to “jackknife” or “leave-one-out” sample-splitting, then there is no residual randomness and the left-hand-side of (ref) collapses to zero. Plugging the bound (ref) into the expression (ref), we find that \begin{equation} m^\star \approx \frac{1}{n} \left(\frac{z_{1-\beta/2}}{\xi}\right)^2 , \end{equation} In other words, (ref) implies that---in contrast to the error tolerance $\xi$---the number of cross-folds $k$ does not affect the required computation. On the other hand, all else equal, less aggregation is needed to remove residual randomness from settings with larger sample sizes $n$. (ref) plays out empirically. (ref) displays measurements of the conditional variance $v_{1,k}(D)$ in dark blue, as $k$ varies, for our applications to casey2021experiment and chakravorty2024can. For the application to casey2021experiment, we use the cross-validated estimate of the mean-squared error at the value of $\lambda$ that minimizes the curves displayed in Panel A of (ref). For the application to chakravorty2024can, we use the $p$-value (ref), associated with the cross-fit estimate of the average treatment effect. If the variance bound (ref) is accurate, then the approximation \begin{equation} v_{1,k}(D) \approx \left(\frac{k^\prime}{k}\right) v_{1,k^\prime}(D) \end{equation} should hold for each pair $k,k^\prime$. To test this, we display estimates of right-hand-side of (ref) in light green, where $k^\prime$ is the largest value of $k$ considered in each sub-figure. In each case, the approximation is remarkably accurate.\footnote{It is worth emphasizing that neither (ref) nor (ref) holds for the application to chakravorty2024can. That is, in this case, nuisance parameters are estimated with random forests, which are random conditional on the data, and the $p$-value (ref) does not admit an exact representation of the form (ref). Rather, in this setting, both of these assumptions are good approximations (so long as the number of trees used to construct the random forests is sufficiently large), and so the predictions of (ref) hold.} \begin{figure}[t] \begin{centering} \caption{Concentration} \begin{tabular}{c} \textit{Panel A: casey2021experiment}\tabularnewline \tabularnewline \textit{Panel B: chakravorty2024can}\tabularnewline \tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) illustrates the concentration of the residual randomness of various cross-split statistics with the number of cross-splits $k$, using data from casey2021experiment and chakravorty2024can. For the application to casey2021experiment, we use the cross-validated estimate of the mean-squared error at the value of $\lambda$ that minimizes the curves displayed in Panel A of (ref). For the application to chakravorty2024can, we use the $p$-value of the form (ref), associated with the cross-fit estimate of the average treatment effect. The $x$-axes give $1/k$. The $y$-axes give measurements of the conditional variance of each statistic. Both the axes are displayed on a logarithmic scale, base 10. The theoretical prediction (ref), based on (ref), is given in light green.} \hrulefill \end{figure} \subsection{Reproducibility} (ref) confers a guarantee. In particular, statistics aggregated with (ref) should be interpreted as being reproducible, up to an error tolerance $\xi$, with probability greater than $\beta$. How accurate is this guarantee? In particular, how does the accuracy of the nominal reproducibility error $\beta$ depend on the choice parameters $\xi$ and $k$? These questions are answered by the following Berry-Esseen type bound on the accuracy of the nominal reproducibility of (ref). \begin{theorem} Suppose that the collections $\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^\prime_{\hat{g}^\prime,k}$ are independently obtained using (ref). If (ref) hold, the conditional variance $v_{1,k}(D) = \Var(a(\mathsf{r}, D)\mid D)$ is strictly positive, almost surely, the data $D$ are independent and identically distributed, and the eighth-order split stability $\zeta^{(8)}$ is finite, then for all sufficiently small $\xi$, the inequality \begin{align} & \bigg\vert P\bigg\{\big\vert a(\mathsf{R}_{\hat{g},k}, D) - a(\mathsf{R}^\prime_{\hat{g}^\prime,k}, D) \big\vert \geq \xi \mid D\bigg\} - \beta \bigg\vert \nonumber \\ &\quad\quad\quad\quad \lesssim \frac{1}{\delta^{3/4}} \frac{1}{k n} \left(\frac{1}{v_{1,k}(D)}\right)^{5/4} \left( \frac{\xi}{z_{1-\beta/2} } \right)^{1/2} , \end{align} holds with probability greater than $1-\delta$ as $D$ varies, where in writing (ref), we have omitted a multiplicative term that converges to zero logarithmically as $\xi$ decreases to zero. \end{theorem} (ref) describes the settings under which (ref) is accurate. To unpack this result, observe that the variance bound (ref) gives \begin{equation} \frac{1}{k n} \left(\frac{1}{v_{1,k}(D)}\right)^{5/4} \left( \frac{\xi}{z_{1-\beta/2} } \right)^{1/2} \gtrsim \frac{k^{1/4}}{n} \left( \frac{\xi}{z_{1-\beta/2} } \right)^{1/2} . \end{equation} This suggests that the performance of (ref) improves as $k$ and $\xi$ decrease and as $n$ increases.\footnote{The dependence of the bound (ref) on $\xi$ is sharp, at least up to the logarithmic factor. This follows from general results concerning randomly stopped sums given in landers1976exact, landers1988sharp.} These predictions hold empirically. Panel A of (ref) displays estimates of the reproducibility error in the application to casey2021experiment as $k$ and $\xi$ vary. An analogous figure for the application to chakravorty2024can is displayed in (ref). As predicted, over most of the range of $\xi$, the reproducibility error is increasing as $k$ increases.\footnote{Note, however, that for large values of $\xi$, the reproducibility error decreases as $k$ increases. This is due to early stopping. That is, if the variance $v_{1,k}(D)$ is large relative to $\xi$, then there is an increased chance that $\hat{g}$ stops immediately after the burn-in period, i.e., $\hat{g}=g_{\mathsf{init}}$.} \begin{figure}[t] \begin{centering} \caption{Performance in Application to casey2021experiment} \begin{tabular}{cc} \multicolumn{2}{c}{\textit{Panel A: Reproducibility Error}}\tabularnewline \multicolumn{2}{c}\tabularnewline \multicolumn{2}{c}{\textit{Panel B: Average Replications, $\hat{g}$}}\tabularnewline \multicolumn{2}{c}\tabularnewline \multicolumn{2}{c}{\textit{Panel C: Discrepancy from Oracle Stopping Time, $\hat{g}/g^{\star} - 1$}}\tabularnewline \multicolumn{2}{c}\tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) displays measurements of the performance of (ref) on the data from casey2021experiment. Panel A displays measurements of the reproducibility error, $P\{\vert a(\mathsf{R}_{\hat{g},k}, D) - a(\mathsf{R}^\prime_{\hat{g}^\prime,k}, D) \vert \geq \xi \mid D\}$, as $\xi$ and $k$ vary. A solid horizontal line is displayed at the nominal error rate $\beta = 0.05$. Panel B displays measurements of the average number of replications $\hat{g}$ s $\xi$ and $k$ vary. The $y$-axis is displayed with a log scale, base 10. Solid horizontal lines are placed at each exponential factor of 10. Panel C displays measurements of the 5th and 95th quantiles of the discrepancy $\hat{g}/g^{\star} - 1$ as $k$ and $\xi$ vary. Further details on the construction of this figure are given in (ref).} \end{figure} Likewise, Panel B displays estimates of the average number of cross-splits $\hat{g}$ used in (ref), at each value of $\xi$ and $k$. Again, as predicted, the total number of splits used by the procedure stays constant as $k$ varies. To see this, observe that the average value of $\hat{g}$ is roughly 10 times smaller for $k = 20$ than for $k=2$, as predicted by the approximation (ref). Comparing Panels A and B, observe that, once the average value of $\hat{g}$ is larger than approximately 500, the reproducibility error is close to the nominal error rate $\beta$. If the total number of splits required to ensure reproducibility at a given error tolerance $\xi$ is constant as $k$ varies, then why does the performance of (ref) decrease with $k$? As part of the proof of (ref), we show that, for each $\varepsilon > 0$, the bound \begin{equation} P\left\{ \bigg\vert \frac{\hat{g}}{g^\star} - 1\bigg\vert \lesssim \frac{1}{n} \frac{1}{k} \left(\frac{1}{v_{1,k}(D)}\right)^{3/2} \frac{\xi}{z_{1-\beta/2}} \sqrt{\frac{\log(\varepsilon^{-1})}{\delta}} \mid D \right\} \geq 1 - \varepsilon \end{equation} holds with probability greater than $1-\delta$, as $D$ varies. Again the variance bound (ref) gives \begin{equation} \frac{1}{n} \frac{1}{k} \left(\frac{1}{v_{1,k}(D)}\right)^{3/2} \frac{\xi}{z_{1-\beta/2}} \gtrsim \frac{k^{1/2}}{n} \frac{\xi}{z_{1-\beta/2}} . \end{equation} Written differently, as $k$ increases the discrepancy between the realized and oracle number of cross-splits, $\vert\hat{g}/g^{\star} - 1\vert$, increases, reducing the accuracy of the nominal error rate. Roughly speaking, this happens because the quality of the estimator $\hat{v}_{g,k}(D)$ for the conditional variance $v_{g,k}(D)$ depends only on the number of cross-splits $g$.\footnote{Similarly, in (ref), we show that the quality of a normal approximation to $a(\mathsf{R}_{g,k},D)$ depends only on $g$, although this discrepancy is not the leading term in the reproducibility error. The close approximation exhibited in (ref) suggests that is may be reasonable to estimate $v_{1,k}(D)$ by taking the sample variance both across and within cross-splits, i.e., computing the sample variance across all $m$ sample-splits. Although this worth further consideration, asymptotic validity, i.e., (ref), would not hold at the same level of generality.} That is, for a given number of sample-splits $\hat{m}=k\cdot \hat{g} $, (ref) is most accurate when $\hat{g}$ is large, as the variance estimate $\hat{v}_{\hat{g},k}(D)$ is more precise. As before, these predictions play out in practice. Panel C of (ref) displays estimates of the 5th and 95th quantiles of the distribution of the discrepancy $\hat{g}/g^{\star} - 1$ as $k$ and $\xi$ vary. Over most of the range of $\xi$ the discrepancy between $\hat{g}$ and $g^\star$ is increasing in $k$. At small values of $k$ and large values of $\xi$, there is an increased chance that $\hat{g}$ stops immediately after the burn-in period, i.e., $\hat{g}=g_{\mathsf{init}}$. It is worth pausing to note that these results do not support eschewing cross-splitting altogether, i.e., setting $k$ equal to one and aggregating over independent splits, as we have restricted attention to the case that $n=k\cdot b$. In (ref), we give analogous results that relax this assumption and show that cross-splitting, i.e., setting $n=k\cdot b$, reduces residual randomness at a faster rate than independent splitting, i.e., setting $k$ equal to one. In other words, all else equal, (ref) performs best, in the sense that the nominal reproducibility error is most accurate, when $k$ is equal to 2 and $b$ is equal to $b/2$. To summarize, the computation needed to achieve a desired bound on residual randomness is highly sensitive to the error tolerance $\xi$, but is not affected by the number of cross-folds $k$. On the other hand, the accuracy of the nominal error rate of (ref) decreases as the number of cross-folds $k$ increases. If the error bound $\xi$ is chosen to be suitably small, such that the realized number of cross-splits $\hat{g}$ is greater than roughly 500, then the nominal reproducibility probability tends to be quite accurate, and insensitive to changes in the number of cross-splits. \section{Recommendations for Practice} Sample-splitting is a helpful tool for simplifying many widely encountered problems in applied econometrics. This simplification comes at the cost of the introduction of residual randomness. We have shown, in several applications, that this residual randomness is large enough to substantively affect results. To address this, we have proposed a simple procedure, summarized in (ref), for removing the auxiliary randomness from sample-split statistics. The procedure takes as input a bound and an error rate. We have shown that, if the procedure were run twice, the chance that the results differ by more than the bound is well-approximated by the error rate. We conclude, in this section, by detailing several recommendations for how to best implement (ref) in practice. The most important choice to make is the specification of the sample-split statistic of interest. For example, suppose that we are interested in estimating an average treatment effect, or other causal contrast, using an estimator based on sample-splitting. Denote this quantity by $\mathsf{est}(\mathsf{r}, D)$ and let $\mathsf{se}(\mathsf{r}, D)$ denote an associated, potentially sample-split, standard error estimate. There are several reasonable choices that one might make. In this case, we recommend applying (ref) to sequentially aggregate the $p$-value \begin{equation} a(\mathsf{r}, D) = 1 - \Phi\left(\frac{\mathsf{est}(\mathsf{r}, D)}{\mathsf{se}(\mathsf{r}, D)}\right) \end{equation} associated with a test that the contrast of interest is greater than zero, where $\Phi(\cdot)$ denotes the standard normal c.d.f. This approach has the benefit of benchmarking residual randomness with an estimate of the sampling error. In some settings, it may be of interest to report the scope of residual randomness for alternative quantities associated with statistics that have aggregated with (ref). For example, suppose that we have applied (ref) to stabilize the $p$-value (ref), and obtain $a(\mathsf{R}_{\hat{g}, k}, D)$. We may also wish to report the associated estimate $\mathsf{est}(\mathsf{R}_{\hat{g}, k}, D)$. In this case, we recommend reporting and interpreting the standard error \begin{equation} \sqrt{ \frac{1}{\hat{g}} \frac{1}{\hat{g}-1} \sum^{\hat{g}}_{i=1} (\mathsf{est}(\mathsf{r}_i, D) - \mathsf{est}(\mathsf{R}_{\hat{g}, k}, D))^2 } \end{equation} in the usual way. That is, the standard error (ref) estimates the residual randomness of the estimate $\mathsf{est}(\mathsf{R}_{\hat{g}, k}, D)$, conditional on the data. If this is too large, the error tolerance $\xi$ should be decreased. The optimal choice of statistic in settings where cross-validation is used for model selection and risk estimation is less immediate. In this case, we recommend applying (ref) to each element of the statistic (ref), i.e., the relative values of the risk estimates. The second most important choice to make is the choice of the error tolerance $\xi$. In some cases, an appropriate choice is clear. For example, it may be desirable to ensure that a $p$-value is reproducible at the level relevant for the determination of statistical significance at level 0.10 or 0.025, e.g., setting $\xi = 0.01$ or $\xi = 0.001$. In other cases, like in applications to cross-validated model selection, this choice is less clear. We propose a procedure for making this choice in that setting in (ref). The error tolerance $\xi$ should be set to a value that is sufficiently small so that the nominal reproducibility error is accurate. On the other hand, $\xi$ should not be set so small that the computation associated with implementing (ref) becomes infeasible. To ensure that the chosen value of $\xi$ satisfies these constraints, we recommend computing the statistic of interest for each of a small, initial sample of cross-splits (e.g., $g$ equal to 20 or 30). Let $\hat{v}_{1,k}(D)$ denote an estimate of the conditional variance $v_{1,k}(D)$ computed using this sample. An estimate of the total number of cross-splits needed to implement (ref) at a specified level of $\xi$ can then be obtained by \begin{equation} g(\xi) = 2 \hat{v}_{1,k}(D)\left(\frac{z_{1-\beta/2}}{\xi}\right)^2 . \end{equation} Motivated by the numerical results reported in (ref), we recommend choosing a value of $\xi$ such that this estimate is greater than 500. In some cases, the estimate (ref) may be too large at practically relevant values of $\xi$ to be computationally feasible. Unfortunately, as we have shown in (ref), changing the number of cross-folds $k$ will not address this issue, and it may be best to consider alternative choices of nuisance parameter estimators. In this situation, one might consider using “leave-one-out” or “jackknife” sample-splitting, i.e., setting $k = n$, which exhibits no residual randomness. This is not an omnibus fix, however. For example, proofs of the asymptotic normality of DML estimates of average treatment effects require that $k$ is small relative to $n$ chernozhukov2018double. Similarly, leave-one-out cross-validation is not necessarily optimal for model selection (see e.g., shao1993linear for early discussion of this point). chetverikov2024tuning gives a review of various, alternative, approaches to tuning parameter selection. If the number of cross-folds $k$ is set too high, the realized number of cross-splits $\hat{g}$ may be too small and (ref) may perform poorly. But, so long as the user ensures that the error tolerance $\xi$ is sufficiently small such that the number of cross-splits $\hat{g}$ tends to be large (e.g., greater than 500, say), small changes in the number of folds $k$ should not have adverse effects. Thus, this choice should be made on substantive grounds, that will depend on the application. In practice, the conventional choices of $k$ equal to 2, 5, or 10 should work well in most settings. We have had relatively little to say about the choice of the nominal reproducibility error $\beta$. We have found that setting $\beta = 0.05$ is suitable for most applications. Finally, although it has not played a central role in our theoretical analysis, it is important to choose the burn-in period $g_{\mathsf{init}}$ to be suitably large to ensure that the variance estimate $\hat{v}_{g_{\mathsf{init},k}}(D)$ is reasonably accurate, otherwise (ref) will tend to stop too early. We recommend setting $g_{\mathsf{init}}$ to be equal to at least 10.

spacing{1.2}

appendix\begin{center} {\it Supplemental Appendix to:} \begin{spacing}{1} {Reproducible Aggregation of Sample-Split Statistics\daggerfootnote{Date: \today}} \end{spacing} \end{center} \begin{center} \begin{tabular}[t]{c@{\extracolsep{4em}}c} {David M. Ritzwoller} & {Joseph P. Romano}\\ {Stanford University} & {Stanford University} \\ \end{tabular} \end{center} \begin{spacing}{1.13} \DoToC \end{spacing} \thispagestyle{empty} \setcounter{page}{0} \setcounter{figure}{0} \begin{spacing}{1.3}