# The Projected Covariance Measure for assumption-lean variable significance
testing

Anton Rask Lundborg†⋄ †Statistical Laboratory, University of Cambridge, United
Kingdom
⋄Department of Mathematical Sciences, University of Copenhagen, Denmark
⋆Department of Statistics and Data Science, Department of Applied Statistics,
Yonsei University, South Korea Ilmun Kim†⋆ †Statistical Laboratory, University
of Cambridge, United Kingdom
⋄Department of Mathematical Sciences, University of Copenhagen, Denmark
⋆Department of Statistics and Data Science, Department of Applied Statistics,
Yonsei University, South Korea Rajen D. Shah† and Richard J. Samworth†
†Statistical Laboratory, University of Cambridge, United Kingdom
⋄Department of Mathematical Sciences, University of Copenhagen, Denmark
⋆Department of Statistics and Data Science, Department of Applied Statistics,
Yonsei University, South Korea

###### Abstract

Testing the significance of a variable or group of variables $X$ for
predicting a response $Y$, given additional covariates $Z$, is a ubiquitous
task in statistics. A simple but common approach is to specify a linear model,
and then test whether the regression coefficient for $X$ is non-zero. However,
when the model is misspecified, the test may have poor power, for example when
$X$ is involved in complex interactions, or lead to many false rejections. In
this work we study the problem of testing the model-free null of conditional
mean independence, i.e. that the conditional mean of $Y$ given $X$ and $Z$
does not depend on $X$. We propose a simple and general framework that can
leverage flexible nonparametric or machine learning methods, such as additive
models or random forests, to yield both robust error control and high power.
The procedure involves using these methods to perform regressions, first to
estimate a form of projection of $Y$ on $X$ and $Z$ using one half of the
data, and then to estimate the expected conditional covariance between this
projection and $Y$ on the remaining half of the data. While the approach is
general, we show that a version of our procedure using spline regression
achieves what we show is the minimax optimal rate in this nonparametric
testing problem. Numerical experiments demonstrate the effectiveness of our
approach both in terms of maintaining Type I error control, and power,
compared to several existing approaches.

##  1 Introduction

Understanding the relationship between a response and associated predictors is
one of the most common problems faced by data analysts across many diverse
areas of science and industry. Often an important step in this task is to
determine which variables or groups of variables are important in this
relationship. To fix ideas, consider data formed of independent copies of a
triple $(X,Y,Z)$, where $Y\in{\mathbb{R}}$ is our response, and we wish to
assess the significance of a group of predictors $X\in{\mathbb{R}}^{d_{X}}$
after adjusting for confounding variables $Z\in{\mathbb{R}}^{d_{Z}}$; we will
consider a more general setting later in this paper where $X$ and $Z$ can be
potentially non-Euclidean. One simple but popular way of addressing this
problem is to fit a linear regression model
$Y=X^{\top}\beta+Z^{\top}\gamma+\varepsilon$, where we assume that the random
error $\varepsilon$ satisfies ${\mathbb{E}}(\varepsilon\,|\,X,Z)=0$, and
perform an $F$-test for the significance of $X$ (i.e. test the null hypothesis
that $\beta=0$). However, in the case that the linear model is not a
sufficiently good approximation of the ground truth, this can result in
wrongly declaring $X$ to be important or unimportant, and other significance
tests based on parametric models suffer from similar issues. The fact that
regressions based on parametric models are typically greatly outperformed by
modern machine learning methods such as deep learning (Goodfellow et al.,,
2016) and random forests (Breiman,, 2001) in regression competitions such as
those hosted by Kaggle (Bojer and Meldgaard,, 2021), suggests that such
parametric models giving poor approximations to the truth is the norm rather
than the exception, at least in contemporary datasets of interest.

In this work we consider the model-free null hypothesis of conditional mean
independence, that is ${\mathbb{E}}(Y\,|\,X,Z)={\mathbb{E}}(Y\,|\,Z)$; in
words, $X$ does not feature in the regression function of $Y$ on $X$ and $Z$.
It is interesting to compare this to the conditional independence null
$Y\mbox{${}\perp\mkern-11.0mu\perp{}$}X\,|\,Z$, which has attracted much
attention in recent years. The latter asks not just for the regression
function to be expressed as a function of $Z$ alone, but in fact for the
entire conditional distribution of $Y$ given $(X,Z)$ to equal the conditional
distribution of $Y$ given $Z$. Any valid test of conditional mean independence
may be used as a test for conditional independence as its size is no larger
than its size over the larger null hypothesis of conditional mean
independence. The two nulls in fact coincide when $Y$ is binary, but more
generally there are important differences. One attractive property of the
conditional mean independence null is that the alternative of conditional mean
dependence may be characterised by the property that $X$ can improve the
prediction of $Y$ in a mean-squared error sense, given knowledge of $Z$. For
example, consider the setting where $X$ is a binary treatment variable, $Z$
contains all pre-treatment confounders and $Y$ is the observed outcome. Under
assumptions (including the absence of unmeasured confounders) that are
standard in the causal inference literature (Neyman,, 1923; Rubin,, 1974),
conditional mean dependence is equivalent to the existence of a subgroup
average treatment effect, that is a (measurable) subset
$\mathcal{A}\subseteq{\mathbb{R}}^{d_{Z}}$ where
${\mathbb{E}}\\{{\mathbb{E}}(Y\,|\,Z,X=1)\,|\,Z\in\mathcal{A}\\}>{\mathbb{E}}\\{{\mathbb{E}}(Y\,|\,Z,X=0)\,|\,Z\in\mathcal{A}\\}$.
On the other hand, rejection of the conditional independence null does not in
general have an immediate interpretation in terms of its predictive
implications.

Despite the attractions of conditional mean independence, an important issue
is that this property is not testable without further restrictions on the null
hypothesis: if $(X,Y,Z)$ have a density that is absolutely continuous with
respect to Lebesgue measure, then the power of any test at any alternative is
at most its size. This comes as a direct consequence of the untestability of
the smaller conditional independence null (Shah and Peters,, 2020). The
conclusion is that in order to test conditional mean independence, one must
further constrain the null hypothesis in some way.

Given the success of machine learning methods in prediction problems, a
natural and convenient way to specify these constraints is based on
restricting the set of nulls to those where user-chosen regression methods can
estimate certain conditional expectations sufficiently well. One strategy, as
adopted in the _Generalised Covariance Measure_ (GCM) of Shah and Peters,
(2020), involves, in the case where $X$ is univariate, regressing each of $X$
and $Y$ on $Z$, computing the covariance between the resulting residuals and
estimating a normalised version of
${\mathbb{E}}\\{{\mathrm{Cov}}(X,Y\,|\,Z)\\}$, a quantity that is zero under
conditional independence. A drawback of this approach, however, is that it has
no power against alternatives to conditional mean independence where
${\mathbb{E}}\\{{\mathrm{Cov}}(X,Y\,|\,Z)\\}=0$.

To gain greater power, Shah and Peters, (2020) suggest to apply the above with
$X$ replaced by each component of
$\bigl{(}\phi_{1}(X,Z),\ldots,\phi_{m}(X,Z)\bigr{)}$, where
$\phi_{1},\ldots,\phi_{m}:{\mathbb{R}}^{d_{X}\times d_{Z}}\to{\mathbb{R}}$ are
a fixed user-chosen collection of transformations of the data. One may then
base a final test on the maximum absolute value of the resulting test
statistics. It is however not clear how one should choose these
transformations, and if $m$ is large, or indeed $d_{X}$ is large and we use
the strategy above but with the $\phi_{j}$ simply extracting the $j$th
component of $X$, then performing all the regressions involved may be
impractical. A related approach to improve the power properties of the GCM is
introduced by Scheidegger et al., (2021), who propose a carefully-weighted
version of the GCM that, under conditions, can have power against alternatives
where we do not have ${\mathrm{Cov}}(X,Y\,|\,Z)=0$ almost surely; see also
Fernández and Rivera, (2022). Nevertheless, it is perfectly possible to have
${\mathrm{Cov}}(X,Y\,|\,Z)=0$ under conditional mean dependence, and here even
the weighted GCM would be powerless: for example, consider the simple setting
where $(X,Z,\varepsilon)\sim N_{3}(0,\boldsymbol{I})$ and
$Y=X^{2}+\varepsilon$. In this case,
${\mathrm{Cov}}(X,Y\,|\,Z)={\mathrm{Cov}}(X,Y)=0$ despite $X$ clearly being
important for the prediction of $Y$. It is therefore of great interest to
develop methods for testing conditional mean independence whose validity, as
in the case of the GCM and its weighted version, relies primarily on the
predictive properties of user-chosen regression methods, but which have power
against much wider classes of alternatives.

While there has a great deal of research effort on the problem of conditional
independence testing in recent years (we review the contributions most
relevant to our work here in Section 1.2), there has been comparatively little
on testing conditional mean independence. One compelling approach is based on
an equivalent way of stating the null hypothesis: defining

$\displaystyle\tau$
$\displaystyle:={\mathbb{E}}\bigl{[}\bigl{\\{}{\mathbb{E}}(Y\,|\,X,Z)-{\mathbb{E}}(Y\,|\,Z)\bigr{\\}}^{2}\bigr{]}$
(1)

we have that $\tau=0$ if and only if $Y$ is conditionally mean independent of
$X$ given $Z$. This suggests a potential strategy for assessing conditional
mean independence via the estimation of $\tau$. Such an approach was adopted
by Williamson et al., (2021), who employed a plug-in estimator of $\tau$, and
showed that, under conditions, it yields a semiparametric efficient estimator,
provided that $\tau>0$. However, as highlighted by Williamson et al., (2021),
under the null where $\tau=0$, semiparametric approaches such as this face a
fundamental difficulty as the influence function is identically zero, and as a
consequence the test statistic has a degenerate distribution.

To circumvent this issue, Williamson et al., (2022) and independently Dai et
al., (2022), utilise an alternative representation of the target parameter
$\tau={\mathbb{E}}[\\{Y-{\mathbb{E}}(Y\,|\,Z)\\}^{2}]-{\mathbb{E}}[\\{Y-{\mathbb{E}}(Y\,|\,X,Z)\\}^{2}]$
and propose a testing procedure via sample splitting, where estimation of
${\mathbb{E}}[\\{Y-{\mathbb{E}}(Y\,|\,Z)\\}^{2}]$ and
${\mathbb{E}}[\\{Y-{\mathbb{E}}(Y\,|\,X,Z)\\}^{2}]$ is done on independent
splits of the data. This restores the asymptotic normality of the test
statistic under the null, but comes with a significant power loss. In
particular, the resulting test becomes asymptotically powerless if
$\sqrt{n}\tau\rightarrow 0$, even for a parametric linear model where the
optimal testing rate is known to be of order $n^{-1}$. Moreover, the
asymptotic normality fails when $Y$ is (close to) independent of $(X,Z)$,
which raises concerns about uniform validity of the test; see Appendix B.2 for
details.

###  1.1 Outline of our approach and contributions

In view of the considerations above, the goal of this paper is to propose a
new framework for testing conditional mean independence that has the following
properties:

  * •

Flexible Type I error control. The user should be able to leverage flexible
regression methods to ensure validity of the test uniformly over classes of
distributions where these methods perform sufficiently well.

  * •

Rate-optimal power in diverse settings. The test should have minimax rate-
optimal power in both simple parametric models, as well as challenging
nonparametric settings, when used with appropriate regression methods.

  * •

Computationally practical. The test should involve performing only a small
number of regressions.

Our approach is based on the following alternative characterisation of
conditional mean independence: $Y$ is conditionally mean independent of $X$
given $Z$ if and only if

${\mathbb{E}}\bigl{[}\bigl{\\{}Y-{\mathbb{E}}(Y\,|\,Z)\bigr{\\}}f(X,Z)\bigr{]}={\mathbb{E}}\bigl{[}{\mathrm{Cov}}\bigl{(}Y,f(X,Z)\,|\,Z\bigr{)}\bigr{]}=0$
(2)

for all functions $f$ such that
${\mathbb{E}}\bigl{(}f(X,Z)^{2}\bigr{)}<\infty$. In words, the residuals
$Y-{\mathbb{E}}(Y\,|\,Z)$ from regressing $Y$ on $Z$ alone are uncorrelated
with any square-integrable function of $X$ and $Z$. On the other hand, under
an alternative, these residuals should not be pure noise but contain some
‘signal’ that can be exposed via an appropriate $f$ such that the left-hand
side of (2) is strictly positive.

To motivate our specific strategy, consider an oracular test statistic that
uses knowledge of the conditional expectation ${\mathbb{E}}(Y\,|\,Z)$: given
independent copies $(X_{i},Y_{i},Z_{i})_{i=1}^{n}$ of $(X,Y,Z)$ and a function
$f$, the random variables
$L_{i}^{*}:=\\{Y_{i}-{\mathbb{E}}(Y_{i}\,|\,Z_{i})\\}f(X_{i},Z_{i})$ for
$i=1,\ldots,n$ are independent and identically distributed, with zero mean
under the null. Writing
$\widetilde{L}_{i}^{*}:=\\{Y_{i}-{\mathbb{E}}(Y_{i}\,|\,X_{i},Z_{i})\\}f(X_{i},Z_{i})$,
we have that under regularity conditions, the studentised statistic

$\displaystyle
T^{*}:=\frac{\frac{1}{\sqrt{n}}\sum_{i=1}^{n}L_{i}^{*}}{\sqrt{\frac{1}{n}\sum_{i=1}^{n}\widetilde{L}_{i}^{*2}}}$
(3)

converges to a standard normal distribution under the null, and may thus form
the basis of a test. Note that since $\widetilde{L}_{i}^{*}=L_{i}^{*}$ under
the null, we may alternatively studentise the test statistic using the
empirical standard deviation of the $L_{i}^{*}$; however this version
simplifies the derivation to follow.

Different choices of $f$ would lead to different power properties under an
alternative. Ideally, we want to maximise the value of the test statistic
under an alternative, so we would like
${\mathbb{E}}(L^{*}_{i})/\sqrt{{\mathrm{Var}}(\widetilde{L}_{i}^{*})}$ to be
as large as possible. It may be shown (see Proposition 19 in Section B.2 of
the appendix) that this is uniquely maximised, up to an arbitrary positive
scaling, by choosing $f(X,Z)=h(X,Z)/v(X,Z)$, where
$h(X,Z):={\mathbb{E}}(Y\,|\,X,Z)-{\mathbb{E}}(Y\,|\,Z)$ and
$v(X,Z):={\mathrm{Var}}(Y\,|\,X,Z)$. We therefore see that the optimal $f$ is
a version of the projection $h$ of $Y$ onto the space of square-integrable
functions of $(X,Z)$ that are orthogonal to functions of $Z$, inversely
weighted by the conditional variance $v$.

The considerations above suggest the following approach: use one portion of
the data to obtain an estimate $\widehat{f}$ of the projection $f$, and then
use the remaining data to evaluate a test statistic of the form (3), with the
unknown conditional expectations there replaced with appropriate regression
estimates. This forms the basis of our proposed test statistic, which we call
the _Projected Covariance Measure_ (PCM). In fact, it turns out to be
advantageous to modify somewhat the basic blueprint described above, for
instance by subtracting from $\widehat{f}(X,Z)$ an estimate of its conditional
expectation given $Z$, to reduce bias; a complete description of our
methodology is given in Section 2.

One important issue to be addressed is the fact that under the null, $h$ is
the zero function, and as a consequence, both the numerator and denominator of
$T^{*}$ are zero. This is not immediately problematic for the oracular
statistic $T^{*}$, as one can always decide not to reject the null when the
numerator is precisely $0$. However, it might appear to be potentially
disastrous for an empirical version of $T^{*}$, where any bias terms in the
numerator could be inflated by division with a denominator that is close to
zero. One of our main contributions in this work is to show that by
formulating our PCM test statistic appropriately, it has an asymptotic
standard Gaussian limit in settings ranging from low- and high-dimensional
linear models to fully nonparametric settings. Moreover, we demonstrate
empirically that this limiting behaviour can be expected to hold more
generally when machine learning methods such as random forests (Breiman,,
2001) are used for the regressions involved.

The rest of the paper is organised as follows. After reviewing some related
literature in Section 1.2, we present a full description of our PCM
methodology in Section 2. In Section 3, we examine the simplest instantiation
of our general framework and study testing in the context of low-dimensional
linear models. An important revelation of this analysis is that in contrast to
the equally general testing frameworks of Williamson et al., (2022) and Dai et
al., (2022), our approach has power against local alternatives where $\tau$ is
of order $n^{-1}$. We go on to show that, under conditions, the PCM maintains
Type I error control in high-dimensional linear models, even when using an
essentially arbitrary machine learning method to estimate the projection $f$.
We present a general theory of the PCM in Section 4, giving conditions
involving prediction errors of the user-chosen regression procedures used in
the PCM that result in Type I error control, and also study the power of the
procedure. In Section 5, we show how our general conditions for Type I error
control may be satisfied in a fully nonparametric regression setting when
using series estimators for the relevant regressions. We also introduce a
slight variant of our approach involving additional sample splitting that
enjoys what we show to be minimax rate optimal power over classes of
alternatives for which $\tau$ in (1) satisfies a lower bound.

In Section 6, we conduct several simulation experiments that demonstrate the
effectiveness of the PCM when used with generalised additive model-based
regressions (Wood,, 2017) and random forests, in terms of both Type I error
control and power. We conclude with a discussion in Section 7 outlining
potential future research directions suggested by our work.

In Sections A and B of the appendix, we include the proofs of all our main
results and related auxiliary lemmas. Section C provides a self-contained
description of spline regression and related results that we use for our
analysis in Section 5. In Section D, we give a more detailed analysis of our
results for linear projections in Section 3; in particular we derive an exact
asymptotic power function of our test. Section E contains the results from
additional numerical experiments beyond those included in Section 6.

###  1.2 Literature review

There is a relatively small body of literature that is explicitly concerned
with conditional mean independence. Early developments on this topic include
the work of Fan and Li, (1996), Lavergne and Vuong, (2000) and Aït-Sahalia et
al., (2001) from the econometrics community. Jin et al., (2018) propose an
approach for testing conditional mean independence in cases where
$\mathbb{E}(Y\,|\,Z)$ is a linear function of $Z$, based on the martingale
difference divergence proposed by Shao and Zhang, (2014).

Recent years have witnessed an increasing use of machine learning (ML) tools
for statistical inference. For example, Chernozhukov et al., (2018) introduce
an ML-driven approach for estimating causal parameters in the presence of
complex nuisance parameters. Shah and Bühlmann, (2018) and Janková et al.,
(2020) propose methods for goodness-of-fit testing in high-dimensional
(generalised) linear models that involve detecting remaining signal in
residuals using ML methods. More closely related to this work, Williamson et
al., (2022), and Dai et al., (2022), propose model-free methods for assessing
conditional mean independence that can take advantage of existing ML
algorithms. Williamson et al., (2022) derive a semiparametrically efficient
estimator $\tau$, but recognise the difficulty of testing the null hypothesis
that $\tau=0$ caused by the fact that the efficient influence function is
identically zero under the null. This means that their sample-splitting
approach lacks validity when $(X,Y,Z)$ are independent, and moreover it turns
out that the test may require larger values of $\tau$ than necessary in order
to achieve power; see Section 3.1 for a more detailed discussion. Dai et al.,
(2022) alleviate the Type I error issue by adding noise to their test
statistic, but this comes at a further price in terms of power, as pointed out
by Verdinelli and Wasserman, (2021). Cai et al., (2022) also propose model-
free tests of conditional mean independence; one of their test statistics has
a similar form to the one in Williamson et al., (2021) that compares the
predictive performance of two regression models, and they use the permutation
approach to calculate a $p$-value. Another related work is that of Zhang and
Janson, (2020), who provide a method of constructing confidence intervals for
$\tau$ in the case where the conditional distribution of $X$ given $Z$ is
(almost) known; see also Candès et al., (2018) and Berrett et al., (2020), who
employ similar assumptions in the context of testing conditional independence.

Many existing tests, including ours, determine their critical values based on
asymptotic theory derived under the null. However, most work (implicitly)
targets pointwise Type I error control that holds only each fixed null. This
type of pointwise analysis leaves room for the existence of a sequence of null
distributions for which the Type I error can be made arbitrarily large. A
classical example is the fact that the $t$-test that has pointwise asymptotic
size $\alpha$ for the class of distributions with finite variance, has uniform
asymptotic size 1 for the same class of distributions (Romano,, 2004). While
it is straightforward to introduce, for instance, moment conditions to restore
uniform size control in that problem, we argue that the issue is even more
pertinent in the context of testing conditional (mean) independence as there
are no canonical choices of restrictions to the null that can yield this form
of error control. In this work,

we therefore put great emphasis on uniform Type I error control over classes
of distributions in order to present more practically-relevant error
guarantees. This uniform analysis is in line with recent work on conditional
independence testing such as Shah and Peters, (2020), Petersen and Hansen,
(2021), Lundborg et al., (2022), Scheidegger et al., (2021) and Neykov et al.,
(2021).

Our work builds on a classical technique, namely sample splitting, that
involves partitioning the data into disjoint subsamples for different
purposes: roughly speaking, a portion of the data is used for seeking a good
direction that potentially contains a high signal and the other portion is
used for conducting a test based on the data projected along the given
direction. Cox, (1975) is one of the earliest papers that applies sample
splitting to testing problems. Since then, many inference procedures have been
developed by leveraging a similar technique to perform variable selection in
high-dimensional models (Wasserman and Roeder,, 2009; Meinshausen et al.,,
2009; Meinshausen and Bühlmann,, 2010; Shah and Samworth,, 2013), inference
after model selection (Rinaldo et al.,, 2019), changepoint detection (Wang and
Samworth,, 2018) and inference based on maximum likelihood estimators
(Wasserman et al.,, 2020), to name just a few. In a similar vein, Kim and
Ramdas, (2020) introduce splitting-based procedures that address an issue of
degenerate $U$-statistics for high-dimensional inference. While our main focus
is on testing, sample splitting has also been considered for estimation
problems, where it typically works as a device to reduce a bias and thus help
to obtain a fast (often optimal) convergence rate (Chernozhukov et al.,, 2018;
Newey and Robins,, 2018; Wang and Shah,, 2020). Some parts of our work are
motivated by Newey and Robins, (2018), who propose cross-fit estimators of
functionals involving conditional expectations.

###  1.3 Preliminaries and notation

Throughout this paper, we adopt the convention that $0/0:=0$ and let
$\mathrm{sgn}(\cdot)$ denote the sign function on $\mathbb{R}$, with the
convention that $\mathrm{sgn}(0):=0$. We let $x\wedge y:=\min(x,y)$ and for
$n\in\mathbb{N}$, let $[n]:=\\{1,\dots,n\\}$. For two sequences $(a_{n})$ and
$(b_{n})$, we write $a_{n}\asymp b_{n}$ and $a_{n}\lesssim b_{n}$ to mean that
there exist $c,C>0$ such that $0<c\leq|a_{n}/b_{n}|\leq C<\infty$ for all $n$,
and $a_{n}\leq Cb_{n}$ for all $n$ respectively. For a vector
$x\in\mathbb{R}^{n}$ and $p\in[1,\infty]$, we denote its $\ell_{p}$ norm by
$\|x\|_{p}$. The operator norm of a matrix
$\boldsymbol{A}\in\mathbb{R}^{n\times m}$ is denoted by
$\|\boldsymbol{A}\|_{\mathrm{op}}$ and the maximum and minimum eigenvalues of
a symmetric matrix $\boldsymbol{B}\in\mathbb{R}^{n\times n}$ are denoted by
$\lambda_{\max}(\boldsymbol{B})$ and $\lambda_{\min}(\boldsymbol{B})$,
respectively. We use the notation $z_{\alpha}$ to denote the $\alpha$th
quantile of the standard normal distribution, whose distribution function is
denoted by $\Phi$.

In order to present our uniform results on testing, we require some
conventions for probabilistic notation used in what follows. Let
$(\Omega,\mathcal{F})$ be a measurable space equipped with a family of
probability measures $(\mathbb{P}_{P})_{P\in\mathcal{P}}$ where $\mathcal{P}$
is a collection of distributions on a Euclidean space. We will permit the
family $\mathcal{P}$ to depend on $n$, to allow for settings where the number
of parameters grows with $n$, but will typically suppress this in the
notation. Given a family of sequences of random variables
$(X_{P,n})_{P\in\mathcal{P},n\in\mathbb{N}}$ on ($\Omega,\mathcal{F}$) whose
distributions are determined by $P\in\mathcal{P}$, we write
$X_{P,n}=o_{\mathcal{P}}(1)$ if
$\sup_{P\in\mathcal{P}}{\mathbb{P}}_{P}(|X_{P,n}|>\epsilon)\rightarrow 0$ for
every $\epsilon>0$. Similarly, we write $X_{P,n}=O_{\mathcal{P}}(1)$ if, for
any $\epsilon>0$, there exist $M_{\epsilon},N_{\epsilon}>0$ such that
$\sup_{n\geq
N_{\epsilon}}\sup_{P\in\mathcal{P}}{\mathbb{P}}_{P}(|X_{P,n}|>M_{\epsilon})<\epsilon$.
In addition, for another family of sequences of random variables
$(Y_{P,n})_{P\in\mathcal{P},n\in\mathbb{N}}$, we write
$X_{P,n}=o_{\mathcal{P}}(Y_{{P,n}})$ if there exists $R_{P,n}$ with
$X_{P,n}=Y_{{P,n}}R_{P,n}$; likewise, we write
$X_{P,n}=O_{\mathcal{P}}(Y_{P,n})$ if $R_{P,n}=O_{\mathcal{P}}(1)$ in this
representation. We say that $(X_{P,n})_{P\in\mathcal{P},n\in\mathbb{N}}$
converges uniformly in distribution to random variable $X$ with distribution
function $F$ if for all continuity points $x$ of $F$, we have

$\lim_{n\to\infty}\sup_{P\in\mathcal{P}}\bigl{|}{\mathbb{P}}_{P}(X_{P,n}\leq
x)-F(x)\bigr{|}=0.$

We will denote different independent datasets by
$\mathcal{D}_{1},\mathcal{D}_{2},\ldots$, each containing $n$ independent
observations. We will frequently abuse notation and write conditional
expectations conditioning on a random function, e.g.
${\mathbb{E}}_{P}\bigl{(}\widehat{f}(X,Z)\,|\,\widehat{f},Z\bigr{)}$ where
$\widehat{f}$ is a function produced by some regression estimator. By this we
mean formally that we condition on the sample used to construct the regression
estimator and any additional randomness involved in the computation of the
regression function. We let $(X,Y,Z)$ be random variables in
$\mathcal{X}\times\mathbb{R}\times\mathcal{Z}$, where $\mathcal{X}$ and
$\mathcal{Z}$ are measurable spaces, although we will at times think of
$\mathcal{X}$ and $\mathcal{Z}$ being specific $d_{X}$\- and
$d_{Z}$-dimensional Euclidean spaces, respectively.

##  2 Projected covariance measure

In this section, we outline our PCM methodology in detail. We first provide
further motivation in Section 2.1, before presenting our final algorithm in
Section 2.2. Given that our approach involves sample splitting, it is
convenient to assume here and also throughout Sections 3 and 4 that we have
$2n$ independent and identically distributed observations
$(X_{i},Y_{i},Z_{i})_{i=1}^{2n}$ rather than the conventional $n$
observations.

###  2.1 Motivation

Recall that the approach sketched in Section 1.1 involves first computing an
estimate $\widehat{f}$ of the weighted projection

$f(X,Z)=\frac{h(X,Z)}{v(X,Z)}=\frac{{\mathbb{E}}(Y\,|\,X,Z)-{\mathbb{E}}(Y\,|\,Z)}{{\mathrm{Var}}(Y\,|\,X,Z)}$

using one portion of the data, say
$\mathcal{D}_{2}:=(X_{i},Y_{i},Z_{i})_{i=n+1}^{2n}$. We discuss how to
construct the estimate $\widehat{f}$ in Section 2.2. Next, given an estimate
$\widehat{m}(\cdot)$ of $m(\cdot):={\mathbb{E}}(Y\,|\,Z=\cdot)$, the oracular
test statistic (3) suggests a numerator of our test statistic of the form

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\\{Y_{i}-\widehat{m}(Z_{i})\\}\widehat{f}(X_{i},Z_{i}).$
(4)

We would like this to have mean close to zero under the null; however it is
well-known (Chernozhukov et al.,, 2018) that when using a nonparametric
estimator $\widehat{m}$, the quantity above may carry a substantial bias and
we should instead consider an orthogonalised version of the form

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}L_{i}\qquad\text{with}\qquad
L_{i}:=\\{Y_{i}-\widehat{m}(Z_{i})\\}\\{\widehat{f}(X_{i},Z_{i})-\widehat{m}_{\widehat{f}}(X_{i},Z_{i})\\},$

where $\widehat{m}_{\widehat{f}}$ is an estimate of
$m_{\widehat{f}}(\cdot):={\mathbb{E}}(\widehat{f}(X,Z)\,|\,Z=\cdot,\widehat{f})$.
Importantly, the bias term then involves a product of the mean squared
prediction error (MSPE) of $\widehat{m}$,

$\frac{1}{n}\sum_{i=1}^{n}\\{m(Z_{i})-\widehat{m}(Z_{i})\\}^{2},$ (5)

and that of $\widehat{f}$, a quantity that may be substantially smaller than
the MSPE of $\widehat{m}$ alone (which would drive the bias in (4)).

Turning to the denominator of our test statistic, instead of studentising by a
quantity requiring an estimate of ${\mathbb{E}}(Y\,|\,X,Z)$ as suggested by
(3), it is practically more convenient to normalise using the empirical
standard deviation of $L_{1},\ldots,L_{n}$ as this does not involve performing
an additional regression. Thus we propose to take as our test statistic

$T:=\frac{\frac{1}{\sqrt{n}}\sum_{i=1}^{n}L_{i}}{\sqrt{\frac{1}{n}\sum_{i=1}^{n}L_{i}^{2}-\bigl{(}\frac{1}{n}\sum_{i=1}^{n}L_{i}^{2}\bigr{)}}}.$
(6)

For local alternatives, both versions are near-identical and so any
differences in power properties should be very slight, as we have also
observed empirically. We choose in practice to train $\widehat{m}$ and
$\widehat{m}_{\widehat{f}}$ on $\mathcal{D}_{1}$ rather than
$\mathcal{D}_{2}$. The errors such as (5) that are required to be controlled
are then _in-sample errors_ , that is the regression methods are trained on
the same data on which they are evaluated, so the regression methods need not
extrapolate to unseen data points, for example. While from a theoretical
perspective in-sample errors and out-of-sample errors are often thought of
similarly, in finite samples, these can behave differently: for example Figure
1 demonstrates that when using additive models (computed using the R package
mgcv (Wood,, 2017)) to estimate $m$ in a setup considered in Section 6.1, out-
of-sample errors can be appreciably larger with non-negligible probability.

Figure 1: In-sample and out-of-sample errors for the model where
$(Z_{1},\dots,Z_{7})\sim N_{7}(0,\boldsymbol{I})$, $Y=\sin(2\pi
Z_{1})+\varepsilon$ with $\varepsilon\sim N(0,1)$ independently of
$Z_{1},\dots,Z_{7}$, and regressions are performed using mgcv; see Section 6.1
for more details on this setup.

As the PCM may be thought of as the GCM applied to a transformed $X$, we would
hope to obtain a standard Gaussian limit for $T$ as in the case of the regular
GCM test statistic. Given that the transformation is designed to result in
large values of $T$ under an alternative, we would perform a one-sided test by
rejecting when $T$ exceeds the appropriate normal quantile. Unfortunately
however, the theory that guarantees asymptotic validity of

the GCM test statistic does not apply in our case:

it would require
${\mathrm{Var}}(\\{Y-m(Z)\\}\\{\widehat{f}(X,Z)-m_{\widehat{f}}(Z)\\}\,|\,\widehat{f})$,
i.e. the (square of the) target of the denominator to be bounded away from
zero under the null. But $f$ is identically $0$ under the null, so
$\widehat{f}$ and hence the above variance, and also both numerator and
denominator of our test statistic, should all converge to $0$.

To see why we can expect a standard Gaussian limit for our test statistic
despite this apparent degeneracy, consider a linear model setting where
$(X,Y,Z)\in{\mathbb{R}}\times{\mathbb{R}}\times{\mathbb{R}}^{d}$ are related
through

$Y=\beta X+Z^{\top}\boldsymbol{\gamma}+\varepsilon\qquad\text{and}\qquad
X=Z^{\top}\boldsymbol{\eta}+\xi,$ (7)

with $\beta=0$ and
${\mathbb{E}}(\varepsilon\,|\,Z)={\mathbb{E}}(\xi\,|\,Z)=0$. If we form
estimates $\widehat{h}$ and $\widehat{m}$ using ordinary least squares, and
for simplicity set $\widehat{v}\equiv 1$ when forming $\widehat{f}$, then
$\widehat{f}(x,z)=\widehat{h}(x,z)$ takes the form
$\widehat{\beta}x+z^{\top}\widetilde{\boldsymbol{\delta}}$ for some
$(\widehat{\beta},\widetilde{\boldsymbol{\delta}})\in{\mathbb{R}}\times{\mathbb{R}}^{d}$.
Note that both $\widehat{\beta}$ and $\|\widetilde{\boldsymbol{\delta}}\|_{2}$
are of stochastic order $1/\sqrt{n}$.

Let us write $\widehat{\boldsymbol{\gamma}}$ and $\widehat{\boldsymbol{\eta}}$
for the regression estimates of $\boldsymbol{\gamma}$ and $\boldsymbol{\eta}$
respectively. The next step of our procedure involves regressing each of
$(Y_{i})_{i=1}^{n}$ and $\bigl{(}\widehat{f}(X_{i},Z_{i})\bigr{)}_{i=1}^{n}$
onto $(Z_{i})_{i=1}^{n}$. The residuals from the latter regression take the
form
$\widehat{\beta}\\{Z_{i}^{\top}(\boldsymbol{\eta}-\widehat{\boldsymbol{\eta}})+\xi_{i}\\}$,
so in our case

$L_{i}=\widehat{\beta}\\{Z_{i}^{\top}(\boldsymbol{\gamma}-\widehat{\boldsymbol{\gamma}})+\varepsilon_{i}\\}\\{Z_{i}^{\top}(\boldsymbol{\eta}-\widehat{\boldsymbol{\eta}})+\xi_{i}\\}.$

Thus, although $L_{i}$ and hence its standard deviation would be of order
$1/\sqrt{n}$ due to the factor of $\widehat{\beta}$, writing
$L_{i}^{\prime}:=L_{i}/|\widehat{\beta}|$, we see that our test statistic is
of the form $\mathrm{sgn}(\widehat{\beta})T^{\prime}$, where $T^{\prime}$ is a
version of $T$ in (6) with $L_{i}$ replaced by $L_{i}^{\prime}$. But
$L_{i}^{\prime}$ is an order 1 quantity (in contrast of $L_{i}$), so under
mild conditions $n^{-1/2}\sum_{i=1}^{n}L_{i}^{\prime}$ will have a non-
degenerate Gaussian limit, yielding a standard Gaussian limit for
$T^{\prime}$. As $\widehat{\beta}$ is independent of $T^{\prime}$, having been
constructed on $\mathcal{D}_{2}$, the final test statistic $T$ will also
converge to a standard Gaussian.

While this argument provides a heuristic justification for the asymptotic
validity of our proposed test under a simple linear model, there remain
challenges in extending the basic intuition of this example to more general
settings. In the above, it was possible to isolate the randomness from
$\widehat{f}$ simply via the sign of $\widehat{\beta}$, which helps bypass the
$0/0$ issue. However, it is by no means straightforward to deal with the
limits of the form $0/0$ in a nonparametric setting where $\widehat{f}$ is
entangled with other sources of randomness in a complicated way. Moreover, in
nonparametric settings one needs to put more effort into ensuring that the
convergence rates of $\widehat{f}$, $\widehat{m}$ and
$\widehat{m}_{\widehat{f}}$ are fast enough that the bias term is
asymptotically negligible. In this process, we are obliged to handle a nested
regression problem that has rarely been touched in the literature with a few
exceptions (e.g. Kennedy,, 2020).

###  2.2 PCM algorithm

Our PCM approach developed in Section 2.1 is set out in Algorithm 1, with some
recommendations for the constructions of $\widehat{h}$ and $\widehat{v}$ that
we discuss in Sections 2.2.1 and 2.2.2 below. In Section 2.2.3, we then put
forward a version of the PCM using multiple sample splits that we recommend
using in practice.

Algorithm 1 Projected Covariance Measure: single sample split

Input: Data $(X_{i},Y_{i},Z_{i})_{i=1}^{2n}$, significance level
$\alpha\in(0,1)$, partition of $[2n]=\mathcal{I}_{1}\cup\mathcal{I}_{2}$ into
index sets $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$, each of size $n$.
Options: Regression methods for each of the regressions.
Define: $\mathcal{D}_{j}=(X_{i},Y_{i},Z_{i})_{i\in\mathcal{I}_{j}}$ for
$j\in[2]$.

1:Form $\widehat{h}$.

  1. (i)

Regress $Y$ onto $(X,Z)$ using $\mathcal{D}_{2}$ to give fitted regression
function $\widehat{g}$.

  2. (ii)

If $\widehat{g}$ can be modified so that all components involving only $Z$ are
set to $0$, let $\widetilde{g}$ be this modified version of $\widehat{g}$.
Alternatively, set $\widetilde{g}:=\widehat{g}$.

  3. (iii)

Regress $\widetilde{g}(X,Z)$ onto $Z$ using $\mathcal{D}_{2}$ to give fitted
regression function $\widetilde{m}$, and then set
$\widetilde{h}(x,z):=\widetilde{g}(x,z)-\widetilde{m}(z)$.

  4. (iv)

Compute

$\widehat{\rho}:=\frac{1}{n}\sum_{i\in\mathcal{I}_{2}}\bigl{\\{}Y_{i}-\widehat{g}(X_{i},Z_{i})+\widetilde{g}(X_{i},Z_{i})-\widetilde{m}(Z_{i})\bigr{\\}}\widetilde{h}(X_{i},Z_{i}),$

and set $\widehat{h}(x,z):=\mathrm{sgn}(\widehat{\rho})\widetilde{h}(x,z)$,

2:Form $\widehat{v}$.

  1. (i)

Regress $\\{Y-\widehat{g}(X,Z)\\}^{2}$ onto $(X,Z)$ using $\mathcal{D}_{2}$ to
give $\widetilde{v}$.

  2. (ii)

Define $a:[0,\infty)\to[0,\infty]$ by

$a(c):=\frac{1}{n}\sum_{i\in\mathcal{I}_{2}}\frac{\\{Y_{i}-\widehat{g}(X_{i},Z_{i})\\}^{2}}{\max\\{\widetilde{v}(X_{i},Z_{i}),0\\}+c}.$

If $a(0)\leq 1$, set $\hat{c}:=0$; otherwise find $\hat{c}$ by solving
$a(c)=1$. Set $\widehat{v}(x,z):=\max\\{\widetilde{v}(x,z),0\\}+\hat{c}$.

3:Compute test statistic.

  1. (i)

Set $\widehat{f}(x,z):=\widehat{h}(x,z)/\widehat{v}(x,z)$ and regress
$\widehat{f}(X,Z)$ onto $Z$ using $\mathcal{D}_{1}$, giving
$\widehat{m}_{\widehat{f}}$.

  2. (ii)

Regress $Y$ onto $Z$ using $\mathcal{D}_{1}$ to give $\widehat{m}$.

  3. (iii)

For $i\in\mathcal{I}_{1}$ set
$L_{i}:=\\{Y_{i}-\widehat{m}(Z_{i})\\}\\{\widehat{f}(X_{i},Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})\\}$
and let

$T:=\frac{\frac{1}{\sqrt{n}}\sum_{i\in\mathcal{I}_{1}}L_{i}}{\sqrt{\frac{1}{n}\sum_{i\in\mathcal{I}_{1}}L_{i}^{2}-\bigl{(}\frac{1}{n}\sum_{i\in\mathcal{I}_{1}}L_{i}^{2}\bigr{)}}}.$

4:Reject $H_{0}$ if $T>z_{1-\alpha}$.

####  2.2.1 Choice of $\widehat{h}$

We would like $\widehat{h}(X,Z)$ to be close to
$h(X,Z)={\mathbb{E}}(Y\,|\,X,Z)-{\mathbb{E}}(Y\,|\,Z)$ in order to maximise
the power of the procedure. There are several ways of estimating $h$, with
perhaps the most obvious being simply to take the difference of the estimated
regression functions $\widehat{g}$ and $\check{m}$ from regressing $Y$ on each
of $(X,Z)$ and $Z$. An alternative approach is based on observing that
$h(X,Z)=g(X,Z)-{\mathbb{E}}\bigl{(}g(X,Z)\,|\,Z\bigr{)}$ where
$g(X,Z):={\mathbb{E}}(Y\,|\,X,Z)$. This suggests subtracting not $\check{m}$
but the output of regressing $\widehat{g}(X,Z)$ onto $Z$. An advantage of this
latter approach is that we are free to subtract any function $r$ of $Z$ from
$\widehat{g}(X,Z)$ prior to this second regression onto $Z$, as we also have
$h(X,Z)=g(X,Z)-r(Z)-{\mathbb{E}}\bigl{(}g(X,Z)-r(Z)\,|\,Z\bigr{)}$. Thus for
example if $\widehat{g}(x,z)=\widehat{g}_{x}(x)+\widehat{g}_{z}(z)$, then we
may form an estimate of $h(X,Z)$ as the residuals from regressing
$\widehat{g}_{x}(X)$ onto $Z$. This second regression can then focus on
removing any $Z$ signal in $\widehat{g}_{x}(X)$, rather than also having to
cancel out $\widehat{g}_{z}(Z)$. We do not make the claim that this always
makes a large improvement on the first approach, and indeed for certain
regression methods such as ordinary least squares (OLS), both approaches are
identical and the ‘cancellation’ is automatic. Nevertheless, we find the
approach of Step 1 of Algorithm 1 to be a sensible default choice.

In Step 1(iv) we make a final modification to the estimate thus constructed by
potentially flipping its sign. The rationale for this is as follows: under an
alternative, we have that
${\mathbb{E}}[\\{Y-{\mathbb{E}}(Y\,|\,Z)\\}h(X,Z)]=\tau>0$. As a basic check
then, we can see if an empirical version of this inequality, with
$\widehat{h}$ taking place of $h$ and an estimate of ${\mathbb{E}}(Y\,|\,Z)$
replacing the population quantity, holds; if not, we can at least flip the
sign of $\widehat{h}$. This does not require performing any further
regressions to estimate ${\mathbb{E}}(Y\,|\,Z)$: noting the identity
${\mathbb{E}}(Y\,|\,Z)={\mathbb{E}}\\{{\mathbb{E}}(Y\,|\,X,Z)-r(Z)\,|\,Z\\}+r(Z)$,
observe that $\widetilde{m}$ in Step 1(iii) is an estimate of the first of
these quantities with $r(Z)=\widehat{g}(X,Z)-\widetilde{g}(X,Z)$, where
$\widetilde{g}$ is defined in Step 1(ii). When using OLS for each of the
regressions, $\widehat{\rho}$ is guaranteed to be non-negative, so no sign
flip is performed.

In high-dimensional settings, we would typically use a sparsity-inducing
regression method such as the Lasso (Tibshirani,, 1996). Considering the
simple case where $X$ is univariate, this can result in the coefficient for
$X$ being set exactly to zero, and so the recommended construction of
$\widehat{h}$ given above would simply produce the zero function. While not a
problem for Type I error control, as our convention (see Section 1.3) is not
to reject the null when $L_{i}=0$ for all $i$, it is wasteful in terms of
power and a better approach here would be to leave the coefficient for $X$
unpenalised. More generally for multivariate $X$, we can additionally regress
on the first principal component of $X$ for example, and leave this
unpenalised.

####  2.2.2 Choice of $\widehat{v}$

A natural way to form $\widehat{v}$ is to regress the square of the residuals
from regressing $Y$ onto $(X,Z)$, and this is what we recommend in Step 2(i)
of Algorithm 1 to produce $\widetilde{v}$. An issue is that while $v$ is
clearly non-negative, and expected to be positive everywhere, $\widetilde{v}$
may in fact be negative. Equally problematic is the possibility that
$\widetilde{v}$ is very close to $0$ at some $(X_{i},Z_{i})$, as then taking
$\widehat{v}=\widetilde{v}$, we would have $\widehat{f}(X_{i},Z_{i})$ very
large and hence $\widehat{f}(X_{i},Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})$
and $L_{i}$ may be greatly inflated and dominate the test statistic. To
mitigate these problems, we modify $\widetilde{v}$ by taking the positive part
of our initial estimate, and then adding a non-negative constant $\hat{c}$.
This constant is chosen such that $a(\hat{c})$ (see Step 2(ii) of Algorithm 1)
is at most $1$, the rationale coming from the population level identity
${\mathbb{E}}[\\{Y-{\mathbb{E}}(Y\,|\,X,Z)\\}^{2}/v(X,Z)]=1$. We also note
that estimation of the conditional variance $v$ is not critical for good power
properties. For example, in Section 5 we show that simply setting
$\widehat{v}\equiv 1$ delivers minimax rate optimal power in a fully
nonparametric setting; however the power properties may improve empirically by
a constant factor; see Section 3.1.

####  2.2.3 Multiple sample splitting

The single sample split in Algorithm 1 crucially ensures independence between
$\widehat{f}$ and the remaining data $\mathcal{D}_{1}$, but has the
consequence of introducing unwanted additional randomness into the test
statistic. This means that two practitioners with the same data may reach
different conclusions about whether or not they should reject the null if they
use different randomisation seeds. To circumvent this issue, we advocate
applying the single split PCM to multiple splits of the data, and averaging
the resulting test statistics, as summarised in Algorithm 2. An alternative to
working with the averaged test statistic would be to combine the $p$-values of
the individual tests, for which several methods are available, ranging from
twice the average or median of $p$-values to the Bonferroni method (e.g. Vovk
and Wang,, 2020; DiCiccio et al.,, 2020; Meinshausen et al.,, 2009, and
references therein). However, our experience is that these approaches tend to
be overly conservative, and typically lose power compared to considering a
single test. Instead, we propose to compare the averaged test statistic
$\bar{T}$ to a standard Gaussian quantile, as with the single split test
statistic $T$. We expect this to be conservative, as by Jensen’s inequality,
$\bar{T}$ is less than or equal to $T$ in the convex ordering, so for example
${\mathrm{Var}}(\bar{T})\leq{\mathrm{Var}}(T)$. However, in practice it does
tend to improve slightly on the power of a single-split test, while at the
same time having the important benefit of derandomising it.

Algorithm 2 Projected Covariance Measure: multiple sample splits

Input: Data $(X_{i},Y_{i},Z_{i})_{i=1}^{2n}$, significance level
$\alpha\in(0,1)$, number of splits $B$.
Options: Regression methods for each of the regressions.

1:Form complementary pairs of index sets
$\\{(\mathcal{I}_{1}^{(b)},\mathcal{I}_{2}^{(b)}):b\in[B]\\}$ each of size
$n$, where $\mathcal{I}_{1}^{(b)}\cup\mathcal{I}_{2}^{(b)}=[2n]$.

2:For $b\in[B]$, apply Algorithm 1 with index sets
$\mathcal{I}_{1}^{(b)},\mathcal{I}_{2}^{(b)}$ to produce test statistic
$T^{(b)}$.

3:Define $\bar{T}:=\sum_{b=1}^{B}T^{(b)}/B$ and reject $H_{0}$ if
$\bar{T}>z_{1-\alpha}$.

##  3 Linear models

In this section we study our PCM methodology in the context of a linear model
for $Y$ on $X$ and $Z$. We begin with the simplest version of this setup,
where we assume that $g(x,z):={\mathbb{E}}(Y\,|\,X=x,Z=z)$ is a linear
function that we estimate using ordinary least squares. This is not the sort
of challenging setting where we would envision applying the PCM in practice,
as clearly a $t$-test (modified to account for potential heteroscedasticity)
would suffice to test for the significance of $X$. We nevertheless present it
to show that in contrast to the general methodologies put forward by
Williamson et al., (2022) and Dai et al., (2022), here our method has power
against $n^{-1/2}$-alternatives. In Section 3.2 below, we show that for both
low- and high-dimensional $Z$, we retain Type I error control even under an
arbitrary model for $X$ and when using an essentially arbitrary estimated
projection $\widehat{f}$.

###  3.1 Linear projection function

We consider a family of $\mathcal{P}$ of joint distributions $P$ of
$(X,Y,Z)\in\mathbb{R}\times\mathbb{R}\times\mathbb{R}^{d}$ satisfying the
linear model

$\displaystyle Y=\beta_{P}X+\boldsymbol{\gamma}_{P}^{\top}Z+\zeta_{P},$ (8)

where $\beta_{P}\in\mathbb{R}$ and $\boldsymbol{\gamma}_{P}\in\mathbb{R}^{d}$
are regression coefficients and $\zeta_{P}$ is a random noise term with
${\mathbb{E}}_{P}(\zeta_{P}\,|\,X,Z)=0$. We further impose the following
moment conditions on $\mathcal{P}$.

######  Assumption 1.

  1. (a)

There exist $C,\delta>0$ such that

$\sup_{P\in\mathcal{P}}\max\bigl{\\{}{\mathbb{E}}_{P}(\|Z\|_{\infty}^{4+\delta}),{\mathbb{E}}_{P}(|Y|^{4+\delta}),{\mathbb{E}}_{P}(|X|^{4+\delta})\bigr{\\}}\leq
C.$

  2. (b)

${\mathbb{E}}_{P}(ZZ^{\top})\in\mathbb{R}^{d\times d}$ is invertible, and
writing
$\boldsymbol{\eta}_{P}:={\mathbb{E}}_{P}(ZZ^{\top})^{-1}{\mathbb{E}}_{P}(XZ)$,
$\xi_{P}:=X-\boldsymbol{\eta}_{P}^{\top}Z$,
$\boldsymbol{\Theta}_{P}:={\mathbb{E}}_{P}(ZZ^{\top}\xi_{P}^{2})$, and
$\boldsymbol{\Sigma}_{P}^{XZ}:={\mathbb{E}}_{P}(WW^{\top})$ where
$W:=(X,Z)\in\mathbb{R}\times\mathbb{R}^{d}$, there exists $c>0$ such that

$\inf_{P\in\mathcal{P}}\min\\{{\mathrm{Var}}_{P}(\zeta_{P}\,|\,X,Z),{\mathrm{Var}}_{P}(\xi_{P}),\lambda_{\min}(\boldsymbol{\Sigma}_{P}^{XZ}),\lambda_{\min}(\boldsymbol{\Theta}_{P})\\}\geq
c.$

######  Proposition 1.

Consider a version of the PCM setting $\widehat{v}\equiv 1$ and using OLS for
each of the regressions involved for a family of distributions $\mathcal{P}$
satisfying Assumption 1 and (8). Let
$\mathcal{P}_{1}(\kappa):=\\{P\in\mathcal{P}:|\beta_{P}|\geq\kappa/\sqrt{n}\\}$.
Given any $\alpha\in(0,1)$, we have

$\lim_{\kappa\rightarrow\infty}\lim_{n\to\infty}\inf_{P\in\mathcal{P}_{1}(\kappa)}{\mathbb{P}}_{P}(T>z_{1-\alpha})=1.$

Proposition 1 gives the reassuring conclusion that in the simplest of
settings, our general PCM framework, when used with appropriately chosen
regression methods, can match up to a constant the power properties of a
$t$-test tailored to this setting. In fact it turns out that the context is
simple enough for us to derive an asymptotic power expression for our test. We
present such an analysis in Section D of the appendix for a version of our
test that uses $n_{1}$ and $n_{2}$ (with $n_{1}+n_{2}=2n$) observations in
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ respectively, rather than the equal
split that we consider here. This shows that the optimal splitting ratio
depends on the unknown signal strength, and therefore supports a default
choice of $n_{1}=n_{2}=n$ for simplicity. We also provide a simulation study
in Section E.1 of the appendix where we compare the local power properties of
the PCM, the Williamson et al., (2022) test and the $F$-test with a robust
variance estimator.

###  3.2 A general estimated projection

We next consider a situation where the model is unspecified under the
alternative, whereas $Y$ has a linear relationship with $Z$ under the null of
conditional mean independence. In this case, it is reasonable to employ a
flexible regression method, such as neural networks or random forests, to
estimate the projection $f$. Our goal here is to identify conditions on
estimators, including $\widehat{f}$, under which the proposed test controls
the Type I error. It turns out that, given a specified null model, the problem
of testing whether $X$ is significant is closely connected to goodness-of-fit
testing for the null model, and we are able to exploit this connection to
study the asymptotic Type I error of the proposed test.

Consider first the case of low-dimensional $Z$. Let $\mathcal{P}_{0}$ denote a
family of distributions of $(X,Y,Z)$ under the null where $Z\in\mathbb{R}^{d}$
has an arbitrary distribution and suppose that
$Y=\boldsymbol{\gamma}_{P}^{\top}Z+\varepsilon_{P}$, where
${\mathbb{E}}_{P}(\varepsilon_{P}\,|\,X,Z)=0$. Then
$m_{P}(z)=\boldsymbol{\gamma}_{P}^{\top}Z$ and it is reasonable to use a
linear regression model for $\widehat{m}$. We will suppose that the
regressions yielding $\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ in
Algorithm 1 are performed using OLS, whereas we will leave the regression
choices involved in the construction of $\widehat{f}$ arbitrary. We define
$T_{\mathrm{OLS}}$ as the resulting test statistic, and make the following
assumptions on $\mathcal{P}_{0}$ to ensure uniform asymptotic normality of the
test statistic.

######  Assumption 2.

  1. (a)

There exist $\delta\in(0,2]$, $c,C>0$ such that
${\mathbb{E}}_{P}(\varepsilon_{P}^{2}\,|\,X,Z)\geq c$ and
${{\mathbb{E}}_{P}(|\varepsilon_{P}|^{2+\delta}\,|\,X,Z)\leq C}$ for all
$P\in\mathcal{P}_{0}$.

  2. (b)

For $i\in[n]$, let
$u_{n,i}:=\widehat{f}(X_{i},Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})$ and
$\upsilon_{n,i}:=u_{n,i}/\bigl{(}\sum_{i^{\prime}=1}^{n}u_{n,i^{\prime}}^{2}\bigr{)}^{1/2}$.
Assume that $\max_{i\in[n]}|\upsilon_{n,i}|=o_{\mathcal{P}_{0}}(1)$ and
$\sum_{i=1}^{n}\upsilon_{n,i}^{2}=1+o_{\mathcal{P}_{0}}(1)$.

  3. (c)

Letting $\boldsymbol{\widehat{\gamma}}$ denote the coefficient from the
$\widehat{m}$ regression, assume that
$\max_{i\in[n]}\|Z_{i}\|_{\infty}\cdot\|\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}_{P}\|_{1}=o_{\mathcal{P}_{0}}(1)$.

Part (a) of Assumption 2 concerns conditional moments of $\varepsilon$, and is
used to establish the asymptotic normality of a suitably normalised version of
the numerator of $T_{\mathrm{OLS}}$. In contrast to prior work on goodness-of-
fit testing, e.g. Janková et al., (2020), we do not assume that the
conditional variance of $\varepsilon$ is constant. Assumption 2(b) asks for no
individual $|\upsilon_{n,i}|$ to be significantly larger than the others, and,
for large enough $n$, that at least one of $\\{u_{n,i}:i\in[n]\\}$ is non-zero
for all $P\in\mathcal{P}_{0}$, so $\widehat{f}(X,Z)$ is not constant in $X$.
The latter condition is important for establishing the asymptotic normality of
our test statistic, but is not crucial for Type I error control. Indeed, when
$u_{n,i}=0$ for all $i\in[n]$, the test statistic is zero, and we do not
reject the null. Finally, in settings where, for example each $Z_{i}$ has a
uniformly bounded $(2+\delta)$th moment for some $\delta>0$, we have
$\max_{i\in[n]}\|Z_{i}\|_{\infty}\lesssim n^{1/(2+\delta)}$ and
$\|\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}_{P}\|_{1}\lesssim
n^{-1/2}$ with high probability, and in that case Part (c) of Assumption 2 is
satisfied.

######  Proposition 2 (Low-dimensional $Z$).

Suppose in the above setting that Assumption 2 holds and that all OLS
estimators exist almost surely. Then the test statistic $T_{\mathrm{OLS}}$
converges to $N(0,1)$ uniformly over $\mathcal{P}_{0}$; i.e.,

$\displaystyle\sup_{P\in\mathcal{P}_{0}}\sup_{t\in\mathbb{R}}\big{|}{\mathbb{P}}_{P}(T_{\mathrm{OLS}}\leq
t)-\Phi(t)\big{|}\rightarrow 0.$

Under the conditions of Proposition 2, the test that rejects the null when
$T_{\mathrm{OLS}}>z_{1-\alpha}$ is uniformly asymptotically of size $\alpha$.
We also note that the only requirement imposed on the projection $\widehat{f}$
is that it satisfies Assumption 2(b). Janková et al., (2020) also consider
this condition, providing supporting empirical evidence in general, and
introducing a specific procedure to guarantee that the condition holds.

We now extend these ideas and the setting described above Assumption 2 to the
case where the dimension of $Z$ is potentially larger than the sample size.
Here, the least squares estimator is not necessarily well-defined, so we
construct $\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ using the Lasso or one
of its variants. Letting $\boldsymbol{\widehat{\gamma}}$ denote the
coefficients from the $\widehat{m}$ regression, the motivation for this comes
from the decomposition

$\displaystyle\sum_{i=1}^{n}(Y_{i}-\boldsymbol{\widehat{\gamma}}^{\top}Z_{i})u_{n,i}=\sum_{i=1}^{n}\varepsilon_{P,i}u_{n,i}-\delta_{\mathrm{bias}},$
(9)

where
$\delta_{\mathrm{bias}}:=\sum_{i=1}^{n}(\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}_{P})^{\top}Z_{i}u_{n,i}$.
This bias term is no longer exactly zero as for the least squares estimators
considered in Proposition 2, but Hölder’s inequality will nevertheless
guarantee that it is sufficiently small for our purposes as long as

$\displaystyle\|\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}_{P}\|_{1}\cdot\bigg{\|}\sum_{i=1}^{n}Z_{i}\upsilon_{n,i}\bigg{\|}_{\infty}=o_{\mathcal{P}_{0}}(1).$
(10)

We denote the test statistic as described in Algorithm 1 in this context by
$T_{\mathrm{Lasso}}$. The next proposition is the analogue of Proposition 2
for $T_{\mathrm{Lasso}}$.

######  Proposition 3 (High-dimensional $Z$).

Suppose in the above setting that Assumption 2 and condition (10) hold. Then

$\displaystyle\sup_{P\in\mathcal{P}_{0}}\sup_{t\in\mathbb{R}}\big{|}{\mathbb{P}}_{P}(T_{\mathrm{Lasso}}\leq
t)-\Phi(t)\big{|}\rightarrow 0.$

In order to ensure that condition (10) holds, one can use the square-root
Lasso (Belloni et al.,, 2011; Sun and Zhang,, 2012), as suggested by Janková
et al., (2020). In particular, for $\lambda_{\mathrm{sq}}>0$, we set
$\widehat{m}_{\widehat{f}}(z)=\boldsymbol{\widehat{\eta}}_{\mathrm{sq}}^{\top}z$
where

$\displaystyle\boldsymbol{\widehat{\eta}}_{\mathrm{sq}}:=\operatorname*{argmin}_{\boldsymbol{\eta}\in\mathbb{R}^{d}}\bigg{\\{}\biggl{(}\frac{1}{n}\sum_{i=1}^{n}\bigl{(}\widehat{f}(X_{i},Z_{i})-\boldsymbol{\eta}^{\top}Z_{i}\bigr{)}^{2}\biggr{)}^{1/2}+\lambda_{\mathrm{sq}}\|\boldsymbol{\eta}\|_{1}\bigg{\\}}.$

With this choice of $\boldsymbol{\widehat{\eta}}_{\mathrm{sq}}$ and by letting
$\lambda_{\mathrm{sq}}=c_{\mathrm{sq}}\sqrt{(\log d_{Z})/n}$ for some constant
$c_{\mathrm{sq}}>0$, the Karush–Kuhn–Tucker conditions for the square-root
Lasso guarantee that $\|\sum_{i=1}^{n}Z_{i}\upsilon_{n,i}\|_{\infty}\leq
c_{\mathrm{sq}}\sqrt{\log d_{Z}}$. Furthermore, under appropriate conditions,
the Lasso estimator $\boldsymbol{\widehat{\gamma}}$ has an error bound
$\|\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}_{P}\|_{1}\lesssim
s\sqrt{(\log d_{Z})/n}$ with high probability, where $s$ denotes the number of
non-zero coefficients of $\boldsymbol{\gamma}_{P}$ (e.g. Corollary 6.2 of
Bühlmann and van de Geer,, 2011). Therefore, in this setting, condition (10)
is satisfied provided that $s(\log d_{Z})/\sqrt{n}\rightarrow 0$.

##  4 General theory

In this section, we present general conditions ensuring uniform asymptotic
validity and power of the test, primarily by imposing assumptions on the
performance of the regressions involved. To facilitate our analysis, it is
helpful to study a slight modification of the test as presented in Algorithm
1, where we form $\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ on an
independent auxiliary sample. In principle, we may accomplish this by further
splitting $\mathcal{D}_{1}$ into two, and using one part to train
$\widehat{m}$ and $\widehat{m}_{\widehat{f}}$, and the other to compute the
test statistic. Moreover, we can exchange the roles of the two parts and
average the resulting test statistics, a process known as cross-fitting
(Chernozhukov et al.,, 2018), which then guarantees no loss in efficiency from
this additional sample split. However, for the reasons discussed in Section
2.1 we do not recommend performing this further sample split in practice.

The following quantities relating to the performances of the regression
methods $\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ will play a key role in
our results. Let us introduce

$\varepsilon_{P,i}:=Y_{i}-m_{P}(Z_{i}),\qquad\xi_{P,i}:=\widehat{f}(X_{i},Z_{i})-m_{P,\widehat{f}}(Z_{i}),$
(11)

for $i\in[n]$, and an analogous version of (11) without a subscript $i$.
Further, define

$\sigma_{P}^{2}:={\mathrm{Var}}_{P}\bigl{(}\xi_{P}\,|\,\widehat{f}\bigr{)},$
(12)

as well as

$\mathcal{E}_{P,1}:=\frac{1}{n}\sum_{i=1}^{n}\\{m_{P}(Z_{i})-\widehat{m}(Z_{i})\\}^{2}\quad\text{and}\quad\mathcal{E}_{P,2}:=\frac{1}{n\sigma_{P}^{2}}\sum_{i=1}^{n}\\{m_{P,\widehat{f}}(Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})\\}^{2}.$
(13)

The second MSPE $\mathcal{E}_{P,2}$ in the display above is normalised by the
variance of the errors $\xi_{P,i}$ featuring in the corresponding regression.
Under the null, we expect this variance to be small as $\widehat{f}$ is then
estimating a zero function, and consequently $\mathcal{E}_{P,2}$ may be
inflated relative to the unnormalised version of this quantity. On the other
hand, as $\widehat{f}$ is small, we can expect that the unnormalised MSPE is
particularly small. For example, writing $\mathcal{P}$ for the simple null
linear model considered in (7), we would have

$\frac{1}{n}\sum_{i=1}^{n}\\{m_{P,\widehat{f}}(Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})\\}^{2}=O_{\mathcal{P}}(n^{-2})\qquad\text{and}\qquad
1/\sigma_{P}^{2}=O_{\mathcal{P}}(n),$

giving $\mathcal{E}_{P,2}=O_{\mathcal{P}}(n^{-1})$.

###  4.1 Type I error control

We consider the following assumption regarding general Type I error control.

######  Assumption 3.

Consider a class of null distributions $\mathcal{P}_{0}$ of $(X,Y,Z)$ on
$\mathcal{X}\times\mathbb{R}\times\mathcal{Z}$ with
${\mathbb{E}}_{P}(Y\,|\,X,Z)={\mathbb{E}}_{P}(Y\,|\,Z)$ for which there exists
$c>0$ such that
$\inf_{P\in\mathcal{P}_{0}}{\mathbb{E}}_{P}(\varepsilon_{P}^{2}\,|\,X,Z)\geq
c,$ and the following hold:

  1. (a)

$\sup_{P\in\mathcal{P}_{0}}{\mathbb{P}}_{P}(\sigma_{P}^{2}=0)=o(1)$.

  2. (b)

The product of the MSPEs satisfies
$\mathcal{E}_{P,1}\mathcal{E}_{P,2}=o_{\mathcal{P}_{0}}(n^{-1})$.

  3. (c)

The weighted MSPEs scaled by $\sigma_{P}^{2}$ satisfy

$\displaystyle\frac{1}{n\sigma_{P}^{2}}\sum_{i=1}^{n}\\{m_{P}(Z_{i})-\widehat{m}(Z_{i})\\}^{2}\xi_{P,i}^{2}$
$\displaystyle=o_{\mathcal{P}_{0}}(1),$
$\displaystyle\frac{1}{n\sigma_{P}^{2}}\sum_{i=1}^{n}\\{m_{P,\widehat{f}}(Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})\\}^{2}\varepsilon_{P,i}^{2}$
$\displaystyle=o_{\mathcal{P}_{0}}(1).$

  4. (d)

There exists $\delta\in(0,2]$ such that
${\mathbb{E}}_{P}\bigl{(}|\varepsilon_{P}\xi_{P}|^{2+\delta}\,|\,\widehat{f}\bigr{)}/\sigma_{P}^{2+\delta}=o_{\mathcal{P}_{0}}(n^{\delta/2})$.

Part (a) of Assumption 3 asks that $\widehat{f}(X,Z)$ is not exactly constant
in $X$ (even though we may expect it not to vary too much with $X$, for the
reasons explained in the discussion in Section 2.1). Part (b) should be
regarded as the primary restriction on $\mathcal{P}_{0}$, and along with (c),
relates directly to the performance of the user-chosen regression methods
involved in the construction of the PCM. As alluded to above, in a simple
linear model setting, we can expect
$\mathcal{E}_{P,1}\mathcal{E}_{P,2}=O_{\mathcal{P}_{0}}(n^{-2})$, which
certainly satisfies the condition. The rate requirement on the product of
MSPEs is however sufficiently slow to accommodate nonparametric models; see
Section 5. We note that the deterministic condition
$\sup_{P\in\mathcal{P}_{0}}\\{{\mathbb{E}}_{P}(\mathcal{E}_{P,1}){\mathbb{E}}_{P}(\mathcal{E}_{P,2})\\}=o(n^{-1})$
is sufficient to guarantee part (b), as can be verified via Markov’s
inequality and the Cauchy–Schwarz inequality. If in addition there exists
$C>0$ such that ${\mathrm{Var}}_{P}(\xi_{P}\,|\,Z,\widehat{f})\leq
C\sigma_{P}^{2}$ and ${\mathrm{Var}}_{P}(\varepsilon_{P})\leq C$, then (c) is
guaranteed when $\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ satisfy the
simple consistency property
$\max\bigl{\\{}{\mathbb{E}}_{P}(\mathcal{E}_{P,1}),\,{\mathbb{E}}(\mathcal{E}_{P,2})\bigr{\\}}=o(1)$.
Part (d) is a conditional Lyapunov condition, and is used to apply the central
limit theorem for triangular arrays.

######  Theorem 4 (Asymptotic normality under the null of a general
procedure).

Suppose that Assumption 3 holds over a class of null distributions
$\mathcal{P}_{0}$. Then

$\displaystyle\sup_{P\in\mathcal{P}_{0}}\sup_{t\in\mathbb{R}}\bigl{|}{\mathbb{P}}_{P}(T\leq
t)-\Phi(t)\bigr{|}\rightarrow 0.$

The proof of Theorem 4 can be found in Appendix A.4, which formalises the
brief explanation of asymptotic normality laid down in Section 2.1. The above
result indicates that the asymptotic normality of $T$ (hence the validity of
the PCM test) is largely determined by the predictive performance of
regression models used in construction of the test statistic. Although we have
stated Theorem 4 for the variant of our PCM procedure where $\widehat{m}$ and
$\widehat{m}_{\widehat{f}}$ are formed on an auxiliary sample, it turns out
that the conclusion holds more generally. Indeed, it follows from the proof
that under Assumption 3, the abstract condition (22) on our procedure
suffices. In addition to being trivially satisfied when $\widehat{m}$ is
formed on an auxiliary sample, (22) holds for the practical version of our
test where $\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ are formed on
$\mathcal{D}_{1}$ provided either that
$X\mbox{${}\perp\mkern-11.0mu\perp{}$}Y\,|\,Z$ or that $\widehat{m}$ is a
linear smoother. See Proposition 20 in Appendix B.2 for more details.

###  4.2 Power properties

When studying the power properties of our test, we restrict attention to a
subset of alternatives that are separated from null distributions
characterised by Assumption 4 below.

######  Assumption 4.

Given a positive sequence $(\epsilon_{n})_{n\in\mathbb{N}}$, let
$\bigl{(}\mathcal{P}_{1}(\epsilon_{n})\bigr{)}_{n\in\mathbb{N}}$ denote a
sequence of collections of alternative distributions such that

$\inf_{P\in\mathcal{P}_{1}(\epsilon_{n})}\tau_{P}\geq\epsilon_{n}.$

Further, suppose that there exists $C>0$ with

$\limsup_{n\to\infty}\sup_{P\in\mathcal{P}_{1}(\epsilon_{n})}\max\\{{\mathrm{Var}}_{P}(Y\,|\,X,Z),h_{P}(X,Z)\\}\leq
C,$

and that the following conditions are satisfied:

  1. (a)

There exists $\beta_{1}>0$ such that
$\mathcal{E}_{P,1}=O_{\mathcal{P}_{1}(\epsilon_{n})}(n^{-\beta_{1}})$.

  2. (b)

There exists $\beta_{2}>0$ such that
$\mathcal{E}_{P,2}=O_{\mathcal{P}_{1}(\epsilon_{n})}(n^{-\beta_{2}})$.

  3. (c)

There exists $\rho>0$ such that

$\sup_{P\in\mathcal{P}_{1}(\epsilon_{n})}{\mathbb{P}}_{P}\Bigl{(}{\mathrm{Corr}}_{P}\bigl{(}h_{P}(X,Z),\xi_{P}\,|\,\widehat{f}\bigr{)}\leq\rho\Bigr{)}=o(1).$

In addition to the rate requirements on MSPEs in (a) and (b), condition (c)
requires $\xi_{P}$, the population residuals from regressing our estimated
projection $\widehat{f}$ onto $Z$, to be positively correlated with $h_{P}$
with high probability (there is no need for the correlation to approach $1$).
To interpret (c), it is helpful to consider a stronger version of this
condition with $\widehat{f}(X,Z)$ replacing $\xi_{P}$; to see that this
results in a stronger condition, note that
${\mathbb{E}}_{P}\bigl{(}h_{P}(X,Z)\xi_{P}\,|\,\widehat{f}\bigr{)}={\mathbb{E}}_{P}\bigl{(}h_{P}(X,Z)\widehat{f}(X,Z)\,|\,\widehat{f}\bigr{)}$
and
${\mathbb{E}}(\xi_{P}^{2}\,|\,\widehat{f})={\mathbb{E}}_{P}\bigl{[}{\mathrm{Var}}_{P}\bigl{(}\widehat{f}(X,Z)\,|\,Z,\widehat{f}\bigr{)}\,|\,\widehat{f}\bigr{]}\leq{\mathbb{E}}\bigl{(}\widehat{f}(X,Z)^{2}\,|\,\widehat{f}\bigr{)}$.
This stronger assumption still permits $\widehat{f}$ to be an inconsistent
estimator of the true $f_{P}(X,Z)=h_{P}(X,Z)/{\mathrm{Var}}_{P}(Y\,|\,X,Z)$ in
that it only requires them to be positively correlated, with probability
approaching one.

The flexibility afforded by this assumption relies on using regression method
$\widehat{m}_{\widehat{f}}$ being scale equivariant in the sense that

$\widehat{m}_{a\cdot\widehat{f}}(Z)=a\cdot\widehat{m}_{\widehat{f}}(Z)$ (14)

for all $a>0$; this is a mild condition, however, that is satisfied by many
regression methods.

We can now state the main result of this subsection.

######  Theorem 5.

Assume that $\widehat{m}_{\widehat{f}}$ is scale equivariant in the sense that
it satisfies (14) and consider the sequence of classes of distributions in
Assumption 4, where

$\epsilon_{n}\cdot n^{\min\\{1,\beta_{1}+\beta_{2}\\}}\rightarrow\infty.$ (15)

Then for any $\alpha\in(0,1)$,

$\inf_{P\in\mathcal{P}_{1}(\epsilon_{n})}{\mathbb{P}}_{P}(T>z_{1-\alpha})\rightarrow
1.$

Theorem 5 shows that if $\beta_{1}+\beta_{2}\geq 1$, as may be expected if $m$
and $m_{\widehat{f}}$ are sufficiently smooth, then the rate $\epsilon_{n}$ is
going to be driven by Assumption 4(c). It is possible to derive a version of
Theorem 5 for the test as described in Algorithm 1 that does not employ the
additional sample splitting we are considering here. The only change is that
(15) becomes $\epsilon_{n}\cdot
n^{\min\\{1,\beta_{1},\beta_{2}\\}}\rightarrow\infty$; however we believe the
version of Theorem 5 above is more in line with the behaviour to be expected
in practice, and empirically we find the version of the test in Algorithm 1 to
provide better discrimination between null and alternatives in finite samples.

##  5 Series estimators

Following the theory in the previous section for a general regression method,
we will now provide more concrete results for a specific class, namely spline
estimators. In particular, our interest is to identify conditions under which
our test is uniformly asymptotically valid and attains near-optimal power in a
nonparametric setting. A formal power analysis, however, is complicated by the
fact that $\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ are computed on the
same subsample as our test statistic. We therefore consider a slightly
modified test statistic that leverages ideas from the literature on cross-
fitting. Throughout this section we assume that
$(X,Z)\in[0,1]^{d_{X}}\times[0,1]^{d_{Z}}$ and set $d:=d_{X}+d_{Z}$.

Due to our additional sample splitting, we will require two additional
independent samples of size $n$, so that we have
$\mathcal{D}_{1},\ldots,\mathcal{D}_{4}$ in total. In Appendix C, we give a
self-contained description of spline spaces and their tensor product B-spline
bases, containing all the results that we require for our analysis. Given a
spline order $r\in\mathbb{N}$ (i.e. degree $r-1$) and $N\in\mathbb{N}_{0}$
equi-spaced interior knots in each dimension, we denote by
$\mathcal{S}_{r,N}^{d_{Z}}$ the corresponding spline space on $[0,1]^{d_{Z}}$,
and by $\boldsymbol{\phi}^{Z}$ its $d_{Z}$-tensor B-spline basis, which
consists of $K_{Z}:=(N+r)^{d_{Z}}$ basis functions. Writing
$\mathcal{S}_{r,N}^{d_{X}}$ for the corresponding spline space on
$[0,1]^{d_{X}}$ with $d_{X}$-tensor B-spline basis $\boldsymbol{\phi}^{X}$,
having $K_{X}:=(N+r)^{d_{X}}$ basis functions, we can define the $d$-tensor
product basis
$\boldsymbol{\phi}(x,z):=\boldsymbol{\phi}^{X}(x)\otimes\boldsymbol{\phi}^{Z}(z)$
for $\mathcal{S}_{r,N}^{d}$, where
$\boldsymbol{u}\otimes\boldsymbol{v}:=\mathrm{vec}(\boldsymbol{u}\boldsymbol{v}^{\top})$,
having $K_{XZ}:=K_{X}K_{Z}$ basis functions. Further, we let
$\boldsymbol{\psi}$ denote the tensor product B-spline basis for
$\mathcal{S}_{2r-1,N}^{d_{Z}}$, and write
$\widetilde{K}_{Z}:=(N+2r-1)^{d_{Z}}$; the higher order of the spline basis
functions that make up $\boldsymbol{\psi}$ affords better approximation
properties that turn out to be useful for our theory.

The following description of the test statistic fixes notation and follows
Algorithm 1 except that we fit $\widehat{m}_{\widehat{f}}$ on
$\mathcal{D}_{3}$, $\widehat{m}$ on $\mathcal{D}_{4}$ and set
$\widehat{v}\equiv 1$ for simplicity. We will also omit discussion of the sign
correction step (Algorithm 1 1(iv)), since $\widehat{\rho}$ is always non-
negative for the estimators considered below. We first regress $Y$ onto
$\boldsymbol{\phi}(X,Z)$ using ordinary least squares (OLS) on
$\mathcal{D}_{2}$, yielding an estimator
$\boldsymbol{\widehat{\beta}}_{XZ}\in\mathbb{R}^{K_{XZ}}$, and set
$\widehat{g}(x,z):=\boldsymbol{\widehat{\beta}}_{XZ}^{\top}\boldsymbol{\phi}(x,z)$.
We then regress $\widehat{g}(X,Z)$ onto $\boldsymbol{\phi}(X,Z)$ using OLS on
$\mathcal{D}_{2}$ again, to obtain an estimator
$\boldsymbol{\widehat{\beta}}_{Z}\in\mathbb{R}^{K_{Z}}$, and set
$\widetilde{m}(z):=\boldsymbol{\widehat{\beta}}_{Z}^{\top}\boldsymbol{\phi}^{Z}(z)$.
Note that this is equivalent to regressing $Y$ onto
$\boldsymbol{\phi}^{Z}(Z)$. We then define the projection
$\widehat{f}(x,z):=\widehat{g}(x,z)-\widetilde{m}(z)$. Using the fact that
$\boldsymbol{\phi}^{X}$ forms a partition of unity (Proposition 28(a) in the
appendix), it follows that if we write
$\widehat{\boldsymbol{\beta}}:=\boldsymbol{\widehat{\beta}}_{XZ}-\boldsymbol{1}\otimes\boldsymbol{\widehat{\beta}}_{Z}$,
where $\boldsymbol{1}\in\mathbb{R}^{K_{X}}$ denotes a vector of ones, then
$\widehat{f}(x,z)=\widehat{\boldsymbol{\beta}}^{\top}\boldsymbol{\phi}(x,z)$.

To estimate $m_{P,\widehat{f}}$, we regress $\widehat{f}(X,Z)$ onto
$\boldsymbol{\psi}(Z)$ on $\mathcal{D}_{3}$ using OLS, giving
$\widehat{m}_{\widehat{f}}$. Similarly, we estimate $m_{P}$ by regressing $Y$
onto $\boldsymbol{\psi}(Z)$ on $\mathcal{D}_{4}$ using OLS to give
$\widehat{m}$. Given $\widehat{f}$, $\widehat{m}$ and
$\widehat{m}_{\widehat{f}}$ as defined above, we compute the test statistic as
in Algorithm 1 (on $\mathcal{D}_{1}$) and denote it by $T_{\mathrm{Spline}}$.
In Theorem 6 in Section 5.1 below, we demonstrate that $T_{\mathrm{Spline}}$
enjoys uniform asymptotic Type I error control under appropriate regularity
conditions, while Theorem 7 and Proposition 8 in Section 5.2 reveal that
$T_{\mathrm{Spline}}$ can achieve the optimal testing rate for this problem.

###  5.1 Type I error control

We start by stating our main distributional assumptions, which rely on the
definitions of Hölder spaces $\mathcal{H}_{s}^{d}$ and Hölder norms
$\|\cdot\|_{\mathcal{H}_{s}}$ given in Definition 29.

######  Assumption 5.

Let $\mathcal{P}$ be a class of distributions of $(X,Y,Z)$ on
$[0,1]^{d_{X}}\times\mathbb{R}\times[0,1]^{d_{Z}}$, and for $P\in\mathcal{P}$,
let $g_{P}(x,z):={\mathbb{E}}_{P}(Y\,|\,X=x,Z=z)$. Assume that there exist
$C\geq 1$ and $c\in(0,1]$ with the following properties:

  1. (a)

For each $P\in\mathcal{P}$, we have
${\mathbb{E}}_{P}(\varepsilon_{P}^{2}\,|\,X,Z)\geq c$ and there exists
$\delta\in(0,2]$ such that
${\mathbb{E}}_{P}(|\varepsilon_{P}|^{2+\delta}\,|\,X,Z)\leq C$.

  2. (b)

For each $P\in\mathcal{P}$, the marginal distribution of $(X,Z)$ is absolutely
continuous with respect to Lebesgue measure on $[0,1]^{d}$, with density
$p_{P}$ satisfying $\sup_{(x,z)\in[0,1]^{d}}p_{P}(x,z)\leq C$ and
$\inf_{(x,z)\in[0,1]^{d}}p_{P}(x,z)\geq c$.

  3. (c)

Let $s\in(0,r]$ and let $p_{X|Z,P}(\cdot\,|\,z)$ denote the conditional
density of $X$ given $Z=z$. Assume that for all $P\in\mathcal{P}$, we have
$p_{X|Z,P}(x\,|\,\cdot)\in\mathcal{H}^{d_{Z}}_{s}$ for all
$x\in[0,1]^{d_{X}}$, and that $m_{P}\in\mathcal{H}^{d_{Z}}_{s}$ and
$g_{P}\in\mathcal{H}^{d}_{s}$, with

$\max\biggl{\\{}\sup_{x\in[0,1]^{d_{X}}}\|p_{X|Z,P}(x,\cdot)\|_{\mathcal{H}_{s}},\|m_{P}\|_{\mathcal{H}_{s}},\|g_{P}\|_{\mathcal{H}_{s}}\biggr{\\}}\leq
C.$

Assumption 5 is closely related to other assumptions commonly used in spline
regression (e.g. Belloni et al.,, 2015; Ichimura and Newey,, 2015; Newey and
Robins,, 2018). In order to state our Type I error control result for
$T_{\mathrm{Spline}}$, it will be convenient to define the projection
$\boldsymbol{\Pi}:\mathbb{R}^{K_{XZ}}\rightarrow\mathbb{R}^{K_{XZ}}$ by
$\boldsymbol{\Pi}(\boldsymbol{x})\equiv\boldsymbol{\Pi}(x_{1},\ldots,x_{K_{XZ}}):=\boldsymbol{x}-\boldsymbol{1}\otimes\bar{\boldsymbol{x}}$,
with $\bar{\boldsymbol{x}}=(\bar{x}_{1},\ldots,\bar{x}_{K_{Z}})$ given by
$\bar{x}_{k}:=K_{X}^{-1}\sum_{\ell=1}^{K_{X}}x_{(k-1)K_{X}+\ell}$ for
$k\in[K_{Z}]$.

######  Theorem 6 (Asymptotic normality of $T_{\mathrm{Spline}}$).

Suppose that Assumption 5 holds for a class of null distributions
$\mathcal{P}_{0}$, i.e. a class of distributions that also satisfies
${\mathbb{E}}_{P}(Y\,|\,X,Z)={\mathbb{E}}_{P}(Y\,|\,Z)$ for every
$P\in\mathcal{P}_{0}$. Assume that
$\sup_{P\in\mathcal{P}_{0}}{\mathbb{P}}_{P}(\|\boldsymbol{\Pi}\boldsymbol{\widehat{\beta}}\|_{\infty}=0)=o(1)$
and that
$\boldsymbol{\Lambda}_{P}:={\mathbb{E}}_{P}\bigl{\\{}{\mathrm{Cov}}_{P}\bigl{(}\boldsymbol{\phi}(X,Z)\,|\,Z\bigr{)}\bigr{\\}}$
satisfies

$\tilde{\lambda}_{\min}(\boldsymbol{\Lambda}_{P}):=\min_{\boldsymbol{x}\in\mathbb{R}^{K_{XZ}}:\boldsymbol{\Pi}\boldsymbol{x}=\boldsymbol{x},\|\boldsymbol{x}\|_{2}=1}\boldsymbol{x}^{\top}\boldsymbol{\Lambda}_{P}\boldsymbol{x}\geq\frac{c}{K_{XZ}},$
(16)

for each $P\in\mathcal{P}_{0}$, where $c\in(0,1]$ is taken from Assumption 5.
Finally, suppose that

$nK_{XZ}\biggl{\\{}\widetilde{K}_{Z}^{-2s/d_{Z}}+\frac{\widetilde{K}_{Z}}{n}\biggr{\\}}^{2}\rightarrow
0$ (17)

and

$\frac{K_{XZ}^{1+2/\delta}}{n}\rightarrow 0$ (18)

where $\delta$ is taken from Assumption 5. Then

$\displaystyle\sup_{P\in\mathcal{P}_{0}}\sup_{t\in\mathbb{R}}\bigl{|}{\mathbb{P}}_{P}(T_{\mathrm{Spline}}\leq
t)-\Phi(t)\bigr{|}\rightarrow 0.$

The proof of Theorem 6 amounts to the verification of Assumption 3, which then
allows us to apply our general Type I error control result, namely Theorem 4.
In addition to Assumption 5, Theorem 6 imposes several additional conditions.
The assumption that
$\sup_{P\in\mathcal{P}_{0}}{\mathbb{P}}_{P}(\|\boldsymbol{\Pi}\boldsymbol{\widehat{\beta}}\|_{\infty}=0)=o(1)$
simply avoids degeneracy of the test statistic and is used to show that
Assumption 3(a) is satisfied. If this condition is not satisfied, then since
we defined $0/0:=0$ in the definition of our test statistic, it can be shown
that

$\liminf_{n\rightarrow\infty}\inf_{P\in\mathcal{P}_{0}}{\mathbb{P}}_{P}(|T_{\mathrm{Spline}}|\leq
t)\geq\Phi(t)-\Phi(-t)$

for all $t\geq 0$ (i.e. $|T_{\mathrm{Spline}}|$ is asymptotically
stochastically dominated by the absolute value of a standard Gaussian random
variable), so the test retains uniform asymptotic Type I error control
provided that $\alpha\leq 1/2$.

Condition (16) can be regarded as a restricted minimum eigenvalue condition;
for $\boldsymbol{x}\in\mathbb{R}^{K_{XZ}}$ with
$\boldsymbol{\Pi}\boldsymbol{x}=0$, we have that
$\boldsymbol{x}^{\top}\boldsymbol{\Lambda}_{P}\boldsymbol{x}=0$, but it turns
out that we are able to restrict attention to the orthogonal complement of
this subspace. Motivation for the form of this condition is provided by the
fact that, writing
$\boldsymbol{\Sigma}_{P}:={\mathbb{E}}_{P}\bigl{(}\boldsymbol{\phi}(X,Z)\boldsymbol{\phi}(X,Z)^{\top}\bigr{)}\in\mathbb{R}^{K_{XZ}\times
K_{XZ}}$, we have

$\tilde{\lambda}_{\min}(\boldsymbol{\Lambda}_{P})\leq\tilde{\lambda}_{\min}(\boldsymbol{\Sigma}_{P})\leq\lambda_{\max}(\boldsymbol{\Sigma}_{P})\leq
C2^{d}K_{XZ}^{-1}$

by Proposition 28(d). Moreover, Lemma 22 in Section B.2 of the appendix shows
that the assumption holds when $X$ and $Z$ are independent.

Condition (17) is used to show that parts (b) and (c) of Assumption 3 are
satisfied while Condition (18) is used to show that part (d) of the assumption
is satisfied. These conditions control the interplay between the growth rate
of the number of basis functions, the smoothness $s$ of the regression
functions and conditional densities and $\delta$. When choosing the knot
spacing to minimise the mean-squared error of the involved regressions, we
would choose $\widetilde{K}_{Z}$ and $K_{Z}$ of order $n^{d_{Z}/(2s+d_{Z})}$
and $K_{X}$ of order $n^{d_{X}/(2s+d_{Z})}$. Thus for (17) to hold, we need
$s>d_{Z}+d_{X}/2$ and for (18) to hold, we need
$\delta>2(d_{X}+d_{Z})/(2s-d_{X})$. Both of conditions could be weakened, at
the expense of additional notational complexity, by choosing different knot
spacings $N_{X}$ and $N_{Z}$ for the $d_{X}$\- and $d_{Z}$-tensor B-spline
bases $\boldsymbol{\phi}^{X}$ and $\boldsymbol{\phi}^{Z}$ for our spline
spaces $\mathcal{S}_{r,N_{X}}^{d_{X}}$ and $\mathcal{S}_{r,N_{Z}}^{d_{Z}}$.
Indeed, by taking $N_{X}$, and hence $K_{X}$, to be of constant order, while
retaining the original choices of $K_{Z}$ and $\widetilde{K}_{Z}$, we see that
(17) holds when $s>d_{Z}$ and (18) holds when $\delta>d_{Z}/s$ (so it would
suffice for Assumption 5(a) to hold with $\delta=1$, provided again that
$s>d_{Z}$).

###  5.2 Power and minimax lower bound

As mentioned at the beginning of this section, we employ additional sample
splitting in the construction of $T_{\mathrm{Spline}}$. This turns out to be
helpful in demonstrating the optimality of our test. To provide insight into
the benefits of sample splitting in this context, consider two generic spline
estimators $\widehat{g}_{1}$ and $\widehat{g}_{2}$ of unknown functions
$g_{1}$ and $g_{2}$, respectively. Suppose that we would like to choose
$\widehat{g}_{1}$ and $\widehat{g}_{2}$ to minimise the empirical cross-
product error

$\displaystyle\widehat{\theta}_{\mathrm{cross}}:=\frac{1}{n}\sum_{i=1}^{n}\\{\widehat{g}_{1}(Z_{i})-g_{1}(Z_{i})\\}\\{\widehat{g}_{2}(Z_{i})-g_{2}(Z_{i})\\}.$
(19)

A naive way of approaching this problem is to construct $\widehat{g}_{1}$ and
$\widehat{g}_{2}$ on the same dataset and to choose the number of spline
functions so as to minimise the mean-squared error of each of
$\widehat{g}_{1}$ and $\widehat{g}_{2}$. The Cauchy–Schwarz inequality then
guarantees that the cross-product error is small as long as the mean-squared
errors are small. However, this indirect approach returns a potentially
suboptimal rate of convergence due to its “own observation” bias, which arises
from using the same datasets to form $\widehat{g}_{1}$ and $\widehat{g}_{2}$.
When employing auxiliary samples to construct $\widehat{g}_{1}$ and
$\widehat{g}_{2}$ we can eliminate this bias; thus a more refined analysis of
terms like $\widehat{\theta}_{\mathrm{cross}}$ that does not directly employ
the Cauchy–Schwarz inequality can result in faster convergence rates; see for
instance Proposition 40 in the appendix. Our main result in this section is as
follows:

######  Theorem 7.

Let $\mathcal{P}$ be a class of distributions satisfying Assumption 5, and let
$\mathcal{P}_{1}(\epsilon_{n}):=\\{P\in\mathcal{P}:\tau_{P}\geq\epsilon_{n}\\}$,
where

$\displaystyle\epsilon_{n}\cdot n^{\frac{4s}{4s+d}}\rightarrow\infty.$ (20)

Further, assume that the tuning parameters are chosen such that $K_{X}\asymp
n^{\frac{2d_{X}}{4s+d}}$ and $K_{Z}\asymp\widetilde{K}_{Z}\asymp
n^{\frac{2d_{Z}}{4s+d}}$ and that $r\geq s\geq 3d/4$. Then

$\displaystyle\inf_{P\in\mathcal{P}_{1}(\epsilon_{n})}{\mathbb{P}}_{P}(T_{\mathrm{Spline}}>z_{1-\alpha})\rightarrow
1.$

Theorem 7 reveals that the test based on $T_{\mathrm{Spline}}$ has uniform
asymptotic power 1 over a class of alternatives that are sufficiently
separated from the null, as defined by $\mathcal{P}_{1}(\epsilon_{n})$.

We remark that in Theorem 7, we have operated in the context of a known
smoothness parameter $s$ for theoretical purposes. It is possible to construct
more involved tests that adapt to unknown smoothness levels following the
strategy of Lepskiǐ, (1991) and Ingster, (2000), but we do not pursue this
direction further.

The separation rate (20) cannot be improved further from a minimax
perspective, as illustrated by the following lower bound result.

######  Proposition 8.

Consider a class of distributions, denoted by $\mathcal{P}$, that satisfy
Assumption 5, and let
$\mathcal{P}_{1}(\epsilon_{n}):=\\{P\in\mathcal{P}:\tau_{P}\geq\epsilon_{n}\\}$.
Then, for a fixed level $\alpha\in(0,1/2)$, there exists $c>0$ such that if
$\limsup_{n\rightarrow\infty}\epsilon_{n}\cdot n^{\frac{4s}{4s+d}}<c$, then
any test $\phi$ having uniform asymptotic size at most $\alpha$ satisfies

$\displaystyle\limsup_{n\rightarrow\infty}\inf_{P\in\mathcal{P}_{1}(\epsilon_{n})}{\mathbb{P}}_{P}(\phi=1)\leq\alpha+1/2.$

Proposition 8 complements Theorem 7 by showing that when $\tau_{P}$ is a small
constant multiple of $n^{-\frac{4s}{4s+d}}$, no test can achieve uniform
consistency under Hölder smoothness. The proof of Proposition 8, which can be
found in Appendix A.8, follows a fairly standard argument (e.g. Ingster,,
1987; Arias-Castro et al.,, 2018) that bounds the $\chi^{2}$-divergence from a
fixed null distribution to a mixture of distributions in the alternative class
$\mathcal{P}_{1}(\epsilon_{n})$.

##  6 Numerical experiments

In this section, we present the results of several simulation experiments that
investigate the empirical performances of both the recommended multiple sample
splitting version of the PCM (see Algorithm 2) with $B=6$ splits, denoted by
pcm, and the single split version (see Algorithm 1) denoted by pcm_ss. We
compare our tests to various conditional (mean) independence tests in the
literature listed below.

  * gam

The test based on the default $p$-value for a smooth when fitting a
generalised additive model (GAM) using the mgcv-package in R (Wood,, 2013,
2017).

  * williamson

The test resulting from applying the approach described in Williamson et al.,
(2021) and employing sample splitting as implemented in the vim function from
the vimp-package in R (Williamson et al.,, 2022).

  * kci

The _kernel conditional independence test_ (Zhang et al.,, 2011) as
implemented in the KCI function of the CondIndTests R package (Heinze-Deml et
al.,, 2018); we use the Bayesian hyperparameter tuning option for sample sizes
of at most $500$.

  * gcm

The _Generalised Covariance Measure_ (GCM) as described in Shah and Peters,
(2020).

  * wgcm.fix

The ‘fixed weight function’ variant of the _Weighted Generalised Covariance
Measure_ (wGCM) (Scheidegger et al.,, 2021) as implemented in the wgcm.fix
function of the weightedGCM R package; we use $\texttt{weight.num}=7$ as in
the simulations of the original paper.

  * wgcm.est

The ‘estimated weight function’ variant of the wGCM as implemented in the
wgcm.est function of the weightedGCM R-package.

In all of our numerical simulations, rejection rates were estimated based on
100 repetitions.

###  6.1 Additive models

We first investigate Type I error control and power in setting where both
${\mathbb{E}}(Y\,|\,Z)$ and ${\mathbb{E}}(X\,|\,Z)$ are additive functions,
and $Z\sim N_{7}(0,\boldsymbol{I})$. For the methods, including the PCM,
requiring choices of regression procedures, we use an additive model fitted
using mgcv. We consider null settings consisting of $n\in\\{250,500,1000\\}$
independent and identically distributed copies of $(X,Y,Z)$ where

$X=\sin(2\pi Z_{1})+0.1\xi\quad\text{and}\quad Y=\sin(2\pi
Z_{1})+\varepsilon.$

and errors $\varepsilon$ and $\xi$ are independent $N(0,1)$ random variables,
independent of $Z$. Such a setup is challenging for Type I error control as
$X$ and $Y$ are highly correlated yet are conditionally independent given $Z$.
Indeed we see from the left panel of Figure 2 that several of the tests are
anti-conservative, most notably kci and gam, which we omit from further
comparisons as their power properties would be hard to interpret given the
high rejection rates under the null. Both the williamson and wgcm.est tests
are also somewhat anti-conservative, but considerably less so. In contrast,
the pcm is conservative here. This is to be expected as the calibration
following the multiple sample splits involved in its construction (Section
2.2.3) is typically conservative; the single split version pcm_ss appears to
have rejection rates close to the 5% mark as suggested by our theory.

We investigate the power properties of the PCM in the following settings,
where as before, $\varepsilon$ and $\xi$ are independent and independent of
$Z$, and moreover $\varepsilon\sim N(0,1)$.

  1. 1.

$\xi\sim N(0,1)$, $X=\sin(2\pi Z_{1})+\xi$ and $Y=\sin(2\pi
Z_{1})+0.2X^{2}+\varepsilon$.

  2. 2.

$\xi+1\sim\textrm{Exp}(1)$, $X=\sin(2\pi Z_{1})-\sin(2\pi Z_{1})\xi$ and
$Y=\sin(2\pi Z_{1})+0.4X^{2}+\varepsilon$.

  3. 3.

$\xi\sim N(0,1)$, $X=\sin(2\pi Z_{1})+\xi$ and $Y=\sin(2\pi
Z_{1})+0.4X^{2}Z_{2}+\varepsilon$.

The settings are chosen such that in setting 1:
${\mathbb{E}}\bigl{(}{\mathrm{Cov}}(X,Y\,|\,Z)\bigr{)}=0$ but
${\mathrm{Cov}}(X,Y\,|\,Z)\neq 0$, in setting 2: ${\mathrm{Cov}}(X,Y\,|\,Z)=0$
but $\tau\neq 0$ and in setting 3: there is only an interaction effect. We
cannot expect this interaction effect to be picked up by methods that fit
additive models, but we nevertheless include this setting to emphasise the
fact that the success of the PCM and related methods is contingent on an
appropriate choice of regression method; see also Section 6.2.

From the right-hand panels of Figure 2, we see that the pcm and williamson
exhibit good power in settings 1 and 2, with the latter also rejecting in
setting 3; however the rejection rates for williamson should be interpreted
carefully given the miscalibration in the null settings. wgcm.est also shows
appreciable power in setting 1, though as expected has little power in setting
2 where ${\mathrm{Cov}}(X,Y\,|\,Z)=0$. For the reasons explained above, the
PCM has no power in setting 3; in the next section we investigate the
performance of the PCM when used in conjunction with a regression method
capable of fitting to such regression functions.

Figure 2: Rejection rates in the various settings considered in Section 6.1
for nominal 5%-level tests. Note that Alternative 3 has only an interaction
effect, so we cannot expect methods that fit additive models to have power.

###  6.2 Non-additive models

In this section, we consider settings where the regression functions are non-
additive and involve complex interactions. We use random forests (Breiman,,
2001) implemented in the ranger R package (Wright and Ziegler,, 2017) as our
regression procedure for the methods considered.

We consider null settings consisting of $n\in\\{10^{4},2\cdot 10^{4},4\cdot
10^{4}\\}$ independent and identically distributed copies of $(X,Y,Z)$ where
$Z\sim N_{7}(0,\boldsymbol{I})$ as before,

$X=\sin(2\pi Z_{1})(1+Z_{3})+\xi\quad\text{and}\quad\sin(2\pi
Z_{1})(1+Z_{3})+v(X)\varepsilon$

with $\varepsilon$ and $\xi$ independent $N(0,1)$ random variables independent
of $Z$, and $v(X):=0.5+\mathbbm{1}_{\\{X>0\\}}$ giving heteroscedastic errors
for the $Y$ regression model. The larger sample sizes considered here reflect
the difficulty of estimating the more complicated regression functions in
these examples. Note that here we do not have
$X\mbox{${}\perp\mkern-11.0mu\perp{}$}Y\,|\,Z$, but the conditional mean
independence ${\mathbb{E}}(Y\,|\,X,Z)={\mathbb{E}}(Y\,|\,Z)$ does hold. The
results are presented in Figure 3. We see that the multiple sample splitting
version of the PCM maintains Type I error control, and is in fact slightly
conservative. All other approaches considered appear to be anti-conservative
to varying degrees: the williamson approach is most clearly miscalibrated
here, and we omit it from our alternative settings described below; wgcm.est
and gcm are also fairly anti-conservative here but the rejection rates appear
to be improving for increasing $n$.

We consider the following alternative settings, where as in Section 6.1,
setting 2 has
${\mathrm{Cov}}(X,Y\,|\,Z)=0$, and setting 3 involves a pure interaction
effect:

  1. 1.

$\xi\sim N(0,1)$, $X=\sin(2\pi Z_{1})(1+Z_{3})+\xi$ and $Y=\sin(2\pi
Z_{1})(1+Z_{3})+0.04X^{2}+v(X)\varepsilon$;

  2. 2.

$\xi\sim\textrm{Exp}(1)$, $X=\sin(2\pi Z_{1})(1+Z_{3})-\sin(2\pi
Z_{1})(\xi-1)$ and $Y=\sin(2\pi Z_{1})(1+Z_{3})+0.04X^{2}+v(X)\varepsilon$;

  3. 3.

$\xi\sim N(0,1)$, $X=\sin(2\pi Z_{1})(1+Z_{3})+\xi$ and $Y=\sin(2\pi
Z_{1})(1+Z_{3})+0.04X^{2}Z_{2}+v(X)\varepsilon$.

Among the methods considered, here only the PCM appear to have good power
across the settings considered. The wgcm.est has reasonable power in setting
1, though this should be interpreted with some care given that Type I error is
not very well controlled in the null settings. However in setting 2, wgcm.est
is powerless as expected.

Figure 3: Rejection rates in the various settings considered in Section 6.2
for nominal 5%-level tests.

##  7 Conclusion

In this work we have introduced a general test statistic called the PCM for
testing conditional mean independence that: (a) can leverage machine learning
methods to yield provable uniform Type I error control across a class of null
distributions where these methods have sufficiently good predictive ability;
and (b) when used in conjunction with appropriate regression methods attains
rate-optimal power in both the parametric setting of the linear model and
fully nonparametric settings. We believe the PCM fills an important gap in the
data analyst’s range of existing tools, which are unable simultaneously to
achieve these desiderata. However, our work also offers several avenues for
further work, some of which we mention below.

##### Verifying the general assumptions for other regression methods

We have verified Assumption 3 for linear regression in linear model settings,
and nonparametric series estimators in fully nonparametric settings. Since we
used the penalised regression splines of mgcv in several of our numerical
experiments, it would be interesting to see for what classes $\mathcal{P}$ of
distributions Assumption 3 is satisfied in that context. Similarly, it would
be very interesting to ask the same question of random forests, which perform
very well in our simulations; however this is likely to be challenging given
the complex nature of the random forest procedure.

##### Aggregation of test statistics from multiple sample splits

While our proposal (Algorithm 2) to average test statistics from multiple
sample splits and compare this to a standard Gaussian quantile works well in
practice, determining interpretable conditions under which this guarantees
Type I error control is an open problem. Moreover, our numerical experiments
suggest that our proposal is somewhat conservative, and it is certainly of
interest, both in our setting and other contexts involving multiple sample
splitting, to develop an aggregation procedure that makes more efficient use
of the sample to attain greater power.

##### Conditional independence testing

Although the problem of testing conditional independence has been studied more
intensively than that of testing conditional mean independence, there do not
exist many practical conditional independence tests that achieve the two
desiderata mentioned at the beginning of this section. One starting point for
constructing a such a test may be the fact that the conditional independence
null $Y\mbox{${}\perp\mkern-11.0mu\perp{}$}X\,|\,Z$ may be viewed as the
intersection of conditional mean independence nulls
${\mathbb{E}}\bigl{(}w(Y)\,|\,X,Z\bigr{)}={\mathbb{E}}\bigl{(}w(Y)\,|\,Z\bigr{)}$
where function $w$ ranges over all monotone functions, for example. It might
therefore be interesting to investigate procedures that seek two
‘projections’: mappings $(X,Z)\to\widehat{f}(X,Z)$ and also
$Y\mapsto\widehat{w}(Y)$, after which one may apply the GCM.

##### Confidence intervals

We have focused on the problem of testing conditional mean independence, but
the problem of deriving confidence intervals for a parameter such as $\tau$
that is $0$ under our null is equally interesting. The pioneering work of
Williamson et al., (2021) proposes an asymptotically optimal approach for this
in the case where $\tau$ is bounded away from $0$. It would be interesting if
the PCM could be used in conjunction with the proposal of Williamson et al.,
(2021) to extend the latter to yield confidence intervals with uniform
coverage for all $\tau$.

## References

  * Aït-Sahalia et al., (2001) Aït-Sahalia, Y., Bickel, P. J., and Stoker, T. M. (2001).  Goodness-of-fit tests for kernel regression with an application to option implied volatilities.  Journal of Econometrics, 105(2):363–412.
  * Arias-Castro et al., (2018) Arias-Castro, E., Pelletier, B., and Saligrama, V. (2018).  Remember the curse of dimensionality: The case of goodness-of-fit testing in arbitrary dimension.  Journal of Nonparametric Statistics, 30(2):448–471.
  * Belloni et al., (2015) Belloni, A., Chernozhukov, V., Chetverikov, D., and Kato, K. (2015).  Some new asymptotic theory for least squares series: Pointwise and uniform results.  Journal of Econometrics, 186(2):345–366.
  * Belloni et al., (2011) Belloni, A., Chernozhukov, V., and Wang, L. (2011).  Square-root lasso: pivotal recovery of sparse signals via conic programming.  Biometrika, 98(4):791–806.
  * Bengs and Holzmann, (2019) Bengs, V. and Holzmann, H. (2019).  Uniform approximation in classical weak convergence theory.  arXiv preprint arXiv:1903.09864.
  * Berrett et al., (2020) Berrett, T. B., Wang, Y., Barber, R. F., and Samworth, R. J. (2020).  The conditional permutation test for independence while controlling for confounders.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(1):175–197.
  * Bojer and Meldgaard, (2021) Bojer, C. S. and Meldgaard, J. P. (2021).  Kaggle forecasting competitions: An overlooked learning opportunity.  International Journal of Forecasting, 37(2):587–603.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine Learning, 45(1):5–32.
  * Bühlmann and van de Geer, (2011) Bühlmann, P. and van de Geer, S. (2011).  Statistics for High-Dimensional Data: Methods, Theory and Applications.  Springer Publishing Company, Incorporated, 1st edition.
  * Cai et al., (2022) Cai, Z., Lei, J., and Roeder, K. (2022).  Model-free prediction test with application to genomics data.  bioRxiv preprint bioRxiv.2022.03.28.486116.
  * Candès et al., (2018) Candès, E., Fan, Y., Janson, L., and Lv, J. (2018).  Panning for gold: ‘model-X’ knockoffs for high dimensional controlled variable selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(3):551–577.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1).
  * Cox, (1975) Cox, D. R. (1975).  A note on data-splitting for the evaluation of significance levels.  Biometrika, 62(2):441–444.
  * Dai et al., (2022) Dai, B., Shen, X., and Pan, W. (2022+).  Significance tests of feature relevance for a blackbox learner.  IEEE Transactions on Neural Networks and Learning Systems, to appear.
  * de Boor, (1976) de Boor, C. (1976).  Splines as linear combinations of b-splines. a survey.
  * DiCiccio et al., (2020) DiCiccio, C. J., DiCiccio, T. J., and Romano, J. P. (2020).  Exact tests via multiple data splitting.  Statistics & Probability Letters, 166:108865.
  * Durrett, (2019) Durrett, R. (2019).  Probability: Theory and Examples.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 5th edition.
  * Fan and Li, (1996) Fan, Y. and Li, Q. (1996).  Consistent model specification tests: omitted variables and semiparametric functional forms.  Econometrica, 64(4):865–890.
  * Fernández and Rivera, (2022) Fernández, T. and Rivera, N. (2022).  A general framework for the analysis of kernel-based tests.  arXiv preprint arXiv:2209.00124.
  * Goodfellow et al., (2016) Goodfellow, I., Bengio, Y., and Courville, A. (2016).  Deep Learning.  MIT press.
  * Gut, (2013) Gut, A. (2013).  Probability: A Graduate Course.  Springer New York.
  * Heinze-Deml et al., (2018) Heinze-Deml, C., Peters, J., and Meinshausen, N. (2018).  Invariant causal prediction for nonlinear models.  Journal of Causal Inference, 6(2):20170016.
  * Huang, (2003) Huang, J. Z. (2003).  Local asymptotics for polynomial spline regression.  The Annals of Statistics, 31(5):1600 – 1635.
  * Ichimura and Newey, (2015) Ichimura, H. and Newey, W. K. (2015).  The influence function of semiparametric estimators.  arXiv preprint arXiv:1508.01378.
  * Ingster, (1987) Ingster, Y. I. (1987).  Minimax testing of nonparametric hypotheses on a distribution density in the $L_{p}$ metrics.  Theory of Probability & Its Applications, 31(2):333–337.
  * Ingster, (2000) Ingster, Y. I. (2000).  Adaptive chi-square tests.  Journal of Mathematical Sciences, 99(2):1110–1119.
  * Janková et al., (2020) Janková, J., Shah, R. D., Bühlmann, P., and Samworth, R. J. (2020).  Goodness-of-fit testing in high dimensional generalized linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(3):773–795.
  * Jin et al., (2018) Jin, Z., Yan, X., and Matteson, D. S. (2018).  Testing for conditional mean independence with covariates through martingale difference divergence.  In Proceedings of the Thirty-Fourth Conference on Uncertainty in Artificial Intelligence, UAI 2018, Monterey, California, USA, August 6-10, 2018, pages 1–12. AUAI Press.
  * Kasy, (2019) Kasy, M. (2019).  Uniformity and the delta method.  Journal of Econometric Methods, 8(1).
  * Kennedy, (2020) Kennedy, E. H. (2020).  Optimal doubly robust estimation of heterogeneous causal effects.  arXiv preprint arXiv:2004.14497.
  * Kim and Ramdas, (2020) Kim, I. and Ramdas, A. (2020).  Dimension-agnostic inference using cross $U$-statistics.  arXiv preprint arXiv:2011.05068.
  * Klenke, (2020) Klenke, A. (2020).  Probability Theory.  Springer International Publishing.
  * Lavergne and Vuong, (2000) Lavergne, P. and Vuong, Q. (2000).  Nonparametric significance testing.  Econometric Theory, 16(4):576–601.
  * Lepskiǐ, (1991) Lepskiǐ, O. V. (1991).  Asymptotically minimax adaptive estimation i: Upper bounds. optimally adaptive estimates.  Theory Probab. Appl., 36:682–697.
  * Lundborg et al., (2022) Lundborg, A. R., Shah, R. D., and Peters, J. (2022+).  Conditional independence testing in Hilbert spaces with applications to functional data analysis.  J. Roy. Statist. Soc., Ser. B., to appear.
  * Meinshausen and Bühlmann, (2010) Meinshausen, N. and Bühlmann, P. (2010).  Stability selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(4):417–473.
  * Meinshausen et al., (2009) Meinshausen, N., Meier, L., and Bühlmann, P. (2009).  P-values for high-dimensional regression.  Journal of the American Statistical Association, 104(488):1671–1681.
  * Newey and Robins, (2018) Newey, W. K. and Robins, J. R. (2018).  Cross-fitting and fast remainder rates for semiparametric estimation.  arXiv preprint arXiv:1801.09138.
  * Neykov et al., (2021) Neykov, M., Balakrishnan, S., and Wasserman, L. (2021).  Minimax optimal conditional independence testing.  The Annals of Statistics, 49(4):2151–2177.
  * Neyman, (1923) Neyman, J. (1923).  Sur les applications de la théorie des probabilités aux experiences agricoles: Essai des principes.  Roczniki Nauk Rolniczych, 10:1–51.
  * Petersen and Hansen, (2021) Petersen, L. and Hansen, N. R. (2021).  Testing conditional independence via quantile regression based partial copulas.  Journal of Machine Learning Research, 22(70):1–47.
  * Powell, (1981) Powell, M. J. D. (1981).  Approximation Theory and Methods.  Cambridge University Press.
  * Rinaldo et al., (2019) Rinaldo, A., Wasserman, L., and G’Sell, M. (2019).  Bootstrapping and sample splitting for high-dimensional, assumption-lean inference.  The Annals of Statistics, 47(6):3438–3469.
  * Romano, (2004) Romano, J. P. (2004).  On Non-parametric Testing, the Uniform Behaviour of the $t$-test, and Related Problems.  Scandinavian Journal of Statistics, 31(4):567–584.
  * Rubin, (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66(5):688.
  * Scheidegger et al., (2021) Scheidegger, C., Hörrmann, J., and Bühlmann, P. (2021).  The Weighted Generalised Covariance Measure.  arXiv preprint arXiv:2111.04361.
  * Schumaker, (2007) Schumaker, L. (2007).  Spline Functions: Basic Theory.  Cambridge Mathematical Library. Cambridge University Press, 3 edition.
  * Shah and Bühlmann, (2018) Shah, R. D. and Bühlmann, P. (2018).  Goodness-of-fit tests for high dimensional linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(1):113–135.
  * Shah and Peters, (2020) Shah, R. D. and Peters, J. (2020).  The hardness of conditional independence testing and the generalised covariance measure.  The Annals of Statistics, 48(3):1514–1538.
  * Shah and Samworth, (2013) Shah, R. D. and Samworth, R. J. (2013).  Variable selection with error control: another look at stability selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(1):55–80.
  * Shao and Zhang, (2014) Shao, X. and Zhang, J. (2014).  Martingale difference correlation and its use in high-dimensional variable screening.  Journal of the American Statistical Association, 109(507):1302–1318.
  * Sun and Zhang, (2012) Sun, T. and Zhang, C.-H. (2012).  Scaled sparse linear regression.  Biometrika, 99(4):879–898.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288.
  * van der Vaart, (1998) van der Vaart, A. W. (1998).  Asymptotic Statistics.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Verdinelli and Wasserman, (2021) Verdinelli, I. and Wasserman, L. (2021).  Decorrelated Variable Importance.  arXiv preprint arXiv:2111.10853.
  * Vovk and Wang, (2020) Vovk, V. and Wang, R. (2020).  Combining $p$-values via averaging.  Biometrika, 107(4):791–808.
  * Wang and Samworth, (2018) Wang, T. and Samworth, R. J. (2018).  High dimensional change point estimation via sparse projection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(1):57–83.
  * Wang and Shah, (2020) Wang, Y. and Shah, R. D. (2020).  Debiased Inverse Propensity Score Weighting for Estimation of Average Treatment Effects with High-Dimensional Confounders.  arXiv preprint arXiv:2011.08661.
  * Wasserman et al., (2020) Wasserman, L., Ramdas, A., and Balakrishnan, S. (2020).  Universal inference.  Proceedings of the National Academy of Sciences, 117(29):16880–16890.
  * Wasserman and Roeder, (2009) Wasserman, L. and Roeder, K. (2009).  High dimensional variable selection.  The Annals of Statistics, 37(5A):2178–2201.
  * White, (1980) White, H. (1980).  A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity.  Econometrica, 48(4):817–838.
  * Williamson et al., (2021) Williamson, B. D., Gilbert, P. B., Carone, M., and Simon, N. (2021).  Nonparametric variable importance assessment using machine learning techniques.  Biometrics, 77:9–22.
  * Williamson et al., (2022) Williamson, B. D., Gilbert, P. B., Simon, N. R., and Carone, M. (2022+).  A general framework for inference on algorithm-agnostic variable importance.  Journal of the American Statistical Association, to appear.
  * Wood, (2013) Wood, S. N. (2013).  On p-values for smooth components of an extended generalized additive model.  Biometrika, 100(1):221–228.
  * Wood, (2017) Wood, S. N. (2017).  Generalized Additive Models.  Chapman and Hall/CRC.
  * Wright and Ziegler, (2017) Wright, M. N. and Ziegler, A. (2017).  ranger: A fast implementation of random forests for high dimensional data in C++ and R.  Journal of Statistical Software, 77(1):1–17.
  * Zeileis and Hothorn, (2002) Zeileis, A. and Hothorn, T. (2002).  Diagnostic checking in regression relationships.  R News, 2(3):7–10.
  * Zhang et al., (2011) Zhang, K., Peters, J., Janzing, D., and Schölkopf, B. (2011).  Kernel-Based Conditional Independence Test and Application in Causal Discovery.  In Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence, UAI’11, pages 804–813, Arlington, Virginia, USA. AUAI Press.
  * Zhang and Janson, (2020) Zhang, L. and Janson, L. (2020).  Floodgate: inference for model-free variable importance.  arXiv preprint arXiv:2007.01283.

In Sections A and B of the appendix, we include the proofs of all of our main
results and related auxiliary lemmas. In Section C, we give a self-contained
description of spline regression and related results that we use for our
analysis in Section 5. Section D contains an analysis of the linear
projections in Section 3 under more general assumptions. We also derive an
exact asymptotic power function of our test in this setting. Section E
contains the results from additional numerical experiments beyond those
included in Section 6.

Throughout the appendix, for a positive semi-definite matrix $\boldsymbol{A}$,
we write $\boldsymbol{A}^{-1}$ for its generalised inverse (i.e. its
Moore–Penrose pseudo-inverse).

##  Appendix A Proofs

In our proofs we often suppress the dependence of quantities on $P$ for ease
of notation.

###  A.1 Proof of Proposition 1

There is no loss of generality in assuming that $\alpha\leq 1/2$, because for
$\alpha>1/2$, we have
$\mathbb{P}(T>z_{1-\alpha})\geq\mathbb{P}(T>z_{1/2})=\mathbb{P}(T>0)$. We
start by checking the assumptions of Lemma 18 for the regressions of $Y$ on
$X$ and $Z$ (of which $\widehat{\beta}$ is one component), $Y$ on $Z$
(yielding $\boldsymbol{\widehat{\theta}}$) and $X$ on $Z$ (yielding
$\boldsymbol{\widehat{\eta}}$). Let $\delta^{\prime}:=\delta/2$.

Recalling that $W=(X,Z)$, we see that condition (i) of Lemma 18 is satisfied
for the $Y$ on $X$ and $Z$ regression by our assumption on
$\boldsymbol{\Sigma}^{XZ}$ and the fact that

$\lambda_{\min}\bigl{(}{\mathbb{E}}(WW^{\top}\zeta^{2})\bigr{)}\geq
c\lambda_{\min}(\boldsymbol{\Sigma}^{XZ})$

by our assumption that ${\mathrm{Var}}(Y\,|\,X,Z)\geq c$. Condition (ii) is
satisfied with $\delta=\delta^{\prime}$ by the Cauchy–Schwarz inequality and
Jensen’s inequality.

By Assumption 1, condition (i) of Lemma 18 is satisfied for the $X$ on $Z$
regression. To see that condition (ii) is satisfied with
$\delta=\delta^{\prime}$, we note that by the Cauchy–Schwarz inequality it
suffices to check that ${\mathbb{E}}(|\xi|^{4+\delta})$ is bounded over
$\mathcal{P}$. Letting $\boldsymbol{\Sigma}:={\mathbb{E}}(ZZ^{\top})$, we have

${\mathbb{E}}(|\xi|^{4+\delta})\leq
2^{3+\delta}\Bigl{(}{\mathbb{E}}(|X|^{4+\delta})+\lambda_{\min}(\boldsymbol{\Sigma})^{-(4+\delta)}\|{\mathbb{E}}(XZ)\|_{2}^{4+\delta}{\mathbb{E}}(\|Z\|_{2}^{4+\delta})\Bigr{)}$

which is bounded under Assumption 1.

To see that condition (i) of Lemma 18 is satisfied for the $Y$ on $Z$
regression, define
$\boldsymbol{\theta}:={\mathbb{E}}(ZZ^{\top})^{-1}{\mathbb{E}}(YZ)\in{\mathbb{R}}^{d}$
and note that

$\displaystyle{\mathbb{E}}\bigl{(}ZZ^{\top}(Y-\boldsymbol{\theta}^{\top}Z)^{2}\bigr{)}$
$\displaystyle={\mathbb{E}}\bigl{(}ZZ^{\top}(\zeta+\beta
X-\beta\boldsymbol{\eta}^{\top}Z)^{2}\bigr{)}$
$\displaystyle={\mathbb{E}}(ZZ^{\top}\zeta^{2})+\beta^{2}{\mathbb{E}}\bigl{(}ZZ^{\top}(X-\boldsymbol{\eta}^{\top}Z)^{2}\bigr{)},$

so the minimum eigenvalue of
${\mathbb{E}}\bigl{(}ZZ^{\top}(Y-\boldsymbol{\theta}^{\top}Z)^{2}\bigr{)}$ is
bounded below by $c\lambda_{\min}(\boldsymbol{\Sigma})$. Condition (ii)
follows by similar arguments as those for the $X$ on $Z$ regression. We
therefore deduce from Lemma 18 that (69) holds with $\sigma_{\beta}^{2}$ given
by the $(1,1)$th entry of the matrix
$(\boldsymbol{\Sigma}^{XZ})^{-1}\mathbb{E}(WW^{\top}\zeta^{2})(\boldsymbol{\Sigma}^{XZ})^{-1}$,
and that

$\sqrt{n}\|\boldsymbol{\widehat{\eta}}-\boldsymbol{\eta}\|_{2}=O_{\mathcal{P}}(1)\quad\text{and}\quad\sqrt{n}\|\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta}\|_{2}=O_{\mathcal{P}}(1).$
(21)

We now verify that the remaining parts of Assumption 6 are satisfied. First,

$\sup_{P\in\mathcal{P}}{\mathbb{E}}\biggl{(}\Bigl{\|}\frac{1}{n}\sum_{i=1}^{n}Z_{i}\xi_{i}\Bigr{\|}_{2}^{2}\biggr{)}=\frac{1}{n}\sup_{P\in\mathcal{P}}{\mathbb{E}}\bigl{(}\|Z\xi\|_{2}^{2}\bigr{)}\leq\frac{1}{n}\sup_{P\in\mathcal{P}}{\mathbb{E}}(\|Z\|_{2}^{2})\cdot\sup_{P\in\mathcal{P}}{\mathbb{E}}(\xi^{2})\to
0,$

so by Lemma 14 and (21) we have that (70) holds. Similar arguments show that
(71) is satisfied.

Next, (21) shows that

$\sqrt{n}\|\boldsymbol{\widehat{\eta}}-\boldsymbol{\eta}\|_{2}\|\boldsymbol{\widehat{\theta}}-\boldsymbol{\theta}\|_{2}=O_{\mathcal{P}}(n^{-1/2})=o_{\mathcal{P}}(1).$

Moreover, by (51) in the proof of Lemma 18 and Assumption 1, we have

$\Bigl{\|}\frac{1}{n}\sum_{i=1}^{n}Z_{i}Z_{i}^{\top}\Bigr{\|}_{\mathrm{op}}\leq\Bigl{\|}\frac{1}{n}\sum_{i=1}^{n}Z_{i}Z_{i}^{\top}-\boldsymbol{\Sigma}\Bigr{\|}_{\mathrm{op}}+\|\boldsymbol{\Sigma}\|_{\mathrm{op}}=O_{\mathcal{P}}(1),$

so (72) is satisfied. The remaining conditions hold by similar arguments using
the moment bounds established earlier, (21) and Lemma 17.

To verify the remaining conditions of Proposition 42, we note that

${\mathrm{Var}}(\varepsilon\xi)\geq{\mathbb{E}}\bigl{\\{}{\mathrm{Var}}(\varepsilon\xi\,|\,X,Z)\bigr{\\}}={\mathbb{E}}\bigl{\\{}\xi^{2}{\mathrm{Var}}(Y\,|\,X,Z)\bigr{\\}}\geq
c^{2}.$

Finally, the moment bound condition in Proposition 42 follows by
Cauchy–Schwarz and the arguments above. The conclusion now follows from
Proposition 42 together with the fact that $\psi_{\alpha,n}$ is an increasing
function of $|\beta|$, so that

$\psi_{\alpha,n}\geq\Phi\biggl{(}\frac{\kappa}{\sigma_{\beta}}\biggr{)}\cdot\Phi\biggl{(}z_{\alpha}+\frac{\kappa\sigma_{\xi}^{2}}{\sigma_{\varepsilon\xi}}\biggr{)}+\Phi\biggl{(}-\frac{\kappa}{\sigma_{\beta}}\biggr{)}\cdot\Phi\biggl{(}z_{\alpha}-\frac{\kappa\sigma_{\xi}^{2}}{\sigma_{\varepsilon\xi}}\biggr{)}\rightarrow
1$

as $\kappa\rightarrow\infty$, as required.

###  A.2 Proof of Proposition 2

Throughout this proof we work on the event that at least one
$u_{n,1},\ldots,u_{n,n}$ is non-zero which is a set of uniform asymptotic
probability $1$ by Assumption 2(b). Let
$\boldsymbol{Z}:=(Z_{1}^{\top},\ldots,Z_{n}^{\top})^{\top}$,
$\boldsymbol{P}:=\boldsymbol{Z}(\boldsymbol{Z}^{\top}\boldsymbol{Z})^{-1}\boldsymbol{Z}^{\top}$,
$\boldsymbol{Y}:=(Y_{1},\ldots,Y_{n})^{\top}$,
$\boldsymbol{\varepsilon}:=(\varepsilon_{1},\ldots,\varepsilon_{n})^{\top}$,
$\boldsymbol{\widehat{f}}:=\bigl{(}\widehat{f}(X_{1},Z_{1}),\dots,\widehat{f}(X_{n},Z_{n})\bigr{)}^{\top}$
and $\boldsymbol{I}$ denote the $d\times d$ identity matrix. Since
$\boldsymbol{P}$ is a matrix representing an orthogonal projection such that
$\boldsymbol{Z}^{\top}(\boldsymbol{I}-\boldsymbol{P})$ is a zero vector, we
have

$\displaystyle\sum_{i=1}^{n}\\{Y_{i}-\boldsymbol{\widehat{\gamma}}^{\top}Z_{i}\\}\\{\widehat{f}(X_{i},Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})\\}$
$\displaystyle=\boldsymbol{\widehat{f}}^{\top}(\boldsymbol{I}-\boldsymbol{P})\boldsymbol{Y}=\boldsymbol{\widehat{f}}^{\top}(\boldsymbol{I}-\boldsymbol{P})\boldsymbol{\varepsilon}$
$\displaystyle=\sum_{i=1}^{n}\varepsilon_{i}\\{\widehat{f}(X_{i},Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})\\}.$

Based on the above identity, we have that

$\displaystyle
T_{\mathrm{OLS}}=\frac{\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}\varepsilon_{i}u_{n,i}}{\sqrt{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}(Y_{i}-\boldsymbol{\widehat{\gamma}}^{\top}Z_{i})^{2}u_{n,i}^{2}-\bigl{(}\frac{1}{n\nu}\sum_{i=1}^{n}\varepsilon_{i}u_{n,i}\bigr{)}^{2}}},$

where

$\nu:=\sqrt{\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}(\varepsilon_{i}^{2}\,|\,X_{i},Z_{i})u_{n,i}^{2}}>0.$

Let $\mathcal{F}_{n}$ denote the $\sigma$-algebra generated by $\widehat{f}$
and $\bigl{(}(X_{i},Z_{i})\bigr{)}_{i=1}^{n}$. Form the triangular array

$W_{n,i}:=\frac{\varepsilon_{i}u_{n,i}}{\nu}$

for $n\in\mathbb{N}$ and $i\in[n]$, and note that this satisfies the first
three assumptions of Lemma 16, since $u_{n,i}$ is measurable with respect to
$\mathcal{F}_{n}$. Finally, the fourth assumption of this lemma is also
satisfied, because

$\displaystyle\frac{1}{n^{1+\delta/2}}$
$\displaystyle\sum_{i=1}^{n}{\mathbb{E}}(|W_{n,i}|^{2+\delta}\,|\,\mathcal{F}_{n})=\frac{1}{n^{1+\delta/2}\nu^{2+\delta}}\sum_{i=1}^{n}{\mathbb{E}}(|\varepsilon_{i}|^{2+\delta}\,|\,X_{i},Z_{i})|u_{n,i}|^{2+\delta}$
$\displaystyle\leq\frac{C}{c^{1+\delta/2}}\sum_{i=1}^{n}|\upsilon_{n,i}|^{2+\delta}\leq\frac{C}{c^{1+\delta/2}}\biggl{(}\sum_{i=1}^{n}|\upsilon_{n,i}|^{2}\biggr{)}\max_{i\in[n]}|\upsilon_{n,i}|^{\delta}=o_{\mathcal{P}_{0}}(1)$

by Assumptions 2(a) and (b) and Lemmas 11 and 15. Lemma 16 thus yields that
the numerator of $T_{\mathrm{OLS}}$ is uniformly asymptotically standard
Gaussian.

For the denominator of $T_{\mathrm{OLS}}$, the uniform version of Slutsky’s
theorem (Bengs and Holzmann,, 2019, Theorem 6.3) yields that
$\frac{1}{n\nu}\sum_{i=1}^{n}\varepsilon_{i}u_{n,i}=o_{\mathcal{P}_{0}}(1)$.
For the first term in the denominator of $T_{\mathrm{OLS}}$, we consider the
decomposition

$\displaystyle\frac{1}{n\nu^{2}}\sum_{i=1}^{n}(Y_{i}-\boldsymbol{\widehat{\gamma}}^{\top}Z_{i})^{2}u_{n,i}^{2}=\underbrace{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\varepsilon_{i}^{2}u_{n,i}^{2}}_{\textup{I}_{n}}$
$\displaystyle+\underbrace{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\\{(\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma})^{\top}Z_{i}\\}^{2}u_{n,i}^{2}}_{\textup{II}_{n}}$
$\displaystyle-\underbrace{\frac{2}{n\nu^{2}}\sum_{i=1}^{n}(\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma})^{\top}Z_{i}\varepsilon_{i}u_{n,i}^{2}}_{\textup{III}_{n}}.$

Similarly to our previous argument, define the triangular array
$\tilde{W}_{n,i}:=W_{n,i}^{2}$ for $n\in\mathbb{N}$ and $i\in[n]$, and note
that $(\tilde{W}_{n,i})_{n\in\mathbb{N},i\in[n]}$ satisfies the conditions of
Lemma 17 with $\mu_{n}=1$ in that result, so
$\textup{I}_{n}=1+o_{\mathcal{P}_{0}}(1)$. Now, by Hölder’s inequality,

$\displaystyle|\textup{II}_{n}|\leq\frac{1}{n\nu^{2}}\|\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}\|_{1}^{2}\sum_{i=1}^{n}\|Z_{i}\|_{\infty}^{2}u_{n,i}^{2}\leq\frac{1}{c}\max_{i\in[n]}\|Z_{i}\|_{\infty}^{2}\|\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}\|_{1}^{2}\sum_{i=1}^{n}\upsilon_{n,i}^{2}=o_{\mathcal{P}}(1),$

by Assumption 2 and Lemma 11. Finally, the Cauchy–Schwarz inequality yields
that

$\displaystyle|\textup{III}_{n}|$ $\displaystyle\leq
2\sqrt{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\varepsilon_{i}^{2}u_{n,i}^{2}}\cdot\sqrt{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\\{(\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma})^{\top}Z_{i}\\}^{2}u_{n,i}^{2}}=2\sqrt{\textup{I}_{n}}\cdot\sqrt{\textup{II}_{n}}=o_{\mathcal{P}}(1)$

by Lemmas 11 and 15. The result follows by the uniform version of Slutsky’s
theorem.

###  A.3 Proof of Proposition 3

As in the proof of Proposition 2, we work on the event that at least one
$u_{n,1},\ldots,u_{n,n}$ is non-zero, which is a set of uniform asymptotic
probability $1$ by Assumption 2(b). Recall the definitions of $\nu$ from the
proof of Proposition 2, and $\delta_{\mathrm{bias}}$ from just after (9). Our
test statistic can be written as

$\displaystyle
T_{\mathrm{Lasso}}=\frac{\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}\varepsilon_{i}u_{n,i}-\frac{1}{\sqrt{n}\nu}\delta_{\mathrm{bias}}}{\sqrt{\frac{1}{ns_{n}^{2}}\sum_{i=1}^{n}(Y_{i}-\boldsymbol{\widehat{\gamma}}^{\top}Z_{i})^{2}u_{n,i}^{2}-\bigl{(}\frac{1}{n\nu}\sum_{i=1}^{n}\varepsilon_{i}u_{n,i}+\frac{1}{n\nu}\delta_{\mathrm{bias}}\bigr{)}^{2}}}.$

Following the proof of Proposition 2, we know that
$\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}\varepsilon_{i}u_{n,i}$ converges
uniformly to $N(0,1)$. Further, by Assumption 2(a) and Hölder’s inequality, we
have

$\biggl{|}\frac{1}{\sqrt{n}\nu}\delta_{\mathrm{bias}}\biggr{|}\leq\frac{1}{c^{1/2}}\biggl{|}\sum_{i=1}^{n}(\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma})^{\top}Z_{i}\upsilon_{n,i}\biggr{|}\leq\frac{1}{c^{1/2}}\|\boldsymbol{\widehat{\gamma}}-\boldsymbol{\gamma}\|_{1}\biggl{\|}\sum_{i=1}^{n}Z_{i}\upsilon_{n,i}\biggr{\|}_{\infty}=o_{\mathcal{P}_{0}}(1)$

under condition (10). A uniform version of Slutsky’s theorem (Bengs and
Holzmann,, 2019, Theorem 6.3) now yields that the numerator of
$T_{\mathrm{Lasso}}$ is uniformly asymptotically standard Gaussian. We can
repeat the arguments of Proposition 2 to show that the denominator is
$1+o_{\mathcal{P}_{0}}(1)$ under Assumption 2, so the uniform version of
Slutsky’s theorem yields the desired result.

###  A.4 Proof of Theorem 4

We prove the result under the given assumptions but instead of assuming that
$\widehat{m}$ and $\widehat{m}_{\widehat{f}}$ are formed on an auxiliary
sample, we let
$R_{ij}:={\mathbb{E}}\bigl{(}M_{i}M_{j}\,|\,(X_{i^{\prime}},Z_{i^{\prime}})_{i^{\prime}=1}^{n}\bigr{)}-{\mathbb{E}}\bigl{(}M_{i}M_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n}\bigr{)}$
for $i\neq j$, where $M_{i}:=m(Z_{i})-\widehat{m}(Z_{i})$, and assume that

$\frac{1}{n\sigma_{n}^{2}}\sum_{i\neq
j}\bigl{|}{\mathbb{E}}\bigl{(}R_{ij}\xi_{i}\xi_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n},\widehat{f}\bigr{)}\bigr{|}=o_{\mathcal{P}_{0}}(1).$
(22)

In Proposition 20 we show that this condition is satisfied if $\widehat{m}$ is
formed out of sample as in the statement of Theorem 4. We also show in that
result that (22) is satisfied if $\widehat{m}$ is a linear smoother or if
$X\mbox{${}\perp\mkern-11.0mu\perp{}$}Y\,|\,Z$. Define
$\nu^{2}:={\mathrm{Var}}(\varepsilon\xi)$ and note that $\nu^{2}\geq
c\sigma^{2}$. Throughout this proof we work on the event
$\Omega_{0}:=\\{\sigma\neq 0\\}$, which satisfies
$\mathbb{P}(\Omega_{0}^{c})=o_{\mathcal{P}_{0}}(1)$ by Assumption 3(a). Define

$T^{(\mathrm{N})}:=\frac{\frac{1}{\sqrt{n}}\sum_{i=1}^{n}L_{i}}{\nu}\quad\text{and}\quad
T^{(\mathrm{D})}:=\frac{\sqrt{\frac{1}{n}\sum_{i=1}^{n}L_{i}^{2}-\bigl{(}\frac{1}{n}\sum_{i=1}^{n}L_{i}\bigr{)}^{2}}}{\nu},$

so that $T=T^{(\mathrm{N})}/T^{(\mathrm{D})}$. We will show that
$T^{(\mathrm{N})}$ converges uniformly in distribution to $N(0,1)$ and
$|T^{(\mathrm{D})}-1|=o_{\mathcal{P}_{0}}(1)$, which yields the desired result
by combining Lemma 15 and the uniform version of Slutsky’s lemma (Bengs and
Holzmann,, 2019, Theorem 6.3).

Define
$\widetilde{M}_{i}:=m_{\widehat{f}}(Z_{i})-\widehat{m}_{\widehat{f}}(Z_{i})$
for $i\in[n]$ and note that

$T^{(\mathrm{N})}=\underbrace{\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}M_{i}\widetilde{M}_{i}}_{a_{n}}+\underbrace{\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}\widetilde{M}_{i}\varepsilon_{i}}_{b_{n}}+\underbrace{\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}M_{i}\xi_{i}}_{c_{n}}+\underbrace{\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}\varepsilon_{i}\xi_{i}}_{U_{n}}.$
(23)

By the Cauchy–Schwarz inequality,

$a_{n}\leq\sqrt{\frac{1}{cn}\biggl{(}\sum_{i=1}^{n}M_{i}^{2}\biggr{)}\biggl{(}\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}^{2}\biggr{)}}=o_{\mathcal{P}_{0}}(1),$

by Assumption 3(b).

To see that $b_{n}=o_{\mathcal{P}_{0}}(1)$, we note that

$b_{n}^{2}=\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}^{2}\varepsilon_{i}^{2}+\frac{1}{n\nu^{2}}\sum_{i\neq
j}\widetilde{M}_{i}\widetilde{M}_{j}\varepsilon_{i}\varepsilon_{j}$ (24)

and the first term is $o_{\mathcal{P}_{0}}(1)$ by Assumption 3(c). Now, for
$i\neq j$,

$\displaystyle{\mathbb{E}}(Y_{i}Y_{j}\,|\,X_{i},X_{j},Z_{i},Z_{j})={\mathbb{E}}(Y_{i}\,|\,X_{i},Z_{i}){\mathbb{E}}(Y_{j}\,|\,X_{j},Z_{j})=m(Z_{i})m(Z_{j}),$

using the fact that $m(Z)={\mathbb{E}}(Y\,|\,Z)={\mathbb{E}}(Y\,|\,X,Z)$ under
$\mathcal{P}_{0}$. Hence,

$\displaystyle{\mathbb{E}}(\varepsilon_{i}\varepsilon_{j}\,|\,X_{i},X_{j},Z_{i},Z_{j})$
$\displaystyle={\mathbb{E}}\bigl{\\{}\bigl{(}Y_{i}-m(Z_{i})\bigr{)}\bigl{(}Y_{j}-m(Z_{j})\bigr{)}\,|\,X_{i},X_{j},Z_{i},Z_{j}\bigr{\\}}$
$\displaystyle={\mathbb{E}}(Y_{i}Y_{j}\,|\,X_{i},X_{j},Z_{i},Z_{j})-m(Z_{i})m(Z_{j})=0.$

It follows that

$\displaystyle{\mathbb{E}}\Bigl{(}\frac{1}{n\nu^{2}}\sum_{i\neq
j}\widetilde{M}_{i}\widetilde{M}_{j}\varepsilon_{i}\varepsilon_{j}$
$\displaystyle\biggm{|}(X_{i^{\prime}},Z_{i^{\prime}})_{i^{\prime}=1}^{n},\widehat{f},\widehat{m}_{\widehat{f}}\Bigr{)}$
$\displaystyle=\frac{1}{n\nu^{2}}\sum_{i\neq
j}\widetilde{M}_{i}\widetilde{M}_{j}{\mathbb{E}}(\varepsilon_{i}\varepsilon_{j}\,|\,X_{i},X_{j},Z_{i},Z_{j})=0,$

and we deduce by Lemmas 14 and 15 that $b_{n}=o_{\mathcal{P}_{0}}(1)$. To see
that $c_{n}=o_{\mathcal{P}_{0}}(1)$, we proceed as above and write

$\displaystyle
c_{n}^{2}=\frac{1}{n\nu^{2}}\sum_{i=1}^{n}M_{i}^{2}\xi_{i}^{2}+\frac{1}{n\nu^{2}}\sum_{i\neq
j}M_{i}M_{j}\xi_{i}\xi_{j},$

where we again note that the first term is $o_{\mathcal{P}_{0}}(1)$ by
Assumption 3(c). Moreover,

$\displaystyle{\mathbb{E}}\Bigl{(}\frac{1}{n\nu^{2}}\sum_{i\neq
j}M_{i}M_{j}\xi_{i}\xi_{j}$
$\displaystyle\biggm{|}(Z_{i^{\prime}})_{i^{\prime}=1}^{n},\widehat{f}\Bigr{)}$
$\displaystyle=\frac{1}{n\nu^{2}}\sum_{i\neq
j}{\mathbb{E}}\bigl{\\{}{\mathbb{E}}\bigl{(}M_{i}M_{j}\,|\,(X_{i^{\prime}},Z_{i^{\prime}})_{i^{\prime}=1}^{n}\bigr{)}\xi_{i}\xi_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n},\widehat{f}\bigr{\\}}$
$\displaystyle=\frac{1}{n\nu^{2}}\sum_{i\neq
j}{\mathbb{E}}(R_{ij}\xi_{i}\xi_{j}\,|\,(Z_{i^{\prime}})_{i=1}^{n},\widehat{f}),$

where the last equality holds since

$\displaystyle{\mathbb{E}}\bigl{\\{}{\mathbb{E}}\bigl{(}M_{i}M_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n}\bigr{)}\xi_{i}\xi_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n},\widehat{f}\bigr{\\}}$
$\displaystyle={\mathbb{E}}\bigl{(}M_{i}M_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n}\bigr{)}{\mathbb{E}}(\xi_{i}\,|\,Z_{i},\widehat{f}){\mathbb{E}}(\xi_{j}\,|\,Z_{j},\widehat{f})$
$\displaystyle=0.$

Continuing, we have by (22) that

$\displaystyle\frac{1}{n\nu^{2}}\sum_{i\neq
j}{\mathbb{E}}\bigl{(}R_{ij}\xi_{i}\xi_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n},\widehat{f}\bigr{)}\leq$
$\displaystyle\frac{1}{cn\sigma^{2}}\sum_{i\neq
j}\big{|}{\mathbb{E}}(R_{ij}\xi_{n,i}\xi_{j}\,|\,(Z_{i^{\prime}})_{i^{\prime}=1}^{n},\widehat{f})\big{|}=o_{\mathcal{P}_{0}}(1).$

Therefore, by Lemmas 14 and 15 we conclude that $c_{n}=o_{\mathcal{P}_{0}}(1)$
as desired.

To deal with the final term, we define the triangular array
$W_{n,i}:=\nu^{-1}\varepsilon_{i}\xi_{i}$ for $n\in\mathbb{N}$ and $i\in[n]$,
and note that $W_{n,i}$ satisfies all the conditions of Lemma 16 by
Assumptions 3(a) and (d) (here we condition on $\widehat{f}$ in applying this
result). Hence, $U_{n}=n^{-1/2}\sum_{i=1}W_{n,i}$, and therefore
$T^{(\mathrm{N})}$, converges uniformly in distribution to $N(0,1)$.

We now show that $|(T^{(\mathrm{D})})^{2}-1|=o_{\mathcal{P}_{0}}(1)$, from
which the desired result follows from Lemma 15. Note that

$(T^{(\mathrm{D})})^{2}=\underbrace{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}L_{i}^{2}}_{p_{n}}-\biggl{(}\underbrace{\frac{1}{\sqrt{n}\nu}\sum_{i=1}^{n}L_{i}}_{q_{n}}\biggr{)}^{2}$
(25)

and that $q_{n}=\frac{1}{\sqrt{n}}T^{(\mathrm{N})}$. We have just shown that
$T^{(\mathrm{N})}=O_{\mathcal{P}_{0}}(1)$, so $q_{n}=o_{\mathcal{P}_{0}}(1)$
and we are therefore done if we can show that
$|p_{n}-1|=o_{\mathcal{P}_{0}}(1)$. Now

$\displaystyle p_{n}$
$\displaystyle=\underbrace{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\varepsilon_{i}^{2}\xi_{i}^{2}}_{\textup{I}_{n}}+\underbrace{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}M_{i}^{2}\widetilde{M}_{i}^{2}}_{\textup{II}_{n}}+\underbrace{\frac{4}{n\nu^{2}}\sum_{i=1}^{n}M_{i}\widetilde{M}_{i}\varepsilon_{i}\xi_{i}}_{\textup{III}_{n}}+\underbrace{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}^{2}\varepsilon_{i}^{2}}_{\textup{IV}_{n}^{\varepsilon}}$
(26) $\displaystyle\hskip
56.9055pt+\underbrace{\frac{1}{n\nu^{2}}\sum_{i=1}^{n}M_{i}^{2}\xi_{i}^{2}}_{\textup{IV}_{n}^{\xi}}+\underbrace{\frac{2}{n\nu^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}^{2}M_{i}\varepsilon_{i}}_{\textup{V}_{n}^{\varepsilon}}+\underbrace{\frac{2}{n\nu^{2}}\sum_{i=1}^{n}M_{i}^{2}\widetilde{M}_{i}\xi_{i}}_{\textup{V}_{n}^{\xi}}$
$\displaystyle\hskip
56.9055pt+\underbrace{\frac{2}{n\nu^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}\xi_{i}\varepsilon_{i}^{2}}_{\textup{VI}_{n}^{\varepsilon}}+\underbrace{\frac{2}{n\nu^{2}}\sum_{i=1}^{n}M_{i}\varepsilon_{i}\xi_{i}^{2}}_{\textup{VI}_{n}^{\xi}}.$

Consider the triangular array $\tilde{W}_{n,i}:=W_{n,i}^{2}$ for
$n\in\mathbb{N}$ and $i\in[n]$, and note that it satisfies all the conditions
of Lemma 17 by Assumptions 3(a) and (d) with $\mu_{n}=1$ (again conditioning
on $\widehat{f}$ in that result), so
$|\textup{I}_{n}-1|=o_{\mathcal{P}_{0}}(1)$. It remains to show that the
remaining terms are $o_{\mathcal{P}_{0}}(1)$. Now

$0\leq\textup{II}_{n}\leq\frac{1}{cn}\biggl{(}\sum_{i=1}^{n}M_{i}^{2}\biggr{)}\biggl{(}\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}^{2}\biggr{)}=o_{\mathcal{P}_{0}}(1)$

by Assumption 3(b). By the Cauchy–Schwarz inequality,

$|\textup{III}_{n}|\leq
4\biggl{(}\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\varepsilon_{i}^{2}\xi_{i}^{2}\biggr{)}^{1/2}\biggl{(}\frac{1}{n\nu^{2}}\sum_{i=1}^{n}M_{i}^{2}\widetilde{M}_{i}^{2}\biggr{)}^{1/2}=4\textup{I}_{n}^{1/2}\textup{II}_{n}^{1/2}=o_{\mathcal{P}_{0}}(1)$

by the above and Lemma 15. Now

$|\textup{IV}_{n}^{\varepsilon}|\leq\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}^{2}\varepsilon_{i}^{2}=o_{\mathcal{P}_{0}}(1)$

by Assumption 3(c). A similar argument shows that
$\textup{IV}_{n}^{\xi}=o_{\mathcal{P}_{0}}(1)$. By the triangle inequality and
the Cauchy–Schwarz inequality,

$|\textup{V}_{n}^{\varepsilon}|\leq
2\biggl{(}\frac{1}{n\nu^{2}}\sum_{i=1}^{n}M_{i}^{2}\widetilde{M}_{i}^{2}\biggr{)}^{1/2}\biggl{(}\frac{1}{n\nu^{2}}\sum_{i=1}^{n}\widetilde{M}_{i}^{2}\varepsilon_{i}^{2}\biggr{)}^{1/2}=2\textup{II}_{n}^{1/2}(\textup{IV}_{n}^{\varepsilon})^{1/2}=o_{\mathcal{P}_{0}}(1)$

by the above and Lemma 15. A similar argument can be used for
$\textup{V}_{n}^{\xi}$. Finally, again by the triangle inequality and the
Cauchy–Schwarz inequality,

$|\textup{VI}_{n}^{\varepsilon}|\leq