- CHANGES IN APE VERSION 2.6-4
+ CHANGES IN APE VERSION 2.7
NEW FEATURES
nucleotide sequence alignment by calling the external programs
of the same names.
+ o Four new functions, diversity.contrast.test, mcconwaysims.test,
+ richness.yule.test, and slowinskiguyer.test, implement various
+ tests of diversification shifts using sister-clade comparisons.
+
o base.freq() gains an option 'all' to count all the possible bases
including the ambiguous ones (defaults to FALSE).
Package: ape
-Version: 2.6-4
-Date: 2011-03-16
+Version: 2.7
+Date: 2011-03-17
Title: Analyses of Phylogenetics and Evolution
Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Klaus Schliep, Korbinian Strimmer, Damien de Vienne
Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
--- /dev/null
+## SlowinskiGuyer.R (2011-03-10)
+
+## Tests of Diversification Shifts with Sister-Clades
+
+## Copyright 2011 Emmanuel Paradis
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+slowinskiguyer.test <- function(x, detail = FALSE)
+{
+ r <- x[, 1]
+ n <- x[, 1] + x[, 2]
+ pp <- (n - r)/(n - 1)
+ chi <- -2 * sum(log(pp))
+ df <- as.integer(2 * length(pp))
+ pval <- pchisq(chi, df, lower.tail = FALSE)
+ res <- data.frame("chisq" = chi, "df" = df, "P.val" = pval, row.names = "")
+ if (detail)
+ res <- list(res, individual_Pvalues = pp)
+ res
+}
+
+mcconwaysims.test <- function(x)
+{
+ LRTp <- function(x) {
+ f <- function(x) ifelse(x == 0, 0, x * log(x))
+ n1 <- x[1]
+ n2 <- x[2]
+ 1.629*(f(n1 - 1) - f(n1) + f(n2 - 1) - f(n2) - f(2) -
+ f(n1 + n2 - 2) + f(n1 + n2))
+ }
+ chi <- sum(apply(x, 1, LRTp))
+ pval <- pchisq(chi, df <- nrow(x), lower.tail = FALSE)
+ data.frame("chisq" = chi, "df" = df, "P.val" = pval, row.names = "")
+}
+
+richness.yule.test <- function(x, t)
+{
+ n1 <- x[, 1]
+ n2 <- x[, 2]
+ n <- c(n1, n2)
+ tb <- c(t, t)
+
+ ## taken from /home/paradis/data/Projets/FISHLOSS/BDfit.R:
+ .PrNt.Yule <- function(N, age, birth) {
+ tmp <- exp(-birth * age)
+ tmp * (1 - tmp)^(N - 1)
+ }
+
+ minusloglik0 <- function(l)
+ -sum(log(.PrNt.Yule(n, tb, l)))
+
+ minusloglika <- function(l)
+ -sum(log(.PrNt.Yule(n1, t, l[1]))) - sum(log(.PrNt.Yule(n2, t, l[2])))
+
+ out0 <- nlminb(0.01, minusloglik0, lower = 0, upper = 1)
+ outa <- nlminb(rep(out0$par, 2), minusloglika,
+ lower = c(0, 0), upper = c(1, 1))
+ chi <- 2 * (out0$objective - outa$objective)
+ pval <- pchisq(chi, 1, lower.tail = FALSE)
+ data.frame("chisq" = chi, "df" = 1, "P.val" = pval, row.names = "")
+}
+
+diversity.contrast.test <-
+ function(x, method = "ratiolog", alternative = "two.sided", nrep = 0, ...)
+{
+ method <- match.arg(method, c("ratiolog", "proportion", "difference"))
+ alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
+
+ minmax <- t(apply(x, 1, sort)) # sort all rows
+ DIFF <- x[, 1] - x[, 2]
+ SIGN <- sign(DIFF)
+
+ CONTRAST <- switch(method,
+ "ratiolog" = {
+ if (any(minmax == 1))
+ minmax <- minmax + 1 # prevent division by 0
+ ## Note: if min = max, no need to set the contrast
+ ## to zero since this is done with sign()
+ log(minmax[, 2]) / log(minmax[, 1])
+ },
+ "proportion" = minmax[, 2] / (minmax[, 2] + minmax[, 1]),
+ "difference" = abs(DIFF))
+
+ y <- SIGN * CONTRAST # the signed contrasts
+
+ if (nrep) {
+ n <- length(SIGN)
+ RND <-
+ replicate(nrep,
+ sum(sample(c(-1, 1), size = n, replace = TRUE) * CONTRAST))
+ cases <- switch(alternative,
+ "two.sided" = sum(abs(RND) > sum(y)),
+ "less" = sum(RND < sum(y)),
+ "greater" = sum(RND > sum(y)))
+ cases/nrep
+ } else wilcox.test(x = y, alternative = alternative, ...)$p.value
+}
-require(ape, quietly = TRUE, save = FALSE)
+require(ape, quietly = TRUE)
bird.families <- read.tree("bird.families.tre")
-require(ape, quietly = TRUE, save = FALSE)
+require(ape, quietly = TRUE)
bird.orders <- read.tree("bird.orders.tre")
-require(ape, quietly = TRUE, save = FALSE)
+require(ape, quietly = TRUE)
cynipids <- read.nexus.data("cynipids.txt")
-require(ape, quietly = TRUE, save = FALSE)
+require(ape, quietly = TRUE)
woodmouse <- read.dna("woodmouse.txt", format = "sequential")
--- /dev/null
+\name{diversity.contrast.test}
+\alias{diversity.contrast.test}
+\title{Diversity Contrast Test}
+\description{
+ This function performs the diversity contrast test comparing pairs of
+ sister-clades.
+}
+\usage{
+diversity.contrast.test(x, method = "ratiolog",
+ alternative = "two.sided", nrep = 0, ...)
+}
+\arguments{
+ \item{x}{a matrix or a data frame with at least two columns: the first
+ one gives the number of species in clades with a trait supposed to
+ increase or decrease diversification rate, and the second one the number of
+ species in the sister-clades without the trait. Each
+ row represents a pair of sister-clades.}
+ \item{method}{a character string specifying the kind of test: \code{"ratiolog"} (default),
+ \code{"proportion"}, \code{"difference"}, or any unambiguous
+ abbreviation of these.}
+ \item{alternative}{a character string defining the alternative
+ hypothesis: \code{"two.sided"} (default), \code{"less"},
+ \code{"greater"}, or any unambiguous abbreviation of these.}
+ \item{nrep}{the number of replications of the randomization test; by
+ default, a Wilcoxon test is done.}
+ \item{\dots}{arguments passed to the function \code{\link[stats]{wilcox.test}}.}
+}
+\details{
+ If \code{method = "ratiolog"}, the test described in Barraclough et
+ al. (1996) is performed. If \code{method = "proportion"}, the version
+ in Barraclough et al. (1995) is used. If \code{method = "difference"},
+ then this is Wiegmann et al.'s (1993) version. Vamosi and Vamosi (2005)
+ gave a detailed account of these three tests which are essentially
+ different versions of the same test.
+
+ If \code{nrep = 0}, a Wilcoxon test is done on the species diversity
+ contrasts with the null hypothesis is that they are distributed around
+ zero. If \code{nrep > 0}, a randomization procedure is done where the
+ signs of the diversity contrasts are randomly chosen. This is used to
+ create a distribution of the test statistic which is compared with the
+ observed value (the sum of the diversity contrasts).
+}
+\value{
+ a single numeric value with the \emph{P}-value.
+}
+\references{
+ Barraclough, T. G., Harvey, P. H. and Nee, S. (1995) Sexual
+ selection and taxonomic diversity in passerine birds.
+ \emph{Proceedings of the Royal Society of London. Series B. Biological
+ Sciences}, \bold{259}, 211--215.
+
+ Barraclough, T. G., Harvey, P. H., and Nee, S. (1996) Rate of
+ \emph{rbc}L gene sequence evolution and species diversification in
+ flowering plants (angiosperms). \emph{Proceedings of the Royal Society
+ of London. Series B. Biological Sciences}, \bold{263}, 589--591.
+
+ Vamosi, S. M. and Vamosi, J. C. (2005) Endless tests: guidelines for
+ analysing non-nested sister-group comparisons. \emph{Evolutionary
+ Ecology Research}, \bold{7}, 567--579.
+
+ Wiegmann, B., Mitter, C. and Farrell, B. 1993. Diversification of
+ carnivorous parasitic insects: extraordinary radiation or specialized
+ dead end? \emph{American Naturalist}, \bold{142}, 737--754.
+}
+\author{Emmanuel Paradis}
+\seealso{
+ \code{\link{slowinskiguyer.test}}, \code{\link{mcconwaysims.test}}
+ \code{\link{richness.yule.test}}
+}
+\examples{
+### data from Vamosi & Vamosi (2005):
+fleshy <- c(1, 1, 1, 1, 1, 3, 3, 5, 9, 16, 33, 40, 50, 100, 216, 393, 850, 947,1700)
+dry <- c(2, 64, 300, 89, 67, 4, 34, 10, 150, 35, 2, 60, 81, 1, 3, 1, 11, 1, 18)
+x <- cbind(fleshy, dry)
+diversity.contrast.test(x)
+diversity.contrast.test(x, alt = "g")
+diversity.contrast.test(x, alt = "g", nrep = 1e4)
+slowinskiguyer.test(x)
+mcconwaysims.test(x)
+}
+\keyword{htest}
--- /dev/null
+\name{mcconwaysims.test}
+\alias{mcconwaysims.test}
+\title{McConway-Sims Test of Homogeneous Diversification}
+\description{
+ This function performs the McConway--Sims test that a trait or
+ variable does not affect diversification rate.
+}
+\usage{
+mcconwaysims.test(x)
+}
+\arguments{
+ \item{x}{a matrix or a data frame with at least two columns: the first
+ one gives the number of species in clades with a trait supposed to
+ increase or decrease diversification rate, and the second one the number of
+ species in the sister-clades without the trait. Each
+ row represents a pair of sister-clades.}
+}
+\details{
+ The McConway--Sims test compares a series of sister-clades where one
+ of the two is characterized by a trait supposed to affect
+ diversification rate. The null hypothesis is that the trait does not
+ affect diversification. The alternative hypothesis is that
+ diversification rate is increased or decreased by the trait (by
+ contrast to the Slowinski--Guyer test). The test is a likelihood-ratio
+ of a null Yule model and an alternative model with two parameters.
+}
+\value{
+ a data frame with the \eqn{\chi^2}{chi2}, the number of degrees of
+ freedom, and the \emph{P}-value.
+}
+\references{
+ McConway, K. J. and Sims, H. J. (2004) A likelihood-based method for
+ testing for nonstochastic variation of diversification rates in
+ phylogenies. \emph{Evolution}, \bold{58}, 12--23.
+}
+\author{Emmanuel Paradis}
+\seealso{
+ \code{\link{balance}}, \code{\link{slowinskiguyer.test}},
+ \code{\link[geiger]{rc}} in \pkg{geiger},
+ \code{\link[apTreeshape]{shift.test}} in \pkg{apTreeshape}
+}
+\examples{
+### simulate 10 clades with lambda = 0.1 and mu = 0.09:
+n0 <- replicate(10, balance(rbdtree(.1, .09, Tmax = 35))[1])
+### simulate 10 clades with lambda = 0.15 and mu = 0.1:
+n1 <- replicate(10, balance(rbdtree(.15, .1, Tmax = 35))[1])
+x <- cbind(n1, n0)
+mcconwaysims.test(x)
+slowinskiguyer.test(x)
+richness.yule.test(x, 35)
+}
+\keyword{htest}
--- /dev/null
+\name{richness.yule.test}
+\alias{richness.yule.test}
+\title{}
+\description{
+ This function performs a test of shift in diversification rate using
+ probabilities from the Yule process.
+}
+\usage{
+richness.yule.test(x, t)
+}
+\arguments{
+ \item{x}{a matrix or a data frame with at least two columns: the first
+ one gives the number of species in clades with a trait supposed to
+ increase or decrease diversification rate, and the second one the number of
+ species in the sister-clades without the trait. Each
+ row represents a pair of sister-clades.}
+ \item{t}{a numeric vector giving the divergence times of each pair of
+ clades in \code{x}.}
+}
+\details{
+}
+\value{
+ a data frame with the \eqn{\chi^2}{chi2}, the number of degrees of
+ freedom (= 1), and the \emph{P}-value.
+}
+\references{
+ Paradis, E. Shift in diversification in sister-clade comparisons: a
+ more powerful test. (manuscript submitted)
+}
+\author{Emmanuel Paradis}
+\seealso{
+ \code{\link{slowinskiguyer.test}}, \code{\link{mcconwaysims.test}}
+ \code{\link{diversity.contrast.test}}
+}
+\examples{
+### see examples(mcconwaysims.test)
+}
+\keyword{htest}
--- /dev/null
+\name{slowinskiguyer.test}
+\alias{slowinskiguyer.test}
+\title{Slowinski-Guyer Test of Homogeneous Diversification}
+\description{
+ This function performs the Slowinski--Guyer test that a trait or
+ variable does not increase diversification rate.
+}
+\usage{
+slowinskiguyer.test(x, detail = FALSE)
+}
+\arguments{
+ \item{x}{a matrix or a data frame with at least two columns: the first
+ one gives the number of species in clades with a trait supposed to
+ increase diversification rate, and the second one the number of
+ species in the corresponding sister-clade without the trait. Each
+ row represents a pair of sister-clades.}
+ \item{detail}{if \code{TRUE}, the individual P-values are appended.}
+}
+\details{
+ The Slowinski--Guyer test compares a series of sister-clades where one
+ of the two is characterized by a trait supposed to increase
+ diversification rate. The null hypothesis is that the trait does not
+ affect diversification. If the trait decreased diversification rate,
+ then the null hypothesis cannot be rejected.
+
+ The present function has mainly a historical interest. The
+ Slowinski--Guyer test generally performs poorly: see McConway and Sims
+ (2004) for an alternative and the functions cited below.
+}
+\value{
+ a data frame with the \eqn{\chi^2}{chi2}, the number of degrees of
+ freedom, and the \emph{P}-value. If \code{detail = TRUE}, a list is
+ returned with the data frame and a vector of individual
+ \emph{P}-values for each pair of sister-clades.
+}
+\references{
+ McConway, K. J. and Sims, H. J. (2004) A likelihood-based method for
+ testing for nonstochastic variation of diversification rates in
+ phylogenies. \emph{Evolution}, \bold{58}, 12--23.
+
+ Slowinski, J. B. and Guyer, C. (1993) Testing whether certain traits
+ have caused amplified diversification: an improved method based on a
+ model of random speciation and extinction. \emph{American Naturalist},
+ \bold{142}, 1019--1024.
+}
+\author{Emmanuel Paradis}
+\seealso{
+ \code{\link{balance}}, \code{\link{mcconwaysims.test}},
+ \code{\link{diversity.contrast.test}},
+ \code{\link{richness.yule.test}},
+ \code{\link[geiger]{rc}} in \pkg{geiger},
+ \code{\link[apTreeshape]{shift.test}} in \pkg{apTreeshape}
+}
+\examples{
+### from Table 1 in Slowinski and Guyer(1993):
+viviparous <- c(98, 8, 193, 36, 7, 128, 2, 3, 23, 70)
+oviparous <- c(234, 17, 100, 4, 1, 12, 6, 1, 481, 11)
+x <- data.frame(viviparous, oviparous)
+slowinskiguyer.test(x, TRUE) # 'P ~ 0.32' in the paper
+xalt <- x
+xalt[3, 2] <- 1
+slowinskiguyer.test(xalt)
+}
+\keyword{htest}
projects / ape.git / commitdiff