New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
A general framework for multiple testing dependence

Communicated by Burton H. Singer, Princeton University, Princeton, NJ, September 4, 2008 (received for review May 8, 2008)
Abstract
We develop a general framework for performing largescale significance testing in the presence of arbitrarily strong dependence. We derive a lowdimensional set of random vectors, called a dependence kernel, that fully captures the dependence structure in an observed highdimensional dataset. This result shows a surprising reversal of the “curse of dimensionality” in the highdimensional hypothesis testing setting. We show theoretically that conditioning on a dependence kernel is sufficient to render statistical tests independent regardless of the level of dependence in the observed data. This framework for multiple testing dependence has implications in a variety of common multiple testing problems, such as in gene expression studies, brain imaging, and spatial epidemiology.
In many areas of science, there has been a rapid increase in the amount of data collected in any given study. This increase is due in part to the ability to computationally handle large datasets and the introduction of various highthroughput technologies. Analyzing data from such highdimensional studies is often carried out by performing simultaneous hypothesis tests for some behavior of interest, on each of thousands or more measured variables. Largescale multiple testing has been applied in fields such as genomics (1–3), astrophysics (4, 5), brain imaging (6–8), and spatial epidemiology (9). By their very definition, highdimensional studies rarely involve the analysis of independent variables, rather, many related variables are analyzed simultaneously. However, most statistical methods for performing multiple testing rely on independence, or some form of weak dependence, among the data corresponding to the variables being tested. Ignoring the dependence among hypothesis tests can result in both highly variable significance measures and bias caused by the confounding of dependent noise and the signal of interest.
Here, we develop an approach for addressing arbitrarily strong multiple testing dependence at the level of the original data collected in a highdimensional study, before test statistics or P values have been calculated. We derive a lowdimensional set of random vectors that fully captures multiple testing dependence in any fixed dataset. By including this lowdimensional set of vectors in the modelfitting process, one may remove arbitrarily strong dependence resulting in independent parameter estimates, test statistics, and P values. This result represents a surprising reversal of the “curse of dimensionality” (10), because of the relatively small sample size in relation to the large number of tests being performed. Essentially, we show that the manifestation of the dependence cannot be too complex and must exist in a lowdimensional subspace of the data, driven by the sample size rather than by the number of hypothesis tests. This approach provides a sharp contrast to currently available approaches to this problem, such as the estimation of a problematically large covariance matrix, the conservative adjustment of P values, or the empirical warping of the test statistics' null distribution.
The main contributions of this article can be summarized as follows. We provide a precise definition of multiple testing dependence in terms of the original data, rather than in terms of P values or test statistics. We also state and prove a theoretical result showing how to account for arbitrarily strong dependence among multiple tests; no assumptions about a restricted dependence structure are required. By exploiting the dimensionality of the problem, we are able to account for dependence on each specific dataset, rather than relying on a populationlevel solution. We introduce a model that, when fit, makes the tests independent for all subsequent inference steps. Utilizing our framework allows all existing multiple testing procedures requiring independence to be extended so that they now provide strong control in the presence of general dependence. Our general characterization of multiple testing dependence directly shows that latent structure in highdimensional datasets, such as population genetic substructure (11) or expression heterogeneity (12), is a special case of multiple testing dependence. We propose and demonstrate an estimation technique for implementing our framework in practice, which is applicable to a large class of problems considered here.
Notation and Assumptions
We assume that m related hypothesis tests are simultaneously performed, each based on an n−vector of data sampled from a common probability space on ℝ^{n}. The data corresponding to hypothesis test i are x_{i} = (x_{i1}, x_{i2}, …, x_{in}), for i = 1, 2, …, m. The overall data can be arranged into an m × n matrix X where the ith row is composed of x_{i}. We assume that there are “primary variables” Y = (y_{1},…,y_{n}) collected, describing the study design or experimental outcomes of interest, and any other covariates that will be employed. Primary variables are those that are both measured and included in the model used to test the hypotheses.
We assume that the goal is to perform a hypothesis test on E[x_{i}Y]. We will also assume that E[x_{i}Y] can be modeled with a standard basisfunction model, which would include linear models, nonparametric smoothers, longitudinal models, and others. To this end, we write E[x_{i}Y] = b_{i}S(Y), where b_{i} is a 1 × d−vector and S(Y) is a d × n matrix of basis functions evaluated at Y d < n. When there is no ambiguity, we will write S = S(Y) to simplify notation. Note that Y can be composed of variables such as time, a treatment, experimental conditions, and demographic variables. The basis S can be arbitrarily flexible to incorporate most of the models commonly used in statistics for continuous data.
The residuals of the model are then e_{i} = x_{i} − E[x_{i}Y] = x_{i} − b_{i}S. Analogously, we let E be the m × n matrix, where the ith row is e_{i}. We make no assumptions about what distribution the residuals follow, although by construction E[e_{i}S(Y)] = 0. We allow for arbitrary dependence across the tests, i.e., dependence across the rows of E. We assume that the marginal model for each e_{i} is known or approximated sufficiently when performing the hypothesis tests. That is, we assume that the marginal null model for each test is correctly specified.
In matrix form, the model can be written as The goal is then to test m hypotheses of the form: where the null and alternative hypothesis tests are identically defined for each of the tests. This setup encompasses what is typically employed in practice, such as in gene expression studies and other applications of microarrays, brain imaging, spatial epidemiology, astrophysics, and environmental modeling (4, 7, 9, 13, 14).
Two Open Problems
The classical approach to testing multiple hypotheses is to first perform each test individually. This involves calculating a 1dimensional statistic for each test, usually as some comparison of the model fit under the constraint of the null hypothesis to that under no constraints. By utilizing the observed test statistics and their null distributions, we calculate a P value for each test (15). An algorithm or point estimate is then applied to the set of P values to determine a significance threshold that controls a specific error measure at a specific level (16), such as the false discovery rate (FDR) at 10% (17, 18). Variations on this approach have been suggested, such as estimating a qvalue for each test (19) or a posterior error probability (20). Regardless of the approach, the validity and accuracy of these procedures are essentially determined by whether the null distributions are correctly specified (or conservatively specified) and whether the data are independent (or weakly dependent) across tests (21, 22).
Two open problems in multiple testing have received a lot of recent attention. The first is concerned with controlling multiple testing error measures, such as the FDR, in the presence of dependence among the P values (23, 24). This dependence is usually formulated as being present in the “noise” component of the models used to obtain the P values. The second open problem is concerned with the fact that latent structure among the tests can distort what would usually be the correct null distribution of the test statistics (11,25–27). The approach proposed here shows that both problems actually stem from sources of variation that are common among tests, which we show is multiple testing dependence, and both problems can be simultaneously resolved through one framework.
The current paradigm for addressing these two problems can be seen in Fig. 1, where the steps taken to get from the original data X to a set of significant tests are shown. It can be seen that existing approaches are applied far downstream in the process. Specifically, adjustments are performed after 1dimensional summaries of each test have been formed, either to the test statistics or P values. As we show below, the information about noise dependence and latent structure is found in the original data X by modeling common sources of variation among tests. Our proposed approach addresses multiple testing dependence (from either noise dependence or latent structure) at the early modelfitting stage (Fig. 1), at which point the tests have been made stochastically independent and the null distribution is no longer distorted.
Proposed Framework
Definition of Multiple Testing Dependence.
Multiple testing dependence has typically been defined in terms of P values or test statistics resulting from multiple tests (21, 24, 26, 28, 29). Here, we form populationlevel and estimationlevel definitions that apply directly to the full dataset, X. The estimationlevel definition also explicitly involves the model assumption and fit utilized in the significance analysis. When fitting model 1, we denote the estimate of B by
Definition: We say that populationlevel multiple testing dependence exists when it is the case that: We say that estimationlevel multiple testing dependence exists when it is the case that:
Multiple testing dependence at the population level is therefore any probabilistic dependence among the x_{i}, after conditioning on Y. In terms of model 1, this is equivalent to the existence of dependence across the rows of E; i.e., dependence among the e_{1}, e_{2}, …, e_{m}. Estimationlevel dependence is equivalent to dependence among the rows of the residual matrix
A General Decomposition of Dependence.
Dependence among the rows of E and among the rows of
A key feature is that, in the multiple testing scenario, the dimension along which the sampling occurs is different than the dimension along which the multivariate inference occurs. In terms of our notation, the sampling occurs with respect to the columns of X, whereas the multiple tests occur across the rows of X. This samplingtoinference structure requires one to develop a specialized approach to multivariate dependence that is different from the classical scenarios. For example, the classical construction and interpretation of a P value threshold is such that a true null test is called significant with P value ≤ α at a rate of α over many independent replications of the study. However, in the multiple testing scenario, the P values that we utilize are not P values corresponding to a single hypothesis test over m independent replications of the study. Rather, the P values result from m related variables that have all been observed in a single study from a single sample of size n. The “sampling variation” that forms the backbone of most statistical thinking is different in our case: we observe one instance of sampling variation among the variables being tested. Therefore, even if each hypothesis test's P value behaves as expected over repeated studies, the set of P values from multiple tests in a single study will not necessarily exhibit the same behavior. Whereas this phenomenon prevents us from invoking wellestablished statistical principles, such as the classical interpretation of a P value, the fact that we have measured thousands of related variables from this single instance of sampling variation allows us to capture and model the common sources of variation across all tests. Multiple testing dependence is variation that is common among hypothesis tests.
Thus, rather than proposing a populationlevel approach to this problem (which includes the population of all hypothetical studies that could take place in terms of sampling of the columns of X), we directly model the random manifestation of dependence in the observed data from a given study, by aggregating the common sampling variation across all tests' data. Including this information in the model during subsequent significance analyses removes the dependence within the study. Therefore dependence is removed across all studies, providing studyspecific and populationlevel solutions. To directly model the random manifestation of dependence in the observed data, we do the following: (i) additively partition E into dependent and independent components, (ii) take the singular value decomposition of the dependent component, and (iii) treat the right singular values as covariates in the model fitting and subsequent hypothesis testing. To this end, we provide the following result, which shows that any dependence can be additively decomposed into a dependent component and an independent component. It is important to note that this is both for an arbitrary distribution for E and an arbitrary (up to degeneracy) level of dependence across the rows of E.
Proposition 1. Let the data corresponding to multiple hypothesis tests be modeled according to Eq. 1. Suppose that for each e_{i}, there is no Borel measurable function g such that e_{i} = g(e_{1}, …, e_{i−1}, e_{i+1}, …, e_{m}) almost surely. Then, there exist matrices Γ_{m×r}, G_{r×n} (r ≤ n), and U_{m×n} such that where the rows of U are jointly independent random vectors so that Also, for all i = 1, 2, …, m, u_{i} ≠ 0 and u_{i} = h_{i}(e_{i}) for a nonrandom Borel measurable function h_{i}.
A formal proof of Proposition 1 and all subsequent theoretical results can be found in the supporting information (SI) Appendix. Note that if we let r = n and then set U = 0 or set U equal to an arbitrary m × n matrix of independently distributed random variables, then the independence of the rows of U is trivially satisfied. However, our added assumption regarding e_{i} allows us to show that a nontrivial U exists where u_{i} ≠ 0 and u_{i} = h_{i}(e_{i}) for a deterministic function h_{i}. In other words, u_{i} is a function of e_{i} in a nondegenerate fashion, which means that U truly represents a rowindependent component of E. The intuition behind these properties is that our assumption guarantees that e_{i} does indeed contain some variation that is independent from the other tests. For hypothesis tests where there does exist a Borel measurable g such that e_{i} = g(e_{1}, …, e_{i−1}, e_{i+1}, …, e_{m}), then the variation of e_{i} is completely dependent with that of the other tests' data. In this case, one can set u_{i} = 0 and the above decomposition is still meaningful.
The decomposition of Proposition 1 immediately indicates one direction to take in solving the multiple testing dependence problem, namely to account for the ΓG component, thereby removing dependence. To this end, we now define a “dependence kernel” for the data X.
Definition: An r × n matrix G forms a dependence kernel for the highdimensional data X, if the following equality holds: where the rows of U are jointly independent as in Proposition 1.
In practice, one would be interested in minimal dependence kernels, which are those satisfying the above definition and having the smallest number of rows, r. Proposition 1 shows that at least one such G exists with r ≤ n rows. As we discuss below in Scientific Applications, the manner in which one incorporates additional information beyond the original observations to estimate and utilize Γ and G is context specific. In the SI Appendix, we provide explicit descriptions for two scientific applications, latent structure as encountered in genomics and spatial dependence as encountered brain imaging. We propose a new algorithm for estimating G in the genomics application and demonstrate that it has favorable operating characteristics.
Dependence Kernel Accounts for Dependence.
An important question arises from Proposition 1. Is including G, in addition to S(Y), in the model used to perform the hypothesis tests sufficient to remove the dependence from the tests? If this is the case, then only an r × n matrix must be known to fully capture the dependence. This is in contrast to the m(m−1)/2 parameters that must be known for a covariance matrix among tests, for example. To put this into context, consider a microarray experiment with 1,000 genes and 20 arrays. In this case, the covariance has ∼500,000 unknown parameters, whereas G has, at most, 400 unknown values. The following two results show that including G in addition to S(Y) in the modeling is sufficient to remove all multiple hypothesis testing dependence.
Corollary 1. Under the assumptions of Proposition 1, all populationlevel multiple testing dependence is removed when conditioning on both Y and a dependence kernel G. That is,
If instead of fitting model 1, suppose that we instead fit the decomposition from Proposition 1, where we assume that S and G are known: It follows that estimationlevel multiple testing independence may then be achieved.
Proposition 2. Assume the data for multiple tests follow model 1, and let G be any valid dependence kernel. Suppose that model 3 is fit by least squares, resulting in residuals
Since G will be unknown in practice, the practical implication of this proposition is that we have to estimate only the relatively small r × n matrix G well in order to account for all of the dependence, while the simple leastsquares solution to Γ suffices. When the row space jointly spanned by S and G has dimension equal to n, then the above proposition becomes trivially true. However, if we assume that S, G, and Γ are known, then the analogous estimationlevel independence holds. In this case, we have to estimate Γ and G well in order to account for dependence. These (m + r)n parameters are still far smaller than the unknown m(m  1)/2 parameters of a covariance matrix, for example.
Strong Control of Multiple Testing Error Rates.
Many methods exist for strongly controlling the familywise error rate (FWER) or FDR (16, 18, 19, 21, 24, 32). These methods are applied to the P values calculated from multiple hypothesis tests. Most of these methods require the P values corresponding to true null hypotheses to be independent in order for the procedure to provide strong control. For example, finitesample strong control of several FDR procedures (21, 24) and the conservative point estimation of the FDR (19) all require the true null P values to be independent. Several methods exist for controlling FWER or FDR when dependence is present. However, these either tend to be quite conservative or require special restrictions on the dependence structure (21, 24).
When utilizing model 3, the statistics formed for testing the hypothesis should be based on a function of the model fits and residuals. When this is the case, we achieve the desired independence of P values.
Corollary 2. Suppose that the assumptions of Proposition 2 hold, model 3 is utilized to perform multiple hypothesis tests, and G is a known dependence kernel. If P values are calculated from test statistics based on a function of the model fits and residuals, then the resulting P values and test statistics are independent across tests.
In other words, Corollary 2 extends all existing multiple testing procedures that have been shown to provide strong control when the null P values are independent to the general dependence case. Instead of deriving new multiple testing procedures for dependence at the level of P values, we can use the existing ones by including G into the model fitting and inference carried out to get the P values themselves.
Scientific Applications
Two causes for multiple testing dependence can be directly derived from scientific problems of interest. In each case, the dependence kernel G has a practical scientific intepretation.
Spatial Dependence.
Spatial dependence usually arises as dependence in the noise because of a structural relationship among the tests. In this case, we will consider the e_{i} of model 1 to simply represent “noise,” an example being the spatial dependence for noise that is typically assumed for brainimaging data (6–8). In this setting, the activity levels of thousands of points in the brain are simultaneously measured, where the goal is to identify regions of the brain that are active. A common model for the measured intensities is a Gaussian random field (6). It is assumed that the Gaussian noise among neighboring points in the brain are dependent, where the covariance between two points in the brain is usually a function of their distance.
In Fig. 2 A and B, we show two datasets generated from a simplified 2dimensional version of this model. It can be seen that the manifestation of dependence changes notably between the two studies, even though they come from the same data generating distribution. Using model 3 for each dataset, we removed the ΓG term. In both cases, the noise among points in the 2dimensional space becomes independent and the P value distributions of points corresponding to true null hypotheses follow the Uniform distribution. It has been shown that null P values following the Uniform(0,1) distribution is the property that confirms that the assumed null distribution is correct (22). Additionally, it can be seen that the null P values from the unadjusted data fluctuate substantially between the two studies, and neither follows the Uniform(0,1) null distribution. This is due to varying levels of correlation between S and G from model 3. In one case, S and G are correlated producing spurious signal among the true null hypotheses; this would lead to a major inflation of significance. In the other case, they are uncorrelated leading to a major loss of power. By accounting for the ΓG term, we have resolved these issues.
Latent Structure.
A second source of multiple testing dependence is a latent structure due to relevant factors not being included in the model. It is possible for there to be unmodeled factors that are common among the multiple tests but that are not included in S. Suppose there exists unmodeled factors Z such that E(x_{i}Y) ≠ E(x_{i}Y, Z) for more than one test. If we utilize model 1 when performing the significance analysis, there will be dependence across the rows of E induced by the common factor Z, causing populationlevel multiple testing dependence. Likewise, there will be dependence across the rows of R causing estimationlevel multiple testing dependence. A similar case can arise when the model for x_{i} in terms of Y is incorrect. For example, it could be the case that E[x_{i}Y] = b_{i}S*(Y), where the differences between S and S* are nontrivial among multiple tests. Here, there will be dependence across the rows of R induced by the variation common to multiple tests due to S* but not captured by S, which would cause estimationlevel multiple testing dependence. Failing to include all relevant factors is a common issue in genomics leading to latent structure (11, 12). The adverse effects of latent structure due to unmodeled factors on differential expression significance analyses has only recently been recognized (12).
Fig. 2 C and D shows independently simulated microarray studies in this scenario, where we have simulated a treatment effect plus effects from several unmodeled variables. The unmodeled factors were simulated as being independently distributed with respect to the treatment, which is equivalent to a study in which the treatment is randomized. As in Fig. 2 A and B, it can be seen that the P values corresponding to true null hypotheses (i.e., genes not differentially expressed with respect to the treatment) are not Uniformly distributed. When utilizing model 3 for these data and subtracting the term ΓG, the residuals are now made independent and the null P values are Uniform(0,1) distributed.
Estimating G in Practice
There are a number of scenarios where estimating G is feasible in practice. One scenario is when nothing is known about the dependence structure, but it is also the case that d + r < n, where d and r are the number of rows of the model S and dependence kernel G, respectively. This is likely when the dependence is driven by latent variables, such as in gene expression heterogeneity (12). In the SI Appendix, we present an algorithm for estimating G in this scenario. It is shown that the proposed algorithm, called iteratively reweighted surrogate variable analysis (IRWSVA), exhibits favorable operating characteristics. We provide evidence for this over a broad range of simulations. Another scenario is when the dependence structure is well characterized at the population level. Here, it may even be the case that d + r ≈ n. This scenario is common in brain imaging (6, 7) and other spatial dependence problems (9), as discussed above. The fact that Γ is largely determined by the known spatial structure allows us to overcome the fact that d + r ≈ n (SI Appendix).
Discussion
We have described a general framework for multiple testing dependence in highdimensional studies. Our framework defines multiple testing dependence as stochastic dependence among tests that remains when conditioning on the model is used in the significance analysis. We presented an approach for addressing the problem of multiple testing dependence based on estimating the dependence kernel, a lowdimensional set of vectors that completely defines the dependence in any highthroughput dataset. We have shown that if the dependence kernel is known and included in the model, then the hypothesis tests can be made stochastically independent. This work extends existing results regarding error rate control under independence to the case of general dependence. An additional advantage of our approach is that we can not only estimate dependence at the level of the data, which is intuitively more appealing than estimating dependence at the level of P values or test statistics, but we can also directly adjust for that dependence in each specific study. We presented an algorithm with favorable operating characteristics for estimating the dependence kernel for one of the main two scientific areas of interest that we discussed. We anticipate that well behaved estimates of the dependence kernel in other scientific areas are feasible.
One important implication of this work is that multiple testing dependence is tractable at the level of the original data. Downstream approaches to dealing with multiple testing dependence are not able to directly capture general dependence structure (Fig. 1). Another implication of this work is that, for a fixed complexity, the stronger the dependence is among tests, the more feasible it is to properly estimate and model it. It has also been shown that the weaker multiple testing dependence is, the more appropriate it is to utilize methods that are designed for the independence case (21). Therefore, there is promise that the full range of multiple testing dependence levels is tractable for a large class of relevant scientific problems.
Acknowledgments
We thank the editor and several anonymous referees for helpful comments. This research was supported in part by National Institutes of Health Grant R01 HG002913.
Footnotes
 ^{1}To whom correspondence should be addressed. Email: jstorey{at}princeton.edu

Author contributions: J.T.L. and J.D.S. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/cgi/content/full/0808709105/DCSupplemental.

Freely available online through the PNAS open access option.
 © 2008 by the National Academy of Sciences of the USA
References
 ↵
 Storey JD,
 Tibshirani R
 ↵
 ↵
 Tusher VG,
 Tibshirani R,
 Chu G
 ↵
 ↵
 Starck JL,
 Pires S,
 Refregier A
 ↵
 ↵
 ↵
 Schwartzman A,
 Dougherty RF,
 Taylor J
 ↵
 Elliott P,
 Wakefield J,
 Best N,
 Briggs D
 ↵
 Bellman R
 ↵
 ↵
 ↵
 Green PJ,
 Silverman BW
 ↵
 Worsley KJ
 ↵
 Rice JA
 ↵
 ↵
 ↵
 Benjamini Y,
 Hochberg Y
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Klebanov L,
 Jordan C,
 Yakovlev
 ↵
 Qui X,
 Klebanov L,
 Yakovlev AY
 ↵
 Mardia KV,
 Kent JT,
 Bibby JM
 ↵
 Owen A
 ↵
 Holm S