+ 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
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 <Emmanuel.Paradis@ird.fr>
-## birthdeath.R (2010-10-17)
+## birthdeath.R (2010-11-04)
## Estimation of Speciation and Extinction Rates
## with Birth-Death Models
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
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")
-## compar.ou.R (2010-03-15)
+## compar.ou.R (2010-11-04)
## Ornstein--Uhlenbeck Model for Continuous Characters
## 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))
-## summary.phylo.R (2010-05-25)
+## summary.phylo.R (2010-11-03)
## Print Summary of a Phylogeny and "multiPhylo" operators
if (details)
for (i in 1:N)
cat("tree", i, ":", length(x[[i]]$tip.label), "tips\n")
- cat("\n")
}
"[[.multiPhylo" <- function(x, i)
\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
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{
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}
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))
### ... 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}
}
\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
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.
\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}}
projects / ape.git / commitdiff