Package: ape
Version: 3.0
-Date: 2012-02-03
+Date: 2012-02-10
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, Andrei-Alin Popescu, Klaus Schliep, Korbinian Strimmer, Damien de Vienne
Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
o boot.phylo() and prop.clades() have a new option rooted =
FALSE. Note that the behaviour of prop.part() is unchanged.
+ o edgelabels() has a new option 'date' to place labels on edges of
+ chronograms using the time scale (suggestion by Rob Lanfear).
+
BUG FIXES
- o mantel.test() could return P-values > 1 with the default
+ o In chronopl(), the code setting the initial dates failed in
+ complicated settings (several dates known within intervals).
+ This has been generally improved and should result in faster
+ and more efficient convergence even in simple settings.
+
+ o mantel.test() sometimes returned P-values > 1 with the default
two-tailed test.
o extract.clade() sometimes shuffled some tip labels (thanks to
o The C code of base.freq() has been improved and is now nearly 8
times faster.
+ o The option 'original.data' of write.nexus() is now deprecated and
+ will be removed in a future release.
+
+ o The code of is.ultrametric() has been improved and is now 3 times
+ faster.
+
+ o The code of vcv.phylo() has been improved and is now 10 or 30
+ times faster for 100 or 1000 tips, respectively. Consequently,
+ fitting models with PGLS will be faster overall.
+
CHANGES IN APE VERSION 2.8
OTHER CHANGES
- o identify.phylo() now returns NULL if the user right-(instead of
- left-)clicks (an error was returned previously).
+ o identify.phylo() now returns NULL if the user right- (instead of
+ left-) clicks (an error was returned previously).
o read.nexus() should be robust to commands inserted in the TREES
block.
-## PGLS.R (2008-05-08)
+## PGLS.R (2012-02-09)
## Phylogenetic Generalized Least Squares
-## Copyright 2004 Julien Dutheil, and 2006-2008 Emmanuel Paradis
+## Copyright 2004 Julien Dutheil, and 2006-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-corBrownian <- function(value = 1, phy, form=~1)
+corBrownian <- function(value = 1, phy, form = ~1)
{
- if (!("phylo" %in% class(phy)))
- stop("object \"phy\" is not of class \"phylo\"")
- attr(value, "formula") <- form
- attr(value, "fixed") <- TRUE
- attr(value, "tree") <- phy
- class(value) <- c("corBrownian", "corPhyl", "corStruct")
- return(value)
+ if (!inherits(phy, "phylo"))
+ stop('object "phy" is not of class "phylo"')
+ attr(value, "formula") <- form
+ attr(value, "fixed") <- TRUE
+ attr(value, "tree") <- phy
+ class(value) <- c("corBrownian", "corPhyl", "corStruct")
+ value
}
-corMartins <- function(value, phy, form=~1, fixed=FALSE)
+corMartins <- function(value, phy, form = ~1, fixed = FALSE)
{
- if(length(value) > 1)
- stop("only one parameter is allowed in corPGLS structure.")
- if(value < 0) stop("parameter alpha must be positive.")
- if (!("phylo" %in% class(phy)))
- stop("object \"phy\" is not of class \"phylo\"")
- attr(value, "formula") <- form
- attr(value, "fixed") <- fixed
- attr(value, "tree") <- phy
- class(value) <- c("corMartins", "corPhyl", "corStruct")
- return(value)
+ if (length(value) > 1)
+ stop("only one parameter is allowed")
+ if (value < 0) stop("the parameter alpha must be positive")
+ if (!inherits(phy, "phylo"))
+ stop('object "phy" is not of class "phylo"')
+ attr(value, "formula") <- form
+ attr(value, "fixed") <- fixed
+ attr(value, "tree") <- phy
+ class(value) <- c("corMartins", "corPhyl", "corStruct")
+ value
}
-corGrafen <- function(value, phy, form=~1, fixed=FALSE)
+corGrafen <- function(value, phy, form = ~1, fixed = FALSE)
{
- if(length(value) > 1)
- stop("only one parameter is allowed in corGrafen structure.")
- if(value < 0) stop("parameter rho must be positive.")
- value <- log(value) # Optimization under constraint, use exponential transform.
- if (!("phylo" %in% class(phy)))
- stop("object \"phy\" is not of class \"phylo\"")
- attr(value, "formula") <- form
- attr(value, "fixed") <- fixed
- attr(value, "tree") <- phy
- class(value) <- c("corGrafen", "corPhyl", "corStruct")
- return(value)
+ if (length(value) > 1)
+ stop("only one parameter is allowed")
+ if (value < 0) stop("parameter rho must be positive")
+ value <- log(value) # Optimization under constraint, use exponential transform.
+ if (!inherits(phy, "phylo"))
+ stop('object "phy" is not of class "phylo"')
+ attr(value, "formula") <- form
+ attr(value, "fixed") <- fixed
+ attr(value, "tree") <- phy
+ class(value) <- c("corGrafen", "corPhyl", "corStruct")
+ value
}
Initialize.corPhyl <- function(object, data, ...)
{
- # The same as in Initialize corStruct:
- form <- formula(object)
- ## Obtaining the group information, if any
- if(!is.null(getGroupsFormula(form))) {
- attr(object, "groups") <- getGroups(object, form, data = data)
- attr(object, "Dim") <- Dim(object, attr(object, "groups"))
- } else { # no groups
- attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
- }
- ## Obtaining the covariate(s)
- attr(object, "covariate") <- getCovariate(object, data = data)
+ ## The same as in Initialize corStruct:
+ form <- formula(object)
+ ## Obtaining the group information, if any
+ if (!is.null(getGroupsFormula(form))) {
+ attr(object, "groups") <- getGroups(object, form, data = data)
+ attr(object, "Dim") <- Dim(object, attr(object, "groups"))
+ } else { # no groups
+ attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
+ }
+ ## Obtaining the covariate(s)
+ attr(object, "covariate") <- getCovariate(object, data = data)
- # Specific to corPhyl:
- phy <- attr(object, "tree")
- if (is.null(data))
- data <- parent.frame()
- ## Added by EP 29 May 2006:
- if (nrow(data) != length(phy$tip.label))
- stop("number of observations and number of tips in the tree are not equal.")
- ## END
- if(is.null(rownames(data))) {
- warning("No row names supplied in dataframe, data taken to be in the same order as in tree.")
- attr(object, "index") <- 1:dim(data)[1]
- } else {
- index <- match(rownames(data), phy$tip.label)
- if(any(is.na(index))) {
- warning("Row names in dataframe do not match tree tip names. data taken to be in the same order as in tree.")
- attr(object, "index") <- 1:dim(data)[1]
+ ## Specific to corPhyl:
+ phy <- attr(object, "tree")
+ if (is.null(data)) data <- parent.frame()
+ ## Added by EP 29 May 2006:
+ if (nrow(data) != length(phy$tip.label))
+ stop("number of observations and number of tips in the tree are not equal.")
+ ## END
+ if (is.null(rownames(data))) {
+ warning("No rownames supplied in data frame, data taken to be in the same order than in tree")
+ attr(object, "index") <- 1:dim(data)[1]
} else {
- attr(object, "index") <- index
+ index <- match(rownames(data), phy$tip.label)
+ if (any(is.na(index))) {
+ warning("Rownames in data frame do not match tree tip names; data taken to be in the same order as in tree")
+ attr(object, "index") <- 1:dim(data)[1]
+ } else {
+ attr(object, "index") <- index
+ }
}
- }
- return(object)
+ object
}
-corMatrix.corBrownian <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
+corMatrix.corBrownian <-
+ function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
- if (!("corBrownian" %in% class(object)))
- stop("object is not of class \"corBrownian\".")
- if(!any(attr(object, "index")))
- stop("object have not been initialized.")
- tree <- attr(object, "tree")
- mat <- vcv.phylo(tree, cor = corr)
- n <- dim(mat)[1]
- # reorder matrix:
- matr <- matrix(nrow=n, ncol=n)
- index <- attr(object, "index")
- for(i in 1:n)
- for(j in i:n)
- matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
- return(matr)
+ if (!("corBrownian" %in% class(object)))
+ stop('object is not of class "corBrownian"')
+ if (!any(attr(object, "index")))
+ stop("object has not been initialized.")
+ tree <- attr(object, "tree")
+ mat <- vcv.phylo(tree, corr = corr)
+ n <- dim(mat)[1]
+ ## reorder matrix:
+ index <- attr(object, "index")
+ mat[index, index]
}
corMatrix.corMartins <-
function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
if (!("corMartins" %in% class(object)))
- stop("object is not of class \"corMartins\".")
- if(!any(attr(object, "index")))
- stop("object have not been initialized.")
+ stop('object is not of class "corMartins"')
+ if (!any(attr(object, "index")))
+ stop("object has not been initialized.")
tree <- attr(object, "tree")
dist <- cophenetic.phylo(tree)
mat <- exp(-object[1] * dist)
if (corr) mat <- cov2cor(mat)
n <- dim(mat)[1]
## reorder matrix:
- matr <- matrix(nrow=n, ncol=n)
index <- attr(object, "index")
- for(i in 1:n)
- for(j in i:n)
- matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
- matr
+ mat[index, index]
}
-corMatrix.corGrafen <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
+corMatrix.corGrafen <-
+ function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
- if (!("corGrafen" %in% class(object)))
- stop("object is not of class \"corGrafen\".")
- if (!any(attr(object, "index")))
- stop("object have not been initialized.")
- tree <- compute.brlen(attr(object, "tree"),
- method = "Grafen", power = exp(object[1]))
- mat <- vcv.phylo(tree, cor = corr)
- n <- dim(mat)[1]
- # reorder matrix:
- matr <- matrix(nrow=n, ncol=n)
- index <- attr(object, "index")
- for(i in 1:n)
- for(j in i:n)
- matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
- return(matr)
+ if (!("corGrafen" %in% class(object)))
+ stop('object is not of class "corGrafen"')
+ if (!any(attr(object, "index")))
+ stop("object has not been initialized.")
+ tree <- compute.brlen(attr(object, "tree"),
+ method = "Grafen", power = exp(object[1]))
+ mat <- vcv.phylo(tree, corr = corr)
+ n <- dim(mat)[1]
+ ## reorder matrix:
+ index <- attr(object, "index")
+ mat[index, index]
}
coef.corBrownian <- function(object, unconstrained = TRUE, ...)
{
- if (!("corBrownian" %in% class(object)))
- stop("object is not of class \"corBrownian\".")
- return(numeric(0))
+ if (!("corBrownian" %in% class(object)))
+ stop('object is not of class "corBrownian"')
+ numeric(0)
}
coef.corMartins <- function(object, unconstrained = TRUE, ...)
{
- if (!("corMartins" %in% class(object)))
- stop("object is not of class \"corMartins\".")
- if(unconstrained) {
- if(attr(object, "fixed")) {
- return(numeric(0))
- } else {
- return(as.vector(object))
+ if (!("corMartins" %in% class(object)))
+ stop('object is not of class "corMartins"')
+ if (unconstrained) {
+ if (attr(object, "fixed")) {
+ return(numeric(0))
+ } else {
+ return(as.vector(object))
+ }
}
- }
- aux <- as.vector(object)
- names(aux) <- "alpha"
- return(aux)
+ aux <- as.vector(object)
+ names(aux) <- "alpha"
+ aux
}
coef.corGrafen <- function(object, unconstrained = TRUE, ...)
{
- if (!("corGrafen" %in% class(object)))
- stop("object is not of class \"corGrafen\".")
- if(unconstrained) {
- if(attr(object, "fixed")) {
- return(numeric(0))
- } else {
- return(as.vector(object))
+ if (!("corGrafen" %in% class(object)))
+ stop('object is not of class "corGrafen"')
+ if (unconstrained) {
+ if (attr(object, "fixed")) {
+ return(numeric(0))
+ } else {
+ return(as.vector(object))
+ }
}
- }
- aux <- exp(as.vector(object))
- names(aux) <- "rho"
- return(aux)
+ aux <- exp(as.vector(object))
+ names(aux) <- "rho"
+ aux
}
### removed node.sons() and node.leafnumber() (2006-10-12)
{
if (value < 0 || value > 1)
stop("the value of lambda must be between 0 and 1.")
- if (!("phylo" %in% class(phy)))
+ if (!inherits(phy, "phylo"))
stop('object "phy" is not of class "phylo"')
attr(value, "formula") <- form
attr(value, "fixed") <- fixed
{
if (!any(attr(object, "index")))
stop("object has not been initialized")
- mat <- vcv.phylo(attr(object, "tree"), cor = corr)
+ mat <- vcv.phylo(attr(object, "tree"), corr = corr)
index <- attr(object, "index")
mat <- mat[index, index]
tmp <- diag(mat)
{
if (value <= 0)
stop("the value of g must be greater than 0.")
- if (!("phylo" %in% class(phy)))
+ if (!inherits(phy, "phylo"))
stop('object "phy" is not of class "phylo"')
attr(value, "formula") <- form
attr(value, "fixed") <- fixed
phy <- attr(object, "tree")
d <- (dist.nodes(phy)[length(phy$tip.label) + 1, ])^(1/object[1])
phy$edge.length <- d[phy$edge[, 2]] - d[phy$edge[, 1]]
- mat <- vcv.phylo(phy, cor = corr)
+ mat <- vcv.phylo(phy, corr = corr)
mat[index, index]
}
scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1)
{
if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo".')
+ stop('object "phy" is not of class "phylo"')
if (is.null(phy$edge.length))
stop("tree has no branch lengths")
type <- match.arg(type, c("continuous", "discrete"))
-## chronopl.R (2011-07-04)
+## chronopl.R (2012-02-09)
## Molecular Dating With Penalized Likelihood
-## Copyright 2005-2011 Emmanuel Paradis
+## Copyright 2005-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
CV = FALSE, eval.max = 500, iter.max = 500, ...)
{
n <- length(phy$tip.label)
- ROOT <- n + 1
- if (length(node) == 1 && node == "root") node <- ROOT
+ ROOT <- n + 1L
+ if (identical(node, "root")) node <- ROOT
if (any(node <= n))
stop("node numbers should be greater than the number of tips")
zerobl <- which(phy$edge.length <= 0)
}
m <- phy$Nnode
el <- phy$edge.length
- e <- phy$edge
- N <- dim(e)[1]
+ e1 <- phy$edge[, 1L]
+ e2 <- phy$edge[, 2L]
+ N <- length(e1)
TIPS <- 1:n
EDGES <- 1:N
## `basal' contains the indices of the basal edges
## (ie, linked to the root):
- basal <- which(e[, 1] == ROOT)
+ basal <- which(e1 == ROOT)
Nbasal <- length(basal)
## `ind' contains in its 1st column the index of all nonbasal
ind <- matrix(0L, N - Nbasal, 2)
ind[, 1] <- EDGES[-basal]
- ind[, 2] <- match(e[EDGES[-basal], 1], e[, 2])
+ ind[, 2] <- match(e1[EDGES[-basal]], e2)
age <- numeric(n + m)
- ##ini.time <- node.depth(phy)[-TIPS] - 1
- ini.time <- node.depth(phy) - 1
+#############################################################################
+### This bit sets 'ini.time' and should result in no negative branch lengths
+
+ seq.nod <- .Call("seq_root2tip", phy$edge, n, phy$Nnode, PACKAGE = "ape")
+
+ ini.time <- age
+ ini.time[ROOT:(n + m)] <- NA
+ ini.time[node] <- if (is.null(age.max)) age.min else (age.min + age.max) / 2
+
+ ## if no age given for the root, find one approximately:
+ if (is.na(ini.time[ROOT]))
+ ini.time[ROOT] <- if (is.null(age.max)) 3 * max(age.min) else 3 * max(age.max)
+
+ ISnotNA.ALL <- unlist(lapply(seq.nod, function(x) sum(!is.na(ini.time[x]))))
+ o <- order(ISnotNA.ALL, decreasing = TRUE)
+
+ for (y in seq.nod[o]) {
+ ISNA <- is.na(ini.time[y])
+ if (any(ISNA)) {
+ i <- 2L # we know the 1st value is not NA, so we start at the 2nd one
+ while (i <= length(y)) {
+ if (ISNA[i]) { # we stop at the next NA
+ j <- i + 1L
+ while (ISNA[j]) j <- j + 1L # look for the next non-NA
+ nb.val <- j - i
+ by <- (ini.time[y[i - 1L]] - ini.time[y[j]]) / (nb.val + 1)
+ ini.time[y[i:(j - 1L)]] <- ini.time[y[i - 1L]] - by * seq_len(nb.val)
+ i <- j + 1L
+ } else i <- i + 1L
+ }
+ }
+ }
- ## first, rescale all times with respect to
- ## the age of the 1st element of `node':
- ratio <- age.min[1]/ini.time[node[1]]
- ini.time <- ini.time*ratio
+ real.edge.length <- ini.time[e1] - ini.time[e2]
- ## because if (!is.null(age.max)), 'node' is modified, so we copy it in case CV = TRUE:
+ if (any(real.edge.length <= 0))
+ stop("some initial branch lengths are zero or negative;
+ maybe you need to adjust the given dates -- see '?chronopl' for details")
+#############################################################################
+
+ ## because if (!is.null(age.max)), 'node' is modified,
+ ## so we copy it in case CV = TRUE:
node.bak <- node
- if (length(node) > 1) {
- ini.time[node] <- age.min
- real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
- while (any(real.edge.length <= 0)) {
- for (i in EDGES) {
- if (real.edge.length[i] > 0) next
- if (e[i, 1] %in% node) {
- ini.time[e[i, 2]] <- ini.time[e[1, 2]] - 2 * real.edge.length[i]
- next
- }
- if (e[i, 2] %in% node) {
- ini.time[e[i, 1]] <- ini.time[e[1, 1]] + 2 * real.edge.length[i]
- next
- }
- ##browser()
- ini.time[e[i, 2]] <- ini.time[e[1, 2]] - real.edge.length[i]
- ini.time[e[i, 1]] <- ini.time[e[1, 1]] + real.edge.length[i]
- }
- real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
- ##print(min(real.edge.length))
- }
- }
## `unknown.ages' will contain the index of the nodes of unknown age:
unknown.ages <- n + 1:m
if (!is.null(age.max)) { # are some nodes known within some intervals?
lower[node - n] <- age.min
upper[node - n] <- age.max
+ ## find nodes known within an interval:
interv <- which(age.min != age.max)
+ ## drop them from the 'node' since they will be estimated:
node <- node[-interv]
- if (length(node)) age[node] <- age.min[-interv]
+ if (length(node)) age[node] <- age.min[-interv] # update 'age'
} else age[node] <- age.min
if (length(node)) {
- unknown.ages <- unknown.ages[n - node]
+ unknown.ages <- unknown.ages[n - node] # 'n - node' is simplification for '-(node - n)'
lower <- lower[n - node]
upper <- upper[n - node]
}
minusploglik.gr <- function(rate, node.time) {
grad <- numeric(N + length(unknown.ages))
age[unknown.ages] <- node.time
- real.edge.length <- age[e[, 1]] - age[e[, 2]]
+ real.edge.length <- age[e1] - age[e2]
if (any(real.edge.length < 0)) {
grad[] <- 0
return(grad)
}
for (i in EDGES) {
- ii <- c(which(e[, 2] == e[i, 1]), which(e[, 1] == e[i, 2]))
+ ii <- c(which(e2 == e1[i]), which(e1 == e2[i]))
if (!length(ii)) next
grad[i] <- grad[i] + lambda*(2*length(ii)*rate[i] - 2*sum(rate[ii]))
}
## gradient for the 'node times'
for (i in 1:length(unknown.ages)) {
nd <- unknown.ages[i]
- ii <- which(e[, 1] == nd)
+ ii <- which(e1 == nd)
grad[i + N] <-
sum(rate[ii] - el[ii]/real.edge.length[ii])#, na.rm = TRUE)
if (nd != ROOT) {
- ii <- which(e[, 2] == nd)
+ ii <- which(e2 == nd)
grad[i + N] <- grad[i + N] -
rate[ii] + el[ii]/real.edge.length[ii]
}
minusploglik <- function(rate, node.time) {
age[unknown.ages] <- node.time
- real.edge.length <- age[e[, 1]] - age[e[, 2]]
+ real.edge.length <- age[e1] - age[e2]
if (any(real.edge.length < 0)) return(1e50)
B <- rate*real.edge.length
loglik <- sum(-B + el*log(B) - lfactorial(el))
attr(phy, "message") <- out$message
age[unknown.ages] <- out$par[-EDGES]
if (CV) ophy <- phy
- phy$edge.length <- age[e[, 1]] - age[e[, 2]]
+ phy$edge.length <- age[e1] - age[e2]
if (CV) attr(phy, "D2") <-
chronopl.cv(ophy, lambda, age.min, age.max, node.bak,
n, S, tol, eval.max, iter.max, ...)
-## is.ultrametric.R (2009-05-10)
+## is.ultrametric.R (2012-02-09)
## Test if a Tree is Ultrametric
-## Copyright 2003-2009 Emmanuel Paradis
+## Copyright 2003-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
is.ultrametric <- function(phy, tol = .Machine$double.eps^0.5)
{
if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo".')
+ stop('object "phy" is not of class "phylo"')
if (is.null(phy$edge.length))
- stop("the tree has no branch lengths.")
+ stop("the tree has no branch lengths")
+
phy <- reorder(phy)
n <- length(phy$tip.label)
- n.node <- phy$Nnode
+ e1 <- phy$edge[, 1]
+ e2 <- phy$edge[, 2]
+ EL <- phy$edge.length
+
+ ## xx: vecteur donnant la distance d'un noeud
+ ## ou d'un tip à partir de la racine
+ xx <- numeric(n + phy$Nnode)
- ## xx: vecteur donnant la distance d'un
- ## noeud ou tip à partir de la racine
- xx <- numeric(n + n.node)
+ ## the following must start at the root and follow the
+ ## edges contiguously; so the tree must be either in cladewise
+ ## order (or in pruningwise but the for loop must start from
+ ## the bottom of the edge matrix)
- for (i in 1:dim(phy$edge)[1])
- xx[phy$edge[i, 2]] <- xx[phy$edge[i, 1]] + phy$edge.length[i]
+ for (i in seq_len(length(e1)))
+ xx[e2[i]] <- xx[e1[i]] + EL[i]
- if (identical(all.equal.numeric(var(xx[1:n]),
- 0, tolerance = tol), TRUE)) TRUE
- else FALSE
+ isTRUE(all.equal.numeric(var(xx[1:n]), 0, tolerance = tol))
}
-## nodelabels.R (2010-07-17)
+## nodelabels.R (2012-02-10)
## Labelling Trees
-## Copyright 2004-2010 Emmanuel Paradis, 2006 Ben Bolker, and 2006 Jim Lemon
+## Copyright 2004-2012 Emmanuel Paradis, 2006 Ben Bolker, and 2006 Jim Lemon
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
function(text, edge, adj = c(0.5, 0.5), frame = "rect",
pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
col = "black", bg = "lightgreen", horiz = FALSE,
- width = NULL, height = NULL, ...)
+ width = NULL, height = NULL, date = NULL, ...)
{
lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
if (missing(edge)) {
XX <- (lastPP$xx[subedge[, 1]] + lastPP$xx[subedge[, 2]]) / 2
YY <- (lastPP$yy[subedge[, 1]] + lastPP$yy[subedge[, 2]]) / 2
}
+
+ ## suggestion by Rob Lanfear:
+ if (!is.null(date))
+ XX[] <- max(lastPP$xx) - date
+
BOTHlabels(text, sel, XX, YY, adj, frame, pch, thermo,
pie, piecol, col, bg, horiz, width, height, ...)
}
-## rTrait.R (2011-06-15)
+## rTrait.R (2012-02-09)
## Trait Evolution
-## Copyright 2010-2011 Emmanuel Paradis
+## Copyright 2010-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
diag(Q) <- -rowSums(Q)
for (i in N:1) {
p <- matexpo(Q * el[i])[x[anc[i]], ]
- x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
+ x[des[i]] <- sample.int(k, size = 1, FALSE, prob = p)
}
}
-## read.nexus.R (2011-03-26)
+## read.nexus.R (2012-02-09)
## Read Tree File in Nexus Format
-## Copyright 2003-2011 Emmanuel Paradis and 2010 Klaus Schliep
+## Copyright 2003-2012 Emmanuel Paradis and 2010 Klaus Schliep
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
class(trees) <- "multiPhylo"
if (!all(nms.trees == "")) names(trees) <- nms.trees
}
- if (length(grep("[\\/]", file)) == 1)
- if (!file.exists(file)) # suggestion by Francois Michonneau
- file <- paste(getwd(), file, sep = "/")
- attr(trees, "origin") <- file
trees
}
-## rtree.R (2010-03-09)
+## rtree.R (2012-02-09)
## Generates Trees
-## Copyright 2004-2010 Emmanuel Paradis
+## Copyright 2004-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
{
foo <- function(n, pos) {
- n1 <- .Internal(sample(n - 1, 1, FALSE, NULL))
+ n1 <- sample.int(n - 1, 1, FALSE, NULL)
n2 <- n - n1
po2 <- pos + 2*n1 - 1
edge[c(pos, po2), 1] <<- nod
i <- which(is.na(edge[, 2]))
edge[i, 2] <- 1:n
} else { # n > 4
- n1 <- .Internal(sample(n - 2, 1, FALSE, NULL))
+ n1 <- sample.int(n - 2, 1, FALSE, NULL)
if (n1 == n - 2) {
n2 <- n3 <- 1
} else {
- n2 <- .Internal(sample(n - n1 - 1, 1, FALSE, NULL))
+ n2 <- sample.int(n - n1 - 1, 1, FALSE, NULL)
n3 <- n - n1 - n2
}
po2 <- 2*n1
tip.label <- paste("t", 1:n, sep = "")
phy <- list(edge = edge, tip.label = tip.label, Nnode = m)
class(phy) <- "phylo"
- attr(phy, "order" <- "cladewise")
+ attr(phy, "order") <- "cladewise"
phy
}
-## vcv.phylo.R (2009-11-19)
+## vcv.phylo.R (2012-02-09)
## Phylogenetic Variance-Covariance or Correlation Matrix
-## Copyright 2002-2009 Emmanuel Paradis
+## Copyright 2002-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
vcv.phylo <- function(phy, model = "Brownian", corr = FALSE, ...)
{
if (is.null(phy$edge.length))
- stop("the tree has no branch lengths")
-
- foo <- function(node, var, endofclade) {
- ## First, get the extent of clade descending
- ## from `node' in the matrix `edge':
- from <- which(phy$edge[, 1] == node)
- to <- c(from[-1] - 1, endofclade)
- ## Get the #'s of the descendants of `node':
- desc <- phy$edge[from, 2]
- ## The variance of each of these is easy:
- vcv[desc, desc] <<- var + phy$edge.length[from]
- nd <- length(desc)
- ## The variance of `node' is equal to the covariance of
- ## each possible pair among its descendant clades.
- for (i in 1:(nd - 1))
- for (j in (i + 1):nd)
- for (k in phy$edge[from[i]:to[i], 2])
- for (l in phy$edge[from[j]:to[j], 2])
- vcv[k, l] <<- vcv[l, k] <<- var
- for (i in 1:nd) {
- if (desc[i] <= n) next
- foo(desc[i], vcv[desc[i], desc[i]], to[i])
+ stop("the tree has no branch lengths")
+
+ pp <- prop.part(phy)
+
+ phy <- reorder(phy, "pruningwise")
+
+ n <- length(phy$tip.label)
+ e1 <- phy$edge[, 1]
+ e2 <- phy$edge[, 2]
+ EL <- phy$edge.length
+
+ ## xx: vecteur donnant la distance d'un noeud
+ ## ou d'un tip à partir de la racine
+ ## (same than in is.ultrametric)
+ xx <- numeric(n + phy$Nnode)
+
+ vcv <- matrix(0, n, n)
+
+ ## the loop below starts from the bottom of the edge matrix, so
+ ## from the root
+
+ for (i in length(e1):1) {
+ var.cur.node <- xx[e1[i]]
+ xx[e2[i]] <- var.cur.node + EL[i] # sets the variance
+ j <- i - 1L
+ while (e1[j] == e1[i] && j > 0) {
+ left <- if (e2[j] > n) pp[[e2[j] - n]] else e2[j]
+ right <- if (e2[i] > n) pp[[e2[i] - n]] else e2[i]
+ vcv[left, right] <- vcv[right, left] <- var.cur.node # sets the covariance
+ j <- j - 1L
}
}
- phy <- reorder(phy)
- n <- length(phy$tip.label)
- n.node <- phy$Nnode
- N <- n.node + n
- vcv <- matrix(0, N, N)
- foo(n + 1, 0, dim(phy$edge)[1])
+ diag(vcv) <- xx[1:n]
- vcv <- vcv[1:n, 1:n]
if (corr) {
## This is inspired from the code of `cov2cor' (2005-09-08):
M <- vcv
vcv[] <- Is * M * rep(Is, each = n)
vcv[1 + 0:(n - 1)*(n + 1)] <- 1
}
- rownames(vcv) <- colnames(vcv) <- phy$tip.label
+
+ dimnames(vcv)[1:2] <- list(phy$tip.label)
+
vcv
}
-## write.nexus.R (2011-03-26)
+## write.nexus.R (2012-02-09)
## Write Tree File in Nexus Format
-## Copyright 2003-2011 Emmanuel Paradis
+## Copyright 2003-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-write.nexus <- function(..., file = "", translate = TRUE, original.data = TRUE)
+write.nexus <- function(..., file = "", translate = TRUE, original.data = NULL)
{
obj <- list(...)
## We insure that all trees are in a list, even if there is a single one:
cat("#NEXUS\n", file = file)
cat(paste("[R-package APE, ", date(), "]\n\n", sep = ""),
file = file, append = TRUE)
- if (original.data) {
- if (!is.null(attr(obj[[1]], "origin"))) {
- if (!file.exists(attr(obj[[1]], "origin"))) {
- warning(paste("the file", attr(obj[[1]], "origin"),
- "cannot be found,
-the original data won't be written with the tree."))
- original.data <- FALSE
- }
- else {
- ORI <- scan(file = attr(obj[[1]], "origin"), what = character(),
- sep = "\n", skip = 1)
- start <- grep("BEGIN TAXA;", ORI)
- ORI <- ORI[-(1:(start - 1))]
- ORI <- gsub("ENDBLOCK;", "END;", ORI)
- endblock <- grep("END;", ORI)
- start <- grep("BEGIN TREES;", ORI)
- end <- endblock[endblock > start][1]
- cat(ORI[1:(start - 1)], file = file, append = TRUE, sep = "\n")
- ORI <- ORI[-(1:end)]
- }
- }
- else original.data <- FALSE
- }
+
+ if (!is.null(original.data))
+ warning("the option 'original.data' is deprecated and will be removed soon. Please update your code.")
+
N <- length(obj[[1]]$tip.label)
- if (!original.data) {
- cat("BEGIN TAXA;\n", file = file, append = TRUE)
- cat(paste("\tDIMENSIONS NTAX = ", N, ";\n", sep = ""),
- file = file, append = TRUE)
- cat("\tTAXLABELS\n", file = file, append = TRUE)
- cat(paste("\t\t", obj[[1]]$tip.label, sep = ""),
- sep = "\n", file = file, append = TRUE)
- cat("\t;\n", file = file, append = TRUE)
- cat("END;\n", file = file, append = TRUE)
- }
+
+ cat("BEGIN TAXA;\n", file = file, append = TRUE)
+ cat(paste("\tDIMENSIONS NTAX = ", N, ";\n", sep = ""),
+ file = file, append = TRUE)
+ cat("\tTAXLABELS\n", file = file, append = TRUE)
+ cat(paste("\t\t", obj[[1]]$tip.label, sep = ""),
+ sep = "\n", file = file, append = TRUE)
+ cat("\t;\n", file = file, append = TRUE)
+ cat("END;\n", file = file, append = TRUE)
+
cat("BEGIN TREES;\n", file = file, append = TRUE)
if (translate) {
cat("\tTRANSLATE\n", file = file, append = TRUE)
"\n", sep = "", file = file, append = TRUE)
}
cat("END;\n", file = file, append = TRUE)
- if(original.data) cat(ORI, file = file, append = TRUE, sep = "\n")
}
\name{chronopl}
\alias{chronopl}
\title{Molecular Dating With Penalized Likelihood}
+\description{
+ This function estimates the node ages of a tree using a
+ semi-parametric method based on penalized likelihood (Sanderson
+ 2002). The branch lengths of the input tree are interpreted as mean
+ numbers of substitutions (i.e., per site).
+}
\usage{
chronopl(phy, lambda, age.min = 1, age.max = NULL,
node = "root", S = 1, tol = 1e-8,
algorithm.}
\item{\dots}{further arguments passed to control \code{nlminb}.}
}
-\description{
- This function estimates the node ages of a tree using a
- semi-parametric method based on penalized likelihood (Sanderson
- 2002). The branch lengths of the input tree are interpreted as mean
- numbers of substitutions (i.e., per site).
-}
\details{
The idea of this method is to use a trade-off between a parametric
formulation where each branch has its own rate, and a nonparametric
c(10, 6), node = c(15, 18)} means that the age of node 15 is 10
units of time, and the age of node 18 is between 5 and 6).
+ If two nodes are linked (i.e., one is the ancestor of the other) and
+ have the same values of \code{age.min} and \code{age.max} (say, 10 and
+ 15) this will result in an error because the medians of these values
+ are used as initial times (here 12.5) giving initial branch length(s)
+ equal to zero. The easiest way to solve this is to change slightly the
+ given values, for instance use \code{age.max = 14.9} for the youngest
+ node, or \code{age.max = 15.1} for the oldest one (or similarly for
+ \code{age.min}).
+
The input tree may have multichotomies. If some internal branches are
of zero-length, they are collapsed (with a warning), and the returned
tree will have less nodes than the input one. The presence of
edgelabels(text, edge, adj = c(0.5, 0.5), frame = "rect",
pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
col = "black", bg = "lightgreen", horiz = FALSE,
- width = NULL, height = NULL, ...)
+ width = NULL, height = NULL, date = NULL, ...)
}
\arguments{
\item{horiz, width, height}{parameters controlling the aspect of
thermometers; by default, their width and height are determined
automatically.}
+ \item{date}{specifies the positions of labels on edges of chronograms
+ with respect to the time scale.}
}
\details{
These three functions have the same optional arguments and the same
\examples{
# Oribatid mite data from Borcard and Legendre (1994)
+\dontrun{
if (require(vegan)) {
data(mite) # Community composition data, 70 peat cores, 35 species
biplot(res, mite.log.st, dir.axis1=-1, dir.axis2=-1)
}
}
+}
\keyword{ multivariate }
### the two followings are the same:
rTraitDisc(bird.orders)
rTraitDisc(bird.orders, model = matrix(c(0, 0.1, 0.1, 0), 2))
+
### two-state model with irreversibility:
rTraitDisc(bird.orders, model = matrix(c(0, 0, 0.1, 0), 2))
+
+### simple two-state model:
+tr <- rcoal(n <- 40, br = runif)
+x <- rTraitDisc(tr, ancestor = TRUE)
+plot(tr, show.tip.label = FALSE)
+nodelabels(pch = 19, col = x[-(1:n)])
+tiplabels(pch = 19, col = x[1:n])
+
### an imaginary model with stasis 0.5 time unit after a node, then
### random evolution:
foo <- function(x, l) {
If several trees are read in the file, the returned object is of class
\code{"multiPhylo"}, and is a list of objects of class \code{"phylo"}.
-
- An attribute \code{"origin"} is further given to the returned object
- which gives the name of the source file (with its path). This is used
- to write a tree in a NEXUS file where all the original data must be
- written (not only the tree) in accordance to the specifications of
- Maddison et al. (1997).
}
\references{
Maddison, D. R., Swofford, D. L. and Maddison, W. P. (1997) NEXUS: an
\alias{write.nexus}
\title{Write Tree File in Nexus Format}
\usage{
-write.nexus(..., file = "", translate = TRUE, original.data = TRUE)
+write.nexus(..., file = "", translate = TRUE, original.data = NULL)
}
\arguments{
\item{\dots}{either (i) a single object of class \code{"phylo"}, (ii) a
\item{translate}{a logical, if \code{TRUE} (the default) a translation
of the tip labels is done which are replaced in the parenthetic
representation with tokens.}
- \item{original.data}{a logical, if \code{TRUE} (the default) the
- data in the original NEXUS file are eventually written in
- \code{"file"} (see below).}
+ \item{original.data}{deprecated; will be removed soon.}
}
\description{
This function writes trees in a file with the NEXUS format.
}
\details{
- If \code{original.data = TRUE}, the file as specified by the attribute
- \code{"origin"} of the first tree is read and its data (except the
- trees) are written in \code{file}.
-
If several trees are given, they must have all the same tip labels.
If among the objects given some are not trees of class \code{"phylo"},
-/* BIONJ.c 2008-01-14 */
+/* BIONJ.c 2012-02-09 */
/* Copyright 2007-2008 Olivier Gascuel, Hoa Sien Cuong,
R port by Vincent Lefort, bionj() below modified by Emmanuel Paradis */
name=(WORD *)calloc(1,sizeof(WORD)); /* taxon name is */
if(name == NULL) /* put in trees */
{
- Rprintf("Out of memories !!");
- exit(0);
+ error("out of memory");
}
else
{
delta[i]=(float *)calloc(n+1, sizeof(float));
if(delta[i] == NULL)
{
- Rprintf("Out of memories!!");
- exit(0);
+ error("out of memory");
}
}
trees=(POINTERS *)calloc(n+1,sizeof(POINTERS));
if(trees == NULL)
{
- Rprintf("Out of memories!!");
- exit(0);
+ error("out of memory");
}
/* initialise and symmetrize the running delta matrix */
r=n;
bran=(WORD *)calloc(1,sizeof(WORD));
if(bran == NULL)
{
- Rprintf("Out of memories");
- exit(0);
+ error("out of memory");
}
else
{
-/* me.c 2008-01-14 */
+/* me.c 2012-02-09 */
/* Copyright 2007-2008 Olivier Gascuel, Rick Desper,
R port by Vincent Lefort, me_*() below modified by Emmanuel Paradis */
return(T);
if (NULL != v->parentEdge)
{
- Rprintf ("Error: root %s is poorly rooted.\n",v->label);
- exit(0);
+ error("root %s is poorly rooted.", v->label);
}
for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
{
}
if (-1 == v->index2)
{
- Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
- exit(0);
+ error("leaf %s in tree not in distance matrix.", v->label);
}
e = depthFirstTraverse(T,NULL);
while (NULL != e)
v = e->head;
if ((leaf(v)) && (-1 == v->index2))
{
- Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
- exit(0);
+ error("leaf %s in tree not in distance matrix.", v->label);
}
e = depthFirstTraverse(T,e);
}
for(X = S; NULL != X; X = X->secondNode)
if (X->firstNode->index2 > -1)
{
- Rprintf("Error node %s in matrix but not a leaf in tree.\n",X->firstNode->label);
- exit(0);
+ error("node %s in matrix but not a leaf in tree.", X->firstNode->label);
}
return;
}
for(i=0;i<d;i++)
{
for(j=0;j<d;j++)
- printf("%lf ", A[i][j]);
- printf("\n");
+ Rprintf("%lf ", A[i][j]);
+ Rprintf("\n");
}
}
-/* newick.c 2010-11-23 */
+/* newick.c 2012-02-09 */
/* Copyright 2007-2008 Vincent Lefort */
case(ReadOpenParenthesis):
if('(' != treeString[0])
{
- Rprintf("Error reading subtree.\n");
- exit(0);
+ error("error reading subtree");
}
i++;
state = ReadSubTree;
centerNode->middleEdge = thisEdge;
else
{
- Rprintf("Error: node %s has too many (>3) children.\n",centerNode->label);
- exit(0);
+ error("node %s has too many (>3) children.", centerNode->label);
}
//sprintf(thisEdge->label,"E%d",edgeCount++);
//snprintf(thisEdge->label,MAX_LABEL_LENGTH,"E%d",edgeCount++);
if ('(' != str[0])
{
- Rprintf("Error reading generated tree - does not start with '('.\n");
- exit(0);
+ error("generated tree does not start with '('");
}
inputLength = strlen (str)+1;
for(i = 0; i < inputLength; i++)
}
else if (parCount < 0)
{
- Rprintf("Error reading generated tree. Too many right parentheses.\n");
- exit(0);
+ error("generated tree has too many right parentheses");
}
}
centerNode = decodeNewickSubtree (str, T, &uCount);
projects / ape.git / commitdiff