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