## mantel test for all data in the Ian dataset
## Andrew Hipp, 18 june 2009
## ahipp@mortonarb.org

ianMantel <- function(phy, traits, perms = 1000, tail = 1, method = "REGRESSION") {
  require(vegan)
  traits <- normalizeData(traits)
  if(class(phy) == "phylo") phy <- list(phy)
  phyD <- lapply(phy, cophenetic)
  numTaxa <- dim(traits)[1]
  traitD <- as.matrix(vegdist(traits, 'euclidean'))
  if(!identical(row.names(phyD[[1]]), row.names(traitD))) stop("Taxon names differ at least in order between matrices")
  traitD <- as.dist(traitD)
  stat <- numeric(length(phy))
  p.value <- numeric(length(phy))
  if(method == "REGRESSION" && tail == 2) {
    tail <- 1
    warning("1-tailed test imposed, as two-tailed test is inappropriate with an r.squared test statistic")
    }
  for(i in seq(length(phy))) {
    print(paste('Doing tree', i))
    if(method == "REGRESSION") stat[i] <- summary(lm(as.dist(phyD[[i]]) ~ traitD))$r.squared
    else stat[i] <- cor.test(as.dist(phyD[[i]]), traitD, method = "pearson")$estimate
    statDist <- numeric(perms)
    for (j in seq(perms)) {
      index <- sample(seq(numTaxa))
      phyDperm <- as.dist(phyD[[i]][index,index])
      if (method == "REGRESSION") statDist[j] <- summary(lm(phyDperm ~ traitD))$r.squared
      else statDist[j] <- cor.test(phyDperm, traitD, method = "pearson")$estimate
      }
    if (tail == 2) p.value[i] <- (2/perms) * min(sum(statDist >= stat[i]), sum(statDist <= stat[i]))
    else p.value[i] <- sum(statDist >= stat[i]) / perms
    }
  mantel.r <- paste(mean(stat), "+/-", (sd(stat) / length(phy)))
  mantel.p <- paste(mean(p.value), "+/-", (sd(p.value) / length(phy)))
  outdata <- list(stat = mantel.r, p = mantel.p)
  }

ianCor.test <- function(phy, traits, native = "lobata", perms = 1000, tail = 1, method = "REGRESSION") {
  traits <- normalizeData(traits)
  noNat <- row.names(traits)[row.names(traits) != native]
  if(class(phy) == "phylo") phy <- list(phy)
  phyD <- lapply(phy, function(x, nat) {cophenetic(x)[nat, noNat]}, nat = native)
  traitD <- as.matrix(vegdist(traits, 'euclidean'))[native, noNat]
  if(!identical(names(phyD[[1]]), names(traitD))) stop("Taxon names differ at least in order between vectors")
  stat <- numeric(length(phy))
  p.value <- numeric(length(phy))
  if(method == "REGRESSION" && tail == 2) {
    tail <- 1
    warning("1-tailed test imposed, as two-tailed test is inappropriate with an r.squared test statistic")
    }
  else if(tail == 1 && method != "REGRESSION") warning("On a correlation test, a 2-tailed p-value is advised")
  for(i in seq(length(phy))) {
    print(paste('Doing tree', i))
    if(method == "REGRESSION") stat[i] <- summary(lm(phyD[[i]] ~ traitD))$r.squared
    else stat[i] <- cor.test(phyD[[i]], traitD, method = "pearson")$estimate
    statDist <- numeric(perms)
    for(j in seq(perms)) {
      phyDperm <- sample(phyD[[i]])
      if (method == "REGRESSION") statDist[j] <- summary(lm(phyDperm ~ traitD))$r.squared
      else statDist[j] <- cor.test(phyDperm, traitD, method = "pearson")$estimate
      }
    if (tail == 2) p.value[i] <- (2/perms) * min(sum(statDist >= stat[i]), sum(statDist <= stat[i]))
    else p.value[i] <- sum(statDist >= stat[i]) / perms
    }
  cor.r <- paste(mean(stat), "+/-", (sd(stat) / length(phy)))
  cor.p <- paste(mean(p.value), "+/-", (sd(p.value) / length(phy)))
  outdata <- list(stat = cor.r, p = cor.p)
  }

ianRegTest <- function(phy, traits, native = "lobata", perms = 1000) {
  traits <- normalizeData(traits)
  if(class(phy) == "phylo") phy <- list(phy)
  phyD <- lapply(phy, function(x, nat) {cophenetic(x)[nat, ]}, nat = native)
  traitD <- as.matrix(vegdist(traits, 'euclidean'))[native, ]
  if(!identical(names(phyD[[1]]), names(traitD))) stop("Taxon names differ at least in order between vectors")
  stat <- numeric(length(phy))
  p.value <- numeric(length(phy))
  for(i in seq(length(phy))) {
    print(paste('Doing tree', i))
    stat[i] <- summary(lm(phyD[[i]] ~ traitD))$r.squared
    statDist <- numeric(perms)
    for(j in seq(perms)) {
      phyDperm <- sample(phyD[[i]])
      statDist[j] <- summary(lm(phyDperm ~ traitD))$r.squared
      }
    p.value[i] <- sum(statDist >= stat[i]) / perms
    }
  cor.r <- paste(mean(stat), "+/-", (sd(stat) / length(phy)))
  cor.p <- paste(mean(p.value), "+/-", (sd(p.value) / length(phy)))
  outdata <- list(stat = cor.r, p = cor.p)
  }