Computing the skewness of the phylogenetic mean pairwise distance in linear time

Background The phylogenetic Mean Pairwise Distance (MPD) is one of the most popular measures for computing the phylogenetic distance between a given group of species. More specifically, for a phylogenetic tree and for a set of species R represented by a subset of the leaf nodes of , the MPD of R is equal to the average cost of all possible simple paths in that connect pairs of nodes in R. Among other phylogenetic measures, the MPD is used as a tool for deciding if the species of a given group R are closely related. To do this, it is important to compute not only the value of the MPD for this group but also the expectation, the variance, and the skewness of this metric. Although efficient algorithms have been developed for computing the expectation and the variance the MPD, there has been no approach so far for computing the skewness of this measure. Results In the present work we describe how to compute the skewness of the MPD on a tree optimally, in Θ(n) time; here n is the size of the tree . So far this is the first result that leads to an exact, let alone efficient, computation of the skewness for any popular phylogenetic distance measure. Moreover, we show how we can compute in Θ(n) time several interesting quantities in , that can be possibly used as building blocks for computing efficiently the skewness of other phylogenetic measures. Conclusions The optimal computation of the skewness of the MPD that is outlined in this work provides one more tool for studying the phylogenetic relatedness of species in large phylogenetic trees. Until now this has been infeasible, given that traditional techniques for computing the skewness are inefficient and based on inexact resampling.


Background
Communities of co-occuring species may be described as "clustered" if species in the community tend to be close phylogenetic relatives of one another, or "overdispersed" if they are distant relatives [1]. To define these terms we need a function that measures the phylogenetic relatedness of a set of species, and also a point of reference for how this function should behave in the absence of ecological and evolutionary processes. One such function is the mean pairwise distance (MPD); given a phylogenetic tree T and a subset of species R that are represented by leaf nodes of T , the MPD of the species in R is equal to average cost of all possible simple paths that connect pairs of nodes in R.
To decide if the value of the MPD for a specific set of species R is large or small, we need to know the average value (expectation) of the MPD for all sets of species in T that consist of exactly r = |R| species. To judge how much larger or smaller is this value from the average, we also need to know the standard deviation of the MPD for all possible sets of r species in T . Putting all these values together, we get the following index that expresses how clustered are the species in R [1]: where MPD(T , R) is the value of the MPD for R in T , and expec(T ) and sd MPD (T , r) are the expected value and the standard deviation respectively of the MPD calculated over all subsets of r species in T . In a previous paper we presented optimal algorithms for computing the expectation and the standard deviation of the MPD of a phylogenetic tree T in O(n) time, where n http://www.almob.org/content/9/1/15 is the number of the edges of T [2]. This enabled exact computations of these statistical moments of the MPD on large trees, which were previously infeasible using traditional slow and inexact resampling techniques. However, an important problem remained unsolved; quantifying our degree of confidence that the NRI value observed in a community reflects non-random ecological and evolutionary processes.
This degree of confidence can be expressed as a statistical P value, that is the probability that we would observe an NRI value as extreme or more so if the community was randomly assembled. Traditionally, estimating P is accomplished by ranking the observed MPD against the distribution of randomized MPD values [3]. If the MPD falls far enough into one of the tails of the distribution (generally below the 2.5 percentile or above the 97.5 percentile, yielding P < 0.05), the community is said to be significantly overdispersed or significantly clustered. However, this approach relies on sampling a large number of random subsets of species in T , and recomputing the MPD for each random subset. Therefore, this method is slow and imprecise. This problem is exacerbated when it is necessary to consider multiple trees at once, arising for example from a Bayesian posterior sample of trees [4,5]. In such cases, sufficient resampling from all trees in the sample can be computationally limiting.
We can approximate the P value of an observed NRI by assuming a particular distribution of the possible MPD values and evaluating its cumulative distribution function at the observed MPD. Because the NRI measures the difference between the observed values and expectation in units of standard deviations, this yields a very simple rule if we assume that possible MPD values are normally distributed: any NRI value larger than 1.96 or smaller than −1.96 is significant. Unfortunately, the distribution of MPD values is often skewed, such that this simple rule will lead to incorrect P value estimates [6,7]. Of particular concern, this skewness introduces a bias towards detecting either significant clustering or significant overdispersion [8]. If the distribution of MPD values for a particular tree can be reasonably approximated using a skew-normal distribution, calculating the skewness analytically would enable us to remove this bias and improve the accuracy of P value estimates. In the last part of the paper, we describe experiments on large randomly generated trees, supporting this argument. Further, when a large sample of trees should be considered, the full distribution of MPD values can be considered as a mixture of skew-normal distributions [9,10], greatly simplifying and speeding up the process of calculating P values across the entire set of trees.
However, so far there has been no result in the related literature that shows how to compute the needed skewness measure efficiently. Hence, given a phylogenetic tree T and an integer r there is the need to design an efficient and exact algorithm that can compute the skewness of the MPD for r species in T . This would provide the last critical piece required for the adoption of a fully analytical and efficient approach for analysing ecological communities using the MPD and the NRI.

Our results
In the present work we show how we can compute efficiently the skewness of the MPD. More specifically, given a tree T that consists of n edges and a positive integer r, we prove that we can compute the skewness of of the MPD over all subsets of r leaf nodes in T optimally, in (n) time. For the calculation of this skewness value we consider that every subset of exactly r species in T is picked uniform with probability out of all possible subsets that have r species. The main contribution of this paper is a constructive proof that leads straightforwardly to an algorithm that computes the skewness of the MPD in (n) time. This is clearly optimal, and it outperforms even the best algorithms that are known so far for computing lower-order statistics for other phylogenetic measures; for example the best known algorithm for computing the variance of the popular Phylogenetic Diversity (PD) runs in O(n 2 ) time [2].
More than that, we prove how we can compute in (n) time several quantities that are related with groups of paths in the given tree; these quantities can be possibly used as building blocks for computing efficiently the skewness (and other statistical moments) of phylogenetic measures that are similar to the MPD. Such an example is the measure which is the equivalent of the MPD for computing the distance between two subsets of species in T [11].
The rest of this paper is, almost in its entirety, an elaborate proof for computing the skewness of the MPD on a tree T in (n) time. In the next section we define the problem that we want to tackle, and we present a group of quantities that we use as building blocks for computing the skewness of the MPD. We prove that all of these quantities can be computed in linear time with respect to the size of the input tree. Then, we provide the main proof of this paper; there we show how we can express the value of the skewness of the MPD in terms of the quantities that we introduced earlier. The proof implies a straightforward linear time algorithm for the computation of the skewness as well. In the last section we provide experimental results that indicate that computing the skewness of the MPD can be a useful tool for improving the estimation of P values when a skew-normal distibution is assumed. There we describe experiments that we conducted on large randomly generated trees to compare two different methods for estimating P values; one method is based on random sampling of a large number of tip sets, and the http://www.almob.org/content/9/1/15 other method relies in calculating the mean, variance, and skewness of the MPD for the given tree.

Description of the problem and basic concepts
Definitions and notation Let T be a phylogenetic tree, and let E be the set of its edges. We denote the number of the edges in T by n, that is n = |E|. For an edge e ∈ E, we use w e to indicate the weight of this edge. We use S to denote the set of the leaf nodes of T . We call these nodes the tips of the tree, and we use s to denote the number of these nodes.
Since a phylogenetic tree is a rooted tree, for any edge e ∈ E we distinguish the two nodes adjacent to e into a parent node and a child node; among these two, the parent node of e is the one for which the simple path from this node to the root does not contain e. We use Ch(e) to indicate the set of edges whose parent node is the child node of e, which of course implies that e / ∈ Ch(e). We indicate the edge whose child node is the parent node of e by parent(e). For any edge e ∈ E, tree T (e) is the subtree of T whose root is the child node of edge e. We denote the set of tips that appear in T (e) as S(e), and we denote the number of these tips by s(e).
Given any edge e ∈ E, we partition the edges of T into three subsets. The first subset consists of all the edges that appear in the subtree of e. We denote this set by Off(e). The second subset consists of all edges e ∈ E for which e appears in the subtree of e . We use Anc(e) to indicate this subset. For the rest of this paper, we define that e ∈ Anc(e), and that e / ∈ Off(e). The third subset contains all the tree edges that do not appear neither in Off(e), nor in Anc(e); we indicate this subset by Ind(e).
For any two tips u, v ∈ S, we use p(u, v) to indicate the simple path in T between these nodes. Of course, the path p(u, v) is unique since T is a tree. We use cost (u, v) to denote the cost of this path, that is the sum of the weights of all the edges that appear on the path. Let u be a tip in S and let e be an edge in E. We use cost(u, e) to represent the cost of the shortest simple path between u and the child node of e. Therefore, if u ∈ S(e) this path does not include e, otherwise it does. For any subset R ⊆ S of the tips of the tree T , we denote the set of all pairs of elements in R, that is the set of all combinations that consist of two distinct tips in R, by (R). Given a phylogenetic tree T and a subset of its tips R ⊆ S, we denote the Mean Pairwise Distance of R in T by MPD(T , R). Let r = |R|. This measure is equal to: cost (u, v) .

Aggregating the costs of paths
Let T be a phylogenetic tree that consists of n edges and s tips, and let r be a positive integer such that r ≤ s. We use sk(T , r) to denote the skewness of the MPD on T when we pick a subset of r tips of this tree with uniform probability. In the rest of this paper we describe in detail how we can compute sk(T , r) in O(n) time, by scanning T only a constant number of times. Based on the formal definition of skewness, the value of sk(T , r) is equal to: where expec(T , r) and var(T , r) are the expectation and the variance of the MPD for subsets of exactly r tips in T , and E R∈Sub (S,r) [ ·] denotes the function of the expectation over all subsets of exactly r tips in S. In a previous paper, we showed how we can compute the expectation and the variance of the MPD on T in O(n) time [2]. Therefore, in the rest of this work we focus on analysing the value E R∈Sub(S,r) [ MPD 3 (T , R)] and expressing this quantity in a way that can be computed efficiently, in linear time with respect to the size of T .
To make things more simple, we break the description of our approach into two parts; in the first part, we define several quantities that come from adding and multiplying the costs of specific subsets of paths between tips of the tree. We also present how we can compute all these quantities in O(n) time in total by scanning T a constant number of times. Then, in the next section, we show how we can express the skewness of the MPD on T based on these quantities, and hence compute the skewness in O(n) time as well. Next we provide the quantities that we want to consider in our analysis; these quantities are described in Table 1. In this table but also in the rest of this work, for any tip u ∈ S, we consider that SQ(u) = SQ(e), and TC(u) = TC(e), such that e is the edge whose child node is u.
We provide now the following lemma.

Lemma 1.
Given a phylogenetic tree T that consists of n edges, we can compute all the quantities that are presented in Table 1 in O(n) time in total.
Proof. Each of the quantities (I)-(X) in Table 1 can be computed by scanning a constant number of times the input tree T , either bottom-up or top-to-bottom. For computing quantity (XI) we follow a more involved divideand-conquer approach.
We showed in a previous paper how we can compute quantity (I) and the quantities in (III) for all e ∈ E in O(n) time in total [2]. http://www.almob.org/content/9/1/15 Table 1 The quantities that we use for expressing the skewness of the MPD For an edge e ∈ E, the quantity in (VII) can be written as: We can compute this quantity for every e ∈ E in linear time as follows; in the first scan we compute for every edge e the number of leaves s(e) in T (e). This can be done in O(n) time by computing in a bottom-up manner s(e) as the sum of the numbers of tips s(e ), ∀e ∈ Ch(e). Then, we can compute TC sub (e) by scanning bottom-up the tree using the following formula: For quantity (VIII), for any e ∈ E we have that: Then SQ sub (e) can be computed for every edge e ∈ E by scanning T bottom up and evaluating the formula: For every edge e in T , quantity (IV) can be written as: In the last expression, value NumPath(e, l, k) is equal to the number of simple paths that connect two tips in T and which also contain all three edges e, l and k. The quantity NumPath(e, l) is equal to the number of simple paths that connect two tips in T and which also contain both edges e and l. Therefore, for any e ∈ E we have: l∈Ind(e) w l · s(l) (2) http://www.almob.org/content/9/1/15 We explain now how we can compute the six quantities in (2) in O(n) time, assuming that we have already computed TC sub (e) and s(e) for every e ∈ E. To make the description simpler, we show in detail how we can compute the second and fourth quantities that appear in the last expression; it is easy to show that the rest of the quantities in (2) can be calculated in a similar manner.
For any e ∈ E, we denote the second quantity as follows: We also define the following quantities: We can calculate SUM 1 (e) for every edge e by traversing the tree top-to-bottom and evaluating the following expressions: .
To compute the fourth quantity in (2), we use the following quantity: This quantity can be evaluated in O(n) time for every e ∈ E with a bottom-up scan of the tree. We also consider the following value which we can precompute in O(n) time: For every edge e ∈ E we calculate in a top-to-bottom manner the formula: SUM 3 (e) = w e (2·TC sub (e)+w e ·s(e))+SUM 3 (parent(e)) .
Then for each tree edge e, the fourth quantity in (2) can be computed in constant time as follows: The remaining quantities in (2) can be computed in a quite similar manner as the two quantities that we already described. Table 1 is equal to:

Quantity (II) in
We have already presented how to compute SQ(e) for every edge e in T in O(n) time in total, hence we can also compute CB(T ) in O(n) time by simply summing up the values w e · SQ(e) for every edge e in the tree. For quantity (V) it holds that: Since we have already computed TC(v) for every tip v ∈ S, we can trivially evaluate v∈S TC(v) in O(n) time. Hence, to compute quantity (V) it remains now to calculate the values SUM 4 (e) = u∈S(e) TC(u) for every edge e ∈ E. We can do this in O(n) time as follows: at each tip u ∈ S we store the value TC(u) that we have already computed. Then we scan T bottom-up and we calculate SUM 4 (e) by summing up the values SUM 4 (l) for all edges l ∈ Ch(e).
Let u be a tip in S, and let e be the edge which is adjacent to u. Then, quantity (VI) is equal to: In the last expression, value v∈S TC(v) can be computed in O(n) time, given that we have already computed TC(v) for every v ∈ S. Value l∈E w l v∈S(l) TC(v) and values x∈S(l) TC(x) for any l ∈ E can be calculated with a bottom-up scan of T in a similar way as we computed TC sub (e) for all e ∈ E. The remaining sums that involve edges in Anc(e) can be computed in linear time for every edge e with a similar mechanism as with SUM 3 (e) that http://www.almob.org/content/9/1/15 we described earlier in this proof. For any edge e ∈ E, quantities PC(e) and PSQ(e) in Table 1 From the two last expressions, and given the description that we provided for other similar quantities in Table 1, it easy to conclude that PC(e) can be evaluated for every edge e in O(n) time by scanning T a constant number of times. Having computed PC(e) for all edges e ∈ E, the quantity PSQ(e) can be computed in a similar manner.
Next we describe a divide-and-conquer approach for computing in (n) time quantity (XI) in Table 1 for every e ∈ E . Before we start our description, we define one more quantity that will help us simplify the rest of this proof. For an edge e ∈ E and a tip u ∈ S(e) we define that TC e (u) is equal to: For any edge e ∈ E it is easy to show that: Therefore, according to (3) we can compute the sum u∈S(e) TC e (u) for all edges e ∈ E in linear time in total, given that we have already computed TC(e) for every e ∈ E, and TC(u) for every u ∈ S.
Next we continue our description for computing QD(e) using a divide-and-conquer approach. We start with the base case; for every edge tree e that is adjacent to a leaf node we have: The first term of the sum in (5) can be expressed as: and can be computed in (|Ch(e)|) time, given that we have already computed SQ sub (l), ∀l ∈ Ch(e). http://www.almob.org/content/9/1/15 The next two parts of the sum in (5) are equal to:

l∈Ch(e) u∈S(l) v∈S(e)\S(l) x∈S(e)\S(l)
cost(u, l) · cost(v, l) The last expression can be computed in (|Ch(e)|) time as well, if we have already computed the sum k∈Ch(e) w k · s(k) and the quantity TC sub (e) for every edge e in the tree. We can rewrite the remaining term in (5) The last expression can be computed in (|Ch(e)|) time in a similar way as the previous terms of the sum in (5).
We left for the end the description of the calculation of the second sum in (4). We can express this sum as follows: We start with the second sum in (9). For this sum we get: (w k · s(k)) − w l · s(l) which takes (|Ch(e)|) time to be computed for each edge e.
To compute the first sum in (9) efficiently, we need to precompute for every edge l ∈ E the following quantity: u∈S(e) cost(u, e) · TC e (u) . http://www.almob.org/content/9/1/15 To do this, we follow again a divide-and-conquer approach. We get the base case for this computation for the edges of T that are adjacent to tips. For any such edge e we have: For any other edge e ∈ E we can compute this quantity based on the respective quantities of the edges in Ch(e). In particular, we have that: The first two sums in the last expression can be computed in (|Ch(e)|) time, given that we have computed already for every l ∈ Ch(e) the quantity TC(l) and the sum u∈S(l) TC(u) (can be done with a single bottom-up scan of the tree). The last sum in (10) can be expressed as: The two last sums in (11) are identical with the quantities that we analysed in (6) and in (7). Finally, the first sum in (11) which can also be computed in (|Ch(e)|) time.
All the sums that we analysed from (4) up to (12) can be computed in (|Ch(e)|) time for every edge e in the tree. From this we conclude that for every edge e ∈ E we can evaluate QD(e) in (4) in (|Ch(e)|) time from the respective values of the edges in Ch(e). Since e∈E |Ch(e)| = (|E|), we prove that we can compute QD(e) for all the edges in T in (n).

Computing the skewness of the MPD
In the previous section we defined the problem of computing the skewness of the MPD for a given phylogenetic tree T . Given a positive integer r ≤ s, we showed that to solve this problem efficiently it remains to find an efficient algorithm for computing E R∈Sub(S,r) [ MPD 3 (T , R)]; this is the mean value of the cube of the MPD among all possible subsets of tips in T that consist of exactly r elements. To compute this efficiently, we introduced in Table 1 eleven different quantities which we want to use in order to express this mean value. In Lemma 1 we proved that these quantities can be computed in O(n) time, where n is the size of T .
Next we prove how we can calculate the value for the mean of the cube of the MPD based on the quantities in Table 1. In particular, in the proof of the following lemma we show how the value E R∈Sub(S,r) [ MPD 3 (T , R)] can be written analytically as an expression that contains the quantities in Table 1. This expression can then be straightforwardly evaluated in O(n) time, given that we have already computed the aforementioned quantities. Because the full form of this expression is very long (it consists of a large number of terms), we have chosen not to include it in the definition of the following lemma. We chose to do so because we considered that including the http://www.almob.org/content/9/1/15 entire expression would not make this work more readable. In any case, the full expression can be easily infered from the proof of the lemma.

Lemma 2. For any given natural r ≤ s, we can compute
Proof. The expectation of the cube of the MPD is equal to: From the last expression we get: where AP R (u, v, x, y, c, d) is a random variable whose value is equal to one in the case that u, v, x, y, c, d ∈ R, otherwise it is equal to zero. For any six tips u, v, x, y, c, d ∈ S, which may not be all of them distinct, we use θ (u, v, x, y, c, d) to denote the number of distinct elements among these tips. Let t be an integer, and let (t) k denote the k-th falling factorial power of t, which means that (t) k = t(t − 1) . . . (t − k + 1). For the expectation of the random variables that appear in the last expression it holds that: Notice that in (14) we have 2 ≤ θ(u, v, x, y, c, d) ≤ 6. The value of the function θ(·) cannot be smaller than two in the above case because we have that u = v, x = y, and c = d. Thus, we can rewrite (13) as: (u,v,x,y,c,d) · cost(u, v) Hence, our goal now is to compute a sum whose elements are the product of costs of triples of paths. Recall that for each of these paths, the end-nodes of the path are a pair of distinct tips in the tree. Although the endnodes of each path are distinct, in a given triple the paths may share one or more end-nodes with each other. Therefore, the distinct tips in any triple of paths may vary from two up to six tips. Indeed, in (15) we get a sum where the triples of paths in the sum are partitioned in five groups; a triple of paths is assigned to a group depending on the number of distinct tips in this triple. In (15) the sum for each group of triples is multiplied by the same factor (r) θ(u,v,x,y,c,d) /(s) θ(u,v,x,y,c,d) , hence we have to calculate the sum for each group of triples separately.
However, when we try to calculate the sum for each of these groups of triples we see that this calculation is more involved; some of these groups of triples are divided into smaller subgroups, depending on which end-nodes of the paths in each triple are the same. To explain this better, we can represent a triple of paths schematically as a graph; let {u, v}, {x, y}, {c, d} ∈ (S) be three pairs of tips in T . As mentioned already, the tips within each pair are distinct, but tips between different pairs can be the same.
We represent the similarity between tips of these three pairs as a graph of six vertices. Each vertex in the graph corresponds to a tip of these three pairs. Also, there exists an edge in this graph between two vertices if the corresponding tips are the same. Thus, this graph is tripartite; no vertices that correspond to tips of the same pair can be connected to each other with an edge. Hence, we have a tripartite graph where each partite set of vertices consists of two vertices-see Figure 1 for an example.
For any triple of pairs of tips {u, v}, {x, y}, {c, d} ∈ (S) we denote the tripartite graph that corresponds to this triple by G [u, v, x, y, c, d]. We call this graph the similarity graph of this triple. Based on the way that similarities may occur between tips in a triple of paths, we can partition the five groups of triples in (15) into smaller subgroups. Each of these subgroups contains triples whose similarity graphs are isomorphic. For a tripartite graph that consists of three partite sets of two vertices each, there can be eight different isomorphism classes. Therefore, the five a b groups of triples in (15) are partitioned into eight subgroups. Figure 2 illustrates the eight isomorphism classes that exist for the specific kind of tripartite graphs that we consider. Since we refer to isomorphism classes, each of the graphs in Figure 2 represents the combinatorial structure of the similarities between three pairs of tips, and it does not correspond to a particular planar embedding, or ordering of the tips. Let X be any isomorphism class that is illustrated in Figure 2. We denote the set of all triples of pairs in (S) whose similarity graphs belong to this class by B X . More formally, the set B X can be defined as follows : v}, {x, y}, {c, d}} : {u, v}, {x, y}, {c, d} ∈ (S) and G [ u, v, x, y, c, d] belongs to class X in Figure 2} .
We introduce also the following quantity: Hence, we can rewrite (15) as follows: 3 · TRS(C) 5 · TRS(G) + 6 · (r) 6 (s) 6 · TRS(H) Notice that some of the terms (r) i (s) i · TRS(X) in (16) are multiplied with an extra constant factor. This happens for the following reason; the sum in TRS(X) counts each triple once for every different combination of three pairs of tips. However, in the triple sum in (15) some triples appear more than once. For example, every triple that belongs in class B appears three times in (15), hence there is an extra factor three in front of TRS(B) in (16).
To compute efficiently E R∈Sub(S,r) [ MPD 3 (T , R)], it remains to compute efficiently each value TRS(X) for every isomorphism class X that is presented in Figure 2. Next we show in detail how we can do that by expressing each quantity TRS(X) as a function of the quantities that appear in Table 1.
For the triples that correspond to the isomorphism class A we have: For TRS(B) we get: The quantity TRS(C) is equal to: The eight isomorphism classes of a tripartite graph of 3 × 2 vertices that represent schematically the eight possible cases of similarities between tips that we can have when we consider three paths between pairs of tips in a tree T . http://www.almob.org/content/9/1/15 The first of the two sums in (18) can be written as: According to Lemma 2, we can compute SM(u) and SQ(u) for all tips u ∈ S in linear time with respect to the size of T . Given these values, we can compute Combining the analyses that we did from (17) up to (26) we get: The value of TRS(D) can be expressed as: For TRS(E) we get: We can rewrite TRS(F) as follows: For the value of TRS(G) we have: We now break the sum in (27) into five pieces and express each piece of this sum in terms of the quantities in Table 1. The first piece of the sum is equal to: The second piece that we take from the sum in (27) can be expressed as: The next piece that we select from (27) is equal to: (29) http://www.almob.org/content/9/1/15 For the fourth piece of the sum in (27) we get: The last piece of the sum in (27) can be expressed as: Combining our analyses from (27) up to (31) we get: We can express TRS(H) using the values of the other isomorphism classes: We get the value of E R∈Sub(S,r) [ MPD 3 (T , R)] by plugging into (16) the values that we got for all eight isomorphism classes of triples. For any isomorphism class X we showed that the value TRS(X) can be computed by using the quantities in Table 1. The lemma follows from the fact that each quantity that appears in this table is used a constant number of times for computing value TRS(X) for any class X, and since we showed that we can precompute all these quantities in (n) time in total. Theorem 3. Let T be a phylogenetic tree that contains s tips, and let r be a natural number with r ≤ s. The skewness of the mean pairwise distance on T among all subsets of exactly r tips of T can be computed in (n) time.  Proof. According to the definition of skewness, as it is also presented in (1), we need to prove that we can compute in (n) time the expectation and the variance of the MPD, and the value of the expression E R∈Sub(S,r) [ MPD 3 (T , R)]. In a previous paper we showed that the expectation and the variance of the MPD can be computed in (n) time. By combining this with Lemma 2 we get the proof of the theorem.

Experiments: improved P value estimation incorporating skewness
Earlier in this paper, we mentioned that distributions of MPD values are often found to be skewed, suggesting that it is necessary to incorporate this skewness into analytical P value estimation. However, it is unclear whether good P value estimates are possible with only the first three moments of the distribution, or if more detailed distributional information is required.
We investigate this question here by considering random phylogenetic trees produced by a pure birth process [12], though results were qualitatively identical when using trees generated by a combined birth-death process (and skewness did not vary as a function of the death rate). We took two approaches for estimating the position of the 2.5 and 97.5 percentile of MPD distribution given a particular tree instance. For any tree T that we Figure 4 Comparison of approximation methods. Errors in P value approximation using different resampling replicates (indicated by the coloured lines), compared to that obtained by assuming a skew-normal distribution of MPD values (indicated as SN). Errors were strongly influenced by tip set size r, and weakly by tree size; on the left side appear the results for a 500 tip tree, and on the right for a 2000 tip tree). In most cases, P value approximation based on the skew-normal distribution performed better than the most commonly-used standard of 1000 set resamplings (blue line), and the relative performance of the skew-normal approach improved with increasing tip set size. http://www.almob.org/content/9/1/15 constructed, we first calculated the distribution of the MPD values using as a point of reference extensive sampling of sets of tips (much more extensive than is usually employed in practice). In particular, for specific values of r we sampled from T a large number of sets that consist of exactly r tips (see Table 2 for the values of r and numbers of sets that we sampled). We simply calculated the percentiles of these distributions, and call these the reference values, recognizing that they neveretheless contain some error, being incomplete samples from the tree. Complete sampling from large trees is computationally infeasible, but we estimate that the error in the calculated percentiles was less than 0.05 distance units in all cases (corresponding to an error of approximately 0.01% relative to the mean MPD-see Figure 3).
The two approaches that we used to estimate the percentile positions reflect two alternatives that might be employed by practising researchers. In the first approach, for each value r that we considered, we sampled again several sets of tips, yet much fewer than the ones we used to calculate the reference values (100, 500, 1000 or 5000 sets). We then compare the absolute difference between percentiles estimated in this manner and the reference values. We refer to this difference as the error between the estimated percentile values and the reference values. The second approach uses the mean, variance and skewness of the MPD distribution to determine the position of the 2.5 and 97.5 percentile of the skew-normal distribution with these moments [13]. The mean, variance and skewness were computed in this case based on all the MPD values that we used to calculate reference percentiles. Although we have implemented algorithms for computing the exact values of the mean and variance of the MPD, we have not implemented so far the algorithm that computes the skewness of the MPD; that is the algorithm outlined in the previous sections of this paper. As with the previous approach, the error of this approximation method was calculated by taking the absolute difference between each estimated percentile position and the corresponding reference value.
The experiment described above was repeated across 100 replicate trees of each of two sizes (500 and 2000 tips), and across a range of tip set sizes (10, 20, 40, 80, 160 and 320). Errors were weakly related to tree size but decreased strongly with tip set size-see Figure 4. This decrease was more pronounced for estimates based on skew-normal approximation than resampling. Notably, the skew-normal approximation yielded smaller errors than the most commonly used standard of 1000 resamplings for all but the smallest tip set sizes.
Thus, we conclude that the errors introduced by assuming a skew-normal distribution of MPD values appear to be comparable to or smaller than those introduced by standard resampling procedures, while also showing better scaling with increased tip sample size. Finally, the computation of P values using skew-normal approximation is typically faster than with resampling, particularly in cases involving large samples of trees.

Conclusions
Given a rooted tree T and a non-negative integer r, we proved that we can compute the skewness of the MPD among all subsets of r leaves in T in O(n) time. An interesting problem for future research would be to implement the algorithm that is outlined by our proof, and show its efficiency in practice. Also, it would be interesting to derive a similar result for the so-called Community Distance measure; this is the equivalent of the MPD when distances between two sets of species are considered [11].