Computing the Skewness of the Phylogenetic
Mean Pairwise Distance in Linear Time

arXiv:1307.7825v1 [q-bio.QM] 30 Jul 2013

Constantinos Tsirogiannis and Brody Sandel
MADALGO? and Department of Bioscience, Aarhus University, Denmark
{constant, brody.sandel}@cs.au.dk

Abstract. 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
T and for a set of species R represented by a subset of the leaf nodes of
T , the MPD of R is equal to the average cost of all possible simple paths
in T 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.
In the present work we describe how to compute the skewness of the
MPD on a tree T optimally, in Θ(n) time; here n is the size of the
tree T . 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 T that can be possibly used as building blocks
for computing efficiently the skewness of other phylogenetic measures.

1

Introduction

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 [7]. 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
?

Center for Massive Data Algorithmics, a Center of the Danish National Research
Foundation.

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 [7]:
NRI =

MPD(T , R) − expecMPD (T , r)
,
sdMPD (T , r)

where MPD(T , R) is the value of the MPD for R in T , and expec(T ) and
sdMPD (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
Θ(n) time, where n is the number of the edges of T [5]. 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, one 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 is a statistical P value, that is the probability that
we would observe an NRI value as extreme or more so if the community were
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.
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 [1,6]. Of particular concern,
this skewness introduces a bias towards detecting either significant clustering or
significant overdispersion [2]. Calculating this skewness analytically would enable
us to remove this bias and improve the accuracy of P value estimates obtained
analytically. However, so far there has been no result in the related literature
that shows how to compute this skewness value.
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 tackle the problem of computing 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 uniformly 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 very efficient, especially if we consider that it outperforms
the best algorithms that are known so far for computing lower-order statistics
for other phylogenetic measures; for example the most efficient known algorithm
for computing the variance of the popular Phylogenetic Distance (PD) runs in
O(n2 ) time [5].
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 [4].
The rest of this paper is, 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. In Section 3 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.

2

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 we 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 e0 ∈ E for which e
appears in the subtree of e0 . 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:
MPD(T , R) =

2.1

2
r(r − 1)

X

cost(u, v) .

{u,v}∈∆(R)

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 the concept of skewness, the value of sk(T , r) is equal to:
"
sk(T , r) = ER∈Sub(S,r)
=

MPD(T , R) − expec(T , r)
var(T , r)

3 #

ER∈Sub(S,r) [MPD3 (T , R)] − 3 · var(T , r)2 − expec(T , r)3
,
var(T , r)3

(1)

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 ER∈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 [5]. Therefore, in the rest of this work we focus on analysing

the value ER∈Sub(S,r) [MPD3 (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 Section 3, 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.
Table 1. The quantities that we use for expressing the skewness of the MPD.
I) TC(T ) =

X

cost(u, v)

II) CB(T ) =

{u,v}∈∆(S)

X

cost3 (u, v)

{u,v}∈∆(S)

X

III) ∀e ∈ E, TC(e) =

{u,v}∈∆(S)
e∈p(u,v)

X

V) ∀e ∈ E, Mult(e) =

X

IV) ∀e ∈ E, SQ(e) =

cost(u, v)

{u,v}∈∆(S)
e∈p(u,v)

TC(u) · TC(v)

VI) ∀u ∈ S, SM(u) =

{u,v}∈∆(S)
e∈p(u,v)

X

VII) ∀e ∈ E, TCsub (e) =

X

X

VIII) ∀e ∈ E, SQsub (e) =

cost(u, e)

X) ∀e ∈ E, PSQ(e) =

cost(u, e)

u∈S

X

cost2 (u, e)

u∈S(e)
2

cost (u, e)

2
X


u∈S(e)

X

X

u∈S


XI) ∀e ∈ E, QD(e) =

cost(u, v) · TC(v)

v∈S\{u}

u∈S(e)

IX) ∀e ∈ E, PC(e) =

cost2 (u, v)

cost(u, v)

v∈S(e)\{u}

For any tip u ∈ S, we define that SQ(u) = SQ(e), and TC(u) = TC(e),
where e is the edge whose child node is u. The proof of the following lemma is
provided in the full version of this paper.
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.

3

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 ER∈Sub(S,r) [MPD3 (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 ten 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 ER∈Sub(S,r) [MPD3 (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 1 .
Lemma 2. For any given natural r ≤ s, we can compute ER∈Sub(S,r) [MPD3 (T , R)]
in Θ(n) time.
Proof. The expectation of the cube of the MPD is equal to:
8
·
ER∈Sub(S,r) [MPD3 (T , R)] = 3
r (r − 1)3

X
X
X
ER∈Sub(S,r) 


cost(u, v) · cost(x, y) · cost(c, d) .

{u,v}∈∆(R) {x,y}∈∆(R) {c,d}∈∆(R)

From the last expression we get:
"
X
X
ER∈Sub(S,r)

#
X

cost(u, v) · cost(x, y) · cost(c, d)

{u,v}∈∆(R) {x,y}∈∆(R) {c,d}∈∆(R)

=

X

X

X

cost(u,v)·cost(x,y)·cost(c, d) ·

{u,v}∈∆(S) {x,y}∈∆(S) {c,d}∈∆(S)

ER∈Sub(S,r) [APR (u,v, x, y, c, d)] ,

(2)

where APR (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:
ER∈Sub(S,r) [APR (u, v, x, y, c, d)] =

(r)θ(u,v,x,y,c,d)
(s)θ(u,v,x,y,c,d)

(3)

Notice that in (3) 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 6= v,
x 6= y, and c 6= d. Thus, we can rewrite (2) as:
X
X
X
(r)θ(u,v,x,y,c,d)
· cost(u,v) · cost(x, y) · cost(c, d)
(s)θ(u,v,x,y,c,d)
{u,v}∈∆(S) {x,y}∈∆(S) {c,d}∈∆(S)

(4)
1

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

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 (4) 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 (4) 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

α β γ

δ



(a)

ζ

η
(b)

Fig. 1. (a) A phylogenetic tree T and (b) an example of the tripartite graph induced
by the triplet of its tip pairs {α, γ}, {δ, γ}, {, δ}, , where {α, γ, δ, } ⊂ S. The dashed
lines in the graph distinguish the partite subsets of vertices; the vertices of each partite
subset correspond to tips of the same pair.

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 Fig. 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 (4) 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 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 Fig. 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.

A:

B:

C:

D:

E:

F:

G:

H:

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

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 BX . More formally, the set BX can be defined as follows :
BX = { {{u, 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:
X
TRS(X) =
cost(u, v) · cost(x, y) · cost(c, d) .
{{u,v},{x,y},{c,d}}∈BX

Hence, we can rewrite (4) as follows:
(r)2
(r)3
(r)3
(r)4
· TRS(A) + 3 ·
· TRS(B) + 6 ·
· TRS(C) + 6 ·
· TRS(D)
(s)2
(s)3
(s)3
(s)4
+3·

(r)4
(r)4
(r)5
(r)6
· TRS(E) + 6 ·
· TRS(F ) + 6 ·
· TRS(G) + 6 ·
· TRS(H) (5)
(s)4
(s)4
(s)5
(s)6

i
Notice that some of the terms (r)
(s)i · TRS(X) in (5) 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 (4) some triples appear more than once. For example, every triple that belongs in class B appears three times in (4), hence there
is an extra factor three in front of TRS(B) in (5).
To compute efficiently ER∈Sub(S,r) [MPD3 (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:
X
TRS(A) =
cost3 (u, v) = CB(T ) .
{u,v}∈∆(S)

For TRS(B) we get:

X

TRS(B) =

cost2 (u, v) 

{u,v}∈∆(S)

cost(u, x) +

x∈S\{u}

X

=


X

X

cost(v, y) − 2 · cost(u, v)

y∈S\{v}

cost2 (u, v) (TC(u) + TC(v) − 2 · cost(u, v))

{u,v}∈∆(S)

X
SQ(u) · TC(u) − 2 · CB(T ) .
=
u∈S

The quantity TRS(C) is equal to:
1X
6

X

u∈S v∈S\{u}

=

cost(u, x) · cost(x, v)

x∈S\{u,v}

X
1X
we
6
e∈E

X

cost(u, v)
X

X

cost(u, x) · cost(x, v) .

(6)

u∈S(e) v∈S−S(e) x∈S\{u,v}

For any e ∈ E we have that:
X
X

X

cost(u, x) · cost(x, v)

u∈S(e) v∈S−S(e) x∈S\{u,v}

X

=

X

X

cost(u, x) · cost(x, v)

(7)

cost(u, x) · cost(x, v) .

(7b)

u∈S(e) v∈S\{u} x∈S\{u,v}

X

−2

X

{u,v}∈∆(S(e)) x∈S\{u,v}

The first of the two sums in (7) can be written as:
X
X
X
cost(u, x) · cost(x, v)
u∈S(e) v∈S\{u} x∈S\{u,v}

=

X

X

X

cost(u, v) · cost(v, x)

u∈S(e) v∈S\{u} x∈S\{u,v}

=

X

X

(cost(u, v) · TC(v) − cost2 (u, v))

u∈S(e) v∈S\{u}

=

X
u∈S(e)

SM(u) − SQ(u) .

(8)

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
P
u∈S(e) SM(u) − SQ(u) for every edge e ∈ E in T with a single bottom-up scan
of the tree. For any edge e in E, the second sum in (7b) is equal to:
X
X
cost(u, x) · cost(x, v)
{u,v}∈∆(S(e)) x∈S\{u,v}

X

=

X

cost(u, x) · cost(x, v)

(9)

{u,v}∈∆(S(e)) x∈S(e)\{u,v}

X

X

+

cost(u, x) · cost(x, v) .

(9b)

{u,v}∈∆(S(e)) x∈S\S(e)

We can express the first sum in (9) as:
X
X
cost(u, x) · cost(x, v)
{u,v}∈∆(S(e)) x∈S(e)\{u,v}


1 X 
=
2
u∈S(e)

2
X

v∈S(e)\{u}

1
1 X
= QD(e) −
2
2

1 X
cost(u, v) −
2

X

cost2 (u, v)

u∈S(e) v∈S(e)\{u}

X

cost2 (u, v) .

(10)

u∈S(e) v∈S(e)\{u}

The last sum in (10) is equal to:
X
X
X
cost2 (u, v) =
SQ(u) − SQ(e). .
u∈S(e) v∈S(e)\{u}

(11)

u∈S(e)

P
The value of the sum u∈S(e) SQ(u) can be computed for every edge 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 (9b) we get:
X
X
cost(u, x) · cost(x, v)
{u,v}∈∆(S(e)) x∈S\S(e)

=

X

X

(cost(u, e) + cost(x, e)) (cost(v, e) + cost(x, e))

{u,v}∈∆(S(e)) x∈S\S(e)

=

X

X

cost(u, e) · cost(v, e)

{u,v}∈∆(S(e)) x∈S\S(e)

+

X

X

cost(x, e) · (cost(u, e) + cost(v, e))

{u,v}∈∆(S(e)) x∈S\S(e)

+

X

X

{u,v}∈∆(S(e)) x∈S\S(e)

cost2 (x, e) .

(12)

The first sum in (12) is equal to:
X
X

cost(u, e) · cost(v, e) = (s − s(e)) TCsub 2 (e) − SQsub (e) .
{u,v}∈∆(S(e)) x∈S\S(e)

(13)
For the second sum in (12) we have:
X

X

cost(x, e) · cost(u, e) + cost(v, e)



{u,v}∈∆(S(e)) x∈S\S(e)

 X


= s(e) − 1
cost(x, e) · TCsub (e) = s(e) − 1 · TCsub (e) · PC(e) − TCsub (e) .
x∈S\S(e)

(14)

The last sum in (12) can be written as:
X

X

cost2 (x, e) =

{u,v}∈∆(S(e)) x∈S\S(e)

s(e)(s(e) − 1)
(PSQ(e) − SQsub (e)) . (15)
2

Combining the analyses that we did from (6) up to (15) we get:

X 
3
1
1
1X
we
SM(u) − SQ(u) + · QD(e) + · SQ(e)
TRS(C) =
6
2
2
2
e∈E

u∈S(e)

2s − 2 · s(e) + s(e)(s(e) − 1)
· SQsub (e)
2
s(e)(s(e) − 1)
· PSQ(e) .
+ (s(e) − 1) · TCsub (e) · PC(e) +
2
+ (s − 2s(e) + 1) · TCsub 2 (e) −

The value of TRS(D) can be expressed as:
X
X
cost(u, v) · cost(u, x) · cost(u, y)
u∈S

v,x,y∈S\{u}
v,x,y are distinct

1
=
6
=

X



3

TC (u) − 2 · TRS(A) − 3 · TRS(B)

u∈S

1 X
2
1
·
TC3 (u) + · CB(T ) − · SQ(u) · TC(u) .
6
3
2
u∈S

For TRS(E) we get:
X

X

cost2 (u, v) · cost(x, y)

{u,v}∈∆(S) {x,y}∈∆(S\{u,v})

=

X

cost2 (u, v)(TC(T ) − TC(u) − TC(v) + cost(u, v))

{u,v}∈∆(S)

= TC(T )

X
e∈E

we · TC(e) −

X
u∈S

(SQ(u) · TC(u)) + CB(T ) .

We can rewrite TRS(F ) as follows:

!
X

cost(u, v) TC(u) · TC(v) − cost2 (u, v) −

{u,v}∈∆(S)

=

X

cost(u, x) · cost(x, v)

x∈S\{u,v}

X

cost(u, v) · TC(u) · TC(v) − CB(T ) − 3 · TRS(C)

{u,v}∈∆(S)

=

X

we · Mult(e) − CB(T ) − 3 · TRS(C) .

e∈E

For the value of TRS(G) we have:

TRS(G) =

1
2

X

X

cost(u, v)

{u,v}∈∆(S)

cost(u, x) + cost(v, x)



TC(T )

x∈S\{u,v}

!
− TC(u) − TC(v) − TC(x) + cost(u, v) + cost(u, x) + cost(v, x)

.
(16)

We now break the sum in (16) 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:

1
2

X

cost(u, v)

{u,v}∈∆(S)

=

(cost(u, x) + cost(v, x)) · TC(T )

x∈S\{u,v}

X
1
· TC(T )
TC2 (u) −
2
u∈S

=

X

X

cost2 (u, v)

{u,v}∈∆(S)

X
X
1
· TC(T )
TC2 (u) −
we · TC(e) .
2
u∈S

e∈E

The second piece that we take from the sum in (16) can be expressed as:
−
=−
=−

1
2
1
2
1
2

X

cost(u, v)

{u,v}∈∆(S)

X

X

(cost(u, x) + cost(v, x)) (TC(u) + TC(v))

x∈S\{u,v}

cost(u, v) (TC(u) + TC(v) − 2 · cost(u, v)) (TC(u) + TC(v))

{u,v}∈∆(S)


cost(u, v) TC2 (u) + TC2 (u) + 2 · TC(u) · TC(v)

X
{u,v}∈∆(S)

− 2 · cost(u, v) · TC(u) + TC(v)



X
1X
cost(v, x) · TC(v) · TC(x)
TC3 (u) −
2
u∈S
{v,x}∈∆(S)
X

+
cost2 (y, z) TC(y) + TC(z)

=−

{y,z}∈∆(S)

=−

X
X
1X
TC3 (u) −
we · Mult(e) +
SQ(u) · TC(u) .
2
u∈S

e∈E

(17)

u∈S

The next piece that we select from (16) is equal to:
−

1
2

X

X

cost(u, v)

{u,v}∈∆(S)

1X
=−
2

X

1X
2

X

x∈S\{u,v}

X

cost(u, v)

u∈S v∈S\{u}

+


cost(u, x) + cost(v, x) · TC(x)
cost(u, x) · TC(x)

x∈S\{u}

cost2 (u, v) · TC(v)

u∈S v∈S\{u}

=−

1X
1
TC(u) · SM(u) +
2
2
u∈S

X

cost2 (u, v) (TC(u) + TC(v))

{u,v}∈∆(S)


1X
=
TC(u) SQ(u) − SM(u) .
2

(18)

u∈S

For the fourth piece of the sum in (16) we get:
1
2

X

cost2 (u, v)

{u,v}∈∆(S)

X

cost(u, x) + cost(v, x)

(19)

x∈S\{u,v}

1
1X
= · TRS(B) =
SQ(u) · TC(u) − CB(T ) .
2
2
u∈S

(20)

The last piece of the sum in (16) can be expressed as:
X
X
1
2
cost(u, v)
(cost(u, x) + cost(v, x))
2
{u,v}∈∆(S)

=

=

=

1
2
1
2

X

x∈S\{u,v}

cost(u, v)

{u,v}∈∆(S)

X

X


cost2 (u, x) + cost2 (v, x) + 3 · TRS(C)

x∈S\{u,v}


cost(u, v) SQ(u) + SQ(v) − 2 · cost2 (u, v) + 3 · TRS(C)

{u,v}∈∆(S)

1X
SQ(u) · TC(u) − CB(T ) + 3 · TRS(C) .
2

(21)

u∈S

Combining our analyses from (16) up to (21) we get:
X
X
1
TRS(G) = · TC(T )
TC2 (u) −
we (TC(e) + Mult(e))
2
u∈S

+

e∈E


1X
TC(u)· 5 · SQ(u)−SM(u) −TC2 (u) − 2 ·CB(T ) +3·TRS(C) .
2
u∈S

We can express TRS(H) using the values of the other isomorphism classes:
X
X
X
1
cost(u, v) · cost(x, y) · cost(c, d)
TRS(H) =
6
{u,v}∈∆(S) {x,y}∈∆(S) {c,d}∈∆(S)

− TRS(A) − 3 · TRS(B) − 6 ·TRS(C) − 6 · TRS(D)
− 3 · TRS(E) − 6 · TRS(F ) − 6 · TRS(G)
=

1
1
1
· TC3 (T ) − · TRS(A) − · TRS(B) − TRS(C) − TRS(D)
6
6
2

−

1
· TRS(E) − TRS(F ) − TRS(G) .
2

We get the value of ER∈Sub(S,r) [MPD3 (T , R)] by plugging into (5) 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.
t
u
Theorem 1. 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 ER∈Sub(S,r) [MPD3 (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.
t
u

References
1. Cooper,N., Rodrı́guez,J., Purvis,A.: A Common Tendency for Phylogenetic
Overdispersion in Mammalian Assemblages. In: Proceedings of the Royal Society
B: Biological Sciences, vol. 275, pp. 2031–2037 (2008)
2. Harmon-Threatt,A.N., AckerlyD.D.: Filtering Across Spatial Scales: Phylogeny,
Biogeography and Community Structure in Bumble Bees. PLoS ONE 8:e60446,
(2013)
3. PontarpM., Canbäck,B., Tunlid,A., Lundberg,P.: Phylogenetic Analysis Suggests
that Habitat Filtering Is Structuring Marine Bacterial Communities Across the
Globe. Microbial Ecology 64, 8–17 (2012)
4. Swenson,N.G.: Phylogenetic Beta Diversity Metrics. Trait Evolution and Inferring
the Functional Beta Diversity of Communities. PLoS ONE 6(6):e21264, (2011)
5. Tsirogiannis,C., Sandel,B. Cheliotis,D.: Efficient Computation of Popular Phylogenetic Tree Measures. In: Proc. Workshop on Algorithmms in Bioinformatics
(WABI), pp. 30–43 (2012)
6. Vamosi,J.C., VamosiS.M.: Body Size, Rarity, and Phylogenetic Community Structure: Insights from Diving Beetle Assemblages of Alberta. Diversity and Distributions 13, 1–10 (2007)
7. Webb,C.O., Ackerly,D.D., McPeek,M.A., Donoghue,M.J.: Phylogenies and Community Ecology. Annual Review of Ecology and Systematics 33, 475–505 (2002)

