From b0548f80b3ac1d2035ddc360d3366eab4f08d247 Mon Sep 17 00:00:00 2001 From: paradis Date: Thu, 17 Mar 2011 03:00:02 +0000 Subject: [PATCH] final packing for ape 2.7 git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@151 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 6 ++- DESCRIPTION | 4 +- R/SlowinskiGuyer.R | 99 ++++++++++++++++++++++++++++++++++ data/bird.families.R | 2 +- data/bird.orders.R | 2 +- data/cynipids.R | 2 +- data/woodmouse.R | 2 +- man/diversity.contrast.test.Rd | 81 ++++++++++++++++++++++++++++ man/mcconwaysims.test.Rd | 52 ++++++++++++++++++ man/richness.yule.test.Rd | 38 +++++++++++++ man/slowinskiguyer.test.Rd | 64 ++++++++++++++++++++++ 11 files changed, 345 insertions(+), 7 deletions(-) create mode 100644 R/SlowinskiGuyer.R create mode 100644 man/diversity.contrast.test.Rd create mode 100644 man/mcconwaysims.test.Rd create mode 100644 man/richness.yule.test.Rd create mode 100644 man/slowinskiguyer.test.Rd diff --git a/ChangeLog b/ChangeLog index a635615..eb386db 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ - CHANGES IN APE VERSION 2.6-4 + CHANGES IN APE VERSION 2.7 NEW FEATURES @@ -14,6 +14,10 @@ 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). diff --git a/DESCRIPTION b/DESCRIPTION index 97b7a0b..3646a52 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ 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 diff --git a/R/SlowinskiGuyer.R b/R/SlowinskiGuyer.R new file mode 100644 index 0000000..1563f02 --- /dev/null +++ b/R/SlowinskiGuyer.R @@ -0,0 +1,99 @@ +## 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 +} diff --git a/data/bird.families.R b/data/bird.families.R index a806ea3..5ec39df 100644 --- a/data/bird.families.R +++ b/data/bird.families.R @@ -1,2 +1,2 @@ -require(ape, quietly = TRUE, save = FALSE) +require(ape, quietly = TRUE) bird.families <- read.tree("bird.families.tre") diff --git a/data/bird.orders.R b/data/bird.orders.R index f1640e1..260a5e7 100644 --- a/data/bird.orders.R +++ b/data/bird.orders.R @@ -1,2 +1,2 @@ -require(ape, quietly = TRUE, save = FALSE) +require(ape, quietly = TRUE) bird.orders <- read.tree("bird.orders.tre") diff --git a/data/cynipids.R b/data/cynipids.R index 649f9a9..f5691b5 100644 --- a/data/cynipids.R +++ b/data/cynipids.R @@ -1,3 +1,3 @@ -require(ape, quietly = TRUE, save = FALSE) +require(ape, quietly = TRUE) cynipids <- read.nexus.data("cynipids.txt") diff --git a/data/woodmouse.R b/data/woodmouse.R index 716be9a..9f4c397 100644 --- a/data/woodmouse.R +++ b/data/woodmouse.R @@ -1,2 +1,2 @@ -require(ape, quietly = TRUE, save = FALSE) +require(ape, quietly = TRUE) woodmouse <- read.dna("woodmouse.txt", format = "sequential") diff --git a/man/diversity.contrast.test.Rd b/man/diversity.contrast.test.Rd new file mode 100644 index 0000000..b24d42b --- /dev/null +++ b/man/diversity.contrast.test.Rd @@ -0,0 +1,81 @@ +\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} diff --git a/man/mcconwaysims.test.Rd b/man/mcconwaysims.test.Rd new file mode 100644 index 0000000..1fb3c6b --- /dev/null +++ b/man/mcconwaysims.test.Rd @@ -0,0 +1,52 @@ +\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} diff --git a/man/richness.yule.test.Rd b/man/richness.yule.test.Rd new file mode 100644 index 0000000..f71200d --- /dev/null +++ b/man/richness.yule.test.Rd @@ -0,0 +1,38 @@ +\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} diff --git a/man/slowinskiguyer.test.Rd b/man/slowinskiguyer.test.Rd new file mode 100644 index 0000000..67877e2 --- /dev/null +++ b/man/slowinskiguyer.test.Rd @@ -0,0 +1,64 @@ +\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} -- 2.39.2