R/cherry.R

   1 ## cherry.R (2009-05-10)
   2
   3 ##     Number of Cherries and Null Models of Trees
   4
   5 ## Copyright 2002-2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 cherry <- function(phy)
  11 {
  12     if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"')
  13     n <- length(phy$tip.label)
  14     nb.node <- phy$Nnode
  15     if (nb.node != n - 1) stop('"phy" is not fully dichotomous')
  16     if (n < 4) stop("not enough tips in your phylogeny for this analysis")
  17     cherry <- sum(tabulate(phy$edge[, 1][phy$edge[, 2] <= n]) == 2)
  18     small.n <- n < 20
  19     if (small.n) {
  20         P.yule <- f.cherry.yule(n, cherry)
  21         P.uniform <- f.cherry.uniform(n, cherry)
  22     }
  23     else {
  24         P.yule <- 2*(1 - pnorm(abs(cherry - n/3)/sqrt(2*n/45)))
  25         mu.unif <- n*(n - 1)/(2*(2*n - 5))
  26         sigma2.unif <- n*(n - 1)*(n - 4)*(n - 5)/(2*(2*n - 5)^2 * (2*n -7))
  27         P.uniform <- 2*(1 - pnorm(abs(cherry - mu.unif)/sqrt(sigma2.unif)))
  28     }
  29     cat("\nAnalysis of the Number of Cherries in a Tree\n\n")
  30     cat("Phylogenetic tree:", deparse(substitute(phy)), "\n")
  31     cat("Number of tips:", n, "\n")
  32     cat("Number of cherries:", cherry, "\n\n")
  33     cat("Null hypothesis: Yule model\n")
  34     cat("    P-value =", round(P.yule, 4), "\n\n")
  35     cat("Null hypothesis: uniform model\n")
  36     cat("    P-value =", round(P.uniform, 4), "\n\n")
  37     if (!small.n) cat("(P-values were computed using normal approximations)\n")
  38 }
  39
  40 f.cherry.yule <- function(n, k)
  41 {
  42     if (k == 0 || k > floor(n/2)) 0 else if (n == 4) if (k == 1) 2/3 else if (k == 2) 1/3 else 0
  43     else (1 - 2*(k - 1)/(n - 1))*f.cherry.yule(n - 1, k - 1) +
  44         2*k/(n - 1)*f.cherry.yule(n - 1, k)
  45 }
  46
  47 f.cherry.uniform <- function(n, k)
  48 {
  49     if (k == 0 || k > floor(n/2)) 0 else if (n == 4) if (k == 1) 4/5 else if (k == 2) 1/5 else 0
  50     else if (k == 1) 0 else (gamma(n + 1)*gamma(n - 2 + 1)*gamma(n - 4 + 1) * 2^(n-2*k)) /
  51         (gamma(n - 2*k + 1)*gamma(2*n - 4 + 1)*gamma(k + 1)*gamma(k - 2 + 1))
  52 }