Mathematical notation↩︎

Notation	Meaning
\(P_i\)	p-value corresponding to hypothesis \(\mathcal{H}_i\)
\(P_{(i)}\)	\(i\)-th sorted p-value
\(Q_i\)	p-value corresponding to the \(i\)-th true hypothesis
\(\mathcal{N}\)	the set of indices of true hypotheses
\(\mathcal{R}\)	the set of indices of rejected hypotheses
\(V\)	the number of false positives, \(|\mathcal{N}\cap \mathcal{R}|\)
\(\mathcal{C}\)	the set of all nonempty subsets of \(\{1,...,m\}\)
\(\mathcal{H}_{\mathcal{I}}\)	intersection hypothesis corresponding to the indices in \(\mathcal{I}\).
\(\mathcal{X}\)	set of all \(\mathcal{I}\in \mathcal{C}\) for which \(\mathcal{H}_{\mathcal{I}}\) is rejected by the CTP
\(\mathbb{R}\)	set of real numbers, \((-\infty,\infty)\)
\(\log(a)\)	natural logarithm of \(a\)
\(\mathbb{1}\)	indicator function
\(a\wedge b\)	minimum of numbers \(a\) and \(b\)
\(a\vee b\)	maximum of numbers \(a\) and \(b\)
\(\lceil a \rceil\)	smallest integer that is at least \(a\)
\(\lfloor a \rfloor\)	largest integer that is at most \(a\)
\(|a|\)	absolute value of \(a\) (if \(a\) is a number)
\(|\mathcal{A}|\)	cardinality of \(\mathcal{A}\) (if \(\mathcal{A}\) is a set)
\(\mathcal{A}\cup \mathcal{B}\)	union of sets \(\mathcal{A}\) and \(\mathcal{B}\)
\(\mathcal{A}\cap \mathcal{B}\)	intersection of sets \(\mathcal{A}\) and \(\mathcal{B}\)
\(\mathcal{A}\setminus \mathcal{B}\)	set \(\mathcal{A}\) minus set \(\mathcal{B}\)
\(\mathcal{A}^c\)	\(\{1,....,m\}\setminus\mathcal{A}\) (if \(\mathcal{A}\subseteq \{1,...,m\}\))
\(\mathcal{A}\subseteq \mathcal{B}\)	\(\mathcal{A}\) is a subset of \(\mathcal{B}\)
\(\emptyset\)	empty set
\(\bigcup_{...}...\)	union over ... of ...
\(\bigcap_{...}...\)	intersection over ... of ...
\(\exists\)	there exists
\(\forall\)	for all

Abbreviations↩︎

Notation	Meaning
i.i.d.	independent and identically distributed
CDF	cumulative distribution function
FWER	familywise error rate
k-FWER	k-familywise error rate
FDX	false discovery exceedance
FDP	false discovery proportion
FDR	false discovery rate
MTP	multiple testing procedure
CTP	closed testing procedure

1 Introduction↩︎

When we test a single hypothesis, we are interested in the probability of committing a type I error, as well as the power. When we test multiple hypotheses, it likewise makes sense to try to know the risks of incurring false positives. At the same time, we want to have a good chance of rejecting hypotheses that are false. These are the aims of multiple hypothesis testing procedures (MTPs).

There are many MTPs. One reason is that there a various kinds of demands that we could make regarding type I error control. For example, some methods ensure that the probability of incurring any false positives is small. This is called familywise error rate (FWER) control. Other methods are more lenient. Another reason is that there are different ways to measure power, and given some definition of power, which method is most powerful may depend on the situation. Further, different methods require different assumptions. Finally, knowledge on multiple testing has improved throughout the years, and sometimes an older method is “uniformly improved” by a new method.

This text covers some of the most important and widely used MTPs. We start off with a section (§2) about testing a single hypothesis, focusing on ideas that will return in later sections. §3 covers some basic building blocks of multiple testing theory, which will return in later sections as well. For example, in §3 we cover various global tests, which will be important in e.g. §4. In that section we discuss closed testing for FWER control. This is a fundamental theory that can be used for designing and understanding FWER controlling methods. §5 covers powerful resampling-based methods for FWER control. This section extends resampling-based tests from §2 to multiple testing settings. §6 discusses false discovery exceedance control, which is less conservative than FWER control. The Romano-Wolf method discussed in §6 relies on so-called k-FWER methods, which are covered in earlier sections.

§7 introduces an approach to multiple testing that is quite different from those in other sections. Traditional MTPs tell the user which hypotheses to reject, and the user has limited freedom to investigate other hypotheses afterwards using the same data. However, the methods in §7, which connect to §4, provide information about all possible subsets of hypotheses simultaneously: for each subset of hypotheses, we obtain a lower bound for the number of false hypotheses in the set. With large probability, all these bounds are simultaneously correct, so that the user can freely explore various subsets of hypotheses in a post hoc manner.

§8 and §9 are about false discovery rate (FDR) control, which means ensuring that the expected proportion of false positives is small. The theory in these sections is quite different from that of §§4-7, so that §8 and §9 can mostly be read independently from those sections. Two proofs in §8 and §9 rely on some limited martingale theory, which is covered in the appendices.

The end of each section offers various exercises that help improve and expand knowledge of the topics. Very nonstandard exercises are indicated with a star. There are many relevant methods which this text does not cover. For example, there is a huge literature on FDR controlling procedures. This text only covers a few such methods, which are relatively prominent. Other relevant topics that are not covered are e.g. graph-based procedures [1], Bayesian approaches [2], e-value-based approaches [3] and selective inference. Regarding the topics that are covered, only a selection of the available theory is provided, and some theory is only developed in exercises.

2 Testing a single hypothesis↩︎

This section recalls core ideas on hypothesis testing. Nearly all ideas introduced here will return in later sections. Moreover, this section defines the meanings of various terms, whose interpretations may differ slightly between textbooks.

2.1 Hypotheses and p-values↩︎

A statistical model is a set of distributions in which we assume that the true distribution of our data lies. More formally, a model is a set of probability measures (see §10.1). For example, the model may be that the data \(X\) is an \(n\)-dimensional vector with i.i.d. normally distributed variables, with unknown mean \(\mu\in \mathbb{R}\) and variance \(\sigma^2>0\). Often we are interested in testing whether the true distribution of our data lies in some particular set of distributions, which is a subset of our model. More precisely, we want to see whether we can statistically prove that the true distribution does not lie in this subset. This subset of our model is then called the null hypothesis \(\mathcal{H}_0\). For example, we may be interested in testing whether \(\mu=0\), i.e., we want to see whether we can prove that \(\mu\neq0\). The null hypothesis is then the set of all distributions in our model for which \(\mu=0\). The null hypothesis could also be \(\mathcal{H}_0: \mu\in [-0.1,0.1]\). This is an example of a composite null hypothesis, because \(\mathcal{H}_0\) corresponds to multiple parameter values. Even if the null hypothesis is \(\mathcal{H}_0: \mu=0\), we will call it composite, since \(\mathcal{H}_0\) contains multiple distributions, namely all \(n\)-vectors of i.i.d. normal variables with mean 0 and positive variance. A hypothesis that contains a single distribution is called simple.

It can be useful to formulate an alternative hypothesis \(\mathcal{H}_{\text{a}}\) as well, which is a subset of our model that does not intersect with the null hypothesis. \(\mathcal{H}_{\text{a}}\) is the set of distributions under which we wish to have a high probability of rejecting \(\mathcal{H}_0\). \(\mathcal{H}_{\text{a}}\) may influence which test statistic we use (see below). Specifying an alternative hypothesis can be relevant when considering power properties of a test, i.e., when studying the probability of rejecting \(\mathcal{H}_0\) given that the data follow \(\mathcal{H}_{\text{a}}\).

Hypothesis testing is related to the concept of a p-value. We will call a random variable \(P\) a valid p-value if it satisfies the following property under every distribution in \(\mathcal{H}_0\): \[\label{eqlunif} \mathbb{P}(P\leq c)\leq c \qquad \text{for all }c\in[0,1].\tag{1}\] This property means that \(P\) is stochastically larger than standard uniform, meaning its CDF is smaller than (or equal to) the CDF of the \(U[0,1]\) distribution.

In order to compute a p-value, we typically compute a test statistic \(T=T(X)\), which is a function of the data \(X\), say. We try to define \(T\) in such a way that it is likely to reject (see below) \(\mathcal{H}_0\) when \(\mathcal{H}_0\) is false, or similarly, when a particular alternative hypothesis is true. A test statistic is a random variable. When we consider a specific realized value of the test statistic, we will write \(T_{\text{obs}}\) here. In case \(\mathcal{H}_0\) is a simple hypothesis and large values of \(T_{\text{obs}}\) are evidence against \(\mathcal{H}_0\), we often compute a p-value as follows: \[p = \mathbb{P}_{\mathcal{H}_0}(T\geq T_{\text{obs}}).\] Here \(\mathbb{P}_{\mathcal{H}_0}(...)\) means the probability of “...” under \(\mathcal{H}_0\), i.e., if \(\mathcal{H}_0\) is true.

If we are testing a hypothesis \(\mathcal{H}_0:\theta=\theta_0\) about a one-dimensional parameter \(\theta\in \mathbb{R}\) and we want to have power against left-sided alternatives, then usually, very low values of \(T_{\text{obs}}\) are evidence against \(\mathcal{H}_0\), so that we compute the p-value as \[p = \mathbb{P}_{\mathcal{H}_0}(T\leq T_{\text{obs}}).\] If \(\theta\in \mathbb{R}\) and we want power in both directions, we can compute a left-sided and a right-sided p-value and take \(p\) to be 2 times their minimum. We wrote “\(p\)” as a lowercase letter, because we considered its realized value. Of course a p-value is a random variable and we will often write “\(P\)” to emphasize that.

When \(\mathcal{H}_0\) is a composite hypothesis (see above), we often compute the p-value as follows: \[p = \sup_{\mathbb{P}\in\mathcal{H}_0}\mathbb{P}(T\geq T_{\text{obs}}).\] To understand this notation, realize that \(\mathcal{H}_0\) is a set of distributions, i.e., it is a set of probability measures \(\mathbb{P}\) on the sample space (see §10.1). By taking the supremum, we ensure that \(p\) is larger than \(\mathbb{P}(T\geq T_{\text{obs}})\) for every \(\mathbb{P}\in\mathcal{H}_0\).

If \(\mathcal{H}_0\) is true and \(P\) is computed in one of the ways indicated, then 1 is satisfied. When \(\mathcal{H}_0\) is simple and the null distribution of \(T\) is continuous and \(P\) is computed in one of the above-mentioned ways, then \(P\) is exactly standard uniform on \([0,1]\) under \(\mathcal{H}_0\). It is equally valid to say that \(P\) is exactly standard uniform on \((0,1]\) or \((0,1)\), since a standard uniform variable is strictly larger than 0 and smaller than 1 almost surely (see Appendix 10.1). Note that if \(P\) is valid, then \(\mathbb{P}_{\mathcal{H}_0}(P=0)=0\). This means that under \(\mathcal{H}_0\), a valid p-value cannot be 0 in practice.

2.2 Hypothesis testing↩︎

One typically rejects \(\mathcal{H}_0\) when \(P\leq \alpha\), where \(\alpha\in(0,1)\) is some prespecified significance threshold. This ensures that the rejection probability is at most \(\alpha\) if \(\mathcal{H}_0\) is true. Note that if \(p\) is computed as \(\mathbb{P}_{\mathcal{H}_0}(T\geq T_{\text{obs}})\), then it is equivalent to rejecting when \(T\) exceeds its \((1-\alpha)\)-quantile. When \(\mathcal{H}_0\) is rejected, it is common to say that the test or p-value is “significant”. Otherwise we do not reject (so we “retain” or “accept” \(\mathcal{H}_0\)). That does not mean that we conclude that \(\mathcal{H}_0\) is true; when we do not reject \(\mathcal{H}_0\), we simply do not make a claim at all. Since we can reject or not, and \(\mathcal{H}_0\) can be true or false, there are four possible outcomes of a hypothesis test, see Table 1.

A true positive is often called a “true discovery” and a false positive is sometimes called a “false discovery”. Of course, this is because scientific discoveries are often based on rejecting a null hypothesis. For example, if the parameter \(\theta\in \mathbb{R}\) indicates the effect of some genetic mutation on a disease or the effect of some governmental policy, then rejecting \(\mathcal{H}_0:\theta=0\) means “discovering” that there is some effect.

If we consider some distribution in \(\mathcal{H}_0\), then we may consider the type I error rate, which is also called the level and means the probability of rejecting \(\mathcal{H}_0\), when the data follow that null distribution. If we consider some distribution in the alternative hypothesis \(\mathcal{H}_{\text{a}}\), then the power is the probability of rejecting \(\mathcal{H}_0\) if the data follow that distribution. The type II error rate is then \(1\) minus the power.

If we are testing a hypothesis \(\mathcal{H}_0: \theta=\theta_0\) about a one-dimensional parameter \(\theta\in \mathbb{R}\), then often we desire to have power in both directions, as mentioned above. Different approaches are then possible, after we have computed a statistic \(T\) that is 0-centered under \(\mathcal{H}_0\):

1. compute a valid left-sided and right-sided p-value and take \(p\) to be 2 times their minimum. This gives a valid p-value for \(\mathcal{H}_0\) and we can reject if \(p\leq \alpha\);

2. reject if \(T\) lies below its \(\alpha/2\) or above its \((1-\alpha/2)\)-quantile;

3. compute the absolute value \(|T|\) and compare it with its \((1-\alpha)\)-quantile.

If the test statistic is symmetric about 0 under \(\mathcal{H}_0\), then these approaches are equivalent. If the test statistic is not symmetric, then the third approach is not equivalent to the first two.

The value \(\alpha\) usually indicates the desired or maximally allowed type I error rate, specified by the user. It is sometimes called the nominal type I error rate. When the actual type I error rate is (potentially) larger than \(\alpha\), then we call the test liberal or anti-conservative. We say that a test controls the type I error rate (or is “valid”) if \(\mathbb{P}(\text{reject})\leq \alpha\) (for every distribution \(\mathbb{P}\) in the null hypothesis, or for some specified set of null distributions). When the test has a type I error rate strictly below \(\alpha\) (for every \(\mathbb{P}\in \mathcal{H}_0\) or for some specified set of distributions in \(\mathcal{H}_0\)), we call the test conservative. The terms size and level are related and sometimes used interchangeably. We define the size of a test to be \(\sup_{\mathbb{P}\in\mathcal{H}_0}\mathbb{P}(\text{reject})\). A test is called exact when \(\mathbb{P}(\text{reject})=\alpha\) (under some specified set of null distributions). The definitions of these terms can differ slightly between textbooks.

A concept closely related to hypothesis tests, are confidence intervals. If we consider a one-dimensional parameter \(\theta\) and a test for hypotheses of the form \(\mathcal{H}_0:\theta=\theta_0\), then a confidence interval for \(\theta\) can be constructed as the set of all \(\theta_0\) for which the corresponding hypothesis test does not reject. This is called inverting a test.

Note that when a null hypothesis \(\mathcal{H}_0\) is rejected (a “discovery”), this often leads researchers to publish the result (if the estimated effect is considered substantial). On the other hand, if \(\mathcal{H}_0\) is not rejected, then this usually not a reason to publish the result. The result may still be published, but even then it will usually receive less attention than rejected hypotheses. This is one of the reasons why type I errors are often considered to be more problematic than type II errors. Also note that a type II error is in a sense not an error at all, because we do not make any claim when we fail to reject \(\mathcal{H}_0\).

In recent years, many articles have been published that warn against an overly dichotomous approach to hypothesis testing, where a finding is treated as important when \(p\leq \alpha\) and unimportant otherwise. In line with this, some researchers argue that “significance testing” should be abandoned in their field. They often propose to instead look at how small the p-value is and to look at the size of the estimated effect — or at the confidence interval. For example, when the sample size is large, then the p-value may be very small, but this does not mean that the estimated effect is large; the estimated effect may be very small, so that it is not very relevant. Note that some people have not only suggested to ban significance testing, but to even stop reporting p-values. However, to many statisticians that goes too far, since a p-value provides valuable information in addition to an effect estimate. Indeed, a large estimated effect is highly dubious when the sample size is too small. Note that confidence intervals are also highly useful in this context, although they can be less practical when the parameter of interest is multi-dimensional.

2.3 Examples of parametric tests↩︎

A parametric statistical model is a model in which finitely many parameters fully describe the possible distributions. An example is the model mentioned earlier, where we assume that the data \((X_1,...,X_n)\) consists of i.i.d. \(N(\mu,\sigma^2)\) variables. That model consists of infinitely many distributions, but each distribution is fully specified by the two parameters \(\mu\in \mathbb{R}\) and \(\sigma>0\).

If we only assume that \((X_1,...,X_n)\) consists of i.i.d. variables, and make no assumption about the distributional shape of the variables, then our model is nonparametric. Such a model is useful if we are not comfortable with assuming that the observations are e.g. normally or exponentially distributed. Another advantage of nonparametric approaches is that they can be very powerful in a multiple testing context. Nonparametric models will be considered later. We first discuss some simple examples of parametric models.

2.3.1 Example 1: Testing for Correlation↩︎

Suppose we observe i.i.d. pairs \((X_1,Y_1),\ldots,(X_n,Y_n)\), where the variables \((X_i,Y_i)\) are assumed to be jointly normally distributed. Note that this model is fully parametric, because if we know the mean vector and covariance matrix of \((X_i,Y_i)\), we know the distribution of the data.

A common question is whether the two variables are correlated. Let \(\rho\) denote the true Pearson correlation between \(X\) and \(Y\). We wish to test \[H_0: \rho = 0 \qquad \text{versus} \qquad H_{\text{a}}: \rho \neq 0.\]

The test statistic is based on the sample Pearson correlation coefficient \[r = \frac{\sum_{i=1}^n (X_i - \bar X)(Y_i - \bar Y)}{\sqrt{\sum_{i=1}^n (X_i - \bar X)^2}\sqrt{\sum_{i=1}^n (Y_i - \bar Y)^2}}.\] If \(\mathcal{H}_0\) is true, the statistic \[T = r \sqrt{\frac{n-2}{1-r^2}}\] follows a \(t\)-distribution with \(n-2\) degrees of freedom [4, p. 218].

The p-value is obtained by comparing the observed value of \(T\) to this reference distribution. We can for example obtain the two-sided p-value by computing 2 times the minimum of the left- and the right-sided p-value. If the model is correct, then \(P\) is exactly standard uniformly distributed under \(\mathcal{H}_0\). This will usually not be the case if the multivariate normality assumption fails. However, when \(n\) is large, the test is still roughly exact, even when that assumption fails.

In Table 2 two columns of some rows of the dataset mtcars from the dataset package in R are shown. This dataset contains characteristics of various car models; here we look at the continuous variables mpg (miles per US gallon) and weight (in 1000 pounds). Applying the correlation test described above (using the R function cor.test()) leads to the p-value \(p= 0.2412\).¹ The empirical Pearson correlation is \(-0.6437\).

2.3.2 Example 2: Non-inferiority Testing↩︎

Suppose we compare a new treatment to a standard treatment. For example, we may want to compare the effect of a new policy with the effect of the standard policy, or the effect of a new drug with the effect of the standard drug. Let \(\mu_N\) and \(\mu_S\) denote the mean outcomes under the new and standard treatments, respectively. Suppose higher values are better.

We are willing to accept the new treatment as adequate if it is not worse than the standard treatment by more than a pre-specified margin \(\Delta > 0\). This leads to the hypotheses \[\mathcal{H}_0: \mu_N - \mu_S \le -\Delta \qquad \text{versus} \qquad \mathcal{H}_1: \mu_N - \mu_S > -\Delta.\]

Here, the null hypothesis represents the case where the new treatment is inferior by more than the acceptable margin. Rejecting \(\mathcal{H}_0\) provides evidence that the new treatment is non-inferior to the standard treatment.

Note that \(\mathcal{H}_0\) is a composite null hypothesis, since it corresponds to all values of \(\mu_N-\mu_S\) below \(-\Delta\). Assume that the differences \[D_1=X^N_1-X^S_1,\ldots,D_n=X^N_n-X^S_n\] are i.i.d. and normally distributed with mean \(\mu_N-\mu_S\).

To test \(\mathcal{H}_0\), we compute the t-statistic

\[T=\frac{\bar{D}-(-\Delta)}{\hat{\sigma}/\sqrt{n}},\] where \[\hat{\sigma} = \sqrt{\frac{1}{n-1} \sum_{i=1}^n( D_i - \bar{D})^2}.\]

The p-value is \[\sup_{\mathbb{P}\in\mathcal{H}_0}\{\mathbb{P}(T\geq T_{\text{obs}})\}=\] \[\mathbb{P}_{\mu_N-\mu_S=-\Delta}(T\geq T_{\text{obs}}).\] This probability can be computed by using that if \(\mu_N-\mu_S=-\Delta\), then \(T\) has a t-distribution with \(n-1\) degrees of freedom.

2.3.3 Example 3: Equivalence Testing↩︎

In equivalence testing, the goal is to show that two variables have similar means, up to a small and practically irrelevant difference. For example, in the field of food safety testing, the aim may be to show that the concentration of a certain molecule in a new variety of potato is similar to the concentration in an existing variety.

Let \(\mu_A\) and \(\mu_B\) denote the means of the two variables, and let \(\Delta > 0\) be a pre-specified equivalence margin. The hypotheses are \[\mathcal{H}_0: |\mu_A - \mu_B| \ge \Delta \qquad \text{versus} \qquad \mathcal{H}_1: |\mu_A - \mu_B| < \Delta.\]

Equivalently, the null hypothesis can be written as \[\mathcal{H}_0 = \mathcal{H}_0^-\cup \mathcal{H}_0^+,\] where \[\mathcal{H}_0^-: \mu_A - \mu_B \le -\Delta \;\;, \qquad \;\; \mathcal{H}_0^+: \mu_A - \mu_B \ge \Delta.\]

Rejecting \(\mathcal{H}_0\) implies that the difference between the means is small enough to be considered practically negligible, and we conclude that the treatments are equivalent. In practice, equivalence testing is often implemented using Two One-Sided Tests (TOST). Thus, we test \(\mathcal{H}_0^-\) and \(\mathcal{H}_0^+\) separately, e.g. using the test from§2.3.2. We reject \(\mathcal{H}_0\) if and only if both \(\mathcal{H}_0^-\) and \(\mathcal{H}_0^+\) are rejected. Note that this is a simple example of multiple hypothesis testing, since we test two hypotheses, \(\mathcal{H}_0^-\) and \(\mathcal{H}_0^+\). However, we do not need a multiple testing correction here. The reason is that at most one of the hypotheses \(\mathcal{H}_0^-\) and \(\mathcal{H}_0^+\) can be true, so that the probability of rejecting a true hypothesis stays below \(\alpha\) with the TOST procedure (see also Example 4).

2.4 Permutation tests and related nonparametric tests↩︎

In §2.4.1 we indicate what nonparametric tests are and why they are useful. In the subsequent sections we give some examples of permutation tests and related tests, to give an impression of how these tests work. Then we provide somewhat general theory on permutation tests. Finally, we discuss bootstrap tests.

2.4.1 Nonparametric models and tests↩︎

Earlier we mentioned the distinction between parametric and nonparametric models: parametric models involve a finite number of parameters, that completely define the possible distributions under this model. An example is a negative binomial model, which is completely define by the coefficients and the overdispersion parameter.

Nonparametric models are models that are not defined in terms of any parameters. For example, the model may be that we have i.i.d. observations \(X_1,...,X_n\), with no known distributional shape.

Semiparametric models are partly defined by parameters, but not completely. An example is a model that is the same as a classical Gaussian linear model, except that we do not assume that the residuals are normally distributed or have any other known shape. That model is sometimes called the general linear model — not be confused with a generalized linear model.

A nonparametric test is a test that can be applied when consider a nonparametric model. An example of a nonparametric test is a permutation test. A permutation test can also be semiparametric, in the sense that a semiparametric model is assumed. For example, there are permutation tests for testing hypotheses about coefficients in the non-Gaussian, general linear model mentioned above.

One advantage of nonparametric tests such as permutation tests is well-known: they require fewer distributional assumptions than parametric tests in order to provide reliable inference. In the context of multiple testing however, there is another major potential advantage: permutation-based (and related) multiple testing methods can be much more powerful than other multiple testing methods. The reason is that they can adapt to the dependence structure in the data, as will be discussed later.

2.4.2 Example 1: a permutation test↩︎

We again consider the dataset from Table 2, which is a subset of the dataset mtcars. Suppose we wish to know whether the variables mpg (\(X\), say) and weight (\(Y\), say) are independent of each other, so \(\mathcal{H}_{\text{indep}}: X \mathrel{\perp\!\!\!\perp}Y\). If we proceed as in §2.3.1 and assume that the pairs \((X_1,Y_1),\ldots(X_5,Y_5)\) are i.i.d. and multivariate Gaussian, then we can obtain a valid p-value for \(\mathcal{H}_0\) by using the parametric correlation test. Indeed, if those assumptions are valid and \(\mathcal{H}_{\text{indep}}\) is true, then the zero-correlation hypothesis \(\mathcal{H}_0\) from §2.3.1 is also true. Thus the probability of rejecting \(\mathcal{H}_0\) will be at most \(\alpha\) and likewise \(\mathbb{P}(P\leq c)\leq c\) for every \(c\in[0,1]\).

However, suppose that the model assumed in §2.3.1 is not valid. More precisely, suppose that \((X_1,Y_1),\ldots(X_5,Y_5)\) are i.i.d. but not multivariate Gaussian. Then the properties of the parametric correlation test break down. We can then instead perform a permutation test, which does not require parametric assumptions such as normality. Besides requiring fewer or no paramatric assumptions, permutation or bootstrap methods also have an additional advantage when we are testing multiple hypotheses; namely, permutation or bootstrap methods then often have superior power (see e.g. §5).

All we need for the permutation test to control the type I error rate, is that under \(\mathcal{H}_{\text{indep}}\), the observations \(Y_1,...,Y_n\) are i.i.d. conditional on \(X=(X_1,...,X_n)\) (or the other way around). The permutation test proceeds as follows. First it computes a test statistic based on the original data. A typical choice is to use the empirical correlation \(r=r(X,Y)\), which is \(-0.6437\) for these data. Note that there are \(n!=5!\) ways to shuffle (permute) the values \(Y=(Y_1,...,Y_n)\), if we include the original ordering. The test permutes \(Y\) in all possible ways, thus creating \(5!-1\) “fake” datasets. For each of these permuted datasets, the test recomputes the test statistic \(r(X,\pi(Y))\), where \(\pi(Y)\) is a permuted version of \(Y\) and \(r(X,\pi(Y))\) is the corresponding empirical Pearson correlation. A histogram of these test statistics is shown in Figure 4. This is a reference distribution — often termed a permutation distribution — with which we compare our original test statistic \(-0.6437\). This original test statistic is somewhat in the left tail of the reference distribution, but not very extreme.

Note that each permutation map \(\pi()\) is a function from \(\mathbb{R}^5\rightarrow \mathbb{R}^5\). Let \(\mathcal{G}\) be the set containing all these 120 permutation maps. It turns out that we can compute a valid right-sided permutation p-value as \[P_r := \frac{|\{\pi\in \mathcal{G}: r(X,\pi(Y))\geq r(X,Y)\}|}{|\mathcal{G}|} =\frac{107}{120}.\] and valid left-sided permutation p-value as \[P_l := \frac{|\{\pi\in \mathcal{G}: r(X,\pi(Y))\leq r(X,Y)\}|}{|\mathcal{G}|} =\frac{14}{120}.\] A valid two-sided p-value is then \[P_{lr}:=2\cdot\min\{P_r,P_l\}=28/120\approx 0.2333.\]

An alternative two-sided p-value — which is also valid but may be somewhat different from the two-sided p-value above — is \[P_{lr}':= \frac{|\{\pi\in \mathcal{G}: |r(X,\pi(Y))|\geq |r(X,Y)|\}|}{|\mathcal{G}|} =\frac{29}{120}.\] In this example, \(P_{lr}'=29/120\approx 0.2417,\) which is slightly different from \(P_{lr}\). In general, these p-values are close when the permutation distribution (shown in Figure 4) is approximately symmetric. As usual, we reject \(\mathcal{H}_{\text{indep}}\) when the p-value is at most \(\alpha\). This ensures the type I error rate is at most \(\alpha\).

It turns out that if the data are continuous (so that with probability 1 there are no ties among the test statistics), then both these p-values are exact. More precisely, if \(\mathcal{H}_{\text{indep}}\) is true, then \(P_{lr}\) is uniformly distributed on \[\{2/120, 4/120, 6/120,...,120/120\}\] and \(P_{lr}'\) is uniform on \[\{1/120, 2/120, 3/120,...,120/120\}.\] Note that an advantage of using \(P_{lr}'\) is that its smallest possible value is \(1/120\), whereas the smallest possible value of \(P_{lr}\) is 1/60.

Recall that the p-value computed with the parametric correlation test (§2.3.1) was \(p= 0.2412\). This is close to the permutation-based p-values. This does not have to be case in general, in particular when \(n\) is small and the data are not Gaussian.

2.4.3 Example 2: a nonparametric test based on sign-flipping↩︎

In [5, Ch. IV] Fisher first proposed a certain nonparametric (or perhaps semiparametric) test based on sign-flipping. This test is somewhat related to the permutation test above. Fisher applied the test to a dataset collected by Charles Darwin. The dataset contains paired observations of heights of maize plants (crossed vs. self-fertilized plants). The first few columns of the dataset are shown in Table 3. For simplicity we will only use these observations as our dataset.

The observations are paired, which means that the rows are independent of each other, but the two observations in each row may be dependent of each other. Denote the differences between the first and second observation in each row by \(d=(d_1,...,d_6)\). A classical parametric way to test the null hypothesis \(\mathcal{H}_0\) that the \(d_i\) have mean 0, is to assume the \(d_i\) are i.i.d. and normally distributed and to perform a t-test. However, we then require the following (if \(\mathcal{H}_0\) is true):

• The \(d_i\) are independent of each other;

• The \(d_i\) are normally distributed;

• The \(d_i\) all have the same variance.

Fisher’s sign-flipping test requires fewer assumptions. It only requires the following:

• The \(d_i\) are independent of each other;

• The \(d_i\) are symmetric about their means.

The test works as follows: we consider all \(2^n=2^6\) sign-flipped versions of the data, i.e., all vectors of the form \((s_1 d_1,...,s_6 d_6)\), where \((s_1,...,s_6)\in\{-1,1\}^6\). Let \(\mathcal{G}_{\pm}\) be the set of all \(2^6\) transformations of the form \((x_1,...,x_6)\mapsto (s_1x_1,...,s_6x_6)\), where \((s_1,...,s_6)\in\{-1,1\}^6\). Thus, the sign-flipped versions of the data are all vectors \(g(d)\), \(g\in \mathcal{G}_{\pm}\).

For each of these \(2^n\) transformed vectors \(g(d)\), we compute the sum of the entries: \[T(g(d))=\sum_{i=1}^n (g(d))_i.\] A histogram of the resulting sums (test statistics) is shown in Figure 5. Note that it is symmetric (why?). To compute a p-value, we can proceed as in §2.4.2, with the difference that we use different transformations (sign-flips instead of permutations) and different test statistics (sum instead of empirical correlation). For example, a valid two-sided p-value is \[P_{lr}:= \frac{|\{g\in \mathcal{G}_{\pm}: |T(g(d))|\geq |T(d)|\}|}{|\mathcal{G}_{\pm}|} =\frac{44}{64}.\] If we would compute a p-value using a t-test, we would likely get a similar result. Indeed, there is no sign that these particular data are very non-Gaussian or heteroscedastic.

Darwin’s full dataset has 15 rows. If we would apply the test to the full dataset, this would mean applying \(2^{15}=32768\) different sign-flipping transformations. In Fisher’s time, that was was very time-consuming, although he managed to perform the test for this particular dataset, see [5, Sec. 21]. Nowadays it is easily doable with a computer.

2.4.4 General theory on permutation and group invariance tests↩︎

In the examples above, we considered permutation maps and sign-flipping maps. It turns out that these sets of transformations are groups, in the algebraic sense. A set \(\mathcal{G}\) of transformations \(g:\mathcal{S}\rightarrow\mathcal{S}\) from a sample space \(\mathcal{S}\) to itself is called a group when it has the following properties:

1. every element \(g\in \mathcal{G}\) has an inverse \(g^{-1}\) in \(\mathcal{G}\);

2. if \(g\), \(h\in \mathcal{G}\), then \(g\circ h\in \mathcal{G}\);

3. the identity map \(id\) is contained in \(\mathcal{S}\);

Here \(g\circ h\), or \(gh\) for short, indicates the composition of the maps \(g\) and \(h\), so that \(g\circ h(x)=g(h(x))\). Note that the third property is in fact implied by the first two. An important property of groups that we will use, is that \[\label{GgisG} \forall g\in \mathcal{G}: \quad \mathcal{G}g=\mathcal{G},\tag{2}\] where \(\mathcal{G}g\) means \(\{h\circ g: h\in \mathcal{G}\}\). Because tests such as those in §2.4.2 and §2.4.3 involve groups of transformations, they are sometimes called group invariance tests. Thus, a permutation test for example, is a special case of a group invariance test.

We will call the data \(X\). This could e.g. be a random vector or matrix. A sufficient property that makes the tests valid (and often exact) is that under \(\mathcal{H}_0\), the distribution of the data is invariant under the transformations, i.e., \[\label{eqperminv} \forall g\in \mathcal{G}: \quad X \,{\buildrel d \over =}\, gX,\tag{3}\] where \(gX\) is short for \(g(X)\). In fact, it is sufficient that the following assumption (which is implied by property 3 ) is satisfied. We will assume that \(\mathcal{G}\) is finite, although this can be generalized.

The test statistic \(T(\cdot)\) can be anything, for example a left-sided test statistic or a two-sided test statistic (e.g., the absolute value of a t-statistic). In order to have an asymptotically exact permutation test, it is not necessary for Assumption 1 to be exactly satisfied [6]–[8].

In the examples in §2.4.2 and §2.4.3, we computed a test statistic for each transformed version \(gX\) of the data. In general a group invariance test considers \(|\mathcal{G}|\) test statistics \(T(gX)\), \(g\in \mathcal{G}\) and bases the test on these test statistics. The null hypothesis \(\mathcal{H}_0\) is rejected if the original statistic \(T(X)\) is extreme compared to the reference statistics \(T(gX)\), \(g\in \mathcal{G}\). The p-value can be computed as \[\label{eq:groupinvpv} P= \frac{|\{g\in \mathcal{G}: T(gX)\geq T(X)\}|}{|\mathcal{G}|}.\tag{5}\]

A basic group invariance test test rejects \(\mathcal{H}_0\) when \(T(X)>T^{(k)}(X)\), where \[T^{(1)}(X)\leq...\leq T^{(|\mathcal{G}|)}(X)\] are the sorted test statistics \(T(gX)\), \(g\in \mathcal{G}\), and \(k=\lceil (1-\alpha)|\mathcal{G}|\rceil\), where \(\lceil\cdot\rceil\) means rounding upwards to an integer. Note that the p-value 5 is at most \(\alpha\) if and only if the group invariance test rejects. As is stated in the following theorem, the test has level at most \(\alpha\).

The above proof and a different proof can be found in e.g. [9]. The different proof considers the orbit of the data \(X\), which is defined as \(O_X=\{gX:g\in G\}\) and shows that under \(\mathcal{H}_0\), conditional on the orbit \(O_X\), the probability of rejecting \(\mathcal{H}_0\) is at most \(\alpha\).

Any proof of the validity of the group invariance test, strongly relies on the group structure of the set \(\mathcal{G}\) of transformations. Indeed, the group invariance test is based on the property that under \(\mathcal{H}_0\), for each \(g\in \mathcal{G}\), the probability \(\mathbb{P}\big\{T(gX)>T^{(k)}(X)\big\}\) is the same. That follows from the fact that for every \(g\in \mathcal{G}\), the joint distribution of \((gX,\mathcal{G}X)\) is the same. That holds because if \(g, g'\in \mathcal{G}\), then under \(\mathcal{H}_0\) we have \[(gX,\mathcal{G}X) = (gX,\mathcal{G}g X) \,{\buildrel d \over =}\, (X,\mathcal{G}X) \,{\buildrel d \over =}\, (g'X,\mathcal{G}g' X) =(g'X,\mathcal{G}X),\] where the property \(\mathcal{G}=\mathcal{G}g\) holds due to the group structure.

If we use a set of transformations that is not a group, then the group invariance test can be very conservative or anti-conservative. An example is when we use a subset (that is not a subgroup) of the group of permutation maps [10]. [9, p. 822] show however that we can validly use set of transformations that is not a group, if we modify the test by applying a certain additional random transformation. This idea is extended in [11].

In §2.4.2 and §2.4.3 we saw two examples of groups of transformations, namely permutation maps and sign-flips. These are perhaps the most common examples. Another example are rotations [12], [13]. Tests based on rotations require more assumptions than tests based on sign flips, so they are used less often in practice. It can be useful to e.g. permute multiple columns of a matrix (i.e. applying the same permutation map to every column). The transformations then also form a group.

In the example from §2.4.2, the group \(\mathcal{G}\) contained \(n!=5!\) elements. Note that the size of the group quickly grows with the sample size. For example, if \(n=10\), then \(n!=3628800\) and if \(n=20\), then \(n!\approx 2.4329\cdot 10^{18}\). As a solution, researchers then do not use the full group, but use a random collection of transformations from \(\mathcal{G}\) that are uniformly sampled. The sampling can be done with or without replacement. For example, if we sample with replacement, we can proceed as follows: let \(M\) be large number (e.g. 1000) and sample \(M-1\) independent, random (uniformly distributed) transformations from \(\mathcal{G}\), then add the identity map to this collection of transformations. Then apply the group invariance test as usual, but using these \(M\) transformations rather than the full group \(\mathcal{G}\). If we sample the random transformations in this way, then (on average) the group invariance test still controls the type I error rate [9]. Tests based on random permutations or other transformations, are usually considered to be part of the class of resampling methods, just like bootstrap tests, which are discussed in §2.4.6.

2.4.5 Permutation tests in more complex models↩︎

Consider a linear model with responses \(Y\in \mathbb{R}^n\) and covariates \(X^1,Z^1,Z^2,...,Z^p\in\mathbb{R}^n\). The responses satisfy \(Y= X \beta +Z\gamma +\epsilon\), where \(Z=(Z^1,...,Z^p)\) and the residuals \(\epsilon_i\) are independent with mean 0. Suppose we are interested in testing \(\mathcal{H}_0:\beta=0\). If the classical Gaussian model holds, we can test \(\mathcal{H}_0\) in an exact way. Suppose however that the residuals are not Gaussian or are not homoscedastic. Then a permutation-type test can be a good alternative, in the sense that it may have type I error rate closer to \(\alpha\) and may have more robust power. A well-known permutation method in this situation is the Freedman-Lane method ([14], and see [15] for a review of related methods). If there are multiple response variable that we are interested in, then besides robustness, a potential advantage of Freedman-Lane is that it can be combined with powerful resampling-based multiple testing methods (see §§5, 6.3, 6.5). This means that if the residuals are perfectly Gaussian and homoscedastic, it can nevertheless be useful to apply Freedman-Lane, since it can be combined with those powerful methods.

The Freedman-Lane method is defined as follows.

1. Compute some (pivotal) test statistic \(T^1\) for testing \(\mathcal{H}_0\), e.g. the (absolute value of the) Wald statisic \(\hat{\beta}/se(\hat{\beta})\).

2. Regress \(Y\) on \(Z\), obtaining \(\hat{Y}= Z(Z'Z)^{-1}Z'Y\) and residuals \(e= Y-\hat{Y}\).

3. (Randomly) permute (or sign-flip) the residuals, obtaining new responses \(Y^*=\hat{Y}+g(e)\), where \(g(e)\) are the permuted or sign-flipped residuals. Compute a test statistic \(T^2\) in the same way as before, but using the changed responses \(Y^*\).

4. Repeat the previous step many times, obtaining test statistics \(T^3,...,T^w\).

5. The (right-sided) p-value is \(P= w^{-1}|\{1\leq j \leq w: T^j\geq T^1\}|\).

The intuition behind Freedman-Lane is as follows. Note that if \(\mathcal{H}_0\) is true and the residuals \(\epsilon_i\) are i.i.d., then \(Y^*\) has roughly the same distribution as \(Y\), so \(T^2, T^3,...\) have roughly the same distribution as \(T^1\), so that \(P\) is roughly uniformly distributed on \(\{1/w,2/w,...,w/w\}\). If \(\mathcal{H}_0\) is false, then \(T^1\) tends to be large, but after permuting the residuals, the conditional dependence between \(Y^*\) and \(X\) will be small, so \(T^2, T^3,...\) do not tend to be large. (For more theory, see [16].)

Permutation-type tests for generalized linear models are available as well, although the literature is small. See [17] for a permutation test for logistic regression. See [18] and [19] for sign-flipping-based tests for generalized linear models. These are implemented in the R package flipscores on CRAN. [20] illustrate how the tests from the latter two papers can be combined with the resampling-based multiple testing method maxT (defined in §5).

2.4.6 Bootstrap tests↩︎

The term “bootstrap method” is often used in contexts where confidence intervals are constructed based on resampling with replacement from the data. By bootstrap testing, however, we mean a type of test. Sometimes such a bootstrap test is simply equivalent to making a bootstrap confidence interval and checking whether \(\theta_0\) (e.g. 0) is within the confidence interval. However, when we combine bootstrap tests with multiple testing method such as maxT and Romano-Wolf (discussed) then we really need a bootstrap test in the strict sense, which does not involve making a confidence interval. Like bootstrapping for making confidence intervals, bootstrap tests are also based on sampling with replacement (except so-called “wild bootstrap tests”).

Like permutation testing, bootstrap testing is manner of nonparametric or semiparametric testing. Bootstrap tests are related to permutation tests in the sense that very roughly speaking, both approaches repeatedly “simulate samples from the null distribution”. This description is indeed very rough; for example, usually there is no such thing as the null distribution, since the null hypothesis typically contains many distributions.

The class of bootstrap tests is quite varied, and there is no simple overarching theory. Some bootstrap tests are so-called wild bootstrap tests [21], [22]. These are not bootstrap tests in the strictest possible sense, but are related to the permutation test for linear models which was discussed in §2.4.5. Below, we focus on “non-wild” bootstrap tests.

The differences between permutation tests (or, more generally, group invariance tests) and (non-wild) bootstrap tests are that, unlike permutation tests,

• bootstrap tests start with adjusting the data by removing the signal;

• bootstrap tests sample with replacement from this adjusted dataset;

• bootstrap tests do not e.g. permute the data.

Those are differences between the two types of tests if we see them purely as procedures. Another important difference lies in the types of hypotheses that can be tested: (non-wild) bootstrap tests are strictly speaking not used for testing a hypothesis of exchangeability, but for testing a hypothesis about a parameter. For example, if our data are \((X_1,...,X_n)\) and we want to test that \(X_1,...,X_n\) are i.i.d. (i.e. exchangeable) then we would use a permutation test. On the other hand, if we assume that \(X_1,...,X_{n/2}\) have mean \(\mu_1\) and \(X_{n/2+1},...,X_{n}\) have mean \(\mu_2\) and want to test \(\mathcal{H}_0:\mu_1=\mu_2\) then in some cases we might prefer a bootstrap test over a permutation test. Another important difference is that bootstrap tests are only asymptotically exact (i.e., the level converges to \(\alpha\) if \(n\rightarrow\infty\)), while permutation tests are sometimes exact for finite samples (as shown in §2.4.4).

Denote the data by \(X\) and suppose the null hypothesis is of the form \(\mathcal{H}_0:\theta=\theta_0\). Generally speaking, a (non-wild) bootstrap test, based on \(b\) bootstrap samples, proceeds as follows:

1. Compute a test statistic \(T(X)\), extreme values of which are evidence against \(\mathcal{H}_0\).

2. Modify the data such that the distribution of the modified data \(X^*\) “falls under \(\mathcal{H}_0\)”.

3. For all \(1\leq j \leq b\), do the following. Create a dataset \(X^j\) by sampling with replacement from \(X^*\), and compute \(T(X^j)\).

4. Reject \(\mathcal{H}_0\) if the original test statistic \(T(X)\) is extreme compared to the reference values \(T(X^1),...,T(X^b)\).

In the above, “extreme” typically means one of the following things:

• \(T(X)\) exceeds the \((1-\alpha)100\%\)-quantile of the reference distribution;

• \(T(X)\) lies below the \((\alpha/2)100\%\)-quantile or above the \((1-\alpha/2)100\%\)-quantile of the reference distribution.

As mentioned, the level of a bootstrap test is not guaranteed to be exactly \(\alpha\), but in various cases the level is \(\alpha\) asymptotically. Finally, note that we can also compute a p-value, by computing the number of reference statistics that are more extreme than the original statistic \(T(X)\) and dividing by \(b\). This p-value can be 0. We can avoid that by including the original statistic among the reference statistics. We now give two examples of bootstrap tests.

One-sample bootstrap test. Consider an i.i.d. sample \(X_1,...,X_n\), denote the mean of the observations by \(\mu\), and consider \(\mathcal{H}_0:\mu=\mu_0\). Permuting the data is not going to help us in this case. A test based on sign-flipping (see §2.4.3) can be a good choice, although is only exact if the variables have symmetric distributions. Another approach is to use the following bootstrap test.

1. Let \(T(X)\) be the usual one-sample t-statistic corresponding to \(\mathcal{H}_0:\mu=\mu_0\).

2. Let \(\hat{\mu}\) be the sample mean \(n^{-1}\sum_{i=1}^n X_i\).

3. Construct the modified dataset \(X^*:=(X_1,...,X_n) - \hat{\mu} + \mu_0\).

4. For all \(1\leq j \leq b\), do the following. Consider a bootstrap sample \(X^j\) by sampling with replacement from \(X^*\). Compute the one-sample t-statistic \(T(X^j)\) corresponding to the hypothesis that the mean is \(\mu_0\).

5. Reject \(\mathcal{H}_0\) if the original test statistic \(T(X)\) is extreme compared to the reference values \(T(X^1),...,T(X^b)\).

A different but equivalent algorithm is the following.

1. Let \(T(X)\) be the usual one-sample t-statistic corresponding to \(\mathcal{H}_0:\mu=\mu_0\).

2. Let \(\hat{\mu}\) be the sample mean \(n^{-1}\sum_{i=1}^n X_i\).

3. For all \(1\leq j \leq b\), do the following. Consider a bootstrap sample \(X^j\) by sampling with replacement from \(X\). Compute the one-sample t-statistic \(T(X^j)\) corresponding to the hypothesis that the mean is \(\hat{\mu}\).

4. Reject \(\mathcal{H}_0\) if the original test statistic \(T(X)\) is extreme compared to the reference values \(T(X^1),...,T(X^b)\).

To see that the two algorithms are equivalent, note that resampling from \(X^*\) and computing a test statistic corresponding to the hypothesis that the mean is \(\mu_0\), is equivalent to resampling directly from \(X\) and computing a test statistic corresponding to the hypothesis that the mean is \(\hat{\mu}\).

Even if \(X_1,...,X_n\) are nicely i.i.d., this bootstrap test is only asymptotically exact. If \(X_1,...,X_n\) are not i.d.d. but have different variances or shapes, then this bootstrap test can in some cases still be asymptotically exact.

Two-sample bootstrap test. Consider two i.i.d. samples \(Y_1,...,Y_{n_1}\) and \(Z_1,...,Z_{n_2}\), with means \(\mu_1\) and \(\mu_2\) respectively. Write \(X=(Y_1,...,Y_{n_1}, Z_1,...,Z_{n_2})\). Consider the null hypothesis \(\mathcal{H}_0: \mu_1=\mu_2\). It is well known that under certain assumptions, the two-sample t-test will be exact. If e.g. the normality assumption does not hold, we could consider a permutation test. For such a test to be exact, we would need that all \(n_1+n_2\) variables are identically distributed under \(\mathcal{H}_0\). If the two groups have different variances, then sometimes a permutation test based on Welch statistics [23] is more reliable. In general, if the two groups do not have the same distribution under \(\mathcal{H}_0\), the following bootstrap test may also be a good choice.

1. Let \(T(X)\) be the two-sample Welch t-statistic.

2. Compute the two sample means \(\hat{\mu}_1\) and \(\hat{\mu}_2\).

3. Let \(Y^*=Y -\hat{\mu}_1\) and \(Z^* = Z- \hat{\mu}_2\).

4. For all \(1\leq j \leq b\), do the following. Take a bootstrap sample \(Y^j\) (of length \(n_1\)) by sampling with replacement from \(Y^*\), and take a bootstrap sample \(Z^j\) (of length \(n_2\)) by sampling with replacement from \(Z^*\). Write \(X^j=(Y^j, Z^j)\). Compute the the two-sample Welch t-statistic \(T(X^j)\).

5. Reject \(\mathcal{H}_0\) if the original test statistic \(T(X)\) is extreme compared to the reference values \(T(X^1),...,T(X^b)\).

The above version of the two-sample bootstrap test samples the bootstrap samples \(Y^j\) and \(Z^j\) independently, and uses the Welch t-statistic. Such an approach is recommended in particular when the two samples cannot be assumed to have identical distributions under \(\mathcal{H}_0\). For example, sometimes the variables in the second sample have different variances than the variables in the first sample.

2.5 Exercises↩︎

1. (a) Consider a valid p-value \(P\) for a hypothesis \(\mathcal{H}_0\). As usual, we reject \(\mathcal{H}_0\) when \(P\leq \alpha\). This generally ensures that the type I error probability stays below \(\alpha\). Does this still generally hold is we choose \(\alpha\) based on the data?

   (b) Related to that, can you mention an advantage of the fact that in several scientific fields, researchers all use \(\alpha=0.05\) by default?

2. Using the fact that \(\mathcal{G}\) is a group, prove that property 2 indeed holds.

3. At the end of §2.4.4, it is explained that if we use (i.i.d. uniform) random transformations in a group invariance test, we should always add the identity map. Explain that if we do not do that, then the p-value of the group invariance test may potentially be 0 under \(\mathcal{H}_0\). Explain why that would mean that the p-value is not valid.

4. Consider the one-sample bootstrap test from §2.4.6. Suppose that the observed values \(x_1,...,x_n\) are all distinct. Suppose that the we use \(b\) bootstrap samples \(X^1,...,X^b\). We compute the p-value by taking the number of reference statistics that are larger than or equal to the original statistic \(T(X)\) and dividing by \(b\). What is the smallest possible p-value that we may get? And if we take the first bootstrap sample \(X^1\) to be the original data \(X\)?

5. \(\bigstar\) Consider the one-sample bootstrap test from §2.4.6. Suppose that the observed values \(x_1,...,x_n\) are all distinct. Recall that a bootstrap sample is a collection \(\{x_1',...,x_n'\}\) such that \(x_i'\in\{x_1,...,x_n\}\) for every \(1\leq x_i' \leq n\). Compute (in terms of \(n\)) the number of different possible bootstrap samples. To get the answer, make use of the following fact from combinatorics: the number of sequences \((b_1,...,b_n)\) such that the \(b_i\) are nonnegative integers and \(b_1+...+b_n=n\), equals \(\binom{2n-1}{n}.\)

6. In recent years, e-values have been a topic of interest in the literature on hypothesis testing. An e-value for a null hypothesis \(\mathcal{H}_0\) is defined as a nonnegative random variable \(e\) that satisfies \(\mathbb{E}(e)\leq 1\) under \(\mathcal{H}_0\). In the context of e-values, one rejects \(\mathcal{H}_0\) if \(e\geq \alpha^{-1}\). This guarantees that \(\mathbb{P}(\text{reject})\leq\alpha\) under \(\mathcal{H}_0\). Prove that.

   Further, prove that for every e-value \(e\) that takes values in \((0,\infty)\), the variable \(P:=\min\{e^{-1},1\}\) is a valid p-value.

   (e-values have certain attractive properties: they allow for online testing with optional stopping, and allow \(\alpha\) to be chosen after seeing the data while still providing error control in a particular sense.)

2.6 Solutions↩︎

1. (a) No. Suppose for example that we choose \(\alpha\) based on the p-value, namely \(\alpha:=P\). Then we always reject \(\mathcal{H}_0\). Suppose \(\mathcal{H}_0\) is true. Then clearly, conditional on the \(\alpha\) we have selected, the type I error probability is 1, which is larger than \(\alpha\) (if \(\alpha<1\)). (hemerik2026choosing? further discusses consequences of choosing \(\alpha\) based on the data.)

   (b) If everyone in a field uses \(\alpha=0.05\), an advantage is that it is clear that researchers in this field do not choose \(\alpha\) based on the data.

2. Let \(g\in\mathcal{G}\). We will show that \(\mathcal{G}\subseteq \mathcal{G}g\) and \(\mathcal{G}g \subseteq \mathcal{G}\).

   To show that \(\mathcal{G}\subseteq \mathcal{G}g\), we show that if \(h\in \mathcal{G}\) then \(h\in \mathcal{G}g\). Thus, we must show that \(h\) is of the form \(ag\), with \(a\in \mathcal{G}\). This is indeed true, since we can take \(a=hg^{-1}\) (which is contained in \(\mathcal{G}\) since \(h\) and \(g^{-1}\) are in \(\mathcal{G}\) and \(\mathcal{G}\) is a group) and then \(h=ag\).

   Next, we must show that \(\mathcal{G}g \subseteq \mathcal{G}\). This follows from the fact if \(a\in \mathcal{G}\), then \(ag\in \mathcal{G}\), by definition of a group.

3. The p-value for the group invariance test is 5 . If we use random transformations from \(\mathcal{G}\), the formula is the same, except that we replace \(\mathcal{G}\) by a random, smaller collection of transformations from \(\mathcal{G}\). Note that if we use random transformations but add the identity map, then it is guaranteed that the numerator in the formula of the p-value is at least 1, since \(T(id(X))\geq T(X)\). Otherwise the numerator may be 0, since all sampled random transformations \(g\) might satisfy \(T(gX)< T(X)\).

   If a p-value has a strictly positive probability of being 0 under \(\mathcal{H}_0\), then it does not satisfy \(\mathbb{P}_{\mathcal{H}_0}(P\leq c)\leq c\) for all \(c\in[0,1]\), so it is not valid. (A test based on such a p-value does not control the type I error rate for all \(\alpha\in(0,1)\).)

4. The smallest possible p-value is 0, since it may happen that \(T(X)\) is more extreme than the reference statistics. The answer to the second question is \(1/b\), since we always have \(T(X^1)\geq T(X)\) in the second case.

5. For each \(x_i'\) there are \(n\) possibilities, so if we would pay attention to the order of \(x_1',...,x_n'\), then there would be \(n^n\) possibilities. However, we only care about the set \(\{x_1',...,x_n'\}\) and not about the order of \(x_1',...,x_n'\). This might suggest that the answer is \(n^n/n!\). However, then we would ignore the fact that there are usually ties among \(x_1',...,x_n'.\) For example, in the extreme case that \(x_1'=...=x_n'\), there are not \(n!\) ordered samples corresponding to that sample, but only one.

   What we really need to compute, is the number of multisets of size \(n\) from a collection \(\mathcal{A}\) of \(n\) elements. A multiset is a set that can contain elements more than once. Such a multiset can be represented by a vector \((b_1,...,b_n)\), where \(b_i\) indicates the number of times that the multiset contains the \(n\)-th element of \(\mathcal{A}\). Note that such a vector satisfies \(b_1+...+b_n=n\). Hence, we can use the given fact, which means that the answer is \(\binom{2n-1}{n}.\)

   (The given fact follows from a so-called stars-and-bars reasoning: imagine we draw \(2n-1\) bars in a row (e.g. \(|||||||\) in case \(2n-1=7\).). Now turn \(n\) of the bars into stars. There are \(\binom{2n-1}{n}\) possible ways to do this. We then have \(n\) stars among \(n-1\) bars, and each such pattern corresponds 1-to-1 to a way of putting \(n\) (identical) things into \(n\) bins; indeed, we can see the \(n-1\) bars as the borders between the bins. Clearly, the number of possible ways of putting \(n\) marbles into \(n\) bins equals the number of sequences \((b_1,...,b_n)\) such that the \(b_i\) are nonnegative integers and \(b_1+...+b_n=n\).)

6. We first prove the first claim. Suppose \(\mathcal{H}_0\) holds. By Markov’s inequality, \(\mathbb{E}(e)\geq \alpha^{-1}\mathbb{P}(e\geq \alpha^{-1})\), which means that \(\mathbb{P}(\text{reject})\leq \alpha \mathbb{E}(e)\leq \alpha\).

   We now prove the second claim. Consider a valid e-value \(e\) and \(P:=1\wedge e^{-1}\). Let \(c\in(0,1]\) and suppose \(\mathcal{H}_0\) is true. We then have \(\mathbb{P}(P\leq c) = \mathbb{P}(e\geq c^{-1})\). By Markov’s inequality, \(\mathbb{E}(e)\geq c^{-1}\mathbb{P}(e\geq c^{-1})\), and it follows that \(\mathbb{P}(P\leq c)\leq \mathbb{E}(e) c\leq c\). (Clearly, \(\mathbb{P}(P\leq c)\leq c\) also holds for \(c=0\).) Thus, \(P\) is a valid p-value.

3 Multiple testing: basic concepts and methods↩︎

When we test a single hypothesis, we want to keep the probability of a false positive (type I error) small. When we test many hypotheses, we likewise need to be careful not to incur a lot of false positives. This can be done by using a multiple testing method. This section discusses basics from multiple testing theory. Many ideas in this section will be built upon in later sections.

3.1 Motivation for multiple testing methods↩︎

Suppose we wish to test \(m>1\) hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) and we compute corresponding p-values \(P_1,...,P_m\). As a simple example, assume for now that all hypotheses are true and all p-values are independent of each other and uniformly distributed on \([0,1]\). Naively we would simply compare every p-value with \(\alpha\) and reject all hypotheses with p-value below \(\alpha\). However, the probability of one or more false findings is then not below \(\alpha\) anymore. Indeed, we then have \[\mathbb{P}(\text{one or more false positives})= 1-\mathbb{P}(\text{no false positives})\] \[1-\prod_{1\leq i \leq m}\mathbb{P}(\mathcal{H}_i \text{ not rejected}) = 1-(1-\alpha)^m,\] which is larger than \(\alpha\). Indeed, it is close to 1 if \(m\) is large. Here we assumed that the p-values are independent of each other, but if they are dependent, then the probability of false positives will usually also be too high. If we would ignore this issue, then testing many hypotheses would mean that we have a high chance of rejecting one or more hypotheses, even if all hypotheses are true. This would flood the scientific literature, companies and institutions with false discoveries. To avoid this, we need multiple testing methods. Some such methods are also known as “multiple testing corrections”.

3.2 Global tests↩︎

We start with global testing, which can be seen as a simple (often too simple) approach to multiple testing, which is nevertheless insightful and serves as an important fundament for later theory. Indeed, global tests will be an essential concept in e.g. §4 and §7.

A global test is a test of the so-called global null hypothesis, which is the hypothesis that all our \(m\) hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) are true. Thus, the global hypothesis, which we will denote by \(\mathcal{H}_{\{1,...,m\}}\), is the intersection of our \(m\) hypotheses: \[\mathcal{H}_{\{1,...,m\}}= \bigcap_{1\leq i \leq m}\mathcal{H}_i.\] To understand the above notation, remember that a hypothesis is a set (of distributions), so the global null hypothesis is the intersection of sets. The global null hypothesis is the set of all distributions that are contained in each of the sets \(\mathcal{H}_1,...,\mathcal{H}_m\). Thus, the global null hypothesis is smaller (or at least not larger) than each of the individual hypotheses, i.e., for every \(1\leq i \leq m\) we have \(\mathcal{H}_{\{1,...,m\}}\subseteq \mathcal{H}_i\). A global test ensures that the probability of rejecting \(\mathcal{H}_{\{1,...,m\}}\) is at most \(\alpha\) if \(\mathcal{H}_{\{1,...,m\}}\) is true. One could debate whether a global test deserves the title “multiple testing method” or not. Indeed, a global test allows saying something about a family of hypotheses, but strictly speaking it tests a single hypothesis.

In a sense, the global null hypothesis is always easier to reject than any of the individual hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\). Indeed, if any of the hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) is false, then that logically implies that \(\mathcal{H}_{\{1,...,m\}}\) is also false. This means that in some applications, we will be able to reject \(\mathcal{H}_{\{1,...,m\}}\), even if we are not able to reject any individual hypotheses.

3.2.1 Global test by Bonferroni↩︎

The name Bonferroni is well known because of the Bonferroni multiple testing correction discussed later. However, the idea underlying Bonferroni’s method can also be used to define a global test. Bonferroni’s global test is defined as follows. Compute p-values \(P_1,...,P_m\), corresponding to the hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\). We will need that each p-value \(P_i\) is valid, in the sense that if \(\mathcal{H}_i\) is true, then \(P_i\) is uniformly distributed on \([0,1]\) — or stochastically larger than uniform. We define \(\mathcal{N}\) to be the set of the indices of the true hypotheses:

\[\mathcal{N}= \{1\leq i \leq m: \mathcal{H}_i \text{ is true}\}.\] The p-values corresponding to \(\mathcal{N}\) are sometimes called the null p-values.

We make this assumption throughout this text, unless stated otherwise. Thus, by “p-value” we will mean “valid p-value”. Note that Assumption 2 does not contain any assumption about the dependence structure of the p-values. Some methods will require such assumptions, but not Bonferroni.

The Bonferroni global test rejects the global hypothesis \(\mathcal{H}_{\{1,...,m\}}\) if and only if at least one of the p-values is smaller than or equal to \(\alpha/m\), i.e., if and only if \[\min\{P_i: 1\leq i \leq m\}\leq \alpha/ m.\]

Later we will discuss Bonferroni’s multiple testing procedure. There is a simple relationship between that procedure and Bonferroni’s global test: Bonferroni’s global test rejects \(\mathcal{H}_{\{1,...,m\}}\) if and only if Bonferroni’s multiple testing method rejects at least one of the hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\).

3.2.2 Šidák’s global test↩︎

Šidák’s global test assumes that the null p-values are independent of each other. More precisely, it is sufficient that under \(\mathcal{H}_{\{1,...,m\}}\), all p-values \(P_{1},..., P_{m}\) are independent of each other.

Let \(c\in[0,1]\). In Section 3.1 we discussed that if all p-values are independent and standard uniform, then \[\mathbb{P}(\exists 1\leq i \leq m: P_i\leq c)=\] \[1-\mathbb{P}(\forall 1\leq i \leq m: P_i> c)=\] \[1-(1-c)^m.\] It follows that that if we choose \(c\) such that \(1-(1-c)^m=\alpha\), then we have an exact global test if we reject \(\mathcal{H}_{\{1,...,m\}}\) if and only if \(\min\{P_i: 1\leq i \leq m\} \leq c.\) It turns out that this value \(c\) is \(1-(1-\alpha)^{1/m}\).

The advantage of Šidák’s global test compared to Bonferroni is that it is more powerful, since \(1-(1-\alpha)^{1/m} > \alpha/m\). However, a major caveat is that it is not valid for all possible dependence structure of the p-values. In some cases it is reasonable to assume that the p-values are independent, e.g. when the \(m\) p-values are computed based on \(m\) independent datasets. In many situations however, there is some dependence between the p-values and then Šidák’s global test may be too liberal. Roughly speaking, this can happen when there are negative correlations between the p-values.

When there are dependencies between the p-values, Šidák’s global test may still be valid. However, this is often difficult to check, for example because the dependence structure of the p-values if often unknown.

3.2.3 Global test using Fisher combinations↩︎

Another well-known global test which can be more powerful than Bonferroni, is Fisher’s combination test. The test is exact when the p-values are i.i.d. and standard uniform. An important caveat is that this test is often anti-conservative when there are positive dependencies between the p-values.

The test rejects for large values of the test statistic \[T = -2\sum_{i=1}^m \ln(P_i).\] Suppose the p-values are i.i.d. and standard uniform. Then \(T\) has a \(\chi^2_{2m}\) distribution (chi-squared distribution with 2 degrees of freedom).

Roughly speaking, when there are positive dependencies between the p-values, the global test based on Fisher combinations tends to be anti-conservative.

3.2.4 Global test using Simes’ inequality↩︎

We have seen that Bonferroni’s global test is always valid as long as the p-values are valid (Assumption 2), while the other global tests that we have discussed may be invalid when the p-values are dependent. The latter is also true for Simes’ global test, although it is in fact quite hard to find a dependence structure that invalidates Simes’ global test. Simes’ global test rejects \(\mathcal{H}_{\{1,...,m\}}\) when there is an \(1\leq i \leq m\) for which \[P_{(i)}\leq \frac{i\alpha}{m},\] where \(P_{(1)}\leq ...\leq P_{(m)}\) indicate the sorted p-values. In other words, the global test rejects when the event \[\label{eventSimes} \bigcup_{1 \leq i \leq m} \{P_{(i)}\leq \frac{i\alpha}{m} \}\tag{6}\] happens. For example, if the second smallest p-value is at most \((2\alpha)/m\), then this global test always rejects. See Figure 6 for an illustration.

Let \(Q_{1},..., Q_{N}\) be the null p-values, i.e., the p-values corresponding to the true hypotheses (in the original order) and let \(Q_{(1)}\leq...\leq Q_{(N)}\) be the corresponding sorted values. Note that Simes’s global test has type I error rate at most \(\alpha\) when the null p-values satisfy Simes’ inequality, i.e., when \[\label{SimesProbEq} \mathbb{P}\big( \bigcup_{1 \leq i \leq N} \{Q_{(i)}\leq \frac{i\alpha}{N} \}\big)\leq \alpha.\tag{7}\]

It turns out that when \(Q_{(1)},...,Q_{(N)}\) are standard uniform and independent of each other, then the probability in 7 is exactly \(\alpha\) [24]. This means that Simes’ global test then has a type I error rate of exactly \(\alpha\). Note that the same is true for Šidák’s global test and the global test based on Fisher combinations.

3.2.5 Comparison of the global tests↩︎

Here we provide a brief comparison of the global tests by Šidák, Fisher and Simes. We first compare Šidák’s global test and Simes’ global test. Both tests reject with probability exactly \(\alpha\) if the p-values are independent and standard uniform. It is clear that these tests do not always reject at the same time. Hence, it is possible that Šidák rejects and Simes not, and vice versa. Each of these two global tests rejects if one of the sorted p-values is smaller than a corresponding critical value. Simes compares each sorted p-value \(P_{(i)}\) with the critical value \(i\alpha/m\), which depends on \(i\). Šidák simply compares each sorted p-value \(P_{(i)}\) with the critical value \(1-(1-\alpha)^{1/m}\), which does not depend on \(i\). Note that the first critical value of Simes is \(\alpha/m\), which is smaller than Šidák’s critical value. However, for \(i\geq 2\), Simes’ critical value is larger. Roughly speaking, Simes’ global test has better power than Šidák’s when the signals are spread out across multiple hypotheses.

As a final rule of thumb, when many of the hypotheses are false and these signals are roughly equally strong, then the global test based on Fisher combinations can be more powerful than the global test based on Simes. However, Fisher’s global test tends to be anti-conservative when there are positive dependencies between the p-values. Simes’ global test is usually valid in practice.

3.2.6 Example of a global test for ANOVA↩︎

Consider the balanced one-way ANOVA model with \(k\) groups and \(n\) observations per group. (The approach below only works when each group has the same number of observations.) Thus, for each \(1\leq i \leq k\) and \(1\leq j \leq n\), the \(j\)-th observation in the \(i\)-th group satisfies \[Y_{ij} \sim \theta_i +\epsilon_{ij},\] where the \(\epsilon_{ij}\) have variance \(\sigma^2\). Then we can consider the global hypothesis \(\mathcal{H}_0:\theta_1=...=\theta_k\). Note that this hypothesis is the intersection of all the pairwise hypotheses of the form \(\mathcal{H}_{ii'}: \theta_i=\theta_{i'}\). The usual estimator of \(\theta_i-\theta_{i'}\) is \(\bar{Y}_{i} - \bar{Y}_{i'}\), with variance \(2\sigma^2/n\).

To test the global hypothesis \(\mathcal{H}_0\), we could use the F-test. An alternative test statistic is the maximum of the pairwise \(|t|\)-statistics \[|T_{ii'}|= |\bar{Y}_{i} - \bar{Y}_{i'}|/(\hat{\sigma}\sqrt{2/n}),\] where \(\hat{\sigma}\) is the sample standard deviation, i.e., \(\hat{\sigma}^2= \sum_{i=1}^k \sum_{j=1}^{n} (Y_{ij}-\bar{Y}_{i})^2/(nk-k).\) Thus, we consider \[\sqrt{2}\max_{1\leq i < i'\leq k} |T_{ii'}| = \max_{1\leq i < i'\leq k}\frac{|\bar{Y}_i-\bar{Y}_{i'}|}{\hat{\sigma}/\sqrt{n}},\] which equals \[\label{eqttTW} \frac{\max_{1\leq i \leq k} \bar{Y}_i- \min_{1\leq i \leq k}\bar{Y}_{i}}{\hat{\sigma}/\sqrt{n}}.\tag{8}\]

Under \(\mathcal{H}_0\), this statistic is distributed as the range of \(k\) independent standard normal random variables divided by an independent \(\sqrt{\chi^2_{\nu}/{\nu}}\) random variable, where \(\nu=k(n-1)\) [25, p. 31]. This variable is denoted by \(Q_{k,\nu}\) and is called the Studentized range variable with parameter \(k\) and \(\nu\) degrees of freedom for the CDF, see [25, p. 376].

We reject the global null hypothesis when \[\sqrt{2}\max_{1\leq i < i'\leq k} |T_{ii'}| \geq Q_{k,\nu}^{(1-\alpha)},\] where \(Q_{k,\nu}^{(1-\alpha)}\) is the \((1-\alpha)\)-quantile of \(Q_{k,\nu}\). This test has better power than the F-test in some situations with sparse signal, and can be computationally convenient in a multiple testing context, as will be discussed in §4.5.

3.3 Familywise error rate (FWER) control↩︎

Earlier we discussed the probability of incurring one or more false positives, when testing multiple hypotheses. This probability is called the familywise error rate (FWER). Thus, if \(V\) is the number of false positives (incorrect rejections), then \[FWER=\mathbb{P}(V>0).\] Note that \(V\) is the number of true hypotheses that are rejected, so \[V=|\mathcal{N}\cap \mathcal{R}|,\] where \(\mathcal{N}\subseteq\{1,...,m\}\) is the set of indices of true hypotheses and \(\mathcal{R}\subseteq\{1,...,m\}\) is the set of indices of rejected hypotheses. If \(FWER\leq \alpha\), we say that the FWER is controlled at level \(\alpha\). Note that if we use a multiple testing method that controls the FWER, then we know that with probability at least \(1-\alpha\), all rejected hypotheses are true discoveries. (This does not mean that conditional on rejecting something, this probability is at least \(1-\alpha\). For example, suppose that all hypotheses are true. Then conditional on \(|\mathcal{R}|>0\), the probability that \(V>0\) is 1.) In §§3.4, 3.6, 3.7 and several later sections, we discuss examples of FWER controlling methods.

Sometimes a distinction is made between weak FWER control and strong FWER control. The former means that the FWER is controlled when all hypotheses are true, but not necessarily in other cases. Thus, such methods provide no guarantees, except that likely nothing is rejected when all hypotheses are true. Strong FWER control means that the FWER is controlled, regardless of how many hypotheses are false. If we simply say “FWER control”, then this is understood to mean strong FWER control.

When we test a single hypothesis, the power of a test is simply defined as the probability of rejecting the hypothesis. When we test multiple hypotheses, we need a new notion of power. Let \(\mathcal{N}^c=\{1,...,m\}\setminus\mathcal{N}\). By the power of a multiple testing method, we will usually mean \[\mathbb{E}\frac{|\mathcal{N}^c\cap \mathcal{R}|}{|\mathcal{N}^c|},\] which is the expected fraction of all false hypotheses that is rejected. The power might also be defined as the probability of rejecting at least one false hypothesis, but that is less common in the literature. We say that procedure A is uniformly more powerful than procedure B if A always rejects all the hypotheses that B rejects, and possibly more.

3.4 Bonferroni↩︎

3.4.1 FWER control with Bonferroni↩︎

We saw that if all hypotheses are true and we only reject the hypotheses \(\mathcal{H}_i\) with \(P_i\leq \alpha/m\), then the probability of any false positives is at most \(\alpha\). It can be easily seen that this is still true if some of the hypotheses are false. This leads us to the Bonferroni method, which is defined as the procedure that rejects all hypotheses with \(\mathcal{H}_i\) with \[P_i\leq \alpha/m.\]

If \(\pi_0\) is near 1, and the null p-values are standard uniform and independent or negatively dependent, then the power of Bonferroni’s method is not so far from optimal. However, it follows from Proposition 5 that the FWER of Bonferroni is much smaller than \(\alpha\) if \(\pi_0\) is much smaller than 1. This means that Bonferroni is then rather suboptima. Thus, if we would know \(\pi_0\), we could use that knowledge to make the method more powerful (how?); however, in practice we often do not know \(\pi_0\).

When there are many strong positive dependencies between the p-values, then the FWER of Bonferroni is also much smaller than \(\alpha\) (even if all hypotheses are true). The reason is that under strong positive dependence, the distribution of \(\min\{P_i:\in \mathcal{N}\}\) is stochastically larger than when the p-values are independent or negatively dependent. In that case, we would still get valid FWER control if we used a rejection threshold that is much larger than \(\alpha/m\). Thus, in a sense Bonferroni then also has rather suboptimal power. However, in practice we often do not know the dependence structure of the p-values.

We now consider two simple examples and discuss them through the lens of Bonferroni’s method.

In case we know that the null p-values are independent of each other, in Bonferroni’s method we can replace \(\alpha/m\) by \(c\) from Šidák’s global test. We then get a procedure that controls the FWER and is slightly more powerful than Bonferroni.

Bonferroni’s method is quite popular. This is because it is simple and valid regardless of the dependence structure of the p-values. However, Bonferroni’s method is uniformly improved by the slightly more complicated Holm procedure (also known as Bonferroni-Holm), see §3.6.

When there are positive dependencies between the p-values, it may be more powerful to use a resampling-based multiple testing method such as the maxT method, which will be discussed in §5.

3.4.2 Simultaneous confidence intervals↩︎

Suppose we consider \(m\) parameters \(\theta_1,...,\theta_1\in \mathbb{R}\) and we have a method that provides a \((1-\alpha)100\%\)-confidence interval (CI) for each of these parameters, for any chosen \(\alpha\in(0,1)\). Note that each CI contains its true parameter with probability (at least) \(1-\alpha\), but the probability that they all contain their true parameters is usually smaller than \(1-\alpha\), which is often not satisfactory. Instead, we can consider simultaneous CIs, i.e., intervals which are such that with probability (at least) \(1-\alpha\), all intervals contain their parameters. A highly popular way of constructing such CIs is related to Bonferroni’s method: we simply change \(\alpha\) to \(\alpha/m\), i.e., we make \(m\) \((1-\alpha/m)\)-CIs. The probability that all CIs contain their parameters is then \[\mathbb{P}(\text{all CIs contain their parameter})=\] \[1- \mathbb{P}(\text{some CIs don't contain their parameter}) \geq 1- m\cdot \alpha/m=1-\alpha,\] where we used Boole’s inequality again. If the CIs are independent of each other, we even known that \[\mathbb{P}(\text{all CIs contain their parameter}) \geq (1-\alpha/m)^m>1-\alpha.\] This allows us to make slightly smaller CIs (use \(c\) from Šidák’s global test instead of \(\alpha/m\)) and still guarantee that they are are simultaneously valid. However, in many cases we cannot assume that the CIs are independent.

Suppose the hypotheses are of the form \(\mathcal{H}_i: \theta_i=\theta_i^0\), \(1\leq i \leq m\), and we compute \(m\) \((1-\alpha/m)\)-CIs by inverting \(m\) t-tests. Then we can control the FWER by rejecting the hypotheses \(\mathcal{H}_i\) corresponding to the CI that do not contain \(\theta_i^0\). This is the same as applying Bonferroni to the p-values corresponding to the t-tests. In case we use the slightly narrower CIs using Šidák’s \(c\), we get a similar relationship with Sidák’s global test. In these examples, we see that there is a simple correspondence between a multiple testing method and a method for constructing simultaneous confidence intervals. However, for most other multiple testing methods, there is no such correspondence.

3.5 Step-up and step-down methods↩︎

In §3.6, we will consider Holm’s procedure [26], also known as Bonferroni-Holm. This is a FWER controlling procedure that is closely related to Bonferroni. Like Bonferroni, Holm controls the FWER when the p-values are valid (Assumption 2). Holm rejects all hypotheses that Bonferroni rejects, but may reject additional hypotheses on top of that. This means that Holm is a uniform improvement of Bonferroni. We say that procedure A uniformly improves procedure B if A (is still valid and) rejects all hypotheses that B rejects and potentially more. Holm is a sequential rejection procedure, because it is a stepwise procedure, that rejects more and more hypotheses until it stops.

More precisely, Holm belongs to the class of step-down procedures. Step-down procedures are often based on p-values, but some step-down procedures are based on test statistics and do not use p-values. For example, the sequential maxT method from §5 is a step-down method which is usually based on test statistics. For now, we discuss procedures based on p-values. Step-down methods based on p-values first sort the p-values from the smallest to the largest. Throughout, we denote the sorted p-values by \(p_{(1)}\leq ...\leq p_{(m)}\). Starting from the smallest p-value, step-down methods check for each p-value \(p_{(i)}\) whether it is smaller than some critical value \(c_i\). They stop as soon as a p-value \(p_{(i)}\) is reached that satisfies \(p_{(i)}>c_i\); Then, all, hypotheses with p-values strictly smaller than \(p_{(i)}\) are rejected, see Algorithm 7.

There are also step-up procedures, for example the Benjamini-Hochberg method, which is discussed later. Step-up methods find the largest \(i\) for which \(p_{(i)}\leq c_i\) and reject all hypotheses with p-values that are at most \(p_{(i)}\). An efficient way to find this index \(i\) is to start with the largest p-value and work backwards, see Algorithm 8.

The terminology step-down and step-up is potentially confusing. Indeed, a step-down method for example, actually steps up: it starts with the smallest p-value and then steps to the next one. The terminology often makes sense however when we reformulate the methods in terms of test statistics: a step-down method often starts with the largest test statistic and then steps down to the second-largest, etcetera.

3.6 FWER control with Holm’s method↩︎

The first step of the step-down method of Holm is identical to Bonferroni, i.e., all p-values are compared with \(\alpha/m\). If no hypotheses are rejected in the first step, the method stops. If \(r\geq 1\) hypotheses are rejected in the first step, then the method moves on to the second step, where all remaining p-values are compared with the threshold \(\alpha/(m-r)\). We will refer to this as step 2 in the proof below. If this leads to no additional rejections then stop. Otherwise, move to step 3, where all remaining p-values are compared with \(\alpha/(m-r)\), where \(r\) is the total number of hypotheses that have been rejected in the previous steps. The method continues like this until a step is reached where no additional hypotheses are rejected, or when all hypotheses have been rejected.

A different way to formulate Holm’s method is the following.

1. Consider the sorted p-values \(p_{(1)}\leq ...\leq p_{(m)}\) and find the smallest \(1\leq k \leq m\) for which \[p_{(k)}>\frac{\alpha}{m-k+1},\] if there is such a \(k\), otherwise take \(k=m+1\).

2. Then reject the hypotheses corresponding to the \(p_{(i)}\) with \(i<k\).

Note that this method is simply the step-down method as defined in §3.5, with critical values \[c_j=\frac{\alpha}{m-j+1}.\]

Other examples of step-down methods are the maxT method and the Romano-Wolf method for FDX control, which will be discussed later. The method in §3.7 is also a step-down method (if we use that term in a somewhat general sense). Holm is an example of a closed testing method, as will also be explained later. Clearly, Holm is uniformly more powerful than Bonferroni. The gain in power is the largest when \(\pi_0\) is far below 1 and the signals are strong. Indeed, in that case many hypotheses tend to be rejected in the first step of Holm, which leads to a substantial increase in the critical value, compared to the initial critical value \(\alpha/m.\)

We mentioned that if there are many strong positive dependencies between the p-values, then Bonferroni is rather conservative. The same is true about Holm. When there are such dependencies, the FWER of Holm is (much) smaller than \(\alpha\). When there are no dependencies, then Holm has fairly good power — among the class of methods that take a vector of \(m\) p-values as input and use no further information. Of course, in many cases, we do not know a priori what the dependence structure of the p-values is. If we can make no assumptions on the dependence structure, then Holm is often a good choice — among the class of methods that take a vector of \(m\) p-values as input and are valid under general dependence and use no further information. One may wonder whether Holm is admissible. By “admissible” we mean that the procedure cannot be uniformly improved without invalidating the method.² There is no known proof that Holm is admissible but it has been conjectured that it is. (It can be proved that Holm admissible in a somewhat different sense [27], but we will not consider that concept here.)

3.7 FWER control for ANOVA↩︎

Here we consider an important but non-straightforward FWER controlling method for one-way ANOVA models, for example the balanced one-way ANOVA model from §3.2.6. The procedure is often referred to as the Tukey-Welsch procedure (see [25, p. 69], [28]). (We do not prove here that this procedure controls the FWER.) We could also use Holm (applied to the p-values from the pairwise comparisons) for example, but it will have low power in this situation, due to dependencies between the pairwise comparisons.

For each subset \(\mathcal{A}\) of \(\{1,...,k\}\), define \(\mathcal{H}_{\mathcal{A}}\) as the hypothesis that all \(\theta_i\) with \(i\in \mathcal{A}\) are equal. (In later sections that are not about ANOVA, the notation \(\mathcal{H}_{\mathcal{A}}\) will have a different meaning.) The goal here is to test \(\mathcal{H}_{\mathcal{A}}\) for every \(\mathcal{A}\subseteq \{1,...,m\}\) with \(\mathcal{A}\geq 2\).

The procedure proceeds in a particular type of “step-down” manner. We define \(\alpha_p=\alpha\) for \(p\in \{k,k-1\}\), and \(\alpha_p = 1-(1-\alpha)^{p/k}\) for \(2\leq p\leq k-2\). Note that as \(p\) decreases, \((1-\alpha)^{p/k}\) increases, so \(\alpha_p\) decreases. Below we give a specific version of the procedure for the one-way balanced ANOVA model. For certain other ANOVA models and test statistics (e.g. F-statistics) the method below also provides FWER control (with the same definition of \(\alpha_p\), but the critical values should be adjusted to the test statistics used). We now give the steps of the Tukey-Welsch procedure.

1. First test \(\mathcal{H}_{\{1,...,k\}}\) at level \(\alpha\) using the test from §3.2.6. If it is not rejected, stop. Otherwise, continue.

2. Test all hypotheses \(\mathcal{H}_{\mathcal{A}}\) with \(\mathcal{A}\subset\{1,...,k\}\) and \(|\mathcal{A}|=m-1\) using the test statistic 8 from §3.2.6 applied to the data corresponding to \(|\mathcal{A}|\), \[\frac{\max_{i\in \mathcal{A}} \bar{Y}_i- \min_{i\in \mathcal{A}}\bar{Y}_{i'}}{\hat{\sigma}/\sqrt{n}}.\] where \(\hat{\sigma}\) still denotes the sample SD based on the full dataset, i.e., \[\hat{\sigma}^2=\sum_{i=1}^k \sum_{j=1}^{n} (Y_{ij}-\bar{Y}_{i})^2/(nk-k),\] and using critical value \(Q_{|\mathcal{A}|,\nu}^{(1-\alpha_p)}\), with \(\nu=k(n-1)\) and \(p=|\mathcal{A}|\). If they are all not rejected, stop. Otherwise continue.

3. For all \(\mathcal{A}\) that were rejected in the previous step, test all hypotheses \(\mathcal{H}_{\mathcal{B}}\) with \(\mathcal{B}\subset\mathcal{A}\) and \(|\mathcal{B}|=m-2\) using the test from §3.2.6 (applied to the data corresponding to \(|\mathcal{B}|\), but using \(\hat{\sigma}\) based on the full dataset, and using critical value \(Q_{|\mathcal{B}|,\nu}^{(1-\alpha_p)}\), with \(\nu=k(n-1)\) and \(p=|\mathcal{B}|\)). If nothing is rejected in this step, stop. Otherwise continue.

4. Consider all sets \(\mathcal{B}\) that were rejected in the previous step and consider their subsets of size \(m-3\). Keep on going like this, until no further subsets remain to be tested.

The method is implemented in the function tukeyWelsch() of the R package asht on CRAN. Some other methods related to Tukey-Welsch are discussed in §4.5.

3.8 k-FWER control↩︎

When we test many hypotheses, often we are willing to have a somewhat higher risk of false positives if this leads to much better power. (Recall that by power, we usually mean the expected fraction of the false hypotheses that is rejected.) One way to increase the power, is to control the k-familywise error rate (k-FWER), where \(1\leq k\leq m\). The k-FWER is defined as follows: \[k\text{-FWER}=\mathbb{P}(\text{reject at least } k \text{ true hypothess}).\] Formulated differently, \[k\text{-FWER}=\mathbb{P}(V\geq k).\] Note that for \(k=1\), this is the usual FWER, Thus, the k-FWER is a generalization of the FWER. Controlling the k-FWER means ensuring that \(k\text{-FWER}\leq \alpha.\) A simple k-FWER controlling procedure is the following generalization of Bonferroni: reject all hypotheses \(\mathcal{H}_j\) with p-values \(P_j\) that satisfy \[\label{eq:thresholdkFWER} P_j\leq k\alpha/m.\tag{10}\]

In the above proof, Markov’s inequality is used. This inequality may be quite conservative; however, it can be shown that there are settings where (Assumption 2 is satisfied and) the k-FWER is exactly \(\alpha\) with this method [29] Thm 2.1(ii).

A different k-FWER method is defined as follows:

1. Reject all hypotheses that Bonferroni (or Bonferroni-Holm) rejects.

2. Consider the p-values that were not rejected in step one, and additionally reject the hypotheses corresponding to the smallest \(k-1\) of these p-values (or reject all \(m\) hypotheses if there were fewer than \(k-1\) left).

This method also controls the k-FWER, but it potentially rejects hypotheses with large p-values, so this method is not recommended.

Clearly, for \(k=1\), the method from Theorem 7 reduces to Bonferroni. Since Bonferroni is uniformly improved by Holm, one may wonder whether the k-FWER method above can likewise be improved. This is indeed the case. For every \(1\leq i \leq m\), consider the following critical value:

\[\label{cvkFWERstepdown} c_i= \begin{cases} \frac{k\alpha}{m}& i\leq k, \\ \frac{k\alpha}{m+k-i} & i>k, \end{cases}\tag{11}\] These critical values are a generalization of the Holm-critical values. Indeed, for \(k=1\), we get Holm’s critical values.

The downside of k-FWER methods compared to FWER control, is that we no longer have \((1-\alpha)100\%\) confidence that there are no false positives. For example, suppose \(k=4\) and we reject 5 hypotheses. Then we have no guarantee that it is likely that most of these findings are true discoveries.

Note that if we do not require FWER control but k-FWER control with \(k\geq 2\), then we could in principle always simply reject \(k-1\) arbitrarily chosen hypotheses (even if all p-values are very large) and still control the k-FWER. Thus we can uniformly improve the method defined by 10 by rejecting some additional hypotheses if that method rejects fewer than \(k-1\) hypotheses. However, if we reject \(k-1\) (or fewer) hypotheses with a k-FWER method, then that is meaningless.

At the beginning of this section (§3.8), we discussed that k-FWER control can be useful in case we want to have higher power. However, instead of k-FWER control, we can then also use a different approach, such as false discovery exceedance (FDX) control (see §6) or false discovery rate control (see §6.1, §8, §9). Such approaches are often more appealing than k-FWER control. Some FDX controlling methods rely on iteratively using k-FWER mehods (see §6.3).

3.9 Adjusted p-values↩︎

Consider hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) and some multiple testing procedure (MTP) which, given \(\alpha\in(0,1)\), decides for each hypothesis whether it is rejected or retained. For every \(1\leq i \leq m\), the adjusted p-value \(p^{\text{adj}}_i\) is defined to be the smallest value \(\alpha\in[0,1]\) for which the MTP still rejects \(\mathcal{H}_i\). More precisely, \[p^{\text{adj}}_i = \min\{\alpha: \mathcal{H}_i \text{ is rejected by the MTP}\}\wedge 1 ,\] where \(a\wedge b\) denotes the minimum of \(a\) and \(b\). For example, note that for Bonferroni, \(p^{\text{adj}}_i = \min\{p_i\cdot m,1\}\). Clearly, \(p_i\cdot m\) can be larger than 1. This illustrates why authors often take the minimum with 1, although this does not really matter, since we do not reject anyway, if \(\min\{\alpha: \mathcal{H}_i \text{ is rejected by the MTP}\}\) is larger than or equal to 1.

To obtain the adjusted p-values for step-down or step-up methods based on p-values, a general algorithm can be used. Consider a step-up or step-down method with critical values \(c_1,...,c_m\). The critical values are typically of the form \(\alpha/a_i\), where \(a_i\) is an adjustment factor which does not depend on \(\alpha\). Thus, \(a_i\) is the ratio between \(\alpha\) and the \(i\)-th critical value. The following algorithm can be used for obtaining the adjusted p-values [30]:

1. Obtain the sorted p-values \(p_{(1)},...,p_{(m)}.\)

2. Multiply each sorted p-values \(p_{(i)}\) by the adjustment factor \(a_i\).

3. If this multiplication violates the ordering, adjust the values. For step-down methods, take \(p^{\text{adj}}_{(i)} = \max_{1\leq j \leq i}a_jp_{(j)}\). For step-up methods, take \(p^{\text{adj}}_{(i)} = \min_{i\leq j \leq m}a_jp_{(j)}\).

4. Set \(p^{\text{adj}}_{(i)}=\min(p^{\text{adj}}_{(i)},1)\) for all \(i\).

For example, for Holm, which is a step-down method, the critical values are \(c_i=\alpha/(m+1-i)\). Thus, to get the adjusted p-values for Holm, we first compute the values \(p_{(i)} \alpha c_i^{-1} = p_{(i)}(m+1-i)\) and then possibly adjust them upwards to enforce monotonicity.

The R function p.adjust() outputs adjusted p-values for some FWER and FDR controlling methods, including Holm, Hommel and Benjamini-Hochberg (the latter two will be discussed later). If one has computed adjusted p-values, then for any value of \(\alpha\), one can immediately see which hypotheses are rejected by the MTP, namely the hypotheses with adjusted p-values below \(\alpha\). Since the function p.adjust() outputs adjusted p-values, it is not necessary to provide \(\alpha\) as an argument and nevertheless, based on the adjusted p-values, one can easily see which hypotheses are rejected for any value of \(\alpha\).

There are also some important caveats when using adjusted p-values. The first one is that just like with regular p-values, it is required to specify \(\alpha\) independently from the data, if we want FWER control (or e.g. FDR control, which is discussed later). Choosing \(\alpha\) based on the data or based on the adjusted p-values, may invalidate a MTP. Also note that, as is clear from the algorithm above, the adjusted p-value of a hypothesis may depend on the other p-values. Thus, adjusted p-values have no independent meaning.

3.10 Application to property sales data↩︎

To illustrate the methods, we apply them to a dataset on residential property sales in Ames, Iowa from 2006 to 2010 [31]. We use the cleaned data which are available through the R package AmesHousing [32]. The dataset contains 2930 rows, which correspond to sold individual properties. The data contain the sale price and various variables that might affect it, such as the lot area, the number of bathrooms above grade, the build year and the roof style.

We removed variables that were near-constant and standardized all quantitative variables, including the sale price. We fitted a standard least squares linear regression model, where sale price was regressed on all other variables, including the dummy-encoded categorical variables. There were 31 numerical predictors and 248 dummy variables, so in total 279 predictors. For each predictor, we considered the hypothesis \(\mathcal{H}_i\) that this predictor has coefficient 0. We computed the 279 two-sided p-values using the usual t-statistics. The smallest p-values are shown in Figure 9.

We take \(\alpha=0.05\). Recall that Bonferroni rejects all hypotheses with p-values below \(\alpha/m=\alpha/279\approx 0.0001792\). In this case, Bonferroni, rejects 36 hypotheses. Indeed, the 36-th smallest p-value, \(p_{(36)}\), is \(\approx 0.0000997\) and \(p_{(37)}\approx 0.0002432\); clearly \(0.0000997\leq \alpha/m\) and \(0.0002432>\alpha/m\).

We now consider Holm. In the first step Holm rejects the same 36 hypotheses as Bonferroni, by definition. In the next step, Holm increases the rejection threshold to \(\alpha/(m-36)=0.05/243\approx 0.0002057613\). This is still smaller than \(p_{(37)}\), hence Holm rejects no further hypotheses in this case.

Next, we apply the k-FWER methods from §3.8. We take \(k=5\) and apply the single-step method from Theorem 7. This leads to 37 rejected hypotheses. Note that in this particular example, the result that Bonferroni provides is more attractive, since it rejects one hypothesis less but it provides much stronger guarantees. The step-down k-FWER method (from Theorem 8) also rejects 37 hypotheses.

We now focus on the p-values for the 31 numerical predictors, see Figure 10. With \(\alpha=0.05\) as before, Bonferroni and Holm both reject 10 hypotheses. The single-step 5-FWER method rejects 12 hypotheses and the sequential 5-FWER method (from Theorem 8) rejects 14 hypotheses. Note that the 5-FWER methods again provide no additional information compared to Holm, in this case.

3.11 Exercises↩︎

1. Consider 5 hypotheses and corresponding p-values \(0.025\), \(0.045\), \(0.010\), \(0.012\) and \(0.245\). Take \(\alpha=0.05\). How many hypotheses are rejected by Holm’s method? And by the k-FWER method from Theorem 8, with \(k=2\)?

2. Consider Bonferroni’s FWER controlling method applied to \(m\geq 2\) valid p-values, which are independent of each other. Prove that the FWER is then strictly smaller than \(\alpha\), regardless of how many hypotheses are true.

3. Let \(\alpha=0.05\) and (using a calculator) compute the ratio between the critical value of Bonferroni’s global test (\(\alpha\)/m) and of Šidák’s global test (\(1-(1-\alpha)^{1/m}\)) for \(m\in\{2,3,4,5\}\). What (simple) pattern do you notice?

4. \(\bigstar\) In exercise [ex:bonsidak], we considered the ratio between the critical values of Bonferroni and Šidák. Prove that as \(m\rightarrow\infty\), this ratio converges to \(-\alpha/\log(1-\alpha)\).

5. Consider \(m\geq 2\) hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) and \(\alpha\in(0,1)\). Suppose all hypotheses are true and all p-values have the uniform distribution on \((0,1]\). Show that under a certain dependence structure of the p-values, the FWER of Bonferroni is exactly \(\alpha\).

6. Assume the p-values \(P_1,...,P_m\) are valid (Assumption 2). Prove that Bonferroni’s FWER-controlling method ensures that the expected value of the number of false positives is at most \(\alpha\), i.e., that \(\mathbb{E}(V)\leq \alpha\). Show that this does not necessarily hold for Holm’s method.

7. In §3.9 a procedure is provided for computing adjusted p-values for step-down (and step-up) methods based on p-values and fixed critical constants. For the case of step-down methods, prove that this algorithm indeed provides the adjusted p-values \(p_{(i)}=\min\{\alpha\in[0,1]: \mathcal{H}_i \text{ is rejected by the MTP}\}\wedge 1.\)

3.12 Solutions↩︎

1. The sorted p-values are \(0.010\), \(0.012\), \(0.025\), \(0.045\) and \(0.245\) and Holm’s critical values are \(c_1=\alpha/5=0.01\), \(c_2=\alpha/4=0.0125\), \(c_3=\alpha/3=0.0166...\), \(c_4=\alpha/2=0.025\), \(c_5=0.05\). Clearly \(P_{(i)}\leq c_i\) for \(1\leq i \leq 2\), but not for \(i=3\). Hence, 2 hypotheses are rejected.

   The 2-FWER method, which is also a step-down method, always rejects at least as many hypotheses as Holm. For the 2-FWER method, the critical values are

   \[c_i= \begin{cases} \frac{2\alpha}{m}& i\leq 2, \\ \frac{2\alpha}{m+2-i} & i>2, \end{cases}\] Thus, \(c_3=2\alpha/4=0.025\), \(c_4=2\alpha/3=0.0333...\) and \(c_5=\alpha\). Note that \(P_{(3)}=0.025\leq c_3.\) Note that \(P_{(4)}=0.045>c_4\). Thus, the 2-FWER method rejects 3 hypotheses.

2. The FWER is \[\mathbb{P}\big(\bigcup_{i\in \mathcal{N}}\{P_i\leq \alpha/m\}\big).\] Note that this is strictly smaller than \(\sum_{i\in \mathcal{N}} \mathbb{P}(P_i\leq \alpha/m),\) since the events \(\{ P_i\leq \alpha/m \}\), \(i \in \mathcal{N}\) are overlapping, due to the assumed independence. This sum equals \(\alpha(N/m)\leq \alpha\). Thus, \(FWER<\alpha\).

3. Approximately 1, 0.987340, 0.983143, 0.981050 and 0.979795. The relative difference becomes larger as \(m\) increases. (Apparently, as \(m\) increases, Bonferroni “pays” more for allowing for negative dependencies. It can be shown that the ratio decreases to \(-\alpha/\log(1-\alpha)\) as \(m\rightarrow \infty\). This is approximately 0.974786.)

4. Note that the ratio is \[\frac{\alpha/m}{1-(1-\alpha)^{1/m}}=\] \[\frac{\alpha/m}{1-\exp((1/m)\log(1-\alpha))}=\]

   \[\frac{\alpha/m}{ 1 - \big( \frac{(1/m\log(1-\alpha))^0}{0!} + \frac{(1/m\log(1-\alpha))^1}{1!} + \frac{(1/m\log(1-\alpha))^2}{2!}...\big) }=\] \[\frac{\alpha}{ m - \big( m + \frac{\log(1-\alpha)^1}{1!} + \frac{1/m(\log(1-\alpha))^2}{2!}...\big) } =\]

   \[\frac{\alpha}{ - \frac{\log(1-\alpha)^1}{1!} - \frac{1/m(\log(1-\alpha))^2}{2!}... }=\] \[\frac{\alpha}{ - \frac{\log(1-\alpha)^1}{1!} +o(m) }.\] As \(m\rightarrow\infty\), this converges to \(\alpha/(-\log(1-\alpha))\).

5. We provide two possible solutions. (For the second approach, no complete proof is provided here.)

   First solution: Let \(i_1,...,i_m\) be a random permutation of \(1,...,m\) (uniformly distributed on all possible permutations) and for every \(1\leq j \leq m\) let \(P_{i_j}\) be uniformly distributed on \(((j-1)/m, j/m]\). Clearly, conditionally on \(i_1,...,i_m\) these p-values do not have the standard uniform distribution, but marginally (i.e. unconditionally) they do. To see this, consider any \(1\leq i \leq m\) and note that for every \(1\leq j \leq m\), with probability \(1/m\), the conditional distribution of \(P_i\) is uniform on \(((j-1)/m, j/m]\); hence marginally, \(P_i\) is uniform on \((0,1]\).

   We now show that for this dependence structure, the FWER of Bonferroni is \(\alpha\). It suffices to show that conditionally on \(i_1,...,i_m\), the probability that Bonferroni rejects something is \(\alpha\). To see this, note that conditionally on \(i_1,...,i_m\), Bonferroni will reject at most one hypothesis, which is \(\mathcal{H}_{i_1}\), and this happens when \(P_{i_1}\leq \alpha/m\), which holds with probability \(\alpha\).

   Second solution: Let \(P_1\) be uniformly distributed on \([0,1]\) and for every \(2\leq i \leq m\) let \(P_i=P_1+(i-1)/m - \lfloor P_1+(i-1)/m \rfloor\) (where \(\lfloor P_1+(i-1)/m \rfloor\) is the “integer part” of \(P_1+(i-1)/m\)). All the p-values are then standard uniform, and neighbouring p-values are a distance \(1/m\) apart.

   It can be checked that for this dependence structure of the p-values, the minimum of the p-values is uniform on \([0,1/m]\), so it is below \(\alpha/m\) with probability \(\alpha\), so the FWER of Bonferroni is then also exactly \(\alpha\).

6. With Bonferroni’s method, it holds that \[\mathbb{E}(V) = \mathbb{E}(\sum_{i\in \mathcal{N}}\mathbb{1}_{\{P_i\leq \alpha m\}})= \sum_{i\in \mathcal{N}} \mathbb{E}\mathbb{1}_{\{P_i\leq \alpha m\}}=\sum_{i\in \mathcal{N}} \mathbb{P}\{P_i\leq \alpha m\} \leq N\alpha m\leq \alpha.\]

   Note that if all hypotheses are true and all p-values are standard uniform, we get an equality: \(\mathbb{E}(V) = \alpha\). In that setting (i.e. \(m=N\) and all p-values standard uniform) Holm always rejects the hypotheses that Bonferroni rejects, and sometimes more (since Holm uses larger critical values \(c_i\) than Bonferroni). Hence, it then holds that \(\mathbb{E}(V) > \alpha\) for Holm.

7. It suffices to show that if \(\mathcal{H}_{(i)}\) is rejected, then \(p^{\text{adj}}_{(i)}\leq \alpha\), and vice versa. Indeed, that implies that \(p^{\text{adj}}_{(i)}= \alpha\) is the borderline case where \(\mathcal{H}_{(i)}\) is barely rejected — in other words, \(p^{\text{adj}}_{(i)}\) is the smallest \(\alpha\) for which the procedure rejects \(\mathcal{H}_{(i)}\).

   Note that \(\mathcal{H}_{(i)}\) is rejected if and only if \(p_{(j)}\leq c_j\) for all \(j\leq i\), which is equivalent to \(\max_{1\leq j \leq i}c_j^{-1}p_{(j)}\leq 1\), which is equivalent to \(\max_{1\leq j \leq i}a_i p_{(j)}\leq \alpha\), which is equivalent to \(p^{\text{adj}}_{(i)}\leq \alpha.\) This finishes the proof.

4 Closed testing for FWER control↩︎

In §3 we discussed some simple methods that control the FWER. One may wonder how we can construct a FWER controlling method that is as powerful as possible. However, one should realize that given some dataset, it is often not known which method will be most powerful. Moreover, not all methods require the same assumptions. Thus, it often depends on the situation which FWER method can best be used. Hence, it is useful to know several methods, each of which can be a good choice in certain specific situations.

Closed testing is a general recipe for constructing multiple testing procedures that control the FWER or provide confidence on false discovery proportions. Moreover, it turns out that all admissible FWER controlling methods are closed testing procedures (CTPs); recall that “admissible” means that the procedure cannot be uniformly improved without violating FWER control. This follows from the fact that a multiple testing procedure that provides FWER control, is a CTP or can be uniformly improved by one (see additional exercises). In §7, CTPs also play an important role, although there we do not focus on the FWER, but on false discovery proportions.

The ingredients required for creating a CTP are so-called local tests. Once we have defined the local tests, the CTP that controls the FWER follows from them. Hence, we will first discuss local tests and then closed testing.

4.1 Local tests↩︎

Let \(\mathcal{C}\) be the set of all nonempty subsets of \(\{1,...,m\}\). Note that \(|\mathcal{C}| = 2^m-1\). When we consider hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\), we can also consider the corresponding intersection hypotheses, which are all hypotheses \(\mathcal{H}_{\mathcal{I}}\) with \(\mathcal{I}\in \mathcal{C}\), where we define \[\mathcal{H}_{\mathcal{I}} = \bigcap_{i\in \mathcal{I}}\mathcal{H}_i.\] Thus, \(\mathcal{H}_{\mathcal{I}}\) is true if and only if all hypotheses \(\mathcal{H}_i\) with \(i\in \mathcal{I}\) are true. The set of all intersection hypotheses is sometimes called the closure of the set of hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\). This is because this set is closed under intersections in the sense that if \(\mathcal{H}_{\mathcal{I}}\) and \(\mathcal{H}_{\mathcal{J}}\) (with \(\mathcal{I}\), \(\mathcal{J}\in \mathcal{C}\)) are intersection hypotheses, then their intersection is also an intersection hypothesis within this family. The hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) that we “start from” are called the elementary hypotheses. (Sometimes a condition is placed on the elementary hypotheses, e.g. that no elementary hypothesis is implied by another one.) Note that if \(\mathcal{I},\mathcal{J}\in \mathcal{C}\), then perhaps confusingly, \[\mathcal{I}\subseteq\mathcal{J}\Rightarrow \mathcal{H}_{\mathcal{I}} \supseteq \mathcal{H}_{\mathcal{J}}.\] For example, if \(\mathcal{I}=\{1,2\}\) and \(\mathcal{J}=\{1,2,3\}\), then \(\mathcal{H}_{\mathcal{J}}\) corresponds to the darkest part of the Venn diagram in Figure 11, while \(\mathcal{H}_{\mathcal{I}}\) corresponds to a larger area.

We have already seen some examples of intersection hypotheses. For instance, if we consider hypotheses \(\mathcal{H}_1:\theta\leq \theta_0\) and \(\mathcal{H}_2:\theta\geq \theta_0\), then the corresponding intersection hypothesis is \(\mathcal{H}_{\{1,2\}}:\theta= \theta_0\). As another example, we considered the global null hypothesis corresponding to a collection of hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\). This is the hypothesis that all \(m\) hypotheses are true, i.e., it is the hypothesis \(\mathcal{H}_{\{1,...,m\}}\).

A local test is simply a hypotheses test for an intersection hypothesis. We will assume throughout that each local test has level at most \(\alpha\). Of course the global hypothesis \(\mathcal{H}_{\{1,...,m\}}\) is also an intersection hypothesis, so we will call a global test a “local test” as well. Suppose we reject the global hypothesis. Then we know (with \((1-\alpha)100\%\) confidence) that at least one of the \(m\) hypotheses is false. Likewise, suppose we reject the intersection hypothesis \(\mathcal{H}_{\{2,5,6\}}\) for example. Then at least one of the hypotheses \(\mathcal{H}_2\), \(\mathcal{H}_5\), and \(\mathcal{H}_6\) is false (unless we made a false discovery).

In §3.2 we considered various global tests. Note that we can use these global tests as local tests. Indeed, a local test is a global test applied to a subcollection of hypotheses.

For example, suppose that for every hypothesis we have a valid p-value \(P_i\). Consider some \(\mathcal{I}\in \mathcal{C}\). By §3.2.1 a valid local test for \(\mathcal{H}_{\mathcal{I}}\) is then the test that rejects \(\mathcal{H}_{\mathcal{I}}\) if and only if \(\min\{P_i:i\in \mathcal{I}\}\leq \alpha/|\mathcal{I}|\). It will be useful to view a test as an indicator function \(\mathbb{1}(\cdot)\), i.e., a function that takes the values 0 and 1, where 1 indicates rejection and 0 indicates a non-rejection. Thus, we can write these global tests as \(\phi_{\mathcal{I}}\), where \(\phi_{\mathcal{I}}\) is the indicator function \[\label{ltBon} \phi_{\mathcal{I}}=\mathbb{1}\Big(\min\{P_i:i\in \mathcal{I}\}\leq \alpha/|\mathcal{I}|\Big).\tag{12}\]

As another example, suppose that for all nonempty \(\mathcal{I}\subseteq \mathcal{N}\), the \(P_i\) with \(i\in \mathcal{I}\) satisfy Simes’ inequality (see §3.2.4), i.e., \[\label{eq:Si} \mathbb{P}\Big\{\forall i \in \mathcal{I}: P_{(i)}^{\mathcal{I}}> \frac{i\alpha}{|\mathcal{I}|}\Big\}\geq1-\alpha,\tag{13}\] where we define \(P_{(1)}^{\mathcal{I}},..., P_{(|\mathcal{I}|)}^{\mathcal{I}}\) to be the sorted values of the \(P_i\) with \(i\in \mathcal{I}\). Then a valid local test for \(\mathcal{H}_{\mathcal{I}}\) is \[\label{lotestSimes} \phi_{\mathcal{I}}'=\mathbb{1}\Big(\exists i \in \mathcal{I}: P_{(i)}^{\mathcal{I}} \leq \frac{i\alpha}{|\mathcal{I}|} \Big).\tag{14}\] Clearly, the local test \(\phi_{\mathcal{I}}'\) is uniformly more powerful than \(\phi_{\mathcal{I}}\) (but requires the extra assumption that 13 holds if \(\mathcal{I}\subseteq \{1,...,m\}\)). Note that if we consider an elementary hypothesis \(\mathcal{I}=\{i\}\), then \(\phi_{\mathcal{I}}\) and \(\phi_{\mathcal{I}}'\) are both simply \(\mathbb{1}(P_i\leq \alpha)\).

We are now ready to define the closed testing procedure based on the local tests. It turns out that the elementary hypotheses rejected by the closed testing procedure based on the local tests 12 , are the same as those rejected by the Holm method. The closed testing procedure based on the local tests \(\phi_{\mathcal{I}}'\) potentially rejects more.

4.2 The general closed testing method for FWER control↩︎

The closed testing procedure tests all intersection hypotheses \(\mathcal{H}_{\mathcal{I}}\), \(\mathcal{I}\in \mathcal{C}\). Thus, to every intersection hypothesis, the procedure assigns a 1 (“reject”) or a 0 (“do not reject”). This means that the procedure in principle considers \(2^m-1\) tests. However, the procedure is often also powerful for simply testing \(\mathcal{H}_1,...,\mathcal{H}_m\), in case we are not interested in the intersection hypotheses.

As explained, the ingredients of the closed testing procedure are the local tests. Thus, suppose that for every intersection hypothesis \(\mathcal{H}_{\mathcal{I}}\) we have defined some local test \(\phi_{\mathcal{I}}\). Then the closed testing procedure is defined as follows.

Thus, the CTP rejects \(\mathcal{H}_{\mathcal{I}}\) if and only if \(\phi_{\mathcal{J}}=1\) for al \(I\subseteq \mathcal{J}\in \mathcal{C}\) — i.e., if all local tests \(\phi_{\mathcal{J}}\) with \(I\subseteq \mathcal{J}\in \mathcal{C}\) reject their corresponding intersection hypotheses. The \(\psi_{\mathcal{I}}\) are sometimes called the effective local tests, to distinguish them from the local tests. The CTP is the collection of all effective local tests. Figure 12 and its caption provide further explanation. Closed testing ensures that the set of rejected intersection hypotheses has a logical structure: If \(\mathcal{I}\subset \mathcal{J}\in\mathcal{C}\) and \(\mathcal{I}\) is rejected by the CTP, then \(\mathcal{J}\) is also rejected. This makes sense since if \(\mathcal{H}_{\mathcal{I}}\) is false, then one of the \(\mathcal{H}_i\) with \(i\in \mathcal{I}\) is false, which implies that \(\mathcal{H}_{\mathcal{J}}\) is false.

There are \(2^m-1\) local tests, and evaluating them all is computationally undoable if \(m\) is large. However, note that if an intersection hypothesis is not rejected, we know that all hypotheses implied by it (i.e. all hypotheses that can be reached by following the arrows in the closed testing tree) are also not rejected by the CTP. This can speed the computations up, but such a procedure can still be computationally feasible. Fortunately, in some cases, exact computational shortcuts are possible. There will be discussed later.

The reason why closed testing is called this way, is that the family of hypotheses \(\mathcal{H}_{\mathcal{I}}\) with \(\mathcal{I}\subseteq\{1,...,m\}\) is closed under intersections. The CTP tests the hypotheses in the closure of the original hypothesis family \(\mathcal{H}_1,...,\mathcal{H}_m\). Even if we are only interested in \(\mathcal{H}_1,...,\mathcal{H}_m\) and not in the intersection hypotheses, CTPs are often powerful. We now check that the CTP indeed controls the FWER. In case we are interested in all intersection hypotheses, the FWER is \(\mathbb{P}(V^{\text{tot}}>0)\), where \(V^{\text{tot}}\) is the number of false positives among all rejected intersection hypotheses: \[V^{\text{tot}}= |\{\mathcal{I}\subseteq\mathcal{N}: \mathcal{I}\neq\emptyset\text{ and } \psi_{\mathcal{I}}=1\}|.\] In case we are interested in the number of false positives among the rejected elementary hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\), we consider \[V^{\text{elem}}= |i\in \mathcal{N}: \psi_{\{i\}}=1\}|.\]

Clearly \(V^{\text{elem}}\leq V^{\text{tot}}\). Even if the CTP does not reject any elementary hypotheses, it might still reject some other intersection hypotheses. However, it turns out that in several important cases, no intersection hypotheses can be rejected by the CTP in case no elementary hypotheses are rejected by the CTP. A CTP is called consonant if it is always the case that for every \(\mathcal{H}_{\mathcal{I}}\) that is rejected by the CTP, there is an \(i\in \mathcal{I}\) such that \(\mathcal{H}_{i}\) is rejected by the CTP. For such procedures, all interesting information is contained in the rejections (by the CTP) of the elementary hypotheses. Indeed, if we know that a CTP is consonant and we know which elementary hypotheses it rejects, then the intersection hypotheses that it rejects are precisely those that are implied by the elementary hypotheses that are rejected by the CTP. For consonant procedures \(V^{\text{elem}}=0\Leftrightarrow V^{\text{tot}}=0\). If a CTP is nonconsonant, it may reject some intersection hypotheses without rejecting any elementary hypotheses. This idea will be important in §7.

An example of a consonant CTP is the CTP with local tests based on Bonferroni (this will be proved in §7). As already mentioned, the elementary hypotheses rejected by this CTP are the hypotheses that Holm rejects. An example of a nonconsonant (or “dissonant”) CTP is the CTP with local tests based on Simes’ inequality (see exercises). We will return to the topic of nonconsonance in §7.

It follows from the proof of Theorem 9 that for the CTP to ensure that \(\mathbb{P}(V^{\text{tot}}>0)\leq \alpha\), it suffices that the local test \(\phi_{\mathcal{N}}\) has type I error rate at most \(\alpha\). Thus, in this sense, we do not require all local tests to be valid, but only \(\phi_{\mathcal{N}}\). However, we do not know \(\mathcal{N}\) of course; often it could be any subset of \(\{1,...,m\}\), hence in a sense we need every local test to be valid.

4.3 Example of a CTP: Holm↩︎

Recall that a local test is simply a global test applied to an intersection hypothesis \(\mathcal{H}_{\mathcal{I}}\), \(\mathcal{I}\in \mathcal{C}\). Let us first consider the test from §3.2.1 by Bonferroni. This is the test \[\label{localtestHolm} \phi_{\mathcal{I}} = \mathbb{1}\big(\min\{P_i: i\in \mathcal{I}\}\leq \alpha /|\mathcal{I}|\big).\tag{15}\] If we consider the corresponding CTP, then it turns our that it rejects exactly the same elementary hypotheses as Holm’s method. Thus, to check which elementary hypotheses are rejected by the CTP, we do not need to perform \(2^m-1\) local tests, but can simply use Holm’s method, which is very fast. Thus, Holm’s method is a computational shortcut for the full CTP. (As has been mentioned and will be proved in §7, this CTP is consonant, so in a sense only the rejected elementary hypotheses are interesting. Holm tells us which these are.)

4.4 Example of a CTP: using Simes-type local tests↩︎

We will now consider Simes local tests, defined in 14 already: \[\phi_{\mathcal{I}} = \mathbb{1}\Big( \bigcup_{i\in \mathcal{I}} \{P^{\mathcal{I}}_{(i)}\leq \frac{i\alpha}{|\mathcal{I}|} \} \Big),\] where \(P^{\mathcal{I}}_{(1)}\leq ...\leq P^{\mathcal{I}}_{(|\mathcal{I}|)}\) are the sorted p-values \(P_i\) with \(i\in \mathcal{I}\).

As explained, the CTP is valid if the local test \(\phi_{\mathcal{N}}\) rejects with probability at most \(\alpha\). Thus, we must assume that inequality 7 holds. As already discussed, this is a very mild assumption on the dependence structure of the null p-values.

[33] first discussed closed testing based on Simes local tests. When statisticians say “Hommel’s method”, note that they usually do not mean the procedure that tests all intersection hypotheses, but just the method for testing the \(m\) elementary hypotheses. In general, “closed testing” does not necessarily mean that we are interested in all \(2^m-1\) intersection hypotheses: often we just use a closed testing method for testing the elementary hypotheses. (However, it turns out that the CTP based on Simes local tests is nonconsonant, so it can be useful to not merely look at the elementary rejections.)

The CTP based on Simes-type local tests is much more complex than Holm’s method in the sense that it is less straightforward to come up with a computationally efficient algorithm (i.e. a shortcut) for evaluating which of the elementary hypotheses are rejected. If we naively performed the CTP, this would become computationally undoable for moderate values of \(m\). Hommel’s original algorithm has a complexity that is quadratic in \(m\), which is already much better than a completely naive implementation. [34] were the first to come up with an algorithm with a complexity that is linear in \(m\) (after sorting the p-values, which takes \(O(m\log(m))\) time).

[33] showed that \(\mathcal{H}_i\) is rejected by the CTP if and only if \[h(\alpha)p_i\leq \alpha,\] where \[h(\alpha) = \max\big\{i\in\{0,...,m\}: i p_{(m-i+j)}> j\alpha \text{ for all } j\in\{1,...,i\}\big\} .\] [34] compute \(h(\alpha)\) in linear time (after sorting the p-values).

The algorithm of [34] is implemented in the R package hommel on CRAN [35]. Its function of the same name, hommel(), can be used as follows: hommel(p, simes = TRUE), where p contains the p-values. simes = TRUE gives the method described above and simes = FALSE gives a more conservative version of the method that does not require any assumptions on the dependence structure. The output of the function are the Hommel-adjusted p-values. Thus, all hypotheses with adjusted p-values at most \(\alpha\) can be rejected.

4.5 Example: one-way ANOVA model↩︎

In §3.2.6 we considered a global test for the one-way balanced ANOVA model. In the section we explore closed testing for this model. Some of these tests also have versions for non-balanced ANOVA models.

4.5.1 Closed testing method↩︎

The test from §3.2.6 rejects when the statistic 8 , \[\sqrt{2}\max_{1\leq i < i'\leq k} |T_{ii'}| =\frac{\max_{1\leq i \leq k} \bar{Y}_i- \min_{1\leq i \leq k}\bar{Y}_{i'}}{\hat{\sigma}/\sqrt{n}},\] is large. As discussed, under the global null hypothesis, this test statistic is a Studentized range variable with parameter \(k\) and \(\nu\) degrees of freedom, where \(\nu=k(n-1)\). We rejected the global null hypothesis when \[\sqrt{2}\max_{1\leq i < i'\leq k} |T_{ii'}| \geq Q_{k,\nu}^{(1-\alpha)}.\]

To see the problem through the lens of closed testing, we must first define elementary hypotheses. These are the hypotheses \(\mathcal{H}_{\mathcal{A}}\) where \(\mathcal{A}\subseteq\{1,...,k\}\) and \(|\mathcal{A}|=2\). Here we interpret the notation \(\mathcal{H}_{\mathcal{A}}\) not as earlier in this section, but as in §3.7: it is the hypothesis that the \(\theta_i\) with \(i\in \mathcal{A}\) are equal to each other. Clearly, the number of elementary hypotheses is \(\binom{k}{2}\). We will call the number of elementary hypotheses \(m\), like before.

We now consider the intersection hypotheses, i.e., all possible intersection of the \(m\) elementary hypotheses. At first sight, there seem to be \(2^m-1\) intersection hypotheses. However, we are then counting some intersection hypotheses multiple times. For example, consider the elementary hypotheses \(\mathcal{H}_{\{1,2\}}\), \(\mathcal{H}_{\{2,3\}}\) and \(\mathcal{H}_{\{1,3\}}\). The intersections \(\mathcal{H}_{\{1,2\}}\cap \mathcal{H}_{\{2,3\}}\), \(\mathcal{H}_{\{2,3\}}\cap\mathcal{H}_{\{1,3\}}\) and \(\mathcal{H}_{\{1,2\}}\cap\mathcal{H}_{\{1,3\}}\) are in fact the same: they are all equal to \(\mathcal{H}_{\{1,2,3\}}\). Further, note that all intersection hypotheses can be written as intersections of the form \[\mathcal{H}_{\mathcal{A}_1}\cap\mathcal{H}_{\mathcal{A}_2}\cap...\] with \(\mathcal{A}_1, \mathcal{A}_2,...\) disjoint.

4.5.2 Shortcuts↩︎

As always, naive closed testing is computationally infeasible when \(m\) is large, so it is useful to look for shortcuts.

• Tukey-Welsch. One important shortcut is the Tukey-Welsch procedure from §3.7. The advantage of Tukey-Welsch is that most of the intersection hypotheses do not need to be individually considered, since we only need to consider some of the hypotheses \(\mathcal{H}_{\mathcal{A}}\). Tukey-Welsch is considered to have good power, although its power is not equal to that of the full CTP (the full CTP is a uniform improvement of Tukey-Welsch).

• Newman-Keuls and Begun-Gabriel. In case studentized range test statistics are used, a small modification of the Tukey-Welsch procedure gives a method that is even faster and is actually uniformly more powerful than the CTP. This method is referred to as the Newman-Keuls procedure [25, p. 66]. This method is not known to control the FWER. However, since it is uniformly more powerful than the CTP, and since the CTP is uniformly more powerful than Tukey-Welsch, we can use Newman-Keuls as an upper bound for the CTP’s decision and Tukey-Welsch as a lower bound for the CTP’s decision. This can be helpful since we know that all hypotheses that are rejected by Tukey-Welsch will also be rejected by the CTP and all hypotheses that are not-rejected by Newman-Keuls will not be rejected by the CTP. Hence, after running Tukey-Welsch and Newman-Keuls, there is only a limited set of hypotheses left that is contentious, i.e., for which it should be further checked whether they can be rejected or not. (This is exploited by a method by Begun and Gabriel, see [25, p. 122].)

  The Newman-Keuls method is defined as folows. We let \(\alpha_p:=\alpha\) for every \(2\leq p \leq k\). We then re-name the parameters such that \(\hat{\theta}_{1}\leq ... \leq \hat{\theta}_{k}\). Then we consider Tukey-Welsch, but with this simpler definition of \(\alpha_p\). Note that if this method rejects \(\mathcal{H}_{\mathcal{A}}\) with \(\min(\mathcal{A})=i\) and \(\max(\mathcal{A})=j\), then it follows that all hypotheses implied by \(\mathcal{H}_{\mathcal{A}}\) and containing \(i\) and \(j\) are rejected, if studentized range statistics are used. Thus, when we reject \(\mathcal{H}_{\mathcal{A}}\), we can automatically reject all subsets of size \(|\mathcal{A}|-1\) as well, except \(|\mathcal{A}|\setminus \{i\}\) and \(|\mathcal{A}|\setminus \{j\}\). Thus, if we reject \(\mathcal{H}_{\mathcal{A}}\), then we can proceed by checking \(\mathcal{H}_{|\mathcal{A}|\setminus \{i\}}\) and \(\mathcal{H}_{|\mathcal{A}|\setminus \{j\}}\), etcetera. In this way, Newman-Keuls saves time compared to Tukey-Welsch.

• Adaptation of Tukey-Welsch. Finally, we can consider the following procedure, which is somewhat more conservative than Tukey-Welsch but faster [25, p. 115]. This adaptation can be used when studentized range test statistics are used. First note that Newman-Keuls exploits the fact that \(Q_{p,\nu}^{(1-\alpha_p)}\) decreases as \(p\) decreases, when \(\alpha_p:=\alpha\) for every \(2\leq p \leq k\). This is not guaranteed is \(\alpha_p\) is defined as in Tukey-Welsch, i.e., we may then sometimes have \(Q_{p,\nu}^{(1-\alpha_{p-1})} > Q_{p,\nu}^{(1-\alpha_p)}\). To resolve this, we can enforce montonicity by making each \(Q_{p,\nu}^{(1-\alpha_{p+1})}\) at least as large as the previous critical value \(Q_{p,\nu}^{(1-\alpha_{p})}\). We then obtain a method that is somewhat more conservative than Tukey-Welsch, but for which we can use an algorithm analogous to that of Newman-Keuls.

• Other methods for ANOVA models. There are also other methods for testing in ANOVA models. A well-known method is Fisher’s LSD (“Least Significant Difference”) procedure, but it only provides weak FWER control (as defined in §3.3). It proceeds by first running a global F-test. If the global null is rejected, the method then tests all pairwise hypotheses without any further multiple testing correction.

  Faster but less powerful alternatives for Tukey-Welsch are Tukey’s HSD (“Honestly Significant Difference”) method and Tukey-Kramer [28]. These do provide strong FWER control.

4.6 Exercises↩︎

1. Remove and add some crosses in the “closed testing tree” from Figure 12 (the crosses represent the hypotheses rejected by the their local tests) and decide which hypotheses are rejected by the closed testing procedure. Additionally, draw the closed testing tree in case there are 4 elementary hypotheses.

2. A multiple testing procedure — which decides for each of the hypotheses in some family whether it is rejected or not — is called coherent if for each pair of hypotheses \(\mathcal{H}^a\) and \(\mathcal{H}^b\) in the family the following always holds: if \(\mathcal{H}^a\subseteq \mathcal{H}^b\)⁴ and \(\mathcal{H}^b\) is rejected, then \(\mathcal{H}^a\) is rejected.

   Consider hypotheses \(\mathcal{H}_1,...,\mathcal{H}_k\) (these need not be elementary hypotheses) and consider some incoherent MTP for these \(k\) hypotheses (i.e. a procedure that “accepts” or rejects each hypothesis). Assume the procedure controls the FWER at level \(\alpha\). Now consider a new MTP, which rejects each hypothesis \(\mathcal{H}_i\) (\(1\leq i \leq k\)) that implies a hypothesis that is rejected by the original procedure. Thus, the new procedure rejects the same hypotheses as the old one, and possibly more. Prove that the new procedure is coherent and that it also controls the FWER at level \(\alpha\) (the proof is short).

3. Consider some family \(\mathcal{H}_1,...,\mathcal{H}_k\) of hypotheses (not necessarily elementary hypotheses). Recall that a multiple testing procedure for this family is called coherent if for each pair of hypotheses \(\mathcal{H}_i\) and \(\mathcal{H}_i\) in the family the following holds: if \(\mathcal{H}_i\subseteq \mathcal{H}_j\) and \(\mathcal{H}_j\) is rejected, then \(\mathcal{H}_i\) is rejected.

   Consider a coherent MTP for the family \(\mathcal{H}_1,...,\mathcal{H}_k\) and suppose the procedure has FWER at most \(\alpha\). We will now prove that this MTP is the result of applying closed testing based on a collection of local tests with size at most \(\alpha\). Thus, every coherent procedure can be seen as a closed testing procedure.

   Consider the closure (under intersections) of the family \(\mathcal{H}_1,...,\mathcal{H}_k\) and define a local test for each hypothesis \(\mathcal{H}\) in this closure as follows. If \(\mathcal{H}\) is in the original family, reject \(\mathcal{H}\) locally if and only the original MTP rejects it. Otherwise, reject \(\mathcal{H}\) locally if and only if at least one of the locally rejected hypotheses (by the new tests) among \(\mathcal{H}_1,...,\mathcal{H}_k\) is implied by \(\mathcal{H}\).

   (a) Now consider the CTP based on these local tests. Prove that this CTP has FWER (over all intersection hypotheses) at most \(\alpha\).

   (b) Prove that among \(\mathcal{H}_1,...,\mathcal{H}_k\), the hypotheses rejected by this CTP are the same as those rejected by the original MTP.

4. Let \(\alpha\in(0,1)\). Consider \(m\geq 3\) elementary hypotheses with corresponding p-values \(P_1,...,P_m\). Use a counterexample to show that the corresponding CTP based on Simes local tests is not necessarily consonant. (Hint: Consider e.g. \(p_1=...=p_{m-1}=\alpha\frac{m-1}{m}\) and \(p_m>\alpha\).)

5. Prove that if \(m=2\), the CTP based on Simes local tests is consonant.

6. \(\bigstar\) Explain that in the specific context of §4.5, the number of distinct intersections of elementary hypotheses equals the number of partitions of the set \(\{1,...,k\}\), minus 1. (A partition is a collection of subsets of \(\{1,...,k\}\), such that each element of \(\{1,...,k\}\) is contained in exactly one of these subsets.)

4.7 Solutions↩︎

1. In case there are 4 hypotheses, the closed testing tree consists of 4 layers. The top layer consists of 1 hypothesis, \(\mathcal{H}_1\cap....\cap\mathcal{H}_4\). The second layer consists of the four hypotheses \(\mathcal{H}_1\cap \mathcal{H}_2\cap\mathcal{H}_3\), \(\mathcal{H}_1\cap \mathcal{H}_2\cap\mathcal{H}_4\), \(\mathcal{H}_1\cap \mathcal{H}_3\cap\mathcal{H}_4\), \(\mathcal{H}_2\cap \mathcal{H}_3\cap\mathcal{H}_4\). The third layer consists of \(\binom{4}{2}=6\) hypotheses and the fourth layer consists of the four elementary hypotheses \(\mathcal{H}_1,....,\mathcal{H}_4\). From each node, you should draw an arrow to the nodes directly below it corresponding to the implied hypotheses. For example, from \(\mathcal{H}_1\cap \mathcal{H}_2\cap\mathcal{H}_3\) there should be arrows to \(\mathcal{H}_1\cap \mathcal{H}_2\), \(\mathcal{H}_1\cap\mathcal{H}_3\) and \(\mathcal{H}_2\cap\mathcal{H}_3\).

2. Consider one of the hypotheses. By construction, all hypotheses that imply it are rejected by the new procedure. Hence, the new procedure is coherent. Let \(\mathcal{E}\) be the event that the old procedure makes no false discoveries. Thus, all hypotheses rejected by the original procedure are false. The new procedure rejects those hypotheses, plus hypotheses which imply them, which must als be false (If \(\mathcal{H}_i\) implies \(\mathcal{H}_j\) and \(\mathcal{H}_j\) is false, then \(\mathcal{H}_i\) must be false.)

3. (a)

   With probability at least \(1-\alpha\), the original MTP makes no false discoveries. Suppose it makes no false discoveries. Thus, for the new procedure, all local rejections of hypotheses among \(\mathcal{H}_1,...,\mathcal{H}_k\) are correct. The new procedure then also rejects every hypotheses \(\mathcal{H}\) that implies one of the rejected hypotheses among \(\mathcal{H}_1,...,\mathcal{H}_k\). Clearly, if \(\mathcal{H}_i\) is correctly rejected, then it is also correct to reject all hypotheses that imply \(\mathcal{H}_i\). Thus, with probability at least \(1-\alpha\), none of the local tests of the new procedure make a false discovery. Consequently, the new CTP has FWER at most \(\alpha\).

   (b) Consider \(1\leq j \leq k\) and suppose \(\mathcal{H}_j\) is rejected by the original procedure. We will show that it is rejected by the defined closed testing procedure. Note that the local test rejects \(\mathcal{H}_j\). Moreover, by coherence, all hypotheses among \(\mathcal{H}_1,...,\mathcal{H}_k\) that imply \(\mathcal{H}_j\), are rejected by the original procedure. Consider an intersection hypothesis \(\mathcal{H}_{\mathcal{I}}\) with \(j\in \mathcal{I}\). \(\mathcal{H}_j\) is rejected locally and is implied by \(\mathcal{H}_{\mathcal{I}}\), so by definition of the local tests, \(\mathcal{H}_{\mathcal{I}}\) is rejected locally. Thus, all \(\mathcal{H}_{\mathcal{I}}\) with \(j\in \mathcal{I}\) are rejected locally, so \(\mathcal{H}_j\) is rejected by the defined CTP.

   Conversely, suppose \(\mathcal{H}_j\) is rejected by the defined closed testing procedure. Then \(\mathcal{H}_j\) is rejected locally, and by construction this means the original procedure rejects it. This finishes the proof.

4. Suppose \(p_1=...=p_{m-1}=\alpha\frac{m-1}{m}\) and \(p_m>\alpha\). Note that \(\mathcal{H}_{\{1,...,m\}}\) is rejected by its local test, since \(p_{m-1}\leq\alpha\frac{m-1}{m}\). Since \(\mathcal{H}_{\{1,...,m\}}\) is at the top of the closed testing tree, it directly follows that it is rejected by the CTP.

   Now consider any \(\mathcal{I}\subseteq \{1,...,m\}\) with \(|\mathcal{I}|=m-1\) and \(m\in \mathcal{I}\). Let \(p_{(1)}^{\mathcal{I}}\leq...\leq p_{(m-1)}^{\mathcal{I}}\) be the sorted p-values corresponding to \(\mathcal{I}\). \(\mathcal{H}_{\mathcal{I}}\) is rejected by its local test if and only if there is an \(1\leq i \leq m-1\) such that \(p_{(i)}^{\mathcal{I}}\leq c^\mathcal{I}_i\), where \(c^\mathcal{I}_i =\frac{\alpha i}{|\mathcal{I}|}\). Note that \(\mathcal{H}_{\mathcal{I}}\) is not rejected by its local test, since \(p^{\mathcal{I}}_{|\mathcal{I}|-1}= \alpha\frac{ m-1}{m}\), which exceeds the critical value \(c^{\mathcal{I}}_{|\mathcal{I}|-1}=\alpha\frac{ |\mathcal{I}|-1}{|\mathcal{I}|}=\alpha\frac{ m-2}{m-1}\), and likewise the other p-values corresponding to \(\mathcal{I}\) exceed their critical values.

   We have seen that none of the hypotheses \(\mathcal{H}_{\mathcal{I}}\) with \(\mathcal{I}\subseteq \{1,...,m\}\) and \(|\mathcal{I}|=m-1\) and \(m\in \mathcal{I}\) are rejected by the their local tests. For every \(1\leq i \leq m\) there is such an \(\mathcal{I}\) which contains \(i\). It follows that none of the elementary hypotheses are rejected by the CTP. However, \(\mathcal{H}_{\{1,...,m\}}\) is rejected by the CTP. Hence, the procedure is nonconsonant.

5. Apart from the elementary hypotheses \(\mathcal{H}_1\) and \(\mathcal{H}_2\), say, there is only one intersection hypothesis, \(\mathcal{H}_{12}\). Let \(p_1, p_2\) be the p-values corresponding to \(\mathcal{H}_1\) and \(\mathcal{H}_2\). Suppose \(\mathcal{H}_{12}\) is rejected by the CTP. That means that either \(p_{(1)}\leq \alpha\frac{1}{2}\) or \(p_{(2)}\leq \alpha\) (or both). Thus, at least one of the p-values \(p_1, p_2\) is at most \(\alpha\). For \(i\in\{1,2\}\), the local test for \(\mathcal{H}_i\) simply rejects when \(p_i\leq \alpha\). Hence, at least one of the elementary hypotheses is rejected locally, and hence also by the CTP. Thus, if \(\mathcal{H}_{12}\) is rejected by the CTP, then at least one of the elementary hypotheses is rejected. We conclude that the method is consonant when \(m=2\).

6. Let \(\mathcal{P}\) be the set of all possible partitions of \(\{1,...,k\}\) exluding the partition that consists of \(m\) singletons. It suffices to show that there is a one-to-one mapping (bijection) from \(\mathcal{P}\) to the set of all intersection hypotheses, as follows. We define this mapping as follows: we map every partition in \(\mathcal{P}\) to the intersection hypothesis \(\mathcal{H}_{\mathcal{A}_1}\cap \mathcal{H}_{\mathcal{A}_2}\cap ....\), where \(\mathcal{A}_1,\mathcal{A}_2,...\) are the elements of the partition that contain at least 2 elements. To see that this is a one-to-one mapping, consider a specific intersection hypothesis and note that it can be uniquely written as \(\mathcal{H}_{\mathcal{A}_1}\cap \mathcal{H}_{\mathcal{A}_2}\cap ....\), where \(\mathcal{A}_1,\mathcal{A}_2,...\) are disjoint and have size at least 2. Clearly, there is only one element in \(\mathcal{P}\) that maps to this intersection hypothesis — namely, the partition that contains \(\mathcal{A}_1,\mathcal{A}_2,...\) and that contains all elements of \(\{1,...,k\}\) that are not in \(\mathcal{A}_1\cup\mathcal{A}_2,...\) as singletons.

5 Resampling-based FWER control↩︎

In §2.4.4 and §2.4.6 in particular, we discussed resampling-based nonparametric and semiparametric testing (e.g. permutation and bootstrap testing) in case there is a single hypothesis. There, we discussed an advantage of such tests: they often require fewer assumptions than fully parametric tests. Among multiple testing methods, there are also some important resampling-based methods. Besides the mentioned robustness, they have another advantage: they take into account the dependence structure of the test statistics, without requiring strong assumptions on this dependence structure. This allows these methods to often have much better power than methods not based on resampling. An example of a method that does not account for the dependence structure, is Holm. That method is very conservative when there are strong positive dependencies between the p-values. The resampling-based maxT method discussed below, is then potentially a much more powerful alternative.

5.1 Definition of the maxT method↩︎

For FWER control, one resampling-based method is by far the most popular: this is the maxT method by [36]. A related method is the minP method, which we do not cover here [30], [37]. MaxT can be based on permutation testing or bootstrap testing. In case of bootstrap testing, the method is sometimes referred to as the Romano-Wolf method [38]. Here we will refer to the approach as maxT. In a general sense, the maxT method need not be based on resampling, and e.g. Holm could then be seen as a special case of maxT. However, maxT tends to be most powerful when combined with resampling, and here we will focus on resampling-based versions of maxT.

Both in §2.4.4 and §2.4.6, we considered “resampled test statistics”:

• In §2.4.4 we considered the statistics \((T(gX):g\in\mathcal{G})\), where \(\mathcal{G}\) is a group of transformations (e.g. permutation maps), or we use a collection of random transformations from \(\mathcal{G}\).

• In §2.4.6 we considered the test statistics \((T(X^j): 1\leq j \leq b)\), where \(X^1,...,X^b\) are bootstrap samples.

Let \(w\) denote the number of transformations used (in case we use group invariance testing, e.g. permutation testing) or the number of bootsrap samples (in case bootstrap testing is used). The maxT method creates resampled test statistics as in §2.4.4 and §2.4.6, except that it does this for every hypothesis, and it applies simultaneous permutation or simultaneous bootstrapping (further explained below). Since the method resamples test statistics for every hypothesis, it does not merely end up with a vector of length \(w\), but with an \(m\)-by-\(w\) matrix, for which we use the following notation: \[\label{eq:matrixT} \boldsymbol{T} = (T_i^j)_{1\leq i \leq m, 1\leq j \leq w}= \begin{bmatrix} T_1^1 & T_1^2 & \ldots & T_1^w \\ T_2^1 & T_2^2 & \ldots & T_2^w \\ \vdots & \vdots & \ddots & \vdots \\ T_m^1 & T_m^2 & \ldots & T_m^w \end{bmatrix} ,\tag{16}\] where we will assume that the first column corresponds to the original data. Thus, if we use random permutations, we use the identity map plus \(w-1\) random permutations. When constructing this matrix, the key point is the simultaneous resampling. In case permutations are used this means that for every \(1\leq j \leq w\), \(T_1^j,...,T_m^j\) are obtained using the same permutation map. In case bootstrapping is used, it means that for every \(1\leq j \leq w\), \(T_1^j,...,T_m^j\) are obtained using the same bootstrap sample of the data. This means that in some particular sense, these transformations preserve the dependence structure of the test statistics. More precisely, let the \(N\)-by \(w\) matrix \(\boldsymbol{T}_{\mathcal{N}}\) contain the rows of \(\boldsymbol{T}\) with indices in \(\mathcal{N}\). A key idea behind the method is that all rows of \(\boldsymbol{T}_{\mathcal{N}}\) have the same joint distribution (for finite \(n\) or asymptoticaly).

There is a single-step maxT procedure and a more powerful sequential, “step-down” version. The single-step maxT method is simply the first step of the sequential procedure — just like Bonferroni is the first step of Holm. The (sequential maxT procedure is in fact equivalent to a CTP, as discussed later. The sequential maxT procedure takes the matrix \(\boldsymbol{T}\) as input and is defined as follows:

1. For each \(1\leq j \leq w\), compute the maximum \(M^j:=\max\{T_i^j: 1\leq i \leq m\}\) and let \(q\) be the \((1-\alpha)\)-quantile of these maxima, i.e., \[q= \min\{x\in \mathbb{R}: \frac{|\{1\leq j \leq w: M^j\leq x \}|}{w}\geq 1-\alpha \}.\] Reject all hypotheses \(\mathcal{H}_i\) with test statistics \(T_i>q\). (Here \(T_i\) means \(T_i^1\).)

2. If nothing was rejected in the previous step, stop. Otherwise, delete all rows from \(\boldsymbol{T}\) corresponding to the rejected hypotheses and repeat the above. This leads to a potentially smaller threshold, \(q_2\) say, which may lead to additional rejections.

3. Continue like this until a step is reached where no additional hypotheses are rejected.

MaxT tends to be more powerful than Bonferroni or Holm when there are positive dependencies between the tested variables. The reason is that the stronger the positive dependence is between the test statistics in a column of \(\boldsymbol{T}\), the smaller the maximum of these test statistics tends to be, and hence the smaller the rejection threshold (\(q\)) tends to be.

As we will discuss further in §5.4, the sequential maxT method is related to Holm’s method. In the same way as Holm is the sequential version of Bonferroni, maxT is the sequential version of single-step maxT. Essentially the only difference between Bonferroni and single-step maxT, is that Bonferroni uses a “worst-case” rejection threshold \(\alpha/m\) that uses no knowledge on the dependence structure, while maxT uses a potentially less restrictive threshold. Another difference is that maxT uses test statistics rather than p-values, but that is not a fundamental difference; for example, we could use \(-P_1,...,-P_m\) as test statistics within maxT.

5.2 Example data analysis↩︎

To illustrate the maxT method, we revisit the dataset from §2.3.1 about car models. An analysis of a larger dataset with maxT is in §6.5. For illustration purposes, we consider the same limited set of car models as in §2.3.1 (\(n=5\)), but now we consider all the variables, see Figure 13 (a). We define the mpg variable to be \(X\) and we define the other 10 variables to be \(Y^1,...,Y^{10}\). Suppose we want to test the hypotheses \(\mathcal{H}_i: X \mathrel{\perp\!\!\!\perp}Y^i\), \(1\leq i \leq m=10\). We discussed in §2.4.2 how we can perform a permutation test to test such a hypothesis. Now we proceed in a similar way, except that we simultaneously permute the columns corresponding to the \(Y^i\). There are \(5!=120\) ways in which this can be done. Permuting the columns corresponding to the \(Y^i\) will be equivalent to simply permuting the first column as done in Figure 13.

As test statistics we will use the absolute values of the correlations, i.e., \(T_i = |\rho(X,Y^i)|\), \(1\leq i \leq m\). For each of the 120 permuted versions of the data we compute the test statistics, and we collect the test statistics in the matrix \(\boldsymbol{T} = (T_i^j)_{i,j}\) as in 16 . For each column of \(\boldsymbol{T}\) we compute the maximum of the statistics. A histogram of the 120 resulting maxima is shown in Figure 14. We take \(\alpha=0.1\) for illustration. The \((1-\alpha)\)-quantile of the maximum-statistics \(M^j\) was \(q\approx 0.931\). Two of the test statistics exceeded \(q\). By coincidence, these were the first two test statistics: \(T_1\approx 0.953\) and \(T_2\approx0.941\). This finishes the first step of the maxT method. For the second step, we remove the columns corresponding to \(Y^1\) and \(Y^2\) and essentially apply the same procedure again. The \((1-\alpha)\)-quantile of the maximum-statistics is now \(0.929\). None of the remaining test statistics \(T_3,...,T_{10}\) exceed this value, so we get no further rejections and stop the procedure. Thus the sequential maxT was able to reject 2 hypotheses in this example.

We have seen that maxT rejects 2 hypotheses. Alternatively, we could compute individual p-values for each of the 10 hypotheses as in §2.3.1, and then apply Bonferroni (or Holm). The smallest (non-adjusted) p-value that we then obtain is \(2/120=1/60\). The smallest Bonferroni-adjusted p-values is thus \((1/60)\cdot m=1/6\), which is larger than \(\alpha\). Thus we see that if we perform classical permutation tests and use Bonferroni (or Holm), we do not reject anything.

If more permutations are available, so that p-values can be very tiny, maxT often still has much better power, when there are strong dependencies among the test statistics. This will be illustrated with a different dataset in §6.5. There Holm rejects 53 hypotheses, while maxT rejects 74 hypotheses. When there are many strong positive dependencies between the test statistics, the relative difference between Holm and maxT can be even much bigger.

5.3 Theory for the maxT method↩︎

Here we make more precise in which case the maxT method controls the FWER, for finite samples or asymptotically. Asymptotic FWER control means that \(\liminf_{n\rightarrow\infty}FWER_n\leq \alpha\), where the substript \(n\) indicates the sample size. We first briefly discuss using bootstrapping and then discuss using permutations (or other groups of transformations).

Bootstrap sampling. Recall that we defined \(\boldsymbol{T}_{\mathcal{N}}\) to be the matrix containing the rows of \(\boldsymbol{T}\) with indices in \(\mathcal{N}\). Recall that we assume in this section that the first row of \(\boldsymbol{T}\) is based on the original data. In case of bootstrapping, a sufficient condition for asymptotic FWER control is that asymptotically, the rows of \(\boldsymbol{T}_{\mathcal{N}}\) are i.i.d.. This will often be the case if we proceed analogously to §2.4.6, where now instead of bootstrap-sampling individual test statistics, we simultaneously bootstrap-sample whole columns of \(\boldsymbol{T}\) (i.e., all values in the column are based on the same bootstrap sample of the data). For finite samples, bootstrap methods are typically not exact.

Permutations. In case we use permutations, we may even have finite-sample FWER control. Consider the following assumption. Let \(X\) be the full dataset and let \(\mathcal{G}\) be the group of transformations (e.g. permutation maps) applied to the data.

This assumption is satisfied in particular if the joint distribution of the part of the data corresponding to \(\mathcal{N}\) is unchanged by the transformations. In the example from §5.2, a sufficient condition for the assumption to be satisfied, is that permuting the first column of the dataset does not change the joint distribution of the first column and the columns corresponding to the true hypotheses — or, equivalently, that permuting column 2 to 11 does not change that joint distribution.

In case random permutations are used (resampled with or without replacement, but such that the first transformation is fixed to be the identity) then maxT also controls the FWER (in a marginal sense). In practice, often random transformations are used, to limit the computational burden.

5.4 MaxT is a closed testing procedure↩︎

Note that maxT is similar to Holm, in the sense that it starts with a first step, and then if there are any rejections, it ignores the hypotheses that have been rejected and repeats the first step, etcetera. We saw that Holm is an exact shortcut for a CTP. Likewise, maxT is an exact shortcut for a (consonant) CTP, namely the CTP based on the following local tests for the intersection hypotheses \(\mathcal{H}_I\): \[\label{eqltmaxT} \phi_{\mathcal{I}}= \mathbb{1}\big( \max\{T_i:i\in\mathcal{I}\} > q^{\mathcal{I}} \big),\tag{17}\] where \(q^{\mathcal{I}}\) is the \((1-\alpha)\)-quantile of the maxima \(\max\{T_i(gX):i\in \mathcal{I}\}\), \(g\in \mathcal{G}\). To show this, we can reason in a way that is similar to how we showed that Holm is a CTP (Theorem 10).

5.5 Software implementations of maxT↩︎

One can easily program maxT oneself, although implementations may vary in their computational efficiencies. Some efficient implementations of the maxT method are contained in the R package multtest [39], [40]. This package also provides maxT-adjusted p-values (see [41] for a formula). Recall that the adjusted p-value for a hypothesis is the smallest \(\alpha\) for which the hypothesis is rejected by the multiple testing procedure.

As mentioned, in case bootstrapping is used, the maxT method is sometimes called the Romano-Wolf method. The method is implemented under that name in Stata [42].

5.6 Extension to k-FWER control↩︎

Just like the methods by Bonferroni and Holm, the maxT method can be generalized to k-FWER control. There is a single-step version and a sequential version. The single-step method is again simply the first step of the sequential version. Here we only provide the single-step version. For the sequential version, see Algorithm 2.1 of [43]. The sequential version (for \(k>1\)) tends to be much more computationally intensive than sequential maxT. This has to do with the fact that after performing e.g. the first step, we cannot be confident that the rejected hypotheses are all false (since we allow for \(k-1\) false positives).

For any collection of numbers, define the k-maximum as the \(k\)-th largest number, i.e., the smallest real number \(x\) such that at most \(k-1\) of the numbers exceed \(x\). The single-step k-FWER version of maxT is exactly the same as maxT, except that we now use the \(k\)-maxima where before we would take the maxima. Thus, \(q\) is now the \((1-\alpha)\)-quantile of the \(k\)-maxima \[k\text{-}\max\{T_i^j: 1\leq i \leq m\},\quad 1\leq j \leq w.\] Then, under Assumption 3, we obtain k-FWER control, and in case of bootstrapping we can obtain asymptotic k-FWER control.

5.7 Exercises↩︎

1. Show that maxT cannot reject any hypotheses if \(w<\alpha^{-1}\).

2. In practice one would take \(w\) large, but here we take \(w=5\) for simplicity. Suppose we are interested in testing \(m=3\) hypotheses and have computed the follow matrix \(\boldsymbol{T}\): \[\boldsymbol{T} = \begin{bmatrix} 6.2 & 6.0 & 4.1 & 5.7 & 4.4 \\ 6.3 & 4.0 & 6.2& 4.0 & 3.7 \\ 3.5 & 5.4 & 4.4 & 2.9 & 3.5\\ \end{bmatrix} ,\] As usual, the first column contains the test statistics corresponding to the original data. Take \(\alpha=0.2\). Which hypotheses do the single-step and the sequential version of maxT reject in this case? How many hypotheses are rejected when \(\alpha<0.2\)?

3. \(\bigstar\) Consider data \(X\), hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) and corresponding test statistics \(T_1(\cdot),...,T_m(\cdot)\). Suppose Assumption 3 holds. Consider some rejection threshold \(t\in \mathbb{R}\). Thus, we reject the hypotheses with indices in \(\mathcal{R}(X):=\{1\leq i \leq m: T_i(X)>t\}\).

   For every \(g\in \mathcal{G}\), let \[R^g(X):=|\{1\leq i \leq m: T_i(gX)>t\}|.\] Let \[R^{(1)}\leq ...\leq R^{(|\mathcal{G}|)}\] be the sorted values \(R^g\), \(g\in \mathcal{G}\).

   Let \(\alpha\in(0,1)\) and let \(\bar{V}(X):= R^{(k)}\), where \(k=\lceil (1-\alpha)|\mathcal{G}|\rceil\).

   Prove that \[\mathbb{P}(V(X)\leq \bar{V}(X))\geq 1-\alpha,\] i.e., \(\bar{V}(X)\) is a \((1-\alpha)100\%\)-confidence upper bound for the number of false positives \(V(X)=|\mathcal{N}\cap\mathcal{R}(X)|\). (This approach is discussed and extended in [44].)

5.8 Solutions↩︎

1. For each \(1\leq j \leq w\), single-step maxT computes the maximum \(M^j:=\max\{T_i^j: 1\leq i \leq m\}\) and lets \(q\) be the \((1-\alpha)\)-quantile of these maxima, i.e., \(q= \min\{x\in \mathbb{R}: \frac{|\{1\leq j \leq w: M^j\leq x \}|}{w}\geq 1-\alpha \}.\) Here \(\frac{|\{1\leq j \leq w: M^j\leq x \}|}{w}\) takes values in \(\{1/w, 2/w, .... w/w\}\). If \(w<\alpha^{-1}\), then \((w-1)/w<1-\alpha\). It follows that \(q\) is actually the maximum of the values \(M^j\), \(1\leq j \leq w\). Thus, \(q\) is the maximum of all values in \(\boldsymbol{T}\). The single-step maxT method only rejects hypotheses \(\mathcal{H}_i\) with \(T_i>q\). It follows that single-step maxT rejects no hypotheses. Consequently, the sequential maxT method rejects no hypotheses either.

2. We start with the first step of the maxT procedure. For each permutation or bootstrap sample, we must compute the maximum of the \(m\) teststatistics. These maxima are 6.3, 6.0, 6.2, 5.7 and 4.4. We must then find the empirical \((1-\alpha)100\%=80\%\)-quantile of these maxima, which is 6.2. We can then reject all hypotheses with test statistics strictly larger than 6.2. There is only one such hypothesis, namely \(\mathcal{H}_2\), which has a statistic of 6.3. Thus, in the first step, we reject \(\mathcal{H}_2\).

   We now perform the second step. After removing the second row, the maxima are 6.2, 6.0, 4.4, 5.7 and 4.4. The \(80\%\)-quantile of these maxima is 6.2. We can then reject all hypotheses with test statistics strictly larger than 6.2. This leads to no additional rejections. Thus, maxT only rejects \(\mathcal{H}_2\).

   In case \(\alpha<0.2\), the \((1-\alpha)100\%\)-quantile of the maxima is simply the maximum of the maxima. No test statistic can exceed the maximum of the maxima, so no hypotheses are rejected when \(\alpha<0.2\).

3. Let \[V^{(1)}\leq ...\leq V^{(|\mathcal{G}|)}\] be the sorted values \(V(gX)\), \(g\in \mathcal{G}\). Here \(V(gX)= |\mathcal{N}\cap\mathcal{R}(gX)|\), where \[\mathcal{R}^g(X):=\{1\leq i \leq m: T_i(gX)>t\}.\]

   Note that if we see the values \(V(gX)\) as the “test statistics”, then Assumption 1 holds, i.e., the joint distribution of the values \(V(gX)\), \(g\in \mathcal{G}\), is invariant under all transformations in \(\mathcal{G}\) of \(X\). Hence, by Theorem 1, \[\mathbb{P}(V(X)\leq V^{(k)}(X))\geq 1-\alpha.\] Finally note that \(V^{(k)}(X)\leq \bar{V}(X)\), so that \[\mathbb{P}(V(X)\leq \bar{V}(X))\geq 1-\alpha,\] as was to be shown.

6 False discovery exceedance control↩︎

6.1 Multiple testing criteria involving false discovery proportions↩︎

When we control the familywise error rate, we are very strict in the sense that we require that with large probability, there are no false positives at all. In many cases however, we would be willing to incur a small percentage of false positives if this leads to more power. In that case it can be useful to consider multiple testing methods that somehow control the false discovery proportion (FDP), which is defined as the fraction of true hypotheses among all the rejections. We define the FDP to be 0 when there are no rejections.

When we apply a multiple testing method, we will usually denote by \(\mathcal{R}\subseteq\{1,...,m\}\) the set containing the indices of the rejected hypotheses and we write \(R=|\mathcal{R}|\). As usual, we denote the number of false positives by \(V=|\mathcal{R}\cap \mathcal{N}|\). the FDP then satisfies \[FDP=\frac{V}{R\vee 1},\] where “\(\vee 1\)” is added to avoid dividing by 0.

There are several possible approaches to “keeping the FDP small”. A very popular approach is to consider the false discovery rate (FDR), which is the expected value of the FDP, \[FDR=\mathbb{E}(FDP).\] Controlling the FDR means ensuring that the FDR is below some value, which we typically denote by \(\alpha\in(0,1)\). Typical choices are \(\alpha=0.05\) and \(\alpha=0.1\). Likewise we could control the median of the FDP, which means ensuring that \[\mathbb{P}(FDP\leq \gamma)\geq 0.5\] for some small value \(\gamma\in[0,1)\).

If we control the mean or median of the FDP, then we control a central tendency of the FDP, but we do not guarantee that with e.g. \(95\%\) confidence, the FDP lies below some small value \(\gamma\). Ensuring that with large probability the FDP is at most \(\gamma\), is called false discovery exceedance (FDX) control (or sometimes simply “FDP control”). Thus, we then ensure that \[\mathbb{P}(FDP\leq \gamma)\geq 1-\alpha\] for some small value \(\alpha\). Note that if we take \(\alpha=0.5\), we control the median of the FDP. Note that a more “positive formulation” of FDX control, is that we ensure that with probability at least \(1-\alpha\), at least \((1-\gamma)100\%\) of the findings are true discoveries: \[\mathbb{P}(TDP\geq 1-\gamma)\geq 1-\alpha,\] where \(TDP:=1-FDP\) is the true discovery proportion.

FDX control is the topic of this section. Providing confidence statements on many FDPs simultaneously, is discussed in §7. There we will allow picking any set of hypotheses post hoc and obtaining a confidence statement on the corresponding FDP (or equivalently, on the TDP). FDR controlling methods are discussed in §8 and §9.

6.2 Lehmann-Romano and related methods↩︎

We will discuss several FDX methods, but we will omit proofs, which can be found in the cited papers. We first define an FDX method by 1. Suppose that at least one of the following assumptions holds:

• the null p-values satisfy \(\mathbb{P}(P_i\leq c)\leq c\) for every \(c\in[0,1]\), conditional on the non-null p-values.

• the null p-values satisfy Simes’ inequality.

Note that the first assumption is satisfied if the \((P_i:i\in \mathcal{N})\) are independent of the \((P_i:i\in \mathcal{N}^c)\). If at least one of the above two assumptions holds, then the p-value-based step-down method (see §3.5) based on critical values \[\label{eq:LRpowerful} c_i = \frac{(\lfloor \gamma i \rfloor +1 )\alpha}{m+\lfloor \gamma i \rfloor +1-i }\tag{18}\] provides FDX control [29, pp. 1148, 1150]. Recall that in §3.5 we defined what we mean by a step-down method: we need to find the smallest sorted p-value that exceeds its critical value, and then reject all hypotheses with smaller p-values. The Lehmann-Romano procedure that we just defined is also the Lehmann-Romano procedure studied in [45], who make some improvements for the case that the data are discrete. The procedure of [29, pp. 1148, 1150] coincides with the FDX method of [43] (see §6.3) if within that procedure the step-down k-FWER method from [29, p. 1143] is used. The procedure of [29, pp. 1148, 1150] is implemented in the function continuous.LR() from the R package FDX on CRAN. The improvement by [45] for discrete data is available through the function discrete.LR(). Take note that the argument alpha required by those R functions is our \(\gamma\), and their argument zeta is our \(\gamma\). ⁵

[29, p. 1152] also provide a method that does not require dependence assumptions at all. This procedure is the same as their other step-down method just discussed, except that we should decrease the critical values \(c_i\) by dividing them by \(C_{(\lceil \gamma m \rceil +1)}\), where \(C_j := \sum_{i=1}^j (1/i)\). Thus the critical values are \[\label{eq:LRgeneral} c_i' = \frac{(\lfloor \gamma i \rfloor +1 )\alpha}{(m+\lfloor \gamma i \rfloor +1-i)C_{(\lceil \gamma m \rceil +1)} }\tag{19}\] For example, for \(j=10\), it holds that \(C_j\approx 2.92\) and for \(j=100\), it holds that \(C_j\approx 5.187\).

We have seen that there are two “Lehmann-Romano (2005)” methods. The second one that we mentioned (with ciritical values \(c_i'\)) was later uniformly improved in [46], Theorem 3.4. That method makes no assumptions on the dependence either, and works by decreasing the critical values 18 not by \(C_{(\lceil \gamma s \rceil +1)}\), but by a number \(D(\gamma, m)\geq 1\), which depends on \(\gamma\) and \(m\). Thus, the critical values are \[\label{eq:RS} c_i'' = \frac{(\lfloor \gamma i \rfloor +1 )\alpha}{(m+\lfloor \gamma i \rfloor +1-i)D(\gamma, m)}.\tag{20}\] For example, for \(\gamma=0.1\) and \(m=1000\), it holds that \(D(\gamma, m)\approx 3.4179.\) [46, p. 42] contains a table that shows \(D(\gamma, m)\) for several other values of \(\gamma\) and \(m\). Clearly, the method from [46], Theorem 3.4 (which makes no assumptions on the dependence) is still less powerful than [29, pp. 1148, 1150] (which does require assumptions on the dependence).

The method from [46], Theorem 3.4 is not implemented in the R package FDX. However, note that if we multiply all p-values by \(D(\gamma, m)\) and then apply the function continuous.LR(), then this is equivalent to applying the method from [46], Theorem 3.4. This is because multiplying all p-values by \(D(\gamma, m)\) has the same effect as dividing the critical constants by \(D(\gamma, m)\).

6.3 Romano-Wolf↩︎

The FDX controlling method of [43] is based on iteratively applying k-FWER methods: the method starts with applying a 1-FWER method, then applies a 2-FWER method, etcetera until a step is reached where some criterion is not met. Romano-Wolf is thus a kind of step-down approach. The k-FWER method that is applied within Romano-Wolf could e.g. be one of the methods discussed in §3.8 and §5.6. Thus, Romano-Wolf is a general method, in the sense that there are multiple options for the k-FWER method that is used within Romano-Wolf.

Given some k-FWER method (for \(k=1,2,...\)), let \(\mathcal{R}^k\) be the set containing the indices of the hypotheses rejected by the k-FWER method, and let \(\mathcal{R}^0=\emptyset\). The Romano-Wolf method is defined as follows:

1. Let \(k^*\) be the smallest \(k\geq 1\) for which \(k/(|\mathcal{R}^k|+1)>\gamma\) (if there is no such \(k\), reject all hypotheses ).

2. Reject all hypotheses with indices in \(\mathcal{R}^{k^*}\).

An equivalent formulation is in Algorithm 15. The intuition behind the method is as follows: if we reject the top \(|\mathcal{R}^k|+1\) hypotheses, then (with probability at least \(1-\alpha\)) there are at most \(k\) false positives, so the FDP is at most \(k/(|\mathcal{R}^k|+1)\). If this value is larger than \(\gamma\), we stop, because we want the FDP to be at most \(\gamma\).

It turns out that if in Romano-Wolf we use the k-FWER method from Theorem 8, then we obtain the Lehmann-Romano method discussed above (the most powerful one, with critical values defined in 18 ).

In [43] it is proved that Romano-Wolf provides FDX control in an asymptotic sense. Further theoretical results are discussed in [47]. In practice, the Romano-Wolf method tends to be valid if the underlying k-FWER method is valid.

6.4 Application to property sales data↩︎

We continue our analysis from §3.10, where we computed 31 p-values for the numerical predictors of the linear regression model predicting sales prices, see Figure 10. We took \(\alpha=0.05\) there and saw that Bonferroni and Holm both rejected 10 hypotheses. If we keep \(\alpha\) and let \(\gamma=0.1\) and apply the Lehmann-Romano method with critical values defined in 18 , we also reject 10 hypotheses. This is because we require the FDP to be at most \(\gamma\), which means that we allow almost no false hypotheses, given the fact that we are merely testing 31 hypotheses. If we increase \(\gamma\) to 0.2, we reject 12 hypotheses. This is not a real improvement, given that we already rejected 10 hypotheses with Bonferroni.

We see that the results based on Lehmann-Romano are somewhat underwhelming. Part of the reason is “bad luck”: in §6.5 the improvement with Lehmann-Romano is more impressive. Another reason is that in sense Lehmann-Romano is quite strict, since it needs to account for worst-case scenarios regarding the dependence structure of the p-values. (It makes some assumptions on the dependence, but these are limited.) This is somewhat comparable to the way in which Bonferroni and Holm tend to strict compared to maxT. The Romano-Wolf method can be potentially much more powerful than Lehmann-Romano, if a resampling-based k-FWER method is used within Romano-Wolf. Based on the existing literature however, it is not fully clear how such resampling-based methods can be used for testing coefficients in a model. Hence, we now move to a different dataset, where we can use a resampling-based version of Romano-Wolf.

6.5 Application to riboflavin production data↩︎

We analyze a dataset about riboflavin (vitamin B2) production with bacteria called B. subtilis. This dataset is freely available [48]. It contains normalized measurements of expression rates of 4088 genes from \(n = 71\) samples. Further, the dataset contains the corresponding 71 measurements of the logarithm of the riboflavin production rate. For each \(1 \leq i \leq 4088\), we are interested in the hypothesis \(\mathcal{H}_i\) that the riboflavin production rate was stochastically independent of the expression level of gene \(i\). For the parametric methods, we computed p-values using two-sided correlation tests as in §2.3.1.

We took \(\alpha=0.05\) again. Bonferroni and Holm rejected 53 hypotheses, while Hommel rejected 54 hypotheses. Next, we kept \(\alpha=0.05\) and took \(\gamma=0.1\) and applied the FDX method of Lehmann Romano (based on critical values 18 ). This method resulted in 103 rejections. Since \(\gamma\) was 0.1, we know that with probability at least \(1-\alpha=0.95\), at most 10 of these 103 rejections are false positives. Thus with probability at least \(1-\alpha=0.95\), at least 93 of these 103 rejections are true positives.

We can have even more power than Lehmann-Romano if we use a permutation-based multiple testing method. Indeed, we can permute the observations of the roboflavin production rate and construct a matrix of resampled test statistics as in §5.1. We can then perform permutation-based methods such as maxT and the k-FWER verision of maxT from §5.6. The test statistics that were used were the absolute values of the empirical correlations between the production rates and gene expressions. The maxT method rejected 74 hypotheses. As expected, this is more than what Holm rejects. We can also use the FDX method of Romano and Wolf, where as the k-FWER method we take the k-FWER version of maxT from §5.6. We took \(w=1000\), i.e, we used 999 random permutations plus the original data. The Romano-Wolf method based on the single-step k-FWER method from §5.6 rejected 201 hypotheses. Note that it was quite expected that the resampling-based Romano-Wolf method would reject more than Lehmann-Romano: the k-FWER method used within Romano-Wolf is relatively powerful, because it takes into account the dependence structure of the gene expression levels. From Romano-Wolf we conclude that with probability at least \(1-\alpha=0.95\), among the 201 hypotheses with the largest test statistics, at least \(90\%\), so at least 181 hypotheses, are true discoveries.

6.6 Exercises↩︎

§8 contains some exercises linking FDX control and FDR control.

1. Provide a simple expression for the last critical constant, \(c_m\), of the most powerful method among the two FDX methods from

2. Suppose \(m=7\) and \(\alpha=\gamma=0.2\). Suppose the sorted p-values are \(0.01\), \(0.02\), \(0.04\), \(0.06\), \(0.14\), \(0.36\) and \(0.74\). Determine how many hypotheses are rejected by the most powerful method of the two FDX methods from 1.

3. Note that the FDP is a random variable, which takes values in a subset of \([0,1]\). Assume for simplicity that the FDP is symmetric about its mean (in practice it is often right-skewed). Let \(\gamma\in(0,1)\) and \(\alpha\in(0,1/2]\). Prove that \(\mathbb{P}(FDP\leq \gamma)\geq 1-\alpha\) implies that \(\mathbb{E}(FDP)\leq \gamma\). Discuss informally that that is not generally the case the other way around. Further, discuss informally that the first statement is not generally true either when the FDP is right-skewed and \(\alpha\) is close to 0.

4. Suppose that all p-values are independent of each other. Suppose that the p-values corresponding to the true hypotheses all have the standard uniform distribution. Suppose that the p-values corresponding to the false hypotheses all have the same distribution as well (not necessarily uniform). Let \(t\in (0,1)\) and suppose we simply reject all hypotheses with p-values at most \(t\). Now suppose we keep adding hypotheses to the problem (so that \(m\rightarrow\infty\)), and that the fraction \(\pi_0\) of true hypotheses converges to some constant in \((0,1)\) as \(m\rightarrow\infty\). Write \(FDP=FDP(m)\) to ephasize the dependence of the FDP on \(m\).

   Prove that as \(m\rightarrow\infty\), \(FDP(m)\) converges to an (unknown) constant \(c\in[0,1]\). Discuss that this means that if we consider \(\alpha, \gamma\in(0,1)\) and choose \(t\) such that \(\lim_{m\rightarrow\infty} \mathbb{E}(FDP(m))< \gamma\), then we have \(\lim_{m\rightarrow\infty} \mathbb{P}(FDP(m)\leq \gamma)\geq 1-\alpha\). (Recall that ensuring that \(\mathbb{E}(FDP)\) is bounded, is called false discovery rate control.) Discuss informally whether these statements are still necessarily true if the p-values are dependent on each other.

6.7 Solutions↩︎

1. We have \(c_m = \frac{(\lfloor\gamma m\rfloor +1 )\alpha}{m+ \lfloor\gamma m\rfloor +1-m } = \frac{(\lfloor\gamma m\rfloor +1 )\alpha}{\lfloor\gamma m\rfloor +1} =\alpha.\) Thus, the last critical value is simply \(\alpha\), just like in the methods of Bonferroni and Holm and in Simes’ global test.

2. For every \(1\leq i \leq 7\), the critical value \(c_i\) satisfies \[c_i = \frac{(\lfloor \gamma i \rfloor +1 )\alpha}{7+\lfloor \gamma i \rfloor +1-i }.\] For \(1\leq i \leq 4\), \(\lfloor \gamma i \rfloor=0\), and for \(5\leq i \leq 7\), \(\lfloor \gamma i \rfloor=1\). Thus, we have \[\begin{align} c_1=& \frac{0.2}{7+1-1} = 0.2/7 =1/35,\\ c_2=& \frac{0.2}{7+1-2} = 0.2/6 =1/30,\\ c_3=& \frac{0.2}{7+1-3} = 0.2/5 =1/25,\\ c_4=& \frac{0.2}{7+1-4} = 0.2/4 =1/20,\\ c_5=& \frac{0.4}{7+1+1-5} = 0.4/4 =1/10,\\ c_6=& \frac{0.4}{7+1+1-6} = 0.4/3 =2/15,\\ c_7=& \frac{0.4}{7+1+1-7} = 0.4/2 =1/5. \end{align}\] Note that \(p_{(i)}\leq c_i\) holds for all \(1\leq i \leq 3\), but not for \(i=4\). Thus, 3 hypotheses are rejected, since the method is a step-down procedure. (Note that we could have stopped computing the critical values after computing \(c_4\).)

3. \(\mathbb{P}(FDP\leq \gamma)\geq 1-\alpha \geq 0.5\) means that the median of the FDP is at most \(\gamma\). By symmetry, \(\mathbb{E}(FDP)\) must also be at most \(\gamma\).

   If \(\mathbb{E}(FDP)\leq \gamma\), there can be a \(c>\alpha\) such that with probability at least \(c\), the FDP is strictly larger than \(\gamma\). (When the FDP is right-skewed but \(\alpha\) is small, this is often true as well.) Then \(\mathbb{P}(FDP\leq \gamma)\leq 1-c< 1-\alpha.\)

4. We have \(FDP(m)=V(m)/(V(m)+S(m))\), where \(S(m)\) is the number of correct rejections \(|\{1\leq i \leq m: \mathcal{H}_i\text{ is false and } P_i\leq t \}|\).

   Because of the independence of the p-values and because \(\pi_0\) is asymptotically proportional to \(m\), \(V(m)/m\) converges to a constant \(c_1\), say, and likewise \(S(m)/m\) converges to some constant \(c_2\). We have \[FDP(m) = \frac{V(m)/m}{V(m)/m+S(m)/m}.\] As \(m\rightarrow\infty\), the numerator converges to \(c_1\) and the denominator converges to \(c_1+c_2\), so that \(\lim_{m\rightarrow\infty}FDP(m) = c_1/(c_1+c_2)\), which meas that \(FDP(m)\) converges to a constant \(c\).

   Since \(\lim_{m\rightarrow\infty} \mathbb{E}(FDP(m))< \gamma\), it follows that \(c<\gamma\). Hence, \(\lim_{m\rightarrow\infty} \mathbb{P}(FDP(m)< \gamma)=1 \geq 1-\alpha\).

   When there is dependence among, say, the null p-values, then the variance of \(V\) does not necessarily vanish relative to its mean, so \(V(m)/m\) does not converge to a constant. In general, when there is dependence among the p-values, \(FDP(m)\) may not converge to a constant but may vary substantially about its mean regardless of \(m\). In that case, if \(\alpha\) is small, requiring that the tail probability \(\mathbb{P}(FDP>\gamma)\) is at most \(\alpha\) is often a stricter requirement than requiring that the mean \(\mathbb{E}(FDP)\) is at most \(\gamma\) (also for large \(m\)).

7 Simultaneous confidence on all FDPs↩︎

7.1 Motivation↩︎

In §6.1 we discussed a few multiple testing criteria that involve false discovery proportions (FDPs). FDX and FDR methods (the latter are discussed in §8) have the limitation that the user must make all data analysis decisions before seeing the data, and these methods only provide a statement on a single set of hypotheses. FDX and FDR methods provide the user with a fixed set of rejected hypotheses, and the user has no freedom to explore other sets of hypotheses after that. For example, if an FDX method rejects 15 hypotheses, then it provides information on these hypotheses only. The method then provides no information on the top 5 hypotheses (i.e., the 5 hypotheses with the most extreme test statistics or p-values) or the top 20 hypotheses — apart from any information that might be deduced from the method’s statement about the 15 rejected hypotheses.

For example, suppose \(\alpha=0.1\) and \(\gamma=0.2\) and the FDX method rejects 15 hypotheses. Then we know with \(90\%\) confidence that among the top 15 rejections, there are at most 3 false positives. This suggest that it is likely that among the top 5 rejections there are strictly fewer than 3 false positives. However, the FDX method does not provide any statement on that. Likewise, if we look at the top 20 hypotheses instead of the rejected 15 hypotheses, then the method tells us nothing about these 5 additional hypotheses. Thus, we can only say that with \(90\%\) confidence, there are at most 3+5 true hypotheses among the top 20 hypotheses. As another example, suppose that among the 15 rejections, 4 of them belong together in some way, or that for some other reason we are interested in these 4. Then we cannot say anything about these 4 hypotheses except that with large probability at most 3 of them are true. As a final example, suppose that we are interested in some subset of 3 hypotheses among the rejected 15. Then we cannot say anything about the number of false hypotheses among these 3. The general problem here is that FDX methods only provide a confidence statements on a single set of hypotheses.

Recall the notation \(\mathcal{C}=\{\text{nonempty subsets of }\{1,...,m\}\}\). In this chapter, we define a general approach that provides confidence statements on the FDPs of all \(|\mathcal{C}|=2^m-1\) sets of hypotheses. By the FDP of a set \(\mathcal{I}\in \mathcal{C}\) we will simply mean the fraction of true hypotheses in \(|\mathcal{I}|\), i.e., \(|\mathcal{N}\cap\mathcal{I}|/|\mathcal{I}|\). Thus, we view the FDP as a property of a set of hypotheses; there are thus \(2^m-1\) FDPs in total. The methods in this section provide a \((1-\alpha)\)-confidence upper bound \(t_{\alpha}(\mathcal{I})\) for the the number of true hypotheses \(|\mathcal{N}\cap\mathcal{I}|\) in every set \(\mathcal{I}\). Dividing by \(|\mathcal{I}|\) gives a confidence bound for \(|\mathcal{N}\cap\mathcal{I}|/|\mathcal{I}|\), the FDP for the set \(\mathcal{I}\). In fact, we will guarantee much more, because the bounds \(t_{\alpha}(\mathcal{I})\) are all simultaneously valid with probability at least \(1-\alpha\), i.e., \[\label{GSresult} \mathbb{P} \Bigg[\bigcap_{\mathcal{I}\in\mathcal{C}}\Big\{ |\mathcal{N}\cap \mathcal{I}| \leq t_{\alpha}(\mathcal{I}) \Big\} \Bigg]\geq 1-\alpha.\tag{21}\]

Because with large probability the bounds \(t_{\alpha}(\mathcal{I})\) are all valid, we can in fact pick any \(\mathcal{I}\) — even after seeing the data — and obtain a valid confidence statement on the number of true hypotheses in \(\mathcal{I}\). This allows users the freedom to keep looking at various sets of hypotheses after looking the data, without invalidating the confidence statements. Thus, users then have much more flexibility than with FDX and FDR methods, which tell the user which hypotheses to look at.

Note that since \(t_{\alpha}(\mathcal{I})\) is a \((1-\alpha)\)-confidence upper bound for \(|\mathcal{N}\cap\mathcal{I}|\), it follows that \[\label{eq:dalpha} d_{\alpha}(\mathcal{I}):=|\mathcal{I}|- t_{\alpha}(\mathcal{I})\tag{22}\] is a \((1-\alpha)\)-confidence lower bound for \(|\mathcal{N}^c\cap \mathcal{I}|\), the number of false hypotheses (true discoveries) among \(\mathcal{I}\). (\(t_{\alpha}\) stands for “true” and \(d_{\alpha}\) stands for “discoveries”.) Thus, the following is an equivalent formulation of 21 : \[\label{GSresultTDP} \mathbb{P} \Bigg[\bigcap_{\mathcal{I}\in\mathcal{C}}\Big\{ |\mathcal{N}^c\cap \mathcal{I}| \geq d_{\alpha}(\mathcal{I}) \Big\} \Bigg]\geq 1-\alpha,\tag{23}\] Such a formulation in terms \(d_{\alpha}(\mathcal{I})\) is sometimes preferred since it gives a more positive perspective: we ensure that there are at least a certain number of true discoveries among \(\mathcal{I}\). Thus, instead of focussing on what goes wrong, we then focus on the correct rejections. Nevertheless, this formulation is completely equivalent to 21 .

7.2 General construction of the simultaneous bounds↩︎

We now define bounds \(t_{\alpha}(\mathcal{I})\) that have the property 21 . The construction is closely related to the general closed testing principle that we covered in §4.2; in fact, we will start by considering a closed testing procedure. Suppose that for every \(\mathcal{I}\in\mathcal{C}\) we have defined a valid local test \(\phi_{\mathcal{I}}\). Recall that the corresponding closed testing procedure (CTP) rejects all hypotheses \(\mathcal{H}_{\mathcal{I}}\) for which \(\psi_{\mathcal{I}}=1\), where \[\begin{align} \psi_{\mathcal{I}} = & \min\{\phi_{\mathcal{J}}: \mathcal{I}\subseteq \mathcal{J}\in \mathcal{C}\}\\ = & \mathbb{1}\big(\forall \mathcal{I}\subseteq \mathcal{J}\in \mathcal{C}: \phi_{\mathcal{J}}=1 \big). \end{align}\] As a shorthand, write \[\mathcal{X}:= \{\mathcal{I}\in \mathcal{C}: \mathcal{H}_{\mathcal{I}} \text{ is rejected by the CTP}\} = \{\mathcal{I}\in \mathcal{C}: \psi_{\mathcal{I}}=1\}.\] For every \(\mathcal{I}\in \mathcal{C}\) we now define \(t_{\alpha}(\mathcal{I})\) as follows [49]: \[\label{eq:talpha} t_{\alpha}(\mathcal{I}):=\max\{|\mathcal{J}|: \mathcal{J}\subseteq \mathcal{I}\text{ and } \mathcal{J}\not \in \mathcal{X}\}.\tag{24}\] If there is no nonempty \(\mathcal{J}\subseteq \mathcal{I}\) with \(\mathcal{J}\not \in \mathcal{X}\), we consider this maximum to be 0.

We see that the proof is short, just like the proof from §4.2 of the fact that a closed testing procedure controls the FWER. Both proofs are completely based on the fact with probability \(1-\alpha\), \(\mathcal{H}_{\mathcal{N}}\) is not rejected by its local test. As long as that is the case, all statements of the closed testing procedure and all the bounds \(t_{\alpha}(\mathcal{I})\) are correct.

The following formulation of \(t_{\alpha}(\mathcal{I})\) (by [50]) can also be useful.

\[\label{eq:talpha2} t_{\alpha}(\mathcal{I})=\max\{|\mathcal{I}\cap \mathcal{K}|: \mathcal{K}\in\mathcal{C}\text{ and }\phi_{\mathcal{K}}=0\}.\tag{25}\]

Which formulation of \(t_{\alpha}(\mathcal{I})\) is more useful or insightful, might depend on the situation. The first formulation, 24 , connects the method to the large literature on closed testing. The second formulation, 25 , is in a sense the simplest, since it does not refer to the closed testing procedure. However, note that if we use the second formulation, we need to check for each \(\mathcal{K}\in \mathcal{C}\) whether \(\phi_{\mathcal{K}}=0\), which is very computationally expensive. In the first formulation, we ‘only’ need to check something for the \(\mathcal{J}\) with \(\mathcal{J}\subseteq \mathcal{I}\). This might be useful in case we have a fast algorithm for checking whether \(\mathcal{J}\in \mathcal{X}\).

7.3 Consonant procedures, e.g. Bonferroni-based local tests↩︎

As a first example, suppose the local tests \(\phi_{\mathcal{I}}\) are Bonferroni tests, i.e., \[\label{localtBonf} \phi_{\mathcal{I}} = \mathbb{1}( \min_{i\in \mathcal{I}}P_i\leq \alpha/|\mathcal{I}|).\tag{26}\] We saw in §4.3 that the CTP based on these local tests is the Holm method; more precisely, they reject the same elementary hypotheses (and since the CTP is consonant, it follows which intersection hypotheses are rejected by the CTP). Now, let us see what simultaneous FDP bounds \(t_{\alpha}(\mathcal{I})\) we get based on these local tests.

We thus see that in case Bonferroni-based local tests are used, the bounds \(t_{\alpha}(\mathcal{I})\) follow directly from knowing which elementary hypotheses are rejected. It turns out that this is always the case if the closed testing procedure is consonant. Recall that the CTP is called consonant if for every \(\mathcal{I}\in \mathcal{X}\), there is an \(i\in \mathcal{I}\) with \(\{i\}\in \mathcal{X}\). Note that for every CTP, if \(\{i\}\in \mathcal{X}\) then all \(\mathcal{I}\) that contain \(i\) are in \(\mathcal{X}\). Thus, if the CTP is consonant, all intersection hypotheses that are rejected, are those that logically have to be false if the rejected elementary hypotheses are false. It may not be surprising that \(t_{\alpha}(\mathcal{I})\) is then simply the number of \(i\in\mathcal{I}\) with \(i\in \mathcal{X}\). This is true for all consonant procedures. We now prove that. Then, we show that the CTP based on Bonferroni-based local tests is consonant. These two results together imply Proposition 15.

Another example of a consonant CTP, is the CTP that corresponds to the maxT method (see §5.4). In the next section, we will look at a CTP that is not consonant.

7.4 Nonconsonant procedures, e.g. Simes local tests↩︎

For nonconsonant (“dissonant”) methods, the bounds \(t_{\alpha}(\mathcal{I})\) do not generally follow trivially from the set of rejected elementary hypotheses. An example of a nonconsonant CTP (see exercises) is the one based on Simes local tests \[\label{localtSimes} \phi_{\mathcal{I}} = \mathbb{1}\Big( \bigcup_{i\in \mathcal{I}} \{P^{\mathcal{I}}_{(i)}\leq \frac{i\alpha}{|\mathcal{I}|} \} \Big).\tag{29}\] We mentioned in §4.4 that the CTP based on these local tests is Hommel’s method — more precisely, Hommel finds the elementary hypotheses that are rejected by this CTP. Since the CTP is not consonant, we can potentially make statements that do not trivially follow from Hommel’s rejections of elementary hypotheses. For example, it may happen that Hommel rejects nothing, but we can still infer that \(|\mathcal{N}^c\cap\mathcal{I}|>0\) for some \(\mathcal{I}\in \mathcal{C}\).

Consider the CTP based on the Simes local tests. Naively computing the \(t_{\alpha}(\mathcal{I})\) for all \(\mathcal{I}\in \mathcal{C}\) is computationally infeasible for moderate or large \(m\). A fast, exact computational shortcut for this CTP is presented in [51]. The method is implemented in the R package Hommel on CRAN [35]. The R package can be used as follows. First create an object of class hommel using \[\verb|hom <- hommel(p)|\] where p contains the \(m\) p-values. Then run \[\verb|discoveries(hom, I, alpha=0.05)|\] where I is the set of (indices of) hypotheses \(\mathcal{I}\) for which we want to compute the bound \(d_{\alpha}(\mathcal{I})=|\mathcal{I}|-t_{\alpha}(\mathcal{I}).\) We then obtain \(d_{\alpha}(\mathcal{I})\).

If we apply this method (with \(\alpha=0.05\) again) to the 4088 p-values from §6.5, we clearly obtain valuable information on top of what Hommel’s FWER method provided. Indeed, we saw that Hommel’s method rejected 54 hypotheses. However, if we compute \(d_{\alpha}(\mathcal{I})\) with \(\mathcal{I}=\{1,...,m\}\), we get 104 rejections. Thus, we know that with probability at least \(0.95\), there are at least 104 false hypotheses. Further, if we compute \(d_{\alpha}(\mathcal{I})\) where \(\mathcal{I}\) contains the hypotheses with the 100 smallest p-values, then we get \(d_{\alpha}(\mathcal{I})=90\); thus, among the 100 hypotheses with the smallest p-values, at least 90 are false. Importantly, all such statements are simultaneously true with probability at least \(0.95\).

If we consider the p-values from §3.10, then unfortunately the bounds \(d_{\alpha}(\mathcal{I})\) provide no additional information in addition to the set of elementary hypotheses rejected by Hommel. Thus, we see that the information that we get on top of the rejected elementary hypotheses can sometimes be a lot, and sometimes nothing.

Other examples of nonconsonant CTPs are procedures with local tests based on sums of test statistics, including Fisher combinations. As is the case for FWER and FDX control, sometimes we can gain more powerful simultaneous FDP statements by using approaches based on permutations or bootstrapping. Examples are the methods in [52]–[55].

7.5 Exercises↩︎

1. Consider Figure 12, where crosses indicate the intersection hyptheses that are rejected by their local tests. Suppose the local tests have significance level \(\alpha=0.05\).

   (a) What can we say with \(95\%\) confidence about the total number of false hypotheses among the three?

   Now suppose that \(\mathcal{H}_2\cap\mathcal{H}_3\) is also rejected by its local test, in addition to the hypotheses indicated with a cross.

   (b) What can we then say about the total number of false hypotheses in \(\mathcal{H}_2\cap\mathcal{H}_3\)?

   (c) What can we then say about the total number of false hypotheses?

2. Suppose we take \(\alpha=0.01\) and are interested in four hypotheses \(\mathcal{H}_1,...,\mathcal{H}_4\) and the \(2^4-1\) intersections. Consider some CTP for these hypotheses. Suppose that all intersection hypothes are rejected by the CTP, except \(\mathcal{H}_{12}\), \(\mathcal{H}_1\) and \(\mathcal{H}_2\).

   (a) What can we say with \(99\%\) confidence about the total number of false hypotheses among \(\mathcal{H}_1,...,\mathcal{H}_4\)?

   (b) What can we say with \(99\%\) confidence about the total number of false hypotheses among \(\mathcal{H}_3\) and \(\mathcal{H}_4\)?

3. Consider elementary hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) and a corresponding CTP. Prove that the following statements are equivalent.

   • The procedure is consonant.

   • For every \(\mathcal{I}\in\mathcal{C}\), \(d_{\alpha}(\mathcal{I})>0\) implies there exists an \(i\in\mathcal{I}\) for which \(d_{\alpha}(\{i\})=1\).

7.6 Solutions↩︎

1. (a) For the first question, we must find \(t_{\alpha}(\{1,2,3\})\), which is defined by \[\max\{|\mathcal{J}|: \mathcal{J}\subseteq \{1,2,3\} \text{ and } \mathcal{J}\not\in \mathcal{X}\}.\] Note that \(\{1,2,3\}\in \mathcal{X}\), so this maximum is smaller than 3. Note that \(\{2,3\}\not\in \mathcal{X}\), so this maximum is 2. Thus, with probability at least \(0.95\) the number of true hypotheses among the three is at most 2, i.e., there is at least one false hypothesis among these three.

   (b) Note that in this case \(\mathcal{H}_1\cap\mathcal{H}_2\), \(\mathcal{H}_1\cap\mathcal{H}_3\) and \(\mathcal{H}_2\cap\mathcal{H}_3\) are all rejected by the CTP. It follows that \(t_{\alpha}(\{1,2,3\})\) is at most 1. It is not 0, since \(\{3\}\not\in \mathcal{X}\). Hence with probability at least \(0.95\) the number of true hypotheses among the three is at most 1, i.e., there is are least 2 false hypotheses among these three.

   (c) For the third question, we look at \(t_{\alpha}(\{2,3\})\), which is 1. Hence, with probability at least \(0.95\), the number of true hypotheses among these two is also at most 1.

2. (a) We must compute \[d_{\alpha}(\{1,...,4\})= 4- t_{\alpha}(\{1,...,m\})=\] \[4-\max\{|\mathcal{J}|: \mathcal{J}\subseteq \{1,...,4\} \text{ and } \mathcal{J}\not \in \mathcal{X}\}=\] \[4-|\{1,2\}|=2.\] We see that with probability at least \(0.99\), there are at least two false elementary hypotheses.

   (b) Note that \[t_{\alpha}(\{3,4\})=\] \[2-\max\{|\mathcal{J}|: \mathcal{J}\subseteq \{3,4\} \text{ and } \mathcal{J}\not \in \mathcal{X}\}.\] There is no \(\mathcal{J}\subseteq \{3,4\}\) for which \(\mathcal{J}\not \in \mathcal{X}\). Hence we know that with probability at least \(0.99\), \(\mathcal{H}_3\) and \(\mathcal{H}_4\) are both false.

3. We first show that the first statement implies the second one. \(d_{\alpha}(\mathcal{I})>0\) means \(t_{\alpha}(\mathcal{I})<|\mathcal{I}|\). Hence, apparently, \(\mathcal{I}\in \mathcal{X}\), since otherwise we would have \(t_{\alpha}(\mathcal{I})=|\mathcal{I}|\). Since the procedure is consonant, this means that there must be an elementary hypothesis \(i\in \mathcal{I}\) that is rejected by the CTP, i.e., \(\{i\}\in \mathcal{X}\). This means that \(t_{\alpha}(\{i\})=0\), so \(d_{\alpha}(\{i\})=1-0=1.\)

   We now show that the second statement implies the first one. Consider \(\mathcal{I}\in \mathcal{X}\). We must show there is an \(i\in \mathcal{I}\) that is rejected by the CTP. Note that \(t_{\alpha}(\mathcal{I})<|\mathcal{I}|\), so \(d_{\alpha}(\mathcal{I})>0\), so there is an \(i\in \mathcal{I}\) with \(d_{\alpha}(\{i\})=1\), so \(t_{\alpha}(\{i\})=0\), which implies that \(\{i\}\in \mathcal{X}\). We conclude that the procedure is consonant.

8 False discovery rate control with Benjamini-Hochberg↩︎

Recall that the false discovery rate (FDR) is defined as \(\mathbb{E}(FDP)\). The by far most popular and well-known FDR controlling method is the Benjamini-Hochberg (BH) method. The method was proposed in [56]. There the authors prove that the method controls the FDR when the null p-values are independent of each other. In fact, they show that the FDR is then at most \(\pi_0\alpha\), where \(\pi_0:=|\mathcal{N}|/m\) is the fraction of true hypotheses among all hypotheses. (Unfortunately, \(\pi_0\) is usually not exactly known.) In [57] it is shown that this still holds under certain dependence structures of the p-values (see below). That paper also provides a variant of the method, which controls the FDR regardless of the dependence structure, but is less powerful.

BH is defined as follows. Let \(k\) be the largest number \(1\leq i \leq m\) for which \(P_{(i)}\leq \frac{i}{m}\alpha\), if there is such an \(i\), and reject the hypotheses corresponding to \(P_{(1)},...,P_{(k)}\). Otherwise, reject nothing. Thus, BH is simply a step-up method with critical values \(c_i=\frac{i}{m}\alpha\), \(1\leq i \leq m\). This means that BH-adjusted p-values can be computed as explained in §3.9.

8.1 Theory↩︎

When \(m\) is large and a substantial proportion of the hypotheses are false, BH can be much more powerful than FWER methods such as Bonferroni, Holm and Hommel. For example, in BH, the \((m/10)\)-th sorted p-value is compared with \(\alpha/10\) (assuming \(m\) is a multiple of 10), while Bonferroni compares all p-values with \(\alpha/m\). Of course, ensuring that \(FDR\leq\alpha\) is a much weaker guarantee than ensuring that \(FWER\leq\alpha\). Also, ensuring that the FDR is below \(0.05\) is weaker than ensuring that the FDP is below \(0.05\) with high probability (FDX control), see some of the exercises.

BH often controls the FDR in practice, especially when two-sided p-values are used. BH is known to control the FDR when the PRDS assumption is satisfied. PRDS stands for positive regression dependence on the subset \(\mathcal{N}\). The PRDS assumption roughly says that all p-values depend on the null p-values in a positive way. To define the PRDS assumption, we need the following concept. A set \(\mathcal{A}\subseteq\mathbb{R}^m\) is called increasing if \(a\in \mathcal{A}\) implies that \(b\in \mathcal{A}\) for all \(b\in \mathbb{R}^m\) satisfying \(b\geq a\) (meaning \(b_i\geq a_i\) for all \(1\leq i \leq m\)).

We will only prove the first part of the theorem here.

The proof is by [59] and is also provided in [58]. ([57] prove a slightly different result.) Note that the critical values of BH are identical to the critical values of Simes’ global test. From the fact that BH provides weak FWER control under the PRDS assumption (see exercises) it follows that Simes’ inequality 7 is valid under the PRDS assumption. The reverse is not always true, i.e., Simes may hold while PRDS does not hold. When there are strong positive dependencies in the data, BH and Simes’s global test tend to be conservative.

BH is implemented in the standard R function p.adjust(). The mentioned conservative variant that controls the FDR regardless of the dependence structure [57] is also available in that function. Applying BH (with \(\alpha=0.05\)) to the 279 p-values from §3.10 results in 58 rejections (Holm rejected 36). Applying it to the 31 p-values corresponding to the numerical covariates results in 17 rejections (Holm rejected 10).

8.2 Exercises↩︎

1. Suppose \(m=5\), \(\alpha=0.05\) and the p-values corresponding to \(\mathcal{H}_1,...,\mathcal{H}_5\) are \(0.200\), \(0.015\), \(0.025\), \(0.030\), and \(0.720\) respectively.

   (a) Which hypotheses are rejected by Benjamini-Hochberg?

   (b) Compute the Benjamini-Hochberg-adjusted p-value for each hypothesis. See §3.9 on how to compute adjusted p-values for step-up methods.

2. Consider \(m\) hypotheses. Below, “hypotheses” always refers to these \(m\) hypotheses and not to any of their intersections.

   For each of the following general statements, indicate whether they are true. Explain your claims.

   • Statement 1: Whenever Benjamini-Hochberg rejects 0 hypotheses, Hommel also rejects 0 of the hypotheses.

   • Statement 2: Whenever Hommel rejects 0 hypotheses, Benjamini-Hochberg also rejects 0 hypotheses.

3. Consider a multiple testing method that controls the FDR at level \(\alpha\), i.e., \(FDR\leq\alpha\), where \(\alpha\in(0,1)\). Prove that the method provides weak FWER control at level \(\alpha\).

4. Consider a multiple testing method that controls the FDR, i.e., \(FDR\leq\alpha\), where \(\alpha\in(0,1)\). Let \(\gamma\in (0,1)\). Prove that \(\mathbb{P}(FDP\geq \gamma) \leq \alpha/ \gamma\), i.e., the method provides FDX control.

5. Suppose a method provides FDX control, i.e., \(\mathbb{P}(FDP\leq \gamma) \geq 1-\alpha\). Show that then \(FDR\leq (1-\alpha)\cdot \gamma + \alpha\).

6. Let \(\mathcal{E}\) be the event that Benjamini-Hochberg rejects at least 1 hypothesis. Consider the following statement: “If the p-values are valid and satisfy the PRDS assumption, BH controls the FDR conditionally on \(\mathcal{E}\), i.e., \(\mathbb{E}(FDP|\mathcal{E})\leq \alpha\)”. Is this statement generally true?

8.3 Solutions↩︎

1. (a) The critical constants are \(0.01, 0.02,...,0.05\). Note that \(p_{(3)}\leq 0.03\), so at least 3 hypotheses are rejected, since Benjamini-Hochberg is a step-up method. The p-values \(p_{(4)}\) and \(p_{(5)}\) do not exceed their critical values. Hence, exactly three hypotheses are rejected. The hypotheses with the smallest p-values are \(\mathcal{H}_2, \mathcal{H}_3, \mathcal{H}_4\), so these hypotheses are rejected.

   (b) First we must compute \(p_{(i)}\alpha c_i^{-1} =p_{(i)}\alpha ((\alpha i/m))^{-1} = p_{(i)}m/i\) for every \(1\leq i \leq 5\). This gives the values \(0.075, 0.062, 0.050, 0.250, 0.720\). Finally, we must enforce monotonicity. Since Benjamini-Hochberg is a step-up method, this means decreasing the values 0.075 and 0.062 to 0.050. Thus, the sorted adjusted p-values are \(0.050, 0.050, 0.050, 0.250, 0.720\). Thus, for \(\mathcal{H}_1,...,\mathcal{H}_5\), the respective adjusted p-values are \(0.250, 0.050, 0.050, 0.050, 0.720\).

2. The critical values of the step-up method of Benjamini and Hochberg are identical to the critical values of Simes’ global test. Consequently, Benjamini-Hochberg rejects something if and only if Simes’ global test rejects.

   Hommel is the CTP based on Simes local tests. A CTP will reject none of the elementary hypotheses if the global hypothesis \(\mathcal{H}_1\cap...\cap \mathcal{H}_m\) is not rejected by its local test. Hence, Hommel will reject nothing if Simes’ global test does not reject. Thus, if Benjamini-Hochberg rejects nothing, then Hommel rejects nothing, so Statement 1 is true.

   If Simes’ global test rejects something, then that does not imply that Hommel rejects one of the elementary hypotheses. This is because the CTP based on Simes local tests is not consonant. Thus, it may happen that Benjamini-Hochberg rejects something, but Hommel rejects none of the \(m\) hypotheses. Thus, Statement 2 is false.

3. Suppose all hypotheses are true. We must show that \(\mathbb{P}(V>0)\leq \alpha\). Note that if all hypotheses are true, then \(V>0\Leftrightarrow FDP=1\) and \(V=0\Leftrightarrow FDP=0\). Hence, \(\mathbb{P}(V>0) =\mathbb{P}(FDP=1) = \mathbb{E}(FDP)\leq\alpha\).

4. Since \(FDP\) is a nonnegative random variable, we can use Markov’s inequality, which says that \(\mathbb{E}(FDP)\geq \gamma \mathbb{P}(FDP\geq \gamma)\). Hence \[\alpha \geq FDR=\mathbb{E}(FDP)\geq \gamma\mathbb{P}(FDP\geq \gamma)\] and the result follows.

5. Note that by the law of total probability, \(\mathbb{E}(FDP)\) equals \[\mathbb{P}(FDP\leq \gamma)\mathbb{E}(FDP|FDP\leq \gamma) + \mathbb{P}(FDP> \gamma)\mathbb{E}(FDP|FDP>\gamma)\] \[\label{eqproofFDRFDX} \leq \mathbb{P}(FDP\leq \gamma) \cdot \gamma + \mathbb{P}(FDP> \gamma) \cdot 1.\tag{31}\] Here we know that \(\mathbb{P}(FDP> \gamma)\) lies in the interval \([0,\alpha]\). 31 is maximized when \(\mathbb{P}(FDP\leq \gamma)\) is minimized and \(\mathbb{P}(FDP> \gamma)\) is maximized, i.e., when \(\mathbb{P}(FDP> \gamma)=\alpha\). Thus, 31 cannot exceed \[(1-\alpha)\cdot \gamma + \alpha,\] as was to be shown.

   Alternative solution: Note that 31 equals \[(1-\mathbb{P}(FDP> \gamma)) \cdot \gamma + \mathbb{P}(FDP> \gamma)\] \[=\mathbb{P}(FDP> \gamma)(1-\gamma)+\gamma\] \[\leq \alpha (1-\gamma)+\gamma\] \[= (1-\alpha)\cdot \gamma + \alpha.\]

6. No. Suppose for example that all hypotheses are true. Then conditional on \(\mathcal{E}\), the FDP is always 1. Thus, conditional on finding something, the error rate is not generally controlled. This is also true for many other multiple testing methods, and also for tests of a single hypothesis.

9 FDR control with knockoffs↩︎

9.1 Motivation and setting↩︎

Consider covariates \(X_1,...,X_m\) and a response variable \(Y\). Often we want to know for every \(1\leq j \leq m\) if \(Y\) depends on \(X_j\) conditional on the other \(m-1\) covariates. For every \(1\leq j \leq m\), let \(\mathcal{H}_j\) be the hypothesis that \(Y\) is conditionally independent on \(X_i\). If \(\mathcal{H}_j\) is true, then that means that given the other covariates, \(X_j\) does not provide any information on \(Y\); in short, \(Y\mathrel{\perp\!\!\!\perp}X_j|X_{-j}\).⁷

To test \(\mathcal{H}_j\), we could for instance assume a generalized linear model where \(\mathbb{E}(Y)= g^{-1}(\beta_1X_1+...+\beta_mX_m)\), and test whether the coefficient of \(X_i\) is 0. Of course, we then make assumptions, since we assume a certain distributional shape for \(Y\) and we assume that \(g(\mathbb{E}(Y))\) depends linearly on \(\beta_1X_1+...+\beta_mX_m\). We could make the model potentially more realistic by adding interactions, but still some substantial assumptions remain, and we should also be careful that the number of parameters does not become too large compared to the sample size.

The knockoffs methodology from [60] provides an alternative method that gets rid of such assumptions, replacing it by assumptions on the joint dependence structure of \((X_1,...,X_m)\). The method tests the conditional-independence hypotheses \(\mathcal{H}_1,...,\mathcal{H}_m\) in such a way that the FDR is controlled (for the knockoffs+ method, to be precise). The method puts no restrictions on \(m\): it may be larger than \(n\). The only assumption that the method makes, is that the joint dependence structure of \((X_1,...,X_m)\) is known. This is of course a big assumption, but in some settings this dependence structure is (approximately) known. [60, pp. 554–555] give these examples where the dependence structure is (approximately) known:

• Sometimes the covariate distribution is known because we control it, e.g. in some experiments.

• Sometimes we very accurately know the covariate distribution, due to the availability of a large amount of additional unlabeled data.

• A third example is when we have pre-existing knowledge about aspects of the joint distribution of the covariates.

Moreover, if the covariate distribution is unknown, it can be estimated based on the data at hand; there are some limited results on knockoffs methods for the situation where the dependence structure of \((X_1,...,X_m)\) is unknown a priori.

The knockoffs method is a multiple testing method, but it can of course be used as a variable selection method. We could therefore denote the set of hypotheses that the method selects by \(\mathcal{S}\) as in [60], but we will use the notation \(\mathcal{R}\) in line with our usual notation. Thus, the (knockoffs+) method guarantees that \(FDR=\mathbb{E}(|\mathcal{N}\cap \mathcal{R}|)/(|\mathcal{R}|\vee 1)\leq \alpha\).

9.2 Definition of knockoffs and test statistics↩︎

Note that property [it:XindY] is automatically satisfied if \(\tilde{X}\) is constructed based on \(X\) without using \(Y\). The vector \((X,\tilde{X})_{\text{swap}(\mathcal{A})}\) is obtained from \((X,\tilde{X})\) by swapping the \(X_i\) with \(i\in \mathcal{A}\) with the corresponding \(\tilde{X}_i\). For example, if \(m=3\) and \(\mathcal{A}=\{2,3\}\), we have \[(X_1,X_2,X_3,\tilde{X}_1,\tilde{X}_2,\tilde{X}_3)_{\text{swap}(\mathcal{A})} = (X_1,\tilde{X}_2,\tilde{X}_3,\tilde{X}_1,X_2,X_3).\] By property [it:swapinv], the distribution of \((X,\tilde{X})\) is invariant under these swapping operations. Thus the knockoff variables are interchangeable with the corresponding original covariates, except that the knockoff variables are conditionally independent of \(Y\). We will discuss the construction of knockoff variables in §9.4.

We now define the data in detail. We consider i.i.d. random vectors \((X_{i1},...,X_{im}, Y_i)\in \mathbb{R}^m\times\mathbb{R}\), which we collect in an \(n\times m\)-matrix \(\mathbf{X}\) and a column vector \(y\) of length \(n\). The knockoffs matrix \(\tilde{\mathbf{X}}\) is constructed in such a way that for each \(1\leq i \leq n\), \((\tilde{X}_{i1},...,\tilde{X}_{im})\) is a knockoff of \((X_{i1},...,X_{im})\) as explained above. We define \((\mathbf{X}, \tilde{\mathbf{X}})_{\text{swap}(\mathcal{A})}\) to be the matrix for which the swapping operation defined above is applied to each row. We need the following lemma.

The knockoffs method takes certain statistics \(W_j\), \(1\leq j \leq m\), as input. We now explain the construction of these \(W_j\), which are a function (\(w_j(\cdot)\), say) of the data and the knockoffs: \[W_j=w_j\big\{(\mathbf{X}, \tilde{\mathbf{X}}),y\big\}.\] In order for the knockoffs method to work, we need to define the test statistics \(W_j\) in such a way that they have the sign flip property, which means that swapping the \(j\)-th variable with its knockoff changes the sign of \(W_j\): \[\label{eq:signflip} w_j\big\{(\mathbf{X}, \tilde{\mathbf{X}})_{\text{swap}(\mathcal{A})},y\big\} = \begin{cases} w_j\big\{(\mathbf{X}, \tilde{\mathbf{X}}),y\big\}, & j\not\in \mathcal{A}, \\ -w_j\big\{(\mathbf{X}, \tilde{\mathbf{X}}),y\big\}, & j\in \mathcal{A}, \end{cases}\tag{32}\] for every \(\mathcal{A}\subseteq \{1,...,m\}\).

Apart from requiring the sign flip property, how should we define these statistics \(W_j\) to have good power? The key idea is that we define them in such a way that if \(\mathcal{H}_j\) is false, then \(w_j\big\{(\mathbf{X}, \tilde{\mathbf{X}}),y\big\}\) tends to be (much) larger than 0. The reason why we want that, is that the procedure will reject all hypotheses with test statistics exceeding some value. For example [60, p. 561], consider a linear regression model with \(2m\) coefficients, corresponding to our \(m\) covariates and their knockoffs — where \(Y_i\) depends on the \(i\)-th row of \((\mathbf{X}, \tilde{\mathbf{X}})\). (Note that do not require this model to be valid. We will only need the assumptions stated above to achieve FDR control.) For this model, consider the lasso,

\[\min_{b\in \mathbb{R}^{2m}} \frac{1}{2} \lVert y - (\mathbf{X}, \tilde{\mathbf{X}}) b\rVert^2_2 +\lambda \lVert b \rVert_1\] and let \(\hat{b}\) be the solution. Here \(\lambda>0\) can be prespecified or based on the data, as long as we define it in such a way that permuting the columns of \((\mathbf{X}, \tilde{\mathbf{X}})\) does not change its value. Then we can define the test statistics \(W_j\) based on the lasso-estimated coefficients in the following way: \[W_j = |\hat{b}_j| - |\hat{b}_{j+m}|.\] With this definition in place, the sign flip property 32 is satisfied. If instead of the lasso we use ridge or elastic net to get \(\hat{b}\), then likewise property 32 is satisfied. These are common ways to define the test statistics, but other definitions satisfying property 32 are also possible. For example, the test statistics may be based on a generalized linear model instead of a linear model. They may also be based on variable importance measures from some machine learning model.

The following lemma [60, p. 561] states the key property that we will use to prove FDR control. Note that a Rademacher variable takes the values \(-1\) and \(1\), each with probability \(0.5\).

9.3 FDR control↩︎

Consider some threshold \(t>0\). By Lemma 2, \[|\{1\leq j \leq m: W_j\leq -t\}|\geq |\{j\in \mathcal{N}: W_j\leq-t\}| \,{\buildrel d \over =}\, |\{j\in \mathcal{N}: W_j\geq t\}|.\] Note that the above equality in distribution holds because conditionally on the indices of “extreme” null statistics \(\{j\in \mathcal{N}: |W_j|\geq t\}\), the signs of these extreme null statistics are i.i.d. Rademacher variables. Thus, if we reject all hypotheses \(\mathcal{H}_j\) with \(W_j\geq t\), then we see that \(|\{1\leq j \leq m: W_j\leq -t\}|\) is an upward biased estimate of \(|\{j\in \mathcal{N}: W_j\geq t\}|\), which is the number false discoveries. Thus, we can view \[\widehat{FDP}(t):= \frac{|\{1\leq j \leq m: W_j\leq-t\}|}{|\{1\leq j \leq m: W_j\geq t\}|}\] as an estimate of the FDP when we use threshold \(t\). The idea of the knockoffs procedure is to choose the rejection threshold \(t\) in a data-dependent way. Namely, we take \(t\) to be the smallest threshold for which the estimated \(FDP\) is still below \(\alpha\).

The following theorem says that if we use the estimate \(\widehat{FDP}(t)\) defined above then we do not control the FDR but something which is similar (if \(m\) is large). To ensure real FDR control, we have to add \(+1\) in the numerator of \(\widehat{FDP}(t)\), which makes the estimate more conservative [60, p. 562]. The proof only makes use of the coin-flip property from Lemma 2 and is given in [61]. Here we give a shortened version of the proof, sometimes referring to [61] for details.

9.4 Constructing knockoffs↩︎

How can we construct knockoff variables, i.e. variables \(\tilde{X}\) such that the properties in Definiton 3 hold?

9.4.1 The case that \(X\) is multivariate normal↩︎

First consider the case that \(X\) has a known multivariate normal distribution \(\mathcal{N}(0,\mathbf{\Sigma})\). Then an example of a joint distribution for \((X, \tilde{X})\) that satisfies property [it:swapinv] is \(\mathcal{N}(0,\mathbf{G})\), where \[\label{refG} \mathbf{G} = \begin{pmatrix} \mathbf{\Sigma}& \mathbf{\Sigma}-\text{diag}(s) \\ \mathbf{\Sigma}-\text{diag}(s) & \mathbf{\Sigma}, \end{pmatrix}\tag{36}\] where \(\text{diag}(s)\) is the diagonal matrix with diagonal \(s=(s_1,...,s_m)'\) and \(s\succeq 0\) (meaning all elements of \(s\) are nonnegative) is chosen in such a way that the matrix is positive semidefinite. Note that by defining \(\mathbf{G}\) in this way, we ensure that for every \(1\leq j \leq m\), \(X_j\) and \(\tilde{X}_j\) have the same covariance with all the other variables⁹, while their covariance with each other is smaller than their variance. That is essential, because if we would have \(\mathop{\mathrm{Cov}}(X_j,\tilde{X}_j) = \mathop{\mathrm{Var}}(X_j)=\mathop{\mathrm{Var}}(\tilde{X}_j)\), then \(X_j\) and \(\tilde{X}_j\) would be perfectly correlated and \(X_j\) and \(\tilde{X}_j\) would be (a.s.) identical. That would mean that the procedure usually has no power for rejecting \(\mathcal{H}_j\).

To sample \(\tilde{X}\) such that \((X, \tilde{X})\sim \mathcal{N}(0,\mathbf{G})\), we can use classical formulas for the conditional distribution of a subvector of a multivariate normal distribution¹⁰. These tell us that \[\tilde{X}|X \,{\buildrel d \over =}\, \mathcal{N}(\mu,\mathbf{V}),\] where \[\begin{align} \mu =& \mathbf{\Sigma}^{-1}(\mathbf{\Sigma}-\text{diag}(s)) X \\ =& X-\mathbf{\Sigma}^{-1}\text{diag}(s)X,\\ \mathbf{V} =& \mathbf{\Sigma}- (\mathbf{\Sigma}-\text{diag}(s))\mathbf{\Sigma}^{-1} (\mathbf{\Sigma}-\text{diag}(s)) \\ =& \mathbf{\Sigma}- \mathbf{\Sigma}+ 2\text{diag}(s) -\text{diag}(s)\mathbf{\Sigma}^{-1}\text{diag}(s) \\ =& 2\text{diag}(s) - \text{diag}(s)\mathbf{\Sigma}^{-1}\text{diag}(s). \end{align}\]

It turns out that \(\mathbf{G}\) is positive semidefinite if and only if \(s\succeq 0\) and \(2\mathbf{\Sigma}\preceq \text{diag}(s)\) [61, pp. 2057, 2062]. (\(\mathbf{A} \succeq \mathbf{B}\) is defined to mean that \(\mathbf{A} - \mathbf{B}\) is positive semidefinite.) How can we choose \(s\) such that this is satisfied and the method has good power? In order to have good power, we want \(\mathop{\mathrm{Cov}}(X_j,\tilde{X}_j)\) to be small, which means we want \(s_j\) to be large. Thus we wish to find large \(s_j\)’s in such a way that the constraints \(s\succeq 0\) and \(2\mathbf{\Sigma}\preceq \text{diag}(s)\) are satisfied. Especially when \(m\) is large, this quicky becomes computationally demanding: in practice we will need to choose a tradeoff between statistical power and computational efficiency. Details are in [60, pp. 564–565].

9.4.2 The general case↩︎

Above we assumed that \(X\) has a known multivariate normal distribution \(\mathcal{N}(0,\mathbf{\Sigma})\). Now suppose we do not have such knowledge. One possible situation is that we know that \(X\) is multivariate normal with mean 0, but we do not know \(\mathbf{\Sigma}\). Then we can estimate \(\mathbf{\Sigma}\) and proceed as above, which leads to approximate FDR control. Another possible situation is that we do not even know whether \(X\) is multivariate normal. Then a proposal by [60] is that we do not require \((X,\tilde{X})_{\text{swap}(\mathcal{A})}\) and \((X,\tilde{X})\) to have exactly the same distribution, but require them to have the same first two moments. (If we have to estimate the covariance matrix \(\mathbf{\Sigma}\), then this moment matching can of course only be done approximately.)

Requiring that \((X,\tilde{X})_{\text{swap}(\mathcal{A})}\) and \((X,\tilde{X})\) have the same covariance matrix, means requiring that \[\mathop{\mathrm{Cov}}(X,\tilde{X})=\mathbf{G},\] with \(\mathbf{G}\) as in 36 . Hence, [60, pp. 564–565] propose to use the same approaches for finding \(s\) as in the case where \(X\) is multivariate normal (see §9.4.1).

Note that when \(\mathbf{\Sigma}\) is unknown and \(m\) is large, estimating \(\mathbf{\Sigma}\) is no small task. This means that the knockoffs(+) method can be quite inaccurate, in the sense that the FDR may substantially exceed \(\alpha\). However, it is worth adding that when \(m\) is large and \(\mathbf{\Sigma}\) is unknown, we might be able to use some pre-existing knowledge on \(\mathbf{\Sigma}\), for example that entries far from the diagonal are 0. And again, sometimes we do know \(\mathbf{\Sigma}\) quite accurately a priori. For example, [60] study a genetic dataset with \(m=377749\) predictor variables (single nucleotide polymorphisms, SNPs). There they use existing knowledge on the correlation structure of such variables. For example, SNPs on different chromosomes are roughly independent of each other (depending on the type of population the sample is from). Further, scientists have some knowledge of the correlations between the SNPs within a chromosome. Alternative methods for sampling approximate knockoffs have also been proposed [62], [63].

To conclude, knockoff methods do not always work well (i.e., they do not always lead to accurate FDR control), but that is mainly because conditional independence testing with many variables is often nearly impossible. Knockoff methods can potentially be a good approach in this context, although research is ongoing. The discussed methods are implemented in the R package knockoff, available on CRAN [64].

9.5 Exercises↩︎

1. Suppose \(m=6\), \(\alpha=0.2\) and \((W_1,...,W_6)=(0.7, 1.1 ,-0.8 ,0.5 ,1.2,0.6)\).

   (a) Which hypotheses are rejected by the knockoffs method as defined in Theorem 19? Which hypotheses are rejected by the knockoffs+ method?

   (b) For each of these two methods, briefly discuss what it guarantees in this situation.

9.6 Solutions↩︎

1. (a) We first consider the knockoffs method. For \(t=0.5\), we have \[\frac{|\{j: W_j\leq -t\}|}{|\{j: W_j\geq t\}|} = 1/5 \leq \alpha.\] For smaller \(t\) the fraction becomes \(\infty\). Hence, \(\tau=0.5\), so all hypotheses except the third one are rejected.

   We now consider the knockoffs+ method. For \(t\in [0.5,0.8]\), we have \[\frac{1+ |\{j: W_j\leq -t\}|}{|\{j: W_j\geq t\}|} \geq 2/5 > \alpha.\] For smaller \(t\) the fraction becomes \(\infty\). For \(t>0.8\), we have \[\frac{1+ |\{j: W_j\leq -t\}|}{|\{j: W_j\geq t\}|} = \frac{1}{|\{j: W_j\geq t\}|} \geq 1/2 > \alpha.\] Hence, \(\tau+=\infty\) and the knockoffs+ method rejects no hypotheses.

   (b) The knockoffs method guarantees that

   \[\mathbb{E}(\frac{|\mathcal{R}\cap\mathcal{N}|}{|\mathcal{R}|+1/\alpha}) =\mathbb{E}(\frac{|\mathcal{R}\cap\mathcal{N}|}{|\mathcal{R}|+5}) \leq 0.2.\] This does not have any interpretation that is similar to FDR control or one of the other error rates that we have seen. Since \(m\) is small and \(\alpha\) is not so small, the method does not even provide something close to FDR control. The knockoffs+ method is stricter, but guarantees that the \(FDR\), i.e. \(\mathbb{E}(\frac{|\mathcal{R}\cap\mathcal{N}|}{|\mathcal{R}|})\), is at most 0.2.

10 Measure theoretic probability↩︎

An extensive treatment of measure theoretic probability and martingale theory is contained in the freely available lecture notes of [65]. In §10.1 and §10.2 we discuss a few concepts from this theory, because theory from §10.2 is used in a few proofs in §8 and §9. In §10.1 we provide a formal definition of random variables, which is required in §10.2. §10.2 contains theory on martingales and stopping times.

10.1 Measure theory and random variables↩︎

In probability theory (and in particular martingale theory), the concept of a measurable random variable can be important. To define this concept, we must first introduce the notion of a \(\sigma\)-algebra [65].

Let \(\mathcal{S}\) be a non-empty set. As usual, let \(2^\mathcal{S}\) denote the set of all subsets of \(\mathcal{S}\). A collection \(\Sigma\subseteq 2^{\mathcal{S}}\) is called a \(\sigma\)-algebra on \(\mathcal{S}\) if

1. \(\mathcal{S}\in \Sigma,\)

2. \(\mathcal{E}\in \Sigma \Rightarrow \mathcal{E}^c\in \Sigma,\)

3. \(\mathcal{E}, \mathcal{F}\in \Sigma \Rightarrow \mathcal{E}\cup \mathcal{F}\in \Sigma\),

4. if \(\mathcal{E}_n\in \Sigma\) (\(n=1,2,...\)) then \(\bigcup_{i=1}^{\infty} \mathcal{E}_n \in \Sigma\).

Note that \(\emptyset\in \Sigma\) and if \(\mathcal{E}, \mathcal{F}\in \Sigma\), then \(\mathcal{E} \cap \mathcal{F}\in \Sigma\), since \(\mathcal{E} \cap \mathcal{F} = (\mathcal{E}^c \cup \mathcal{F}^c)^c.\)

If \(\Sigma\) is a \(\sigma\)-algebra on \(\mathcal{S}\), then \((\mathcal{S},\Sigma)\) is called a measurable space and the elements of \(\Sigma\) are called measurable sets.

If \(\mathcal{D}\) is a collection of subsets of \(\mathcal{S}\), then by \(\sigma(\mathcal{D})\) we denote the smallest \(\sigma\)-algebra containing \(\mathcal{D}\). If \(\Sigma = \sigma(\mathcal{D})\), we say that \(\mathcal{D}\) generates \(\Sigma\).

If \(\mathcal{S}=\mathbb{R}\) is the set of real numbers, then the \(\sigma\)-algebra on \(\mathbb{R}\) that is usually considered is the Borel \(\sigma\)-algebra, which is the \(\sigma\)-algebra generated by all open subsets of \(\mathbb{R}\). This is the same as the \(\sigma\)-algebra generated by all open intervals \((a,b)\). This is also identical to the \(\sigma\)-algebra generated by all closed intervals \([a,b]\) for example. This is because an open interval is a union of closed intervals. For example, \((0,5)\) is the union of all intervals of the form \([0+1/n,5-1/n]\), where \(n\in\{1,2,...\}\).

A map \(\mu: \Sigma \rightarrow [0,\infty]\) is called a measure if it is countable addititive, which means that

• \(\mu(\emptyset)=0,\)

• \(\mu(\cup_i\mathcal{E}_i) = \sum_i \mu(\mathcal{E}_i)\) for every countable collection \(\mathcal{E}_1, \mathcal{E}_2,...\) of disjoint elements of \(\Sigma\).

Further, if \(\mu(\mathcal{S})=1\), then \(\mu\) is called a probability measure. Such a measure is often denoted by \(\mathbb{P}\) and \((\mathcal{S}, \Sigma, \mathbb{P})\) is then called a probability space. It is then also common to use the symbol \(\Omega\) for \(\mathcal{S}\). We say that a subset of \(\mathcal{S}\) is measurable if it is contained in \(\Sigma\).

We will also need the notion of a measurable function. Let \(\mathcal{B}\) denote the Borel \(\sigma\)-algebra. A function \(X: \mathcal{S} \rightarrow \mathbb{R}\) is called (\(\Sigma\)-)measurable if for every \(\mathcal{A}\in \mathcal{B}\), \(X^{-1}(\mathcal{A})\in \Sigma\). \(X\) being (\(\Sigma\)-)measurable means that if we know for each element of \(\Sigma\) whether it contains \(\omega\), then we know \(X(\omega)\).

A random variable is a measurable function from a probability space \(\Omega\) to \(\mathbb{R}\) [65]. More precisely, it is a function \(X: \Omega \rightarrow \mathbb{R}\) such that for every \(\mathcal{A}\in \mathcal{B}\), \(X^{-1}( \mathcal{A})\in \Sigma\). We could also interpret the term “random variable” more broadly as a measurable map from a probability space to some other measurable space, e.g. \(\mathbb{R}^n\) with the Borel \(\sigma\)-algebra on \(\mathbb{R}^n\), or a space of functions, etc.. In case of \(\mathbb{R}^n\), we often speak of a random vector. In general, the space where \(X\) takes its values is called the sample space. The interpretation of the function \(X\) is as follows: “the probability that the random variable \(X\) takes a value in \(\mathcal{A}\)” is \(\mathbb{P}(X^{-1}( \mathcal{A}))\). This quantity is well-defined since we know that \(X^{-1}( \mathcal{A}) \in \Sigma\), since \(X\) is measurable. Note that when we write \(\mathbb{P}(X\in \mathcal{A})\), this is interpreted as \(\mathbb{P}(\omega\in \Omega: X(\omega)\in \mathcal{A})\), which is the probability that \(X\) is in \(\mathcal{A}\).

The \(\sigma\)-algebra generated by a real-valued variable \(X\) is the set \[\sigma(X):=\sigma(X^{-1}(\mathcal{A}): \mathcal{A}\in \mathcal{B}),\] i.e., the \(\sigma\)-algebra generated by \(\{X^{-1}(\mathcal{A}): \mathcal{A}\in \mathcal{B}\}\). \(\sigma(X)\) is the smallest \(\sigma\)-algebra with respect to which \(X\) is measurable. Conditioning on \(\sigma(X)\) means conditioning on knowing \(X\). Indeed, if we know which elements of \(\sigma(X)\) contain \(\omega\), then we know \(X(\omega)\), and vice versa.

The set \[\{X\in \mathcal{A}\}=\{\omega\in \Omega: X(\omega)\in \mathcal{A}\} = \{X^{-1}( \mathcal{A})\}\] is called the event that \(X\) is in \(A\). Formally, an event is a measurable subset of the measurable space \((\Omega, \Sigma)\). We say that an event \(\mathcal{A}\) holds almost surely (“a.s.”) if \(\mathbb{P}(\mathcal{A})=1,\) i.e., there may be some \(\omega\)’s outside \(\mathcal{A}\), but their measure \(\mathbb{P}(\Omega\setminus\mathcal{A})\) is \(0\).

In practice we often do not explicitly define random variables or vectors as functions. Instead, we often just define a probability measure on the sample space — e.g. by defining a joint PDF on the sample space \(\mathbb{R}^n\), which implies a probability measure on \(\mathbb{R}^n\).

The probability measure \(\mathbb{P}\) of the probability space \((\Omega, \Sigma, \mathbb{P})\) and the variable \(X\) imply a probability measure \(\mu_X\) on the outcome space (with \(\sigma\)-algebra \(\mathcal{B}\), say), if for every \(\mathcal{A}\in \mathcal{B}\) we define \(\mu_X(\mathcal{A})=\mathbb{P}(X^{-1}(\mathcal{A}))\). Note that \(\mu_X(\mathcal{A})\) then has the interpretation of being the probability that the random variable \(X\) takes a value in \(\mathcal{A}\).

The expected value of a random variable \(X\) is defined as the (Lebesgue) integral of \(X(\cdot)\) with respect to \(\mathbb{P}\), \[\mathbb{E}(X) =\int_{\Omega} X(\omega) d\mathbb{P}(\omega),\] see [65] for details on Lebesgue integration, which is a generalization of Riemann integration. In many basic texts on probability theory, the expectation of a continuous random variable is defined in terms of the PDF of the variable, and the expectation of a discrete random variable is defined as a sum. Such definitions have limitations. For example, a random variable might be neither fully continuous nor fully discrete. The above definition is more general.

10.2 Martingales and optional stopping↩︎

A stochastic process (in \(\mathbb{R}\)) is a sequence of random variables (measurable functions from a measurable space \((\Omega, \mathcal{F})\) to \((\mathbb{R}, \mathcal{B})\)) indexed by a set \(T\) that indicates time. That set could e.g. be \(T=\{0,1,2,...\}\) (a discrete set of time points), \(T=[0,\infty)\) or \(T=[0,1]\) (a continuous set of time points). We denote the process by \((X_t)_{t\in T}\). Just like a random variable, the process \((X_t)_{t\in T}\) is a mapping from \(\Omega\). A realization of the process corresponds to some \(\omega\in \Omega\), so it can be written as \((X_t)_{t\in T}(\omega)\).

In order to define martingales and to give the optional stopping theorem, we need the concept of a filtration \(\mathbb{F}\), which is a sequence \((\mathcal{F}_t)_{t\in T}\) of sub-\(\sigma\)-algebras on \(\mathcal{F}\), such that for all \(s,t\in T\) with \(s\leq t\) we have \(\mathcal{F}_s\subseteq \mathcal{F}_t\subseteq \mathcal{F}\).

A stochastic process \((X_t)_{t\in T}\) is called adapted to \(\mathbb{F}\) if for every \(t\in T\), \(X_t\) is measurable with respect to \(\mathcal{F}_t\). Note that \(X_t\) being measurable with respect to \(\mathcal{F}_t\) means that if we know for every element of \(\mathcal{F}_t\) whether it contains \(\omega\), then we know \(X(\omega)\).

Conditioning on a \(\sigma\)-algebra has the interpretation of knowing for every element of the \(\sigma\)-algebra whether it contains \(\omega\). The interpretation of a filtration \(\mathbb{F}\) is often that \(\mathcal{F}_t\) represents the information available by time \(t\), in the following sense: if we have observed the process \((X_t)_{t\in T}(\omega)\) up to time \(s<\max(T)\), then we have observed part of the process but not the whole process, so we only have limited information about which \(\omega\in \Omega\) our process corresponds to. However, we do know for each element of \(\mathcal{F}_s\) whether it contains \(\omega\).

The interpretation is as follows: conditioning on \(\mathcal{F}_s\) means that we know in which subsets of \(\mathcal{F}_s\) our realization \(\omega\) is contained, which means that we know \((X_t)_{t\leq s}\). The martingale property says that given our information at time \(s\), the expected future values of the process are exactly \(X_s\). A supermartingale is defined in the same way, except that the equality in 37 becomes “\(\leq\)”.

For example, suppose we are playing a gambling game with multiple rounds, and in every round we win a euro or lose a euro, each with probability 0.5. For simplicity, suppose that it is possible to have negative capital, so that we can always keep playing. Then if we are currently in timepoint \(s\) and have 40 euros, then conditional on everything that has happened until now, our expected capital at any future timepoint is also 40 euros.

Often we are interested in evaluating or stopping a process at some timepoint that we pick based on everything that has happened before it. For example, we might want to stop the gambling game as soon as our capital exceeds 1000. This relates to the concept of stopping times.

Thus, \(\tau: \Omega \rightarrow T\cup\{\infty\}\) is a stopping time if the decision to stop at time \(t\) depends on the information available at time \(t\).

Then we have the following result, which is known a version of the optional stopping theorem.

Thus, in the mentioned gambling game (where we can have negative capital), if we start with 100 euros, and we stop once we have won 3 times in a row, then our expected capital at the stopping time is 100 euros. (Note that if we stop when we have 200 euros, then Theorem 20 does not apply since \(\tau\) is not bounded then. Note that if we would use that strategy, with probability 1, we would at some point reach 200 euros. However, it turns out that the expected time until we reach the value 200 is infinite. This strategy is not useful in practice, because we do not get expected gains after finite time.)

For super-martingales the above optional stopping theorem also holds, except that the equality becomes “\(\leq\)”. In general, the theorem is called the optional stopping theorem because the stopping decision is based on observed information, i.e., you do not need to stop at some fixed time.

References↩︎

---

1. If we compute the two-sided p-value for the hypothesis that \(\beta_1=0\) in a simple Gaussian regression model (using lm() in R), then we will obtain the same p-value as above. This is related to the fact that if the \((X_i,Y_i)\) are multivariate normal and \(\rho=0\), then conditional on \(X_1,...,X_n\), the slope \(\beta_1\) in the regression model is 0. Note however that the traditional regression model does not assume that the \((X_i,Y_i)\) are multivariate normal — in fact, it does not require the \(X_i\) to be stochastic. A consequence is that for the traditional regression model, \(\beta_1\neq0\) does not necessarily mean that \(\rho\neq0\).↩︎

2. Admissibility is only defined given the assumptions: if we make additional assumptions, then we can potentially uniformly improve a method, even if we could not under the original assumptions.↩︎

3. An equivalent formulation is that the CTP rejects \(\mathcal{H}_{\mathcal{I}}\) if and only if for all \(I\subseteq \mathcal{J}\in \mathcal{C}\), \(\phi_{\mathcal{J}}=1.\)↩︎

4. Recall that a hypothesis is a set of distributions, so \(\mathcal{H}^a\subseteq \mathcal{H}^b\) means that if the true distribution is in \(\mathcal{H}^a\), it is also in \(\mathcal{H}^b\). In other words, if \(\mathcal{H}^a\) is true, then \(\mathcal{H}^b\) is true, i.e., \(\mathcal{H}^a\) implies \(\mathcal{H}^b\).↩︎

5. The reason why our \(\gamma\) is called alpha in that package, is that the literature sometimes uses \(\alpha\) to denote the “target FDP”.↩︎

6. Here we use that \(\mathcal{H}_i\) is rejected by BH if and only if there is a \(1\leq j\leq m\) with \(P_i\leq P_j\) and \(P_j\leq \frac{R(P_j)}{m}\alpha\).↩︎

7. This means in particular that conditional on \((X_j:j\in \mathcal{N}^c)\), \(Y\) is independent of each individual \(X_j\) with \(j\in \mathcal{N}\). However, it does not always mean that conditional on \((X_j:j\in \mathcal{N}^c)\), \(Y\) is independent of the whole vector \((X_j:j\in \mathcal{N})\). [60, pp. 557–558] give an example to illustrate that. However, if we exclude some pathological cases, these two notions are equivalent.↩︎

8. In fact it is sufficient to consider the case \(c=0.5\), for proving FDR control of knockoffs+. Then \((1-c)/c\) is simply \(1\).↩︎

9. For example, for \(i\neq j\), \(\mathop{\mathrm{Cov}}(X_j,X_i)=\Sigma_{ji}=(\Sigma-\text{diag}(s))_{ji} = \mathop{\mathrm{Cov}}(\tilde{X}_j,X_i)\).↩︎

10. The well known result that we use is the following. If a random vector \(Z\) has a \(N(\mu_z,\mathbf{G}_z)\) distribution with \[Z=\begin{bmatrix} {Z_1}\\ {Z_2} \end{bmatrix},\quad \mu_z=\begin{bmatrix} {\mu_1}\\ {\mu_2} \end{bmatrix},\quad \mathbf{G}_z= \begin{bmatrix} \mathbf{G}_{11}& \mathbf{G}_{12} \\ \mathbf{G}_{21} &\mathbf{G}_{22} \end{bmatrix} ,\] then the conditional distribution of \(Z_2\) given \(Z_1=z_1\) is \[N\big(\mu_2+\mathbf{G}_{21}\mathbf{G}_{11}^{-1}(Z_1-\mu_1), \mathbf{G}_{22}- \mathbf{G}_{21}\mathbf{G}_{11}^{-1}\mathbf{G}_{12}\big).\]↩︎