o There is now a c() method for lists of class "DNAbin".
+ o yule.cov() now fits the null model, and its help page has been
+ corrected with respect to this change.
+
BUG FIXES
o boot.phylo() now treats correctly data frames.
+ o del.gaps() did not copy the rownames of a matrix.
+
+ o A small bug was fixed in CDAM.global().
+
OTHER CHANGES
o Instances of the form class(phy) == "phylo" have been replaced
by inherits(phy, "phylo").
- o rcoal() should be faster.
+ o rcoal() is now faster.
Package: ape
Version: 2.3-1
-Date: 2009-05-19
+Date: 2009-06-08
Title: Analyses of Phylogenetics and Evolution
Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne
Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
Rmat <- apply(Y,2,rank)
## Correction factors for tied ranks (eq. 3.3)
- t.ranks <- apply(Rmat, 2, function(x) summary(as.factor(x)))
+ t.ranks <- apply(Rmat, 2, function(x) summary(as.factor(x), maxsum=nd))
TT <- sum(unlist(lapply(t.ranks, function(x) sum((x^3)-x))))
# if(!silent) cat("TT = ",TT,'\n')
n <- dim(x)[1]
y <- vector("list", n)
for (i in 1:n) y[[i]] <- x[i, ]
+ names(y) <- rownames(x)
x <- y
rm(y)
}
-## yule.R (2007-10-18)
+## yule.R (2009-06-08)
## Fits Yule Model to a Phylogenetic Tree
## yule: standard Yule model (constant birth rate)
## yule.cov: Yule model with covariates
-## Copyright 2003-2007 Emmanuel Paradis
+## Copyright 2003-2009 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
bt <- rev(bt) # branching times from past to present
ni <- cumsum(rev(table(bt))) + 1
X <- model.matrix(formula, data)
- Xi <- X[phy$edge[, 1], ]
- Xj <- X[phy$edge[, 2], ]
+ Xi <- X[phy$edge[, 1], , drop = FALSE]
+ Xj <- X[phy$edge[, 2], , drop = FALSE]
dev <- function(b) {
2 * sum(((1/(1 + exp(-(Xi %*% b)))) +
(1/(1 + exp(-(Xj %*% b)))))
* phy$edge.length/2) -
2 * (sum(log(ni[-length(ni)])) +
- sum(log((1/(1 + exp(-(X[-(1:(n + 1)), ] %*% b)))))))
+ sum(log((1/(1 + exp(-(X[-(1:(n + 1)), , drop = FALSE] %*% b)))))))
}
out <- nlm(function(p) dev(p), p = c(rep(0, ncol(X) - 1), -1),
hessian = TRUE)
para <- matrix(NA, ncol(X), 2)
para[, 1] <- out$estimate
if (any(out$gradient == 0))
- warning("The likelihood gradient seems flat in at least one dimension (null gradient):\ncannot compute the standard-errors of the transition rates.\n")
+ warning("The likelihood gradient seems flat in at least one dimension (null gradient):\ncannot compute the standard-errors of the parameters.\n")
else para[, 2] <- sqrt(diag(solve(out$hessian)))
rownames(para) <- colnames(X)
colnames(para) <- c("Estimate", "StdErr")
+ ## fit the intercept-only model:
+ X <- model.matrix(~ 1, data = data.frame(X))
+ Xi <- X[phy$edge[, 1], , drop = FALSE]
+ Xj <- X[phy$edge[, 2], , drop = FALSE]
+ Dev.null <- nlm(function(p) dev(p), p = -1)$minimum
cat("\n---- Yule Model with Covariates ----\n\n")
cat(" Phylogenetic tree:", deparse(substitute(phy)), "\n")
cat(" Number of tips:", n, "\n")
cat(" Parameter estimates:\n")
print(para)
cat("\n")
+ cat("Null Deviance:", Dev.null, "\n")
+ cat(" Test of the fitted model: ")
+ chi <- Dev.null - Dev
+ df <- nrow(para) - 1
+ cat("chi^2 =", round(chi, 3), " df =", df,
+ " P =", round(1 - pchisq(chi, df), 3), "\n")
}
The three methods are described below:
\itemize{
- \item{``number''}{The labels are created with 1, 2, \dots suffixed
+ \item{``number''}{The labels are created with 1, 2, \dots prefixed
with the argument \code{prefix}; thus the default is to have
Node1, Node2, \dots Set \code{prefix = ""} to have only numbers.}
\item{``md5sum''}{For each node, the labels of the tips descendant
The user must obtain the values for the nodes separately.
- Note that the method in its present implementation assumes that the
- change in a species trait is more or less continuous between two nodes
- or between a node and a tip. Thus reconstructing the ancestral values
- with a Brownian motion model may be consistent with the present
- method. This can be done with the function \code{\link{ace}}.
+Note that the method in its present implementation assumes that the
+change in a species trait is more or less continuous between two nodes
+or between a node and a tip. Thus reconstructing the ancestral values
+with a Brownian motion model may be consistent with the present
+method. This can be done with the function \code{\link{ace}}.
}
\value{
- A NULL value is returned, the results are simply printed.
+ A NULL value is returned, the results are simply printed. The output
+ includes the deviance of the null (intercept-only) model and a
+ likelihood-ratio test of the fitted model against the null model.
+ Note that the deviance of the null model is different from the one
+ returned by \code{\link{yule}} because of the different parametrizations.
}
\references{
Paradis, E. (2005) Statistical analysis of diversification with
### the standard-error for x should be as large as
### the estimated parameter
yule.cov(bird.orders, ~ x)
-### compare with the simple Yule model, eventually
-### with a likelihood ratio test
-yule(bird.orders)
### another example with a tree that has a multichotomy
data(bird.families)
y <- rnorm(272) # 137 tips + 135 nodes
yule.cov(bird.families, ~ y)
-yule(multi2di(bird.families))
}
\keyword{models}
projects / ape.git / commitdiff