R/gammaStat.R

   1 ## gammaStat.R (2009-05-10)
   2
   3 ##   Gamma-Statistic of Pybus and Harvey
   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 gammaStat <- function(phy)
  11 {
  12     if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"')
  13     N <- length(phy$tip.label)
  14     bt <- sort(branching.times(phy))
  15     g <- rev(c(bt[1], diff(bt))) # internode intervals are from past to present
  16     ST <- sum((2:N) * g)
  17     stat <- sum(cumsum((2:(N - 1)) * g[-(N - 1)]))/(N - 2)
  18     m <- ST/2
  19     s <- ST * sqrt(1/(12 * (N - 2)))
  20     (stat - m)/s
  21 }