1 ## SlowinskiGuyer.R (2011-05-05)
3 ## Tests of Diversification Shifts with Sister-Clades
5 ## Copyright 2011 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 slowinskiguyer.test <- function(x, detail = FALSE)
15 chi <- -2 * sum(log(pp))
16 df <- as.integer(2 * length(pp))
17 pval <- pchisq(chi, df, lower.tail = FALSE)
18 res <- data.frame("chisq" = chi, "df" = df, "P.val" = pval, row.names = "")
20 res <- list(res, individual_Pvalues = pp)
24 mcconwaysims.test <- function(x)
27 f <- function(x) ifelse(x == 0, 0, x * log(x))
30 1.629*(f(n1 - 1) - f(n1) + f(n2 - 1) - f(n2) - f(2) -
31 f(n1 + n2 - 2) + f(n1 + n2))
33 chi <- sum(apply(x, 1, LRTp))
34 pval <- pchisq(chi, df <- nrow(x), lower.tail = FALSE)
35 data.frame("chisq" = chi, "df" = df, "P.val" = pval, row.names = "")
38 richness.yule.test <- function(x, t)
45 ## taken from /home/paradis/data/Projets/FISHLOSS/BDfit.R:
46 .PrNt.Yule <- function(N, age, birth) {
47 tmp <- exp(-birth * age)
48 tmp * (1 - tmp)^(N - 1)
51 minusloglik0 <- function(l)
52 -sum(log(.PrNt.Yule(n, tb, l)))
54 minusloglika <- function(l)
55 -sum(log(.PrNt.Yule(n1, t, l[1]))) - sum(log(.PrNt.Yule(n2, t, l[2])))
57 out0 <- nlminb(0.01, minusloglik0, lower = 0, upper = 1)
58 outa <- nlminb(rep(out0$par, 2), minusloglika,
59 lower = c(0, 0), upper = c(1, 1))
60 chi <- 2 * (out0$objective - outa$objective)
61 pval <- pchisq(chi, 1, lower.tail = FALSE)
62 data.frame("chisq" = chi, "df" = 1, "P.val" = pval, row.names = "")
65 diversity.contrast.test <-
66 function(x, method = "ratiolog", alternative = "two.sided", nrep = 0, ...)
68 method <- match.arg(method, c("ratiolog", "proportion",
69 "difference", "logratio"))
70 alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
72 minmax <- t(apply(x, 1, sort)) # sort all rows
73 DIFF <- x[, 1] - x[, 2]
76 CONTRAST <- switch(method,
79 minmax <- minmax + 1 # prevent division by 0
80 ## Note: if min = max, no need to set the contrast
81 ## to zero since this is done with sign()
82 log(minmax[, 2]) / log(minmax[, 1])
84 "proportion" = minmax[, 2] / (minmax[, 2] + minmax[, 1]),
85 "difference" = abs(DIFF),
86 "logratio" = log(minmax[, 1] / minmax[, 2]))
88 y <- SIGN * CONTRAST # the signed contrasts
94 sum(sample(c(-1, 1), size = n, replace = TRUE) * CONTRAST))
95 cases <- switch(alternative,
96 "two.sided" = sum(abs(RND) > sum(y)),
97 "less" = sum(RND < sum(y)),
98 "greater" = sum(RND > sum(y)))
100 } else wilcox.test(x = y, alternative = alternative, ...)$p.value