From: paradis Date: Mon, 23 Mar 2009 07:49:41 +0000 (+0000) Subject: change nlm to nlminb in ace() + new makeNodeLabel + fixed drop.tip X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=f5793d2af1299540f8ff7ad9ec8c1a6de1fee4cf change nlm to nlminb in ace() + new makeNodeLabel + fixed drop.tip git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@66 6e262413-ae40-0410-9e79-b911bd7a66b7 --- diff --git a/ChangeLog b/ChangeLog index b6bffe8..d659113 100644 --- a/ChangeLog +++ b/ChangeLog @@ -10,13 +10,22 @@ NEW FEATURES o The new function yule.time fits a user-defined time-dependent Yule model by maximum likelihood. + o The new function makeNodeLabel creates and/or modifies node + labels in a flexible way. + BUG FIXES o drop.tip() shuffled tip labels in some cases. + o drop.tip() did not handle node.label correctly. + o is.ultrametric() now checks the ordering of the edge matrix. + o ace() sometimes returned negative values of likelihoods of + ancestral states (thanks to Dan Rabosky for solving this long + lasting bug). + OTHER CHANGES diff --git a/DESCRIPTION b/DESCRIPTION index 1cc4f4a..51c4cfc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 2.3 -Date: 2009-03-10 +Date: 2009-03-22 Title: Analyses of Phylogenetics and Evolution Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, diff --git a/R/ace.R b/R/ace.R index c38faeb..0439600 100644 --- a/R/ace.R +++ b/R/ace.R @@ -1,4 +1,4 @@ -## ace.R (2009-01-19) +## ace.R (2009-03-22) ## Ancestral Character Estimation @@ -170,13 +170,18 @@ as the number of categories in `x'") if (output.liks) return(liks[-(1:nb.tip), ]) - 2 * log(sum(liks[nb.tip + 1, ])) } - out <- nlm(function(p) dev(p), p = rep(ip, length.out = np), - hessian = TRUE) - obj$loglik <- -out$minimum / 2 - obj$rates <- out$estimate - if (any(out$gradient == 0)) + out <- nlminb(rep(ip, length.out = np), function(p) dev(p), + lower = rep(0, np), upper = rep(Inf, np)) + obj$loglik <- -out$objective/2 + obj$rates <- out$par + oldwarn <- options("warn") + options(warn = -1) + h <- nlm(function(p) dev(p), p = obj$rates, iterlim = 1, + stepmax = 0, hessian = TRUE)$hessian + options(oldwarn) + if (any(h == 0)) warning("The likelihood gradient seems flat in at least one dimension (gradient null):\ncannot compute the standard-errors of the transition rates.\n") - else obj$se <- sqrt(diag(solve(out$hessian))) + else obj$se <- sqrt(diag(solve(h))) obj$index.matrix <- index.matrix if (CI) { lik.anc <- dev(obj$rates, TRUE) diff --git a/R/drop.tip.R b/R/drop.tip.R index 69c9576..c0bda8b 100644 --- a/R/drop.tip.R +++ b/R/drop.tip.R @@ -1,4 +1,4 @@ -## drop.tip.R (2009-03-04) +## drop.tip.R (2009-03-22) ## Remove Tips in a Phylogenetic Tree @@ -76,6 +76,7 @@ drop.tip <- { if (class(phy) != "phylo") stop('object "phy" is not of class "phylo"') + phy <- reorder(phy) Ntip <- length(phy$tip.label) NEWROOT <- ROOT <- Ntip + 1 Nnode <- phy$Nnode @@ -145,7 +146,10 @@ drop.tip <- TIPS <- phy$edge[, 2] <= Ntip ## keep the ordering so no need to reorder tip.label: phy$edge[TIPS, 2] <- rank(phy$edge[TIPS, 2]) - Ntip <- length(phy$tip.label) # update Ntip + ## 3) update node.label if needed + if (!is.null(phy$node.label)) + phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip] + Ntip <- length(phy$tip.label) # 4) update Ntip ## make new tip labels if necessary if (subtree || !trim.internal) { diff --git a/R/makeNodeLabel.R b/R/makeNodeLabel.R new file mode 100644 index 0000000..d679f5e --- /dev/null +++ b/R/makeNodeLabel.R @@ -0,0 +1,61 @@ +## makeNodeLabel.R (2009-03-22) + +## Makes Node Labels + +## Copyright 2009 Emmanuel Paradis + +## This file is part of the R-package `ape'. +## See the file ../COPYING for licensing issues. + +makeNodeLabel <- function(phy, method = "number", prefix = "Node", + nodeList = list(), ...) +{ + method <- sapply(method, match.arg, c("number", "md5sum", "user"), + USE.NAMES = FALSE) + + if ("number" %in% method) + phy$node.label <- paste(prefix, 1:phy$Nnode, sep = "") + + if ("md5sum" %in% method) { + nl <- character(phy$Nnode) + pp <- prop.part(phy, check.labels = FALSE) + labs <- attr(pp, "labels") + fl <- tempfile() + for (i in seq_len(phy$Nnode)) { + cat(sort(labs[pp[[i]]]), sep = "\n", file = fl) + nl[i] <- tools::md5sum(fl) + } + unlink(fl) + phy$node.label <- nl + } + + if ("user" %in% method) { + if (is.null(phy$node.label)) + phy$node.label <- character(phy$Nnode) + nl <- names(nodeList) + if (is.null(nl)) stop("argument 'nodeList' has no names") + Ntip <- length(phy$tip.label) + seq.nod <- .Call("seq_root2tip", phy$edge, Ntip, phy$Nnode, + PACKAGE = "ape") + ## a local version to avoid the above call many times: + .getMRCA <- function(seq.nod, tip) { + sn <- seq.nod[tip] + MRCA <- Ntip + 1 + i <- 2 + repeat { + x <- unique(unlist(lapply(sn, "[", i))) + if (length(x) != 1) break + MRCA <- x + i <- i + 1 + } + MRCA + } + for (i in seq_along(nodeList)) { + tips <- sapply(nodeList[[i]], grep, phy$tip.label, ..., + USE.NAMES = FALSE) + j <- .getMRCA(seq.nod, unique(unlist(tips))) + phy$node.label[j - Ntip] <- nl[i] + } + } + phy +} diff --git a/Thanks b/Thanks index 8f93403..57e3cb2 100644 --- a/Thanks +++ b/Thanks @@ -6,8 +6,8 @@ comments, or bug reports: thanks to all of you! Significant bug fixes were provided by Cécile Ané, James Bullard, Éric Durand, Olivier François, Bret Larget, Nick Matzke, -Elizabeth Purdom, Klaus Schliep, Li-San Wang, Yan Wong, and Peter -Wragg. Contact me if I forgot someone. +Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Li-San Wang, Yan Wong, +and Peter Wragg. Contact me if I forgot someone. Kurt Hornik, of the R Core Team, helped in several occasions to fix some problems and bugs. diff --git a/man/ace.Rd b/man/ace.Rd index 85cfb21..73934df 100644 --- a/man/ace.Rd +++ b/man/ace.Rd @@ -41,7 +41,7 @@ ace(x, phy, type = "continuous", method = "ML", CI = TRUE, \item{k}{a numeric value giving the penalty per estimated parameter; the default is \code{k = 2} which is the classical Akaike information criterion.} - \item{...}{further arguments passed to or from other methods.} + \item{\dots}{further arguments passed to or from other methods.} } \description{ This function estimates ancestral character states, and the associated diff --git a/man/makeNodeLabel.Rd b/man/makeNodeLabel.Rd new file mode 100644 index 0000000..501f405 --- /dev/null +++ b/man/makeNodeLabel.Rd @@ -0,0 +1,69 @@ +\name{makeNodeLabel} +\alias{makeNodeLabel} +\title{Makes Node Labels} +\description{ + This function makes node labels in a tree in a flexible way. +} + +\usage{ +makeNodeLabel(phy, method = "number", prefix = "Node", nodeList = list(), ...) +} +\arguments{ + \item{phy}{an object of class \code{"phylo"}.} + \item{method}{a character string giving the method used to create the + labels. Three choices are possible: \code{"number"} (the default), + \code{"md5sum"}, and \code{"user"}, or any unambiguous abbreviation + of these.} + \item{prefix}{the prefix used if \code{method = "number"}.} + \item{nodeList}{a named list specifying how nodes are names if + \code{method = "user"} (see details and examples).} + \item{\dots}{further arguments passed to \code{grep}.} +} +\details{ + The three methods are described below: + + \itemize{ + \item{``number''}{The labels are created with 1, 2, \dots suffixed + 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 + from this node are extracted, sorted alphabetically, and written + into a temporary file, then the md5sum of this file is extracted + and used as label. This results in a 32-character string which is + unique (even accross trees) for a given set of tip labels.} + \item{``user''}{the argument \code{nodeList} must be a list with + names, the latter will be used as node labels. For each element of + \code{nodeList}, the tip labels of the tree are searched for + patterns present in this element: this is done using + \code{\link[base]{grep}}. Then the most recent common ancestor of + the matching tips is given the corresponding names as labels. This + is repeated for each element of \code{nodeList}.} + } + + The method \code{"user"} can be used in combination with either of the + two others (see examples). Note that this method only modifies the + specified node labels (so that if the other nodes have already labels + they are not modified) while the two others change all labels. +} + +\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}} +\seealso{ + \code{\link{makeLabel}}, \code{\link[base]{grep}} +} +\examples{ +tr <- +"((Pan_paniscus,Pan_troglodytes),((Homo_sapiens,Hom_erectus),Homo_abilis));" +tr <- read.tree(text = tr) +tr <- makeNodeLabel(tr, "u", nodeList = list(Pan = "Pan", Homo = "Homo")) +plot(tr, show.node.label = TRUE) +### does not erase the previous node labels: +tr <- makeNodeLabel(tr, "u", nodeList = list(Hominid = c("Pan","Homo"))) +plot(tr, show.node.label = TRUE) +### the two previous commands could be combined: +L <- list(Pan = "Pan", Homo = "Homo", Hominid = c("Pan","Homo")) +tr <- makeNodeLabel(tr, "u", nodeList = L) +### combining different methods: +tr <- makeNodeLabel(tr, c("n", "u"), prefix = "#", nodeList = list(Hominid = c("Pan","Homo"))) +plot(tr, show.node.label = TRUE) +} +\keyword{manip}