diff -Nru r-cran-phytools-0.7-47/data/datalist r-cran-phytools-0.7-70/data/datalist --- r-cran-phytools-0.7-47/data/datalist 2020-06-01 05:35:45.000000000 +0000 +++ r-cran-phytools-0.7-70/data/datalist 2020-09-19 14:37:45.000000000 +0000 @@ -1,9 +1,13 @@ anole.data anoletree +flatworm.data +flatworm.tree mammal.data mammal.tree salamanders sunfish.data sunfish.tree +vertebrate.data +vertebrate.tree wasp.data wasp.trees Binary files /tmp/tmpq4s5d8/HczhHbu3Zm/r-cran-phytools-0.7-47/data/flatworm.data.rda and /tmp/tmpq4s5d8/Fo10e83ML8/r-cran-phytools-0.7-70/data/flatworm.data.rda differ Binary files /tmp/tmpq4s5d8/HczhHbu3Zm/r-cran-phytools-0.7-47/data/flatworm.tree.rda and /tmp/tmpq4s5d8/Fo10e83ML8/r-cran-phytools-0.7-70/data/flatworm.tree.rda differ Binary files /tmp/tmpq4s5d8/HczhHbu3Zm/r-cran-phytools-0.7-47/data/mammal.data.rda and /tmp/tmpq4s5d8/Fo10e83ML8/r-cran-phytools-0.7-70/data/mammal.data.rda differ Binary files /tmp/tmpq4s5d8/HczhHbu3Zm/r-cran-phytools-0.7-47/data/mammal.tree.rda and /tmp/tmpq4s5d8/Fo10e83ML8/r-cran-phytools-0.7-70/data/mammal.tree.rda differ Binary files /tmp/tmpq4s5d8/HczhHbu3Zm/r-cran-phytools-0.7-47/data/vertebrate.data.rda and /tmp/tmpq4s5d8/Fo10e83ML8/r-cran-phytools-0.7-70/data/vertebrate.data.rda differ Binary files /tmp/tmpq4s5d8/HczhHbu3Zm/r-cran-phytools-0.7-47/data/vertebrate.tree.rda and /tmp/tmpq4s5d8/Fo10e83ML8/r-cran-phytools-0.7-70/data/vertebrate.tree.rda differ diff -Nru r-cran-phytools-0.7-47/debian/changelog r-cran-phytools-0.7-70/debian/changelog --- r-cran-phytools-0.7-47/debian/changelog 2020-07-02 12:50:25.000000000 +0000 +++ r-cran-phytools-0.7-70/debian/changelog 2020-09-21 09:18:07.000000000 +0000 @@ -1,3 +1,9 @@ +r-cran-phytools (0.7-70-1) unstable; urgency=medium + + * New upstream version + + -- Andreas Tille Mon, 21 Sep 2020 11:18:07 +0200 + r-cran-phytools (0.7-47-1) unstable; urgency=medium * Team upload. diff -Nru r-cran-phytools-0.7-47/debian/control r-cran-phytools-0.7-70/debian/control --- r-cran-phytools-0.7-47/debian/control 2020-07-02 12:50:24.000000000 +0000 +++ r-cran-phytools-0.7-70/debian/control 2020-09-21 09:18:07.000000000 +0000 @@ -9,7 +9,6 @@ r-base-dev, r-cran-ape (>= 4.0), r-cran-maps, - r-cran-animation, r-cran-clustergeneration, r-cran-coda, r-cran-combinat, diff -Nru r-cran-phytools-0.7-47/DESCRIPTION r-cran-phytools-0.7-70/DESCRIPTION --- r-cran-phytools-0.7-47/DESCRIPTION 2020-06-01 21:00:02.000000000 +0000 +++ r-cran-phytools-0.7-70/DESCRIPTION 2020-09-19 22:20:02.000000000 +0000 @@ -1,19 +1,19 @@ Package: phytools -Version: 0.7-47 -Date: 2020-6-01 +Version: 0.7-70 +Date: 2020-9-19 Title: Phylogenetic Tools for Comparative Biology (and Other Things) Author: Liam J. Revell Maintainer: Liam J. Revell Depends: R (>= 3.5.0), ape (>= 4.0), maps -Imports: animation, clusterGeneration, coda, combinat, expm, graphics, - grDevices, gtools, MASS, methods, mnormt, nlme, numDeriv, - phangorn (>= 2.3.1), plotrix, scatterplot3d, stats, utils -Suggests: geiger, RColorBrewer, rgl +Imports: clusterGeneration, coda, combinat, expm, graphics, grDevices, + gtools, MASS, methods, mnormt, nlme, numDeriv, phangorn (>= + 2.3.1), plotrix, scatterplot3d, stats, utils +Suggests: animation, geiger, RColorBrewer, rgl ZipData: no Description: A wide range of functions for phylogenetic analysis. Functionality is concentrated in phylogenetic comparative biology, but also includes numerous methods for visualizing, manipulating, reading or writing, and even inferring phylogenetic trees and data. Included among the functions in phylogenetic comparative biology are various for ancestral state reconstruction, model-fitting, simulation of phylogenies and data, and multivariate analysis. There are a broad range of plotting methods for phylogenies and comparative data which include, but are not restricted to, methods for mapping trait evolution on trees, for projecting trees into phenotypic space or a geographic map, and for visualizing correlated speciation between trees. Finally, there are a number of functions for reading, writing, analyzing, inferring, simulating, and manipulating phylogenetic trees and comparative data not covered by other packages. For instance, there are functions for randomly or non-randomly attaching species or clades to a phylogeny, for estimating supertrees or consensus phylogenies from a set, for simulating trees and phylogenetic data under a range of models, and for a wide variety of other manipulations and analyses that phylogenetic biologists might find useful in their research. License: GPL (>= 2) -URL: http://github.com/liamrevell/phytools -Packaged: 2020-06-01 20:18:09 UTC; liamj +URL: https://github.com/liamrevell/phytools +Packaged: 2020-09-19 18:15:53 UTC; liamj Repository: CRAN -Date/Publication: 2020-06-01 21:00:02 UTC +Date/Publication: 2020-09-19 22:20:02 UTC NeedsCompilation: no diff -Nru r-cran-phytools-0.7-47/man/anoletree.Rd r-cran-phytools-0.7-70/man/anoletree.Rd --- r-cran-phytools-0.7-47/man/anoletree.Rd 2020-06-01 16:23:43.000000000 +0000 +++ r-cran-phytools-0.7-70/man/anoletree.Rd 2020-09-19 16:53:51.000000000 +0000 @@ -1,49 +1,69 @@ \name{anoletree} \alias{anole.data} \alias{anoletree} +\alias{flatworm.tree} +\alias{flatworm.data} \alias{mammal.tree} \alias{mammal.data} \alias{salamanders} \alias{sunfish.tree} \alias{sunfish.data} +\alias{vertebrate.tree} +\alias{vertebrate.data} \alias{wasp.trees} \alias{wasp.data} \title{Phylogenetic datasets} \description{ \code{anoletree} is a phylogeny of Greater Antillean anole species with a mapped discrete character - \emph{ecomorph class}. \code{anole.data} is a data frame of morphological characters. Data and tree are from Mahler et al. (2010). + \code{flatworm.tree} and \code{flatworm.data} are a phylogeny and dataset of habitat preferences for flatworms from Benitez-Alvarez et al. (2000). \code{flatworm.tree} has been made ultrametric using penalized likelihood. + \code{mammal.tree} and \code{mammal.data} are the phylogeny and dataset for mammal body size and home range size from Garland et al. (1992). \code{salamanders} is a phylogeny of \emph{Plethodon} salamanders from Highton and Larson (1979). According to Wikipedia, the genus \emph{Plethodon} contains 55 species in total. \code{sunfish.tree} and \code{sunfish.data} are the phylogeny and dataset for Centrarchidae and buccal morphology (respectively) from Revell and Collar (2009). + \code{vertebrate.tree} is a time-calibrated phylogeny of vertebrates and \code{vertebrate.data} is a dataset of phenotypic traits. The phylogeny is from \url{http://www.timetree.org/} (Hedges et al. 2006). + \code{wasp.trees} and \code{wasp.data} are the phylogeny and host-parasite associations from Lopez-Vaamonde et al. (2001). } \usage{ data(anole.data) data(anoletree) +data(flatworm.tree) +data(flatworm.data) data(mammal.data) data(mammal.tree) data(salamanders) data(sunfish.data) data(sunfish.tree) +data(vertebrate.tree) +data(vertebrate.data) data(wasp.data) data(wasp.trees) } \format{ \code{anoletree} is an object of class \code{"simmap"}. \code{anole.data} is a data frame. + \code{flatworm.tree} is an object of class \code{"phylo"}. \code{flatworm.data} is a data frame. + \code{mammal.tree} is an object of class \code{"phylo"}. \code{mammal.data} is a data frame. \code{salamanders} is an object of class \code{"phylo"}. \code{sunfish.tree} is an object of class \code{"simmap"}. \code{sunfish.data} is a data frame. + \code{vertebrate.tree} is an object of class \code{"phylo"}. \code{vertebrate.data} is a data frame. + \code{wasp.trees} is an object of class \code{"multiPhylo"}. \code{wasp.data} is a data frame. } \source{ + Benitez-Alvarez, L., A. Maria Leal-Zanchet, A. Oceguera-Figueroa, R. Lopes Ferreira, D. de Medeiros Bento, J. Braccini, R. Sluys, and M. Riutort. (2020) Phylogeny and biogeography of the Cavernicola (Platyhelminthes: Tricladida): Relicts of an epigean group sheltering in caves? \emph{Molecular Phylogenetics and Evolution}, \bold{145}, 106709. + Garland, T., Jr., P. H. Harvey, and A. R. Ives. (1992) Procedures for the analysis of comparative data using phylogenetically independent contrasts. \emph{Systematic Biology}, \bold{41}, 18-32. + + Hedges, S. B., J. Dudley, and S. Kumar. (2006) TimeTree: A public knowledgebase of divergence times among organisms. \emph{Bioinformatics}, \bold{22}, 2971-2972. Highton, R., and A. Larson. (1979) The genetic relationships of the salamanders of the genus \emph{Plethodon}. \emph{Systematic Zoology}, \bold{28}, 579-599. diff -Nru r-cran-phytools-0.7-47/man/as.Qmatrix.Rd r-cran-phytools-0.7-70/man/as.Qmatrix.Rd --- r-cran-phytools-0.7-47/man/as.Qmatrix.Rd 2020-05-31 16:39:00.000000000 +0000 +++ r-cran-phytools-0.7-70/man/as.Qmatrix.Rd 2020-06-03 16:46:32.000000000 +0000 @@ -15,10 +15,12 @@ \item{...}{optional arguments.} } \description{ - This function extracts a Q-matrix (in the form of an object of class \code{"Qmatrix"} from a fitted M\emph{k} model. + This function extracts a \bold{Q}-matrix (in the form of an object of class \code{"Qmatrix"} from a fitted M\emph{k} model. } \value{ An object of class \code{"Qmatrix"}. + + \code{plot.Qmatrix} invisibly returns the coordinates of vertices of the plotted \bold{Q}-matrix. } \references{ Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). \emph{Methods Ecol. Evol.}, \bold{3}, 217-223. diff -Nru r-cran-phytools-0.7-47/man/drop.tip.contMap.Rd r-cran-phytools-0.7-70/man/drop.tip.contMap.Rd --- r-cran-phytools-0.7-47/man/drop.tip.contMap.Rd 2015-02-04 14:23:26.000000000 +0000 +++ r-cran-phytools-0.7-70/man/drop.tip.contMap.Rd 2020-08-21 01:35:49.000000000 +0000 @@ -1,14 +1,16 @@ \name{drop.tip.contMap} \alias{drop.tip.contMap} \alias{drop.tip.densityMap} +\alias{keep.tip.contMap} \title{Drop tip or tips from an object of class "contMap" or "densityMap"} \usage{ drop.tip.contMap(x, tip) drop.tip.densityMap(x, tip) +keep.tip.contMap(x, tip) } \arguments{ \item{x}{an object of class \code{"contMap"} or \code{"densityMap"}.} - \item{tip}{name or names of species to be dropped.} + \item{tip}{name or names of species to be dropped or kept.} } \description{ This function drops one or multiple tips from an object of class \code{"contMap"} or \code{"densityMap"}. This function is equivalent to \code{\link{drop.tip}} but for an object of this class. @@ -24,7 +26,7 @@ } \author{Liam Revell \email{liam.revell@umb.edu}} \seealso{ - \code{\link{contMap}}, \code{\link{densityMap}}, \code{\link{drop.tip}}, \code{\link{drop.tip.simmap}} + \code{\link{contMap}}, \code{\link{densityMap}}, \code{\link{drop.tip}}, \code{\link{drop.tip.simmap}}, \code{\link{keep.tip}} } \keyword{phylogenetics} \keyword{utilities} diff -Nru r-cran-phytools-0.7-47/man/drop.tip.multiPhylo.Rd r-cran-phytools-0.7-70/man/drop.tip.multiPhylo.Rd --- r-cran-phytools-0.7-47/man/drop.tip.multiPhylo.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-phytools-0.7-70/man/drop.tip.multiPhylo.Rd 2020-06-26 21:16:43.000000000 +0000 @@ -0,0 +1,31 @@ +\name{drop.tip.multiPhylo} +\alias{drop.tip.multiPhylo} +\alias{drop.tip.multiSimmap} +\title{Drop tip or tips from an object of class "multiPhylo" or "multiSimmap"} +\usage{ +drop.tip.multiPhylo(phy, tip, ...) +drop.tip.multiSimmap(phy, tip) +} +\arguments{ + \item{phy}{an object of class \code{"multiPhylo"} or \code{"multiSimmap"}.} + \item{tip}{name or names of species to be dropped.} + \item{...}{optional arguments to be passed to \code{\link{drop.tip}}. Most optional arguments work, with the exception of \code{interactive=TRUE} which will return an error.} +} +\description{ + This function drops one or multiple tips from all the trees of an object of class \code{"multiPhylo"} or \code{"multiSimmap"}. +} +\details{ + This function merely wraps \code{\link{drop.tip}} and \code{\link{drop.tip.simmap}}. Note that \code{drop.tip.multiSimmap} is merely just an alias of \code{drop.tip.multiPhylo}. +} +\value{ + An object of class \code{"multiPhylo"} or \code{"multiSimmap"}, depending on the input object class. +} +\references{ + Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). \emph{Methods Ecol. Evol.}, \bold{3}, 217-223. +} +\author{Liam Revell \email{liam.revell@umb.edu}} +\seealso{ + \code{\link{drop.tip}}, \code{\link{drop.tip.simmap}} +} +\keyword{phylogenetics} +\keyword{utilities} diff -Nru r-cran-phytools-0.7-47/man/edge.widthMap.Rd r-cran-phytools-0.7-70/man/edge.widthMap.Rd --- r-cran-phytools-0.7-47/man/edge.widthMap.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-phytools-0.7-70/man/edge.widthMap.Rd 2020-07-20 18:09:46.000000000 +0000 @@ -0,0 +1,50 @@ +\name{edge.widthMap} +\alias{edge.widthMap} +\alias{plot.edge.widthMap} +\title{Map continuous trait evolution on the tree} +\usage{ +edge.widthMap(tree, x, ...) +\method{plot}{edge.widthMap}(x, max.width=0.9, legend="trait value", ...) +} +\arguments{ + \item{tree}{object of class \code{"phylo"}.} + \item{x}{a numerical vector of phenotypic trait values for species. \code{names(x)} should contain the species names and match \code{tree$tip.label}. Or, for \code{plot.edge.widthMap}, an object of class \code{"edge.widthMap"}.} + \item{max.width}{maximum edge width in plot units.} + \item{legend}{label for the plot legend.} + \item{...}{optional arguments - especially for the \code{plot} method. Perhaps the most important of these is \code{min.width}, which defaults to \code{0} but could probably be increased for many datasets and graphical devices. Other arguments are passed internally to \code{\link{plotTree}}.} +} +\description{ + Function maps a discrete character onto the edges of the tree using variable edge widths. +} +\value{ + \code{edge.widthMap} returns an object of class \code{"edge.widthMap"}. + + \code{plot.edge.widthMap} can be used to plot this object. +} +\references{ + Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). \emph{Methods Ecol. Evol.}, \bold{3}, 217-223. +} +\author{Liam Revell \email{liam.revell@umb.edu}} +\seealso{ + \code{\link{contMap}}, \code{\link{fastAnc}} +} +\examples{ +## load data from Garland et al. (1992) +data(mammal.tree) +data(mammal.data) +## extract character of interest +ln.bodyMass<-log(setNames(mammal.data$bodyMass, + rownames(mammal.data))) +## create "edge.widthMap" object +mammal.ewMap<-edge.widthMap(mammal.tree,ln.bodyMass, + min.width=0.05) +## plot it +plot(mammal.ewMap,legend="log(body mass)") +par(mar=c(5.1,4.1,4.1,2.1)) ## reset margins to default +} +\keyword{ancestral states} +\keyword{phylogenetics} +\keyword{plotting} +\keyword{comparative method} +\keyword{continuous character} +\keyword{maximum likelihood} diff -Nru r-cran-phytools-0.7-47/man/fitMk.Rd r-cran-phytools-0.7-70/man/fitMk.Rd --- r-cran-phytools-0.7-47/man/fitMk.Rd 2020-06-01 02:16:24.000000000 +0000 +++ r-cran-phytools-0.7-70/man/fitMk.Rd 2020-09-19 18:14:54.000000000 +0000 @@ -4,11 +4,14 @@ \alias{plot.gfit} \alias{fitmultiMk} \alias{fitpolyMk} +\alias{graph.polyMk} \alias{plot.fitpolyMk} \alias{mcmcMk} \alias{plot.mcmcMk} \alias{density.mcmcMk} \alias{plot.density.mcmcMk} +\alias{fitHRM} +\alias{plot.fitHRM} \title{Fits Mk model} \usage{ fitMk(tree, x, model="SYM", fixedQ=NULL, ...) @@ -16,11 +19,14 @@ \method{plot}{gfit}(x, ...) fitmultiMk(tree, x, model="ER", ...) fitpolyMk(tree, x, model="SYM", ordered=FALSE, ...) +graph.polyMk(k=2, model="SYM", ordered=FALSE, ...) \method{plot}{fitpolyMk}(x, ...) mcmcMk(tree, x, model="ER", ngen=10000, ...) \method{plot}{mcmcMk}(x, ...) \method{density}{mcmcMk}(x, ...) \method{plot}{density.mcmcMk}(x, ...) +fitHRM(tree, x, model="ARD", ncat=2, ...) +\method{plot}{fitHRM}(x, ...) } \arguments{ \item{tree}{an object of class \code{"phylo"}. In the case of \code{fitmultiMk} an object of class \code{"simmap"} with a mapped discrete character.} @@ -28,8 +34,10 @@ \item{model}{model. See \code{make.simmap} or \code{ace} for details.} \item{fixedQ}{fixed value of transition matrix \code{Q}, if one is desired.} \item{ordered}{for \code{fitpolyMk}, a logical value indicating whether or not the character should be treated as ordered. For now the function assumes alphanumerical order (i.e., numbers sorted by their initial and then successive digits followed by characters or character strings in alphabetical order).} + \item{k}{For \code{graph.polyMk}, the number of monomorphic states for the discrete trait.} \item{ngen}{number of generations of MCMC for \code{mcmcMk}.} - \item{...}{optional arguments, including \code{pi}, the prior distribution at the root node (defaults to \code{pi="equal"}). For \code{plot} method optional arguments include (but may not be limited to): \code{signif}, the number of digits for the rates to be plotted; \code{main}, a character vector of length two with the headings for each subplot; \code{cex.main}, \code{cex.traits}, and \code{cex.rates}, font sizes for the various text elements of the plot; and \code{show.zeros}, a logical argument specifying whether or not to plot arrows with the ML estimated transition rate is not different from zero (with tolerance specified by the optional argument \code{tol}).} + \item{ncat}{number of rate categories (per level of the discrete trait) in the hidden-rate model.} + \item{...}{optional arguments, including \code{pi}, the prior distribution at the root node (defaults to \code{pi="equal"}). Other options for \code{pi} include \code{pi="fitzjohn"} (which implements the prior distribution of Fitzjohn et al. 2009), \code{pi="estimated"} (which finds the stationary distribution of state frequencies and sets that as the prior), or an arbitrary prior distribution specified by the user. For \code{plot} method optional arguments include (but may not be limited to): \code{signif}, the number of digits for the rates to be plotted; \code{main}, a character vector of length two with the headings for each subplot; \code{cex.main}, \code{cex.traits}, and \code{cex.rates}, font sizes for the various text elements of the plot; and \code{show.zeros}, a logical argument specifying whether or not to plot arrows with the ML estimated transition rate is not different from zero (with tolerance specified by the optional argument \code{tol}). Finally, for \code{fitpolyMk}, \code{max.poly} can be set for the \code{ordered=TRUE} model. \code{max.poly} defaults to the highest level of polymorphism observed in the data.} } \description{ The function \code{fitMk} fits a so-called M\emph{k} model for discrete character evolution. @@ -40,15 +48,19 @@ The function \code{fitpolyMk} fits an M\emph{k} model to data for a discrete character with intraspecific polymorphism. Polymorphic species should be coded with the name of the two or more states recorded for the species separated by a space (e.g., \code{A+B} would indicate that both states \code{A} and \code{B} are found in the corresponding taxon). Invariably it's assumed that transitions between states must occur through a polymorphic condition, whereas transitions \emph{cannot} occur directly between two incompatible polymorphic conditions. For instance, a transition between \code{A+B} and \code{B+C} would have to occur through the monomorphic state \code{B}. At time of writing, this function permits the models \code{"ER"} (equal rates for all permitted transitions), \code{"SYM"} (symmetric backward & forward rates for all permitted transitions), \code{"ARD"} (all-rates-different for permitted transitions), and a new model called \code{"transient"} in which the acquisition of polymorphism (e.g., \code{A -> A+B}) is assumed to occur at a different rate than its loss (e.g., \code{A+B -> B}). The method \code{plot.fitpolyMk} plots the fitted M\emph{k} model with intraspecific polymorphism. - Finally, the function \code{mcmcMk} runs a Bayesian MCMC version of \code{fitMk}. The shape of the prior distribution of the transition rates is \eqn{\Gamma}, with \eqn{\alpha} and \eqn{\beta} via the argument \code{prior}, which takes the form of a list. The default value of \eqn{\alpha} is 0.1, and \eqn{\beta} defaults to a value such tha \eqn{\alpha/\beta} is equal to the parsimony score for \code{x} divided by the sum of the edge lengths of the tree. The shape of the proposal distribution is normal, with mean zero and a variance that can be controlled by the user via the optional argument \code{prior.var}. The argument \code{auto.tune}, if \code{TRUE} or \code{FALSE}, indicates whether or not to 'tune' the proposal variance up or down to target a particular acceptance rate (defaults to 0.5). \code{auto.tune} can also be a numeric value between 0 and 1, in which case this value will be the target acceptance ratio. The argument \code{plot} indicates whether the progress of the MCMC should be plotted (defaults to \code{TRUE}, but runs much faster when set to \code{FALSE}). + The function \code{mcmcMk} runs a Bayesian MCMC version of \code{fitMk}. The shape of the prior distribution of the transition rates is \eqn{\Gamma}, with \eqn{\alpha} and \eqn{\beta} via the argument \code{prior}, which takes the form of a list. The default value of \eqn{\alpha} is 0.1, and \eqn{\beta} defaults to a value such tha \eqn{\alpha/\beta} is equal to the parsimony score for \code{x} divided by the sum of the edge lengths of the tree. The shape of the proposal distribution is normal, with mean zero and a variance that can be controlled by the user via the optional argument \code{prior.var}. The argument \code{auto.tune}, if \code{TRUE} or \code{FALSE}, indicates whether or not to 'tune' the proposal variance up or down to target a particular acceptance rate (defaults to 0.5). \code{auto.tune} can also be a numeric value between 0 and 1, in which case this value will be the target acceptance ratio. The argument \code{plot} indicates whether the progress of the MCMC should be plotted (defaults to \code{TRUE}, but runs much faster when set to \code{FALSE}). The method \code{plot.mcmcMk} plots a log-likelihood trace and a trace of the rate parameters from the MCMC. (This the samem graph that is created by setting \code{plot=TRUE} in \code{mcmcMk}.) The method \code{density.mcmcMk} computes a posterior density on the transition rates in the model from the posterior sample obtained in the MCMC, will import the package \emph{coda} if it is available, and returns an object of class \code{"density.mcmcMk"}. Finally, the method \code{plot.density.mcmcMk} creates a plot of the posterior density (or a set of plots) for the transition rates between states. + + Finally, the function \code{fitHRM} fits a hidden-rate M\emph{k} model following Beaulieu et al. (2013). For the hidden-rate model we need to specify a number of rate categories for each level of the trait - and this can be a vector of different values for each trait. We can also choose a model (\code{"ER"}, \code{"SYM"}, or \code{"ARD"}), as well as whether or not to treat the character as a 'threshold' trait (\code{umbral=TRUE}, defaults to \code{FALSE}). This latter model is basically one that allows absorbing conditions for some hidden states. This is a work in progress. } \details{ - Note that both \code{fitMk} and \code{fitmultiMk} recycle code from \code{\link{ace}} in the \emph{ape} package for computing the likelihood. \code{fitpolyMk} and \code{mcmcMk} use \code{fitMk} internally to compute the likelihood. + Note that both \code{fitMk} and \code{fitmultiMk} recycle code from \code{\link{ace}} in the \emph{ape} package for computing the likelihood. \code{fitpolyMk}, \code{mcmcMk}, and \code{fitHRM} use \code{fitMk} internally to compute the likelihood. } \value{ - An object of class \code{"fitMk"}, \code{"fitmultiMk"}, \code{"fitpolyMk"}, or \code{"mcmcMk"}. In the case of \code{density.mcmcMk} an object of class \code{"density.mcmcMk"}. + An object of class \code{"fitMk"}, \code{"fitmultiMk"}, \code{"fitpolyMk"}, \code{"mcmcMk"}, or \code{"fitHRM"}. In the case of \code{density.mcmcMk} an object of class \code{"density.mcmcMk"}. + + \code{plot.fitMk}, \code{plot.gfit}, and \code{plot.HRM} invisibly return the coordinates of vertices of the plotted \bold{Q}-matrix. } \references{ Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). \emph{Methods Ecol. Evol.}, \bold{3}, 217-223. @@ -72,6 +84,31 @@ print(fit.ARD) ## compare the models AIC(fit.ER,fit.ARD) + +## load tree and data from Benitez-Alvarez et al. (2000) +data(flatworm.data) +data(flatworm.tree) +## extract discrete character (habitat) +habitat<-setNames(flatworm.data$Habitat, + rownames(flatworm.data)) +## fit polymorphic models "ER" and "transient" +fitpoly.ER<-fitpolyMk(flatworm.tree,habitat, + model="ER") +fitpoly.transient<-fitpolyMk(flatworm.tree,habitat, + model="transient") +## print fitted models +print(fitpoly.ER) +print(fitpoly.transient) +## compare model +AIC(fitpoly.ER,fitpoly.transient) +## plot models +par(mfrow=c(2,1)) +plot(fitpoly.ER) +mtext("a) ER polymorphic model",adj=0,line=1) +plot(fitpoly.transient) +mtext("b) Transient polymorphic model",adj=0, + line=1) +par(mfrow=c(1,1)) } \keyword{phylogenetics} \keyword{comparative method} diff -Nru r-cran-phytools-0.7-47/man/phylomorphospace.Rd r-cran-phytools-0.7-70/man/phylomorphospace.Rd --- r-cran-phytools-0.7-47/man/phylomorphospace.Rd 2020-05-31 21:28:38.000000000 +0000 +++ r-cran-phytools-0.7-70/man/phylomorphospace.Rd 2020-09-19 14:09:39.000000000 +0000 @@ -47,7 +47,7 @@ node.size=c(0,1.2),xlab="relative buccal length", ylab="relative gape width") title(main="Phylomorphospace of buccal morphology in Centrarchidae", - font.main=3) + font.main=3) } \keyword{ancestral states} \keyword{phylogenetics} diff -Nru r-cran-phytools-0.7-47/man/phytools-package.Rd r-cran-phytools-0.7-70/man/phytools-package.Rd --- r-cran-phytools-0.7-47/man/phytools-package.Rd 2020-06-01 15:12:48.000000000 +0000 +++ r-cran-phytools-0.7-70/man/phytools-package.Rd 2020-09-19 16:41:36.000000000 +0000 @@ -10,7 +10,7 @@ The complete list of functions can be displayed with \code{library(help=phytools)}. - The \emph{phytools} development page is \url{http://github.com/liamrevell/phytools/}. More information on \emph{phytools} can also be found at \url{http://blog.phytools.org} or \url{http://www.phytools.org}. + The \emph{phytools} development page is \url{https://github.com/liamrevell/phytools/}. More information on \emph{phytools} can also be found at \url{http://blog.phytools.org} or \url{http://www.phytools.org}. If you use \emph{phytools} (or other packages that depend on \emph{phytools}) in a publication, please \emph{cite it}. The appropriate citation for \emph{phytools} is listed below or can be obtained using \code{citation("phytools")} with the package installed. } diff -Nru r-cran-phytools-0.7-47/man/plotTree.Rd r-cran-phytools-0.7-70/man/plotTree.Rd --- r-cran-phytools-0.7-47/man/plotTree.Rd 2015-02-04 14:52:07.000000000 +0000 +++ r-cran-phytools-0.7-70/man/plotTree.Rd 2020-09-19 14:30:00.000000000 +0000 @@ -22,8 +22,10 @@ \code{\link{plot.phylo}}, \code{\link{plotSimmap}} } \examples{ -tree<-pbtree(n=25) -plotTree(tree,color="blue",ftype="i") +data(vertebrate.tree) +plotTree(vertebrate.tree,fsize=1.2,ftype="i") +## reset margins +par(mar=c(5.1,4.1,4.1,2.1)) } \keyword{phylogenetics} \keyword{plotting} diff -Nru r-cran-phytools-0.7-47/man/plotTree.wBars.Rd r-cran-phytools-0.7-70/man/plotTree.wBars.Rd --- r-cran-phytools-0.7-47/man/plotTree.wBars.Rd 2020-06-01 05:15:54.000000000 +0000 +++ r-cran-phytools-0.7-70/man/plotTree.wBars.Rd 2020-09-19 14:29:47.000000000 +0000 @@ -56,6 +56,21 @@ plotTree.barplot(anoletree,exp(svl), args.plotTree=list(fsize=0.5), args.barplot=list(xlab="SVL (mm)")) + +## load vertebrate tree and data +data(vertebrate.tree) +data(vertebrate.data) +## plotTree.barplot +options(scipen=3) ## change sci-notation +plotTree.barplot(vertebrate.tree, + setNames(vertebrate.data$Mass, + rownames(vertebrate.data)), + args.barplot=list(log="x", + xlab="mass (kg)", + xlim=c(0.01,500000), + col=palette()[4])) +options(scipen=0) + ## reset par to defaults par(mfrow=c(1,1),mar=c(5.1,4.1,4.1,2.1)) } diff -Nru r-cran-phytools-0.7-47/man/strahlerNumber.Rd r-cran-phytools-0.7-70/man/strahlerNumber.Rd --- r-cran-phytools-0.7-47/man/strahlerNumber.Rd 2015-02-04 14:37:32.000000000 +0000 +++ r-cran-phytools-0.7-70/man/strahlerNumber.Rd 2020-09-19 16:43:41.000000000 +0000 @@ -12,7 +12,7 @@ \item{plot}{logical value indicating whether to plot the tree with Strahler numbers for node labels.} } \description{ - The function \code{strahlerNumber} computes the Strahler number of all nodes and tips in the tree. For more information about Strahler numbers see \url{http://en.wikipedia.org/wiki/Strahler_number}. The function \code{extract.strahlerNumber} extracts all of the most inclusive clades of Strahler number \code{i}. + The function \code{strahlerNumber} computes the Strahler number of all nodes and tips in the tree. For more information about Strahler numbers see \url{https://en.wikipedia.org/wiki/Strahler_number}. The function \code{extract.strahlerNumber} extracts all of the most inclusive clades of Strahler number \code{i}. } \value{ Either a vector with the Strahler number for each tip and internal node; or (for \code{extract.strahlerNumber} the set of (most inclusive) subtrees with Strahler number \code{i} as an object of class \code{"multiPhylo"}. diff -Nru r-cran-phytools-0.7-47/MD5 r-cran-phytools-0.7-70/MD5 --- r-cran-phytools-0.7-47/MD5 2020-06-01 21:00:02.000000000 +0000 +++ r-cran-phytools-0.7-70/MD5 2020-09-19 22:20:02.000000000 +0000 @@ -1,5 +1,5 @@ -b7b062158fbbe9975b3a4dc4f17d1e24 *DESCRIPTION -cb2681944760521b438b3f20e8960549 *NAMESPACE +cb3cb1f758e68aed77578440c5ccdcbc *DESCRIPTION +314d7193aefe01038124f81631f06d17 *NAMESPACE ec1dee2fc5dce861a4a896e13a10a9a5 *R/Dtest.R ee8f279c7051fd15e5f3b69a78b13d42 *R/add.everywhere.R 00125e455cd6c2a9a162f414c82b31d3 *R/add.random.R @@ -8,24 +8,24 @@ 24a634dd6c26e4e706ba8afef7591b40 *R/anc.Bayes.R 717587314843beab5cc2e340d878fb54 *R/anc.ML.R 227c5163b6303989b36dbe9d50a7f53b *R/anc.trend.R -aec489dfdb755d6cff97e439e29157a5 *R/ancThresh.R +96a051e7661aa1e6e5bb86e2fc88b055 *R/ancThresh.R a47da3dc00dd6df56555a968d2902b01 *R/backbonePhylo.R 3f839fe6b127ee44026ad7f20c32e030 *R/bd.R 9d8aa3bd858bdd62de5df03497e96d2c *R/bmPlot.R -040ca7767b407d159e39b66e90f18237 *R/branching.diffusion.R -3202d9af6f88d617a44909986a210cff *R/brownie.lite.R +894d5f9b5dbf164c92714fd15d41dbf9 *R/branching.diffusion.R +2ce09cebde349b458acdba93c6bc85b9 *R/brownie.lite.R d81faab1ac079af1ce2016b7e004bcce *R/brownieREML.R ac74214ded5662595cc95cbf8c4f49f5 *R/collapseTree.R 1001c3a0359677faf03cd3912d5aacd5 *R/compare.chronograms.R be284dbd922c8e640b06f22aeea30eb6 *R/consensus.edges.R -e44c25cfd1b0c7d1f07400d96396c6c8 *R/contMap.R +de76451e2805b839c1a7120305bc6921 *R/contMap.R ca84810482c72fff11b7cc05c1ca994a *R/cophylo.R c6caa20ceaa72a8a1dc2e71efa9ea19e *R/cospeciation.R -5a714a5ff9ab89fc5d41ca1dbd38e3ce *R/ctt.R +f7ffd39dc1f6c31378ee311508708bfb *R/ctt.R deb83911c92974c2eb3b653cbcbbc6df *R/densityMap.R a651a11526c62a939580958d664aea3d *R/densityTree.R 4d0db3cfa22ee367a93e56c4459c0db4 *R/dotTree.R -70a8db5e7ea6c926c19e792168050491 *R/drop.tip.simmap.R +51ac2df4f93f187de8e0e71a589b4fc2 *R/drop.tip.simmap.R bef780c00ce923e2306295f2dcd7bf8b *R/estDiversity.R 016534b15b0de0b6183169eb2358b998 *R/evol.rate.mcmc.R 2606d5683b834c8576792d3a107553d3 *R/evol.vcv.R @@ -33,21 +33,22 @@ 164c5aba763c4823fae1a2747cbbdde8 *R/exhaustiveMP.R 0c94431c709b011b1a19a5238b7a3352 *R/export.as.xml.R 59c8fa4c521018fcc8a9e8bbe8e8e718 *R/fancyTree.R -a1ebd3392c1313f685add9284121415c *R/fastAnc.R +10c2c1db23b6325cd65f17e4d3549aba *R/fastAnc.R 1fb690b492fa00b02e86eaf49b5e48b6 *R/fastBM.R fd4043fe4370b493b19344fc752bcdf5 *R/fitBayes.R 2180796336a37244d12cd384282ae33b *R/fitDiversityModel.R -af013ec1fb38db66642353cbaabe727a *R/fitMk.R -31916a41cf94b3475d3737ff0ffd6ee2 *R/fitPagel.R -f7c288b33bb7cc270f77eff1af16bfba *R/fitmultiMk.R -297860c8a4d1639ea312e384d08a270a *R/fitpolyMk.R +c4529c3b9702c79aa051b338023f40ab *R/fitHRM.R +0e5dd4f025704ecef036080adafb6fb8 *R/fitMk.R +ef4e25e489eae81b918f39d9aeaf74c8 *R/fitPagel.R +2313c5868bb1a48ceb89a3f12af424b5 *R/fitmultiMk.R +4eac4a53f81fab9432c1748b756da06c *R/fitpolyMk.R 14efc2d387e2360f9f37721f7a6b77a6 *R/locate.fossil.R ae5cd392abeb954e238903bbbfc01239 *R/locate.yeti.R e41cd28e97c6101b3e0c7c4b955abb1a *R/ls.consensus.R -b8c4eb43939eed6b2c6c2b8353743d70 *R/ltt.R +d634251656b1ef6cd2605223d5647898 *R/ltt.R d1669270abede61af74c637d742da68c *R/ltt95.R 07c37a4c09941dfa8088c4eb91ef5c4a *R/make.era.map.R -2b878acd43e9afca86fc10febba22942 *R/make.simmap.R +e689da3727825702e9193f0a9ac85da1 *R/make.simmap.R b535cd4a22f407601477a0c787a8f869 *R/map.overlap.R c292a07e2f329ba3550ab11a573f9e2e *R/map.to.singleton.R bf6c23b046ed30ee9e9183582aceb070 *R/mcmcBM.R @@ -61,56 +62,60 @@ 4251fbfa33b24cedd655a575878bee3a *R/paintSubTree.R 86f26b836e23ca7099d522bc2c8ab79c *R/pbtree.R a189790284945deed14a053d6e0c4c38 *R/pgls.Ives.R -0a3f71e1a3714c0fc5ed5fce53b9aa54 *R/phenogram.R +db533074aef3ced57817fd11bf0a26c5 *R/phenogram.R 93cdbb83a29c009a6381af932cb78b87 *R/phyl.RMA.R 3cbec3ddb43718bf824e1b4b5b1bb92b *R/phyl.cca.R 2e0314f7d261000550c22b529bfb6908 *R/phyl.pairedttest.R -2550b66a5f43a0e214c6b6ae38efe40d *R/phyl.pca.R +9ea7e788ffc6c18279f2e01dffc1a641 *R/phyl.pca.R 17d0e8f613fafc8cf7e0a719120dd9cd *R/phyl.resid.R 9265547de9fc901a0e0af7d58706b142 *R/phylANOVA.R -fc18886487d639e3c2c400b48f34648e *R/phylo.heatmap.R +82754e445e944f49ccacc078a185f36c *R/phylo.heatmap.R 3fd0311c0d039889d2cca7347e3aaee8 *R/phylo.impute.R -57c437adb54a1e56d5144cea6ef1ab87 *R/phylo.to.map.R -cf5760ebfd99d61a222daca28a054bd1 *R/phylomorphospace.R +e6c26951fcaebf42542b6fbc1affb620 *R/phylo.to.map.R +963b45797954ac24ffab58bc0a91551b *R/phylomorphospace.R ba04cc646aeb4d77be7af8bf3b1f21ac *R/phylomorphospace3d.R -4c68b9ff8ed3e59069e9d286624b463a *R/phylosig.R -cd926e56e60606371ed57490b57ab79e *R/plotBranchbyTrait.R +a321ac68ac46f667e1fdd45cd53c3410 *R/phylosig.R +3a5045f97da2d4f697eea1e86e46edb7 *R/plotBranchbyTrait.R 7b978ab147f37fa38a8d22412e6969d7 *R/plotSimmap.R -c9f1bff41eba65d96f74ad36910f0ca4 *R/plotTree.datamatrix.R +d6a34cd3cc500aafcfaab6113e5a29c6 *R/plotTree.datamatrix.R a58e6ab65dca2e69dde3266c9d05e51e *R/plotTree.errorbars.R 9f36a04d3449c0f6a4e59b0154ee802a *R/plotTree.wBars.R bc0699095386afd623fa2385c68eb218 *R/project.phylomorphospace.R de8c966ad29e2ef0b7d6838e34bd8f81 *R/ratebystate.R c393517b7ce4f9ccd5a44dfa4d541583 *R/ratebytree.R -5f808f6878e031af61d4042b0dcbcb5e *R/rateshift.R +47895a08678f00f93229807ffc536c57 *R/rateshift.R 10aa4276f62187281eaa23d509df6f58 *R/read.newick.R 613a9311d192d89fefaf7f5c6e039de3 *R/read.simmap.R -782540c3e7967df671440da8834fbb5d *R/rerootingMethod.R +2658841a3af53c82b86513482e4f0bd2 *R/rerootingMethod.R ba26eefe861f14105c46ee13652c5b66 *R/resolveNodes.R 9bb4ea1e15eaeba8b3ea0f4ffc3c9255 *R/roundPhylogram.R d8a6bb753f6779a1dd49e769b18775a0 *R/sim.corrs.R a02ce76ea678143f57284dae55f34531 *R/sim.history.R f83f44a29df59f5ff489e0120c4c328b *R/sim.rates.R -9deeee566b523eaedbbabf8586533c2e *R/simBMphylo.R +de02f79e452899513b9e72530461e288 *R/simBMphylo.R 5ae6a2a091d203ecc04ad6081e11293d *R/skewers.R 606ad5a902660277dce5f3ef7fed7461 *R/splitplotTree.R c8bfa07dd3238d97fcaafe5a369e5e89 *R/starTree.R 0567456c325cbac2c305158ffa05df91 *R/strahlerNumber.R -7f60fdf45387dfe3ae42629306a202a1 *R/threshBayes.R +87cc055735bb22fe3784662ef6f24abb *R/threshBayes.R 5a8a89c2fe9a2e54e4d7ef961e241198 *R/tree.grow.R 7e547545cd50de19581af1e1c179e967 *R/treeSlice.R -01efb8d2e615551bd3a7b3aea009cb6b *R/utilities.R +08dcee63ed285b0babccc79abe56bea3 *R/utilities.R 4c418f19c184c846f28f696969142bd0 *R/write.simmap.R 67b049f3e083513ff41ab9ca2bd940b5 *R/writeAncestors.R 65bf87a538f44806311e8914317b242c *R/writeNexus.R 8e78dfc6a432989f6cf7c6d9459b88d9 *data/anole.data.rda 0cc04b5f3865c524c8ec36dc84eccef6 *data/anoletree.rda -f23addc91d7cec7da27faeb56a937b80 *data/datalist -3cea6deec311a975254717398b8dac68 *data/mammal.data.rda -18633b2cb0a80c8db5213fcaf7df04a9 *data/mammal.tree.rda +9e309b71d0240676b06fc8ff09f56f38 *data/datalist +0416c11b7757a8a1b7d4ac5b1e4b4fb5 *data/flatworm.data.rda +0a5fc0e193972b61d818e4fdf6b733d1 *data/flatworm.tree.rda +29bd6b204223041383a082d85e235e9a *data/mammal.data.rda +114305e2f421cba146fd300da570ab92 *data/mammal.tree.rda 5807bfb760e5b9d95822d9efe4524639 *data/salamanders.rda c7367501b639064f0b47dbc088213730 *data/sunfish.data.rda 5e76e37fc53b3002673868b074ca89ff *data/sunfish.tree.rda +e3023f4223c0c3863596df1581c4a468 *data/vertebrate.data.rda +f5e1db184a75019f551868c00ed2293a *data/vertebrate.tree.rda b7fbd42a08d4aca73426dbb42de518c3 *data/wasp.data.rda 5f4d2752a6176f1bc9ec31cab308b7ce *data/wasp.trees.rda 0eb3bef3466bd6d264c669571621c2d1 *inst/CITATION @@ -127,9 +132,9 @@ 8dc3065b1e42478de13d946eb134e67e *man/anc.ML.Rd 0d5d3d8a34b5631dcc4b36414678d290 *man/anc.trend.Rd 7b21c6c1faa82f81ad95d2584761a3c9 *man/ancThresh.Rd -333cfa037f2a47e0eaafd9bebcde1193 *man/anoletree.Rd +dc0c204ce220f15c233b29e9e735c03e *man/anoletree.Rd ba7274925fef31f711d97b43ea309e57 *man/applyBranchLengths.Rd -f03624a71c7986344d6016b084100b5a *man/as.Qmatrix.Rd +a6d1ef15876e7dc57e9d7dd6750bcd33 *man/as.Qmatrix.Rd 9438d11405f55849b8b13173e0aa677f *man/as.multiPhylo.Rd f92f0ea6ad70e2e1a10dbc4c47b4d79e *man/ave.rates.Rd 09f837558c285442823a52ef98631230 *man/averageTree.Rd @@ -158,8 +163,10 @@ 76ba9e500f6c1f963cbd878a39ed0752 *man/dotTree.Rd c92235efc1cf6791a2871a0b59bf94b1 *man/drop.clade.Rd 139cf35fd81a07386ea5489a3c8a549e *man/drop.leaves.Rd -4dbacec8bf2b24739708948c5c6f2251 *man/drop.tip.contMap.Rd +6a5ce468aa63a264cf9f755f1af20edf *man/drop.tip.contMap.Rd +622f7e000f871efa589f42fecb153c22 *man/drop.tip.multiPhylo.Rd 0f4c4d6f030921fe8f033e54a5ba2226 *man/drop.tip.simmap.Rd +847828b538311f8b3c10e46b402b28a6 *man/edge.widthMap.Rd a856f0f113e72db407bb0f4d3d3816ff *man/edgeProbs.Rd 5fe47b2846351ede4b8190aa483fed51 *man/estDiversity.Rd 029397a04141e0ac95659f6c80a367c3 *man/evol.rate.mcmc.Rd @@ -176,7 +183,7 @@ f403634da6aaa52e261930730fbbc669 *man/fit.bd.Rd 8c9a8bf66fa03d03fc4d30902bba058e *man/fitBayes.Rd 50d9e4e565e23081519af79600a7fc0b *man/fitDiversityModel.Rd -4e877a6c5099c3720d065430367061ae *man/fitMk.Rd +2423d3313f375f7a1f00615bcb85f78e *man/fitMk.Rd bec587cf8339fb2c411708eb27b2389e *man/fitPagel.Rd ba9a3f05a5e21858e0b73853b71893f0 *man/force.ultrametric.Rd 1b8c690f7a0c981bdeb6023ceea6986d *man/gammatest.Rd @@ -235,18 +242,18 @@ 6b3088a2e0a7f3afa415a1a3567c3c68 *man/phylo.to.map.Rd 95aa8b498df77e1b9ea962512e42ce2f *man/phylo.toBackbone.Rd 3a2d310bb460a7a129fc6fbffc3692ad *man/phyloDesign.Rd -19c19789a7aa60a16e00dd704b46cb84 *man/phylomorphospace.Rd +204134b0c17cb8aabdded5a994dc544c *man/phylomorphospace.Rd 0860a04032f550a8c1d32a53523a9c57 *man/phylomorphospace3d.Rd 4b20335dd81ce22fb269f0f594e12580 *man/phylosig.Rd -1803382cd66705cc6014798294212a73 *man/phytools-package.Rd +9b9e628b355be7cf665647cbd922f4c2 *man/phytools-package.Rd aecc4ccca42174f15c6280dbb5aa13fe *man/plot.backbonePhylo.Rd c414a71b143bca6924c8ff87c8808f96 *man/plotBranchbyTrait.Rd a37f0674203b4beb950a5839bdbef9af *man/plotSimmap.Rd 29b91327ad06627986ba73959ff1afaa *man/plotThresh.Rd -7abaec7aa366e41647df0bf3073a36ef *man/plotTree.Rd +b7c741a12fcc5c028d021260351d925a *man/plotTree.Rd d6e4d7fa3d7e518b7c75318a1438fcd1 *man/plotTree.datamatrix.Rd 028fb0ca902c0cf9ffcdc47ad85da56e *man/plotTree.errorbars.Rd -2835404ddf5c5f78559c43c91ad46d83 *man/plotTree.wBars.Rd +f77f83cffdb74d5d1d7dc60781beea93 *man/plotTree.wBars.Rd 1b2c88a32aa808cc49ebc2e769a85b9f *man/posterior.evolrate.Rd 9f409f4328ac4802d1700eec2233361b *man/posthoc.Rd 4bfa72a7ac47c3e185164f8c35dd5e13 *man/print.backbonePhylo.Rd @@ -278,7 +285,7 @@ f0b7c55fe94d87a63fcc473557ba356c *man/splitTree.Rd c4c25d51cb5afc0a3ca9e8c50aa4f8a5 *man/splitplotTree.Rd e3ebd622b81458023c71dc8b846bc8bc *man/starTree.Rd -e3cb564ff8ac467e211b63893e5bd8ab *man/strahlerNumber.Rd +e66a9557f578c6102a325e45413fa26d *man/strahlerNumber.Rd 0ffebdf089a196fac35debbf50ac3add *man/threshBayes.Rd 2066958f8984e50838e5da011b83c841 *man/threshDIC.Rd 4c42ea4cbb70dd27915a35eb6759a176 *man/threshState.Rd diff -Nru r-cran-phytools-0.7-47/NAMESPACE r-cran-phytools-0.7-70/NAMESPACE --- r-cran-phytools-0.7-47/NAMESPACE 2020-06-01 04:06:19.000000000 +0000 +++ r-cran-phytools-0.7-70/NAMESPACE 2020-08-29 00:09:22.000000000 +0000 @@ -8,16 +8,18 @@ export(consensus.edges, contMap, cophylo, countSimmap, compute.mr, cospeciation, ctt) export(density.anc.Bayes, density.multiSimmap, densityMap, densityTree, describe.simmap) export(di2multi.simmap, dot.legend, dotTree, drop.clade, drop.leaves, drop.tip.contMap) -export(drop.tip.densityMap, drop.tip.simmap, drop.tip.singleton, Dtest) -export(edgelabels.cophylo, edgeProbs, errorbar.contMap, estDiversity, evol.rate.mcmc) -export(evol.vcv, evolvcv.lite, exhaustiveMP, expand.clade, export.as.xml) +export(drop.tip.densityMap, drop.tip.multiPhylo, drop.tip.multiSimmap, drop.tip.simmap) +export(drop.tip.singleton, Dtest) +export(edge.widthMap, edgelabels.cophylo, edgeProbs, errorbar.contMap, estDiversity) +export(evol.rate.mcmc, evol.vcv, evolvcv.lite, exhaustiveMP, expand.clade, export.as.xml) export(extract.clade.simmap, extract.strahlerNumber) export(fancyTree, fastAnc, fastBM, fastDist, fastHeight, fastMRCA, findMRCA, fit.bd) -export(fit.yule, fitBayes, fitMk, fitmultiMk, fitpolyMk, fitDiversityModel, fitPagel) -export(force.ultrametric) +export(fit.yule, fitBayes, fitHRM, fitMk, fitmultiMk, fitpolyMk, fitDiversityModel) +export(fitPagel, force.ultrametric) export(gammatest, genus.to.species.tree, genSeq, geo.palette, geo.legend, get.treepos) export(getCladesofSize, getDescendants, getExtant, getExtinct, getnode, getParent) -export(getSisters, getStates, gtt) +export(getSisters, getStates, graph.polyMk, gtt) +export(keep.tip.contMap) export(labelnodes, ladderize.simmap, lambda.transform, lik.bd, likMlambda) export(likSurface.rateshift, linklabels, locate.fossil, locate.yeti) export(ls.consensus, ls.tree, ltt, ltt95) @@ -32,17 +34,18 @@ export(phyl.vcv, phylo.heatmap, phylo.impute, phylo.toBackbone, phylo.to.map, phylANOVA) export(phyloDesign, phylomorphospace, phylomorphospace3d, phyloScattergram, phylosig) export(plot.anc.Bayes, plot.changesMap, plot.contMap, plot.cophylo, plot.densityMap) -export(plot.fitMk, plot.fitPagel, plot.gfit, plot.phyl.RMA, plot.phylo.to.map, plot.Qmatrix) -export(plot.simBMphylo, plotBranchbyTrait, plot.phylosig, plotSimmap, plotThresh, plotTree) -export(plotTree.barplot, plotTree.boxplot, plotTree.datamatrix, plotTree.errorbars) -export(plotTree.singletons, plotTree.splits, plotTree.wBars, posterior.evolrate, posthoc) -export(posthoc.ratebytree, print.Qmatrix, project.phylomorphospace) +export(plot.edge.widthMap, plot.fitMk, plot.fitPagel, plot.gfit, plot.phyl.RMA) +export(plot.phylo.to.map, plot.Qmatrix, plot.simBMphylo, plotBranchbyTrait, plot.phylosig) +export(plotSimmap, plotThresh, plotTree, plotTree.barplot, plotTree.boxplot) +export(plotTree.datamatrix, plotTree.errorbars, plotTree.singletons, plotTree.splits) +export(plotTree.wBars, posterior.evolrate, posthoc, posthoc.ratebytree, print.Qmatrix) +export(project.phylomorphospace) export(ratebystate, ratebytree, rateshift, read.newick, read.simmap, reorder.backbonePhylo) export(reorderSimmap, rep.multiPhylo, rep.phylo, repPhylo, reroot, rerootingMethod) export(rescaleSimmap, resolveAllNodes, resolveNode, rootedge.to.singleton, rotate.multi) export(rotateNodes, roundBranches, roundPhylogram, rstate) -export(sampleFrom, scores, scores.phyl.pca, setMap, sim.corrs, sim.ctt, sim.history, sim.Mk) -export(sim.multiCtt, sim.multiMk, sim.ratebystate, sim.rates, simBMphylo, skewers) +export(sampleFrom, scores, scores.phyl.pca, setMap, sim.corrs, sim.ctt, sim.history) +export(sim.Mk, sim.multiCtt, sim.multiMk, sim.ratebystate, sim.rates, simBMphylo, skewers) export(splitEdgeColor, splitplotTree, splitTree, starTree, strahlerNumber) export(threshBayes, threshDIC, threshState, tiplabels.cophylo, tipRotate, to.matrix) export(tree.grow, treeSlice) @@ -203,30 +206,46 @@ S3method(summary, anc.Bayes) S3method(plot, Qmatrix) S3method(logLik, anc.ML) +S3method(print, edge.widthMap) +S3method(plot, edge.widthMap) +S3method(plot, fitHRM) +S3method(print, fitHRM) +S3method(as.Qmatrix, corhmm) -importFrom(animation, ani.options, ani.record, ani.replay, saveVideo) -importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin, as.phylo, bind.tree, branching.times, collapse.singles) -importFrom(ape, consensus, compute.brlen, cophenetic.phylo, corBrownian, di2multi, dist.dna, dist.nodes, drop.tip, edgelabels, extract.clade) -importFrom(ape, getMRCA, is.binary) -importFrom(ape, is.monophyletic, is.rooted, is.ultrametric, ladderize, matexpo, mrca, multi2di) -importFrom(ape, nodelabels, Ntip, pic, plot.phylo, prop.part, read.tree, reorder.phylo, root.phylo, rotate, rtree, stree, tiplabels) -importFrom(ape, unroot, vcv, vcv.phylo, write.tree) +importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin) +importFrom(ape, as.phylo, bind.tree, branching.times, collapse.singles, consensus) +importFrom(ape, compute.brlen, cophenetic.phylo, corBrownian, di2multi, dist.dna) +importFrom(ape, dist.nodes, drop.tip, edgelabels, extract.clade, getMRCA, is.binary) +importFrom(ape, is.monophyletic, is.rooted, is.ultrametric, ladderize, matexpo, mrca) +importFrom(ape, multi2di, nodelabels, Ntip, pic, plot.phylo, prop.part, read.tree) +importFrom(ape, reorder.phylo, root.phylo, rotate, rtree, stree, tiplabels, unroot, vcv) +importFrom(ape, vcv.phylo, write.tree) +importFrom(coda, as.mcmc, HPDinterval) +importFrom(combinat, permn) importFrom(clusterGeneration, genPositiveDefMat) +importFrom(expm, expm) +importFrom(graphics, strwidth, par, segments, locator, lines, text, strheight, symbols) +importFrom(graphics, plot, layout, plot.new, title, axis, points, plot.window, polygon) +importFrom(graphics, rect, curve, image, mtext, arrows, barplot, boxplot, contour, hist) +importFrom(graphics, abline, legend, grid) +importFrom(grDevices, palette, colors, dev.hold, dev.flush, rainbow, heat.colors, gray) +importFrom(grDevices, colorRampPalette, dev.new, rgb, pdf, dev.off, col2rgb) +importFrom(gtools, combinations) importFrom(maps, map) +importFrom(MASS, ginv) +importFrom(methods, hasArg) importFrom(mnormt, dmnorm, pd.solve) +importFrom(nlme, gls, varFixed) importFrom(numDeriv, hessian) -importFrom(phangorn, allTrees, Ancestors, as.phyDat, dist.hamming, midpoint, NJ, nni, nnls.tree, optim.parsimony, parsimony, phyDat, pratchet, treedist, RF.dist, KF.dist, path.dist, SPR.dist, Children, Descendants, threshStateC) +importFrom(phangorn, allTrees, Ancestors, as.phyDat, dist.hamming, midpoint, NJ, nni) +importFrom(phangorn, nnls.tree, optim.parsimony, parsimony, phyDat, pratchet, treedist) +importFrom(phangorn, RF.dist, KF.dist, path.dist, SPR.dist, Children, Descendants) +importFrom(phangorn, threshStateC) importFrom(plotrix, arctext, draw.arc, draw.circle, draw.ellipse, textbox) importFrom(scatterplot3d, scatterplot3d) -importFrom(stats, cophenetic, reorder, runif, setNames, optim, dexp, dnorm, rnorm, var, nlminb, biplot, pchisq, rexp, optimize, logLik, as.dist, pnorm, median, cor, aggregate, rgamma, dgamma, lm, rbinom, rgeom, pt, anova, p.adjust, screeplot, rchisq, optimHess, rmultinom, sd, dunif, dist, plogis, density, coef, cmdscale) -importFrom(methods, hasArg) -importFrom(graphics, strwidth, par, segments, locator, lines, text, strheight, symbols, plot, layout, plot.new, title, axis, points, plot.window, polygon, rect, curve, image, mtext, arrows, barplot, boxplot, contour, hist, abline, legend, grid) +importFrom(stats, cophenetic, reorder, runif, setNames, optim, dexp, dnorm, rnorm, var) +importFrom(stats, nlminb, biplot, pchisq, rexp, optimize, logLik, as.dist, pnorm, median) +importFrom(stats, cor, aggregate, rgamma, dgamma, lm, rbinom, rgeom, pt, anova, p.adjust) +importFrom(stats, screeplot, rchisq, optimHess, rmultinom, sd, dunif, dist, plogis) +importFrom(stats, density, coef, cmdscale) importFrom(utils, flush.console, head, installed.packages, str, capture.output) -importFrom(grDevices, palette, colors, dev.hold, dev.flush, rainbow, heat.colors, gray, colorRampPalette, dev.new, rgb, -pdf, dev.off, col2rgb) -importFrom(combinat, permn) -importFrom(coda, as.mcmc, HPDinterval) -importFrom(nlme, gls, varFixed) -importFrom(MASS, ginv) -importFrom(expm, expm) -importFrom(gtools, combinations) diff -Nru r-cran-phytools-0.7-47/R/ancThresh.R r-cran-phytools-0.7-70/R/ancThresh.R --- r-cran-phytools-0.7-47/R/ancThresh.R 2020-05-30 19:25:56.000000000 +0000 +++ r-cran-phytools-0.7-70/R/ancThresh.R 2020-07-17 15:04:09.000000000 +0000 @@ -4,10 +4,20 @@ ancThresh<-function(tree,x,ngen=100000,sequence=NULL,method="mcmc", model=c("BM","OU","lambda"),control=list(),...){ - if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".") + if(!inherits(tree,"phylo")) + stop("tree should be an object of class \"phylo\".") + + if(hasArg(auto.tune)) auto.tune<-list(...)$auto.tune + else auto.tune<-TRUE + if(is.logical(auto.tune)) if(auto.tune==TRUE) auto.tune<-0.234 + else if(is.numeric(auto.tune)) if(auto.tune>=1.0||auto.tune<=0.0){ + cat("value for auto.tune outside allowable range. Resetting....\n") + auto.tune<-0.234 + } # check method - if(method!="mcmc") stop(paste(c("do not recognize method =",method,",quitting"))) + if(method!="mcmc") stop(paste(c("do not recognize method =",method, + ",quitting"))) # get model for the evolution of liability model<-model[1] @@ -34,9 +44,22 @@ X<-X[,seq] # order columns by seq # ok, now set starting thresholds - th<-c(1:length(seq))-1; names(th)<-seq + th<-c(1:length(seq))-1 + names(th)<-seq x<-to.vector(X) - l<-sapply(x,function(x) runif(n=1,min=th[x]-1,max=th[x])) # set plausible starting liability + ## set plausible starting liabilities + MU<-mean(th) + SD<-2*sd(th) + # now change the upper limit of th to Inf + th[length(th)]<-Inf + l<-setNames(vector(mode="numeric",length=length(x)),names(x)) + for(i in 1:length(l)){ + l[i]<-rnorm(n=1,mean=MU,sd=SD) + while(threshState(l[i],thresholds=th)!=x[i]) l[i]<-rnorm(n=1,mean=MU, + sd=SD) + } + ## set plausible starting liability + ## l<-sapply(x,function(x) runif(n=1,min=th[x]-1,max=th[x])) if(model=="OU") alpha<-0.1*max(nodeHeights(tree)) if(model=="lambda") lambda<-1.0 @@ -77,10 +100,13 @@ # now set ancestral liabilities, by first picking ancestral states from their prior temp<-apply(con$pr.anc,1,rstate) # assign random liabilities consistent with the starting thresholds - a<-sapply(temp,function(x) runif(n=1,min=th[x]-1,max=th[x])) - - # now change the upper limit of th to Inf - th[length(th)]<-Inf + ## a<-sapply(temp,function(x) runif(n=1,min=th[x]-1,max=th[x])) + a<-setNames(vector(mode="numeric",length=length(temp)),names(temp)) + for(i in 1:length(a)){ + a[i]<-rnorm(n=1,mean=MU,sd=SD) + while(threshState(a[i],thresholds=th)!=temp[i]) a[i]<-rnorm(n=1,mean=MU, + sd=SD) + } # compute some matrices & values V<-if(model=="BM") vcvPhylo(tree) @@ -94,14 +120,21 @@ lik1<-likLiab(l,a,V,invV,detV)+log(probMatch(X,l,th,seq)) # store - A<-matrix(NA,ngen/con$sample+1,tree$Nnode,dimnames=list(NULL,n+1:tree$Nnode)) - B<-if(model=="BM") matrix(NA,ngen/con$sample+1,m+2,dimnames=list(NULL,c("gen",names(th),"logLik"))) - else if(model=="OU") matrix(NA,ngen/con$sample+1,m+3,dimnames=list(NULL,c("gen",names(th),"alpha","logLik"))) - else if(model=="lambda") matrix(NA,ngen/con$sample+1,m+3,dimnames=list(NULL,c("gen",names(th),"lambda","logLik"))) + A<-matrix(NA,ngen/con$sample+1,tree$Nnode,dimnames=list(NULL, + n+1:tree$Nnode)) + B<-if(model=="BM") matrix(NA,ngen/con$sample+1,m+2, + dimnames=list(NULL,c("gen",names(th),"logLik"))) + else if(model=="OU") matrix(NA,ngen/con$sample+1,m+3, + dimnames=list(NULL,c("gen",names(th),"alpha","logLik"))) + else if(model=="lambda") matrix(NA,ngen/con$sample+1,m+3, + dimnames=list(NULL,c("gen",names(th),"lambda","logLik"))) - C<-matrix(NA,ngen/con$sample+1,tree$Nnode+n,dimnames=list(NULL,c(tree$tip.label,1:tree$Nnode+n))) + C<-matrix(NA,ngen/con$sample+1,tree$Nnode+n,dimnames=list(NULL, + c(tree$tip.label,1:tree$Nnode+n))) A[1,]<-threshState(a,thresholds=th) - B[1,]<-if(model=="BM") c(0,th,lik1) else if(model=="OU") c(0,th,alpha,lik1) else if(model=="lambda") c(0,th,lambda,lik1) + B[1,]<-if(model=="BM") c(0,th,lik1) else + if(model=="OU") c(0,th,alpha,lik1) else + if(model=="lambda") c(0,th,lambda,lik1) C[1,]<-c(l[tree$tip.label],a[as.character(1:tree$Nnode+n)]) # run MCMC @@ -143,7 +176,8 @@ } } lik2<-likLiab(lp,ap,Vp,invVp,detVp)+log(probMatch(X,lp,thp,seq)) - p.odds<-min(c(1,exp(lik2+logPrior(threshState(ap,thresholds=thp),thp,con)-lik1-logPrior(threshState(a,thresholds=th),th,con)))) + p.odds<-min(c(1,exp(lik2+logPrior(threshState(ap,thresholds=thp),thp,con)- + lik1-logPrior(threshState(a,thresholds=th),th,con)))) if(p.odds>runif(n=1)){ a<-ap; l<-lp; th<-thp @@ -154,7 +188,9 @@ } else logL<-lik1 if(i%%con$sample==0){ A[i/con$sample+1,]<-threshState(a,thresholds=th) - B[i/con$sample+1,]<-if(model=="BM") c(i,th[colnames(B)[1+1:m]],logL) else if(model=="OU") c(i,th[colnames(B)[1+1:m]],alpha,logL) else if(model=="lambda") c(i,th[colnames(B)[1+1:m]],lambda,logL) + B[i/con$sample+1,]<-if(model=="BM") c(i,th[colnames(B)[1+1:m]],logL) else + if(model=="OU") c(i,th[colnames(B)[1+1:m]],alpha,logL) else + if(model=="lambda") c(i,th[colnames(B)[1+1:m]],lambda,logL) C[i/con$sample+1,]<-c(l[tree$tip.label],a[as.character(1:tree$Nnode+n)]) } } @@ -181,12 +217,14 @@ ## some S3 methods (added in 2017) print.ancThresh<-function(x,...){ - cat("\nObject containing the results from an MCMC analysis\nof the threshold model using ancThresh.\n\n") + cat("\nObject containing the results from an MCMC analysis\n") + cat("of the threshold model using ancThresh.\n\n") cat("List with the following components:\n") cat(paste("ace:\tmatrix with posterior probabilities assuming",x$burnin, "\n\tburn-in generations.\n")) cat("mcmc:\tposterior sample of liabilities at tips & internal\n") - cat(paste("\tnodes (a matrix with",nrow(x$mcmc),"rows &",ncol(x$mcmc),"columns).\n")) + cat(paste("\tnodes (a matrix with",nrow(x$mcmc),"rows &",ncol(x$mcmc), + "columns).\n")) cat("par:\tposterior sample of the relative positions of the\n") cat(paste("\tthresholds, the log-likelihoods, and any other\n", "\tmodel variables (a matrix with",nrow(x$par),"rows).\n\n")) @@ -202,7 +240,10 @@ else burnin<-x$burnin args<-list(...) if(is.null(args$lwd)) args$lwd<-1 - if(is.null(args$ylim)) args$ylim<-c(-0.1*Ntip(x$tree),Ntip(x$tree)) + if(is.null(args$type)) args$type<-"phylogram" + if(is.null(args$ylim)) + if(args$type=="phylogram") + args$ylim<-c(-0.1*Ntip(x$tree),Ntip(x$tree)) if(is.null(args$offset)) args$offset<-0.5 if(is.null(args$ftype)) args$ftype="i" args$tree<-x$tree @@ -212,8 +253,8 @@ THRESH<-as.matrix(x$par)[ii:nrow(x$par),1:length(x$seq)+1] STATES<-matrix(NA,nrow(LIAB),ncol(LIAB),dimnames=dimnames(LIAB)) for(i in 1:nrow(LIAB)) STATES[i,]<-threshState(LIAB[i,],THRESH[i,]) - PP<-t(apply(STATES,2,function(x,levs) summary(factor(x,levels=levs))/length(x), - levs=x$seq)) + PP<-t(apply(STATES,2,function(x,levs) + summary(factor(x,levels=levs))/length(x),levs=x$seq)) if(hasArg(piecol)) piecol<-list(...)$piecol else piecol<-setNames(colorRampPalette(c("blue", "yellow"))(length(x$seq)),x$seq) @@ -225,7 +266,7 @@ else tip.cex<-0.4 tiplabels(pie=PP[x$tree$tip.label,],piecol=piecol, cex=tip.cex) - legend(x=par()$usr[1],y=par()$usr[1],legend=x$seq,pch=21,pt.bg=piecol, + legend(x="bottomleft",legend=x$seq,pch=21,pt.bg=piecol, pt.cex=2.2,bty="n") invisible(PP) } @@ -245,11 +286,14 @@ # plot tree par(lend=2) - plotTree(tree,ftype="i",lwd=1,ylim=if(legend) c(-0.1*length(tree$tip.label),length(tree$tip.label)) else NULL,...) + plotTree(tree,ftype="i",lwd=1,ylim=if(legend) c(-0.1*length(tree$tip.label), + length(tree$tip.label)) else NULL,...) if(legend){ zz<-par()$cex; par(cex=0.6) for(i in 1:length(piecol)) - add.simmap.legend(leg=leg[i],colors=piecol[i],prompt=FALSE,x=0.02*max(nodeHeights(tree)),y=-0.1*length(tree$tip.label),vertical=FALSE,shape="square",fsize=1) + add.simmap.legend(leg=leg[i],colors=piecol[i],prompt=FALSE, + x=0.02*max(nodeHeights(tree)),y=-0.1*length(tree$tip.label),vertical=FALSE, + shape="square",fsize=1) par(cex=zz) } # pull matrices from mcmc @@ -277,8 +321,10 @@ # plot tip labels if(tipcol=="input") tiplabels(pie=X,piecol=piecol[colnames(X)],cex=0.6) else if(tipcol=="estimated") { - XX<-matrix(NA,nrow(liab),length(tree$tip),dimnames=list(rownames(liab),colnames(liab)[1:length(tree$tip)])) - for(i in 1:nrow(liab)) XX[i,]<-threshState(liab[i,1:length(tree$tip)],thresholds=param[i,1:ncol(X)+1]) + XX<-matrix(NA,nrow(liab),length(tree$tip),dimnames=list(rownames(liab), + colnames(liab)[1:length(tree$tip)])) + for(i in 1:nrow(liab)) XX[i,]<-threshState(liab[i,1:length(tree$tip)], + thresholds=param[i,1:ncol(X)+1]) X<-t(apply(XX,2,function(x) summary(factor(x,levels=colnames(X))))) tiplabels(pie=X/rowSums(X),piecol=piecol[colnames(X)],cex=0.6) } @@ -319,9 +365,14 @@ thBar<-colMeans(mcmc$par[start:nrow(mcmc$par),2:(ncol(mcmc$par)-k)]) liabBar<-colMeans(mcmc$liab[start:nrow(mcmc$liab),]) if(model=="BM") V<-vcvPhylo(tree) - else if(model=="OU") V<-vcvPhylo(tree,model="OU",alpha=mean(mcmc$par[start:nrow(mcmc$par),"alpha"])) - else if(model=="lambda") V<-vcvPhylo(tree,model="lambda",lambda=mean(mcmc$par[start:nrow(mcmc$par),"lambda"])) - Dtheta<--2*(likLiab(liabBar[tree$tip.label],liabBar[as.character(1:tree$Nnode+length(tree$tip))],V,solve(V),determinant(V,logarithm=TRUE)$modulus[1])+log(probMatch(X[tree$tip.label,],liabBar[tree$tip.label],thBar,seq))) + else if(model=="OU") V<-vcvPhylo(tree,model="OU", + alpha=mean(mcmc$par[start:nrow(mcmc$par),"alpha"])) + else if(model=="lambda") V<-vcvPhylo(tree,model="lambda", + lambda=mean(mcmc$par[start:nrow(mcmc$par),"lambda"])) + Dtheta<--2*(likLiab(liabBar[tree$tip.label], + liabBar[as.character(1:tree$Nnode+length(tree$tip))],V,solve(V), + determinant(V,logarithm=TRUE)$modulus[1])+log(probMatch(X[tree$tip.label,], + liabBar[tree$tip.label],thBar,seq))) D<--2*mcmc$par[start:nrow(mcmc$par),"logLik"] Dbar<-mean(D) if(method=="pD"){ diff -Nru r-cran-phytools-0.7-47/R/branching.diffusion.R r-cran-phytools-0.7-70/R/branching.diffusion.R --- r-cran-phytools-0.7-47/R/branching.diffusion.R 2020-06-01 00:25:50.000000000 +0000 +++ r-cran-phytools-0.7-70/R/branching.diffusion.R 2020-06-06 21:56:44.000000000 +0000 @@ -16,6 +16,10 @@ cat(" record != NULL requires the package \"animation\"\n") cat(" Animation will play but not record\n\n") record<-NULL + ani.options<-function(...) NULL + ani.record<-function(...) NULL + ani.replay<-function(...) NULL + saveVideo<-function(...) NULL } if(!is.null(record)){ if(is.null(path)) path="C:/Program Files/ffmpeg/bin/ffmpeg.exe" diff -Nru r-cran-phytools-0.7-47/R/brownie.lite.R r-cran-phytools-0.7-70/R/brownie.lite.R --- r-cran-phytools-0.7-47/R/brownie.lite.R 2019-09-05 23:40:34.000000000 +0000 +++ r-cran-phytools-0.7-70/R/brownie.lite.R 2020-09-15 16:05:25.000000000 +0000 @@ -42,8 +42,9 @@ method="L-BFGS-B",lower=l) logL2<--model2$value while(logL20) print(fits[[i]]) + logL<-sapply(fits,logLik) + if(!quiet){ + cat(paste("log-likelihood from current iteration:", + logLik(fits[[i]]),"\n")) + cat(paste(" --- Best log-likelihood so far:",max(logL), + "---\n")) + flush.console() + } + } + obj<-fits[[which(logL==max(logL))[1]]] + obj$ncat<-ncat + obj$model<-MODEL + obj$umbral<-umbral + obj$all.fits<-fits + obj$data<-X + class(obj)<-c("fitHRM","fitMk") + obj +} + +## print method for objects of class "fitHRM" +print.fitHRM<-function(x,digits=6,...){ + cat("Object of class \"fitHRM\".\n\n") + ss<-unique(sapply(x$states,function(x) strsplit(x,"*",fixed=TRUE)[[1]][1])) + cat(paste("Observed states: [ ",paste(ss,collapse=", ")," ]\n", + sep="")) + cat(paste("Number of rate categories per state: [ ", + paste(x$ncat,collapse=", ")," ]\n\n",sep="")) + cat("Fitted (or set) value of Q:\n") + Q<-matrix(NA,length(x$states),length(x$states)) + Q[]<-c(0,x$rates)[x$index.matrix+1] + diag(Q)<-0 + diag(Q)<--rowSums(Q) + colnames(Q)<-rownames(Q)<-x$states + print(round(Q,digits)) + cat("\nFitted (or set) value of pi:\n") + print(round(x$pi,digits)) + cat(paste("due to treating the root prior as (a) ",x$root.prior,".\n", + sep="")) + cat(paste("\nLog-likelihood:",round(x$logLik,digits),"\n")) + cat(paste("\nOptimization method used was \"",x$method,"\"\n\n", + sep="")) +} + +isOdd<-function(x) (x%%2)==1 + +## S3 plot method for objects of class "fitHRM" +plot.fitHRM<-function(x,...){ + if(hasArg(tol)) tol<-list(...)$tol + else tol<-1e-6 + umbral<-x$umbral + ncat<-x$ncat + Q<-as.Qmatrix(x) + II<-x$index.matrix + for(i in 1:nrow(II)) for(j in 1:ncol(II)) + if(Q[i,j]1) sum(x$ncat[1:(i-1)]) else 0 + mm[1:x$ncat[i]+cs]<-if(isOdd(i)) nn[1:x$ncat[i]+cs] else nn[x$ncat[i]:1+cs] + } + if(isOdd(length(x$ncat))) mm<-nn + Q<-Q[mm,mm] + class(Q)<-"Qmatrix" + args.list<-list(...) + args.list$show.zeros<-FALSE + plot(Q,show.zeros=FALSE,umbral=umbral,ncat=ncat,...) +} + +as.Qmatrix.corhmm<-function(x,...){ + Q<-x$solution + Q[is.na(Q)]<-0 + diag(Q)<--rowSums(Q) + class(Q)<-"Qmatrix" + Q +} diff -Nru r-cran-phytools-0.7-47/R/fitMk.R r-cran-phytools-0.7-70/R/fitMk.R --- r-cran-phytools-0.7-47/R/fitMk.R 2020-06-01 04:04:54.000000000 +0000 +++ r-cran-phytools-0.7-70/R/fitMk.R 2020-08-28 15:44:44.000000000 +0000 @@ -11,6 +11,10 @@ else opt.method<-"nlminb" if(hasArg(min.q)) min.q<-list(...)$min.q else min.q<-1e-12 + if(hasArg(max.q)) max.q<-list(...)$max.q + else max.q<-max(nodeHeights(tree))*100 + if(hasArg(logscale)) logscale<-list(...)$logscale + else logscale<-FALSE N<-Ntip(tree) M<-tree$Nnode if(is.matrix(x)){ @@ -25,14 +29,24 @@ } if(hasArg(pi)) pi<-list(...)$pi else pi<-"equal" - if(pi[1]=="equal") pi<-setNames(rep(1/m,m),states) - else if(pi[1]=="estimated"){ - pi<-if(!is.null(fixedQ)) statdist(fixedQ) else statdist(summary(fitMk(tree,x,model),quiet=TRUE)$Q) - cat("Using pi estimated from the stationary distribution of Q assuming a flat prior.\npi =\n") + if(is.numeric(pi)) root.prior<-"given" + if(pi[1]=="equal"){ + pi<-setNames(rep(1/m,m),states) + root.prior<-"flat" + } else if(pi[1]=="estimated"){ + pi<-if(!is.null(fixedQ)) statdist(fixedQ) else + statdist(summary(fitMk(tree,x,model),quiet=TRUE)$Q) + cat(paste("Using pi estimated from the stationary", + "distribution of Q assuming a flat prior.\npi =\n")) print(round(pi,6)) cat("\n") - } - else pi<-pi/sum(pi) + root.prior<-"stationary" + } else if(pi[1]=="fitzjohn") root.prior<-"nuisance" + if(is.numeric(pi)){ + pi<-pi/sum(pi) + if(is.null(names(pi))) pi<-setNames(pi,states) + pi<-pi[states] + } if(is.null(fixedQ)){ if(is.character(model)){ rate<-matrix(NA,m,m) @@ -70,7 +84,10 @@ rate[rate==0]<-k+1 liks<-rbind(x,matrix(0,M,m,dimnames=list(1:M+N,states))) pw<-reorder(tree,"pruningwise") - lik<-function(Q,output.liks=FALSE,pi){ + lik<-function(Q,output.liks=FALSE,pi,...){ + if(hasArg(output.pi)) output.pi<-list(...)$output.pi + else output.pi<-FALSE + if(is.Qmatrix(Q)) Q<-unclass(Q) if(any(is.nan(Q))||any(is.infinite(Q))) return(1e50) comp<-vector(length=N+M,mode="numeric") parents<-unique(pw$edge[,1]) @@ -84,45 +101,73 @@ for(j in 1:length(v)){ v[[j]]<-EXPM(Q*el[j])%*%liks[desc[j],] } - vv<-if(anc==root) Reduce('*',v)[,1]*pi else Reduce('*',v)[,1] + if(anc==root){ + if(is.numeric(pi)) vv<-Reduce('*',v)[,1]*pi + else if(pi[1]=="fitzjohn"){ + D<-Reduce('*',v)[,1] + pi<-D/sum(D) + vv<-D*D/sum(D) + } + } else vv<-Reduce('*',v)[,1] + ## vv<-if(anc==root) Reduce('*',v)[,1]*pi else Reduce('*',v)[,1] comp[anc]<-sum(vv) liks[anc,]<-vv/comp[anc] } - if(output.liks)return(liks[1:M+N,,drop=FALSE]) - logL<--sum(log(comp[1:M+N])) - return(if(is.na(logL)) Inf else logL) + if(output.liks) return(liks[1:M+N,,drop=FALSE]) + else if(output.pi) return(pi) + else { + logL<--sum(log(comp[1:M+N])) + if(is.na(logL)) logL<-Inf + return(logL) + } } if(is.null(fixedQ)){ if(length(q.init)!=k) q.init<-rep(q.init[1],k) - if(opt.method=="optim") - fit<-optim(q.init,function(p) lik(makeQ(m,p,index.matrix),pi=pi), - method="L-BFGS-B",lower=rep(min.q,k)) - else if(opt.method=="none") + q.init<-if(logscale) log(q.init) else q.init + if(opt.method=="optim"){ + fit<-if(logscale) + optim(q.init,function(p) lik(makeQ(m,exp(p),index.matrix),pi=pi), + method="L-BFGS-B",lower=rep(log(min.q),k),upper=rep(log(max.q),k)) else + optim(q.init,function(p) lik(makeQ(m,p,index.matrix),pi=pi), + method="L-BFGS-B",lower=rep(min.q,k),upper=rep(max.q,k)) + } else if(opt.method=="none"){ fit<-list(objective=lik(makeQ(m,q.init,index.matrix),pi=pi), par=q.init) - else - fit<-nlminb(q.init,function(p) lik(makeQ(m,p,index.matrix),pi=pi), - lower=rep(0,k),upper=rep(1e50,k)) + } else { + fit<-if(logscale) + nlminb(q.init,function(p) lik(makeQ(m,exp(p),index.matrix),pi=pi), + lower=rep(log(min.q),k),upper=rep(log(max.q),k)) + else nlminb(q.init,function(p) lik(makeQ(m,p,index.matrix), + pi=pi),lower=rep(0,k),upper=rep(max.q,k)) + } + if(logscale) fit$par<-exp(fit$par) + if(pi[1]=="fitzjohn") pi<-setNames( + lik(makeQ(m,fit$par,index.matrix),FALSE,pi=pi,output.pi=TRUE), + states) obj<-list(logLik= if(opt.method=="optim") -fit$value else -fit$objective, rates=fit$par, index.matrix=index.matrix, states=states, pi=pi, - method=opt.method) + method=opt.method, + root.prior=root.prior) if(output.liks) obj$lik.anc<-lik(makeQ(m,obj$rates,index.matrix),TRUE, pi=pi) } else { fit<-lik(Q,pi=pi) + if(pi[1]=="fitzjohn") pi<-setNames(lik(Q,FALSE,pi=pi,output.pi=TRUE),states) obj<-list(logLik=-fit, rates=Q[sapply(1:k,function(x,y) which(x==y),index.matrix)], index.matrix=index.matrix, states=states, - pi=pi) + pi=pi, + root.prior=root.prior) if(output.liks) obj$lik.anc<-lik(makeQ(m,obj$rates,index.matrix),TRUE, pi=pi) } - lik.f<-function(q) -lik(q,output.liks=FALSE,pi=pi) + lik.f<-function(q) -lik(q,output.liks=FALSE, + pi=if(root.prior=="nuisance") "fitzjohn" else pi) obj$lik<-lik.f class(obj)<-"fitMk" return(obj) @@ -147,9 +192,12 @@ colnames(Q)<-rownames(Q)<-x$states print(round(Q,digits)) cat("\nFitted (or set) value of pi:\n") - print(x$pi) + print(round(x$pi,digits)) + cat(paste("due to treating the root prior as (a) ",x$root.prior,".\n", + sep="")) cat(paste("\nLog-likelihood:",round(x$logLik,digits),"\n")) - cat(paste("\nOptimization method used was \"",x$method,"\"\n\n",sep="")) + cat(paste("\nOptimization method used was \"",x$method,"\"\n\n", + sep="")) } ## summary method for objects of class "fitMk" @@ -186,9 +234,10 @@ plot.gfit<-function(x,...){ if("mkn"%in%class(x$lik)==FALSE){ stop("Sorry. No plot method presently available for objects of this type.") + object<-NULL } else { chk<-.check.pkg("geiger") - if(chk) plot(as.Qmatrix(x),...) + if(chk) object<-plot(as.Qmatrix(x),...) else { obj<-list() QQ<-.Qmatrix.from.gfit(x) @@ -199,9 +248,10 @@ obj$index.matrix[col(obj$index.matrix)!=row(obj$index.matrix)]<-1:k obj$rates<-QQ[sapply(1:k,function(x,y) which(x==y),obj$index.matrix)] class(obj)<-"fitMk" - plot(obj,...) + object<-plot(obj,...) } } + invisible(object) } ## S3 method for "Qmatrix" object class @@ -225,18 +275,34 @@ else mar<-c(1.1,1.1,3.1,1.1) if(hasArg(lwd)) lwd<-list(...)$lwd else lwd<-1 - spacer<-0.1 + if(hasArg(umbral)) umbral<-list(...)$umbral + else umbral<-FALSE + if(hasArg(ncat)) ncat<-list(...)$ncat + else ncat<-NULL + if(hasArg(spacer)) spacer<-list(...)$spacer + else spacer<-0.1 plot.new() par(mar=mar) xylim<-c(-1.2,1.2) plot.window(xlim=xylim,ylim=xylim,asp=1) if(!is.null(main)) title(main=main,cex.main=cex.main) nstates<-nrow(Q) - step<-360/nstates - angles<-seq(0,360-step,by=step)/180*pi - if(nstates==2) angles<-angles+pi/2 - v.x<-cos(angles) - v.y<-sin(angles) + if(!umbral||is.null(ncat)){ + step<-360/nstates + angles<-seq(0,360-step,by=step)/180*pi + if(nstates==2) angles<-angles+pi/2 + v.x<-cos(angles) + v.y<-sin(angles) + } else { + v.x<-v.y<-vector() + for(i in 1:length(ncat)){ + Q<-Q[sort(rownames(Q)),sort(colnames(Q))] + xp<--1+2*(i-1)/(length(ncat)-1) + v.x<-c(v.x,rep(xp,ncat[i])) + yp<-seq(1,-1,length.out=max(ncat))[1:ncat[i]] + v.y<-c(v.y,yp) + } + } for(i in 1:nstates) for(j in 1:nstates) if(if(!isSymmetric(Q)) i!=j else i>j){ dx<-v.x[j]-v.x[i] @@ -271,6 +337,8 @@ } text(v.x,v.y,rownames(Q),cex=cex.traits, col=make.transparent("black",0.7)) + object<-data.frame(states=rownames(Q),x=v.x,y=v.y) + invisible(object) } ## wraps around expm @@ -368,5 +436,8 @@ print.Qmatrix<-function(x,...){ cat("Estimated Q matrix:\n") - print(unclass(x)) + print(unclass(x),...) } + +is.Qmatrix<-function(x) "Qmatrix" %in% class(x) + diff -Nru r-cran-phytools-0.7-47/R/fitmultiMk.R r-cran-phytools-0.7-70/R/fitmultiMk.R --- r-cran-phytools-0.7-47/R/fitmultiMk.R 2020-06-01 04:05:56.000000000 +0000 +++ r-cran-phytools-0.7-70/R/fitmultiMk.R 2020-08-24 18:33:29.000000000 +0000 @@ -1,5 +1,5 @@ ## new function to fit multi-regime Mk model -## written by Liam J. Revell 2017 +## written by Liam J. Revell 2017, 2020 ## adapted from ape::ace by E. Paradis et al. 2013 fitmultiMk<-function(tree,x,model="ER",...){ @@ -103,6 +103,8 @@ regimes=regimes, pi=pi, method=opt.method) + lik.f<-function(q) -lik(q,output.liks=FALSE,pi=pi) + obj$lik<-lik.f class(obj)<-"fitmultiMk" return(obj) } diff -Nru r-cran-phytools-0.7-47/R/fitPagel.R r-cran-phytools-0.7-70/R/fitPagel.R --- r-cran-phytools-0.7-47/R/fitPagel.R 2020-05-28 18:13:18.000000000 +0000 +++ r-cran-phytools-0.7-70/R/fitPagel.R 2020-08-29 15:05:17.000000000 +0000 @@ -121,25 +121,27 @@ ## print method for objects of class "fitPagel" ## written by Liam J. Revell 2014, 2016 print.fitPagel<-function(x,...){ + if(hasArg(digits)) digits<-list(...)$digits + else digits<-6 cat("\nPagel's binary character correlation test:\n") cat(paste("\nAssumes \"",x$model, "\" substitution model for both characters\n",sep="")) cat("\nIndependent model rate matrix:\n") - print(x$independent.Q) + print(round(x$independent.Q,digits)) tmp<-if(x$dep.var=="xy") "x & y" else if(x$dep.var=="x") "x only" else if(x$dep.var=="y") "y only" cat(paste("\nDependent (",tmp,") model rate matrix:\n",sep="")) - print(x$dependent.Q) + print(round(x$dependent.Q,digits)) cat("\nModel fit:\n") obj<-matrix(c(x$independent.logL,x$dependent.logL, x$independent.AIC,x$dependent.AIC),2,2) rownames(obj)<-c("independent","dependent") colnames(obj)<-c("log-likelihood","AIC") - print(obj) + print(round(obj,digits)) cat("\nHypothesis test result:\n") - cat(paste(" likelihood-ratio: ",signif(x$lik.ratio,7),"\n")) - cat(paste(" p-value: ",signif(x$P,7),"\n")) + cat(paste(" likelihood-ratio: ",signif(x$lik.ratio,digits),"\n")) + cat(paste(" p-value: ",signif(x$P,digits),"\n")) cat(paste("\nModel fitting method used was",x$method,"\n\n")) } @@ -197,6 +199,8 @@ else cex.traits<-0.9 if(hasArg(cex.rates)) cex.rates<-list(...)$cex.rates else cex.rates<-0.8 + if(hasArg(lwd)) lwd<-list(...)$lwd + else lwd<-2 ## only used if lwd.by.rate=FALSE if(hasArg(lwd.by.rate)) lwd.by.rate<-list(...)$lwd.by.rate else lwd.by.rate<-FALSE if(lwd.by.rate){ @@ -207,7 +211,7 @@ else max.lwd<-10 LWD.ind[LWD.ind>max.lwd]<-max.lwd LWD.dep[LWD.dep>max.lwd]<-max.lwd - } else LWD.ind<-LWD.dep<-matrix(2,nrow(x$dependent.Q), + } else LWD.ind<-LWD.dep<-matrix(lwd,nrow(x$dependent.Q), ncol(x$dependent.Q)) par(mfrow=c(2,1)) ## INDEPENDENT MODEL diff -Nru r-cran-phytools-0.7-47/R/fitpolyMk.R r-cran-phytools-0.7-70/R/fitpolyMk.R --- r-cran-phytools-0.7-47/R/fitpolyMk.R 2020-06-01 04:05:34.000000000 +0000 +++ r-cran-phytools-0.7-70/R/fitpolyMk.R 2020-08-27 04:17:35.000000000 +0000 @@ -1,5 +1,6 @@ ## fitpolyMk -## written by Liam J. Revell 2019 +## fits several polymorphic discrete character evolution models +## written by Liam J. Revell 2019, 2020 fitpolyMk<-function(tree,x,model="SYM",ordered=FALSE,...){ if(hasArg(quiet)) quiet<-list(...)$quiet @@ -7,6 +8,13 @@ if(is.factor(x)) x<-setNames(as.character(x),names(x)) X<-strsplit(x,"+",fixed=TRUE) ns<-sapply(X,length) + if(ordered){ + if(hasArg(max.poly)) max.poly<-list(...)$max.poly + else if(hasArg(max.states)){ + max.states<-list(...)$max.states + max.poly<-max.states + } else max.poly<-max(ns) + } if(all(ns==1)){ cat("No polymorphic species found. Use fitMk.\n\n") object<-NULL @@ -17,9 +25,12 @@ states<-sort(unique(unlist(X))) if(ordered){ ss<-vector() - for(i in 1:(length(states)-1)) - ss<-c(ss,states[i],paste(states[i],states[i+1],sep="+")) - ss<-c(ss,states[i+1]) + for(i in 1:(length(states)-1)){ + ss<-c(ss,states[i]) + for(j in (i+1):length(states)) + if((j-i)j){ - dx<-v.x[j]-v.x[i] - dy<-v.y[j]-v.y[i] - slope<-abs(dy/dx) - shift.x<-0.02*sin(atan(dy/dx))*sign(j-i)*if(dy/dx>0) 1 else -1 - shift.y<-0.02*cos(atan(dy/dx))*sign(j-i)*if(dy/dx>0) -1 else 1 - s<-c(v.x[i]+spacer*cos(atan(slope))*sign(dx)+ - if(isSymmetric(Q)) 0 else shift.x, - v.y[i]+spacer*sin(atan(slope))*sign(dy)+ - if(isSymmetric(Q)) 0 else shift.y) - e<-c(v.x[j]+spacer*cos(atan(slope))*sign(-dx)+ - if(isSymmetric(Q)) 0 else shift.x, - v.y[j]+spacer*sin(atan(slope))*sign(-dy)+ - if(isSymmetric(Q)) 0 else shift.y) - if(x$index.matrix[i,j]!=0){ - if(abs(diff(c(i,j)))==1||abs(diff(c(i,j)))==(nstates-1)) + if(attr(x$ordered,"max.poly")==2){ + step<-360/nstates + angles<-seq(-floor(nstates/2)*step,360-ceiling(nstates/2)*step,by=step)/180*pi + if(nstates==2) angles<-angles+pi/2 + v.x<-sin(angles) + v.y<-cos(angles) + for(i in 1:nstates) for(j in 1:nstates){ + if(if(!isSymmetric(Q)) i!=j else i>j){ + dx<-v.x[j]-v.x[i] + dy<-v.y[j]-v.y[i] + slope<-abs(dy/dx) + shift.x<-0.02*sin(atan(dy/dx))*sign(j-i)*if(dy/dx>0) 1 else -1 + shift.y<-0.02*cos(atan(dy/dx))*sign(j-i)*if(dy/dx>0) -1 else 1 + s<-c(v.x[i]+spacer*cos(atan(slope))*sign(dx)+ + if(isSymmetric(Q)) 0 else shift.x, + v.y[i]+spacer*sin(atan(slope))*sign(dy)+ + if(isSymmetric(Q)) 0 else shift.y) + e<-c(v.x[j]+spacer*cos(atan(slope))*sign(-dx)+ + if(isSymmetric(Q)) 0 else shift.x, + v.y[j]+spacer*sin(atan(slope))*sign(-dy)+ + if(isSymmetric(Q)) 0 else shift.y) + if(x$index.matrix[i,j]!=0){ + if(abs(diff(c(i,j)))==1||abs(diff(c(i,j)))==(nstates-1)) + text(mean(c(s[1],e[1]))+1.5*shift.x, + mean(c(s[2],e[2]))+1.5*shift.y, + round(Q[i,j],signif),cex=cex.rates, + srt=atan(dy/dx)*180/pi) + else + text(mean(c(s[1],e[1]))+0.3*diff(c(s[1],e[1]))+ + 1.5*shift.x, + mean(c(s[2],e[2]))+0.3*diff(c(s[2],e[2]))+ + 1.5*shift.y, + round(Q[i,j],signif),cex=cex.rates, + srt=atan(dy/dx)*180/pi) + arrows(s[1],s[2],e[1],e[2],length=0.05, + code=if(isSymmetric(Q)) 3 else 2,lwd=lwd) + } + } + } + text(v.x,v.y,x$states,cex=cex.traits, + col=make.transparent(par()$fg,0.7)) + } else { + nlevs<-attr(x$ordered,"max.poly") + Ns<-inv.ncombn2(nstates,nlevs) + v.x<-v.y<-vector() + xx<-seq(-1,1,length.out=Ns) + for(i in 1:Ns){ + for(j in 1:min(nlevs,Ns-i+1)){ + v.x<-c(v.x,mean(xx[i:(i+j-1)])) + v.y<-c(v.y,1-2*(j-1)/(nlevs-1)) + } + } + for(i in 1:nstates) for(j in 1:nstates){ + if(if(!isSymmetric(Q)) i!=j else i>j){ + dx<-v.x[j]-v.x[i] + dy<-v.y[j]-v.y[i] + slope<-abs(dy/dx) + shift.x<-0.02*sin(atan(dy/dx))*sign(j-i)*if(dy/dx>0) 1 else -1 + shift.y<-0.02*cos(atan(dy/dx))*sign(j-i)*if(dy/dx>0) -1 else 1 + s<-c(v.x[i]+spacer*cos(atan(slope))*sign(dx)+ + if(isSymmetric(Q)) 0 else shift.x, + v.y[i]+spacer*sin(atan(slope))*sign(dy)+ + if(isSymmetric(Q)) 0 else shift.y) + e<-c(v.x[j]+spacer*cos(atan(slope))*sign(-dx)+ + if(isSymmetric(Q)) 0 else shift.x, + v.y[j]+spacer*sin(atan(slope))*sign(-dy)+ + if(isSymmetric(Q)) 0 else shift.y) + if(x$index.matrix[i,j]!=0){ text(mean(c(s[1],e[1]))+1.5*shift.x, mean(c(s[2],e[2]))+1.5*shift.y, round(Q[i,j],signif),cex=cex.rates, srt=atan(dy/dx)*180/pi) - else - text(mean(c(s[1],e[1]))+0.3*diff(c(s[1],e[1]))+ - 1.5*shift.x, - mean(c(s[2],e[2]))+0.3*diff(c(s[2],e[2]))+ - 1.5*shift.y, - round(Q[i,j],signif),cex=cex.rates, - srt=atan(dy/dx)*180/pi) - arrows(s[1],s[2],e[1],e[2],length=0.05, - code=if(isSymmetric(Q)) 3 else 2,lwd=lwd) + arrows(s[1],s[2],e[1],e[2],length=0.05, + code=if(isSymmetric(Q)) 3 else 2,lwd=lwd) + } } } - text(v.x,v.y,x$states,cex=cex.traits, - col=make.transparent("black",0.7)) + text(v.x,v.y,x$states,cex=cex.traits, + col=make.transparent(par()$fg,0.7)) + } } else { Ns<-inv.ncombn(nstates) step.y<-2/(Ns-1) @@ -215,7 +271,7 @@ } } text(v.x,v.y,x$states,cex=cex.traits, - col=make.transparent("black",0.7)) + col=make.transparent(par()$fg,0.7)) } } @@ -233,3 +289,196 @@ } return(n-1) } +inv.ncombn2<-function(N,m){ + n<-2 + Nc<-0 + while(Nc!=N){ + Nc<-sum((n:1)[1:min(m,n)]) + n<-n+1 + } + return(n-1) +} + +## maybe this should be plot.mkModel or something? +## then you create a model class & plot it? + +graph.polyMk<-function(k=2,model="SYM",ordered=FALSE,...){ + if(hasArg(states)) states<-list(...)$states + else states<-0:(k-1) + if(hasArg(max.poly)) max.poly<-list(...)$max.poly + else max.poly<-k + if(hasArg(main)) main<-list(...)$main + else main<-NULL + if(hasArg(cex.main)) cex.main<-list(...)$cex.main + else cex.main<-1.2 + if(hasArg(cex.traits)) cex.traits<-list(...)$cex.traits + else cex.traits<-1 + if(hasArg(mar)) mar<-list(...)$mar + else mar<-c(1.1,1.1,3.1,1.1) + if(hasArg(lwd)) lwd<-list(...)$lwd + else lwd<-1 + if(hasArg(quiet)) quiet<-list(...)$quiet + else quiet<-FALSE + if(hasArg(xlim)) xlim<-list(...)$xlim + else xlim<-c(-1.2,1.2) + if(hasArg(ylim)) ylim<-list(...)$ylim + else ylim<-c(-1.2,1.2) + if(hasArg(plot)) plot<-list(...)$plot + else plot<-TRUE + if(hasArg(asp)) asp<-list(...)$asp + else asp<-1 + if(ordered){ + ss<-vector() + for(i in 1:(length(states)-1)){ + ss<-c(ss,states[i]) + for(j in (i+1):length(states)) + if((j-i)0 + &&(0%in%c(length(SDij),length(SDji))) + &&(1%in%c(length(SDij),length(SDji)))){ + if(model=="ER"){ + tmodel[i,j]<-tmodel[j,i]<-1 + } else if(model%in%c("ARD","SYM")){ + index<-index+1 + tmodel[i,j]<-index + if(model=="SYM") tmodel[j,i]<-index + else { + tmodel[j,i]<-index+1 + index<-index+1 + } + } else if(model=="transient"){ + if(length(poly[[i]])>length(poly[[j]])){ + tmodel[i,j]<-1 + tmodel[j,i]<-2 + } else { + tmodel[i,j]<-2 + tmodel[j,i]<-1 + } + } + } + } + } + if(plot){ + spacer<-0.1 + plot.new() + par(mar=mar) + plot.window(xlim=xlim,ylim=ylim,asp=asp) + if(!is.null(main)) title(main=main,cex.main=cex.main) + if(ordered){ + if(max.poly==2){ + step<-360/nstates + angles<-seq(-floor(nstates/2)*step,360-ceiling(nstates/2)*step,by=step)/180*pi + if(k==2) angles<-angles+pi/2 + v.x<-sin(angles) + v.y<-cos(angles) + for(i in 1:nstates) for(j in 1:nstates){ + if(if(!isSymmetric(tmodel)) i!=j else i>j){ + dx<-v.x[j]-v.x[i] + dy<-v.y[j]-v.y[i] + slope<-abs(dy/dx) + shift.x<-0.02*sin(atan(dy/dx))*sign(j-i)*if(dy/dx>0) 1 else -1 + shift.y<-0.02*cos(atan(dy/dx))*sign(j-i)*if(dy/dx>0) -1 else 1 + s<-c(v.x[i]+spacer*cos(atan(slope))*sign(dx)+ + if(isSymmetric(tmodel)) 0 else shift.x, + v.y[i]+spacer*sin(atan(slope))*sign(dy)+ + if(isSymmetric(tmodel)) 0 else shift.y) + e<-c(v.x[j]+spacer*cos(atan(slope))*sign(-dx)+ + if(isSymmetric(tmodel)) 0 else shift.x, + v.y[j]+spacer*sin(atan(slope))*sign(-dy)+ + if(isSymmetric(tmodel)) 0 else shift.y) + if(tmodel[i,j]!=0){ + arrows(s[1],s[2],e[1],e[2],length=0.05, + code=if(isSymmetric(tmodel)) 3 else 2,lwd=lwd) + } + } + } + text(v.x,v.y,colnames(tmodel),cex=cex.traits, + col=make.transparent(par()$fg,0.7)) + } else { + nlevs<-max.poly + Ns<-inv.ncombn2(nstates,nlevs) + v.x<-v.y<-vector() + xx<-seq(-1,1,length.out=Ns) + for(i in 1:k){ + for(j in 1:min(nlevs,Ns-i+1)){ + v.x<-c(v.x,mean(xx[i:(i+j-1)])) + v.y<-c(v.y,1-2*(j-1)/(nlevs-1)) + } + } + for(i in 1:nstates) for(j in 1:nstates){ + if(if(!isSymmetric(tmodel)) i!=j else i>j){ + dx<-v.x[j]-v.x[i] + dy<-v.y[j]-v.y[i] + slope<-abs(dy/dx) + shift.x<-0.02*sin(atan(dy/dx))*sign(j-i)*if(dy/dx>0) 1 else -1 + shift.y<-0.02*cos(atan(dy/dx))*sign(j-i)*if(dy/dx>0) -1 else 1 + s<-c(v.x[i]+spacer*cos(atan(slope))*sign(dx)+ + if(isSymmetric(tmodel)) 0 else shift.x, + v.y[i]+spacer*sin(atan(slope))*sign(dy)+ + if(isSymmetric(tmodel)) 0 else shift.y) + e<-c(v.x[j]+spacer*cos(atan(slope))*sign(-dx)+ + if(isSymmetric(tmodel)) 0 else shift.x, + v.y[j]+spacer*sin(atan(slope))*sign(-dy)+ + if(isSymmetric(tmodel)) 0 else shift.y) + if(tmodel[i,j]!=0){ + arrows(s[1],s[2],e[1],e[2],length=0.05, + code=if(isSymmetric(tmodel)) 3 else 2,lwd=lwd) + } + } + } + text(v.x,v.y,rownames(tmodel),cex=cex.traits, + col=make.transparent(par()$fg,0.7)) + } + } else { + step.y<-2/(k-1) + v.x<-v.y<-vector() + for(i in 1:k){ + nc<-ncombn(k,i) + v.x<-c(v.x,if(nc>1) seq(-1,1,by=2/(nc-1)) else 0) + v.y<-c(v.y,rep(1-rep((i-1)*step.y,nc))) + } + for(i in 1:ncol(tmodel)) for(j in 1:ncol(tmodel)){ + if(if(!isSymmetric(tmodel)) i!=j else i>j){ + dx<-v.x[j]-v.x[i] + dy<-v.y[j]-v.y[i] + slope<-abs(dy/dx) + shift.x<-0.02*sin(atan(dy/dx))*sign(j-i)*if(dy/dx>0) 1 else -1 + shift.y<-0.02*cos(atan(dy/dx))*sign(j-i)*if(dy/dx>0) -1 else 1 + s<-c(v.x[i]+spacer*cos(atan(slope))*sign(dx)+ + if(isSymmetric(tmodel)) 0 else shift.x, + v.y[i]+spacer*sin(atan(slope))*sign(dy)+ + if(isSymmetric(tmodel)) 0 else shift.y) + e<-c(v.x[j]+spacer*cos(atan(slope))*sign(-dx)+ + if(isSymmetric(tmodel)) 0 else shift.x, + v.y[j]+spacer*sin(atan(slope))*sign(-dy)+ + if(isSymmetric(tmodel)) 0 else shift.y) + if(tmodel[i,j]!=0){ + arrows(s[1],s[2],e[1],e[2],length=0.05, + code=if(isSymmetric(tmodel)) 3 else 2,lwd=lwd) + } + } + } + text(v.x,v.y,rownames(tmodel),cex=cex.traits, + col=make.transparent(par()$fg,0.7)) + } + } + invisible(tmodel) +} \ No newline at end of file diff -Nru r-cran-phytools-0.7-47/R/ltt.R r-cran-phytools-0.7-70/R/ltt.R --- r-cran-phytools-0.7-47/R/ltt.R 2020-06-01 05:31:09.000000000 +0000 +++ r-cran-phytools-0.7-70/R/ltt.R 2020-08-20 19:59:05.000000000 +0000 @@ -303,6 +303,12 @@ if(hasArg(main)) main<-list(...)$main else main=expression(paste("null distribution of ", gamma)) + if(hasArg(las)) las<-list(...)$las + else las<-par()$las + if(hasArg(cex.axis)) cex.axis<-list(...)$cex.axis + else cex.axis<-par()$cex.axis + if(hasArg(cex.lab)) cex.lab<-list(...)$cex.lab + else cex.lab<-par()$cex.lab hh<-hist(x$null.gamma,breaks=min(c(max(12, round(length(x$null.gamma)/10)),20)), plot=FALSE) @@ -310,7 +316,7 @@ else ylim<-c(0,1.15*max(hh$counts)) plot(hh,xlim=range(c(x$gamma,x$null.gamma)), main=main,xlab=expression(gamma),col="lightgrey", - ylim=ylim) + ylim=ylim,las=las,cex.axis=cex.axis,cex.lab=cex.lab) arrows(x0=x$gamma,y0=par()$usr[4],y1=0,length=0.12, col=make.transparent("blue",0.5),lwd=2) text(x$gamma,0.96*par()$usr[4], diff -Nru r-cran-phytools-0.7-47/R/make.simmap.R r-cran-phytools-0.7-70/R/make.simmap.R --- r-cran-phytools-0.7-47/R/make.simmap.R 2020-05-31 21:00:23.000000000 +0000 +++ r-cran-phytools-0.7-70/R/make.simmap.R 2020-09-18 02:10:22.000000000 +0000 @@ -100,7 +100,7 @@ # written by Liam J. Revell 2013 printmessage<-function(Q,pi,method){ if(method=="empirical"||method=="fixed") - cat("make.simmap is sampling character histories conditioned on the transition matrix\n\nQ =\n") + cat("make.simmap is sampling character histories conditioned on\nthe transition matrix\n\nQ =\n") else if(method=="mcmc"){ cat("make.simmap is simulating with a sample of Q from\nthe posterior distribution\n") cat("\nMean Q from the posterior is\nQ =\n") @@ -303,7 +303,7 @@ if(N>printlen) cat(paste("\t",paste(x$tip.label[1:printlen],collapse=", "),", ...\n",sep="")) else print(x$tip.label) ss<-sort(unique(c(getStates(x,"tips"),getStates(x,"nodes")))) - cat(paste("\nThe tree includes a mapped, ",length(ss),"-state discrete character with states:\n", + cat(paste("\nThe tree includes a mapped, ",length(ss),"-state discrete character\nwith states:\n", sep="")) if(length(ss)>printlen) cat(paste("\t",paste(ss[1:printlen],collapse=", "),", ...\n",sep="")) else cat(paste("\t",paste(ss,collapse=", "),"\n",sep="")) diff -Nru r-cran-phytools-0.7-47/R/phenogram.R r-cran-phytools-0.7-70/R/phenogram.R --- r-cran-phytools-0.7-47/R/phenogram.R 2017-09-19 16:57:24.000000000 +0000 +++ r-cran-phytools-0.7-70/R/phenogram.R 2020-08-21 00:38:42.000000000 +0000 @@ -1,5 +1,5 @@ ## function creates a phenogram (i.e., 'traitgram') -## written by Liam J. Revell 2011, 2012, 2013, 2014, 2015, 2016 +## written by Liam J. Revell 2011, 2012, 2013, 2014, 2015, 2016, 2020 phenogram<-function(tree,x,fsize=1.0,ftype="reg",colors=NULL,axes=list(),add=FALSE,...){ ## get optional arguments @@ -48,6 +48,12 @@ else quiet<-FALSE if(hasArg(label.pos)) label.pos<-list(...)$label.pos else label.pos<-NULL + if(hasArg(cex.axis)) cex.axis<-list(...)$cex.axis + else cex.axis<-par()$cex.axis + if(hasArg(cex.lab)) cex.lab<-list(...)$cex.lab + else cex.lab<-par()$cex.lab + if(hasArg(las)) las<-list(...)$las + else las<-par()$las ## end optional arguments # check tree if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".") @@ -68,8 +74,10 @@ if(!is.null(axes$time)) xlim<-axes$time if(!add&&is.null(xlim)){ pp<-par("pin")[1] - sw<-fsize*(max(strwidth(tree$tip.label,units="inches")))+offsetFudge*offset*fsize*strwidth("W",units="inches") - alp<-optimize(function(a,H,link,sw,pp) (a*1.04*(max(H)+link)+sw-pp)^2,H=H,link=link,sw=sw,pp=pp,interval=c(0,1e6))$minimum + sw<-fsize*(max(strwidth(tree$tip.label,units="inches")))+ + offsetFudge*offset*fsize*strwidth("W",units="inches") + alp<-optimize(function(a,H,link,sw,pp) (a*1.04*(max(H)+link)+sw-pp)^2,H=H, + link=link,sw=sw,pp=pp,interval=c(0,1e6))$minimum xlim<-c(min(H),max(H)+link+sw/alp) } if(!quiet&&Ntip(tree)>=40&&spread.labels){ @@ -82,13 +90,18 @@ if(is.null(tree$maps)){ if(is.null(colors)) colors<-"black" if(!add){ - plot(H[1,],X[1,],type=type,lwd=lwd,lty=lty,col=colors,xlim=xlim,ylim=ylim,log=log,asp=asp,xlab="",ylab="",frame=FALSE, axes=FALSE) - if(spread.labels) tt<-spreadlabels(tree,x,fsize=fsize,cost=spread.cost,range=spread.range,label.pos=label.pos) else tt<-x[1:length(tree$tip)] + plot(H[1,],X[1,],type=type,lwd=lwd,lty=lty,col=colors,xlim=xlim,ylim=ylim, + log=log,asp=asp,xlab="",ylab="",frame=FALSE,axes=FALSE) + if(spread.labels) tt<-spreadlabels(tree,x,fsize=fsize,cost=spread.cost, + range=spread.range,label.pos=label.pos) else tt<-x[1:length(tree$tip)] if(tree$edge[1,2]<=length(tree$tip)){ if(fsize&&!add){ - text(gsub("_"," ",tree$tip.label[tree$edge[1,2]]),x=H[1,2]+link,y=tt[tree$edge[1,2]],cex=fsize,font=ftype,pos=4,offset=offset) - tip.coords[tree$tip.label[tree$edge[1,2]],]<-c(H[1,2]+link,tt[tree$edge[1,2]]) - if(link>0) lines(x=c(H[1,2],H[1,2]+link),y=c(X[1,2],tt[tree$edge[1,2]]),lty=3) + text(gsub("_"," ",tree$tip.label[tree$edge[1,2]]),x=H[1,2]+link, + y=tt[tree$edge[1,2]],cex=fsize,font=ftype,pos=4,offset=offset) + tip.coords[tree$tip.label[tree$edge[1,2]],]<-c(H[1,2]+link, + tt[tree$edge[1,2]]) + if(link>0) lines(x=c(H[1,2],H[1,2]+link),y=c(X[1,2], + tt[tree$edge[1,2]]),lty=3) } } s<-2 @@ -97,9 +110,12 @@ lines(H[i,],X[i,],type=type,lwd=lwd,lty=lty,col=colors) if(tree$edge[i,2]<=length(tree$tip)){ if(fsize&&!add){ - text(gsub("_"," ",tree$tip.label[tree$edge[i,2]]),x=H[i,2]+link,y=tt[tree$edge[i,2]],cex=fsize,font=ftype,pos=4,offset=offset) - tip.coords[tree$tip.label[tree$edge[i,2]],]<-c(H[i,2]+link,tt[tree$edge[i,2]]) - if(link>0) lines(x=c(H[i,2],H[i,2]+link),y=c(X[i,2],tt[tree$edge[i,2]]),lty=3) + text(gsub("_"," ",tree$tip.label[tree$edge[i,2]]),x=H[i,2]+link, + y=tt[tree$edge[i,2]],cex=fsize,font=ftype,pos=4,offset=offset) + tip.coords[tree$tip.label[tree$edge[i,2]],]<-c(H[i,2]+link, + tt[tree$edge[i,2]]) + if(link>0) lines(x=c(H[i,2],H[i,2]+link),y=c(X[i,2],tt[tree$edge[i,2]]), + lty=3) } } } @@ -115,24 +131,33 @@ a<-c(y,y+tree$maps[[i]][j]) b<-m*(a-H[i,1])+X[i,1] if(i==1&&j==1&&!add) { - plot(a,b,col=colors[names(tree$maps[[i]])[j]],type=type,lwd=lwd,lty=lty,xlim=xlim,ylim=ylim,log=log,asp=asp,axes=FALSE,xlab="",ylab="") - if(spread.labels) tt<-spreadlabels(tree,x[1:length(tree$tip)],fsize=fsize,cost=spread.cost,range=spread.range) else tt<-x[1:length(tree$tip)] - print(tt) - } else lines(a,b,col=colors[names(tree$maps[[i]])[j]],lwd=lwd,lty=lty,type=type) + plot(a,b,col=colors[names(tree$maps[[i]])[j]],type=type,lwd=lwd, + lty=lty,xlim=xlim,ylim=ylim,log=log,asp=asp,axes=FALSE,xlab="", + ylab="") + if(spread.labels) tt<-spreadlabels(tree,x[1:length(tree$tip)], + fsize=fsize,cost=spread.cost,range=spread.range) else + tt<-x[1:length(tree$tip)] + } else lines(a,b,col=colors[names(tree$maps[[i]])[j]],lwd=lwd,lty=lty, + type=type) y<-a[2] } if(tree$edge[i,2]<=length(tree$tip)){ if(fsize&&!add){ - text(gsub("_"," ",tree$tip.label[tree$edge[i,2]]),x=H[i,2]+link,y=tt[tree$edge[i,2]],cex=fsize,font=ftype,pos=4,offset=offset) - tip.coords[tree$tip.label[tree$edge[i,2]],]<-c(H[i,2]+link,tt[tree$edge[i,2]]) - if(link>0) lines(x=c(H[i,2],H[i,2]+link),y=c(X[i,2],tt[tree$edge[i,2]]),lty=3) + text(gsub("_"," ",tree$tip.label[tree$edge[i,2]]),x=H[i,2]+link, + y=tt[tree$edge[i,2]],cex=fsize,font=ftype,pos=4,offset=offset) + tip.coords[tree$tip.label[tree$edge[i,2]],]<-c(H[i,2]+link, + tt[tree$edge[i,2]]) + if(link>0) lines(x=c(H[i,2],H[i,2]+link),y=c(X[i,2], + tt[tree$edge[i,2]]),lty=3) } } } } if(!add){ at<-round(0:(nticks-1)*max(H)/(nticks-1),digits) - axis(1,at=at); axis(2); title(xlab=xlab,ylab=ylab,main=main,sub=sub) + axis(1,at=at,cex.axis=cex.axis,cex.lab=cex.lab,las=las) + axis(2,cex.axis=cex.axis,cex.lab=cex.lab,las=las) + title(xlab=xlab,ylab=ylab,main=main,sub=sub) } if(hold) null<-dev.flush() xx<-setNames(c(H[1,1],H[,2]),c(tree$edge[1,1],tree$edge[,2])) @@ -156,7 +181,8 @@ else { if(is.null(range)) range<-range(x) yy<-x[1:Ntip(tree)] - zz<-setNames((rank(yy,ties.method="random")-1)/(length(yy)-1)*diff(range(yy))+range(yy)[1],names(yy)) + zz<-setNames((rank(yy,ties.method="random")-1)/(length(yy)-1)*diff(range(yy))+ + range(yy)[1],names(yy)) mm<-max(fsize*strheight(tree$tip.label)) ff<-function(zz,yy,cost,mo=1,ms=1){ ZZ<-cbind(zz-mm/2,zz+mm/2) @@ -173,7 +199,8 @@ ms<-ff(yy,zz,cost=c(0,1)) if(mo==0&&ms==0) return(yy) else { - rr<-optim(zz,ff,yy=yy,mo=mo,ms=ms,cost=cost,method="L-BFGS-B",lower=rep(range[1],length(yy)),upper=rep(range[2],length(yy))) + rr<-optim(zz,ff,yy=yy,mo=mo,ms=ms,cost=cost,method="L-BFGS-B", + lower=rep(range[1],length(yy)),upper=rep(range[2],length(yy))) return(rr$par) } } diff -Nru r-cran-phytools-0.7-47/R/phylo.heatmap.R r-cran-phytools-0.7-70/R/phylo.heatmap.R --- r-cran-phytools-0.7-47/R/phylo.heatmap.R 2020-06-01 01:29:38.000000000 +0000 +++ r-cran-phytools-0.7-70/R/phylo.heatmap.R 2020-07-09 19:02:50.000000000 +0000 @@ -48,12 +48,17 @@ segments(x, y[1], x, y[length(y)]) segments(x[1], y, x[length(x)], y) } - if(legend) add.color.bar(leg=END-START,cols=colors,lims=range(X,na.rm=TRUE), - title=if(standardize) "standardized value" else "value", - subtitle=if(standardize) "SD units" else "",prompt=FALSE,x=START, - y=-1/(2*(Ntip(cw)-1))-3*fsize[3]*strheight("W"), - digits=if(max(abs(X),na.rm=TRUE)<1) round(log10(1/max(abs(X),na.rm=TRUE)))+1 - else 2,fsize=fsize[3]) + if(legend){ + if(hasArg(leg.title)) leg.title<-list(...)$leg.title + else leg.title<-if(standardize) "standardized value" else "value" + if(hasArg(leg.subtitle)) leg.subtitle<-list(...)$leg.subtitle + else leg.subtitle<-if(standardize) "SD units" else "" + add.color.bar(leg=END-START,cols=colors,lims=range(X,na.rm=TRUE), + title=leg.title,subtitle=leg.subtitle,prompt=FALSE,x=START, + y=-1/(2*(Ntip(cw)-1))-3*fsize[3]*strheight("W"), + digits=if(max(abs(X),na.rm=TRUE)<1) round(log10(1/max(abs(X),na.rm=TRUE)))+1 + else 2,fsize=fsize[3]) + } if(labels) text(x=seq(START,END,by=(END-START)/(ncol(X)-1)), y=rep(1+1/(2*(Ntip(cw)-1))+0.4*fsize[2]*strwidth("I"),ncol(X)), colnames(X),srt=70,adj=c(0,0.5),cex=fsize[2]) diff -Nru r-cran-phytools-0.7-47/R/phylomorphospace.R r-cran-phytools-0.7-70/R/phylomorphospace.R --- r-cran-phytools-0.7-47/R/phylomorphospace.R 2020-05-31 21:27:56.000000000 +0000 +++ r-cran-phytools-0.7-70/R/phylomorphospace.R 2020-08-20 14:28:09.000000000 +0000 @@ -1,5 +1,5 @@ ## this funciton creates a phylomorphospace plot (Sidlauskas 2006) -## written by Liam J. Revell 2010-13, 2015, 2018 +## written by Liam J. Revell 2010-13, 2015, 2018, 2020 phylomorphospace<-function(tree,X,A=NULL,label=c("radial","horizontal","off"),control=list(),...){ # some minor error checking @@ -59,6 +59,8 @@ else pch<-21 if(hasArg(bty)) bty<-list(...)$bty else bty<-"o" + if(hasArg(las)) las<-list(...)$las + else las<-par()$las # deprecate to logical label argument label<-label[1] if(label==TRUE||label==FALSE) @@ -72,7 +74,7 @@ YY<-matrix(bb[as.character(tree$edge)],nrow(tree$edge),2) # plot projection if(!add) plot(x=A[1,1],y=A[1,2],xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,pch=16,cex=0.1, - col="white",axes=axes,frame.plot=TRUE,bty=bty) + col="white",axes=axes,frame.plot=TRUE,bty=bty,las=las) if(is.null(tree$maps)){ for(i in 1:nrow(XX)) lines(XX[i,],YY[i,],col=con$col.edge[as.character(tree$edge[i,2])], lwd=lwd) diff -Nru r-cran-phytools-0.7-47/R/phylosig.R r-cran-phytools-0.7-70/R/phylosig.R --- r-cran-phytools-0.7-47/R/phylosig.R 2019-12-13 14:31:56.000000000 +0000 +++ r-cran-phytools-0.7-70/R/phylosig.R 2020-09-14 18:49:14.000000000 +0000 @@ -1,5 +1,6 @@ -## function for computing phylogenetic signal by the lambda (Pagel 1999) of K (Blomberg et al. 2003) methods -## written by Liam J. Revell 2011/2012, 2019 +## function for computing phylogenetic signal by the lambda (Pagel 1999) +## or K (Blomberg et al. 2003) methods +## written by Liam J. Revell 2011/2012, 2019, 2020 phylosig<-function(tree,x,method="K",test=FALSE,nsim=1000,se=NULL,start=NULL,control=list()){ # some minor error checking @@ -229,11 +230,19 @@ if(attr(x,"method")=="K"&&attr(x,"test")) "K" else "sig2" if(hasArg(res)) res<-list(...)$res else res<-100 + if(hasArg(las)) las<-list(...)$las + else las<-par()$las + if(hasArg(cex.lab)) cex.lab<-list(...)$cex.lab*par()$cex + else cex.lab<-par()$cex.lab + if(hasArg(cex.axis)) cex.axis<-list(...)$cex.axis*par()$cex + else cex.axis<-par()$cex.axis + if(hasArg(bty)) bty<-list(...)$bty + else bty<-par()$bty if(attr(x,"method")=="lambda"){ lambda<-seq(0,max(c(1,x$lambda)),length.out=res) logL<-sapply(lambda,x$lik) plot(lambda,logL,xlab=expression(lambda),ylab="log(L)", - type="l",bty="l") + type="l",bty=bty,las=las,cex.lab=cex.lab,cex.axis=cex.axis) lines(rep(x$lambda,2),c(par()$usr[3],x$logL),lty="dotted") text(x=x$lambda+0.01*diff(par()$usr[1:2]), par()$usr[3]+0.5*diff(par()$usr[3:4]), @@ -256,12 +265,13 @@ cat("Sorry. This is not a valid plotting option for your object.\n\n") else { hist(x$sim.K,breaks=min(c(max(12,round(length(x$sim.K)/10)), - 20)),bty="l",col="lightgrey",border="lightgrey", - main="",xlab="K",ylab="null distribution of K") + 20)),bty=bty,col="lightgrey",border="lightgrey", + main="",xlab="K",ylab="null distribution of K", + las=las,cex.lab=cex.lab,cex.axis=cex.axis) arrows(x0=x$K,y0=par()$usr[4],y1=0,length=0.12, col=make.transparent("blue",0.5),lwd=2) text(x$K,0.95*par()$usr[4],"observed value of K", - pos=if(x$K>mean(x$sim.K)) 2 else 4) + pos=if(x$K>mean(range(x$sim.K))) 2 else 4) } } else if(what=="sig2"){ if(attr(x,"se")==FALSE) @@ -270,7 +280,8 @@ sig2<-seq(0.25*x$sig2,1.75*x$sig2,length.out=res) logL<-sapply(sig2,x$lik) plot(sig2,logL,xlab=expression(sigma^2),ylab="log(L)", - type="l",bty="l") + type="l",bty=bty,las=las,cex.lab=cex.lab, + cex.axis=cex.axis) lines(rep(x$sig2,2),c(par()$usr[3],x$logL),lty="dotted") text(x=x$sig2+0.01*diff(par()$usr[1:2]), par()$usr[3]+0.5*diff(par()$usr[3:4]), diff -Nru r-cran-phytools-0.7-47/R/phylo.to.map.R r-cran-phytools-0.7-70/R/phylo.to.map.R --- r-cran-phytools-0.7-47/R/phylo.to.map.R 2019-05-21 18:11:46.000000000 +0000 +++ r-cran-phytools-0.7-70/R/phylo.to.map.R 2020-07-02 16:51:25.000000000 +0000 @@ -32,7 +32,7 @@ } ## S3 method to plot object of class "phylo.to.map" -## written by Liam J. Revell 2013, 2014, 2016, 2019 +## written by Liam J. Revell 2013, 2014, 2016, 2019, 2020 plot.phylo.to.map<-function(x,type=c("phylogram","direct"),...){ type<-type[1] @@ -89,6 +89,8 @@ else lty<-"dashed" if(hasArg(pts)) pts<-list(...)$pts else pts<-TRUE + if(hasArg(col.edge)) col.edge<-list(...)$col.edge + else col.edge<-rep(par()$fg,nrow(x$tree$edge)) if(type=="phylogram"){ if(x$direction=="downwards"&&direction=="rightwards"){ cat("\"phylo.to.map\" direction is \"downwards\" but plot direction has been given as \"rightwards\".\n") @@ -163,9 +165,10 @@ for(i in 1:nrow(Y)) lines(rep(x[cw$edge[i,2]],2),Y[i,], lwd=lwd[1],lend=2) # plot horizontal relationships - for(i in 1:cw$Nnode+n) - lines(range(x[cw$edge[which(cw$edge[,1]==i),2]]), - Y[which(cw$edge[,1]==i),1],lwd=lwd[1],lend=2) + for(i in 1:cw$Nnode+n){ + ii<-which(cw$edge[,1]==i) + if(length(ii)>1) lines(range(x[cw$edge[ii,2]]),Y[ii,1],lwd=lwd[1],lend=2) + } # plot tip labels for(i in 1:n){ if(ftype) text(x[i],Y[which(cw$edge[,2]==i),2], @@ -211,8 +214,10 @@ } points(coords,pch=pch,cex=cex.points[2],bg=colors[,2]) for(i in 1:nrow(X)) lines(X[i,],rep(y[cw$edge[i,2]],2),lwd=lwd[1],lend=2) - for(i in 1:cw$Nnode+n) lines(X[which(cw$edge[,1]==i),1], - range(y[cw$edge[which(cw$edge[,1]==i),2]]),lwd=lwd[1],lend=2) + for(i in 1:cw$Nnode+n){ + ii<-which(cw$edge[,1]==i) + if(length(ii)>1) lines(X[ii,1],range(y[cw$edge[ii,2]]),lwd=lwd[1],lend=2) + } for(i in 1:n){ if(ftype) text(X[which(cw$edge[,2]==i),2],y[i], paste(" ",sub("_"," ",cw$tip.label[i]),sep=""), @@ -240,7 +245,7 @@ ## rotates all nodes to try and match tip an ordering -## written by Liam J. Revell 2013, 2015 +## written by Liam J. Revell 2013, 2015, 2020 minRotate<-function(tree,x,...){ if(hasArg(print)) print<-list(...)$print else print<-TRUE @@ -248,7 +253,7 @@ nn<-1:tree$Nnode+Ntip(tree) x<-x[tree$tip.label] for(i in 1:tree$Nnode){ - tt<-read.tree(text=write.tree(rotate(tree,nn[i]))) + tt<-read.tree(text=write.tree(if(isSingleton(tree,nn[i])) tree else rotate(tree,nn[i]))) oo<-sum(abs(order(x[tree$tip.label])-1:length(tree$tip.label))) pp<-sum(abs(order(x[tt$tip.label])-1:length(tt$tip.label))) if(oo>pp) tree<-tt @@ -258,6 +263,8 @@ return(tree) } +isSingleton<-function(tree,node) if(sum(tree$edge[,1]==node)<=1) TRUE else FALSE + print.phylo.to.map<-function(x,...){ cat("Object of class \"phylo.to.map\" containing:\n\n") cat(paste("(1) A phylogenetic tree with",Ntip(x$tree),"tips and",x$tree$Nnode,"internal nodes.\n\n",sep=" ")) diff -Nru r-cran-phytools-0.7-47/R/phyl.pca.R r-cran-phytools-0.7-70/R/phyl.pca.R --- r-cran-phytools-0.7-47/R/phyl.pca.R 2019-08-21 19:42:38.000000000 +0000 +++ r-cran-phytools-0.7-70/R/phyl.pca.R 2020-08-26 14:37:22.000000000 +0000 @@ -1,7 +1,8 @@ ## function to perform phylogenetic principal components analysis ## multiple morphological traits in Y ## also can use lambda transformation in which lambda is optimized by ML or REML (in progress) -## written by Liam Revell 2010, 2011, 2013, 2015, 2016, 2017, 2019 ref. Revell (2009; Evolution) +## written by Liam Revell 2010, 2011, 2013, 2015, 2016, 2017, 2019, 2020 +## ref. Revell (2009; Evolution) phyl.pca<-function(tree,Y,method="BM",mode="cov",...){ ## get optional argument @@ -195,10 +196,8 @@ else newdata<-NULL if(!is.null(newdata)){ if(!is.matrix(newdata)) newdata<-as.matrix(newdata) - if(ncol(newdata)!=ncol(object$Evec)){ - if(nrow(newdata)==ncol(object$Evec)) newdata<-t(newdata) - else stop("Dimensions of newdata incorrect.") - } + if(ncol(newdata)!=nrow(object$Evec)) + stop("Dimensions of newdata incorrect.") n<-nrow(newdata) m<-ncol(newdata) A<-matrix(rep(object$a,n),n,m,byrow=TRUE) diff -Nru r-cran-phytools-0.7-47/R/plotBranchbyTrait.R r-cran-phytools-0.7-70/R/plotBranchbyTrait.R --- r-cran-phytools-0.7-47/R/plotBranchbyTrait.R 2020-02-29 18:05:10.000000000 +0000 +++ r-cran-phytools-0.7-70/R/plotBranchbyTrait.R 2020-07-20 18:15:28.000000000 +0000 @@ -166,3 +166,167 @@ subtitle,pos=1,cex=fsize,srt=90) } } + +## function to create "edge.widthMap" object +edge.widthMap<-function(tree,x,...){ + if(!inherits(tree,"phylo")) + stop("tree should be an object of class \"phylo\".") + tree<-as.phylo(tree) + a<-fastAnc(tree,x) + node.values<-c(x[tree$tip.label],a) + edge.values<-apply(tree$edge,1,function(e,nv) + mean(nv[e]),nv=node.values) + edge.widths<-edge.values + object<-list(tree=tree,edge.widths=edge.widths, + node.values=node.values) + class(object)<-"edge.widthMap" + object +} + +## print method +print.edge.widthMap<-function(x,...){ + cat("Object of class \"edge.widthMap\" containing:\n") + cat(paste("(1) Phylogenetic tree with",Ntip(x$tree), + "tips and",x$tree$Nnode,"internal nodes.\n")) + cat("(2) Vector of node values for a mapped quantitative\n") + cat(" trait.\n\n") +} + +## plot method +plot.edge.widthMap<-function(x,max.width=0.9,legend="trait value",...){ + if(hasArg(min.width)) min.width<-list(...)$min.width + else min.width<-0 + if(hasArg(vertical.as.polygon)) + vertical.as.polygon<-list(...)$vertical.as.polygon + else vertical.as.polygon<-TRUE + if(hasArg(lwd)) lwd<-list(...)$lwd + else lwd<-2 + h<-max(nodeHeights(x$tree)) + node.values<-x$node.values-min(x$node.values) + node.values<-node.values*((max.width-min.width)/ + max(node.values))+min.width + args.list<-list(...) + args.list$tree<-x$tree + args.list$type<-"phylogram" + if(!is.null(args.list$direction)){ + if(!args.list$direction%in%c("leftwards","rightwards")) + args.list$direction<-"rightwards" + } else args.list$direction<-"rightwards" + if(is.null(args.list$ylim)) + args.list$ylim<-c(1,Ntip(x$tree)+Ntip(x$tree)/25) + if(is.null(args.list$ftype)) args.list$ftype<-"i" + if(is.null(args.list$fsize)) + args.list$fsize<-36*par()$pin[2]/par()$pin[1]/ + Ntip(x$tree) + if(is.null(args.list$color)){ + args.list$color<-"transparent" + color<-"gray62" + } else { + color<-args.list$color + args.list$color<-"transparent" + } + if(is.null(args.list$border)){ + border<-color + } else { + border<-args.list$border + args.list$border<-NULL + } + do.call(plotTree,args.list) + obj<-get("last_plot.phylo",envir=.PlotPhyloEnv) + asp<-par()$pin[2]/par()$pin[1]/ + (diff(par()$usr[4:3])/diff(par()$usr[2:1])) + for(i in 1:nrow(x$tree$edge)){ + if(vertical.as.polygon){ + xx<-obj$xx[c(x$tree$edge[i,1], + x$tree$edge[i,1:2], + x$tree$edge[i,2:1], + x$tree$edge[i,1])]+c(0, + asp*node.values[x$tree$edge[i,1]]/2, + 0,0,asp*node.values[x$tree$edge[i,1]]/2,0) + yy<-rep(obj$yy[x$tree$edge[i,2]],6)+ + c(node.values[x$tree$edge[i,1]], + node.values[x$tree$edge[i,1:2]], + -node.values[x$tree$edge[i,2:1]], + -node.values[x$tree$edge[i,1]])/2 + } + else { + xx<-obj$xx[c(x$tree$edge[i,1:2], + x$tree$edge[i,2:1])] + yy<-rep(obj$yy[x$tree$edge[i,2]],4)+ + c(node.values[x$tree$edge[i,1:2]], + -node.values[x$tree$edge[i,2:1]])/2 + } + polygon(x=crop.to.h(xx,h),y=yy, + border=border,col=color,lwd=lwd) + } + for(i in 1:x$tree$Nnode+Ntip(x$tree)){ + nn<-x$tree$edge[which(x$tree$edge[,1]==i),2] + yy<-range(obj$yy[nn]) + if(vertical.as.polygon){ + xx<-rep(obj$xx[i],4)+ + asp*c(-node.values[i]/2,node.values[i]/2, + node.values[i]/2,-node.values[i]/2) + polygon(x=crop.to.h(xx,h), + y=c(rep(yy[1],2),rep(yy[2],2))+ + c(-rep(node.values[i],2), + rep(node.values[i],2))/2, + border=border,col=color,lwd=lwd) + } else { + lines(rep(obj$xx[i],2),yy+c(-node.values[i], + node.values[i])/2,lend=2,col=border, + lwd=lwd) + } + } + if(border!=color&&vertical.as.polygon){ + for(i in 1:nrow(x$tree$edge)){ + xx<-obj$xx[c(x$tree$edge[i,1], + x$tree$edge[i,1:2], + x$tree$edge[i,2:1], + x$tree$edge[i,1])]+c(0, + asp*node.values[x$tree$edge[i,1]]/2, + 0,0,asp*node.values[x$tree$edge[i,1]]/2,0) + yy<-rep(obj$yy[x$tree$edge[i,2]],6)+ + c(node.values[x$tree$edge[i,1]], + node.values[x$tree$edge[i,1:2]], + -node.values[x$tree$edge[i,2:1]], + -node.values[x$tree$edge[i,1]])/2 + polygon(x=crop.to.h(xx,h),y=yy, + border=FALSE,col=color) + } + for(i in 1:x$tree$Nnode+Ntip(x$tree)){ + nn<-x$tree$edge[which(x$tree$edge[,1]==i),2] + yy<-range(obj$yy[nn]) + xx<-rep(obj$xx[i],4)+ + asp*c(-node.values[i]/2,node.values[i]/2, + node.values[i]/2,-node.values[i]/2) + polygon(x=crop.to.h(xx,h), + y=c(rep(yy[1],2),rep(yy[2],2))+ + c(-rep(node.values[i],2), + rep(node.values[i],2))/2, + border=FALSE,col=color) + } + } + leg.length<-0.4*h + x.adj<-if(obj$direction=="rightwards") 0 else obj$x.lim[2]-leg.length + polygon(x=c(0,0,leg.length,leg.length)+x.adj, + y=Ntip(x$tree)+Ntip(x$tree)/25+ + c(-min.width/2,min.width/2,max(node.values)/2, + -max(node.values)/2), + border=border,col=color,lwd=lwd) + if(border!=color) + polygon(x=c(0,0,leg.length,leg.length)+x.adj, + y=Ntip(x$tree)+Ntip(x$tree)/25+ + c(-min.width/2,min.width/2,max(node.values)/2, + -max(node.values)/2), + border=FALSE,col=color) + text(0+x.adj,Ntip(x$tree)+Ntip(x$tree)/25-0.2*max.width, + round(min(x$node.values),2),pos=1, + cex=0.8) + text(leg.length+x.adj,Ntip(x$tree)+Ntip(x$tree)/25-0.2*max.width, + round(max(x$node.values),2),pos=1,cex=0.8) + text(leg.length/2+x.adj,Ntip(x$tree)+Ntip(x$tree)/25-0.2*max.width, + legend,pos=1,cex=0.8) + +} + +crop.to.h<-function(x,h) sapply(x,function(x,h) if(x<=h) x else h,h=h) \ No newline at end of file diff -Nru r-cran-phytools-0.7-47/R/plotTree.datamatrix.R r-cran-phytools-0.7-70/R/plotTree.datamatrix.R --- r-cran-phytools-0.7-47/R/plotTree.datamatrix.R 2018-02-14 00:36:59.000000000 +0000 +++ r-cran-phytools-0.7-70/R/plotTree.datamatrix.R 2020-08-24 19:25:53.000000000 +0000 @@ -1,9 +1,9 @@ ## function to plot a grid of discrete character data next to the tips of a tree -## written by Liam J. Revell 2018 +## written by Liam J. Revell 2018, 2020 plotTree.datamatrix<-function(tree,X,...){ N<-Ntip(tree) - ss<-sapply(X,function(x) levels(x)) + ss<-lapply(X,function(x) levels(x)) k<-sapply(ss,length) if(hasArg(fsize)) fsize<-list(...)$fsize else fsize<-40*par()$pin[2]/par()$pin[1]/Ntip(tree) @@ -21,7 +21,8 @@ while(length(palettes)1) make.era.map(tree,c(0,fit$par)) else make.era.map(tree,0) obj<-fn(mtree,x) - H<-optimHess(c(obj$sig2.multiple,fit$par),likHess,tree=tree,y=x,nrates=nrates,tol=tol,maxh=h) - vcv<-if(nrates>1) solve(H) else 1/H + if(compute.se){ + H<-optimHess(c(obj$sig2.multiple,fit$par),likHess,tree=tree,y=x,nrates=nrates,tol=tol,maxh=h) + vcv<-if(nrates>1) solve(H) else 1/H + } else vcv<-matrix(-1,2*nrates-1,2*nrates-1) if(nrates>1) rownames(vcv)<-colnames(vcv)<-c(paste("sig2(",1:nrates,")",sep=""),paste(1:(nrates-1),"<->",2:nrates,sep="")) else rownames(vcv)<-colnames(vcv)<-"sig2(1)" @@ -119,7 +123,7 @@ } ## S3 print method for object of class "rateshift" -## written by Liam J. Revell 2013 +## written by Liam J. Revell 2013, 2020 print.rateshift<-function(x,...){ sqroot<-function(x){ if(length(x)==1) if(x>=0) sqrt(x) else NaN @@ -129,28 +133,47 @@ else digits<-4 x<-lapply(x,function(a,b) if(is.numeric(a)) round(a,b) else a,b=digits) cat(paste("ML ",length(x$sig2),"-rate model:\n",sep="")) - cat(paste(c("",paste("s^2(",names(x$sig2),")","\tse(",names(x$sig2),")",sep=""), - "k","logL","\n"),collapse="\t")) - cat(paste(paste(c("value",paste(x$sig2,round(sqroot(diag(x$vcv)[1:length(x$sig2)]),digits), - sep="\t"),2*length(x$sig2),x$logL),collapse="\t"),"\n\n",sep="")) + tmp<-list() + nn<-vector() + for(i in 1:length(x$sig2)){ + tmp[2*i-1]<-x$sig2[i] + tmp[2*i]<-sqroot(diag(x$vcv)[i]) + nn[2*i-1]<-paste("s^2(",names(x$sig2)[i],")",sep="") + nn[2*i]<-paste("se(",names(x$sig2)[i],")",sep="") + } + tmp[2*i+1]<-2*length(x$sig2) + tmp[2*i+2]<-x$logL + nn[2*i+1]<-"k" + nn[2*i+2]<-"logL" + tmp<-as.data.frame(tmp) + colnames(tmp)<-nn + rownames(tmp)<-"value" + print(tmp,digits=digits) if(!is.null(x$shift)){ - cat("Shift point(s) between regimes (height above root):\n") - nn<-sapply(strsplit(names(x$shift),"<->"),paste,collapse="|") - cat(paste(c("",paste(nn,paste("se(",nn,")",sep=""),sep="\t"),"\n"), - collapse="\t")) - cat(paste(paste(c("value",paste(x$shift, - round(sqroot(diag(x$vcv)[1:length(x$shift)+length(x$sig2)]),digits), - sep="\t")),collapse="\t"),"\n\n",sep="")) - } else cat("This is a one-rate model.\n\n") - if(x$method=="ML") cat("Model fit using ML.\n\n") - else if(x$method=="REML") cat("Model fit using REML.\n\n") + cat("\nShift point(s) between regimes (height above root):\n") + tmp<-list() + nn<-vector() + for(i in 1:length(x$shift)){ + tmp[2*i-1]<-x$shift[i] + tmp[2*i]<-sqroot(diag(x$vcv)[i+length(x$sig2)]) + nn[2*i-1]<-paste(strsplit(names(x$shift[i]),"<->")[[1]],collapse="|") + nn[2*i]<-paste("se(",paste(strsplit(names(x$shift[i]),"<->")[[1]], + collapse="|"),")",sep="") + } + tmp<-as.data.frame(tmp) + colnames(tmp)<-nn + rownames(tmp)<-"value" + print(tmp,digits=digits) + } else cat("\nThis is a one-rate model.\n") + if(x$method=="ML") cat("\nModel fit using ML.\n\n") + else if(x$method=="REML") cat("\nModel fit using REML.\n\n") cat(paste("Frequency of best fit:",x$frequency.best,"\n\n")) if (x$convergence==0) cat(paste("R thinks it has found the",x$method,"solution.\n\n")) else cat("Optimization may not have converged.\n\n") } ## S3 logLik method for object of class "rateshift" -## written by Liam J. Revell 2013 +## written by Liam J. Revell 2013, 2020 logLik.rateshift<-function(object,...){ logLik<-object$logL class(logLik)<-"logLik" @@ -159,7 +182,7 @@ } ## S3 plot method for object of class "rateshift" -## written by Liam J. Revell 2015 +## written by Liam J. Revell 2015, 2020 plot.rateshift<-function(x,...){ if(length(x$sig2)>1){ cols<-colorRampPalette(c("blue","purple","red"))(101) @@ -179,10 +202,9 @@ colors<-setNames("blue",1) plot(x$tree,ylim=c(-0.1*Ntip(x$tree),Ntip(x$tree)), colors=colors,...) - txt<-as.character(round(x$sig2,3)) - add.simmap.legend(leg=expression(paste(sigma^2," = ",sep="")), - colors="blue",prompt=FALSE,x=0,y=-0.05*Ntip(x$tree)) - text(x=5.5*strwidth("W"),y=-0.05*Ntip(x$tree),round(x$sig2,3)) + legend(x=0,y=0, + legend=bquote(sigma^2 == .(round(x$sig2,3))), + pch=15,col="blue",pt.cex=2,bty="n") } } diff -Nru r-cran-phytools-0.7-47/R/rerootingMethod.R r-cran-phytools-0.7-70/R/rerootingMethod.R --- r-cran-phytools-0.7-47/R/rerootingMethod.R 2019-05-04 15:52:29.000000000 +0000 +++ r-cran-phytools-0.7-70/R/rerootingMethod.R 2020-07-06 17:40:13.000000000 +0000 @@ -1,5 +1,6 @@ -## function to compute the marginal posterior probabilities for nodes using the rerooting method -## written by Liam J. Revell 2013, 2015, 2017, 2018 +## function to compute the marginal posterior probabilities for nodes +## using the rerooting method +## written by Liam J. Revell 2013, 2015, 2017, 2018, 2020 rerootingMethod<-function(tree,x,model=c("ER","SYM"),...){ if(!inherits(tree,"phylo")) @@ -22,14 +23,17 @@ Q<-matrix(c(0,YY$rates)[YY$index.matrix+1],length(YY$states), length(YY$states),dimnames=list(YY$states,YY$states)) diag(Q)<--colSums(Q,na.rm=TRUE) - nn<-if(tips) c(1:n,2:tree$Nnode+n) else 2:tree$Nnode+n + nn<-if(tips) c(1:n,if(tree$Nnode>1) 2:tree$Nnode+n) else { + if(tree$Nnode>1) 2:tree$Nnode+n else vector() } ff<-function(nn){ - tt<-if(nn>Ntip(tree)) ape::root.phylo(tree,node=nn) else reroot(tree,nn,tree$edge.length[which(tree$edge[,2]==nn)]) + tt<-if(nn>Ntip(tree)) ape::root.phylo(tree,node=nn) else + reroot(tree,nn,tree$edge.length[which(tree$edge[,2]==nn)]) fitMk(tt,yy,model=model,fixedQ=Q,output.liks=TRUE)$lik.anc[1,] } XX<-t(sapply(nn,ff)) - if(tips) XX<-rbind(XX[1:n,],YY$lik.anc[1,],XX[(n+1):nrow(XX),]) - else XX<-rbind(yy,YY$lik.anc[1,],XX) + if(tips) XX<-rbind(XX[1:n,],YY$lik.anc[1,],if(tree$Nnode>1) + XX[(n+1):nrow(XX),]) + else XX<-rbind(yy,YY$lik.anc[1,],if(tree$Nnode>1) XX) rownames(XX)<-1:(tree$Nnode+n) if(tips) rownames(XX)[1:n]<-tree$tip.label XX<-if(tips) XX else XX[1:tree$Nnode+n,] @@ -39,21 +43,24 @@ } print.rerootingMethod<-function(x,digits=6,printlen=NULL,...){ - cat("Ancestral character estimates using re-rooting method\nof Yang et al. (1995):\n") + cat("Ancestral character estimates using re-rooting method") + cat("\nof Yang et al. (1995):\n") if(is.null(printlen)) print(round(x$marginal.anc,digits)) else { print(round(x$marginal.anc[1:printlen,],digits)) cat("...\n") } cat("\nEstimated transition matrix,\nQ =\n") print(round(x$Q,digits)) - cat("\n**Note that if Q is not symmetric the marginal\nreconstructions may be invalid.\n") + cat("\n**Note that if Q is not symmetric the marginal") + cat("\nreconstructions may be invalid.\n") cat(paste("\nLog-likelihood =",round(x$loglik,digits),"\n\n")) } plot.rerootingMethod<-function(x, ...){ args<-list(...) if(is.null(args$lwd)) args$lwd<-1 - if(is.null(args$ylim)) args$ylim<-c(-0.1*Ntip(x$tree),Ntip(x$tree)) + if(is.null(args$ylim)) args$ylim<-c(-0.1*Ntip(x$tree), + Ntip(x$tree)) if(is.null(args$offset)) args$offset<-0.5 if(is.null(args$ftype)) args$ftype="i" args$tree<-x$tree diff -Nru r-cran-phytools-0.7-47/R/simBMphylo.R r-cran-phytools-0.7-70/R/simBMphylo.R --- r-cran-phytools-0.7-47/R/simBMphylo.R 2019-12-14 17:55:21.000000000 +0000 +++ r-cran-phytools-0.7-70/R/simBMphylo.R 2020-08-20 15:41:58.000000000 +0000 @@ -38,6 +38,8 @@ else cex.axis<-0.9 if(hasArg(cex.lab)) cex.lab=list(...)$cex.lab else cex.lab<-1 + if(hasArg(las)) las<-list(...)$las + else las<-par()$las plotTree(tree,mar=c(0.1,4.1,2.1,1.1), xlim=c(0,1.05*max(H)), ylim=c(1-2*(Ntip(tree)-1)*0.04, @@ -48,8 +50,8 @@ plot.new() par(mar=c(5.1,4.1,1.1,1.1)) plot.window(xlim=c(0,1.05*max(H)),ylim=range(xx)) - axis(1,cex.axis=cex.axis) - axis(2,cex.axis=cex.axis) + axis(1,cex.axis=cex.axis,las=las) + axis(2,cex.axis=cex.axis,las=las) for(i in 1:length(xx)) lines(H[i,1]:H[i,2],xx[[i]]) for(i in 1:Ntip(tree)){ diff -Nru r-cran-phytools-0.7-47/R/threshBayes.R r-cran-phytools-0.7-70/R/threshBayes.R --- r-cran-phytools-0.7-47/R/threshBayes.R 2020-04-12 18:28:11.000000000 +0000 +++ r-cran-phytools-0.7-70/R/threshBayes.R 2020-08-25 01:46:59.000000000 +0000 @@ -100,7 +100,8 @@ cont<-which(types=="cont") for(i in 1:length(disc)) Y[Y[,disc[i]]==0,disc[i]]<--1 } - Y[,disc]<-Y[,disc]*abs(rnorm(n=length(as.vector(Y[,disc])))) + Y[,disc]<-Y[,disc]*abs(rnorm(n=length(as.vector(Y[,disc])), + sd=sqrt(max(nodeHeights(tree))))) npar<-npar-n*(m-length(disc))-length(disc) if(is.null(colnames(X))) colnames(Y)<-paste("V",1:m,sep="") @@ -164,7 +165,7 @@ } d<-i%%npar if(ngen>=con$print.interval) if(i%%con$print.interval==0) if(!con$quiet){ - cat(paste("genearation: ",i,"; mean acceptance rate: ",round(mean(accept.rate),2),"\n",sep="")) + cat(paste("generation: ",i,"; mean acceptance rate: ",round(mean(accept.rate),2),"\n",sep="")) flush.console() } Yp[]<-Y @@ -308,14 +309,20 @@ else xlim<-c(-1,1) if(hasArg(ylim)) ylim<-list(...)$ylim else ylim<-c(0,1.2*max(d$y)) + if(hasArg(bty)) bty<-list(...)$bty + else bty<-"n" + if(hasArg(cex.lab)) cex.lab<-list(...)$cex.lab + else cex.lab<-1 + if(hasArg(cex.axis)) cex.axis<-list(...)$cex.axis + else cex.axis<-1 plot(d,xlim=xlim,ylim=ylim,col="blue",xlab="Posterior sample of r", - ylab="Density",main="",bty="l") + ylab="Density",main="",bty=bty,cex.lab=cex.lab,cex.axis=cex.axis) polygon(x=c(min(d$x),d$x,max(d$x)),y=c(0,d$y,0), col=make.transparent("blue",0.2)) - lines(rep(r,2),c(0,par()$usr[4]),col="blue",lty="dashed", + lines(rep(r,2),c(0,max(d$y)),col="blue",lty="dashed", lwd=2) - text(r,0.95*par()$usr[4],"mean post-burnin\nvalue of r",cex=0.7, - pos=if(r>0) 2 else 4) + text(r,max(d$y),"mean post-burnin\nvalue of r",cex=0.7, + pos=if(r>0) 2 else 4,font=3) } density.threshBayes<-function(x,...){ @@ -343,24 +350,30 @@ plot.threshBayes<-function(x,...){ if(hasArg(bw)) bw<-list(...)$bw else bw<-floor(length(x$par$gen)/100) + if(hasArg(bty)) bty<-list(...)$bty + else bty<-"n" + if(hasArg(las)) las<-list(...)$las + else las<-1 + if(hasArg(cex.main)) cex.main<-list(...)$cex.main + else cex.main<-1 par(mfrow=c(3,1)) par(mar=c(5.1,4.1,2.1,1.1)) - plot(x$par$gen,x$par$logL,type="l",bty="l",col=make.transparent("grey",0.5), - xlab="generation",ylab="log(L)") - mtext("a) log-likelihood trace",side=3,line=0,cex=1,at=0,outer=FALSE, - adj=0) + plot(x$par$gen,x$par$logL,type="l",bty=bty,col=make.transparent("grey",0.5), + xlab="generation",ylab="log(L)",las=las) + mtext("a) log-likelihood trace",side=3,line=0.5,cex=cex.main,at=0, + outer=FALSE,adj=0) par(mar=c(5.1,4.1,2.1,1.1)) accept<-vector() for(i in 1:length(x$par$gen)) accept[i]<-mean(x$par$accept_rate[max(c(1,i-bw)):i]) - plot(x$par$gen,accept,type="l",bty="l",col=make.transparent("red",0.5), - xlab="generation",ylab="mean acceptance rate") + plot(x$par$gen,accept,type="l",bty=bty,col=make.transparent("red",0.5), + xlab="generation",ylab="mean acceptance rate",las=las) if(is.numeric(attr(x,"auto.tune"))) lines(c(par()$usr[1],max(x$par$gen)), rep(attr(x,"auto.tune"),2),lty="dotted") mtext(paste("b) mean acceptance rate (sliding window: bw=",bw,")",sep=""), - side=3,line=0,cex=1,at=0,outer=FALSE,adj=0) - plot(x$par$gen,x$par$r,type="l",bty="l",col=make.transparent("blue",0.5), - xlab="generation",ylab="r") - mtext("c) trace of the correlation coefficient, r",side=3,line=0,cex=1,at=0, - outer=FALSE,adj=0) + side=3,line=0.5,cex=cex.main,at=0,outer=FALSE,adj=0) + plot(x$par$gen,x$par$r,type="l",bty=bty,col=make.transparent("blue",0.5), + xlab="generation",ylab="r",las=las) + mtext("c) trace of the correlation coefficient, r",side=3,line=0.5, + cex=cex.main,at=0,outer=FALSE,adj=0) } diff -Nru r-cran-phytools-0.7-47/R/utilities.R r-cran-phytools-0.7-70/R/utilities.R --- r-cran-phytools-0.7-47/R/utilities.R 2020-05-29 16:29:28.000000000 +0000 +++ r-cran-phytools-0.7-70/R/utilities.R 2020-08-31 22:40:42.000000000 +0000 @@ -444,7 +444,7 @@ } ## function mostly to interactively label nodes by clicking -## written by Liam J. Revell 2017 +## written by Liam J. Revell 2017, 2020 labelnodes<-function(text,node=NULL,interactive=TRUE, shape=c("circle","ellipse","rect"),...){ shape<-shape[1] @@ -454,6 +454,8 @@ else rect.exp<-1.6 if(hasArg(cex)) cex<-list(...)$cex else cex<-1 + if(hasArg(bg)) bg<-list(...)$bg + else bg<-"white" obj<-get("last_plot.phylo",envir=.PlotPhyloEnv) h<-cex*strheight("A") w<-cex*strwidth(text) @@ -474,15 +476,15 @@ node[i]<-ii } else ii<-node[i] if(shape=="circle") - draw.circle(obj$xx[ii],obj$yy[ii],rad,col="white") + draw.circle(obj$xx[ii],obj$yy[ii],rad,col=bg) else if(shape=="ellipse") draw.ellipse(obj$xx[ii],obj$yy[ii],0.8*w[i],h, - col="white") + col=bg) else if(shape=="rect") rect(xleft=obj$xx[ii]-0.5*rect.exp*w[i], ybottom=obj$yy[ii]-0.5*rect.exp*h, xright=obj$xx[ii]+0.5*rect.exp*w[i], - ytop=obj$yy[ii]+0.5*rect.exp*h,col="white", + ytop=obj$yy[ii]+0.5*rect.exp*h,col=bg, ljoin=1) text(obj$xx[ii],obj$yy[ii],label=text[i],cex=cex) } @@ -2128,8 +2130,8 @@ return(type) } -# function 'untangles' (or attempts to untangle) a tree with crossing branches -# written by Liam J. Revell 2013, 2015 +## function 'untangles' (or attempts to untangle) a tree with crossing branches +## written by Liam J. Revell 2013, 2015, 2020 untangle<-function(tree,method=c("reorder","read.tree")){ if(inherits(tree,"multiPhylo")){ tree<-lapply(tree,untangle,method=method) @@ -2140,8 +2142,11 @@ method<-method[1] if(method=="reorder") tree<-reorder(reorder(tree,"pruningwise")) else if(method=="read.tree"){ + tip.label<-tree$tip.label + tree$tip.label<-1:Ntip(tree) if(inherits(tree,"simmap")) tree<-read.simmap(text=write.simmap(tree)) else tree<-if(Ntip(tree)>1) read.tree(text=write.tree(tree)) else read.newick(text=write.tree(tree)) + tree$tip.label<-tip.label[as.numeric(tree$tip.label)] } ii<-!names(obj)%in%names(attributes(tree)) attributes(tree)<-c(attributes(tree),obj[ii])