library(adephylo) data(lizards) head(lizards$traits) table.value(lizards$traits) table.value(data.frame(scale(lizards$traits))) ## simple pca pca1 <- dudi.pca(lizards$traits) summary(pca1) scatter(pca1) ## remove size effect newdat <- data.frame(lapply(log(lizards$traits), function(v) residuals(lm(v~log(lizards$traits$mean.L))))) rownames(newdat) <- rownames(lizards$traits) newdat <- newdat[,-1] pca1 <- dudi.pca(newdat) ## build a tree library(ape) library(phylobase) lizards$hprA liz.tre <- read.tree(text = lizards$hprA) plot(liz.tre) summary(liz.tre) ## associate data and phylogeny liz.4d0 <- phylo4d(liz.tre, log(lizards$traits)) par(mar=rep(.1,4)) table.phylo4d(liz.4d0,show.node=FALSE, cex.lab=1.2) liz.4d <- phylo4d(liz.tre, newdat) table.phylo4d(liz.4d,show.node=FALSE, cex.lab=1.2) ## Remove duplicated populations liz.4d <- prune(liz.4d, c(7,14)) table.phylo4d(liz.4d) ## get/change some information tipLabels(liz.4d) tipData(liz.4d) ## from phylogeny to distances d1 <- distTips(liz.4d) d2 <- distTips(liz.4d, method = "nNodes") d3 <- distTips(liz.4d, method = "Abouheif") d4 <- distTips(liz.4d, method = "sumDD") obj <- liz.4d tipData(obj) <- as.matrix(d1) table.phylo4d(obj, center = FALSE, scale = FALSE, cex.sym = 1) ## from phylogeny to proximities prox1 <- proxTips(liz.4d) prox2 <- proxTips(liz.4d, method = "nNodes") prox3 <- proxTips(liz.4d, method = "Abouheif") prox4 <- proxTips(liz.4d, method = "sumDD") obj <- liz.4d tipData(obj) <- as.data.frame(prox3) table.phylo4d(obj, center = FALSE, scale = FALSE, cex.sym = 1) ## compute Moran's Index of phylogenetic autocorrelation moran.idx(tipData(liz.4d)[,2], prox3, addInfo=TRUE) abouheif.moran(tipData(liz.4d)[,1:2], prox3) ## orthonormal basis based on partition dummy <- phylo4d(as(liz.4d, "phylo4"), treePart(liz.4d, result="dummy")) table.phylo4d(dummy, scale = FALSE, center = FALSE) dummy.ortho <- phylo4d(as(liz.4d, "phylo4"), treePart(liz.4d, result="orthobasis")) table.phylo4d(dummy.ortho, scale = FALSE, center = FALSE) ## compute Moran's eigenvectors mev <- me.phylo(prox = prox3) obj <- liz.4d tipData(obj) <- (mev) table.phylo4d(obj) ## decompose the variance of the traits corTab <- cor(tipData(liz.4d), mev)^2 barplot(corTab[1,], ylab = expression(r^2)) barplot(corTab, beside=TRUE, col=grey(seq(1, 0, length = 7)), las=3, ylab=expression(r^2)) legend("topright", fill=grey(seq(1, 0, length = 7)), leg=names(tipData(liz.4d)), title="traits") ## test the orthogram ortho4 <- orthogram(tipData(liz.4d)[,4], orthobas = mev) ## multivariate data pca2 <- dudi.pca(tipData(liz.4d), scannf = FALSE) abouheif.moran(pca2$li[,1], prox3) orthogram(pca2$li[,1], orthobas = mev) ## pPCA liz.ppca <- ppca(liz.4d, prox = prox3) liz.ppca barplot(liz.ppca$eig) summary(liz.ppca) scatter(liz.ppca) par(mar=rep(.1,4)) plot(liz.ppca,ratio.tree=.7) dotchart(liz.ppca$c1[,1],lab=rownames(liz.ppca$c1), main="Global principal component 1") abline(v=0,lty=2) dotchart(liz.ppca$c1[,2],lab=rownames(liz.ppca$c1), main="Local principal component 1") abline(v=0,lty=2) s.arrow(liz.ppca$c1,xlim=c(-1,1),clab=1.3,cgrid=1.3)