# Original R script provided by the authors of the following paper as supplementary material # Harmon, L.J & Glor, R.E. (2010) # Poor statistical performance of the Mantel test in phylogenetic comparative analyses # Evolution 64(7): 2173-2178 ################################################### ## Modified by O.Missa to suit different needs ## ## (working from a genetic distance Matrix ## ## rather than a phylogenetic tree) ## ################################################### ############################################################ ## Calculates probabilities for phylogenetic permutations ## ## from Lapointe and Garland 2001 ## phyloProb<-function(gen.dist, k=1) { pd <- gen.dist pdr<-pd/max(pd) s<-k-pdr p<-s/rowSums(s) p } ## (probs <- phyloProb(virus.dist)) ## just to test it works OK ######################################################### ## Permutes species according to phylogenetic distance ## ## Returns new permuted order ## phyloPermute<-function(gen.dist, k=1) { p<-phyloProb(gen.dist, k) tt<-rownames(p) nsp<-length(tt) order<-sample(1:nsp, replace=F) ttnew<-character(nsp) cpm<-p for(j in order[-nsp]) { cpm<-cpm/rowSums(cpm) rr<-which(rownames(cpm)==tt[j]) pp<-cpm[rr,] s2<-sample(names(pp), size=1, replace=T, prob=pp) slot<-which(tt==s2) rc<-which(colnames(cpm)==s2) ttnew[slot]<-tt[j] cpm<-cpm[-rr,-rc] } ttnew[which(ttnew=="")]<-tt[order[nsp]] ttnew trans <- match(ttnew, tt) trans } ##################################################### ## Adapting the mantel function from vegan package ## pp.mantel.f <- function (xdis, ydis, phylodist, k=1, method = "pearson", permutations = 9999) { xdis <- as.dist(xdis) ydis <- as.vector(as.dist(ydis)) tmp <- cor.test(as.vector(xdis), ydis, method = method) statistic <- as.numeric(tmp$estimate) variant <- tmp$method if (permutations) { N <- attributes(xdis)$Size perm <- rep(0, permutations) for (i in 1:permutations) { take <- phyloPermute(as.matrix(phylodist),k) permvec <- as.vector(as.dist(as.matrix(xdis)[take, take])) perm[i] <- cor(permvec, ydis, method = method) } signif <- (sum(perm >= statistic) + 1)/(permutations + 1) } else { signif <- NA perm <- NULL } res <- list(call = match.call(), method = variant, statistic = statistic, signif = signif, perm = perm, permutations = permutations) class(res) <- "mantel" res } ############################################################# ## Adapting the mantel.partial function from vegan package ## pp.mantel.partial.f <- function (xdis, ydis, zdis, phylodist, k=1, method = "pearson", permutations = 999) { part.cor <- function(rxy, rxz, ryz) { (rxy - rxz * ryz)/sqrt(1 - rxz * rxz)/sqrt(1 - ryz * ryz) } xdis <- as.dist(xdis) ydis <- as.vector(as.dist(ydis)) zdis <- as.vector(as.dist(zdis)) rxy <- cor.test(as.vector(xdis), ydis, method = method) rxz <- cor(as.vector(xdis), zdis, method = method) ryz <- cor(ydis, zdis, method = method) variant <- rxy$method rxy <- rxy$estimate statistic <- part.cor(rxy, rxz, ryz) if (permutations) { N <- attributes(xdis)$Size perm <- rep(0, permutations) for (i in 1:permutations) { take <- phyloPermute(as.matrix(phylodist),k) permvec <- as.vector(as.dist(as.matrix(xdis)[take, take])) rxy <- cor(permvec, ydis, method = method) rxz <- cor(permvec, zdis, method = method) perm[i] <- part.cor(rxy, rxz, ryz) } signif <- (sum(perm >= statistic) + 1)/(permutations + 1) } else { signif <- NA perm <- NULL } res <- list(call = match.call(), method = variant, statistic = statistic, signif = signif, perm = perm, permutations = permutations) class(res) <- c("mantel.partial", "mantel") res }