Statistical inference is often simplified by sample-splitting. This
simplification comes at the cost of the introduction of randomness that is not
native to the data. We propose a simple procedure for sequentially aggregating
statistics constructed with multiple splits of the same sample. The user
specifies a bound and a nominal error rate. If the procedure is implemented
twice on the same data, the nominal error rate approximates the chance that
the results differ by more than the bound. We provide a non-asymptotic
analysis of the accuracy of the nominal error rate and illustrate the
application of the procedure to several widely applied statistical methods.

Keywords: Sample-splitting, Cross-Validation, Replicability, Exchangeable
Pairs, Stability
MSC 2020: 62E17, 62F40, 62L12

##  1\. Introduction

Sample-splitting is ubiquitous in modern statistical theory. Routine
statistical tasks—model selection, dimension reduction, nuisance parameter
estimation—can be implemented on a randomly selected subsample of a data set,
without contaminating a statistical inference produced on its complement. This
principle underlies the widely applied practices of cross-validation for
predictive risk estimation (Stone,, 1974; Arlot and Celisse,, 2010) and cross-
splitting for adaptive estimation of semiparametric models (Bickel,, 1982;
Klaassen,, 1987; Chernozhukov et al.,, 2018), among many other applications.

For a fixed data set, statistics constructed with sample-splitting are not
deterministic. Two researchers can compute the same statistic on the same data
and obtain different values. If this residual variability is large,
substantive results based on sample-split statistics may not be credible.
Researchers are incentivized to report significant results. If there is scope
to materially alter the statistics that they report through the choice of the
split of their sample, should this choice be left to chance?

In this paper, we propose a method for sequentially aggregating sample-split
statistics over multiple splits of the same sample. Our method is based on the
fixed-length sequential confidence intervals of Anscombe, (1952) and Chow and
Robbins, (1965). The procedure takes as input a bound and an error rate. The
statistic of interest is iteratively computed on different splits of the
sample and aggregated in a running average. The procedure is stopped after an
estimate of the residual variability of the running average falls below a
fixed threshold. If the procedure were run twice, we show that the chance that
the outputs differ by more the bound is well-approximated by the error rate.
That is, by setting the bound and the error rate to be sufficiently small,
sample-split statistics aggregated with the procedure are reproducible.

To fix ideas, consider the setting of Banerjee et al., (2015), who study data
collected from a randomized evaluation of a poverty alleviation program
implemented in several developing countries.111Further details on our
treatment of the Banerjee et al., (2015) data are given in Appendix A. We
focus on their evaluation of the effect of the program on the level of
household consumption in Bangladesh. For each household in their sample, they
observe a vector $D_{i}=(Y_{i},W_{i},X_{i})$, where $Y_{i}$ is a measurement
of consumption three years after the implementation of the program, $W_{i}$ is
an indicator denoting assignment to the program, and $X_{i}$ is a vector of
pretreatment covariates. Let $D=(D_{i})_{i=1}^{n}$ collect the data in their
sample.

The augmented inverse propensity score weighted (AIPW) estimator of Robins et
al., (1994) is a standard approach to estimating average treatment effects in
randomized experiments. This estimator is formed by randomly splitting the
sample into two parts. Let $\mathsf{s}$ and $\tilde{\mathsf{s}}$ denote sets
containing the indices of the households in each part. Using data from the
households $i$ in $\tilde{\mathsf{s}}$, nonparametric estimates of the
propensity score and conditional outcome regression

$\pi(x)=P\\{W_{i}=1\mid
X_{i}=x\\}\quad\text{and}\quad\mu_{w}(x)=\mathbb{E}[Y_{i}\mid
X_{i}=x,W_{i}=w]$

are formed. For example, conditional outcome regressions are often estimated
with the Lasso estimator of Tibshirani, (1996), through

$\hat{\mu}_{w}(X_{i})=\hat{\beta}_{w}(\lambda)^{\top}X_{i},\quad\text{where}\quad\hat{\beta}_{w}(\lambda)=\underset{\beta}{\arg\min}\left\\{\sum_{i\in\tilde{\mathsf{s}}}\mathbb{I}\\{W_{i}=w\\}(Y_{i}-\beta^{\top}X_{i})^{2}+\lambda\|\beta\|_{1}\right\\}$
(1.1)

for some penalization parameter $\lambda$, with $\mathbb{I}\\{\cdot\\}$
denoting the indicator function.222For the purposes of this section, we choose
the penalization parameter with a single implementation of $10$-fold cross-
validation of the mean-squared error. Similarly, we estimate the propensity
score with a linear probability model constructed with $10$-fold cross-
validated lasso. Further details on design of our simulation are given in
Appendix A. Collect these estimates into
$\hat{\eta}(D_{\tilde{\mathsf{s}}})=(\hat{\pi},\hat{\mu}_{w})$, where
$D_{\tilde{\mathsf{s}}}$ collects the data with indices in
$\tilde{\mathsf{s}}$. The AIPW estimator is given by

$\displaystyle T(\mathsf{s},D)$
$\displaystyle=\frac{1}{|\mathsf{s}|}\sum_{i\in\mathsf{s}}\psi(D_{i},\hat{\eta}(D_{\tilde{\mathsf{s}}})),\quad\text{where}$
(1.2) $\displaystyle\psi(D_{i},\hat{\eta}(D_{\tilde{\mathsf{s}}}))$
$\displaystyle=\hat{\mu}_{1}(X_{i})-\hat{\mu}_{0}(X_{i})+\frac{W_{i}(Y_{i}-\hat{\mu}_{1}(X_{i}))}{\hat{\pi}(X_{i})}-\frac{(1-W_{i})(Y_{i}-\hat{\mu}_{0}(X_{i}))}{1-\hat{\pi}(X_{i})}.$

Let $\mathsf{s}$ and $\mathsf{s}^{\prime}$ denote independently drawn subsets
of $[n]=\\{1,\ldots,n\\}$ of size $n/2$. Panel A of Figure 1 compares the
densities of the single-split, double-split, and two-fold cross-split
estimators

$T(\mathsf{s},D),\quad\frac{1}{2}(T(\mathsf{s},D)+T(\mathsf{s}^{\prime},D)),\quad\text{and}\quad\frac{1}{2}(T(\mathsf{s},D)+T(\tilde{\mathsf{s}},D)),$
(1.3)

over repeated samples of $(\mathsf{s},\mathsf{s}^{\prime})$ for the Banerjee
et al., (2015) data.

Figure 1. Residual Randomness

Panel A: One Split, Two Splits | Panel B: $10$-Fold Cross Split
---|---
|
Panel C: $10$-Fold Cross Validation

Notes: Panel A of Figure 1 displays kernel estimates of the density of the
single split, double split, and two-fold cross split estimators defined in
(1.3) over repeated samples of $(\mathsf{s},\mathsf{s}^{\prime})$ for the
Banerjee et al., (2015) data. Panel B displays a discretized heat map of the
joint distribution of the 10-fold cross split estimator $a(\mathsf{r}_{10},D)$
defined in (1.4) and the critical value estimator
$\mathsf{CV}_{\alpha}(\mathsf{r}_{10},D)$ defined in (1.5). The black line
denotes the threshold where the estimator is equal to the estimated critical
value. Panel C displays quantiles of the $10$-fold cross-validated estimate of
the mean-squared error (1.6) of the Lasso regression (1.1) over a grid of
values of $\lambda$. The probabilities that each value of $\lambda$ minimizes
the cross-validated risk estimate are displayed with large grey dots.

The variability in the estimator across sample-splits is substantial and
sufficient to switch the sign. A moderate increase in the number of splits
does not address the problem. The $k$-fold cross-split estimator is given by

$a(\mathsf{r}_{k},D)=\frac{1}{k}\sum_{j=1}^{k}T(\mathsf{s}_{j},D)\leavevmode\nobreak\
,$ (1.4)

where $\mathsf{r}_{k}=(\mathsf{s}_{i})_{i=1}^{k}$ is a collection of $k$
equally sized sets that form a partition of $[n]$.333If $\pi$ and $\mu_{w}$
are estimated sufficiently quickly, this estimator is semiparametrically
efficient for the average effect of the program on household consumption
(Chernozhukov et al.,, 2018; Hahn,, 1998). An asymptotically exact test of the
null hypothesis that the average effect of the program is less than or equal
to zero can be constructed by comparing (1.4) with the critical value

$\mathsf{CV}_{\alpha}(\mathsf{r}_{k},D)=\frac{z_{1-\alpha}}{n}\left(\sum_{j=1}^{k}\sum_{i\in\mathsf{s}_{j}}\left(\psi(Y_{i},\hat{\eta}(D_{\tilde{\mathsf{s_{j}}}}))-a(\mathsf{r}_{k},D)\right)^{2}\right)^{1/2}\leavevmode\nobreak\
,$ (1.5)

where $z_{1-\alpha}$ is the $1-\alpha$ quantile of the standard normal
distribution (Chernozhukov et al.,, 2018). Panel B of Figure 1 displays a heat
map of the joint distribution of the 10-fold cross-split statistic
$a(\mathsf{r}_{10},D)$ and the critical value
$\mathsf{CV}_{0.05}(\mathsf{r}_{10},D)$ over random draws of the partition
$\mathsf{r}_{10}$. The residual variability in the estimator has been greatly
reduced, but is still large relative to estimates of the sampling variability.
The choice of the partition $\mathsf{r}_{10}$ determines whether a test that
the treatment effect of the program is negative is rejected.

The problem of residual variability is even more severe when choosing the
penalization parameter $\lambda$ in the Lasso regression (1.1) with cross-
validation. Here, the sample-split statistic of interest is the out-of-sample
mean-squared error estimate

$T(\mathsf{s},D)=\frac{1}{\sum_{i\in\mathsf{s}}\mathbb{I}\\{W_{i}=w\\}}\sum_{i\in\mathsf{s}}\mathbb{I}\\{W_{i}=w\\}(Y_{i}-\hat{\beta}_{w}(\lambda)^{\top}X_{i})^{2}\leavevmode\nobreak\
.$ (1.6)

Often, the penalization parameter $\lambda$ is chosen as the value that
minimizes the cross-fit statistic (1.4). Panel C of Figure 1 displays the
distribution of $10$-fold cross-validated estimates of the mean-squared error
of the regression (1.1), with $w$ equal to one, over a grid of values of
$\lambda$.444We consider the values of $\lambda$ reported by the “glmnet” R
software package (Friedman et al.,, 2021). The probabilities that each value
of $\lambda$ minimizes the cross-validated risk estimate are displayed with
large grey dots. There is substantial variability in the penalization
parameter chosen by $10$-fold cross-validation.

The primary methodological contribution of this paper is the introduction of a
computationally efficient procedure for aggregating sample-split statistics
across multiple splits to ensure that residual variability is small. The
proposal is based on the sequential fixed-width confidence intervals developed
by Anscombe, (1952) and Chow and Robbins, (1965). In Section 2, we detail the
proposed algorithm and give very general conditions for its validity, in a
particular asymptotic sense.

A user must make several choices when implementing the procedure. These
include specifying a suitable bound and error rate as well as choosing the
size and joint distribution of the random splits. The primary theoretical
contribution of this paper is a non-asymptotic analysis that motivates
approaches for making these choices. In particular, in Section 3, we restrict
attention to linearly separable statistics that are symmetric in each part of
each split. We give two main results. First, we derive concentration
inequalities for averages of independently drawn cross-split statistics of the
form (1.4) around their conditional mean. Our bounds characterize the dramatic
difference between the residual variability of the double split and two-fold
cross split statistics displayed in Panel A of Figure 1. Second, we provide a
Berry-Esseen type bound on the accuracy of the nominal error rate for the
procedure. As the procedure is sequential, this calculation is nonstandard.
This result illustrates a trade-off between the computational efficiency of
the procedure and the accuracy of the nominal error rate.

The main theoretical challenge posed by this analysis is the accommodation of
the dependence between statistics computed on cross-splits of a sample. That
is, conditioned on the data, the summands, over $j$ from $1$ to $k$, of the
aggregate statistic (1.4) will be dependent. We address this through an
application of the method of exchangeable pairs (see e.g., Stein,, 1986;
Ross,, 2011; Chen et al.,, 2011). To construct an appropriate exchangeable
pair for our problem, we develop a novel application of a coupling argument
due to Chatterjee, (2005). See also Paulin et al., (2013, 2016). The various
concentration and moment inequalities we develop are then obtained by applying
results due to Chatterjee, (2005, 2007). These arguments are outlined in
Section 5.

In Section 4, we measure the performance of our procedure in the applications
considered in Figure 1. Section 6 concludes. Details concerning the data and
simulations considered in this paper are given in Appendix A. Unless otherwise
specified, proofs for all results stated in the main text are given in
Appendix B. Proofs supporting Lemmas used in Appendix B are given in Appendix
C. Appendix D gives additional results that will be introduced at appropriate
points throughout the paper.

###  1.1. Related Literature

The procedure studied in this paper is based on the sequential fixed-width
confidence intervals of Anscombe, (1952) and Chow and Robbins, (1965). Our
non-asymptotic analysis is related to the Berry-Esseen bounds for randomly
indexed sums given in Landers and Rogge, (1976, 1988). Applications of these
results to fixed-width confidence intervals have been given in Csenki, (1980),
Mukhopadhyay, (1981), and Callaert and Janssen, (1981). A more recent non-
asymptotic consideration of the application of fixed width confidence
intervals to Monte Carlo techniques is given in Hickernell et al., (2013). Our
more specialized analysis differs from these results by accommodating and
quantifying the dependence induced by cross-splitting.

Our setting is related to a large literature that studies methods for
constructing confidence intervals for cross-validated estimates of
generalization error. Several recent examples include Dietterich, (1998),
Nadeau and Bengio, (1999), Lei, (2020), Bayle et al., (2020), Austern and
Zhou, (2020), and Bates et al., (2023). In contrast, we are interested in the
randomness conditional on the data. Formally, our results are most similar to
the Berry-Esseen bounds given in Austern and Zhou, (2020), who study the
unconditional normal approximation of statistics similar to the those
considered in Section 3. Our bounds apply under strictly weaker conditions and
are substantially simpler.

Some of our arguments build on a literature studying the role of algorithmic
stability in the accuracy of cross-validation (Kale et al.,, 2011; Kumar et
al.,, 2013). These papers are related to a broader literature that derives
generalization bounds for stable algorithms, originating with Bousquet and
Elisseeff, (2002). The stability of several specific estimators is studied in
Elisseeff et al., (2005), Hardt et al., (2016), Celisse and Guedj, (2016),
Chen et al., (2022), and Du et al., (2023). Some of the concentration
inequalities we derive can be compared to the results of Cornec, (2010) and
Abou-Moustafa and Szepesvári, (2019). Again, the setting we study is different
and our conditions are substantially simpler.

Chen et al., (2022) show that sample-splitting is unnecessary for the
asymptotic normality of estimators of the form (1.2) if the nuisance parameter
estimator $\hat{\eta}$ satisfies a stability condition related to some of the
conditions studied in this paper. The intention of our paper is not to argue
that sample-splitting should, necessarily, be used for any particular
application. Rather, our premise is that sample-splitting is widely applied in
practice, irrespective of its optimality.

The procedure studied in this paper is applicable to a large variety of
statistical methods constructed with sample-splitting. General methods for
testing statistical hypothesis with sample-splitting are studied in Guo and
Romano, (2017), DiCiccio and Romano, (2019), DiCiccio et al., (2020),
Wasserman et al., (2020), and Tse and Davison, (2022). Rinaldo et al., (2019)
study a general procedure for selective inference based on sample-splitting.
Meinshausen et al., (2009) and Meinshausen and Bühlmann, (2010) study methods
for inference in high-dimensional linear models that use multiple data
splitting.

Although our emphasis is on statistics constructed with sample-splitting, the
algorithmic and formal methods studied in this paper are potentially
applicable to randomized algorithms more generally. See Beran and Millar,
(1987) for an analysis of the asymptotics of randomized tests and estimators.
Guo and Shah, (2023) and Zhang et al., (2023) propose methods for conducting
inference with randomized algorithms through subsampling (Politis et al.,,
1999).

###  1.2. Notation

We let $\tilde{\mathsf{s}}$ denote the complement of the set $\mathsf{s}$ in
$[n]=\\{1,\ldots,n\\}$. The set $\mathcal{S}_{n,b}$ consists of all subsets of
$\left[n\right]$ of size $b$. The set $\mathcal{R}_{n,k,b}$ contains all
collections of $k$ mutually exclusive elements of $\mathcal{S}_{n,b}$. We
refer to elements of $\mathcal{R}_{n,k,b}$ as cross-splits. That is, if
$n=k\cdot b$, then $\mathcal{R}_{n,k,b}$ is the set of all partitions of $[n]$
into $k$ mutually exclusive subsets of size $b$. Likewise,
$\mathcal{R}_{n,1,b}$ is the collection of sets containing one element of
$\mathcal{S}_{n,b}$. We say $a\lesssim b$ if there exists a universal constant
$C$ such that $a\leq Cb$. We let $z_{\alpha}$ denote the $\alpha$ quantile of
the standard normal distribution and $\lfloor x\rfloor$ denote the largest
integer smaller than or equal to $x$.

##  2\. Reproducible Aggregation

Consider a sample $D=(D_{i})_{i=1}^{n}$. Let $\mathsf{s}$ be some subset of
$\left[n\right]$ with complement $\tilde{\mathsf{s}}$. We are interested in
reporting a real-valued sample-split statistic of the form

$T\left(\mathsf{s},D\right)=\Psi\left(D_{\mathsf{s}},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right)\right),$
(2.1)

where the collection $D_{\mathsf{s}}=(D_{i})_{i\in\mathsf{s}}$ contains data
with indices in $\mathsf{s}$, $\hat{\eta}$ is an estimator of an unknown
nuisance parameter $\eta$, and $\Psi\left(\cdot\right)$ is some function. The
structure encoded in statistics of the form (2.1) is quite general and
encompasses most applications of sample-splitting encountered in practice,
including out-of-sample evaluation of prediction error, model selection, and
nuisance parameter estimation.

We study methods for aggregating statistics of the form (2.1) across multiple
splits of the same sample. Our objective is to ensure that the auxiliary
randomness introduced by sample-splitting is small. In particular, for each
collection $\mathsf{R}_{g,k}=(\mathsf{r}_{i})_{i=1}^{g}$ of $g$ elements of
$\mathcal{R}_{n,k,b}$, where we write
$\mathsf{r}_{i}=(\mathsf{s}_{i,j})_{j=1}^{k}$, we consider aggregate
statistics of the form

$a(\mathsf{R}_{g,k},D)=\frac{1}{g}\frac{1}{k}\sum_{i=1}^{g}\sum_{j=1}^{k}T(\mathsf{s}_{i,j},D)\leavevmode\nobreak\
.$ (2.2)

The statistic (2.2) generalizes several standard methods for aggregating
sample-split statistics. If $k=n/b$ and $g=1$, then (2.2) is a $k$-fold cross-
split statistic (also known as, e.g., $k$-fold cross-validation or $k$-fold
cross-fitting). If $k=1$, then (2.2) averages over $g$ independently drawn
sample-splits. In this case, if $\hat{\eta}$ is constant, then (2.2) is an
incomplete $U$-statistic of order $b$.555A special case of the class of
statistics with $k=1$ are subsampled statistics (Politis et al.,, 1999). We
entertain the intermediate cases, i.e., where $k$ is between $1$ and $n/b$, so
that the theoretical results to follow interpolate between independent
splitting and cross-splitting. Throughout, we denote the proportion $b/n$ by
$\varphi$.

Our task is to formulate a method for choosing the number of collections of
mutually exclusive splits $g$ to ensure that the residual variability of the
aggregate statistic (2.2) is small. This objective is formalized in the
following criterion.

######  Definition 2.1 (Reproducibility).

Suppose that the integers $\hat{g}$ and $\hat{g}^{\prime}$ and the collections
$\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^{\prime}_{\hat{g}^{\prime},k}$, are
independent and identically distributed, conditional on the data $D$. We say
that $\mathsf{R}_{\hat{g},k}$ is $(\xi,\beta)$-reproducible if

$P\bigg{\\{}\big{|}a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}^{\prime}_{\hat{g}^{\prime},k},D)\big{|}\leq\xi\mid
D\bigg{\\}}\geq 1-\beta$ (2.3)

almost surely.

######  Remark 2.1.

Suppose that the data $D$ are given to two researchers. Each researcher is
tasked with producing an estimate of the form (2.2). They generate the
collections of splits $\mathsf{R}_{\hat{g},k}$ and
$\mathsf{R}^{\prime}_{\hat{g}^{\prime},k}$ independently using the same
procedure. That is, the two collections of splits used by the two researchers
are independent and identically distributed conditional on the data $D$. If
this procedure is $(\xi,\beta)$-reproducible, then the probability that the
two researchers’ estimates differ by more than $\xi$ is less than $\beta$.
$\blacksquare$

We propose a sequential method for constructing a reproducible collection
$\mathsf{R}_{\hat{g},k}$. Our proposal is based on the fixed-length sequential
confidence intervals of Anscombe, (1952) and Chow and Robbins, (1965). Define
the variance estimator

$\displaystyle\hat{v}\left(\mathsf{R}_{g,k},D\right)$
$\displaystyle=\frac{1}{g\left(g-1\right)k^{2}}\sum_{i=1}^{g}\sum_{j,j^{\prime}=1}^{k}\left(T(\mathsf{s}_{i,j},D)-a\left(\mathsf{R}_{g,k},D\right)\right)\left(T(\mathsf{s}_{i,j^{\prime}},D)-a\left(\mathsf{R}_{g,k},D\right)\right).$
(2.4)

Observe that (2.4) is simply the sample variance of the summands in (2.2)
across the $g$ collections of splits. Our procedure iteratively draws a
collection uniformly at random and computes the variance (2.4) until the
condition

$\hat{v}\left(\mathsf{R}_{g,k},D\right)\leq\mathsf{cv}_{\xi,\beta}=\frac{1}{2}\left(\frac{\xi}{z_{1-\beta/2}}\right)^{2}$
(2.5)

is satisfied. We let $\hat{g}$ denote the smallest value of $g$ greater than
or equal to $g_{\mathsf{init}}$ such that (2.5) is satisfied, where
$g_{\mathsf{init}}\geq 2$ is an integer chosen by the user. This procedure is
summarized in Algorithm 1.

Input: Data $D$, tolerance $\xi$, error rate $\beta$, collection size $k$,
split size $b$, initialization $g_{\mathsf{init}}$

1

2Set $g\leftarrow g_{\mathsf{init}}$

3Draw $\mathsf{r}_{1}$, …, $\mathsf{r}_{g_{\mathsf{init}}}$ independently and
uniformly from $\mathcal{R}_{n,k,b}$. Collect
$\mathsf{R}_{g_{\mathsf{init}},k}=(\mathsf{r}_{j})_{j=1}^{g_{\mathsf{init}}}$.

4while  _$\hat{v}(\mathsf{R}_{g,k},D) >\mathsf{cv}_{\xi,\beta}$_ do

5

6 Set $g\leftarrow g+1$

7 Draw $\mathsf{r}_{g}$ uniformly from $\mathcal{R}_{n,k,b}$. Collect
$\mathsf{R}_{g,k}=(\mathsf{R}_{g-1},\mathsf{r}_{g})$.

8 end while

9Set $\hat{g}\leftarrow g$

10return $\mathsf{R}_{\hat{g},k}$

Algorithm 1 Anscombe-Chow-Robbins Aggregation

The following theorem gives an asymptotic sense in which the collection of
sample-splits chosen with Algorithm 1 are $(\xi,\beta)$-reproducible. The
proof is closely related to the arguments of Anscombe, (1952) and Chow and
Robbins, (1965). See e.g., Theorem 3.1 of Gut, (2009) for a textbook
treatment.

###### Theorem 2.1.

Suppose that the real-valued statistic $a(\mathsf{R}_{1,k},D)$ has a non-zero
variance conditional on $D$ almost surely. If the collections
$\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^{\prime}_{\hat{g}^{\prime},k}$ are
independently computed with Algorithm 1, then

$P\bigg{\\{}\big{|}a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}^{\prime}_{\hat{g}^{\prime},k},D)\big{|}\leq\xi\mid
D\bigg{\\}}\to 1-\beta\quad\text{as}\quad\xi\to 0\leavevmode\nobreak\ .$ (2.6)

######  Remark 2.2.

Theorem 2.1 applies to a fixed data set $D$. Asymptotics are taken as $\xi\to
0$ with all other quantities, including the sample size $n$, fixed. Roughly,
Theorem 2.1 is established by showing that

$\hat{g}/g^{\star}\overset{\text{a.s}}{\to}1\quad\text{as}\quad\xi\to
0,\quad\text{where}\quad
g^{\star}=\min\left\\{g:\operatorname{Var}\left(a(\mathsf{R}_{g,k},D)\mid
D\right)\leq\mathsf{cv}_{\xi,\beta}\right\\}\leavevmode\nobreak\ ,$ (2.7)

and using the fact that a central limit theorem will hold for the sequence
$a(\mathsf{R}_{g^{\star},k},D)$ as $\xi\to 0$, conditional on the data. The
discrepancy between $a(\mathsf{R}_{\hat{g},k},D)$ and
$a(\mathsf{R}_{g^{\star},k},D)$ can be handled with an appropriate maximal
inequality, an idea due to Rényi, (1957). $\blacksquare$

######  Remark 2.3.

The proof of Theorem 2.1 uses only the fact that the $g$ collections
$\mathsf{r}_{i}$ in $\mathsf{R}_{g,k}$ are drawn independently and
identically. However, Figure 1 suggests that we should expect the statistic
$a(\mathsf{R}_{g,k},D)$ to concentrate as $k$ increases. Some insight into
this phenomenon can be gained without any assumptions. In particular, let
$\mathsf{s}$ and $\mathsf{s}^{\prime}$ each be uniformly distributed on
$\mathcal{S}_{n,b}$ such that their intersection is empty with probability
one. Define the conditional variance and covariance

$\phi_{n,b}(D)=\operatorname{Var}(T(\mathsf{s},D)\mid
D)\quad\text{and}\quad\gamma_{n,b}(D)=\operatorname{Cov}(T(\mathsf{s},D),T(\mathsf{s}^{\prime},D)\mid
D)$ (2.8)

respectively. We can evaluate

$v_{g,k}(D)=\operatorname{Var}(a(\mathsf{R}_{g,k},D)\mid
D)=\frac{1}{g}\left(\frac{1}{k}\phi_{n,b}(D)+\frac{k-1}{k}\gamma_{n,b}(D)\right).$
(2.9)

Observe that $\gamma_{n,b}(D)\leq\phi_{n,b}(D)$ by the Cauchy-Schwarz
inequality. Thus, for a fixed total number of cross-splits $g$ and split size
$b$, the conditional variance $v_{g,k}(D)$ decreases as $k$ increases. That
is, cross-splitting (i.e., setting $k>1$) always has a smaller conditional
variance than independent splitting (i.e., setting $k=1$) for a fixed value of
$g$. On the other hand, for a fixed number of total splits $m=gk$ and split
size $b$, the conditional variance $v_{g,k}(D)$ is smaller for cross-splitting
than for independent splitting if and only if the conditional covariance
$\gamma_{n,b}(D)$ is negative.

In many standard applications, the conditional covariance $\gamma_{n,b}(D)$ is
negative.666This will not always be the case. In Section D.1, we give examples
where $\gamma_{n,b}(D)$ obtains both endpoints of the Cauchy-Schwarz bound
$-\phi_{n,b}(D)\leq\gamma_{n,b}(D)\leq\phi_{n,b}(D)$. Consider again the
Banerjee et al., (2015) data. Figure 2 compares the variance components

$\frac{1}{k}\phi_{n,b}(D)\quad\text{and}\quad-\frac{k-1}{k}\gamma_{n,b}(D)$
(2.10)

under the restriction that $n=kb$ for both the treatment effect and risk
estimation applications introduced in Section 1. In both cases, the two terms
(2.10) are nearly equal, i.e., the two terms in the variance (2.9) nearly
cancel each other out. In other words, cross-splitting is able to account for
almost all of the residual variability introduced by independent splitting. By
placing restrictions on the statistic (2.1), the subsequent section provides a
non-asymptotic account of the error in the approximation (2.6) that quantifies
and compares the concentration in $a(\mathsf{R}_{g,k},D)$ with both $g$ and
$k$. $\blacksquare$

Figure 2. Covariance Across Cross-Splits

Panel A: Treatment Effect Estimation | Panel B: Risk Estimation
---|---
|

Notes: Figure 2 displays measurements of the quantities (2.10) for the
treatment effect and risk estimation applications to the Banerjee et al.,
(2015) data introduced in Section 1. For the risk estimation application, we
study the regression (1.1) for the treated outcome, i.e., $w$ equal to one,
for a fixed choice of $\lambda$. In each case, the parameters $k$ and $b$ are
chosen to satisfy $n=kb$. The $x$-axes are displayed in a log-scale base 10.
Further details on this figure are given in Appendix A.

##  3\. Non-asymptotic Theory

Theorem 2.1 is very general. However, interpretation of the asymptotic
approximation with $\xi\to 0$ is somewhat nebulous, or at least unfamiliar.
What would be a reasonable value of $\xi$ to choose such that the
approximation (2.6) is accurate? Or, how quickly should we expect the error in
(2.6) to diminish as $\xi$ shrinks?

We give a non-asymptotic analysis of the quality of the nominal error rate
$\beta$ in Algorithm 1. This is accomplished by placing restrictions on the
statistics of interest (2.1). We emphasize that Theorem 2.1 should provide
assurance that Algorithm 1 can be applied widely. The more specialized
analysis that follows is aimed at (i) articulating why cross-split statistics
often exhibit significantly lower residual variability than independently
split statistics and (ii) providing intuition for how the error in the
approximation (2.6) depends on the choices of $k$ and $\xi$. Proofs for
results stated in this section are deferred to Section 5.

###  3.1. Symmetry, Linearity, and Stability

We impose two simplifying restrictions on the statistic of interest. First, we
assume that the statistic is symmetric and deterministic in each part of the
split sample.

###### Assumption 3.1 (Symmetry and Determinism).

For all sets $\mathsf{s}\subseteq[n]$ and data $D$, the statistic
$T\left(\mathsf{s},D\right)$ is deterministic and invariant to permutations of
the data with indices in $\mathsf{s}$ and of the data with indices in
$\tilde{\mathsf{s}}$.

######  Remark 3.1.

The intention of Assumption 3.1 is to restrict the residual randomness under
consideration to the randomness introduced by sample-splitting. Assumption 3.1
rules out procedures where the estimator $\hat{\eta}$ is random conditional on
the data. This includes settings where, e.g., $\hat{\eta}$ is estimated with
stochastic gradient descent, bagging or subsampling, or is itself constructed
with data splitting. It does not, however, exclude cases where $\hat{\eta}$ is
deterministic, e.g., when $\hat{\eta}$ is a coefficient vector determined by a
regularized regression or in the applications to hypothesis testing considered
by DiCiccio et al., (2020) or Wasserman et al., (2020). $\blacksquare$

Second, we assume that the statistic (2.1) is linearly separable in the first
part of the split sample.

###### Assumption 3.2 (Linearity).

For all sets $\mathsf{s}\subseteq[n]$ and data $D$, the statistic
$T\left(\mathsf{s},D\right)$ can be written

$T\left(\mathsf{s},D\right)=\Psi\left(D_{\mathsf{s}},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right)\right)=\frac{1}{|\mathsf{s}|}\sum_{i\in\mathsf{s}}\psi(D_{i},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right))$
(3.1)

for some function $\psi(\cdot,\cdot)$.

######  Remark 3.2.

Assumption 3.2 is satisfied in the applications to treatment effect and risk
estimation studied in Section 1 in addition to the applications to hypothesis
testing considered in DiCiccio et al., (2020) and Wasserman et al., (2020).
Extending our analysis to cases where (2.1) satisfies a component-wise
Lipschitz or bounded differences condition is worth further consideration.
$\blacksquare$

Both the remaining assumptions, and the ensuing bounds, specified in this
section are expressed in terms of two objects that measure the sensitivity of
the statistics of interest (2.1) to perturbations of the data and of the
splits, respectively. We refer to these objects as stabilities. Let
$D^{\prime}$ denote an independent and identical copy of the data $D$. For
each $\mathsf{q}\subseteq[n]$, let $\tilde{D}^{(\mathsf{q})}$ be constructed
by replacing $D_{i}$ with $D^{\prime}_{i}$ in $D$ for each $i$ in
$\mathsf{q}$.

######  Definition 3.1 (Sample Stability).

Fix a set $\mathsf{s}\subseteq\mathcal{S}_{n,b}$. Let $I$ be an index drawn
uniformly from $\mathsf{s}$. Let $\mathsf{q}$ be a randomly selected subset of
$\tilde{\mathsf{s}}$ of cardinality $q$. We refer to the quantities

$\displaystyle\sigma_{\mathsf{va\ell id}}^{(r)}$
$\displaystyle=\mathbb{E}\left[\big{|}\psi(D_{I},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right))-\psi(D^{\prime}_{I},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right))\big{|}^{r}\right]\quad\text{and}$
$\displaystyle\sigma_{\mathsf{train}}^{(r,q)}$
$\displaystyle=\mathbb{E}\left[\big{|}\psi(D_{I},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right))-\psi(D_{I},\hat{\eta}(\tilde{D}^{(\mathsf{q})}_{\tilde{\mathsf{s}}}))\big{|}^{r}\right]$

as the $r$th-order validation and $(r,q)$th-order training sample stabilities,
respectively. Similarly, we refer to the quantity

$\sigma^{(r)}_{\mathsf{max}}=\max\left\\{\sigma_{\mathsf{va\ell
id}}^{(r)},\sigma_{\mathsf{train}}^{(r,1)}\right\\}$ (3.2)

as the $r$th-order full sample stability.

######  Definition 3.2 (Split Stability).

We refer to the quantity

$\zeta^{(r)}=\mathbb{E}\left[\max_{\mathsf{s},\mathsf{s}^{\prime}\in\mathsf{S}_{n,b}}\left(T\left(\mathsf{s},D\right)-T\left(\mathsf{s}^{\prime},D\right)\right)^{r}\right].$

as the $r$th-order split stability.

At several points in the analysis that follows, we focus attention on
statistics whose training sample stabilities shrink at an appropriate rate as
$n$ increases.

######  Definition 3.3.

A statistic $T\left(\mathsf{s},D\right)$ that satisfies Assumption 3.2 is
sample stable if

$\sigma_{\mathsf{train}}^{(r,q)}\lesssim\left(\frac{q}{n-b}\right)^{r}$ (3.3)

uniformly for all even values of the parameter $r$.

######  Remark 3.3.

The $(2,1)$th-order training sample stability
$\sigma_{\mathsf{train}}^{(2,1)}$ is a widely studied object and is referred
to as mean-square stability in the statistical learning literature (see e.g.,
Kale et al.,, 2011; Kumar et al.,, 2013).777Connections between sample
stability and generalization are studied in Bousquet and Elisseeff, (2002) and
Celisse and Guedj, (2016). In Section D.2, we give simple sufficient
conditions for the sample stability of statistics that satisfy Assumption 3.2
when $\hat{\eta}$ is an empirical risk minimizer of a potentially regularized
strictly convex loss. The argument is closely related to a similar result
given in Austern and Zhou, (2020). Analogous bounds have been obtained for
many estimators that do not necessarily satisfy Assumption 3.1, such as bagged
or subsampled estimators (Chen et al.,, 2022), ensemble methods (Elisseeff et
al.,, 2005), and estimators computed with stochastic gradient descent (Hardt
et al.,, 2016). It most settings, it is reasonable to expect that the
validation sample stability $\sigma_{\mathsf{va\ell id}}^{(r)}$ is bounded for
small values of $r$. The split stability $\zeta^{(r)}$ is a less frequently
studied object, and we will only require that it is finite for $r$ equal to 4
or 8, depending on the setting. This is a weak restriction that will hold, for
example, if $T\left(\mathsf{s},D\right)$ is bounded. $\blacksquare$

###  3.2. Concentration and Normal Approximation

We begin by studying the concentration of the aggregate statistic
$a(\mathsf{R}_{g,k},D)$ about its conditional mean

$\overline{a}\left(D\right)=\mathbb{E}\left[a(\mathsf{R}_{g,k},D)\mid
D\right]=\mathbb{E}\left[T(\mathsf{s}_{i,j},D)\mid
D\right]\leavevmode\nobreak\ .$ (3.4)

The nonstandard aspect of our analysis is that we account for the dependence
in the summands in $a(\mathsf{R}_{g,k},D)$ across cross-splits, i.e., the
dependence induced by cross-fitting. This is accomplished by constructing a
suitable exchangeable pair with a coupling argument due to Chatterjee, (2005)
and applying a method for deriving concentration inequalities with
exchangeable pairs, also due to Chatterjee, (2005, 2007). Aspects of our
argument are closely related to the analysis of Paulin et al., (2013, 2016).

###### Theorem 3.1.

Suppose that Assumptions 3.1 and 3.2 hold, that the data $D$ are independently
and identically distributed, and that the statistic $a(\mathsf{R}_{1,k},D)$
has a non-zero variance conditional on $D$ almost surely. If the fourth-order
split stability $\zeta^{(4)}$ is finite, then the concentration inequality

$\displaystyle\log
P\left\\{a(\mathsf{R}_{g,k},D)-\overline{a}\left(D\right)\geq t\mid
D\right\\}$ $\displaystyle\leq-\frac{g}{2^{4}(2-\varphi
k-\varphi)^{2}}\frac{\delta t^{2}}{\Gamma_{k,\varphi,b}}$ (3.5)

holds with probability greater than $1-\delta$ as $D$ varies, where

$\displaystyle\Gamma_{k,\varphi,b}$
$\displaystyle=\left(\frac{1-k\varphi}{1-\varphi}\right)4\sigma^{(2)}_{\mathsf{max}}+\left(\frac{k\varphi-\varphi}{1-\varphi}\right)\sigma^{(2,b-1)}_{\mathsf{train}}$
(3.6)

and $\varphi=b/n$.

######  Remark 3.4.

In the inequality (3.5), a conditional probability that depends on the data
$D$ is bounded by a term that depends on the fixed parameters $g$, $k$,
$\varphi$, $\delta$, and $t$ in addition to the unconditional quantity
$\Gamma_{k,\varphi,b}$. In other words, if $\mathcal{F}$ is the event that the
inequality (3.5) holds, then $P\\{\mathcal{F}\\}>1-\delta$ unconditionally.
This transference from a conditional quantity to an unconditional quantity is
obtained by a Markov type bound at one point in the proof. In the sequel, we
give many conditional bounds of a similar form. This strategy—bounding the
conditional quantities of interest with unconditional quantities—is helpful
because the resultant unconditional objects are more tractable and are widely
studied.

######  Remark 3.5.

The quantity $\Gamma_{k,\varphi,b}$ interpolates between
$\sigma^{(2)}_{\mathsf{max}}$ and $\sigma^{(2,b-1)}_{\mathsf{train}}$ as $k$
varies between its bounds $1$ and $1/\varphi$.888The pre-factor $(2-\varphi
k-\varphi)$ also reduces by a factor of 2 as $k$ increases from $1$ to
$1/\varphi$. In the case of independent splitting, the rate of concentration
for the statistic $a(\mathsf{R}_{g,k},D)$ is driven by the full sample
stability $\sigma^{(2)}_{\mathsf{max}}$. On the other hand, in the case of
cross-splitting, concentration depends only on the training sample stability
$\sigma^{(2,b-1)}_{\mathsf{train}}$. At the extreme, if $b$ is set equal to
one and $k$ is set equal to $n$, known as leave-one-out or jackknife cross-
fitting, $\Gamma_{k,\varphi,b}$ will equal zero and there will be no residual
randomness (that is, the assumption that $a(\mathsf{R}_{1,k},D)$ has a non-
zero conditional variance is violated).

To provide a quantitative sense of this tradeoff, Figure 3 displays estimates
of $\sigma^{(2)}_{\mathsf{va\ell id}}$ and $\sigma^{(2,b-1)}_{\mathsf{train}}$
for the risk and treatment effect estimation applications to the Banerjee et
al., (2015) data, where we have approximated the unknown data generating
distribution with the empirical distribution of the data. The training
stability is dramatically smaller than the validation stability and decreases
rapidly as $b$ decreases.

Figure 3. Sample Stability

Panel A: Treatment Effect Estimation | Panel B: Risk Estimation
---|---
|

Notes: Figure 3 displays estimates of the validation and training sample
stabilities, defined in 3.1, for the treatment effect and risk estimation
applications to the Banerjee et al., (2015) data introduced in Section 1. For
the risk estimation application, we study the regression (1.1) for the treated
outcome, i.e., $w$ equal to one, for a fixed choice of $\lambda$. We
approximate the unknown true data generating distribution with the empirical
distribution of the data. The $y$-axes are displayed in a log-scale base 10.
Further details on this figure are given in Appendix A.

An interesting structure emerges under the assumption that
$T\left(\mathsf{s},D\right)$ is sample stable. In particular, restricting
attention to the case of cross-fitting $b=n/k$, and plugging (3.3) into the
bound (3.5) yields

$\frac{\delta
gt^{2}}{2^{4}(1-\varphi)^{2}}\left(\frac{n-b}{b-1}\right)^{2}\geq\frac{\delta}{2^{4}}gk^{2}t^{2}$

In other words, in the case of cross-splitting, large deviation probabilities
converge to zero with an exponential rate that is quadratic in $k$. By
contrast, the exponential rate of concentration is always linear in $g$. That
is, for a fixed total number of splits $m=gk$, the exponential rates of
concentration for sample stable statistics aggregated with independent
splitting and with cross-splitting are $m$ and $m\cdot k$, respectively.
$\blacksquare$

Next, we derive bounds for the centered conditional moments of
$a\left(\mathsf{R}_{g,k},D\right)$ through an argument closely related to the
proof of Theorem 3.1. Bounds of this form are known as Burkholder-Davis-Gundy
inequalities (Burkholder,, 1973). By combining these bounds with a standard
Berry-Esseen inequality (Shevtsova,, 2011), we obtain a non-asymptotic central
limit theorem as a corollary.

###### Theorem 3.2.

Suppose that Assumptions 3.1 and 3.2 hold and that the data $D$ are
independent and identically distributed. If $r=2^{c-1}$ for some positive
integer $c$ and the $4r$th-order split stability $\zeta^{(4r)}$ is finite,
then the inequality

$\displaystyle\mathbb{E}\left[\left(a\left(\mathsf{R}_{g,k},D\right)-\overline{a}\left(D\right)\right)^{2r}\right]$
$\displaystyle\leq(2r-1)^{r}\left(\frac{2^{4}(2-\varphi
k-\varphi)^{2}}{g}\right)^{r}\Gamma_{k,\varphi,b}^{(r)}\leavevmode\nobreak\ ,$
(3.7)

holds, where

$\displaystyle\Gamma_{k,\varphi,b}^{(r)}$
$\displaystyle=\left(\frac{1-k\varphi}{1-\varphi}\right)2^{2r}\sigma_{\mathsf{max}}^{(2r,1)}+\left(\frac{k\varphi-\varphi}{1-\varphi}\right)\sigma_{\mathsf{va\ell
id}}^{(2r,b-1)}$ (3.8)

and $\varphi=b/n$.

######  Remark 3.6.

By setting $r=1$, the inequality (3.7) gives a variance bound. In this case,
the right-hand side of the inequality (3.7) is equal to negative one times the
inverse of the right-hand side of the inequality (3.5). In this sense, the
concentration inequality (3.5) can be thought of as an exponential Efron-Stein
inequality adapted to the dependence inherent in our problem (Boucheron et
al.,, 2003).999Abou-Moustafa and Szepesvári, (2019) give related exponential
Efron-Stein inequality in a different setting. For values of $r$ greater than
1, the unconditional quantity $\Gamma_{k,\varphi,b}^{(r)}$ again interpolates
between the $2r$th order full sample stability
$\sigma_{\mathsf{max}}^{(2r,1)}$ and the $(2r,b-1)$th order training sample
stability $\sigma_{\mathsf{va\ell id}}^{(2r,b-1)}$ as $k$ increases from $1$
to $1/\varphi$. For sample stable statistics under cross-fitting, the bound
(3.7) can be re-written

$\mathbb{E}\left[\left(a\left(\mathsf{R}_{g,k},D\right)-\overline{a}\left(D\right)\right)^{2r}\right]\lesssim\left(\frac{1}{gk^{2}}\right)^{r}$
(3.9)

by omitting constants that depend only on $r$. The rate of convergence of
central moments of cross-split sample stable statistics is much faster with
$k$ than it is with $g$.

The Kolmogorov distance between two measures $P$ and $Q$ on the real line is
given by

$d_{\text{K}}\left(P,Q\right)=\sup_{w\in\mathbb{R}}\left\\{|P\left\\{(-\infty,w]\right\\}-Q\left\\{(-\infty,w]\right\\}|\right\\}.$
(3.10)

With a mild abuse of notation, we let $d_{\text{K}}\left(X,Y\right)$ denote
the Kolmogorov distance between the probability measures of the real-valued
random variables $X$ and $Y$.

###### Corollary 3.1.

Let $W$ denote a standard normal random variable. Suppose that Assumptions 3.1
and 3.2 hold, that the data $D$ are independently and identically distributed,
and that the statistic $a(\mathsf{R}_{1,k},D)$ has a non-zero variance
conditional on $D$ almost surely. If the eighth-order split stability
$\zeta^{(8)}$ is finite, then for all $g$, the conditional Berry-Esseen
inequality

$\displaystyle
d_{\text{K}}\left(\frac{a\left(\mathsf{R}_{g,k},D\right)-\overline{a}\left(D\right)}{\sqrt{v_{g,k}\left(D\right)}},W\bigg{|}D\right)$
$\displaystyle\leq\frac{2^{6}3^{2}}{\delta}\frac{\left(2-\varphi
k-\varphi\right)^{3}}{g^{1/2}}\left(\frac{(\Gamma_{k,\varphi,b}^{(2)})^{1/2}}{v_{1,k}(D)}\right)^{3/2}$
(3.11)

is satisfied with probability greater than $1-\delta$ as $D$ varies.101010To
the best of our knowledge, the only Stein’s method central limit theorem in
the Kolmogorov distance, applicable to our setting, is given in Zhang, (2022),
which generalizes the argument of Shao and Zhang, (2019). In Section D.3, we
show that this central limit theorem, in conjunction with the bounds used to
derive Theorem 3.2, implies a bound that does not reduce as either $g$ or $k$
increase.

######  Remark 3.7.

Consider again the case of cross-splitting, where $k=n/b=1/\varphi$. By
Theorem 3.2, the Berry-Esseen bound (3.11) is at least

$\frac{(1-\varphi)^{3}}{g^{1/2}}\left(\frac{(\sigma^{(4,b-1)}_{\mathsf{train}})^{1/2}}{v_{1,k}(D)}\right)^{3/2}\gtrsim\frac{1}{g^{1/2}}\left(\frac{(\sigma^{(4,b-1)}_{\mathsf{train}})^{1/2}}{\sigma^{(2,b-1)}_{\mathsf{train}}}\right)^{3/2}\leavevmode\nobreak\
.$ (3.12)

Suppose that the statistic is sample stable. As long as
$\sigma^{(4,b-1)}_{\mathsf{train}}$ and $\sigma^{(2,b-1)}_{\mathsf{train}}$
are well-approximated by their upper bounds (3.3), the quantity (3.12) is of
order $g^{-1/2}$ and does not decrease as $k$ increases.

This argument suggests that, despite the fast decay in the residual randomness
demonstrated in Theorem 3.2, the quality of a normal approximation to
$a\left(\mathsf{R}_{g,k},D\right)$ does not increase as $k$ increases. This is
again born out in our running examples. Figure 4 displays the Kolmogorov
distance between a standard normal random variable and the studentized
distribution of the cross-fit statistic $a\left(\mathsf{R}_{1,k},D\right)$
conditioned on $D$ over a range of values of $k$ for our two applications to
the Banerjee et al., (2015) data. That is, Figure 4 plots the left-hand side
of the inequality (3.11). In both cases, despite a substantial initial
improvement, the quality of a normal approximation does not increase for large
$k$. In fact, in the risk estimation application, the quality of a normal
approximation decays for values of $k$ larger than 20. On the other hand, it
is worth noting that a normal approximation in the treatment effect
application is remarkably accurate for values of $k$ larger than 20.
$\blacksquare$

Figure 4. Quality of Normal Approximation

Panel A: Treatment Effect Estimation | Panel B: Risk Estimation
---|---
|

Notes: Figure 4 displays measurements of the Kolmogorov distance (3.10)
between a standard normal random variable and the studentized distribution of
the cross-fit statistic $a\left(\mathsf{R}_{1,k},D\right)$ conditioned on $D$
over a range of values of $k$ for our two applications to the Banerjee et al.,
(2015) data. This quantity is defined in the left hand side of the inequality
(3.11). The $y$-axis of Panel A is displayed in a log-scale base 10. For the
risk estimation application, we study the regression (1.1) for the treated
outcome, i.e., $w$ equal to one, for a fixed choice of $\lambda$. Further
details on this figure are given in Appendix A.

###  3.3. Reproducibility

We now return our attention to the sequential procedure for aggregating
sample-split statistics formulated in Section 2. Our main result is a Berry-
Esseen type bound on the accuracy of the nominal error rate $\beta$ for the
Anscombe-Chow-Robbins method specified in Algorithm 1.

###### Theorem 3.3.

Suppose that the statistic $a(\mathsf{R}_{1,k},D)$ has a non-zero and
variance, conditional on $D$, almost surely and that the collections
$\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^{\prime}_{\hat{g}^{\prime},k}$ are
independently computed with Algorithm 1. If Assumptions 3.1 and 3.2 hold, the
data $D$ are independent and identically distributed, and the eighth-order
split stability $\zeta^{(8)}$ is finite, then the Berry-Esseen inequality

$\displaystyle\bigg{|}P\bigg{\\{}\big{|}a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}^{\prime}_{\hat{g}^{\prime},k},D)\big{|}\leq\xi\mid
D\bigg{\\}}-(1-\beta)\bigg{|}$
$\displaystyle\lesssim\rho_{k,\varphi,b}(\xi,\beta\mid
D)+\lambda_{k,\varphi,b}(\xi,\beta\mid D)\leavevmode\nobreak\ ,$ (3.13)

holds with probability greater than $1-\delta$ as $D$ varies, for

$\displaystyle\rho_{k,\varphi,b}(\xi,\beta\mid D)$
$\displaystyle=\frac{\xi}{z_{1-\beta/2}}\frac{(2-\varphi
k-\varphi)^{3}(\Gamma^{(2)}_{k,\varphi,b})^{3/4}}{\delta\left(v_{1,k}\left(D\right)\right)^{2}}\quad\text{and}$
(3.14) $\displaystyle\lambda_{k,\varphi,b}(\xi,\beta\mid D)$
$\displaystyle=\left(\frac{\xi}{z_{1-\beta/2}}\frac{(2-\varphi
k-\varphi)^{4}\Gamma^{(1)}_{k,\varphi,b}(\Gamma_{k,\varphi,b}^{(2)})^{1/2}}{\delta^{3/2}(v_{1,k}\left(D\right))^{5/2}}\right)^{1/2}\leavevmode\nobreak\
,$ (3.15)

where a multiplicative factor that grows logarithmically as $\xi$ shrinks has
been omitted in writing (3.15) and an additional additive term that shrinks
exponentially as $\xi$ shrinks has been omitted in writing (3.13).

######  Remark 3.8.

Recall the definition of the oracle stopping time $g^{\star}$ given in (2.7).
The bound (3.13) is obtained by decomposing the approximation error into two
terms. The first term involves the error in a normal approximation to the
difference
$a(\mathsf{R}_{g^{\star},k},D)-a(\mathsf{R}^{\prime}_{g^{\star},k},D)$ and
contributes the term $\rho_{\varphi,k}(\xi,\beta\mid D)$. As $g^{\star}$ is
deterministic conditional on $D$, a bound on this term follows from an
argument very similar to the proof of 3.1. The second term involves the
differences $a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}_{g^{\star},k},D)$ and
$a(\mathsf{R}^{\prime}_{\hat{g}^{\prime},k},D)-a(\mathsf{R}^{\prime}_{g^{\star},k},D)$
and contributes the quantity $\lambda_{\varphi,k}(\xi,\beta\mid D)$. The key
step in bounding these quantities involves deriving a high probability bound
for the difference $\hat{g}-g^{\star}$. This follows by combining a
concentration inequality analogous to Theorem 3.1 for the conditional variance
estimator $\hat{v}_{g,k}(D)$ with a Bernstein-type maximal inequality, due to
Steiger, (1970). $\blacksquare$

######  Remark 3.9.

The dependence of the bound (3.13) on $\xi$ is sharp, at least up to the
logarithmic factor. This follows from the examples given in Landers and Rogge,
(1976, 1988). Berry-Esseen type bounds on the coverage error of Anscombe-Chow-
Robbins fixed-width sequential confidence intervals have been given in Csenki,
(1980), Mukhopadhyay, (1981), and Callaert and Janssen, (1981). Our result
differs in several respects. Most importantly for our application, the
dependence of our bound on the quantities $k$, $\varphi$, and $v_{1,k}(D)$ is
explicit. We explore the structure of the bound through these terms in the
following remark. Secondarily, we give an explicit characterization of the
logarithmic factors in the bound by using Bernstein-type concentration
inequalities, e.g., the bound given in Theorem 3.1, rather than moment bounds.
$\blacksquare$

######  Remark 3.10.

Restrict attention again to the case of cross-fitting and assume that the
statistic of interest is sample stable. We will refer to the procedure that
uses the oracle stopping time $g^{\star}$ in the place of $\hat{g}$ as the
oracle procedure. Observe that the oracle stopping time $g^{\star}$ is the
unique integer that satisfies

$g^{\star}-1\leq 2v_{1,k}(D)\left(\frac{z_{1-\beta/2}}{\xi}\right)^{2}\leq
g^{\star}\leavevmode\nobreak\ .$ (3.16)

Replacing the variance $v_{1,k}(D)$ by the bound derived in Theorem 3.2, we
get that $g^{\star}$ is proportional to

$\frac{2}{k^{2}}\left(\frac{z_{1-\beta/2}}{\xi}\right)^{2}\leavevmode\nobreak\
.$ (3.17)

Thus, if $m^{\star}=kg^{\star}$ is the total number of splits used by the
oracle procedure, then $m^{\star}$ is proportional to $k^{-1}$ as $k$
increases with all other quantities held fixed. In other words, the total
number of splits used by the oracle procedure decreases rapidly as $k$
increases. In the proof of Theorem 3.3, we give high probability bounds for
$|g^{\star}-\hat{g}|$. Consequently, similar intuition will hold for Algorithm
1. That is, with high probability, the computational cost scales like $k^{-1}$
as $k$ increases.

On the other hand, by Theorem 3.2, the leading term
$\lambda_{k,\varphi,b}(\xi,\beta\mid D)$ in the bound (3.13) is at least

$\left(\frac{\xi}{z_{1-\beta/2}}\frac{1}{\delta^{3/2}(1-\varphi)}\frac{(\sigma_{\mathsf{train}}^{(4,b-1)})^{1/2}}{(\sigma_{\mathsf{train}}^{(2,b-1)})^{3/2}}\right)^{1/2}\leavevmode\nobreak\
.$ (3.18)

If $\sigma_{\mathsf{train}}^{(4,b-1)}$ and $\sigma_{\mathsf{train}}^{(2,b-1)}$
are proportional to their upper bounds (3.3), then the quantity (3.18) is of
order

$\left(\frac{\xi}{z_{1-\beta/2}}\frac{k}{\delta^{3/2}}\right)^{1/2}\leavevmode\nobreak\
,$ (3.19)

which is in turn of order $(g^{\star})^{-1/4}$ by (3.17).

Comparison of the quantities (3.16) and (3.19) illustrates a fundamental
trade-off in the performance of Algorithm 1 over the choice of $k$. Increased
values of $k$ reduce the conditional variance $v_{g,k}(D)$, and thereby,
reduce the oracle sample size $g^{\star}$, easing computational cost. But as a
consequence, the quality of the accuracy of the nominal error rate may
deteriorate. In particular, at smaller values of $g$, the accuracy of the
estimator $\hat{v}_{g,k}(D)$ for the conditional variance $v_{g,k}(D)$ will be
lower, increasing the risk that $\hat{g}$ is smaller than $g^{\star}$. The
increased risk of early stopping decreases the probability that a result is
reproducible. It is desirable, then, to choose values for $k$ and $\xi$ such
that the variance $v_{1,k}(D)$ is small enough to make Algorithm 1
computationally efficient, but not so small that the effective sample size
$g^{\star}$ is too small for the nominal error rate $\beta$ to be accurate.

##  4\. Performance

We now study the performance of Algorithm 1 in the applications to the
Banerjee et al., (2015) data introduced in Section 1. Figure 5 displays
measurements of reproducibility error

$P\bigg{\\{}\big{|}a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}^{\prime}_{\hat{g}^{\prime},k},D)\big{|}\geq\xi\mid
D\bigg{\\}}$ (4.1)

when $\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^{\prime}_{\hat{g}^{\prime},k}$
are obtained independently with Algorithm 1 for the risk estimation
application. Figure A.2, given in Section A.3, displays analogous measurements
for the treatment effect estimation application. We set the nominal error rate
to $\beta=0.1$ and vary the values of $k$ and $\xi$.

Figure 5. Reproducibility Error in Risk Estimation

Panel A: Early Stopping with $g_{\mathsf{init}}=2$
---
Panel B: Reproducibility with $g_{\mathsf{init}}=10$

Notes: Figure 5 displays measurements of reproducibility error (4.1) for
Algorithm 1 in the risk estimation application to the Banerjee et al., (2015)
data. A dashed horizontal line is displayed at the nominal error rate
$\beta=0.1$. The $x$-axes vary the bound $\xi$ and are displayed in a log-
scale base 10. We study the regression (1.1) for the treated outcome, i.e.,
$w$ equal to one, for a fixed choice of $\lambda$. To ease computation, we
only display reproducibility error estimates for a $(\xi,k)$ pair if
$g^{\star}$ is less than $10^{4}$. Further details on this figure are given in
Appendix A.

Panel A displays measurements of the reproducibility error with the
initialization $g_{\mathsf{init}}$ set equal to 2, the minimal feasible value.
A transition is evident. For large values of $\xi$, the reproducibility error
is monotonically decreasing in $k$. As $\xi$ decreases, this pattern reverses
and the reproducibility probability approaches the nominal error rate. This
phenomenon is caused by severe early stopping at large values of $\xi$. In
particular, at large values of $\xi$, there is a substantial probability that
$\hat{g}$ is chosen to be a very small value, e.g., 2 or 3, often several
orders of magnitude less than $g^{\star}$.111111In the risk estimation
application, with $k=10$, $g^{\star}$ is 54.4 when $\xi$ is 10.4 and is 570.5
when $\xi$ is 3.2. We document this behavior in Figure A.3, given in Section
A.3.

This issue can be addressed through a very small increase in
$g_{\mathsf{init}}$. Panel B displays measurements of the reproducibility
error with the initialization $g_{\mathsf{init}}$ increased to 10. The nominal
error rate is now substantially more accurate over the full support of values
for $\xi$ considered. The moderate inaccuracy at larger values of $\xi$ has
two sources. First, there is still some risk of early stopping at values of
$g$ not much greater than $10$. Second, if $\xi$ is large relative to the
conditional variance $v_{1,k}(D)$, ensuring that $\hat{g}$ is at least 10 can
cause the reproducibility error to be close to 0. At moderate values of $\xi$,
the trade off documented in Section 3 is apparent. As $k$ increases, the
effective sample size $g^{\star}$ decreases, reducing the accuracy of the
nominal error rate. Figure A.3, given in Section A.3, demonstrates that this
inaccuracy can be largely explained by the error in the approximation of
$g^{\star}$ with $\hat{g}$. For small values of $\xi$, at the level of
precision likely to be of practical interest, the nominal error rate is very
accurate.

##  5\. A Stein Representer, Concentration, and Reproducibility

In this section, we outline the proofs for results stated in Section 3. The
primary difficulty is accommodating the dependence in the summands of
$a(\mathsf{R}_{g,k},D)$ across elements of the same cross-split. We tackle
this problem through the method of exchangeable pairs (Stein,, 1986).

Suppose that we are interested in studying the statistic $f\left(X\right)$,
where $X$ is random variable valued on the separable metric space
$\mathcal{X}$. The method of exchangeable pairs has two ingredients. First, we
need to construct a random variable $X^{\prime}$ such that
$\left(X,X^{\prime}\right)$ is an exchangeable pair. Second, we need to
construct an antisymmetric function $F\left(X,X^{\prime}\right)$ such that

$f\left(X\right)=\mathbb{E}\left[F\left(X,X^{\prime}\right)\mid X\right]$
(5.1)

almost surely. We will refer to the function $F\left(X,X^{\prime}\right)$ as a
“Stein representer” for $f\left(X\right)$.

###  5.1. Constructing a Stein Representer

In many applications, the Stein representer (5.1) induced by a suitable
exchangeable pair $\left(X,X^{\prime}\right)$ can be derived in closed form.
See e.g Chen et al., (2011) and Ross, (2011). A closed form derivation is more
challenging in our setting, where the collection $\mathsf{R}_{g,k}$ takes the
place of the random variable $X$. To address this issue, we apply a method due
to Chatterjee, (2005) for constructing Stein representers through a pair of
coupled Markov chains induced by an appropriately chosen exchangeable pair.

Chatterjee’s construction is founded on a observation, due to Stein, (1986),
that an exchangeable pair $\left(X,X^{\prime}\right)$ induces a reversible
Markov kernel $K$ through

$Kg(X)=\mathbb{E}\left[g(X^{\prime})\mid X=x\right]\leavevmode\nobreak\ ,$

where $g$ is any function satisfying $\mathbb{E}\left[|g(X)|\right]<\infty$.
Suppose that $\\{X_{m}\\}_{m\geq 0}$ and $\\{X^{\prime}_{m}\\}_{m\geq 0}$ are
two Markov chains constructed with the kernel induced by
$\left(X,X^{\prime}\right)$ and coupled in a such way that the marginal
distributions of $X_{m}$ and $X^{\prime}_{m}$ depend only on the initial
conditions $X_{0}$ and $X^{\prime}_{0}$, respectively. Chatterjee makes the
following observation. If there exists a constant $C$ such that

$\sum_{m=0}^{\infty}|\mathbb{E}\left[f(X_{m})-f(X_{m}^{\prime})\mid
X_{0}=x,X_{0}^{\prime}=y\right]|\leq C\leavevmode\nobreak\ ,$ (5.2)

for each $x$ and $y$ in $\mathcal{X}$, then the function

$F(x,y)=\sum_{m=0}^{\infty}\mathbb{E}\left[f(X_{m})-f(X_{m}^{\prime})\mid
X_{0}=x,X_{0}^{\prime}=y\right]$

is a Stein representer for $f\left(X\right)$. See Paulin et al., (2013, 2016)
for applications of this idea to the derivation of matrix concentration
inequalities. To apply this construction to our setting, we have two tasks.
First, we need to specify a suitable exchangeable pair. Second, we need to
construct a pair of coupled Markov chains induced by this pair that satisfy
the finiteness condition (5.2). Throughout, for a vector $x=(x_{i})_{i=1}^{n}$
we let $(x_{-\ell},y)$ denote the vector formed by replacing the $\ell$th
component of $x$ with $y$.

Our construction is premised on the observation that the random collection of
splits $\mathsf{R}_{g,k}$ can be generated by a random collection of
permutations. To see this, let $\mathcal{P}_{n}$ denote the set of
permutations of the set $\left[n\right]$, treating each
$\pi\in\mathcal{P}_{n}$ as a bijection from $\left[n\right]$ to
$\left[n\right]$. Observe that each permutation $\pi\in\mathcal{P}_{n}$ can be
associated with an element of $\mathcal{R}_{n,k,b}$, denoted by
$\mathsf{\bm{r}}_{k}\left(\pi\right)=\left(\mathsf{s}_{1}\left(\pi\right),\ldots,\mathsf{s}_{k}\left(\pi\right)\right)$,
through

$\mathsf{s}_{i}\left(\pi\right)=\left\\{\pi\left(k\cdot\left(i-1\right)+1\right),\ldots,\pi\left(k\cdot\left(i-1\right)+b\right)\right\\}.$

If $\bm{\pi}=\left(\pi_{i}\right)_{i=1}^{g}$denotes a collection of
permutations drawn independently and uniformly at random from
$\mathcal{P}_{n}$, then the collection
$\mathsf{R}_{g,k}\left(\bm{\pi}\right)=\left(\mathsf{r}_{k}\left(\pi_{i}\right)\right)_{i=1}^{g}$
is equidistributed with the collection $\mathsf{R}_{g,k}$ defined in Section
2.

Now, we construct an exchangeable pair $(\bm{\pi},\bm{\pi}^{\prime})$, keeping
in mind that our aim is to verify a condition of the form (5.2). For any
permutation $\pi\in\mathcal{P}_{n}$ and indices $i,j\in\left[n\right]$, define
the updated permutation

$\hat{\pi}\left(i,j\right)\left(x\right)=\begin{cases}j,&x=i,\\\
\pi\left(i\right),&x=\pi^{-1}\left(j\right),\\\
\pi\left(x\right),&\text{otherwise.}\end{cases}$

In other words, $\hat{\pi}\left(i,j\right)$ is identical to $\pi$, except that
$i$ maps to $j$ and $\pi^{-1}\left(j\right)$ maps to $\pi\left(i\right)$. Let
$L$ be distributed uniformly on $\left[g\right]$ and let $I$ and $J$ be
independently and uniformly distributed on $\left[n\right]$. Define the
modified collection

$\displaystyle\bm{\pi}^{\prime}$
$\displaystyle=(\pi_{-L},\hat{\pi}_{L}(I,J)).$

and observe that $(\bm{\pi},\bm{\pi}^{\prime})$ is an exchangeable pair.

With the exchangeable pair $(\bm{\pi},\bm{\pi}^{\prime})$ in place, we choose
a coupled pair of Markov chains that it induces. For each $m\geq 1$, let
$L_{m}$ be distributed uniformly on $\left[g\right]$ and let $I_{m}$ and
$J_{m}$ be distributed uniformly on $\left[n\right]$. Construct
$\left(\bm{\pi}_{m},\bm{\pi}_{m}^{\prime}\right)$ from the pair
$\left(\bm{\pi}_{0},\bm{\pi}^{\prime}_{0}\right)=\left(\bm{\pi},\bm{\pi}^{\prime}\right)$
by setting

$\displaystyle\bm{\pi}_{m}$
$\displaystyle=(\pi_{m-1,-L_{m}},\hat{\pi}_{m-1,L_{m}}(I_{m},J_{m}))\quad\text{and}\quad\bm{\pi}_{m}^{\prime}=(\pi_{m-1,-L_{m}}^{\prime},\hat{\pi}_{m-1,L_{m}}^{\prime}(I_{m},J_{m}))$
(5.3)

recursively. In other words, the Markov chain $\bm{\pi}_{0}^{\prime}$ is
initialized by choosing one permutation in $\bm{\pi}_{0}$ and swapping one
pair of indices. In the $m$th iteration of the Markov chain
$(\bm{\pi}_{m},\bm{\pi}_{m}^{\prime})_{m\geq 0}$, the permutations
$\pi_{m-1,L_{m}}$ and $\pi^{\prime}_{m-1,L_{m}}$ are selected and updated so
that $I_{m}$ now maps to $J_{m}$. As the iterations proceed, $I_{m}$ will
continue to map to the same index in the $L_{m}$th permutation in both
collections. That is, in each iteration that a new permutation $L$ and index
$I$ are selected, the two collections become more similar. Eventually, every
$L$ and $I$ will have been selected, the two collections $\bm{\pi}_{m}$ and
$\bm{\pi}_{m}^{\prime}$ will be the same, and so
$a(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D)-a(\mathsf{R}_{g,k}(\bm{\pi}_{m}^{\prime}),D)=0$.
This will imply the finiteness (5.2). This argument is formalized in the proof
of the following Lemma.

###### Lemma 5.1.

The function

$A\left(\bm{\pi},\bm{\pi}^{\prime}\mid
D\right)=\sum_{m=0}^{\infty}\mathbb{E}\left[\left(a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}^{\prime}),D\right)\right)\mid\bm{\pi}_{0}=\bm{\pi},\bm{\pi}_{0}^{\prime}=\bm{\pi}^{\prime}\right],$

is finite, antisymmetric, and satisfies the equality

$\mathbb{E}\left[A\left(\bm{\pi},\bm{\pi}^{\prime}\mid
D\right)\mid\bm{\pi},D\right]=a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)$
(5.4)

almost surely.

###  5.2. Proofs for Theorems 3.1, 3.2, and 3.1

To obtain the concentration inequality (3.5) and Burkholder-Davis-Gundy
inequality (3.7), we apply the following result due to Chatterjee, (2005,
2007), which we have augmented to be applicable under a set of assumptions
considered in Paulin et al., (2016).

###### Lemma 5.2.

Let $\mathcal{X}$ be a separable metric space and suppose that
$\left(X,X^{\prime}\right)$ is an exchangeable pair of $\mathcal{X}$-valued
random variables. Let $f:\mathcal{X}\to\mathbb{R}$ be a square integrable
function and let $F:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be a Stein
representer for $f$. For each positive integer $(r)$, define the quantities

$U_{f}^{(r)}\left(X\right)=\frac{1}{2}\mathbb{E}\left[\left(f\left(X\right)-f\left(X^{\prime}\right)\right)^{2r}\mid
X\right]\quad\text{and}\quad
U_{F}^{(r)}\left(X\right)=\frac{1}{2}\mathbb{E}\left[F\left(X,X^{\prime}\right)^{2r}\mid
X\right].$ (5.5)

If there exist nonnegative constants $u$ and $s$ such that

$U^{(1)}_{f}\left(X\right)\leq s^{-1}u\quad\text{and}\quad
U^{(1)}_{F}\left(X\right)\leq su,$ (5.6)

then the concentration inequality

$P\left\\{|f\left(X\right)|\geq\delta\right\\}\leq
2\exp\left(-t^{2}/2u\right)$ (5.7)

holds for all $t\geq 0$. Moreover, the moment inequality

$\mathbb{E}\left[f\left(X\right)^{2r}\right]\leq\left(2r-1\right)^{r}\left(s\mathbb{E}\left[U_{F}^{(r)}(X)\right]+s^{-1}\mathbb{E}\left[U_{f}^{(r)}(X)\right]\right)$
(5.8)

holds for all positive integers $r$ for any positive constant $s$.

Throughout, to ease notation, we use the short hand

$Z=(\mathsf{R}_{g,k}(\bm{\pi}),D)\quad\text{and}\quad
Z^{\prime}=(\mathsf{R}_{g,k}(\bm{\pi}^{\prime}),D)\leavevmode\nobreak\ .$

The Markov chains $(Z_{m})_{m\geq 0}$ and $(Z^{\prime}_{m})_{m\geq 0}$ are
defined accordingly and we simplify $A\left(\bm{\pi},\bm{\pi}^{\prime}\mid
D\right)$ to $A(Z,Z^{\prime})$. To apply Lemma 5.2, we are required to develop
bounds for the objects

$\displaystyle U_{a}^{(r)}\left(Z\right)$
$\displaystyle=\frac{1}{2}\mathbb{E}\left[\left(a\left(Z\right)-a\left(Z^{\prime}\right)\right)^{2r}\mid
Z\right]\quad\text{and}\quad
U_{A}^{(r)}\left(Z\right)=\frac{1}{2}\mathbb{E}\left[A\left(Z,Z^{\prime}\right)^{2r}\mid
Z\right].$

Our approach is based on the following Lemma, which combines a generalization
of an idea due to Lemma 10.4 of Paulin et al., (2016) with a Markov type
bound.

###### Lemma 5.3.

Let $r=2^{c}$ for some positive integer $c$. Let $f:\mathcal{Z}\to\mathbb{R}$
be a function such that there exists a square integrable random variable $W$
with

$\left(\sum_{m=0}^{\infty}\mathbb{E}\left[f\left(Z_{m}\right)-f\left(Z_{m}^{\prime}\right)\mid
Z_{0}=Z,Z_{0}^{\prime}=Z^{\prime}\right]\right)^{r}\leq W$ (5.9)

almost surely. If the inequality

$\mathbb{E}\left[\mathbb{E}\left[f(Z_{m})-f(Z_{m}^{\prime})\mid
Z_{0}=Z,Z_{0}^{\prime}=Z^{\prime}\right]^{r}\mid\bm{\pi}\right]\leq h_{m}^{r}$
(5.10)

holds for each $m\geq 0$ and each collection $\bm{\pi}$, where
$\left(h_{m}\right)_{m\geq 0}$ is a deterministic sequence of nonnegative
numbers, then the inequalities

$\displaystyle\frac{1}{2}\mathbb{E}\left[\left(f\left(Z\right)-f\left(Z^{\prime}\right)\right)^{r}\mid
Z\right]$ $\displaystyle\leq\frac{h_{0}^{r}}{\delta}\quad\text{and}$ (5.11)
$\displaystyle\frac{1}{2}\mathbb{E}\left[\left(\sum_{m=0}^{\infty}\mathbb{E}\left[f\left(Z_{m}\right)-f\left(Z_{m}^{\prime}\right)\mid
Z_{0}=Z,Z_{0}^{\prime}=Z^{\prime}\right]\right)^{r}\mid Z\right]$
$\displaystyle\leq\frac{1}{\delta}\left(\sum_{m=0}^{\infty}h_{m}\right)^{r}$
(5.12)

both hold with probability greater than $1-\delta$ as $D$ varies

To apply this Lemma, we begin by noting that the inequality

$\displaystyle\left(\sum_{m=0}^{\infty}\mathbb{E}\left[a\left(Z_{m}\right)-a\left(Z_{m}^{\prime}\right)\mid
Z_{0}=Z,Z_{0}^{\prime}=Z^{\prime}\right]\right)^{2}$ $\displaystyle\lesssim
g^{2}n^{4}\max_{\mathsf{s},\mathsf{s}^{\prime}\in\mathsf{S}_{n,b}}\left(T\left(\mathsf{s},U,D\right)-T\left(\mathsf{s}^{\prime},U,D\right)\right)^{2}$
(5.13)

follows from Lemma B.1, stated in Section B.2. Moreover, the right hand side
of (5.13) is square integrable, as the fourth order split-stability
$\zeta^{(4)}$ is finite by assumption. Deterministic bounds of the form (5.10)
are obtained through the following Lemma.

###### Lemma 5.4.

Under Assumptions 3.1 and 3.2, for all integers $m\geq 0$ and $r\geq 1$, the
inequality

$\displaystyle\mathbb{E}\left[\mathbb{E}\left[a(Z_{m})-a(Z_{m}^{\prime})\mid
Z_{0}=Z,Z_{0}^{\prime}=Z^{\prime}\right]^{2r}\mid\bm{\pi}\right]$
$\displaystyle\leq
2^{4r}\left(1-\frac{2}{gn^{2}}\right)^{2mr}\left(\frac{2n-bk-b}{gn^{2}}\right)^{2r}\Gamma_{k,\varphi,b}^{(r)}$

holds almost surely, where $\Gamma_{k,\varphi,b}^{(r)}$ is defined in the
statement of Theorem 3.2.

Combining Lemmas 5.3 and 5.4, we have that

$\displaystyle U_{A}^{(r)}(Z)$
$\displaystyle\leq\frac{1}{\delta}\left(\frac{gn^{2}}{2}\right)^{r}\left(\frac{2^{4}\left(2n-bk-b\right)^{2}}{gn^{2}}\right)^{r}\Gamma_{k,\varphi,b}^{(r)}$

and

$\displaystyle U_{a}^{(r)}(Z)$
$\displaystyle\leq\frac{1}{\delta}\left(\frac{2}{gn^{2}}\right)^{r}\left(\frac{2^{4}\left(2n-bk-b\right)^{2}}{gn^{2}}\right)^{r}\Gamma_{k,\varphi,b}^{(r)}$

with probability $1-\delta$. Theorem 3.1 is obtained by setting $r=1$ and
applying (5.7) of Lemma 5.2 with $s=gn^{2}/2$ and

$u=\frac{1}{\delta}\frac{2^{4}\left(2n-bk-b\right)^{2}}{gn^{2}}\Gamma_{n,k,b}^{(1)}.$

Similarly, Theorem 3.2 is obtained by applying (5.8) of Lemma 5.2 with
$s=(gn^{2}/2)^{r}$.

The normal approximation error bound (3.11) is a corollary of Theorem 3.2. To
see this, consider the centered statistic

$a\left(Z\right)-\overline{a}\left(D\right)=\frac{1}{g}\sum_{\ell=1}^{g}\overline{a}\left(\mathsf{r}_{\ell},D\right),\quad\text{where}\quad\overline{a}\left(\mathsf{r}_{\ell},D\right)=\frac{1}{k}\sum_{i=1}^{k}\left(T(\mathsf{s}_{\ell,i},D)-\overline{a}\left(D\right)\right)$
(5.14)

for each $\ell$ in $\left[g\right]$. Conditional on the data $D$, the
statistics $\overline{a}\left(\mathsf{r}_{\ell},D\right)$ are independent,
identically distributed, and mean zero. Moreover, we have that

$\displaystyle\mathbb{E}\left[|\overline{a}\left(\mathsf{r}_{\ell},D\right)|^{3}\right]$
$\displaystyle\leq 3^{2}2^{6}\left(2-\varphi
k-\varphi\right)^{3}(\Gamma_{k,\varphi,b}^{(2)})^{3/4}$ (5.15)

by Hölder’s inequality and Theorem 3.2. Hence, we find that

$\displaystyle\frac{\mathbb{E}\left[|\overline{a}\left(\mathsf{r}_{\ell},D\right)|^{3}\mid
D\right]}{g^{1/2}\left(v_{1,k}(D)\right)^{3/2}}$
$\displaystyle\leq\frac{3^{2}2^{6}}{\delta}\frac{\left(2-\varphi
k-\varphi\right)^{3}}{\left(v_{1,k}(D)\right)^{3/2}}\frac{(\Gamma_{k,\varphi,b}^{(2)})^{3/4}}{g^{1/2}}$

holds with probability greater than $1-\delta,$ by combining (5.15) with the
Markov inequality. The proof then follows by the standard Berry-Esseen
inequality. See e.g., Corollary 1 of Shevtsova, (2011).

###  5.3. Proof of Theorem 3.3

The first step for verifying Theorem 3.3 is deriving the exponential rate of
concentration for the variance estimator
$\hat{v}\left(\mathsf{R}_{g,k},D\right)$ defined in (2.4). This is obtained in
the following Lemma, which follows from an argument very similar to the proof
of Theorem 3.1.

###### Lemma 5.5.

Suppose that Assumptions 3.1 and 3.2 hold and that the data $D$ are
independent and identically distributed. If the eighth-order split-stability
$\zeta^{(8)}$ is finite, then the conditional concentration inequality

$\displaystyle\log\frac{1}{4}P\left\\{\bigg{|}\frac{\hat{v}(\mathsf{R}_{g,k},D)}{v_{g,k}\left(D\right)}-1\bigg{|}\geq
t\mid
D\right\\}\lesssim-\frac{\delta\left(v_{1,k}(D)\right)^{2}}{\left(2-\varphi
k-\varphi\right)^{4}}\frac{gt^{2}}{\Gamma_{k,\varphi,b}^{(2)}}$

holds for all $t>0$ with probability greater than $1-\delta$ as $D$ varies.

We focus our analysis on the error

$|P\left\\{a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\leq\xi\mid
D\right\\}-\left(1-\beta/2\right)|.$ (5.16)

An analogous argument will yield the same bound for the lower tail. We begin
by bounding (5.16) with quantities that will be easier to handle in isolation.
Define the objects

$\displaystyle U\left(\mathsf{R}_{g,k},D\right)$
$\displaystyle=\frac{1}{g^{\star}}\left(\sum_{i=1}^{g}\overline{a}\left(\mathsf{r}_{i,k},D\right)-\sum_{i=1}^{g^{\star}}\overline{a}\left(\mathsf{r}_{i,k},D\right)\right)\quad\text{and}$
(5.17) $\displaystyle Q\left(\mathsf{R}_{g,k},D\right)$
$\displaystyle=\left(1-\frac{g}{g^{\star}}\right)\left(a(\mathsf{R}_{g,k},D)-\overline{a}\left(D\right)\right)\leavevmode\nobreak\
.$ (5.18)

The following Lemma bounds the error (5.16) in terms of the error in the
normal approximation to the quantity
$a(\mathsf{R}_{g^{\star},k},D)-a(\mathsf{R}_{g^{\star},k}^{\prime},D)$ and
generic high probability bounds on (5.17) and (5.18).

###### Lemma 5.6.

Let $W$ denote a standard normal random variable. Define the events

$\displaystyle\mathcal{U}_{k,\lambda}\left(D\right)$
$\displaystyle=\left\\{|U(\mathsf{R}_{\hat{g},k},D)-U(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)|\leq\lambda\sqrt{2v_{g^{\star},k}\left(D\right)}\right\\}\quad\text{and}$
$\displaystyle\mathcal{Q}_{k,\lambda}(D)$
$\displaystyle=\left\\{\Big{|}Q(\mathsf{R}_{g,k},D)-Q(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\Big{|}\leq\lambda\sqrt{2v_{g^{\star},k}\left(D\right)}\right\\}.$

The quantity (5.16) is bounded above by

$d_{K}\left(\frac{a(\mathsf{R}_{g^{\star},k},D)-a(\mathsf{R}_{g^{\star},k}^{\prime},D)}{\sqrt{2v_{g^{\star},k}\left(D\right)}},W\mid
D\right)+2\lambda+\left(1-P\left\\{\mathcal{U}_{k,\lambda}\left(D\right)\cap\mathcal{Q}_{k,\lambda}(D)\mid
D\right\\}\right)\leavevmode\nobreak\ ,$ (5.19)

where the Kolmogorov distance $d_{K}\left(\cdot,\cdot\right)$ is defined in
(3.10).

Thus, it remains to give suitable bounds for the objects in (5.19). These are
obtained in the following Lemma.

###### Lemma 5.7.

Let $W$ denote a standard normal random variable. Recall the objects
$\rho_{\varphi,k}\left(\xi,\beta\mid D\right)$ and
$\lambda_{\varphi,k}\left(\xi,\beta\mid D\right)$ defined in the statement of
Theorem 3.3. Suppose that Assumptions 3.1 and 3.2 hold, that the data $D$ are
independent and identically distributed, and that the eighth-order split
stability $\zeta^{(8)}$ is finite.

(i) The Berry-Esseen inequality

$\displaystyle
d_{K}\left(\frac{a(\mathsf{R}_{g^{\star},k},D)-a(\mathsf{R}_{g^{\star},k}^{\prime},D)}{\sqrt{2v_{g^{\star},k}\left(D\right)}},W\mid
D\right)$ $\displaystyle\lesssim\rho_{k,\varphi,b}\left(\xi,\beta\mid
D\right)$

is satisfied with probability greater than $1-\delta$ as $D$ varies.

(ii) The conditional concentration inequality

$\displaystyle
P\left\\{\frac{|U(\mathsf{R}_{\hat{g},k},D)-U(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)|}{\sqrt{2v_{g^{\star},k}\left(D\right)}}\geq\lambda_{k,\varphi,b}\left(\xi,\beta\mid
D\right)\tilde{\lambda}_{k,\varphi,b}\left(D\right)\mid
D\right\\}\lesssim\rho_{k,\varphi,b}\left(\xi,\beta\mid D\right)$

holds with probability greater than $1-\delta$ as $D$ varies, where
$\tilde{\lambda}_{\varphi,k}\left(D\right)$ is a term that grows
logarithmically as $\xi$ shrinks and is characterized explicitly in the proof.

(iii) The conditional concentration inequality

$P\left\\{\frac{|Q(\mathsf{R}_{\hat{g},k},D)-Q(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)|}{\sqrt{2v_{g^{\star},k}\left(D\right)}}\geq\lambda_{k,\varphi,b}\left(\xi,\beta\mid
D\right)\right\\}\lesssim\rho_{k,\varphi,b}\left(\xi,\beta\mid
D\right)+\tilde{\rho}_{k,\varphi,b}\left(\xi,\beta\mid D\right)$

holds with probability greater than $1-\delta$ as $D$ varies, where
$\tilde{\rho}_{\varphi,k}\left(\xi,\beta\mid D\right)$ is a term that shrinks
exponentially as $\xi$ shrinks and is characterized explicitly in the proof.

The proof is complete by substituting the bounds specified in Lemma 5.7 into
Lemma 5.6.

##  6\. Conclusion

We propose a method for sequentially aggregating randomized statistics to
ensure that residual randomness is small. The method is applicable under very
general conditions. By restricting the set of statistics under consideration,
we give a non-asymptotic analysis of the performance of the method. We have
two main findings. First, aggregation of sample-split statistics with cross-
splitting reduces residual randomness at a much faster rate than for
independent splitting. Second, cross-splitting does not necessarily improve
the quality of the nominal error rate of the procedure. Thus, users navigate a
tradeoff. Finer cross-splitting reduces the conditional variance of the
aggregate statistic, reducing the computation needed to achieve a given bound
on the residual randomness. But as a consequence, the quality of the nominal
error rate may have also been reduced, limiting the user’s ability to give a
high quality estimate of the probability that their results are reproducible.

Our computational analysis has highlighted the substantial practical
importance of choosing the initialization parameter $g_{\mathsf{init}}$ to be
sufficiently large. By contrast, our theoretical analysis is agnostic to this
issue. In the context of fixed-length sequential confidence intervals,
Mukhopadhyay and Datta, (1996) study the choice of $g_{\mathsf{init}}$. They
show that if $g_{\mathsf{init}}$ is chosen appropriately, a non-uniform
asymptotic approximation to the coverage error, based on Edgeworth expansions
due to Aras and Woodroofe, (1993), has a leading term of order $O(\xi^{2})$.
This is a substantial improvement on the $O(\xi^{1/2})$ rate given in Theorem
3.3. Further research should give a non-asymptotic treatment of the choice of
the parameter $g_{\mathsf{init}}$, aimed at assessing whether the
$O(\xi^{1/2})$ rate can be improved through further assumptions and a
principled choice of $g_{\mathsf{init}}$.

## References

  * Abou-Moustafa and Szepesvári, (2019) Abou-Moustafa, K. and Szepesvári, C. (2019).  An exponential efron-stein inequality for $\ell_{q}$ stable learning rules.  In Algorithmic Learning Theory, pages 31–63. PMLR.
  * Anscombe, (1952) Anscombe, F. J. (1952).  Large-sample theory of sequential estimation.  In Mathematical Proceedings of the Cambridge Philosophical Society, volume 48, pages 600–607. Cambridge University Press.
  * Aras and Woodroofe, (1993) Aras, G. and Woodroofe, M. (1993).  Asymptotic expansions for the moments of a randomly stopped average.  The annals of Statistics, pages 503–519.
  * Arlot and Celisse, (2010) Arlot, S. and Celisse, A. (2010).  A survey of cross-validation procedures for model selection.  Statistics Surveys, 4:40 – 79.
  * Austern and Zhou, (2020) Austern, M. and Zhou, W. (2020).  Asymptotics of cross-validation.  arXiv preprint arXiv:2001.11111.
  * Banerjee et al., (2015) Banerjee, A., Duflo, E., Goldberg, N., Karlan, D., Osei, R., Parienté, W., Shapiro, J., Thuysbaert, B., and Udry, C. (2015).  A multifaceted program causes lasting progress for the very poor: Evidence from six countries.  Science, 348(6236):1260799.
  * Bates et al., (2023) Bates, S., Hastie, T., and Tibshirani, R. (2023).  Cross-validation: what does it estimate and how well does it do it?  Journal of the American Statistical Association, pages 1–12.
  * Bayle et al., (2020) Bayle, P., Bayle, A., Janson, L., and Mackey, L. (2020).  Cross-validation confidence intervals for test error.  Advances in Neural Information Processing Systems, 33:16339–16350.
  * Beran and Millar, (1987) Beran, R. and Millar, P. W. (1987).  Stochastic estimation and testing.  The Annals of Statistics, 15(3):1131–1154.
  * Bickel, (1982) Bickel, P. J. (1982).  On adaptive estimation.  The Annals of Statistics, pages 647–671.
  * Boucheron et al., (2003) Boucheron, S., Lugosi, G., and Massart, P. (2003).  Concentration inequalities using the entropy method.  The Annals of Probability, 31(3):1583–1614.
  * Bousquet and Elisseeff, (2002) Bousquet, O. and Elisseeff, A. (2002).  Stability and generalization.  The Journal of Machine Learning Research, 2:499–526.
  * Burkholder, (1973) Burkholder, D. L. (1973).  Distribution function inequalities for martingales.  The Annals of Probability, 1(1):19–42.
  * Callaert and Janssen, (1981) Callaert, H. and Janssen, P. (1981).  The convergence rate of fixed-width sequential confidence intervals for the mean.  Sankhyā: The Indian Journal of Statistics, Series A, pages 211–219.
  * Celisse and Guedj, (2016) Celisse, A. and Guedj, B. (2016).  Stability revisited: new generalisation bounds for the leave-one-out.  arXiv preprint arXiv:1608.06412.
  * Chatterjee, (2005) Chatterjee, S. (2005).  Concentration inequalities with exchangeable pairs (Ph. D. thesis).  arXiv preprint math/0507526.
  * Chatterjee, (2007) Chatterjee, S. (2007).  Stein’s method for concentration inequalities.  Probability Theory and Related Fields, 1(138):305–321.
  * Chen and Ritzwoller, (2023) Chen, J. and Ritzwoller, D. M. (2023).  Semiparametric estimation of long-term treatment effects.  Journal of Econometrics, 237(2):105545.
  * Chen et al., (2011) Chen, L. H., Goldstein, L., and Shao, Q.-M. (2011).  Normal approximation by Stein’s method, volume 2.  Springer.
  * Chen et al., (2022) Chen, Q., Syrgkanis, V., and Austern, M. (2022).  Debiased machine learning without sample-splitting for stable estimators.  Advances in Neural Information Processing Systems, 35:3096–3109.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters: Double/debiased machine learning.  The Econometrics Journal, 21(1).
  * Chow and Robbins, (1965) Chow, Y. S. and Robbins, H. (1965).  On the asymptotic theory of fixed-width sequential confidence intervals for the mean.  The Annals of Mathematical Statistics, 36(2):457–462.
  * Cornec, (2010) Cornec, M. (2010).  Concentration inequalities of the cross-validation estimate for stable predictors.  arXiv preprint arXiv:1011.5133.
  * Csenki, (1980) Csenki, A. (1980).  On the convergence rate of fixed-width sequential confidence intervals.  Scandinavian Actuarial Journal, 1980(2):107–111.
  * Dembo, (2021) Dembo, A. (2021).  Probability theory: Stat310/math230.
  * DiCiccio and Romano, (2019) DiCiccio, C. and Romano, J. P. (2019).  Multiple data splitting for testing.  Department of Statistics, Stanford University Stanford, CA.
  * DiCiccio et al., (2020) DiCiccio, C. J., DiCiccio, T. J., and Romano, J. P. (2020).  Exact tests via multiple data splitting.  Statistics & Probability Letters, 166:108865.
  * Dietterich, (1998) Dietterich, T. G. (1998).  Approximate statistical tests for comparing supervised classification learning algorithms.  Neural Computation, 10(7):1895–1923.
  * Du et al., (2023) Du, J.-H., Patil, P., Roeder, K., and Kuchibhotla, A. K. (2023).  Extrapolated cross-validation for randomized ensembles.  arXiv preprint arXiv:2302.13511.
  * Elisseeff et al., (2005) Elisseeff, A., Evgeniou, T., Pontil, M., and Kaelbing, L. P. (2005).  Stability of randomized learning algorithms.  Journal of Machine Learning Research, 6(1).
  * Friedman et al., (2021) Friedman, J., Hastie, T., Tibshirani, R., Narasimhan, B., Tay, K., Simon, N., and Qian, J. (2021).  Package ‘glmnet’.  CRAN R Repositary, 595.
  * Guo and Shah, (2023) Guo, F. R. and Shah, R. D. (2023).  Rank-transformed subsampling: inference for multiple data splitting and exchangeable p-values.  arXiv preprint arXiv:2301.02739.
  * Guo and Romano, (2017) Guo, W. and Romano, J. P. (2017).  Analysis of error control in large scale two-stage multiple hypothesis testing.  arXiv preprint arXiv:1703.06336.
  * Gut, (2009) Gut, A. (2009).  Stopped random walks.  Springer.
  * Hahn, (1998) Hahn, J. (1998).  On the role of the propensity score in efficient semiparametric estimation of average treatment effects.  Econometrica, pages 315–331.
  * Hardt et al., (2016) Hardt, M., Recht, B., and Singer, Y. (2016).  Train faster, generalize better: Stability of stochastic gradient descent.  In International Conference on Machine Learning, pages 1225–1234. PMLR.
  * Hickernell et al., (2013) Hickernell, F. J., Jiang, L., Liu, Y., and Owen, A. B. (2013).  Guaranteed conservative fixed width confidence intervals via monte carlo sampling.  In Monte Carlo and Quasi-Monte Carlo Methods 2012, pages 105–128. Springer.
  * Kale et al., (2011) Kale, S., Kumar, R., and Vassilvitskii, S. (2011).  Cross-validation and mean-square stability.  In ICS, pages 487–495.
  * Klaassen, (1987) Klaassen, C. A. (1987).  Consistent estimation of the influence function of locally asymptotically linear estimators.  The Annals of Statistics, 15(4):1548–1562.
  * Knight and Fu, (2000) Knight, K. and Fu, W. (2000).  Asymptotics for lasso-type estimators.  The Annals of Statistics, 28(5):1356–1378.
  * Kumar et al., (2013) Kumar, R., Lokshtanov, D., Vassilvitskii, S., and Vattani, A. (2013).  Near-optimal bounds for cross-validation via loss stability.  In International Conference on Machine Learning, pages 27–35. PMLR.
  * Landers and Rogge, (1976) Landers, D. and Rogge, L. (1976).  The exact approximation order in the central-limit-theorem for random summation.  Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 36(4):269–283.
  * Landers and Rogge, (1988) Landers, D. and Rogge, L. (1988).  Sharp orders of convergence in the random central limit theorem.  Journal of Approximation Theory, 53(1):86–111.
  * Lei, (2020) Lei, J. (2020).  Cross-validation with confidence.  Journal of the American Statistical Association, 115(532):1978–1997.
  * Levin and Peres, (2017) Levin, D. A. and Peres, Y. (2017).  Markov chains and mixing times, volume 107.  American Mathematical Soc.
  * Meinshausen and Bühlmann, (2010) Meinshausen, N. and Bühlmann, P. (2010).  Stability selection.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 72(4):417–473.
  * Meinshausen et al., (2009) Meinshausen, N., Meier, L., and Bühlmann, P. (2009).  P-values for high-dimensional regression.  Journal of the American Statistical Association, 104(488):1671–1681.
  * Mukhopadhyay, (1981) Mukhopadhyay, N. (1981).  Convergence rates of sequential confidence intervals and tests for the mean of a u-statistic.  Communications in Statistics-Theory and Methods, 10(21):2231–2244.
  * Mukhopadhyay and Datta, (1996) Mukhopadhyay, N. and Datta, S. (1996).  On sequential fixed-width confidence intervals for the mean and second-order expansions of the associated coverage probabilities.  Annals of the Institute of Statistical Mathematics, 48:497–507.
  * Nadeau and Bengio, (1999) Nadeau, C. and Bengio, Y. (1999).  Inference for the generalization error.  Advances in Neural Information Processing Systems, 12.
  * Paulin et al., (2013) Paulin, D., Mackey, L., and Tropp, J. A. (2013).  Deriving matrix concentration inequalities from kernel couplings.  arXiv preprint arXiv:1305.0612.
  * Paulin et al., (2016) Paulin, D., Mackey, L., and Tropp, J. A. (2016).  Efron-Stein inequalities for random matrices.  The Annals of Probability, 44(5):3431 – 3473.
  * Politis et al., (1999) Politis, D. N., Romano, J. P., and Wolf, M. (1999).  Subsampling.  Springer Science & Business Media.
  * Rényi, (1957) Rényi, A. (1957).  On the asymptotic distribution of the sum of a random number of independent random variables.  Acta Math, 8:193–199.
  * Rinaldo et al., (2019) Rinaldo, A., Wasserman, L., and G’Sell, M. (2019).  Bootstrapping and sample splitting for high-dimensional, assumption-lean inference.  The Annals of Statistics, 47(6):3438–3469.
  * Robins et al., (1994) Robins, J. M., Rotnitzky, A., and Zhao, L. P. (1994).  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American statistical Association, 89(427):846–866.
  * Ross, (2011) Ross, N. (2011).  Fundamentals of Stein’s method.  Probability Surveys, 8:210–293.
  * Shao and Zhang, (2019) Shao, Q.-M. and Zhang, Z.-S. (2019).  Berry–Esseen bounds of normal and nonnormal approximation for unbounded exchangeable pair.  The Annals of Probability, 47(1):61–108.
  * Shevtsova, (2011) Shevtsova, I. (2011).  On the absolute constants in the berry-esseen type inequalities for identically distributed summands.  arXiv preprint arXiv:1111.6554.
  * Steiger, (1970) Steiger, W. (1970).  Bernstein’s inequality for martingales.  Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 16(2):104–106.
  * Stein, (1986) Stein, C. (1986).  Approximate computation of expectations.  IMS Lecture Notes—Monograph Series, 7.
  * Stone, (1974) Stone, M. (1974).  Cross-validatory choice and assessment of statistical predictions.  Journal of the royal statistical society: Series B (Methodological), 36(2):111–133.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267–288.
  * Tse and Davison, (2022) Tse, T. and Davison, A. C. (2022).  A note on universal inference.  Stat, 11(1):e501.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-dimensional statistics: A non-asymptotic viewpoint, volume 48.  Cambridge university press.
  * Wasserman et al., (2020) Wasserman, L., Ramdas, A., and Balakrishnan, S. (2020).  Universal inference.  Proceedings of the National Academy of Sciences, 117(29):16880–16890.
  * Zhang et al., (2023) Zhang, Z., Lee, S., and Dobriban, E. (2023).  A framework for statistical inference via randomized algorithms.  arXiv preprint arXiv:2307.11255.
  * Zhang, (2022) Zhang, Z.-S. (2022).  Berry–esseen bounds for generalized u-statistics.  Electronic Journal of Probability, 27:1–36.

Supplemental Appendix to: Reproducible Aggregation of Sample-Split
Statistics111Date:  David M. Ritzwoller Joseph P. Romano Stanford University
Stanford University

Contents
appendix.Asubsection.A.1subsection.A.2subsubsection.A.2.1subsubsection.A.2.2subsubsection.A.2.3subsubsection.A.2.4subsubsection.A.2.5subsection.A.3appendix.Bsubsection.B.1subsection.B.2subsection.B.3subsection.B.4subsection.B.5subsection.B.6subsection.B.7subsection.B.8subsection.B.9subsection.B.10appendix.Csubsection.C.1subsection.C.2subsection.C.3subsection.C.4subsection.C.5subsection.C.6subsection.C.7subsection.C.8appendix.Dsubsection.D.1subsection.D.2subsubsection.D.2.1subsection.D.3subsubsection.D.3.1

##  Appendix A Details Concerning Data and Simulations

###  A.1. Data

Banerjee et al., (2015) study randomized evaluations of several similar
poverty-alleviation programs implemented by BRAC, a large non-governmental
organization. The data from Banerjee et al., (2015) were acquired from
https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/NHIXNT
on September 10, 2021. See Appendix D of Chen and Ritzwoller, (2023) for
further details concerning the cleaning of these data. Our discussion in this
subsection uses language similar to the language used to discuss the data
there.

The programs studied by Banerjee et al., (2015) randomly allocated productive
assets (typically livestock) to participating households and measured economic
outcomes two and three years later. We restrict attention to data from the
program evaluation in Pakistan. In the sample considered in this paper, there
are 854 households. Of these households, 446 were randomly assigned to
participate in the program. In both the treatment effect and risk estimation
applications considered throughout the paper, the outcome of interest $Y_{i}$
is the total monthly consumption for each household. The covariate vector
$X_{i}$ collects measurements of 20 pretreatment variables. These variables
can be grouped into five categories: consumption, food security, assets,
finance, and income and revenue. The consumption variables are pretreatment
total monthly consumption and total monthly consumptions on food, non-food,
and durable commodities. The food-security variables are an index for overall
food-security and five binary variables indicating different aspects of food
security (e.g., did a child in the household skip a meal). The assets
variables are value of assets, value of productive assets, and value of the
household’s assets, each aggregated in two ways. The financial variables are
total amount of formal and informal loans outstanding and total value of the
household’s savings. The income and revenue variables are income from
agriculture, income from business, income from paid labor, revenue from
animals, and self-assessed perception of economic status. We refer the reader
to the appendix of Banerjee et al., (2015) for further information on the
construction of these variables.

###  A.2. Simulations

Throughout the paper, we consider three Lasso regressions

$\displaystyle\hat{\beta}_{w}(\lambda)$
$\displaystyle=\underset{\beta}{\arg\min}\left\\{\sum_{i\in\tilde{\mathsf{s}}}\mathbb{I}\\{W_{i}=w\\}(Y_{i}-\beta^{\top}X_{i})^{2}+\lambda\|\beta\|_{1}\right\\},\quad
w\in\\{0,1\\},\quad\text{and}$ (A.1) $\displaystyle\tilde{\beta}(\lambda)$
$\displaystyle=\underset{\beta}{\arg\min}\left\\{\sum_{i\in\tilde{\mathsf{s}}}(W_{i}-\beta^{\top}X_{i})^{2}+\lambda\|\beta\|_{1}\right\\}\
$ (A.2)

which produce the nuisance parameter estimates

$\hat{\mu}_{w}(X_{i})=\hat{\beta}_{w}(\lambda)^{\top}X_{i}\quad\text{and}\quad\hat{\pi}(X_{i})=\tilde{\beta}(\lambda)^{\top}X_{i}\leavevmode\nobreak\
.$

Estimates of the average treatment effect are formed with the estimator (1.2).
Figure A.1 displays average cross-validated risk estimates, over one million
cross-splits, for the regressions (A.1) and (A.2) as $k$ is varied. For all
figures, except for Figure 1, we choose $\lambda$ based off of the estimates
reported in Figure A.1. Specifically, for the two outcome regressions (A.1).
we choose the value of $\lambda$ that minimizes the average $10$-fold cross-
validated mean-squared error.121212The chosen values of $\lambda$ are equal to
approximately 3.38 and 1.85 for the untreated and treated outcome regression,
respectively. For the propensity score regression, the average $k$-fold mean-
squared error is monotonically decreasing in $\lambda$. To encourage some
selection of covariates, we set the value of $\lambda$ sub-optimally, at
approximately 0.02.

Figure A.1. Mean Squared Error for Nuisance Parameter Estimates

Panel A: Control Group Outcome Regression
---
Panel B: Treatment Group Outcome Regression
Panel C: Propensity Score

Notes: Figure A.1 displays the average of cross-validated risk estimate, taken
over one million cross-splits, for the regressions (A.1) and (A.2) implemented
with the Banerjee et al., (2015) data. Each line displays this estimate for a
different value of $k$. We consider the values of $\lambda$ reported by the
“glmnet” R software package (Friedman et al.,, 2021).

####  A.2.1. Figure 1

To emulate the way that cross-splitting is often applied in practice, for the
purpose of Panels A and B of Figure 1, we choose $\lambda$ by implementing
10-fold cross validation once. We use one hundred thousand cross-splits to
construct each Panel.

####  A.2.2. Figure 2

For the purpose of Figure 2, we reduce the sample under consideration so that
the total sample size has many perfect divisors. In particular, for the
treatment effect application, we randomly reduce the sample of 854 households
to a sample of 800 households. For the risk estimation application, we
randomly reduce the sample of 446 treated households to 400 treated
households. We consider the set of values of $b$ that perfectly divide 800 and
400 respectively. For each value of $b$, we generate one million sample
sample-splitsplits. In both applications, to address some numerical
instability in the estimation of $\phi_{n,b}(D)$ and $\gamma_{n,b}(D)$ we drop
the 0.25% of simulation draws with the largest and smallest values of
$T(\mathsf{s},D)$.

####  A.2.3. Figure 3

One million sample-splits are used to estimate $\sigma_{\mathsf{va\ell
id}}^{(2)}$ and $\sigma_{\mathsf{train}}^{(2,b-1)}$ for each value of $b$
considered in Figure 3. In particular, for a given value of $b$, we draw one
million sample-splits $\mathsf{s}$ and compute $T(\mathsf{s},D)$. For each
sample-split $\mathsf{s}$, we draw a random index $I$ in $\mathsf{s}$ and
compute

$\left(\psi(D_{I},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right))-\psi(D^{\prime}_{I},\hat{\eta}\left(D_{\tilde{\mathsf{s}}}\right)\right)^{2}\leavevmode\nobreak\
,$ (A.3)

where $D^{\prime}_{I}$ is a random data point drawn from the uniformly from
the data $D$.131313In the case of the risk estimation application to the
treated outcome regression, the data set $D$ is composed of only the data for
the treated households. The quantity $\sigma_{\mathsf{va\ell id}}^{(2)}$ is
estimated by averaging the quantity (A.3) over the one million sample-splits
$\mathsf{s}$. The estimate of $\sigma_{\mathsf{train}}^{(2,b-1)}$ is
constructed similarly, where now we replace $b-1$ random data points in
$\tilde{\mathsf{s}}$ with $b-1$ observations drawn uniformly from the data
$D$.

####  A.2.4. Figure 4

The total number of simulation draws used in the measurements displayed in
Figure 4 varies by the value of $k$. In the treatment effect estimation
application, one million cross-split estimates are drawn for each value of $k$
less than or equal to 20. For each value of $k$ greater than 20, five hundred
thousand cross-splits are used. Similarly, in the risk estimation application,
one million cross-splits are used for $k$ less than or equal to 20 and five
hundred thousand cross-splits are used for $k$ greater than 20. To address
issues with numerical instability, when computing the variance $v_{1,k}(D)$ in
the left-hand side of the Berry-Esseen inequality (3.11), we drop the 0.1% of
simulation draws with the largest and smallest values of
$a(\mathsf{R}_{1,k},D)$.

####  A.2.5. Figure 5

The measurements displayed in this figure are obtained by first computing one
million replicates of $a(\mathsf{R}_{1,k},D)$ for each value of $k$. For each
value of $k$ and and $\xi$, we implement Algorithm 1 twice fifty thousand
times by sampling with replacement from the replicates of
$a(\mathsf{R}_{1,k},D)$. The measurements of the reproducibility error give
the proportions of time that these two measurements differ by more than $\xi$.

###  A.3. Additional Figures and Discussion

Figure A.2 displays measurements of reproducibility error (4.1) when
$\mathsf{R}_{\hat{g},k}$ and $\mathsf{R}^{\prime}_{\hat{g}^{\prime},k}$ are
obtained independently with Algorithm 1 for the treatment effect application
to the Banerjee et al., (2015) data. Again, we set the nominal error rate to
$\beta=0.1$ and vary the values of $k$ and $\xi$. These measurements are
qualitatively similar to the results for the risk estimation application
displayed in (5).

Figure A.2. Reproducibility Error in Treatment Effect Estimation

Panel A: Early Stopping with $g_{\mathsf{init}}=2$
---
Panel B: Reproducibility with $g_{\mathsf{init}}=10$

Notes: Figure A.2 displays measurements of reproducibility error (4.1) for
Algorithm 1 in the treatment effect estimation application to the Banerjee et
al., (2015) data. A dashed horizontal line is displayed at the nominal error
rate $\beta=0.1$. The $x$-axes vary the bound $\xi$ and are displayed in a
log-scale base 10. To ease computation, we do not display reproducibility
error estimates for a $(\xi,k)$ pair if $k$ is less than 10 and $\xi$ is less
than $0.03$. We reduce the number of simulation replications from 50,000 to
10,000 for all estimates with $k$ less than or equal to $7$.

The details of the simulation are also similar. We first compute one million
replicates of $a(\mathsf{R}_{1,k},D)$ for each value of $k$. However, for
numerical stability, we remove the 0.001% of replicates with the largest and
smallest values of $a(\mathsf{R}_{1,k},D)$. These extreme values are caused by
the division by estimates of the propensity score in the formula (1.2). For
each value of $k$ and and $\xi$, we implement Algorithm 1 twice fifty thousand
times by sampling with replacement from the replicates of
$a(\mathsf{R}_{1,k},D)$. To ease computation, we do not display
reproducibility error estimates for a $(\xi,k)$ pair if $k$ is less than 10
and $\xi$ is less than $0.03$. We reduce the number of replications from fifty
thousand to ten thousand for all estimates with $k$ less than or equal to $7$.

Figure A.3 displays the quantiles of the distribution of the sample size
discrepancy

$\hat{g}/g^{\star}-1$ (A.4)

for the risk estimation application to the Banerjee et al., (2015) data. The
95th and 5th quantiles of the distribution of this difference are displayed
with triangles and dots, respectively. Panel A gives results for
$g_{\mathsf{init}}$ set equal to 2. The high probability of early stopping at
large values of $\xi$ is evident. Interestingly, there appears to be a phase
transition as $\xi$ decreases, wherein early stopping becomes much less
likely. Panel B gives results for $g_{\mathsf{init}}$ set equal to 10. Here,
the features that characterize the measurements in Panel A only occur at very
large values of $\xi$. That is, the problem of early stopping has largely been
resolved. At moderate and small values of $\xi$, the variance of $\hat{g}$
around $g^{\star}$ increases as $k$ increases.

Figure A.3. Sample Size Discrepancy in Risk Estimation

Panel A: Early Stopping with $g_{\mathsf{init}}=2$
---
Panel B: Reproducibility with $g_{\mathsf{init}}=10$

Notes: Figure A.3 displays quantiles of measurements of sample size
discrepancy (A.4) for Algorithm 1 in the risk estimation application to the
Banerjee et al., (2015) data. The 95th and 5th quantiles of the sample size
error are displayed with triangles and circles, respectively. We set the error
rate at $\beta=0.1$. The $x$-axes vary the bound $\xi$ and are displayed in a
log-scale base 10. We study the regression (1.1) for the treated outcome,
i.e., $w$ equal to one, for a fixed choice of $\lambda$. To ease computation,
we only display reproducibility error estimates for a $(\xi,k)$ pair if
$g^{\star}$ is less than $10^{4}$.

##  Appendix B Proofs for Results Stated in the Main Text

###  B.1. Proof of Theorem 2.1

For each collection $\mathsf{r}_{k}=(\mathsf{s}_{j})_{j=1}^{k}$ in
$\mathcal{R}_{n,k,b}$, we write

$\overline{a}(\mathsf{r}_{k},D)=\frac{1}{k}\sum_{j=1}^{k}T(\mathsf{s}_{j},D)-\mathbb{E}\left[T(\mathsf{s}_{j},D)\mid
D\right].$

Define the conditional variance

$v_{g,k}(D)=\operatorname{Var}(a(\mathsf{R}_{g,k},D)\mid D)$

and the oracle stopping time

$g^{\star}=\underset{g\geq
2}{\arg\min}\left\\{v_{g,k}(D)\leq\frac{1}{2}\left(\frac{\xi}{z_{1-\beta/2}}\right)^{2}\right\\}.$

Observe that $v_{g,k}(D)=(1/g)\cdot v_{1,k}(D)$ as the collections of cross-
splits are drawn independently and identically. Similarly, we have that

$\frac{\hat{v}(\mathsf{R}_{g,k},D)}{v_{g,k}(D)}=\frac{g\hat{v}(\mathsf{R}_{g,k},D)}{v_{1,k}(D)}\overset{\text{a.s}}{\to}1$
(B.1)

as $g\to\infty$, by the strong law of large numbers. Observe that

$\frac{\hat{g}\cdot\hat{v}(\mathsf{R}_{\hat{g},k},D)}{v_{1,k}(D)}\leq\frac{\hat{g}}{2v_{1,k}(D)}\left(\frac{\xi}{z_{1-\beta/2}}\right)^{2}\leq\frac{\hat{g}\cdot\hat{v}(\mathsf{R}_{\hat{g}-1,k},D)}{v_{1,k}(D)},$
(B.2)

by definition. Thus, as $\hat{g}\to\infty$ and

$g^{\star}\left(2v_{1,k}(D)\left(\frac{z_{1-\beta/2}}{\xi}\right)^{2}\right)\to
1$

as $\xi\to 0$, we have that

$\hat{g}/g^{\star}\overset{\text{a.s}}{\to}1$ (B.3)

as $\xi\to 0$ by (B.1) and (B.2).

Define the objects

$\displaystyle U\left(\mathsf{R}_{g,k},D\right)$
$\displaystyle=\frac{1}{g^{\star}}\left(\sum_{i=1}^{g}\overline{a}\left(\mathsf{r}_{i,k},D\right)-\sum_{i=1}^{g^{\star}}\overline{a}\left(\mathsf{r}_{i,k},D\right)\right)\quad\text{and}$
(B.4) $\displaystyle Q\left(\mathsf{R}_{g,k},D\right)$
$\displaystyle=\left(1-\frac{g}{g^{\star}}\right)\left(a(\mathsf{R}_{g,k},D)-\mathbb{E}\left[T(\mathsf{s}_{j},D)\mid
D\right]\right)$ (B.5)

Observe that we can decompose

$\displaystyle\left(a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\right)/\sqrt{2v_{g^{\star},k}(D)}$
$\displaystyle=\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(\frac{1}{g^{\star}}\sum_{i=1}^{g^{\star}}\left(\overline{a}\left(\mathsf{r}_{i,k},D\right)-\overline{a}\left(\mathsf{r}_{i,k}^{\prime},D\right)\right)\right)$
$\displaystyle+\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(\frac{1}{\hat{g}}\sum_{i=1}^{\hat{g}}\overline{a}\left(\mathsf{r}_{i,k},D\right)-\frac{1}{g^{\star}}\sum_{i=1}^{g^{\star}}\overline{a}\left(\mathsf{r}_{i,k},D\right)\right)$
$\displaystyle-\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(\frac{1}{\hat{g}^{\prime}}\sum_{i=1}^{\hat{g}^{\prime}}\overline{a}\left(\mathsf{r}_{i,k}^{\prime},D\right)-\frac{1}{g^{\star}}\sum_{i=1}^{g^{\star}}\overline{a}\left(\mathsf{r}_{i,k}^{\prime},D\right)\right)$
$\displaystyle=\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(a(\mathsf{R}_{g^{\star},k},D)-a(\mathsf{R}_{g^{\star},k}^{\prime},D)\right)+\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(U(\mathsf{R}_{\hat{g},k},D)-U(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\right)$
$\displaystyle+\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(Q(\mathsf{R}_{\hat{g},k},D)-Q(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\right).$
(B.6)

First, we have that

$\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(Q(\mathsf{R}_{\hat{g},k},D)-Q(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\right)=o_{p}(1)$

as $\xi\to 0$ by (B.3). To handle the terms involving (B.4), fix some
$\varepsilon>0$ and observe that

$\displaystyle
P\left\\{\bigg{|}\sum_{i=1}^{\hat{g}}a(\mathsf{r}_{i,k},D)-\sum_{i=1}^{g^{\star}}a(\mathsf{r}_{i,k},D)\bigg{|}>\varepsilon\sqrt{g^{\star}}\mid
D\right\\}$ $\displaystyle\leq
P\left\\{\bigg{|}\sum_{i=1}^{\hat{g}}a(\mathsf{r}_{i,k},D)-\sum_{i=1}^{g^{\star}}a(\mathsf{r}_{i,k},D)\bigg{|}>\varepsilon\sqrt{g^{\star}},\hat{g}\in\left[g^{\star}\left(1-\varepsilon^{3}\right),g^{\star}\left(1+\varepsilon^{3}\right)\right]\mid
D\right\\}$
$\displaystyle+P\left\\{\bigg{|}\sum_{i=1}^{\hat{g}}a(\mathsf{r}_{i,k},D)-\sum_{i=1}^{g^{\star}}a(\mathsf{r}_{i,k},D)\bigg{|}>\varepsilon\sqrt{g^{\star}},\hat{g}^{\prime}\not\in\left[g^{\star}\left(1-\varepsilon^{3}\right),g^{\star}\left(1+\varepsilon^{3}\right)\right]\mid
D\right\\}$ $\displaystyle\leq
P\left\\{\max_{g^{\star}\left(1-\varepsilon^{3}\right)\leq g\leq
g^{\star}}\bigg{|}\sum_{i=1}^{g}a(\mathsf{r}_{i,k},D)\bigg{|}>\varepsilon\sqrt{g^{\star}}\mid
D\right\\}$ (B.7) $\displaystyle+P\left\\{\max_{g^{\star}\leq
g\leq\left(1+\varepsilon^{3}\right)g^{\star}}\bigg{|}\sum_{i=1}^{g}a(\mathsf{r}_{i,k},D)\bigg{|}>\varepsilon\sqrt{g^{\star}}\mid
D\right\\}+P\left\\{\hat{g}\not\in\left[g^{\star}\left(1-\varepsilon^{3}\right),g^{\star}\left(1+\varepsilon^{3}\right)\right]\mid
D\right\\}\leavevmode\nobreak\ .$

The third term in (B.7) is smaller than $\varepsilon$ for all sufficiently
small $\xi$ by (B.3). To handle the first two terms, observe that

$\displaystyle P\left\\{\max_{g^{\star}\leq
g\leq\left(1+\varepsilon^{3}\right)g^{\star}}\bigg{|}\sum_{i=1}^{g}a(\mathsf{r}_{i,k},D)\bigg{|}>\varepsilon\sqrt{\hat{g}}\mid
D\right\\}\leq\frac{1}{\varepsilon^{2}}\frac{1}{g^{\star}}\operatorname{Var}\left(\sum_{i=1}^{\lfloor\varepsilon^{3}\hat{g}\rfloor+1}a(\mathsf{r}_{i,k},D)\mid
D\right)\leq\varepsilon,$

where the first inequality follows from Kolmogorov’s maximal inequality (see
e.g., Proposition 2.3.16 of Dembo, (2021)). An analogous bound holds for the
remaining term. Thus, we have that

$\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(U(\mathsf{R}_{\hat{g},k},D)-U(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\right)=o_{p}\left(1\right)$

as $\xi\to 0$. Hence, we have that

$\displaystyle
P\left\\{\big{|}a(\mathsf{R}_{\hat{g},k},D)-a(\mathsf{R}_{\hat{g}^{\prime},k}^{\prime},D)\big{|}>\xi\mid
D\right\\}$
$\displaystyle=P\left\\{\bigg{|}\sqrt{\frac{g^{\star}}{2v_{1,k}(D)}}\left(\frac{1}{g^{\star}}\sum_{i=1}^{g^{\star}}a(\mathsf{r}_{i,k},D)-a(\mathsf{r}_{i,k}^{\prime},D)\right)\bigg{|}>z_{1-\beta/2}\mid
D\right\\}+o\left(1\right)$ $\displaystyle=1-\beta+o\left(1\right)$

as $\xi\to 0$ by the central limit theorem.

###  B.2. Proof of Lemma 5.1

We use the following Lemma in several places.

###### Lemma B.1.

Let $\bm{\psi}$ and $\bm{\psi}^{\prime}$ be two sets, each containing $g$
elements of $\mathcal{P}_{n}$. The inequality

$\displaystyle\sum_{m=0}^{\infty}\bigg{|}\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m})\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}^{\prime}_{m})\right)\mid\bm{\pi}_{0}=\bm{\psi},\bm{\pi}^{\prime}_{0}=\bm{\psi}^{\prime},D\right]\bigg{|}$
$\displaystyle\leq
2gn^{2}\max_{\mathsf{s},\mathsf{s}^{\prime}\in\mathsf{S}_{n,b}}\left(T\left(\mathsf{s},D\right)-T\left(\mathsf{s}^{\prime},D\right)\right)$

holds almost surely.

Observe that the quantity

$\max_{\mathsf{s},\mathsf{s}^{\prime}\in\mathsf{S}_{n,b}}\left(T\left(\mathsf{s},D\right)-T\left(\mathsf{s}^{\prime},D\right)\right)$

is finite almost surely. Thus, the convergence of the series defining
$A\left(\bm{\pi},\bm{\pi}^{\prime}\mid D\right)$ follows from Lemma B.1.
Define the operator

$\displaystyle K:\mathcal{F}$ $\displaystyle\to\mathcal{F}$ $\displaystyle
f\left(\cdot\right)$
$\displaystyle\mapsto\mathbb{E}\left[f\left(\bm{\pi}^{\prime}\right)\mid\bm{\pi}=\cdot\right],$

where $\mathcal{F}$ is the set of all measurable functions supported on the
domain of $\bm{\pi}$. Observe that

$\displaystyle\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}^{\prime}_{m}),D\right)\mid\bm{\pi}_{0}=\bm{\pi},\bm{\pi}^{\prime}_{0}=\bm{\pi}^{\prime},D\right]$
$\displaystyle=\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-\overline{a}(D)\mid\bm{\pi}_{0}=\bm{\pi},D\right]-\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}^{\prime}),D\right)-\overline{a}(D)\mid\bm{\pi}^{\prime}_{0}=\bm{\pi}^{\prime},D\right]$
$\displaystyle=K^{m}\left(a\left(\mathsf{R}_{g,k}(\bm{\pi})\right)-\overline{a}(D)\right)-K^{m+1}\left(a\left(\mathsf{R}_{g,k}(\bm{\pi}^{\prime})\right)-\overline{a}(D)\right).$

Thus, for any $m^{\prime}$, we have that

$\displaystyle\sum_{m=0}^{m^{\prime}}\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}^{\prime}_{m}),D\right)\mid\bm{\pi}_{0}=\bm{\pi},\bm{\pi}^{\prime}_{0}=\bm{\pi}^{\prime},D\right]$
$\displaystyle=a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)-K^{m^{\prime}+1}\left(a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)\right).$
(B.8)

By Lemma B.1, the partial sums (B.8) converge almost everywhere and so the
sequence

$\left(K^{m+1}\left(a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)\right)\right)_{m\geq
0}$ (B.9)

also converges almost everywhere. Lemma B.1 also implies that the limit of
(B.9) depends only on $D$, as

$K^{m}\left(a\left(\left(\mathsf{R}_{g,k}\left(\bm{\psi}\right),D\right)\right)-\overline{a}\left(D\right)\right)-K^{m}a\left(\left(\mathsf{R}_{g,k}\left(\bm{\psi}^{\prime}\right),D\right)-\overline{a}\left(D\right)\right)\to
0$

for any $\bm{\psi}$ and $\bm{\psi}^{\prime}$ each containing $g$ elements of
$\mathcal{P}_{n}$. Therefore, we have that

$\displaystyle\mathbb{E}\left[A\left(\bm{\pi},\bm{\pi}^{\prime}\mid
D\right)\mid D\right]$
$\displaystyle=\mathbb{E}\left[\lim_{n\to\infty}\left(a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)-K^{m+1}a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)\right)\mid
D\right]$
$\displaystyle=\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)\mid
D\right]-b\left(D\right),$

for some quantity

$b\left(D\right)=\lim_{m\to\infty}K^{m}\left(a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)\right)$

that depends only on $D$. Observe that

$\displaystyle\mathbb{E}\left[A\left(\bm{\pi},\bm{\pi}^{\prime}\mid
D\right)\mid D\right]$
$\displaystyle=\mathbb{E}\left[\lim_{m^{\prime}\to\infty}\sum_{m=0}^{m^{\prime}}\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}^{\prime}),D\right)\mid\bm{\pi}_{0}=\bm{\pi},\bm{\pi}^{\prime}_{0}=\bm{\pi}^{\prime},D\right]\mid
D\right]$
$\displaystyle=\lim_{m^{\prime}\to\infty}\sum_{m=0}^{m^{\prime}}\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}^{\prime}),D\right)\mid
D\right]$ (Dominated Conv.) $\displaystyle=0,$ (Exchangeability)

where the applicability of the Dominated Convergence Theorem follows from
Lemma B.1. Thus, as

$\displaystyle\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)\mid
D\right]$ $\displaystyle=0,$

we can conclude that $b\left(D\right)=0$ almost surely. Hence, we find that

$\displaystyle\mathbb{E}\left[A\left(\bm{\pi},\bm{\pi}^{\prime}\mid
D\right)\mid\bm{\pi},D\right]$
$\displaystyle=\mathbb{E}\left[\lim_{m^{\prime}\to\infty}\sum_{m=0}^{m^{\prime}}\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}^{\prime}),D\right)\mid\bm{\pi}_{0}=\bm{\pi},\bm{\pi}^{\prime}_{0}=\bm{\pi}^{\prime},D\right]\mid\bm{\pi},D\right]$
$\displaystyle=\lim_{m^{\prime}\to\infty}\sum_{m=0}^{m^{\prime}}\mathbb{E}\left[a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}),D\right)-a\left(\mathsf{R}_{g,k}(\bm{\pi}_{m}^{\prime}),D\right)\mid\bm{\pi}_{0}=\bm{\pi},\bm{\pi}^{\prime}_{0}=\bm{\pi}^{\prime},D\right]$
(Dominated Conv.)
$\displaystyle-\lim_{m^{\prime}\to\infty}K^{m^{\prime}+1}\left(a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right)\right)$
$\displaystyle=a\left(\mathsf{R}_{g,k}(\bm{\pi}),D\right)-\overline{a}\left(D\right),$

completing the proof.

###  B.3. Proof of Lemma 5.2

First, observe that

$U_{F}\left(X\right)\geq\frac{1}{2}\left(\mathbb{E}\left[F\left(X,X^{\prime}\right)\mid
X\right]\right)^{2}=\frac{1}{2}f\left(X\right)^{2},$

where the inequality follows from Jensen’s inequality and the definition of
the Stein representer $F$. By (5.11), we have that

$su\geq\frac{1}{2}f\left(X\right)^{2}.$

Hence, the random variable $f\left(X\right)$ is bounded almost surely.

Now, suppose $h:\mathcal{X}\to\mathbb{R}$ is any measurable function such that
$\mathbb{E}\left[h\left(X\right)F\left(X,X^{\prime}\right)\right]<\infty$.
Then,
$\mathbb{E}\left[h\left(X\right)f\left(X\right)\right]=\mathbb{E}\left[h\left(X\right)F\left(X,X^{\prime}\right)\right]$.
Using the exchangeability of $X$ and $X^{\prime}$ and the fact that $F$ is
antisymmetric, we have that

$\mathbb{E}\left[h\left(X\right)F\left(X,X^{\prime}\right)\right]=\mathbb{E}\left[h\left(X^{\prime}\right)F\left(X^{\prime},X\right)\right]=-\mathbb{E}\left[h\left(X^{\prime}\right)F\left(X,X^{\prime}\right)\right]$

and that therefore

$\mathbb{E}\left[h\left(X\right)f\left(X\right)\right]=\frac{1}{2}\mathbb{E}\left[\left(h\left(X\right)-h\left(X^{\prime}\right)\right)F\left(X,X^{\prime}\right)\right].$
(B.10)

Let

$m\left(\theta\right)=\mathbb{E}\left[\exp\left(\theta
f\left(X\right)\right)\right]$

denote the moment generating function of $f\left(X\right)$. As
$f\left(X\right)$ is bounded almost surely, we can exchange differentiation
and expectation in the differentiation of $m\left(\theta\right)$. Thus, we
obtain

$\displaystyle m^{\prime}\left(\theta\right)$
$\displaystyle=\mathbb{E}\left[\exp\left(\theta
f\left(X\right)\right)f\left(X\right)\right].$
$\displaystyle=\frac{1}{2}\mathbb{E}\left[\left(\exp\left(\theta
f\left(X\right)\right)-\exp\left(\theta
f\left(X^{\prime}\right)\right)\right)F\left(X,X^{\prime}\right)\right],$
(B.11)

where the second inequality follows from (B.10). To bound
$m^{\prime}\left(\theta\right)$ we apply the following exponential mean-value
inequality, stated in a more general form in Paulin et al., (2016).

###### Lemma B.2.

For all constants $x$, $y$, and $c$ in $\mathbb{R}$ and $s>0$, it holds that

$|\left(e^{x}-e^{y}\right)c|\leq\frac{1}{4}\left(s\left(x-y\right)^{2}+s^{-1}c^{2}\right)\left(e^{x}+e^{y}\right).$

In particular, by (B.11) and Lemma B.2, we obtain the bound

$\displaystyle|m^{\prime}\left(\theta\right)|$
$\displaystyle\leq\frac{1}{2}\mathbb{E}\left[|\left(\exp\left(\theta
f\left(X\right)\right)-\exp\left(\theta
f\left(X^{\prime}\right)\right)\right)F\left(X,X^{\prime}\right)|\right]$
$\displaystyle\leq\frac{1}{8}\inf_{t>0}\mathbb{E}\left[\left(t\left(\theta
f\left(X\right)-\theta
f\left(X^{\prime}\right)\right)^{2}+t^{-1}F\left(X,X^{\prime}\right)^{2}\right)\left(\exp\left(\theta
f\left(X\right)\right)+\exp\left(\theta
f\left(X^{\prime}\right)\right)\right)\right]$
$\displaystyle=\frac{|\theta|}{4}\inf_{t>0}\mathbb{E}\left[\left(t\left(f\left(X\right)-f\left(X^{\prime}\right)\right)^{2}+t^{-1}F\left(X,X^{\prime}\right)^{2}\right)\exp\left(\theta
f\left(X\right)\right)\right]$
$\displaystyle=\frac{|\theta|}{2}\inf_{t>0}\mathbb{E}\left[\left(\frac{t}{2}\mathbb{E}\left[\left(f\left(X\right)-f\left(X^{\prime}\right)\right)^{2}\mid
X\right]+\frac{1}{2t}\mathbb{E}\left[F\left(X,X^{\prime}\right)^{2}\mid
X\right]\right)\mathbb{E}\left[\exp\left(\theta f\left(X\right)\right)\mid
X\right]\right].$
$\displaystyle=\frac{|\theta|}{2}\inf_{t>0}\mathbb{E}\left[\left(tU_{f}\left(X\right)+t^{-1}U_{F}\left(X\right)\right)\mathbb{E}\left[\exp\left(\theta
f\left(X\right)\right)\mid X\right]\right]$
$\displaystyle\leq\frac{|\theta|}{2}\mathbb{E}\left[\left(sU_{f}\left(X\right)+s^{-1}U_{F}\left(X\right)\right)\mathbb{E}\left[\exp\left(\theta
f\left(X\right)\right)\mid X\right]\right]$
$\displaystyle\leq|\theta|v\mathbb{E}\left[\exp\left(\theta
f\left(X\right)\right)\right]$

for all $\theta\in\mathbb{R}$. Thus, we have that

$m^{\prime}\left(\theta\right)\leq u\theta m\left(\theta\right)$

for all $\theta>0$. As $m\left(\cdot\right)$ is a convex function and
$m^{\prime}\left(0\right)=0$, $m^{\prime}\left(\theta\right)$ always has the
same sign as $\theta$, we find that

$\frac{\text{d}}{\text{d}\theta}\log m\left(\theta\right)\leq u\theta.$

As a consequence, and by $m\left(0\right)=1$, we have that

$\log
m\left(\theta\right)\leq\int_{0}^{\theta}ut\text{d}t\leq\frac{u\theta^{2}}{2}.$

By the Chernoff bound (see e.g., Wainwright,, 2019, Equation 2.5, ), we have

$\log P\left\\{f\left(X\right)\geq\delta\right\\}\leq\inf_{\theta\geq
0}\left(\log m\left(\theta\right)-\theta\delta\right)\leq\inf_{\theta\geq
0}\left(\frac{u\theta^{2}}{2}-\theta\delta\right)$

Solving this minimization with $\theta=\delta/u$, we find

$P\left\\{f\left(X\right)\geq
t\right\\}\leq\exp\left(\frac{-\delta^{2}}{2u}\right),$

as required. The analogous lower tail bound follows from an identical
argument, which completes the proof of (5.7). To prove (5.8) we apply the
following result stated in Chatterjee, (2007).

###### Theorem B.1 (Theorem 1.5, (iii), Chatterjee,, 2007).

Reintroduce the notation and assumptions from the statement of Theorem 5.2.
Define

$\Delta\left(X\right)=\frac{1}{2}\mathbb{E}\left[|F\left(X,X^{\prime}\right)\left(f\left(X\right)-f\left(X^{\prime}\right)\right)|\mid
X\right].$

The Burkholder-Davis-Gundy inequality

$\mathbb{E}\left[f\left(X\right)^{2r}\right]\leq\left(2r-1\right)^{r}\mathbb{E}\left[\Delta\left(X\right)^{r}\right]$

holds for any positive integer $r$.

The inequality

$\displaystyle\mathbb{E}\left[\Delta\left(X\right)^{r}\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[|\left(f\left(X\right)-f\left(X^{\prime}\right)\right)F\left(X,X^{\prime}\right)|\mid
X\right]^{r}\right]$
$\displaystyle\leq\mathbb{E}\left[\mathbb{E}\left[|\left(f\left(X\right)-f\left(X^{\prime}\right)\right)F\left(X,X^{\prime}\right)|^{r}\mid
X\right]\right]$ (Jensen)
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[\left(\left(s^{-1}\left(f\left(X\right)-f\left(X^{\prime}\right)\right)^{2r}\right)\left(sF\left(X,X^{\prime}\right)^{2r}\right)\right)^{1/2}\mid
X\right]\right]$
$\displaystyle\leq\mathbb{E}\left[\mathbb{E}\left[s^{-1}\left(f\left(X\right)-f\left(X^{\prime}\right)\right)^{2r}+sF\left(X,X^{\prime}\right)^{2r}\mid
X\right]\right]$ (Young)
$\displaystyle=s^{-1}\mathbb{E}\left[{U}_{f}^{(r)}(X)\right]+s\mathbb{E}\left[{U}_{F}^{(r)}(X)\right]$

then completes the proof.

###  B.4. Proof of Lemma 5.3

We apply the following Lemma.

###### Lemma B.3.

Let $\left(X_{m}\right)_{m\geq 0}$ be a sequence of real-valued random
variables. Suppose that the inequality

$\mathbb{E}\left[X_{m}^{2^{c}}\right]\leq h_{m}^{2^{c}}$ (B.12)

holds for each $m\geq 0$ and positive integer $c$, where
$\left(h_{m}\right)_{m\geq 0}$ is a deterministic sequence of nonnegative real
numbers. If there exists a square integrable random variable $W$ such that

$\left(\sum_{m=0}^{\infty}X_{m}\right)^{2^{c}}\leq W$

almost surely, then the inequality

$\mathbb{E}\left[\left(\sum_{m=0}^{\infty}X_{m}\right)^{2^{c}}\right]\leq\left(\sum_{m=0}^{\infty}h_{m}\right)^{2^{c}}$

holds almost surely.

Define the objects

$\displaystyle
U_{f}^{(r)}\left(Z\right)=\frac{1}{2}\mathbb{E}\left[\left(f\left(Z\right)-f\left(Z^{\prime}\right)\right)^{2r}\mid
Z\right]$ and $\displaystyle
U_{F}^{(r)}\left(Z\right)=\frac{1}{2}\mathbb{E}\left[\left(\sum_{m=0}^{\infty}\mathbb{E}\left[f\left(Z_{m}\right)-f\left(Z_{m}^{\prime}\right)\mid
Z_{0}=Z,Z_{0}^{\prime}=Z^{\prime}\right]\right)^{2r}\mid Z\right].$

By the (5.9), we have that

$\mathbb{E}\left[\sum_{m=0}^{\infty}f(Z_{m})-f(Z_{m}^{\prime})\mid
Z_{0}=Z,Z_{0}^{\prime}=Z^{\prime}\right]^{2^{c}}\leq W.$

Thus, (5.10) guarantees that the conditions of Lemma B.3 are satisfied, and we
have that

$\mathbb{E}\left[U_{F}^{(2^{c-1})}\left(Z\right)\mid\bm{\pi}\right]\leq\frac{1}{2}\left(\sum_{i=0}^{\infty}h_{i}\right)^{2^{c}}\quad\text{and}\quad\mathbb{E}\left[U_{f}^{(2^{c-1})}\left(Z\right)\mid\bm{\pi}\right]\leq\frac{1}{2}h_{0}^{2^{c}}.$

By Markov’s inequality, we obtain

$P\left\\{U_{F}^{(2^{c-1})}\left(Z\right)\geq\frac{2}{\delta}\mathbb{E}\left[U_{F}^{(2^{c-1})}\left(Z\right)\mid\bm{\pi}\right]\quad\text{or}\quad
U_{f}^{(2^{c-1})}\left(Z\right)\geq\frac{2}{\delta}\mathbb{E}\left[U_{f}^{(2^{c-1})}\left(Z\right)\mid\bm{\pi}\right]\mid\bm{\pi}\right\\}\leq\delta.$

Hence, by Lemma B.3 and DeMorgan’s law, the probability that both

$\displaystyle U_{F}^{(2^{c-1})}\left(Z\right)$
$\displaystyle\leq\frac{2}{\delta}\mathbb{E}\left[U_{F}^{(2^{c-1})}\left(Z\right)\mid\bm{\pi}\right]\leq\frac{1}{\delta}\left(\sum_{m=0}^{\infty}h_{m}\right)^{2^{c}}\quad\text{and}$
$\displaystyle U_{f}^{(2^{c-1})}\left(Z\right)$
$\displaystyle\leq\frac{2}{\delta}\mathbb{E}\left[U_{f}^{(2^{c-1})}\left(Z\right)\mid\bm{\pi}\right]\leq\frac{h_{0}^{2^{c}}}{\delta}$

hold is greater than $1-\delta$.

###  B.5. Proof of Lemma 5.4

Recall the definition of the collections
$\left(\bm{\pi}_{m},\bm{\pi}_{m}^{\prime}\right)_{m\geq 0}$ given in (5.3).
Fix $Z_{0}=Z$ and $Z_{0}^{\prime}=Z^{\prime}$ throughout. For any $i$ in
$\left[n\right]$, let $\mathsf{s}_{m,\ell}\left(i\right)$ denote the element
of the collection

$\mathsf{r}(\pi_{m,\ell})=\left(\mathsf{s}_{1}(\pi_{m,\ell}),\ldots,\mathsf{s}_{g}(\pi_{m,\ell})\right)$

that contains the index $i$ and set $\mathsf{s}_{m,\ell}(i)$ equal to
$\varnothing$ if no element of $\mathsf{r}(\pi_{m,\ell})$ contains $i$.

We begin by defining three events that will determine the structure of our
argument. By construction, the collections $\bm{\pi}_{m}$ and
$\bm{\pi}_{m}^{\prime}$ are either identical or differ in exactly two indices
in their $L$th element. Let $\mathcal{E}_{m}$ denote the event that
$\bm{\pi}_{m}$ and $\bm{\pi}_{m}^{\prime}$ differ. On the event
$\mathcal{E}_{m}$, let $i_{1,m}$ and $i_{2,m}$ denote the two indices in which
the $L$th elements of $\bm{\pi}_{m}$ and $\bm{\pi}_{m}^{\prime}$ differ.
Define the random variables $B_{m}$ and $C_{m}$ such that, conditional on
$\mathcal{E}_{m}$, each is uniformly distributed on
$\left\\{i_{1,m},i_{2,m}\right\\}$ such that $B_{m}\neq C_{m}$. On the
complement of $\mathcal{E}_{m}$, set these indices uniformly at random. Let
$\mathcal{F}_{m}$ denote the event that

$\mathsf{s}_{m,L}\left(B_{m}\right)\neq\mathsf{s}_{m,L}\left(C_{m}\right),$

i.e., the event that the indices that differ are not in the same element of
the collection $\mathsf{r}(\pi_{m,\ell})$. Finally, let $\mathcal{G}_{m}$
denote the event that

$\mathsf{s}_{m,L}\left(B_{m}\right)\neq\varnothing\quad\text{and}\quad\mathsf{s}_{m,L}\left(C_{m}\right)\neq\varnothing,$

i.e., the event that the indices that differ are both in the collection
$\mathsf{r}(\pi_{m,\ell})$.

By Assumption 3.1, the event
$\mathcal{H}_{m}=\mathcal{E}_{m}\cap\mathcal{F}_{m}$ is a necessary condition
for $a\left(Z_{m}\right)-a\left(Z_{m}^{\prime}\right)\neq 0$. Thus, we can
compute

$\displaystyle P\left\\{\mathcal{H}_{m}\mid\mathcal{E}_{0}\right\\}$
$\displaystyle=P\left\\{\mathcal{F}_{m}\mid\mathcal{E}_{m},\mathcal{E}_{0}\right\\}P\left\\{\mathcal{E}_{m}\mid\mathcal{E}_{0}\right\\}$
$\displaystyle=\Big{(}P\left\\{\mathsf{s}_{m,L}\left(B_{m}\right)\neq\mathsf{s}_{m,L}\left(C_{m}\right)\mid\mathsf{s}_{m,L}\left(B_{m}\right)\neq\varnothing,\mathcal{E}_{m}\right\\}P\left\\{\mathsf{s}_{m,L}\left(B_{m}\right)\neq\varnothing\mid\mathcal{E}_{m}\right\\}$
$\displaystyle\quad+P\left\\{\mathsf{s}_{m,L}\left(B_{m}\right)\neq\mathsf{s}_{m,L}\left(C_{m}\right)\mid\mathsf{s}_{m,L}\left(B_{m}\right)=\varnothing,\mathcal{E}_{m}\right\\}P\left\\{\mathsf{s}_{m,L}\left(B_{m}\right)=\varnothing\mid\mathcal{E}_{m}\right\\}\Big{)}$
$\displaystyle\quad\quad\cdot
P\left\\{\mathcal{E}_{m}\mid\mathcal{E}_{0}\right\\}$
$\displaystyle=\left(\frac{n-b}{n-1}\frac{kb}{n}+\frac{kb}{n-1}\frac{n-kb}{b}\right)\left(1-\frac{2}{gn^{2}}\right)^{m}$
$\displaystyle=\left(\frac{kb\left(2n-kb-b\right)}{n(n-1)}\right)\left(1-\frac{2}{gn^{2}}\right)^{m}$

for all $m\geq 0$ by the law of total probability. Moreover, we have that

$P\left\\{\mathcal{H}_{m}\mid Z,Z^{\prime}\right\\}\leq
P\left\\{\mathcal{H}_{m}\mid\mathcal{E}_{0}\right\\}$

almost surely. Consequently, we find that

$\displaystyle\mathbb{E}\left[\mathbb{E}\left[a\left(Z_{m}\right)-a\left(Z_{m}^{\prime}\right)\mid
Z,Z^{\prime}\right]^{2r}\mid\bm{\pi}\right]$
$\displaystyle=\mathbb{E}\left[\left(P\left\\{\mathcal{H}_{m}\mid
Z,Z^{\prime}\right\\}\mathbb{E}\left[a\left(Z_{m}\right)-a\left(Z_{m}^{\prime}\right)\mid
Z,Z^{\prime},\mathcal{H}_{m}\right]\right)^{2r}\mid\bm{\pi}\right]$
$\displaystyle\leq\left(P\left\\{\mathcal{H}_{m}\mid\mathcal{E}_{0}\right\\}\right)^{2r}\mathbb{E}\left[\mathbb{E}\left[\left(a\left(Z_{m}\right)-a\left(Z_{m}^{\prime}\right)\right)^{2r}\mid
Z,Z^{\prime},\mathcal{H}_{m}\right]\mid\bm{\pi}\right]$ (Jensen)
$\displaystyle=\left(1-\frac{2}{gn^{2}}\right)^{2mr}\left(\frac{kb\left(2n-kb-b\right)}{n(n-1)}\right)^{2r}\mathbb{E}\left[\mathbb{E}\left[\left(a\left(Z_{m}\right)-a\left(Z_{m}^{\prime}\right)\right)^{2r}\mid\mathcal{H}_{m},\bm{\pi}\right]\mid\bm{\pi}\right]$
$\displaystyle\leq\left(1-\frac{2}{gn^{2}}\right)^{2mr}\left(\frac{2kb\left(2n-kb-b\right)}{n^{2}}\right)^{2r}\mathbb{E}\left[\left(a\left(Z_{m}\right)-a\left(Z_{m}^{\prime}\right)\right)^{2r}\mid\mathcal{H}_{m}\right]\leavevmode\nobreak\
,$ (B.13)

where the final inequality follows from the fact that $\bm{\pi}$ is uniformly
distributed independently of $D$ and the elementary inequality $n/(n-1)\leq
2$. Thus, it remains to bound the expectation in (B.13).

To ease notation, we now drop the dependence on $m$ and $L$. Observe that

$\displaystyle P\left\\{\mathcal{G}\mid\mathcal{H}\right\\}$
$\displaystyle=P\left\\{\mathsf{s}\left(B\right)\neq\varnothing\mid\mathsf{s}\left(C\right)\neq\varnothing\right\\}=\frac{kb-b}{n-b}$

and that therefore

$\displaystyle\mathbb{E}\left[\left(a\left(Z\right)-a\left(Z^{\prime}\right)\right)^{2r}\mid\mathcal{H},\bm{\pi}\right]$
$\displaystyle=\left(\frac{kb-b}{n-b}\right)\mathbb{E}\left[\left(a\left(Z\right)-a\left(Z^{\prime}\right)\right)^{2r}\mid\mathcal{\mathcal{G}}\right]$