o The new function stree generates trees with regular shapes.
- o drop.tip(), extract.clade(), and root() now have an 'interactive'
- option to make the operation on a plotted tree.
+ o It is now possible to bind two trees with x + y (see ?bind.tree for
+ details).
+
+ o drop.tip(), extract.clade(), root(), and bind.tree() now have an
+ 'interactive' option to make the operation on a plotted tree.
BUG FIXES
with NA values. Nothing is drawn now like with 'text' or 'pch'.
The same bug occurred with the 'pie' option.
+ o A bug was fixed in compar.ou() and the help page was clarified.
+
+ o bind.tree() has been rewritten fixing several bugs and making it
+ more efficient.
+
+ o plot.phylo(type = "p") sometimes failed to colour correctly the
+ vertical lines representing the nodes.
+
CHANGES IN APE VERSION 2.5
Package: ape
Version: 2.5-1
-Date: 2010-03-12
+Date: 2010-03-15
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>
-## bind.tree.R (2010-02-12)
+## bind.tree.R (2010-03-15)
## Bind Trees
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-bind.tree <- function(x, y, where = "root", position = 0)
+`+.phylo` <- function(x, y)
{
- nb.tip <- length(x$tip.label)
- nb.node <- x$Nnode
- ROOT <- nb.tip + 1
- if (where == 0 || where == "root")
- where <- ROOT
- if (position < 0) position <- 0
- if (where > nb.tip + nb.node) stop("node number out of range for tree 'x'")
- nb.edge <- dim(x$edge)[1]
- yHasNoRootEdge <- is.null(y$root.edge)
- xHasNoRootEdge <- is.null(x$root.edge)
+ p <- if (is.null(x$root.edge)) 0 else x$root.edge
+ bind.tree(x, y, position = p)
+}
- ## check whether both trees have branch lengths:
- wbl <- TRUE
- noblx <- is.null(x$edge.length)
- nobly <- is.null(y$edge.length)
- if (noblx && nobly) wbl <- FALSE
- if (xor(noblx, nobly)) {
- if (nobly) x$edge.length <- NULL
- else y$edge.length <- NULL
- wbl <- FALSE
- warning("one tree has no branch lengths, they will be ignored")
- }
+bind.tree <- function(x, y, where = "root", position = 0, interactive = FALSE)
+{
+ nx <- length(x$tip.label)
+ mx <- x$Nnode
+ ROOTx <- nx + 1L
+ ny <- length(y$tip.label)
+ my <- y$Nnode
- ## To avoid problems with tips or nodes with indentical
- ## labels we substitute the one where `y' is grafted:
- if (where <= nb.tip) {
- Tip.Label.where <- x$tip.label[where]
- x$tip.label[where] <- "TheTipWhereToGraftY"
+ if (interactive) {
+ lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
+ if (lastPP$type != "phylogram" || lastPP$direction != "rightwards")
+ stop("you must plot tree 'x' as a 'rightward phylogram'")
+ cat("Click where you want to graft tree 'y'...\n")
+ xy <- locator(1)
+ d <- abs(xy$y - lastPP$yy)
+ d[lastPP$xx - xy$x < 0] <- Inf
+ where <- which.min(d)
+ position <- lastPP$xx[where] - xy$x
+ if (position < 0) position <- 0
+ cat("The following parameters are used:\n")
+ cat(" where =", where, " position =", position, "\n")
+ } else {
+ if (where == 0 || where == "root") where <- ROOTx
+ if (position < 0) position <- 0
+ if (where > nx + mx)
+ stop("argument 'where' out of range for tree 'x'")
}
- if (where > ROOT) {
- xHasNoNodeLabel <- TRUE
- if (is.null(x$node.label)) {
- x$node.label <- paste("NODE", 1:nb.node, sep = "")
- x$node.label[where - nb.tip] <- "TheNodeWhereToGraftY"
- } else {
- Node.Label.where <- x$node.label[where - nb.tip]
- x$node.label[where - nb.tip] <- "TheNodeWhereToGraftY"
- xHasNoNodeLabel <- FALSE
- }
+
+ ## check whether both trees have branch lengths:
+ switch(is.null(x$edge.length) + is.null(y$edge.length) + 1L,
+ wbl <- TRUE, {
+ x$edge.length <- y$edge.length <- NULL
+ wbl <- FALSE
+ warning("one tree has no branch lengths, they will be ignored")
+ },
+ wbl <- FALSE)
+
+ yHasNoRootEdge <- is.null(y$root.edge)
+ xHasNoRootEdge <- is.null(x$root.edge)
+
+ ## find the row of 'where' before renumbering
+ if (where == ROOTx) case <- 1 else {
+ i <- which(x$edge[, 2] == where)
+ case <- if (where <= nx) 2 else 3
}
+ ## case = 1 -> y is bound on the root of x
+ ## case = 2 -> y is bound on a tip of x
+ ## case = 3 -> y is bound on a node of x
- ## if we bind `y' under a node or tip of `y', we first
- ## adjust the edge lengths if needed
- if (position && wbl) {
- if (where == ROOT) {
- if (xHasNoRootEdge) stop("tree 'x' has no root edge")
- if (x$root.edge < position)
- stop("argument 'position' is larger than the root edge.")
- x$root.edge <- x$root.edge - position
+ ## check that 'position' is correct
+ if (position) {
+ if (!wbl)
+ stop("'position' is non-null but trees have no branch lengths")
+ if (case == 1) {
+ if (xHasNoRootEdge)
+ stop("tree 'x' has no root edge")
+ if (position > x$root.edge)
+ stop("'position' is larger than x's root edge")
} else {
- i <- which(x$edge[, 2] == where)
if (x$edge.length[i] < position)
- stop("argument 'position' is larger than the specified edge.")
- x$edge.length[i] <- x$edge.length[i] - position
+ stop("'position' is larger than the branch length")
}
- y$root.edge <- if (yHasNoRootEdge) position else y$root.edge + position
}
- if (is.null(y$root.edge) && where > nb.tip) y$root.edge <- 0
+ x <- reorder(x)
+ y <- reorder(y)
- X <- write.tree(x)
- Y <- write.tree(y)
- Y <- substr(Y, 1, nchar(Y) - 1)
+### because in all situations internal nodes need to be renumbered,
+### they are changed to negatives first, and nodes that eventually
+### will be added will be numbered sequentially
- if (where <= nb.tip) {
- if (position)
- X <- gsub("TheTipWhereToGraftY",
- paste("(", "TheTipWhereToGraftY", ",", Y, ")",
- sep = ""), X)
- else
- X <- gsub("TheTipWhereToGraftY", Y, X)
- }
- if (where == ROOT) {
- rmvx <- if (xHasNoRootEdge) "\\);$" else ";$"
- X <- gsub(rmvx, "", X)
- Y <- gsub("^\\(", "", Y)
- if (!xHasNoRootEdge) X <- paste("(", X, sep = "")
- X <- paste(X, ",", Y, ";", sep = "")
+ nodes <- x$edge > nx
+ x$edge[nodes] <- -(x$edge[nodes] - nx) # -1, ..., -mx
+ nodes <- y$edge > ny
+ y$edge[nodes] <- -(y$edge[nodes] - ny + mx) # -(mx+1), ..., -(mx+my)
+ ROOT <- -1L # may change later
+ next.node <- -(mx + my) - 1L
+
+ ## renumber now the tips in y:
+ new.nx <- if (where <= nx && !position) nx - 1L else nx
+ y$edge[!nodes] <- y$edge[!nodes] + new.nx
+
+ ## if 'y' as a root edge, use it:
+ if (!yHasNoRootEdge) {
+ y$edge <- rbind(c(0, y$edge[1]), y$edge)
+ ## ^ will be filled later
+ next.node <- next.node - 1L
+ if (wbl) y$edge.length <- c(y$root.edge, y$edge.length)
}
- if (where > ROOT) {
+
+ switch(case, { # case = 1
+ if (position) {
+ x$root.edge <- x$root.edge - position
+ x$edge <- rbind(c(next.node, x$edge[1]), x$edge)
+ ROOT <- next.node
+ if (wbl) x$edge.length <- c(position, x$edge.length)
+ }
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- ROOT
+ } else y$edge[1] <- ROOT
+ x$edge <- rbind(x$edge, y$edge)
+ if (wbl)
+ x$edge.length <- c(x$edge.length, y$edge.length)
+ }, { # case = 2
if (position) {
- ## find where is the node in `X':
- ## below 19 is: nchar("TheNodeWhereToGraftY") - 1
- for (i in 1:nchar(X)) {
- if ("TheNodeWhereToGraftY" == substr(X, i, i + 19))
- break
- i <- i + 1
+ x$edge[i, 2] <- next.node
+ x$edge <- rbind(x$edge[1:i, ], c(next.node, where), x$edge[-(1:i), ])
+ if (wbl) {
+ x$edge.length[i] <- x$edge.length[i] - position
+ x$edge.length <- c(x$edge.length[1:i], position, x$edge.length[-(1:i)])
}
- ## now go back to find the left matching parentheses
- n.paren <- 1
- i <- i - 2
- while (n.paren > 0) {
- if (substr(X, i, i) == ")") n.paren <- n.paren + 1
- if (substr(X, i, i) == "(") n.paren <- n.paren - 1
- i <- i - 1
+ i <- i + 1L
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- x$edge[i, 1]
+ } else y$edge[1] <- x$edge[i, 1]
+ } else {
+ if (yHasNoRootEdge) x$edge[i, 2] <- y$edge[1]
+ else {
+ ## the root edge of y is fused with the terminal edge of x
+ if (wbl) y$edge.length[1] <- y$edge.length[1] + x$edge.length[i]
+ y$edge[1] <- x$edge[i, 1]
+ ## delete i-th edge in x:
+ x$edge <- x$edge[-i, ]
+ i <- i - 1L
+ if (wbl) x$edge.length <- x$edge.length[-i]
+ }
+ x$tip.label <- x$tip.label[-where]
+ ## renumber the tips that need to:
+ ii <- which(x$edge[, 2] > where & x$edge[, 2] <= nx)
+ x$edge[ii, 2] <- x$edge[ii, 2] - 1L
+ }
+ x$edge <- rbind(x$edge[1:i, ], y$edge, x$edge[-(1:i), ])
+ if (wbl)
+ x$edge.length <- c(x$edge.length[1:i], y$edge.length, x$edge.length[-(1:i)])
+ }, { # case = 3
+ if (position) {
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- next.node
+ } else y$edge[1] <- next.node
+ x$edge <- rbind(x$edge[1:i, ], c(next.node, x$edge[i, 2]), x$edge[-(1:i), ])
+ x$edge[i, 2] <- next.node
+ if (wbl) {
+ x$edge.length[i] <- x$edge.length[i] - position
+ x$edge.length <- c(x$edge.length[1:i], position, x$edge.length[-(1:i)])
}
- ## insert the left parenthesis:
- ## here 21 is: nchar("TheNodeWhereToGraftY") + 1
- X <- paste(substr(X, 1, i - 1), "(",
- substr(X, i, 21), sep = "")
- ## and insert `y':
- X <- gsub("TheNodeWhereToGraftY",
- paste("TheNodeWhereToGraftY", ",", Y,
- sep = ""), X)
+ i <- i + 1L
} else {
- xx <- paste(")", "TheNodeWhereToGraftY", sep = "")
- X <- gsub(xx, paste(",", Y, xx, sep = ""), X)
+ if (yHasNoRootEdge) {
+ j <- which(y$edge[, 1] == y$edge[1])
+ y$edge[j, 1] <- x$edge[i, 2]
+ } else y$edge[1] <- x$edge[i, 2]
}
+ x$edge <- rbind(x$edge[1:i, ], y$edge, x$edge[-(1:i), ])
+ if (wbl)
+ x$edge.length <- c(x$edge.length[1:i], y$edge.length, x$edge.length[-(1:i)])
+ })
+
+ x$tip.label <- c(x$tip.label, y$tip.label)
+
+ if (is.null(x$node.label)) {
+ if (!is.null(y$node.label))
+ x$node.label <- c(rep(NA, mx), y$node.label)
+ } else {
+ x$node.label <-
+ if (is.null(y$node.label)) c(x$node.label, rep(NA, my))
+ else c(x$node.label, y$node.label)
}
- phy <- read.tree(text = X)
- ## restore the labels:
- if (where <= nb.tip)
- phy$tip.label[which(phy$tip.label == "TheTipWhereToGraftY")] <-
- Tip.Label.where
- if (where > ROOT) {
- if (xHasNoNodeLabel) phy$node.label <- NULL
- else
- phy$node.label[which(phy$node.label == "TheNodeWhereToGraftY")] <-
- Node.Label.where
- }
- phy
+
+ n <- length(x$tip.label)
+ x$Nnode <- dim(x$edge)[1] + 1L - n
+
+ ## update the node labels before renumbering (this adds NA for
+ ## the added nodes, and drops the label for those deleted)
+ if (!is.null(x$node.label))
+ x$node.label <- x$node.label[-unique(x$edge[, 1])]
+
+ ## renumber nodes:
+ newNb <- integer(x$Nnode)
+ newNb[-ROOT] <- n + 1L
+ sndcol <- x$edge[, 2] < 0
+ ## executed from right to left, so newNb is modified before x$edge:
+ x$edge[sndcol, 2] <- newNb[-x$edge[sndcol, 2]] <-
+ (n + 2):(n + x$Nnode)
+ x$edge[, 1] <- newNb[-x$edge[, 1]]
+
+ if (!is.null(x$node.label))
+ x$node.label <- x$node.label[newNb[newNb == 0] - n]
+
+ x
}
-## compar.ou.R (2009-05-10)
+## compar.ou.R (2010-03-15)
## Ornstein--Uhlenbeck Model for Continuous Characters
-## Copyright 2005-2009 Emmanuel Paradis
+## Copyright 2005-2010 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
compar.ou <- function(x, phy, node = NULL, alpha = NULL)
{
if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo".')
+ stop('object "phy" is not of class "phylo".')
if (!is.numeric(x)) stop("'x' must be numeric.")
if (!is.null(names(x))) {
if (all(names(x) %in% phy$tip.label)) x <- x[phy$tip.label]
else warning('the names of argument "x" and the tip labels of the tree did not match: the former were ignored in the analysis.')
}
- nb.tip <- length(phy$tip.label)
- root <- nb.tip + 1
+ n <- length(phy$tip.label)
+ root <- n + 1L
if (is.null(node)) node <- numeric(0)
- if (root %in% node) node <- node[node != root]
+ if (is.character(node)) {
+ if (is.null(phy$node.label))
+ stop("argument 'node' is character but 'phy' has no node labels")
+ node <- match(node, phy$node.label) + n
+ phy$node.label <- NULL
+ }
+ if (root %in% node) node <- node[-1]
bt <- branching.times(phy)
Tmax <- bt[1]
- Wend <- matrix(0, nb.tip, length(node) + 1)
- colnames(Wend) <- c(names(sort(bt[node])), as.character(root))
+ Wend <- matrix(0, n, length(node) + 1)
+ colnames(Wend) <- c(names(sort(bt[node - n])), as.character(root))
Wstart <- Wend
Wstart[, ncol(Wstart)] <- Tmax
- root2tip <- .Call("seq_root2tip", phy$edge, nb.tip,
+ root2tip <- .Call("seq_root2tip", phy$edge, n,
phy$Nnode, PACKAGE = "ape")
- for (i in 1:nb.tip) {
+ for (i in 1:n) {
last.change <- names(Tmax)
- for (j in root2tip[[i]]) {#[-1]) {# don't need to look at the root
+ for (j in root2tip[[i]]) {
if (j %in% node) {
jb <- as.character(j)
Wend[i, last.change] <- Wstart[i, jb] <- bt[jb]
}
W <- cophenetic.phylo(phy)
dev <- function(p) {
- ##M <- rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
- ##M <- -rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
- M <- rowSums((exp(-p[1] * Wend) - exp(-p[1] * Wstart)) * p[-(1:2)])
- V <- exp(-p[1]*W) * (1 - exp(-2*p[1]*(Tmax - W/2)))
- nb.tip*log(2*pi*p[2]) + log(det(V)) +
- (t(x - M) %*% chol2inv(V) %*% (x - M)) / p[2]
+ alpha <- p[1]
+ sigma2 <- p[2]
+ if (sigma2 <= 0) return(1e100)
+ theta <- p[-(1:2)]
+ ## 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
}
- if (is.null(alpha))
- out <- nlm(function(p) dev(p),
- p = c(0.1, 1, rep(mean(x), ncol(Wstart))),
- hessian = TRUE)
- else
- out <- nlm(function(p) dev(c(alpha, p)),
- p = c(1, rep(mean(x), ncol(Wstart))),
- hessian = TRUE)
- para <- cbind(out$estimate, sqrt(diag(solve(out$hessian))))
+
+ out <- if (is.null(alpha))
+ nlm(function(p) dev(p),
+ p = c(0.1, 1, rep(mean(x), ncol(Wstart))), hessian = TRUE)
+ else nlm(function(p) dev(c(alpha, p)),
+ p = c(1, rep(mean(x), ncol(Wstart))), hessian = TRUE)
+
+ ## if alpha is estimated it may be that the Hessian matrix has the
+ ## corresponding column and row filled with 0, making solve() fail
+ se <- if (is.null(alpha) && all(out$hessian[1, ] == 0))
+ c(NA, sqrt(diag(solve(out$hessian[-1, -1]))))
+ else sqrt(diag(solve(out$hessian)))
+
+ para <- cbind(out$estimate, se)
nms <- c("sigma2", paste("theta", 1:ncol(Wstart), sep = ""))
if (is.null(alpha)) nms <- c("alpha", nms)
dimnames(para) <- list(nms, c("estimate", "stderr"))
-## plot.phylo.R (2010-01-04)
+## plot.phylo.R (2010-03-15)
## Plot Phylogenies
edge.color <- rep(edge.color, length.out = Nedge)
edge.width <- rep(edge.width, length.out = Nedge)
edge.lty <- rep(edge.lty, length.out = Nedge)
- DF <- data.frame(edge.color, edge.width, edge.lty)
+ DF <- data.frame(edge.color, edge.width, edge.lty, stringsAsFactors = FALSE)
color.v <- rep("black", Nnode)
width.v <- rep(1, Nnode)
lty.v <- rep(1, Nnode)
\alias{bind.tree}
\title{Binds Trees}
\usage{
-bind.tree(x, y, where = "root", position = 0)
+bind.tree(x, y, where = "root", position = 0, interactive = FALSE)
+\method{+}{phylo}(x, y)
}
\arguments{
\item{x}{an object of class \code{"phylo"}.}
\item{position}{a numeric value giving the position from the tip or
node given by \code{node} where the tree \code{y} is binded;
negative values are ignored.}
+ \item{interactive}{if \code{TRUE} the user is asked to choose the tip
+ or node of \code{x} by clicking on the tree which must be plotted.}
}
\description{
This function binds together two phylogenetic trees to give a single
}
\details{
The argument \code{x} can be seen as the receptor tree, whereas
- \code{y} is the donor tree. The root of \code{y} is then sticked on a
+ \code{y} is the donor tree. The root of \code{y} is then grafted on a
location of \code{x} specified by \code{where} and, possibly,
\code{position}. If \code{y} has a root edge, this is added as in
internal branch in the resulting tree.
+
+ \code{x + y} is a shortcut for:
+
+ \preformatted{
+ bind.tree(x, y, position = if (is.null(x$root.edge)) 0 else
+ x$root.edge)
+ }
}
\value{
an object of class \code{"phylo"}.
}
-\note{
- For the moment, this function handles only trees with branch lengths,
- and does not handle node labels.
-
- Further testing/improvements may be needed.
-}
-\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
+\author{Emmanuel Paradis}
\seealso{
\code{\link{drop.tip}}, \code{\link{root}}
}
tree.bird1 <- read.tree("ex1.tre")
tree.bird2 <- read.tree("ex2.tre")
unlink(c("ex1.tre", "ex2.tre")) # clean-up
-birds <- bind.tree(tree.bird1, tree.bird2, where = "root",
- position = tree.bird1$root.edge)
-birds
+(birds <- tree.bird1 + tree.bird2)
layout(matrix(c(1, 2, 3, 3), 2, 2))
plot(tree.bird1)
plot(tree.bird2)
plot(birds)
-layout(matrix(1))
+
+### examples with random trees
+x <- rtree(4, tip.label = LETTERS[1:4])
+y <- rtree(4, tip.label = LETTERS[5:8])
+x <- makeNodeLabel(x, prefix = "x_")
+y <- makeNodeLabel(y, prefix = "y_")
+x$root.edge <- y$root.edge <- .2
+
+z <- bind.tree(x, y, po=.2)
+plot(y, show.node.label = TRUE, font = 1, root.edge = TRUE)
+title("y")
+plot(x, show.node.label = TRUE, font = 1, root.edge = TRUE)
+title("x")
+plot(z, show.node.label = TRUE, font = 1, root.edge = TRUE)
+title("z <- bind.tree(x, y, po=.2)")
+
+z <- bind.tree(x, y, 2, .1)
+plot(y, show.node.label = TRUE, font = 1, root.edge = TRUE)
+title("y")
+plot(x, show.node.label = TRUE, font = 1, root.edge = TRUE)
+title("x")
+plot(z, show.node.label = TRUE, font = 1, root.edge = TRUE)
+title("z <- bind.tree(x, y, 2, .1)")
+
+x <- rtree(100)
+y <- rtree(100)
+x$root.edge <- y$root.edge <- .2
+z <- x + y
+plot(y, show.tip.label = FALSE, root.edge = TRUE); axisPhylo()
+title("y")
+plot(x, show.tip.label = FALSE, root.edge = TRUE); axisPhylo()
+title("x")
+plot(z, show.tip.label = FALSE, root.edge = TRUE); axisPhylo()
+title("z <- x + y")
+layout(1)
}
\keyword{manip}
character.}
\item{phy}{an object of class \code{"phylo"}.}
\item{node}{a vector giving the number(s) of the node(s) where the
- parameter `theta' (the character optimum) is assumed to change. By
- default there is no change (same optimum thoughout lineages).}
+ parameter `theta' (the trait optimum) is assumed to change. The
+ node(s) can be specified with their labels if \code{phy} has node
+ labels. By default there is no change (same optimum thoughout lineages).}
\item{alpha}{the value of \eqn{\alpha}{alpha} to be used when fitting
the model. By default, this parameter is estimated (see details).}
}
Hansen, T. F. (1997) Stabilizing selection and the comparative
analysis of adaptation. \emph{Evolution}, \bold{51}, 1341--1351.
}
-\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
+\author{Emmanuel Paradis}
\seealso{
\code{\link{ace}}, \code{\link{compar.lynch}},
\code{\link{corBrownian}}, \code{\link{corMartins}}, \code{\link{pic}}
}
\examples{
-\dontrun{
data(bird.orders)
### This is likely to give you estimates close to 0, 1, and 0
### for alpha, sigma^2, and theta, respectively:
### for the two clades of birds...
x <- c(rnorm(5, 0), rnorm(18, 5))
### ... the model with two optima:
-compar.ou(x, bird.orders, node = -2, alpha = .1)
+compar.ou(x, bird.orders, node = 25, alpha = .1)
### ... and the model with a single optimum:
compar.ou(x, bird.orders, node = NULL, alpha = .1)
### => Compare both models with the difference in deviances
-## with follows a chi^2 with df = 1.
-}
+## wicth follows a chi^2 with df = 1.
}
\keyword{models}
projects / ape.git / commitdiff