# FDP control in multivariate linear models using the bootstrap

Samuel Davenport, Bertrand Thirion, Pierre Neuvial

###### Abstract

In this article we develop a method for performing post hoc inference of the
False Discovery Proportion (FDP) over multiple contrasts of interest in the
multivariate linear model. To do so we use the bootstrap to simulate from the
distribution of the null contrasts. We combine the bootstrap with the post hoc
inference bounds of Blanchard et al. (2020) and prove that doing so provides
simultaneous asymptotic control of the FDP over all subsets of hypotheses.
This requires us to demonstrate consistency of the multivariate bootstrap in
the linear model, which we do via the Lindeberg Central Limit Theorem,
providing a simpler proof of this result than that of Eck (2018). We
demonstrate, via simulations, that our approach provides simultaneous control
of the FDP over all subsets and is typically more powerful than existing,
state of the art, parametric methods. We illustrate our approach on functional
Magnetic Resonance Imaging data from the Human Connectome project and on a
transcriptomic dataset of chronic obstructive pulmonary disease.

Keywords: FDP control, bootstrap, simultaneous inference, post hoc inference

##  1 Introduction

Statistical analysis of functional Magnetic Resonance Imaging data grounds the
inference of associations between external conditions (such as disease status
and experimental factors) and the signals recorded in brain regions, that is
assumed to reflect brain activity. In particular, practitioners typically aim
to uncover associations between local signals and conditions that are specific
to a given area; such specificity is essential for interpretation purposes.
The most standard framework is that of mass-univariate inference, in which
models are fit separately at each brain location, in order to detect
significant associations. This framework is simple and computationally
efficient but, given mm-scale resolution reached by current imaging setups,
results in a dire multiple comparison issue.

Statistical analysis of genomic data encounters a similar multiple comparison
problem. In particular, this is the case in Genome-Wide Association Studies
that aim to identify Single Nucleotide Polymorphisms that are associated with
one or more phenotypes of interest, and in gene expression studies, the goal
of which is to identify genes where activity is associated with one or more
variables of biological or clinical interest. In this field, the state-of-the-
art framework is also based on univariate tests that are performed for each
genomic marker. While imaging data typically consist of a smooth volume-domain
voxel grid, the dependence structure of genomic data is dictated by the
interdependence between genomic markers, which is mediated by haplotypic
blocks encountered in Genome-Wide Association Studies and by gene networks or
pathways in expression studies.

In both of these scientific fields (and in others), control of the false
discovery rate (FDR) has quickly become a de facto standard, as it yields
image or genome-level error control together with acceptable power (Genovese
et al., 2002; Storey and Tibshirani, 2003). In practice, most researchers
control the FDR using the Benjamini-Hochberg procedure (Benjamini and
Hochberg, 1995), under the assumption of positive regression dependence
(Benjamini and Yekutieli, 2001). This assumption is generally considered
reasonable given the positive correlation that typically exists between voxels
or genomic markers. However, users often interpret FDR-control as a control of
false discovery proportion (FDP), which is incorrect, as the FDR is only the
expected value of the FDP. Overall, this approach can result in unreliable
error control, especially when there is dependence within the data, see Korn
et al. (2004) and Figure 2.1 in Neuvial (2020). It is instead desirable to
provide probabilistic control on the _proportion_ or _number_ of false
discoveries.

Genomic and brain imaging datasets frequently involve the simultaneous test of
several contrasts (Smyth, 2004; Alberton et al., 2020). Such simultaneous
tests are important because they can ground double dissociation (Henson,
2006), ensuring the specificity of discoveries and leading to unambiguous
interpretations of the results. A difficulty arises here as the tests of the
different contrasts that are considered at each feature (voxel/gene) are
typically dependent and it may no longer be reasonable to assume positive
regression dependence. It is thus of interest to consider controlling the FDP
under the null hypothesis for each contrast, without making unwanted
assumptions.

The notion of _post hoc inference_ was introduced by Goeman and Solari (2011),
following earlier works by Genovese and Wasserman (2006); Meinshausen (2006)
on the probabilistic control of the FDP. The idea of post hoc inference is to
provide confidence bounds on the number or proportion of true/false
discoveries among arbitrary and possibly data-driven subsets of variables of
interest. By construction, such guarantees address the issue of circular
inference (Rosenblatt et al., 2018).

Post hoc bounds can be obtained as a by-product of the control of a multiple
testing risk called the joint error rate (JER) by a simple interpolation
argument (Blanchard et al., 2020). Using this construction, state-of-the-art
post hoc bounds (Goeman and Solari, 2011; Rosenblatt et al., 2018) can be
recovered from the Simes inequality, a classical result from the multiple
testing literature (Simes, 1986). The resulting bounds are valid under
positive regression dependence. They have also been shown to be conservative
in genomics and neuroimaging applications (Blanchard et al., 2021).

Since the joint error rate only depends on the joint distribution of the test
statistics under the null hypothesis, joint error rate control can
alternatively be obtained by randomization techniques (Blanchard et al., 2020,
2021; Hemerik et al., 2019). In particular, sharp data-driven joint error rate
control and associated post hoc bounds have been obtained for one-sample tests
using sign-flipping, and for two-sample tests using permutations (Blanchard et
al., 2020, 2021). However, obtaining joint error rate control more generally
in linear models, especially when doing inference on multiple contrasts,
remains an open question to the best of our knowledge.

In order to perform post hoc inference over contrasts we will need to be able
to obtain the joint null distribution of the test-statistics of multiple
contrasts within the framework of the linear model. To do so we use the
bootstrap, adjusting the approach of Westfall (2011) to the multivariate
setting. Justification for bootstrapping the residuals in a one-dimensional
linear model was first provided in Freedman (1981) based on theory proved in
Bickel and Freedman (1981). These results (and their proofs) were extended to
multivariate linear models in Eck (2018). In this work we provide an
alternative, simpler, proof of the validity of the bootstrap in the linear
model that relies on the Lindeberg Central Limit Theorem. We use these results
to show that we can obtain asymptotically valid simultaneous FDP control. This
extends the work of Blanchard et al. (2020) to the setting of the linear
model.

A python package with code to run the methods detailed in this paper is
available at: http://github.com/sjdavenport/pyperm, Jupyter notebooks with
illustrated examples of its use in practice are also available there.
Moreover, code to reproduce the analyses and figures of this paper is
available at: http://github.com/sjdavenport/lmfdp. Proofs and further
theoretical and simulation results are available in the Appendix.

##  2 Notation and general framework

###  2.1 Random Fields on a Lattice

Throughout we will take $(\Omega,\mathcal{F},\mathbb{P})$ to be a probability
space, write $\mathbb{E}$ to denote expectation and will define random
variables with respect to this space. We will also take $\mathbb{N}$ to be the
set of positive integers. We will primarily be working with random fields,
observed at a finite number of points, as our data. These are defined as
follows.

######  Definition 2.1.

Given $D,L\in\mathbb{N}$ and a finite set $\mathcal{V}\subset\mathbb{R}^{D}$,
we define a random field on $\mathcal{V}$ to be a measurable mapping
$f:\Omega\rightarrow\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\mathbb{R}^{L}}\right\\}$.
We will say that $f$ has dimension $L.$

Given $\omega\in\Omega$ and $v\in\mathcal{V}$ we will write
$f(\omega,v)=f(\omega)(v)$ and will typically drop dependence on $\omega$ and
simply refer to the random variable $f(v):\Omega\rightarrow\mathbb{R}^{L}$
when indexing $f$ and say that $f$ is a random field on $\mathcal{V}$. We
define the mean of $f$ to be the function
$\mu:\mathcal{V}\rightarrow\mathbb{R}^{L}$ sending $v\in\mathcal{V}$ to
$\mathbb{E}\mathopen{}\mathclose{{}\left[f(v)}\right]$. To each $f$ we
associate a covariance $\textswab{c}$ and a correlation function $\rho$ which
map $\mathcal{V}\times\mathcal{V}$ to $\mathbb{R}^{L\times L}$ and are defined
as

$\textswab{c}(u,v)=\text{cov}(f(u),f(v))=\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(f(u)-\mu(u)}\right)\mathopen{}\mathclose{{}\left(f(v)-\mu(v)}\right)^{T}}\right]$

and
$\rho(u,v)=\textswab{c}(u,v)\mathopen{}\mathclose{{}\left(\textswab{c}(u,u)\textswab{c}(v,v)}\right)^{-1/2}$
for all $u,v\in\mathcal{V}$.

For $1\leq j\leq L$, we define the random fields
$f_{j}:\Omega\rightarrow\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\mathbb{R}}\right\\}$
which send $\omega\in\Omega$ to
$f_{j}(\omega)(\cdot)=f(\omega)(\cdot)_{j}=f(\cdot)_{j}$. We will call
$f_{1},\dots,f_{L}$ the components of $f$ and will write the combination as
$f=[f_{1},\dots,f_{L}]^{T}.$ Convergence in distribution and probability
(which we will denote by $\overset{d}{\longrightarrow}$ and
$\overset{\mathbb{P}}{\longrightarrow}$) of random fields is well defined via
vectorization, see Section A.1 for a formalization of this. We also define
Gaussian random fields as follows.

######  Definition 2.2.

Given functions $\mu:\mathcal{V}\rightarrow\mathbb{R}^{L}$ and
$\textswab{c}:\mathcal{V}\times\mathcal{V}$ we write
$f\sim\mathcal{G}(\mu,\textswab{c})$ if $f$ is a random field with mean $\mu$
and covariance $\textswab{c}$ and such that the vector
$(f_{j}(v):v\in\mathcal{V},1\leq j\leq L)$ has a multivariate Gaussian
distribution.

###  2.2 Linear Model Framework

Let $\mathcal{V}\subset\mathbb{R}^{D}$ be a finite set of points corresponding
to the domain of interest (this could for instance be the voxels of the brain
or points representing genes). Suppose that we observe random fields
$y_{i}:\mathcal{V}\rightarrow\mathbb{R}$, for $1\leq i\leq n$ and some number
of subjects $n\in\mathbb{N}$. At each point $v\in\mathcal{V}$, we assume that

$Y_{n}(v)=X_{n}\beta(v)+E_{n}(v)$ (1)

where for each $v\in\mathcal{V}$, $Y_{n}(v)=[y_{1}(v),\dots,y_{n}(v)]^{T}$ is
a vector giving the observed data, $\beta(v)\in\mathbb{R}^{p}$ is the vector
of parameters (some $p\in\mathbb{N}$), $X_{n}\in\mathbb{R}^{n\times p}$ is the
design matrix of the covariates (note that this may include nuisance
variables) and $E_{n}=[\epsilon_{1},\dots,\epsilon_{n}]^{T}$ is an
$n$-dimensional random field on $\mathcal{V}$ which represents the unobserved
noise, where $(\epsilon_{n})_{n\in\mathbb{N}}$ is an i.i.d sequence of
1-dimensional random fields on $\mathcal{V}$ taking values in $\mathbb{R}$.
Note that we give the design matrix $X_{n}$ a subscript $n$ as we will allow
it to grow with $n.$

Then given contrasts $c_{1},\dots,c_{L}\in\mathbb{R}^{p}$ for some number of
contrasts $L\in\mathbb{N}$, we are interested in testing the null hypotheses:
$H_{0,l}(v):c_{l}^{T}\beta(v)=0,$ for $1\leq l\leq L$ and each
$v\in\mathcal{V}$. For each $v\in\mathcal{V}$ we can test the intersection
hypothesis

$H_{0}(v):c_{l}^{T}\beta(v)=0\text{ for }1\leq l\leq L$

using an $F$-test at each $v\in\mathcal{V}$ given by

$F_{n}(v)=\frac{(C\hat{\beta}_{n}(v))^{T}(C(X_{n}^{T}X_{n})^{-1}C^{T})^{-1}(C\hat{\beta}_{n}(v))/\text{rank}(C)}{\hat{\sigma}_{n}(v)^{2}}.$
(2)

Here $\hat{\beta}_{n}(v)=(X_{n}^{T}X_{n})^{-1}X_{n}^{T}Y_{n}(v)$ and
$C=(c_{1},\dots,c_{L})^{T}\in\mathbb{R}^{L\times p}$ is the matrix of
contrasts. $\hat{\sigma}_{n}^{2}:\mathcal{V}\rightarrow\mathbb{R}$ is the
estimate of the variance based on the residuals which sends $v\in\mathcal{V}$
to

$\hat{\sigma}_{n}^{2}(v)=\frac{1}{n-r_{n}}\mathopen{}\mathclose{{}\left\lVert
Y_{n}(v)-X_{n}\hat{\beta}_{n}(v)}\right\rVert^{2}.$

where $I_{n}$ is the $n\times n$ identity matrix and $r_{n}$ is the rank of
$X_{n}$. The individual null hypotheses can be tested using test statistics:

$T_{n,l}(v)=\frac{c_{l}^{T}\hat{\beta}_{n}(v)}{\sqrt{\hat{\sigma}_{n}(v)^{2}c_{l}^{T}(X_{n}^{T}X_{n})^{-1}c_{l}}}.$
(3)

Under $H_{0,l}(v)$ and assuming that the noise is Gaussian, conditional on
$X_{n}$, $T_{n,l}(v)$ is distributed as a $t$-statistic with $n-r_{n}$ degrees
of freedom. This allows a $p$-value to be calculated for each contrast $l$ at
each point $v$ as
$p_{n,l}(v)=2(1-\Phi_{n-r_{n}}(\mathopen{}\mathclose{{}\left|T_{n,l}(v)}\right|))$
where $\Phi_{d}$ is the CDF of a $t$-statistic with $d\in\mathbb{N}$ degrees
of freedom. Dropping the Gaussianity assumption, the $p$-values are still
asymptotically valid under reasonable assumptions (see e.g. Theorem B.4).
Moreover, for each $1\leq l\leq L,$ $T_{n,l}$ is a 1-dimensional random field
and we define $T_{n}=[T_{n,1},\dots,T_{n,L}]^{T}$.

###  2.3 Bounds on the False Discovery Proportion

The above framework gives us $m=LV$ different hypothesis tests (where $V$ is
the number of elements of $\mathcal{V}$), and results in a multiple testing
problem, which can be quite severe e.g. if the size of $\mathcal{V}$ is large.
Let $\mathcal{H}=\mathopen{}\mathclose{{}\left\\{(l,v):1\leq l\leq L\text{ and
}v\in\mathcal{V}}\right\\}$ index the hypotheses. For $H\subseteq\mathcal{H}$,
let $\mathopen{}\mathclose{{}\left|H}\right|$ denote the number of elements
within $H$. Finally let $\mathcal{N}\subset\mathcal{H}$ index the true null
hypotheses. Given $0<\alpha<1$ we will seek to provide a function
$V:\mathcal{H}\rightarrow\mathbb{N}$ such that

$\mathbb{P}\mathopen{}\mathclose{{}\left(|H\cap\mathcal{N}|\leq
V(H),\,\text{for all }H\subset\mathcal{H}}\right)\geq 1-\alpha.$ (4)

If (4) holds then simultaneously over all $H\subset\mathcal{H}$, with
probability $1-\alpha$, $V(H)$ provides an upper bound on the number of false
positives within $H.$ Suppose that for some $K\in\mathbb{N}$ we have sets
$R_{1},\dots,R_{K}\subset\mathcal{H}$ (which depend on the data) and constants
$\zeta_{1},\dots,\zeta_{K}\in\mathbb{N}$ and define

$\text{JER}\mathopen{}\mathclose{{}\left((R_{k},\zeta_{k})_{1\leq k\leq
K}}\right):=\mathbb{P}\mathopen{}\mathclose{{}\left(|R_{k}\cap\mathcal{N}|>\zeta_{k},\text{
some }1\leq k\leq K}\right)$ (5)

to be the joint error rate of the collection $(R_{k},\zeta_{k})_{1\leq k\leq
K}$. Blanchard et al. (2020) showed that if
$\text{JER}\mathopen{}\mathclose{{}\left((R_{k},\zeta_{k})_{1\leq k\leq
K}}\right)\leq\alpha$, then the bound
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V:\mathcal{H}\rightarrow\mathbb{R}$,
sending $H\subset\mathcal{H}$ to

$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V(H)=\min_{1\leq k\leq
K}(|H\setminus R_{k}|+\zeta_{k})\wedge|H|,$ (6)

satisfies (4) and thus provides an $\alpha$-level bound over the number of
false positives within each chosen rejection set. If the sets
$R_{1},\dots,R_{K}$ are nested then
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V$ is in fact optimal:
this follows by Blanchard et al. (2020)’s Proposition 2.5. We will follow the
approach of Blanchard et al. (2020) and define the collections
$(R_{k},\zeta_{k})_{1\leq k\leq K}$, that we will use, using template
families. In this case, on top of being statistically optimal, an important
practical feature of the bound
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V(H)$ is that can be
computed in linear time in $|H|$, see Algorithm 2 in Enjalbert-Courrech and
Neuvial (2022).

######  Definition 2.3.

Given $K\in\mathbb{N}$, we say that a family of functions $(t_{k})_{1\leq
k\leq K}$ is a template family if for each $1\leq k\leq K$,
$t_{k}:[0,1]\rightarrow\mathbb{R}$, $t_{k}(0)=0$ and $t_{k}$ is strictly
increasing and continuous. The parameter $K$ is called the size of the
template.

The simplest and most commonly used template family is the linear template
which, for $K\in\mathbb{N}$, is given by $t_{k}(x)=\frac{xk}{m}$ for $1\leq
k\leq K$ and $x\in[0,1]$. Existing post hoc bounds associated with this
template are briefly described in Section 4.3. However other choices are
available and the optimal choice of template may depend on the dataset under
consideration: we refer to Section 7 for further details and discussion of the
choice of template as well as Hemerik et al. (2019), Blanchard et al. (2020)
and Blain et al. (2022). Given a template family and $\lambda\in[0,1]$, for
each $1\leq k\leq K$ and $n\in\mathbb{N}$, we will take
$R_{k}(\lambda)=\mathopen{}\mathclose{{}\left\\{(l,v)\in\mathcal{H}:p_{n,l}(v)\leq
t_{k}(\lambda)}\right\\},$ set $\zeta_{k}=k-1,$ and let
$p^{n}_{(k:\mathcal{N})}$ be the $k$th smallest $p$-value in the set
$\mathopen{}\mathclose{{}\left\\{p_{n,l}(v):(l,v)\in\mathcal{N}}\right\\}$
(setting $p^{n}_{(k:\mathcal{N})}=1$ if $k>|\mathcal{N}|$). We will refer to
the collection $(R_{k}(\lambda),k-1)_{1\leq k\leq K}$ as the canonical
reference family.

######  Lemma 2.4.

For each $\lambda\in[0,1]$,

$\displaystyle\text{JER}\mathopen{}\mathclose{{}\left((R_{k}(\lambda),k-1)_{1\leq
k\leq K}}\right)=\mathbb{P}\mathopen{}\mathclose{{}\left(\min_{1\leq k\leq
K\wedge|\mathcal{H}|}t_{k}^{-1}(p^{n}_{(k:\mathcal{N})})\leq\lambda}\right).$

Thus for a given template family, in order to obtain an upper bound on the
number of false positives we can choose a threshold $\lambda\in[0,1]$ such
that

$\mathbb{P}\mathopen{}\mathclose{{}\left(\min_{1\leq k\leq
K\wedge|\mathcal{H}|}t_{k}^{-1}(p^{n}_{(k:\mathcal{N})})\leq\lambda}\right)\leq\alpha.$
(7)

Then the joint error rate of the family $(R_{k}(\lambda),k-1)_{1\leq k\leq K}$
is controlled to a level $\alpha$ and so the corresponding bound:
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V$, provides a
$(1-\alpha)$-level simultaneous upper bound on the number of false positives.

Blanchard et al. (2020) chose $\lambda$ via permutation testing, using the
fact that under an exchangeability assumption permutation allows the
probability in (7) to be controlled exactly. In the linear model, permutation
of the response does not satisfy the exchangeability assumption when there are
multiple potentially non-zero covariates in the model (see Appendix C.2 for a
discussion of this). In what follows we take a different approach that
proceeds via bootstrapping the data and results in asymptotic control of the
error rate.

For $\alpha\in(0,1)$
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V(H)$ provides an
$(1-\alpha)$-level simultaneous upper bound on the number of false positives
within $H$. From (4) we have

$\mathbb{P}\mathopen{}\mathclose{{}\left(\frac{|H\cap\mathcal{N}|}{|H|}\leq\frac{\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V(H)}{|H|},\,\,\forall
H\subset\mathcal{H}}\right)\geq 1-\alpha.$ (8)

It thus follows that for each $H\in\mathcal{H}$,
$\frac{\overline{V}(H)}{\mathopen{}\mathclose{{}\left|H}\right|}$ provides an
upper bound on the proportion of false positives within $H$ also known as the
false discovery proportion or FDP. Similarly $\frac{|H|-\overline{V}(H)}{|H|}$
provides a $(1-\alpha)$-level simultaneous lower bounds on the true discovery
proportion or TDP.

##  3 Bootstrapping in the Linear Model

###  3.1 Bootstrapping

There are several different ways to bootstrap data in the linear model (see
Freedman (1981)). We shall focus on the residual bootstrap. Given
$n\in\mathbb{N}$, this proceeds by calculating the residuals

$\hat{E}_{n}=Y_{n}-X_{n}\hat{\beta}_{n}=(I_{n}-X_{n}(X_{n}^{T}X_{n})^{-1}X_{n}^{T})E_{n},$
(9)

where $I_{n}$ is the $n\times n$ identity matrix and

$\hat{\beta}_{n}=(X_{n}^{T}X_{n})^{-1}X_{n}^{T}Y_{n}=\beta+(X_{n}^{T}X_{n})^{-1}X_{n}^{T}E_{n}.$
(10)

Given a number of bootstraps to perform: $B\in\mathbb{N}$ for each $1\leq
b\leq B$, conditional on the data, a selection:
$\hat{\epsilon}_{1}^{b},\dots,\hat{\epsilon}_{n}^{b}$ is chosen independently
with replacement from
$\mathopen{}\mathclose{{}\left\\{\hat{E}_{n,1},\dots,\hat{E}_{n,n}}\right\\}$
resulting in a combined random field
$E^{b}_{n}=[\hat{\epsilon}_{1}^{b},\dots,\hat{\epsilon}_{n}^{b}]^{T}.$ Given
this let $Y^{b}_{n}=X_{n}\hat{\beta}_{n}+E^{b}_{n}$ and define bootstrapped
parameter estimates
$\hat{\beta}_{n}^{b}=(X_{n}^{T}X_{n})^{-1}X_{n}^{T}Y^{b}_{n}.$

###  3.2 Convergence Results

Given this set-up we will now prove convergence of the bootstrapped test-
statistics. To do so we will require the following assumption.

######  Assumption 1.

zzzzz
zzzzza) For $n\in\mathbb{N}$, $X_{n}=[x_{1},\dots,x_{n}]^{T}$ for a sequence
of i.i.d vectors
$\mathopen{}\mathclose{{}\left(x_{n}}\right)_{n\in\mathbb{N}}$ in
$\mathbb{R}^{p}$ such that
$\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\lVert
x_{1}}\right\rVert^{5/2}}\right]<\infty$ and whose multivariate density is
bounded above.
zzzzzb) $(\epsilon_{n})_{n\in\mathbb{N}}$ is an i.i.d sequence of
1-dimensional random fields on $\mathcal{V}$ which is independent of
$\mathopen{}\mathclose{{}\left(x_{n}}\right)_{n\in\mathbb{N}}$ and such that
$\max_{v\in\mathcal{V}}\mathbb{E}\mathopen{}\mathclose{{}\left[\epsilon_{1}(v)^{4}}\right]<\infty$
and $\min_{v\in\mathcal{V}}\text{var}(\epsilon_{1}(v))>0.$

######  Theorem 3.1.

(Bootstrap convergence.) Suppose that $(X_{n})_{n\in\mathbb{N}}$ and
$(\epsilon_{n})_{n\in\mathbb{N}}$ satisfy Assumption 1. Then conditional on
$\mathopen{}\mathclose{{}\left(X_{m},Y_{m}}\right)_{m\in\mathbb{N}}$, for
almost all sequences
$\mathopen{}\mathclose{{}\left(X_{m},Y_{m}}\right)_{m\in\mathbb{N}}$, for each
$1\leq b\leq B,$ as $n\rightarrow\infty$,

$\sqrt{n}(\hat{\beta}_{n}^{b}-\hat{\beta}_{n})\overset{d}{\longrightarrow}\mathcal{G}(0,\textswab{c}_{\epsilon}\Sigma_{X}^{-1})$
$\text{ and
}\quad\hat{\sigma}^{b}_{n}\overset{\mathbb{P}}{\longrightarrow}\sigma.$

Using this result we obtain the following theorem.

######  Theorem 3.2.

(Bootstrap test-statistic convergence.) Suppose that
$(X_{n})_{n\in\mathbb{N}}$ and $(\epsilon_{n})_{n\in\mathbb{N}}$ satisfy
Assumption 1 and, for each $1\leq b\leq B,$ let
$T^{b}_{n}:\mathcal{V}\rightarrow\mathbb{R}$ be the $L$-dimensional random
field on $\mathcal{V}$ such that, for $1\leq l\leq L$,

$T_{n,l}^{b}=\frac{c_{l}^{T}(\hat{\beta}_{n}^{b}-\hat{\beta}_{n})}{\hat{\sigma}_{n}^{b}\sqrt{c_{l}^{T}(X_{n}^{T}X_{n})^{-1}c_{l}}}.$

Then conditional on
$\mathopen{}\mathclose{{}\left(X_{m},Y_{m}}\right)_{m\in\mathbb{N}}$, for
almost every sequence
$\mathopen{}\mathclose{{}\left(X_{m},Y_{m}}\right)_{m\in\mathbb{N}}$, for each
$1\leq b\leq B,$

$T_{n}^{b}\overset{d}{\longrightarrow}\mathcal{G}(0,\textswab{c}^{\prime})$

as $n\rightarrow\infty$. Here
$\textswab{c}^{\prime}:\mathcal{V}\times\mathcal{V}\rightarrow\mathbb{R}$
takes $u,v\in\mathcal{V}$ to
$\textswab{c}^{\prime}(u,v)=\rho_{\epsilon}(u,v)AC\Sigma_{X}^{-1}C^{T}A^{T}$
where $A\in\mathbb{R}^{L\times L}$ is a diagonal matrix with
$A_{ll}=(c_{l}^{T}\Sigma_{X}^{-1}c_{l})^{-1/2}$ for $1\leq l\leq L.$

Crucially the limiting distribution in this result is the same as the limiting
distribution of the test-statistics under the global null that $\beta=0$, see
Theorem B.4. It follows that the bootstrap provides consistent estimates of
the quantiles of functionals of the data under the global null, see Theorem
3.3. Our proof of Theorem 3.1 (see Section A) - which uses the Lindeberg
Central Limit Theorem - is substantially simpler than existing proofs that we
are aware of. Notably Eck (2018) proved a version of Theorem 3.1 by extending
the results of Freedman (1981) to multiple dimensions. Their approach, while
interesting, is rather complex and relies on the fact that convergence in
distribution is equivalent to convergence in the Mallows metric (Bickel and
Freedman, 1981).

###  3.3 Consistency of the bootstrap quantile

When bootstrapping we use the bootstrap samples to estimate quantiles of the
test-statistic under the null hypothesis. In what follows we will demonstrate
that as the number of bootstraps and subjects tends to infinity, the derived
quantiles converge to a limit. In order to do so given
$G:\mathbb{R}\rightarrow[0,1]$, define $G^{-}:[0,1]\rightarrow\mathbb{R}$,
which takes $y\in[0,1]$ to
$G^{-}(y)=\inf\mathopen{}\mathclose{{}\left\\{x:G(x)\geq y}\right\\}$, be the
generalized inverse of $G.$ Then we have the following theorem.

######  Theorem 3.3.

Let $(f_{n})_{n\in\mathbb{N}},f$ be functions from
$\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\mathbb{R}^{L}}\right\\}$
to $\mathbb{R}$ such that conditional on $(X_{m},Y_{m})_{m\in\mathbb{N}}$ for
almost all sequences $(X_{m},Y_{m})_{m\in\mathbb{N}}$, for each
$b\in\mathbb{N},$

$f_{n}(T^{b}_{n})\overset{d}{\longrightarrow}f(\mathcal{G}(0,\textswab{c}^{\prime})).$

For each $n,B\in\mathbb{N}$ and $0<\alpha<1$, let

$\lambda^{*}_{\alpha,n,B}=\inf\mathopen{}\mathclose{{}\left\\{\lambda:\frac{1}{B}\sum_{b=1}^{B}1\mathopen{}\mathclose{{}\left[f_{n}(T^{b}_{n})\leq\lambda}\right]\geq\alpha}\right\\}.$

Take $F$ to be the CDF of $f(\mathcal{G}(0,\textswab{c}^{\prime}))$, i.e. for
$\lambda\in[0,1]$,
$F(\lambda)=\mathbb{P}\mathopen{}\mathclose{{}\left(f(\mathcal{G}(0,\textswab{c}^{\prime}))\leq\lambda}\right)$
and assume that $F$ is strictly increasing and continuous. Then, letting
$\lambda_{\alpha}=F^{-1}(\alpha)$,

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\lambda^{*}_{\alpha,n,B}=\lambda_{\alpha}$

almost surely.

We refer to $\lambda^{*}_{\alpha,n,B}$ as the $\alpha$-quantile of the
bootstrap distribution of $f_{n}(T_{n})$ based on $B$ bootstraps. By Theorem
3.1, for suitable choices of $f_{n},f$,
$f_{n}(T_{n}^{b})\overset{d}{\longrightarrow}f(\mathcal{G}(0,\textswab{c}^{\prime}))$
as $n\rightarrow\infty$ and so Theorem 3.3 shows that the bootstrapped
$t$-statistics can be used to provide consistent estimates of the quantiles of
the limiting distribution of $f(\mathcal{G}(0,\textswab{c}^{\prime}))$. The
easiest example of such a suitable sequence of functions is to take $f$ to be
continuous and let $f_{n}=f$ for all $n\in\mathbb{N}$, because of the
Continuous Mapping Theorem. A more general result that provides a sufficient
condition, based on uniform convergence of $f_{n}$ to $f,$ is given in Lemma
C.3.

##  4 Joint error rate control in the linear model

In this section we will state and prove our main results. We will show that
given $0<\alpha<1$, choosing $\lambda$ to be the $\alpha$-quantile of the
bootstrapped distribution results in asymptotic $\alpha$ level control of the
joint error rate and thus results in simultaneous control of the FDP.

###  4.1 Joint error rate control

To set this up, given a test-statistic
$T:\mathcal{V}\rightarrow\mathbb{R}^{L}$, a subset $H\subset\mathcal{H}$, and
$n\in\mathbb{N}$, for $1\leq k\leq\mathopen{}\mathclose{{}\left|H}\right|$,
let $p^{n}_{(k:H)}(T)$ be the $k$th minimum value in the set

$\mathopen{}\mathclose{{}\left\\{2-2\Phi_{n-r_{n}}(\mathopen{}\mathclose{{}\left|T_{l}(v)}\right|):(l,v)\in
H}\right\\}.$

Using the results we have proved so far we can obtain the following theorem.

######  Theorem 4.1.

For $H\subset\mathcal{H}$, let
$f_{n,H}:\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\mathbb{R}^{L}}\right\\}\rightarrow\mathbb{R}$
send

$T\mapsto\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|H}\right|}t_{k}^{-1}(p^{n}_{(k:H)}(T))$

and for $n\in\mathbb{N}$ let $\lambda^{*}_{\alpha,n,B}(H)$ be the
$\alpha$-quantile of the bootstrap distribution (based on $B\in\mathbb{N}$
bootstraps) of $f_{n,H}(T_{n})$ conditional on the observed data. Assume that
the conditions of Assumption 1 hold and that $r_{n}=o(n)$. Then for
$\mathcal{N}\subset H$, conditional on
$\mathopen{}\mathclose{{}\left(X_{m},Y_{m}}\right)_{m\in\mathbb{N}}$,

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(f_{n,\mathcal{N}}(T_{n})\leq\lambda^{*}_{\alpha,n,B}(H)}\right)\leq\alpha.$

The limit holds with equality if $H=\mathcal{N}$. In particular taking
$H=\mathcal{H}$, it follows that

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(\min_{1\leq
k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|}t_{k}^{-1}(p^{n}_{(k:\mathcal{N})}(T_{n}))\leq\lambda^{*}_{\alpha,n,B}(\mathcal{H})}\right)\leq\alpha$

Applying this result and using Claim 2.4 we are thus able to obtain asymptotic
control of the joint error rate of the canonical reference family. Following
the discussion in Section 2.3 this means that we obtain asymptotic post hoc
FDP control. In particular we having the following corollary.

######  Corollary 4.2.

Under the assumptions of Theorem 4.1, for $0<\alpha<1,$ and
$H\subset\mathcal{H}$, let

$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V_{\alpha,n,B}(H)=\min_{1\leq
k\leq K}(|H\setminus
R_{k}(\lambda^{*}_{\alpha,n,B}(\mathcal{H}))|+k-1)\wedge|H|.$ $\text{Then
}\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(|H\cap\mathcal{N}|\leq\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V_{\alpha,n,B}(H),\,\,\forall
H\subset\mathcal{H}}\right)\geq 1-\alpha.$

Thus in order to provide FDP control, given a number of bootstraps
$B\in\mathbb{N}$, we can calculate $\lambda^{*}_{\alpha,n,B}(H)$, the
$\alpha$-quantile of the bootstrap distribution of $f_{n,H}(T_{n})$
conditional on the observed data. Then
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V_{\alpha,n,B}(H)$
provides a $(1-\alpha)$ level simultaneous upper bound on the number of false
positives in $H\subset\mathcal{H}$.

###  4.2 Bootstrap step-down procedure

It is possible to improve on the power of the above procedure by taking a
step-down approach in the spirit of (Romano and Wolf, 2005). This is based on
the idea that joint error rate control implies familywise error rate control,
see Section C.1. As such it is possible to obtain an estimate of the set of
null hypotheses and thereby obtain a tighter bound. The procedure can be
iterated as follows.

Algorithm 1 step-down bootstrap

1:Set $j\leftarrow 0$ and $H_{n}^{(0)}\leftarrow\mathcal{H}$

2:repeat

3: Set $j\leftarrow j+1$,
$\lambda_{n,j}\leftarrow\lambda^{*}_{\alpha,n,B}(H_{n}^{(j-1)})$ and
$H_{n}^{(j)}\leftarrow\mathopen{}\mathclose{{}\left\\{(l,v):p_{n,l}(v)\geq
t_{1}(\lambda_{n,j})}\right\\}$

4:until $H_{n}^{(j)}=H_{n}^{(j-1)}$

5:Set $\hat{H}_{n}\leftarrow H_{n}^{(j)}$ and return $\hat{H}_{n}$

As the follow theorem demonstrates, the step-down approach controls the joint
error rate and therefore provides simultaneous FDP control.

######  Theorem 4.3.

Under the assumptions of Theorem 4.1, for $0<\alpha<1,$ let $\hat{H}_{n}$ be
the set generated by applying Algorithm 1. Then

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(f_{n,\mathcal{N}}(T_{n})<\lambda^{*}_{\alpha,n,B}(\hat{H}_{n})}\right)\leq\alpha.$

Thus, for $H\subset\mathcal{H},$ letting
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V_{\alpha,n,B}(H)=|H|\wedge\min_{1\leq
k\leq K}(|H\setminus R_{k}(\lambda^{*}_{\alpha,n,B}(\hat{H}_{n}))|+k-1),$ it
follows that

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(|H\cap\mathcal{N}|\leq\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V_{\alpha,n,B}(H),\,\,\forall
H\subset\mathcal{H}}\right)\geq 1-\alpha.$

######  Remark 4.4.

The results in this subsection and the one previous have been stated for two-
sided $p$-values however they also hold for one-sided $p$-values without
change. All that is required to show this is to re-define $p^{n}_{(k:H)}(T)$
as the $k$th minimum value in the set

$\mathopen{}\mathclose{{}\left\\{1-\Phi_{n-r_{n}}(T_{l}(v)):(l,v)\in
H}\right\\}.$

In the definition of $\lambda^{*}_{\alpha,n,B}$ we require the computation of
$|\mathcal{H}|$ statistics for each bootstrap each of which is based on a
sample of size $n$. As such the complexity of these algorithms is
$O(nB|\mathcal{H}|)$.

###  4.3 Parametric Approaches

In this section we will discuss two parametric approaches to simultaneous FDP
inference which are based on the Simes inequality (12). Here we use the term
parametric to indicate that dependency assumptions on the data are required in
order for the methods to be valid. The first one is the original Simes post
hoc bound introduced in Goeman and Solari (2011). The second one is the method
of Rosenblatt et al. (2018) and Goeman et al. (2019). It corresponds to an
improvement on the basic Simes bound that is adaptive to the proportion of
true null hypotheses - i.e. it is a step-down version of the Simes bound. This
method has been applied to brain imaging data in Rosenblatt et al. (2018), and
is called ARI which stands for “All Resolutions Inference”. Both methods can
be conveniently formulated in terms of the bound $\overline{V}$ defined in
(6), associated to the linear template family $(t_{k})_{1\leq k\leq m}$, where
$t_{k}(\lambda)=\lambda k/m$, i.e.

$\displaystyle\overline{V}_{\lambda}(S)$ $\displaystyle=\min_{1\leq k\leq
m}\mathopen{}\mathclose{{}\left\\{\sum_{i\in
S}1\mathopen{}\mathclose{{}\left[p_{i}\geq\frac{\lambda
k}{m}}\right]+k-1}\right\\}\,.$ (11)

As noted by Blanchard et al. (2020), the Simes post hoc bound of Goeman and
Solari (2011) is simply
$\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V_{\alpha}$. Moreover,
letting $\overline{\alpha}=\alpha m/h(\alpha)$, where

$h(\alpha)=\max\mathopen{}\mathclose{{}\left\\{i\in\\{1,\dots,m\\},\forall
j\in\\{1,\dots,i\\},p_{(m-i+j)}>\frac{\alpha j}{i}}\right\\},$

the ARI bound of Goeman et al. (2019) is $\overline{V}_{\overline{\alpha}}$.
The quantity $h(\alpha)$ is called the Hommel factor (Hommel, 1988) and can be
interpreted as a $(1-\alpha)$-level upper confidence bound on $|\mathcal{N}|$,
the number of true null hypotheses.

If the null $p$-values satisfy positive regression dependence then both of
these methods result in simultaneous $(1-\alpha)$-level FDP control. This is
shown formally in Goeman and Solari (2011) and Goeman et al. (2019) via closed
testing and can also be shown to hold by combining the Simes inequality with
the joint error rate framework of Section 2.3. To see this note that if the
null $p$-values are positive regression dependent (Sarkar et al., 2008), then
the Simes inequality is satisfied, that is:

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left(\exists 1\leq
k\leq\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|:p^{n}_{(k:\mathcal{N})}\leq\frac{\alpha
k}{|\mathcal{N}|}}\right)\leq\alpha,$ (12)

with equality if the null $p$-values are independent.

In particular taking $\lambda=\alpha$, and noting that
$\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|\leq m$, the Simes
inequality implies that (7) holds (taking $K=m$ and $(t_{k})_{1\leq k\leq m}$
to be the linear reference family). Moreover, Goeman et al. (2019)’s Lemma 2
implies that if the null $p$-values satisfy positive regression dependence,
then

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left(\exists 1\leq
k\leq\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|:p^{n}_{(k:\mathcal{N})}\leq\frac{\alpha
k}{h(\alpha)}}\right)\leq\alpha.$ (13)

Thus taking $\lambda=\overline{\alpha}=\alpha m/h(\alpha)$, it follows that
(7) holds with respect to the linear reference family. In particular the Simes
procedure, which uses $\overline{V}_{\alpha}$ as a bound, and ARI, which uses
$\overline{V}_{\overline{\alpha}}$, provide simultaneous $(1-\alpha)$-level
control of the FDP.

In our results, presented in the following sections, we compare the
performance of the non-parametric bootstrap approach to these parametric
alternatives.

##  5 Simulation Results

###  5.1 Simulation Setup

In order to assess empirically that our method correctly controls the joint
error rate we run numerical simulations. We create the noise in these
simulations by generating 2-dimensional stationary Gaussian random fields on
domains which are $25$ by $25$, $50$ by $50$ and $100$ by $100$ pixels. To do
so we smooth Gaussian white noise with a Gaussian kernel with full width at
half maximum (FWHM) in $\mathopen{}\mathclose{{}\left\\{0,4,8}\right\\}$ (in
pixel units), accounting for edge effects to ensure stationarity (see e.g.
Davenport and Nichols (2022)), and scaled so that the variance is 1
everywhere.

We let the total number of subjects $n$ range from $20$ to $100$. For each
$n$, smoothness level, image size and
$\pi_{0}\in\mathopen{}\mathclose{{}\left\\{0.5,0.8,0.9,1}\right\\}$ , we run
5000 simulations - each with 100 bootstraps - to test the joint error rate.
For each simulation we do the following. First we generate $n$ Gaussian random
fields $\epsilon_{1},\dots,\epsilon_{n}$ as described above and add signal to
them (as detailed in the next paragraph). We then randomly divide these images
into 3 disjoint groups:
$G_{1},G_{2},G_{3}\subset\mathopen{}\mathclose{{}\left\\{1,\dots,n}\right\\}$
\- performing assignment to each group with equal probability (we eliminate
assignments where a given group has no entries). We test for the difference
between the first and the second group and between the second and the third
group - giving us two contrasts to differentiate between. We thus, in total,
test $5000$ hypotheses for the 50 by 50 scenario and $20000$ for the 100 by
100 case.

We vary the amount of signal in the datasets as follows. Given a proportion
$\pi_{0}$ we randomly choose a subset $\mathcal{N}$ of size
$\pi_{0}\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|$ of
$\mathcal{H}=\mathopen{}\mathclose{{}\left\\{(l,v):1\leq l\leq
2,v\in\mathcal{V}}\right\\}$ to be null (which is thus different in each
simulation) and add signal to ensure that the remainder are non-null. To do
so, for $1\leq i\leq n,$ and each $v\in\mathcal{V}$, we set

$Y_{i}(v)=1[i\in G_{2},(1,v)\not\in\mathcal{N}]+1[i\in
G_{3},(1,v)\not\in\mathcal{N}]+1[i\in
G_{3},(2,v)\not\in\mathcal{N}]+\epsilon_{i}(v).$

This ensures that the power of the test to detect a difference at $h$ is equal
for any $h\in\mathcal{N}^{C}$. If $\pi_{0}=1$ then all hypotheses are null. An
example realisation is shown in Figure 10.

In the next subsections we compare our bootstrap procedures, in terms of false
positive control and power, to two parametric alternatives: the Simes
procedure (Goeman and Solari, 2011) and its step-down: all resolutions
inference (ARI, Rosenblatt et al. (2018)) which are described in Section 4.3.

###  5.2 False Positive control

In each simulation setting, for $1\leq j\leq 5000$, we calculate a test
statistic random field $T_{n}^{(j)}$ and obtain $\lambda$ thresholds for the
single-step bootstrap, step-down bootstrap, Simes and ARI methods based on the
data as described in Section 4 , where we have used 100 bootstraps for the
non-parametric procedures. For each method we obtain $\lambda$-thresholds
$\lambda_{1},\dots,\lambda_{5000}$ allowing us to estimate the joint error
rate via the statistic

$\frac{1}{5000}\sum_{j=1}^{5000}1[f_{n,\mathcal{N}}(T_{n}^{(j)})\leq\lambda_{j}]$

which we refer to as the empirical joint error rate. Here $1[\cdot]$ denotes
the indicator function.

The results for the $50$ by $50$ simulations are displayed in Figure 1 and
those for the other domain sizes are shown in Figures 11 and 12. The results
for the bootstrap methods are shown in blue whilst those for the parametric
methods are shown in red. The solid lines indicate the step-down methods (i.e.
ARI and the step-down bootstrap). These plots demonstrate that, given a
reasonable number of subjects, the joint error rate of the bootstrap
procedures converges to the nominal level, in this case 0.1. Empirically the
parametric procedures are valid in all settings considered. However, their
control of the joint error rate is substantially below the nominal level i.e.
when the applied smoothing is non-zero, while the bootstrap approaches
demonstrate tighter control. The step-down procedures provide an improvement
on their single-step counterparts. This difference increases as $\pi_{0}$
decreases. See Section 5.3 for further details on the effect of $\pi_{0}$.

Figure 1: Comparing the empirical joint error rate across methods for the
simulation setting described in Section 5.1 for $\alpha=0.1$ in the $50$ by
$50$ simulation setting. The bootstrap procedures typically provide tighter
control of the joint error rate than the parametric ones, except under
independence. The bootstrap methods are shown in blue whilst the parametric
methods are shown in red. The solid lines indicate the step-down methods. The
thin flat black dashed lines provide 95% marginal confidence bands based on
the normal approximation to the binomial distribution.

###  5.3 Power

In this section we compare the power of the various methods in the simulation
setting described in Section 5.1 in the case where the applied FWHM is 4
pixels. We have chosen to focus on this level of smoothness because it
represents a realistic level of applied smoothness and illustrates the
benefits that can be achieved when using the bootstrap under dependence.

Here we shall use a notion of power originally proposed in Blanchard et al.
(2020) to compare the ability of joint error rate controlling procedures to
detect signal. Given a set $R\subset\mathcal{H}$, define

$\text{Pow}(R):=\mathbb{E}\mathopen{}\mathclose{{}\left[\frac{|R|-\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V(R)}{\mathopen{}\mathclose{{}\left|R\cap(\mathcal{H}\setminus\mathcal{N})}\right|}\middle|\mathopen{}\mathclose{{}\left|R\cap(\mathcal{H}\setminus\mathcal{N})}\right|>0}\right]$
(14)

where $\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V$ is defined as in
equation 6. Here we consider the following choices of $R$ with which we
compare the power (as in Blanchard et al. (2020)). 1) $R=\mathcal{H}$ and 2)
Taking $R$ to be the hypotheses of $\mathcal{H}$ which are rejected by the
Benjamini Hochberg procedure, applied to the $p$-values
$\mathopen{}\mathclose{{}\left\\{p_{n,l}(v):(l,v)\in\mathcal{H}}\right\\}$, at
a level $0.05$. Note that, unlike in Blanchard et al. (2020), no additional
level of randomness in the choice of the sets in 2) is prescribed. We also
consider taking $R=\mathopen{}\mathclose{{}\left\\{(l,v):p_{n,l}(v)\leq
0.05}\right\\}$, see Section E.5, the results for which are similar in nature
to scenario 1 from above. The results for cases 1) and 2) are illustrated
graphically in Figure 13. These are for simulations on the $50$ by $50$
domain.

Figure 2: Plotting the power of the different methods against the number of
subjects. The power for setting 1 (i.e. $R$ = $\mathcal{H}$) is shown in the
top row and the power for setting 2 (i.e. taking $R$ to be the Benjamini-
Hochberg rejection set) is shown in the bottom row.

From these plots we can see that overall the bootstrap based approaches have a
higher power than the parametric ones. The power of ARI only becomes
comparable (or higher) to that of the bootstrap in the extreme scenario
($\pi_{0}=0.5$) given a large enough sample size. Additionally the bootstrap
is not robust at the smallest sample size considered (i.e. $n=10$) where it is
slightly conservative. However it is important to note that in typical high-
dimensional applications (neuroimaging, genetics) $\pi_{0}>0.9$ and $n$ is
often substantially greater than $20.$

The lower the value of $\pi_{0}$, the greater the increase in power that is
obtained by using the step-down algorithms. ARI is always more powerful than
Simes by construction. In the relatively sparse scenarios (i.e. $\pi_{0}\geq
0.8$) they have a very similar power however for $\pi_{0}=0.5$, ARI provides a
marked improvement over Simes. The bootstrap step-down always improves on the
standard bootstrap approach though the difference is not particularly large:
even when $\pi=0.5$ this increase is relatively small. The similarity of the
standard and step-down procedures, for both the parametric and bootstrap
methods, is consistent with the results obtained on real data which are
described in the next subsections.

##  6 Real Data results

###  6.1 Neuroimaging data application

We have 3D functional Magnetic Resonance Imaging data from $n=386$ unrelated
subjects, who performed an $m$-back working memory task, from the Human
Connectome Project. After pre-processing (described in Section D) we obtain a
3-dimensional contrast image for each subject. We fit a linear model to these
images including sex, height, weight, body mass index, two different measures
of blood pressure, handedness and IQ (measured using the PMAT24_A_CR test
score). We consider sex and IQ as a variables of interest. We obtain test-
statistic contrasts for sex and IQ and a $p$-value at each voxel for each
contrast. We form clusters using a cluster defining threshold on the
$p$-values of $p=0.001$, with each cluster being a contiguous set of voxels
above the threshold (clusters are defined separately for each contrast of
interest).

We use our bootstrap framework to provide a lower bound on the proportion of
active voxels within each cluster. This illustrates that multiple clusters, in
different regions of the brain, have a relatively large proportion of active
voxels for the contrast of IQ. For the contrast of sex only a single cluster
has a non-zero lower bound on the number of true positives. The bounds
provided using the step-down bootstrap procedure are the same as the single-
step version in this example.

We compare to the results that are obtained using the parametric methods of
Goeman and Solari (2011) and Rosenblatt et al. (2018) and see that our
bootstrap approach results in higher lower bounds on the number of active
voxels. The bounds obtained using the single-step and step-down parametric
methods are very similar, which is not surprising given the sparsity of the
signal. For the IQ contrast the lower bounds provided by the bootstrap and ARI
for the number of true positives and on the TDP within each cluster are shown
graphically in the upper panel of Figure 8. The corresponding plot for the sex
contrast is shown in Section E.2. Direct comparison of the lower bounds is
shown in Figure 4.

Figure 3: TDP bounds within clusters for the contrast for IQ in the linear
regression model fit to the HCP data. Each cluster is shaded a single colour
which is the lower bound on the TDP. The upper panel gives the TDP bounds
within each cluster provided by the bootstrap procedure. The lower panel gives
the bounds provided by using ARI. The bounds given by the bootstrap are larger
(as indicated by the lighter colours) indicating that the method is more
powerful. Note that these images are 2D slices through the 3D brain and so
voxels that are part of the same cluster are not necessarily connected.

Figure 4: Comparing the TDP and true positive lower bounds across clusters for
the different methods. The bootstrap lower bounds are consistently higher than
the parametric methods. Clusters are organized from left to right in terms of
their size. Only one cluster for the sex contrast is found: this is the 2nd
smallest cluster overall with a TP lower bound of 1 voxel. The sizes and
bounds of the clusters in the IQ contrast are larger. For the largest cluster
we are able to conclude that it contains 908 true positives using the
bootstrap approach.

###  6.2 Transcriptomic data application

In this section, we illustrate the application of our methods to a specific
gene expression data set. Gene expression studies use microarray or sequencing
biotechologies in order to measure the activity (or “expression level”) of a
large number of genes simultaneously. We focus on a study of chronic
obstructive pulmonary disease (COPD), Bahr et al. (2013), whose main goal was
to identify genes whose expression level is significantly associated with lung
function. In order to do this, the authors fit a linear model for this
association for each gene, while controlling for the following covariates:
age, sex, body mass index, parental history of COPD, and two smoking variables
(smoking status and pack-years). The number of subjects is $n=135$ while the
number of genes is $V=12,531$, leading to a large-scale multiple testing
problem. Using the Benjamini-Hochberg method to control the FDR at the $5\%$
level, $1,745$ genes were found to be significantly associated.

We fit this linear model to the data, regressing the gene level data against
the controlled covariates and lung function and considering a single contrast
for lung function. We performed 1000 bootstraps and used these to obtain
$\lambda^{*}_{\alpha,135,1000}=0.22$, where we took $\alpha=0.1.$ This allows
us to provide a $(1-\alpha)$-level simultaneous lower bounds on the number of
true positives within any specified set of genes via (6). In particular it
allows us to conclude (with $90\%$ confidence) that at least 1,354 of the
1,745 genes within the Benjamini-Hochberg significance set are active. The
stepdown bootstrap provides the same bound as the single step version in this
case. Simes and ARI provide lower bounds on the number of true positives in
this set of 917 and 966 respectively, which are substantially less informative
than the bootstrap bounds.

In the absence of prior information on genes, a natural idea is to rank them
by decreasing statistical significance. Our post hoc methods provide upper
confidence curves on the proportion of true positives among the most
significant genes. Such curves are displayed in Figure 5, where the blue lines
correspond to our proposed single step and step-down bootstrap-based methods,
and the red lines correspond to the parametric approaches of Goeman and Solari
(2011) and Rosenblatt et al. (2018). These results are consistent with the
numerical experiments of Section 5. First, the bootstrap method yields post
hoc bounds that are substantially more informative that their parametric
counterpart. Second, the difference between single-step methods and their
step-down or $\pi_{0}$-adaptive counterpart is very small, which is consistent
with the fact that the signal is expected to be small in such genomic data
sets, corresponding to $\pi_{0}$ close to $1$. For the bootstrap there is in
fact no difference between the single step and step-down approach in this
example.

A widely used approach in differential expression studies is to select genes
based on the conjunction of a threshold on the $p$-values and a threshold on
its effect size (Cui and Churchill, 2003). Ebrahimpoor et al. (2020) recently
noted that this type of double selection can lead to inflated numbers of false
discoveries when used in conjunction with FDR-based multiple testing
corrections, whereas post hoc inference is by construction robust against this
issue. The use of our proposed post hoc bounds in this context is illustrated
in the volcano plot in Figure 6 (Cui and Churchill, 2003). In this plot, each
gene is represented in two dimensions by estimates of its effect size ($x$
axis, also kown as “fold change” in genomics) and significance ($y$ axis), in
a logarithmic scale. Figure 6 illustrates a particular selection,
corresponding to the genes whose $p$-value is below $0.001$ and whose effect
size is above $0.5$. Our bootstrap-based bound ensures that with probability
$1-\alpha=90\%$, among these $546$ genes, at least $490$ are true positives,
corresponding to a FDP below $0.1$. Importantly, the $p$-value and effect size
thresholds can be chosen post hoc, and multiple such choices can be made
without compromising the statistical coverage of the associated bound. For
example, the bounds associated to the gene subsets with positive and negative
effect size are also displayed in Figure 6.

Figure 5: False discovery proportion and true positive plots for the
transcriptomic dataset. In the upper panels, for $k=1,\dots,3000$, upper
bounds on the FDP and lower bounds on the number of true positives are
provided by each of the methods for the sets comprised of the hypotheses with
the $k$ smallest $p$-values. The silver vertical line corresponds to the
location of the Benjamini-Hochberg rejection set. The lower panels provide a
zoomed in version of the same plot for for the $1000$ smallest $p$-values. The
bootstrap methods provide substantially better bounds than the parametric
ones. ARI slightly improves on Simes while the step-down bootstrap is
indistinguishable from the single step bootstrap approach in this setting.
Figure 6: A volcano plot for the $p$-values for the transcriptomic data. For
each gene this plots the estimated contrast effect size (labelled as fold
change and corresponding to $c^{T}\hat{\beta}_{135}$ where $c$ is the contrast
vector for COPD) against the $p$-value, where both are measured in log scale.
Two regions (shown via shading) are selected containing the genes whose
$p$-values are less than $10^{-3}$ and for which the absolute fold change is
greater than $10^{0.5}$. Bounds on the true positives (TP) and FDP, overall
and for the shaded regions are provided on the plot.

##  7 Discussion

In this paper we have introduced a bootstrap method which provides
simultaneous control of the FDP over subsets of hypotheses of multiple
contrasts in the linear model. We have proved the asymptotic validity of this
approach and shown, via simulation, that the error rate is controlled to the
correct level given a reasonable number of subjects.

From our simulations and real data examples, we can see that the bootstrap
approach typically provides better bounds than existing, state of the art
parametric methods (i.e. Simes and ARI). This occurs because we are able to
model the dependence within the data. The parametric methods, on the other
hand, rely on the Simes inequality which is only exact under independence.
Moreover the Simes inequality is only valid under positive regression
dependence whereas the non-parametric bootstrap makes relatively few
assumptions other than finite moments of the noise and the design. Moreover in
real data situations there is typically relatively strong dependence within
the data and so we would expect the bootstrap to give better bounds. This is
illustrated in our brain imaging and transcriptomic examples where the
bootstrap bounds provided substantial improvements over the ones derived using
the parametric methods. Overall, these results are consistent with those
obtained by Enjalbert-Courrech and Neuvial (2022) in the specific case of two-
sample tests, where post hoc bounds based on non-parametric joint error rate
control substantially outperformed their parametric counterpart.

The step-down bootstrap approach improves the power whilst maintaining control
of the error rate. However in practice, as illustrated in the real datasets,
the improvement is likely to be small as $\pi_{0}$ will be close to 1. Indeed
in both of our real data examples there is no noticeable improvement. The
improvement of ARI over Simes is typically non-zero but is rather small. These
results demonstrate that the step-down methods, whether parametric or
otherwise, appear to require a relatively small value of $\pi_{0}$ before they
substantial improvement on their single step versions. It is worth noting that
the improvement in the bounds provided by ARI relative to Simes is greater
than that of the step-down bootstrap relative to the single-step version. One
possible reason for this discrepancy is that in the step-down bootstrap
(Algorithm 1), only the first threshold $t_{1}$ is used at each step. This
implies that only part of the information on the true positives is exploited.

It is important to note that for the bootstrap approach, control of the FDP is
asymptotic. In Section 5.2 we showed that given reasonable sample sizes and
smoothness levels (e.g. $n\geq 20$ for FWHM = $4,8$) the joint error rate was
controlled at the correct rate. At low smoothness levels and low sample sizes
the error rate can be slightly inflated. In this scenario this inflation is
counter balanced by a small amount of signal. Moreover at very small sample
sizes e.g. $n=10$ the bootstrap can be conservative (see Figure 13). At the
smoothness levels and sample sizes used in real data analyses, based on our
theoretical and simulation results, we would expect the bootstrap to control
the joint error rate to the desired level.

The template $(t_{k})_{1\leq k\leq K}$ is a free parameter of the proposed
method, and optimising this choice for a particular application could lead to
tighter TDP bounds. For the numerical experiments of this paper, we have only
considered the linear template, for which $t_{k}(\lambda)=\frac{\lambda
k}{m}$, which is the most widely used in the post hoc inference literature
(Goeman et al., 2019) and in particular for neuroimaging applications
(Rosenblatt et al., 2018). Other parametric templates are considered in
Blanchard et al. (2020); Andreella et al. (2020). However, the experiments
reported in Andreella et al. (2020) suggest that the linear template may be
difficult to beat. A natural idea to go beyond parametric templates is to
learn from the data the shape of the template itself. This has been advocated
by Meinshausen (2006), but the proof of the proposed method is invalid as it
suffers from a circularity issue (Blanchard et al., 2020, Remark 5.3).
Recently, Blain et al. (2022) used an independent data set to learn the
optimal template in the case of one- and two-sample testing. Extending this
idea to multivariate linear models as considered in the present paper is an
interesting perspective for future research.

Our results can also be used to provide strong control of the familywise error
rate over multiple contrasts (see Appendix C.1 for a proof, a formal
discussion of this and the results of simulations). This comes about in two
different ways. Firstly it arises as a direct consequence of joint error rate
control when using the canonical reference family and taking $\zeta_{k}=k-1$.
Secondly the familywise error rate can be targeted directly, along the lines
of the approach of Westfall (2011) (i.e. not simultaneously with joint error
rate control), this follows from Theorem 4.1 by taking $K=1$, a result that is
stated formally in Theorem C.1. Recently Alberton et al. (2020) sought to
provide direct familywise error rate control over multiple contrasts. They
resampled their data using Manly permutation (personal correspondence with the
authors) which provides weak (rather than strong) control when there are
multiple covariates in the model which may or may not be non-zero. This occurs
because Manly permutation acts by permuting the $Y$s and thus does not
generate resamples under the full null hypothesis - see our Appendix C.2 for
details. The bootstrap is able to avoid these issues, and thus provide strong
control, because it centres the residuals before resampling. As discussed
below, other forms of permutation testing could be used to provide the desired
error rate control as an alternative to the bootstrap.

An alternative approach to controlling the joint error rate over multiple
contrasts could be developed by considering permutations of the residuals
rather than bootstrapping. There are number of methods which could be used to
permute the data (see Andersen et al. (1999) for a comparison). Given the
widespread using of permutation testing based methods in the linear model
(Winkler et al., 2014) this is an interesting avenue for future research.
Importantly, like the bootstrap, permuting the residuals via these methods
will typically be valid asymptotically. This is because, when dealing with
multiple contrasts in the linear model, exchangeability does not hold and so
permutation is not exact: see Section C.2 for further details.

The choice of the error rate to control in a scenario where many hypotheses
are being tested depends strongly on the goals of the researcher. The bounds
that we have provided on the FDP provide more informative inference than
simply controlling the FDR. As discussed in Neuvial (2020), under dependence
controlling the FDR can lead to non-nonsensical results. Instead bounds on the
FDP allow statements about the number of active voxels with a given set to be
made. Moreover this inference is valid simultaneously over all sets and so
guards against circular inference.

## Acknowledgments

We are grateful to Alexandre Blain at the University of Toulouse for his help
with demonstrating how to use the sanssouci code to get ARI to work and for
checking that the python implementation was consistent with the implementation
in the ARIbrain R package. We are also grateful to François Bachoc at the
University of Toulouse for his help with the proof of Lemma C.4. SD is also
grateful to Fabian Telschow at Humboldt University for useful discussions on
bootstrapping.

SD was supported by NIH grant R01EB026859. SD and PN were supported by the
SansSouci ANR project (ANR-16-CE40-0019). BT was supported by the KARAIB AI
chair (ANR-20-CHIA-0025-01) and the FastBig ANR project (ANR-17-CE23-0011).

Data were provided in part by the Human Connectome Project, WU-Minn Consortium
(Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657)
funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for
Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at
Washington University.

##  Appendix A  Consistency of the bootstrap in the linear model

###  A.1 Further theory for random fields

Given two random fields $f,g:\mathcal{V}\rightarrow\mathbb{R}^{L}$ operations
of addition and subtraction can be performed pointwise and so $f+g$ and $f-g$
are well defined. Moreover if instead $g:\mathcal{V}\rightarrow\mathbb{R}$
then multiplication and division can also be performed pointwise and so, in
that case, $fg$ and $f/g$ are well-defined.

Given $D\in\mathbb{N}$, suppose that
$\mathcal{V}=\mathopen{}\mathclose{{}\left\\{u_{1},\dots,u_{V}}\right\\}$ for
some $V\in\mathbb{N}$ and $u_{1},\dots,u_{V}\in\mathbb{R}^{D}$. For
$L\in\mathbb{N}$, let $f:\mathcal{V}\rightarrow\mathbb{R}^{L}$ be a random
field. Then we define vec$(f)\in\mathbb{R}^{LV}$ to be the vector whose
$((i-1)L+j)$th element is $f_{j}(u_{i})$ for $1\leq i\leq V$ and $1\leq j\leq
L.$ We refer this operation as vectorization. This allows us to easily define
notions of convergence. Given a sequence: $((f_{n})_{n\in\mathbb{N}},f)$ of
random fields from $\mathcal{V}$ to $\mathbb{R}^{L}$ we say that $f_{n}$
converges to $f$ in distribution (resp. probability/almost surely) if
vec$(f_{n})$ converges in distribution (resp. probability/almost surely) to
vec$(f)$. We will write this as $f_{n}\overset{d}{\longrightarrow}f$ (resp.
$f_{n}\overset{\mathbb{P}}{\longrightarrow}f/f_{n}\overset{a.s.}{\longrightarrow}f$).
Given such a sequence we will write $f_{n,j}\,\,(1\leq j\leq L)$ to denote its
components.

######  Definition A.1.

For $L,L^{\prime}\in\mathbb{N}$ let $f:\mathcal{V}\rightarrow\mathbb{R}^{L}$
be a random field then we define the random field $Mf$ which sends
$v\in\mathcal{V}$ to $Mf(v)\in\mathbb{R}^{L^{\prime}}$.

######  Lemma A.2.

For $L,L^{\prime}\in\mathbb{N}$ let $f:\mathcal{V}\rightarrow\mathbb{R}^{L}$
be a random field with covariance $\textswab{c}$ and let
$M\in\mathbb{R}^{L^{\prime}\times L}$, then $Mf$ has covariance

$M\,\textswab{c}\,M^{T}.$

Moreover if $f$ is Gaussian then so is $Mf.$

###  A.2 Lindeberg Central Limit Theorem

In order to proof our main results we require Proposition A.4 (stated below)
which we prove using the Lindeberg CLT (see e.g. Van der Vaart (2000) Chapter
2.8). We will also require the following lemma.

######  Lemma A.3.

Let $X$ and $Y$ be random variables such that
$\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|X}\right|^{2+\eta}}\right]<\infty$
and
$\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|Y}\right|^{K}}\right]<\infty$
for some $K,\eta>0$, then for all $a\in\mathbb{R},$

$\mathbb{E}\mathopen{}\mathclose{{}\left[X^{2}1[a|Y|>\gamma]}\right]\leq\gamma^{-K/q}a^{K/q}\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|X}\right|^{2+\eta}}\right]^{1/(1+\eta/2)}\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|Y}\right|^{K}}\right]^{1/q}$

where $q=1-(1+\eta/2)^{-1}.$

###### Proof.

By Holder’s inequality for $p,q\in\mathbb{R}_{>0}$ such that
$\frac{1}{p}+\frac{1}{q}=1,$

$\displaystyle\mathbb{E}\mathopen{}\mathclose{{}\left[X^{2}1[a|Y|>\gamma]}\right]$
$\displaystyle\leq\mathbb{E}\mathopen{}\mathclose{{}\left[X^{2p}}\right]^{1/p}\mathbb{E}\mathopen{}\mathclose{{}\left[1[a|Y|>\gamma]}\right]^{1/q}=\mathbb{E}\mathopen{}\mathclose{{}\left[X^{2p}}\right]^{1/p}\mathbb{P}\mathopen{}\mathclose{{}\left(a|Y|>\gamma}\right)^{1/q}$
$\displaystyle\leq\mathbb{E}\mathopen{}\mathclose{{}\left[X^{2p}}\right]^{1/p}\mathopen{}\mathclose{{}\left(\frac{\mathbb{E}\mathopen{}\mathclose{{}\left[a^{K}\mathopen{}\mathclose{{}\left|Y}\right|^{K}}\right]}{\gamma^{K}}}\right)^{1/q}=\gamma^{-K/q}a^{K/q}\mathbb{E}\mathopen{}\mathclose{{}\left[X^{2p}}\right]^{1/p}\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|Y}\right|^{K}}\right]^{1/q}$

where the middle inequality holds by Markov’s inequality. Taking $p=1+\eta/2$
and $q=1-\frac{1}{p}$, the result follows. ∎

######  Proposition A.4.

Given a sequence $(k_{n})_{n\in\mathbb{N}}$, let
$\mathopen{}\mathclose{{}\left\\{\xi_{n,i}:n,i\in\mathbb{N},1\leq i\leq
k_{n}}\right\\}$ be a triangular array of mean-zero random fields on
$\mathcal{V}$ which are i.i.d within rows and have finite covariance. Let
$\mathopen{}\mathclose{{}\left\\{a_{ni}:n,i\in\mathbb{N},1\leq i\leq
n}\right\\}$ be a triangular array of $D-$dimensional vectors such that
$\sum_{i=1}^{n}\mathopen{}\mathclose{{}\left\lVert
a_{ni}}\right\rVert^{2+K/q}\rightarrow 0$ as $n\rightarrow\infty$ and
$\sup_{i,n}\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|\xi_{n,i}}\right|^{\max(K,2+\eta)}}\right]<\infty$
for some $K>0$, any $\eta>0$ and $q=1-(1+\eta/2)^{-1}$. Let
$A_{n}=(a_{n1},\dots,a_{nk_{n}})\in\mathbb{R}^{D\times k_{n}}$ and suppose
that $A_{n}^{T}A_{n}\rightarrow\Sigma\in\mathbb{R}^{D\times D}.$ For
$n\in\mathbb{N}$, let $\textswab{c}_{n}$ be the covariance function of
$\xi_{n,1}$ and suppose that as
$n\rightarrow\infty,\textswab{c}_{n}\rightarrow\textswab{c}$ (pointwise) for
some covariance function $\textswab{c}$ on $\mathcal{V}$. Then as
$n\rightarrow\infty,$

$\sum_{i=1}^{k_{n}}a_{ni}\xi_{n,i}\overset{d}{\longrightarrow}\mathcal{G}(0,\textswab{c}\Sigma).$

###### Proof.

The proof is an application of the Lindeberg CLT (see e.g. van der Vaart
(1998) Proposition 2.27) to the vectors $\text{vec}(a_{ni}\xi_{n,i})$. There
are two conditions to verify. The first is to show that the covariance
converges. We can show this blockwise, i.e., for each $u,v\in\mathcal{V}$,

$\displaystyle\sum_{i=1}^{k_{n}}\text{cov}(a_{ni}\xi_{n,i}(u),a_{ni}\xi_{n,i}(v))$
$\displaystyle=\sum_{i=1}^{k_{n}}\mathbb{E}\mathopen{}\mathclose{{}\left[a_{ni}\xi_{n,i}(u)\xi_{n,i}(v)a_{ni}^{T}}\right]$
$\displaystyle=\textswab{c}_{n}(u,v)\sum_{i=1}^{k_{n}}a_{ni}a_{ni}^{T}=\textswab{c}_{n}(u,v)A_{n}^{T}A_{n}.$

which converges to $\textswab{c}(u,v)\Sigma$ as $n\rightarrow\infty.$ For the
second condition we need to show that for all $\gamma>0,$

$\sum_{i=1}^{k_{n}}\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\lVert\text{vec}(a_{ni}\xi_{n,i})}\right\rVert^{2}1[\mathopen{}\mathclose{{}\left\lVert\text{vec}(a_{ni}\xi_{n,i})}\right\rVert>\gamma]}\right]\underset{n\rightarrow\infty}{\longrightarrow}0.$

We can expand the left hand side as

$\displaystyle\sum_{i=1}^{k_{n}}$
$\displaystyle\mathbb{E}\mathopen{}\mathclose{{}\left[\sum_{v\in\mathcal{V}}\mathopen{}\mathclose{{}\left\lVert
a_{ni}\xi_{n,i}(v)}\right\rVert^{2}1\mathopen{}\mathclose{{}\left[\sum_{u\in\mathcal{V}}\mathopen{}\mathclose{{}\left\lVert
a_{ni}\xi_{n,i}(u)}\right\rVert^{2}>\gamma^{2}}\right]}\right]$ (15)
$\displaystyle\hskip
56.9055pt\leq\sum_{i=1}^{k_{n}}\sum_{v\in\mathcal{V}}\mathopen{}\mathclose{{}\left\lVert
a_{ni}}\right\rVert^{2}\mathbb{E}\mathopen{}\mathclose{{}\left[\xi_{n,i}(v)^{2}\sum_{u\in\mathcal{V}}1\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\lVert
a_{ni}}\right\rVert\mathopen{}\mathclose{{}\left|\xi_{n,i}(u)}\right|>\gamma\mathopen{}\mathclose{{}\left|\mathcal{V}}\right|^{-1/2}}\right]}\right]$
(16) $\displaystyle\hskip
56.9055pt=\sum_{i=1}^{k_{n}}\mathopen{}\mathclose{{}\left\lVert
a_{ni}}\right\rVert^{2}\sum_{u,v\in\mathcal{V}}\mathbb{E}\mathopen{}\mathclose{{}\left[\xi_{n,i}(v)^{2}1\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\lVert
a_{ni}}\right\rVert\mathopen{}\mathclose{{}\left|\xi_{n,i}(u)}\right|>\gamma\mathopen{}\mathclose{{}\left|\mathcal{V}}\right|^{-1/2}}\right]}\right]$
(17) $\displaystyle\hskip 56.9055pt\leq
C\sum_{i=1}^{k_{n}}\mathopen{}\mathclose{{}\left\lVert
a_{ni}}\right\rVert^{2+K/q}$ (18)

for some fixed constant $C>0,$ chosen in accordance with Lemma A.3. This bound
converges to zero as $n\rightarrow\infty$. ∎

###  A.3 Proof of Theorem 3.1

Here we prove Theorem 3.1 from the main text.

###### Proof.

Expanding, we have that

$\sqrt{n}(\hat{\beta}_{n}^{b}-\hat{\beta}_{n})=\sqrt{n}(X_{n}^{T}X_{n})^{-1}X_{n}^{T}E^{b}_{n}=\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}x_{i}E^{b}_{n,i}.$

Applying Lemma B.2,
$\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}$
converges a.s. to $\Sigma_{X}^{-1}$. Moreover,
$(E^{b}_{n,i})_{n\in\mathbb{N},1\leq i\leq n}$ is a triangular array which is
mean-zero and i.i.d within rows so if we can show that its covariance
converges the result will follow by applying Proposition A.4 with $\eta=2$ and
$K=1$. To demonstrate this convergence, for each $u,v\in\mathcal{V}$,
conditional on $(X_{m},Y_{m})_{m\in\mathbb{N}}$

$\displaystyle\text{cov}(E^{b}_{n,1}(u),E^{b}_{n,1}(v))$
$\displaystyle=\sum_{j=1}^{n}\frac{1}{n}\mathopen{}\mathclose{{}\left(\hat{E}_{n,j}(u)-\frac{1}{n}\sum_{l=1}^{n}\hat{E}_{n,l}(u)}\right)\hat{E}_{n,j}(v)$
$\displaystyle=\sum_{j=1}^{n}\frac{1}{n}\mathopen{}\mathclose{{}\left(\hat{E}_{n,j}(u)-\frac{1}{n}\sum_{l=1}^{n}\hat{E}_{n,l}(u)}\right)\hat{E}_{n,j}(v)$
$\displaystyle=\frac{1}{n}\hat{E}_{n}(u)^{T}\hat{E}_{n}(v)-\mathopen{}\mathclose{{}\left(\frac{1}{n}\sum_{j=1}^{n}\hat{E}_{n,j}(u)}\right)\mathopen{}\mathclose{{}\left(\frac{1}{n}\sum_{j=1}^{n}\hat{E}_{n,j}(v)}\right)$

Now, letting $P_{n}=X_{n}(X_{n}^{T}X_{n})^{-1}X_{n}$ and letting $I_{n}$ be
the $n\times n$ identity matrix.

$\frac{1}{n}\hat{E}_{n}(u)^{T}\hat{E}_{n}(v)=\frac{1}{n}E_{n}^{T}(u)(I_{n}-P_{n})E_{n}(v)=\frac{1}{n}E_{n}^{T}(u)E_{n}(v)-\frac{1}{n}E_{n}^{T}(u)P_{n}E_{n}(v).$

We can write
$\frac{1}{n}E_{n}^{T}(u)E_{n}(v)=\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}(u)\epsilon_{i}(v)$,
which converges almost surely to $\textswab{c}(u,v)$ by the strong law of
large numbers. Moreover,

$\displaystyle\frac{1}{n}E_{n}^{T}(u)P_{n}E_{n}(v)$
$\displaystyle=\frac{1}{n}E_{n}^{T}(u)X_{n}(X_{n}^{T}X_{n})^{-1}X_{n}^{T}E_{n}(v)$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}E_{n}(u)}{n}}\right)^{T}\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}E_{n}(u)}{n}}\right)$

which converges almost surely to zero as $n\rightarrow\infty.$ Finally,

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}\hat{E}_{n,j}(u)$
$\displaystyle=\frac{1}{n}1_{n}^{T}(I_{n}-P_{n})E_{n}(u)=\frac{1}{n}1_{n}^{T}E_{n}(u)-\frac{1}{n}1_{n}^{T}X_{n}(X_{n}^{T}X_{n})^{-1}X_{n}^{T}E_{n}(u)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}(u)-\mathopen{}\mathclose{{}\left(\frac{1}{n}\sum_{i=1}^{n}x_{i}}\right)\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}\mathopen{}\mathclose{{}\left(\frac{1}{n}\sum_{i=1}^{n}x_{i}\epsilon_{i}(u)}\right)$

which converges almost surely to 0 as $n\rightarrow\infty.$ To show that the
variance converges, note that

$(\hat{\sigma}^{b}_{n})^{2}=\frac{1}{n}\sum_{i=1}^{n}(E^{b}_{n,i})^{2}-\mathopen{}\mathclose{{}\left(\frac{1}{n}\sum_{i=1}^{n}E^{b}_{n,i}}\right)^{2}$

The $E^{b}_{n,i}$ are i.i.d and mean-zero and the covariance of $E^{b}_{n,i}$
converges as shown above. As such by the Lindeberg CLT,
$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}E^{b}_{n,i}$ converges in distribution and,
dividing by $\sqrt{n}$, it follows that $\frac{1}{n}\sum_{i=1}^{n}E^{b}_{n,i}$
converges almost surely to zero as $n\rightarrow\infty.$ For the first term,
note that

$\mathbb{E}(E^{b}_{n,i})^{2}=\sum_{j=1}^{n}\frac{1}{n}\mathopen{}\mathclose{{}\left(\hat{E}_{n,j}-\frac{1}{n}\sum_{l=1}^{n}\hat{E}_{n,l}}\right)^{2}$
(19)

which converges to $\sigma^{2}$ almost surely as $n\rightarrow\infty.$ As such
the result follows by the triangular weak law of large numbers so long as we
can demonstrate that $\sup_{n\in\mathbb{N},1\leq i\leq
n}\mathbb{E}(E^{b}_{n,i})^{4}<\infty$. To show this note that for each
$n\in\mathbb{N}$ and $1\leq i\leq n$ and $1\leq b\leq B,$

$\displaystyle\mathbb{E}(E^{b}_{n,i})^{4}$
$\displaystyle=\sum_{j=1}^{n}\frac{1}{n}\mathopen{}\mathclose{{}\left(\hat{E}_{n,j}-\frac{1}{n}\sum_{l=1}^{n}\hat{E}_{n,l}}\right)^{4}=\frac{1}{n}\sum_{j=1}^{n}\mathopen{}\mathclose{{}\left(\epsilon_{j}-(P_{n}E_{n})_{j}-\frac{1}{n}\sum_{l=1}^{n}\mathopen{}\mathclose{{}\left(\epsilon_{l}-(P_{n}E_{n})_{l}}\right)}\right)^{4}$

Now $\lVert P_{n}E_{n}\rVert$ converges in probability to 0 by Lemma A.5 (see
below) and so

$\max_{1\leq l\leq
n}|(P_{n}E_{n})_{l}|\overset{\mathbb{P}}{\longrightarrow}0,$

since $\max_{1\leq l\leq n}(P_{n}E_{n})_{l}^{2}\leq\lVert
P_{n}E_{n}\rVert^{2}.$ In particular it follows that for $M>0$,

$\lim_{k\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(\max_{n\geq
k}\max_{1\leq l\leq n}|(P_{n}E_{n})_{l}|>M}\right)\rightarrow 0$ (20)

as $k\rightarrow\infty.$ For $k\in\mathbb{N}$, Let
$A_{k}=\mathopen{}\mathclose{{}\left\\{\max_{n\geq k}\max_{1\leq l\leq
n}|(P_{n}E_{n})_{l}|\leq M}\right\\}$, then equation (20) implies that
$\mathbb{P}\mathopen{}\mathclose{{}\left(\cup_{k}A_{k}}\right)=1$ since the
sets are nested. As such for $\omega\in\cup_{k}A_{k}$, $\omega$ is contained
in $A_{K}$ some $K=K(\omega)\in\mathbb{N}.$ It follows that

$\max_{n>K}\max_{1\leq l\leq n}|(P_{n}E_{n})_{l}|\leq M$

almost everywhere which implies that

$\max_{n\in\mathbb{N}}\max_{1\leq l\leq n}|(P_{n}E_{n})_{l}|\leq
M^{\prime}=M+\max_{1\leq n\leq K}\max_{1\leq l\leq n}|(P_{n}E_{n})_{l}|.$

We can thus bound $\mathbb{E}(E^{b}_{n,i})^{4}$ by

$\frac{1}{n}\sum_{j=1}^{n}\sum_{k=0}^{4}\mathopen{}\mathclose{{}\left(\epsilon_{j}-\frac{1}{n}\sum_{l=1}^{n}\epsilon_{l}}\right)^{k}(2M^{\prime})^{4-k}\leq(2M^{\prime})^{4}\frac{1}{n}\sum_{j=1}^{n}\sum_{k=0}^{4}\mathopen{}\mathclose{{}\left(\epsilon_{j}-\frac{1}{n}\sum_{l=1}^{n}\epsilon_{l}}\right)^{k}.$

The right hand side converges almost surely by the strong law of large numbers
to a quantity that is the same for each $i$. It follows that the supremum over
$i,n$ of $\mathbb{E}(E^{b}_{n,i})^{4}$ is bounded, a fact that is true almost
everywhere since $\mathbb{P}(A)=1.$ ∎

######  Lemma A.5.

Under Assumption 1, letting $P_{n}=X_{n}(X_{n}^{T}X_{n})^{-1}X_{n}^{T}$, as
$n\rightarrow\infty$,

$\lVert P_{n}E_{n}\rVert\overset{\mathbb{P}}{\longrightarrow}0.$

###### Proof.

We have,

$P_{n}E_{n}=X_{n}(X_{n}^{T}X_{n})^{-1}X_{n}^{T}E_{n}=\frac{X_{n}}{n^{0.45}}\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}E_{n}}{n^{0.55}}}\right).$

Thus,

$\lVert
P_{n}E_{n}\rVert=\mathopen{}\mathclose{{}\left\lVert\frac{X_{n}}{n^{0.45}}\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}E_{n}}{n^{0.55}}}\right)}\right\rVert\leq\mathopen{}\mathclose{{}\left\lVert\frac{X_{n}}{n^{0.45}}}\right\rVert\mathopen{}\mathclose{{}\left\lVert\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}}\right\rVert\mathopen{}\mathclose{{}\left\lVert\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}E_{n}}{n^{0.55}}}\right)}\right\rVert.$

$\frac{X_{n}^{T}E_{n}}{\sqrt{n}}$ converges in distribution (see e.g. the
proof of Lemma LABEL:lem:lm2) so
$\mathopen{}\mathclose{{}\left\lVert\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}E_{n}}{n^{0.55}}}\right)}\right\rVert\overset{\mathbb{P}}{\longrightarrow}0$
and $\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}$
converges almost surely to $\Sigma_{X}^{-1}$ by Lemma B.2. Applying the
Gershgorin circle theorem and the AM-RM inequality, we have

$\mathopen{}\mathclose{{}\left\lVert X_{n}}\right\rVert\leq\max_{1\leq i\leq
n}\sum_{j=1}^{p}|(X_{n})_{ij}|=\max_{1\leq i\leq
n}\sum_{j=1}^{p}|(x_{i})_{j}|\leq\frac{p}{\sqrt{p}}\max_{1\leq i\leq n}\lVert
x_{i}\rVert.$

$n^{-0.45}\max_{1\leq i\leq n}\lVert
x_{i}\rVert\overset{a.s.}{\longrightarrow}0$ since $\mathbb{E}(\lVert
x_{1}\rVert^{5/2})<\infty$, so in particular $\lVert
n^{-0.45}X_{n}\rVert\overset{a.s.}{\longrightarrow}0.$ Combining these results
and using Slutsky, it follows that $\lVert
P_{n}E_{n}\rVert\overset{\mathbb{P}}{\longrightarrow}0.$ ∎

##  Appendix B  Proofs for the main text

###  B.1 Proofs for Section 2

####  B.1.1 Proof of Claim 2.4

###### Proof.

The event

$\displaystyle\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left|R_{k}(\lambda)\cap\mathcal{N}}\right|>k-1}\right\\}$
$\displaystyle=\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left\\{(l,v)\in\mathcal{N}:p_{n,l}(v)\leq
t_{k}(\lambda)}\right\\}}\right|>k-1}\right\\}$
$\displaystyle=\mathopen{}\mathclose{{}\left\\{p^{n}_{(k:\mathcal{N})}\leq
t_{k}(\lambda)}\right\\}=\mathopen{}\mathclose{{}\left\\{t_{k}^{-1}(p^{n}_{(k:\mathcal{N})})\leq\lambda}\right\\}$

As such,

$\displaystyle\bigcup_{1\leq k\leq
K}\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left|R_{k}(\lambda)\cap\mathcal{N}}\right|>k-1}\right\\}=\mathopen{}\mathclose{{}\left\\{\min_{1\leq
k\leq
K}t_{k}^{-1}(p^{n}_{(k:\mathcal{N})})\leq\lambda}\right\\}=\mathopen{}\mathclose{{}\left\\{\min_{1\leq
k\leq
K\wedge|\mathcal{H}|}t_{k}^{-1}(p^{n}_{(k:\mathcal{N})})\leq\lambda}\right\\}.$

∎

######  Remark B.1.

This claim can be generalized to arbitrary $\zeta_{k}$. The result in that
case is that

$\displaystyle\text{JER}\mathopen{}\mathclose{{}\left((R_{k}(\lambda),\zeta_{k})_{1\leq
k\leq K}}\right)=\mathbb{P}\mathopen{}\mathclose{{}\left(\min_{1\leq k\leq
K\wedge|\mathcal{H}|}t_{k}^{-1}(p^{n}_{(\zeta_{k}+1:\mathcal{N})})\leq\lambda}\right).$

Throughout the main text we take $\zeta_{k}=k-1$, this can be motivated by the
fact that it implies that each individual rejection region $R_{k}(\lambda)$
controls the $k$-familywise error rate. However other choices provide valid
inference, see Blanchard et al. (2020) for a discussion of the different
choices of $\zeta_{k}$. As such the results in Section 4 can trivially be
generalized to arbitrary $\zeta_{k}$.

###  B.2 Proofs for Section 3

We will need the following useful Lemma which is Davenport et al. (2021)’s
Lemma 8.2.

######  Lemma B.2.

Suppose that $(X_{n})_{n\in\mathbb{N}}$ satisfies Assumption 1a and let
$\Sigma_{X}=\mathbb{E}\mathopen{}\mathclose{{}\left[x_{1}x_{1}^{T}}\right]$,
then $\Sigma_{X}$ is invertible and

$\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}\overset{a.s.}{\longrightarrow}\Sigma_{X}^{-1}.$

####  B.2.1 Convergence in the Linear Model

In this section we establish results for asymptotics of coefficients and test-
statistics in the linear model, written in terms of the framework of random
fields.

######  Lemma B.3.

Suppose that $(X_{n})_{n\in\mathbb{N}}$ and $(\epsilon_{n})_{n\in\mathbb{N}}$
satisfy Assumption 1. Then

$\sqrt{n}(\hat{\beta}_{n}-\beta)\overset{d}{\longrightarrow}\mathcal{G}(0,\textswab{c}_{\epsilon}\Sigma_{X}^{-1}).$

###### Proof.

For each $n\in\mathbb{N},$

$\sqrt{n}(\hat{\beta}_{n}-\beta)=\sqrt{n}(X_{n}^{T}X_{n})^{-1}X_{n}^{T}\epsilon_{n}=\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}x_{i}\epsilon_{i}.$

By the Central Limit Theorem,
$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}x_{i}\epsilon_{i}$ converges to a
$p$-dimensional Gaussian random field with covariance

$\displaystyle\text{cov}\mathopen{}\mathclose{{}\left(x_{1}\epsilon_{1}(u),x_{1}\epsilon_{1}(v)}\right)=\mathbb{E}\mathopen{}\mathclose{{}\left[x_{1}\epsilon_{1}(u)\epsilon_{1}(v)x_{1}^{T}}\right]=\mathbb{E}\mathopen{}\mathclose{{}\left[\epsilon_{1}(u)\epsilon_{1}(v)}\right]\mathbb{E}\mathopen{}\mathclose{{}\left[x_{1}x_{1}^{T}}\right]=\textswab{c}_{\epsilon}(u,v)\Sigma_{X}$

for $u,v\in\mathcal{V}$.
$\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}$
converges almost surely to $\Sigma_{X}^{-1}$ by Lemma B.2 and so the result
follows by applying Lemma A.2 and Slutsky as the limiting distribution has
covariance (for each $u,v,\in\mathcal{V}$)

$\Sigma_{X}^{-1}(\textswab{c}_{\epsilon}(u,v)\Sigma_{X})\Sigma_{X}^{-1}=\textswab{c}_{\epsilon}(u,v)\Sigma_{X}^{-1}.$

∎

Let $\textswab{c}^{\prime}:\mathcal{V}\times\mathcal{V}\rightarrow\mathbb{R}$
be the covariance function such that for all $u,v\in\mathcal{V}$

$\textswab{c}^{\prime}(u,v)=\rho_{\epsilon}(u,v)AC\Sigma_{X}^{-1}C^{T}A^{T}$
(21)

where $A\in\mathbb{R}^{L\times L}$ is a diagonal matrix with
$A_{ll}=(c_{l}^{T}\Sigma_{X}^{-1}c_{l})^{-1/2}$ for $1\leq l\leq L.$ Then we
having the following results.

######  Theorem B.4.

For $n\in\mathbb{N}$, let $S_{n}$ be the $L$-dimensional random field on
$\mathcal{V}$ defined by

$S_{n,l}=\frac{c_{l}^{T}(\hat{\beta}_{n}-\beta)}{\hat{\sigma}_{n}\sqrt{c_{l}^{T}(X_{n}^{T}X_{n})^{-1}c_{l}}}.$

for $1\leq l\leq L$. Then, under the conditions of Lemma B.3, as
$n\rightarrow\infty,$

$S_{n}\overset{d}{\longrightarrow}\mathcal{G}(0,\textswab{c}^{\prime})$

and it follows that

$T_{n}|_{\mathcal{N}}\overset{d}{\longrightarrow}\mathcal{G}(0,\textswab{c}^{\prime})|_{\mathcal{N}}.$

###### Proof.

We can write

$S_{n}=\sqrt{n}A_{n}C(\hat{\beta}_{n}-\beta_{n})/\hat{\sigma}_{n}.$

where $A_{n}$ is a diagonal matrix with
$(A_{n})_{ll}=\mathopen{}\mathclose{{}\left(c_{l}^{T}\mathopen{}\mathclose{{}\left(\frac{X_{n}^{T}X_{n}}{n}}\right)^{-1}c_{l}}\right)^{-1/2}$.
$A_{n}\overset{a.s.}{\longrightarrow}A$ by Lemma B.2 and
$\hat{\sigma}_{n}\overset{a.s.}{\longrightarrow}\sigma$ as
$n\rightarrow\infty$. So applying Lemmas B.3 and A.2 and Slutsky, the first
result follows. For $(v,l)\in\mathcal{N}$, $c^{T}_{l}\beta(v)=0$. As such
$S_{n}|_{\mathcal{N}}=T_{n}|_{\mathcal{N}}$ and it follows that

$T_{n}|_{\mathcal{N}}\overset{d}{\longrightarrow}\mathcal{G}|_{\mathcal{N}}.$

∎

####  B.2.2 Proof of Theorem 3.1

###### Proof.

Our proof of this result is available, see Section A.3. What follows here is
an alternative proof using Theorem 1 of Eck (2018).

Applying Eck (2018)’s Theorem 1 (conditioning on $(Y_{n})_{n\in\mathbb{N}}$
and restricting to the probability 1 event that
$(\frac{1}{n}X_{n}^{T}X_{n})^{-1}\longrightarrow\Sigma_{X}^{-1}$), we see that

$\sqrt{n}(\text{vec}(\hat{\beta}_{n}^{b})-\text{vec}(\hat{\beta}_{n}))\longrightarrow
N(0,\Sigma\otimes\Sigma_{X}^{-1}),$

where $\Sigma=\text{cov}(\text{vec}(\epsilon_{1}))$. It follows that
$\sqrt{n}(\hat{\beta}_{n}^{b}-\hat{\beta}_{n})$ converges in distribution to a
Gaussian random field which limiting covariance
$\textswab{c}_{\epsilon}\Sigma^{-1}_{X}$. The form of the covariance in the
statement of the theorem follows as writing
$\mathcal{V}=\mathopen{}\mathclose{{}\left\\{u_{1},\dots,u_{V}}\right\\}$, for
$1\leq l,m\leq L$ and $1\leq j,k\leq V$,

$(\Sigma\otimes\Sigma_{X}^{-1}))_{L(l-1)+j,L(m-1)+k}=\textswab{c}_{\epsilon}(u_{j},u_{k})(\Sigma_{X}^{-1})_{lm}.$
(22)

∎

######  Remark B.5.

Eck (2018)’s theorem needs to be applied with care as they write the model
$Y=\beta X+\epsilon$ rather than via the more standard formulation of
$Y=X\beta+\epsilon$, i.e. they takes $\beta$ to be a row vector rather than a
column vector. Their vec operation is thus the result of stacking a transposed
matrix the resulting distribution in the statement of their Theorem 1 is
$N(0,\Sigma_{X}^{-1}\otimes\Sigma)$ rather than
$N(0,\Sigma\otimes\Sigma_{X}^{-1})$.

######  Remark B.6.

Eck (2018)’s Theorem 1 is stated in terms of fixed design matrices which
converge. Here we assume that the design is random but condition on it which
allows us to apply their Theorem 1 because $(x_{n})_{n\in\mathbb{N}}$ and
$(y_{n})_{n\in\mathbb{N}}$ are independent. Eck (2018) has an alternative
result (their Theorem 2) which applies when $(x_{n},y_{n})_{n\in\mathbb{N}}$
has a joint distribution, however this requires an alternative form of the
bootstrap first introduced in Freedman (1981).

We prove this result an alternative, somewhat simpler way, using the Lindeberg
Central Limit Theorem. See Section A for details.

####  B.2.3 Proof of Theorem 3.2

###### Proof.

We can write

$T_{n}^{b}=\sqrt{n}A_{n}C(\hat{\beta}_{n}^{b}-\hat{\beta}_{n})/\hat{\sigma}^{b}_{n}.$

where $A_{n}$ is defined as in the proof of B.4. Applying Eck (2018)’s Theorem
1b it follows that, as $n\rightarrow\infty$,
$\hat{\sigma}_{n}^{b}\overset{a.s.}{\longrightarrow}\sigma$. Moreover
$A_{n}\overset{a.s.}{\longrightarrow}A$ so applying Theorem 3.1, Lemma A.2 and
Slutsky, the first result holds. The second result immediately follows from
the first. ∎

###  B.3 Proof of Theorem 3.3

###### Proof.

Let $F_{n}:\mathbb{R}\rightarrow[0,1]$ send $\lambda\in\mathbb{R}$ to
$\mathbb{P}\mathopen{}\mathclose{{}\left(f(T^{1}_{n})\leq\lambda\middle|(X_{m},Y_{m})_{m\in\mathbb{N}}}\right)$.
Define a sequence $(\eta_{n})_{n\in\mathbb{N}}\geq 0$ such that
$\alpha\pm\eta_{n}$ are continuity points of $F_{n}^{-}$ and
$\eta_{n}\rightarrow 0$ as $n\rightarrow\infty$. To do so let $\eta_{n}=0$ if
$\alpha$ is a continuity point of $F_{n}^{-}$ and take
$\eta_{n}=\frac{1}{2n^{n}}$ is $\alpha$ otherwise. Note that there are at most
$n^{n}$ distinct values that $f(T_{n}^{1})$ can take, so $F_{n}$ is a step
function with at most $n^{n}$ steps, meaning that the height difference
between steps is at least $\frac{1}{n^{n}}.$ The points of discontinuity of
$F_{n}^{-}$ are the values in the range of $F_{n}$ and so if $\alpha$ is not a
point of continuity of $F_{n}^{-}$ then $\alpha\pm\frac{1}{2n^{n}}$ must be.
Now,

$\lambda^{*}_{\alpha-\eta_{n},n,B}\leq\lambda^{*}_{\alpha,n,B}\leq\lambda^{*}_{\alpha+\eta_{n},n,B}.$
(23)

The values $\alpha\pm\eta_{n}$ are continuity points of $F_{n}^{-}$ and for
$\lambda\in\mathbb{R}$, conditional on $(X_{m},Y_{m})_{m\in\mathbb{N}}$, by
the SLLN,
$\frac{1}{B}\sum_{b=1}^{B}1\mathopen{}\mathclose{{}\left[f(T^{b}_{n})\leq\lambda}\right]$
converges almost surely to $F_{n}(\lambda)$ as $B\rightarrow\infty$. As such,
applying Lemma 1.1.1 from De Haan and Ferreira (2006),
$\lambda^{*}_{\alpha\pm\eta_{n},n,B}\rightarrow F_{n}^{-}(\alpha\pm\eta_{n})$
almost surely as $B\rightarrow\infty.$ Moreover, as $n\rightarrow\infty,$
$F_{n}$ converges pointwise to $F$ (as
$f(T^{1}_{n})|(X_{m},Y_{m})_{m\in\mathbb{N}}\overset{d}{\longrightarrow}f(\mathcal{G}(0,\textswab{c}^{\prime})))$
which is an increasing invertible function with continuous inverse. As such
$F_{n}^{-}(\alpha\pm\eta_{n})\rightarrow\lambda_{\alpha}$ as
$n\rightarrow\infty.$ To see this note that for all $\delta>0,$ there exists
$N\in\mathbb{N}$ such that for all $n\geq N,\eta_{n}<\delta$ and

$F_{n}^{-}(\alpha)\leq F_{n}^{-}(\alpha+\eta_{n})\leq
F_{n}^{-}(\alpha+\delta)$

and $F_{n}^{-}(\alpha+\delta)$ converges to $F^{-1}(\alpha+\delta)$ as
$n\rightarrow\infty$ by applying Lemma 1.1.1 from De Haan and Ferreira (2006)
once again. $F^{-1}$ is continuous and so $F^{-1}(\alpha+\delta)\rightarrow
F^{-1}(\alpha)$ as $\delta\rightarrow 0$. Arguing similarly for the sequence
$\alpha-\eta_{n}$ the result follows. Taking limits and using the bound in
equation (23), it almost surely follows that

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\lambda^{*}_{\alpha,n,B}=\lambda_{\alpha}.$

∎

###  B.4 Proofs for Section 4

####  B.4.1 Setup

In what follows we will require the following Lemma.

######  Lemma B.7.

Let $(F_{n})_{n\in\mathbb{N}},F$ be CDFs such that $F_{n}$ converges to $F$
pointwise and $F$ is continuous. Let
$(\lambda_{n})_{n\in\mathbb{N}}\in\mathbb{R}$ be a sequence such that
$\lambda_{n}\rightarrow\lambda\in\mathbb{R}$ as $n\rightarrow\infty$, then

$F_{n}(\lambda_{n})\rightarrow F(\lambda).$

###### Proof.

We can write

$F_{n}(\lambda_{n})-F(\lambda)=F_{n}(\lambda_{n})-F(\lambda_{n})+F(\lambda_{n})-F(\lambda).$

$F_{n}$ converges uniformly to $F$ (as CDFs which converge pointwise to a
continuous limit do so uniformly) so
$F_{n}(\lambda_{n})-F(\lambda_{n})\rightarrow 0$ as $n\rightarrow 0$ and
$F(\lambda_{n})-F(\lambda)\rightarrow 0$ because $F$ is continuous. ∎

Moreover we will want to restrict random fields to subsets. This is defined
formally as follows.

######  Definition B.8.

Given a a set valued function: $\mathcal{N}$ on $\mathcal{V}$, such that for
each
$v\in\mathcal{V},\mathcal{N}_{v}\subset\mathopen{}\mathclose{{}\left\\{1,\dots,L}\right\\}$,
we define the restriction of $f$ to $\mathcal{N}$ to be the map
$f|_{\mathcal{N}}:\Omega\rightarrow\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\bigcup_{1\leq
j\leq L}\mathbb{R}^{j}}\right\\}$ such that $f|_{\mathcal{N}}(\omega)(v)$ is
the vector
$(f_{k}(v):k\in\mathcal{N}_{v})^{T}\in\mathbb{R}^{\mathopen{}\mathclose{{}\left|\mathcal{N}_{v}}\right|}.$

Given a set function $\mathcal{N}$, defined as in Definition B.8, we can stack
the entries of $f|_{\mathcal{N}}$ to create vec$(f|_{\mathcal{N}})$ and thus
define
$f_{n}|_{\mathcal{N}}\overset{d}{\longrightarrow}f|_{\mathcal{N}}$,$f_{n}|_{\mathcal{N}}\overset{\mathbb{P}}{\longrightarrow}f|_{\mathcal{N}}$
and $f_{n}|_{\mathcal{N}}\overset{a.s.}{\longrightarrow}f|_{\mathcal{N}}$.
Because of the Central Limit Theorem convergence will typically be to a
Gaussian random field which is defined as follows.

######  Definition B.9.

Moreover, for a set function $\mathcal{N}$ as defined above, we shall write
$\mathcal{G}|_{\mathcal{N}}(\mu,\textswab{c})$ to denote the distribution of
the restricted random field. Given

$h:\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\mathbb{R}^{L}}\right\\}\rightarrow\mathbb{R}$

we shall write $X\sim h(\mathcal{G}(\mu,\textswab{c}))$ to indicate that $X$
is a real valued random variable which has the same distribution as $h(G)$
where $G\sim\mathcal{G}(\mu,\textswab{c}).$ Given

$h:\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\bigcup_{1\leq
j\leq L}\mathbb{R}^{j}}\right\\}\rightarrow\mathbb{R}$ (24)

we similarly define the notation
$h(\mathcal{G}|_{\mathcal{N}}(\mu,\textswab{c}))$.

####  B.4.2 Proof of Theorem 4.1

In order to facilitate the proof we will first make some further definitions.
Firstly, given $H\subset\mathcal{H}$ and
$T:\mathcal{V}\rightarrow\mathbb{R}^{L}$ define $p_{(k:H)}(T)$ to be the
minimum value in the set

$\mathopen{}\mathclose{{}\left\\{2-2\Phi(\mathopen{}\mathclose{{}\left|T_{l}(v)}\right|):(l,v)\in
H}\right\\}$ (25)

where $\Phi$ is the CDF of a standard normal distribution. Secondly, given
$H\subset\mathcal{H}$, let
$f_{H}:\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\mathbb{R}^{L}}\right\\}\rightarrow\mathbb{R}$
send
$T\in\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\mathbb{R}^{L}}\right\\}$
to $\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|}t_{k}^{-1}(p_{(k:H)}(T))$.
Thirdly given a function $S$ such that

$S:\mathcal{V}\rightarrow\bigcup_{0\leq j\leq L}\mathbb{R}^{j},$

$n\in\mathbb{N}$ and $1\leq k\leq|\mathcal{H}|$ we shall define $q^{n}_{k}(S)$
to be the $k$th minimum value in the set

$\mathopen{}\mathclose{{}\left\\{2-2\Phi_{n-r}(\mathopen{}\mathclose{{}\left|S_{l}(v)}\right|):v\in\mathcal{V},l\leq\text{dim}(S(v))}\right\\}$

when this is well defined and take $q^{n}_{k}(S)$ to be 1 when it is not (i.e.
when $k$ is larger than the size of the set). Here for $z\in\bigcup_{1\leq
j\leq L}\mathbb{R}^{j}$, $\text{dim}(z)$ denotes the dimension of $z$.
Similarly define $q_{k}(S)$ to be the $k$th minimum value in the set

$\mathopen{}\mathclose{{}\left\\{2-2\Phi(\mathopen{}\mathclose{{}\left|S_{l}(v)}\right|):v\in\mathcal{V},l\leq\text{dim}(S(v))}\right\\}.$

Finally we define functions
$\phi_{n}:\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\bigcup_{0\leq
j\leq L}\mathbb{R}^{j}}\right\\}\rightarrow\mathbb{R}$ which send

$S\mapsto\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|}t_{k}^{-1}(q^{n}_{k}(S))$

and
$\phi:\mathopen{}\mathclose{{}\left\\{g:\mathcal{V}\rightarrow\bigcup_{0\leq
j\leq L}\mathbb{R}^{j}}\right\\}\rightarrow\mathbb{R}$ which sends
$S\mapsto\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|}t_{k}^{-1}(q_{k}(S)).$

With these definitions in mind we are ready to prove Theorem 4.1.

###### Proof.

Defining $\textswab{c}^{\prime}$ as in Section 3.2, $T_{n}|_{\mathcal{N}}$
converges to $\mathcal{G}(0,\textswab{c}^{\prime})|_{\mathcal{N}}$ in
distribution by Theorem B.4. As such, using the fact that $\phi_{n}$ is the
composition of functions which are either continuous or converge uniformly
with range $[0,1]$ (since the minimum is continuous), by Lemma C.3 and the
Continuous Mapping Theorem,

$f_{n,\mathcal{N}}(T_{n})=\phi_{n}(T_{n}|_{\mathcal{N}})\overset{d}{\longrightarrow}\phi(\mathcal{G}(0,\textswab{c}^{\prime})|_{\mathcal{N}})=f_{\mathcal{N}}(\mathcal{G}(0,\textswab{c}^{\prime})).$
(26)

By the same logic, and applying Theorem 3.2, for sets $H$ such that
$\mathcal{N}\subset H\subset\mathcal{H}$,

$f_{n,H}(T_{n}^{b})\overset{d}{\longrightarrow}f_{H}(\mathcal{G}(0,\textswab{c}^{\prime})).$
(27)

This convergence occurs conditional on the data, a fact that we take as
implicit in (27) and in the rest of the proof. As such, applying Theorem 3.3,
it follows almost surely that
$\lambda^{*}_{\alpha,n,B}(H)\rightarrow\lambda_{\alpha}=F^{-1}(\alpha)$, where
$F$ is the CDF of $f_{H}(\mathcal{G}(0,\textswab{c}^{\prime}))$ using the fact
that $F$ is strictly increasing (which follows from the form of $f_{H}$ and
the fact that the density of the multivariate normal distribution is positive
everywhere) and continuous, by Lemma C.4. Letting $F_{n}$ be the CDF of
$f_{n,\mathcal{N}}(T_{n})$ and $F_{0}$ be the CDF of
$f_{\mathcal{N}}(\mathcal{G}(0,\textswab{c}^{\prime}))$, we have
$F_{n}\longrightarrow F_{0}$ pointwise using (26) and the fact that $F_{0}$ is
continuous (which follows from Lemma C.4). As such, applying Lemma B.7 (since
$F_{n}$ and $F_{0}$ are CDFs and $F_{0}$ is continuous), it follows that for
all $\epsilon>0,$

$\displaystyle\lim_{n\rightarrow\infty}$
$\displaystyle\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(f_{n,\mathcal{N}}(T_{n})\leq\lambda^{*}_{\alpha,n,B}(H)}\right)$
$\displaystyle\quad\leq\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(f_{n,\mathcal{N}}(T_{n})\leq\lambda_{\alpha}+\epsilon}\right)+\mathbb{P}(|\lambda^{*}_{\alpha,n,B}(H)-\lambda_{\alpha}|>\epsilon)$
$\displaystyle\quad=\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}F_{n}(\lambda_{\alpha}+\epsilon)+\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}(|\lambda^{*}_{\alpha,n,B}(H)-\lambda_{\alpha}|>\epsilon)$
$\displaystyle\quad=F_{0}(\lambda_{\alpha}+\epsilon)\leq
F(\lambda_{\alpha}+\epsilon)=\alpha+\epsilon.$

Taking $\epsilon$ to zero proves the result. Note that the inequality holds
because

$f_{H}(\mathcal{G}(0,\textswab{c}^{\prime}))=\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|}t_{k}^{-1}(p_{(k:H)}(\mathcal{G}(0,\textswab{c}^{\prime})))\leq\min_{1\leq
k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|}t_{k}^{-1}(p_{(k:\mathcal{N})}(\mathcal{G}(0,\textswab{c}^{\prime})))=f_{\mathcal{N}}(\mathcal{G}(0,\textswab{c}^{\prime}))$

and so

$F_{0}(\lambda_{\alpha})=\mathbb{P}\mathopen{}\mathclose{{}\left(f_{\mathcal{N}}(\mathcal{G}(0,\textswab{c}^{\prime}))\leq\lambda_{\alpha}}\right)\leq\mathbb{P}\mathopen{}\mathclose{{}\left(f_{H}(\mathcal{G}(0,\textswab{c}^{\prime}))\leq\lambda_{\alpha}}\right)=F(\lambda_{\alpha}).$

∎

####  B.4.3 Proof of Corollary 4.2

###### Proof.

For any $\epsilon>0$, and all large enough $n$ and $B$, we have

$\mathbb{P}\mathopen{}\mathclose{{}\left(\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{H}}\right|}t_{k}^{-1}(p_{(k:\mathcal{N})}(T_{n}))\leq\lambda^{*}_{\alpha,n,B}(\mathcal{H})}\right)\leq\alpha+\epsilon$

and so, arguing as in Blanchard et al. (2020),

$\mathbb{P}\mathopen{}\mathclose{{}\left(|H\cap\mathcal{N}|\leq\makebox[0.0pt]{$\phantom{V}\overline{\phantom{V}}$}V_{\alpha,n,B}(H),\,\,\forall
H\subset\mathcal{H}}\right)\leq 1-\alpha-\epsilon.$

The result follows by sending $\epsilon$ to zero. ∎

####  B.4.4 Proof of Theorem 4.3

The proof is similar to that of Proposition 4.5 of Blanchard et al. (2020).

###### Proof.

Let

$\Omega_{n}=\mathopen{}\mathclose{{}\left\\{p^{n}_{(k:\mathcal{N})}(T_{n})\geq
t_{k}(\lambda^{*}_{\alpha,n,B}(\mathcal{N}))\text{ for all }1\leq k\leq
K}\right\\}.$

Then by Theorem 4.1,

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}(\Omega_{n})=1-\alpha.$

We claim that on the event $\Omega_{n}$, $\mathcal{N}\subset\hat{H}_{n}$. We
prove this inductively, using the notation from Algorithm 1.
$\mathcal{N}\subset H^{(0)}$ trivially. Assuming that $\mathcal{N}\subset
H^{(j-1)}$ for some $j\in\mathbb{N},$ it follows that
$p^{n}_{(k:H^{(j-1)})}\leq p^{n}_{(k:\mathcal{N})}$ i.e. that
$f_{n,H^{(j-1)}}\leq f_{n,\mathcal{N}}$. In particular,

$\lambda^{*}_{\alpha,n,B}(H^{(j-1)})\leq\lambda^{*}_{\alpha,n,B}(\mathcal{N})$

and thus (since we are on $\Omega_{n}$),

$p^{n}_{(1:\mathcal{N})}(T_{n})\geq
t_{1}(\lambda^{*}_{\alpha,n,B}(H^{(j-1)}))$

which implies that $\mathcal{N}\subset H^{(j)}$. Thus
$\mathcal{N}\subset\hat{H}_{n}$ and so for all $1\leq k\leq K,$

$p^{n}_{(k:\mathcal{N})}(T_{n})\geq
t_{k}(\lambda^{*}_{\alpha,n,B}(\mathcal{N}))\geq
t_{k}(\lambda^{*}_{\alpha,n,B}(\hat{H}_{n}))$

and so

$f_{n,\mathcal{N}}(T_{n})\geq\lambda^{*}_{\alpha,n,B}(\hat{H}_{n}).$

The post hoc bound result follows as in the proof of Corollary 4.2. ∎

##  Appendix C  Further Theory

###  C.1 FWER inference

FWER inference is commonly used in brain imaging in order to identify areas of
activation in the brain. This corresponds to performing multiple testing
inference on the data and returning a set of active hypotheses
$R\subset\mathcal{H}$ such that the familywise error rate (FWER), defined as

$\text{FWER}=\mathbb{P}\mathopen{}\mathclose{{}\left(R\cap\mathcal{N}}\right)\leq\alpha.$

When a single test is being used (for a single contrast or an $F$-test at each
voxel), brain imaging studies have typically used a permutation based
procedure (Winkler et al., 2014) in order to control these error rates. In the
case of multiple contrasts this approach is not always applicable - see
Section C.2. However the bootstrap approach can be applied. In particular we
have the following theorem which follows as a corollary of Theorem 4.1 by
taking $K=1$ and using the linear template.

######  Theorem C.1.

For $0\leq\alpha\leq 1$ and $n,B\in\mathbb{N}$, let
$\lambda^{\prime}_{\alpha,n,B}$ be the $\alpha$-quantile of the bootstrap
distribution (based on $B$ bootstraps) of

$p_{1:\mathcal{H}}(T_{n})=\min_{(l,v)\in\mathcal{H}}p_{n,l}(v).$

Let
$R_{n,B}=\mathopen{}\mathclose{{}\left\\{(l,v)\in\mathcal{H}:p_{n,l}(v)\leq\lambda^{\prime}_{\alpha,n,B}}\right\\}$.
Then

$\lim_{n\rightarrow\infty}\lim_{B\rightarrow\infty}\mathbb{P}\mathopen{}\mathclose{{}\left(R_{n,B}\cap\mathcal{N}}\right)\leq\alpha.$

So choosing $R_{n,B}$ as the rejection set provides asymptotic control of the
FWER.

This control of the FWER does not occur simultaneously with the control of the
joint error rate. However let $\lambda^{*}_{\alpha,n,B}$ is the
$\alpha-$quantile of the bootstrap distribution of $f_{n,\mathcal{N}}$ (as
defined in the statement of Theorem 4.1). Then FWER is automatically entailed
with control of the joint error rate by using the rejection set
$R=\mathopen{}\mathclose{{}\left\\{(l,v)\in\mathcal{H}:p_{n,l}\leq
t_{1}^{-1}(\lambda^{*}_{\alpha,n,B})}\right\\}$. When $K>1$, typically
$t_{1}^{-1}(\lambda^{*}_{\alpha,n,B})$ will be less than the value of
$\lambda^{\prime}_{\alpha,n,B}$ from Theorem C.1 so this will result in less
power but comes with the advantage of holding jointly with control of the
joint error rate. Which version is to be preferred depends on which error rate
one desires to control.

Figure 7: Direct FWER control using the different methods. Here the parametric
procedures are Bonferroni and its step-down: Holm (1979).

###  C.2 Permutation in the Linear Model

Here we show that under the alternative that $\beta\not=0$ at a given point
(e.g. voxel or gene), permuting the data does not necessarily generate data
under the global null even when the noise is exchangeable under permutation.

######  Claim C.2.

Suppose that the global null is not true, i.e. $\beta\not=0$, then permuting
$Y$ is not equivalent to generating data under the global null (and so cannot
be used to generate under the null and provide strong control over contrasts).

###### Proof.

Let $P$ be a permutation matrix, then

$PY=P(X\beta+\epsilon)=PX\beta+P\epsilon.$

Now

$P\epsilon\sim\epsilon$

by exchangability. However $PX\beta\not=0$ so it is not true that
$PY\sim\epsilon$ which is what we want (because we need to simulate under the
null model, in order to apply the subset pivotality condition to provide
strong control over contrasts). $PX\beta$ is a random variable (due to
randomness in $P$) with a non-zero mean and variance. So regressing $PY$
against $X$ gives linear model coefficients of

$\displaystyle\hat{\beta}$
$\displaystyle=(X^{T}X)^{-1}X^{T}PY=(X^{T}X)^{-1}X^{T}P(X\beta+\epsilon)$
$\displaystyle=(X^{T}X)^{-1}X^{T}PX\beta+(X^{T}X)^{-1}X^{T}P\epsilon.$

Now, under exchangeability,

$(X^{T}X)^{-1}X^{T}P\epsilon\sim(X^{T}X)^{-1}X^{T}\epsilon$

which indeed is the distribution of the linear model estimates under the null,
however

$(X^{T}X)^{-1}X^{T}PX\beta\not=0$

which causes a problem. ∎

###  C.3 Additional Lemmas for the proofs

######  Lemma C.3.

Suppose that $(Z_{n})_{n\in\mathbb{N}},Z$ are $\mathbb{R}^{M}$ valued random
variables, for some $M\in\mathbb{N}$. Let $(f_{n})_{n\in\mathbb{N}},f$ be
functions from $\mathbb{R}^{M}\rightarrow I$ for some compact set
$I\subset\mathbb{R}$. Suppose that $f_{n}$ converges uniformly to $f$, that
$f$ is continuous and that $Z_{n}\overset{d}{\longrightarrow}Z$, then

$f_{n}(Z_{n})\overset{d}{\longrightarrow}f(Z).$

###### Proof.

Given any continuous and bounded function $g:\mathbb{R}\rightarrow\mathbb{R}$,

$\mathopen{}\mathclose{{}\left|\mathbb{E}\mathopen{}\mathclose{{}\left[g(f_{n}(Z_{n}))}\right]-\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z))}\right]}\right|\leq\mathopen{}\mathclose{{}\left|\mathbb{E}\mathopen{}\mathclose{{}\left[g(f_{n}(Z_{n}))}\right]-\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z_{n}))}\right]}\right|+\mathopen{}\mathclose{{}\left|\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z_{n}))}\right]-\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z))}\right]}\right|.$

the functions $(f_{n})_{n\in\mathbb{N}},f$ range values within a compact set
$I$ and so without loss of generality we may assume that $g$ is uniformly
continuous. So for any $\epsilon>0$ there is some $\delta$ such that
$\mathopen{}\mathclose{{}\left|g(x)-g(y)}\right|<\epsilon$ for all
$x,y\in\mathbb{R}$ such that
$\mathopen{}\mathclose{{}\left|x-y}\right|<\delta$. By uniform convergence,
there is some $N\in\mathbb{N}$ such that for all $n>N$,
$\mathopen{}\mathclose{{}\left|f_{n}(z)-f(z)}\right|<\delta$ for all
$z\in\mathbb{R}^{L}$. As such

$\mathopen{}\mathclose{{}\left|\mathbb{E}\mathopen{}\mathclose{{}\left[g(f_{n}(Z_{n}))}\right]-\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z_{n}))}\right]}\right|\leq\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|g(f_{n}(Z_{n}))-g(f(Z_{n}))}\right|}\right]<\mathbb{E}\mathopen{}\mathclose{{}\left[\epsilon}\right]=\epsilon.$

So this term converges to zero as $n\rightarrow\infty$. The second term:
$\mathopen{}\mathclose{{}\left|\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z_{n}))}\right]-\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z))}\right]}\right|$
also converges to zero as $g\circ f$ is a continuous bounded function and
$Z_{n}\overset{d}{\longrightarrow}Z$ as $n\rightarrow\infty$ (by applying the
Portmanteau Theorem). Thus, as $n\rightarrow\infty,$

$\mathbb{E}\mathopen{}\mathclose{{}\left[g(f_{n}(Z_{n}))}\right]\rightarrow\mathbb{E}\mathopen{}\mathclose{{}\left[g(f(Z))}\right].$

Since this holds for any continuous bounded $g$ the result follows by
Portmanteau. ∎

######  Lemma C.4.

Let $F_{0}$ be the CDF of $\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|}t_{k}^{-1}(p_{(k:H)}(T))$
where $T\sim\mathcal{G}(0,\textswab{c}^{\prime})$ and $\textswab{c}^{\prime}$
is defined as in Section 3.2, then $F_{0}$ is continuous.

###### Proof.

It is sufficient to show that for all $\lambda\in\mathbb{R}$,
$\mathbb{P}\mathopen{}\mathclose{{}\left(\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|}t_{k}^{-1}(p_{(k:\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|)}(T))=\lambda}\right)=0.$
To show this, choose $\lambda\in\mathbb{R}$, then

$\displaystyle\mathbb{P}\mathopen{}\mathclose{{}\left(\min_{1\leq k\leq
K\wedge\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|}t_{k}^{-1}(p_{(k:H)}(T))=\lambda}\right)$
$\displaystyle\leq\mathbb{P}\mathopen{}\mathclose{{}\left(\exists 1\leq
k\leq\mathopen{}\mathclose{{}\left|\mathcal{N}}\right|\text{ s.t.
}t_{k}^{-1}(p_{(k:\mathcal{N})}(T))=\lambda}\right)$
$\displaystyle=\mathbb{P}\mathopen{}\mathclose{{}\left(\exists 1\leq k\leq
m\text{ s.t. }p_{(k:\mathcal{N})}(T)=t_{k}(\lambda)}\right)$
$\displaystyle\leq\sum_{k=1}^{m}\mathbb{P}\mathopen{}\mathclose{{}\left(p_{(k:\mathcal{N})}(T)=t_{k}(\lambda)}\right)$
$\displaystyle\leq\sum_{k=1}^{m}\mathbb{P}\mathopen{}\mathclose{{}\left(\exists(l,v)\in\mathcal{N}:2(1-\Phi(\mathopen{}\mathclose{{}\left|T_{l}(v)}\right|))=t_{k}(\lambda)}\right)$
$\displaystyle\leq\sum_{k=1}^{m}\sum_{(l,v)\in\mathcal{H}}\mathbb{P}\mathopen{}\mathclose{{}\left(2(1-\Phi(\mathopen{}\mathclose{{}\left|T_{l}(v)}\right|))=t_{k}(\lambda)}\right).$

Now given $(l,v)\in\mathcal{H}$ and $1\leq k\leq m,$

$\mathbb{P}\mathopen{}\mathclose{{}\left(2(1-\Phi(\mathopen{}\mathclose{{}\left|T_{l}(v)}\right|))=t_{k}(\lambda)}\right)=\mathbb{P}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left|T_{l}(v)}\right|=\Phi^{-1}\mathopen{}\mathclose{{}\left(1-t_{k}(\lambda)/2}\right)}\right)=0$

since $T_{l}(v)$ is a Gaussian random variable. The result follows. ∎

##  Appendix D  fMRI data pre-processing

Participants underwent a working memory task in which they were shown images
asked to remember them. They were reshown them at a subsequent point. This is
known as an $m$-back task when $m\in\mathbb{N}$ is number of intervals between
when each image is shown and then repeated - see Barch et al. (2013) for
further details. The data we have consists of images that give the difference
between the brain scans of participants under the 2-back and 0-back
conditions. The data was pre-processed at the first level using nilearn. the
images were then smoothed using an isotropic Gaussian kernel with an FWHM of
$4/3$ voxels (4 mm).

##  Appendix E  Further figures

###  E.1 Simes vs ARI for the IQ contrast

Figure 8: TDP bounds within clusters for the contrast for IQ in the linear
regression model fit to the HCP data. Each cluster is shaded a single colour
which is the lower bound on the TDP. The upper panel gives the TDP bounds
within each cluster provided by the bootstrap procedure. The lower panels
gives the bounds provided by using ARI and the Simes procedure. The bounds
given by the bootstrap are larger (as indicated by the light colours)
indicating that the method is more powerful. (Note that the step-down
bootstrap gave the same bounds as the bootstrap and so is not shown.) Note
that these images are 2D slices through the 3D brain and so voxels that are
part of the same cluster are not necessarily connected.

###  E.2 The contrast for sex

Much less activation is found for the contrast of sex in the linear model fit
to the HCP data. In this case only a single cluster above the cluster defining
threshold has non-zero lower bound. The bound provided is the same for all the
parametric and bootstrap methods that we consider, in particular they all
conclude that at least one of the 17 voxels within this cluster has non-zero
activation. The cluster (and its TDP) is illustrated in Figure 9.

Figure 9: Illustrating the cluster in the sex contrast with non-zero
activation.

###  E.3 Illustrating the simulation setup

Figure 10: Illustrating the simulation setup for a domain of size [25,25] and
a smoothness of 4 pixels. Left: the signal for the first contrast. Right: a
realisation of one of the subjects in $G_{2}.$

###  E.4 Additional JER control plots

In this section we present the results of the simulations to consider JER
control where the domain of the data in the simulations is $25$ by $25$ or
$100$ by $100$ rather than $50$ by $50$. The results for the $25$ by $25$
simulations are shown in Figure 11 and those for the $100$ by $100$
simulations are shown in Figure 12. The results are similar to those in the
main text.

Figure 11: Empirical joint error rate for the $25$ by $25$ simulations.

Figure 12: Empirical joint error rate for the $100$ by $100$ simulations.

###  E.5 Additional power plots

Figure 13: Plotting the power of the different methods against the numbers of
a subjects for setting 3, i.e. taking
$R=\mathopen{}\mathclose{{}\left\\{(l,v):p_{n,l}(v)\leq 0.05}\right\\}$ in
(14).

## References

  * Alberton et al. (2020) Bianca AV Alberton, Thomas E Nichols, Humberto R Gamba, and Anderson M Winkler.  Multiple testing correction over contrasts for brain imaging.  _NeuroImage_ , 216:116760, 2020.
  * Andersen et al. (1999) Anders H. Andersen, Don M. Gash, and Malcolm J. Avison.  Principal component analysis of the dynamic response measured by fMRI: A generalized linear systems framework.  _Magnetic Resonance Imaging_ , 17(6):795–815, 1999.  ISSN 0730725X.  doi: 10.1016/S0730-725X(99)00028-4.
  * Andreella et al. (2020) Angela Andreella, Jesse Hemerik, Wouter Weeda, Livio Finos, and Jelle Goeman.  Permutation-based true discovery proportions for fmri cluster analysis.  arXiv preprint arXiv:2012.00368, 2020.
  * Bahr et al. (2013) Timothy M Bahr et al.  Peripheral blood mononuclear cell gene expression in chronic obstructive pulmonary disease.  _American journal of respiratory cell and molecular biology_ , 49(2):316–323, 2013.
  * Barch et al. (2013) Deanna M Barch, Gregory C Burgess, Michael P Harms, Steven E Petersen, Bradley L Schlaggar, Maurizio Corbetta, Matthew F Glasser, Sandra Curtiss, Sachin Dixit, Cindy Feldt, et al.  Function in the human connectome: task-fmri and individual differences in behavior.  _Neuroimage_ , 80:169–189, 2013.
  * Benjamini and Hochberg (1995) Yoav Benjamini and Yosef Hochberg.  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  _Journal of the Royal statistical society: series B (Methodological)_ , 57(1):289–300, 1995.
  * Benjamini and Yekutieli (2001) Yoav Benjamini and Daniel Yekutieli.  The control of the false discovery rate in multiple testing under dependency.  _Annals of statistics_ , pages 1165–1188, 2001.
  * Bickel and Freedman (1981) Peter J Bickel and David A Freedman.  Some Asymptotic Theory for the Bootstrap.  _Annals of Statistics_ , 9(6):1196–1217, 1981\.  ISSN 00905364.  doi: 10.1214/aos/1176342871.
  * Blain et al. (2022) Alexandre Blain, Bertrand Thirion, and Pierre Neuvial.  Notip: Non-parametric true discovery proportion control for brain imaging.  _NeuroImage_ , page 119492, 2022.  URL https://hal.archives-ouvertes.fr/hal-03649114.
  * Blanchard et al. (2020) Gilles Blanchard, Pierre Neuvial, Etienne Roquain, et al.  Post hoc confidence bounds on false positives using reference families.  _Annals of Statistics_ , 48(3):1281–1303, 2020\.
  * Blanchard et al. (2021) Gilles Blanchard, Pierre Neuvial, and Etienne Roquain.  On agnostic post hoc approaches to false positive control.  In _Handbook of Multiple Comparisons_ , Handbooks of Modern Statistical Methods. Chapman & Hall/CRC, November 2021.  URL https://hal.archives-ouvertes.fr/hal-02320543.
  * Cui and Churchill (2003) Xiangqin Cui and Gary A Churchill.  Statistical tests for differential expression in cDNA microarray experiments.  _Genome Biol_ , 4(4):210, 2003.
  * Davenport and Nichols (2022) Samuel Davenport and Thomas E. Nichols.  The expected behaviour of random fields in high dimensions: contradictions in the results of [1].  _Magnetic Resonance Imaging_ , 2022.
  * Davenport et al. (2021) Samuel Davenport, Fabian Telschow, Thomas E. Nichols, and Armin Schwarzman.  Confidence regions for the location of peaks of a smooth random field.  2021\.
  * De Haan and Ferreira (2006) Laurens De Haan and Ana Ferreira.  _Extreme value theory: an introduction_ , volume 21.  Springer, 2006.
  * Ebrahimpoor et al. (2020) Mitra Ebrahimpoor, Pietro Spitali, Kristina Hettne, Roula Tsonaka, and Jelle Goeman.  Simultaneous enrichment analysis of all possible gene-sets: unifying self-contained and competitive methods.  _Briefings in bioinformatics_ , 21(4):1302–1312, 2020.
  * Eck (2018) Daniel J Eck.  Bootstrapping for multivariate linear regression models.  _Statistics & Probability Letters_, 134:141–149, 2018\.
  * Enjalbert-Courrech and Neuvial (2022) Nicolas Enjalbert-Courrech and Pierre Neuvial.  Powerful and interpretable control of false discoveries in differential expression studies.  bioRxiv preprint: https://doi.org/10.1101/2022.03.08.483449, 2022.
  * Freedman (1981) David A Freedman.  Bootstrapping regression models.  _The Annals of Statistics_ , 9(6):1218–1228, 1981\.
  * Genovese and Wasserman (2006) Christopher R Genovese and Larry Wasserman.  Exceedance control of the false discovery proportion.  _Journal of the American Statistical Association_ , 101(476):1408–1417, 2006.
  * Genovese et al. (2002) Christopher R Genovese, Nicole A Lazar, and Thomas E Nichols.  Thresholding of statistical maps in functional neuroimaging using the false discovery rate.  _Neuroimage_ , 15(4):870–878, 2002.
  * Goeman and Solari (2011) Jelle J. Goeman and Aldo Solari.  Multiple Testing for Exploratory Research.  _Statistical Science_ , 26(4):584–597, 2011.  ISSN 0883-4237.  doi: 10.1214/11-STS356.
  * Goeman et al. (2019) Jelle J Goeman, Rosa J Meijer, Thijmen JP Krebs, and Aldo Solari.  Simultaneous control of all false discovery proportions in large-scale multiple hypothesis testing.  _Biometrika_ , 106(4):841–856, 2019.
  * Hemerik et al. (2019) Jesse Hemerik, Aldo Solari, and Jelle J Goeman.  Permutation-based simultaneous confidence bounds for the false discovery proportion.  _Biometrika_ , 106(3):635–649, 2019.
  * Henson (2006) R. Henson.  Forward inference using functional neuroimaging: Dissociations versus associations.  _Trends in cognitive sciences_ , 10:64–69, 2006.
  * Holm (1979) Sture Holm.  A simple sequentially rejective multiple test procedure.  _Scand. J. Statist._ , 6(2):65–70, 1979.  ISSN 0303-6898.
  * Hommel (1988) Gerhard Hommel.  A stagewise rejective multiple test procedure based on a modified bonferroni test.  _Biometrika_ , 75(2):383–386, 1988.
  * Korn et al. (2004) Edward L Korn, James F Troendle, Lisa M McShane, and Richard Simon.  Controlling the number of false discoveries: application to high-dimensional genomic data.  _Journal of Statistical Planning and Inference_ , 124(2):379–398, 2004.
  * Meinshausen (2006) Nicolai Meinshausen.  False discovery control for multiple tests of association under general dependence.  _Scandinavian Journal of Statistics_ , 33(2):227–237, 2006.
  * Neuvial (2020) Pierre Neuvial.  _Contributions to statistical inference from genomic data_.  Habilitation thesis, Université Toulouse III (France), 2020.  Available from https://tel.archives-ouvertes.fr/tel-02969229.
  * Romano and Wolf (2005) Joseph P Romano and Michael Wolf.  Exact and Approximate Stepdown Methods for Multiple Hypothesis Testing.  _Journal of the American Statistical Association_ , 100(469):94–108, 2005.  ISSN 0162-1459.  doi: 10.1198/016214504000000539.
  * Rosenblatt et al. (2018) Jonathan D Rosenblatt, Livio Finos, Wouter D Weeda, Aldo Solari, and Jelle J Goeman.  All-resolutions inference for brain imaging.  _Neuroimage_ , 181:786–796, 2018.
  * Sarkar et al. (2008) Sanat K Sarkar et al.  On the simes inequality and its generalization.  In _Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen_ , pages 231–242. Institute of Mathematical Statistics, 2008.
  * Simes (1986) R J Simes.  An improved Bonferroni procedure for multiple tests of significance.  _Biometrika_ , 73(3):751–754, 1986.
  * Smyth (2004) Gordon K. Smyth.  Linear models and empirical bayes methods for assessing differential expression in microarray experiments.  _Statistical Methods in Genetics and Molecular Biology_ , 3(3), 2004.  doi: 10.2202/1544-6115.1027.
  * Storey and Tibshirani (2003) John D Storey and Robert Tibshirani.  Statistical significance for genomewide studies.  _Proceedings of the National Academy of Sciences_ , 100(16):9440–9445, 2003.
  * Van der Vaart (2000) Aad W Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * van der Vaart (1998) A.W. van der Vaart.  _Asymptotic Statistics_.  1998\.
  * Westfall (2011) Peter H. Westfall.  On Using the Bootstrap for Multiple Comparisons.  _Journal of Biopharmaceutical Statistics_ , 21(6):1187–1205, 2011.  ISSN 1054-3406.  doi: 10.1080/10543406.2011.607751.  URL http://www.tandfonline.com/doi/abs/10.1080/10543406.2011.607751.
  * Winkler et al. (2014) Anderson M. Winkler, Gerard R. Ridgway, Matthew A. Webster, Stephen M. Smith, and Thomas E. Nichols.  Permutation inference for the general linear model.  _NeuroImage_ , 92:381–397, 2014.  ISSN 10959572.  doi: 10.1016/j.neuroimage.2014.01.060.

## References

  * Alberton et al. (2020) Bianca AV Alberton, Thomas E Nichols, Humberto R Gamba, and Anderson M Winkler.  Multiple testing correction over contrasts for brain imaging.  _NeuroImage_ , 216:116760, 2020.
  * Andersen et al. (1999) Anders H. Andersen, Don M. Gash, and Malcolm J. Avison.  Principal component analysis of the dynamic response measured by fMRI: A generalized linear systems framework.  _Magnetic Resonance Imaging_ , 17(6):795–815, 1999.  ISSN 0730725X.  doi: 10.1016/S0730-725X(99)00028-4.
  * Andreella et al. (2020) Angela Andreella, Jesse Hemerik, Wouter Weeda, Livio Finos, and Jelle Goeman.  Permutation-based true discovery proportions for fmri cluster analysis.  arXiv preprint arXiv:2012.00368, 2020.
  * Bahr et al. (2013) Timothy M Bahr et al.  Peripheral blood mononuclear cell gene expression in chronic obstructive pulmonary disease.  _American journal of respiratory cell and molecular biology_ , 49(2):316–323, 2013.
  * Barch et al. (2013) Deanna M Barch, Gregory C Burgess, Michael P Harms, Steven E Petersen, Bradley L Schlaggar, Maurizio Corbetta, Matthew F Glasser, Sandra Curtiss, Sachin Dixit, Cindy Feldt, et al.  Function in the human connectome: task-fmri and individual differences in behavior.  _Neuroimage_ , 80:169–189, 2013.
  * Benjamini and Hochberg (1995) Yoav Benjamini and Yosef Hochberg.  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  _Journal of the Royal statistical society: series B (Methodological)_ , 57(1):289–300, 1995.
  * Benjamini and Yekutieli (2001) Yoav Benjamini and Daniel Yekutieli.  The control of the false discovery rate in multiple testing under dependency.  _Annals of statistics_ , pages 1165–1188, 2001.
  * Bickel and Freedman (1981) Peter J Bickel and David A Freedman.  Some Asymptotic Theory for the Bootstrap.  _Annals of Statistics_ , 9(6):1196–1217, 1981\.  ISSN 00905364.  doi: 10.1214/aos/1176342871.
  * Blain et al. (2022) Alexandre Blain, Bertrand Thirion, and Pierre Neuvial.  Notip: Non-parametric true discovery proportion control for brain imaging.  _NeuroImage_ , page 119492, 2022.  URL https://hal.archives-ouvertes.fr/hal-03649114.
  * Blanchard et al. (2020) Gilles Blanchard, Pierre Neuvial, Etienne Roquain, et al.  Post hoc confidence bounds on false positives using reference families.  _Annals of Statistics_ , 48(3):1281–1303, 2020\.
  * Blanchard et al. (2021) Gilles Blanchard, Pierre Neuvial, and Etienne Roquain.  On agnostic post hoc approaches to false positive control.  In _Handbook of Multiple Comparisons_ , Handbooks of Modern Statistical Methods. Chapman & Hall/CRC, November 2021.  URL https://hal.archives-ouvertes.fr/hal-02320543.
  * Cui and Churchill (2003) Xiangqin Cui and Gary A Churchill.  Statistical tests for differential expression in cDNA microarray experiments.  _Genome Biol_ , 4(4):210, 2003.
  * Davenport and Nichols (2022) Samuel Davenport and Thomas E. Nichols.  The expected behaviour of random fields in high dimensions: contradictions in the results of [1].  _Magnetic Resonance Imaging_ , 2022.
  * Davenport et al. (2021) Samuel Davenport, Fabian Telschow, Thomas E. Nichols, and Armin Schwarzman.  Confidence regions for the location of peaks of a smooth random field.  2021\.
  * De Haan and Ferreira (2006) Laurens De Haan and Ana Ferreira.  _Extreme value theory: an introduction_ , volume 21.  Springer, 2006.
  * Ebrahimpoor et al. (2020) Mitra Ebrahimpoor, Pietro Spitali, Kristina Hettne, Roula Tsonaka, and Jelle Goeman.  Simultaneous enrichment analysis of all possible gene-sets: unifying self-contained and competitive methods.  _Briefings in bioinformatics_ , 21(4):1302–1312, 2020.
  * Eck (2018) Daniel J Eck.  Bootstrapping for multivariate linear regression models.  _Statistics & Probability Letters_, 134:141–149, 2018\.
  * Enjalbert-Courrech and Neuvial (2022) Nicolas Enjalbert-Courrech and Pierre Neuvial.  Powerful and interpretable control of false discoveries in differential expression studies.  bioRxiv preprint: https://doi.org/10.1101/2022.03.08.483449, 2022.
  * Freedman (1981) David A Freedman.  Bootstrapping regression models.  _The Annals of Statistics_ , 9(6):1218–1228, 1981\.
  * Genovese and Wasserman (2006) Christopher R Genovese and Larry Wasserman.  Exceedance control of the false discovery proportion.  _Journal of the American Statistical Association_ , 101(476):1408–1417, 2006.
  * Genovese et al. (2002) Christopher R Genovese, Nicole A Lazar, and Thomas E Nichols.  Thresholding of statistical maps in functional neuroimaging using the false discovery rate.  _Neuroimage_ , 15(4):870–878, 2002.
  * Goeman and Solari (2011) Jelle J. Goeman and Aldo Solari.  Multiple Testing for Exploratory Research.  _Statistical Science_ , 26(4):584–597, 2011.  ISSN 0883-4237.  doi: 10.1214/11-STS356.
  * Goeman et al. (2019) Jelle J Goeman, Rosa J Meijer, Thijmen JP Krebs, and Aldo Solari.  Simultaneous control of all false discovery proportions in large-scale multiple hypothesis testing.  _Biometrika_ , 106(4):841–856, 2019.
  * Hemerik et al. (2019) Jesse Hemerik, Aldo Solari, and Jelle J Goeman.  Permutation-based simultaneous confidence bounds for the false discovery proportion.  _Biometrika_ , 106(3):635–649, 2019.
  * Henson (2006) R. Henson.  Forward inference using functional neuroimaging: Dissociations versus associations.  _Trends in cognitive sciences_ , 10:64–69, 2006.
  * Holm (1979) Sture Holm.  A simple sequentially rejective multiple test procedure.  _Scand. J. Statist._ , 6(2):65–70, 1979.  ISSN 0303-6898.
  * Hommel (1988) Gerhard Hommel.  A stagewise rejective multiple test procedure based on a modified bonferroni test.  _Biometrika_ , 75(2):383–386, 1988.
  * Korn et al. (2004) Edward L Korn, James F Troendle, Lisa M McShane, and Richard Simon.  Controlling the number of false discoveries: application to high-dimensional genomic data.  _Journal of Statistical Planning and Inference_ , 124(2):379–398, 2004.
  * Meinshausen (2006) Nicolai Meinshausen.  False discovery control for multiple tests of association under general dependence.  _Scandinavian Journal of Statistics_ , 33(2):227–237, 2006.
  * Neuvial (2020) Pierre Neuvial.  _Contributions to statistical inference from genomic data_.  Habilitation thesis, Université Toulouse III (France), 2020.  Available from https://tel.archives-ouvertes.fr/tel-02969229.
  * Romano and Wolf (2005) Joseph P Romano and Michael Wolf.  Exact and Approximate Stepdown Methods for Multiple Hypothesis Testing.  _Journal of the American Statistical Association_ , 100(469):94–108, 2005.  ISSN 0162-1459.  doi: 10.1198/016214504000000539.
  * Rosenblatt et al. (2018) Jonathan D Rosenblatt, Livio Finos, Wouter D Weeda, Aldo Solari, and Jelle J Goeman.  All-resolutions inference for brain imaging.  _Neuroimage_ , 181:786–796, 2018.
  * Sarkar et al. (2008) Sanat K Sarkar et al.  On the simes inequality and its generalization.  In _Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen_ , pages 231–242. Institute of Mathematical Statistics, 2008.
  * Simes (1986) R J Simes.  An improved Bonferroni procedure for multiple tests of significance.  _Biometrika_ , 73(3):751–754, 1986.
  * Smyth (2004) Gordon K. Smyth.  Linear models and empirical bayes methods for assessing differential expression in microarray experiments.  _Statistical Methods in Genetics and Molecular Biology_ , 3(3), 2004.  doi: 10.2202/1544-6115.1027.
  * Storey and Tibshirani (2003) John D Storey and Robert Tibshirani.  Statistical significance for genomewide studies.  _Proceedings of the National Academy of Sciences_ , 100(16):9440–9445, 2003.
  * Van der Vaart (2000) Aad W Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * van der Vaart (1998) A.W. van der Vaart.  _Asymptotic Statistics_.  1998\.
  * Westfall (2011) Peter H. Westfall.  On Using the Bootstrap for Multiple Comparisons.  _Journal of Biopharmaceutical Statistics_ , 21(6):1187–1205, 2011.  ISSN 1054-3406.  doi: 10.1080/10543406.2011.607751.  URL http://www.tandfonline.com/doi/abs/10.1080/10543406.2011.607751.
  * Winkler et al. (2014) Anderson M. Winkler, Gerard R. Ridgway, Matthew A. Webster, Stephen M. Smith, and Thomas E. Nichols.  Permutation inference for the general linear model.  _NeuroImage_ , 92:381–397, 2014.  ISSN 10959572.  doi: 10.1016/j.neuroimage.2014.01.060.