R/diversi.gof.R

   1 ## diversi.gof.R (2006-10-16)
   2
   3 ##   Tests of Constant Diversification Rates
   4
   5 ## Copyright 2002-2006 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 diversi.gof <- function(x, null = "exponential", z = NULL)
  11 {
  12     n <- length(x)
  13     if (null == "exponential") {
  14         delta <- n/sum(x)
  15         z <- 1 - exp(-delta * sort(x))
  16     }
  17     else {
  18         nmsz <- deparse(substitute(z))
  19         z <- sort(z) # utile ???
  20     }
  21     i <- 1:n
  22     W2 <- sum((z - (2*i - 1)/(2*n))^2) + 1/12*n
  23     A2 <- -sum((2*i - 1)*(log(z) + log(1 - rev(z))))/n - n
  24     if (null == "exponential") {
  25         W2 <- W2*(1 - 0.16/n)
  26         A2 <- A2*(1 + 0.6/n)
  27     }
  28     else W2 <- (W2 - 0.4/n + 0.6/n^2)/(1 + 1/n)
  29     cat("\nTests of Constant Diversification Rates\n\n")
  30     cat("Data:", deparse(substitute(x)), "\n")
  31     cat("Number of branching times:", n, "\n")
  32     cat("Null model: ")
  33     if (null == "exponential") cat("exponential\n\n")
  34     else cat(nmsz, "(user-specified)\n\n")
  35     cat("Cramer-von Mises test: W2 =", round(W2, 3))
  36     if (null == "exponential") {
  37         if (W2 < 0.177) cat("   P > 0.1\n")
  38         if (W2 >= 0.177 && W2 < 0.224) cat("   0.05 < P < 0.1\n")
  39         if (W2 >= 0.224 && W2 < 0.273) cat("   0.025 < P < 0.05\n")
  40         if (W2 >= 0.273 && W2 < 0.337) cat("   0.01 < P < 0.025\n")
  41         if (W2 > 0.337) cat("   P < 0.01\n")
  42     }
  43     else {
  44         if (W2 < 0.347) cat("   P > 0.1\n")
  45         if (W2 >= 0.347 && W2 < 0.461) cat("   0.05 < P < 0.1\n")
  46         if (W2 >= 0.461 && W2 < 0.581) cat("   0.025 < P < 0.05\n")
  47         if (W2 >= 0.581 && W2 < 0.743) cat("   0.01 < P < 0.025\n")
  48         if (W2 > 0.743) cat("   P < 0.01\n")
  49     }
  50     cat("Anderson-Darling test: A2 =", round(A2, 3))
  51     if (null == "exponential") {
  52         if (A2 < 1.078) cat("   P > 0.1\n")
  53         if (A2 >= 1.078 && A2 < 1.341) cat("   0.05 < P < 0.1\n")
  54         if (A2 >= 1.341 && A2 < 1.606) cat("   0.025 < P < 0.05\n")
  55         if (A2 >= 1.606 && A2 < 1.957) cat("   0.01 < P < 0.025\n")
  56         if (A2 > 1.957) cat("   P < 0.01\n")
  57     }
  58     else {
  59         if (A2 < 1.933) cat("   P > 0.1\n")
  60         if (A2 >= 1.933 && A2 < 2.492) cat("   0.05 < P < 0.1\n")
  61         if (A2 >= 2.492 && A2 < 3.070) cat("   0.025 < P < 0.05\n")
  62         if (A2 >= 3.070 && A2 < 3.857) cat("   0.01 < P < 0.025\n")
  63         if (A2 > 3.857) cat("   P < 0.01\n")
  64     }
  65 }