From 071f848fadb1be486490a10f3de8324d3aaa5cb7 Mon Sep 17 00:00:00 2001 From: paradis Date: Thu, 4 Nov 2010 08:38:08 +0000 Subject: [PATCH] a bunch of modifications and fixes git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@135 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 10 +++++++++ DESCRIPTION | 4 ++-- R/birthdeath.R | 54 +++++++++++++++++++++++++++-------------------- R/compar.ou.R | 7 +++--- R/summary.phylo.R | 3 +-- man/bd.ext.Rd | 24 +++++++++++++++++++-- man/compar.ou.Rd | 6 +++--- man/write.tree.Rd | 10 ++++++--- 8 files changed, 80 insertions(+), 38 deletions(-) diff --git a/ChangeLog b/ChangeLog index 50d75d7..8db603f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,13 @@ + CHANGES IN APE VERSION 2.6-2 + + +NEW FEATURES + + o bd.ext() has a new option conditional = TRUE to use probabilities + conditioned on no extinction for the taxonomic data. + + + CHANGES IN APE VERSION 2.6-1 diff --git a/DESCRIPTION b/DESCRIPTION index e38fc9f..01bf1ad 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape -Version: 2.6-1 -Date: 2010-11-01 +Version: 2.6-2 +Date: 2010-11-04 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/birthdeath.R b/R/birthdeath.R index b4a25e1..3aee631 100644 --- a/R/birthdeath.R +++ b/R/birthdeath.R @@ -1,4 +1,4 @@ -## birthdeath.R (2010-10-17) +## birthdeath.R (2010-11-04) ## Estimation of Speciation and Extinction Rates ## with Birth-Death Models @@ -79,9 +79,10 @@ print.birthdeath <- function(x, ...) cat(" b - d: [", x$CI[2, 1], ", ", x$CI[2, 2], "]", "\n\n", sep = "") } -bd.ext <- function(phy, S) +bd.ext <- function(phy, S, conditional = TRUE) { - if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"') + if (!inherits(phy, "phylo")) + stop('object "phy" is not of class "phylo"') if (!is.null(names(S))) { if (all(names(S) %in% phy$tip.label)) S <- S[phy$tip.label] else warning('the names of argument "S" and the tip labels @@ -91,27 +92,34 @@ did not match: the former were ignored.') x <- branching.times(phy) x <- c(x[1], x) trm.br <- phy$edge.length[phy$edge[, 2] <= N] - dev <- function(a, r) - { - -2 * (lfactorial(N - 1) - + (N - 2) * log(r) - + (3 * N) * log(1 - a) - + 2 * r * sum(x[2:N]) - - 2 * sum(log(exp(r * x[2:N]) - a)) - + r * sum(trm.br) - + sum((S - 1) * log(exp(r * trm.br) - 1)) - - sum((S + 1) * log(exp(r * trm.br) - a))) - } - out <- nlm(function(p) dev(p[1], p[2]), c(0, 0.2), hessian = TRUE) - if (out$estimate[1] < 0) { - out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE) - para <- c(0, out$estimate) - se <- c(0, sqrt(diag(solve(out$hessian)))) - } - else { - para <- out$estimate - se <- sqrt(diag(solve(out$hessian))) + if (conditional) { + dev <- function(a, r) { + if (a >= 1 || a < 0 || r <= 0) return(1e50) + ert <- exp(r * trm.br) + zeta <- (ert - 1)/(ert - a) + -2 * (lfactorial(N - 1) + + (N - 2) * log(r) + + N * log(1 - a) + + 2 * r * sum(x[2:N]) + - 2 * sum(log(exp(r * x[2:N]) - a)) + + sum(log(1 - zeta) + (S - 1)*log(zeta))) + } + } else { + dev <- function(a, r) { + if (a >= 1 || a < 0 || r <= 0) return(1e50) + -2 * (lfactorial(N - 1) + + (N - 2) * log(r) + + (3 * N) * log(1 - a) + + 2 * r * sum(x[2:N]) + - 2 * sum(log(exp(r * x[2:N]) - a)) + + r * sum(trm.br) + + sum((S - 1) * log(exp(r * trm.br) - 1)) + - sum((S + 1) * log(exp(r * trm.br) - a))) + } } + out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE) + para <- out$estimate + se <- sqrt(diag(solve(out$hessian))) Dev <- out$minimum cat("\nExtended Version of the Birth-Death Models to\n") cat(" Estimate Speciation and Extinction Rates\n\n") diff --git a/R/compar.ou.R b/R/compar.ou.R index 7732b31..2e1d3a5 100644 --- a/R/compar.ou.R +++ b/R/compar.ou.R @@ -1,4 +1,4 @@ -## compar.ou.R (2010-03-15) +## compar.ou.R (2010-11-04) ## Ornstein--Uhlenbeck Model for Continuous Characters @@ -53,8 +53,9 @@ compar.ou <- function(x, phy, node = NULL, alpha = NULL) ## fixed a bug below: must be '%*% theta' instead of '* theta' (2010-03-15) M <- rowSums((exp(-alpha * Wend) - exp(-alpha * Wstart)) %*% theta) V <- exp(-alpha * W) * (1 - exp(-2 * alpha * (Tmax - W/2))) - n * log(2 * pi * sigma2) + log(det(V)) + - (t(x - M) %*% chol2inv(V) %*% (x - M)) / sigma2 + R <- chol(V) # correction by Cecile Ane (2010-11-04) + n * log(2 * pi * sigma2) + 2 * sum(log(diag(R))) + + (t(x - M) %*% chol2inv(R) %*% (x - M)) / sigma2 } out <- if (is.null(alpha)) diff --git a/R/summary.phylo.R b/R/summary.phylo.R index bf01266..e44d48d 100644 --- a/R/summary.phylo.R +++ b/R/summary.phylo.R @@ -1,4 +1,4 @@ -## summary.phylo.R (2010-05-25) +## summary.phylo.R (2010-11-03) ## Print Summary of a Phylogeny and "multiPhylo" operators @@ -118,7 +118,6 @@ print.multiPhylo <- function(x, details = FALSE, ...) if (details) for (i in 1:N) cat("tree", i, ":", length(x[[i]]$tip.label), "tips\n") - cat("\n") } "[[.multiPhylo" <- function(x, i) diff --git a/man/bd.ext.Rd b/man/bd.ext.Rd index 2f246c9..f35236d 100644 --- a/man/bd.ext.Rd +++ b/man/bd.ext.Rd @@ -3,11 +3,14 @@ \title{Extended Version of the Birth-Death Models to Estimate Speciation and Extinction Rates} \usage{ -bd.ext(phy, S) +bd.ext(phy, S, conditional = TRUE) } \arguments{ \item{phy}{an object of class \code{"phylo"}.} \item{S}{a numeric vector giving the number of species for each tip.} + \item{conditional}{whether probabilities should be conditioned on no + extinction (mainly to compare results with previous analyses; see + details).} } \description{ This function fits by maximum likelihood a birth-death model to the @@ -32,12 +35,28 @@ bd.ext(phy, S) than the tip labels of \code{phy}, and a warning message is issued. Note that the function does not check that the tree is effectively - ultrametric, so if it is not, the returned result may not be meaningful. + ultrametric, so if it is not, the returned result may not be + meaningful. + + If \code{conditional = TRUE}, the probabilities of the taxonomic data + are calculated conditioned on no extinction (Rabosky et al. 2007). In + previous versions of the present function (until ape 2.6-1), + unconditional probabilities were used resulting in underestimated + extinction rate. Though it does not make much sense to use + \code{conditional = FALSE}, this option is provided to compare results + from previous analyses: if the species richnesses are relatively low, + both versions will give similar results (see examples). } \references{ Paradis, E. (2003) Analysis of diversification: combining phylogenetic and taxonomic data. \emph{Proceedings of the Royal Society of London. Series B. Biological Sciences}, \bold{270}, 2499--2505. + + Rabosky, D. L., Donnellan, S. C., Talaba, A. L. and Lovette, + I. J. (2007) Exceptional among-lineage variation in diversification + rates during the radiation of Australia's most diverse vertebrate + clade. \emph{Proceedings of the Royal Society of London. Series + B. Biological Sciences}, \bold{274}, 2915-2923. } \author{Emmanuel Paradis} \seealso{ @@ -53,5 +72,6 @@ data(bird.orders) S <- c(10, 47, 69, 214, 161, 17, 355, 51, 56, 10, 39, 152, 6, 143, 358, 103, 319, 23, 291, 313, 196, 1027, 5712) bd.ext(bird.orders, S) +bd.ext(bird.orders, S, FALSE) # same than older versions } \keyword{models} diff --git a/man/compar.ou.Rd b/man/compar.ou.Rd index 4c1c02b..84be454 100644 --- a/man/compar.ou.Rd +++ b/man/compar.ou.Rd @@ -86,9 +86,9 @@ compar.ou(x, phy, node = NULL, alpha = NULL) data(bird.orders) ### This is likely to give you estimates close to 0, 1, and 0 ### for alpha, sigma^2, and theta, respectively: -compar.ou(rnorm(23), bird.orders) +compar.ou(x <- rnorm(23), bird.orders) ### Much better with a fixed alpha: -compar.ou(rnorm(23), bird.orders, alpha = 0.1) +compar.ou(x, bird.orders, alpha = 0.1) ### Let us 'mimick' the effect of different optima ### for the two clades of birds... x <- c(rnorm(5, 0), rnorm(18, 5)) @@ -97,6 +97,6 @@ compar.ou(x, bird.orders, node = 25, alpha = .1) ### ... and the model with a single optimum: compar.ou(x, bird.orders, node = NULL, alpha = .1) ### => Compare both models with the difference in deviances -## wicth follows a chi^2 with df = 1. +## which follows a chi^2 with df = 1. } \keyword{models} diff --git a/man/write.tree.Rd b/man/write.tree.Rd index 64e0c70..7a51244 100644 --- a/man/write.tree.Rd +++ b/man/write.tree.Rd @@ -26,7 +26,7 @@ write.tree(phy, file = "", append = FALSE, } \value{ a vector of mode character if \code{file = ""}, none (invisible - `NULL') otherwise. + \code{NULL}) otherwise. } \details{ The node labels and the root edge length, if available, are written in @@ -34,7 +34,9 @@ write.tree(phy, file = "", append = FALSE, If \code{tree.names == TRUE} then a variant of the Newick format is written for which the name of a tree precedes the Newick format tree - (parentheses are eventually deleted beforehand). + (parentheses are eventually deleted beforehand). The tree names are + taken from the \code{names} attribute if present (they are ignored if + \code{tree.names} is a character vector). } \references{ Felsenstein, J. The Newick tree format. @@ -44,7 +46,9 @@ write.tree(phy, file = "", append = FALSE, \url{http://evolution.genetics.washington.edu/phylip/newick_doc.html} } -\author{Emmanuel Paradis and Daniel Lawson \email{dan.lawson@bristol.ac.uk}} +\author{Emmanuel Paradis, Daniel Lawson + \email{dan.lawson@bristol.ac.uk}, and Klaus Schliep + \email{kschliep@snv.jussieu.fr}} \seealso{ \code{\link{read.tree}}, \code{\link{read.nexus}}, \code{\link{write.nexus}} -- 2.39.2