adiv (version 2.2.1)

evopca: Evolutionary Principal Component analyses

Description

The functions evopcachord and evopcahellinger perform two adaptations of Principal Component Analysis (PCA) for the analysis of phylogenetic diversity patterns across species communities: the evolutionary PCA based on Chord distance (evoPCAChord) and the evolutionary PCA based on Hellinger distance (evoPCAHellinger) (Pavoine 2016).

Usage

evopcachord(phyl, comm, option = c("centred", "decentred"), 
w = c("evoab", "even", "speciesab"), scannf = TRUE, 
nf = 2, abundance = TRUE)
evopcahellinger(phyl, comm, option = c("centred", "decentred"),
w = c("evoab", "even", "speciesab"), scannf = TRUE, 
nf = 2, abundance = TRUE)

Value

evopcachord and evopcahellinger both return an object of class evopca and of class dudi (see package ade4). Graphical tools are associated with class dudi in packages ade4 and adegraphics (see section "Examples" below).

The returned object contains the following components:

tab

a data frame with n rows and p columns, with communities as rows and nodes of the phylogeny as columns. Internal data frame used by the algorithm;

cw

weights attributed to the nodes of the phylogeny, a vector with p components;

lw

weights attributed to the communities, a vector with n components;

eig

vector of eigenvalues;

rank

integer, number of axes;

nf

integer, number of kept axes;

c1

normed scores for the nodes of the phylogeny, data frame with p rows and nf columns;

l1

normed scores for the communities, data frame with n rows and nf columns;

co

scores for the nodes of the phylogeny, data frame with p rows and nf columns;

li

scores for the communities, data frame with n rows and nf columns;

call

the original call.

If X is an object of class evopca, then attributes(X)$phy contains the phylogenetic tree (of class phylo) with names for internal nodes.

Arguments

phyl

an object inheriting the class phylo (see package ape), phylo4 (see package phylobase) or hclust. To ease the interpretation of factorial maps, it is advised to use a rooted phylogenetic tree and to add labels to the nodes of the phylogeny. For example, If your phylogeny is named 'tree' and is of class phylo, you can root the phylogeny using tree$root.edge <- 0 and you can name the nodes with the following command: tree$node.label <- paste("n", 1:tree$Nnode, sep=""). Use plot(tree, show.node.label=TRUE) to see the result. If it is rooted and of class phylo4, you can use: nodeLabels(tree) <- paste("n", 1:nNodes(tree), sep="").

comm

a data frame or a matrix typically with communities as rows, species as columns and presence/absence (1/0) or an index of abundance as entries. Species should be labeled as in the phylogenetic tree where they are the tips.

option

a string: either "centered" or "decentered". If several strings are given, only the first one is used. See details.

w

a string: either "evoab", "even", or "speciesab". If several strings are given, only the first one is used. See details.

scannf

a logical value indicating whether the screeplot (eigenvalues) should be displayed for choosing the number of axes to be kept.

nf

if scannf is FALSE, an integer indicating the number of axes to be kept.

abundance

a logical value. If FALSE, only the presences/absences of the evolutionary units in communities are considered.

Author

Sandrine Pavoine

Details

Let \(L_k\) be the length of branch k in the phylogenetic tree (out of K branches); \(a_{jk}\) the sum of abundances, in community j, for all species descending from branch k; \(w_j\) a positive weight attributed to community j (the definition for \(w_j\) is flexible with the only requirement that \(\sum_{j=1}^m w_j=1\)); \(a_{j+}=\sum_{k=1}^K L_k a_{jk}\); \(a_{+k}=\sum_{j=1}^m a_{jk}\); \(a_{++}=\sum_{j=1}^m\sum_{k=1}^K L_k a_{jk}\).

The weights of the communities (argument w) can be "even" (even weights, i.e. relative abundances are considered for evolutionary units), "evoab" (proportional to the summed abundances of all evolutionary units, i.e. absolute abundances are considered for evolutionary units), or "speciesab" (proportional to the summed abundances of all species). Note that if the phylogenetic tree is ultrametric (the distance from any species to the root is constant), then options "evoab" and "speciesab" are equivalent.

In evopcahellinger, as recommended by Rao (1995), the vector used to centre matrix \((\sqrt{a_{jk}/a_{j+}})_{j,k}\) in PCA can be defined as $$\left(\sum_{j=1}^m w_j \sqrt{a_{jk}/a_{j+}}\right)_k$$ (ordinary weighted mean, option "centered") or as $$\left(\sqrt{(\sum_{j=1}^m w_j a_{jk})/(\sum_{k=1}^K L_k \sum_{j=1}^m w_j a_{jk})}\right)_k$$ (option "decentered"); see Pavoine (2016) for an introduction of all ordination approaches.

Similarly, in evopcachord, the vector used to centre matrix \((a_{jk}/\sqrt{\sum_{k=1}^K a_{jk}^2})_{j,k}\) can be defined as $$\left(\sum_{j=1}^m w_j a_{jk}/\sqrt{\sum_{k=1}^K a_{jk}^2}\right)_k$$ (ordinary weighted mean, option "centered") or as $$\left((\sum_{j=1}^m w_j a_{jk})/\sqrt{\sum_{k=1}^K L_k (\sum_{j=1}^m w_j a_{jk})^2}\right)_k$$ (option "decentered").

References

Pavoine, S. (2016) A guide through a family of phylogenetic dissimilarity measures among sites. Oikos, 125, 1719--1732.

Rao, C.R. (1995) A review of canonical coordinates and an alternative to correspondence analysis using Hellinger distance. Questiio, 19, 23--63.

See Also

evoCA, evoNSCA, evodiss

Examples

if (FALSE) {
if(require(ade4) && require(ape) && require(adegraphics)){
O <- adegpar()$plabels$optim
adegpar("plabels.optim" = TRUE)

data(batcomm)
ab <- batcomm$ab
phy <- read.tree(text=batcomm$tre)
plot(phy, show.node=TRUE)

evopcaHbat <- evopcahellinger(phy, ab, scan=FALSE, nf=3)
dist(evopcaHbat$li)
evodiss(phy, ab, "Hellinger")

evopcaHbat$eig/sum(evopcaHbat$eig)
s.label(evopcaHbat$li)
s.label(evopcaHbat$co)
s.arrow(evopcaHbat$co)
inertia.dudi(evopcaHbat, row=TRUE)$row.abs
inertia.dudi(evopcaHbat, col=TRUE)$col.abs

evopcaCbat <- evopcachord(phy, ab, scan=FALSE, nf=3)
dist(evopcaCbat$li)
evodiss(phy, ab, "Chord")

evopcaCbat$eig/sum(evopcaCbat$eig)
s.label(evopcaCbat$li)
s.label(evopcaCbat$co)
s.arrow(evopcaCbat$co)
inertia.dudi(evopcaCbat, row=TRUE)$row.abs
inertia.dudi(evopcaCbat, col=TRUE)$col.abs

adegpar("plabels.optim" = O)
}
}