{'fraction_non_alphanumeric': 0.038150330044715734, 'fraction_numerical': 0.018241181063240824, 'mean_word_length': 4.277954673159768, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
False discovery rate envelope for functional test statistics

False discovery rate (FDR) is a common way to control the number of false discoveries in multiple testing. In this paper, the focus is on functional test statistics which are discretized into m highly correlated hypotheses and thus resampling based methods are investigated. The aim is to find a graphical envelope that detects the outcomes of all individual hypotheses by a simple rule: the hypothesis is rejected if and only if the empirical test statistic is outside of the envelope. Such an envelope offers a straightforward interpretation of the test results similarly as in global envelope testing recently developed for controlling the family-wise error rate. Two different algorithms are developed to fulfill this aim. The proposed algorithms are adaptive single threshold procedures which include the estimation of the true null hypotheses. The new methods are illustrated by two real data examples.


Motivation and overview
Nowadays functional test statistics appear in many applications, for example in spatial statistics, functional regression and neuroimaging. A wealth of literature can be found for multiple comparison testing, in particular for micro-arrays applications. However, new challenges arise with functional applications: (i) the functional test statistics are highly correlated over the functional domain, (ii) the distribution of test statistic can change across the domain, and even (iii) the correlation of the test statistic can change across the functional domain. Further, the distribution of the test statistic is rarely known, therefore nonparametric methods are required. With the functional domain, also the visualization of the test statistic together with its rejection or acceptance region is of practical importance.
Recently, we developed a global envelope test (Myllymäki et al., 2017) which provides a global envelope, i.e. the acceptence region for testing with functional test statistic, under the control of family-wise error rate (FWER). The global envelope is based on the ranks of the empirical and simulated test statistics among each other, which solves the problem of the change of the distribution across the domain. The high correlation of the test statistics across the functional domain is handled by using the resampling based method to obtain the distribution under the null hypothesis. In such resampling methods based on the ranks, similar to the well known p-min procedure, there is the problem of ties due to the discreteness of ranks or p-values, which is solved in our approach by refining the extreme rank measure (p-min) by measures which uniquely order the functions, namely extreme rank length (Myllymäki et al., 2017;Narisetty and Nair, 2016), continuous rank (Hahn, 2015; or area rank . These measures also satisfy the property of intrinsic graphical interpretation (IGI) (Myllymäki and Mrkvička, 2020): the global envelope fulfills the IGI property if any function which lies at least in one point outside of such envelope is more extreme with respect to the chosen measure than all the functions forming the envelope.
In functional tests, it is often important to obtain not only if the test is significant but also estimate the whole domain which is responsible for the rejection. The first attempt in this direction was made by Xu and Reiss (2020a) by defining a step-down procedure on a functional test statistic based on the multiple testing procedure of Westfall and Young (1993). However, for such step-down procedures, the direct visualization is lost.
For the purpose of estimating all alternative hypotheses the false discovery rate (FDR) control was defined in Benjamini and Hochberg (1995). More specifically, the difference between FDR and FWER can be seen from Table 1 of outcomes when testing m hypotheses. The FWER is defined to be P(V ≥ 1), whereas FDR is defined to be E(Q), where Q = V ρ whenever ρ > 0, and Q = 0 if ρ = 0. Clearly, FWER provides more strict control than FDR. Only in the case of all hypothesis being true null (i.e. m 0 = m), the FDR control implies FWER control (Benjamini and Hochberg, 1995).
Since Benjamini and Hochberg (1995) many theoretical achievements have been made for FDR. Benjamini and Hochberg (1995) proved that their linear step-up procedure controls the FDR for independent test statistics. The subsequent article of Benjamini and Yekutieli (2001) established the control of FDR for test statistics with more general dependence structures, such as positive regression dependence. Benjamini and Hochberg (2000) and Benjamini et al. (2006) defined adaptive step-up procedures which include the estimation of  Storey (2002) proposed to control instead the positive FDR (i.e. pFDR defined as E( V ρ |ρ > 0)), which is by definition more powerful than FDR, since when positive findings occur, the interesting quantity E( V ρ |ρ > 0) is controlled at level α/P(ρ > 0) assuming the control of FDR.
Since we are dealing with highly dependent test statistics, we will stick to the resampling based methods to control FDR. The first attempt for doing so was in Yekutieli and Benjamini (1999), where the FDR conditional on S = s (see Table 1) was estimated. This was critisized by Storey and Tibshirani (2001) and they proposed a resampling method for estimating pFDR and FDR based on the estimation of the number of true null hypothesis for a single threshold rejection rule. Ge et al. (2003) then rephrased this procedure. Dudoit et al. (2008) defined a Bayes posterior resampling FDR control, under the assumption of mixture model, which assumes that the distributions of all true null statistic are equal and the distributions of all false null statistic are equal. Romano et al. (2008) constructed a bootstrap procedure to estimate critical values in a step-down FDR procedure; this method requires bootstrapping the whole data set, which makes its applicability limited. Finally, Ge et al. (2008) proposed a resampling procedure based on the step-down procedure of Westfall and Young (1993), which requires the subset pivotality and relies on the assumption that the joint distribution of statistics from the true nulls is independent of the joint distribution of statistics from the false nulls. The last assumption is hardly satisfied for functional test statistic. Also the step-down procedure does not have direct visualization. In the review paper Goeman and Solari (2014), it is claimed that none of the permutation FDR methods comes with user-friendly software.
The following software can be used for the computation of FDR: The R/Bioconductor package multtest (Gilbert et al., 2009) contains the procedures of Benjamini et al. (2006) and Dudoit et al. (2008). The package fdrtool (Strimmer, 2008) offers the unified approach for computation of FDR by a semiparametric method allowing for parametric modelling of tests interdependencies. The FDR has further been modelled by many parametric approaches, see Strimmer (2008) and the references therein. The R/Bioconductor package qvalue (Storey et al., 2020) determines the q-value of Storey (2002) for given raw p-values without resampling. The q-value gives the scientist a hypothesis testing error measure for each observed statistic with respect to pFDR. This corresponds to the multiple testing adjusted p-value in the FWER sense.

Outline of the work
The aim of this work is to introduce an envelope which instead of FWER controls the FDR. Because the data are assumed to be highly dependent, nonparametric resampling methods will be used. We also assume a general model for alternative hypotheses, since e.g. in functional permutation tests the differences between groups of functions will range from zero (null hypothesis) to small (alternative hypothesis) and large ones (alternative hypothesis). Since our aim is to develop an FDR envelope which allows for direct visualization, we will concentrate on a single threshold multiple testing procedure, such as developed by Storey (2002).
We will utilize the resampling algorithm of Storey and Tibshirani (2001). However, we propose a different algorithm for the estimation of probability that a hypothesis is a true null. Further, since we assume that the distribution of the test statistic can change across the domain, we introduce this algorithm in the space of ranks, whereas in Storey and Tibshirani (2001) the test statistic itself was considered. In order to deal with ties that occur because the resampling p-values (ranks) are discrete, we will adopt refinement measures proposed by Myllymäki et al. (2017); Narisetty and Nair (2016); Hahn (2015); . As such we will obtain an envelope which retains the IGI property, which is directly able to identify all the false hypothesis under the control of FDR, which is not sensitive to the change of the distribution of the test statistic across the functional domain, and which has almost no ties. This means that, for a given level α of FDR, we are able to find an envelope which has maximal FDR from all subsequent envelopes given by the chosen measure and its IGI. The FDR of the resulting envelope remains smaller than the given level α.
Section 2 first presents a motivating example for the FDR envelopes. The setup for multiple hypothesis testing and the proposed rejection (and acceptance) regions based on ranks are described in Section 3. Section 4 proposes algorithms for the FDR envelopes. The proposed algorithms are compared in two examples (Sections 5 and 6), and conclusions given in Section 7.
We note that the problem of ties in p-values for non resampling algorithms was also discussed recently in the literature (Chen et al., 2018;Chen, 2020).

A comparison of FWER and FDR global envelopes: a motivational example
To illustrate the merits of the FDR global envelopes, we explore linear trends in annual water temperature curves sampled at the water level of Rimov reservoir in Czech republic every day from 1979 to 2014 (Figure 1, top). Viewing the observations as daily samples of functional data, we fit the model where y i (s) is the water temperature at day s, s = 1, . . . , 365, in year i, i = 1979, . . . , 2014, and β 0 (s) and β 1 (s) are model parameters and e(s) denotes the error. Then 7000 simulations were obtained under the null hypothesis using simple permutation of raw data. The test statistic was taken to be the estimator of the regression coefficient, following the methodology proposed in . We utilized the R packages GET (Myllymäki and Mrkvička, 2020) and pppvalue (Xu and Reiss, 2020b,a). Figure 1 (middle and bottom) shows the output of the global ERL envelope test and the adjusted p-values of the step-down procedure of Xu and Reiss (2020a) on testing the null hypothesis (2). The global ERL test shows the behavior ofβ 1 (s) estimated from the data in comparison to the 95% global envelope constructed from the permutations, and detects significant increase in temperature for 24 days around 120th day of the year. The step-down p-values obtained in this example indicate no significant increases (or decreases).
Both the global ERL test as well as the step-down method control FWER. Thus, the plots show the reasons of rejection, but the tests may not find all the false null hypotheses. The FDR envelopes presented in this work control instead FDR, offering higher power for detecting the false null hypotheses. Figure 2 shows an FDR envelope proposed in this work (the envelope of Algorithm 4.3 presented in Section 4). Thus, the test controlling FDR rejected null hypotheses in spring time, around the 120th day of the year as the FWER envelope test, but additionally also in summer time around the 180th day of the year. The null hypothesis was rejected for 35 days in total.

Multiple hypothesis testing
Assume the functional hypothesis testing problem, with hypothesis H x , x ∈ R d , which is discretized into m hypotheses. Assume now that all H x are simple hypothesis. The discretization can be arbitrary, equidistances are not required. Thus we assume, that we are testing m hypotheses H k , k = 1, . . . , m, simultaneously. Denote H k = 0 if the k-th hypothesis is a true null hypothesis and H i = 1 if the k-th hypothesis is a false null hypothesis. The sets  Figure 1 by an FDR envelope. The grey area represents the 95% FDR envelope (obtained by Algorithm 4.3) and the red points denote the days with significant increase in the temperature.
W . On the other hand, the number of rejected and not rejected true null hypotheses V and U and the number of rejected and not rejected false null hypotheses S and T are not observed (see Table 1).
Every single test k is tested by the statistic T k . The vector of the statistics, T = (T 1 , . . . , T m ), can have an arbitrary distribution. The marginal distributions of T k , k = 1, . . . , m, can be interdependent and the distribution of T k can change across k. Also the interdependence can vary across k. A usual approach is to consider the parametric statistics as t, F, . . ., but since such a parametric statistic can change its marginal distribution we will switch to rank statistics of these parametric statistics. Obviously, when the statistic changes its marginal distribution across k, the simultaneous tests will be blind to alternatives that occur for k with smaller variability of test statistic T; the rank statistics overcome this problem.
When R denotes the rejection region in a single test, P(T ∈ R|H = 0) is called the type I error rate and is controlled by a predetermined level α. On the other hand, P(T / ∈ R|H = 1) is called the type II error rate and is not controlled.
In multiple testing, the problem of controlling the type I error is more complex and it is generalised in several ways. The most common is family-wise error rate (FWER), where the probability of at least one type I error is controlled, i.e. FWER= P(V ≥ 1). But also per comparison error rate can be defined as PCER= E(V )/m or per family error rate PFER= E(V ). In functional setting, Pini and Vantini (2017) recently introduced the interval wise error rate, in order to control the type I error for each subinterval of the functional domain. All the above mentioned controls are designated to control the existence of rejection of true null hypothesis, whereas the famous FDR was introduced in order to determine M 0 and M 1 . The natural way to control FDR would be by E(V /ρ), but the case of division by zero has to be treated. Therefore, the first definition of FDR (Benjamini and Hochberg, 1995) was FDR= E( V ρ |ρ > 0)P(ρ > 0). Then Storey (2002) defined positive false discovery rate by pFDR= E( V ρ |ρ > 0), which controls error rate when positive findings occurred. The rejection region for the multiple hypothesis test can be given by a step-up procedure, subsequently rejecting the hypothesis from the smallest p-values, as proposed by Benjamini and Hochberg (1995). Such rejection region is estimated to fulfill the FDR control. If an opposite approach is taken, it is possible to estimate FDR or pFDR for a general fixed rejection region Γ = (R 1 , . . . , R m ). Storey (2002) proposed to use a single value threshold γ as a natural rejection region for the continuous statistics T k . While we would like to use rank statistics, this approach brings many ties and has to be refined. For example, assume that the null hypothesis is resampled s times, thus the single test statistic can achieve only ranks from 1 to s, thus the single threshold rejection region can achieve also only these values and thus the estimated pFDR would be rather sparse and the control would be achieved only approximately. But it is possible to use information about dependence between neighbouring test statistics and use the rejection region defined by a certain measure of outlyingness. For this purpose we will use as acceptance regions the global envelopes defined in Myllymäki et al. (2017),  and Myllymäki and Mrkvička (2020). Their advantage is that they offer a useful intrinsic graphical interpretation (IGI) and rely on rank statistics.

Global envelopes
Assume now, that we have s resamples of the complete null hypothesis, i.e. m 0 = m, which gives us s vectors T i = (T i1 , . . . , T im ), i = 1, . . . , s. For such a sample, we define now a sequence of global envelopes E 1 , . . . , E s with respect to the ordering ≺ of the vectors T i , i = 1, . . . , s, which is induced by a univariate measure M i . That is, M i ≥ M j iff T i ≺ T j , which means that T i is less extreme or as extreme as T j . That is, the smaller the measure M i , the more extreme is the T i . Denote I r = {i ∈ {1, . . . , s} : M i ≥ M [r] }, where M [r] is the r-th smallest value from M 1 , . . . , M s . Then, when the alternative hypothesis is two-sided, the global envelope E r is bounded by In the one-sided case where large values of the test statistic are considered extreme, the global envelope is defined by The case where small values are significant is defined analogously.
Definition 3.1. The global envelope E r has intrinsic graphical interpretation (IGI) with respect to the ordering ≺ if low k or T ik > T The following four measures induce ordering of vectors for which E r follows the IGI.

Global rank envelope
The global rank envelope is determined by the extreme rank measure R i of the vector T i , which is defined as the minimum of the pointwise ranks of T i , namely where the pointwise rank R ik is the rank of the element T ik among the corresponding elements T 1k , T 2k , . . . , T sk of the s vectors such that the lowest ranks correspond to the most extreme values of the statistics (Myllymäki et al., 2017). How the pointwise ranks are determined, depends on whether a one-sided or a two-sided envelope is to be constructed: Let r 1k , r 2k , . . . , r sk be the raw ranks of T 1k , T 2k , . . . , T sk , such that the smallest T ik has rank 1.
In the case of ties, the raw ranks are averaged. The pointwise ranks are then calculated as for the one-sided case, where small T ik is considered extreme s + 1 − r ik , for the one-sided case, where large T ik is considered extreme min(r ik , s + 1 − r ik ), for the two-sided case.
(6) Since the extreme rank can achieve many ties, it is necessary to have large s for the global rank envelope in order to obtain precise control of FDR. The following three refinements of the extreme rank solve the ties problem.

Global extreme rank length (ERL) envelope
The extreme rank length (ERL) measure (Myllymäki et al., 2017;Narisetty and Nair, 2016) of T i is defined based on the vector of the pointwise ordered ranks where the ranks are arranged from smallest to largest, i.e., R i[k] ≤ R i[k ] whenever k ≤ k . While the extreme rank (5) corresponds to R i = R i[1] , the ERL measure takes all these ranks into account by the reverse lexical ordering. The ERL measure for T i is the normalized rank induced by the lexical ordering, where The number of functions s scales the values between 0 and 1. Even though formally the ERL measure uses all the pointwise ranks, often just the number of k where the extreme rank is reached, i.e. {#k : R ik = R i }, is enough to break the ties among those T i that obtain the same extreme rank.

Global continuous rank envelope
The ties can alternatively be broken by the continuous rank measure (Hahn, 2015; which refines the extreme rank measure R i by considering instead of the (discrete) pointwise ranks R ik continuous pointwise ranks C ik defined by the ratios to the closest other T jk . The continuous rank measure is defined analogous to the extreme rank, namely with normalization by s. The pointwise continuous rank C ik is the rank of the element T ik among the corresponding elements T 1k , T 2k , . . . , T sk of the s vectors such that the lowest ranks correspond to the most extreme values of the statistics; its formal definition is given in Appendix A following  and Myllymäki and Mrkvička (2020).

Global area rank envelope
Another refinement of rank envelope is the area rank measure , This measure breaks the ties in the extreme ranks by the sum (area) of the differences between the extreme rank R i and the pointwise continuous rank C ik from those k = 1, . . . , m where C ik is smaller than R i . These critical k correspond to those k for which R ik = R i . Thus, the larger the number of k where T i reaches the extreme rank R i and the smaller pointwise continuous ranks at these critical k are, the smaller is the area measure and the more extreme is the corresponding test statistic T i .

FDR estimation
This section describes an algorithm for FDR estimation from s resamples of complete null hypothesis. We adopt the Storey (2002) approach of estimation FDR, rather than the Benjamini and Hochberg (1995) approach of estimation of rejection region. This allows us to estimate FDR for a general global envelope E , which represents the acceptance region correspondent to the rejection region Γ.
The advantage of this algorithm is that it makes no assumptions on the correlation structure of the test statistics, and there are further no assumptions on the false null hypothesis. The only assumption here is that the true null hypotheses are same and they are simple.
Denote the observed test vector T = (T 1 , . . . , T m ) and its s resamples from complete null hypothesis by T i = (T i1 , . . . , T im ), i = 1, . . . , s. For a fixed Γ, we can observe only ρ(Γ), the number of positives, and ρ 0 i (Γ), i = 1, . . . , s, the number of positives for resampled complete null hypothesis, together with the corresponding numbers of acceptances, W (Γ) and W 0 i (Γ). The number of false positives V (Γ) is not observable, but it can be estimated. Storey and Tibshirani (2001) showed that where ρ 0 (Γ) is the number of positives under complete null hypothesis. In order to estimate the probability that a hypothesis is true null, π 0 , they proposed to choose a smaller rejection region Γ and computeπ When resampling highly correlated data, this estimator of π 0 has a significant drawback: the W (Γ ) having the correlation structure of some false null hypothesis is compared to W 0 (Γ ) having the correlation structure of all null hypothesis. This feature causes the estimator to fail for any choice of Γ . In non resampling approaches, the choice of Γ attracted many researchers in order to define so called adaptive FDR procedures based on the Storey estimator (e.g. Liang and Nettleton, 2012;Heesen and Janssen, 2016). However, since these procedures are still based on the basic Storey estimator, we have searched for another estimator of π 0 , which would be free of this difficulty. Hwang (2011) compared several estimators from which the lowest slope method developed by Benjamini and Hochberg (2000) and the mean of differences method (Hsueh et al., 2003) seem to be the best in dependent case, but Hwang (2011) did not include in the comparison the two stage approach of Benjamini et al. (2006). The two stage method was compared to the method of Benjamini and Hochberg (2000) in Benjamini et al. (2006) and it was find to be significantly better. Also this two stage approach satisfies the conservative estimation of π 0 under independence of all hypothesis. From these reasons we have chosen the two stage approach of Benjamini et al. (2006), which is defined in the following algorithm. We can now in the two stage procedure replace the linear step-up procedure by the resampling based procedure of Storey and Tibshirani (2001) and obtain the following two stage resampling based procedure.
Let r 1 be the number of rejected hypotheses for Γ * . If r 1 = 0, do not reject any hypothesis, take Γ = Γ * and stop; if r 1 = m reject all m hypotheses, take Γ = Γ * and stop; otherwise continue.

Find the largest rejection region Γ for which
The expectations in the steps 1. and 3. of the algorithm are obtained from the resamples of T 1 , . . . , T s from the complete null hypothesis, i.e.: It seems that it is not possible to prove the conservativeness of the ATSGE procedure, but it is possible to prove it for the following algorithm, where the the probability π 0 is estimated iteratively.
Algorithm 4.3. Iterative adaptive two stage global envelope (IATSGE): 1. Find the largest rejection region Γ * for which Let r 1 be the number of rejected hypotheses for Γ * . If r 1 = 0, do not reject any hypothesis, take Γ = Γ * and stop; otherwise continue.
The expectations in steps 1. and 3. of Algorithm 4.3 are obtained from the resamples of T 1 , . . . , T s from the complete null hypothesis by (10) similarly as in Algorithm 4.2. Note here that the ATSGE method controls the estimate of π 0 by assuming greater significance level in the first step, in the same way as the ATS algorithm of Benjamini et al. (2006). On the other hand, the IATSGE method controls the estimate of π 0 by adding the expected number of false null hypothesis to the estimate.

Application 1: Day-varying year effect in water temperatures
To compare Algorithms 4.1, 4.2 and 4.3 and the different types of envelopes (ERL, continuous or area rank), we applied them on the motivating example of Section 2. In order to allow the smallest p-value to be less or equal to α/m in the linear step-up procedure (Algorithm 4.1), the number of resamples must be at least m/α. Therefore, we first made 7000 resamples (365/(α/(1 − α)) = 6935), and based all the tests to this same set of permutations. Figure 2 shows the FDR envelope of the IATSGE algorithm using the ERL measure. The ATSGE algorithm with the ERL measure led to an equivalent outcome, rejecting 35 hypotheses in total. In fact, both algorithms with the continuous and area measures also led to the same outcome. On the other hand, the ATS algorithm rejected hypotheses both around the 120th and 180th day of the year, but the total number of rejected hypotheses was 33 (missing days 169 and 188 in comparison to the other).
Running the same example with 2000 permutations, all the ATSGE and IATSGE methods rejected the null for 34 days.
6. Application 2: A test of complete spatial randomness of the locations of trees Figure 3 shows the locations of Norwegian spruce trees in a plot of size a 56 × 38 metres in a natural forest stand in Saxonia, Germany. The data are available from the R library spatstat (Baddeley et al., 2015). For a spatial point pattern like this the very basic hypothesis is the complete spatial randomness (CSR) of the points. This hypothesis is commonly tested by the centred L-function, L(r) = (K(r)/π) − r, where K(r) is Ripley's K-function and r is the inter-point distance. Myllymäki et al. (2017) proposed global envelope tests for testing such hypotheses for spatial processes with the control of FWER. When the interest is to find all possible distances r at which the null hypothesis should be rejected, an alternative is to control the FDR instead of FWER.
Again all the algorithms and measures led to similar outcomes: using 2500 permutations, ATSGE and IATSGE tests rejected the null hypothesis at 243 distances out of the total 513 distances considered on the interval from 0 to 9.5 m (see Figure 4). The ATS algorithm led almost to the same result with 241 rejections.

Discussion
In this work, two different algorithms, namely ATSGE 4.2 and IATSGE 4.3, were proposed for constructing FDR envelopes for functional test statistics. In the illustrative examples, the algorithms led to rather similar outcomes as the adaptive two stage procedure of Benjamini et al. (2006) (given in Algorithm 4.1), which is the best FDR controlling method that we have found from the literature. The ATSGE and IATSGE algorithms are additionally provided by graphical interpretation showing the test statistic together with its acceptance/rejection region. Their performance is to be further considered under different scenarios.  Figure 4: Test of complete spatial randomness for the pattern of spruces shown in Figure 3 using the centred L-function. The grey area represents the 95% FDR envelope obtained by the IATSGE algorithm (4.3) and the ERL measure. The red points denote the distances with significant deviations from CSR.