Training and Testing with Multiple Splits: A Central Limit Theorem for Split-Sample Estimators

Training and Testing with Multiple Splits: A Central Limit Theorem for Split-Sample Estimators

Introduction

As predictive algorithms transform empirical economics, policy, and industry, it is now routine to use the same dataset to train and evaluate a new model. For example, when training a machine learning algorithm to predict treatment effects, create a targeted policy rule, or automate consumer credit scoring, it is essential to evaluate the quality of the predictions and assess whether implementation would generate disparate impact across demographic groups. However, standard methods for training and evaluating typically waste part of the data by splitting the sample into training and testing sets. I develop a new inference approach that uses the entire sample for both tasks by combining multiple splits, improving statistical power and reproducibility. I provide valid confidence intervals under weak conditions for model properties such as accuracy and fairness calculated using cross-fitting, repeated sample-splitting or repeated cross-fitting.

Specifically, I study a setting in which an analyst (a researcher, policymaker, or industry practitioner) wishes to use the same dataset to both:

For example, consider a government using machine learning (ML) to target recipients of a poverty alleviation program. Step (i) consists of training a model to predict, for example, families at higher risk of falling below the poverty line, while step (ii) consists of constructing a confidence interval for the out-of-sample mean squared error or rate of correct classifications.

Using the same observations for both steps (i) and (ii) creates a form of statistical dependence that makes inference challenging. For example, standard central limit theorems (CLTs) assume independence, which is violated in this setting. This difficulty is often overcome by randomly splitting the sample into two, one part to train the model (training sample), and the other to evaluate its properties (evaluation sample). Since each task is conducted with separate data, such statistical dependence is not generated, and one can use standard approaches to inference. This procedure, however, has three drawbacks: it uses only part of the data for training the model, only part of the data for evaluating its properties, and different random splits can lead to widely different estimates and potentially affect statistical significance.

My main contribution is an inference approach that averages estimates across multiple sample splits, improving upon a standard 50/50 split by using more data for training, using twice as much data for evaluation, and improving reproducibility. In empirical applications and Monte Carlo experiments, I show that these improvements often reveal statistical significance that would otherwise be missed. The main challenge of using multiple splits is a new form of statistical dependence due to reusing observations in both training and evaluation roles across different splits. I address this challenge by proving a new CLT for a large class of split-sample estimators under weak conditions. My CLT builds on previous literature in two key dimensions: (i) it applies to a large class of estimators and split-sample procedures, and (ii) it imposes no restrictions on model complexity or rates of convergence and stability, only requiring that the estimated model converges to an arbitrary limit at any rate. This generality is crucial for accommodating popular ML algorithms, such as random forests or neural networks. Moreover, I characterize when confidence intervals based on the normal approximation are valid, showing they may fail in important cases such as comparing the performance of two models or learning features of heterogeneous treatment effects. I develop a new inference approach for such cases that explicitly accounts for the dependence across splits, leveraging my CLT. Finally, I develop a reproducibility measure for p-values from split-sample procedures, quantifying whether the number of repetitions is sufficiently large to ensure reproducible inference.

To illustrate the technical challenges and empirical implications of my results, consider the problem of predicting poverty as a simple running example, which is one of my applications. Accurate out-of-sample poverty prediction is central to Development Economics for understanding poverty dynamics and designing targeted interventions. I focus on assessing predictive accuracy as the natural starting point, though my framework applies more broadly. Consider a sample $D = (Y_i, X_i)_{i=1}^n$ of $n$ households, where $X_i$ are covariates and $Y_i$ is a binary indicator equal to one if household $i$ is below the poverty line 13 years after the covariates were measured. The goal is to use the sample to (i) train a model ${\hat{\eta}}(x)$ to predict poverty by estimating $P(Y = 1 | X = x)$, for example using a machine learning algorithm, and (ii) evaluate its accuracy, for example by estimating and calculating a confidence interval (CI) for the out-of-sample mean squared error (MSE) $$\theta_{\hat{\eta}} = \E{\left( Y_{new} - {\hat{\eta}}(X_{new}) \right)^2 | D} = \int \left( y - {\hat{\eta}}(x) \right)^2 dP(y,x),$$ where $(Y_{new},X_{new})$ is an out-of-sample observation drawn from the same distribution as the sample. Note that $\theta_{\hat{\eta}}$ is data-dependent, and is thus different from targeting a parameter $\theta_{\eta_{0}}$ for some fixed ${\eta_{0}}$. In the policy prediction example, the researcher is not interested in the out-of-sample accuracy of an ideal but unknown model ${\eta_{0}}$. Instead, they are interested in the accuracy of the actually estimated model ${\hat{\eta}}$.

In this context, confidence intervals are often constructed using sample-splitting. If the entire sample is used for both tasks, standard CLTs do not apply to the average $$\frac{1}{n} \sum_{i=1}^{n} (Y_i - {\hat{\eta}}(X_i))^2,$$ since the summands are not independent. For example, $Y_1 - {\hat{\eta}}(X_1)$ and $Y_2 - {\hat{\eta}}(X_2)$ are dependent since ${\hat{\eta}}$ is estimated with both $(Y_1, X_1)$ and $(Y_2, X_2)$. A standard approach to handle this dependence is to impose complexity restrictions on how ${\hat{\eta}}$ is estimated, such as Donsker conditions. These restrictions hold for simple procedures like ordinary least squares, but fail for complex machine learning algorithms frequently used in applied problems chernozhukov2018double. Sample-splitting avoids this dependence without strong assumptions: randomly split $\{1, \dots, n\}$ into sets $\mathsf{s}_1$ and $\mathsf{s}_2$ of size for example $n/2$, use data in $\mathsf{s}_1$ to estimate ${\hat{\eta}}_{1}$, and data in $\mathsf{s}_2$ to calculate the average

Since the summands in $\hat{\theta}_{{\hat{\eta}}_{1}}$ are independent conditional on $\mathsf{s}_1$, standard CLTs apply, and the normal approximation gives a valid CI for $\theta_{{\hat{\eta}}_1}$. However, this procedure uses only half of the data for each task, and different random splits can lead to widely different estimates and potentially different conclusions about statistical significance.

Using multiple splits can improve upon these drawbacks but introduces a new challenge. Consider, for example, two-fold cross-fitting, where the roles of samples $\mathsf{s}_1$ and $\mathsf{s}_2$ are reversed and the final estimator averages the split-specific estimates. That is, estimate ${\hat{\eta}}_{1}$ using $\mathsf{s}_1$ and ${\hat{\eta}}_{2}$ using $\mathsf{s}_2$, then calculate the final estimator

where ${\hat{\eta}} = ({\hat{\eta}}_1, {\hat{\eta}}_2)$. The estimand in this case is the MSE of an average model, as discussed in the next paragraph. While this estimator averages over all $n$ observations, standard CLTs do not apply due to a different form of statistical dependence: the first sum is not independent of the second since both use the entire dataset. My first main contribution is a central limit theorem for a large class of estimators that includes $\hat{\theta}_{\hat{\eta}}$, which I use to construct valid CIs. In addition to using the entire sample for evaluation in (ref), which reduces the variance of the asymptotic distribution compared to that of (ref), more data can be used for training by increasing the number of folds. With 3 folds, for example, three models are trained, each using two-thirds of the data, with the remaining third used to evaluate the MSE. Finally, reproducibility is improved by repeating the splitting process multiple times and averaging the estimators over repetitions.

I show that $\sqrt{n} (\hat{\theta}_{\hat{\eta}} - \theta_{\hat{\eta}})$ is asymptotically normal under weak conditions, targeting its out-of-sample expectation $$\theta_{\hat{\eta}} = \E{\frac{1}{2} \left( Y_{new} - {\hat{\eta}}_1(X_{new}) \right)^2 + \frac{1}{2} \left( Y_{new} - {\hat{\eta}}_2(X_{new}) \right)^2 | D}.$$ In the example above, $\theta_{\hat{\eta}}$ is mathematically equivalent to the MSE of the average model, that is, $\theta_{\hat{\eta}} = \theta_{\bar{\eta}}$, where $$\theta_{\bar{\eta}} = \E{\left( Y_{new} - {\bar{\eta}}(X_{new}) \right)^2 | D}, \qquad {\bar{\eta}}(x) = \frac{1}{2} {\hat{\eta}}_1(x) + \frac{1}{2} {\hat{\eta}}_2(x).$$ This happens anytime the outcome is binary, and holds for the MSE, mean absolute deviation, among others, including when averaging over multiple folds and repetitions. In the poverty prediction example, this means that a researcher or policymaker can use model ${\bar{\eta}}$ for out-of-sample predictions, which will have MSE $\theta_{\hat{\eta}}$. For continuous outcomes, the researcher has two options. The first is to use a model ${\tilde{\eta}}(x)$ that predicts a value in $({\hat{\eta}}_1(x), {\hat{\eta}}_2(x))$ at random. This model has an out-of-sample MSE equal to $\theta_{\hat{\eta}}$. Alternatively, one could still use ${\bar{\eta}}$, which has the guarantee to perform better or equal than ${\tilde{\eta}}$ in terms of out-of-sample accuracy due to a risk-contraction property (for details, see (ref)).

I make three main contributions. First, I prove a new central limit theorem for a large class of split-sample estimators under mild conditions. Specifically, I make no restrictions on the complexity of the models ${\hat{\eta}}$, or on their rates of convergence or algorithmic stability. For sample-average estimators, my CLT follows under a standard moments condition and assuming that ${\hat{\eta}}$ converges to an arbitrary limit, at any rate. I show that the normal approximation yields a valid CI in many applications, but may fail to do so in important cases of interest, such as comparing the performance between two models or some instances when ${\hat{\eta}}$ converges to a constant. My second contribution builds on the CLT to develop a new inference approach that covers such cases, explicitly accounting for the dependence across splits. I focus on the case of comparing the performance between two models, and discuss how the arguments apply more broadly to other cases. Finally, I develop a reproducibility measure for p-values obtained from split-sample estimators. It addresses a common concern: another researcher using the same dataset, but different splits, may reach a different conclusion about statistical significance. For a given (large) number of repetitions of sample-splitting/cross-fitting, my measure quantifies p-value reproducibility, assessing whether the number of repetitions is sufficiently large to ensure reproducible inference.

Other contributions include a central limit theorem for split-sample empirical processes, which I use to prove my main central limit theorem, and may be of independent interest. I also apply this CLT to develop a new ensemble method for learning features of heterogeneous treatment effects in randomized experiments, following the framework of chernozhukov2025generic. The ensemble method improves on previous alternatives by using the entire sample for evaluation, more data for training, and combining multiple machine learning predictors, potentially improving power and avoiding issues of multiple hypothesis testing.

I apply my inference approaches to two important problems in development and public economics: predicting poverty and learning heterogeneous treatment effects in randomized experiments. In the first application, using a panel from Ghana ghanapaneldataset and Monte Carlo experiments, I show that repeated cross-fitting outperforms previous alternative approaches in detecting predictive power for being below the poverty line 13 years ahead. In the second application, I revisit the experiment of karlan2007does on charitable giving and conduct Monte Carlo simulations, and show that my ensemble method achieves improved power for detecting heterogeneous treatment effects compared to previous alternatives.

The rest of the paper is structured as follows. (ref) summarizes related work, and (ref) establishes the setup and notation. (ref) establishes a central limit theorem for split-sample Z-estimators, and (ref) develops inference using the normal approximation and for comparing two models. I introduce my measure of reproducibility in (ref). Finally, I implement my inference approaches in two empirical applications: predicting poverty in Ghana in (ref) and heterogeneous treatment effects in charitable giving in (ref). (ref) concludes. Proofs are delayed to (ref).

Related Work

I contribute to the literature on inference using multiple splits of the data. The literature on risk estimation via cross-validation provides related results establishing asymptotic normality for sample-average estimators based on multiple splits. Like my approach, these CLTs target data-dependent parameters, but rely on different types of assumptions and focus on the particular case of sample-averages. dudoit2005asymptotics consider estimators that average over all possible splits or cross-splits of the sample, assume bounded loss function and require ${\hat{\eta}}$ to be loss consistent for a risk minimizing function, whereas I assume that ${\hat{\eta}}$ converges to an arbitrary limit. austern2020asymptotics and bayle2020cross provide CLTs under rate assumptions on the algorithmic stability of ${\hat{\eta}}$. bayle2020cross provide two CLTs using estimators based on a single repetition of cross-fitting, one relying on rate conditions for algorithmic stability, and the second requires a “conditional variance convergence” assumption that they verify using rates for loss stability. My result does not require verifying loss stability conditions, which may not be satisfied in some high-dimensional settings bates2024cross, and my result allows for any ML algorithm as long as ${\hat{\eta}}$ converges to an arbitrary limit at any rate, in the sense established in (ref). ledell2015computationally provide a CLT for the particular case of estimating the area under the curve (AUC) measure via cross-validation, and andrews2022transfer derive confidence intervals for the different but related problem of learning the transfer error of a model across domains. Moreover, I document that asymptotic normality of split-sample estimators does not immediately lead to valid inference for the important problem of comparing the performance between two models, and I construct a new inference approach for this case that explicitly incorporates the dependence across splits.

A different class of related results show asymptotic normality using cross-fitting when targeting parameters that are not data-dependent. These approaches require stronger conditions on ${\hat{\eta}}$ that may not hold in general for nonparametric models with more than a handful of covariates, such as requiring ${\hat{\eta}}$ to converge in probability at some specified rate luedtke2016statistical,belloni2017program,chernozhukov2018double,benkeser2020improved,imai2025statistical. Leveraging the data-dependent parameter of interest, my CLT ((ref)) requires no complexity restrictions, and assumes ${\hat{\eta}}$ converges in probability to any limit at any rate.

In the context of learning features of heterogeneous treatment effects in randomized trials, chernozhukov2025generic proposed taking the median of estimators, confidence intervals and p-values across splits, similarly focusing on a data-dependent parameter, without relying on complexity or rate assumptions. wager2024sequential proposed a modified, sequential approach based on luedtke2016statistical, and chernozhukov2025reply suggested taking the median over repetitions of the sequential approach. In the same framework, imai2025statistical developed inference using cross-fitting, relying on rate assumptions. My results build on this literature in four main dimensions, relying on the mild assumption that the trained models converge to any limit, at any rate. First, my estimator uses all observations on the role of evaluation sample, leading to a smaller variance of its asymptotic distribution. Second, my approach does not exhibit a tradeoff between training and evaluation sample sizes, allowing for more data to be used to train the models. Third, I provide inference for an interpretable estimand under no rate assumptions on the trained models, while chernozhukov2025generic require a rate of concentration condition for coverage of their median estimand, which requires for example the training sample to be large relative to the evaluation sample. Finally, I introduce a new ensemble method that combines predictions from multiple ML algorithms, potentially improving statistical power for detecting HTE, and avoiding issues of multiple hypothesis testing.

The literature on learning features of heterogeneous treatment effects with multiple splits is a subset of a broader literature on aggregating potentially dependent p-values ruger1978maximale,ruschendorf1982random,meng1994posterior,meinshausen2009pvalues,gasparin2025combining. These approaches similarly apply to data-dependent parameters under weak conditions, and typically target size control under the worst data generating process, thus being conservative in general. My confidence intervals are asymptotically exact and improve statistical power.

Finally, my work is complementary to ritzwoller2023reproducible. They provide a stopping algorithm for determining how many times to repeat sample-splitting to ensure reproducibility of averages over split-sample statistics, for example, the average point estimate. I take the number of repetitions as given and focus on inference, providing a measure of reproducibility for p-values calculated using multiple splits. While ritzwoller2023reproducible uses an asymptotic framework that takes the sample size fixed and assumes a small threshold for the variability of the average split-sample statistic, my framework uses a growing sample size and number of repetitions.

Setup

I consider a setup in which an analyst (a researcher, policymaker, or industry practitioner) wishes to use a dataset to both (i) train a new model and (ii) evaluate some of its properties. This is typically the case when one wants to train a new model to automate or assist decision-making, for example using a machine learning algorithm. Since these algorithms, despite their potential, may perform poorly in practice or have disparate performance across groups, it is often important to assess their accuracy and fairness. I use the term fairness as in the algorithmic fairness literature chouldechova2018frontiers,cowgill2020algorithmic,barocas2016big and provide example measures below. I state the analyst's goals, discuss the parameter of interest with examples, and introduce the split-sample procedures I study.

The first goal of the analyst is to train a model ${\hat{\eta}} \in H$ using an algorithm and a dataset $D = (W_i)_{i=1}^n$, where each $W_i \in \mathcal{W}$ is an iid draw from a distribution $P$. I use train to denote the fitting/estimation of ${\hat{\eta}}$ using $D$, training algorithm (or just algorithm) for the procedure that maps $D$ to ${\hat{\eta}}$, and estimated model (or just model) for the realized function ${\hat{\eta}}$. For example, one can use the Random Forests algorithm to train a new model ${\hat{\eta}}$. The sets $H$ and $\mathcal{W}$ are in principle unconstrained, and $H$ depends on the choice of training algorithm. Typically, the analyst will estimate ${\hat{\eta}}$ by minimizing some loss function. My setup, however, is agnostic to the choice of training algorithm, and all results hold for any algorithm as long as ${\hat{\eta}}$ converges to an arbitrary limit at any rate, in the sense defined in (ref).

The second goal of the analyst is to use $D$ to evaluate some performance property of ${\hat{\eta}}$, denoted $\theta_{\hat{\eta}}$. Specifically, the analyst wishes to construct a confidence interval $\widehat{\rm CI}_{\alpha}$ for $\theta_{\hat{\eta}}$ such that, for $\alpha \in (0,1)$,

where the probability accounts for the randomness in both ${\hat{\eta}}$ and $\widehat{\rm CI}_{\alpha}$. $\theta_{\hat{\eta}}$ can be, for example, a measure of accuracy or fairness of ${\hat{\eta}}$. The parameter of interest, $\theta_{\hat{\eta}}$, depends on the data through the estimated model ${\hat{\eta}}$. This differs from the standard semiparametric literature, where the parameter of interest takes the form of $\theta_{\eta_{0}}$ for some nuisance function ${\eta_{0}}$. In the applications I consider, the object of interest is $\theta_{\hat{\eta}}$ since the analyst/policymaker is interested in the accuracy or fairness of the specific estimated model ${\hat{\eta}}$ they can actually implement. This is different from evaluating the performance of an ideal but unknown model ${\eta_{0}}$.

I provide three examples of such parameters of interest, and then discuss related cases in the literature where the parameter of interest is data-dependent.

There are several examples in the literature where the parameter of interest takes the form of a data-dependent $\theta_{\hat{\eta}}$. This occurs anytime the hypothesis of interest is selected only after the data has been used dawid1994selection. An important case is the approach of chernozhukov2025generic to inference on features of heterogeneous effects in randomized trials, which I revisit in (ref). Other examples include evaluating the impacts of data-driven algorithms in policy applications potash2015predictive,kuzmanovic2024causal, measuring welfare gains generated from data-driven rules kitagawa2018should,ida2024dynamic, and the “inference on winners” framework of andrews2024inference.

My setup also applies to some cases where the parameter of interest is not data dependent, but is estimated using split-sample techniques. For example, in fava2025predicting I develop an approach to inference on points of the distribution of treatment effects. Although the parameter of interest, $\theta$, is not data dependent, I incorporate covariate-adjustment terms ${\hat{\eta}}$ that yield bounds $\theta_{{\hat{\eta}},L} \le \theta \le \theta_{{\hat{\eta}},U}$. Inference on $\theta$ can then be derived from the asymptotic distribution of split-sample estimators $(\hat{\theta}_{{\hat{\eta}},L},\hat{\theta}_{{\hat{\eta}},U})$, centered around the bounds $(\theta_{{\hat{\eta}},L},\theta_{{\hat{\eta}},U})$. Other examples where $\theta_{\hat{\eta}}$ is informative about a parameter $\theta$ include learning the mean outcome under an optimal treatment regime shi2020breaking, fischer2024bridging, and averages of intersection bounds ji2024model,semenova2025debiased. Another type of application is when $\theta = \theta_{\hat{\eta}}$ does not depend on ${\hat{\eta}}$, yet estimating $\theta_{\hat{\eta}}$ leads to some better properties. This is the case of adding a covariate-adjustment term for learning the average treatment effect in a randomized trial, as I discuss in Appendix (ref).

I consider four split-sample procedures for attaining the analyst's goals: 1) sample-splitting, 2) cross-fitting, 3) repeated sample-splitting, and 4) repeated cross-fitting. First, I introduce some notation. Let $\left[ n \right] = \left\{ 1, \dots, n \right\}$ and the dataset $D = ( W_i)_{i \in \left[ n \right]}$ be an iid sample of $W \sim P$. I denote the training algorithm by $\mathcal{A} : \mathcal{W}^m \to H$, a function that takes a sample of size $m$ and returns a value $\eta \in H$. The dependence on $m$ is suppressed in the notation of $\mathcal{A}$. For any subsample $\mathsf{s} \subset \left[ n \right]$, $D_{\mathsf{s}} = \{ W_i \}_{i \in \mathsf{s}}$.

Sample-splitting consists of taking a random subsample $\mathsf{s} \subseteq \left[ n \right]$ of size $b$, using its complement $\tilde{\mathsf{s}} = \left[ n \right] \setminus \mathsf{s}$ to train the model ${\hat{\eta}}_{\tilde{\mathsf{s}}} = \mathcal{A}(D_{\tilde{\mathsf{s}}})$, and calculating $\widehat{\rm CI}_{\alpha}$ from $\mathsf{s}$ for the parameter $\theta_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}$. Cross-fitting consists of partitioning $\left[ n \right]$ into $K$ roughly equal-sized subsets (folds) $(\mathsf{s}_k)_{k=1}^K$, at random. For $k=1,\dots,K$, train a model ${\hat{\eta}}_{\tilde{\mathsf{s}}_k} = \mathcal{A}(D_{\tilde{\mathsf{s}}_k})$, that is, using all observations except those in fold $k$. Each model ${\hat{\eta}}_{\tilde{\mathsf{s}}_k}$ is trained from $n (K-1) / K$ observations when $n$ is a multiple of $K$. In (ref), I discuss different ways to aggregate the $K$ models into an estimand $\theta_{\hat{\eta}}$, where ${\hat{\eta}} = ({\hat{\eta}}_{\tilde{\mathsf{s}}_k})_{k=1}^K$, and the construction of a confidence interval $\widehat{\rm CI}_{\alpha}$ for $\theta_{\hat{\eta}}$. I consider $K$ fixed as $n \to \infty$.

Repeated sample-splitting and cross-fitting consist of repeating the procedures above $M$ times. That is, for repeated sample-splitting, take $M$ independent, random subsamples of $\left[ n \right]$ of size $b$, $(\mathsf{s}_{m,1})_{m=1}^M$, and train $M$ models $({\hat{\eta}}_{\tilde{\mathsf{s}}_{m,1}})_{m=1}^M$. For repeated cross-fitting, take $M$ independent, random partitions of $\left[ n \right]$ into $K$ roughly equal-sized folds, $\mathcal{R} = (r_m)_{m=1}^M$, where each $r_m = (\mathsf{s}_{m,k})_{k=1}^K$ forms a partition of $\left[ n \right]$. For each subsample $\mathsf{s}_{m,k}$, train a model ${\hat{\eta}}_{\tilde{\mathsf{s}}_{m,k}}$ using all observations except the ones in $\mathsf{s}_{m,k}$, giving a total of $M K$ models. I discuss different ways to aggregate the multiple splits in (ref).

I give a unified notation to the four split-sample procedures described above. Let $\mathcal{R} = \left( r_m \right)_{m=1}^M$ denote a collection of $M$ random splits of the sample, where each split can be either:

I use $K=1$ to denote sample-splitting for convenience. $K=1$ means that $r_m$ consists of one subsample, of size $b < n$ chosen by the researcher. For cross-fitting, I assume folds are equal-sized if $n$ is a multiple of $K$, and have sizes $\lfloor n/K \rfloor$ and $\lceil n/K \rceil$ otherwise, and define $b = n/K$. Define $\pi = \lim_{n \to \infty} b/n$, $\pi \in (0,1)$. With this notation, $K=1$ denotes sample-splitting and $K>1$ denotes cross-fitting. I allow $M$ to grow as $n$ increases, and denote $$\bar{M} = \lim_{n \to \infty} M \in \mathbb{N} \cup \left\{ +\infty \right\}.$$ This notation unifies the four split-sample procedures described previously, as shown in (ref).

I use the term multiple splits to denote any of the three procedures that use more than one split ($\bar{M} > 1$ and/or $K>1$). In all cases, I assume that the splits are taken at random uniformly over all possible splits or cross-splits. Although the number of possible splits is finite for any given $n$, I consider that the $M$ repetitions are taken independently, with repetition. This assumption reflects common practice, as the computationally feasible number of repetitions is usually much smaller than the total number of possible splits, so that the probability of taking two identical splits is negligible.

I compare the four split-sample procedures in terms of statistical power, modeling power, and reproducibility properties in (ref).

CLT for Split-Sample Z-Estimators

I prove a central limit theorem for split-sample Z-estimators, defined as zeroes of empirical moment equations. Z-estimators are a large class of estimators which include averages, linear regressions, and most M-estimators, since the parameter value that maximizes some objective function is the same that sets its partial derivatives to zero. This CLT can be used off-the-shelf in many applications, including the poverty prediction application in (ref). First, in (ref), I define split-sample Z-estimators and Z-estimands, introduce the assumptions used, and state the CLT. Finally, in (ref), I compare the four split-sample procedures (sample-splitting, cross-fitting, repeated sample-splitting, and cross-fitting).

I provide a more accessible exposition for the particular case of sample average estimators, such as the MSE ((ref)), in (ref). I prove a new CLT for split-sample empirical processes in (ref), which I use to prove my CLT for Z-estimators and may be of independent interest.

Main Result

Since Z-estimators can be nonlinear, unlike the mean squared error ((ref)), different approaches to aggregating multiple splits lead to different estimators and estimands. I discuss three such approaches. Let $\left\lVert \cdot \right\rVert$ be the Euclidean norm, $\psi_{\theta,\eta}: \mathcal{W} \to \mathbb{R}^d$ be measurable functions for $\theta \in \Theta \subseteq \mathbb{R}^d$ and $\eta \in H$ ($H$ is defined as in (ref)), and $d \ge 1$. For $\eta \in H$, let $\Psi_{\eta}(\theta) = P \psi_{\theta, \eta}$, $\hat{\Psi}_{\mathsf{s},\eta}(\theta) = |\mathsf{s}|^{-1} \sum_{i \in \mathsf{s}} \psi_{\theta, \eta}(W_i)$, and $\dot{\Psi}_{\eta}$ be the Jacobian matrix of $\Psi_{\eta}(\theta_\eta)$. As in (ref), let $\mathcal{R}$ denote a collection of splits with $M$ repetitions and $K$ folds, and let ${\hat{\eta}} = {\hat{\eta}}_{\mathcal{R}} = \left( {\hat{\eta}}_{\tilde{\mathsf{s}}_{m,k}} \right)_{m \in \left[ M \right], k \in \left[ K \right]}$.

The first type of estimand is an average across split-specific estimands:

where $\theta_{\eta}^{(1)}$ for $\eta \in H$ is the unique solution for $\theta$ in $\Psi_{\eta}(\theta) = 0$, i.e., $$\Psi_{\eta}(\theta_{\eta}^{(1)}) = 0.$$ (ref) consists of solving the moment condition $\Psi_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}(\theta) = 0$ for each split $\tilde{\mathsf{s}}$, and averaging over the split-specific estimands. The Z-estimator for (ref) is

where $\hat{\theta}_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}^{(1)} \in \arg\min_{\theta \in \Theta} \left\lVert \hat{\Psi}_{\mathsf{s},{\hat{\eta}}_{\tilde{\mathsf{s}}}}(\theta) \right\rVert$. This approach is analogous to the DML1 estimator in chernozhukov2018double.

The second type of estimand solves the average of the moment conditions. That is, $\theta_{{\hat{\eta}}}^{(2)}$ uniquely solves

The associated Z-estimator is given by

This approach is analogous to the DML2 estimator in chernozhukov2018double.

Finally, the third type of estimand is a hybrid of the previous two approaches. It solves the moment condition at each cross-split of the sample, and averages across repetitions. That is,

where $\theta_{{\hat{\eta}}_{r}}^{(2)}$ uniquely solves $$\frac{1}{K} \sum_{\mathsf{s} \in r} \Psi_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}\left( \theta_{{\hat{\eta}}_{r}}^{(2)} \right) = 0.$$ The associated Z-estimator is given by

where

In this approach, each $\hat{\theta}_{{\hat{\eta}}_{r}}^{(3)}$ uses the whole sample both for calculating ${\hat{\eta}}_{r}$ and the average in (ref), and the final estimator $\hat{\theta}_{{\hat{\eta}}}^{(3)}$ is the average of the cross-fitting estimators across repetitions. Note that $\theta_{{\hat{\eta}}}^{(1)} = \theta_{{\hat{\eta}}}^{(3)}$ if $K=1$, $\theta_{{\hat{\eta}}}^{(2)} = \theta_{{\hat{\eta}}}^{(3)}$ if $M = 1$, and $\theta_{{\hat{\eta}}}^{(1)} = \theta_{{\hat{\eta}}}^{(2)} = \theta_{{\hat{\eta}}}^{(3)}$ if $K=M=1$. The estimators are not assumed to be unique, but I assume the estimands and the limit of the estimators to be unique.

For a concrete example, I consider below the particular case of sample-averages, as in the example of calculating the MSE for poverty prediction ((ref)).

$\hat{\theta}_{{\hat{\eta}}}^{(2)}$ can be interpreted as the value of $\theta$ that solves the moment condition for a randomized function that takes value across $\left( {\hat{\eta}}_{\tilde{\mathsf{s}}_{m,k}} \right)_{m \in \left[ M \right], k \in \left[ K \right]}$ uniformly at random. That is, (ref) is equivalent to $$\int \psi_{\hat{\theta}_{{\hat{\eta}}}^{(2)}, {\hat{\eta}}_{\xi}}(w) dP(w, \xi) = 0,$$ where $dP(w, \xi) = dP(w) dP(\xi)$ and $\xi$ takes value in $\left( \tilde{\mathsf{s}}_{m,k} \right)_{m \in \left[ M \right], k \in \left[ K \right]}$ uniformly at random. If, for example, each ${\hat{\eta}}_{\tilde{\mathsf{s}}}$ is a probabilistic classifier as in (ref), $\hat{\theta}_{{\hat{\eta}}}^{(2)}$ can be interpreted as solving the moment condition for a randomized rule ${\bar{\eta}}(x)$ that predicts a positive classification with probability $(MK)^{-1} \sum_{r \in \mathcal{R}} \sum_{\mathsf{s} \in r} {\hat{\eta}}_{\tilde{\mathsf{s}}}(x)$.

I provide a CLT for the three estimators $(\hat{\theta}_{{\hat{\eta}}}^{(1)},\hat{\theta}_{{\hat{\eta}}}^{(2)},\hat{\theta}_{{\hat{\eta}}}^{(3)})$. Below, I establish my main regularity conditions.

(ref) is a standard moments condition for CLTs. (ref) is a mild stability condition on ${\tilde{\eta}}$. Importantly, ${\tilde{\eta}}$ is allowed to converge at any rate and to any limit $\eta^*_P$, which may depend on $P$. It holds, for example, if $$\psi_{\theta, {\tilde{\eta}}}(w) \xrightarrow{P} \psi_{\theta, {\eta^{*}_{P}}}(w)$$ pointwise for every $w \in \mathcal{W}$ and $\theta \in \Theta'$. This condition is more interpretable but stronger than required (see (ref) in (ref)). (ref) differs from the typical approach in the double machine learning literature where faster convergence rates (often $n^{-1/4}$) are required for nuisance functions chernozhukov2018double. The key difference between the two approaches is that I target a different, data-dependent parameter.

My CLT relies on the additional technical regularity conditions (ref), which I delay to (ref). This assumption adapts standard conditions for consistency and asymptotic normality of Z-estimators to the context of split-sample estimators van2000asymptotic,van2023weak. This is a weak assumption that holds in many settings, and it mostly concerns the choice of $\psi_{\theta,\eta}$. First, it assumes that the classes $\mathcal{F}_\eta = \left\{ \psi_{\theta, \eta, j} : \theta \in \Theta' \right\}$ are Donsker, which restricts complexity along $\theta \in \Theta'$ but does not restrict the complexity of $H$. Second, it requires $\hat{\theta}_{{\hat{\eta}}}^{(j)}$ to nearly solve the moment conditions, and $\theta_{{\hat{\eta}}}^{(j)}$ to be unique and well-separated zeroes of the population moment conditions. Finally, it assumes that $\Psi_{\eta}$ is differentiable in $\theta$ for $\eta \in H$, and the Jacobian is continuous in $\eta$ around ${\eta^{*}_{P}}$. (ref) holds, for example, in the case of sample averages ((ref)), or the “fraction in poverty by tercile” estimator in the poverty prediction application in (ref).

(ref) is the first main result of this paper.

The limiting variance $V_{{\eta^{*}_{P}}}$ is the product of two terms, the scalar $V_{\bar{M}, K}$ and a positive semidefinite matrix. The choice of split-sample procedure only affects $V_{{\eta^{*}_{P}}}$ through $V_{\bar{M}, K}$, which acts as a variance-inflating term since $V_{\bar{M}, K} \ge 1$. When using a single split ($K=1, \bar{M}=1$), the asymptotic variance is inflated by $V_{\bar{M}, K} = \pi^{-1}$, where $\pi$ is the fraction of the sample used to evaluate $\hat{\theta}_{{\hat{\eta}}}^{(j)}$ (as opposed to training ${\hat{\eta}}$). This occurs because $\hat{\theta}_{{\hat{\eta}}}^{(j)}$ is calculated from only $b = \pi n$ observations. When using repeated sample-splitting ($K = 1$ and $\bar{M} > 1$), $V_{\bar{M}, K} = \bar{M}^{-1} \pi^{-1} + \bar{M}^{-1} (\bar{M} - 1)$ is an average of $\pi^{-1}$ and $1$ with weights proportional to $1$ and $\bar{M} - 1$. This occurs since each observation is picked a different number of times across splits for calculating $\hat{\theta}_{{\hat{\eta}}}^{(j)}$. A larger number of repetitions leads to more balance in how often each observation is selected, and $V_{\bar{M}, K}$ decreases with larger $\bar{M}$. In fact, if $\bar{M} = \infty$, there is perfect balance in large samples and $V_{\bar{M}, K} = 1$. When using cross-fitting ($K > 1$), all observations are used an equal amount of times, and $V_{\bar{M}, K} = 1$. For intuition on this result, consider the particular case of sample averages ((ref)). In this case, $$\hat{\theta}_{{\hat{\eta}}} = \frac{1}{M} \sum_{r \in \mathcal{R}} \frac{1}{K} \sum_{\mathsf{s} \in r} \frac{1}{b} \sum_{i \in \mathsf{s}} f_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}(W_i)$$ is the same for $j=1,2,3$. If $\bar{M}=K=1$, $\hat{\theta}_{{\hat{\eta}}}$ averages over $b = \pi n$ observations. If $\bar{M} > 1$ and $K=1$, different observations are picked by splits $\mathsf{s}$ a different, random amount of times, and larger $\bar{M}$ leads to more balance. If $K > 1$, $\frac{1}{K} \sum_{\mathsf{s} \in r} \frac{1}{b} \sum_{i \in \mathsf{s}} f_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}(W_i)$ is an average over all observations, the entire sample is used equally, and the variance-inflation term is minimum. Hence, the asymptotic variance is minimized using cross-fitting with any number of folds $K > 1$ and repetitions $\bar{M}$.

(ref) appears to be new. The literature on risk estimation via cross-validation provides related results establishing asymptotic normality for sample average estimators based on multiple splits. Like my approach, these CLTs target data-dependent parameters, though they rely on different types of assumptions, and focus on the specific case of sample-averages. dudoit2005asymptotics consider estimators that average over all possible splits or cross-splits of the sample, assume bounded loss function and requires ${\hat{\eta}}$ to be loss consistent for a risk minimizing function, whereas I assume ${\hat{\eta}}$ converges to any limit. austern2020asymptotics,bayle2020cross provide CLTs under rate assumptions on the algorithmic stability of ${\hat{\eta}}$. bayle2020cross provides two CLTs using estimators based on a single repetition of cross-fitting, one relying on rate condition for algorithmic stability, and the second requires a “conditional variance convergence” assumption that they verify using rates for loss stability. My result does not require verifying a loss stability condition, which may not be satisfied in some high-dimensional settings bates2024cross, and my result allows for any ML algorithm as long as (ref) holds. ledell2015computationally provides a CLT for the particular case of estimating the area under the curve (AUC) measure via cross-validation.

A different class of related results are CLTs with cross-fitting for parameters that are not data-dependent. These approaches require stronger conditions on ${\hat{\eta}}$, such as requiring ${\hat{\eta}}$ to converge in probability at some specified rate, typically $n^{-1/4}$ luedtke2016statistical,chernozhukov2018double,benkeser2020improved,imai2025statistical. (ref) requires no complexity restrictions, and assumes ${\hat{\eta}}$ converges in probability to any limit at any rate.

A central limit theorem for the class of split-sample Z-estimators appears to be new. The characterization of the asymptotic variance, specifically how the variance-inflating term $V_{\bar{M}, K}$ depends on the number of splits $\bar{M}$ when $K=1$, also appears to be new. The proof uses a new CLT for split-sample empirical stated in (ref), which also appears to be new and may be of independent interest.

Comparison of Split-Sample Procedures

I compare the four split-sample procedures (sample-splitting, cross-fitting, repeated sample-splitting, and repeated cross-fitting) in terms of statistical power, modeling power, reproducibility, and computation time.

Cross-fitting and repeated cross-fitting, as well as repeated sample-splitting with $\bar{M} = \infty$, all exhibit the highest statistical power since they all minimize the variance of the asymptotic distribution in (ref). Repeated sample-splitting with $1 < \bar{M} < \infty$ comes second, and sample-splitting yields the largest variance.

I say that an estimator has better modeling power than another if the models in $\left( {\hat{\eta}}_{\tilde{\mathsf{s}}_{m,k}} \right)_{m \in \left[ M \right], k \in \left[ K \right]}$ are trained using larger datasets. Using more data for training typically leads to models with smaller expected loss, as I formalize in (ref). For sample-splitting or repeated sample-splitting, modeling power increases by picking a smaller $b$ (and $\pi$), so that more data is used to train each ${\hat{\eta}}_{\tilde{\mathsf{s}}_{m,k}}$. However, if $\bar{M} < \infty$, a smaller $b$ leads to smaller statistical power, since fewer data are used as evaluation sample at each split. When using cross-fitting, modeling power increases with $K$, since $b=n/K$. In this case, the returns to increasing $K$ are diminishing. For example, if $K=2$, ${\hat{\eta}}_{\tilde{\mathsf{s}}_{m,k}}$ is calculated with $50$% of the sample, and this fraction raises to $90$% with $K=10$. If $K=20$, however, the fraction only raises by another $5$%. Although a large value of $K$ or small value of $\pi$ (when $K=1$) lead to better modelling power, my asymptotic framework takes these quantities as fixed. This means that the quality of the asymptotic approximation may be poor if $K$ is large (or $\pi$ small) relative to the sample size. For example, my asymptotic framework does not accommodate for leave-one-out cross-fitting, that is, $K=n$.

I formalize the fact that increasing $M$ leads to better reproducibility properties in (ref). For example, as $M$ increases, it becomes more likely that two researchers using the same dataset but different random splits will reach the same conclusion about statistical significance of $\theta_{\hat{\eta}}$. Although I make no formal comparison between the cases $K=1$ and $K>1$ in terms of reproducibility, I note that ritzwoller2023reproducible documented the difference in variance between repeated sample-splitting and cross-fitting in an earlier draft.\footnote{This discussion appears in the second version at \url{https://arxiv.org/abs/2311.14204} (dated December 9, 2023).} Comparing cross-fitting with $M$ repetitions to sample-splitting with $KM$ repetitions, they argued that in principle it is possible that split-sample estimators have smaller variance conditional on data when $K=1$ instead of $K>1$, but show empirical evidence that cross-fitting typically leads to better reproducibility.

(ref) summarizes the comparison of the four procedures.

The choices of $M$, $K$, and $\pi$ (when $K=1$) involve tradeoffs. Statistical power is maximized when $K>1$ or $\bar{M} = \infty$ ((ref)), and the reproducibility properties improve with larger $M$ and are ambiguously affected by $K$, despite empirical evidence that $K>1$ usually leads to better properties ritzwoller2023reproducible. For $K=1$ and $\bar{M} < \infty$, there is a tradeoff between statistical and modelling powers, unlike with cross-fitting. A larger $M$ is always beneficial in terms of reproducibility (and statistical power when $K=1$), but this comes at the cost of higher computation time. Hence, I recommend choosing $M$ as large as computationally convenient, and $K>1$ but small, since that provides valid asymptotic inference, maximum statistical power, and likely better reproducibility properties. In (ref), I provide a measure to assess whether a given $M$ is sufficiently large to ensure reproducibility of p-values calculated from split-sample Z-estimators.

Inference on Split-Sample Estimands

The CLT in (ref) can be directly applied to conduct inference on many split-sample estimands. However, confidence intervals based on the normal approximation may fail to cover $\theta_{\hat{\eta}}$ at the nominal level in some important cases of interest. First, in (ref), I consider inference when the normal approximation is asymptotically exact, and discuss why this approximation may not be precise in some contexts. Then, in (ref), I propose a new approach for the particular cases of inference on comparisons between models, which explicitly accounts for the dependence across splits.

I discuss in (ref) that a typical case when the normal approximation CI may have coverage probability smaller than nominal is when the variance of a moment function is allowed to be zero in the limit. I provide a general method for inference that covers this case in (ref), by exploring the faster-than-$\sqrt{n}$ convergence rate of the empirical moment functions and introducing a tuning parameter. I also discuss in (ref) that although (ref) considers the specific case of comparing two models, the arguments developed in that section apply more broadly, covering other cases such as the Generic ML approach of chernozhukov2025generic (see (ref)).

Inference from Normal Approximation

Consider the problem of conducting inference on $h(\theta_{\hat{\eta}})$, where $\theta_{\hat{\eta}}$ is any of the split-sample Z-estimands in (ref), and $h$ is any scalar differentiable function with row-vector of partial derivatives $\dot{h}(\theta_{{\eta^{*}_{P}}}) \neq 0$. This encompasses many cases of interest, for example when $h(\theta_{\hat{\eta}})$ is a subset of the vector $\theta_{\hat{\eta}}$ or a linear combination of its entries, as in the application of (ref). An application of (ref) and the delta-method yields

where $\hat{\theta}_{\hat{\eta}}$ is a Z-estimator as in (ref), and $$\sigma^2_{{\eta^{*}_{P}}} = \dot{h}(\theta_{{\eta^{*}_{P}}}) V^*_{{\eta^{*}_{P}}} \dot{h}(\theta_{{\eta^{*}_{P}}})^T.$$ If $\sigma^2_{{\eta^{*}_{P}}} > 0$, one can calculate the plug-in estimator

where $\hat{V}_{{\hat{\eta}}}$ is given in (ref) in (ref), and the confidence interval

contains $h(\theta_{\hat{\eta}})$ with probability approaching $1 - \alpha$, where $z_\alpha$ is the $\alpha$-th quantile of the standard normal distribution.

(ref) assumes the existence of a consistent estimator of $\dot{\Psi}_{{\eta^{*}_{P}}}$. If $\psi_{\theta, \eta}(w)$ is differentiable in $\theta$, this assumption is satisfied by the plug-in estimator defined in (ref) in (ref) under a uniform integrability condition on this derivative. Otherwise, consistent estimators of $\dot{\Psi}_{{\eta^{*}_{P}}}$ can typically be constructed on a case-by-case basis hansen2022econometrics. Note that the probability in (ref) is taken over both the random estimand $h(\theta_{\hat{\eta}})$ and the CI.

(ref) implies that (ref) contains $h(\theta_{\hat{\eta}})$ with probability approaching $1 - \alpha$ in many settings. However, in some cases, (ref) may not cover $h(\theta_{\hat{\eta}})$ with nominal probability, as illustrated in the two examples below.

(ref) have two features in common: the normal approximation CI may undercover $\theta_{\hat{\eta}}$ only when ${\eta^{*}_{P}}(x)$ is constant, and one of the empirical moment equations evaluated at the true parameter converges to zero at a rate faster than $n^{-1/2}$:

where $\psi_{\theta,\eta,j}$ is the $j$-th entry of the vector $\psi_{\theta,\eta}$. In (ref), I develop an approach that can be used to test whether ${\eta^{*}_{P}}(x)$ is constant, and although I focus on the particular case of comparing the performance between two models, the arguments apply more broadly and could be used to provide a valid CI for the problems in (ref) under the same conditions of (ref). In (ref), I establish a general approach to inference on $\theta_{\hat{\eta}}$ that allows (ref) to happen. The approach explores the faster-than-$\sqrt{n}$ convergence rate to provide an asymptotically uniformly valid CI by introducing a tuning parameter.

Inference on Model Comparisons

In several applications, the goal is not only to create a new model ${\hat{\eta}}$ and assess some property $\theta_{\hat{\eta}}$, but to compare such properties between two models. For example, if $\theta_{\hat{\eta}}$ is a measure of accuracy such as the mean squared error ((ref)), one might want to infer if ${\hat{\eta}}$ has better performance than a baseline model that predicts the sample mean of $Y$ for all observations. This is the case in the application of (ref), where the goal is to assess whether a random forest model has predictive power for poverty, that is, whether it achieves smaller MSE than using the sample average. Alternatively, one might want to compare the performance of using different machine learning algorithms, such as training ${\hat{\eta}}$ with neural networks versus random forests. I show that the CLTs of the previous sections give a valid inference approach when both models do not have similar performances in large samples. However, if the models have similar performance, the asymptotic distribution of the difference in performance is degenerate at the $\sqrt{n}$ rate, and CIs based on the asymptotic approximation may fail to control size. In this section, I build on the CLT of (ref) to develop an inference method that is valid for both cases. Although this section focuses on the particular case of comparing two models, I discuss in the end of (ref) that the arguments developed in this section apply more broadly.

The setting is as follows. $\hat{\theta}_{\hat{\eta}}$ denotes any of the estimators $(\hat{\theta}_{\hat{\eta}}^{(1)},\hat{\theta}_{\hat{\eta}}^{(2)},\hat{\theta}_{\hat{\eta}}^{(3)})$ of (ref), assumed to be a scalar ($d=1$) (alternatively, one could consider a scalar transformation $h(\hat{\theta}_{\hat{\eta}})$ as in (ref)). I refer to the parameter $\theta_{\hat{\eta}}$ (defined analogously) as a performance measure for expositional convenience, though the results apply more generally. I focus on comparing $\theta_{\hat{\eta}}$ to the performance $\theta_{\hat{b}}$ of a baseline model $\hat{b} \in H$ computed using the entire sample, that is, without forms of sample-splitting. $\hat{b}$ is assumed to come from a parametric model, and it can be, for example, the sample average $\hat{b}(x) = n^{-1} \sum_{i=1}^{n} Y_i$ in (ref). Following the notation of (ref), $\theta_{\hat{b}}$ is the unique solution for $\theta$ in $\Psi_{\hat{b}}(\theta) = 0$, i.e., $$\Psi_{\hat{b}}(\theta_{\hat{b}}) = 0.$$ Similarly, the estimator $\hat{\theta}_{\hat{b}}$ is a (near) zero of the empirical moment condition, $$\hat{\theta}_{\hat{b}} \in \arg\min_{\theta \in \Theta} \left\lVert \frac{1}{n} \sum_{i=1}^{n} \psi_{\theta,\hat{b}}(W_i) \right\rVert.$$ In (ref), I discuss how to extend the current setting for comparing $\theta_{\hat{\eta}}$ to the performance of another model ${\hat{\eta}}'$ computed with the same split-sample approach as ${\hat{\eta}}$. Let $$\mathcal{S} = \left( \mathsf{s}_{m,k} \right)_{m \in \left[ M \right], k \in \left[ K \right]}$$ be a collection of $M K$ splits of the sample, that is, a vectorization of $\mathcal{R}$ defined in (ref). Notice that each $\mathsf{s} \in \mathcal{S}$ is associated with a model ${\hat{\eta}}_{\tilde{\mathsf{s}}}$, as in (ref).

To see the challenge of conducting inference based on $\hat{\theta}_{\hat{\eta}} - \hat{\theta}_{\hat{b}}$, consider a simplified setting where each $\hat{\theta}_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}$ (as in (ref)) is a sample average, that is, $\psi_{\theta,\eta}(w) = f_\eta(w) - \theta$ for some $f_\eta$ and $\hat{\theta}_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} = \left| \mathsf{s} \right|^{-1} \sum_{i \in \mathsf{s}} f_{{\hat{\eta}}_{\tilde{\mathsf{s}}}}(W_i)$. The CLT in (ref) gives $$\sqrt{n} \left( \hat{\theta}_{\hat{\eta}} - \theta_{\hat{\eta}} \right) = \frac{1}{\sqrt{n}} \sum_{i=1}^{n} \left( f_{{\eta^{*}_{P}}}(W_i) - Pf_{{\eta^{*}_{P}}} \right) + o_P(1),$$ and the normal approximation gives an asymptotically valid CI for $\theta_{\hat{\eta}}$. Similarly, if $\hat{b}$ converges to some model $b_P \in H$, $$\sqrt{n} \left( \hat{\theta}_{\hat{b}} - \theta_{\hat{b}} \right) = \frac{1}{\sqrt{n}} \sum_{i=1}^{n} \left( f_{b_P}(W_i) - Pf_{b_P} \right) + o_P(1),$$ and these two results can be combined to construct a CI for $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ based on a normal approximation. However, if the baseline model $b_P$ is the same as ${\eta^{*}_{P}}$, both estimators have the same limit, the difference

has a degenerate limit in probability, and the CLT of (ref) does not inform how to compute a CI for $\theta_{{\hat{\eta}}} - \theta_{\hat{b}}$.

First, I develop a test for whether any of the models ${\hat{\eta}}_{\tilde{\mathsf{s}}}$ perform better than $\hat{b}$, then show how this test can be used to construct a CI for $\theta_{{\hat{\eta}}} - \theta_{\hat{b}}$. Both results build on my CLT for Z-estimators.

A Multivariate One-sided Test for Model Differences

From (ref), the asymptotic distribution of $\hat{\theta}_{{\hat{\eta}}} - \hat{\theta}_{\hat{b}}$ centered around the parameter of interest $\theta_{{\hat{\eta}}} - \theta_{\hat{b}}$ is degenerate at the $n^{-1/2}$ rate if $b_P = {\eta^{*}_{P}}$. Yet, for each split $\mathsf{s} \in \mathcal{S}$,

has a non-degenerate limit since the first average does not include observations $i \in \tilde{\mathsf{s}}$. I explore this fact to construct a test of whether any model ${\hat{\eta}}_{\tilde{\mathsf{s}}}$ has better performance than $\hat{b}$, then develop a CI for $\theta_{{\hat{\eta}}} - \theta_{\hat{b}}$ in the following subsection.

Consider the hypothesis test

If $\theta_{\eta}$ is a measure of performance such as the mean squared error, having $\theta_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} - \theta_{\hat{b}} < 0$ means that ${\hat{\eta}}_{\tilde{\mathsf{s}}}$ performs better than $\hat{b}$. The hypotheses $H_{0,{\hat{\eta}}}$ and $H_{A,{\hat{\eta}}}$ depend on ${\hat{\eta}}$ due to the data-dependent parameter of interest $\theta_{\hat{\eta}}$. Testing such hypotheses is analogous to constructing a confidence interval for a data-dependent parameter as in (ref). Let $$\delta_{\hat{\eta}} = \left( \theta_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} - \theta_{\hat{b}} \right)_{\mathsf{s} \in \mathcal{S}},$$ and similarly define $${\hat{\delta}}_{\hat{\eta}} = \left( \hat{\theta}_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} - \hat{\theta}_{\hat{b}} \right)_{\mathsf{s} \in \mathcal{S}}.$$ An application of (ref) gives $$\sqrt{n} \left( {\hat{\delta}}_{\hat{\eta}} - \delta_{\hat{\eta}} \right) \leadsto \mathcal{N} \left( 0, \Sigma \right),$$ for some nonzero $\Sigma$ that can be consistently estimated with $\hat{\Sigma}$ (see equation (ref) in (ref)). Since splits reuse observations, the off-diagonal terms of $\Sigma$ explicitly incorporate the dependence across splits.

Denote by $\hat{\sigma}^2_\mathsf{s}$ the entry of the main diagonal of $\hat{\Sigma}$ associated with $\mathsf{s} \in \mathcal{S}$, that is, with the term $\hat{\theta}_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} - \hat{\theta}_{\hat{b}}$. I propose computing the test-statistic $$T({\hat{\delta}}_{\hat{\eta}}, n^{-1} \hat{\Sigma}) = \sum_{\mathsf{s} \in \mathcal{S}} \left( \min \left\{ \sqrt{n} \frac{\hat{\theta}_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} - \hat{\theta}_{\hat{b}}}{\hat{\sigma}_{\mathsf{s}}}, 0 \right\} \right)^2.$$ This type of test statistic has been considered for example in chernozhukov2007estimation,romano2008inference,andrews2009validity,romano2010inference in the context of moment inequalities. Critical values $\hat{c}_{1 - \alpha}$ can be computed via Monte Carlo: simulate $Z \sim \mathcal{N}(0,\hat{\Sigma})$ and estimate $\hat{c}_{1 - \alpha}$ as the $1-\alpha$ quantile of $T(Z, n^{-1} \hat{\Sigma})$. I note that, alternatively, one could use the likelihood ratio test statistic.

Asymptotic exactness of this test under the least favorable null follows from similar conditions to (ref), established below.

(ref) consists of more technical conditions, which are delayed to the appendix for ease of exposition. For example, they extend the Z-estimator assumptions on ${\eta^{*}_{P}}$ to $b_P$. (ref) requires the limiting variance of $\sqrt{n} (\hat{\theta}_{\hat{\eta}} - \theta_{\hat{\eta}})$ to be positive, and (ref) defines the requirements on the baseline (parametric) estimator $\hat{b}$. It holds, for example, if $\hat{b}$ belongs to a Donsker class with probability approaching one, which typically happens for parametric models such as the sample average $\hat{b}(x) = n^{-1} \sum_{i=1}^{n} Y_i$.

(ref) appears to be new. It establishes size control: the probability of rejecting the null hypothesis, conditional on it being true, does not exceed $\alpha$ in large samples. Note that the probabilities in (ref) are not random objects, they integrate over the distribution of the data conditional on the events $\delta_{\hat{\eta}} \ge 0$ or $\delta_{\hat{\eta}} = 0$. Alternative approaches for testing across multiple splits of the sample typically aggregate p-values or confidence intervals computed separately for each split, without accounting for the dependence structure across splits chernozhukov2025generic,gasparin2025combining. For example, chernozhukov2025generic propose aggregating the median of p-values or CIs across splits. Because these methods do not incorporate the correlation across splits, they are conservative in most data-generating processes, as they guard against the worst-case dependence structure. In contrast, my approach explicitly accounts for the dependence across splits, which enables the test to achieve exactness under the least favorable null $\delta_{\hat{\eta}} = 0$ in a uniform sense across DGPs. The proof is made possible by the decomposition in (ref), which follows from the new CLT in (ref).

The result above requires the probability of the conditioning event to be bounded away from zero using the constant $\bar{c}_2 > 0$. This could lead to an apparent uniformity issue for sequences of DGPs $(P_n)_{n \ge 1}$ with $P_n(\delta_{\hat{\eta}} \ge 0) \to 0$, for example. For such sequences, the probability of rejecting the null conditional on the null being true could be greater than $\alpha$. This is not, however, an issue for empirical practice: for such sequences the probability of being under the null itself converges to zero. Incorrectly rejecting the null is not a concern when the probability of the null being true is zero.

A Confidence Interval for the Average Performance

I construct a new confidence interval for $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ based on two insights from the previous subsections. The first is that a CI based on the normal approximation using (ref) is asymptotically exact if $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ converges in probability to a value different from zero, since in this case the terms in (ref) do not cancel out. The second insight is that the case $\theta_{\hat{\eta}} - \theta_{\hat{b}} \xrightarrow{P} 0$ is closely connected with the null hypothesis of the one-sided test developed in the previous subsection. Hence, my CI consists of using the normal approximation if the one-sided test is rejected, and an extended CI in case it is not.

Define the normal approximation CI $$\widehat{\rm CI}_{\alpha,\mathcal{N}} = \left[ \left( \hat{\theta}_{\hat{\eta}} - \hat{\theta}_{\hat{b}} \right) - z_{1 - \alpha / 2} \frac{\hat{\sigma}_{\hat{\delta}}}{\sqrt{n}}, \left( \hat{\theta}_{\hat{\eta}} - \hat{\theta}_{\hat{b}} \right) + z_{1 - \alpha / 2} \frac{\hat{\sigma}_{\hat{\delta}}}{\sqrt{n}} \right],$$ where $\hat{\sigma}_{\hat{\delta}}$ is a standard error for $\hat{\theta}_{\hat{\eta}} - \hat{\theta}_{\hat{b}}$ (see equation (ref) in (ref)), and an extended CI $$\widehat{\rm CI}_{\alpha,{\rm ext}} = \operatorname{Conv}\left( \widehat{\rm CI}_{\alpha,\mathcal{N}} \cup \{0\} \right),$$ where $\operatorname{Conv}(\cdot)$ denotes the convex hull, that is, $\widehat{\rm CI}_{\alpha,{\rm ext}}$ has all the elements in $\widehat{\rm CI}_{\alpha,\mathcal{N}}$, $0$, and all elements in between. The final CI is given by $$\widehat{\rm CI}_{\alpha} =

$$ $\widehat{\rm CI}_{\alpha}$ is based on a pre-test, using different inference approaches depending on whether the one-sided test is rejected or not. This construction is motivated by the following facts, which are formalized in \Cref{th.diff_pointwise,th.diff_ci1,th.diff_ci1_power}. If $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ converges in probability to a negative value, $P(T({\hat{\delta}}_{\hat{\eta}}, \hat{\Sigma}) > \hat{c}_{1 - \alpha}) \to 1$, and $\widehat{\rm CI}_{\alpha,\mathcal{N}}$ is used, which is asymptotically exact. If $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ converges in probability to a positive value, $P(T({\hat{\delta}}_{\hat{\eta}}, \hat{\Sigma}) > \hat{c}_{1 - \alpha}) \to 0$, $\widehat{\rm CI}_{\alpha,\mathcal{N}}$ is asymptotically exact but $\widehat{\rm CI}_{\alpha,{\rm ext}}$ is used, which is valid since it is wider than $\widehat{\rm CI}_{\alpha,\mathcal{N}}$, although conservative. This asymmetric construction is a choice, which reflects the motivating problem of this section of learning whether the new model ${\hat{\eta}}$ performs better (instead of worse) than the baseline model $\hat{b}$. Finally, if $\theta_{\hat{\eta}} - \theta_{\hat{b}} \xrightarrow{P} 0$, intuitively $P(T({\hat{\delta}}_{\hat{\eta}}, \hat{\Sigma}) > \hat{c}_{1 - \alpha})$ should be close to $\alpha$ given \Cref{th.diff_onesided_test}. If that happens, $\widehat{\rm CI}_{\alpha,{\rm ext}}$ covers $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ with high probability since it includes $0$, the limit of $\theta_{\hat{\eta}} - \theta_{\hat{b}}$. However, this guarantee depends on additional conditions as I discuss next, since $P (\delta_{\hat{\eta}} \ge 0)$ may not converge to one even if $\delta_{\hat{\eta}} \xrightarrow{P} 0$.

First, I show that $\widehat{\rm CI}_{\alpha}$ is valid pointwise in $P \in \mathcal{P}$, assuming that if ${\eta^{*}_{P}} = b_P$, then the parametric model is well-specified in the sense that it minimizes the error $\theta_\eta$ in $\eta$, that is, $\theta_\eta \ge \theta_{b_P}$ for all $\eta \in H$. Then, I establish conditions under which $\widehat{\rm CI}_{\alpha}$ is valid uniformly in $P \in \mathcal{P}$.

Further, I show that $\widehat{\rm CI}_{\alpha}$ is asymptotically valid for most sequences of $\theta_{{\hat{\eta}}} - \theta_{\hat{b}}$, and discuss why it may fail for specific sequences. Then, I establish that the additional condition (ref) is sufficient for $\widehat{\rm CI}_{\alpha}$ to be asymptotically uniformly valid in $P \in \mathcal{P}$. Later, I propose a modification to $\widehat{\rm CI}_{\alpha}$ that gives uniform validity under only (ref).

Under (ref), $\widehat{\rm CI}_{\alpha}$ covers $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ when this difference is positive or “sufficiently” negative, with $\mathcal{P}_{\bar{c}_3,\bar{c}_4}$ only requiring this event to happen with positive probability. If $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ converges to any negative value, coverage is asymptotically exact ((ref)). If it converges to a positive value, similarly, the normal approximation CI is exact, and the extended $\widehat{\rm CI}_{\alpha,{\rm ext}}$ is conservative. Failure of coverage may happen only if $\theta_{\hat{\eta}} - \theta_{\hat{b}} \xrightarrow{P} 0^-$, that is, it converges to zero “from the left”. For such sequences, the components of $\delta_{{\hat{\eta}}}$, $\theta_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} - \theta_{\hat{b}}$, may be enough negative so that the one-sided test rejects the null with high probability, but since they converge to zero, the terms in (ref) cancel out, and the normal approximation CI may undercover. Importantly, $\widehat{\rm CI}_{\alpha}$ is valid in the case of interest $\theta_{\hat{\eta}} - \theta_{\hat{b}} \ge 0$, that is, when ${\hat{\eta}}$ performs equally or worse than the baseline model $\hat{b}$. This is the case, for example, when the parametric model is well-specified, as in (ref), since $\sqrt{n} \left( \theta_{\hat{b}} - \theta_{b_P} \right) \xrightarrow{P} 0$ from (ref). Hence, $\widehat{\rm CI}_{\alpha}$ may overstate the advantage of ${\hat{\eta}}$ when it slightly outperforms $\hat{b}$, but not when their performances are equal or when ${\hat{\eta}}$ performs worse.

(ref) rules out the problematic sequences by ensuring that if $\theta_{\hat{\eta}} - \theta_{\hat{b}} \xrightarrow{P} 0^-$, $\theta_{\hat{\eta}}$ is close enough to $\theta_{\hat{b}}$ in large samples so that the one-sided test does not reject with probability higher than $\alpha$. It is motivated by the fact that machine learning algorithms typically penalize deviations from the mean. If there is little signal to be learned by ${\hat{\eta}}$, that is, $\theta_{\eta^{*}_{P_n}}$ is close to $\theta_{b_{P_n}}$, it may be reasonable to expect that regularization will make the estimates ${\hat{\eta}}$ closer to $\hat{b}$ than to ${\eta^{*}_{P_n}}$. For example, in the case of estimating a linear model with the Lasso, if the true coefficients are very small, penalization leads to estimated coefficients exactly equal to $0$ with high probability zhao2006model,zhang2008sparsity,wuthrich2023omitted. However, this assumption may not lead to a good approximation for the behavior of DGPs where $\theta_{\eta^{*}_{P_n}}$ is sufficiently distant from $\theta_{b_{P_n}}$ and ${\hat{\eta}}$ is estimated with no or little regularization.

Next, I provide an alternative, more conservative CI that gives uniform coverage without relying on (ref). It deals with sequences with $\theta_{\hat{\eta}} - \theta_{\hat{b}} \xrightarrow{P} 0^{-}$ by modifying $\widehat{\rm CI}_{\alpha}$ to be more conservative in the one-sided test. For any $\bar{c}_5 > 0$, consider the modified version of the test in (ref):

$\bar{c}_5$ represents a degree of slackness on how large $-(\theta_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} - \theta_{\hat{b}})$ has to be to reject the null hypothesis. The final CI is given by $$\widehat{\rm CI}_{\alpha}' =

$$ and the critical value $\hat{c}_{1 - \alpha}$ is the same as before. A large $\bar{c}_5$ gives more robustness in finite samples in the sense that $$P \left( \left( \theta_{\hat{\eta}} - \theta_{\hat{b}} \right) \in \widehat{\rm CI}_{\alpha}' \right)\alpha}'}$$ is (weakly) increasing in $\bar{c}_5$. On the other hand, a large $\bar{c}_5$ leads to less power. Importantly, this approach is not necessary if the goal is to test the null $H_{0,{\hat{\eta}}}: \theta_{{\hat{\eta}}} - \theta_{\hat{b}} = 0$, since this cased is covered by \Cref{th.diff_ci1}. The modified confidence interval $\widehat{\rm CI}_{\alpha}'$ is intended for researchers who may want to be careful not to overestimate the magnitude of $\theta_{\hat{\eta}} - \theta_{\hat{b}}$ when it is small but negative.

Reproducibility

The split-sample estimators and estimands defined in (ref) depend not only on the algorithm used to estimate ${\hat{\eta}}$, but also on the specific splits of the sample $\mathcal{R}$. In applications, this may lead to the undesirable phenomenon that different researchers with the same dataset, using different random splits $\mathcal{R}$, may reach different conclusions in terms of statistical significance. Intuitively, by averaging over multiple splits, this phenomenon becomes less likely. In this section, I first formalize this intuition by establishing basic reproducibility properties of split-sample Z-estimators. Then, I develop a measure that quantifies the reproducibility of p-values from hypothesis tests based on Z-estimators for a given number of repetitions $M$.

Basic Reproducibility Properties

I establish two basic reproducibility properties for the three versions of split-sample Z-estimators defined in (ref). The two results formalize the notion that, for fixed $n$, choosing to use a larger number of repetitions $M$ improves reproducibility of the estimators. The results exploit the fact that $\hat{\theta}_{\hat{\eta}}^{(1)}$ and $\hat{\theta}_{\hat{\eta}}^{(3)}$ are averages over $M$ independent repetitions $r \in \mathcal{R}$. For the second estimator, I use the fact that, conditional on the data $D$, $\hat{\theta}_{\hat{\eta}}^{(2)}$ is still a Z-estimator where the “observations” are the splits $r \in \mathcal{R}$ and the target parameter is the value of $\theta$ that solves the moment condition averaged over all possible splits. This allows me to explore large $M$ properties of $\hat{\theta}_{\hat{\eta}}^{(2)}$ using arguments similar to those applied to Z-estimators (e.g., Theorem 5.9 in van2000asymptotic). For $\hat{\theta}_{\hat{\eta}}^{(2)}$, I require an additional technical condition which I delay to (ref). This assumption is analogous to standard conditions for proving consistency of Z-estimators, and holds, for example, if $\Theta'$ is bounded, $\psi_{\theta,\eta}$ is Lipschitz in $\theta$ with a Lipschitz constant that does not depend on $\eta$ or $w$, and if the solution to the moment condition averaged over all possible splits is unique.

The first reproducibility property is that, for fixed $n$, the variance of the Z-estimators conditional on the data converge to zero as $M$ grows. Conditional on the data, the estimators vary only due to the random partitioning. This approximates the behavior of the estimators when the number of repetitions $M$ is chosen to be large. This guarantees that two researchers with the same dataset and different splits will calculate estimators that are arbitrarily close to each other with high probability for large enough $M$.

For the estimators $\hat{\theta}_{{\hat{\eta}}}^{(1)}$ and $\hat{\theta}_{{\hat{\eta}}}^{(3)}$, I show that the conditional variance is strictly decreasing in $M$. This establishes a stronger property than the asymptotic result in (ref): not only does reproducibility improve as $M \to \infty$, but every increase in $M$ strictly reduces variance and thus improves reproducibility.

A Reproducibility Measure

I propose a reproducibility measure for p-values from hypothesis tests based on transformations of split-sample Z-estimators. Specifically, I study reproducibility of the p-value for testing $H_{0,{\hat{\eta}}}: h(\theta_{{\hat{\eta}}}) = \tau$ versus $H_{A,{\hat{\eta}}}: h(\theta_{{\hat{\eta}}}) \neq \tau$ (and its one-sided versions) for $h: \Theta \to \mathbb{R}$ differentiable. The hypotheses $H_{0,{\hat{\eta}}}$ and $H_{A,{\hat{\eta}}}$ depend on ${\hat{\eta}}$ since the parameter of interest, $\theta_{\hat{\eta}}$, depends on ${\hat{\eta}}$. Testing this hypothesis is analogous to constructing a CI for $\theta_{\hat{\eta}}$: in fact, inverting this test for all values of $\tau$ at significance level $\alpha$ gives the confidence interval of (ref).

I begin by defining the reproducibility measure, then describe the asymptotic framework I use and the technical challenges involved. Finally, I establish the limit distribution of the difference of t-statistics constructed from different random splits, and apply this result to construct the reproducibility measure. As in (ref), I consider $M$ repetitions of sample-splitting with $K$ folds ($K=1$ denotes repeated sample-splitting).

The goal of this section is to construct a measure $\hat{\delta}(\beta)$, for $\beta \in (0, 0.5)$, that satisfies $$P \left( p_2 > p_1 + \hat{\delta}(\beta) \Biggm| D \right) = \beta + o_P(1),$$ where $p_1$ and $p_2$ are p-values for $H_{0,{\hat{\eta}}}$ calculated with separate, independent splits. This measure provides the following guarantee: if a researcher calculates a p-value $p_1$ using one set of random splits, then a second researcher using the same dataset, but different splits, will obtain a p-value exceeding $p_1 + \hat{\delta}(\beta)$ with probability approximately $\beta$. This allows researcher 1 to assess whether their result would remain statistically significant without the computational cost of re-running the analysis. For example, if $p_1 < 0.05$ but $p_1 + \hat{\delta}(\beta) > 0.05$ for some small $\beta$, the researcher may need to increase $M$ to guarantee reproducibility of their finding.

I consider an asymptotic regime where both the number of repetitions $M$ and the sample size $n$ grow to infinity, which is the main technical challenge for proving validity of my reproducibility measure. An alternative framework is to consider the data $D$ fixed, let $M \to \infty$, and treat each repetition as an independent observation. Although this alternative regime facilitates statistical analysis, it provides asymptotic guarantees only when $M$ is large relative to $n$. In practice, choosing $M$ much larger than $n$ is often computationally intractable. My asymptotic framework better reflects much of empirical practice by allowing $M$ to grow slower than $n$, so that $M$ can be, for instance, a small fraction of $n$. The proofs of my results under this asymptotic regime rely on the CLT of (ref).

I focus on the estimator $\hat{\theta}_{{\hat{\eta}}} = \hat{\theta}_{{\hat{\eta}}}^{(2)}$ from (ref), and similar results can be extended to $\hat{\theta}_{{\hat{\eta}}}^{(1)}$ and $\hat{\theta}_{{\hat{\eta}}}^{(3)}$ using similar techniques. The $\hat{\theta}_{{\hat{\eta}}}^{(2)}$ case is much more challenging because, unlike $\hat{\theta}_{{\hat{\eta}}}^{(1)}$ and $\hat{\theta}_{{\hat{\eta}}}^{(3)}$, $\hat{\theta}_{{\hat{\eta}}}^{(2)}$ is not an average of $M$ independent terms conditional on the data.

The setting follows (ref). Additionally, let $\mathcal{R}_1$ and $\mathcal{R}_2$ be independent collections of $M$ splits of the data with $K$ folds (uniformly at random). Let ${\hat{\eta}}_1$ and ${\hat{\eta}}_2$ be calculated with $\mathcal{R}_1$ and $\mathcal{R}_2$ respectively, which leads to analogous definitions of $\hat{\theta}_{{\hat{\eta}}_{j}}$, $\theta_{{\hat{\eta}}_{j}}$, and $\hat{\sigma}_{{\hat{\eta}}_j}$ for $j=1,2$. Under the null hypothesis and the conditions of (ref), the t-statistic $$t_{{\hat{\eta}}_{j}} = \frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_{j}}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_{j}}} \leadsto \mathcal{N}(0,1),$$ where $\hat{\sigma}_{{\hat{\eta}}_{j}}$ is given as in (ref), $\dot{h}(\theta)$ is a row vector with the partial derivatives of $h(\theta)$ evaluated at $\theta$, and $\hat{V}_{{\hat{\eta}}}$ is a plug-in estimator for $V_{{\eta^{*}_{P}}}$ defined in (ref). Based on this result, one can calculate p-values $$p_j^{\pm} = 2 \Phi \left( - \left| \frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_j}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_j}} \right| \right),$$ $$p_j^+ = \Phi \left( \frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_j}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_j}} \right), \quad p_j^- = \Phi \left( -\frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_j}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_j}} \right),$$ for $H_{0,{\hat{\eta}}_j}: h(\theta_{{\hat{\eta}}_j}) = \tau$ versus $H_{A,{\hat{\eta}}}: h(\theta_{{\hat{\eta}}}) \neq \tau$ and its one-sided versions.

The asymptotic regime assumes $M^{-1} n \sigma^2_D = O_P(1)$, where $\sigma^2_D$ (defined in (ref) in the appendix) reflects the variance of $t_{{\hat{\eta}}_{1}}$ conditional on the data. Since $h(\hat{\theta}_{{\hat{\eta}}})$ and $\hat{\sigma}_{{\hat{\eta}}}$ converge to non-random quantities $h(\theta_{{\eta^{*}_{P}}})$ and $\sigma_{{\eta^{*}_{P}}}$ respectively, $\sigma^2_D \xrightarrow{P} 0$. Hence, the asymptotic regime requires $M \to \infty$ at a rate slower than $n$. The rate of convergence of $\sigma^2_D$ depends on the rate at which ${\tilde{\eta}} = \mathcal{A}(D)$ converges to ${\eta^{*}_{P}}$, and may be slow especially when ${\tilde{\eta}}$ is estimated nonparametrically. In (ref), I show that a safe guideline for achieving the reproducibility guarantees established below is to choose $M$ of comparable magnitude to $n$.

I characterize below a central limit theorem for the difference of t-statistics constructed using different splits, which is the main ingredient for deriving my reproducibility measure in (ref). Both results rely on the fairly technical (ref), stated in (ref). The key condition is a Donsker-type requirement on $\left\{ \Psi_{{\hat{\eta}}_{\tilde{\mathsf{s}}}} : \mathsf{s} \subseteq \left[ n \right] \right\}$ and $\psi_{\theta, \eta, i} \psi_{\theta, \eta, j}$. This condition holds, for example, if $\Theta'$ and $\psi_{\theta,\eta}$ are bounded and the cross products of the entries of $\psi_{\theta,\eta}$ are Lipschitz. Importantly, (ref) does not restrict the complexity of ${\hat{\eta}}$, it only restricts the complexity of the function classes over $\theta \in \Theta'$, and not over $\eta \in H$.

I introduce my reproducibility measure for each of the three tests (two-sided and both one-sided tests), where $\Phi$ is the standard normal cdf, and formalize their guarantees in (ref). $$\hat{\delta}^\pm(\beta) = 2 \Phi \left( -\left| \frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_1}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_1}} \right| - \frac{\sqrt{n} \hat{\sigma}_{D}}{\sqrt{M}} \Phi^{-1}(\beta/2) \right) - 2 \Phi \left( -\left| \frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_1}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_1}} \right| \right),$$ $$\hat{\delta}^+(\beta) = \Phi \left( \frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_1}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_1}} - \frac{\sqrt{n} \hat{\sigma}_{D}}{\sqrt{M}} \Phi^{-1}(\beta) \right) - \Phi \left( \frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_1}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_1}} \right),$$ $$\hat{\delta}^-(\beta) = \Phi \left( -\frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_1}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_1}} - \frac{\sqrt{n} \hat{\sigma}_{D}}{\sqrt{M}} \Phi^{-1}(\beta) \right) - \Phi \left( -\frac{\sqrt{n} (h(\hat{\theta}_{{\hat{\eta}}_1}) - \tau)}{\hat{\sigma}_{{\hat{\eta}}_1}} \right).$$

(ref) gives a novel measure of reproducibility for p-values based on split-sample Z-estimators. The guarantee of reproducibility in (ref) is inspired by the definition of $(\xi, \beta)$-reproducibility in ritzwoller2023reproducible. They provide an algorithm for deciding how many repetitions $M$ of the sample-splitting procedure are necessary to guarantee reproducibility of the average across split-sample statistics. This covers, for example, the estimators $\hat{\theta}_{{\hat{\eta}}}^{(1)}$ and $\hat{\theta}_{{\hat{\eta}}}^{(3)}$. My approach complements theirs by focusing on reproducibility of inference, examining p-value rather than average statistics. My results hold for $\hat{\theta}_{{\hat{\eta}}}^{(2)}$, and the arguments can easily be extended to $\hat{\theta}_{{\hat{\eta}}}^{(1)}$ and $\hat{\theta}_{{\hat{\eta}}}^{(3)}$. ritzwoller2023reproducible's procedure takes as input the desired level of reproducibility, and outputs the required number of repetitions $M$ that guarantees such reproducibility. My approach takes $M$ as input (assumed “large”), and outputs a measure of how much reproducibility is guaranteed by such $M$. The asymptotic regimes also differ: ritzwoller2023reproducible takes the data as fixed and considers that the desired threshold for the variability of the average split-sample statistic is small, while my framework considers $n$ and $M$ large.

The result in (ref) relies on choosing $M$ such that $M^{-1} n \sigma^2_D = O_P(1)$. In practice, it may be hard to choose $M$ that satisfies this condition since the rate at which $\sigma^2_D \xrightarrow{P} 0$ is in general unknown. I show that if $M$ grows too fast, i.e., if $M^{-1} n \sigma^2_D \xrightarrow{P} 0$, the distribution in (ref) collapses and the guarantees in (ref) hold conservatively. This gives a safe guideline for empirical implementation: choose $M$ to be at least a small fraction of $n$, such as $M= 0.1 n$, and the guarantee in (ref) will hold conservatively.

Application 1: Poverty Prediction in Ghana

Understanding the drivers of poverty is at the root of much of Development Economics. For research, being able to better predict poverty dynamics is of first-order importance to both form hypotheses and then validate theories that explain poverty and poverty dynamics. For policy, accurate predictions of current or future poverty could enable better targeting of interventions (ideally then combined with causal inference on policies and interventions).

Using a sample of 319 households in urban Accra from the ISSER-Northwestern-Yale Long Term Ghana Socioeconomic Panel Survey (GSPS) ghanapaneldataset, I examine how well I can predict which households will be below the poverty line 13 years ahead. The outcome of interest is an indicator for whether a household is below the poverty line in the fourth wave of GSPS (2022/2023), and I use covariates measured in wave 1 (2009/2010), that is, 13 years before. Of the 319 households, 22 were below the poverty line in wave 4 (around 7%). I use predictive covariates including household demographics, parental education, religion, political and traditional leadership experience, asset holdings, and financial indicators (see (ref) for details). Although I focus on the binary indicator of below the poverty line, the approach applies more broadly and could use other outcomes such as level of consumption or assets.

I estimate two quantities: the mean squared error (MSE) and the fraction in poverty by tercile of predicted probability of being below the poverty line. In both cases, I use repeated cross-fitting with $K=3$ and $M=200$, and fit random forest models using the R package ranger implemented through mlr3. Let $i \in \left\{ 1,\dots,319 \right\}$, $Y_i$ denote the indicator of whether household $i$ is below the poverty line in wave 4 of GSPS and $X_i$ the set of covariates measured in wave 1. The estimated MSE is given by $$\hat{\theta}_{{\hat{\eta}}, {\rm MSE}} = \frac{1}{M} \sum_{r \in \mathcal{R}} \frac{1}{n} \sum_{i=1}^{n} \left( Y_i - {\hat{\eta}}_{\tilde{\mathsf{s}}}(X_i) \right)^2.$$ For $j \in \left\{ 1,2 \right\}$, let $\hat{t}_{j,\tilde{\mathsf{s}}}$ be the first and second terciles of $({\hat{\eta}}_{\tilde{\mathsf{s}}}(X_i))_{i=1}^n$, that is, $$\hat{t}_{j,\tilde{\mathsf{s}}} = \inf\left\{t : \frac{1}{|\mathsf{s}|} \sum_{i \in \mathsf{s}} \mathbf{1}\{{\hat{\eta}}_{\tilde{\mathsf{s}}}(X_i) \le t\} \ge \frac{j}{3}\right\},$$ and let $\hat{t}_{0,\tilde{\mathsf{s}}} = -\infty$, $\hat{t}_{3,\tilde{\mathsf{s}}} = \infty$. For $j \in \left\{ 1,2,3 \right\}$, the fraction in poverty in tercile $j$ of predicted probability of being below the poverty line is given by $$\hat{\theta}_{{\hat{\eta}}, {\rm Frac} j} = \frac{1}{M K} \sum_{r \in \mathcal{R}} \sum_{\mathsf{s} \in r} \frac{\sum_{i \in \mathsf{s}} Y_i \mathbb{I}\left(\hat{t}_{j - 1,\tilde{\mathsf{s}}} < {\hat{\eta}}_{\tilde{\mathsf{s}}}(X_i) \le \hat{t}_{j,\tilde{\mathsf{s}}}\right)}{\sum_{i \in \mathsf{s}} \mathbb{I}\left(\hat{t}_{j - 1,\tilde{\mathsf{s}}} < {\hat{\eta}}_{\tilde{\mathsf{s}}}(X_i) \le \hat{t}_{j,\tilde{\mathsf{s}}}\right)}.$$ I show in (ref) that $\hat{\theta}_{{\hat{\eta}}, {\rm Frac} j}$ is a Z-estimator.

I also compare the MSE of the models estimated with random forests to the MSE of using the sample average, as described in (ref). In particular, I report p-values for the test of (ref).

I calculate the MSE estimators and the one-sided test both in the real data and in two Monte Carlo designs, described in (ref). The data generating processes are designed to be similar to the original dataset, preserving the empirical marginals and rank-based dependence structure of the observed data. In the first design, denoted Correlated, the outcome $Y$ is correlated to the covariates $X$. In the second design, denoted Uncorrelated, the outcome is independent of the covariates. I run around 5,000 Monte Carlo iterations for each of the three designs -- real data, “correlated” and “uncorrelated” simulated data --, drawing 200 new random splits of the sample at each Monte Carlo iteration. For the real data, the only source of randomness are the 200 splits, while for the simulation designs I draw a new dataset at each iteration (with 200 splits for each dataset). For each simulated dataset and split, I also calculate the difference between top and bottom terciles, $\hat{\theta}_{{\hat{\eta}}, {\rm Frac} 3} - \hat{\theta}_{{\hat{\eta}}, {\rm Frac} 1}$.

I compare the estimates and p-values of using repeated cross-fitting (RCF) with three alternatives. The first is the standard “twice the median” (TTM) rule ruger1978maximale,gasparin2025combining,chernozhukov2025generic: calculate the p-value (for difference in MSE or “top minus bottom” estimator) separately for each fold, that is, using a third of the data, take the median of the 600 p-values (200 repetitions, 3 folds) and multiply it by 2. The second is the Sequential Aggregation (Seq) approach of luedtke2016statistical and wager2024sequential: train a random forest using only fold 1, compute the t-statistic using fold 2, then train a random forest using folds 1 and 2 and compute the t-statistic in fold 3. The p-value for each repetition of cross-fitting uses as final t-statistic $\sqrt{2}$ times the average of the two t-statistics. Finally, the final p-value for each Monte Carlo iteration is twice the median over the 200 p-values coming from the 200 repetitions, similar to chernozhukov2025reply. The third method is standard sample-splitting (SS): train a random forest using two thirds of the data, calculate p-value in the excluded third, not aggregating across repetitions.

(ref) presents the p-values for whether random forest MSE is lower than sample average MSE, and accuracy ($1 - MSE$) point estimates across Monte Carlo iterations for the two simulation designs as well as for the real data. In the uncorrelated design, all methods exhibit similar accuracy on average, with sample-splitting having larger variance since it does not aggregate across multiple splits. All methods are conservative: the p-values concentrate around 1. For sample-splitting, this happens since the sample average is the best predictor of $Y$ in this design, and the random forests are noisy estimates that have larger MSE. The other methods are conservative for the same reason, and TTM and Seq are more conservative since they take twice the median p-value, which guards against the worst DGP. For the correlated design, all methods correctly give small p-values, with RCF being more concentrated around zero.

In the real dataset, Seq often has the smallest accuracy, and TTM the highest, while RCF stands in between. Seq has smaller accuracy since one the two models that it averages over is trained with only a third of the data. RCF and TTM, on the other hand, always use two thirds of the data for training. The only difference between the two numbers is that RCF averages model performances over 200 repetitions while TTM takes the median. Hence, the higher accuracy of TTM reflects the distribution of model accuracies being left-skewed. Only RCF manages to consistently reject the null, concluding that poverty can be predicted from the observed covariates using a random forest model. TTM and Seq are more conservative, with Seq having larger p-values than TTM due to its lower accuracy.

A comparison similar to (ref) for the top minus bottom estimator $\hat{\theta}_{{\hat{\eta}}, {\rm Frac} 3} - \hat{\theta}_{{\hat{\eta}}, {\rm Frac} 1}$ is presented in (ref).

(ref) shows the point estimates and CIs for the estimators $\hat{\theta}_{{\hat{\eta}}, {\rm Frac} 1}$, $\hat{\theta}_{{\hat{\eta}}, {\rm Frac} 3}$, and their difference using RCF and TTM in the real dataset, as well as p-values for testing whether the difference between top and bottom groups is positive (that is, top tercile has a larger fraction below the poverty line than the bottom tercile). These final estimates aggregate over all the 7,104,600 splits displayed in (ref), averaging for RCF and taking the median for TTM. I do not display the point estimates for Seq since wager2024sequential focuses on testing, but the p-value indicates that the difference between top and bottom groups is not significant. (ref) shows that the difference between top and bottom terciles is statistically significant only for RCF.

Application 2: Heterogeneous Treatment Effects in Charitable Giving

There has been growing interest in the literature for learning features of heterogeneous treatment effects using machine learning (chernozhukov2025generic,wager2024sequential,imai2025statistical; for applications, see, e.g., bryan2024big,athey2025machine,johnson2023improving). I revisit the Generic Machine Learning framework of chernozhukov2025generic (henceforth CDDF), and propose a new ensemble estimator that uses the entire sample for calculating confidence intervals, more data for training machine learning algorithms, and aggregates predictions over multiple ML predictors into an ensemble. I first revisit CDDF's approach, and second introduce my ensemble estimator. Theoretical properties are delayed to (ref). Finally, I compare my estimator to the approaches of CDDF and of wager2024sequential in a Monte Carlo design and in an empirical application using data from karlan2007does. The simulation exercise shows gains in power using the ensemble method, and the ensemble approach is the only to detect statistically significant treatment effect heterogeneity in the empirical application.

The Generic ML Approach of \texorpdfstring{chernozhukov2025generic}{Chernozhukov et al.\ (2025b)}

CDDF proposed a method for learning features of treatment effect heterogeneity in randomized trials. In this section, I focus on their Sorted Group Average Treatment Effects (GATES) estimand. This approach consists of using a machine learning (ML) algorithm and pre-treatment covariates to find groups of individuals with larger and smaller average treatment effects (ATEs). If such groups exist, this means that treatment effect is heterogeneous and that this heterogeneity can be explained at least in part by observable characteristics. Moreover, one can explore how these groups differ in terms of these characteristics. They call this last step Classification Analysis (CLAN), and although I focus on GATES to simplify exposition, my results also hold for CLAN.

First, I define some notation. Let $D = \left( Y_i, T_i, X_i \right)_{i=1}^n$ denote the data, where $Y$ is a scalar outcome, $T$ is the treatment assignment indicator, and $X$ is a vector of pre-treatment covariates. I assume that $(Y_i, T_i, X_i)$ are drawn i.i.d. from a distribution $P \in \mathcal{P}$. Let $\mathcal{A}$ denote an ML algorithm, a function that takes a dataset as input, and outputs an estimate of the Conditional Average Treatment Effect (CATE) function, $$\eta_P(x) = \E[P]{Y(1) - Y(0) | X = x}.$$ For example, $\mathcal{A}$ could be Causal Forests wager2018estimation, or based on Random Forests, Neural Networks, or Gradient Boosting.\footnote{For example, one could use any of these three algorithms to estimate separately the functions $\E[P]{Y(1) | X = x}$ and $\E[P]{Y(0) | X = x}$, and use the difference of the two estimated functions as an estimate of the CATE.} For any subsample $\mathsf{s} \subseteq \left\{ 1,\dots,n \right\}$, let $D_{\mathsf{s}} = \left\{ Y_i, T_i, X_i \right\}_{i \in s}$, $\tilde{\mathsf{s}} = \left\{ 1,\dots,n \right\} \setminus \mathsf{s}$, and ${\hat{\eta}}_{\tilde{\mathsf{s}}} = \mathcal{A}(D_{\tilde{\mathsf{s}}})$, that is, ${\hat{\eta}}_{\tilde{\mathsf{s}}}$ is the model trained with algorithm $\mathcal{A}$ using the subsample $D_{\tilde{\mathsf{s}}}$.

The procedure is given as follows. First, take $M$ random subsets of $\left\{ 1,\dots,n \right\}$ of size $\pi n$. For each $m=1,\dots,M$, denote the subsample by $\mathsf{s}_m$, where $\mathsf{s}_m \subseteq \left\{ 1,\dots,n \right\}$ and $\left| \mathsf{s}_m \right| = \pi n$. For each repetition $m$, call $\mathsf{s}_m$ the main sample, and $\tilde{\mathsf{s}}_m = \left\{ 1,\dots,n \right\} \setminus \mathsf{s}_m$ the auxiliary sample. For $m=1,\dots,M$, train the model

using data from the auxiliary sample. In the main sample, calculate predicted individual treatment effects (ITEs) $\hat{\tau}_i = {\hat{\eta}}_{\tilde{\mathsf{s}}_m}(X_i)$. Sort $\left( \hat{\tau}_{i} \right)_{i \in \mathsf{s}}$ into $J$ quantile groups $G_1, \dots, G_J$, where

with $I_j = [\hat{d}_{j-1}, \hat{d}_{j})$, $-\infty = \hat{d}_0 < \hat{d}_1 < \dots < \hat{d}_J = \infty$, and $(\hat{d}_j)_{j=0}^J$ are calculated such that the number of observations in $(G_j)_{j=1}^J$ is balanced or nearly balanced. For example, with $J=4$, $(G_j)_{j=1}^J$ is a partition of the sample into quartiles of $\left( \hat{\tau}_{i} \right)_{i \in s}$. Calculate the split-specific GATES estimator by running the weighted regression

with weights $\omega_i = \left\{ p(X_i) \left[ 1 - p(X_i) \right] \right\}^{-1}$, where $p(x) = P(T = 1 | X = x)$ is the (known) propensity score. These weights guarantee correct identification of ATEs when the propensity score is not constant, that is, it ensures $$\gamma_{j}^{(m)} = \E[P]{Y_i(1) - Y_i(0) | i \in G_j}.$$ Denote the estimates by $(\hat{\gamma}_{j}^{(m)})_{j=1}^J$. A frequent parameter of interest is $$\delta^{(m)} = \gamma_{J}^{(m)} - \gamma_{1}^{(m)},$$ the difference in ATEs between the top and bottom groups of predicted ITEs. This parameter can be estimated with the analogue $$\hat{\delta}^{(m)} = \hat{\gamma}_{J}^{(m)} - \hat{\gamma}_{1}^{(m)},$$ and a CI can be calculated as usual,

where $\hat{\sigma}^{(m)} / \sqrt{\pi n}$ is a heteroscedasticity-robust standard error for $\hat{\delta}^{(m)}$ calculated as usual from the OLS regression (ref), and $z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of the standard normal distribution. Finally, the final estimators and CIs are given by $$\hat{\delta} = {\rm Med}(\hat{\delta}^{(m)})$$ and $$(L, U) = \left( {\rm Med}(L^{(m)}), {\rm Med}(U^{(m)}) \right),$$ where ${\rm Med}$ denotes the median across repetitions $m$. Conditions for the validity of this CI are established in Theorem 4.3 of CDDF.

This approach carries a tradeoff that's not present in my method, and it considers a single ML algorithm $\mathcal{A}$. The tradeoff regards the choice of $\pi$: a larger $\pi$ means more data is used to estimate the regression (ref), leading to narrower CIs in (ref); but fewer data are used to train the ML model in (ref), likely yielding a worse estimate of the CATE. Moreover, regularity condition R3 in CDDF requires $\pi$ to be relatively small to guarantee that the CI $[L, U]$ covers the median of $\delta^{(m)}$ across all possible splits. My ensemble approach presented next avoids this tradeoff since it uses the entire sample for estimation and a larger sample for training. The ensemble estimator also incorporates more than one ML algorithm, which is important if one does not want to commit beforehand to any specific algorithm. Although CDDF's approach can be repeated with different algorithms, that comes with potential issues of multiple hypothesis testing.

In the next subsection I propose a new GATES estimator that (i) uses the entire sample to calculate $(\hat{\gamma}_{j})_{j=1}^J$ in (ref), and (ii) combines predictions from multiple ML algorithms to form an ensemble, eliminating the need for algorithm selection.

An Ensemble Estimator

Before defining my ensemble estimator, I introduce some additional notation. Theoretical properties are delayed to (ref). Let $A$ denote the number of machine learning algorithms that will be used for predicting ITEs. For $a=1,\dots,A$, let $\mathcal{A}_a$ denote an ML algorithm, that is, a function that takes a dataset as input, and outputs an estimate of the CATE. For example, one could choose $\mathcal{A}_1$ to use Random Forests, $\mathcal{A}_2$ Neural Nets, and $\mathcal{A}_3$ Gradient Boosting. For $\mathsf{s} \subseteq \left[ 1,\dots,n \right]$ and $a=1,\dots,A$, let $${\hat{\eta}}_{\mathsf{s},a} = \mathcal{A}_a(D_{\mathsf{s}}),$$ that is, ${\hat{\eta}}_{\mathsf{s},a}$ is the model trained with algorithm $\mathcal{A}_a$ using the subsample $D_\mathsf{s}$.

The ensemble approach is summarized in (ref). The first difference is that instead of splitting the sample into two sets, I split it into $K$ roughly equal-sized folds $(\mathsf{s}_k)_{k=1}^K$, again repeating the process $M$ times. I calculate $A$ predicted ITEs for each individual using the $A$ ML algorithms, trained using all folds except the one that contains observation $i$. I denote the predicted ITEs by $\hat{\tau}_{i,a} = {\hat{\eta}}_{\tilde{\mathsf{s}}_{k(i)}}(X_i)$, where $k(i)$ is such that $i \in \mathsf{s}_{k(i)}$. Then, to calibrate the weights for combining the multiple ML predictions into one, I split the sample again into $L$ different folds, for each repetition $m=1,\dots,M$. Let $\{\mathsf{s}'_{\ell}\}_{\ell=1}^L$ denote the $L$ folds ($m$ is not incorporated in the notation to simplify exposition). For $\ell=1,\dots,L$, estimate the weighted regression

with weights $\omega_i = \left\{ p(X_i) \left[ 1 - p(X_i) \right] \right\}^{-1}$. In (ref), $\bar{\tau}_{a} = \frac{1}{n - \left| \mathsf{s}'_{\ell} \right|} \sum_{i \not\in \mathsf{s}'_{\ell}} \hat{\tau}_{i,a}$, $p(X_i)$ is the propensity score, and $Z_i$ is a vector of functions of $X_i$, for example $Z_i = (X_{1,i}, p(X_i))'$, where $X_{1,i}$ is a subset of $X_i$. The role of $Z_i$ is only reducing noise in estimation, so this term can be omitted if desired. Denote the estimates of $(\beta_{\ell, a})_{a=1}^A$ by $({\hat{\beta}}_{\ell, a})_{a=1}^A$. The final predicted ITE is then given by $$\hat{\tau}_{i} = \sum_{a = 1}^A {\hat{\beta}}_{\ell,a} \hat{\tau}_{i,a}, i \in \mathsf{s}'_\ell.$$ Repeating this process for $\ell=1,\dots,L$ gives $\hat{\tau}_{i}$ for every observation. I sort $\left( \hat{\tau}_{i} \right)_{i \in \mathsf{s}}$ into groups separately by fold. That is, for $k=1,\dots,K$,

with $I_{j,k} = [\hat{d}_{j-1,k}, \hat{d}_{j,k})$, $-\infty = \hat{d}_{0,k} < \hat{d}_{1,k} < \dots < \hat{d}_{J,k} = \infty$, and $(\hat{d}_{j,k})_{j=0}^J$ are calculated such that the number of observations in $(G_{j,k})_{j=1}^J$ is balanced or nearly balanced. Finally, the split-specific GATES estimator uses the whole sample, defining

and running the weighted regression

with weights $\omega_i = \left\{ p(X_i) \left[ 1 - p(X_i) \right] \right\}^{-1}$.

(ref) is very close to the Best Linear Predictor (BLP) regression of CDDF, except that it uses the $A$ predicted ITEs instead of just one. The intuition behind (ref) is that $(\beta_a)_{a=1}^A$ are the best linear predictor coefficients of a regression where the true CATE $\eta_P(X_i)$ is the response variable, and $(\hat{\tau}_{i,a})_{a=1}^A$ are the independent variables (see Theorem 3.1 of CDDF). Hence, $\sum_{a = 1}^A \beta_a \hat{\tau}_{i,a}$ is the best linear approximation of $\eta_P(X_i)$ given $(\hat{\tau}_{i,a})_{a=1}^A$.

The final estimator averages over repetitions, $$\hat{\delta}_{\hat{\eta}} = \hat{\delta} = \frac{1}{M} \sum_{m=1}^M \hat{\delta}^{(m)},$$ where, as before, $\hat{\delta}^{(m)} = \hat{\gamma}_{J}^{(m)} - \hat{\gamma}_{1}^{(m)}$, with $\hat{\gamma}_{J}^{(m)}$ and $\hat{\gamma}_{1}^{(m)}$ being the estimates from (ref). The final standard error is

where $\hat{\sigma}^{(m)} / \sqrt{n}$ is a heteroscedasticity-robust standard error for $\hat{\delta}^{(m)}$ calculated as usual from the OLS regression (ref). The parameter of interest is $$\delta_{\hat{\eta}} = \delta = \frac{1}{M} \sum_{m=1}^M \gamma_{J}^{(m)} - \gamma_{1}^{(m)},$$ where $\gamma_{J}^{(m)}$ and $\gamma_{1}^{(m)}$ are defined in (ref).

Application to Charitable Giving and Monte Carlo Experiments

I compare my new ensemble approach to two alternative methods in an empirical application and in Monte Carlo experiments. I revisit karlan2007does, which sent fundraising letters to prior donors of a liberal nonprofit organization in the United States, randomizing the match ratio offered (1:1, 2:1, or 3:1) versus no match for a control group. I pool all match treatments into a single treatment group, focusing on the binary treatment of receiving any match offer versus none. The outcome of interest is the amount donated in dollars. The predictive covariates I use include individual donation history (frequency, recency, amount), gender, state-level political variables (Bush vote share, count of court cases in which the organization was either a party to or filed a brief), and zip code-level demographics and economics (race, age, household size, income, homeownership, education, urbanization) (see (ref) for details). I focus on the subset of 6,419 donors who donated within the last two months, as they were more responsive to the solicitation and the smaller sample facilitates computation of the Monte Carlo experiments.

I compare the ensemble with CDDF's approach, described in (ref), and the sequential aggregation approach of luedtke2016statistical, wager2024sequential, and chernozhukov2025reply. Sequential aggregation (Seq) consists of splitting the sample into $K$ folds, for $k=2,\dots,K$ train an ML model using folds $1$ through $k-1$, and compute GATES in the $K$-th fold. The final estimator is the average over the $K-1$ estimates, and the p-value uses the final t-statistic equal to $\sqrt{K-1}$ times the average of the fold-specific t-statistics. This approach uses more data for calculating GATES and p-values ($n (K-1)/K$ observations), but trains some ML models using fewer data (the first model uses $n/K$ observations). I aggregate the final estimates and p-values taking the median over $M$ repetitions as in chernozhukov2025reply.

I compute the three approaches across four designs: (i) using the real data (real), (ii) using the real data but shuffling the treatment assignment indicator at random (so there is no treatment effect heterogeneity) (real-shuffled), (iii) drawing from a DGP where treatment effect is partially predictable using covariates (mc-hte), (iv) drawing from a DGP where treatment effect heterogeneity is independent of covariates (mc-nohte). The two DGPs are meant to be similar to the real data, preserving the marginal distributions of covariates and rank-correlation structure, as described in (ref). Across all methods and datasets, at each Monte Carlo iteration I use 100 repetitions of sample-splitting, take random samples (without replacement) of sizes $n=500,1000,2000,6419$ (entire dataset), and compare the number of folds $K=2,3,5,10$ (for CDDF, the ML is trained with $n(K-1)/K$ observations and GATES calculated in the remaining sample). For Ensemble, I draw at random between 1 and 4 ML algorithms among 10 popular algorithms available in R's mlr3verse: XGBoost (xgboost), Random Forest (ranger), Neural Networks (nnet), Elastic Net (glmnet), k-Nearest Neighbors (\texttt{kknn}), Linear Regression (\texttt{lm}), Decision Trees (\texttt{rpart}), Fast Nearest Neighbors (\texttt{fnn}), Multivariate Adaptive Regression Splines (\texttt{earth}), and Gradient Boosting (\texttt{gbm}). For CDDF and Seq, I draw one of the same ten algorithms at random, for each Monte Carlo iteration. I show the number of iterations used for each specification in (ref) in the appendix.

(ref) shows the gains in power of using the ensemble method in the real dataset. It displays boxplots of one-sided p-values for testing whether the top tercile of predicted treatment effects has a larger ATE than the bottom tercile. A small p-value means rejecting the null hypothesis of no detectable treatment effect heterogeneity. With $n=6419$ (the entire dataset), Ensemble with 4 algorithms detects treatment effect heterogeneity at the 10% level in more than 75% of the iterations. Seq and CDDF give p-values above 10% in most iterations. None of the methods are powered enough to reject the null consistently with $n = 2000$.

(ref) is similar to (ref), except that it uses the synthetic DGP where there is no detectable heterogeneity. It shows that all methods correctly fail to reject the null in most iterations. Similar figures for designs real-shuffled and mc-hte are presented in (ref).

(ref) shows the rejection probabilities at the 5% significance level, that is, the percentage of iterations with p-value below 5%. For the two datasets with no detectable heterogeneity, real-shuffled and mc-nohte, all methods are conservative when $K=2$ or $K=3$, they yield rejection probabilities below the nominal level. In the real-shuffle design with $n=6419$ and $K=5$ or $K=10$, the ensemble methods reject the null with probability slightly higher than nominal, but smaller than 10%. With $n=2000$, only Ensemble 4 rejects the null with probability higher than nominal with $K \ge 5$ in the real-shuffled design. In the real dataset, CDDF almost never detects HTE, and Seq detects in less than 20% of iterations with $K=10$ and $n=6419$. The ensemble methods have higher power especially in the specifications using the entire dataset. For example, Ensemble 2 detects heterogeneity in around 50% of iterations with $K=3$ folds. In the synthetic dataset where there is detectable heterogeneity, mc-hte, as well as in the real data, Ensemble 2 and 4 have higher power across all specifications.

As I discuss in (ref), the rejection probability under the null of no detectable heterogeneity could in principle be above the nominal level when using the normal approximation CI. In (ref), I propose an alternative CI that controls size under the null, at the expense of being more conservative and requiring more computational time. However, I note that extensive simulation experiments, including but not limited to the design of (ref), suggest that Ensemble 4 is conservative for relatively small values of $K$. Hence, my recommendation for empirical practice is to use the normal approximation CI with Ensemble 4 and $K=3$.

Conclusion

As predictive algorithms become increasingly popular, using the same dataset to both train and test a new model has become routine across research, policy, and industry. I derived a new inference approach on model properties that averages across several splits of the sample, where at each split one part is used to train a model and the remaining to evaluate it. Compared to a standard 50-50 sample-splitting, my approach improves statistical and modeling power by using more data for training and evaluating, and improves reproducibility, so two researchers using different splits are more likely to reach the same conclusion about statistical significance. Although the practice of averaging over multiple splits is not new, the confidence intervals and establishing their validity appears to be new.

I addressed the main technical challenge, the dependence created by reusing observations across splits, by proving a central limit theorem for the large class of split-sample Z-estimators. Leveraging the data-dependent parameter of interest, my CLT does not require restricting the complexity of the model or its convergence rate, unlike in the classic semiparametrics problem that used cross-fitting and focused on a different parameter that is not data-dependent. This generality is important as it allows the model to be learned with potentially complex machine learning algorithms, as is commonly done across research, policy, and industry.

Using the CLT, I constructed CIs based on the normal approximation that are valid in a large class of problems, and documented cases where this approximation may fail to cover the parameter of interest at nominal rate. I provided a new approach to inference for such problems, focusing on the particular case of inference when comparing the performance between two models. The approach builds on my CLT, and I discussed how the arguments can be extended to other problems. I also provided a general approach that allows the moment functions to have zero limit variance in (ref), by exploring the faster-than-$\sqrt{n}$ convergence of the empirical moment equations and a tuning parameter.

In (ref), I derived a new reproducibility measure for p-values calculated with split-sample Z-estimators. This measure is especially useful when computational resources are limited, quantifying whether a given number of split-sample repetitions suffices for two researchers using different splits to reach similar conclusions about statistical significance with high probability.

Finally, I illustrated the empirical implications of my results by revisiting two important problems in development and public economics: predicting poverty and learning heterogeneous treatment effects in randomized experiments. Using a panel from Ghana ghanapaneldataset and Monte Carlo experiments, repeated cross-fitting performed better than previous alternatives in detecting predictive power for being below the poverty line 13 years ahead. For the heterogeneous treatment effects application, I developed a new ensemble method that uses the entire sample for evaluation, more data for training, and combines multiple machine learning predictors. I revisited karlan2007does's experiment on charitable giving and conducted Monte Carlo simulations. In both cases, my ensemble method achieved improved power for detecting heterogeneous treatment effects compared to previous alternatives.