From 91cbce9b55b05cef1f7167f646bc30b3e568ebf9 Mon Sep 17 00:00:00 2001 From: paradis Date: Tue, 9 Jun 2009 06:51:23 +0000 Subject: [PATCH] fixed null model in yule.cov + various things git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@75 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 9 ++++++++- DESCRIPTION | 2 +- R/CDAM.global.R | 2 +- R/DNA.R | 1 + R/yule.R | 23 +++++++++++++++++------ man/makeNodeLabel.Rd | 2 +- man/yule.cov.Rd | 20 ++++++++++---------- 7 files changed, 39 insertions(+), 20 deletions(-) diff --git a/ChangeLog b/ChangeLog index accd346..8ad799a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -5,6 +5,9 @@ NEW FEATURES 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 @@ -15,6 +18,10 @@ 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 @@ -23,7 +30,7 @@ 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. diff --git a/DESCRIPTION b/DESCRIPTION index 7f69fe2..fac51bd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ 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 diff --git a/R/CDAM.global.R b/R/CDAM.global.R index 184bffa..5644ca5 100644 --- a/R/CDAM.global.R +++ b/R/CDAM.global.R @@ -143,7 +143,7 @@ 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') diff --git a/R/DNA.R b/R/DNA.R index b73d795..d74d188 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -19,6 +19,7 @@ del.gaps <- function(x) 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) } diff --git a/R/yule.R b/R/yule.R index bb3d226..041bda5 100644 --- a/R/yule.R +++ b/R/yule.R @@ -1,11 +1,11 @@ -## 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. @@ -40,14 +40,14 @@ yule.cov <- function(phy, formula, data = NULL) 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) @@ -55,10 +55,15 @@ yule.cov <- function(phy, formula, data = NULL) 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") @@ -68,4 +73,10 @@ yule.cov <- function(phy, formula, data = NULL) 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") } diff --git a/man/makeNodeLabel.Rd b/man/makeNodeLabel.Rd index 501f405..6fda0e6 100644 --- a/man/makeNodeLabel.Rd +++ b/man/makeNodeLabel.Rd @@ -23,7 +23,7 @@ makeNodeLabel(phy, method = "number", prefix = "Node", nodeList = list(), ...) 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 diff --git a/man/yule.cov.Rd b/man/yule.cov.Rd index b749aa5..7e3b1a2 100644 --- a/man/yule.cov.Rd +++ b/man/yule.cov.Rd @@ -61,14 +61,18 @@ yule.cov(phy, formula, data = NULL) 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 @@ -87,13 +91,9 @@ x <- rnorm(45) # the tree has 23 tips and 22 nodes ### 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} -- 2.39.2