From: paradis Date: Tue, 16 Mar 2010 09:54:18 +0000 (+0000) Subject: fixes in compar.ou and plot.phylo + new bind.tree X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=90f18c75d642f56b020bc6e0cdd0c5949c1d9a1d fixes in compar.ou and plot.phylo + new bind.tree git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@114 6e262413-ae40-0410-9e79-b911bd7a66b7 --- diff --git a/ChangeLog b/ChangeLog index f59984b..c551ded 100644 --- a/ChangeLog +++ b/ChangeLog @@ -5,8 +5,11 @@ NEW FEATURES 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 @@ -18,6 +21,14 @@ 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 diff --git a/DESCRIPTION b/DESCRIPTION index ce0536b..c23ac10 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ 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 diff --git a/R/bind.tree.R b/R/bind.tree.R index 8810f41..4202e2e 100644 --- a/R/bind.tree.R +++ b/R/bind.tree.R @@ -1,4 +1,4 @@ -## bind.tree.R (2010-02-12) +## bind.tree.R (2010-03-15) ## Bind Trees @@ -7,127 +7,202 @@ ## 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 } diff --git a/R/compar.ou.R b/R/compar.ou.R index ada36d7..7732b31 100644 --- a/R/compar.ou.R +++ b/R/compar.ou.R @@ -1,8 +1,8 @@ -## 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. @@ -10,27 +10,33 @@ 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] @@ -40,22 +46,30 @@ compar.ou <- function(x, phy, node = NULL, alpha = NULL) } 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")) diff --git a/R/plot.phylo.R b/R/plot.phylo.R index e02ba30..de7278e 100644 --- a/R/plot.phylo.R +++ b/R/plot.phylo.R @@ -1,4 +1,4 @@ -## plot.phylo.R (2010-01-04) +## plot.phylo.R (2010-03-15) ## Plot Phylogenies @@ -431,7 +431,7 @@ phylogram.plot <- function(edge, Ntip, Nnode, xx, yy, horizontal, 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) diff --git a/man/bind.tree.Rd b/man/bind.tree.Rd index ccc5714..ff8cafe 100644 --- a/man/bind.tree.Rd +++ b/man/bind.tree.Rd @@ -2,7 +2,8 @@ \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"}.} @@ -13,6 +14,8 @@ bind.tree(x, y, where = "root", position = 0) \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 @@ -20,21 +23,22 @@ bind.tree(x, y, where = "root", position = 0) } \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}} } @@ -55,13 +59,45 @@ cat("(Turniciformes:27.0,(Piciformes:26.3,((Galbuliformes:24.4,", 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} diff --git a/man/compar.ou.Rd b/man/compar.ou.Rd index b6f4b5c..4c1c02b 100644 --- a/man/compar.ou.Rd +++ b/man/compar.ou.Rd @@ -9,8 +9,9 @@ compar.ou(x, phy, node = NULL, alpha = NULL) 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).} } @@ -76,13 +77,12 @@ compar.ou(x, phy, node = NULL, alpha = NULL) 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: @@ -93,11 +93,10 @@ compar.ou(rnorm(23), bird.orders, alpha = 0.1) ### 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}