R/SlowinskiGuyer.R

   1 ## SlowinskiGuyer.R (2011-05-05)
   2
   3 ##   Tests of Diversification Shifts with Sister-Clades
   4
   5 ## Copyright 2011 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 slowinskiguyer.test <- function(x, detail = FALSE)
  11 {
  12     r <- x[, 1]
  13     n <- x[, 1] + x[, 2]
  14     pp <- (n - r)/(n - 1)
  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 = "")
  19     if (detail)
  20         res <- list(res, individual_Pvalues = pp)
  21     res
  22 }
  23
  24 mcconwaysims.test <- function(x)
  25 {
  26     LRTp <- function(x) {
  27         f <- function(x) ifelse(x == 0, 0, x * log(x))
  28         n1 <- x[1]
  29         n2 <- x[2]
  30         1.629*(f(n1 - 1) - f(n1) + f(n2 - 1) - f(n2) - f(2) -
  31                f(n1 + n2 - 2) + f(n1 + n2))
  32     }
  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 = "")
  36 }
  37
  38 richness.yule.test <- function(x, t)
  39 {
  40     n1 <- x[, 1]
  41     n2 <- x[, 2]
  42     n <- c(n1, n2)
  43     tb <- c(t, t)
  44
  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)
  49     }
  50
  51     minusloglik0 <- function(l)
  52         -sum(log(.PrNt.Yule(n, tb, l)))
  53
  54     minusloglika <- function(l)
  55         -sum(log(.PrNt.Yule(n1, t, l[1]))) - sum(log(.PrNt.Yule(n2, t, l[2])))
  56
  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 = "")
  63 }
  64
  65 diversity.contrast.test <-
  66     function(x, method = "ratiolog", alternative = "two.sided", nrep = 0, ...)
  67 {
  68     method <- match.arg(method, c("ratiolog", "proportion",
  69                                   "difference", "logratio"))
  70     alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
  71
  72     minmax <- t(apply(x, 1, sort)) # sort all rows
  73     DIFF <- x[, 1] - x[, 2]
  74     SIGN <- sign(DIFF)
  75
  76     CONTRAST <- switch(method,
  77                        "ratiolog" = {
  78                            if (any(minmax == 1))
  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])
  83                        },
  84                        "proportion" = minmax[, 2] / (minmax[, 2] + minmax[, 1]),
  85                        "difference" = abs(DIFF),
  86                        "logratio" = log(minmax[, 1] / minmax[, 2]))
  87
  88     y <- SIGN * CONTRAST # the signed contrasts
  89
  90     if (nrep) {
  91         n <- length(SIGN)
  92         RND <-
  93             replicate(nrep,
  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)))
  99         cases/nrep
 100     } else wilcox.test(x = y, alternative = alternative, ...)$p.value
 101 }