R/sh.test.R

   1 ## sh.test.R (2006-07-06)
   2
   3 ##   Shimodaira-Hasegawa Test
   4
   5 ## Copyright 2006 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 sh.test <- function(..., x, model = DNAmodel(), B = 100)
  11 {
  12     ## Prepare the list of trees:
  13     phy <- list(...)
  14     if (length(phy) == 1 && class(phy[[1]]) != "phylo")
  15       phy <- unlist(phy, recursive = FALSE)
  16     ntree <- length(phy)
  17
  18     ## Arrange the sequences as a matrix:
  19     if (is.list(x)) {
  20         nm <- names(x)
  21         n <- length(x)
  22         x <- unlist(x)
  23         nL <- length(x)
  24         x <- matrix(x, n, nL/n, byrow = TRUE)
  25         rownames(x) <- nm
  26     }
  27
  28     ## Step 1:
  29     foo <- function(PHY)
  30       attr(mlphylo(model, x, PHY, search.tree = FALSE, quiet = TRUE), "loglik")
  31     Talpha <- sapply(phy, foo)
  32     Talpha <- max(Talpha) - Talpha
  33
  34     ## Do the bootstrap resampling (Step 2):
  35     M <- matrix(NA, ntree, B)
  36     for (i in 1:B) {
  37         boot.samp <- x[, sample(ncol(x), replace = TRUE)]
  38         for (j in 1:ntree)
  39           M[j, i] <- attr(mlphylo(model, boot.samp, phy[[j]],
  40                                   search.tree = FALSE, quiet = TRUE),
  41                           "loglik")
  42     }
  43     M <- M - rowMeans(M) # Step 3
  44     ## Step 4: <FIXME> This can greatly simplified </FIXME>
  45     for (i in 1:B)
  46       for (j in 1:ntree)
  47         M[j, i] <- max(M[j, i] - M[, i])
  48     ## Step 5:
  49     count <- numeric(ntree)
  50     for (j in 1:ntree)
  51       count[j] <- sum(M[j, ] > Talpha[j])
  52     count <- count/B
  53     names(count) <- names(phy)
  54     count
  55 }