diff --git a/.Rbuildignore b/.Rbuildignore index bd4c0f2..7dfa56a 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -16,6 +16,7 @@ lib library check www +examples ^(.+?)\.pdf ^(.+?)\.tar\.gz ^(.+?)\.tgz diff --git a/DESCRIPTION b/DESCRIPTION index e7b5996..d31f922 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,13 +2,37 @@ Package: ouch Type: Package Title: Ornstein-Uhlenbeck Models for Phylogenetic Comparative Hypotheses Version: 2.15-1 -Date: 2021-02-16 -Author: Aaron A. King and Marguerite A. Butler +Date: 2021-02-17 +Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa@umich.edu"), + person(given=c("Marguerite","A."),family="Butler",role=c("ctb"))) Maintainer: Aaron A. King Description: Fit and compare Ornstein-Uhlenbeck models for evolution along a phylogenetic tree. -Depends: R(>= 3.6), methods, stats, graphics, subplex +Depends: R(>= 3.6) +Imports: methods, stats, graphics, grDevices, utils, subplex Suggests: ape, geiger URL: https://kingaa.github.io/ouch/ License: GPL-3 LazyLoad: true -Collate: ouchtree.R ape2ouch.R glssoln.R rmvnorm.R brown.R hansen.R methods.R paint.R plot.R +LazyData: true +RoxygenNote: 7.1.1 +BugReports: https://github.com/kingaa/ouch/issues/ +Encoding: UTF-8 +Collate: + 'anolis.R' + 'package.R' + 'ouchtree.R' + 'ape2ouch.R' + 'bimac.R' + 'rmvnorm.R' + 'glssoln.R' + 'brown.R' + 'simulate.R' + 'hansen.R' + 'bootstrap.R' + 'coef.R' + 'loglik.R' + 'paint.R' + 'plot.R' + 'print.R' + 'summary.R' + 'update.R' diff --git a/Makefile b/Makefile index 480b05f..337dfb0 100644 --- a/Makefile +++ b/Makefile @@ -2,6 +2,7 @@ REXE = R --vanilla RCMD = $(REXE) CMD RCMD_ALT = R --no-save --no-restore CMD RSCRIPT = Rscript --vanilla +REPODIR = ../www PDFLATEX = pdflatex BIBTEX = bibtex @@ -18,6 +19,7 @@ PKGVERS = $(PKG)_$(VERSION) SOURCE=$(shell ls R/*R src/*.c src/*.h data/*) CSOURCE=$(shell ls src/*.c) TESTS=$(shell ls tests/*R) +REVDEPS=surface mvSLOUCH pmc default: @echo $(PKGVERS) @@ -25,24 +27,19 @@ default: .PHONY: clean win wind tests check dist manual vignettes: export R_QPDF=qpdf -dist manual vignettes: export R_GSCMD=gs -dist manual vignettes: export GS_QUALITY=ebook -dist manual vignettes: export R_HOME=$(shell $(REXE) RHOME) +roxy dist manual vignettes: export R_HOME=$(shell $(REXE) RHOME) check xcheck xxcheck: export FULL_TESTS=yes -xcheck tests: export R_PROFILE_USER=$(CURDIR)/.Rprofile -session htmldocs vignettes data tests manual: export R_LIBS=$(CURDIR)/library -session: export R_DEFAULT_PACKAGES=datasets,utils,grDevices,graphics,stats,methods,ouch -xxcheck: export R_LIBS=$(CURDIR)/check - -includes: - -headers: - -inst/include/%.h: src/%.h - $(CP) $^ $@ +dist revdeps session tests check xcheck xxcheck: export R_KEEP_PKG_SOURCE=yes +revdeps xcheck tests: export R_PROFILE_USER=$(CURDIR)/.Rprofile +revdeps session xxcheck htmldocs vignettes data tests manual: export R_LIBS=$(CURDIR)/library +session: export R_DEFAULT_PACKAGES=datasets,utils,grDevices,graphics,stats,methods,tidyverse,subplex,ouch htmldocs: inst/doc/*.html +htmlhelp: install + rsync -avz library/ouch/html/ www/manual + (cd www/manual; (cat links.ed && echo w ) | ed - 00Index.html) + vignettes: manual install $(MAKE) -C www/vignettes @@ -59,13 +56,17 @@ www/NEWS.html: inst/NEWS.Rd session: install exec $(REXE) +revdeps: install + mkdir -p library check + $(REXE) -e "pkgs <- strsplit('$(REVDEPS)',' ')[[1]]; download.packages(pkgs,destdir='library',repos='https://mirrors.nics.utk.edu/cran/')" + $(RCMD) check --library=library -o check library/*.tar.gz + roxy: $(SOURCE) -## $(REXE) -e "pkgload::load_all(compile=FALSE); devtools::document(roclets=c('rd','collate','namespace'))" $(REXE) -e "pkgbuild::compile_dll(); devtools::document(roclets=c('rd','collate','namespace'))" dist: NEWS $(PKGVERS).tar.gz -$(PKGVERS).tar.gz: $(SOURCE) $(TESTS) includes headers +$(PKGVERS).tar.gz: $(SOURCE) $(TESTS) $(RCMD) build --force --no-manual --resave-data --compact-vignettes=both --md5 . binary: dist @@ -74,9 +75,9 @@ binary: dist rm -rf plib publish: dist manual news - $(RSCRIPT) -e 'drat::insertPackage("$(PKGVERS).tar.gz",repodir="../www",action="prune")' - -$(RSCRIPT) -e 'drat::insertPackage("$(PKGVERS).tgz",repodir="../www",action="prune")' - -$(RSCRIPT) -e 'drat::insertPackage("$(PKGVERS).zip",repodir="../www",action="prune")' + $(RSCRIPT) -e 'drat::insertPackage("$(PKGVERS).tar.gz",repodir="$(REPODIR)",action="prune")' + -$(RSCRIPT) -e 'drat::insertPackage("$(PKGVERS).tgz",repodir="$(REPODIR)",action="prune")' + -$(RSCRIPT) -e 'drat::insertPackage("$(PKGVERS).zip",repodir="$(REPODIR)",action="prune")' $(CP) $(PKG).pdf ../www/manuals win: dist @@ -98,10 +99,10 @@ qqcheck: dist $(RCMD) check --library=check -o check --no-codoc --no-examples --no-vignettes --no-manual --no-tests $(PKGVERS).tar.gz xcheck: dist - mkdir -p check + mkdir -p check library $(RCMD_ALT) check --no-stop-on-test-error --as-cran --library=library -o check $(PKGVERS).tar.gz -xxcheck: xcheck +xxcheck: install xcheck mkdir -p check $(REXE) -d "valgrind --tool=memcheck --track-origins=yes --leak-check=full" < check/$(PKG).Rcheck/$(PKG)-Ex.R 2>&1 | tee $(PKG)-Ex.Rout @@ -123,10 +124,15 @@ install: library/$(PKG) library/$(PKG): dist mkdir -p library - $(RCMD) INSTALL --library=library $(PKGVERS).tar.gz + $(RCMD) INSTALL --html --library=library $(PKGVERS).tar.gz remove: - -$(RCMD) REMOVE --library=library $(PKG) + if [ -d library ]; then \ + $(RCMD) REMOVE --library=library $(PKG); \ + rmdir library; \ + fi + +fresh: clean remove inst/doc/*.html: install diff --git a/NAMESPACE b/NAMESPACE index 118e34f..747c9be 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,21 +1,32 @@ -useDynLib(ouch,.registration=TRUE) +# Generated by roxygen2: do not edit by hand + +export(ape2ouch) +export(brown) +export(hansen) +export(ouchtree) +export(paint) +exportMethods(bootstrap) +exportMethods(coef) +exportMethods(logLik) +exportMethods(plot) +exportMethods(print) +exportMethods(show) +exportMethods(simulate) +exportMethods(summary) +exportMethods(update) import(methods) -importFrom(subplex,subplex) -importFrom(stats,deviance,simulate,update,coef,logLik,optim, - rnorm,runif,setNames) -importFrom(utils,head,tail) -importFrom(graphics,plot,par,text) importFrom(grDevices,rainbow) -exportClasses('ouchtree','browntree','hansentree') -exportMethods( - 'plot','logLik','summary','coef','coerce', - 'print','show', - 'simulate','update','bootstrap' -) -export( - ouchtree, - brown, - hansen, - ape2ouch, - paint -) +importFrom(graphics,par) +importFrom(graphics,text) +importFrom(stats,coef) +importFrom(stats,logLik) +importFrom(stats,optim) +importFrom(stats,rnorm) +importFrom(stats,runif) +importFrom(stats,setNames) +importFrom(stats,simulate) +importFrom(stats,update) +importFrom(subplex,subplex) +importFrom(utils,head) +importFrom(utils,tail) +useDynLib(ouch, .registration = TRUE) diff --git a/R/anolis.R b/R/anolis.R new file mode 100644 index 0000000..6292236 --- /dev/null +++ b/R/anolis.R @@ -0,0 +1,43 @@ +#' Greater Antillean anolis lizard sexual size dimorphism data +#' +#' The dataset consists of sexual size-dimorphism data for 38 species of anoles from Cuba, Hispaniola, Jamaica, and Puerto Rico (Butler, Schoener, and Losos 2000). +#' Each of these species belongs to one of six microhabitat types, or \dQuote{ecomorphs} (sensu Williams, 1972): +#' trunk-ground, grass-bush, trunk, trunk-crown, twig, and crown-giant. +#' The data were used to demonstrate an evolutionary association between habitat type and degree of sexual size dimorphism. +#' +#' Size dimorphism was calcuated as the log-ratio of male snout-to-vent length to female snout-to-vent length (males are larger). +#' +#' In this example, we tested three models of evolution: +#' Brownian motion, Ornstein-Uhlenbeck with one global optimum, and Ornstein-Uhlenbeck with 7 optima (one for each ecomorph type plus an additional one for an \dQuote{unknown} type). +#' +#' For the 7-optima model, we assigned each terminal branch to an optimum according to the ecomorph type of the extant species. +#' Because we had no information to help guide hypotheses about internal branches, we assigned +#' internal branches to the \dQuote{unknown} selective regime. +#' The phylogeny of these species is consistent with and adaptive radiation, with a burst of speciation events early in the evolutionary history of this clade (see phylogeny in Butler & King (2004) or example below). +#' +#' @name anolis.ssd +#' @rdname anolis_ssd +#' @docType data +#' @format +#' A data frame with 38 observations on the following 6 variables. +#' \describe{ +#' \item{node}{Labels for the nodes.} +#' \item{species}{Names of extant species.} +#' \item{log.SSD}{Log sexual size dimorphism of extant species.} +#' \item{ancestor}{Ancestor node.} +#' \item{time}{Time of node.} +#' \item{OU.1}{a factor with levels \code{ns}} +#' \item{OU.7}{a factor with levels corresponding to ecomorph +#' (\code{tg} \code{tc} \code{gb} \code{cg} \code{tw} \code{tr} \code{anc})} +#' } +#' @author Marguerite A. Butler, Aaron A. King +#' @references +#' \Butler2000 +#' +#' \Williams1972 +#' +#' @source \Butler2004 +#' @keywords models +#' @example examples/anolis.R +#' +NULL diff --git a/R/ape2ouch.R b/R/ape2ouch.R index dbdd9ca..eb59da7 100644 --- a/R/ape2ouch.R +++ b/R/ape2ouch.R @@ -1,23 +1,37 @@ -# This function takes a tree file in the phylo format (from -# ape/read.tree) and creates an ouchtree object. -# The option exists whereby users can change branch lengths -# while keeping the same topology. -# -# t = object of type 'phylo', as returned by ape/read.tree -# branch.lengths = optional branch length vector in same -# order as t$edge.length; default is to use t$edge.length -# -# D. Ackerly, dackerly at berkeley dot edu, July 18, 2006. -# Modified by A.A. King, kingaa at umich dot edu, 5/15/2007. - +#' Convert an "ape" tree to an "ouch" tree. +#' +#' \code{ape2ouch} translates \pkg{ape}'s \code{phylo} representation of a phylogenetic tree into \pkg{ouch}'s \code{ouchtree} representation. +#' The user can change the branch lengths while preserving the topology. +#' +#' @param tree a tree of class \code{phylo} created in package \pkg{ape}. +#' @param scale if \code{scale=TRUE}, the tree's depth will be scaled to 1. +#' If \code{scale} is a number, then the branch lengths will be scaled by this number. +#' @param branch.lengths optional vector of branch lengths. +#' @author A. A. King, D. Ackerly +#' @keywords models +#' +#' @rdname ape2ouch +#' @include ouchtree.R +#' @export ape2ouch ape2ouch <- function (tree, scale = TRUE, branch.lengths = tree$edge.length) { + ## This function takes a tree file in the phylo format (from + ## ape/read.tree) and creates an ouchtree object. + ## The option exists whereby users can change branch lengths + ## while keeping the same topology. + ## + ## t = object of type 'phylo', as returned by ape/read.tree + ## branch.lengths = optional branch length vector in same + ## order as t$edge.length; default is to use t$edge.length + ## + ## D. Ackerly, July 18, 2006. + ## Modified by A.A. King, 5/15/2007. if (!inherits(tree,'phylo')) stop(sQuote("tree")," must be of class ",sQuote("phylo")) nnodes <- nrow(tree$edge)+1 # number of nodes n.term <- length(tree$tip.label) # number of terminal nodes - n.int <- nnodes-n.term # internal nodes + n.int <- nnodes-n.term # internal nodes tmp <- matrix(NA,nnodes,2) tmp[-1,1:2] <- tree$edge @@ -33,9 +47,9 @@ ape2ouch <- function (tree, scale = TRUE, branch.lengths = tree$edge.length) { for (n in 2:nnodes) { anc <- which(tmp[,2]==tmp[n,1]) if (length(anc)>1) stop("invalid tree") - if (length(anc)>0) { # the node has a non-root ancestor + if (length(anc)>0) { # the node has a non-root ancestor ancestor[n] <- node[anc] - } else { # the node has the root as an ancestor + } else { # the node has the root as an ancestor ancestor[n] <- node[1] } } @@ -57,11 +71,11 @@ ape2ouch <- function (tree, scale = TRUE, branch.lengths = tree$edge.length) { } ouchtree( - nodes=node, - ancestors=ancestor, - times=times, - labels=species - ) + nodes=node, + ancestors=ancestor, + times=times, + labels=species + ) } branch.height <- function (node, anc, bl, k) { diff --git a/R/bimac.R b/R/bimac.R new file mode 100644 index 0000000..2ec489e --- /dev/null +++ b/R/bimac.R @@ -0,0 +1,57 @@ +#' Anolis bimaculatus lizard size data. +#' +#' This is the \emph{Anolis bimaculatus} dataset used in Butler & King (2004). +#' It is used to test a hypothesis of character displacement using an interspecific dataset of body sizes and current data on sympatry/allopatry. +#' The data frame has the following columns: +#' \code{species} which are species names, +#' \code{size} which is the phenotypic data, +#' and the variables \code{ancestor} and \code{time} which specify the topology of the phylogeny and the location of the nodes in time, respectively. +#' The columns \code{OU.1}, \code{OU.3}, \code{OU.4}, and \code{OU.LP} specify four hypothetical arrangements of selective regimes. +#' Explanations of the data follow: +#' \describe{ +#' \item{Body size.}{ +#' We use the phenotypic data and phylogeny of Losos (1990), which employed the head lengths (of males) as a proxy for body size. +#' In this group of lizards, head length correlates very strongly with snout-to-vent length and the cube root of mass, which are standard measures of body size. +#' The data are head lengths in mm; note that we use the log of this value in analyses. +#' } +#' \item{Tree topology}{ +#' The tree topology is encoded via two vectors: \code{ancestor} and \code{time}. +#' Each node of the' phylogenetic tree has a corresponding row in the data frame, numbered from 1 to 45. +#' The columns \code{ancestor} and \code{time} specify the phylogeny. +#' The \code{ancestor} variable specifies the topology: it is a list indicating the ancestor of each node. +#' The root node has ancestor 0. +#' The variable \code{time} specifies the temporal location of each node, with the root node being at time 0. +#' } +#' \item{Specifications of selective regimes.}{ +#' (Columns \code{OU.1}, \code{OU.3}, \code{OU.4}, \code{OU.LP}). +#' These columns are factors, the levels of which correspond to the \dQuote{paintings} of the respective adaptive regime hypotheses onto the phylogeny. +#' Each selective regime is named (small, medium, large, etc.). +#' Each column corresponds to a different painting of the selective regimes, and thus to a different hypothesis. +#' In this example, there are 3 alternative models (see Butler & King 2004): \code{OU.4} is 4-regime model, \code{OU.3} is 3-regime model (all ancestors are medium), \code{OU.LP} is the linear parsimony model. +#' } +#' } +#' +#' @name bimac +#' @docType data +#' @format A data frame with 45 observations on the following 8 variables. +#' \describe{ +#' \item{node}{Labels for the nodes.} +#' \item{species}{Species names for extant species.} +#' \item{size}{Body size (head length in mm) of extant species.} +#' \item{ancestor}{Ancestral node.} +#' \item{time}{Time of node.} +#' \item{OU.1}{a factor with levels \code{ns}} +#' \item{OU.3}{a factor with levels \code{small}, \code{medium}, \code{large}} +#' \item{OU.4}{a factor with levels \code{small}, \code{medium}, \code{large}, \code{anc}} +#' \item{OU.LP}{a factor with levels \code{small}, \code{medium}, \code{large}} +#' } +#' @author Marguerite A. Butler and Aaron A. King +#' @references +#' \Lazell1972 +#' +#' \Losos1990 +#' @source \Butler2004 +#' @keywords models +#' @example examples/bimac.R +#' +NULL diff --git a/R/bootstrap.R b/R/bootstrap.R new file mode 100644 index 0000000..e0b9d78 --- /dev/null +++ b/R/bootstrap.R @@ -0,0 +1,75 @@ +#' Bootstrapping for uncertainty quantification +#' +#' Generic bootstrapping for \pkg{ouch} models. +#' +#' \code{bootstrap} performs a parametric bootstrap for estimation of confidence intervals. +#' +#' @rdname bootstrap +#' @name bootstrap +#' +#' @include brown.R hansen.R +#' +#' @param object A fitted model object. +#' @param nboot integer; number of bootstrap replicates. +#' @param seed integer; setting \code{seed} to a non-\code{NULL} value allows one to fix the random seed (see \code{\link[ouch:simulate]{simulate}}). +#' @param ... Additional arguments are passed to \code{\link[ouch:update]{update}}. +NULL + +setGeneric( + "bootstrap", + function (object, ...) { + standardGeneric("bootstrap") + } +) + +#' @rdname bootstrap +#' @export +setMethod( + "bootstrap", + signature=signature(object="browntree"), + function (object, nboot = 200, seed = NULL, ...) { + simdata <- simulate(object,nsim=nboot,seed=seed) + results <- vector(mode='list',length=nboot) + toshow <- c("sigma.squared","theta","loglik","aic","aic.c","sic","dof") + for (b in seq_len(nboot)) { + results[[b]] <- summary(update(object,data=simdata[[b]],...)) + } + as.data.frame(t(sapply(results,function(x)unlist(x[toshow])))) + } +) + +#' @rdname bootstrap +#' @export +setMethod( + "bootstrap", + signature=signature(object="hansentree"), + function (object, nboot = 200, seed = NULL, ...) { + simdata <- simulate(object,nsim=nboot,seed=seed) + results <- vector(mode='list',length=nboot) + toshow <- c("alpha","sigma.squared","optima","loglik","aic","aic.c","sic","dof") + for (b in seq_len(nboot)) { + results[[b]] <- summary(update(object,data=simdata[[b]],...)) + } + as.data.frame(t(sapply(results,function(x)unlist(x[toshow])))) + } +) + +#' @rdname bootstrap +#' @export +setMethod( + "bootstrap", + signature=signature(object="missing"), + definition=function (...) { + reqd_arg("bootstrap","object") + } +) + +#' @rdname bootstrap +#' @export +setMethod( + "bootstrap", + signature=signature(object="ANY"), + definition=function (object, ...) { + undef_method("bootstrap",object) + } +) diff --git a/R/brown.R b/R/brown.R index 80327bd..c6b9a88 100644 --- a/R/brown.R +++ b/R/brown.R @@ -1,3 +1,19 @@ +#' Phylogenetic Brownian motion models +#' +#' The function \code{brown} creates a \code{browntree} object by fitting a +#' Brownian-motion model to data. +#' +#' @name brown +#' @rdname brown +#' @author Aaron A. King +#' @seealso \code{\link{ouchtree}}, \code{\link{hansen}}, \code{\link{bimac}}, \code{\link{anolis.ssd}} +#' @keywords models +#' @references +#' \Butler2004 +#' @keywords models +#' +NULL + setClass( 'browntree', contains='ouchtree', @@ -11,6 +27,19 @@ setClass( ) ) +#' @rdname brown +#' @include ouchtree.R glssoln.R rmvnorm.R +#' +#' @param data Phenotypic data for extant species, i.e., at the terminal ends +#' of the phylogenetic tree. This can either be a numeric vector or a list. +#' If it is a numeric vector, there must be one entry for every node. If it is +#' a list, it must consist entirely of numeric vectors, each of which has one +#' entry per node. A data-frame is coerced to a list. +#' @param tree A phylogenetic tree, specified as an \code{ouchtree} object. +#' +#' @return \code{brown} returns an object of class \code{browntree}. +#' +#' @export brown <- function (data, tree) { if (!is(tree,'ouchtree')) @@ -93,9 +122,14 @@ brown.deviate <- function(n = 1, object) { apply(X,3,as.data.frame) } +#' @rdname simulate +#' @include simulate.R +#' @importFrom stats runif +#' @param object fitted model object +#' @export setMethod( 'simulate', - 'browntree', + signature=signature(object='browntree'), function (object, nsim = 1, seed = NULL, ...) { if (!is.null(seed)) { if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1) @@ -110,33 +144,117 @@ setMethod( } ) +#' @rdname update +#' @importFrom stats update +#' @export setMethod( 'update', - 'browntree', - function (object, data) { + signature=signature(object='browntree'), + function (object, data, ...) { if (missing(data)) data <- object@data brown( data=data, - tree=object + tree=object, + ... ) } ) -bootstrap <- function (object, nboot = 200, seed = NULL, ...) { - stop("function ",sQuote("bootstrap")," is undefined for objects of class ",sQuote(class(object))) -} -setGeneric('bootstrap') +#' @rdname coef +#' @importFrom stats coef +#' @return \code{coef} applied to a \code{browntree} object extracts a list with three elements: +#' \describe{ +#' \item{\code{sigma}}{the coefficients of the sigma matrix.} +#' \item{\code{theta}}{a list of the estimated optima, one per character.} +#' \item{\code{sigma..sq.matrix}}{the sigma-squared matrix itself.} +#' } +#' @export +setMethod( + 'coef', + signature=signature(object="browntree"), + function (object, ...) { + list( + sigma=object@sigma, + theta=object@theta, + sigma.sq.matrix=sym.par(object@sigma) + ) + } +) +#' @rdname print +#' @export setMethod( - "bootstrap", - "browntree", - function (object, nboot = 200, seed = NULL, ...) { - simdata <- simulate(object,nsim=nboot,seed=seed) - results <- vector(mode='list',length=nboot) - toshow <- c("sigma.squared","theta","loglik","aic","aic.c","sic","dof") - for (b in seq_len(nboot)) { - results[[b]] <- summary(update(object,data=simdata[[b]],...)) - } - as.data.frame(t(sapply(results,function(x)unlist(x[toshow])))) + 'show', + signature=signature(object="browntree"), + function (object) { + print(as(object,'browntree')) + invisible(NULL) + } +) + +#' @rdname print +#' @return \code{print} displays the tree as a table, along with the coefficients of the fitted model and diagnostic information. +#' @export +setMethod( + 'print', + signature=signature(x='browntree'), + function (x, ...) { + cat("\ncall:\n") + print(x@call) + print(as(x,'data.frame'),...) + sm <- summary(x) + cat('\nsigma squared:\n') + print(sm$sigma.squared) + cat('\ntheta:\n') + print(sm$optima) + print(unlist(sm[c("loglik","deviance","aic","aic.c","sic","dof")])) + invisible(x) + } +) + +#' @rdname loglik +#' @importFrom stats logLik +#' @export +setMethod( + "logLik", + signature=signature(object="browntree"), + function(object)object@loglik +) + +#' @rdname summary +#' @return \code{summary} applied to a \code{browntree} object returns information about the fitted model, including parameter estimates and quantities describing the goodness of fit. +#' @export +setMethod( + "summary", + signature=signature(object="browntree"), + function (object, ...) { + cf <- coef(object) + dof <- object@nchar+length(object@sigma) + deviance=-2*logLik(object) + aic <- deviance+2*dof + aic.c <- aic+2*dof*(dof+1)/(object@nterm*object@nchar-dof-1) + sic <- deviance+log(object@nterm*object@nchar)*dof + list( + call=object@call, + sigma.squared=cf$sigma.sq.matrix, + theta=cf$theta, + loglik=logLik(object), + deviance=deviance, + aic=aic, + aic.c=aic.c, + sic=sic, + dof=dof + ) + } +) + +setAs( + from='browntree', + to='data.frame', + def = function (from) { + cbind( + as(as(from,'ouchtree'),'data.frame'), + as.data.frame(from@data) + ) } ) diff --git a/R/coef.R b/R/coef.R new file mode 100644 index 0000000..f3ac78f --- /dev/null +++ b/R/coef.R @@ -0,0 +1,10 @@ +#' Model coefficients +#' +#' \code{coef} extracts the parameters from a fitted model object. +#' +#' @name coef +#' @rdname coef +#' @include brown.R hansen.R +#' @param object fitted model object. +#' @param ... additional arguments, ignored. +NULL diff --git a/R/glssoln.R b/R/glssoln.R index 9c52a07..a8750b6 100644 --- a/R/glssoln.R +++ b/R/glssoln.R @@ -1,14 +1,14 @@ glssoln <- function (a, x, v, tol = sqrt(.Machine$double.eps)) { n <- length(x) vh <- try( - chol(v), - silent=FALSE - ) + chol(v), + silent=FALSE + ) if (inherits(vh,'try-error')) { warning( - "glssoln: Choleski decomposition of variance-covariance matrix fails", - call.=FALSE - ) + "glssoln: Choleski decomposition of variance-covariance matrix fails", + call.=FALSE + ) y <- rep(NA,ncol(a)) e <- rep(NA,n) dim(y) <- ncol(a) @@ -20,15 +20,15 @@ glssoln <- function (a, x, v, tol = sqrt(.Machine$double.eps)) { r <- length(svals) svals <- diag(1/svals,nrow=r,ncol=r) y <- (s$v[,inds,drop=FALSE]%*% - (svals %*% - t(s$u[,inds,drop=FALSE])))%*% - forwardsolve(vh,x,upper.tri=TRUE,transpose=TRUE) + (svals %*% + t(s$u[,inds,drop=FALSE])))%*% + forwardsolve(vh,x,upper.tri=TRUE,transpose=TRUE) e <- a%*%y-x dim(y) <- dim(y)[1] dim(e) <- n } list( - coeff=y, - residuals=e - ) + coeff=y, + residuals=e + ) } diff --git a/R/hansen.R b/R/hansen.R index fb7dd9f..501b199 100644 --- a/R/hansen.R +++ b/R/hansen.R @@ -1,3 +1,46 @@ +#' Ornstein-Uhlenbeck models of trait evolution +#' +#' The function \code{hansen} fits an Ornstein-Uhlenbeck model to data. +#' The fitting is done using \code{optim} or \code{subplex}. +#' +#' The Hansen model for the evolution of a multivariate trait \eqn{X} along a lineage can be written as a stochastic differential equation (Ito diffusion) +#' \deqn{dX=\alpha(\theta(t)-X(t))dt+\sigma dB(t),}{dX = alpha (theta(t)-X(t)) dt + sigma dB(t),} +#' where \eqn{t} is time along the lineage, +#' \eqn{\theta(t)}{theta(t)} is the optimum trait value, \eqn{B(t)} is a standard Wiener process (Brownian motion), +#' and \eqn{\alpha}{alpha} and \eqn{\sigma}{sigma} are matrices +#' quantifying, respectively, the strength of selection and random drift. +#' Without loss of generality, one can assume \eqn{\sigma}{sigma} is lower-triangular. +#' This is because only the infinitesimal variance-covariance matrix +#' \eqn{\sigma^2=\sigma\sigma^T}{sigma^2 = sigma\%*\%transpose(sigma)} +#' is identifiable, and for any admissible variance-covariance matrix, we can choose \eqn{\sigma}{sigma} to be lower-triangular. +#' Moreover, if we view the basic model as describing evolution on a fitness landscape, then \eqn{\alpha}{alpha} will be symmetric. +#' If we further restrict ourselves to the case of stabilizing selection, \eqn{\alpha}{alpha} will be positive definite as well. +#' We make these assumptions and therefore can assume that the matrix \eqn{\alpha}{alpha} has a lower-triangular square root. +#' +#' The \code{hansen} code uses unconstrained numerical optimization to maximize the likelihood. +#' To do this, it parameterizes the \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices in a special way: +#' each matrix is parameterized by \code{nchar*(nchar+1)/2} parameters, where \code{nchar} is the number of quantitative characters. +#' Specifically, the parameters initialized by the \code{sqrt.alpha} argument of \code{hansen} are used +#' to fill the nonzero entries of a lower-triangular matrix (in column-major order), +#' which is then multiplied by its transpose to give the selection-strength matrix. +#' The parameters specified in \code{sigma} fill the nonzero entries in the lower triangular \eqn{\sigma}{sigma} matrix. +#' When \code{hansen} is executed, the numerical optimizer maximizes the likelihood over these parameters. +#' +#' @name hansen +#' @rdname hansen +#' @author Aaron A. King +#' @seealso +#' \code{\link{ouchtree}}, \code{\link{brown}}, \code{\link[stats:optim]{optim}}, \code{\link[subplex:subplex]{subplex}}, \code{\link{bimac}}, \code{\link{anolis.ssd}} +#' @references +#' \Butler2004 +#' +#' \Cressler2015 +#' +#' @keywords models +#' @example examples/hansen.R +#' +NULL + setClass( 'hansentree', contains='ouchtree', @@ -16,6 +59,43 @@ setClass( ) ) +NULL + +#' @rdname hansen +#' @include ouchtree.R glssoln.R rmvnorm.R +#' @importFrom stats optim +#' @importFrom subplex subplex +#' +#' @param data Phenotypic data for extant species, i.e., species at the terminal twigs of the phylogenetic tree. +#' This can either be a single named numeric vector, a list of \code{nchar} named vectors, or a data frame containing \code{nchar} data variables. +#' There must be an entry per variable for every node in the tree; use \code{NA} to represent missing data. +#' If the +#' data are supplied as one or more named vectors, the names attributes are taken to correspond to the node names specified when the \code{ouchtree} was constructed (see \code{\link{ouchtree}}). +#' If the data are supplied as a +#' data-frame, the rownames serve that purpose. +#' @param tree A phylogenetic tree, specified as an \code{ouchtree} object. +#' @param regimes A vector of codes, one for each node in the tree, specifying the selective regimes hypothesized to have been operative. +#' Corresponding to each node, enter the code of the regime hypothesized for the branch segment terminating in that node. +#' For the root node, because it has no branch segment terminating on it, the regime specification is irrelevant. +#' If there are \code{nchar} quantitative characters, then one can specify a single set of \code{regimes} for all characters or a list of \code{nchar} regime specifications, one for each character. +#' @param sqrt.alpha,sigma These are used to initialize the optimization algorithm. +#' The selection strength matrix \eqn{\alpha}{alpha} and the random drift variance-covariance matrix \eqn{\sigma^2}{sigma^2} are parameterized by their matrix square roots. +#' Specifically, these initial guesses are each packed into lower-triangular matrices (column by column). +#' The product of this matrix with its transpose is the \eqn{\alpha}{alpha} or \eqn{\sigma^2}{sigma^2} matrix. +#' See Details. +#' @param fit If \code{fit=TRUE}, then the likelihood will be maximized. +#' If \code{fit=FALSE}, the likelihood will be evaluated at the specified values of \code{sqrt.alpha} and \code{sigma}; +#' the optima \code{theta} will be returned as well. +#' @param method The method to be used by the optimization algorithm. +#' See \code{\link[subplex]{subplex}} and \code{\link{optim}} for information on the available options. +#' @param hessian If \code{hessian=TRUE}, then the Hessian matrix will be computed by \code{optim}. +#' @param \dots Additional arguments will be passed as \code{control} options to \code{optim} or \code{subplex}. +#' See \code{\link{optim}} and \code{\link[subplex:subplex]{subplex}} for information on the available options. +#' +#' @return \code{hansen} returns an object of class \code{hansentree}. +#' For details on the methods of that class, see \code{\link{hansen}}. +#' +#' @export hansen hansen <- function (data, tree, regimes, sqrt.alpha, sigma, fit = TRUE, method = c("Nelder-Mead","subplex","BFGS","L-BFGS-B"), @@ -172,13 +252,13 @@ hansen <- function (data, tree, regimes, sqrt.alpha, sigma, if (hessian) { hs <- opt$hessian - # se <- sqrt(diag(solve(0.5*hs))) - # se.alpha <- se[seq(nalpha)] - # se.sigma <- se[nalpha+seq(nsigma)] + ## se <- sqrt(diag(solve(0.5*hs))) + ## se.alpha <- se[seq(nalpha)] + ## se.sigma <- se[nalpha+seq(nsigma)] } else { hs <- matrix(NA,0,0) - # se.alpha <- rep(NA,nalpha) - # se.sigma <- rep(NA,nalpha) + ## se.alpha <- rep(NA,nalpha) + ## se.sigma <- rep(NA,nalpha) } sol <- ou.lik.fn( @@ -331,9 +411,13 @@ hansen.deviate <- function (n = 1, object) { apply(X,3,as.data.frame) } +#' @rdname simulate +#' @include simulate.R +#' @importFrom stats runif +#' @export setMethod( 'simulate', - 'hansentree', + signature=signature(object='hansentree'), function (object, nsim = 1, seed = NULL, ...) { if (!is.null(seed)) { if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1) @@ -348,9 +432,13 @@ setMethod( } ) +#' @rdname update +#' @importFrom stats update +#' @inheritParams hansen +#' @export setMethod( 'update', - 'hansentree', + signature=signature(object='hansentree'), function (object, data, regimes, sqrt.alpha, sigma, ...) { if (missing(sqrt.alpha)) sqrt.alpha <- object@sqrt.alpha if (missing(sigma)) sigma <- object@sigma @@ -365,17 +453,114 @@ setMethod( } ) +#' @rdname coef +#' @importFrom stats coef +#' @return \code{coef} applied to a \code{hansentree} object returns a named list containing the estimated \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices(given as the \code{alpha.matrix} and \code{sigma.sq.matrix} elements, respectively) but also the MLE returned by the optimizer +#' (as \code{sqrt.alpha} and \code{sigma}, respectively). +#' \strong{The latter elements should not be interpreted, but can be used to restart the algorithm, etc.} +#' @export +setMethod( + 'coef', + signature=signature(object='hansentree'), + function (object, ...) { + list( + sqrt.alpha=object@sqrt.alpha, + sigma=object@sigma, + theta=object@theta, + alpha.matrix=sym.par(object@sqrt.alpha), + sigma.sq.matrix=sym.par(object@sigma) + ) + } +) + +#' @rdname loglik +#' @importFrom stats logLik +#' @export +setMethod( + "logLik", + signature=signature(object='hansentree'), + function (object) object@loglik +) + +#' @rdname summary +#' @return \code{summary} applied to a \code{hansentree} method displays the estimated \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices as well as various quantities describing the goodness of model fit. +#' @export +setMethod( + "summary", + signature=signature(object='hansentree'), + function (object, ...) { + cf <- coef(object) + ## if (length(object@hessian)>0) + ## se <- sqrt(diag(solve(0.5*object@hessian))) + dof <- length(object@sqrt.alpha)+length(object@sigma)+sum(sapply(object@theta,length)) + deviance=-2*logLik(object) + aic <- deviance+2*dof + aic.c <- aic+2*dof*(dof+1)/(object@nterm*object@nchar-dof-1) + sic <- deviance+log(object@nterm*object@nchar)*dof + list( + call=object@call, + conv.code=object@optim.diagn$convergence, + optimizer.message=object@optim.diagn$message, + alpha=cf$alpha.matrix, + sigma.squared=cf$sigma.sq.matrix, + optima=cf$theta, + loglik=logLik(object), + deviance=deviance, + aic=aic, + aic.c=aic.c, + sic=sic, + dof=dof + ) + } +) +#' @rdname print +#' @return \code{print} displays the tree as a table, along with the coefficients of the fitted model and diagnostic information. +#' @export setMethod( - "bootstrap", - "hansentree", - function (object, nboot = 200, seed = NULL, ...) { - simdata <- simulate(object,nsim=nboot,seed=seed) - results <- vector(mode='list',length=nboot) - toshow <- c("alpha","sigma.squared","optima","loglik","aic","aic.c","sic","dof") - for (b in seq_len(nboot)) { - results[[b]] <- summary(update(object,data=simdata[[b]],...)) + 'print', + signature=signature(x='hansentree'), + function (x, ...) { + cat("\ncall:\n") + print(x@call) + print(as(x,'data.frame'),...) + if (length(x@optim.diagn)>0) { + if (x@optim.diagn$convergence!=0) + cat("\n",sQuote("optim")," convergence code: ",x@optim.diagn$convergence) + if (!is.null(x@optim.diagn$message)) + cat("\n",sQuote("optim")," diagnostic message: ",x@optim.diagn$message) } - as.data.frame(t(sapply(results,function(x)unlist(x[toshow])))) + sm <- summary(x) + cat('\nalpha:\n') + print(sm$alpha) + cat('\nsigma squared:\n') + print(sm$sigma.squared) + cat('\ntheta:\n') + print(sm$optima) + print(unlist(sm[c("loglik","deviance","aic","aic.c","sic","dof")])) + invisible(x) + } +) + +#' @rdname print +#' @export +setMethod( + 'show', + signature=signature(object='hansentree'), + function (object) { + print(as(object,'hansentree')) + invisible(NULL) + } +) + +setAs( + from='hansentree', + to='data.frame', + def = function (from) { + cbind( + as(as(from,'ouchtree'),'data.frame'), + as.data.frame(from@regimes), + as.data.frame(from@data) + ) } ) diff --git a/R/loglik.R b/R/loglik.R new file mode 100644 index 0000000..56a63bb --- /dev/null +++ b/R/loglik.R @@ -0,0 +1,11 @@ +#' Log likelihood of a fitted model +#' +#' \code{logLik} extracts the log likelihood from a fitted model object. +#' +#' @name loglik +#' @rdname loglik +#' @include brown.R hansen.R +#' @param object fitted model object +#' @return \code{logLik} returns a numeric value. +#' +NULL diff --git a/R/methods.R b/R/methods.R deleted file mode 100644 index f8742a4..0000000 --- a/R/methods.R +++ /dev/null @@ -1,201 +0,0 @@ -setMethod( - 'coef', - 'browntree', - function (object, ...) { - list( - sigma=object@sigma, - theta=object@theta, - sigma.sq.matrix=sym.par(object@sigma) - ) - } -) - -setMethod( - 'coef', - 'hansentree', - function (object, ...) { - list( - sqrt.alpha=object@sqrt.alpha, - sigma=object@sigma, - theta=object@theta, - alpha.matrix=sym.par(object@sqrt.alpha), - sigma.sq.matrix=sym.par(object@sigma) - ) - } -) - - # coerce methods allow for coercion to a data.frame -setAs( - from='ouchtree', - to='data.frame', - def = function (from) { - df <- data.frame( - nodes=from@nodes, - ancestors=from@ancestors, - times=from@times, - labels=from@nodelabels, - row.names=from@nodes - ) - rownames(df) <- from@nodes - df - } -) - -setAs( - from='browntree', - to='data.frame', - def = function (from) { - cbind( - as(as(from,'ouchtree'),'data.frame'), - as.data.frame(from@data) - ) - } -) - -setAs( - from='hansentree', - to='data.frame', - def = function (from) { - cbind( - as(as(from,'ouchtree'),'data.frame'), - as.data.frame(from@regimes), - as.data.frame(from@data) - ) - } -) - - -setMethod('logLik','browntree',function(object)object@loglik) - -setMethod( - "summary", - "browntree", - function (object, ...) { - cf <- coef(object) - dof <- object@nchar+length(object@sigma) - deviance=-2*logLik(object) - aic <- deviance+2*dof - aic.c <- aic+2*dof*(dof+1)/(object@nterm*object@nchar-dof-1) - sic <- deviance+log(object@nterm*object@nchar)*dof - list( - call=object@call, - sigma.squared=cf$sigma.sq.matrix, - theta=cf$theta, - loglik=logLik(object), - deviance=deviance, - aic=aic, - aic.c=aic.c, - sic=sic, - dof=dof - ) - } -) - -setMethod('logLik','hansentree',function(object)object@loglik) - -setMethod( - "summary", - "hansentree", - function (object, ...) { - cf <- coef(object) - ## if (length(object@hessian)>0) - ## se <- sqrt(diag(solve(0.5*object@hessian))) - dof <- length(object@sqrt.alpha)+length(object@sigma)+sum(sapply(object@theta,length)) - deviance=-2*logLik(object) - aic <- deviance+2*dof - aic.c <- aic+2*dof*(dof+1)/(object@nterm*object@nchar-dof-1) - sic <- deviance+log(object@nterm*object@nchar)*dof - list( - call=object@call, - conv.code=object@optim.diagn$convergence, - optimizer.message=object@optim.diagn$message, - alpha=cf$alpha.matrix, - sigma.squared=cf$sigma.sq.matrix, - optima=cf$theta, - loglik=logLik(object), - deviance=deviance, - aic=aic, - aic.c=aic.c, - sic=sic, - dof=dof - ) - } -) - -setMethod( - 'print', - 'ouchtree', - function (x, ...) { - print(as(x,'data.frame'),...) - invisible(x) - } -) - -setMethod( - 'print', - 'browntree', - function (x, ...) { - cat("\ncall:\n") - print(x@call) - print(as(x,'data.frame'),...) - sm <- summary(x) - cat('\nsigma squared:\n') - print(sm$sigma.squared) - cat('\ntheta:\n') - print(sm$optima) - print(unlist(sm[c("loglik","deviance","aic","aic.c","sic","dof")])) - invisible(x) - } -) - -setMethod( - 'print', - 'hansentree', - function (x, ...) { - cat("\ncall:\n") - print(x@call) - print(as(x,'data.frame'),...) - if (length(x@optim.diagn)>0) { - if (x@optim.diagn$convergence!=0) - cat("\n",sQuote("optim")," convergence code: ",x@optim.diagn$convergence) - if (!is.null(x@optim.diagn$message)) - cat("\n",sQuote("optim")," diagnostic message: ",x@optim.diagn$message) - } - sm <- summary(x) - cat('\nalpha:\n') - print(sm$alpha) - cat('\nsigma squared:\n') - print(sm$sigma.squared) - cat('\ntheta:\n') - print(sm$optima) - print(unlist(sm[c("loglik","deviance","aic","aic.c","sic","dof")])) - invisible(x) - } -) - -setMethod( - 'show', - 'ouchtree', - function (object) { - print(as(object,'ouchtree')) - invisible(NULL) - } -) - -setMethod( - 'show', - 'browntree', - function (object) { - print(as(object,'browntree')) - invisible(NULL) - } -) - -setMethod( - 'show', - 'hansentree', - function (object) { - print(as(object,'hansentree')) - invisible(NULL) - } -) diff --git a/R/ouchtree.R b/R/ouchtree.R index 93e745c..da19f4f 100644 --- a/R/ouchtree.R +++ b/R/ouchtree.R @@ -1,3 +1,22 @@ +#' Phylogenetic tree object in 'ouch' format. +#' +#' An object containing a phylogenetic tree in a form suitable for using \pkg{ouch} methods. +#' +#' \code{ouchtree} creates an \code{ouchtree} object. +#' This contains the topology, branch times, and epochs. +#' It also (optionally) holds names of taxa for display purposes. +#' +#' @name ouchtree +#' @rdname ouchtree +#' @aliases ouchtree ouchtree-class +#' +#' @author Aaron A. King +#' @seealso \code{ouchtree}, \code{ape2ouch}, \code{brown}, \code{hansen} +#' @keywords models +#' @example examples/bimac1.R +#' +NULL + setClass( 'ouchtree', representation=representation( @@ -17,6 +36,24 @@ setClass( ) ) +#' @rdname ouchtree +#' @param nodes A character vector giving the name of each node. +#' These are used internally and must be unique. +#' @param ancestors Specification of the topology of the phylogenetic tree. +#' This is in the form of a character vector specifying the name +#' (as given in the \code{nodes} argument) +#' of the immediate ancestor of each node. +#' In particular, the i-th name is that of the ancestor of the i-th node. +#' The root node is distinguished by having no ancestor (i.e., \code{NA}). +#' @param times A vector of nonnegative numbers, one per node in the tree, +#' specifying the time at which each node is located. +#' Time should be increasing from the root node to the terminal twigs. +#' @param labels Optional vector of node labels. +#' These will be used in plots to label nodes. +#' It is not necessary that these be unique. +#' +#' @include package.R +#' @export ouchtree ouchtree <- function (nodes, ancestors, times, labels = as.character(nodes)) { nodes <- as.character(nodes) @@ -160,3 +197,40 @@ is.root.node <- function (anc) { is.na(anc) } +#' @rdname print +#' @export +setMethod( + 'print', + signature=signature(x='ouchtree'), + function (x, ...) { + print(as(x,'data.frame'),...) + invisible(x) + } +) + +#' @rdname print +#' @export +setMethod( + 'show', + signature=signature(object='ouchtree'), + function (object) { + print(as(object,'ouchtree')) + invisible(NULL) + } +) + +setAs( + from='ouchtree', + to='data.frame', + def = function (from) { + df <- data.frame( + nodes=from@nodes, + ancestors=from@ancestors, + times=from@times, + labels=from@nodelabels, + row.names=from@nodes + ) + rownames(df) <- from@nodes + df + } +) diff --git a/R/package.R b/R/package.R new file mode 100644 index 0000000..3742cba --- /dev/null +++ b/R/package.R @@ -0,0 +1,48 @@ +#' Ornstein-Uhlenbeck methods for comparative phylogenetic hypotheses +#' +#' The \pkg{ouch} package provides facilities for phylogenetic comparative analysis based on Ornstein-Uhlenbeck models of trait evolution along a phylogeny. +#' Multivariate data and complex adaptive hypotheses are supported. +#' +#' @name ouch-package +#' @rdname package +#' @docType package +#' @section Classes: +#' The basic class, \code{ouchtree}, is provided to encode a +#' phylogenetic tree. Plot and print methods are provided. +#' +#' The class \code{browntree} derives from class \code{ouchtree} and encodes +#' the results of fitting a Brownian Motion model to data. +#' +#' The class \code{hansentree} derives from class \code{ouchtree} and encodes +#' the results of fitting a Hansen model to data. +#' @section Detailed Documentation: +#' \describe{ +#' \item{Phylogenies in \pkg{ouch} format}{\code{\link{ouchtree}}, \code{\link{ape2ouch}}} +#' \item{Brownian motion models}{\code{\link{brown}}} +#' \item{Ornstein-Uhlenbeck models}{\code{\link{hansen}}} +#' \item{Simulation of models}{\code{\link[ouch:simulate]{simulate}}} +#' \item{Display of data}{\code{\link[ouch:plot]{plot}}} +#' } +#' @author Aaron A. King +#' @references +#' \Butler2004 +#' +#' \Cressler2015 +#' +#' @useDynLib ouch, .registration = TRUE +#' @import methods +#' @keywords models +NULL + +undef_method <- function (method, object) { + o <- deparse(substitute(object)) + stop(sQuote(method)," is undefined for ",sQuote(o)," of class ", + sQuote(class(object)),".",call.=FALSE) +} + +reqd_arg <- function (method, object) { + if (is.null(method) || length(method)==0) + stop(sQuote(object)," is a required argument.",call.=FALSE) + else + stop(method,sQuote(object)," is a required argument.",call.=FALSE) +} diff --git a/R/paint.R b/R/paint.R index d696a5a..df09387 100644 --- a/R/paint.R +++ b/R/paint.R @@ -1,3 +1,35 @@ +#' Painting regimes on a phylogenetic tree +#' +#' Function to paint selective regimes on a phylogenetic tree. +#' +#' The names of \code{subtree} and \code{branch} must be the names of nodes of \code{tree}. +#' The painting proceeds in a particular order: +#' one can overpaint a branch. +#' The subtrees indicated by the elements of \code{subtree} are painted first, in order. +#' Then the branches indicated by \code{branch} are painted. +#' If \code{tree} is a simple \code{ouchtree} object, then \code{paint} begins with a blank canvas, +#' i.e., a tree painted with the single regime \dQuote{nonspec}. +#' If \code{tree} is of class \code{hansentree}, then \code{paint} begins with the regimes specified in the \code{regimes} slot of \code{tree}. +#' Note that, if \code{tree} is a multivariate \code{hansentree}, then there are multiple regime specifications contained in \code{tree}. +#' In this case, the argument \code{which} lets you pick which one you wish to begin with; +#' by default, the first is used. +#' +#' @param tree An object of class \code{ouchtree}. +#' @param subtree An optional named vector specifying the root nodes of subtrees. +#' Each branch that descends from this node will be painted with the +#' specified regime. +#' @param branch An optional named vector specifying the end nodes of branches. +#' The unique branch that terminates at the named node will be painted with the specified regime. +#' @param which integer; +#' if \code{tree} is a \code{hansentree}, start not with a blank canvas but with the regime specifications \code{tree} contains for the character indicated by \code{which}. +#' @return A vector of class \sQuote{factor} with names corresponding to the nodes in \code{tree}, specifying selective regimes. +#' @author Aaron A. King +#' @seealso \code{ouchtree}, \code{hansen} +#' @keywords models +#' @example examples/paint.R +#' @include ouchtree.R +#' @importFrom utils head tail +#' @export paint paint <- function (tree, subtree, branch, which = 1) { if (!is(tree,'ouchtree')) stop(sQuote("tree")," must be of class ",sQuote("ouchtree")) @@ -48,4 +80,3 @@ paint <- function (tree, subtree, branch, which = 1) { } as.factor(regimes) } - diff --git a/R/plot.R b/R/plot.R index 5968c56..095ce43 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,3 +1,76 @@ +#' Plotting functions. +#' +#' Plot phylogenies. +#' +#' @name plot +#' @rdname plot +#' @param x object to plot. +#' @param y ignored. +#' @param regimes factor or character; a vector of regime paintings. +#' @param node.names logical; should node names be displayed? +#' @param legend logical; display a legend? +#' @param labels character; taxon labels. +#' @param ... additional arguments, passed to \code{text}. +#' +NULL + +#' @rdname plot +#' @importFrom grDevices rainbow +#' @importFrom graphics par text +#' @importFrom stats setNames +#' @export +setMethod( + "plot", + signature=signature(x="ouchtree"), + function (x, y, regimes = NULL, node.names = FALSE, legend = TRUE, ..., labels) { + if (!missing(y)) warning(sQuote("y")," is ignored.") + labels <- x@nodelabels + if (node.names) { + lbld <- !is.na(labels) + labels[lbld] <- paste(x@nodes[lbld],labels[lbld]) + labels[!lbld] <- x@nodes[!lbld] + } + if (is.data.frame(regimes)) { + nm <- rownames(regimes) + regimes <- lapply(as.list(regimes),function(x){names(x)<-nm;x}) + } + if (is.list(regimes)) { + if (any(sapply(regimes,length)!=x@nnodes)) + stop("each element in ",sQuote("regimes")," must be a vector with one entry per node of the tree") + } else if (!is.null(regimes)) { + if (length(regimes)!=x@nnodes) + stop("there must be one entry in ",sQuote("regimes")," per node of the tree") + nm <- deparse(substitute(regimes))[1] + regimes <- list(regimes) + names(regimes) <- nm + } + if (is.null(regimes)) { + invisible(tree.plot.internal(x,regimes=NULL,labels=labels,legend=legend,...)) + } else { + oldpar <- par(mfrow=c(1,length(regimes))) + on.exit(par(oldpar)) + retval <- lapply( + regimes, + function (r) tree.plot.internal(x,regimes=r,labels=labels,...) + ) + invisible(retval) + } + } +) + +#' @rdname plot +#' @export +setMethod( + "plot", + signature=signature(x="hansentree"), + function (x, y, regimes, ...) { + if (!missing(y)) warning(sQuote("y")," is ignored.") + if (missing(regimes)) regimes <- x@regimes + f <- getMethod("plot","ouchtree") + f(x,regimes=regimes,...) + } +) + tree.plot.internal <- function (x, regimes = NULL, labels = x@nodelabels, legend = TRUE, ...) { rx <- range(x@times,na.rm=T) rxd <- 0.1*diff(rx) @@ -61,51 +134,3 @@ arrange.tree <- function (root, anc, term) { y[oo] } } - -setMethod( - 'plot', - 'ouchtree', - function (x, y, regimes = NULL, node.names = FALSE, legend = TRUE, ..., labels) { - labels <- x@nodelabels - if (node.names) { - lbld <- !is.na(labels) - labels[lbld] <- paste(x@nodes[lbld],labels[lbld]) - labels[!lbld] <- x@nodes[!lbld] - } - if (is.data.frame(regimes)) { - nm <- rownames(regimes) - regimes <- lapply(as.list(regimes),function(x){names(x)<-nm;x}) - } - if (is.list(regimes)) { - if (any(sapply(regimes,length)!=x@nnodes)) - stop("each element in ",sQuote("regimes")," must be a vector with one entry per node of the tree") - } else if (!is.null(regimes)) { - if (length(regimes)!=x@nnodes) - stop("there must be one entry in ",sQuote("regimes")," per node of the tree") - nm <- deparse(substitute(regimes))[1] - regimes <- list(regimes) - names(regimes) <- nm - } - if (is.null(regimes)) { - invisible(tree.plot.internal(x,regimes=NULL,labels=labels,legend=legend,...)) - } else { - oldpar <- par(mfrow=c(1,length(regimes))) - on.exit(par(oldpar)) - retval <- lapply( - regimes, - function (r) tree.plot.internal(x,regimes=r,labels=labels,...) - ) - invisible(retval) - } - } -) - -setMethod( - 'plot', - 'hansentree', - function (x, y, regimes, ...) { - if (missing(regimes)) regimes <- x@regimes - f <- getMethod("plot","ouchtree") - f(x=x,y=y,regimes=regimes,...) - } -) diff --git a/R/print.R b/R/print.R new file mode 100644 index 0000000..9c1ce7e --- /dev/null +++ b/R/print.R @@ -0,0 +1,10 @@ +#' Print methods +#' +#' @name print +#' @rdname print +#' @include ouchtree.R brown.R hansen.R +#' @param object,x object to display. +#' @param ... additional arguments, ignored. +#' @return \code{print} displays the tree as a table, with (possibly) other information. +#' +NULL diff --git a/R/rmvnorm.R b/R/rmvnorm.R index 203bcf9..94f9d3d 100644 --- a/R/rmvnorm.R +++ b/R/rmvnorm.R @@ -1,3 +1,5 @@ +#' @importFrom stats rnorm + rmvnorm <- function (n = 1, mean, var) { p <- length(mean) if (!all(dim(var)==c(p,p))) @@ -5,3 +7,4 @@ rmvnorm <- function (n = 1, mean, var) { cf <- t(chol(var)) matrix(mean,p,n)+cf%*%matrix(rnorm(p*n),p,n) } + diff --git a/R/simulate.R b/R/simulate.R new file mode 100644 index 0000000..9835639 --- /dev/null +++ b/R/simulate.R @@ -0,0 +1,20 @@ +#' Simulations of a phylogenetic trait model. +#' +#' \code{simulate} generates random deviates from a fitted model. +#' +#' @name simulate +#' @docType methods +#' @rdname simulate +#' @include brown.R hansen.R +#' @importFrom stats simulate +#' +#' @return +#' \code{simulate} returns a list of data-frames, each comparable to the original data. +#' +#' @param object fitted model object to simulate. +#' @param nsim integer; number of independent simulations. +#' @param seed integer; if non-\code{NULL}, the RNG will be initialized with this seed for the simulations. +#' The RNG will be reset to its pre-existing state when \code{simulate} returns. +#' @param ... additional arguments, ignored. +NULL + diff --git a/R/summary.R b/R/summary.R new file mode 100644 index 0000000..89966be --- /dev/null +++ b/R/summary.R @@ -0,0 +1,9 @@ +#' Summary methods +#' +#' @name summary +#' @rdname summary +#' @include ouchtree.R brown.R hansen.R +#' @param object fitted model object. +#' @param ... additional arguments, ignored. +#' +NULL diff --git a/R/update.R b/R/update.R new file mode 100644 index 0000000..aa5f69f --- /dev/null +++ b/R/update.R @@ -0,0 +1,18 @@ +#' Update and refit a model. +#' +#' \code{update} will update a model and re-fit. +#' This allows one to change the data and/or parameters. +#' +#' @name update +#' @docType methods +#' @rdname update +#' @include brown.R hansen.R +#' @importFrom stats update +#' +#' @return +#' @return \code{update} returns a new fitted-model object of the same class as \code{object}. +#' @param object fitted model object. +#' @param data data that replace those used in the original fit. +#' @param ... Additional arguments replace the corresponding arguments in the original call. +NULL + diff --git a/examples/anolis.R b/examples/anolis.R new file mode 100644 index 0000000..59a6abb --- /dev/null +++ b/examples/anolis.R @@ -0,0 +1,15 @@ +tree <- with(anolis.ssd,ouchtree(node,ancestor,time/max(time),species)) +plot(tree,node.names=TRUE) + +h1 <- brown(anolis.ssd['log.SSD'],tree) +h1 +plot(h1) + +h2 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.1'],sqrt.alpha=1,sigma=1) +h2 +plot(h2) + +h3 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.7'],sqrt.alpha=1,sigma=1) +h3 +plot(h3) + diff --git a/examples/bimac.R b/examples/bimac.R new file mode 100644 index 0000000..d32afd1 --- /dev/null +++ b/examples/bimac.R @@ -0,0 +1,28 @@ +tree <- with(bimac,ouchtree(node,ancestor,time/max(time),species)) +plot(tree,node.names=TRUE) + +h1 <- brown(log(bimac['size']),tree) +h1 +plot(h1) + +h2 <- hansen(log(bimac['size']),tree,bimac['OU.1'],sqrt.alpha=1,sigma=1) +h2 +plot(h2) + +h3 <- hansen(log(bimac['size']),tree,bimac['OU.3'],sqrt.alpha=1,sigma=1) +h3 +plot(h3) + +h4 <- hansen(log(bimac['size']),tree,bimac['OU.4'],sqrt.alpha=1,sigma=1) +h4 +plot(h4) + +h5 <- hansen(log(bimac['size']),tree,bimac['OU.LP'],sqrt.alpha=1,sigma=1,reltol=1e-5) +h5 <- update(h5,method='subplex',reltol=1e-11,parscale=c(0.1,0.1),hessian=TRUE) +h5 + +simdat <- simulate(h5,nsim=10) +hsim <- update(h5,data=simdat[[1]]) +summary(hsim) +bsim <- update(h1,data=simdat[[1]]) +summary(bsim) diff --git a/examples/bimac1.R b/examples/bimac1.R new file mode 100644 index 0000000..cf574f1 --- /dev/null +++ b/examples/bimac1.R @@ -0,0 +1,8 @@ +tree <- with( + bimac, + ouchtree(nodes=node,ancestors=ancestor,times=time,labels=species) +) +tree + +plot(tree) +plot(tree,node.names=TRUE) diff --git a/examples/hansen.R b/examples/hansen.R new file mode 100644 index 0000000..5ce0d19 --- /dev/null +++ b/examples/hansen.R @@ -0,0 +1,48 @@ +\dontrun{ + if (library(geiger,logical.return=TRUE)) { + +### an example data set (Darwin's finches) + data(geospiza) + str(geospiza) + sapply(geospiza,class) + +### check the correspondence between data and tree tips: + print(nc <- with(geospiza,name.check(geospiza.tree,geospiza.data))) +### looks like one of the terminal twigs has no data associated +### drop that tip: + tree <- with(geospiza,drop.tip(geospiza.tree,nc$tree_not_data)) + dat <- as.data.frame(geospiza$dat) + + +### make an ouchtree out of the phy-format tree + ot <- ape2ouch(tree) + +### merge data with tree info + otd <- as(ot,"data.frame") +### in these data, it so happens that the rownames correspond to node names +### we will exploit this correspondence in the 'merge' operation: + dat$labels <- rownames(dat) + otd <- merge(otd,dat,by="labels",all=TRUE) + rownames(otd) <- otd$nodes + print(otd) +### this data-frame now contains the data as well as the tree geometry + +### now remake the ouch tree + ot <- with(otd,ouchtree(nodes=nodes,ancestors=ancestors,times=times,labels=labels)) + + b1 <- brown(tree=ot,data=otd[c("tarsusL","beakD")]) + summary(b1) + +### evaluate an OU model with a single, global selective regime + otd$regimes <- as.factor("global") + h1 <- hansen( + tree=ot, + data=otd[c("tarsusL","beakD")], + regimes=otd["regimes"], + sqrt.alpha=c(1,0,1), + sigma=c(1,0,1), + maxit=10000 + ) + summary(h1) + } +} diff --git a/examples/paint.R b/examples/paint.R new file mode 100644 index 0000000..d716f84 --- /dev/null +++ b/examples/paint.R @@ -0,0 +1,14 @@ +x <- with( + bimac, + ouchtree(nodes=node,times=time/max(time),ancestors=ancestor,labels=species) +) + +r <- paint(x,subtree=c("1"="medium","9"="large","2"="small"), + branch=c("38"="large","2"="medium")) +plot(x,regimes=r,node.names=TRUE) + +## compare to bimac['OU.LP'] +h5 <- hansen(data=log(bimac['size']),tree=x,regimes=bimac['OU.LP'], + sqrt.alpha=1,sigma=1,reltol=1e-5) +r <- paint(h5,branch=c("18"="large"),subtree=c("9"="small")) +plot(x,regimes=r,node.names=TRUE) diff --git a/inst/NEWS b/inst/NEWS index 2f73c13..7c89dd2 100644 --- a/inst/NEWS +++ b/inst/NEWS @@ -1,5 +1,11 @@ _N_e_w_s _f_o_r _P_a_c_k_a_g_e '_o_u_c_h' +_C_h_a_n_g_e_s _i_n _o_u_c_h _v_e_r_s_i_o_n _2._1_5: + + • The ‘plot’ function now spaces the terminal taxa evenly. + + • The documentation has received an overhaul. + _C_h_a_n_g_e_s _i_n _o_u_c_h _v_e_r_s_i_o_n _2._1_4: • Depends on R>=3.6. This change necessitated by modifications diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 306a83b..2983bb4 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -1,5 +1,11 @@ \name{NEWS} \title{News for Package 'ouch'} +\section{Changes in ouch version 2.15}{ + \itemize{ + \item The \code{plot} function now spaces the terminal taxa evenly. + \item The documentation has received an overhaul. + } +} \section{Changes in ouch version 2.14}{ \itemize{ \item Depends on R>=3.6. diff --git a/man/anolis.ssd.Rd b/man/anolis.ssd.Rd deleted file mode 100644 index 612d8e9..0000000 --- a/man/anolis.ssd.Rd +++ /dev/null @@ -1,63 +0,0 @@ -\name{anolis.ssd} -\alias{anolis.ssd} -\docType{data} -\title{Greater Antillean anolis lizard sexual size dimorphism data.} -\description{ - The dataset consists of sexual size-dimorphism data for 38 species of anoles from Cuba, Hispaniola, Jamaica, and Puerto Rico (Butler, Schoener, and Losos 2000). - Each of these species belongs to one of six microhabitat types, or ``ecomorphs'' (sensu Williams, 1972): - trunk-ground, grass-bush, trunk, trunk-crown, twig, and crown-giant. - The data were used to demonstrate an evolutionary association between habitat type and degree of sexual size dimorphism. -} -\usage{data(anolis.ssd)} -\format{ - A data frame with 38 observations on the following 6 variables. - \describe{ - \item{node}{Labels for the nodes.} - \item{species}{Names of extant species.} - \item{log.SSD}{Log sexual size dimorphism of extant species.} - \item{ancestor}{Ancestor node.} - \item{time}{Time of node.} - \item{OU.1}{a factor with levels \code{ns}} - \item{OU.7}{a factor with levels corresponding to ecomorph - (\code{tg} \code{tc} \code{gb} \code{cg} \code{tw} \code{tr} - \code{anc})} - } -} -\details{ - Size dimorphism was calcuated as the log-ratio of male snout-to-vent length to female snout-to-vent length (males are larger). - - In this example, we tested three models of evolution: - Brownian motion, Ornstein-Uhlenbeck with one global optimum, and Ornstein-Uhlenbeck with 7 optima (one for each ecomorph type plus an additional one for an ``unknown'' type). - - For the 7-optima model, we assigned each terminal branch to an optimum according to the ecomorph type of the extant species. - Because we had no information to help guide hypotheses about internal branches, we assigned internal branches to the ``unknown'' selective regime. - The phylogeny of these species is consistent with and adaptive radiation, with a burst of speciation events early in the evolutionary history of this clade (see phylogeny in Butler & King (2004) or example below. -} -\source{ - Butler, M.A. and A.A. King. 2004. - Phylogenetic comparative analysis: a modeling approach for adaptive evolution. - American Naturalist 164:683-695. -} -\references{ - Butler, M. A., T. W. Schoener, and J. B. Losos. 2000. - The relationship between sexual size dimorphism and habitat use in Greater Antillean Anolis lizards. - Evolution, 54:259-272. - - Williams, E. E. 1972. - The origin of faunas. Evolution of lizard congeners in a complex island fauna: a trial analysis. - Evol. Biol., 6:47-89. -} -\examples{ -data(anolis.ssd) -tree <- with(anolis.ssd,ouchtree(node,ancestor,time/max(time),species)) -plot(tree,node.names=TRUE) -print(h1 <- brown(anolis.ssd['log.SSD'],tree)) -plot(h1) -print(h2 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.1'],sqrt.alpha=1,sigma=1)) -plot(h2) -print(h3 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.7'],sqrt.alpha=1,sigma=1)) -plot(h3) -} -\author{Marguerite A. Butler and Aaron A. King } -\keyword{models} - diff --git a/man/anolis_ssd.Rd b/man/anolis_ssd.Rd new file mode 100644 index 0000000..a6945ee --- /dev/null +++ b/man/anolis_ssd.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/anolis.R +\docType{data} +\name{anolis.ssd} +\alias{anolis.ssd} +\title{Greater Antillean anolis lizard sexual size dimorphism data} +\format{ +A data frame with 38 observations on the following 6 variables. +\describe{ + \item{node}{Labels for the nodes.} + \item{species}{Names of extant species.} + \item{log.SSD}{Log sexual size dimorphism of extant species.} + \item{ancestor}{Ancestor node.} + \item{time}{Time of node.} + \item{OU.1}{a factor with levels \code{ns}} + \item{OU.7}{a factor with levels corresponding to ecomorph + (\code{tg} \code{tc} \code{gb} \code{cg} \code{tw} \code{tr} \code{anc})} +} +} +\source{ +\Butler2004 +} +\description{ +The dataset consists of sexual size-dimorphism data for 38 species of anoles from Cuba, Hispaniola, Jamaica, and Puerto Rico (Butler, Schoener, and Losos 2000). +Each of these species belongs to one of six microhabitat types, or \dQuote{ecomorphs} (sensu Williams, 1972): +trunk-ground, grass-bush, trunk, trunk-crown, twig, and crown-giant. +The data were used to demonstrate an evolutionary association between habitat type and degree of sexual size dimorphism. +} +\details{ +Size dimorphism was calcuated as the log-ratio of male snout-to-vent length to female snout-to-vent length (males are larger). + +In this example, we tested three models of evolution: +Brownian motion, Ornstein-Uhlenbeck with one global optimum, and Ornstein-Uhlenbeck with 7 optima (one for each ecomorph type plus an additional one for an \dQuote{unknown} type). + +For the 7-optima model, we assigned each terminal branch to an optimum according to the ecomorph type of the extant species. +Because we had no information to help guide hypotheses about internal branches, we assigned +internal branches to the \dQuote{unknown} selective regime. +The phylogeny of these species is consistent with and adaptive radiation, with a burst of speciation events early in the evolutionary history of this clade (see phylogeny in Butler & King (2004) or example below). +} +\examples{ +tree <- with(anolis.ssd,ouchtree(node,ancestor,time/max(time),species)) +plot(tree,node.names=TRUE) + +h1 <- brown(anolis.ssd['log.SSD'],tree) +h1 +plot(h1) + +h2 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.1'],sqrt.alpha=1,sigma=1) +h2 +plot(h2) + +h3 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.7'],sqrt.alpha=1,sigma=1) +h3 +plot(h3) + +} +\references{ +\Butler2000 + +\Williams1972 +} +\author{ +Marguerite A. Butler, Aaron A. King +} +\keyword{models} diff --git a/man/ape2ouch.Rd b/man/ape2ouch.Rd index 91a1fcf..d6c1ad1 100644 --- a/man/ape2ouch.Rd +++ b/man/ape2ouch.Rd @@ -1,24 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ape2ouch.R \name{ape2ouch} \alias{ape2ouch} \title{Convert an "ape" tree to an "ouch" tree.} -\description{ - \code{ape2ouch} translates \pkg{ape}'s \code{phylo} representation of a phylogenetic tree into \pkg{ouch}'s \code{ouchtree} representation. - The user can change the branch lengths while preserving the topology. -} \usage{ ape2ouch(tree, scale = TRUE, branch.lengths = tree$edge.length) } \arguments{ - \item{tree}{ - a tree of class \code{phylo} created in package \pkg{ape}. - } - \item{scale}{ - if \code{scale=TRUE}, the tree's depth will be scaled to 1. - If \code{scale} is a number, then the branch lengths will be scaled by this number. - } - \item{branch.lengths}{ - optional vector of branch lengths. - } +\item{tree}{a tree of class \code{phylo} created in package \pkg{ape}.} + +\item{scale}{if \code{scale=TRUE}, the tree's depth will be scaled to 1. +If \code{scale} is a number, then the branch lengths will be scaled by this number.} + +\item{branch.lengths}{optional vector of branch lengths.} +} +\description{ +\code{ape2ouch} translates \pkg{ape}'s \code{phylo} representation of a phylogenetic tree into \pkg{ouch}'s \code{ouchtree} representation. +The user can change the branch lengths while preserving the topology. +} +\author{ +A. A. King, D. Ackerly } -\author{Aaron A. King } \keyword{models} diff --git a/man/bimac.Rd b/man/bimac.Rd index 5ba23cd..88f0858 100644 --- a/man/bimac.Rd +++ b/man/bimac.Rd @@ -1,93 +1,94 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bimac.R +\docType{data} \name{bimac} \alias{bimac} -\docType{data} \title{Anolis bimaculatus lizard size data.} -\description{ - This is the Anolis bimaculatus dataset used in Butler & King (2004). - It is used to test a hypothesis of character displacement using an interspecific dataset of body sizes and current data on sympatry/allopatry. - The data frame consists of the following columns: - \code{species} which are species names, - \code{size} which is the phenotypic data, - and the variables \code{ancestor} and \code{time} which specify the topology of the phylogeny and the location of the nodes in time, respectively. - The columns \code{OU.1}, \code{OU.3}, \code{OU.4}, and \code{OU.LP} specify four hypothetical arrangements of selective regimes. - Explanations of the data are given below. -} -\usage{data(bimac)} \format{ - A data frame with 45 observations on the following 8 variables. - \describe{ - \item{node}{Labels for the nodes.} - \item{species}{Species names for extant species.} - \item{size}{Body size (head length in mm) of extant species.} - \item{ancestor}{Ancestral node.} - \item{time}{Time of node.} - \item{OU.1}{a factor with levels \code{ns}} - \item{OU.3}{a factor with levels \code{small} \code{medium} \code{large}} - \item{OU.4}{a factor with levels \code{small} \code{medium} \code{large} \code{anc}} - \item{OU.LP}{a factor with levels \code{small} \code{medium} \code{large}} - } +A data frame with 45 observations on the following 8 variables. +\describe{ + \item{node}{Labels for the nodes.} + \item{species}{Species names for extant species.} + \item{size}{Body size (head length in mm) of extant species.} + \item{ancestor}{Ancestral node.} + \item{time}{Time of node.} + \item{OU.1}{a factor with levels \code{ns}} + \item{OU.3}{a factor with levels \code{small}, \code{medium}, \code{large}} + \item{OU.4}{a factor with levels \code{small}, \code{medium}, \code{large}, \code{anc}} + \item{OU.LP}{a factor with levels \code{small}, \code{medium}, \code{large}} } -\details{ - \describe{ - \item{Body size.}{ - We use the phenotypic data and phylogeny of Losos (1990), which employed the head lengths (of males) as a proxy for body size. - In this group of lizards, head length correlates very strongly with snout-to-vent length and the cube root of mass, which are standard measures of body size. - The data are head lengths in mm, note that we use the log of this value in analyses. - } - \item{Tree topology}{ - The tree topology is encoded via two vectors: \code{ancestor} and \code{time}. - Each node of the phylogenetic tree has a corresponding row in the data frame, numbered from 1 to 45. - The columns \code{ancestor} and \code{time} specify the phylogeny. - The \code{ancestor} variable specifies the topology: it is a list indicating the ancestor of each node. - The root node has ancestor 0. - The variable \code{time} specifies the temporal location of each node, with the root node being at time 0. - } - \item{Specifications of selective regimes.}{ - (Columns \code{OU.1}, \code{OU.3}, \code{OU.4}, \code{OU.LP}). - These columns are factors, the levels of which correspond to the ``paintings'' of the respective adaptive regime hypotheses onto the phylogeny. - Each selective regime is named (small, medium, large, etc.). - Put the corresponding name on each branch segment to indicate which selective regime it belongs to. - Each column corresponds to a different painting of the selective regimes, and thus to a different hypothesis. - In this example, there are 3 alternative models (see Butler & King 2004): - \code{OU.4} is 4-regime model, - \code{OU.3} is 3-regime model (all ancestors are medium), - \code{OU.LP} is linear parsimony model. - } - } } \source{ - Butler, M.A. and A.A. King. 2004. - Phylogenetic comparative analysis: a modeling approach for adaptive evolution. - American Naturalist 164:683-695. +\Butler2004 +} +\description{ +This is the \emph{Anolis bimaculatus} dataset used in Butler & King (2004). +It is used to test a hypothesis of character displacement using an interspecific dataset of body sizes and current data on sympatry/allopatry. +The data frame has the following columns: +\code{species} which are species names, +\code{size} which is the phenotypic data, +and the variables \code{ancestor} and \code{time} which specify the topology of the phylogeny and the location of the nodes in time, respectively. +The columns \code{OU.1}, \code{OU.3}, \code{OU.4}, and \code{OU.LP} specify four hypothetical arrangements of selective regimes. +Explanations of the data follow: +\describe{ + \item{Body size.}{ + We use the phenotypic data and phylogeny of Losos (1990), which employed the head lengths (of males) as a proxy for body size. + In this group of lizards, head length correlates very strongly with snout-to-vent length and the cube root of mass, which are standard measures of body size. + The data are head lengths in mm; note that we use the log of this value in analyses. + } + \item{Tree topology}{ + The tree topology is encoded via two vectors: \code{ancestor} and \code{time}. + Each node of the' phylogenetic tree has a corresponding row in the data frame, numbered from 1 to 45. + The columns \code{ancestor} and \code{time} specify the phylogeny. + The \code{ancestor} variable specifies the topology: it is a list indicating the ancestor of each node. + The root node has ancestor 0. + The variable \code{time} specifies the temporal location of each node, with the root node being at time 0. + } + \item{Specifications of selective regimes.}{ + (Columns \code{OU.1}, \code{OU.3}, \code{OU.4}, \code{OU.LP}). + These columns are factors, the levels of which correspond to the \dQuote{paintings} of the respective adaptive regime hypotheses onto the phylogeny. + Each selective regime is named (small, medium, large, etc.). + Each column corresponds to a different painting of the selective regimes, and thus to a different hypothesis. + In this example, there are 3 alternative models (see Butler & King 2004): \code{OU.4} is 4-regime model, \code{OU.3} is 3-regime model (all ancestors are medium), \code{OU.LP} is the linear parsimony model. + } } -\references{ - Lazell, J. D. 1972. - The anoles (Sauria: Iguanidae) of the Lesser Antilles. - Bull. Mus. Comp. Zool., 143:1-115. - - Losos, J. B. 1990. - A phylogenetic analysis of character displacement in Caribbean Anolis lizards. - Evolution, 44:558-569. } \examples{ -data(bimac) tree <- with(bimac,ouchtree(node,ancestor,time/max(time),species)) plot(tree,node.names=TRUE) -print(h1 <- brown(log(bimac['size']),tree)) + +h1 <- brown(log(bimac['size']),tree) +h1 plot(h1) -print(h2 <- hansen(log(bimac['size']),tree,bimac['OU.1'],sqrt.alpha=1,sigma=1)) + +h2 <- hansen(log(bimac['size']),tree,bimac['OU.1'],sqrt.alpha=1,sigma=1) +h2 plot(h2) -print(h3 <- hansen(log(bimac['size']),tree,bimac['OU.3'],sqrt.alpha=1,sigma=1)) + +h3 <- hansen(log(bimac['size']),tree,bimac['OU.3'],sqrt.alpha=1,sigma=1) +h3 plot(h3) -print(h4 <- hansen(log(bimac['size']),tree,bimac['OU.4'],sqrt.alpha=1,sigma=1)) + +h4 <- hansen(log(bimac['size']),tree,bimac['OU.4'],sqrt.alpha=1,sigma=1) +h4 plot(h4) + h5 <- hansen(log(bimac['size']),tree,bimac['OU.LP'],sqrt.alpha=1,sigma=1,reltol=1e-5) -print(h5 <- update(h5,method='subplex',reltol=1e-11,parscale=c(0.1,0.1),hessian=TRUE)) +h5 <- update(h5,method='subplex',reltol=1e-11,parscale=c(0.1,0.1),hessian=TRUE) +h5 + simdat <- simulate(h5,nsim=10) hsim <- update(h5,data=simdat[[1]]) -print(summary(hsim)) +summary(hsim) bsim <- update(h1,data=simdat[[1]]) -print(summary(bsim)) +summary(bsim) +} +\references{ +\Lazell1972 + +\Losos1990 +} +\author{ +Marguerite A. Butler and Aaron A. King } -\author{Marguerite A. Butler and Aaron A. King } \keyword{models} diff --git a/man/bootstrap.Rd b/man/bootstrap.Rd new file mode 100644 index 0000000..4b86a13 --- /dev/null +++ b/man/bootstrap.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bootstrap.R +\name{bootstrap} +\alias{bootstrap} +\alias{bootstrap,browntree-method} +\alias{bootstrap,hansentree-method} +\alias{bootstrap,missing-method} +\alias{bootstrap,ANY-method} +\title{Bootstrapping for uncertainty quantification} +\usage{ +\S4method{bootstrap}{browntree}(object, nboot = 200, seed = NULL, ...) + +\S4method{bootstrap}{hansentree}(object, nboot = 200, seed = NULL, ...) + +\S4method{bootstrap}{missing}(object, ...) + +\S4method{bootstrap}{ANY}(object, ...) +} +\arguments{ +\item{object}{A fitted model object.} + +\item{nboot}{integer; number of bootstrap replicates.} + +\item{seed}{integer; setting \code{seed} to a non-\code{NULL} value allows one to fix the random seed (see \code{\link[ouch:simulate]{simulate}}).} + +\item{...}{Additional arguments are passed to \code{\link[ouch:update]{update}}.} +} +\description{ +Generic bootstrapping for \pkg{ouch} models. +} +\details{ +\code{bootstrap} performs a parametric bootstrap for estimation of confidence intervals. +} diff --git a/man/brown.Rd b/man/brown.Rd index 4921746..d3e42c8 100644 --- a/man/brown.Rd +++ b/man/brown.Rd @@ -1,32 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brown.R \name{brown} \alias{brown} -\title{Fit a Brownian-motion model of evolution along a phylogenetic tree} -\description{ - \code{brown} fits a Brownian motion model to the given data. -} +\title{Phylogenetic Brownian motion models} \usage{ brown(data, tree) } \arguments{ - \item{data}{Phenotypic data for extant species, i.e., at the terminal - ends of the phylogenetic tree. This can either be a numeric vector - or a list. If it is a numeric vector, there must be one entry for - every node. If it is a list, it must consist entirely of numeric - vectors, each of which has one entry per node. A data-frame is - coerced to a list.} - \item{tree}{A phylogenetic tree, specified as an \code{ouchtree} - object.} +\item{data}{Phenotypic data for extant species, i.e., at the terminal ends +of the phylogenetic tree. This can either be a numeric vector or a list. +If it is a numeric vector, there must be one entry for every node. If it is +a list, it must consist entirely of numeric vectors, each of which has one +entry per node. A data-frame is coerced to a list.} + +\item{tree}{A phylogenetic tree, specified as an \code{ouchtree} object.} } \value{ - \code{brown} returns an object of class \code{browntree}. See - \code{\link{browntree-class}} for information on the methods of this - class. +\code{brown} returns an object of class \code{browntree}. +} +\description{ +The function \code{brown} creates a \code{browntree} object by fitting a +Brownian-motion model to data. } -\seealso{\code{\link{ouchtree}}, \code{\link{browntree}}} \references{ - Butler, M.A. and A.A. King (2004) Phylogenetic comparative analysis: a - modeling approach for adaptive evolution. American Naturalist - 164:683-695. +\Butler2004 +} +\seealso{ +\code{\link{ouchtree}}, \code{\link{hansen}}, \code{\link{bimac}}, \code{\link{anolis.ssd}} +} +\author{ +Aaron A. King } -\author{Aaron A. King } \keyword{models} diff --git a/man/browntree-class.Rd b/man/browntree-class.Rd deleted file mode 100644 index 7d9fc38..0000000 --- a/man/browntree-class.Rd +++ /dev/null @@ -1,76 +0,0 @@ -\name{browntree} -\docType{class} -\alias{browntree} -\alias{browntree-class} -\alias{logLik,browntree-method} -\alias{logLik-browntree} -\alias{summary,browntree-method} -\alias{summary-browntree} -\alias{print,browntree-method} -\alias{print-browntree} -\alias{show,browntree-method} -\alias{show-browntree} -\alias{plot,browntree-method} -\alias{plot-browntree} -\alias{simulate,browntree-method} -\alias{simulate-browntree} -\alias{update,browntree-method} -\alias{update-browntree} -\alias{bootstrap,browntree-method} -\alias{bootstrap-browntree} -\alias{coef,browntree-method} -\alias{coef-browntree} -\alias{coerce,browntree,data.frame-method} -\title{Fitted phylogenetic Brownian motion model} -\description{ - A fitted phylogenetic Brownian-motion model object. -} -\details{ - The function \code{brown} creates a \code{browntree} object by fitting - a Brownian-motion model to data. -} -\section{Methods}{ - \describe{ - \item{\code{plot()}}{plots the tree.} - \item{\code{print()}}{ - prints the tree as a table, along with the coefficients of the fitted model and diagnostic information. - } - \item{\code{show()}}{displays the fitted \code{browntree} object.} - \item{\code{summary()}}{ - displays information on the call, the fitted coefficients, and model selection statistics. - } - \item{coerce}{A \code{browntree} object can be coerced to a data-frame - via \code{as(object,"data.frame")}.} - \item{\code{coef(object,\dots)}}{extracts the coefficients of the - fitted model. This is a list with three elements: - \describe{ - \item{\code{sigma}:}{the coefficients of the sigma matrix.} - \item{\code{theta}:}{a list of the estimated optima, one per character.} - \item{\code{sigma..sq.matrix}:}{the sigma-squared matrix itself.} - } - } - \item{\code{logLik(object,\dots)}}{ - extracts the log likelihood of the fitted model. - } - \item{\code{update(object, \dots)}}{ - refits the model. - \code{object} is the \code{browntree} object. - Additional arguments (in \code{\dots}) replace the corresponding arguments in the original call. - } - \item{\code{bootstrap(object, nboot = 200, seed = NULL, \dots)}}{ - performs a parametric bootstrap for estimation of confidence intervals. - \code{object} is the \code{browntree} object. - \code{nboot} is the number of bootstraps. - \code{seed} allows one to fix the random seed (see \code{simulate} below). - Additional arguments (in \code{\dots}) are passed to \code{update}. - } - \item{\code{simulate(object, nsim = 1, seed = NULL, \dots)}}{ - generates random deviates from the fitted model. - \code{object} is the \code{browntree} object, \code{nsim} is the desired number of replicates, and \code{seed} is (optionally) the random seed to use. - \code{simulate} returns a list of data-frames, each comparable to the original data. - } - } -} -\author{Aaron A. King kingaa at umich dot edu} -\seealso{\code{\link{ouchtree}}, \code{\link{brown}}} -\keyword{models} diff --git a/man/coef.Rd b/man/coef.Rd new file mode 100644 index 0000000..10a7e40 --- /dev/null +++ b/man/coef.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brown.R, R/hansen.R, R/coef.R +\name{coef,browntree-method} +\alias{coef,browntree-method} +\alias{coef,hansentree-method} +\alias{coef} +\title{Model coefficients} +\usage{ +\S4method{coef}{browntree}(object, ...) + +\S4method{coef}{hansentree}(object, ...) +} +\arguments{ +\item{object}{fitted model object.} + +\item{...}{additional arguments, ignored.} +} +\value{ +\code{coef} applied to a \code{browntree} object extracts a list with three elements: +\describe{ + \item{\code{sigma}}{the coefficients of the sigma matrix.} + \item{\code{theta}}{a list of the estimated optima, one per character.} + \item{\code{sigma..sq.matrix}}{the sigma-squared matrix itself.} +} + +\code{coef} applied to a \code{hansentree} object returns a named list containing the estimated \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices(given as the \code{alpha.matrix} and \code{sigma.sq.matrix} elements, respectively) but also the MLE returned by the optimizer +(as \code{sqrt.alpha} and \code{sigma}, respectively). +\strong{The latter elements should not be interpreted, but can be used to restart the algorithm, etc.} +} +\description{ +\code{coef} extracts the parameters from a fitted model object. +} diff --git a/man/hansen.Rd b/man/hansen.Rd index 1ad9841..3536bbe 100644 --- a/man/hansen.Rd +++ b/man/hansen.Rd @@ -1,135 +1,146 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hansen.R \name{hansen} \alias{hansen} -\title{Hansen model of evolution along a phylogenetic tree} -\description{ - Fits the Ornstein-Uhlenbeck-based Hansen model to data. - The fitting is done using \code{optim} or \code{subplex}. -} +\title{Ornstein-Uhlenbeck models of trait evolution} \usage{ -hansen(data, tree, regimes, sqrt.alpha, sigma, - fit = TRUE, - method = c("Nelder-Mead","subplex","BFGS","L-BFGS-B"), - hessian = FALSE, \dots) +hansen( + data, + tree, + regimes, + sqrt.alpha, + sigma, + fit = TRUE, + method = c("Nelder-Mead", "subplex", "BFGS", "L-BFGS-B"), + hessian = FALSE, + ... +) } \arguments{ - \item{data}{ - Phenotypic data for extant species, i.e., species at the terminal twigs of the phylogenetic tree. - This can either be a single named numeric vector, a list of \code{nchar} named vectors, or a data-frame containing \code{nchar} data variables. - There must be an entry per variable for every node in the tree; use \code{NA} to represent missing data. - If the data are supplied as one or more named vectors, the names attributes are taken to correspond to the node names specified when the \code{ouchtree} was constructed (see \code{\link{ouchtree}}). - If the data are supplied as a data-frame, the rownames serve that purpose. - } - \item{tree}{ - A phylogenetic tree, specified as an \code{ouchtree} object. - } - \item{regimes}{ - A vector of codes, one for each node in the tree, specifying the selective regimes hypothesized to have been operative. - Corresponding to each node, enter the code of the regime hypothesized for the branch segment terminating in that node. - For the root node, because it has no branch segment terminating on it, the regime specification is irrelevant. - If there are \code{nchar} quantitative characters, then one can specify a single set of \code{regimes} for all characters or a list of \code{nchar} regime specifications, one for each character. - } - \item{sqrt.alpha, sigma}{ - These are used to initialize the optimization algorithm. - The selection strength matrix \eqn{\alpha}{alpha} and the random drift variance-covariance matrix \eqn{\sigma^2}{sigma^2} are parameterized by their matrix square roots. - Specifically, these initial guesses are each packed into lower-triangular matrices (column by column). - The product of this matrix with its transpose is the \eqn{\alpha}{alpha} or \eqn{\sigma^2}{sigma^2} matrix. - See Details, below. - } - \item{fit}{ - If \code{fit=TRUE}, then the likelihood will be maximized. - If \code{fit=FALSE}, the likelihood will be evaluated at the specified values of \code{sqrt.alpha} and \code{sigma}; - the optima \code{theta} will be returned as well. - } - \item{method}{ - The method to be used by the optimization algorithm, \code{optim}. - See \code{\link[subplex]{subplex}} and \code{\link{optim}} for information on the available options. - } - \item{hessian}{ - If \code{hessian=TRUE}, then the Hessian matrix will be computed by \code{optim}. - } - \item{\dots}{ - Additional arguments will be passed as \code{control} options to \code{optim} or \code{subplex}. - See \code{\link{optim}} and \code{\link[subplex]{subplex}} for information on the available options. - } +\item{data}{Phenotypic data for extant species, i.e., species at the terminal twigs of the phylogenetic tree. +This can either be a single named numeric vector, a list of \code{nchar} named vectors, or a data frame containing \code{nchar} data variables. +There must be an entry per variable for every node in the tree; use \code{NA} to represent missing data. +If the +data are supplied as one or more named vectors, the names attributes are taken to correspond to the node names specified when the \code{ouchtree} was constructed (see \code{\link{ouchtree}}). +If the data are supplied as a +data-frame, the rownames serve that purpose.} + +\item{tree}{A phylogenetic tree, specified as an \code{ouchtree} object.} + +\item{regimes}{A vector of codes, one for each node in the tree, specifying the selective regimes hypothesized to have been operative. +Corresponding to each node, enter the code of the regime hypothesized for the branch segment terminating in that node. +For the root node, because it has no branch segment terminating on it, the regime specification is irrelevant. +If there are \code{nchar} quantitative characters, then one can specify a single set of \code{regimes} for all characters or a list of \code{nchar} regime specifications, one for each character.} + +\item{sqrt.alpha, sigma}{These are used to initialize the optimization algorithm. +The selection strength matrix \eqn{\alpha}{alpha} and the random drift variance-covariance matrix \eqn{\sigma^2}{sigma^2} are parameterized by their matrix square roots. +Specifically, these initial guesses are each packed into lower-triangular matrices (column by column). +The product of this matrix with its transpose is the \eqn{\alpha}{alpha} or \eqn{\sigma^2}{sigma^2} matrix. +See Details.} + +\item{fit}{If \code{fit=TRUE}, then the likelihood will be maximized. +If \code{fit=FALSE}, the likelihood will be evaluated at the specified values of \code{sqrt.alpha} and \code{sigma}; +the optima \code{theta} will be returned as well.} + +\item{method}{The method to be used by the optimization algorithm. +See \code{\link[subplex]{subplex}} and \code{\link{optim}} for information on the available options.} + +\item{hessian}{If \code{hessian=TRUE}, then the Hessian matrix will be computed by \code{optim}.} + +\item{\dots}{Additional arguments will be passed as \code{control} options to \code{optim} or \code{subplex}. +See \code{\link{optim}} and \code{\link[subplex:subplex]{subplex}} for information on the available options.} } \value{ - \code{hansen} returns an object of class \code{hansentree}. - For details on the methods of that class, see \code{\link{hansentree}}. +\code{hansen} returns an object of class \code{hansentree}. +For details on the methods of that class, see \code{\link{hansen}}. +} +\description{ +The function \code{hansen} fits an Ornstein-Uhlenbeck model to data. +The fitting is done using \code{optim} or \code{subplex}. } \details{ - The Hansen model for the evolution of a multivariate trait \eqn{X} along a lineage can be written as a stochastic differential equation (Ito diffusion) \deqn{dX=\alpha(\theta(t)-X(t))dt+\sigma dB(t),}{dX = alpha (theta(t)-X(t)) dt + sigma dB(t),} where \eqn{t} is time along the lineage, \eqn{\theta(t)}{theta(t)} is the optimum trait value, \eqn{B(t)} is a standard Wiener process (Brownian motion), and \eqn{\alpha}{alpha} and \eqn{\sigma}{sigma} are matrices quantifying, respectively, the strength of selection and random drift. - Without loss of generality, one can assume \eqn{\sigma}{sigma} is lower-triangular. - This is because only the infinitesimal variance-covariance matrix \eqn{\sigma^2=\sigma\sigma^T}{sigma^2 = sigma\%*\%transpose(sigma)} is identifiable, and for any admissible variance-covariance matrix, we can choose \eqn{\sigma}{sigma} to be lower-triangular. - Moreover, if we view the basic model as describing evolution on a fitness landscape, then \eqn{\alpha}{alpha} will be symmetric and if we further restrict ourselves to the case of stabilizing selection, \eqn{\alpha}{alpha} will be positive definite as well. - We make these assumptions and therefore can assume that the matrix \eqn{\alpha}{alpha} has a lower-triangular square root. - - The \code{hansen} code uses unconstrained numerical optimization to maximize the likelihood. - To do this, it parameterizes the \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices in a special way: - each matrix is parameterized by \code{nchar*(nchar+1)/2} parameters, where \code{nchar} is the number of quantitative characters. - Specifically, the parameters initialized by the \code{sqrt.alpha} argument of \code{hansen} are used to fill the nonzero entries of a lower-triangular matrix (in column-major order), which is then multiplied by its transpose to give the selection-strength matrix. - The parameters specified in \code{sigma} fill the nonzero entries in the lower triangular \eqn{\sigma}{sigma} matrix. - When \code{hansen} is executed, the numerical optimizer maximizes the likelihood over these parameters. - The \code{print}, \code{show}, and \code{summary} methods for the resulting \code{hansentree} object display (among other things) the estimated \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices. - The \code{coef} method extracts a named list containing not only these matrices (given as the \code{alpha.matrix} and \code{sigma.sq.matrix} elements) but also the MLEs returned by the optimizer (as \code{sqrt.alpha} and \code{sigma}, respectively). - \strong{The latter elements should not be interpreted, but can be used to restart the algorithm, etc.} +The Hansen model for the evolution of a multivariate trait \eqn{X} along a lineage can be written as a stochastic differential equation (Ito diffusion) +\deqn{dX=\alpha(\theta(t)-X(t))dt+\sigma dB(t),}{dX = alpha (theta(t)-X(t)) dt + sigma dB(t),} +where \eqn{t} is time along the lineage, +\eqn{\theta(t)}{theta(t)} is the optimum trait value, \eqn{B(t)} is a standard Wiener process (Brownian motion), +and \eqn{\alpha}{alpha} and \eqn{\sigma}{sigma} are matrices +quantifying, respectively, the strength of selection and random drift. +Without loss of generality, one can assume \eqn{\sigma}{sigma} is lower-triangular. +This is because only the infinitesimal variance-covariance matrix +\eqn{\sigma^2=\sigma\sigma^T}{sigma^2 = sigma\%*\%transpose(sigma)} +is identifiable, and for any admissible variance-covariance matrix, we can choose \eqn{\sigma}{sigma} to be lower-triangular. +Moreover, if we view the basic model as describing evolution on a fitness landscape, then \eqn{\alpha}{alpha} will be symmetric. +If we further restrict ourselves to the case of stabilizing selection, \eqn{\alpha}{alpha} will be positive definite as well. +We make these assumptions and therefore can assume that the matrix \eqn{\alpha}{alpha} has a lower-triangular square root. + +The \code{hansen} code uses unconstrained numerical optimization to maximize the likelihood. +To do this, it parameterizes the \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices in a special way: +each matrix is parameterized by \code{nchar*(nchar+1)/2} parameters, where \code{nchar} is the number of quantitative characters. +Specifically, the parameters initialized by the \code{sqrt.alpha} argument of \code{hansen} are used +to fill the nonzero entries of a lower-triangular matrix (in column-major order), +which is then multiplied by its transpose to give the selection-strength matrix. +The parameters specified in \code{sigma} fill the nonzero entries in the lower triangular \eqn{\sigma}{sigma} matrix. +When \code{hansen} is executed, the numerical optimizer maximizes the likelihood over these parameters. } \examples{ \dontrun{ -if (library(geiger,logical.return=TRUE)) { - + if (library(geiger,logical.return=TRUE)) { + ### an example data set (Darwin's finches) -data(geospiza) -str(geospiza) -sapply(geospiza,class) + data(geospiza) + str(geospiza) + sapply(geospiza,class) ### check the correspondence between data and tree tips: -print(nc <- with(geospiza,name.check(geospiza.tree,geospiza.data))) + print(nc <- with(geospiza,name.check(geospiza.tree,geospiza.data))) ### looks like one of the terminal twigs has no data associated ### drop that tip: -tree <- with(geospiza,drop.tip(geospiza.tree,nc$tree_not_data)) -dat <- as.data.frame(geospiza$dat) + tree <- with(geospiza,drop.tip(geospiza.tree,nc$tree_not_data)) + dat <- as.data.frame(geospiza$dat) ### make an ouchtree out of the phy-format tree -ot <- ape2ouch(tree) + ot <- ape2ouch(tree) ### merge data with tree info -otd <- as(ot,"data.frame") + otd <- as(ot,"data.frame") ### in these data, it so happens that the rownames correspond to node names ### we will exploit this correspondence in the 'merge' operation: -dat$labels <- rownames(dat) -otd <- merge(otd,dat,by="labels",all=TRUE) -rownames(otd) <- otd$nodes -print(otd) + dat$labels <- rownames(dat) + otd <- merge(otd,dat,by="labels",all=TRUE) + rownames(otd) <- otd$nodes + print(otd) ### this data-frame now contains the data as well as the tree geometry ### now remake the ouch tree -ot <- with(otd,ouchtree(nodes=nodes,ancestors=ancestors,times=times,labels=labels)) + ot <- with(otd,ouchtree(nodes=nodes,ancestors=ancestors,times=times,labels=labels)) -b1 <- brown(tree=ot,data=otd[c("tarsusL","beakD")]) -summary(b1) + b1 <- brown(tree=ot,data=otd[c("tarsusL","beakD")]) + summary(b1) ### evaluate an OU model with a single, global selective regime -otd$regimes <- as.factor("global") -h1 <- hansen( - tree=ot, - data=otd[c("tarsusL","beakD")], - regimes=otd["regimes"], - sqrt.alpha=c(1,0,1), - sigma=c(1,0,1), - maxit=10000 - ) -summary(h1) + otd$regimes <- as.factor("global") + h1 <- hansen( + tree=ot, + data=otd[c("tarsusL","beakD")], + regimes=otd["regimes"], + sqrt.alpha=c(1,0,1), + sigma=c(1,0,1), + maxit=10000 + ) + summary(h1) + } } } +\references{ +\Butler2004 + +\Cressler2015 } \seealso{ - \code{\link{ouchtree}}, \code{\link{hansentree}}, \code{\link{optim}}, \code{\link{bimac}}, \code{\link{anolis.ssd}} +\code{\link{ouchtree}}, \code{\link{brown}}, \code{\link[stats:optim]{optim}}, \code{\link[subplex:subplex]{subplex}}, \code{\link{bimac}}, \code{\link{anolis.ssd}} } -\references{ - Butler, M.A. and A.A. King (2004) - Phylogenetic comparative analysis: a modeling approach for adaptive evolution. - American Naturalist 164:683-695. +\author{ +Aaron A. King } -\author{Aaron A. King } \keyword{models} diff --git a/man/hansentree-class.Rd b/man/hansentree-class.Rd deleted file mode 100644 index 64f4167..0000000 --- a/man/hansentree-class.Rd +++ /dev/null @@ -1,17 +0,0 @@ -\name{hansentree} -\docType{class} -\title{hansen tree object} -\alias{hansentree} -\alias{hansentree-class} -\description{ - A fitted phylogenetic Hansen-model object. -} -\details{ - The function \code{hansen} creates a \code{hansentree} object by fitting a Hansen model to data. -} -\section{Methods}{ - See \link{hansentree-methods}. -} -\author{Aaron A. King kingaa at umich dot edu} -\seealso{\code{\link{ouchtree}}, \code{\link{hansen}}, \link{hansentree-methods}} -\keyword{models} diff --git a/man/hansentree-methods.Rd b/man/hansentree-methods.Rd deleted file mode 100644 index e3b3f5f..0000000 --- a/man/hansentree-methods.Rd +++ /dev/null @@ -1,115 +0,0 @@ -\name{hansentree-methods} -\docType{methods} -\alias{hansentree-methods} -\alias{logLik,hansentree-method} -\alias{logLik-hansentree} -\alias{summary,hansentree-method} -\alias{summary-hansentree} -\alias{print,hansentree-method} -\alias{print-hansentree} -\alias{show,hansentree-method} -\alias{show-hansentree} -\alias{plot,hansentree-method} -\alias{plot-hansentree} -\alias{update,hansentree-method} -\alias{update-hansentree} -\alias{bootstrap} -\alias{bootstrap,hansentree-method} -\alias{bootstrap-hansentree} -\alias{simulate,hansentree-method} -\alias{simulate-hansentree} -\alias{coef,hansentree-method} -\alias{coef-hansentree} -\alias{as,hansentree-method} -\alias{coerce,hansentree,data.frame-method} -\title{Methods of the "hansentree" class} -\description{Methods of the "hansentree" class.} -\usage{ -\S4method{logLik}{hansentree}(object) -\S4method{coef}{hansentree}(object, \dots) -\S4method{summary}{hansentree}(object, \dots) -\S4method{show}{hansentree}(object) -\S4method{print}{hansentree}(x, \dots) -\S4method{plot}{hansentree}(x, y, regimes, \dots) -\S4method{simulate}{hansentree}(object, nsim = 1, seed = NULL, \dots) -\S4method{update}{hansentree}(object, data, regimes, sqrt.alpha, sigma, \dots) -\S4method{bootstrap}{hansentree}(object, nboot = 200, seed = NULL, \dots) -\S4method{as}{hansentree}(object, class) -\S4method{coerce}{hansentree,data.frame}(from, to = "data.frame", strict = TRUE) -} -\arguments{ - \item{object}{The \code{hansentree} object.} - \item{x}{the \code{hansentree} object.} - \item{class}{ - character; - name of the class to which \code{object} should be coerced. - } - \item{from, to}{ - the classes betwen which coercion should be performed. - } - \item{nsim}{ - The number of simulations to perform. - } - \item{nboot}{ - The number of boostraps to perform. - } - \item{seed}{The random seed to use in simulations.} - \item{regimes, sqrt.alpha, sigma}{See \code{\link{hansen}}.} - \item{data}{see \code{\link{hansen}}.} - \item{y, strict}{Ignored.} - \item{\dots}{ - Further arguments (either ignored or passed to underlying functions). - In the case of \code{update}, these replace the corresponding arguments in the original call. - } -} -\section{Methods}{ - \describe{ - \item{\code{plot()}}{ - plots the tree, with branches colored according to the selective regimes. - See \link{plot-ouchtree} for more details. - } - \item{\code{print()}}{ - prints the tree as a table, along with the coefficients of the fitted model and diagnostic information. - } - \item{\code{show()}}{ - displays the fitted \code{hansentree} object. - } - \item{\code{summary()}}{ - displays information on the call, the fitted coefficients, and model selection statistics. - } - \item{coerce}{ - A \code{hansentree} object can be coerced to a data-frame via \code{as(object,"data.frame")}. - } - \item{\code{coef()}}{ - extracts the coefficients of the fitted model. - This is a list with five elements: - \describe{ - \item{\code{sqrt.alpha}:}{the coefficients that parameterize the alpha matrix.} - \item{\code{sigma}:}{the coefficients that parameterize the sigma matrix.} - \item{\code{theta}:}{ - a list of the estimated optima, one per character. - Each element of the list is a vector containing one optimal value per regime. - } - \item{\code{alpha.matrix}:}{the alpha matrix itself.} - \item{\code{sigma.sq.matrix}:}{the sigma-squared matrix itself.} - } - } - \item{\code{logLik()}}{ - extracts the log likelihood of the fitted model. - } - \item{\code{update()}}{ - refines the model fit. - } - \item{\code{bootstrap()}}{ - performs a parametric bootstrap for confidence intervals. - } - \item{\code{simulate()}}{ - generates random deviates from the fitted model. - \code{object} is the \code{hansentree} object, \code{nsim} is the desired number of replicates, and \code{seed} is (optionally) the random seed to use. - \code{simulate} returns a list of data-frames, each comparable to the original data. - } - } -} -\author{Aaron A. King kingaa at umich dot edu} -\seealso{\code{\link{ouchtree}}, \code{\link{hansen}}} -\keyword{models} diff --git a/man/loglik.Rd b/man/loglik.Rd new file mode 100644 index 0000000..0cf1bb1 --- /dev/null +++ b/man/loglik.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brown.R, R/hansen.R, R/loglik.R +\name{logLik,browntree-method} +\alias{logLik,browntree-method} +\alias{logLik,hansentree-method} +\alias{loglik} +\title{Log likelihood of a fitted model} +\usage{ +\S4method{logLik}{browntree}(object) + +\S4method{logLik}{hansentree}(object) +} +\arguments{ +\item{object}{fitted model object} +} +\value{ +\code{logLik} returns a numeric value. +} +\description{ +\code{logLik} extracts the log likelihood from a fitted model object. +} diff --git a/man/macros/citations.Rd b/man/macros/citations.Rd new file mode 100644 index 0000000..0e9b3ba --- /dev/null +++ b/man/macros/citations.Rd @@ -0,0 +1,13 @@ +% unlike in latex, new lines within \newcommand end the command (even using %) + +\newcommand{\Butler2000}{Butler, M. A., T. W. Schoener, and J. B. Losos. 2000. The relationship between sexual size dimorphism and habitat use in Greater Antillean Anolis lizards. Evolution, 54:259--272.} + +\newcommand{\Williams1972}{Williams, E. E. 1972. The origin of faunas. Evolution of lizard congeners in a complex island fauna: a trial analysis. Evol. Biol., 6:47--89.} + +\newcommand{\Butler2004}{Butler, M.A. and A.A. King. 2004. Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am. Nat. 164:683--695.} + +\newcommand{\Lazell1972}{Lazell, J. D. 1972. The anoles (Sauria: Iguanidae) of the Lesser Antilles. Bull. Mus. Comp. Zool., 143:1--115.} + +\newcommand{\Losos1990}{Losos, J. B. 1990. A phylogenetic analysis of character displacement in Caribbean Anolis lizards. Evolution, 44:558--569.} + +\newcommand{\Cressler2015}{Cressler, C. E., Butler, M. A., and King, A. A. 2015. Detecting adaptive evolution in phylogenetic comparative analysis using the Ornstein-Uhlenbeck model. Syst. Biol., 64:953--968.} diff --git a/man/ouch-package.Rd b/man/ouch-package.Rd deleted file mode 100644 index ca2c2f2..0000000 --- a/man/ouch-package.Rd +++ /dev/null @@ -1,45 +0,0 @@ -\name{ouch-package} -\docType{package} -\alias{ouch-package} -\title{Ornstein-Uhlenbeck methods for comparative phylogenetic hypotheses} -\description{The \pkg{ouch} package provides facilities for Ornstein-Uhlenbeck based methods of phylogenetic comparative analysis.} -\section{Classes}{ - The basic class, \code{ouchtree}, is provided to encode a phylogenetic tree. - Plot and print methods are provided. - - The class \code{browntree} derives from class \code{ouchtree} and encodes the results of fitting a Brownian Motion model to data. - - The class \code{hansentree} derives from class \code{ouchtree} and encodes the results of fitting a Hansen model to data. -} -\section{Detailed Documentation}{ - \describe{ - \item{ouchtree}{\code{\link{ouchtree}}} - \item{Brownian motion methods}{ - \code{\link{brown}}, - \code{\link{browntree-class}} - } - \item{Ornstein-Uhlenbeck methods}{ - \code{\link{hansen}}, - \code{\link{hansentree-class}} - } - \item{simulate methods}{ - \code{\link{simulate-browntree}}, - \code{\link{simulate-hansentree}} - } - \item{plot methods}{ - \code{\link{plot-browntree}}, - \code{\link{plot-hansentree}} - } - \item{\code{ape2ouch}}{ - Convert a tree in \pkg{ape} format to \pkg{ouch} format. - See \code{\link{ape2ouch}} - } - } -} -\references{ - Butler, M.A. and A.A. King (2004) - Phylogenetic comparative analysis: a modeling approach for adaptive evolution. - American Naturalist 164:683--695. -} -\author{Aaron A. King (kingaa at umich dot edu)} -\keyword{models} diff --git a/man/ouchtree-class.Rd b/man/ouchtree-class.Rd deleted file mode 100644 index 7df5286..0000000 --- a/man/ouchtree-class.Rd +++ /dev/null @@ -1,69 +0,0 @@ -\name{ouchtree} -\docType{class} -\alias{ouchtree} -\alias{ouchtree-class} -\alias{show,ouchtree-method} -\alias{show-ouchtree} -\alias{plot,ouchtree-method} -\alias{plot-ouchtree} -\alias{print,ouchtree-method} -\alias{print-ouchtree} -\alias{coerce,ouchtree,data.frame-method} -\title{Phylogenetic tree object in 'ouch' format.} -\description{ - An object containing a phylogenetic tree in a form suitable for using \pkg{ouch} methods. -} -\usage{ -ouchtree(nodes, ancestors, times, labels = as.character(nodes)) -} -\arguments{ - \item{nodes}{ - A character vector giving the name of each node. These are used internally and must be unique. - } - \item{ancestors}{ - Specification of the topology of the phylogenetic tree. - This is in the form of a character vector specifying the name (as given in the \code{nodes} argument) of the immediate ancestor of each node. - In particular, the i-th name is that of the ancestor of the i-th node. - The root node is distinguished by having no ancestor (i.e., NA). - } - \item{times}{ - A vector of nonnegative numbers, one per node in the tree, specifying the time at which each node is located. - Time should be increasing from the root node to the terminal twigs. - } - \item{labels}{ - Optional vector of node labels. - These will be used in plots to label nodes. - It is not necessary that these be unique. - } -} -\details{ - \code{ouchtree} creates an \code{ouchtree} object. - This contains the topology, branch times, and epochs. - It also holds (optionally) names of taxa for display purposes. -} -\section{Methods}{ - \describe{ - \item{\code{plot(tree,regimes=NULL,node.names=FALSE,legend=TRUE,...)}}{ - displays the phylogenetic tree graphically. - } - \item{\code{print()}}{ - displays the tree in table form. - } - \item{\code{show()}}{ - displays the tree. - } - \item{coerce}{ - An \code{ouchtree} object can be coerced to a data-frame via \code{as(object,"data.frame")}. - } - } -} -\author{Aaron A. King kingaa at umich dot edu} -\seealso{\code{ouchtree}, \code{ape2ouch}, \code{brown}, \code{hansen}} -\examples{ -data(bimac) -tree <- with(bimac,ouchtree(nodes=node,ancestors=ancestor,times=time,labels=species)) -plot(tree) -plot(tree,node.names=TRUE) -print(tree) -} -\keyword{models} diff --git a/man/ouchtree.Rd b/man/ouchtree.Rd new file mode 100644 index 0000000..7e1ed3c --- /dev/null +++ b/man/ouchtree.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ouchtree.R +\name{ouchtree} +\alias{ouchtree} +\alias{ouchtree-class} +\title{Phylogenetic tree object in 'ouch' format.} +\usage{ +ouchtree(nodes, ancestors, times, labels = as.character(nodes)) +} +\arguments{ +\item{nodes}{A character vector giving the name of each node. +These are used internally and must be unique.} + +\item{ancestors}{Specification of the topology of the phylogenetic tree. +This is in the form of a character vector specifying the name +(as given in the \code{nodes} argument) +of the immediate ancestor of each node. +In particular, the i-th name is that of the ancestor of the i-th node. +The root node is distinguished by having no ancestor (i.e., \code{NA}).} + +\item{times}{A vector of nonnegative numbers, one per node in the tree, +specifying the time at which each node is located. +Time should be increasing from the root node to the terminal twigs.} + +\item{labels}{Optional vector of node labels. +These will be used in plots to label nodes. +It is not necessary that these be unique.} +} +\description{ +An object containing a phylogenetic tree in a form suitable for using \pkg{ouch} methods. +} +\details{ +\code{ouchtree} creates an \code{ouchtree} object. +This contains the topology, branch times, and epochs. +It also (optionally) holds names of taxa for display purposes. +} +\examples{ +tree <- with( + bimac, + ouchtree(nodes=node,ancestors=ancestor,times=time,labels=species) +) +tree + +plot(tree) +plot(tree,node.names=TRUE) +} +\seealso{ +\code{ouchtree}, \code{ape2ouch}, \code{brown}, \code{hansen} +} +\author{ +Aaron A. King +} +\keyword{models} diff --git a/man/package.Rd b/man/package.Rd new file mode 100644 index 0000000..00f7409 --- /dev/null +++ b/man/package.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/package.R +\docType{package} +\name{ouch-package} +\alias{ouch-package} +\title{Ornstein-Uhlenbeck methods for comparative phylogenetic hypotheses} +\description{ +The \pkg{ouch} package provides facilities for phylogenetic comparative analysis based on Ornstein-Uhlenbeck models of trait evolution along a phylogeny. +Multivariate data and complex adaptive hypotheses are supported. +} +\section{Classes}{ + +The basic class, \code{ouchtree}, is provided to encode a +phylogenetic tree. Plot and print methods are provided. + +The class \code{browntree} derives from class \code{ouchtree} and encodes +the results of fitting a Brownian Motion model to data. + +The class \code{hansentree} derives from class \code{ouchtree} and encodes +the results of fitting a Hansen model to data. +} + +\section{Detailed Documentation}{ + +\describe{ + \item{Phylogenies in \pkg{ouch} format}{\code{\link{ouchtree}}, \code{\link{ape2ouch}}} + \item{Brownian motion models}{\code{\link{brown}}} + \item{Ornstein-Uhlenbeck models}{\code{\link{hansen}}} + \item{Simulation of models}{\code{\link[ouch:simulate]{simulate}}} + \item{Display of data}{\code{\link[ouch:plot]{plot}}} +} +} + +\references{ +\Butler2004 + +\Cressler2015 +} +\author{ +Aaron A. King +} +\keyword{models} diff --git a/man/paint.Rd b/man/paint.Rd index 6f9621c..64f0d89 100644 --- a/man/paint.Rd +++ b/man/paint.Rd @@ -1,53 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/paint.R \name{paint} \alias{paint} -\title{Painting selective regimes on a phylogenetic tree.} -\description{ - Function to paint selective regimes on a phylogenetic tree. -} +\title{Painting regimes on a phylogenetic tree} \usage{ paint(tree, subtree, branch, which = 1) } \arguments{ - \item{tree}{ - An object of class \code{ouchtree}. - } - \item{subtree}{ - An optional named vector specifying the root nodes of subtrees. - Each branch that descends from this node will be painted with the specified regime. - } - \item{branch}{ - An optional named vector specifying the end nodes of branches. - The unique branch that terminates at the named node will be painted with the specified regime. - } - \item{which}{ - integer; - if \code{tree} is a \code{hansentree}, start not with a blank canvas but with the regime specifications \code{tree} contains for the character indicated by \code{which}. - } -} -\details{ - The names of \code{subtree} and \code{branch} must be the names of nodes of \code{tree}. - The painting proceeds in a particular order: one can overpaint a branch. - The subtrees indicated by the elements of \code{subtree} are painted first, in order. - Then the branches indicated by \code{branch} are painted. - If \code{tree} is a simple \code{ouchtree} object, then \code{paint} begins with a blank canvas, i.e., a tree painted with the single regime "nonspec". - If \code{tree} inherits class \code{hansentree}, then \code{paint} begins with the regimes specified in the \code{regimes} slot of \code{tree}. - Note that, if \code{tree} is a multivariate \code{hansentree}, then there are multiple regime specifications contained in \code{tree}. - In this case, the argument \code{which} lets you pick which one you wish to begin with; - by default, the first is used. +\item{tree}{An object of class \code{ouchtree}.} + +\item{subtree}{An optional named vector specifying the root nodes of subtrees. +Each branch that descends from this node will be painted with the +specified regime.} + +\item{branch}{An optional named vector specifying the end nodes of branches. +The unique branch that terminates at the named node will be painted with the specified regime.} + +\item{which}{integer; +if \code{tree} is a \code{hansentree}, start not with a blank canvas but with the regime specifications \code{tree} contains for the character indicated by \code{which}.} } \value{ - A vector of class 'factor' with names corresponding to the nodes in \code{tree}, specifying selective regimes. +A vector of class \sQuote{factor} with names corresponding to the nodes in \code{tree}, specifying selective regimes. +} +\description{ +Function to paint selective regimes on a phylogenetic tree. +} +\details{ +The names of \code{subtree} and \code{branch} must be the names of nodes of \code{tree}. +The painting proceeds in a particular order: +one can overpaint a branch. +The subtrees indicated by the elements of \code{subtree} are painted first, in order. +Then the branches indicated by \code{branch} are painted. +If \code{tree} is a simple \code{ouchtree} object, then \code{paint} begins with a blank canvas, +i.e., a tree painted with the single regime \dQuote{nonspec}. +If \code{tree} is of class \code{hansentree}, then \code{paint} begins with the regimes specified in the \code{regimes} slot of \code{tree}. +Note that, if \code{tree} is a multivariate \code{hansentree}, then there are multiple regime specifications contained in \code{tree}. +In this case, the argument \code{which} lets you pick which one you wish to begin with; +by default, the first is used. } -\author{Aaron A. King kingaa at umich dot edu} -\seealso{\code{ouchtree}, \code{hansen}} \examples{ -data(bimac) -x <- with(bimac,ouchtree(nodes=node,times=time/max(time),ancestors=ancestor,labels=species)) -r <- paint(x,subtree=c("1"="medium","9"="large","2"="small"),branch=c("38"="large","2"="medium")) +x <- with( + bimac, + ouchtree(nodes=node,times=time/max(time),ancestors=ancestor,labels=species) +) + +r <- paint(x,subtree=c("1"="medium","9"="large","2"="small"), + branch=c("38"="large","2"="medium")) plot(x,regimes=r,node.names=TRUE) -# compare to bimac['OU.LP'] -h5 <- hansen(data=log(bimac['size']),tree=x,regimes=bimac['OU.LP'],sqrt.alpha=1,sigma=1,reltol=1e-5) + +## compare to bimac['OU.LP'] +h5 <- hansen(data=log(bimac['size']),tree=x,regimes=bimac['OU.LP'], + sqrt.alpha=1,sigma=1,reltol=1e-5) r <- paint(h5,branch=c("18"="large"),subtree=c("9"="small")) plot(x,regimes=r,node.names=TRUE) } +\seealso{ +\code{ouchtree}, \code{hansen} +} +\author{ +Aaron A. King +} \keyword{models} diff --git a/man/plot.Rd b/man/plot.Rd new file mode 100644 index 0000000..0c8ca2f --- /dev/null +++ b/man/plot.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.R +\name{plot} +\alias{plot} +\alias{plot,ouchtree-method} +\alias{plot,hansentree-method} +\title{Plotting functions.} +\usage{ +\S4method{plot}{ouchtree}(x, y, regimes = NULL, node.names = FALSE, legend = TRUE, ..., labels) + +\S4method{plot}{hansentree}(x, y, regimes, ...) +} +\arguments{ +\item{x}{object to plot.} + +\item{y}{ignored.} + +\item{regimes}{factor or character; a vector of regime paintings.} + +\item{node.names}{logical; should node names be displayed?} + +\item{legend}{logical; display a legend?} + +\item{...}{additional arguments, passed to \code{text}.} + +\item{labels}{character; taxon labels.} +} +\description{ +Plot phylogenies. +} diff --git a/man/print.Rd b/man/print.Rd new file mode 100644 index 0000000..4e199eb --- /dev/null +++ b/man/print.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ouchtree.R, R/brown.R, R/hansen.R, R/print.R +\name{print,ouchtree-method} +\alias{print,ouchtree-method} +\alias{show,ouchtree-method} +\alias{show,browntree-method} +\alias{print,browntree-method} +\alias{print,hansentree-method} +\alias{show,hansentree-method} +\alias{print} +\title{Print methods} +\usage{ +\S4method{print}{ouchtree}(x, ...) + +\S4method{show}{ouchtree}(object) + +\S4method{show}{browntree}(object) + +\S4method{print}{browntree}(x, ...) + +\S4method{print}{hansentree}(x, ...) + +\S4method{show}{hansentree}(object) +} +\arguments{ +\item{...}{additional arguments, ignored.} + +\item{object, x}{object to display.} +} +\value{ +\code{print} displays the tree as a table, along with the coefficients of the fitted model and diagnostic information. + +\code{print} displays the tree as a table, along with the coefficients of the fitted model and diagnostic information. + +\code{print} displays the tree as a table, with (possibly) other information. +} +\description{ +Print methods +} diff --git a/man/simulate.Rd b/man/simulate.Rd new file mode 100644 index 0000000..ee7ae1b --- /dev/null +++ b/man/simulate.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brown.R, R/simulate.R, R/hansen.R +\docType{methods} +\name{simulate,browntree-method} +\alias{simulate,browntree-method} +\alias{simulate} +\alias{simulate,hansentree-method} +\title{Simulations of a phylogenetic trait model.} +\usage{ +\S4method{simulate}{browntree}(object, nsim = 1, seed = NULL, ...) + +\S4method{simulate}{hansentree}(object, nsim = 1, seed = NULL, ...) +} +\arguments{ +\item{object}{fitted model object to simulate.} + +\item{nsim}{integer; number of independent simulations.} + +\item{seed}{integer; if non-\code{NULL}, the RNG will be initialized with this seed for the simulations. +The RNG will be reset to its pre-existing state when \code{simulate} returns.} + +\item{...}{additional arguments, ignored.} +} +\value{ +\code{simulate} returns a list of data-frames, each comparable to the original data. +} +\description{ +\code{simulate} generates random deviates from a fitted model. +} diff --git a/man/summary.Rd b/man/summary.Rd new file mode 100644 index 0000000..52d18f0 --- /dev/null +++ b/man/summary.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brown.R, R/hansen.R, R/summary.R +\name{summary,browntree-method} +\alias{summary,browntree-method} +\alias{summary,hansentree-method} +\alias{summary} +\title{Summary methods} +\usage{ +\S4method{summary}{browntree}(object, ...) + +\S4method{summary}{hansentree}(object, ...) +} +\arguments{ +\item{object}{fitted model object.} + +\item{...}{additional arguments, ignored.} +} +\value{ +\code{summary} applied to a \code{browntree} object returns information about the fitted model, including parameter estimates and quantities describing the goodness of fit. + +\code{summary} applied to a \code{hansentree} method displays the estimated \eqn{\alpha}{alpha} and \eqn{\sigma^2}{sigma^2} matrices as well as various quantities describing the goodness of model fit. +} +\description{ +Summary methods +} diff --git a/man/update.Rd b/man/update.Rd new file mode 100644 index 0000000..0c4fc5b --- /dev/null +++ b/man/update.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brown.R, R/hansen.R, R/update.R +\docType{methods} +\name{update,browntree-method} +\alias{update,browntree-method} +\alias{update,hansentree-method} +\alias{update} +\title{Update and refit a model.} +\usage{ +\S4method{update}{browntree}(object, data, ...) + +\S4method{update}{hansentree}(object, data, regimes, sqrt.alpha, sigma, ...) +} +\arguments{ +\item{object}{fitted model object.} + +\item{data}{data that replace those used in the original fit.} + +\item{...}{Additional arguments replace the corresponding arguments in the original call.} + +\item{regimes}{A vector of codes, one for each node in the tree, specifying the selective regimes hypothesized to have been operative. +Corresponding to each node, enter the code of the regime hypothesized for the branch segment terminating in that node. +For the root node, because it has no branch segment terminating on it, the regime specification is irrelevant. +If there are \code{nchar} quantitative characters, then one can specify a single set of \code{regimes} for all characters or a list of \code{nchar} regime specifications, one for each character.} + +\item{sqrt.alpha}{These are used to initialize the optimization algorithm. +The selection strength matrix \eqn{\alpha}{alpha} and the random drift variance-covariance matrix \eqn{\sigma^2}{sigma^2} are parameterized by their matrix square roots. +Specifically, these initial guesses are each packed into lower-triangular matrices (column by column). +The product of this matrix with its transpose is the \eqn{\alpha}{alpha} or \eqn{\sigma^2}{sigma^2} matrix. +See Details.} + +\item{sigma}{These are used to initialize the optimization algorithm. +The selection strength matrix \eqn{\alpha}{alpha} and the random drift variance-covariance matrix \eqn{\sigma^2}{sigma^2} are parameterized by their matrix square roots. +Specifically, these initial guesses are each packed into lower-triangular matrices (column by column). +The product of this matrix with its transpose is the \eqn{\alpha}{alpha} or \eqn{\sigma^2}{sigma^2} matrix. +See Details.} +} +\value{ + + +\code{update} returns a new fitted-model object of the same class as \code{object}. +} +\description{ +\code{update} will update a model and re-fit. +This allows one to change the data and/or parameters. +} diff --git a/tests/boot.R b/tests/boot.R index 593e47f..be08c8c 100644 --- a/tests/boot.R +++ b/tests/boot.R @@ -1,5 +1,4 @@ library(ouch) -data(bimac) tree <- with(bimac,ouchtree(node,ancestor,time/max(time),species)) h1 <- brown( data=log(bimac['size']), diff --git a/tests/boot.Rout.save b/tests/boot.Rout.save index 4eeda5c..99146cc 100644 --- a/tests/boot.Rout.save +++ b/tests/boot.Rout.save @@ -1,6 +1,6 @@ -R version 3.6.0 (2019-04-26) -- "Planting of a Tree" -Copyright (C) 2019 The R Foundation for Statistical Computing +R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -18,8 +18,6 @@ Type 'demo()' for some demos, 'help()' for on-line help, or Type 'q()' to quit R. > library(ouch) -Loading required package: subplex -> data(bimac) > tree <- with(bimac,ouchtree(node,ancestor,time/max(time),species)) > h1 <- brown( + data=log(bimac['size']), @@ -42,7 +40,7 @@ Loading required package: subplex > b5 <- sapply(simdat,function(x)summary(update(h5,data=x))$aic.c) > toc <- Sys.time() > print(toc-tic) -Time difference of 6.676345 secs +Time difference of 7.572737 secs > cat("approximate 95% AIC.c cutoff",signif(quantile(b1-b5,0.95),digits=3),"\n") approximate 95% AIC.c cutoff 2.58 > diff --git a/tests/exacttree.Rout.save b/tests/exacttree.Rout.save index 47714a1..b5a9915 100644 --- a/tests/exacttree.Rout.save +++ b/tests/exacttree.Rout.save @@ -1,6 +1,6 @@ -R version 3.6.0 (2019-04-26) -- "Planting of a Tree" -Copyright (C) 2019 The R Foundation for Statistical Computing +R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -18,7 +18,6 @@ Type 'demo()' for some demos, 'help()' for on-line help, or Type 'q()' to quit R. > library(ouch) -Loading required package: subplex > > treedat <- data.frame(node=1:7,anc=NA,time=NA) > treedat$anc[-1] <- c(1,2,2,1,5,5) diff --git a/tests/simtree.Rout.save b/tests/simtree.Rout.save index f2be2ee..89e996f 100644 --- a/tests/simtree.Rout.save +++ b/tests/simtree.Rout.save @@ -1,6 +1,6 @@ -R version 3.6.0 (2019-04-26) -- "Planting of a Tree" -Copyright (C) 2019 The R Foundation for Statistical Computing +R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -18,7 +18,6 @@ Type 'demo()' for some demos, 'help()' for on-line help, or Type 'q()' to quit R. > library(ouch) -Loading required package: subplex > > oubranch <- function (x0, t0, t1, alpha, sigma, theta, nstep = 10) { + x <- x0 @@ -93,7 +92,7 @@ Loading required package: subplex > tree <- with(x,ouchtree(node,ancestor,time)) > toc <- Sys.time() > print(toc-tic) -Time difference of 1.985994 secs +Time difference of 1.638595 secs > > bfit <- brown(data=x[c("A","B")],tree) > print(summary(bfit)) @@ -136,7 +135,7 @@ $dof > hfit <- hansen(data=x[c("A","B")],tree=tree,regimes=x['reg'],sqrt.alpha=a,sigma=s,fit=F) > toc <- Sys.time() > print(toc-tic) -Time difference of 0.1315975 secs +Time difference of 0.1367433 secs > print(summary(hfit)) $call hansen(data = x[c("A", "B")], tree = tree, regimes = x["reg"], @@ -204,7 +203,7 @@ $dof > tree <- with(x,ouchtree(node,ancestor,time)) > toc <- Sys.time() > print(toc-tic) -Time difference of 1.880172 secs +Time difference of 1.591859 secs > > bfit <- brown(data=x[c("A","B")],tree) > print(summary(bfit)) @@ -284,7 +283,7 @@ $dof > hfit <- hansen(data=x[c("A","B")],tree=tree,regimes=x['reg'],sqrt.alpha=a,sigma=s,fit=F) > toc <- Sys.time() > print(toc-tic) -Time difference of 0.2379658 secs +Time difference of 0.1803641 secs > print(summary(hfit)) $call hansen(data = x[c("A", "B")], tree = tree, regimes = x["reg"],