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

- Constantinos Tsirogiannis
^{1, 2}Email author and### Affiliated with

- Brody Sandel
^{1, 2}### Affiliated with

**9**:15

**DOI: **10.1186/1748-7188-9-15

© Tsirogiannis and Sandel; licensee BioMed Central Ltd. 2014

**Received: **9 December 2013

**Accepted: **27 May 2014

**Published: **14 June 2014

## Abstract

### 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.

### Keywords

Algorithms for phylogenetic trees Mean pairwise distance Skewness## 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
and a subset of species *R* that are represented by leaf nodes of
, the MPD of the species in *R* is equal to average cost of all possible simple paths that connect pairs of nodes in *R*.

*R*is large or small, we need to know the average value (expectation) of the MPD for all sets of species in 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 . Putting all these values together, we get the following index that expresses how clustered are the species in

*R*[1]:

where
is the value of the MPD for *R* in
, and
and
are the expected value and the standard deviation respectively of the MPD calculated over all subsets of *r* species in
.

In a previous paper we presented optimal algorithms for computing the expectation and the standard deviation of the MPD of a phylogenetic tree
in *O*(*n*) time, where *n* is the number of the edges of
[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
, 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
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
. 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
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
optimally, in *Θ*(*n*) time. For the calculation of this skewness value we consider that every subset of exactly *r* species in
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
[11].

The rest of this paper is, almost in its entirety, an elaborate proof for computing the skewness of the MPD on a tree
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 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
be a phylogenetic tree, and let *E* be the set of its edges. We denote the number of the edges in
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
. 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
is the subtree of
whose root is the child node of edge *e*. We denote the set of tips that appear in
as *S*(*e*), and we denote the number of these tips by *s*(*e*).

Given any edge *e*∈*E*, we partition the edges of
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*).

*u*,

*v*∈

*S*, we use

*p*(

*u*,

*v*) to indicate the simple path in between these nodes. Of course, the path

*p*(

*u*,

*v*) is unique since is a tree. We use

*c*

*o*

*s*

*t*(

*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

*c*

*o*

*s*

*t*(

*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 , 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 and a subset of its tips

*R*⊆

*S*, we denote the Mean Pairwise Distance of

*R*in by . Let

*r*=|

*R*|. This measure is equal to:

### Aggregating the costs of paths

*n*edges and

*s*tips, and let

*r*be a positive integer such that

*r*≤

*s*. We use to denote the skewness of the MPD on 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 in

*O*(

*n*) time, by scanning only a constant number of times. Based on the formal definition of skewness, the value of is equal to:

where
and
are the expectation and the variance of the MPD for subsets of exactly *r* tips in
, 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
in *O*(*n*) time [2]. Therefore, in the rest of this work we focus on analysing the value
and expressing this quantity in a way that can be computed efficiently, in linear time with respect to the size of
.

*O*(

*n*) time in total by scanning a constant number of times. Then, in the next section, we show how we can express the skewness of the MPD on 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*.

**The quantities that we use for expressing the skewness of the MPD**

I) | II) |
---|---|

III) | IV) |

V) | VI) |

VII) | VIII) |

IX) | X) |

XI) |

We provide now the following lemma.

#### Lemma 1

Given a phylogenetic tree
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 , either bottom-up or top-to-bottom. For computing quantity (XI) we follow a more involved divide-and-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].

*e*∈

*E*, the quantity in (VII) can be written as:

*e*∈

*E*in linear time as follows; in the first scan we compute for every edge

*e*the number of leaves

*s*(

*e*) in . 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

*T*

*C*

_{ s u b }(

*e*) by scanning bottom-up the tree using the following formula:

*e*∈

*E*we have that:

*S*

*Q*

_{ s u b }(

*e*) can be computed for every edge

*e*∈

*E*by scanning bottom up and evaluating the formula:

*e*in , quantity (IV) can be written as:

*e*,

*l*,

*k*) is equal to the number of simple paths that connect two tips in 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 and which also contain both edges

*e*and

*l*. Therefore, for any

*e*∈

*E*we have:

We explain now how we can compute the six quantities in (2) in *O*(*n*) time, assuming that we have already computed T*C*
_{
s
u
b
}(*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.

*e*∈

*E*, we denote the second quantity as follows:

_{1}(

*e*) for every edge

*e*by traversing the tree top-to-bottom and evaluating the following expressions:

*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:

*e*∈

*E*we calculate in a top-to-bottom manner the formula:

*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.

*e*) for every edge

*e*in in

*O*(

*n*) time in total, hence we can also compute 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
in *O*(*n*) time. Hence, to compute quantity (V) it remains now to calculate the values
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
bottom-up and we calculate SUM_{4}(*e*) by summing up the values SUM_{4}(*l*) for all edges *l*∈Ch(*e*).

*u*be a tip in

*S*, and let

*e*be the edge which is adjacent to

*u*. Then, quantity (VI) is equal to:

*O*(

*n*) time, given that we have already computed TC(

*v*) for every

*v*∈

*S*. Value and values for any

*l*∈

*E*can be calculated with a bottom-up scan of in a similar way as we computed

*T*

*C*

_{ s u b }(

*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 we described earlier in this proof. For any edge

*e*∈

*E*, quantities PC(

*e*) and PSQ(

*e*) in Table 1 are equal to:

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
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.

*Θ*(

*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:

*e*∈

*E*it is easy to show that:

Therefore, according to (3) we can compute the sum
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*.

*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:

*e*∈

*E*that is not adjacent to a leaf node, we can calculate QD(

*e*) using the values of the respective quantities of the edges in Ch(

*e*):

*Θ*(|Ch(

*e*)|) time for each edge

*e*, given that we have already computed the values QD(

*l*) for every

*l*∈Ch(

*e*). We leave the description for calculating the second sum in (4) for the end of this proof. The third sum in this expression is equal to:

and can be computed in *Θ*(|Ch(*e*)|) time, given that we have already computed *S*
*Q*
_{
s
u
b
}(*l*),∀*l*∈Ch(*e*).

*Θ*(|Ch(

*e*)|) time as well, if we have already computed the sum and the quantity

*T*

*C*

_{ s u b }(

*e*) for every edge

*e*in the tree. We can rewrite the remaining term in (5) as:

The last expression can be computed in *Θ*(|Ch(*e*)|) time in a similar way as the previous terms of the sum in (5).

which takes *Θ*(|Ch(*e*)|) time to be computed for each edge *e*.

*l*∈

*E*the following quantity:

*e*we have:

*e*∈

*E*we can compute this quantity based on the respective quantities of the edges in Ch(

*e*). In particular, we have that:

*Θ*(|Ch(

*e*)|) time, given that we have computed already for every

*l*∈Ch(

*e*) the quantity TC(

*l*) and the sum (can be done with a single bottom-up scan of the tree). The last sum in (10) can be expressed as:

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
, we prove that we can compute QD(*e*) for all the edges in
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
. Given a positive integer *r*≤*s*, we showed that to solve this problem efficiently it remains to find an efficient algorithm for computing
; this is the mean value of the cube of the MPD among all possible subsets of tips in
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
.

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
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 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*
*in Θ (n) time.*

### Proof.

*A*

*P*

_{ 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:

*θ*(

*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:

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 end-nodes 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
. As mentioned already, the tips within each pair are distinct, but tips between different pairs can be the same.

*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 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.

*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 . More formally, the set can be defined as follows :

Notice that some of the terms
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
, 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.

*A*we have:

*B*) we get:

*C*) is equal to:

*e*∈

*E*we have that:

*u*) and SQ(

*u*) for all tips

*u*∈

*S*in linear time with respect to the size of . Given these values, we can compute for every edge

*e*∈

*E*in with a single bottom-up scan of the tree. For any edge

*e*in

*E*, the second sum in (18) is equal to:

*e*in

*Θ*(

*n*) time in total as follows; for every tip

*u*∈

*S*we store SQ(

*u*) together with this tip, and then scan bottom-up the tree adding those values that are in the subtree of each edge. For the remaining part of (20) we get:

*D*) can be expressed as:

*E*) we get:

*F*) as follows:

*G*) we have:

*H*) using the values of the other isomorphism classes:

We get the value of
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*
*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*
*among all subsets of exactly r tips of*
*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
. 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.

*r*we sampled from 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 sizes of tip samples that we considered for our experiments, together with the number of sets that we sampled for each tip size in order to derive the “true” values**

Size of each tip sample | Number of sampled sets |
---|---|

10 | 10 |

20 | 10 |

40 | 5·10 |

80 | 3·10 |

160 | 2·10 |

320 | 10 |

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.

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
and a non-negative integer *r*, we proved that we can compute the skewness of the MPD among all subsets of *r* leaves in
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].

## Declarations

### Acknowledgements

We want to thank Dr. Erick Matsen for his insightful comments on the paper and the productive discussion that we had regarding the motivation for deriving an algorithm for the MPD skewness.

## Authors’ Affiliations

## References

- Webb CO, Ackerly DD, McPeek MA, Donoghue MJ:
**Phylogenies and community ecology****.***Annu Rev Ecol Systemat*2002,**33:**475–505. 10.1146/annurev.ecolsys.33.010802.150448View Article - Tsirogiannis C, Sandel B, Cheliotis D:
**Efficient computation of popular phylogenetic tree measures****.**In*Proceedings of 12th International Workshop on Algorithmms in Bioinformatics (WABI)*. Edited by: Raphael B, Tang J. Springer-Verlag Berlin Heidelberg; 2012:30–43.View Article - Pontarp M, Canbäck B, Tunlid A, Lunberg P:
**Phylogenetic analysis suggests that habitat filtering is structuring marine bacterial communities across the globe****.***Microb Ecol*2012,**64:**8–17. 10.1007/s00248-011-0005-7View ArticlePubMed CentralPubMed - Jetz W, Thomas GH, Joy JB, Hartmann K, Mooers AO:
**The global diversity of birds in space and time****.***Nature*2012,**491:**444–448. 10.1038/nature11631View ArticlePubMed - Barnagaud J-Y, Daniel Kissling W, Sandel B, Eiserhardt WL, Şekercioğlu CH, Enquist BJ, Tsirogiannis C, Svenning J-C:
**Ecological traits influence the phylogenetic structure of bird species co-occurrences worldwide****.***Ecol Lett*2014. To appear - Cooper N, Rodríguez J, Purvis A:
**A common tendency for phylogenetic overdispersion in mammalian assemblages****.***Proc Biol Sci*2008,**275:**2031–2037. 10.1098/rspb.2008.0420View ArticlePubMed CentralPubMed - Vamosi JC, Vamosi SM:
**Body size, rarity, and phylogenetic community structure: Insights from diving beetle assemblages of alberta****.***Divers Distributions*2007,**13:**1–10. - Harmon-Threat AN, Ackerly DD:
**Filtering across spatial scales: phylogeny, biogeography and community structure in bumble bees****.***PLoS ONE*2013,**8:**60446. 10.1371/journal.pone.0060446View Article - Lin TI, Lee JC, Yen SY:
**Finite mixture modelling using the skew normal distribution****.***Statistica Sinica*2007,**17:**209–227. - Lee SX, McLachlan GJ:
**On mixtures of skew normal and skew t-distributions****.***Adv Data Anal Classif*2013,**7:**241–266. 10.1007/s11634-013-0132-8View Article - Swenson NG:
**Phylogenetic beta diversity metrics. trait evolution and inferring the functional beta diversity of communities****.***PLoS ONE*2011,**6**(6)**:**21264. 10.1371/journal.pone.0021264View Article - Harmon LJ, Weir JT, Brock CD, Glor RE, Challenger W:
**Geiger: investigating evolutionary radiations****.***Bioinformatics*2008,**24:**129–131. 10.1093/bioinformatics/btm538View ArticlePubMed - Azzalini A:
**The r ‘sn’ package: the skew-normal and skew-t distributions (version 1.0–0)****.**[http://cran.r-project.org/web/packages/sn/index.html] 2014. Accessed 2014–05–7

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.