False Discovery Rate for Functional Data

TL;DR

This paper extends the false discovery rate concept to functional data analysis, introducing a continuous BH procedure and demonstrating its effectiveness through simulations and real-world satellite temperature data.

Contribution

It proposes a novel continuous FDR control method for functional data and provides theoretical conditions for its validity.

Findings

The functional BH procedure controls FDR under specified conditions.

Simulation studies show the method's effectiveness compared to existing approaches.

Application to satellite data identifies significant temperature increase regions.

Abstract

Since Benjamini and Hochberg introduced false discovery rate (FDR) in their seminal paper, this has become a very popular approach to the multiple comparisons problem. An increasingly popular topic within functional data analysis is local inference, i.e., the continuous statistical testing of a null hypothesis along the domain. The principal issue in this topic is the infinite amount of tested hypotheses, which can be seen as an extreme case of the multiple comparisons problem. In this paper we define and discuss the notion of false discovery rate in a very general functional data setting. Moreover, a continuous version of the Benjamini-Hochberg procedure is introduced along with a definition of adjusted p-value function. Some general conditions are stated, under which the functional Benjamini-Hochberg procedure provides control of the functional FDR. Two different simulation studies are presented; the first study has a one-dimensional domain and a comparison with another state of the art method, and the second study has a planar two-dimensional domain. Finally, the proposed method is applied to satellite measurements of Earth temperature. In detail, we aim at identifying the regions of the planet where temperature has significantly increased in the last decades. After adjustment, large areas are still significant.

Tables

Table 1. Table 1: Results from simulation at 5% threshold with 255x255 sample points. The varying number of replications is due to the upper limit of 2500 signals in total

$\begin{array}{c}\text{Setup}\\ \text{no.}\end{array}$	$\begin{array}{c}\text{Signal}\\ \text{size}\end{array}$	$\begin{array}{c}\text{\# of samples}\\ \text{per test}\end{array}$	Replications	Sensitivity	FPR	FDR	$\begin{array}{c}\text{FDR,}\\ \text{unadjusted}\end{array}$
1	2.0	20	125	0.498	0.00785	0.0252	0.113
2	2.0	10	250	0.226	0.00553	0.0257	0.135
3	2.0	40	62	0.654	0.00810	0.0242	0.117
4	1.0	20	125	0.116	0.00467	0.0249	0.155
5	0.5	20	125	0.0065	0.00295	0.0258	0.235

Table 2. Table 2: Areas of significance of the correction methods at various significance levels as percentage of Earth total

Significance level	Unadjusted p value	fBH-adjusted p value
0.10	0.410	0.229
0.05	0.324	0.150
0.01	0.178	0.062
0.001	0.074	0.023

Equations

$$\mathrm{E}[Q]=\mathrm{E}\left[\frac{\#\{i:G_i\text{ is true but rejected}\}}{\#\{i:G_i\text{ is rejected}\}}\right]$$

$$\mathrm{E}[Q]=\mathrm{E}\left[\frac{\sum _{i:G_i\text{ is true but rejected}}{w}_i}{\sum _{i:G_i\text{ is rejected}}{w}_i}\right]$$

$$k=\arg \underset{i}{\max }\left[p_{(i)}\le \frac{i}{m}\alpha \right]$$

$$k=\arg \underset{i}{\max }\left[p_{(i)}\le \sum _{j:p_j\le p(i)}{w}_j\alpha \right]$$

$${\tilde{p}}_{(i)}=\min (1,\frac{m}{i}{p}_{(i)},\dots ,\frac{m}{m-1}{p}_{(m-1)},p_{(m)})$$

$${\tilde{p}}_{(i)}=\min \left\{1,\frac{p_{(i)}}{{\sum }_{j:p_j\le p_{(i)}}{w}_j},\frac{p_{(i+1)}}{{\sum }_{j:p_j\le p_{(i+1)}}{w}_j},\dots ,\frac{p_{(m-1)}}{{\sum }_{j:p_j\le p_{(m-1)}}{w}_j},p_{(m)}\right\},i=1,\dots m$$

$$x\le y\Rightarrow P(\mathbf{X}\in D|X_i=x)\le P(\mathbf{X}\in D|X_i=y)$$

$$\mathrm{E}[Q]\le \frac{m_0}{m}\alpha \le \alpha$$

$$\mathrm{fFDR}=\mathrm{E}[Q]=\mathrm{E}\left[\frac{\nu (V)}{\nu (V\cup S)}{1}_{\nu (V\cup S)>0}\right]$$

$$p(t)\le {\alpha }^{\ast }\text{where}{\alpha }^{\ast }=\arg \underset{r}{\max }\frac{\nu (\{s:p(s)\le r\})}{\nu (\mathbf{T})}\ge {\alpha }^{-1}r$$

$$\tilde{p}(t)=\underset{s\ge p(t)}{\min }\left\{1,\frac{\nu (\mathbf{T})s}{\nu (r:p(r)\le s)}\right\},t\in \mathbf{T}$$

$$p(t)\le {\alpha }^{\ast }⟺\tilde{p}(t)\le \alpha ,t\in \mathbf{T}$$

$$\begin{gathered} p(t)\le {\alpha }^{\ast }\iff \exists s\ge p(t):\frac{\nu (\{r:p(r)\le s\})}{\nu (\mathbf{T})}\ge {\alpha }^{-1}s⟺ \\ \exists s\ge p(t):\frac{\nu (\mathbf{T})s}{\nu (\{r:p(r)\le s\})}\le \alpha \iff \tilde{p}(t)\le \alpha \end{gathered}$$

$$P\left[\underset{k\to \infty }{\lim }\underset{r}{\sup }\frac{\#(S_k\cap \{s:p(s)\le r\})}{\#S_k}-\mu \{s:p(s)\le r\}\to 0\right]=1$$

$$P\left[\underset{k\to \infty }{\lim }\underset{r}{\sup }\frac{\#(S_k\cap \{s:p(s)\le r\}\cap U)}{\#S_k}-\mu (\{s:p(s)\le r\}\cap U)\to 0\right]=1$$

$$\mu \{s:p(s)=t\}=0\forall t\in \mathbf{T}$$

$$Q_k=\frac{\#\{t\in S_k:p(t)\le b_k\}\cap U}{\#\{t\in S_k:p(t)\le b_k\}},b_k=\arg \underset{r}{\max }\frac{\#\{s\in S_k:p(s)\le r\}}{\#S_k}\ge {\alpha }^{-1}r$$

$$\underset{k\to \infty }{\lim }{Q}_k=Q.$$

$$\underset{k\to \infty }{\mathrm{limsup}}\mathrm{E}[Q_k]\le \alpha \frac{\mu (U)}{\mu (\mathbf{T})}\le \alpha$$

$$P\left[\underset{k\to \infty }{\lim }\underset{r}{\sup }\frac{{\sum }_{i\in S_k\cap \{s:p(s)\le r\}}f(i)}{\#S_k}-{\int }_{\{s:p(s)\le t\}}f(x)\mathrm{d}x\to 0\right]=1$$

$$P\left[\underset{k\to \infty }{\lim }\underset{r}{\sup }\frac{{\sum }_{i\in S_k\cap \{s:p(s)\le r\}\cap U}f(i)}{\#S_k}-{\int }_{\{s:p(s)\le t\}\cap U}f(x)\mathrm{d}x\to 0\right]=1$$

$$\nu \{s:p(s)=t\}=0\forall t\in \mathbf{T}$$

$$Q_k=\frac{{\sum }_{t\in S_k\cap U:p(t)\le b_k}f(t)}{{\sum }_{t\in S_k:p(t)\le b_k}f(t)},b_k=\arg \underset{r}{\max }\frac{{\sum }_{\{s\in S_k:p(s)\le r\}}f(s)}{{\sum }_{\{s\in S_k\}}f(s)}\ge {\alpha }^{-1}r$$

$$\underset{k\to \infty }{\lim }{Q}_k=Q.$$

$$\underset{k\to \infty }{\mathrm{limsup}}\mathrm{E}[Q_k]\le \alpha \frac{\nu (U)}{\nu (\mathbf{T})}\le \alpha$$

$$\nu \{s:p(s)=t\}=0\forall t\in \mathbf{T}$$

$$\mathrm{E}[Q]$$

$$=\underset{k\to \infty }{\mathrm{limsup}}\mathrm{E}[Q_k]\le \underset{k\to \infty }{\mathrm{limsup}}\frac{\#(S_k\cap U)}{\#S_k}\alpha =\alpha \frac{\nu (U)}{\nu (\mathbf{T})}\le \alpha ,$$

$$\#\{s\in [0,D]|p(s)=t\}\le N_C\forall t\in (0,1]$$

$$y_i(t)=\beta (t)x_i+{\epsilon }_i(t)i=1,\dots ,n,t\in (0,1)$$

$$H_t^0:\beta (t)=0;H_t^1:\beta (t)\ne 0.$$

Full text

In this paper we define and discuss the notion of false discovery rate in a very general functional data setting. Moreover, a continuous version of the Benjamini-Hochberg procedure is introduced along with a definition of adjusted p-value function. Some general conditions are stated, under which the functional Benjamini-Hochberg procedure provides control of the functional FDR. Two different simulation studies are presented; the first study has a one-dimensional domain and a comparison with another state of the art method, and the second study has a planar two-dimensional domain.

Finally, the proposed method is applied to satellite measurements of Earth temperature. In detail, we aim at identifying the regions of the planet where temperature has significantly increased in the last decades. After adjustment, large areas are still significant.

1 Department of Mathematical Sciences, University of Copenhagen

2 Department of Statistical Sciences, Università Cattolica di Sacro Cuore

3 Department of Mathematics, Politecnico di Milano

1 Introduction

Statistical inference, and in particular hypothesis testing, is central in the field of statistics. Functional data analysis (FDA) often deals with “summary statistics” such as identifying mean curves/trajectories and principal modes of variation, including clustering and classification. Yet hypothesis testing is still a key part of FDA, where it is often done in conjunction with functional regression. Indeed, it is important to assess whether there is a significant effect of a covariate on the response, where covariates and response variables can be functions or scalars, depending on the setup.

In the common case of a functional response and scalar covariates, functional regression is usually modelled through a linear model on the form $y_{i}=Ax_{i}+\text{noise}$, where $y_{i}$ is a curve belonging to a suitable function space e.g. $L^{2}[0,1]$, $x_{i}$ are covariate(s), and $A\in L$ is a linear operator which ought to be estimated. In the simple case of one covariate, the null hypothesis would be asking whether $A=0$, and more generally the null hypothesis would be asking if $A\in U$ for a given subspace $U\subset L$.

This is a hard question to ask for functional data and is not that frequently encountered. It is absent in some textbooks [13, 5] but has extensive treatment in [11]. There are various approaches to this issue, which importantly also depend on the scope of the test. We will distinguish between two kinds of tests, (i) global tests: does covariate $x$ have ‘influence’ on curve $\theta$ in at least one part of the domain of $\theta$, and (ii) local tests or domain selections: if covariate $x$ has an influence, which part(s) of the domain of $\theta$ are significantly affected?

Global tests

Global tests have been studied by various authors, for references see e.g. [11]. However, a crucial feature is that many of the available methods rely on (strong) parametric assumptions on the data distribution, such as Gaussianity111Note: This includes independence of PC scores (among other things), which may be valid asymptotically but can be problematic for the usually small sample sizes and infinite dimensions that are characteristic of functional data analysis.

Non-parametric approaches such as permutation tests are popular alternatives to parametric tests. Curves are (randomly) permuted wrt. likelihood-independent transformations, and some suitable test statistic (e.g. deviation from the mean) is evaluated for each permutation. However this is computationally expensive, and permutation tests are only asymptotically exact in the presence of several covariates.

Local tests

Local tests have not been studied to the same extent as global tests. In functional data analysis performing local inference carries several issues. The most important one is how to define, and how to control the probability of committing a type 1 error globally over the whole domain. One recent framework for doing local testing on functional data is the Interval-wise testing (IWT) introduced by [12] and extended to linear models by [1]. These procedures perform non-parametric inference based on permutation schemes and provide (asymptotical) control of the family-wise error rate (FWER) on each sub-interval of the domain. In detail, the probability of falsely selecting at least part of an interval where the null hypothesis is not violated is controlled. Another procedure proposed in the literature, and focusing on the FWER control is the Fmax-procedure [10, 17]. The Fmax-procedure is a method that provides strong control of the FWER, and like the IWT-framework it is based on permutation tests. The procedure is multivariate in nature, but it can also be applied to functional data discretized on a fine grid.

An alternative approach to local testing in functional or spatial data is to use discrete features of the observed functions such as local maxima or zero values. That is, instead of assessing a continuum of tests, one selects a finite number of data features for testing. Using discrete features has some challenge with respect to interpretation, as one has to specify when two different observations can be considered instances of the same discrete feature, and there are no obvious definitions of domain selection. [6, 14] present some interesting methods with this, where they also proof control of the false discovery rate (FDR).

In this work we likewise focus on the control of the FDR. In particular, we focus on defining a continuum of statistical tests over the domain of functional data, and secondly we propose a generalization of the Benjamini-Hochberg procedure [2] for adjusting such tests in order to control the FDR.

The remainder of the paper is organised as follows: Section 2 describes false discovery rate in the multivariate case and reviews related work. Section 3 presents the novel work of functional false discovery rate and the functional Benjamini-Hochberg procedure. Two different simulation studies are presented in Section 4, and in Section 5 we apply the proposed procedure to a data set on climate change. Finally in Section 6 we highlight and discuss important points of this article. Proof of the main theorems are provided in the appendix.

2 False discovery rate for multivariate data

Background

Multiple testing is a central topic in statistics. Observed within virtually every area of statistics, it is fundamental in many statistical applications and multiple testing is recognised as an important statistical issue within many sciences.

Many ways to deal with multiple testing have been proposed with various advantages and disadvantages – one popular approach is to use the False Discovery Rate (FDR) [2]. The false discovery rate is the expected proportion of false rejections (“discoveries”) among all among rejected hypotheses. The procedures controlling the FDR are generally more powerful than the ones controlling the FWER.

FDR is often applied in cases when a single or comparatively few false positives is not considered a serious issue, as long as their rate among all discoveries can be controlled. In [2] the Benjamini-Hochberg (BH) procedure for controlling FDR is introduced. In the succeeding literature, a number of other procedures for controlling FDR have been proposed (see [9] for a discussion). The paper [4] is important, as it introduces a modification of the BH procedure that controls FDR without specifying any dependency assumptions. More importantly, they show that the original procedure introduced in [2] controls FDR under a weaker assumption than independence, namely positive regression dependence on the subset of true null hypotheses (PDRS).

Other closely related quantities for assessing the errors within the paradigm of multiple testing have been proposed such as the weighted false discovery rate (WDFR) [3], (see below) the positive false discovery rate [15], and the local false discovery rate [7].

However, in this paper we will only focus on the BH procedure and FDR which – due to its simple interpretation and definition – is still the most popular method for multiplicity correction. False Discovery Rate is only defined for finite numbers of hypotheses, and one must be careful when defining FDR on infinite sets of hypotheses.

False discovery rate

Assume we are given a set of $m$ null hypotheses, $G_{1},\dots,G_{m}$, each of which can either be true or false, and can either be accepted or rejected by a statistical test. Furthermore, let $w_{1},\dots,w_{m}$ be strictly positive weights with $\sum w_{i}=1$, which we assume are a priori known. This will be used in the case of weighted false discovery rate. These weights can e.g. be interpreted as how important the different tests are, where the ”usual” false discovery rate corresponds to the case of equal weights.

Definition 2.1** (False discovery rate, unweighted case).**

The false discovery rate is defined as:

[TABLE]

with $Q:=0$ whenever the denominator is zero.

Definition 2.2** (False discovery rate, weighted case).**

The false discovery rate in the weighted case is defined as:

[TABLE]

with $Q:=0$ whenever the denominator is zero.

The Benjamini-Hochberg procedure

Let $\{p_{(i)}\}_{i=1}^{m}$ be the $p$-values sorted in increasing order. Let $\{G_{(i)}\}_{i=1}^{m}$ be the corresponding ordering of the hypotheses and $\{w_{(i)}\}_{i=1}^{m}$ the corresponding ordering of weights. The very popular and easily applicable Benjamini-Hochberg (BH) procedure for multiple comparison adjustment [2, 3] is defined as follows:

Definition 2.3** (Benjamini-Hochberg procedure, unweighted case).**

Define

[TABLE]

The Benjamini-Hochberg procedure is: reject hypotheses $G_{(1)},\dots,G_{(k)}$ corresponding to the $k$ smallest $p$-values and accept the rest.

Definition 2.4** (Benjamini-Hochberg procedure, weighted case).**

Define

[TABLE]

The weighted Benjamini-Hochberg procedure is: reject hypotheses $G_{(1)},\dots,G_{(k)}$ corresponding to the $k$ smallest $p$-values and accept the rest.

Unlike most adjustment procedures introduced prior to this, the BH procedure is scalable: if data is duplicated such that one has twice the amount of hypotheses, the inference by using the BH procedure will still be the same.

Adjusted p-values

A concept often used in context of multiple testing is adjusted (or corrected) p-values. Informally, the adjusted p-values for a multiple testing procedure are defined as corrections $\{\tilde{p}_{i}\}$ of the original $p$-values such that a null hypothesis $G_{i}$ can be rejected at level $\alpha$ if $\tilde{p}_{i}\leq\alpha$.

The adjusted p-values for the unweighted Benjamini-Hochberg procedure are

[TABLE]

where $p_{(\cdot)}$ and $\tilde{p}_{(\cdot)}$ are the order statistics of $p$ and $\tilde{p}$, respectively. By construction, the ordering is the same.

We can likewise define adjusted p-values for the weighted Benjamini-Hochberg procedure:

[TABLE]

PRDS and control of false discovery rate

Benjamini and Hochberg showed in their seminal paper [2] that if test statistics for different hypotheses are independent, then the false discovery rate is controlled by $\frac{m_{0}}{m}\alpha$ where $m_{0}$ is the total number of correct null hypotheses. The independence assumption was later relaxed by [4] to positive regression dependency on one (PRDS).

Below we define the PRDS property and extend it to the infinite-dimensional case, which will be needed later.

Definition 2.5** (Positive regression dependency on one (PRDS)).**

Let ’$\leq$’ be the usual ordering on $\mathbb{R}^{l}$. An increasing set $D\subseteq\mathbb{R}^{l}$ is a set satisfying $x\in D\wedge y\geq x\Rightarrow y\in D$.

A random variable $\mathbf{X}$ on $\mathbb{R}^{l}$ is said to be PRDS on $I_{0}$, where $I_{0}$ is a subset of $\{1,\dots,l\}$, if it for any increasing set $D$ and $i\in I_{0}$ holds that

[TABLE]

Let $\mathbf{Z}$ be an infinite-dimensional random variable, where instances of $\mathbf{Z}$ are functions $T\rightarrow\mathbb{R}$. We say that $\mathbf{Z}$ is PRDS on $U\subseteq T$ if all finite-dimensional distributions of $\mathbf{Z}$ are PRDS. That is, for all finite subsets $I=\{i_{1},\dots,i_{l}\}\subseteq T$, it holds that $Z(i_{1}),\dots,Z(i_{l})$ is PDRS on $I\cap U$.

We refer to [4] for a discussion on the PRDS property and how it relates to other types of dependency.

Theorem 2.6**.**

Given a set of hypotheses $\{H_{1},\dots,H_{m}\}$ and corresponding p-values $(p_{1},\dots,p_{m})$, let $I_{0}=\{i_{1},\dots,i_{m_{0}}\}\subseteq\{1,\dots,m\}$ be the index set corresponding to true null hypotheses $\{H_{i_{1}},\dots,H_{i_{m_{0}}}\}$.

If the joint distribution of the p-values $(p_{1},\dots,p_{m})$ is PDRS on $I_{0}$, the BH procedure controls the FDR at level $\frac{m_{0}}{m}\alpha$ in the sense that

[TABLE]

where $Q$ is the proportion of false discoveries.

Proof.

See [4, Theorem 1.2] ∎

3 False discovery rate for functional data

In this section we define the false discovery rate for functional data and propose a functional extension of the Benjamini-Hochberg procedure. These are the functional versions of the classical false discovery rate and the BH procedure used in the multivariate cases.

3.1 Definition of functional false discovery rate

Our definition of false discovery rate is related to that of [16], although [16] uses a seemingly arbitrary lower bound for the measure of rejection region, not present in other papers on FDR, and avoids references to p-values which are natural in the context of multiple testing.

Having defined functional false discovery rate and introduced the functional Benjamini-Hochberg procedure, we proceed to define the adjusted p-value function, an important tool for practical applications, in Section 3.3. Section 3.4 contains a discussion and results about finite approximations of the fBH procedure, which is easily adopted into an algorithm. Finally we prove control of the false discovery rate under regularity conditions in Section 3.5.

For the remainder of this section we assume that we have $N$ functional samples, $y_{1},\dots,y_{N}$, with $y_{i}:\mathbf{T}\rightarrow\mathbb{R}$, where $\mathbf{T}\subset R^{d}$ is an open and bounded subset of $\mathbb{R}^{d}$.

Suppose that for each point $t\in\mathbf{T}$, we have a null hypothesis $H^{0}_{t}$ together with an alternative hypothesis $H^{A}_{t}$, that we are interesting in testing. Furthermore, suppose that by pointwise application of some statistical test we obtain p-values $p(t)$ for every $t$ with the property that $(H^{0}_{t}$ true$)\implies p(t)\sim U(0,1)$, where $U(0,1)$ is the uniform distribution on $(0,1)$. The $p$-values together make up a function $p:\mathrm{\mathbf{T}}\rightarrow[0,1]$, the unadjusted $p$-value function [12].

Let $U$ be the set of the domain where the null hypothesis is true, ie. $U=\{t\in\mathrm{\mathbf{T}}:H_{0}^{t}$ is true $\}$. Let $\nu$ be a bounded measure on $\mathrm{\mathbf{T}}$ that is absolutely continuous wrt. the Lebesgue measure, which we denote by $\mu$. By absolute continuity we have $\nu=f\cdot\mu$ for some measurable function $f:\mathrm{\mathbf{T}}\rightarrow[0,\infty)$. The function $f$ can be interpreted as a weight function assigning more weight to some regions of $\mathrm{\mathbf{T}}$ than others, and is the functional counterpart of the weights used in section 2.

Given $U$ and an instance of $p(t)$, let $V=\{t:H^{0}_{t}$ is true and $H^{0}_{t}$ is rejected$\}$ be the region where the null hypothesis is wrongly rejected, and let $S=\{t:H^{0}_{t}$ is false and $H^{0}_{t}$ is rejected$\}$ be the region where the null hypothesis is correctly rejected. The set $V$ corresponds to committing type I errors, and in a given research situation, it is desired that $V$ is as small as possible and $S$ is as large as possible. Specifically, the functional false discovery rate can be defined as follows.

Definition 3.1** (Functional false discovery rate).**

Define the functional false discovery rate (fFDR) as

[TABLE]

where $Q$ is the proportion of false discoveries.

Remark 3.2** (False discovery rate for other manifolds).**

In this paper we define false discovery rate for functional data defined on open subsets of $\mathbb{R}^{d}$. However, many smooth manifolds can be diffeomorphically mapped into open and bounded subsets of $\mathbb{R}^{d}$. The mapping gives naturally rise to a measure on this set, which can be used as measure for functional false discovery rate.

3.2 The functional Benjamini-Hochberg procedure: the adjusted threshold

Analogous to the multivariate case, we can define the fBH procedure, that is the Benjamini-Hochberg procedure for functional data. The counts and sums used in Defintions 2.3 and 2.4 are replaced by the measure $\nu$, that incorporates any weighting of the domain $\mathrm{\mathbf{T}}$.

The main theoretical result of this article is that the fBH procedure can be approximated by the multivariate BH procedure on the pointwise evaluations of functional data, and that it controls the fFDR by $\alpha\nu(U)/\nu(\mathrm{\mathbf{T}})$ under regularity conditions.

Definition 3.3** (Functional Benjamini-Hochberg procedure).**

Let $\alpha\in(0,1)$ be a desired significance level for the tests. The functional Benjamini-Hochberg (fBH) procedure is:

Reject hypotheses $H^{0}_{t}$ that satisfy**

[TABLE]

We will refer to $\alpha^{*}$ as the adjusted threshold of the procedure, and the function $a(r)=\nu(\{s:p(s)\leq r\})$ as the cumulated p-value function.

Two examples of the functional BH procedure are shown in Figure 1.

3.3 The functional Benjamini-Hochberg procedure: the adjusted p-value function

In the application of FDR-controlling procedures, it is often of interest to have the possibility of changing the threshold. For the functional case, this amounts to asking for which $\alpha$ that $H_{0}^{t}$ is rejected after adjustment by the fBH procedure. This naturally leads to the notion of the fFDR-adjusted p-value function:

[TABLE]

This is analogous to adjusted p-values, which applies in the discrete case, and plays a similar role, ie. if $\tilde{p}(t)\leq\alpha$ this means that $H^{0}_{t}$ will be rejected when the (weighted) fBH procedure is applied with threshold level $\alpha$.

Adjusted $p$-values allow us to quantify significance after adjustment and to simultaneously compute rejection areas for all values of $\alpha$. By 3.4 below, adjustment using either the adjusted threshold or adjusted p-value function gives rise to the same inference.

Proposition 3.4**.**

Let $\alpha$ be a given threshold, and define the adjusted threshold $\alpha^{*}$ according to 3.3. Then adjustment using the adjusted threshold or the adjusted p-value function gives the identical results, ie.

[TABLE]

Proof.

Using the definitions of the adjusted threshold and the adjusted p-value function, respectively, the result is true since

[TABLE]

∎

3.4 The functional Benjamini-Hochberg procedure: Finite approximation and computational cost

In Section 3.5 (see Theorem 3.9), we will prove that the fBH procedure controls the fFDR. Before doing so, we start defining a finite approximation of the fBH procedure, that is proven to converge to the theoretical one. Indeed, a crucial issue regarding the fBH procedure is that it makes use of an infinite amount of hypotheses and p-values, which in principle is computationally unattainable. We present here a simple algorithm for approximating the fBH procedure. In Propositions 3.5 and 3.7 we show that under regularity conditions, this algorithm approaches the fBH procedure in the limit.

We present the algorithm and proposition in two versions: a notationally simpler version with equal weights, which is a special case of the more general, weighted version.

3.4.1 Unweighted case

In many applications there is no a priori reason for assigning different weights to different parts of the domain, $\mathrm{\mathbf{T}}$. Thus we will assume that $\nu=\mu$ corresponding to the weight function $f$ being constant.

Algorithm** (Unweighted case).**

*Let $S$ be a dense, uniform grid in $\mathrm{\mathbf{T}}$. The fBH procedure can be approximated by evaluation of the ‘usual’ BH procedure:

Apply the BH procedure to $\{p(t):t\in S\}$.*

Obviously, not all grids are suitable, and what a ’dense uniform grid’ constitutes, is stated in Proposition 3.5. The typical choice of a grid would be a lattice, and by decreasing the mesh of the lattice, the approximation becomes better and better.

Proposition 3.5**.**

Let $\{S_{k}\}_{k=1}^{\infty},S_{1}\subset S_{2}\subset\dots$ be a dense, uniform grid in $\mathrm{\mathbf{T}}$ in sense that $S_{k}$ uniformly approximates all level sets of $p$ and $p|_{U}$ with probability one:

[TABLE]

and

[TABLE]

Furthermore, assume that the assumptions about the $p$-value function below hold true with probability one:

1. (a1)

All level sets of $p$ have zero measure,

[TABLE] 2. (a2)

$\alpha^{*}\in(0,\alpha]\Rightarrow$: for any open neighbourhood $O$ around $\alpha^{*}$ there exists $s_{1},s_{2}\in O$ s.t. $a(s_{1})>\alpha^{-1}s_{1},a(s_{2})<\alpha^{-1}s_{2}$, where $a$ is the cumulated p-value function (Definition 3.3). 3. (a3)

$[\alpha^{*}=0]\Rightarrow\min p(t)>0$.

Now define the k’th step false discovery proportion $Q_{k}$ by applying the (usual) BH procedure at level $\alpha$ to $p$ evaluated in $S_{k}$.

Mathematically, this can be defined by

[TABLE]

Then $Q_{k}$ behaves asymptotically as the functional false discovery proportion:

[TABLE]

where $Q$ is defined as in (3.1). Furthermore, if $p$ is PDRS wrt. the set of true null hypotheses with probability one, then the false discovery rate $\mathrm{E}[Q_{k}]$ is controlled by $\alpha\mu(U)/\mu(\mathrm{\mathbf{T}})$:

[TABLE]

Proof.

See appendix. Also a special case of 3.7. ∎

From the proof of the theorem, we have the following important corollary which states that as the grid $S_{k}$ becomes tighter and tighter, hypotheses are eventually rejected or accepted:

Corollary 3.6**.**

For $t\in\cup_{m=1}^{\infty}S_{m}$ and $k\geq 1$, define $H_{t,k}=(t\in S_{k})\wedge(p(t)\leq b_{k})$ where $b_{k}$ is is given by (3.4). That is, $H_{t,k}$ is true if the adjusted threshold at step $k$ is larger than $p(t)$. Assume $p(t)\neq\alpha$. Eventually, $H_{t,k}$ is either rejected or accepted.

Proof.

Proposition A.3 ∎

3.4.2 Weighted case

Allowing the weighting function, $f$, to vary allows for more general measures, such as those which arise from Remark 3.2. However, the notation is more tedious, so we have separated unweighted case as a special case.

Let $\nu=f\cdot\mu$, and assume that $f$ is a bounded and strictly positive density function $\mathrm{\mathbf{T}}\rightarrow\mathbb{R}$.

Algorithm** (Weighted case).**

*Let $S$ be a dense, uniform grid in $\mathrm{\mathbf{T}}$. The fBH procedure can be approximated by evaluation of the ’usual’ BH procedure:

Apply the BH procedure to $\{p(t):t\in S\}$ with weights $(f(t):t\in S)$ where the weights have been normalised to one.*

Under almost similar assumptions to the unweighted case, we are able to control the false discovery rate at level $\alpha$:

Proposition 3.7**.**

Let $\{S_{k}\}_{k=1}^{\infty},S_{1}\subset S_{2}\subset\dots$ be a dense, uniform grid in $\mathrm{\mathbf{T}}$ in sense that $S_{k}$ weighted by $f$ uniformly approximates all level sets of $p$ and $p|_{U}$ with probability one:

[TABLE]

and

[TABLE]

Furthermore, assume that the assumptions about the $p$-value function below hold true with probability one:

1. (a1)

All level sets of $p$ have zero measure,

[TABLE] 2. (a2)

$\alpha^{*}\in(0,\alpha]\Rightarrow$: for any open neighbourhood $O$ around $\alpha^{*}$ there exists $s_{1},s_{2}\in O$ s.t. $a(s_{1})>\alpha^{-1}s_{1},a(s_{2})<\alpha^{-1}s_{2}$, where $a$ is the cumulated p-value function (Definition 3.3). 3. (a3)

$[\alpha^{*}=0]\Rightarrow\min p(t)>0$.

Now define the $k$th step false discovery proportion $Q_{k}$ by applying the (usual) BH procedure at level $\alpha$ to $p$ evaluated in $S_{k}$, weighted by $f$ evaluated in $S_{k}$.

Mathematically, this can be defined by

[TABLE]

Then $Q_{k}$ behaves asymptotically as the functional false discovery proportion:

[TABLE]

where $Q$ is defined as in (3.1). Furthermore, if $p$ is PDRS wrt. the set of true null hypotheses with probability one, then the false discovery rate $\mathrm{E}[Q_{k}]$ is controlled by $\alpha\nu(U)/\nu(\mathrm{\mathbf{T}})$:

[TABLE]

Proof.

See appendix. ∎

Analogous to 3.6 we have the following corollary which states that as the grid becomes tighter and tighter, hypotheses are eventually rejected or accepted:

Corollary 3.8**.**

For $t\in\cup_{m=1}^{\infty}S_{m}$ and $k\geq 1$, define $H_{t,k}=(t\in S_{k})\wedge(p(t)\leq b_{k})$ where $b_{k}$ is given by (3.7). That is, $H_{t}$ is true if the adjusted threshold at step $k$ is larger than $p(t)$. Assume $p(t)\neq\alpha$. Eventually, $H_{t,k}$ is either rejected or accepted.

The results in this section are important as they outline how to do inference using the fBH procedure. The implementation is very easy and requires no additional tools besides an implementation of the multivariate BH procedure.

Computational cost

The fBH procedure adds $O(n\log n)$ computational cost to the calculation of adjusted p-values, where $n=\#S_{k}$ is the number of approximation points,. This extra computational cost is due to the sorting of $p$-values, which has $O(n\log n)$ computational cost. As a default, calculation of pointwise $p$-values has complexity $O(n)$, thus the total computational cost is $O(n\log n)$.

However, sorting algorithms on modern computer are very fast, whereas calculating pointwise p-values can be comparatively slow, in particular if permutation tests are used, such as in Section 4.1. Thus, the computational cost from fBH adjustment procedure is expected to be negligible in practice. This was also the case for the simulation studies and the case study presented in Sections 4 and 5.

3.5 Control of false discovery rate for functional data

With 3.7 in place, it is easy to prove that the functional Benjamini-Hochberg procedure controls FDR under PRDS and regularity assumptions.

Theorem 3.9**.**

Assume that there exists a uniform grid $\{S_{k}\}_{k=1}^{\infty},S_{1}\subset S_{2}\subset\dots$ dense in $\mathrm{\mathbf{T}}$, such that with probability one assumptions (3.5) and (3.6) are met. Furthermore, assume that $p(t)$ is PRDS wrt. the set of true null hypotheses with probability one, and that the assumptions about $p$-value function below hold true with probability one:

1. (a1)

All level sets of $p$ have zero measure,

[TABLE] 2. (a2)

$\alpha^{*}\in(0,\alpha]\Rightarrow$: for any open neighbourhood $O$ around $\alpha^{*}$ there exists $s_{1},s_{2}\in O$ s.t. $a(s_{1})>\alpha^{-1}s_{1},a(s_{2})<\alpha^{-1}s_{2}$, where $a$ is the cumulated p-value function (Definition 3.3). 3. (a3)

$[\alpha^{*}=0]\Rightarrow\min p(t)>0$.

Then the functional BH procedure controls FDR at level $\alpha\nu(U)/\nu(\mathrm{\mathbf{T}})$, ie. $\mathrm{E}[Q]\leq\alpha\nu(U)/\nu(\mathrm{\mathbf{T}})\leq\alpha$, when applying the functional BH procedure at level $\alpha$.

Note that the assumptions (3.5) and (3.6) are much simplified in the equal-weight case where $f\equiv 1$.

Proof.

For the ease of presentation, we will only show the proof in the equal-weight case, ie. when $f$ is constant (and thus can be disregarded).

From 3.7 and the notation of that proposition, it shows that $Q_{k}\rightarrow Q$ almost surely.

Since the PDRS assumption is fullfilled for all finite sets, we can apply Benjamini and Yekuteli’s original proposition to $S_{k}$ (Theorem 2.6), from which it holds that $E[Q_{k}]\leq\frac{\#(S_{k}\cap U)}{\#S_{k}}\alpha$.

As $0\leq Q_{k}\leq 1$ for all $k$, it is now a simple application of the dominated convergence theorem to show that $\mathrm{E}[Q]\leq\alpha\nu(U)/\nu(\mathrm{\mathbf{T}})$:

[TABLE]

where the last equation follows from assumption (3.6) (cf. its unweigthed form (3.3)). ∎

As remarked in [4] the PDRS assumption is sufficient but not necessary, and (2.1) needs only to be true for certain sets defined in relation to the order statistics of $p$. The details are quite technical, and we refer to [4, Remark 4.2] and the general discussion of that paper.

Sufficient criteria for the one-dimensional case

The assumptions (a1)-(a3) and equations (3.5) and (3.6) are consequences of the more simple criteria in the one-dimensional case, which must be true with probability one:

1. (d1)

$p$ is continuous. 2. (d2)

There exists a maximal number of crossings, $N_{C}$, i.e.:

[TABLE] 3. (d3)

$U$ is a finite union of disjoint intervals.

These criteria will generally be true for smooth curves as we typically use to model functional data.

4 Simulations

In this section two different simulation studies are performed, which differ both in scope and setting – in the first simulation study, a functional-on-scalar regression is performed using permutation tests, while in the second simulation we test for mean equal to zero using one-sided t-tests.

The first simulation study has a more theoretical flavour and compares our proposed method to the Fmax method. The second study is intended to simulate a more realistic scenario with a number of disjoint peaks, and the performance of our method for various levels of $\alpha$ is analysed.

4.1 1D-simulation

Description of simulation

In this section we wanted to numerically assess the performance of our fBH procedure, and to compare it with the Fmax-procedure [10, 17]. The Fmax-procedure is a method that provides strong control of the family-wise error rate in a multivariate high-dimensional setting. It can be extended to functional data by applying it to the discrete point-wise evaluations of the functional data.

We simulated functional data according to the following functional-on-scalar linear model:

[TABLE]

where $n=10$, $x_{i}=\frac{i-1}{n-1}$, and $\beta(t)=d\cdot f(t)$, with $d$ ranging from 0 to 5. We modelled the function $f(t)$ with a cubic B-spline expansion with 40 basis functions and equally-spaced knots. The first $h$ coefficients of the expansion were set to one, and the last $40-h$ coefficients were set to zero. The resulting function assumed the value 1 in the first part of the domain, 0 in the second part of the domain, with a smooth transition. We explored three values of the parameter $h$: $h\in\{10,20,30\}$.

The error functions $\varepsilon_{i}(t)$ were obtained by simulating the coefficients of the same cubic B-spline expansion. The 40 coefficients were sampled independently from a standard normal distribution. The three panels in the first column of Figure 2 show an instance of the simulated functional data with $d=5$, and $h=10,20,30$, respectively. The functional data are coloured in grayscale that is proportional to the value of $x_{i}$. The study is much similar to the simulation study presented in [1], but here it is of interest to vary the domain where the null hypothesis is true.

The fBH and Fmax procedures are applied to test the following hypotheses:

[TABLE]

The unadjusted $p$-value at point $t$ was computed with a permutation test based on the Freedman and Lane method [8].

With $d\in\{0,1,\ldots,5\}$ and $h\in\{10,20,30\}$, we had 18 scenarios in total, but only 16 different ones, as the scenarios were identical for $d=0$.

Simulation results

We measured the performances of the two methods by evaluating the FWER, FDR, false positivity rate (i.e., the measure of the incorrectly rejected part of the domain over the measure of the domain where the null hypothesis is true) and sensitivity (i.e., the measure of the correctly rejected part of the domain relative to the total measure of the domain where the null hypothesis is false). We performed the tests at nominal level $\alpha=0.05$.

Figure 2 reports the results of the simulation obtained averaging over 1000 instances. Each row of the figure report the results with a different value of $h$. Each panel on columns 2-5 reports one of the measures discussed before as a function of the parameter $d$ for the unadjusted $p$-value (black), the fBH procedure (dark grey) and the Fmax procedure (light grey).

As expected by theory, the unadjusted $p$-value function does not control the FWER nor the FDR. It controls instead the false positive rate. The Fmax method controls the FWER, and by consequence, also the FDR. Finally, the fBH method controls the FWER only weakly (i.e., when $d=0$ and by consequence $H_{t}^{0}$ is true for all $t$). It controls instead the FDR in all scenarios. We notice a trade-off between the type of control and the sensitivity of the methods. The Fmax, being provided with a stronger control, is also the less sensible to deviations from the null hypothesis. Instead, the fBH is provided with a less strong control, but it is more sensible.

4.2 2D-simulation

Description of simulation

The base signal $\theta$ for this simulation consisted of 9 conical spikes with height $h=1$ and diameter $d=0.2$ arranged on the grid $\{0.25,0.50,0.75\}^{2}$ with the unit square $\mathbf{T}=(0,1)^{2}$ as doma in. Five spikes had positives values and the remaining four spikes had negative values. On top we added an error signal generated as a smooth Matérn field $x_{i}$ with varying scale parameter. The observed signal is $y_{i}=\theta+x_{i}$ for $i=1,\dots,N$. See Figure 3 for an illustration.

The simulation was inspired by [6], who also studies FDR in a setting of random fields, albeit with a quite different scope. We tested against the pointwise null hypothesis $H_{0}(t):\theta(t)=0$ using a pointwise, one-sided t-test over a fine lattice.

The Benjamini-Hochberg procedure was applied, and false positive rate (FPR), false discovery rate and sensitivity values were evaluated, along with false discovery rate for the unadjusted p-values for comparison. We simulated 2500 samples of the error process at grid $255\times 255$. We assessed performance of the test in five different setups with varying strengths of the base signal and varying numbers of samples per t-test. Samples of the error process were recycled and used for all experimental setups, this is also the reason for the varying number of replications in table 1.

We remark that the observed data will show a ”bend” at the edges of the spikes, making these easily detectable by other methods. However, the use of pointwise tests will ignore such features.

Simulation results

We tested at various thresholds, $\alpha\in\{0.001,0.01,0.02,0.03,0.04,0.05,0.10\}$ for five different experimental setups described in Table 1; results are shown in Figure 4. The sensitivity values and false postive rate (FPR) are defined as proportions of rejected/accepted hypotheses to the total number of false and correct hypotheses, respectively, after p-value adjustment. $\mu(U)=71.7\%$ of the signal was zero, thus we would expect FDR to be controlled by $0.717\alpha$; for the 5% threshold of table 1, this is roughly $0.036$.

Unsurprisingly, sensitivity (power level), FPR and FDR increase with $\alpha$. Experiment 5 has almost no power, even at $\alpha=0.10$, but still a comparatively large FDR, indicating that the Benjamini-Hochberg procedure is not too conservative in this setup. The false discovery rate is remarkably stable across experimental setups, and shows a linear tendency. As expected, FDR is well below $\alpha$ in all instances, and also below $\mu(U)\alpha$. In comparison, the FDR for the unadjusted p-values exceeds $\alpha$ in all setups except one, and varies considerably with the experimental setting.

5 Application: Analysis of Climate Data

Climate change is a huge issue, both politically and scientifically. The main issue are increasing temperatures with many adverse effects on weather and climate. Knowing that temperature has increased significantly on a global scale, we wanted to test where on Earth temperature has increased.

5.1 Data and model

Data consists of yearly averages of temperatures, starting in 1983 and ending in 2007, for each $1^{\circ}\times 1^{\circ}$ tile on Earth, using standard latitudes and longitudes. Temperatures are satellite measurements collected by NASA. These data were obtained from the NASA Langley Research Center Atmospheric Science Data Center Surface meteorological and Solar Energy (SSE) web portal supported by the NASA LaRC POWER Project. 222http://eosweb.larc.nasa.gov

One crucial feature is that data was more densely sampled closer to the poles than close to Equator. Naturally, data also exhibits behaviour depending on local geography. Rather than viewing data as truly areal, we considered data to sit on the midpoints of the tiles, e.g. $(73.5^{\circ}S,24.5^{\circ}W)$ corresponds to the tile $(74^{\circ}S,73^{\circ}S)\times(24^{\circ}W,25^{\circ}W)$.

We applied the linear regression model $y_{st}=a_{s}+b_{s}\text{year}_{t}+\epsilon_{st},s\in S^{2},t\in\{1983,\dots,2007\}$, testing for positive trend, i.e. $H^{0}_{s}:b_{s}=0$ with alternative hypothesis $H^{A}_{s}:b_{s}>0$. Figure 5 displays the temperature changes, values of the pointwise t-statistics and corresponding (unadjusted) p-values. Observe how much the test statistic varies across the globe, and how differences between land and ocean are more visible in the lower the plot compared to the upper plot.

To perform the BH procedure, we mapped the sphere into $\mathbf{T}=(-\pi,\pi)\times(-\pi/2,\pi/2)$ by (scaled) polar coordinates, ie. longitude and latitude. This mapping gives rise to a measure $\nu=f\cdot\mu$ on $\mathbf{T}$ where $f$ is proportional to $\cos($latitude) cf. Remark 3.2. This measure gives uniform weights to all points on Earth, assuming Earth to be a perfect sphere. One-sided t-tests were used for obtaining unadjusted p-values. We used the same grid as the observations for approximating the BH procedure, in total $180\times 360=64800$ grid points.

5.2 Results

The coverage areas at various significance levels are provided in Table 2.

Although more conservative by construction than unadjusted p-values, the fBH-adjusted p-values still retained large significant areas, indicating that the temperature increase observed in these areas is very unlikely to be a coincidence in but a small fraction of these areas. If we take look at the map, the North Atlantic Ocean and northern China stands out; it is evident that these regions have experienced temperatures far above the normal in the latest years with the adverse weather effects this may cause.

6 Discussion

In this paper we defined the false discovery rate for functional data defined on generic subdomains of $\mathbb{R}^{k}$, and by remark 3.2, this is easily extended to non-euclidean domains such as the sphere used in the data application. Furthermore, we devised a correction method for controlling FDR, the functional BH procedure, along with an adjusted p-value function. We showed a finite dimensional approximation of the fBH procedure, and proved that it converges to its theoretical infinite-dimensional definition. The fBH procedure was successfully applied in two simulations and on a data set on climate change.

The two simulation studies allowed us to study the fBH procedure in a setting where the distributions and true null hypotheses were known. The false discovery rate was controlled by FBH procedure in all instances, unlike the unadjusted p-value functions. The signal-to-noise ratio proved to be important: the sensitivity of the fBH method increased with signal strength, but the false discovery rate was remarkably constant.

By means of the climate change analysis, we demonstrated the applicability of the method and also gained insight into which regions of Earth that have seen temperature increases due to global warming in a recent time span. More advanced models and tests may further increase our understanding of local temperature changes in connection to warming, but we leave this as future work.

The fBH procedure is easily applicable in practice thanks to its finite dimensional approximation. However, it should be noted that as generally is the case in functional data analysis, the performances of the procedure depends on how well the discrete grid approximates the functional data on one hand, and how well the functions are smoothed from data on the other hand.

We would like to stress the minimal assumptions required of the fBH approach: the dependence structure of functional data such as the Earth climate data can be complex and difficult to model. However, our approach does not require specific modelling of the covariance structure of the data; we merely require a certain degree of positive association among $p$-values.

Functional FDR and the fBH procedure suffer from some of the same issues as their corresponding discrete versions, see. e.g. [15]. As noted in Section 2, many authors have proposed methods or quantities to deal with multiple testing. Given the success of formulating FDR and the BH procedure in a functional framework, it is likely that some of these other methods or quantities can be expanded to the functional case as well.

Due to its simple applicability, general setting and easy understanding, we expect fFDR to have a great potential as a tool for local inference in functional data analysis. The functional BH procedure require only little computational power, and should at most be a minor issue in applications.

Appendix A Proofs

Proof of Proposition 3.5 and Proposition 3.7

A.1 Assumptions

We begin by repeating the assumptions of Proposition 3.5 and Proposition 3.7.

Definition** (PDRS).**

Let ’$\leq$’ be the partial/usual ordering on $\mathbb{R}^{l}$. An increasing set $D\subseteq\mathbb{R}^{l}$ is a set satisfying $x\in D\wedge y\geq x\Rightarrow y\in D$.

A random variable $\mathbf{X}$ on $\mathbb{R}^{l}$ is said to be PDRS on $I_{0}$, where $I_{0}$ is a subset of $\{1,\dots,l\}$, if it for any increasing set $D$ and $i\in I_{0}$ holds that

[TABLE]

Let $\mathbf{Z}$ be an infinite-dimensional random variable, where instances of $\mathbf{Z}$ are functions $T\rightarrow\mathbb{R}$. We say that $\mathbf{Z}$ is PDRS on $U\subseteq T$ if all finite-dimensional distributions of $\mathbf{Z}$ are PDRS. That is, for all finite subsets $I=\{i_{1},\dots,i_{l}\}\subseteq T$, it holds that $Z(i_{1}),\dots,Z(i_{l})$ is PDRS on $I\cap U$.

Let $\{S_{k}\}_{k=1}^{\infty},S_{1}\subset S_{2}\subset\dots$ be a dense, uniform grid in $\mathrm{\mathbf{T}}$ in sense that $S_{k}$ uniformly approximates all level sets of $p$ and $p|_{U}$ with probability one.

For Theorem 3.5 that amounts to,

[TABLE]

and

[TABLE]

whereas for Theorem 3.7, we need the density function $f$:

[TABLE]

and

[TABLE]

Furthermore, assume that $p$ is PDRS wrt. the set of true null hypotheses with probability one, and that the assumptions about $p$-value function below hold true with probability one:

1. (a1)

All level sets of $p$ have zero measure,

[TABLE] 2. (a2)

$\alpha^{*}\in(0,\alpha]\Rightarrow$: for any open neighbourhood $O$ around $\alpha^{*}$ there exists $s_{1},s_{2}\in O$ s.t. $a(s_{1})>\alpha^{-1}s_{1},a(s_{2})<\alpha^{-1}s_{2}$, where $a$ is the cumulated p-value function (Definition 3.3). 3. (a3)

$[\alpha^{*}=0]\Rightarrow\min p(t)>0$.

A.2 Proof details

For the ease of presentation, we will only consider Theorem 3.5 and furthermore assume that $\mu(\mathrm{\mathbf{T}})=1$; the latter can be done without loss of generalisation. The proof of Proposition 3.7 is analogous but notionally tedious, as the counts are replaced by sums and the measures by integrals.

Let $a_{k}$ be the cumulated p-value function for the $k$’th iteration of the BH procedure:

[TABLE]

and define the k’th step false discovery proportion $Q_{k}$ by applying the (usual) BH procedure at level $q$ to $p$ evaluated in $S_{k}$:

[TABLE]

or equivalently $b_{k}=\arg\max_{t}a_{k}(t)\geq\alpha^{-1}r$.

Lemma A.1**.**

$a_{k}$* converges to $a$ uniformly as $k\rightarrow\infty$.*

Proof: Follows from assumption (A.1) and definitions of $a_{k}$ and $a$.

Lemma A.2**.**

$b_{k}$* converges to $\alpha^{*}$ as $k\rightarrow\infty$*

Proof.

By Lemma A.1, $a_{k}$ converges uniformly to $a$. There are two cases: $\alpha^{*}=0$ and $\alpha^{*}\in(0,\alpha]$.

Case 1, $\alpha^{*}=0$

Let $O$ be any open neighbourhood around zero. $O^{C}$ (where the complement is wrt. [0,1]) is a closed set that satisfies $a(t)<\alpha^{-1}t$. By continuity of $a$, there exists an $\epsilon>0$ s.t. $a(t)<\alpha^{-1}t-\epsilon$ for all $t\in O^{C}$. As $a_{k}$ converges uniformly to $a$, eventually for large enough $k$, $a_{k}(t)<\alpha^{-1}t$ for $t\in T\backslash O$, and thus $b_{k}\in O$ eventually. This was true for any $O$, and we conclude that $b_{k}\rightarrow 0$.

Case 2, $\alpha^{*}\in(0,\alpha]$

By assumption, for any open neighbourhood $O\ni\alpha^{*}$, there exist $s_{1},s_{2}\in O$ s.t. $a(s_{1})>\alpha^{-1}s_{1}$, $a(s_{2})<\alpha^{-1}s_{2}$.

For $t>\alpha^{*},t\notin O$, we have that $\alpha^{-1}t-a(t)>\epsilon$ for some $\epsilon>0$ by continuity of $a$. Hence by uniform convergence, it must hold that for $k$ sufficiently large we have $a_{k}(t)<\alpha^{-1}t$ for $t>\alpha^{*},t\notin O$. This was true for any $O$, and we conclude $\limsup b_{k}\leq\alpha^{*}$.

Conversely, we can show that $\liminf b_{k}\geq\alpha^{*}$, and thus $\lim b_{k}=\alpha^{*}$.

∎

Define $A_{k}=\#\{t\in S_{k}:p(t)\leq b_{k}\}$, and define $Q_{k}$ as the false discovery proportion for the $k$’th iteration:

[TABLE]

Rejection areas

Now we intend to prove that $H_{t,k}$ converges eventually. Note that $p(t)$ is independent of $k$, and that $H_{t,k}=(t\in S_{k})\cap(p(t)\leq b_{k})$, i.e. the event that the BH threshold at step $k$ is larger than $p(t)$.

Proposition A.3**.**

For all $t$ that satisfies $p(t)\neq\alpha^{*}$, $H_{t,k}$ converges eventually.

Proof.

First note that if $t\notin S_{k}$ for all $k$, then $H_{t,k}$ is trivially zero for all $k$. So assume $t\in S_{k_{0}}$ for some $k_{0}$. As $k\rightarrow\infty$, $b_{k}\rightarrow\alpha^{*}$, and by assumption $p(t)\neq\alpha^{*}$. Eventually, as $k\rightarrow\infty$, $p(t)$ is either strictly larger or strictly smaller than $b_{k}$, proving the result. ∎

Convergence of $Q_{k}$

Finally we need to show that $Q_{k}\rightarrow Q$. We show this by proving convergence of the nominator and denominator, and arguing that $Q=0$ implies that $Q_{k}$ is 0 eventually.

Define $H^{0}=\{t|p(t)>\alpha^{*}\}$, ie. the acceptance region, and $H^{1}=T\backslash H^{0}$, the rejection region. Note that $\mu(H^{1})=a(\alpha^{*})=\alpha^{-1}\alpha^{*}$. Also note that $H^{1}=V\cup S=\{t:p(t)\leq\alpha^{*}\}$ and $H^{1}\cap U=V$.

Proposition A.4**.**

$N_{k}\#A_{k}\rightarrow\mu(H^{1})$* and $N_{k}\#(A_{k}\cap U)\rightarrow\mu(H^{1}\cap U)$.*

Proof.

For $k$, define $J_{k}=\{t:p(t)\leq b_{k}\}$. Note that $A_{k}=J_{k}\cap S_{k}$. Observe that by the assumption about uniform convergence on levels sets (Eq. (A.1)):

[TABLE]

Next observe that due to (1) continuity of $a$, (2) $b_{k}\rightarrow\alpha^{*}$ and (3) the fact that we are considering sets on the form $\{t:p(t)\leq x\}$, we are able to conclude that

[TABLE]

and we conclude $N_{k}\#A_{k}\rightarrow\mu(H^{1})$.

For the second part, observe that from Eq. (A.2)) it follows that

[TABLE]

We just argued that $\mu(J_{k}\triangle H^{1})\rightarrow 0$. It remains true when ”conditioning” on a measurable set, in this case $U$:

[TABLE]

and we conclude $N_{k}\#(A_{k}\cap U)\rightarrow\mu(H^{1})$. ∎

For $\alpha^{*}=0$ we have the following stronger result:

Lemma A.5**.**

If $\alpha^{*}=0$, then $\#A_{k}=0$ eventually.

From this lemma it follows that $N_{k}\#A_{k}=0$ (and thus $Q_{k}$ as well) eventually.

Proof.

Since $\alpha^{*}=0$, $a(t)<\alpha^{-1}t$ for all $t>0$. By assumption, $\min p(t)>0$, and thus $a_{k}(s)=0$ for $s<\min p(t)$ and all $k$.

By continuity of $a$, it follows that there exists $\epsilon>0$ s.t. $\alpha^{-1}t-a(t)>\epsilon$ on the interval $[\min p(t),1]$, and by uniform convergence of $a_{k}$ we get that for large enough $k$, $a_{k}(t)<\alpha^{-1}t$ for all $t\geq\min p(t)$.

Combining with $a_{k}(t)=0$ for $t<\min p(t)$, we get that eventually $a_{k}(t)<\alpha^{-1}t$ for every $t>0$ and thus $b_{k}=0$. From this (remember $\min p(t)>0$) we conclude that all hypotheses are rejected eventually, ie. $\#A_{k}=0$ for $k$ sufficiently large. ∎

Theorem A.6**.**

$Q_{k}$* converges to $Q$ almost surely, and $\limsup_{k\rightarrow\infty}\mathrm{E}[Q_{k}]\leq\alpha\mu(U)$.*

Proof.

By Lemma A.5, $Q_{k}$ converges to $Q$ when $\alpha^{*}=0$, and by Proposition A.4 $Q_{k}$ converges to $Q$ when $\alpha^{*}>0$ since $\mu(H^{1})=\alpha^{*}/\alpha>0$.

Applying Benjamini and Yekuteli’s original proposition, Theorem 2.6, (now we use the PDRS assumption), we have $\mathrm{E}[Q_{k}]\leq\alpha N_{k}\#(S_{k}\cap U)$ for all $k$. By setting $r=1$ it follows from (A.2) that $\lim_{k\rightarrow\infty}N_{k}\#(S_{k}\cap U)=\mu(U)$, and hence $\limsup_{k\rightarrow\infty}\mathrm{E}[Q_{k}]\leq\alpha\mu(U)$.

∎