# Rank-transformed subsampling:
Inference for multiple data splitting and exchangeable p-values

F. Richard Guo<EMAIL_ADDRESS>Statistical Laboratory, University of
Cambridge, UK  Rajen D. Shah<EMAIL_ADDRESS>Statistical Laboratory,
University of Cambridge, UK

###### Abstract

Many testing problems are readily amenable to randomised tests such as those
employing data splitting, which divide the data into disjoint parts for
separate purposes. However despite their usefulness in principle, randomised
tests have obvious drawbacks. Firstly, two analyses of the same dataset may
lead to different results. Secondly, the test typically loses power because it
does not fully utilise the entire sample. As a remedy to these drawbacks, we
study how to combine the test statistics or p-values resulting from multiple
random realisations such as through random data splits. We introduce rank-
transformed subsampling as a general method for delivering large sample
inference about the combined statistic or p-value under mild assumptions. We
apply our methodology to a range of problems, including testing unimodality in
high-dimensional data, testing goodness-of-fit of parametric quantile
regression models, testing no direct effect in a sequentially randomised trial
and calibrating cross-fit double machine learning confidence intervals. For
the latter, our method improves coverage in finite samples and for the testing
problems, our method is able to derandomise and improve power. Moreover, in
contrast to existing p-value aggregation schemes that can be highly
conservative, our method enjoys type-I error control that asymptotically
approaches the nominal level.

Keywords—Data-splitting; Cross-fitting; Goodness-of-fit; Rejection sampling;
Subsampling; Unimodality; Verma constraint.

##  1 Introduction

Many modern statistical procedures are randomised in the sense that the output
is a random function of the data. A prominent class of such procedures are
hypothesis testing methods that involve splitting the dataset into independent
parts (Moran, 1973; Cox, 1975). These procedures randomly divide the data into
two non-overlapping sample sets, A and B, and then perform two steps which can
be described as “hunt and test”: first, sample A is used to choose one from
among a number of test statistics; next the chosen statistic is applied to
sample B to produce the final test statistic. This approach is attractive
because in the first stage an arbitrarily complicated procedure may be
employed to “hunt” for an appropriate test. Clearly, were we to ignore the
fact that our test statistic was selected from data and simply apply it to
sample A, we would fail to control the type-I error, a phenomenon sometimes
referred to as “double dipping” or “data snooping”. However, because the data
in A and B are independent, we can, in the second stage, effectively forget
that the test statistic was chosen in a data-driven way, which permits
straightforward calibration even when using a complicated “hunting” procedure.
This strategy is particularly useful in settings with complex alternatives, as
the test may be chosen to target the particular alternative from which the
data appear to have arisen. This approach has been used for a variety of
problems, such as testing the location of multiple samples (Cox, 1975),
constructing conformal prediction intervals (Lei et al., 2018; Solari and
Djordjilović, 2022), goodness-of-fit testing (Janková et al., 2020),
conditional (mean) independence testing (Scheidegger et al., 2021; Lundborg et
al., 2022), and conducting inference that is agnostic to the asymptotic regime
(Kim and Ramdas, 2020), to list just a few. As we show in this work (see
Section 4.1.2), a hunt and test approach can also be used to test for a
clustering structure, i.e., for testing the null of unimodality, in high-
dimensional data.

Another use of data splitting is related to nonparametric or semiparametric
methods where the estimator for the parameter of interest depends on nuisance
parameters that must also be estimated. To ensure proper asymptotic behaviour
of the final estimator, the bias from nuisance parameter estimation needs to
be controlled, and this may be achieved by employing a form of sample
splitting known as cross-fitting. Here, the data are first split into folds of
roughly equal size, and then estimators are computed on each fold using
nuisance parameters estimated from out-of-fold data. The per-fold estimates
are then combined to produce the final estimate. The independence afforded by
data splitting permits the use of flexible machine learning methods to
estimate these parameters, as adopted by targeted estimation or
double/debiased machine learning methods. This second use of sample splitting
has recently become popular in practice, although the idea has a long history;
see also Newey and Robins (2018); Chernozhukov et al. (2018) and references
therein.

Randomised procedures also arise naturally in settings where a null hypothesis
one wishes to test may more easily be stated in terms of a re-weighted
distribution. Thams et al. (2021) show how a wide range of problems may be
cast in this framework, including testing properties of a new policy in a
contextual bandits setting, model selection after covariate shift and testing
so-called generalised conditional independencies (also known as dormant
independencies, or Verma constraints; see Section 4.2 for further discussion).
In this work we focus on the latter, which after appropriate reweighting may
be reduced to simpler independence testing problems. Thams et al. (2021)
further propose to use resampling or rejection sampling to select from the
original set of observations, a random subset that behaves like a sample from
the reweighted distribution, to which an off-the-shelf testing procedure may
then be applied.

Despite its simplicity and broad applicability, as pointed out by Cox (1975),
randomised procedures have obvious drawbacks. Firstly, the extra randomness
hinders replication of the analysis, an issue that is particularly concerning
in view of today’s “replication crisis” in many scientific disciplines
(Ioannidis, 2005; Open Science Collaboration, 2015; Baker, 2016). Although one
may insist that the random seed used in a randomised algorithm should be part
of the replication, when the substantive result of an analysis hinges on a
particular seed, considering this as replicable is questionable. Moreover,
although it is sometimes argued that replicability is less of an issue when
the sample size is large, as we show below, it remains a problem when the
effect size is moderate.

A second major issue in the context of hypothesis testing is that when
sampling or data-splitting is used in the construction of a test, we may
expect a loss of power due to the sample not being fully utilised. In the
context of debiased machine learning, cross-fitting can sometimes be applied
to alleviate the two drawbacks, to a certain extent. Under conditions, the
per-fold estimates are asymptotically independent and jointly Gaussian, thus
giving an approximately Gaussian final averaged statistic. However, in finite
samples there can be visible correlations between these estimates resulting in
under-coverage of standard confidence intervals; we discuss this further in
Section 4.3.

To illustrate the two main drawbacks mentioned above, consider the following
toy example.

######  Example 1 (Gaussian location experiment).

Let $T^{(1)},T^{(2)},\dots$ be test statistics resulting from repeatedly
applying a randomised (e.g., data splitting) procedure to a given dataset.
Because $T^{(1)},T^{(2)},\dots$ are iid conditional on the data,
_unconditionally_ they are exchangeable in the sense that the joint
distribution is invariant under any finite permutation (Kallenberg, 2005,
§1.1). Suppose the statistic is constructed such that the location of
$T^{(1)}$ captures the signal. In particular, suppose
$T^{(1)}\sim\mathcal{N}(\mu,1)$ marginally and we are interested in testing
$H_{0}:\mu=0$ against $H_{1}:\mu>0$. The corresponding $\alpha$-level test is
to reject $H_{0}$ when $T^{(1)}>z_{\alpha}:=\Phi^{-1}(1-\alpha)$, which we
refer to as the single-split test. We compare this to a test that properly
aggregates the values from $L$ realizations (e.g., random splits). In this
case, it is natural to aggregate by taking the average
$\bar{T}_{1:L}:=(T^{(1)}+\dots+T^{(L)})/L$. The aggregated test rejects
$H_{0}$ when $\bar{T}_{1:L}>\bar{t}_{\alpha}$, where the critical value
$\bar{t}_{\alpha}$ is determined by the null distribution of $\bar{T}_{1:L}$.

In order to permit closed form expressions for the distribution of
$\bar{T}_{1:L}$ and hence $\bar{t}_{\alpha}$, let us assume
$T^{(1)},T^{(2)},\dots$ are jointly normal, i.e., follow a Gaussian process.
One can show that the likelihood ratio test against any fixed $\mu>0$ is
monotone in $\bar{T}_{1:L}$. Due to the exchangeability, the distribution is
parametrised by $(\mu,\rho)$, where $\rho\in[0,1)$ is the pairwise correlation
between $T^{(i)}$ and $T^{(j)}$ for every $i\neq j$. In practice, we expect
$\rho>0$.

Let $\phi_{T^{(l)}}:=\mathbb{I}\\{T^{(l)}>z_{\alpha}\\}$ be the test
corresponding to the $l$-th randomised test statistic. For each $\mu$ and
$\rho$, we may compute the probability that the test cannot be replicated by
another application of the same procedure on the same dataset:
$\operatorname{\mathbb{P}}(\phi_{T^{(1)}}\neq\phi_{T^{(2)}})$. Similarly, the
probability that the aggregated test does not replicate is
$\operatorname{\mathbb{P}}(\phi_{\bar{T}_{1:L}}\neq\phi_{\bar{T}_{(L+1):2L}})$,
where $\phi_{\bar{T}_{1:L}}:=\mathbb{I}\\{\bar{T}_{1:L}>\bar{t}_{\alpha}\\}$.
The two probabilities (see Appendix A for their expressions) are compared in
Fig. 1.1. Note that the probability of non-replication for $\phi_{T^{(1)}}$
can be quite high when the effect size (relative to the sample size for
constructing $\phi_{T^{(1)}}$) is neither too weak nor too strong. The
probability is significantly reduced by aggregating $L=200$ realizations and
approaches zero as $L$ tends to infinity.

Moreover, the aggregation also boosts the power as evidenced by Fig. 1.1. One
can show that the power of the aggregated test is given by

$\operatorname{\mathbb{E}}\phi_{\bar{T}_{1:L}}=\Phi\left(\frac{\mu}{\sqrt{1/L+\rho(L-1)/L}}-z_{\alpha}\right),$

where $L=1$ corresponds to the power of the single-split test. For a small
$\mu$ and a large $L$, we have

$\operatorname{\mathbb{E}}\phi_{\bar{T}^{(1)}}\approx\alpha+\phi(z_{\alpha})\mu,\quad\operatorname{\mathbb{E}}\phi_{\bar{T}_{1:L}}\approx\alpha+\phi(z_{\alpha})\mu/\sqrt{\rho},$

so the slope of local power is improved by a factor of $1/\sqrt{\rho}$. This
improvement can be particularly significant when $\rho$ takes a small positive
value, which is common for the settings considered in the paper. For example,
statistics resulting from two random splits are expected to be weakly
correlated.

Figure 1.1: Probability of non-replication (left) and power (right) for the
Gaussian location experiment in Example 1 ($\alpha=0.05$, $L=200$). The
conservative aggregation rule rejects $H_{0}$ when $\bar{T}_{1:L}>2z_{\alpha}$
(see Appendix B).

From the example above, one may be tempted to conclude that to perform the
aggregated test, we only need to estimate the correlation $\rho$ (e.g.,
through the empirical variance computed from $T^{(1)},T^{(2)},\dots$).
However, the joint normality assumption on the statistics
$T^{(1)},T^{(2)},\dots$, which further implies the normality of
$\bar{T}_{1:L}$, need not hold in practice. In other words, even if every
$T^{(l)}$ is marginally normal or asymptotically normal, the _dependence_
among the statistics need not be a normal copula even in large samples.
Consider the following example due to Kim and Ramdas (2020).

######  Example 2.

Let $X_{1},\dots,X_{n}$ be iid random vectors in $\mathbb{R}^{p}$ with mean
$\mu$ and covariance $\Sigma$. We are interested in testing $H_{0}:\mu=\bm{0}$
against $H_{1}:\mu\neq\bm{0}$. Let $I_{n}$ be a random subset of
$\\{1,\dots,n\\}$ of size $n_{1}$ and let $n_{2}:=n-n_{1}$. Let
$\widehat{\mu}_{1}:=n_{1}^{-1}\sum_{i\in I_{n}}X_{i}$ and
$\widehat{\mu}_{2}:=n_{2}^{-1}\sum_{i\in I_{n}^{c}}X_{i}$. We expect
$\widehat{\mu}_{1}^{\top}\widehat{\mu}_{2}>0$ if $\mu$ deviates from zero.
Treating $\widehat{\mu}_{1}$ as a fixed vector we see that
$\widehat{\mu}_{1}^{\top}\widehat{\mu}_{2}=n_{2}^{-1}\sum_{i\in
I_{n}^{c}}X_{i}^{\top}\widehat{\mu}_{1}$ is simply a sample mean computed from
$I_{n}^{c}$. By studentising this quantity, we may obtain via a central limit
theorem (under appropriate conditions), a test statistic

$T_{n}(X_{1},\dots,X_{n};I_{n}):=\frac{\sqrt{n_{2}}\widehat{\mu}_{1}^{\top}\widehat{\mu}_{2}}{\widehat{\mu}_{1}^{\top}\widehat{\Sigma}_{2}\widehat{\mu}_{1}}\rightarrow_{d}\mathcal{N}(\|\mu\|^{2},1)$

as $n\rightarrow\infty$ and $n_{2}/n\rightarrow q\in(0,1)$, where
$\widehat{\Sigma}_{2}$ is the empirical covariance computed from $I^{c}$. We
reject $H_{0}$ when $T_{n}$ is large compared to $\mathcal{N}(0,1)$. An
advantage of using such a test statistic is that unlike for example a norm of
the empirical mean from the whole sample, the limit distribution does not
depend delicately on the asymptotic limit of $p/n$.

As indicated earlier, a disadvantage is that the approach does not fully
utilise the information in the sample. One might hope that this can be
alleviated by considering the cross-fitted statistic
$[T_{n}(X_{1},\dots,X_{n};I_{n})+T_{n}(X_{1},\dots,X_{n};I_{n}^{c})]/2$.
However, Kim and Ramdas (2020, Proposition A.1) showed that this does not have
a normal limit, making calibration challenging . Further, let $I_{n}^{(l)}$
for $l=1,\dots,L$ independent random subsets of size $\lfloor qn\rfloor$ and
define $T_{n}^{(l)}:=T_{n}(X_{1},\dots,X_{n};I_{n}^{(l)})$ and the aggregated
statistic $S_{n}:=\sum_{l=1}^{L}T_{n}^{(l)}/L$. We see in Fig. 1.2 that its
sampling distribution under
$X_{1},\dots,X_{n}\overset{\textsf{iid}}{\sim}\mathcal{N}(\mathbf{0},\Sigma)$
is clearly non-normal even in large sample.

Figure 1.2: Non-normal null distribution of the aggregated statistic $S_{n}$
in Example 2 under
$X_{1},\dots,X_{n}\overset{\textsf{iid}}{\sim}\mathcal{N}(\mathbf{0},\Sigma)$
with $L=200$, $p=3$ and $\Sigma_{ij}=2^{-|i-j|}$.

Using an aggregate $S_{n}$ of the exchangeable test statistics
$T_{n}^{(1)},\dots,T_{n}^{(L)}$ is more sensitive to departures from the null
as it makes better use of the full data. In addition, because the conditional
variability in $S_{n}$ given the data decreases in $L$ (e.g., like $1/L$ when
$S_{n}$ is the average), by taking a relatively large $L$, $S_{n}$ is
effectively derandomised; see Sections 4.1.1 and D.1 for the improvements in
Example 2. However it should be clear from this very simple example that the
main challenge in general lies in calibrating $S_{n}$. In particular, we need
to handle the unknown, potentially complicated dependence among
$T_{n}^{(1)},\dots,T_{n}^{(L)}$, which typically causes $S_{n}$ to not follow
any textbook distribution. To solve this problem, we introduce a data-driven
calibration scheme based on subsampling to obtain tests with sizes that under
mild assumptions will approach their nominal levels. Our approach is
applicable to all of the randomised tests mentioned above, and as we will see,
by inverting particular aggregated hypothesis tests, we may also obtain
confidence intervals in double machine learning settings that can deliver
improved coverage in finite samples. Before introducing our method, we briefly
summarise some existing proposals in the literature.

###  1.1 Existing proposals and related literature

There is a long line of work that has considered the problem of combining
multiple test statistics. Most approaches either aim to control the type-I
error under an arbitrary dependence of the test statistics, or make some
specific assumptions about their dependence, such as a Gaussian copula.

Most approaches of the first type consider aggregating p-values, i.e., when
$T^{(1)}$ is marginally (super-)uniformly distributed under the null.
Classical examples include Bonferroni correction, i.e., taking the minimum
p-value and multiplying by the total number of p-values, and the average of
the p-values multiplied by $2$ (Rüschendorf, 1982; Meng, 1994). Related
results have been shown for quantiles (Meinshausen et al., 2009; DiCiccio et
al., 2020) and generalised means (Vovk and Wang, 2020). Other aggregation
rules include taking a weighted sum of Cauchy transformations (Liu and Xie,
2020), those developed through concentration inequalities (DiCiccio et al.,
2020) and those involving converting p-values to the so-called e-values (Vovk
and Wang, 2021), to name a few.

While these methods have the attractive guarantee of giving finite sample
valid p-values, as they necessarily must cater for the worst-case dependence,
this benefit comes with the downside of the combination methods typically
being conservative. Indeed, as pointed out by DiCiccio et al. (2020), when
used with p-values produced through sample-splitting, these methods can
sometimes be outperformed by using a single-split test; see Section D.3 for
several such numerical examples from applications considered in this paper.
For combining $Z$-statistics, the right panel of Fig. 1.1 illustrates a
similar phenomenon in the context of Example 1 with the ‘conservative’ rule
that rejects when $\bar{T}_{1:L}>2z_{\alpha}$, where the extra factor of $2$
guarantees its validity; see Theorem B.1 in the Appendix. Indeed, as we show
in Proposition B.1, directly comparing $\bar{T}_{1:L}$ to a standard normal
can lead to a size of $2\alpha$, even when enforcing exchangeability of the
underlying $Z$-statistics. This may come as a surprise given that by Jensen’s
inequality, $\bar{T}_{1:L}$ has variance at most $1$ for example; see Appendix
B for further results on a scaling constant $c\geq 1$ that ensures comparing
$c\bar{T}_{1:L}$ to $z_{\alpha}$ yields a valid test. An analogous negative
result for averaging p-values, which shows that the doubling rule cannot be
improved under exchangeability, is proved in Choi and Kim (2022).

In view of this, other work has considered approaches relying on asymptotic
joint normality of the underlying test statistics
$T_{n}^{(1)},T_{n}^{(2)},\dots$; see for example Romano and DiCiccio (2019,
Theorem 4.1), DiCiccio (2018, Theorem 3.2), Tian et al. (2021) and Liu et al.
(2022). However, as evidenced by the case of Example 2, such a dependence
structure among the test statistics is unlikely to hold in practice for data
splitting.

Our work also connects more broadly to a body of literature on subsampling;
see Politis et al. (1999) for a monograph on the topic. Berg et al. (2010);
McMurry et al. (2012) consider subsampling for hypothesis testing and use data
to centre the estimated null distribution to improve power; this is in similar
spirit to our rank transform introduced in Section 2.2 which would enforce the
mean and all other moments of the null. Subsampling has also been used to
reduce the variance of an estimator through what is known as ‘subagging’
(Bühlmann and Yu, 2002). Stability selection procedures exploit this for
variable selection (Meinshausen and Bühlmann, 2010; Shah and Samworth, 2013),
and can provide a form of finite-sample error control for a user-chosen
variable selection method.

A large number of randomised algorithms have been developed for various
testing problems to which our method would be applicable. Some of these have
been mentioned in the introduction under the umbrellas of hunt and test and
reweighting. Other methods that do not explicitly fall within these categories
include approximate co-sufficient sampling (Barber and Janson, 2022) for
goodness-of-fit testing and several approaches for assessing variable
importance nonparametrically (Cai et al., 2022; Williamson et al., 2022+; Dai
et al., 2022+; Tansey et al., 2022).

###  1.2 Main contributions and organisation of the paper

We develop a general framework for hypothesis testing with an aggregated
statistic $S_{n}:=S(T_{n}^{(1)},\dots,T_{n}^{(L)})$, where $S(\cdot)$ is a
user-specified, symmetric aggregation function such as the arithmetic mean. To
fully exploit the signal in $S_{n}$, it is essential that the aggregated test
is not conservative. Therefore, we construct a test from $S_{n}$ whose type-I
error asymptotically approaches the nominal level $\alpha$, and which
outperforms existing conservative aggregation rules.

To achieve this, we employ subsampling (Politis and Romano, 1994; Politis et
al., 1999), a generic tool for approximating sampling distributions. However,
using subsampling alone is not enough, because it not only approximates the
sampling distribution under the null, but also mimics the sampling
distribution under local alternatives. Hence, naively comparing $S_{n}$ to its
subsampling critical value leads to a test powerless against local
alternatives. To fix this issue, we exploit the fact that the asymptotic null
distribution of $T_{n}^{(1)}$ is typically known, e.g., $\mathrm{unif}(0,1)$
for a p-value or $\mathcal{N}(0,1)$ for a $Z$-statistic. We introduce a rank
transform that we apply to the subsampled statistics to enforce the known null
marginal distribution of $(T_{n}^{(1)},\dots,T_{n}^{(L)})$. In other words, to
approximate the null distribution of $(T_{n}^{(1)},\dots,T_{n}^{(L)})$ and
hence the aggregated $S_{n}$, we effectively combine the _known marginal null
distribution_ with the _unknown copula_ estimated from subsampling. We
demonstrate favourable performance of our method with three types of
applications.

  1. (i)

“Hunt-and-test” procedures: we use our framework to develop a new tests for
the goodness of fit of parametric quantile regression models, and for testing
unimodality in high-dimensional data. We illustrate the latter on cancer gene
expression data to detect the presence of subtypes.

  2. (ii)

Testing hypotheses under reweighting or “distributional shift” with resampling
or rejection sampling: specifically we study testing the sharp null of no
direct treatment effect in the context of a sequentially randomised trial.

  3. (iii)

Calibrating confidence intervals for cross-fitted, double/debiased machine
learning estimators: in our simulations we look in particular at confidence
intervals for partially linear models, though the methodology we develop is
applicable more broadly.

The rest of this paper is organised as follows. In Section 2, we introduce our
rank-transformed subsampling procedure. We present an aggregated, multiple-
split test (Algorithm 2) and a variant (Algorithm 3) that allows for several
user-specified aggregation functions and adapts to the best one. In Section 3,
we study the theoretical properties of our method. We show that these
algorithms give tests that asymptotically have size equal to a given
significance level $\alpha$. Further, we show that if the test statistic and
the aggregated statistic converge uniformly under the null, our procedures
inherit such uniformity in terms of error control. In terms of power, we show
that, if the copula among the single-split statistics converge locally
uniformly, then under any sequence of local alternatives, our test is as
powerful as an oracle test that has access to the null distribution of the
aggregated statistic. In Section 4, we study applications of our method and
present additional theory permitting the use of estimated weights in
randomised tests based on rejection sampling, and further methodology and
theory specific to the construction of confidence intervals for
double/debiased machine learning estimators. Finally, we conclude with a
discussion in Section 5 outlining possible directions for future research. The
appendix contains all proofs, additional simulation results and further
details on the numerical studies. An R package MultiSplit implementing our
method and scripts for reproducing numerical results are available from
https://github.com/richardkwo/MultiSplit.

##  2 Method

###  2.1 Setup

Let $X_{1},\dots,X_{n}\in\mathcal{X}$ be data points drawn iid from an
underlying distribution $P$. We are interested in testing

$H_{0}:P\in\mathcal{P}_{0}\quad\text{vs}\quad
H_{1}:P\in\mathcal{P}\setminus\mathcal{P}_{0},$

where $\mathcal{P}$ is the set of relevant laws of $X$ that includes both the
null and the alternative. Let $T_{n}^{(1)},\dots,T_{n}^{(L)}$ be test
statistics that can be computed from sample $X:=(X_{1},\dots,X_{n})$ and a
piece of external randomness $\Omega$ generated by the analyst. Throughout, we
assume that the random vector

$(T_{n}^{(1)},\dots,T_{n}^{(L)})\quad\text{is exchangeable}$ (1)

under $(X,\Omega)\sim P^{n}\times P_{\Omega}$ for every $P\in\mathcal{P}$,
where $P_{\Omega}$ is the distribution of $\Omega$.

Often, such $T_{n}^{(1)},\dots,T_{n}^{(L)}$ are obtained by applying the same
randomised procedure $L$ times on $X$. That is, we have

$T_{n}^{(l)}=T_{n}(X_{1},\dots,X_{n};\Omega^{(l)}),\quad l=1,\dots,L,$ (2)

where, without loss of generality, we assume
$\Omega^{(l)}\overset{\textsf{iid}}{\sim}\mathrm{unif}(0,1)$ independently
from $X$. Formally, for $n=1,2,\dots$, the randomised test is a measurable map
$T_{n}:\mathcal{X}^{n}\times[0,1]\rightarrow\mathbb{R}$. For example,
$\Omega^{(l)}$ can be used by $T_{n}$ to realise a random data split or a
sequence of $\mathrm{unif}(0,1)$ random variables for acceptance-rejection
sampling. In this case, $\Omega=(\Omega^{(1)},\dots,\Omega^{(L)})$.

Alternatively, it can be that every $T_{n}^{(l)}$ is a deterministic function
of $X_{1},\dots,X_{n}$. This can happen, for example, when there are $L$ pre-
specified ways of splitting the full sample. We will study this case in the
context of cross-fitting in Section 4.3.

Throughout, we require that under the null, $T_{n}^{(1)}$ converges to a
known, continuous distribution $F_{0}$, such as $\mathrm{unif}(0,1)$ or
$\mathcal{N}(0,1)$. Without loss of generality, we assume $H_{0}$ is rejected
for large values of $T_{n}^{(1)}$; other cases can be handled by redefining
$T_{n}^{(1)}$, e.g., replacing $T_{n}^{(1)}$ with $|T_{n}^{(1)}|$ for a two-
sided test, or with $1-T_{n}^{(1)}$ for a p-value. To abuse the term slightly,
we call the test that rejects when $T_{n}^{(1)}>F_{0}^{-1}(1-\alpha)$ the
“single-split” test, even though $T_{n}^{(1)}$ may not be constructed with
data splitting. Consider the aggregated or “multiple-split” statistic

$S_{n}:=S(T_{n}^{(1)},\dots,T_{n}^{(L)})$

constructed with a symmetric, continuous aggregation function
$S:\mathbb{R}^{L}\rightarrow\mathbb{R}$, such as the arithmetic mean or the
maximum. By taking $L$ reasonably large, we can expect that
$\operatorname{\mathrm{var}}[S_{n}\mid X_{1},\dots,X_{n}]$ is small enough
such that the aggregated test statistic is effectively derandomised. Note the
restriction that $S$ is symmetric is rather reasonable: it follows from the
Neyman–Pearson lemma that a most powerful test (for a simple null against a
simple alternative) necessarily combines the exchangeable test statistics
symmetrically in some fashion; see Proposition C.1 in the Appendix.

We will make the mild assumption (see Assumption 1 and following discussion)
that $S_{n}$ converges to some distribution $G_{P}$ under the null; in
practice the limit $G_{P}$ is typically an unknown, non-Gaussian continuous
distribution that depends on $P\in\mathcal{P}_{0}$. Our aggregated test
rejects $H_{0}$ for large values of $S_{n}$, and under the null aims to mimic
an oracle procedure that rejects whenever $S_{n}$ exceeds the unknown upper
$\alpha$ quantile of $G_{P}$. To do this, we use subsampling to compute
$\widetilde{G}_{n}$, an approximation to $G_{P}$, and use its quantile to
determine the critical values for $S_{n}$. As mentioned earlier in Section
1.2, however, $\widetilde{G}_{n}$ must be sufficiently different from the
sampling distribution of $S_{n}$ under local alternatives, so that the test
has non-trivial power. For example, when $T_{n}^{(1)}\sim\mathcal{N}(0,1)$
under $H_{0}$ and $S_{n}:=(T_{n}^{(1)}+\dots+T_{n}^{(L)})/L$, under local
alternatives, $\widetilde{G}_{n}$, as expected from an oracle procedure,
should maintain zero mean even when $S_{n}$ takes a positive mean. This is
achieved by the rank transform introduced below.

###  2.2 Rank-transformed subsampling

In this section, we describe our procedure when using a single aggregation
function $S$. We first introduce some notation relating to distribution
functions and then set out our subsampling scheme.

##### Notation

Given a set of points $\\{x_{i}\\}$ on the real line, we use
$\mathbb{F}_{\\{x_{i}\\}}$ to denote their empirical distribution function.
For a real-valued function $F$, $\|F\|_{\infty}:=\sup_{x}|F(x)|$. For a
distribution function $F$, its upper $\alpha$ quantile is defined as
$F^{-1}(1-\alpha):=\inf\\{x:F(x)\geq 1-\alpha\\}$.

##### Subsampling

Let $m<n$ be a user-chosen subsample size. Our theory will require that
$m\rightarrow\infty$, $m/n\rightarrow 0$; for the description of the
algorithms and all the numerical experiments in this paper, we use $m=\lfloor
n/\log n\rfloor$. We randomly select a total of $B$ sets of indices, each of
size $m$, such that there is a sufficiently low degree of overlap among the
sets. To do this, we first choose a positive integer (e.g., $J=100$) and let
$B:=J\lfloor n/m\rfloor$. Then our collection of sets of indices
$\mathcal{B}:=\\{(i_{1,b},\dots,i_{m,b}):b=1,\dots,B\\}$ is formed using
Algorithm 1 below.

Algorithm 1 Generate ordered tuples

Input: Sample size $n$, subsample size $m$, positive integer $J$.

1:$\mathcal{B}\leftarrow\emptyset$.

2:for $j=1,\dots,J$ do

3: $\pi\leftarrow$ a random permutation of $\\{1,\dots,n\\}$.

4:
$\mathcal{B}\leftarrow\mathcal{B}\cup\left\\{(\pi_{1},\dots,\pi_{m}),\,(\pi_{m+1},\dots,\pi_{2m}),\,\dots,\,(\pi_{(\lfloor
n/m\rfloor-1)m+1},\dots,\pi_{\lfloor n/m\rfloor m})\right\\}$.

5:end for

6:return $\mathcal{B}$

Note that the construction guarantees that $\mathcal{B}$ contains $J$
collections of $\lfloor n/m\rfloor$ sets that are non-overlapping, and so
statistics computed on these subsamples will be independent. This will allow
us to obtain guarantees for subsampling that do not rely on approximating a
scheme (e.g., Politis et al., 1999, §2.4) where statistics on every possible
subsample of size $m$ are evaluated.

Let $\widehat{\mathbf{H}}=(\widehat{H}_{b,l})$ be a $B\times L$ matrix
consisting of rows

$\widehat{\mathbf{H}}_{b,\cdot}:=\left(T_{m}^{(1)}(X_{i_{1,b}},\dots,X_{i_{m,b}}),\quad\dots\quad,T_{m}^{(L)}(X_{i_{1,b}},\dots,X_{i_{m,b}})\right),\quad
b=1,\dots,B,$ (3)

i.e., by computing $(T_{m}^{(1)},\dots,T_{m}^{(L)})$ on each subsample listed
in $\mathcal{B}$. When the statistic is a randomised test in the form of Eq.
2, the external random number is regenerated for every entry of
$\widehat{\mathbf{H}}$; that is,

$\widehat{H}_{b,l}=T_{m}(X_{i_{1,b}},\dots,X_{i_{m,b}};\Omega^{(b,l)}),\quad\Omega^{(b,l)}\overset{\textsf{iid}}{\sim}\mathrm{unif}(0,1),\quad
b=1,\dots,B,\>l=1,\dots,L.$

Note that although we have arranged the subsampled test statistics into a
matrix, entries in the same column but different rows do not correspond
directly to one another, that is $\widehat{H}_{b,l}$ is no more related to
$\widehat{H}_{b^{\prime},l}$ for $b^{\prime}\neq b$ than
$\widehat{H}_{b^{\prime},l^{\prime}}$, where $l^{\prime}\neq l$.

Now if we were to apply the aggregation function $S$ to each row of
$\mathbf{H}$, we would obtain

$\widehat{S}_{b}:=S(\widehat{H}_{b,1},\dots,\widehat{H}_{b,L}),\quad
b=1,\dots,B,$

whose empirical distribution function
$\widehat{G}_{n}(x):=\mathbb{F}_{\\{\widehat{S}_{b}\\}}(x)$ is the natural
subsampling estimate for $G_{P}(x)$. By the standard consistency result of
subsampling (Politis et al., 1999; see also Section C.1), we have
$\|\widehat{G}_{n}-G_{P}\|_{\infty}\rightarrow_{p}0$ under
$P\in\mathcal{P}_{0}$. However, directly using $\widehat{G}_{n}$ to construct
the test is suboptimal because $\widehat{G}_{n}$ would continue to approximate
the asymptotic distribution of $S_{n}$ (were it to exist) under a sequence of
local alternatives, resulting in low power; see the bottom left plot in Fig.
2.1. Therefore, instead of using $\widehat{G}_{n}$ to calibrate our test
statistic $S_{n}$, we perform the rank transform introduced below.

##### Rank transform

Let $\mathbb{F}_{\widehat{\mathbf{H}}}$ denote the empirical distribution
function based on entries of $\widehat{\mathbf{H}}$. With this, we form a
rank-transformed version of $\widehat{\mathbf{H}}$ denoted
$\widetilde{\mathbf{H}}=(\widetilde{H}_{b,l})$ filled with entries

$\begin{split}\widetilde{H}_{b,l}&:=F_{0}^{-1}\left(\mathbb{F}_{\widehat{\mathbf{H}}}(\widehat{H}_{b,l})-1/(2B\,L)\right)\\\
&=F_{0}^{-1}\left(\frac{(\text{rank of $\widehat{H}_{b,l}$ among entries in
$\widehat{\mathbf{H}}$})-1/2}{B\,L}\right).\end{split}$ (4)

We then compute the aggregated statistics

$\widetilde{S}_{b}:=S(\widetilde{H}_{b,1},\dots,\widetilde{H}_{b,L}),\quad
b=1,\dots,B$

and their resulting empirical distribution function
$\widetilde{G}_{n}:=\mathbb{F}_{\\{\widetilde{S}_{b}\\}}(x)$, which we then
use to determine the critical value $\widetilde{G}_{n}^{-1}(1-\alpha)$ for
$S_{n}$. The full procedure is given in Algorithm 2.

Figure 2.1: Illustration of the rank transform. Here $L=2$ and the rows of
$\widehat{\mathbf{H}}$ and $\widetilde{\mathbf{H}}$ are plotted as points in
two dimensions in the left and right panels respectively. Arrows indicate for
certain points their image after the rank transform. We consider aggregation
function $S$ as the arithmetic mean which may be visualised as projection onto
the red dashed line. Top (red) histograms: distributions $\widehat{G}_{n}$ and
$\widetilde{G}_{n}$ (curve: null density of $S_{n}$); side (grey) histograms:
marginal distributions $\mathbb{F}_{\widehat{\mathbf{H}}}$ and
$\mathbb{F}_{\widetilde{\mathbf{H}}}$ ($F_{0}=\mathrm{unif}(0,1)$).

Consider first the null case where $P\in\mathcal{P}_{0}$. As
$n\rightarrow\infty$ (and hence $n/m\rightarrow\infty$, $B\rightarrow\infty$),
by consistency of subsampling and exchangeability of
$(T_{m}^{(1)},\dots,T_{m}^{(L)})$, we expect
$\mathbb{F}_{\widehat{\mathbf{H}}}(\cdot)-1/(2BL)\approx F_{0}(\cdot)$ in Eq.
4. Therefore, $\widetilde{\mathbf{H}}\approx\widehat{\mathbf{H}}$ under the
null. Because $\widetilde{G}_{n}$ is computed from $\widetilde{\mathbf{H}}$ in
the same way as $\widehat{G}_{n}$ is computed from $\widehat{\mathbf{H}}$, we
can expect that $\widetilde{G}_{n}\approx\widehat{G}_{n}\approx G_{P}$ with
high probability under the null. This is formalised in Theorem 3.1 and
illustrated in the top panel of Fig. 2.1, where we see that under the null the
rank-transform leaves the points largely unchanged, and both $\widehat{G}_{n}$
and $\widetilde{G}_{n}$ well-approximate the sampling distribution of $S_{n}$.
Note that the subtraction of $1/2$ from the ranks in Eq. 4 is simply a finite
sample correction to prevent infinity being produced when applying
$F_{0}^{-1}$.

Under an alternative $P\in\mathcal{P}\setminus\mathcal{P}_{0}$, we expect that
$\mathbb{F}_{\widehat{\mathbf{H}}}(\cdot)-1/(2BL)\approx F_{P}(\cdot)$, with
$F_{P}$ the distribution function of the test statistic $T_{m}^{(1)}$ on the
subsample. Thus

$\displaystyle
F_{0}^{-1}\left(\mathbb{F}_{\widehat{\mathbf{H}}}(\widehat{H}_{b,l})-1/(2B\,L)\right)$
$\displaystyle\approx F_{0}^{-1}\left(F_{P}(\widehat{H}_{b,l})\right)\approx
F_{0}^{-1}\left(U_{b,l}\right)$

for some $U_{b,l}\sim\mathrm{unif}(0,1)$. Importantly, this suggests that the
marginal distribution of the rank transformed $\widetilde{H}_{b,l}$ should in
fact continue to be close to $F_{0}$, the asymptotic null distribution of
$\widehat{H}_{b,l}$, just as in the null case. We cannot be sure that joint
distribution of $\widetilde{H}_{b,1},\ldots,\widetilde{H}_{b,L}$ is close to
that of $\widehat{H}_{b,1},\ldots,\widehat{H}_{b,L}$ under the null, and
indeed typically the distribution of the latter would depend on where in the
null $\mathcal{P}_{0}$ one was, so this is certainly too much to ask for.
Nevertheless, using $\widetilde{G}_{n}^{-1}(1-\alpha)$ as a critical value
amounts to comparing $S_{n}$ to the aggregation function $S$ applied to a
version of the random vector $(T_{n}^{(1)},\ldots,T_{n}^{(L)})$ transformed
such that its margins coincide with those expected under any null. This
indicates that using such a calibration should give non-trivial power under
local alternatives, a fact which we formalise in Theorem 3.3. The
effectiveness of the rank transform in restoring null-like behaviour under the
alternative is shown in the bottom panel of Fig. 2.1.

Algorithm 2 Aggregated multiple-split test

Input: Data $(X_{1},\dots,X_{n})$, exchangeable single-split test statistics
$(T_{n}^{(1)},\dots,T_{n}^{(L)})$, asymptotic null distribution function
$F_{0}$, aggregation function $S$, significance level $\alpha\in(0,1)$,
positive integer $J$.

1:$m\leftarrow\lfloor n/\log n\rfloor$, $B\leftarrow J\lfloor n/m\rfloor$.

2:Run Algorithm 1 to obtain $\mathcal{B}=\\{(i_{1,b},\dots,i_{m,b}):b=1,\dots
B\\}$.

3:Initialise $B\times L$ matrices $\widehat{\mathbf{H}}$,
$\widetilde{\mathbf{H}}$ and $B$-dimensional vector $\widetilde{S}$.

4:for $b=1,\dots,B$ do

5:
$\widehat{\mathbf{H}}_{b,\cdot}\leftarrow(T_{m}^{(1)}(X_{i_{1,b}},\dots,X_{i_{m,b}}),\dots,T_{m}^{(L)}(X_{i_{1,b}},\dots,X_{i_{m,b}}))$

6:end for

7:for $b=1,\dots,B$ do

8: for $l=1,\dots,L$ do

9: $\widetilde{H}_{b,l}\leftarrow F_{0}^{-1}\left((\\{\text{rank of
$\widehat{H}_{b,l}$ in $\widehat{\mathbf{H}}$}\\}-1/2)/B\,L\right)$

10: end for

11: $\widetilde{S}_{b}\leftarrow
S(\widetilde{H}_{b,1},\dots,\widetilde{H}_{b,L})$

12:end for

13:$\widetilde{G}_{n}\leftarrow\mathbb{F}_{\\{\widetilde{S}_{b}\\}}$

14:Compute $S_{n}\leftarrow S(T_{n}^{(1)},\dots,T_{n}^{(L)})$ from
$X_{1},\dots,X_{n}$.

15:Reject $H_{0}$ if $S_{n}>\widetilde{G}^{-1}_{n}(1-\alpha)$ and report
p-value $1-\widetilde{G}_{n}(S_{n})$.

###  2.3 Adapting to the best aggregation function

We can further improve power by choosing a good aggregation function $S$. The
performance of an aggregation function, however, depends on the joint
behaviour of the single-split tests under the alternative of interest, which
is usually unknown. For example, we expect
$S=(T^{(1)}_{n}+\dots+T^{(L)}_{n})/L$ to work particularly well if most of the
single-split statistics are large under the alternative; in contrast,
$S=\max(T^{(1)}_{n},\dots,T^{(L)}_{n})$ should perform better if only a few of
them are large under the alternative. This motivates us to allow the user to
specify multiple candidate aggregation functions $S^{1},\dots,S^{W}$, which
are expected to accommodate different cases. In Algorithm 3, we present a
variant of our procedure that aims to adapt to the best aggregation function
among $S^{1},\dots,S^{W}$.

The algorithm rejects for large values of

$R_{n}:=\max\left(\widetilde{G}_{n}^{1}(S_{n}^{1}),\dots,\widetilde{G}_{n}^{W}(S_{n}^{W})\right),$
(5)

where $\widetilde{G}_{n}^{w}$ and $S_{n}^{w}$ respectively are the
counterparts of $\widetilde{G}_{n}$ and $S_{n}$ in Algorithm 2 but relate to
the $w$-th aggregation function. The quantity $R_{n}$ in Eq. 5 is therefore
one minus the minimum p-value corresponding to each of the aggregation
functions. Thus if any one of the aggregation functions yields good power, we
should expect $R_{n}$ to be large: in this way, the test statistic aims to
achieve power close to that of the best $S^{w}$ under consideration.

We could calibrate $R_{n}$ using a Bonferroni correction, but this would give
a conservative test potentially sacrificing any power we may have gained in
using multiple aggregation functions. Instead, we can _reuse_ our subsampling
aggregate statistics $\widetilde{S}_{b}^{w}$ to approximate the sampling
distribution of $R_{n}$ under the null; the subsampled versions of $R_{n}$
used for this are computed in line 10 of Algorithm 3. The advantage of this
approach is that it properly takes account of the dependence among
$\widetilde{G}_{n}^{1}(S_{n}^{1}),\dots,\widetilde{G}_{n}^{W}(S_{n}^{W})$
involved in the construction of $R_{n}$. As a consequence, the resulting test
has asymptotic size equal to its prescribed level (Theorem 3.2), and so in
this sense no power is lost.

Algorithm 3 Aggregated multiple-split test that adapts to the best aggregation
function

Input: As in Algorithm 2 but with the single aggregation function $S$ replaced
by a collection $S^{1},\dots,S^{W}$.

1:Run up to line 10 of Algorithm 2 to obtain $\widetilde{H}_{b,l}$ for
$b=1,\ldots,B$ and $l=1,\ldots,L$.

2:Initialise $B$-dimensional vectors $\widetilde{R}$,
$\widetilde{S}^{1},\dots,\widetilde{S}^{W}$.

3:for $w=1,\dots,W$ do

4: for $b=1,\dots,B$ do

5: $\widetilde{S}_{b}^{w}\leftarrow
S^{w}(\widetilde{H}_{b,1},\dots,\widetilde{H}_{b,L})$

6: end for

7:
$\widetilde{G}_{n}^{w}\leftarrow\mathbb{F}_{\\{\widetilde{S}^{w}_{b}:\,b=1,\dots,B\\}}$

8:end for

9:for $b=1,\dots,B$ do

10:
$\widetilde{R}_{b}\leftarrow\max\left(\widetilde{G}_{n}^{1}(\widetilde{S}^{1}_{b}),\dots,\widetilde{G}_{n}^{W}(\widetilde{S}^{W}_{b})\right)$

11:end for

12:$\widetilde{Q}_{n}\leftarrow\mathbb{F}_{\\{\widetilde{R}_{b}\\}}$

13:for $w=1,\dots,W$ do

14: Compute $S_{n}^{w}\leftarrow S^{w}(T_{n}^{(1)},\dots,T_{n}^{(L)})$ from
$X_{1},\dots,X_{n}$.

15:end for

16:$R_{n}\leftarrow\max\left(\widetilde{G}_{n}^{1}(S_{n}^{1}),\dots,\widetilde{G}_{n}^{W}(S_{n}^{W})\right)$

17:Reject $H_{0}$ if $R_{n}>\widetilde{Q}_{n}^{-1}(1-\alpha)$ and report
p-value $1-\widetilde{Q}_{n}(R_{n})$.

##  3 Theory

###  3.1 Behaviour under the null

In this section, we establish that our algorithms lead to asymptotically level
$\alpha$ tests under a set of mild assumptions, which ensure the consistency
of subsampling and the validity of rank transform.

######  Condition 1 (Asymptotic pivotal null).

For every $P\in\mathcal{P}_{0}$, under $(X,\Omega)\sim P^{n}\times
P_{\Omega}$, as $n\rightarrow\infty$ it holds that

$T_{n}^{(1)}\rightarrow_{d}F_{0},$

where $F_{0}$ is a known, continuous distribution function.

Note that when $T_{n}^{(l)}$ is a randomised test of the form in Eq. 2, a
sufficient condition for the above is that for every fixed $\omega\in(0,1)$
and every $P\in\mathcal{P}_{0}$, it holds that
$T_{n}(X_{1},\dots,X_{n};\omega)\rightarrow_{d}F_{0}$.

We study the tests constructed from rank-transformed subsampling under the
null for two leading cases, when (i) $F_{0}$ is $\mathrm{unif}(0,1)$ (e.g.,
one minus p-value), (ii) $F_{0}$ is standard normal (i.e., $Z$-statistic). For
other null distributions, probability integral transform and its inverse can
be applied to convert the statistic to one of these cases.

First, we show that Algorithm 2 is an asymptotic level $\alpha$ test under
mild assumptions. Further, when $T_{n}^{(1)}$ and $S_{n}$ converge uniformly
to their limiting distributions over the null, under a mild finite density
condition, we show that the test also controls size below $\alpha$
_uniformly_. It can be argued that uniform asymptotic size control, as opposed
to pointwise asymptotic size control, is more relevant to practice because in
contrast to the latter, it ensures that the sample size required to control
the actual type-I error below, say 0.051, does not depend on the underlying
$P\in\mathcal{P}_{0}$; see Lehmann and Romano (2005, §11.1). As uniform size
control involves consideration of the behaviour of random variables when they
are governed by different $P$ (rather than a single $P$ in pointwise
asymptotics), in the below, we will use a subscript
$\operatorname{\mathbb{P}}_{P}(\cdot)$ to denote that $(X,\Omega)\sim
P^{n}\times P_{\Omega}$.

Recall that Algorithm 2 rejects $H_{0}$ whenever $S_{n}$ exceeds the upper
$\alpha$ quantile of the rank-transformed subsampling distribution
$\widetilde{G}_{n}$. Control of the size is therefore intimately linked to the
behaviour of $\widetilde{G}_{n}$, for which we will require the following.
Note that in the below, all densities are with respect to the Lebesgue
measure.

######  Condition 2 (Lipschitz aggregation).

The aggregation function $S$ is $1$-Lipschitz in $\|\cdot\|_{\infty}$.

Examples include $S_{n}=(T_{n}^{(1)}+\dots+T_{n}^{(L)})/L$ and
$S_{n}=\max(T_{n}^{(1)},\dots,T_{n}^{(L)})$. Note that given any Lipschitz
$S$, by scaling, the Lipschitz constant of $1$ in 2 is not a restriction;
however such a scaling would affect the value of $g_{P,\max}$ in Assumption 1
below, which we require to be finite.

######  Assumption 1 (Stability of $S_{n}$).

For every $P\in\mathcal{P}_{0}$, under $(X,\Omega)\sim P^{n}\times
P_{\Omega}$, it holds that

$S_{n}\rightarrow_{d}G_{P},$

where $G_{P}$ is a continuous distribution that can depend on $P$. Further,
$G_{P}$ has density function $g_{P}$ such that
$g_{P,\max}:=\sup_{x}g_{P}(x)<\infty$.

If $(T_{n}^{(1)},\dots,T_{n}^{(L)})$ has a limiting joint distribution under
every $P\in\mathcal{P}_{0}$, then the first part of Assumption 1 holds by
definition of $S_{n}$. Given that $(T_{n}^{(1)},\dots,T_{n}^{(L)})$ is
exchangeable with a limit marginal law (1), we typically expect the joint
distribution to be stable as well (see Fig. 1.2 for an example). In fact, for
a given $P$, the sequence $(S_{n})_{n=1}^{\infty}$ is uniformly tight, and so
by Prohorov’s theorem (van der Vaart, 2000, Thm. 2.4), there always exists a
subsequence that converges in distribution; see Proposition C.3 in Appendix C.
Thus the only way the stability assumption may fail is that the copula
$(F_{n,P}(T_{n}^{(1)}),\dots,F_{n,p}(T_{n}^{(L)}))$ in some sense “oscillates”
as $n\to\infty$, where $F_{n,P}$ is the distribution function of
$T_{n}^{(1)}$.

Below we present results on both pointwise and uniform asymptotic size control
for the test in Algorithm 2.

######  Theorem 3.1 (Validity of Algorithm 2).

Let $(X,\Omega)\sim P^{n}\times P_{\Omega}$ for $P\in\mathcal{P}_{0}$. Suppose
$T_{n}^{(1)},\dots,T_{n}^{(L)}$ are exchangeable and 1 holds with
$F_{0}=\mathrm{unif}(0,1)$ or $F_{0}=\mathcal{N}(0,1)$. Suppose $S$ is chosen
such that 2 holds. Then for all $\alpha\in(0,1)$, the following hold:

  1. (i)

Under Assumption 1, the test in Algorithm 2 is pointwise asymptotically level
$\alpha$.

  2. (ii)

Suppose $T_{n}^{(1)}$ and $S_{n}$ converge to their respective limit
distributions uniformly over the null, i.e., for every $x\in\mathbb{R}$

$\sup_{P\in\mathcal{P}_{0}}|F_{n,P}(x)-F_{0}(x)|\rightarrow
0,\quad\sup_{P\in\mathcal{P}_{0}}|G_{n,P}(x)-G_{P}(x)|\rightarrow 0,$

where $F_{n,P}$ and $G_{n,P}$ respectively denote the distribution function of
$T_{n}^{(1)}$ and $S_{n}$. Also, suppose
$\sup_{P\in\mathcal{P}_{0}}g_{\max,P}<\infty$. Then, the test in Algorithm 2
is uniformly asymptotically level $\alpha$, i.e.,

$\sup_{P\in\mathcal{P}_{0}}\left|\mathbb{P}_{P}\left\\{S_{n}>\widetilde{G}_{n}^{-1}(1-\alpha)\right\\}-\alpha\right|\rightarrow
0.$

Now we establish similar results for the adaptive test in Algorithm 3, under
the following joint stability assumption on the chosen aggregation functions
$(S^{1},\dots,S^{W})$.

######  Assumption 2 (Joint stability of multiple aggregation functions).

For every $P\in\mathcal{P}_{0}$, under $(X,\Omega)\sim P^{n}\times
P_{\Omega}$, it holds that

$(S_{n}^{1},\dots,S_{n}^{W})\rightarrow_{d}(S^{1},\dots,S^{W}),$

where every $S^{w}$ has a continuous distribution function $G_{P}^{w}$ that
can depend on $P$. Further, suppose $G_{P}^{w}$ permits a density $g_{P}^{w}$
such that $g_{P,\max}^{w}:=\sup_{x}g_{P}^{w}(x)<\infty$ for $w=1,\dots,W$.

######  Theorem 3.2 (Validity of Algorithm 3).

Let $(X,\Omega)\sim P^{n}\times P_{\Omega}$ for $P\in\mathcal{P}_{0}$. Suppose
$T_{n}^{(1)},\dots,T_{n}^{(L)}$ are exchangeable and 1 holds with
$F_{0}=\mathrm{unif}(0,1)$ or $F_{0}=\mathcal{N}(0,1)$. Also, suppose 2 holds
for every $S^{w}$ ($w=1,\dots,W)$. Then for every $\alpha\in(0,1)$ the
following hold:

  1. (i)

Under Assumption 2, the test in Algorithm 3 has pointwise asymptotic level
$\alpha$.

  2. (ii)

Suppose $T_{n}^{(1)},S_{n}^{1},\dots,S_{n}^{W}$ and
$\max\\{G_{P}^{1}(S_{n}^{1}),\dots,G_{P}^{W}(S_{n}^{W})\\}$ converge to their
respective limit distributions uniformly over the null, i.e., for every
$x\in\mathbb{R}$ and each $w=1,\ldots,W$,

$\displaystyle\sup_{P\in\mathcal{P}_{0}}|F_{n,P}(x)-F_{0}(x)|\rightarrow
0,\quad\sup_{P\in\mathcal{P}_{0}}|G^{w}_{n,P}(x)-G^{w}_{P}(x)|\rightarrow
0,\quad\sup_{P\in\mathcal{P}_{0}}|Q_{n,P}(x)-Q_{P}(x)|\rightarrow 0,$

where $Q_{n,P}$ is the distribution function of
$\max\\{G_{P}^{1}(S_{n}^{1}),\dots,G_{P}^{W}(S_{n}^{W})\\}$, and $G^{w}_{n,P}$
is the distribution function of $S_{n}^{w}$. Also, suppose
$\max_{w}\sup_{P\in\mathcal{P}_{0}}g^{w}_{\max,P}<\infty$.Then, the test in
Algorithm 3 is uniformly asymptotically level $\alpha$, i.e.,

$\sup_{P\in\mathcal{P}_{0}}\left|\mathbb{P}_{P}\left\\{R_{n}>\widetilde{Q}_{n}^{-1}(1-\alpha)\right\\}-\alpha\right|\to
0.$

###  3.2 Power

Our next result shows that for a fixed $P_{0}\in\mathcal{P}_{0}$, if the
copula of $(T_{n}^{(1)},\dots,T_{n}^{(L)})$ converges in a locally uniform
fashion in the neighbourhood of $P_{0}$, then under alternatives local to
$P_{0}$, the test in Algorithm 2 asymptotically achieves the same power as an
oracle test that has access to the null distribution of $S_{n}$, i.e. its
asymptotic distribution under $P_{0}$.

######  Theorem 3.3 (Behaviour under local uniform convergence of copula).

Let $(T_{n}^{(1)},\dots,T_{n}^{(L)})$ be exchangeable under both the null and
the alternative. Suppose 1 holds with $F_{0}=\mathrm{unif}(0,1)$ or
$F_{0}=\mathcal{N}(0,1)$ and $S$ is chosen to satisfy 2.

Fix $P_{0}\in\mathcal{P}_{0}$ and suppose Assumption 1 holds. Suppose that for
every sequence $P_{n}\in\mathcal{P}$ such that $P_{n}\rightarrow_{d}P_{0}$,
under $(X,\Omega)\sim P_{n}^{n}\times P_{\Omega}$, it holds that

$\left(F_{n,P_{n}}(T_{n}^{(1)}),\dots,F_{n,P_{n}}(T_{n}^{(L)})\right)\rightarrow_{d}(C^{(1)},\dots,C^{(L)})$
(6)

for some $(C^{(1)},\dots,C^{(L)})$ whose distribution does not depend on the
sequence $P_{n}$. Then, for every sequence $P_{n}\in\mathcal{P}$ such that
$P_{n}\rightarrow_{d}P_{0}$, under $(X,\Omega)\sim P_{n}^{n}\times
P_{\Omega}$, we have that the critical value of our test in Algorithm 2
converges to that of the oracle test, i.e., for $\alpha\in(0,1)$,

$\widetilde{G}^{-1}_{n}(1-\alpha)\rightarrow_{p}G_{P_{0}}^{-1}(1-\alpha),$

where $\widetilde{G}_{n}$ is computed using Algorithm 2 and $G_{P_{0}}$ is the
asymptotic distribution of $S_{n}$ under $(X,\Omega)\sim P_{0}^{n}\times
P_{\Omega}$.

As a consequence of the convergence of critical values above, at Lebesgue-
almost every $\alpha$, the difference in power between our test and the oracle
test tends to zero. Note that the assumption of local uniform convergence of
the copula in Eq. 6 is guaranteed by Le Cam’s third lemma when
$T_{n}^{(1)},\dots,T_{n}^{(L)}$ are asymptotically linear and the log-
likelihood ratio
$\log(\mathop{}\\!\mathrm{d}P_{n}^{n}/\mathop{}\\!\mathrm{d}P_{0}^{n})$ is
also asymptotically linear; see van der Vaart (2000, Example 6.7).

##  4 Applications

We illustrate our method with three types of applications. First, we study
data-splitting, “hunt-and-test” procedures: specifically we look at the
approach of Kim and Ramdas (2020) for testing whether a high-dimensional
random vector has mean zero; we develop a new test for unimodality in high
dimensions; and introduce a simple, flexible approach for goodness-of-fit
testing of parametric regression models such as parametric quantile
regression. Next, we consider using the data from a distribution $P$ to test a
property of a different distribution $Q$, where $Q$ is related to $P$ through
reweighting. We study this in the context of causal inference, where $P$ is
the observational distribution and $Q$ is an intervened distribution. Finally,
we study the inference of cross-fitted, double/debiased machine learning
estimators. We show that the cross-fold dependence in these estimators, though
often argued to be asymptotically negligible, can in finite sample leads to
under-coverage of confidence intervals. We present an alternative construction
of confidence intervals using our method that captures such dependence and
restores the desired coverage.

For each application, we present numerical results to illustrate and benchmark
the new methods we develop. Additional numerical results, including the
performance of conservative aggregation rules, can be found in Appendix D.

###  4.1 Hunt and test

In this section, we consider testing a hypothesis that can be expressed as a
conjunction of simpler hypotheses

$H_{0}=\bigcap_{\delta\in\mathcal{D}}H_{0}(\delta),$ (7)

where we already have an off-the-shelf test for each $H_{0}(\delta)$. As
explained in the introduction, such null hypotheses are amenable to a hunt and
test, data splitting approach where one part of the data is used to find an
appropriate $\widehat{\delta}$ and the remaining data is used to test
$H_{0}(\widehat{\delta})$. Here is another perspective due to Moran (1973).
Consider testing

$H_{0}:\theta\in\Theta_{0}\quad\text{vs}\quad H_{1}:\theta\in\Theta_{1},$ (8)

where $\Theta_{1}$ does not contain $\theta_{0}$. The alternative parameter
space $\Theta_{1}$ might be so large or heterogeneous that a reasonable test
for $H_{0}$ only has power against certain alternatives in $\Theta_{1}$.
Again, we can split our data and perform hunt-and-test: use the first part to
estimate $\widehat{\theta}_{1}\in\Theta_{1}$ and then use the second part to
test $\theta\in\Theta_{0}$ versus $\theta=\widehat{\theta}_{1}$.

####  4.1.1 Revisiting Example 2

In Example 2, we considered testing $H_{0}:\mu=\mathbf{0}$ versus
$H_{1}:\mu\neq\mathbf{0}$ given data consisting of iid random vectors, using a
hunt and test approach. Clearly, the problem is an instance of Eq. 8; it can
also be viewed as an instance of Eq. 7, namely
$H_{0}=\cap_{\delta\in\mathbb{R}^{p}}\\{\mu:\mu^{\top}\delta=0\\}$. Fig. 4.1
compares the performance of various tests based on aggregated
$S_{n}=(T_{n}^{(1)}+\dots+T_{n}^{(L)})/L$ with $L=200$ and the single-split
test based on $T_{n}^{(1)}$ alone. In our simulation, we draw
$X_{1},\dots,X_{n}\overset{\textsf{iid}}{\sim}\mathcal{N}(\mu,\Sigma)$ where
$\Sigma\in\mathbb{R}^{3\times 3}$ has entries given by
$\Sigma_{ij}=2^{-|i-j|}$ and $\mu=\tau n^{-1/2}v_{1}$, where $v_{1}$ is the
normalised principal eigenvector of $\Sigma$. When $\tau=0$, we see that our
method (Algorithm 2) controls the type-I error at the nominal level. Further,
as $\tau$ grows, its power clearly dominates the single-split test and is
insensitive to the split ratio $q$. In fact, its performance is very close to
the performance of an oracle test that has access to the null distribution of
$S_{n}$. Further, our aggregated test significantly reduces the chance of non-
replication. For example, when $q=0.5$, for a random dataset, there is less
than 5% chance that two applications of Algorithm 2 will give contradicting
results, while the probability can be as large as 30% for the single-split
test; see Section D.1.

Figure 4.1: Testing $\mu=\mathbf{0}$ in Example 2: power (95% CI) at level
$\alpha=0.05$ (dashed horizontal). Location is $\mu=\tau n^{-1/2}v_{1}$ with
$v_{1}$ the normalised principal eigenvector of $\Sigma$. We compare the
single-split test (single) with tests based on $S_{n}$ being the arithmetic
mean (S=avg): Algorithm 2 is our method; oracle is an oracle test that
compares $S_{n}$ to its sampling distribution under the null; no.rank uses the
subsampling distribution $\widehat{G}_{n}$ without the rank transform to
determine the critical value; conservative compares $S_{n}/2$ to a standard
normal, where division by two ensures validity when the statistics are not
jointly normal as evidenced by Fig. 1.2 (see also Appendix B).

####  4.1.2 Testing multivariate unimodality

Testing non-trivial clustering structure of a high-dimensional dataset is a
long-standing problem. The problem cannot be directly answered by clustering
algorithms, because typically these (e.g., k-means or hierarchical clustering)
always return clusters even when the data comes from a homogenous population
(Huang et al., 2015a). This problem is closely connected to selecting the
number of clusters in that our problem may be framed as testing that the true
number is one. Here, we work with Euclidean data and we take the perspective
that there is only one cluster if the population distribution is unimodal.

A univariate distribution is called unimodal if there is a point $a$ such that
the distribution function is convex on $(-\infty,a)$ and concave on
$(a,+\infty)$ (Khintchine, 1938). While there are different notions of
multivariate unimodality, we take _linear unimodality_ as our definition: we
say a random vector $(X_{1},\dots,X_{p})$ is unimodal if $\sum_{i}a_{i}X_{i}$
is unimodal for every non-zero coefficient vector $a=(a_{1},\dots,a_{p})$.
That is,

$H_{0}:\bigcap_{a\neq\bm{0}}\left\\{\sum_{i=1}^{p}a_{i}X_{i}\text{ is
unimodal}\right\\}.$ (9)

Linear unimodality is implied by several other notions related to multivariate
unimodality such as log-concavity (Dharmadhikari and Joag-Dev, 1988, Lemma
2.1); see also Dharmadhikari and Joag-Dev (1988, Theorem 2.15). The
formulation in Eq. 9 naturally leads to the following “hunt-and-test”
procedure after randomly splitting the data into two parts A and B:

  1. 1.

Identify a direction $\widehat{a}$ using any suitable clustering algorithm on
part A of the data;

  2. 2.

Test univariate unimodality of $\sum_{i}\widehat{a}_{i}X_{i}$ on part B.

We note that the idea of reducing to a univariate test is not new, and is for
example used by Ahmed and Walther (2012) in projecting data onto its principal
curve, and in a likelihood ratio test for log-concavity using random
projection and data-splitting (Dunn et al., 2021).

##### Dip hunting test

To identify a good direction $\widehat{a}$, on part A we run a 2-means
algorithm (initialised with k-means++, Vassilvitskii and Arthur, 2006) and
choose $\widehat{a}$ as the normalised vector connecting the two cluster
centres. Then to test for unimodality, we use a test based on the dip
statistic due to Hartigan and Hartigan (1985), which we describe below; an
alternative would be the test of Silverman (1981). Fig. 4.2 shows a schematic
of our procedure, which may be described as “dip hunting” (by analogy with the
bump hunting procedure of Good and Gaskins (1980)).

Let $\mathbb{F}_{n}$ be the empirical distribution function of $Y$ and let
$\mathcal{U}$ be the set of unimodal univariate distributions. The dip
statistic is defined as

$\rho_{n}:=\inf_{Q\in\mathcal{U}}\|\mathbb{F}_{n}-Q\|_{\infty},$

and may be computed efficiently using the R package diptest (Maechler, 2021).
Hartigan and Hartigan (1985) recommend comparing $\rho_{n}$ to the dip
statistic of a sample drawn from $\mathrm{unif}(0,1)$, which serves as the
least favourable null distribution. However, this approach typically results
in very conservative p-values. An alternative approach to calibration that
avoids this problem is offered by Cheng and Hall (1998) who show that when the
density $f$ of $Y$ is unimodal, under mild regularity conditions, we have
$2n^{3/5}\rho_{n}\rightarrow_{d}c\mathcal{Z}$, where $\mathcal{Z}$ is a
particular function of a standard Wiener process. The constant $c$ depends on
the density $f$ and is given by

$c=\left\\{f(x_{0})^{3}/|f^{\prime\prime}(x_{0})|\right\\}^{1/5},$

where $x_{0}$ is the unique mode of $f$. The only unknown quantity $c$ can be
estimated with $\widehat{c}$, which is a plugin from kernel density estimates
$\widehat{f}$ and $\widehat{f^{\prime\prime}}$ evaluated at
$\widehat{x}_{0}=\operatorname*{arg\,max}\widehat{f}$. We use R package kedd
(Guidoum, 2015) to estimate $\widehat{f}$ and $\widehat{f^{\prime\prime}}$,
for which the respective bandwidths are selected with maximum likelihood cross
validation (function h.mlcv, Habbema et al., 1974; Duin, 1976). From our
experience, for $\widehat{c}$ to behave properly, it is essential to centre
and rescale $Y$ so that its value lies between 0 and 1, which does not affect
$c$. Because the asymptotic distribution of the dip statistic only depends on
$c$, Cheng and Hall (1998) suggest the following approach to obtain the
corresponding p-value, which we adopt here. Given our observed dip statistic,
we compare this to the distribution of $\rho_{n}$ based on samples drawn from
a known distribution whose $c$ equals $\widehat{c}$; three families of such
distributions covering the range of $c$ are provided by Cheng and Hall (1998).

Figure 4.2: Dip hunting test for multivariate unimodality

##### Simulations

We consider the following settings.

  1. 1.

Mixture of unit balls. Let $\mathfrak{B}(p)$ be the unit $p$-dimensional ball
centred at the origin. Consider the following density with bounded support:

$f(x)=1/2\,\mathrm{unif}_{\mathfrak{B}(p)}(x)+1/2\,\mathrm{unif}_{\mathfrak{B}(p)}(x-x_{0}).$

Here, $\mathrm{unif}_{\mathfrak{B}(p)}$ denotes the uniform density on the
unit ball $\mathfrak{B}(p)$ and $\|x_{0}\|$ is the Euclidean distance between
the centres of the two balls. We set $\|x_{0}\|=2\tau/\sqrt{2+p}$ so that the
density in the direction connecting the two ball centres111The one-dimensional
projection of $\mathrm{unif}_{\mathfrak{B}(p)}(x)$ has density
$f(x)=\Gamma(p/2+1)/[\sqrt{\pi}\Gamma(p/2+1/2)](1-x^{2})^{(p-1)/2}$ for
$x\in[-1,1]$, which leads to variance $1/(2+p)$. becomes bimodal roughly when
$\tau\geq 1$.

  2. 2.

Mixture of multivariate t’s. We consider a heavy-tailed setting

$f(x)=1/2\,t_{4}(x;\Sigma)+1/2\,t_{4}(x-x_{0};\Sigma),$

where $t_{4}$ is the density of the $p$-dimensional multivariate
t-distribution with 4 degrees of freedom, mean zero and scale matrix
$\Sigma\in\mathbb{R}^{p\times p}$ with entries $\Sigma_{ij}=2^{-|i-j|}$. We
set $x_{0}=\tau v_{2}\sqrt{p}$, where $v_{2}$ is the second normalised
eigenvector of $\Sigma$.

Fig. 4.3 shows the results based on sample size $n=1000$. We compare the
single-split dip hunting test (single) with the multiple-split versions that
aggregate $L=50$ dip hunting p-values, including $S=\text{avg}$ and $S=\min$
(Algorithm 2), as well as the adaptive test (Algorithm 3) with
$(S^{1}=\text{avg},S^{2}=\min)$. We compare dip hunting to SigClust (Liu et
al., 2008; Huang et al., 2015b, 2022), a widely-used clustering significance
testing method based on a Gaussian mixture model. Other methods include a
nonparametric bootstrap approach suitable for ellipsoidal clusters (Maitra et
al., 2012), and an approach based on simulating from an estimated Gaussian
copula model (Helgeson et al., 2021); see also the review paper Adolfsson et
al. (2019) and references therein. We see that at the null ($\tau=0$),
SigClust incurs a large type-I error for the multivariate t settings, while
all the dip hunting tests maintain the correct level in both settings; note
that this is to be expected as SigClust assumes a Gaussian distribution under
the null. For the unit ball setting, SigClust loses power as $p$ increases. We
can see that $S=\text{avg}$ is more powerful than $S=\min$ for the unit ball
setting and conversely for the multivariate t setting. Nevertheless, our
adaptive test is able to achieve the better performance between the two and
also shows significant power improvement over the single-split test.

Figure 4.3: Detecting a mixture of two $p$-dimensional unimodal components
whose centres are separated $\propto\tau$ away: power (95% CI) at level
$\alpha=0.05$ (dashed horizontal). The multivariate density is linearly
unimodal when $\tau=0$. We compare dip hunting with clustering significance
testing method SigClust. Test single is the single-split dip hunting test;
S=avg, S=min and S=adaptive are multiple-split tests that aggregate $L=50$ dip
hunting p-values. The S=adaptive test is Algorithm 3 with
$(S^{1}=\text{avg},S^{2}=\min)$, which is able to adapt to the better
performance between the two. See also Fig. D.3 for the performance of
conservative p-value aggregation rules.

##### Gene expression of cancer subtypes

We apply our test to gene expression data on renal cell carcinoma (RCC), which
mainly consists of three subtypes: clear cell (ccRCC), papillary (PRCC) and
chromophobe (ChRCC). We use the ICGC/TCGA Pan-Cancer dataset (ICGC/TCGA Pan-
Cancer Analysis of Whole Genomes Consortium, 2020, available from
https://bit.ly/3r65AUw), which contains mRNA expression levels (FPKM-UQ
normalised) from 111 kidney samples, including 37 cc, 31 P and 43 Ch RCC
cases. We use the expression levels of 1,000 genes that are most relevant to
RCC by selecting those with the highest $|\mu_{1}-\mu_{0}|/\sigma_{0}$, where
$\mu_{1}$ and $\mu_{0}$ are case and control means, $\sigma_{0}$ is the
control standard deviation. We apply both the single-split and the aggregated
($S$ is the arithmetic mean, $L=6000$) dip hunting tests to every subtype,
every mixture of two subtypes and the whole sample. Fig. 4.4 shows the
distribution of p-values. The aggregated test produces stable p-values, which
indicate clear separation between subtypes and relative homogeneity of each
subtype. In contrast, it is more difficult to disentangle subtypes ccRCC and
PRCC from the single-split dip hunting test.

Figure 4.4: Testing homogeneity of gene expression levels of three subtypes of
renal cell carcinoma: clear cell (cc), papillary (P) and chromophobe (Ch). The
plot shows p-values resulting from repeatedly applying single-split and
aggregated ($S=\text{avg}$, $L=6000$) dip hunting tests to each subtype, every
mixture of two subtypes and the whole sample.

####  4.1.3 Goodness-of-fit testing for parametric regression models

As another application, we use hunt and test to construct flexible goodness-
of-fit tests for parametric regression models. The approach we take here is
closely related to the generalised residual prediction (RP) test (Janková et
al., 2020; see also Shah and Bühlmann (2018)) for assessing the goodness of
fit of (potentially high-dimensional) generalised linear models. Generalised
RP tests also employs sample splitting and could equally well benefit from our
rank transform p-value aggregation scheme. However our construction here is
applicable more broadly to test the specification of models of the form

$h(X)=\beta^{\top}X,$ (10)

where $h(X)$ is the conditional mean or a conditional quantile of outcome
$Y\in\mathbb{R}$ given $X\in\mathbb{R}^{p}$.

The starting point of our approach is the simple observation that
$h(X)=\beta^{\top}X$ is equivalent to having $\gamma=0$ in
$h(X)=\beta^{\top}X+\gamma g(X)$ where $g:\mathbb{R}^{p}\to\mathbb{R}$ is any
nonlinear (measurable) function. For a given $g$, we may test for whether
$\gamma=0$ by regressing $Y$ on $(X,g(X))$ and utilising existing inference
tools for the model at hand to assess the significance of $g(X)$. To obtain
good power under an alternative, we would like to pick an appropriate $g$ to
expose the particular lack of fit present. This suggests a hunt and test
procedure where we randomly divide our data into parts A and B, and use part A
to hunt for a suitable $g$, and part B to assess the significance of our
artificially constructed additional covariate $g(X)$.

To find an appropriate $g$, we take inspiration from gradient boosting
(Friedman, 2001) and proceed as follows. Let
$(X_{i},Y_{i})_{i=1}^{n^{\prime}}$ be iid covariate–response pairs in part A
of our data. Suppose we have an M-estimator $\widehat{\beta}$ for estimating
$\beta$ in Eq. 10 that minimises $\sum_{i}\ell(Y_{i}-\beta^{\top}X_{i})$ for
some loss function $\ell:\mathbb{R}\to\mathbb{R}$. Defining residual
$r_{i}:=Y_{i}-\widehat{\beta}^{\top}X_{i}$, upon introducing a potential new
covariate $g(X)$, the loss can be locally approximated as

$\sum_{i}\ell(Y_{i}-\widehat{\beta}^{\top}X_{i}-\gamma
g(X))\approx\sum_{i}\ell(r_{i})-\gamma\sum_{i}g(X_{i})\,\ell^{\prime}(r_{i}),$

where $\ell^{\prime}$ is the derivative of the loss function. Hence, when
$\gamma>0$ and fixing $\left(\sum_{i}g(X_{i})^{2}\right)^{1/2}$, to locally
decrease the loss by the greatest amount, we should attempt to choose $g$ such
that approximately
$(g(X_{i}))_{i=1}^{n^{\prime}}\propto(\ell^{\prime}(r_{i}))_{i=1}^{n^{\prime}}$.
To achieve this, we regress $(\ell^{\prime}(r_{i}))_{i=1}^{n^{\prime}}$ onto
the covariates using any flexible regression or machine learning method, and
take the fitted regression function to be $g$. Since we expect that under an
alternative the resulting $g(X)$ should have a positive coefficient, we take
as our single-split statistic $T_{n}$ a $Z$-statistic for the significance of
$g(X)$ computed on part B. We expect $T_{n}$ to be large and positive under an
alternative.

We demonstrate the effectiveness of this approach for quantile regression.
Consider a quantile regression model specified as
$q_{\tau}(X)=\beta_{0}+\beta^{\top}X$ for a fixed $\tau\in(0,1)$, where
$q_{\tau}(X)$ is the $\tau$-th conditional quantile of $Y$ given
$p$-dimensional covariates $X$. The construction of goodness-of-fit tests, or
more commonly called lack-of-fit tests in the related literature, have largely
relied on asymptotic properties of certain statistics or processes concerning
the residual; see, e.g., Horowitz and Spokoiny (2002); He and Zhu (2003);
Escanciano and Velasco (2010); Escanciano and Goh (2014). These tests tend to
have difficulty scaling up to more than a handful of covariates (Conde-Amboage
et al., 2015). Recently, Dong et al. (2019) recast the goodness-of-fit problem
as a two-sample test problem and developed a different, highly competitive
method that can handle moderate or large $p$. It is worth mentioning that,
unlike our approach that directly repurposes existing parameter inference for
quantile regression (Koenker, 2005, Chap. 3), these aforementioned methods
rely on asymptotic results that can require substantial development.

For quantile regression, we have $\ell_{\tau}(r)=r(\tau-\mathbb{I}_{r<0})$ and
$\ell^{\prime}_{\tau}(r)=\tau-\mathbb{I}_{r<0}$ for $r\neq 0$. Therefore, we
use part A to train a classifier (e.g., random forest, Breiman, 2001) that
predicts the sign of the residual from $X$; we take
$g:\mathbb{R}^{p}\to\\{-1,1\\}$ to be the resulting prediction function. To
improve numerical stability, we also partial out $X$ from $g(X)$ before adding
$g(X)$ to the covariates. Define the single-split test statistic as
$T_{n}:=\sqrt{n/2}\widehat{\beta}^{\prime}/\widehat{\sigma}$, where
$\widehat{\sigma}$ is estimated from bootstrap. The statistic can be readily
computed using R package quantreg (Koenker, 2022).

Fig. 4.5 shows results from a simulation study under sample size $n=1000$.
Covariate vector $X$ is drawn from $p$-dimensional Gaussian with covariance
$\Sigma_{ij}=2^{-|i-j|}$. We fix $\tau=0.5$ and consider two specifications

$\begin{split}\text{(i)}\quad&Y=1+\beta_{0}^{\top}X+(v/\sqrt{n})\eta(X)+(1+X_{2}+X_{3})\varepsilon,\\\
\text{(ii)}\quad&Y=1+\beta_{0}^{\top}X+\left[1+X_{2}+X_{3}+(v/\sqrt{n})\eta(X)\right]\varepsilon,\end{split}$
(11)

where we use non-linear function $\eta(X)=4\sqrt{X_{1}^{2}+X_{2}^{2}}$ to
introduce misspecification. We set $\beta_{0}^{\top}=(-1,2,0,-1,2,0,\dots)$.
Error $\varepsilon$ is drawn from $\\{\mathsf{Exp}(1),t_{3}\\}$ for
specification (i) and from $\mathsf{Exp}(1)$ for specification (ii). Because
$\mathsf{Exp}(1)$ has a non-zero median, observe that the quantile regression
model $q_{\tau}(X)=\beta_{0}+\beta^{\top}X$ is well-specified if and only if
$v=0$. We choose $g$ as a random forest trained with R package ranger (Wright
and Ziegler, 2017). We run Algorithm 2 with $S$ being the arithmetic mean. Our
approach is competitive with the state-of-the-art method of Dong et al.
(2019); see Section D.2 for results in another setting.

Figure 4.5: Testing goodness-of-fit of a quantile regression model
$q_{0.5}(X)=\beta_{0}+\beta^{\top}X$: power (95% CI) at level $\alpha=0.05$
(dashed horizontal) under $n=1000$. The model is well-specified if and only if
$v=0$. The non-linear function in Eq. 11 is
$\eta(X)=4\sqrt{X_{1}^{2}+X_{2}^{2}}$. See also Figs. D.2 and D.4 in the
Appendix.

###  4.2 Testing generalised conditional independence

In this section, we consider the use of randomised procedures in causal
inference. To draw causal conclusions from data, we are often faced with the
challenge that the quantity of interest is defined with respect to an
“intervened” distribution, which is related to but different from the data-
generating distribution. This is also known as a “distributional shift” in
machine learning (Shimodaira, 2000).

For example, consider a two-stage sequentially randomised trial represented by
graph $\mathcal{G}$ in Fig. 4.6: $A_{1}$ is the first treatment, $L$ is the
first outcome, $A_{2}$ is the second treatment and $Y$ is the final outcome.
The first treatment $A_{1}$ is completely randomised; the second treatment
$A_{2}$ is randomised according to $A_{1}$ and $L$. Variable $U$ represents,
for example, the underlying health status that affects both outcomes and is
unobserved. Graph $\mathcal{G}$, typically with extra base covariates that we
omit here, is often employed to represent observational or follow-up studies
where the treatment is time-varying and is affected by previous outcomes. For
example, in HIV studies, we have $A_{t}=1$ ($t=1,2$) if the individual
receives antiretroviral therapy at time $t$ and $A_{t}=0$ otherwise. Outcomes
$L$ and $Y$ denote CD4 cell counts that measure the effectiveness of the
therapy. See, e.g., Hernán and Robins (2020, Chap. 19) for more background.

Suppose we are interested in the first treatment’s _direct effect_ $\tau$ on
the final outcome (i.e., not through the second treatment), represented by the
dashed edge in $\mathcal{G}$. Because $A_{1}$ also affects $Y$ through
$A_{2}$, we cannot learn $\tau$ by regressing $Y$ on $A_{1}$ alone. Moreover,
because of the latent $U$ that affects both $L$ and $Y$ (such variables are
called “phantoms” by Bates et al., 2022), nor can we learn $\tau$ by
additionally adjusting for $L$ or $A_{2}$ (or both) in the regression.

$A_{1}$$L$$A_{2}$$Y$$U$$\tau$$\mathcal{G}$$A_{1}$$L$$A_{2}$$Y$$U$$\tau$$\mathcal{G}^{\prime}$
Figure 4.6: Graph $\mathcal{G}$ depicts a sequentially randomised trial:
$A_{1}$ is the first treatment, $L$ is the first outcome, $A_{2}$ is the
second treatment, and $Y$ is the final outcome. Latent variable $U$ represents
the underlying health status that affects both outcomes. The dashed edge
represents the direct effect from $A_{1}$ on $Y$. Graph $\mathcal{G}^{\prime}$
represents the population where both $A_{1}$ and $A_{2}$ are completely
randomised. In this case, if $A_{1}$ has no direct effect on $Y$ (i.e., the
dashed edge is absent), then we can observe
$A_{1}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}Y$.

To learn $\tau$, it is useful to imagine another trial, drawn as
$\mathcal{G}^{\prime}$ in Fig. 4.6, where both treatments $A_{1}$ and $A_{2}$
are completely randomised so that $A_{1}$ only affects $Y$ directly. There, we
can easily learn $\tau$ by regressing $Y$ on $A_{1}$ in the
$\mathcal{G}^{\prime}$ distribution. Although we did not carry out the
$\mathcal{G}^{\prime}$ trial, its data distribution can be approximated by
reweighting our data obtained from the $\mathcal{G}$ trial according to the
inverse propensity of $A_{2}$ given $A_{1}$ and $L$. Further, we can even
artificially _simulate_ data from the $\mathcal{G}^{\prime}$ trial by
resampling (e.g., importance resampling or rejection sampling) our data from
the $\mathcal{G}$ trial. In other words, reweighting or resampling provides
access to our distribution of interest. This idea underlies a general approach
known as the g-methods (Robins, 1986; Naimi et al., 2017).

Specifically, let us consider testing the sharp null hypothesis

$H_{0}:\text{$A_{1}$ has no individual direct effect on $Y$},$ (12)

which is represented by the dashed edge from $A_{1}$ to $Y$ being absent from
$\mathcal{G}$. As explained earlier, this amounts to
$A_{1}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}Y$ in the
population $Q$ represented by $\mathcal{G}^{\prime}$, given by

$dQ/dP=q(A_{2})/p(A_{2}\mid A_{1},L),$ (13)

where $q$ is an arbitrary positive distribution on $A_{2}$. Constraints as
such, which prescribe independence or conditional independence in a reweighted
distribution, are called generalised conditional independence or Verma
constraints in the literature; see Richardson et al. (2017) and the references
therein.

Because the independence holds under $Q$ instead of $P$, the usual permutation
test is not applicable; nor is it applicable through simple reweighting as
employed by Berrett et al. (2020) for testing regular conditional
independence. Instead, the standard approach in causal inference is through
inverse probability weighting (IPW). Here, we consider an alternative
approach: we resample our data to represent $Q$ according to Eq. 13 and then
use the resampled data to test
$A_{1}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}Y$, e.g., by
a permutation test or any off-the-shelf test for independence.

More generally, as demonstrated by Thams et al. (2021), such a test-after-
subsampling procedure is applicable to testing any property under
“distributional shift”. Suppose we observe iid sample $X_{1},\dots X_{n}\sim
P$ but we are interested in testing a property of a different target
distribution $Q$. Distribution $Q$ is related to $P$ through a density ratio
$r=\mathop{}\\!\mathrm{d}Q/\mathop{}\\!\mathrm{d}P$, which is either known or
can be estimated from $P$. Suppose we already have a suitable test for the
property based on a test statistic
$T_{n}:=T_{n}(\widetilde{X}_{1},\dots,\widetilde{X}_{n})$ for iid sample
$\widetilde{X}_{1},\dots,\widetilde{X}_{n}\sim Q$. Assume $r$ is bounded from
above by a constant $C$. Then, we can test the property in two steps:

  1. 1.

Obtain a sample from $Q$ through rejection sampling;

  2. 2.

Test the property with the test statistic computed from the accepted sample.

The next result shows that the test statistic computed as such inherits the
desired asymptotic distribution while permitting the use of an estimated
density ratio $\widehat{r}_{n}$.

######  Proposition 4.1 (Test after rejection sampling).

Let $T_{n}(X_{1},\dots,X_{n})$ be a test statistic. Define $T_{0}:=0$. With
$X:=(X_{1},\dots,X_{n})$, suppose $T_{n}(X_{1},\dots,X_{n})\rightarrow_{d}T$
under $X\sim Q^{n}$. Let $P$ be a distribution and $C$ be a positive constant,
such that $Q$ is absolute continuous with respect to $P$ and
$r:=\mathop{}\\!\mathrm{d}Q/\mathop{}\\!\mathrm{d}P$ satisfies

$r(x)\leq C,\quad\text{$P$-almost every $x$}.$

Let $\widehat{r}_{n}$ be a random sequence of approximate density ratios such
that almost surely,

$\int\widehat{r}_{n}\mathop{}\\!\mathrm{d}P=1,\quad
0\leq\widehat{r}_{n}(x)\leq C\text{ for $P$-almost every $x$}.$

Also, suppose

$\operatorname{\mathbb{E}}\int|\widehat{r}_{n}-r|\mathop{}\\!\mathrm{d}P=o(n^{-1}).$
(14)

With $U:=(U_{1},\dots,U_{n})$, consider $(X,U)\sim
P^{n}\times\mathrm{unif}^{n}(0,1)$. For $i=1,\dots,n$, let the $i$-th sample
be accepted if $U_{i}<\widehat{r}_{n}(X_{i})/C$. Let the set of accepted
sample be denoted as
$\widetilde{X}^{n}_{1},\dots,\widetilde{X}^{n}_{\Gamma_{n}}$, where
$0\leq\Gamma_{n}\leq n$ is the number of acceptances. Then it holds that

$T_{\Gamma_{n}}(\widetilde{X}^{n}_{1},\dots,\widetilde{X}^{n}_{\Gamma_{n}})\rightarrow_{d}T.$

While the rate in Eq. 14 seems to demand a large separate sample for
estimating $r$, we expect that such rate can be relaxed under additional
conditions on the density and the test statistic. For example, when the
density ratio is parametrically specified and the test statistic is simple,
see Proposition C.2 and the following remark in the Appendix.

When the density ratio $r$ is unbounded, the rejection sampling can be
replaced by the distinct replacement sampling (DRPL) introduced by Thams et
al. (2021), which selects a resampled dataset
$(\widetilde{X}_{1}=X_{i_{1}},\dots,\widetilde{X}_{k}=X_{i_{k}})$ of given
size $k$ by sampling distinct indices $(i_{1},\dots,i_{k})$ with probability
proportional to $\prod_{l=1}^{k}r(X_{i_{l}})$. It is shown that, under a
moment condition on $r$, by choosing $k=o(\sqrt{n})$, the statistic $T_{n}$
evaluated on the DRPL resampled data converges to its asymptotic distribution
under $Q$ (Thams et al., 2021, Theorem 1). However, DRPL produces an
$o(\sqrt{n})$ sample while rejection sampling produces an $O(n)$ sample with
high probability.

Both rejection sampling and DRPL lead to tests that are randomised and
potentially of low power due to a reduced sample size. However, as we
demonstrate here in the case of testing Eq. 12 in Fig. 4.6, our aggregation
scheme is able to restore highly competitive levels of power. For simplicity,
we choose the post-rejection/DRPL-sampling test statistic $T_{n}$ as the
permutation p-value for $|\operatorname{\mathrm{cov}}_{Q}(A_{1},Y)|$, which is
uniformly distributed when
$A_{1}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}Y$ holds
under $Q$. As a benchmark, we also consider testing Eq. 12 with the standard
IPW version of the statistic Robins (1999, p. 10), which directly uses the
sample under $P$. Let

$Z_{i}:=\frac{Y_{i}(A_{1,i}-\operatorname{\mathbb{E}}A_{1})}{p(A_{2,i}\mid
L_{i},A_{1,i})},$

where $\operatorname{\mathbb{E}}A_{1}$ and $p(A_{2}\mid L,A_{1})$ can be
replaced by their consistent, asymptotically linear estimators. Under $H_{0}$,
we have

$\chi_{n}:=\frac{\sum_{i}Z_{i}}{\sqrt{\sum_{i}Z_{i}^{2}}}\rightarrow_{d}\mathcal{N}(0,1),$
(15)

which leads to a two-sided test; see also Proposition C.4 in the Appendix.

We consider the following data generating mechanism. We have binary treatments
$A_{1},A_{2}$ and real-valued outcomes $L,Y$. The latent $U$ is a
4-dimensional random vector. Let

$A_{1}\sim\mathsf{Ber}(1/2),\quad
A_{2}\sim\mathsf{Ber}(\text{expit}(2A_{1}-L+2)),$

and

$U\sim\mathcal{N}(0,\Sigma_{ij}=2^{-|i-j|}),\quad
L=A_{1}+\beta_{U,L}^{\top}U+\varepsilon_{L},\quad Y=\tau
A_{1}-A_{2}+\beta_{U,Y}^{\top}U+\varepsilon_{Y},$

where $\beta_{U,L}=(1,1,-2,2)^{\top}$, $\beta_{U,Y}=(2,-1,3,-10)^{\top}$ and
$\varepsilon_{L},\varepsilon_{Y}\sim\mathcal{N}(0,1)$ independently. Parameter
$\tau$ controls the effect size: $H_{0}$ in Eq. 12 holds if and only if
$\tau=0$.

Fig. 4.7 shows the result for sample size $n=1000$ under $P$. For rejection
sampling, $q(A_{2})$ in Eq. 13 is chosen to maximise the acceptance rate; for
DRPL, we choose $q(A_{2})=p(A_{2})$ and $k=\lfloor\sqrt{n}\rfloor$. We use
Algorithm 2 to aggregate $L=20$ post-rejection/DRPL-sampling permutation
p-values ($S$=avg). Presumably because rejection sampling generates a larger
sample than DRPL, rejection sampling tests are more powerful than DRPL tests
here. For both rejection sampling and DRPL, ‘avg’ considerably boosts the
power while maintaining the correct level at $\tau=0$. Perhaps surprisingly,
the ‘avg’ rejection sampling test outperforms the IPW test in Eq. 15, even
though they target the same population quantity
$\operatorname{\mathrm{cov}}_{Q}(A_{1},Y)$.

Figure 4.7: Testing the absence of the dashed edge in graph $\mathcal{G}$ of
Fig. 4.6: power (95% CI) based on $n=1000$ sample at level 0.05 (dashed
horizontal). The null hypothesis corresponds to $\tau=0$. Rej: post-rejection-
sampling permutation test for $|\operatorname{\mathrm{cov}}_{Q}(A_{1},Y)|$,
DRPL: post-DRPL-sampling (Thams et al., 2021) permutation test for
$|\operatorname{\mathrm{cov}}_{Q}(A_{1},Y)|$, IPW: inverse probability
weighted test based on $\operatorname{\mathrm{cov}}_{Q}(A_{1},Y)=0$ given by
Eq. 15. Rej and DRPL are randomised tests (‘single’). Their ‘$S$=avg’ version
is the aggregated test by taking the average of 20 p-values, calibrated by
Algorithm 2. Left panel: the true propensity $p(A_{2}\mid A_{1},L)$ is used
for sampling; right panel: the propensity used for sampling is estimated using
logistic regression. See also Fig. D.5 in the Appendix for comparison with
conservatively merged p-values.

##### More challenging setting

Our strategy is particularly useful when it is difficult to detect the
dependence between $A_{1}$ and $Y$ under $Q$ with simple, “mean-like”
statistics such as covariance and a more complex statistic must be employed.
In such cases, it is not always clear how to reweight and calibrate the
statistic under $P$ using IPW. In contrast, post-resampling test works out of
the box whenever a test under $Q$ can be constructed.

For example, consider changing the specification of $Y$ to

$Y=-A_{2}+\beta_{H,Y}^{\top}H+(1-A_{1})\varepsilon_{Y}+A_{1}\xi_{\nu},\quad\nu\in(0,\infty],$
(16)

where $\xi_{\nu}$ is an independent $t_{\nu}$-distributed random variable. The
null hypothesis holds when $\nu=\infty$ and hence
$\xi_{\nu}=_{d}\varepsilon_{Y}$. Because
$\operatorname{\mathrm{cov}}_{Q}(A_{1},Y)=0$ under the alternative, using it
as the statistic leads to trivial power. Instead, we employ the maximum mean
discrepancy (MMD) statistic (Gretton et al., 2012) between $Y\mid A_{1}=0$ and
$Y\mid A_{1}=1$ under $Q$, with the Gaussian kernel and bandwidth chosen by
the median heuristic (Gretton et al., 2012, §8), which is able to detect
distributional discrepancy beyond the mean. Under the null hypothesis, the
asymptotic distribution of the statistic depends on unknown parameters
(Gretton et al., 2012, Theorem 12) that are very difficult to handle under
IPW. However aggregated post-resampling permutation tests work as expected;
see Fig. 4.8.

Figure 4.8: Testing the absence of the dashed edge in graph $\mathcal{G}$ of
Fig. 4.6: power (95% CI) under the more difficult setting Eq. 16, where tests
based on $|\operatorname{\mathrm{cov}}_{Q}(A_{1},Y)|$ cannot detect the
effect. The null hypothesis corresponds to $\tau=0$. Here ‘single’ refers to
post-rejection-sampling (Rej) and post-DRPL (DRPL) permutation test of the MMD
statistic; ‘avg’ refers to the corresponding aggregated test based on the
average of 20 p-values, calibrated by Algorithm 2. Sample size is $n=1000$ and
level is 0.05 (dashed horizontal). See also Fig. D.6 in the appendix.

###  4.3 Calibrating cross-fitted double machine learning

Double/debiased machine learning (DML) methods, also known as doubly robust
targeted estimation, employ flexible, machine learning methods to estimate
low-dimensional target parameters. These methods are widely used for a variety
of problems; see Díaz (2020); Kennedy (2022) for recent reviews. In this
section, we follow the setup of Chernozhukov et al. (2018) and focus on
estimation with Neyman orthogonal scores, but our method can be extended to
other settings such as those that are based on influence functions with a
mixed bias property (Rotnitzky et al., 2021).

Suppose we have an estimating equation $\psi(W;\theta,\eta)$ for an unknown,
real-valued target parameter $\theta_{0}$ such that
$\operatorname{\mathbb{E}}\psi(W;\theta_{0},\eta_{0})=0$. Here $W$ denotes the
data point and $\eta_{0}$ denotes two or more unknown nuisance parameters.
Suppose $\psi(W;\theta,\eta)$ satisfies the Neyman orthogonality condition
(Chernozhukov et al., 2018, Definition 2.1) and takes the form

$\psi(W;\theta,\eta)=\psi^{a}(W;\eta)\,\theta+\psi^{b}(W;\eta).$

The DML approach proposes to estimate $\eta_{0}$ using flexible machine
learning methods, and employs sample-splitting to control bias that may
otherwise be introduced through overfitting. Suppose iid data
$W_{1},\dots,W_{n}$ are split into $L$ equal-sized folds. For $l=1,\dots,L$,
let $\widehat{\eta}^{(-l)}$ be the nuisance parameters learned from all the
data except the $l$-th fold. Then, let $\widehat{\theta}^{(l)}$ be the
solution to

$\mathbb{P}_{n}^{(l)}\psi(W;\theta,\eta^{(-l)})=0,\quad l=1,\dots,L,$

where $\mathbb{P}_{n}^{(l)}$ denotes the empirical measure of data in the
$l$-th fold. Under regularity conditions and a sufficiently fast rate of
nuisance estimation (Chernozhukov et al., 2018, Assumptions 3.1 & 3.2), it can
be shown that the single-split statistic

$T_{n}^{(l)}:=\frac{\sqrt{n/L}(\widehat{\theta}_{n}^{(l)}-\theta_{0})}{\sigma}\rightarrow_{d}\mathcal{N}(0,1),\quad
l=1,\dots,L,$ (17)

where the asymptotic variance is

$\sigma^{2}=\operatorname{\mathbb{E}}\psi^{2}(W;\theta_{0},\eta_{0})/\\{\operatorname{\mathbb{E}}\psi^{a}(W;\eta_{0})\\}^{2}.$
(18)

The so-called cross-fitted DML-1 estimator is simply the average
$\widehat{\theta}^{\text{dml}}:=(\widehat{\theta}^{(1)}+\dots+\widehat{\theta}^{(L)})/L$.
Under the same conditions as above, it holds that

$\sqrt{L}S_{n}=\frac{1}{\sqrt{L}}\sum_{l=1}^{T}T_{n}^{(l)}=\frac{\sqrt{n}(\widehat{\theta}_{n}^{\text{dml}}-\theta)}{\sigma}\rightarrow_{d}\mathcal{N}(0,1),$
(19)

where $S_{n}$ is the average of exchangeable statistics
$T_{n}^{(1)},\dots,T_{n}^{(L)}$. Provided with a consistent estimator of
$\sigma$, standard large-sample DML confidence intervals are constructed from
Eq. 19.

In view of Eq. 17, Eq. 19 effectively ignores the between-fold correlations
among statistics $T_{n}^{(1)},\dots,T_{n}^{(L)}$ because their contribution is
asymptotically negligible under the assumptions. However, in practice, such
correlations can be visible and hence undermine the coverage of normal
confidence intervals constructed from Eq. 19. When the between-fold
correlation is $\rho$, the variance of
$\sqrt{n}(\widehat{\theta}_{n}^{\text{dml}}-\theta)$ is roughly
$\sigma^{2}(1+\rho(L-1))$ instead of $\sigma^{2}$. We propose to use Algorithm
2 to account for such correlation and construct confidence intervals with
better finite-sample coverage. In this context, $\widehat{H}_{b,l}$ in
Algorithm 2 is the subsample counterpart of Eq. 17, which involves unknown
$\sigma$ and $\theta_{0}$. Nevertheless, by monotonicity, the rank-transformed
$\widetilde{\mathbf{H}}$ can be computed from the ranks of the subsampling
counterparts of $\widehat{\theta}_{n}^{(l)}$ directly as

$\widetilde{H}_{b,l}=\Phi^{-1}\left(\mathbb{F}_{\\{\widehat{\theta}_{m,b}^{(l)}\\}}(\widehat{\theta}_{m,b}^{(l)})-1/(2BL)\right),\quad
b=1,\dots,B,\quad l=1,\dots,L.$

Here $\widehat{\theta}_{m,b}^{(l)}$ is the $l$-th single-split estimator
computed from the $b$-th subsample. Further, because
$\widetilde{H}_{b,l}\approx\widehat{H}_{b,l}=\sqrt{m/L}(\widehat{\theta}_{m,b}^{(l)}-\theta_{0})/\sigma$
with high probability, we get the following estimator for $\sigma$ as a by-
product:

$\widehat{\sigma}_{\text{ls}}:=\sqrt{m/L}/\widehat{\beta}_{1}.$

Here $\widehat{\beta}_{1}$ the slope from least squares
$\widetilde{H}_{b,l}\sim\beta_{1}\widehat{\theta}_{m,b}^{(l)}+\beta_{0}$. We
perform this least squares regression only using those points for which
$\mathbb{F}_{\\{\widehat{\theta}_{m,b}^{(l)}\\}}(\widehat{\theta}_{m,b}^{(l)})\in(\varepsilon/2,1-\varepsilon/2)$
for some small $\varepsilon>0$ to avoid high-leverage points influencing the
fit too much; in our numerical experiments we use $\varepsilon=0.1$.

######  Proposition 4.2.

Suppose $\widehat{\theta}_{n}^{(1)},\dots\widehat{\theta}_{n}^{(L)}$ are
exchangeable and Eq. 17 holds. Then,
$\widehat{\sigma}_{\text{ls}}\rightarrow_{p}\sigma$.

Consequently, our rank-transformed subsampling confidence interval for
$\theta_{0}$ is

$\left[\widehat{\theta}_{n}^{\text{dml}}-\sqrt{L/n}\,\widehat{\sigma}_{\text{ls}}\,\widetilde{G}_{n}^{-1}(1-\alpha/2),\>\widehat{\theta}_{n}^{\text{dml}}-\sqrt{L/n}\,\widehat{\sigma}_{\text{ls}}\,\widetilde{G}_{n}^{-1}(\alpha/2)\right].$
(20)

As a simple example, consider inferring $\theta_{0}$ in a partially linear
model

$\displaystyle D=m_{0}(X)+V,\quad$
$\displaystyle\operatorname{\mathbb{E}}[V\mid X]=0,$ $\displaystyle
Y=D\theta_{0}+g_{0}(X)+\xi,\quad$
$\displaystyle\operatorname{\mathbb{E}}[\xi\mid D,X]=0$

using $\psi$ given by Robinson’s score function (Chernozhukov et al., 2018,
§4.1)

$\psi(W;\theta,\eta)=\\{Y-l(X)-\theta(D-m(X))\\}(D-m(X)),\quad
l(X):=\operatorname{\mathbb{E}}[Y\mid
X],\>m(X)=\operatorname{\mathbb{E}}[D\mid X].$

In our numerical example, we choose errors $V$ and $\xi$ to be heteroscedastic
(see Section D.4 for details) and fit nuisance parameters $\eta:=(l,m)$ with
random forests. Table 4.1 compares the coverage of Eq. 20 with the coverage of
standard DML confidence intervals based on Eq. 19, where $\sigma$ is estimated
by an empirical version of Eq. 18. While DML confidence intervals tend to
undercover for smaller samples, where $\rho(L-1)$ is larger, our confidence
intervals seem well calibrated.

Table 4.1: Coverage of nominal $95\%$ confidence intervals (brackets: median width of intervals; $\rho$ is the between-fold correlation of $T_{n}^{(l)}$) |  | $n=500$ |  | $n=1000$ |  | $n=2000$
---|---|---|---|---|---|---
method |  | $L=2$ | $L=5$ |  | $L=2$ | $L=5$ |  | $L=2$ | $L=5$
$\rho(L-1)$ |  | 0.46 | 0.31 |  | 0.36 | 0.18 |  | 0.25 | 0.14
Eq. 20 |  | 0.94 [0.21] | 0.93 [0.19] |  | 0.95 [0.15] | 0.95 [0.13] |  | 0.96 [0.10] | 0.95 [0.09]
DML |  | 0.86 [0.15] | 0.88 [0.15] |  | 0.88 [0.11] | 0.92 [0.11] |  | 0.91 [0.08] | 0.92 [0.08]

##  5 Discussion

Rank-transformed subsampling provides a framework that properly aggregates
results from multiple applications of a randomised statistical test. Backed by
this framework, we are free to design explicitly randomised tests that employ
data splitting, resampling or any other random processing to solve difficult
problems. In particular, data splitting and resampling can often reduce a
complex hypothesis to a simpler one, for which an off-the-shelf test can be
repurposed in a “plug-and-play” fashion. Even though the resulting “single-
split” randomised procedure may have high variability or lose power, as we
have demonstrated, these aspects can be significantly improved by our
aggregation algorithms.

There are several aspects worth further investigation. Firstly, it is of
interest to go beyond Algorithm 3 and study how to learn an optimal
aggregation function. Secondly, with a large $L$ and $B$, constructing
$\widehat{\mathbf{H}}$ can be computationally intensive, although it is
embarrassingly parallelisable across both rows and columns. It is desirable to
develop faster approximations to $\widehat{\mathbf{H}}$. Thirdly, we find
$m=\lfloor n/\log n\rfloor$ works surprisingly well in our numerical studies,
in spite of only guaranteeing around 7 independent subsamples for $n=1000$.
Essentially rank transformed subsampling is able to approximate the null
sampling distribution far better than the theory might suggest, and it would
be interesting to investigate why this might be the case, and whether there is
a better choice of subsample size.

##### Acknowledgements

The authors thank Richard Moulange and Sach Mukherjee for help with the TCGA
data, and Richard Samworth and Qingyuan Zhao for helpful comments. Both
authors were funded by EPSRC grant EP/N031938/1.

## References

  * Adolfsson et al. (2019) Andreas Adolfsson, Margareta Ackerman, and Naomi C Brownstein.  To cluster, or not to cluster: An analysis of clusterability methods.  _Pattern Recognition_ , 88:13–26, 2019.
  * Ahmed and Walther (2012) Murat O Ahmed and Guenther Walther.  Investigating the multimodality of multivariate data with principal curves.  _Computational Statistics & Data Analysis_, 56(12):4462–4469, 2012.
  * Baker (2016) Monya Baker.  1,500 scientists lift the lid on reproducibility.  _Nature_ , 533(7604), 2016.
  * Barber and Janson (2022) Rina Foygel Barber and Lucas Janson.  Testing goodness-of-fit and conditional independence with approximate co-sufficient sampling.  _The Annals of Statistics_ , 50(5):2514–2544, 2022.
  * Bates et al. (2022) Stephen Bates, Edward Kennedy, Robert Tibshirani, Valerie Ventura, and Larry Wasserman.  Causal inference with orthogonalized regression adjustment: Taming the phantom.  _arXiv preprint arXiv:2201.13451_ , 2022.
  * Berg et al. (2010) Arthur Berg, Timothy L McMurry, and Dimitris N Politis.  Subsampling p-values.  _Statistics & Probability Letters_, 80(17-18):1358–1364, 2010.
  * Berrett et al. (2020) Thomas B Berrett, Yi Wang, Rina Foygel Barber, and Richard J Samworth.  The conditional permutation test for independence while controlling for confounders.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82(1):175–197, 2020.
  * Breiman (2001) Leo Breiman.  Random forests.  _Machine Learning_ , 45(1):5–32, 2001.
  * Bühlmann and Yu (2002) Peter Bühlmann and Bin Yu.  Analyzing bagging.  _The Annals of Statistics_ , 30(4):927–961, 2002\.
  * Cai et al. (2022) Zhanrui Cai, Jing Lei, and Kathryn Roeder.  Model-free prediction test with application to genomics data.  _Proceedings of the National Academy of Sciences_ , 119(34):e2205518119, 2022.
  * Cheng and Hall (1998) M-Y Cheng and Peter Hall.  Calibrating the excess mass and dip tests of modality.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 60(3):579–589, 1998.
  * Chernozhukov et al. (2018) Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters.  _The Econometrics Journal_ , 21(1):C1–C68, 01 2018.  ISSN 1368-4221.  doi: 10.1111/ectj.12097.
  * Choi and Kim (2022) Woohyun Choi and Ilmun Kim.  Averaging p-values under exchangeability.  _Statistics & Probability Letters_, page 109748, 2022.  doi: https://doi.org/10.1016/j.spl.2022.109748.
  * Conde-Amboage et al. (2015) Mercedes Conde-Amboage, César Sánchez-Sellero, and Wenceslao González-Manteiga.  A lack-of-fit test for quantile regression models with high-dimensional covariates.  _Computational Statistics & Data Analysis_, 88:128–138, 2015.
  * Cox (1975) David R Cox.  A note on data-splitting for the evaluation of significance levels.  _Biometrika_ , 62(2):441–444, 1975.
  * Dai et al. (2022+) Ben Dai, Xiaotong Shen, and Wei Pan.  Significance tests of feature relevance for a blackbox learner.  _IEEE Transactions on Neural Networks and Learning Systems, to appear_ , 2022+.
  * Dharmadhikari and Joag-Dev (1988) Sudhakar Dharmadhikari and Kumar Joag-Dev.  _Unimodality, Convexity, and Applications_.  Academic Press, San Diego, 1988.
  * Díaz (2020) Iván Díaz.  Machine learning in the estimation of causal effects: targeted minimum loss-based estimation and double/debiased machine learning.  _Biostatistics_ , 21(2):353–358, 2020.
  * DiCiccio (2018) Cyrus J DiCiccio.  _Hypothesis Testing Using Multiple Data Splitting_.  Stanford University, 2018.
  * DiCiccio et al. (2020) Cyrus J DiCiccio, Thomas J DiCiccio, and Joseph P Romano.  Exact tests via multiple data splitting.  _Statistics & Probability Letters_, 166:108865, 2020.
  * Dong et al. (2019) Chen Dong, Guodong Li, and Xingdong Feng.  Lack-of-fit tests for quantile regression models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 81(3):629–648, 2019.
  * Duin (1976) R. P. W Duin.  On the choice of smoothing parameters for Parzen estimators of probability density functions.  _IEEE Transactions on Computers_ , 25(11):1175–1179, 1976.
  * Dunn et al. (2021) Robin Dunn, Larry Wasserman, and Aaditya Ramdas.  Universal inference meets random projections: a scalable test for log-concavity.  _arXiv preprint arXiv:2111.09254_ , 2021.
  * Dvoretzky et al. (1956) Aryeh Dvoretzky, Jack Kiefer, and Jacob Wolfowitz.  Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator.  _The Annals of Mathematical Statistics_ , pages 642–669, 1956.
  * Escanciano and Goh (2014) Juan Carlos Escanciano and Sze-Chuan Goh.  Specification analysis of linear quantile models.  _Journal of Econometrics_ , 178:495–507, 2014.
  * Escanciano and Velasco (2010) Juan Carlos Escanciano and Carlos Velasco.  Specification tests of parametric dynamic conditional quantiles.  _Journal of Econometrics_ , 159(1):209–221, 2010\.
  * Friedman (2001) Jerome H Friedman.  Greedy function approximation: a gradient boosting machine.  _Annals of Statistics_ , pages 1189–1232, 2001.
  * Good and Gaskins (1980) IJ Good and RA Gaskins.  Density estimation and bump-hunting by the penalized likelihood method exemplified by scattering and meteorite data.  _Journal of the American Statistical Association_ , 75(369):42–56, 1980.
  * Gretton et al. (2012) Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _Journal of Machine Learning Research_ , 13(25):723–773, 2012.  URL http://jmlr.org/papers/v13/gretton12a.html.
  * Guidoum (2015) Arsalane Chouaib Guidoum.  _kedd: Kernel estimator and bandwidth selection for density and its derivatives._ , 2015.  URL http://CRAN.R-project.org/package=kedd.  R package version 1.0.3.
  * Habbema et al. (1974) JDF Habbema, J. Hermans, and K van den Broek.  A stepwise discriminant analysis program using density estimation.  In _Compstat 1974: Proceedings in Computational Statistics_ , Vienna, 1974. Physica Verlag.
  * Hartigan and Hartigan (1985) John A Hartigan and Pamela M Hartigan.  The dip test of unimodality.  _The Annals of Statistics_ , pages 70–84, 1985.
  * He and Zhu (2003) Xuming He and Li-Xing Zhu.  A lack-of-fit test for quantile regression.  _Journal of the American Statistical Association_ , 98(464):1013–1022, 2003.
  * Helgeson et al. (2021) Erika S Helgeson, David M Vock, and Eric Bair.  Nonparametric cluster significance testing with reference to a unimodal null distribution.  _Biometrics_ , 77(4):1215–1226, 2021.
  * Hernán and Robins (2020) M. A. Hernán and J. M. Robins.  _Causal Inference: What If_.  Chapman & Hall/CRC, Boca Raton, 2020.
  * Horowitz and Spokoiny (2002) Joel L Horowitz and Vladimir G Spokoiny.  An adaptive, rate-optimal test of linearity for median regression models.  _Journal of the American Statistical Association_ , 97(459):822–835, 2002.
  * Huang et al. (2015a) Hanwen Huang, Yufeng Liu, David Neil Hayes, Andrew Nobel, J. S Marron, and Christian Hennig.  Significance testing in clustering.  _Handbook of cluster analysis_ , pages 315–335, 2015a.
  * Huang et al. (2015b) Hanwen Huang, Yufeng Liu, Ming Yuan, and J. S Marron.  Statistical significance of clustering using soft thresholding.  _Journal of Computational and Graphical Statistics_ , 24(4):975–993, 2015b.
  * Huang et al. (2022) Hanwen Huang, Yufeng Liu, and J. S Marron.  _SigClust: Statistical Significance of Clustering_ , 2022.  URL https://CRAN.R-project.org/package=sigclust.  R package version 1.1.0.1.
  * ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium (2020) ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium.  Pan-cancer analysis of whole genomes.  _Nature_ , 578(7793):82–93, 2020.  doi: 10.1038/s41586-020-1969-6.  URL https://doi.org/10.1038/s41586-020-1969-6.
  * Ioannidis (2005) John PA Ioannidis.  Why most published research findings are false.  _PLoS Medicine_ , 2(8):e124, 2005.
  * Janková et al. (2020) Jana Janková, Rajen D Shah, Peter Bühlmann, and Richard J Samworth.  Goodness-of-fit testing in high dimensional generalized linear models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82(3):773–795, 2020.
  * Kallenberg (2005) Olav Kallenberg.  _Probabilistic Symmetries and Invariance Principles_.  Springer, New York, 2005.
  * Kennedy (2022) Edward H Kennedy.  Semiparametric doubly robust targeted double machine learning: a review.  _arXiv preprint arXiv:2203.06469_ , 2022.
  * Khintchine (1938) A Ya Khintchine.  On unimodal distributions.  _Izvestiya Nauchno-Issledovatel’skogo Instituta Matematiki i Mekhaniki_ , 2(2):1–7, 1938.
  * Kim and Ramdas (2020) Ilmun Kim and Aaditya Ramdas.  Dimension-agnostic inference using cross U-statistics.  _arXiv preprint arXiv:2011.05068_ , 2020.
  * Koenker (2005) Roger Koenker.  _Quantile Regression_.  Econometric Society Monographs. Cambridge University Press, 2005.
  * Koenker (2022) Roger Koenker.  _quantreg: Quantile Regression_ , 2022.  URL https://CRAN.R-project.org/package=quantreg.  R package version 5.94.
  * Lehmann and Romano (2005) E. L. Lehmann and Joseph P. Romano.  _Testing Statistical Hypotheses_.  Springer Texts in Statistics. Springer, New York, third edition, 2005\.  ISBN 0-387-98864-5.
  * Lei et al. (2018) Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman.  Distribution-free predictive inference for regression.  _Journal of the American Statistical Association_ , 113(523):1094–1111, 2018.
  * Liu et al. (2022) Wanjun Liu, Xiufan Yu, and Runze Li.  Multiple-splitting projection test for high-dimensional mean vectors.  _Journal of Machine Learning Research_ , 23(71):1–27, 2022.
  * Liu and Xie (2020) Yaowu Liu and Jun Xie.  Cauchy combination test: a powerful test with analytic p-value calculation under arbitrary dependency structures.  _Journal of the American Statistical Association_ , 115(529):393–402, 2020.
  * Liu et al. (2008) Yufeng Liu, David Neil Hayes, Andrew Nobel, and J. S Marron.  Statistical significance of clustering for high-dimension, low–sample size data.  _Journal of the American Statistical Association_ , 103(483):1281–1293, 2008.  URL https://doi.org/10.1198/016214508000000454.
  * Lundborg et al. (2022) Anton Rask Lundborg, Ilmun Kim, Rajen D Shah, and Richard J Samworth.  The projected covariance measure for assumption-lean variable significance testing.  _arXiv preprint arXiv:2211.02039_ , 2022.
  * Maechler (2021) Martin Maechler.  _diptest: Hartigan’s Dip Test Statistic for Unimodality - Corrected_ , 2021.  URL https://CRAN.R-project.org/package=diptest.  R package version 0.76-0.
  * Maitra et al. (2012) Ranjan Maitra, Volodymyr Melnykov, and Soumendra N Lahiri.  Bootstrapping for significance of compact clusters in multidimensional datasets.  _Journal of the American Statistical Association_ , 107(497):378–392, 2012.
  * Massart (1990) Pascal Massart.  The tight constant in the Dvoretzky–Kiefer–Wolfowitz inequality.  _The Annals of Probability_ , pages 1269–1283, 1990.
  * McMurry et al. (2012) Timothy L McMurry, Dimitris N Politis, and Joseph P Romano.  Subsampling inference with k populations and a non-standard behrens–fisher problem.  _International Statistical Review_ , 80(1):149–175, 2012.
  * Meinshausen and Bühlmann (2010) Nicolai Meinshausen and Peter Bühlmann.  Stability selection.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 72(4):417–473, 2010.
  * Meinshausen et al. (2009) Nicolai Meinshausen, Lukas Meier, and Peter Bühlmann.  P-values for high-dimensional regression.  _Journal of the American Statistical Association_ , 104(488):1671–1681, 2009.
  * Meng (1994) Xiao-Li Meng.  Posterior Predictive $p$-Values.  _The Annals of Statistics_ , 22(3):1142 – 1160, 1994.
  * Moran (1973) Patrick AP Moran.  Dividing a sample into two parts a statistical dilemma.  _Sankhyā: The Indian Journal of Statistics, Series A_ , pages 329–333, 1973.
  * Müller and Stoyan (2002) Alfred Müller and Dietrich Stoyan.  _Comparison Methods for Stochastic Models and Risks_.  Wiley Series in Probability and Statistics. John Wiley & Sons Ltd., Chichester, 2002.
  * Naimi et al. (2017) Ashley I Naimi, Stephen R Cole, and Edward H Kennedy.  An introduction to g methods.  _International Journal of Epidemiology_ , 46(2):756–762, 2017.
  * Newey and Robins (2018) Whitney K Newey and James R Robins.  Cross-fitting and fast remainder rates for semiparametric estimation.  _arXiv preprint arXiv:1801.09138_ , 2018.
  * Open Science Collaboration (2015) Open Science Collaboration.  Estimating the reproducibility of psychological science.  _Science_ , 349(6251):aac4716, 2015.
  * Politis and Romano (1994) Dimitris N Politis and Joseph P Romano.  Large sample confidence regions based on subsamples under minimal assumptions.  _The Annals of Statistics_ , pages 2031–2050, 1994.
  * Politis et al. (1999) Dimitris N Politis, Joseph P Romano, and Michael Wolf.  _Subsampling_.  Springer Science & Business Media, 1999.
  * Richardson et al. (2017) Thomas S Richardson, Robin J Evans, James M Robins, and Ilya Shpitser.  Nested markov properties for acyclic directed mixed graphs.  _arXiv preprint arXiv:1701.06686_ , 2017.
  * Robins (1986) James Robins.  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  _Mathematical Modelling_ , 7(9-12):1393–1512, 1986.
  * Robins (1999) James M Robins.  Testing and estimation of direct effects by reparameterizing directed acyclic graphs with structural nested models.  _Computation, Causation, and Discovery_ , pages 349–405, 1999.
  * Romano and DiCiccio (2019) Joseph P Romano and Cyrus DiCiccio.  Multiple data splitting for testing.  Technical report, Department of Statistics, Stanford University, 2019\.  URL https://purl.stanford.edu/fb041jg0790.
  * Rotnitzky et al. (2021) Andrea Rotnitzky, Ezequiel Smucler, and James M Robins.  Characterization of parameters with a mixed bias property.  _Biometrika_ , 108(1):231–238, 2021.
  * Rüschendorf (1982) Ludger Rüschendorf.  Random variables with maximum sums.  _Advances in Applied Probability_ , 14(3):623–632, 1982.
  * Scheidegger et al. (2021) Cyrill Scheidegger, Julia Hörrmann, and Peter Bühlmann.  The Weighted Generalised Covariance Measure.  _arXiv preprint arXiv:2111.04361_ , 2021.
  * Shah and Bühlmann (2018) Rajen D Shah and Peter Bühlmann.  Goodness-of-fit tests for high dimensional linear models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 80(1):113–135, 2018.
  * Shah and Samworth (2013) Rajen D Shah and Richard J Samworth.  Variable selection with error control: another look at stability selection.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 75(1):55–80, 2013.
  * Shimodaira (2000) Hidetoshi Shimodaira.  Improving predictive inference under covariate shift by weighting the log-likelihood function.  _Journal of Statistical Planning and Inference_ , 90(2):227–244, 2000.
  * Silverman (1981) Bernard W Silverman.  Using kernel density estimates to investigate multimodality.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 43(1):97–99, 1981.
  * Solari and Djordjilović (2022) Aldo Solari and Vera Djordjilović.  Multi split conformal prediction.  _Statistics & Probability Letters_, 184:109395, 2022.
  * Tansey et al. (2022) Wesley Tansey, Victor Veitch, Haoran Zhang, Raul Rabadan, and David M Blei.  The holdout randomization test for feature selection in black box models.  _Journal of Computational and Graphical Statistics_ , 31(1):151–162, 2022.
  * Thams et al. (2021) Nikolaj Thams, Sorawit Saengkyongam, Niklas Pfister, and Jonas Peters.  Statistical testing under distributional shifts.  _arXiv preprint arXiv:2105.10821_ , 2021.
  * Tian et al. (2021) Jinjin Tian, Xu Chen, Eugene Katsevich, Jelle Goeman, and Aaditya Ramdas.  Large-scale simultaneous inference under dependence.  _arXiv preprint arXiv:2102.11253_ , 2021.
  * van der Vaart (2000) Aad van der Vaart.  _Asymptotic Statistics_.  Cambridge University Press, 2000.
  * Vassilvitskii and Arthur (2006) Sergei Vassilvitskii and David Arthur.  k-means++: The advantages of careful seeding.  In _Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms_ , pages 1027–1035, 2006.
  * Vershynin (2018) Roman Vershynin.  _High-Dimensional Probability: an Introduction with Applications in Data Science_ , volume 47.  Cambridge University Press, 2018.
  * von Neumann (1951) John von Neumann.  Various techniques used in connection with random digits.  _Applied Math Series_ , 12(36-38):3, 1951.
  * Vovk and Wang (2020) Vladimir Vovk and Ruodu Wang.  Combining p-values via averaging.  _Biometrika_ , 107(4):791–808, 2020.
  * Vovk and Wang (2021) Vladimir Vovk and Ruodu Wang.  E-values: Calibration, combination and applications.  _The Annals of Statistics_ , 49(3):1736–1754, 2021.
  * Wang and Shah (2020) Yuhao Wang and Rajen D Shah.  Debiased inverse propensity score weighting for estimation of average treatment effects with high-dimensional confounders.  _arXiv preprint arXiv:2011.08661_ , 2020.
  * Williamson et al. (2022+) Brian D. Williamson, Peter B. Gilbert, Noah R. Simon, and Marco Carone.  A general framework for inference on algorithm-agnostic variable importance.  _Journal of the American Statistical Association, to appear_ , 2022+.  doi: 10.1080/01621459.2021.2003200.
  * Wright and Ziegler (2017) Marvin N. Wright and Andreas Ziegler.  ranger: A fast implementation of random forests for high dimensional data in C++ and R.  _Journal of Statistical Software_ , 77(1):1–17, 2017.  doi: 10.18637/jss.v077.i01.

##  Appendix A Gaussian location experiment

Here we derive expressions for the power and probability of replication
relating to Example 1. The sequence $T^{(1)},T^{(2)},\dots$ follows a Gaussian
process with mean and covariance

$\mu_{i}=\mu,\quad\Sigma_{\rho}=\begin{cases}1,&\quad i=j\\\ \rho,&\quad i\neq
j\end{cases}.$

It follows that

$\begin{pmatrix}\bar{T}_{1:L}\\\
\bar{T}_{(1+L):2L}\end{pmatrix}:=\begin{pmatrix}L^{-1}\sum_{l=1}^{L}T^{(l)}\\\
L^{-1}\sum_{l=L+1}^{2L}T^{(l)}\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix}\mu\\\
\mu\end{pmatrix},\begin{pmatrix}1/L+\rho(L-1)/L&\rho\\\
\rho&1/L+\rho(L-1)/L\end{pmatrix}\right).$

The critical value is

$\bar{t}_{\alpha}=\sqrt{1/L+\rho(L-1)/L}\,z_{\alpha},$

which reduces to $z_{\alpha}$ when $L=1$. The power function of the aggregated
test given by

$\displaystyle\operatorname{\mathbb{E}}\phi_{\bar{T}_{1:L}}=\operatorname{\mathbb{P}}(\bar{T}_{1:L}>\bar{t}_{\alpha})$
$\displaystyle=\Phi\left(\frac{\mu}{\sqrt{1/L+\rho(L-1)/L}}-z_{\alpha}\right)$
$\displaystyle=\alpha+\phi(z_{\alpha})\frac{\mu}{\sqrt{\frac{1}{L}+\frac{L-1}{L}\rho}}+o(\mu)$

as $\mu\searrow 0$. To derive the expression for
$\operatorname{\mathbb{P}}(\phi_{\bar{T}_{1:L}}\neq\phi_{\bar{T}_{(L+1):2L}})$,
consider the representation

$\begin{pmatrix}\bar{T}_{1:L}\\\
\bar{T}_{(1+L):2L}\end{pmatrix}=_{d}\begin{pmatrix}\mu+\sqrt{\rho}Z_{0}+\sqrt{(1-\rho)/L}Z_{1}\\\
\mu+\sqrt{\rho}Z_{0}+\sqrt{(1-\rho)/L}Z_{2}\end{pmatrix},$

where $Z_{0},Z_{1},Z_{2}$ are independent standard normals. Then we can write

$\displaystyle\quad\operatorname{\mathbb{P}}(\bar{T}_{1:L}>\bar{t}_{\alpha},\bar{T}_{(L+1):2L}>\bar{t}_{\alpha})$
$\displaystyle=\operatorname{\mathbb{E}}\operatorname{\mathbb{P}}(\sqrt{(1-\rho)/L}Z_{1}>\bar{t}_{\alpha}-\mu-\sqrt{\rho}Z_{0},\sqrt{(1-\rho)/L}Z_{2}>\bar{t}_{\alpha}-\mu-\sqrt{\rho}Z_{0}\mid
Z_{0})$
$\displaystyle=\operatorname{\mathbb{E}}\left\\{1-\Phi\left(\frac{\bar{t}_{\alpha}-\mu-\sqrt{\rho}Z_{0}}{\sqrt{(1-\rho)/L}}\right)\right\\}^{2}=\operatorname{\mathbb{E}}\Phi^{2}\left(\frac{\mu+\sqrt{\rho}Z_{0}-\bar{t}_{\alpha}}{\sqrt{(1-\rho)/L}}\right),$

and similarly

$\operatorname{\mathbb{P}}(\bar{T}_{1:L}\leq\bar{t}_{\alpha},\bar{T}_{(L+1):2L}\leq\bar{t}_{\alpha})=\operatorname{\mathbb{E}}\Phi^{2}\left(-\frac{\mu+\sqrt{\rho}Z_{0}-\bar{t}_{\alpha}}{\sqrt{(1-\rho)/L}}\right).$

Thus

$\displaystyle\operatorname{\mathbb{P}}\left(\phi_{\bar{T}_{1:L}}\neq\phi_{\bar{T}_{(L+1):2L}}\right)=1-\int\left\\{\Phi^{2}(y)+\Phi^{2}(-y)\right\\}\phi(x)\mathop{}\\!\mathrm{d}x,$
$\displaystyle
y=\frac{\mu+\sqrt{\rho}z-\sqrt{1/L+\rho(L-1)/L}\,z_{\alpha}}{\sqrt{(1-\rho)/L}},$

where $\phi(x)$ is the standard Gaussian density evaluated and $x$.
$\operatorname{\mathbb{P}}(\phi_{T^{(1)}}\neq\phi_{T^{(2)}})$ is given by the
same expression with $L=1$.

##  Appendix B Tail bound for the average statistic

Let $T^{(1)},\dots,T^{(L)}$ be exchangeable statistics with $T^{(1)}\sim
F_{0}$ under the null. In this Appendix, we consider aggregation rules of the
form

$S:=c\,\bar{T}=c\,(T^{(1)}+\dots+T^{(L)})/L,$ (21)

that guarantee finite sample type-I error control when a test is constructed
by comparing $S$ to the $\alpha$ quantile of $F_{0}$. No assumption is made on
the copula of $(T^{(1)},\dots,T^{(L)})$. We specifically consider two cases:
$F_{0}=\mathrm{unif}(0,1)$ and $F_{0}=\mathcal{N}(0,1)$. We summarise results
from the literature on the former case, and also develop new results on the
latter case.

###  B.1 The $F_{0}=\mathrm{unif}(0,1)$ case

It is well-known that when each $T^{(l)}$ is a valid p-value, Eq. 21 with
$c=2$ is also a valid p-value; see Rüschendorf (1982); Meng (1994); Vovk and
Wang (2021). It has been recently shown that the constant $c=2$ cannot be
improved under an exchangeability condition on $(T^{(1)},\dots,T^{(L)})$ (Choi
and Kim, 2022).

###  B.2 The $F_{0}=\mathcal{N}(0,1)$ case

Without loss of generality, suppose the null hypothesis is rejected whenever
$S>\Phi(1-\alpha)$.

######  Theorem B.1 (Averaging for $Z$-statistics).

Let $T^{(1)},\dots,T^{(L)}$ be exchangeable statistics with
$T^{(1)}\sim\mathcal{N}(0,1)$ under the null. Then, for $0<\alpha\leq
1-\Phi(1)\approx 0.159$, rejecting the null whenever
$S=\bar{T}/2>\Phi(1-\alpha)$ controls the type-I error below $\alpha$.

We will show Theorem B.1 using properties of convex ordering on random
variables.

######  Definition 1 (Convex order).

For random variables $X$ and $Y$, we say $X\leq_{\mathrm{cx}}Y$ if
$\operatorname{\mathbb{E}}f(X)\leq\operatorname{\mathbb{E}}f(Y)$ for every
convex $f$.

For a random variable $X$, let its integrated survival function be defined as

$\pi_{X}(t)=\int_{t}^{\infty}\bar{F}_{X}(s)\mathop{}\\!\mathrm{d}s=\int_{t}^{\infty}(1-F_{X}(s))\mathop{}\\!\mathrm{d}s=\operatorname{\mathbb{E}}(X-t)_{+}.$

It can be shown that $\pi_{X}(t)$ is decreasing and convex.

######  Lemma B.1 (Müller and Stoyan (2002, §1.5)).

$X\leq_{\mathrm{cx}}Y$ if and only if $\pi_{X}(t)\leq\pi_{Y}(t)$ for every $t$
and $\operatorname{\mathbb{E}}X=\operatorname{\mathbb{E}}Y$.

######  Lemma B.2 (Mills’ ratio).

For $\lambda>0$,
$\frac{\lambda}{\lambda^{2}+1}\phi(\lambda)<1-\Phi(\lambda)<\frac{1}{\lambda}\phi(\lambda)$.

The next lemma derives a tail bound for a random variable dominated by a
standard Gaussian in the convex order; see Meng (1994, Lemma 1) for a similar
result on convex order with respect to a $\mathrm{unif}(0,1)$ random variable.

######  Lemma B.3.

Suppose $\bar{Z}\leq_{\mathrm{cx}}Z$ for $Z\sim\mathcal{N}(0,1)$. Then we have

$\mathbb{P}(\bar{Z}>s)=\bar{F}_{\bar{Z}}(s)<\frac{1-\Phi(\beta)}{\beta(s-\beta)},\quad
0<\beta<s.$

  * Proof.

Because $\pi_{\bar{Z}}$ is convex, we have

$\pi_{\bar{Z}}(\beta)\geq\pi(s)+\pi_{\bar{Z}}^{\prime}(s)(\beta-s)=\pi_{\bar{Z}}(s)+\bar{F}_{\bar{Z}}(s)(s-\beta),$

which gives

$\bar{F}_{\bar{Z}}(s)\leq\frac{\pi_{\bar{Z}}(\beta)-\pi_{\bar{Z}}(s)}{s-\beta},\quad\beta<s.$

By $\bar{Z}\leq_{\mathrm{cx}}Z$ and Lemma B.1, we have
$\pi_{\bar{Z}}(\beta)\leq\pi_{Z}(\beta)$ and it follows that

$\bar{F}_{\bar{Z}}(s)\leq\frac{\pi_{Z}(\beta)-\pi_{\bar{Z}}(s)}{s-\beta},\quad\beta<s.$

Now we bound $\pi_{Z}(\beta)$ for $\beta>0$. Using Lemma B.2, we have

$\begin{split}\pi_{Z}(\beta)<\int_{\beta}^{+\infty}t\phi(t)\mathop{}\\!\mathrm{d}t=\int_{\beta}^{+\infty}t\mathop{}\\!\mathrm{d}\Phi(t)&=[t\Phi(t)]|_{\beta}^{+\infty}-\int_{\beta}^{+\infty}\Phi(t)\mathop{}\\!\mathrm{d}(1/t)\\\
&=-\Phi(\beta)/\beta+\int_{\beta}^{+\infty}\Phi(t)t^{-2}\mathop{}\\!\mathrm{d}t\\\
&\leq-\Phi(\beta)/\beta+\int_{\beta}^{+\infty}t^{-2}\mathop{}\\!\mathrm{d}t=\frac{1-\Phi(\beta)}{\beta}.\end{split}$

Hence, for $0<\beta<s$, we have

$\bar{F}_{\bar{Z}}(s)<\frac{(1-\Phi(\beta))/\beta-\pi_{\bar{Z}}(s)}{s-\beta}\leq\frac{1-\Phi(\beta)}{\beta(s-\beta)}.$

∎

  * Proof of Theorem B.1.

We first claim that $\bar{T}\leq_{\mathrm{cx}}T^{(1)}\sim\mathcal{N}(0,1)$.
Indeed, given any convex function $f$, we have

$\operatorname{\mathbb{E}}f(\bar{T})\leq\operatorname{\mathbb{E}}L^{-1}\sum_{l=1}^{L}f(T^{(l)})=\operatorname{\mathbb{E}}f(T^{(1)})$

if $\operatorname{\mathbb{E}}f(T^{(1)})<\infty$, where we used Jensen’s
inequality and exchangeability. Now, take $\beta=s/2$ for $s\geq 2$ in Lemma
B.3. We have

$\mathbb{P}(\bar{T}/2>s/2)=\mathbb{P}(\bar{T}>s)<\frac{1-\Phi(s/2)}{(s/2)^{2}}\leq
1-\Phi(s/2).$

∎

From our experience, the aggregation rule with $c=1/2$ is often quite
conservative; see, for example, Fig. 4.1. It is sometimes argued that using
$c=1$ seems to work in practice (Wang and Shah, 2020; Lundborg et al., 2022).
Proposition B.1 below shows that in certain cases using $c=1$ can lead to
type-I error inflated by a factor of two. Further, Proposition B.2 shows that,
when $\alpha$ is small, the actual type-I error is at most $e\alpha$. However,
determining the optimal $c$, which lies between $1/2$ and $1$, remains an open
question.

######  Proposition B.1.

Under the assumption of Theorem B.1, using $S=\bar{T}$ can lead to size
$2\alpha$ when $L$ is even.

  * Proof.

Let $0\leq\alpha\leq 1/2$ and consider Fig. B.1. Let
$U_{1}\sim\mathrm{unif}(0,1)$ and let

$U_{2}=\begin{cases}U_{1},\quad&U_{1}\leq 1-2\alpha\\\ 2-2\alpha-
U_{1},\quad&U_{1}>1-2\alpha\end{cases}.$

Pair $U_{1},U_{2}$ are exchangeable with uniform margin. Then, let
$T^{(1)}=\Phi^{-1}(U_{1})$ and $T^{(2)}=\Phi^{-1}(U_{2})$. It follows that
$T^{(1)},T^{(2)}$ are exchangeable with standard normal margin. Further, with
probability $2\alpha$, the graph of $(T^{(1)},T^{(2)})$ is the red curve, on
which the minimum value of $\bar{T}$ is $\Phi^{-1}(1-\alpha)$, when
$U_{1}=U_{2}=1-\alpha$. Hence,

$\operatorname{\mathbb{P}}(\bar{T}>\Phi^{-1}(1-\alpha))=2\alpha,$

or equivalently,

$\operatorname{\mathbb{P}}(\bar{T}>t)=2(1-\Phi(t)).$

The construction is specific to $t$ for $t>0$.

It is easy to extend this construction for $L=2k$, $k\geq 1$. To construct
$T^{(1)},\dots,T^{(2k)}$, randomly choose half of them to be equal to
$\Phi^{-1}(U_{1})$ and the rest to be equal to $\Phi^{-1}(U_{2})$. Observe
that $T^{(1)},\dots,T^{(2k)}$ are exchangeable with standard normal margin.
Further, $\bar{T}=(\Phi^{-1}(U_{1})+\Phi^{-1}(U_{2}))/2$ so the bound remains
valid. ∎

######  Proposition B.2 (Tight tail bound up to a constant).

Given $\bar{Z}\leq_{\mathrm{cx}}\mathcal{N}(0,1)$, it holds that

$\operatorname{\mathbb{P}}(\bar{Z}>s)<\frac{1-\Phi(s-1/s)}{1-1/s^{2}},\quad
s>1,$

where $\mathrm{RHS}\sim e(1-\Phi(s))$ as $s\rightarrow\infty$.

  * Proof.

For $s>1$, choose $\beta=s-1/s>0$ in Lemma B.3 and we get

$\operatorname{\mathbb{P}}(\bar{Z}>s)<\frac{1-\Phi(s-1/s)}{1-1/s^{2}}.$

By L’Hôpital’s rule, it is easy to see that

$\lim_{s\rightarrow+\infty}\frac{1-\Phi(s-1/s)}{1-1/s^{2}}/(1-\Phi(s))=\lim_{s\rightarrow+\infty}\frac{\phi(s-1/s)}{\phi(s)}=e.$

∎

######  Remark B.1.

This means for small $\alpha$, comparing $\bar{Z}$ to a standard normal
inflates size by $e$.

Figure B.1: Construction for exchangeable $(T^{(1)},T^{(2)})$ with standard
normal margin such that $\operatorname{\mathbb{P}}(\bar{T}>t)=2(1-\Phi(t))$.

##  Appendix C Proofs of theoretical results

In this Appendix, we prove results for a generic $m=m(n)$ sequence such that
$1<m<n$, $m\nearrow\infty$ and $n/m\rightarrow\infty$, although we choose
$m=\lfloor n/\log n\rfloor$ throughout the paper. We write $B_{n}$ in lieu of
$B$ to highlight its dependence on $n$.

###  C.1 Finite-sample bounds of subsampling

The standard consistency result is given by Politis et al. (1999, Theorem
2.2.1). The following is an extension that allows multiple exchangeable
statistics. Throughout, let
$\widehat{F}_{n}:=\mathbb{F}_{\widehat{\mathbf{H}}}$, where
$\widehat{\mathbf{H}}$ is given by Eq. 3. Let $F_{n,P}$ be the distribution
function of $T_{n}^{(1)}$.

######  Lemma C.1 (Uniform Hoeffding bound).

Let $T_{n}^{(1)},\dots,T_{n}^{(L)}$ be exchangeable statistics under
$(X,\Omega)\sim P^{n}\times P_{\Omega}$. We have
$\operatorname{\mathbb{E}}\widehat{F}_{n}(x)=F_{m,P}(x)$ for every $x$ and

$\operatorname{\mathbb{P}}(\|\widehat{F}_{n}-F_{m,P}\|_{\infty}>t)\leq
2\exp(-2\lfloor n/m\rfloor t^{2}/c_{0}^{2}),\quad t\geq 0,$

where $c_{0}$ is a universal constant.

  * Proof.

By the construction of $\mathcal{B}$ via Algorithm 1 and noting that
$|\mathcal{B}|=B=J\lfloor n/m\rfloor$, we have

$\widehat{F}_{n}(x)-F_{m,P}(x)=\frac{1}{JL}\sum_{l=1}^{L}\sum_{\pi\in\\{\pi_{1},\dots,\pi_{J}\\}}\left[\widehat{F}_{\pi,n}^{(l)}(x)-F_{m,P}(x)\right],$
(22)

where

$\widehat{F}_{\pi,n}^{(l)}(x):=\frac{1}{\lfloor
n/m\rfloor}\biggl{[}\mathbb{I}\left\\{T_{m}^{(l)}(X_{\pi(1)},\dots,X_{\pi(m)})\leq
x\right\\}+\\\
\mathbb{I}\left\\{T_{m}^{(l)}(X_{\pi(m+1)},\dots,X_{\pi(2m)})\leq
x\right\\}+\dots+\mathbb{I}\left\\{T_{m}^{(l)}(X_{\pi((\lfloor
n/m\rfloor-1)m+1)},\dots,X_{\pi(\lfloor n/m\rfloor m)})\leq
x\right\\}\biggr{]}.$

The summands in the square brackets correspond to $\lfloor n/m\rfloor$
disjoint subsamples of data

$(X_{\pi(1)},\dots,X_{\pi(m)}),\quad(X_{\pi(m+1)},\dots,X_{\pi(2m)}),\quad\dots\quad,(X_{\pi((\lfloor
n/m\rfloor-1)m+1)},\dots,X_{\pi(\lfloor n/m\rfloor m)})$

and they are therefore independent.

Fix $\pi$ and $l$. By definition, we know
$\operatorname{\mathbb{E}}\widehat{F}_{\pi,n}^{(l)}(x)=F_{m,P}(x)$ and hence
$\operatorname{\mathbb{E}}\widehat{F}_{n}(x)=F_{m,P}(x)$. In view of
$\widehat{F}_{\pi,n}^{(l)}$ as an empirical distribution function, the
inequality due to Dvoretzky et al. (1956) and Massart (1990) states that

$\operatorname{\mathbb{P}}\left(\|\widehat{F}_{\pi,n}^{(l)}-F_{m,P}\|_{\infty}>t\right)\leq
2\exp\left(-2\lfloor n/m\rfloor t^{2}\right),\quad t\geq 0.$

This inequality does not require $F_{m,P}$ to be continuous; see Massart
(1990, Comment 2(iii)). The inequality implies that
$\|\widehat{F}_{\pi,n}^{(l)}-F_{m,P}\|_{\infty}$ is $c_{0}/(2\sqrt{\lfloor
n/m\rfloor})$-sub-Gaussian for a universal constant $c_{0}>0$ (Vershynin,
2018, Proposition 2.5.2). From Eq. 22, it follows that

$\|\widehat{F}_{n}-F_{m,P}\|_{\infty}\leq\frac{1}{JL}\sum_{l=1}^{L}\sum_{\pi\in\\{\pi_{1},\dots,\pi_{J}\\}}\|\widehat{F}_{\pi,n}^{(l)}-F_{m,P}\|_{\infty},$

where we observe that every $\|\widehat{F}_{\pi,n}^{(l)}-F_{m,P}\|_{\infty}$
is identically distributed. By Jensen’s inequality, we have that the upper
bound above is also $c_{0}/(2\sqrt{\lfloor n/m\rfloor})$-sub-Gaussian and
hence so is $\|\widehat{F}_{n}-F_{m,P}\|_{\infty}$, from which the result
follows. ∎

Recall that $\widehat{G}_{n}(x):=\mathbb{F}_{\\{\widehat{S}_{b}\\}}(x)$ is the
natural subsampling estimation of $G_{P}$, the asymptotic distribution
function of $S_{n}$. For $b=1,\dots,B$, $\widehat{S}_{b}$ is the result of
applying $S$ to the $b$-th row of $\widehat{\mathbf{H}}$. Let $G_{n,P}$ be the
distribution function of $S_{n}$. We define

$K_{n}:=\sqrt{\frac{\lfloor n/m\rfloor}{\log(n/m)}},$ (23)

which is a positive sequence that tends to infinity.

######  Lemma C.2 (High probability bound on subsampling approximation).

Let $T_{n}^{(1)},\dots,T_{n}^{(L)}$ be exchangeable statistics under
$(X,\Omega)\sim P^{n}\times P_{\Omega}$. Given any distribution functions $F$
and $G$, with probability at least $1-4(n/m)^{-2/c_{0}^{2}}$, it holds that

$\displaystyle\|\widehat{F}_{n}-F\|_{\infty}$ $\displaystyle\leq
1/K_{n}+\|F_{m,P}-F\|_{\infty},$ $\displaystyle\|\widehat{G}_{n}-G\|_{\infty}$
$\displaystyle\leq 1/K_{n}+\|G_{m,P}-G\|_{\infty}.$

  * Proof.

Taking $t=1/K_{n}$ in Lemma C.1, by the triangle inequality, the first
inequality above holds with probability at least $1-2(n/m)^{-2/c_{0}^{2}}$. By
the definition of $\widehat{G}_{n}$, a tail bound of the same type as Lemma
C.1 holds for $\|\widehat{G}_{n}-G_{m,P}\|_{\infty}$. Hence, the second
inequality above also holds with probability at least
$1-2(n/m)^{-2/c_{0}^{2}}$. The result follows from a union bound. ∎

######  Corollary C.1.