{'fraction_non_alphanumeric': 0.041099541223277684, 'fraction_numerical': 0.014410538539148312, 'mean_word_length': 4.406751749691231, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Distance-based analysis of variance: approximate inference and an application to genome-wide association studies

In several modern applications, ranging from genetics to genomics and neuroimaging, there is a need to compare observations across different populations, such as groups of healthy and diseased individuals. The interest is in detecting a group effect. When the observations are vectorial, real-valued and follow a multivariate Normal distribution, multivariate analysis of variance (MANOVA) tests are routinely applied. However, such traditional procedures are not suitable when dealing with more complex data structures such as functional (e.g. curves) or graph-structured (e.g. trees and networks) objects, where the required distributional assumptions may be violated. In this paper we discuss a distance-based MANOVA-like approach, the DBF test, for detecting differences between groups for a wider range of data types. The test statistic, analogously to other distance-based statistics, only relies on a suitably chosen distance measure that captures the pairwise dissimilarity among all available samples. An approximate null probability distribution of the DBF statistic is proposed thus allowing inferences to be drawn without the need for costly permutation procedures. Through extensive simulations we provide evidence that the proposed methodology works well for a range of data types and distances, and generalizes the traditional MANOVA tests. We also report on an application of the proposed methodology for the analysis of a multi-locus genome-wide association study of Alzheimer's disease, which has been carried out using several genetic distance measures.


Introduction
Multivariate analysis of variance (MANOVA) techniques are routinely applied to compare realvalued multivariate observations collected in two or more populations and test whether the mean vectors in those groups have been drawn from the same underlying sampling distribution (see Anderson (1984) and Krzanowski (2000), for example). These methods are applicable in a wide range of scientific areas. For instance, in genomics, microarrays allow the researcher to investigate the behaviour of thousands of genes simultaneously under various conditions. A common task consists of observing the expression levels of a group of genes sharing the same biological function, and assessing whether the mean expression levels in these gene sets differ across experimental conditions. Compared to a univariate approach involving only one gene expression measurement at a time, a MANOVA test may capture potential gene interactions and hence provide a more powerful and biologically meaningful way of detecting subsets of differentially expressed genes (Szabo et al., 2003;Tsai and Chen, 2009;Xu and Cui, 2008;Shen et al., 2011).
In several applications, however, MANOVA tests may be inappropriate for at least two main reasons. Firstly, when the observations are multivariate, the multivariate normality assumption may not necessarily hold, e.g., if the observations are heavily skewed or are discrete-valued. Secondly, the observations may not necessarily be represented by vectorial data structures. In an increasingly large number of applications, the data being compared across groups are functional (e.g. curves) and graph-structured (e.g. trees and networks) objects. This is indeed the case with several applications we have a particular interest in, of which the following are three representative examples: • In genome-wide association (GWA) studies, the allele frequencies observed at multiple genetic markers (single nucleotide polymorphisms, or SNPs) in healthy and diseased subjects are compared in order to identify genetic variants that are associated with disease risk. Although GWAs have emerged as popular methodologies for gene mapping, traditional analyses that rely on testing one individual marker at a time are believed to have low power to detect small genetic effects. They have also been shown to have limited reproducibility, and are unable to detect effects which are due to gene-gene interactions. For these reasons, sets of multiple markers (i.e. markers observed at multiple loci) can be used instead, giving rise to discrete-valued vectors to be compared across groups (Wu et al., 2010).
• In longitudinal microarray experiments, repeated measurements of a gene's expression level are recorded at different time-points, and the resulting time course is often modelled as a functional object or curve, i.e. an infinite-dimensional data structure (Storey et al., 2005) (see Figure 1 for an example of such a functional dataset). The resulting gene curves observed under two or more experimental conditions are then compared using a test statistic that accounts for the functional nature of the data (Berk et al., 2012;Berk and Montana, 2009).
• In neuroimaging studies, there is increasing interest in comparing individual brain connectivity networks inferred in subjects that belong to two or more clinically relevant groups, e.g. diseased and healthy individuals. A network is a collection of nodes (vertices) and links (edges) between pairs of nodes. Such brain connectivity networks describe the set of connections in the neural system, or connectome, in which nodes could be neurons or cortical areas, and edges could be axons or fibre tracts. Thus, edges could refer to the structural connectivity of a neural network, or could signify correlations between the activity patterns of nodes forming functional connectivity (Rubinov and Sporns, 2010). Although univariate test statistics can be used to compare particular summary features extracted from the graphs, such as centrality measures, there is a need to compare individual graphs directly.
Since standard MANOVA techniques assume real-valued vectorial data representations to compare mean vectors, they are not always suitable. Moreover, in many of these applications, a very large number of tests are needed to be simultaneously performed within the same experiment.
For instance, in genomic experiments such as those involving gene sets or time courses, as well as GWA studies, the number of comparisons can range from many tens of thousands to millions.
Computationally efficient inferential methods are therefore required.
In this article we consider the problem of comparing two or more independent random samples, representative of some underlying populations, by using a distance-based analysis of variance approach. As in traditional MANOVA, the total variability observed in the random samples is decomposed into the sum of between-group and within-group variability. The decomposition only relies on a suitable distance measure defining how dissimilar any two samples are, and does not require the computation of means in each group; both attractive properties when dealing with non-vectorially structured objects. For a given application, a suitable distance measure is defined depending on the nature of the data (whether vectorial or not, for example), and on the specific objectives of the study, since each distance measure captures a different feature of the data. We discuss particular distances for genetic and functional data later in this article. Pairwise distances between all available samples are computed and stored in a square, symmetric distance matrix, with samples deemed more similar if the distances between them are small.
A number of distance-based tests for no group effect have been proposed in the literature as generalizations of classical MANOVA procedures, and are applicable when the MANOVA assumptions are violated, such as with data of the type described above. These include the Mantel test (Mantel, 1967), the MRPP test (Mielke and Berry, 2007), and the more recent DBF test (Minas et al., 2011). Each statistic utilizes some notion of within-and between-group variability in terms of distances: Mantel considers between-group variability, MRPP considers within-group variability, and DBF considers the ratio of between-to within-group variability. An important, and in many cases limiting, feature of such distance-based tests is that exact distribution theory of the statistics under the null hypothesis of no group effect is unavailable. This is due to the wide variety of distance measures available for different data types leading to different distributional characteristics of the distances (Mantel, 1967). Due to this limitation, non-parametric and computationally intensive statistical procedures based on permutations are typically applied, since they do not rely on any distributional assumptions. The only requirement is that samples are exchangeable across groups under the null hypothesis of no group effect (Pesarin and Salmaso, 2010). Statistical significance of an observed statistic is assessed by comparing the estimate against the discrete sampling distribution obtained by permuting the samples across groups and recomputing the statistic each time. For each permutation, this is equivalent to permuting the rows and columns of the corresponding distance matrix simultaneously, and recomputing the statistic with the permuted distance matrix (Mielke and Berry, 2007;Minas et al., 2011). Thus the distance matrix needs only be computed once, and inferences can be drawn based on this matrix.
Despite being routinely adopted, such a permutation approach suffers from some important limitations. Even for low sample sizes, exact p-values require a very large number of permutations to be obtained; for instance, for 12 observations, the total number of permutations required is O(10 8 ). Exact p-values are often impossible to obtain due to computational and time constraints, and Monte Carlo procedures are used instead whereby a smaller number of randomly chosen permutations are used to approximate the p-values. This approach inevitably introduces sampling errors (Berry and Mielke, 1983). Moreover, whereas large p-values can be well approximated by a Monte Carlo approach, smaller ones will be estimated less accurately (Mielke and Berry, 2007;Knijnenburg et al., 2009). For example, in order to obtain a permutation p-value within 10 −5 of the true p-value, it has been shown that O(10 7 ) permutations are required. In several real applications, we have observed that less than O(10 5 ) permutations are typically used, especially when a very large number of statistical tests must be performed simultaneously. In these cases, the total number of permutations per test may be vastly limited even when ad-hoc parallel implementations that run on high-performance computing facilities are being used.
In order to deal with these limitations, in this work we propose a closed-form, continuous approximation to the null sampling distribution of the DBF statistic (Minas et al., 2011), which was originally proposed for the analysis of functional data arising from longitudinal microarray experiments. Our main contribution here is to enable distance-based analysis of variance testing in a wide range of applications, in a computationally efficient manner without the need for permutations. Using extensive simulations, we demonstrate that the proposed methodology works well for several different data structures, including functional and discrete-valued vectorial data, in addition to several distances, and that it provides a good approximation of the permutation p-values. We also show that traditional ANOVA and MANOVA tests are special cases of the proposed approach when their assumptions are satisfied. For these special cases we show that the proposed continuous distribution approximates the true distributions well. The computational cost of analyzing large datasets is thus reduced significantly without compromising the accuracy of the results.
A second contribution of this work consists of applying of the proposed methodology to a multi-locus GWA study of the Alzheimer's Disease Neuroimaging Initiative (ADNI) cohort (see for instance, Vounou et al. (2010Vounou et al. ( , 2012). To the best of our knowledge, this is the first multi-locus case-control genetic study on this cohort. By using a range of genetic distances we show that the DBF test identifies previously known genetic variants associated with Alzheimer's disease, and discuss a comparison with a distance-based method that has been specifically designed for the analysis of such large-scale case-control studies.
This article is organized as follows. In Section 2, we briefly review the classical variance decomposition and MANOVA tests, followed by a detailed account of the proposed distance-based variance decomposition. The corresponding DBF statistic is defined, and theoretical connections with classical MANOVA and ANOVA tests are included. Section 3 describes how the sampling distribution of the DBF statistic exhibits skewed characteristics, and we justify the use of the Pearson type III distribution for its approximation. Approximate distribution theory for the DBF statistic then follows in Section 4. We provide evidence that the approximation works well for a range of distances and data types in Section 5, including empirical comparisons with ANOVA and MANOVA. Finally, in Section 6 we showcase the generality and applicability of the proposed approach by analyzing the ADNI dataset and reporting on the genetic variants we have identified.
Conclusive remarks are given in Section 7.

A distance-based test of equality between populations 2.1 Traditional variance decomposition and MANOVA tests
We first consider the case of N independent vectorial and real-valued observations {y i } N i=1 in R P that belong to one of G populations with means {µ g } G g=1 . Each group is of size N g such that N = G g=1 N g . We are interested in testing the null hypothesis that all the population means are equal, versus the alternative hypothesis that some of the means are not equal.
Several standard statistical tests are based on the multivariate analysis of variance (MANOVA).
The overall sample mean is given byȳ = 1 N N i=1 y i and each within-group sample mean is given byȳ g = 1 Ng N i=1 y i I gi for g = 1, . . . , G where I gi is an indicator variable taking the value 1 if observation y i is in group g and 0 otherwise. The P × P total sum of squares matrix is given by and can be partitioned into the sum of between-and within-group sum of squares matrices, respectively. Existing MANOVA test statistics make use of different elements of this total sum of squares decomposition. For G = 2, the well-known Hotelling's T 2 statistic is given by N 1 N 2 (N − 2)(ȳ 1 − y 2 ) T W −1 (ȳ 1 −ȳ 2 )/N . For G > 2, MANOVA statistics include Wilks' Lamda, det(W )/det(T ), the Pillai trace, tr T −1 B , and the Lawley-Hotelling trace, tr W −1 B (Rencher, 2002;Krzanowski, 2000). Under the assumption that the observations are independent and identically distributed from a multivariate Normal distribution with mean µ and variance-covariance matrix Σ, some distributional results are available. For instance, for G = 2, Hotelling's T 2 statistic (multiplied by a constant depending on N and P ) has an exact F distribution under the null; where F P,N −P −1 denotes the F distribution with degrees of freedom P and N − P − 1. For G > 2, Wilks' Lamda, the Pillai trace and the Lawley-Hotelling trace can all be similarly transformed to statistics which are well-approximated by the F distribution with degrees of freedom dependent on N , G and P (see, for example, Rencher (2002)). When N < P , as is typically the case with genomic datasets, the classical MANOVA tests cannot be applied directly. This is because the T and W matrices are singular, and so have at least one zero-valued eigenvalue and cannot be inverted. Several high-dimensional MANOVA settings have been considered in the literature, and tests of equality between groups have been proposed, some using traditional MANOVA statistics with generalized inverses (see Srivastava (2007) and Schott (2007) for good reviews, and Tsai and Chen (2009) for an application to gene expression data). One of the first tests was proposed by Dempster for G = 2 (Dempster, 1960), where an F-type statistic defined as the ratio of within-to between-group variability was proposed.
This statistic was generalized for G > 2 and named the Dempster trace criterion four decades later by Fujikoshi et al. (2004). They noticed that the trace operator could be applied to the B and W sum of squares matrices to yield equivalent expressions to those proposed by Dempster. Although not stated explicitly, they summarize the total variability exhibited by the P variables by tr(T ). They then partition this into between-and within-group components via tr(T ) = tr(B) + tr(W ), which follows by applying the trace operator to the decomposition T = B + W . The Dempster trace criterion was then defined as tr(B)/tr(W ), and a transformation of this statistic was shown to be asymptotically Gaussian. In the next section we describe a further elaboration of this approach which only relies on pairwise distances among samples.

The distance-based variance decomposition
As before, we consider N independent observations {y i } N i=1 belonging to G groups each of size N g for g = 1, . . . , G. Now, however, we place no restriction on the nature of these observations; they can be of any form, for example, scalar-valued, vector-valued, functional, graph-structured, or images, amongst others. The fundamental assumption is that we are able to define a distance measure d, which may be either semimetric or metric, which quantifies the dissimilarity between any pair of observations in the random sample. All pairwise distances are then arranged in an The choice of distance depends on the type of data and scientific problem at hand. We are then interested in testing the null hypothesis that there is no difference between the observations from different groups with respect to the chosen distance measure d, against the alternative that a difference exists between any two groups. That is, under the null, all observations are assumed to belong to the same population.
First, we introduce a generalization of the variance decomposition approach using pairwise distances. Consider the quantity tr(T ) associated with a set of vectorial observations. It is a measure of spread found by summing the squared Euclidean distance of each observation to the population mean vector. This quantity can be equivalently written using only pairwise Euclidean distances between observations, as follows: where d E denotes the Euclidean distance. Thus, tr(T ) is proportional to the sum of squared interpoint Euclidean distances between all N observations. This well-known connection shows that the total variability of a given set of vectorial observations, traditionally found using the population mean, can be computed using only the inter-point Euclidean distances (Gower and Krzanowski, 1999;Anderson, 2001;Minas et al., 2011). In an analogous manner, the within-and between-group variability quantities tr(W ) and tr(B) can also be written in terms of squared Euclidean distances, as follows: and since tr(B) = tr(T ) − tr(W ), we obtain Generalizations of these quantities can be defined by replacing the Euclidean distance, d E , with any distance d. Thus we can define the total variability of a set of observations {y i } N i=1 with respect to distance d as the within-group variability with respect to d as and the between-group variability with respect to d as The total variability in the data captured by T ∆ can hence be decomposed into the sum of two components quantifying within-and between-group variability. That is, analogously to the decomposition tr(T ) = tr(W ) + tr(B). We can write T ∆ , B ∆ and W ∆ more compactly in matrix form by introducing the N × N Gower's centered inner product matrix (Gower, 1966), defined as where I N is the identity matrix of size N , J N is the square matrix of ones of size N and This matrix contains all the information on the inter-point distances between the N observations, and is such that its trace equals T ∆ ; Therefore we rewrite T ∆ more conveniently as tr(G ∆ ). For W ∆ and B ∆ we define the centered N × N matrix of constants encoding group membership of each observation to one of the G groups as where J a is the square matrix of ones of size a. Since this matrix is centered, we have that and we use this fact in the evaluation of the quantity tr(H c G ∆ ) to derive expressions for W ∆ and

The DBF test statistic
Making use of the distance-based variance decomposition above, a distance-based test statistic has been defined to test the null hypothesis of equality between groups with respect to the chosen distance measure. This statistic, denoted DBF, is of the form and was originally used to compare functional data arising in genomic experiments (Minas et al., 2011). Analogously to the Dempster trace criterion and Lawley-Hotelling trace statistic, this statistic considers a ratio of between-to within-group variability. Larger values of this statistic provide evidence against the null hypothesis, as larger between-group variability and smaller within-group variability suggest that observations in the same group are more similar than observations in different groups. A statistic of similar form was also proposed by Anderson (2001) for application in ecology, but with degrees of freedom divisors G − 1 and N − G in the numerator and denominator, respectively. When the observations are P -dimensional vectors, it has been shown that F ∆ is related monotonically to the Lawley-Hotelling and Pillai trace MANOVA statistics for G > 2 (Minas et al., 2011). For instance, on defining the distance matrices it was shown that Lawley-Hotelling = P F ∆W and Pillai trace = P F ∆T 1 + (1 − P )F ∆T .
For G = 2, it was also shown that F ∆ is monotonically related to Hotelling's T 2 statistic via the equation Expanding on this relationship between T 2 and F ∆T , since we know that (N −P −1)T 2 /((N −2)P ) has an exact F distribution under the null, it follows that That is, a transformation of F ∆T follows the exact F distribution with degrees of freedom P and N − P − 1. As a special case, when the data consists of scalar observations and the Euclidean distance d E , is applied, F ∆E is identical to the classical one-way ANOVA F statistic, ignoring the degrees of freedom divisors G − 1 and N − G in the numerator and denominator, respectively. Thus, Anderson, 2001;Minas et al., 2011).
Given an observed value of the test statistic,F ∆ , computed for any suitably chosen distance measure d, inference can be carried out using a non-parametric approach. That is, the p-value can be found using permutations. Given N π permutations π ∈ Π, where π is a one-to-one map-ping of the set {1, . . . , N } to itself, the set {F ∆π } π∈Π is generated by recalculating B ∆ for each permutation, denotedB ∆π , and using the monotonic relationship T ∆ is fixed for all permutations so that permutated values of F ∆ are monotonically related to permuted values of B ∆ . The p-value is then computed as the proportion of the N π permuted statistics greater than or equal to the observedF ∆ , i.e., Clearly, this is a one-sided test, since only larger values of F ∆ provide evidence against the null.
As an alternative to this expensive permutation-based testing approach, we propose an approximate distribution for the null sampling distribution of F ∆ in Section 3. From this approximate distribution p-values can be well-approximated without the need for permutations. Since F ∆ is related to B ∆ via (3), we first approximate the null distribution of B ∆ .
3 Distance-based between-group variability

Skewed characteristics
For general data structures and distance measures, the null sampling distribution of DBF test statistic (1) is unknown. The reason for this is that the between-group variability quantity, B ∆ , featuring in the statistic will, in general, follow some unknown distribution which depends on the specific distance measure being used (Mantel, 1967). On denoting the (i, j) th element of H c by h ij and recalling that H c is centered, B ∆ can be expressed as the weighted sum of squared distances Thus even if each d 2 (y i , y j ), for i = j, was assumed to be a random variable with known distribution, B ∆ would be a weighted sum of correlated and uncorrelated random variables, whose distribution would be difficult to evaluate. For instance, the problem of evaluating the sum of correlated and uncorrelated Chi-squared and Gamma random variables has been considered extensively (see, for example, Solomon and Stephens (1977), Kourouklis and Moschopoulos (1985)). Although it has been argued that B ∆ has the appearance of a U-statistic which is asymptotically normal (Mantel, 1967;Hoeffding, 1948), in our experience with different vectorial and non-vectorial data types arising in real applications, even for large samples sizes, B ∆ often appears to be skewed to various degrees. Further insights into this problem can be obtained by exploring the empirical permutation distribution of B ∆ for a number of real datasets involving different data structures and distances.
As an illustration, we consider three different examples: (i) Vectorial and real-valued data: the data consists of the P = 50 gene expression measurements observed on N = 103 biological samples from the Novartis multi-tissue dataset described in Monti et al. (2003). In this case, G = 4, corresponding to four different tissues. For this dataset, we considered the Euclidean, Mahalanobis and Manhattan distances.
(ii) Vectorial and discrete-valued data: the data consists of P = 5 randomly selected SNPs observed on N = 254 samples from the ADNI dataset (see Section 6 for further details). (iii) Functional data (curves): the data consists of N = 18 gene expression functional data replicates for a randomly selected gene in a dataset on M.tuberculosis analyzed by Tailleux et al. (2008). In this case, G = 2, corresponding to two different types of cell, and replicate time courses were observed at 4 time-points for each type of cell. These were smoothed via cubic smoothing splines to yield the 18 replicate curves (Minas et al., 2011). Figure 1 shows observed time courses and their fitted curves for two randomly selected genes. The L 2 , Visual L 2 and Curvature distance measures were applied to this dataset. The L 2 measure captures the difference in magnitude between curves, the Visual L 2 measure captures their scale-invariant differences in shape, and the Curvature measure captures their difference in rate of change regardless of direction. See Appendix for further details.  The exact permutation distribution of B ∆ in each case would be given by the set {B ∆π } π∈Π where Π contains all N ! permutations π of the elements of {1, . . . , N }. For large N , due to the computational effort required in enumerating all possible permutations, the exact distribution is generally unavailable. Figure 2 shows the approximate sampling distribution of B ∆ obtained by using 10 6 Monte Carlo permutations for each of the data types and distance measures considered. In almost all cases the distributions are heavily skewed. Remarkably, this is also the case for large sample sizes.

A Person Type III approximation
In order to approximate the exact permutation distribution B ∆ , we propose an approach based on moment matching. Given the skewness often observed in real datasets, we assume that B ∆ follows a Pearson type III distribution. This is a type of Gamma distribution encompassing other distributions as special cases, such as the Exponential and Normal distributions (Josse et al., 2008), and has been successfully used for modeling skewed sampling distributions of various other statistics (Berry and Mielke, 1983;Kazi-Aoual et al., 1995;Josse et al., 2008).
The first three moments of the exact permutation distribution of B ∆ are given by respectively. By evaluating the expressions analytically, closed form manipulations of the three moments have been derived allowing their efficient computation for N > 6 without the need for permutations (Kazi-Aoual et al., 1995). These closed form expressions require only the square, symmetric and centered matrices H c and G ∆ . The mean and variance, for example, are given by The expression for skewness is much more involved; the reader is referred to Kazi-Aoual et al. (1995).
On standardizing B ∆ by subtracting µ B and dividing by σ B , the Pearson type III distribution can be parameterized by only the skewness parameter, γ B . That is, where P T III denotes the Pearson type III distribution. By assumption of this model, the support for γ B > 0 and −2/γ B ≤ b < ∞, where Γ(·) denotes the usual Gamma function, for γ B < 0 and −∞ < b ≤ −2/γ B , and the standard Normal distribution for γ B = 0 (Mielke and Berry, 2007).
In practice γ B is not equal to zero exactly, as exhibited for some real datasets in Figure 2. We see that all sampling distributions, and corresponding approximate PDFs, are skewed to some degree. Thus in the next section we only consider the cases where γ B > 0 and γ B < 0, and do not consider the trivial case of γ B = 0.

Approximate null distribution of the DBF statistic
The DBF statistic F ∆ is related to the standardized distance-based between-group variability quantity B s ∆ via the one-to-one function h : We aim to derive the approximate null distribution of F ∆ in terms of the distribution of B s ∆ via transformation h, which is required to be continuous over the given supports of B s ∆ . However, it is not continuous in the positive plane at β = (T ∆ − µ B )/σ B because T ∆ = tr (G ∆ ) > µ B due to tr(H c ) = 1. Since the boundaries of both supports of B s ∆ depend on the skewness, the position of β may or may not affect the continuity of the distribution of B s ∆ over the given support. We thus consider dealing with the discontinuity separately for both supports of B s ∆ , i.e. for both cases of skewness.
First consider the positive skewness case, where the support of B s ∆ is [−2/γ B , ∞). Since γ B > 0, −2/γ B is negative and the discontinuity shows itself as B s ∆ increases form −2/γ B to ∞. Figure 3 (a) shows how F ∆ behaves as a function of B s ∆ over this support; F ∆ is an increasing function of B s ∆ on both sides of the discontinuity at β. We thus divide the support of B s ∆ into the two regions, namely where the equivalent supports of F ∆ are given by [α, ∞), where satisfies h (α) = −2/γ B , and (−∞, −1). We can show by contradiction that α > −1, since α ≤ −1 implies T ∆ γ B ≤ 0, and by definition both T ∆ and γ B are positive. In these regions we can apply the transformation since there are no discontinuities, and define the CDF of F ∆ in terms of the CDF of B s ∆ by for γ B > 0. The relevant derivations are in the Appendix, including a proof that this is a valid CDF. Now we turn our attention to the negative skewness case, where the support of B s ∆ is (−∞, −2/γ B ]. In this case, γ B < 0 so that −2/γ B is positive. The discontinuity at B s ∆ = β thus only needs to be considered if β is to the left of −2/γ B , otherwise it can be ignored since it is not included in the support of B s ∆ . We consider these two cases separately. First consider the case where β < −2/γ B . We have that Applying this with the equation for α given by (6) yields α < −1. Thus we define the occurrence of this first case when α < −1. Figure 3 (b) shows how F ∆ behaves as a function of B s ∆ over the support (−∞, −2/γ B ] when α < −1. As with the positive skewness case, F ∆ is an increasing function of B s ∆ on both sides of the discontinuity at β. We thus divide the support of B s ∆ into the two regions (−∞, β) and β, −2 γ B , with equivalent supports of F ∆ given by (−1, ∞) and (−∞, α], in which F ∆ and B s ∆ are monotonically related with no discontinuities. Thus we define the CDF of F ∆ as for γ B < 0 and α < −1.
Now consider the case where β > −2/γ B ; this is equivalent to α > −1. In this case there are no discontinuities in the support of B s ∆ , so F ∆ and B s ∆ are monotonically related everywhere. The support for B s ∆ given by (−∞, −2/γ B ] is equivalent to the support for F ∆ of (−1, α]. Thus the CDF of F ∆ is defined as for γ B < 0 and α > −1.
Using these results, the approximate p-value of an observedF ∆ can be readily obtained without permutations. On computing the permutational mean µ B , variance σ 2 B , and skewness γ B , and additionally α if γ B < 0, the p-value is given by 1 is the CDF chosen for the specific case of skewness and α value.
For the given case of skewness and α value, the PDF of where the range of f is given by the selected case of CDF.

Comparison with ANOVA and MANOVA
Given that F ∆ equals the ANOVA F statistic up to a constant as a special case for univariate data, we verify that the distribution of F ∆ approximates that of the ANOVA F statistic well as N and G increase. Also, since F ∆ is related to Hotelling's T 2 as a special case for multivariate data with G = 2, we verify that our proposed distribution, transformed via (2), approximates the distribution of T 2 well as N increases. That is, we aim to show that for the special cases the DBF test is approximately equivalent to the ANOVA F and Hotelling's T 2 tests, respectively.
For the univariate case, data was generated under the null and the DBF statistic using the Euclidean distance measure, d E , and the ANOVA F statistic were computed. P-values were found by comparing against their respective distributions. For N = 40, 100, 500, 1000 and G = 2, 4, 5, the k th Monte Carlo run consisted of simulating y 1 , . . . , y N ∼ N (µ k , σ 2 k ), where µ k ∼ U (−10, 10), σ 2 k ∼ U (0, 10), and U (a, b) denotes the Uniform distribution over [a, b]. The mean and standard deviation of the absolute differences between the p-values obtained for B = 200 Monte Carlo simulations are reported in Table 1. It can be seen that as N and G increase, the absolute difference between the p-values decreases, thus showing that the approximate distribution of the DBF statistic behaves as expected in this case. For the multivariate case, the DBF statistic using the Mahalanobis-like distance measure d T , and the Hotelling's T 2 statistic were computed. P-values were found by comparing against their respective distributions. For N = 40, 100, 500, 1000 and P = 10, the k th Monte Carlo run consisted of simulating y 1 , . . . y N ∼ N P (µ k , Σ k ), where µ k = (µ k1 , . . . , µ P k ) T with µ jk ∼ U (−6, 6) for j = 1, . . . , P , and Σ k a random Wishart matrix of size P × P . The mean and standard deviation of the absolute differences between the p-values obtained for B = 200 Monte Carlo runs are reported in Table 1. As N increases the difference between the p-values decreases, showing that the DBF and Hotelling's T 2 tests are approximately equivalent as N increases.
A further experiment was performed to show that, as N increases, the proposed approximate null distribution of the DBF statistic approximates the true ANOVA F and Hotelling's T 2 distributions, on applying the required transformations, quite well. In particular, we show that it yields a better approximation than a permutation-based CDF, especially when the number of permutations is low.
For P = 1, G = 2, and each of N = 50, 70, one set of univariate observations were generated under the null from a Normal distribution as above. The DBF null CDF, suitably transformed, and the ANOVA F CDF were obtained, and the Kolmogorov-Smirnov (KS) statistic was used to compute the difference between these distributions. This statistic is computed as the maximum distance between two vectors representing the CDFs of interest; we used a vector of 1000 equally spaced points across the range of the approximate DBF distribution. For the given dataset for  Figure 4. We see that for N = 50, using more than 5 × 10 4 permutations yields a permutation distribution which is directly comparable with our approximate distribution. For N = 70, however, the approximate DBF distribution better approximates the true underlying ANOVA F distribution than the permutation distributions typically used in practice; typically not more than 10 5 permutations are used for real data analyses.
For P = 10, G = 2 and N = 50, one set of multivariate observations were generated under the null from a Multivariate Normal distribution as described above. The DBF null CDF, suitably transformed, and the Hotelling's T 2 CDF were obtained. Repeating as above, and using the KS statistic to quantify the difference between the transformed DBF permutation CDF and true Hotelling's T 2 CDF, the results are given in plot (c) of Figure 4. We see that for N = 50 the approximate DBF distribution yields a better approximation of the true distribution than the permutation distributions.

Approximate null distribution for various data types and distances
In this section we illustrate how the approximate distribution of the DBF statistic compares with the Monte Carlo permutation distribution for a number of data types and distances. In our setting, we explore a range of sample sizes and distance measures for datasets simulated to mimic the real datasets introduced in Section 3.1. For vectorial and real-valued data, we consider the Euclidean, Bray-Curtis, Canberra, Manhattan and Maximum distances. For vectorial and discrete-valued genetic data, we consider the identity-by-state (IBS), Simple Matching, Sokal and Sneath, Rogers and Tanimoto I (RTI) and Hamman I distances. For functional data we consider the L 2 , Visual L 2 and Curvature distances (see Appendix for further details on the genetic and functional distance measures). On selection of data type, distance measure and number of samples N , the datasets are simulated as follows: (i) Vectorial and real-valued data: P -dimensional vectors y i = (y i1 , . . . , y iP ) T were simulated such that y ij ∼ N (0, 4) for i = 1, . . . , N and j = 1, . . . , 1000. For the Bray-Curtis and Canberra distances where positive values are required, we take absolute values.
(ii) Vectorial and discrete-valued data: 5-dimensional vectors y i = (y i1 , . . . , y i5 ) T for i = 1, . . . , N were simulated based on the observations of the 153 control subjects from chromosome 1 of the ADNI dataset introduced in Section 6. N control subjects were randomly selected and their minor allele counts at 5 randomly chosen SNPs across the chromosome selected.
(iii) Functional data (curves): N curves {y i (t)} N i=1 were simulated over the range t ∈ [0, 48] by using quadratic Bezier curves (Farin, 1992) and a sampling procedure involving smoothing splines (Ramsay and Silverman, 2005). N Bezier curves were randomly generated, and N 1000-dimensional vectors were sampled from these curves at equally spaced points across [0, 48] with standard Gaussian error. These were then smoothed via cubic smoothing splines to yield the N curves; see Minas et al. (2011) for further details of this curve-generation procedure. This procedure generates random curves similar to those observed in real longitudinal datasets, such as those shown in Figure 1.  i.e., no group effect. For N = 10, all N ! permutations were used to compute the permutation p-value, but for N = 30, 100, a Monte Carlo set of 10 6 permutations was used. The theoretical and permutation p-values were computed, and the mean and standard deviation of the absolute difference between these for each combination of data type, distance measure and N are reported in Table 2. As expected, the absolute difference between the p-values decreases as N increases for each distance measure applied to each data type.

A genome-wide association study of Alzheimer's disease
In traditional case-control association studies, subjects are genotyped for a comprehensive range of genetic markers across the entire genome, and markers associated with disease risk are sought. The objective of such studies is to identify genetic variants in the human genome that are associated with disease risk (see, for example, Altshuler et al. (2008) and Pearson and Manolio (2008) for good overviews of traditional GWA studies).
Genetic variations are captured by observing SNPs. Human SNPs are biallelic genetic markers, and as such, the genotype of an individual at a given SNP is represented by one of three combinations of two alleles; a major allele occurring more commonly in the study cohort, and a minor allele which is less common. The possible combinations of alleles are 'major, major', 'major, minor' and 'minor, minor', and the genotype of the individual at a given SNP is summarized by the minor allele count. The genotype of individual i at a given SNP k, denoted y ik , is thus represented by either a 0, 1 or 2 corresponding to homozygotes for the major allele, heterozygotes and homozygotes for the minor allele, respectively.
A common approach of scanning the genome in search of causal variants is to group SNPs together into subsets where it is plausible that some dependence exists between them, for example, if the SNPs are in the same gene or biological pathway; generally referred to as the multi-locus approach (Mukhopadhyay et al., 2010;Wu et al., 2010;Yang et al., 2009). Sliding windows which partition the genome into overlapping subsets of SNPs of the same length have also been proposed to group SNPs together, and ensure coverage of the entire genome without omitting possibly useful information in intergenic regions.
We describe an application of the DBF test with the approximate inference approach of Section 4 by conducting several multi-locus GWA studies on the Alzheimer Disease Neuroimaging Initiative (ADNI) cohort. The available dataset consists of 254 subjects, 101 cases of AD and 153 controls, all genotyped at 316, 348 SNPs. We assume the genome has been partitioned into subsets of P = 5 SNPs so that, for each subset, each of the N individuals in the given study cohort is represented by the discrete-valued P -dimensional vector {y i = (y i1 , . . . , y iP )} N i=1 . This results in a total number of 316, 260 SNP sets to be compared across the two populations. For each sliding window, five genetic distances (IBS, Simple Matching, Sokal and Sneath, Rogers and Tanimoto I and Hamman I) were applied, and the DBF statistic and corresponding p-value computed using the approach described in Section 4. Due to the multiple testing problem, a genome-wide significance threshold of 10 −7 was set to identify statistically significant SNP sets.
In Figure 5 we provide the corresponding Manhattan plot which depicts the significant SNP subsets across the entire genome for the Sokal and Sneath distance measure, showing the greatest effects around chromosomes 18 and 19. The results of all distance measures were summarized by the unique SNP and gene combinations identified; see Table 3. All significant SNPs are identified in chromosomes 18 and 19. In particular, chromosome 19 contains two genes, APOE4 and TOMM40, which are the major genetic variants found in many studies (see, for example, Braskie et al. (2011) and Shen et al. (2010)). Other previously reported genetic variants that overlap with our findings include APOC2, APOC4, PVRL2 and CLPTM1 (Takei et al., 2009;Yu et al., 2007). The DCC gene has also been previously identified (Bredesen, 2009;Lourenco et al., 2009).
As an illustrative comparison with the permutation approach, in Figure 6 we compare the approximate distribution of the DBF statistic under the null with the null permutation distribution obtained using by 10 6 Monte Carlo permutations. This was done for the first sliding window of chromosome 19 and three genetic distances. Here again it can be clearly seen that the approximate null distribution of the DBF statistic provides a good fit. We have also compared our results with those obtained using a statistical procedure developed specifically for multi-locus case-control studies, the kernel-machine logistic regression approach of Wu et al. (2010), which we denote LKMT. This method makes use of similarities rather than distances via kernel functions K(·, ·); note that kernels are related to distances, for example, the IBS kernel between individuals i and j is just one minus their IBS distance. Since LKMT makes use of an approximate distribution rather than permutations, we are able to compare the two methods without major computational issues. LKMT consists of a logistic kernel-machine regression model where the probability that a given subject is a case subject is modelled as some function of the similarities between that subject and all other subjects. We denote the case-control status of the N subjects by {s i } N i=1 , where s i = 0 for controls and s i = 1 for cases. The model is then given by logit[P (s i = 1)] = β 0 + m(y i ) for i = 1, . . . , N , where β 0 is an intercept, and m(z i ) = N j=1 δ j K(y i , y j ) for some constants {δ j } N j=1 so that m(·) is completely specified by the similarities between the subjects depicted by the kernel function K(·, ·). The null hypothesis of no group effect is stated as H 0 : h(y i ) = 0 for i = 1, . . . , N , and the variance-component score statistic The p-values obtained by LKMT are of the same order as our approach with the IBS distance measure when using the comparable IBS kernel. For instance, for chromosome 19, which is the chromosome containing the smallest p-values, Figure 7 provides a visual comparison of the two methods when using the IBS measure; similar results are obtained for other chromosomes (not shown). Although a rigorous power comparison is beyond the scope of this study, it can be noted that both methods identified the same causal SNPs in Chromosome 19 at the significance level of 10 −7 , i.e., the ones listed in Table 3. Thus the well-known APOE and TOMM40 genes were identified by both approaches.

Conclusion
For large datasets, such as those arising in biological applications, it is desirable to apply distancebased tests for equality between groups without permutations. This requires that the discrete sampling distribution of the chosen statistic under the null is approximated by a suitably chosen continuous distribution. In this article we focused on the DBF statistic due to the intuitive distance-based variance decomposition it uses. We described some results relating the DBF statistic to classical MANOVA and ANOVA statistics. Its connection with ANOVA is well-known, allowing an exact null distribution of the DBF statistic to be obtained when applied in the special case of univariate Normal data with the Euclidean distance measure. Its connection with MANOVA, however, is not so well-known. We presented the result that for multivariate Normal data with G = 2, the DBF statistic applied with a Mahalanobis-like distance measure is related to Hotelling's T 2 statistic. From this connection the exact null distribution of the transformed DBF statistic is attainable.
For more general data types, and corresponding distance measures, where the exact null distribution of the DBF statistic is unknown, we considered its permutation distribution. We showed that it depends on the permutation distribution of the between-group variability component of the distance-based variance decomposition. On presenting the skewed characteristics of the betweengroup variability for real biological datasets, we justified the use of the Pearson type III distribution to model its skewed nature. We then used its monotonic relationship with the DBF statistic to derive an approximate distribution for the DBF statistic. For a range of data types and distance measures we provided evidence to support the performance of the proposed approximation. We also used it to conduct several case-control GWA studies on the ADNI dataset in order to show its applicability to a large dataset. This study identified a number of genes inline with validated genetic risk factors for Alzheimer's disease. We also showed that the proposed approach performs comparably to an existing test statistic that has been specifically designed for case-control association studies. The R software for performing the DBF test with the approximate null distribution is available from http://www2.imperial.ac.uk/∼gmontana.

Appendix A Multi-locus genetic distance measures
Assume N individuals in the given study cohort are observed at P SNPs, yielding the discretevalued P -dimensional vectors {y i = (y i1 , . . . , y iP )} N i=1 . The identity-by-state (IBS) distance measure is commonly used, and measures the proportion of risk alleles (minor alleles) shared between individuals across the SNPs considered (Wessel and Schork, 2006;Wu et al., 2010;Mukhopadhyay et al., 2010). The distance between individuals i and j is defined by where s(y ik , y jk ) = 0 if y ik = 0 and y jk = 2, or if y ik = 2 and y jk = 0, s(y ik , y jk ) = 1 if y ik = 1 and y jk = 1, or if y jk = 1 and y ik = 1, and s(y ik , y jk ) = 2 if y ik = y jk . This distance takes values between 0 and 1. Weighted versions of this distance also exist where a weight is attached to each of the P SNPs depending on properties such as functional significance or frequency of the minor allele (Wessel and Schork, 2006;Li et al., 2009). Genetic distances have also been proposed based on the contingency table between individuals i and j containing the frequency that each combination of y ik and y jk values occur over the range of SNPs considered (Selinski and Ickstadt, 2005); see Table 4 below. The key statistics in this table are the number of complete matches of the minor alleles, m + = 2 k=0 m kk , and the number of mismatches, m − = P − m + , where the total number of comparisons made is P . Based on these quantities, the following 'matching coefficient' distance measures can be defined; the Simple Matching distance the Sokal and Sneath distance d SS (y i , y j ) = 1 − m + m + + 1 2 m − , the Rogers and Tanimoto I distance where s * (y i , y j ) = s HI (y i , y j ) + | min i,j {s HI (y i , y j )}| with s HI (y i , y j ) = (m + − m − )/P . All of these distances take values between 0 and 1.

B Distance measures between curves
Assume a set of N curves {y i (t)} N i=1 defined for t ∈ τ where all curves have been defined over the same range τ . The L 2 distance represents the area between curves, and hence the magnitude of difference between them (Ferraty and Vieu, 2006;Minas et al., 2011;Salem et al., 2010), and is defined by The curvature distance quantifies the difference in the rate of change between two curves (Ferraty and Vieu, 2006;Minas et al., 2011), and is defined by The visual L 2 distance quantifies the difference in the scale-invariant shape between curves, analogously to the difference detected by the human eye (Marron and Tsybakov, 1995;Minas et al., 2011). For this distance, curves y i and y j are scaled both in time and magnitude, so that their values range between 0 and 1 in time period [0, 1]; denote these y s i (t) and y s j (t) where t ∈ [0, 1]. They are then represented as infinite sets of points in the two-dimensional plane, denoted p i = {(t, y s i (t)) | t ∈ [0, 1]} and p j = (t, y s j (t)) | t ∈ [0, 1] . The visual L 2 distance is then defined by where δ ij (t) is the minimum Euclidean distance between the point y s i (t) and all points p j representing y s j ; note that δ ij (t) is not necessarily equal to δ ji (t).
C Derivation of the CDF of F ∆ for γ B > 0 First consider the case where −∞ < f < −1. By inspection of Figure 3 (a), we see that we need only consider the relationship between F ∆ and B s ∆ where B s ∆ > β. We thus have that Now consider the case where α ≤ f < ∞. From Figure 3 (a) we see that we must consider the relationship between F ∆ and B s ∆ for B s ∆ < β, while adding the cumulative component of all values of F ∆ for which B s ∆ > β. That is, and since −2/γ T ≤ B s ∆ < ∞, we have F B s ∆ (−2/γ B ; γ B ) = 0 and F B s ∆ (∞; γ B ) = 1, so that Thus we have that the cumulative distribution function of F ∆ is given by for γ B > 0, as required.
Next we show that this is a valid CDF by showing that the following conditions are satisfied: (i) The limit of F F∆ (f ) as f tends to −∞ from the right is 0, and as f tends to ∞ from the left is 1. That is (ii) F F∆ (f ) is a monotone, non-decreasing function of f . That is, for f 1 < f 2 , F F∆ (f 1 ) ≤ F F∆ (f 2 ).
For −∞ < f < −1 we have that and for α ≤ f < ∞ we have that as required.