# richfitz/diversitree

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
1754 lines (1578 sloc) 63 KB
 \documentclass[11pt]{article} \usepackage{fullpage} \usepackage{natbib} \usepackage{amsmath} \usepackage{xcolor} \definecolor{navy}{rgb}{0,0,0.4} \usepackage[colorlinks,citecolor=navy,linkcolor=navy,urlcolor=navy]{hyperref} \usepackage{dtsweave} % Needed or Sweave adds it in: % \usepackage{Sweave} \SweaveOpts{keep.source=TRUE,eps=FALSE} \SweaveOpts{prefix.string=cache/diversitree-tutorial} <>= library(cacheSweave) ## Remove stars, as they don't look nice when typeset options(show.signif.stars=FALSE, continue=" ") ## This helps add real chi^2 numbers to the text chi <- function(m1, m2, digits=1) round(2*(m1$lnLik - m2$lnLik), digits) @ %\VignetteIndexEntry{Introduction to diversitree} %\VignettePackage{diversitree} \title{Analysing diversification with diversitree} \author{Rich FitzJohn, with GeoSSE by Emma Goldberg} \date{25 March 2012, version 0.9-2} \begin{document} \maketitle \tableofcontents \clearpage \section{Introduction} The diversitree package includes a number of comparative phylogenetic methods, mostly focusing on analysing diversification and character evolution. These methods all share a common set of tools for performing maximum likelihood (ML) parameter estimation, hypothesis testing and model comparison, and Bayesian parameter estimation through Markov chain Monte Carlo (MCMC). Included methods include: \begin{itemize} \item Diversification \begin{itemize} \item Constant rate birth-death models \citep{Nee-1994-305} \end{itemize} \item Character evolution \begin{itemize} \item Discrete trait evolution \citep{Pagel-1994-37} \item Univariate Brownian motion and Ornstein-Uhlenbeck models of continuous trait evolution \end{itemize} \item Joint character evolution and diversification: \begin{itemize} \item BiSSE: Binary trait speciation and extinction \citep{Maddison-2007-701}, and extensions for incompletely resolved phylogenies \citep{FitzJohn-2009-595}. \item MuSSE: Multiple State Speciation and Extinction. \item QuaSSE: Quantitative State Speciation and Extinction \citep{FitzJohn-2010-619}. \end{itemize} \item Joint character evolution and diversification, with speciation-triggered changes in character states: \begin{itemize} \item GeoSSE: Geographic State Speciation and Extinction \citep{Goldberg-2011-451}. \item BiSSEness: BiSSE-Node Enhanced State Shift \citep{Magnuson-Ford-2012}. \item ClaSSE: Cladogenetic State Speciation and Extinction \citep{Goldberg-classe}. \end{itemize} \end{itemize} In addition, variants of these models are available: \begin{itemize} \item Time dependent'': different time epochs have different parameters (implemented for BiSSE, MuSSE), and where rates are arbitrary functions of time (implemented for birth death, BiSSE, and MuSSE). \item MEDUSA-style partitioned analyses, where different regions of the tree have different parameters (implemented for birth-death, BiSSE, MuSSE, QuaSSE, and GeoSSE). \item Marginal ancestral state reconstruction for discrete characters (Pagel94''), BiSSE, and MuSSE. \item Stochastic character mapping for discrete traits \citep{Bollback-2006-88} \end{itemize} In the future, new methods will include \begin{itemize} \item Reflected Brownian motion for bounded traits \item Stochastic character mapping for discrete traits that affect speciation or extinction rates \end{itemize} For all methods, inference can be carried out under maximum likelihood, or in Bayesian analyses via MCMC (Markov Chain Monte Carlo). Phylogenies can also be simulated under most of the models. This tutorial is designed to give an overview of the features in diversitree. It does not aim to be a compete reference to the package, or claim to always follow best practice. The manual follows the structure above. Many of the examples are just taken from the online documentation; further examples can be found there. Most are fairly contrived -- if you have examples you would rather see here, I would welcome data sets. \subsection{Regenerating this file} This file is written in Sweave'', which allows mixing R code and \LaTeX\ markup. If you want to regenerate the file, or run the empirical examples, you will need some data files, available at \url{http://www.zoology.ubc.ca/prog/diversitree/files/}, which must be present in the directory \code{data}. Some of the code chunks take a very long time to run (total processing time is currently about 40 hours on a 2.8~GHz MacPro), and use the \code{cacheSweave} package to speed subsequent runs up. % This is no longer true -- see the makefile... % With Gregor Gorjanc's % \href{http://gregor.gorjanc.googlepages.com/Sweave.sh}{Sweave.sh}, the % file can be easily regenerated with \code{Sweave.sh -c -ld % diversitree-tutorial.Rnw}. All code requires that the diversitree package be loaded. <<>>= library(diversitree) @ \section{Constant-rate birth-death models} The ape package already has some support for constant-rate birth-death models, but diversitree duplicates this for completeness. The major differences are (1) the function is not constrained to positive diversification rates ($\mu$ can exceed $\lambda$), (2) support for both random taxon sampling and unresolved terminal clades (but see ape's \code{bd.ext}), and (3) run both MCMC and MLE fits to birth death trees. % The constant rate birth death model is a special case of the other diversification models implemented in the package, and is the simplest model in diversitree. Start with a simulated phylogeny, with speciation rate $\lambda=0.1$ and extinction rate $\mu=0.03$: <<>>= set.seed(2) phy <- tree.bd(c(.1, .03), max.taxa=100) @ % (plotted in figure \ref{fig:bd-tree}). The first step in any analysis in diversitree is to construct a likelihood function. For constant rate birth death models, this is done with \code{make.bd}: <<>>= lik <- make.bd(phy) @ To see the argument names of the likelihood function, use the \code{argnames} function <<>>= argnames(lik) @ This shows that \code{lik} takes a vector of two parameters $(\lambda, \mu)$. It will return the log likelihood of the parameters, following the calculations in \citep{Nee-1994-305}. <<>>= lik(c(.1, .03)) # -7.74086 @ % Most likelihood functions accept additional arguments; these are documented on their online help pages. The only additional argument accepted for birth-death model is to disable conditioning on survival (by default, the likelihood is conditional on two lineages surviving to the present, following \citep{Nee-1994-305}. <<>>= lik(c(.1, .03), condition.surv=FALSE) # -10.74823 @ \begin{figure} \centering <>= plot(phy, no.margin=TRUE, show.tip.label=FALSE) @ \caption{A birth-death tree, with speciation rate $\lambda=0.1$ and extinction rate $\mu=0.03$, simulated until it has 100 species (the actual simulation runs until the speciation event that would create the 101st species, but does not add this species).} \label{fig:bd-tree} \end{figure} To do a ML model fit, pass the likelihood function and a starting point guess to the \code{find.mle} function: <<>>= fit <- find.mle(lik, c(.1, .03), method="subplex") @ % The final argument here selects the method \code{subplex}'' for the ML search; other methods are available. (The default for birth-death models is \code{nlm}'', which may produce warnings about failure to converge for the example here.) To extract the coefficients from the fitted object, use the \code{coef} function: <<>>= coef(fit) @ and to extract the log-likelihood value at the ML point, use the \code{logLik} function <<>>= logLik(fit) @ which extracts the coefficients with some additional information, or extract the \code{lnLik} element from the list directly: <<>>= fit$lnLik @ Does this model fit much better than a model without extinction (a Yule, or pure birth, model)? You can constrain parameters of likelihood functions using the \code{constrain} function. To specify that the extinction rate,$\mu$should be zero: <<>>= lik.yule <- constrain(lik, mu ~ 0) @ % Argument names here must match those given by \code{argnames}. Run the ML search the same way as above, specifying a single starting parameter (I've used the speciation rate from the full model here): <<>>= fit.yule <- find.mle(lik.yule, coef(fit)[1], method="subplex") @ To perform a likelihood ratio test, use the \code{anova} function\footnote{This is an unfortunate convention in R: many packages use this function as a general model comparison function, and I've taken their lead here -- in a future version, I may add a \code{lrt} function, which should be clearer. No analysis of variance is performed.} The model with the nonzero extinction estimate is preferred, with$\chi^2_1=\Sexpr{chi(fit, fit.yule)}$<<>>= anova(fit, yule=fit.yule) @ Alternatively, we can use Markov chain Monte Carlo (MCMC) to perform a Bayesian analysis. Here, I will use a uniform prior on the interval$[0,\infty)$for both parameters by not specifying any prior. The \code{w} parameter is the tuning parameter. Here, it affects how many function evaluations will be needed per sample, but will not generally affect the rate of mixing (see the online help for \code{mcmc} for more information). <>= set.seed(1) samples <- mcmc(lik, fit$par, nsteps=10000, w=.1, print.every=0) @ The posterior distribution of these parameters, and the code to generate it, is in figure \ref{fig:bd-mcmc-distr}. \begin{figure} \centering <>= samples$r <- samples$lambda - samples$mu col <- c("#eaab00", "#004165", "#618e02") profiles.plot(samples[c("lambda", "mu", "r")], col.line=col, las=1, opacity=.75, legend.pos="topright") abline(v=c(.1, .03, .07), col=col, lty=2) @ \caption{Posterior probability distributions for the parameters of the constant rate birth death model. True values are indicated by the solid vertical lines. The bars at the bottom of the distributions and the shaded areas correspond to the 95\% credibility intervals. These include the two parameter$\lambda$and$\mu$, though the true diversification rate$r$lies above the 95\% credibility interval for that parameter.} \label{fig:bd-mcmc-distr} \end{figure} Analyses can also use trees where only a fraction of species are present in the phylogeny. To demonstrate this, let's drop 25 of the 100 species from the original tree at random: <<>>= set.seed(1) phy.sub <- drop.tip(phy, sample(100, 25)) @ When constructing the likelihood function, pass an argument \code{sampling.f} in, with a value on$(0,1]$representing the fraction of species that are descended from the root node that are included in the phylogeny (here,$75/100$). Then, run a ML analysis with \code{find.mle} as before: <<>>= lik.sub <- make.bd(phy.sub, sampling.f=75/100) fit.sub <- find.mle(lik.sub, c(.1, .03), method="subplex") coef(fit.sub) @ With fewer included species, test to see whether the full model is still preferred over the Yule model: <<>>= lik.sub.yule <- constrain(lik.sub, mu ~ 0) fit.sub.yule <- find.mle(lik.sub.yule, coef(fit.sub)[1], method="subplex") anova(fit.sub, yule=fit.sub.yule) @ In this small tree, the support for the model with extinction is only$\chi^2_1=\Sexpr{chi(fit.sub, fit.sub.yule)}$(compared with$\chi^2_1=\Sexpr{chi(fit, fit.yule)}$when we had all species), which is no longer significant at the 5\% level. \clearpage % Force figure to appear. \section{Markov models of discrete character evolution} As with the birth death models above, diversitree supports both simulating discrete characters and estimating rates of character evolution. The implemented methods are a little idiosyncratic, as they are primarily here for completeness (representing special cases of the joint diversification-character evolution models), but are also useful in their own right. First, we will simulate a binary trait on a birth death tree. The tree above with$100$tips is a bit unwieldy to visualise, so we'll make a smaller$50$tip tree: <<>>= set.seed(2) phy <- tree.bd(c(.1, .03), max.taxa=50) @ % Dammit - this is broken apparently... Then, on this tree, simulate a binary character where the rate of transition from state$0$to state$1$is$0.1$, and the reverse rate is$0.2$. We'll start the tree in state$0$(specified by \code{x0}). The argument \code{model="mk2"} specifies that we will use a Mk2 model: <<>>= set.seed(1) states <- sim.character(phy, c(.1, .2), x0=0, model="mk2") @ The simulation remembers the history at nodes too, which is displayed in figure \ref{fig:mk2-tree} (eventually full history will be recorded). \begin{figure} \centering <>= plot(phy, show.tip.label=FALSE, no.margin=TRUE) col <- c("#004165", "#eaab00") tiplabels(col=col[states+1], pch=19, adj=1) nodelabels(col=col[attr(states, "node.state")+1], pch=19) @ \caption{Character history for simulated trait and tree. Blue is state$0$, yellow is state$1$.} \label{fig:mk2-tree} \end{figure} Next, build a likelihood function with the \code{make.mk2} function, and run a ML analysis with \code{find.mle}, using an initial parameter guess of$(0.1, 0.1)$: <<>>= lik.mk2 <- make.mk2(phy, states) fit.mk2 <- find.mle(lik.mk2, c(.1, .1), method="subplex") coef(fit.mk2) @ In the fit, the$q_{10}$parameter is higher than$q_{01}$, the difference being a little larger than in the true model. See if this difference is statistically justified by running a model where the two$q$values are constrained to be equal: <<>>= lik.mk1 <- constrain(lik.mk2, q10 ~ q01) fit.mk1 <- find.mle(lik.mk1, .1, method="subplex") anova(fit.mk2, mk1=fit.mk1) @ This is significant ($\chi^2_1=\Sexpr{chi(fit.mk2, fit.mk1)}$) so we can conclude that asymmetric model fits better. \subsection{Drawing samples with MCMC} It is possible to run an MCMC analysis. However, care should be taken to choose priors carefully, as while$q_{10}/(q_{01} + q_{10})$is usually well characterised by the data, the overall rate$(q_{01} + q_{10})$is poorly defined. For small trees like this, essentially infinite values of character evolution are consistent with the data, with the tip states just drawn from the stationary distribution of the process. There are two supplied prior functions, but any function that takes a vector of parameters and returns the log prior probability may be used. First, consider an exponential prior with rate$10$, which gives a mean of$1/10$, and assume the same prior distribution for both parameters. <<>>= prior.exp <- make.prior.exponential(10) @ To run the MCMC, we need to specify a starting point (again, I have used$(.1, .1)$, but the ML point might be preferable). I have discarded the first 500 samples (10\%), which is probably overkill for this model, as the autocorrelation between samples is extremely small. <>= set.seed(1) samples <- mcmc(lik.mk2, c(.1, .1), nsteps=5000, prior=prior.exp, w=.1, print.every=0) samples <- subset(samples, i > 500) @ \begin{figure} \centering <>= col <- c("#004165", "#eaab00") profiles.plot(samples[c("q01", "q10")], col.line=col, las=1, legend.pos="topright") abline(v=c(.1, .2), col=col) @ \caption{Posterior probability distributions for the parameters of the Mk2 model. True values are indicated by the solid vertical lines.} \label{fig:mk2-mcmc} \end{figure} The marginal distributions of these parameters are shown in figure \ref{fig:mk2-mcmc}, and overlap substantially. However, that is because the distributions are correlated (increasing$q_{01}$fits best if$q_{10}$is also increased). The$q_{01}$parameter is estimated to be greater than the$q_{10}$parameter only a small fraction of the time: <<>>= mean(samples$q01 > samples$q10) @ % TODO: Meristic and multitrait. \clearpage \section{Binary traits and diversification: BiSSE} \label{sec:bisse} The BiSSE (Binary State Speciation and Extinction) model combines the features of the constant-rates birth-death model with the two-state Markov model. Again, start with a simulated tree. The parameters here are in the order$\lambda_0$,$\lambda_1$,$\mu_0$,$\mu_1$,$q_{01}$,$q_{10}$, so the parameters below correspond to an asymmetry in the speciation rate where state$1$speciates at twice the rate as state$0$. All other parameters are equal between states. <<>>= pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01) set.seed(4) phy <- tree.bisse(pars, max.t=30, x0=0) states <- phy$tip.state head(states) @ This gives a 52 species tree, shown with its true history in figure \ref{fig:bisse-tree}. The character states are now stored in the \code{states} vector. This vector is named, so that each element can be easily associated with a tip in the tree. \begin{figure} \centering <>= par(mar=rep(0, 4)) col <- c("#004165", "#eaab00") plot(history.from.sim.discrete(phy, 0:1), phy, col=col) @ \caption{A BiSSE tree, with parameters $\lambda_0=0.1$, $\lambda_1=0.2$, $\mu_0=\mu_1=0.03$, and $q_{01}=q_{10}=0.01$. Blue is state $0$, yellow is state $1$.} \label{fig:bisse-tree} \end{figure} The \code{make.bisse} takes as its first two arguments a tree and set of character states (these are the only mandatory arguments): <<>>= lik <- make.bisse(phy, states) lik(pars) # -159.71 @ To perform an ML search, we need a starting point. The \code{starting.point.bisse} function produces a basic heuristic guess of a sensible starting point, based on the character-independent birth-death fit. There are no guarantees that this is at all close to the ML point, or that the ML point can be reached from this point while climbing uphill only (which most optimisers assume). <<>>= p <- starting.point.bisse(phy) p @ Start an ML search from this point (this may take some time) <>= fit <- find.mle(lik, p) @ As above, the \code{fit.mle} object has an element \code{lnLik} with the log-likelihood value <<>>= fit$lnLik @ and coefficients may be extracted with \code{coef} (rounded for clarity): <<>>= round(coef(fit), 3) @ Let's test the hypothesis that the speciation rates are different for the different states. We can use \code{constrain} to enforce equal speciation rates to be equal across character states: <<>>= lik.l <- constrain(lik, lambda1 ~ lambda0) @ and then start the ML search again: <>= fit.l <- find.mle(lik.l, p[argnames(lik.l)]) fit.l$lnLik # -158.74 @ % (the statement \code{p[argnames(lik.l)]}'' drops the $\lambda_1$ element from the starting parameter vector). This fit has quite different parameters to the full model (compare $\mu_0$) <<>>= round(rbind(full=coef(fit), equal.l=coef(fit.l, TRUE)), 3) @ (the \code{TRUE} argument forces \code{coef} to return values for constrained parameters). However, the difference in fits is not statistically supported, with $\chi^2_1=\Sexpr{chi(fit, fit.l)}$ <<>>= anova(fit, equal.l=fit.l) @ \subsection{Analysis with MCMC} Because we are fitting six parameters to a tree with only 52 species, priors will be needed so that the posterior distribution is proper. I will use an exponential prior with rate $1/(2r)$, where $r$ is the character independent diversification rate: <<>>= prior <- make.prior.exponential(1 / (2 * (p[1] - p[3]))) @ The MCMC sampler in diversitree uses slice sampling \citep{Neal-2003-705} for parameter updates. The step size'' (argument \code{w}) does not need to be carefully tuned as it does not affect the rate of mixing -- just the number of function evaluations per update. Ideally it will be on the same order as the width of the high probability region''. An easy way of setting this is to run a short chain (say, 100 steps) and use the range of observed samples as a measure of this. <>= set.seed(1) tmp <- mcmc(lik, fit$par, nsteps=100, prior=prior, lower=0, w=rep(1, 6), print.every=0) w <- diff(sapply(tmp[2:7], range)) w @ Run the chain for 10,000 steps (this will take a while) <>= samples <- mcmc(lik, fit$par, nsteps=10000, w=w, lower=0, prior=prior, print.every=0) @ The marginal distributions for the two speciation rates are shown in figure \ref{fig:bisse-mcmc}, which shows the 95\% credibility intervals for $\lambda_0$ completely overlapping those for $\lambda_1$. \begin{figure} \centering <>= col <- c("#004165", "#eaab00") profiles.plot(samples[c("lambda0", "lambda1")], col.line=col, las=1, xlab="Speciation rate", legend="topright") abline(v=c(.1, .2), col=col) @ \caption{Posterior probability distributions for $\lambda_0$ and $\lambda_1$ for a BiSSE model. True values are indicated by the solid vertical lines. The bars at the bottom of the distributions and the shaded areas correspond to the 95\% credibility intervals.} \label{fig:bisse-mcmc} \end{figure} \subsection{Incomplete taxonomic sampling} Not all phylogenies are complete, but the basic BiSSE calculations assume that they are. If given tree contains a random sample of all extant species, the calculations can be corrected. To demonstrate this, we will generate a larger tree (150 taxa), and drop 50 taxa from it. Here, I am using the same parameters as earlier. <<>>= pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01) set.seed(4) phy <- tree.bisse(pars, max.taxa=150, x0=0) states <- phy$tip.state phy.s <- drop.tip(phy, setdiff(seq_len(150), sample(150, 50))) states.s <- states[phy.s$tip.label] @ Calculate what the sampling fraction is for this tree. You can either assume that the sampling fraction is independent of the character state: <<>>= sampling.f <- 50 / 150 @ or you can assume that it varies with character state <<>>= sampling.f <- as.numeric(table(states.s) / table(states)) sampling.f @ Pass this in to \code{make.bisse} and construct a new likelihood function that accounts for the sampling: <<>>= lik.s <- make.bisse(phy.s, phy.s$tip.state, sampling.f=sampling.f) @ This can then be optimised, as before: <>= p <- starting.point.bisse(phy) fit.s <- find.mle(lik.s, p) fit.s[1:2] @ \subsection{Terminally unresolved trees} Another way that phylogenies might be incompletely resolved is that higher level relationships may be known (say, genera), but little or nothing is known about species relationships within these groups. This results in trees where some taxa'' represent a number of species -- terminal clades''. There are a couple of different ways that unresolved clade information may be specified. To demonstrate this, I will use an example of sexual dimorphism in shorebirds; this is the same example as in \citep{FitzJohn-2009-595}. The phylogeny is a supertree constructed by \citet{Thomas-2004-28}, and the data on sexual size dimorphism are derived from \citet{Lislevand-2007-1605}. The required files can be downloaded from \url{http://www.zoology.ubc.ca/prog/diversitree/files/} Read in the phylogenetic tree <<>>= tree <- read.nexus("data/Thomas-tree.nex") tree <- ladderize(tree) @ The tree contains many polytomies; the original tree, and a simplified tree with polytomies converted into clades are shown in figure \ref{fig:shorebird-tree}. \setkeys{Gin}{width=1.0\textwidth} \begin{figure}[p] \centering <>= par(mar=rep(.5, 4), mfrow=c(1, 2)) plot(tree, show.tip.label=FALSE) tree.clade <- clades.from.polytomies(tree) plt <- plot(tree.clade, show.tip.label=FALSE, clade.fill="gray") @ \caption{The shorebird supertree \citep{Thomas-2004-28}. On the left, the original tree with all polytomies and 350 species. On the right, the polytomies have been collapsed, to leave 135 tips, with tips representing from 1 to 47 species.} \label{fig:shorebird-tree} \end{figure} \setkeys{Gin}{width=0.8\textwidth} The character states are stored as the size of the difference of mass between sexes, divided by the mean across sexes \citep[see][]{FitzJohn-2009-595}. <<>>= d <- read.csv("data/Lislevand-states.csv", as.is=TRUE) states <- d$dimorph names(states) <- d$species states <- states[tree$tip.label] names(states) <- tree$tip.label head(states) @ These will need converting to a binary character for use, for example -- to convert this into a binary character where an absolute relative difference of 10\% would be considered dimorphic'': <<>>= head((abs(states) > .1) + 0) @ The simplest way of working with this tree is to use the \code{clades.from.polytomies} function. This collapses all daughters of any polytomy into a clade. (Be careful - if you have a polytomy at the base of your tree, the entire tree will collapse into a single clade!) This can be visualised with \code{plot} functions as normal -- see \code{?plot.clade.tree} for more information. This tree can then be passed into \code{make.bisse}, along with a plain vector of state names. The catch is that every taxon still needs state information, not just those at the tips. Our state vector \code{states} includes states for all 350 species, so we are OK to use this. <<>>= states.15 <- (abs(states) > 0.15) + 0 tree.clade <- clades.from.polytomies(tree) lik <- make.bisse(tree.clade, states.15) @ \setkeys{Gin}{width=1.0\textwidth} \begin{figure}[p] \centering <>= ## Diversitree imports: group.label.tip <- diversitree:::group.label.tip filled.arcs <- diversitree:::filled.arcs add.alpha <- diversitree:::add.alpha ## Determine family for each tip: families <- read.csv("data/shorebird-families.csv", as.is=TRUE) family <- families$family[match(tree.clade$tip.label, families$spp)] ## Summarise trait over clades, and make vector where 1 is grey, 0 is ## black, and values in between are darker or lighter. nn <- sapply(tree.clade$clades, function(x) tabulate(states.15[x]+1, 2)) p <- states.15[tree.clade$tip.label] p[colnames(nn)] <- nn[2,] / colSums(nn) cols <- add.alpha("black", (1 - p)*2/3 + 1/3) names(cols) <- tree.clade$tip.label w <- max(branching.times(tree)) / 30 plt <- plot(tree.clade, show.tip.label=FALSE, no.margin=TRUE, transform=sqrt, clade.fill="gray", label.offset=33, type="fan") group.label.tip(plt, family, "black", "black", 1, 2*w, 3*w, cex=.8) theta <- plt$xy$theta[seq_along(tree.clade$tip.label)] dt <- (plt$n.taxa/2 - .1) / plt$n.spp * 2 * pi filled.arcs(theta - dt, theta + dt, max(plt$xy$x) + w/2, w, cols) @ \caption{Phylogenetic tree of the 350 species of shorebirds (Charadriiformes) and measures of sexual dimorphism, based on \citet{Thomas-2004-28}. Gray triangles indicate unresolved clades, with the height of the triangle being proportional to the square root of the number of species. Character states at the 15\% threshold level are indicated at the tips; gray indicates sexually dimorphic, black indicates sexually monomorphic, and white indicates no data. For unresolved clades, the degree of shading indicates the proportion of species in each state. For clarity, only family or subfamily names are shown.} \label{fig:shorebird-tree-with-traits} \end{figure} \setkeys{Gin}{width=0.8\textwidth} A sensible starting point can still be computed with \code{starting.point.bisse}; this takes into account the clade information automatically. <<>>= p <- starting.point.bisse(tree.clade) @ You can then perform a ML search or MCMC analysis as usual. However, unresolved clades slow down the analysis considerably, and this will take several minutes. <>= fit.full <- find.mle(lik, p) @ In the full model, speciation rates for dimorphic species (state $1$) are greater than those for monomorphic species, but character transition rates away from dimorphism ($q_{10}$) are greater than the reverse transition ($q_{10}$): <<>>= round(coef(fit.full), 3) @ To test whether these differences are significant, we can use a likelihood ratio test. First, construct the reduced models, constraining $\lambda_1 \sim \lambda_0$ or $q_{10} \sim q_{01}$: <<>>= lik.l <- constrain(lik, lambda1 ~ lambda0) lik.q <- constrain(lik, q10 ~ q01) @ Then, rerun the <>= fit.l <- find.mle(lik.l, p[argnames(lik.l)]) fit.q <- find.mle(lik.q, p[argnames(lik.q)]) @ These constrained models are significantly worse fits than the full model ($p \approx 0.04$ for both). <<>>= anova(fit.full, equal.l=fit.l, equal.q=fit.q) @ The analysis can also be run with MCMC. Here, I am using an exponential prior with rate $1 / (2r)$, as earlier. <<>>= prior <- make.prior.exponential(1 / (2*p[1])) @ The MCMC itself takes a very long time to run (approximately 4--5~s/sample) <>= set.seed(1) samples <- mcmc(lik, coef(fit.full), 10000, w=.3, prior=prior, print.every=0) @ <>= .save <- "cache/shorebirds-mcmc-samples.Rdata" if ( file.exists(.save) ) { load(.save) } else { <> save(samples, file=.save) } @ \begin{figure}[p] \centering <>= par(mar=c(4.1, 4.1, .5, .5)) col <- c("#004165", "#618e02") ## col <- c("#004165", "#eaab00") profiles.plot(samples[c("lambda0", "lambda1")], col, las=1, legend.pos="topright", xlab="Speciation rate estimate") @ \caption{Posterior probability distributions for speciation rates for non-dimorphoc species (blue) and dimorphic species (yellow), estimated using MCMC.} \label{fig:shorebird-mcmc-plot} \end{figure} Setting up the likelihood function above assumed that the tree being used contained polytomies, and this is the source of the unresolved clades. However, it is probably more common to have an exemplar'' tree, where the unresolved species were never included in the first place. The tree \code{Thomas-tree-exemplar.nex} (which was derived from the tree above) does not contain any reference the species that are contained within the unresolved clades. <<>>= tree.ex <- read.nexus("data/Thomas-tree-exemplar.nex") states.ex <- states.15[tree.ex$tip.label] names(states.ex) <- tree.ex$tip.label @ We need to define a \code{data.frame} with information about the unresolved clades. The file \code{Thomas-unresolved.csv} contains the information in the correct format: <<>>= unresolved <- read.csv("data/Thomas-unresolved.csv", as.is=TRUE) head(unresolved) @ All the columns here are required: \begin{itemize} \item \code{tip.label}: the tip label within the tree \item \code{Nc}: the total number of species that the tip represents \item \code{n0}: the number of species known to be in state $0$ \item \code{n1}: the number of species known to be in state $1$ \end{itemize} (additional columns are fine and will be silently ignored). Note that \code{Nc} can be greater than \code{n0 + n1}: this allows for species with unknown state. This \code{unresolved} object is passed into the \code{make.bisse} function: <<>>= lik.ex <- make.bisse(tree.ex, states.ex, unresolved=unresolved) @ This likelihood function should be identical to the one created above. <<>>= lik.ex(coef(fit.full)) lik(coef(fit.full)) @ \section{Multiple state characters and diversification: MuSSE} MuSSE (Multiple State Speciation and Extinction) generalises the BiSSE model to allow characters with more than two states. Following \citet{Pagel-1994-37}, this can also allow for multiple characters, each of which might be binary, by recoding the states. To illustrate, we'll start with a simple simulated example with a three-level state. The tree is simulated where character evolution is only possible among neighbouring states (i.e., $1\to 3$ and $3\to 1$ transitions are disallowed). All other transitions are equal, and both speciation and extinction rates increase as the character number increases. For a three state case, the parameter vector is in the order $(\lambda_1, \lambda_2, \lambda_3, \mu_1, \mu_2, \mu_3, q_{12}, q_{13}, q_{21}, q_{23}, q_{31}, q_{32})$. This order can be seen here (sorry -- clunky at the moment) <<>>= diversitree:::default.argnames.musse(3) @ (the order of the $q$ parameters is row-wise through the transition rate matrix, skipping diagonal elements). Simulate a 30 species tree, with the tree starting in state $1$. <<>>= pars <- c(.1, .15, .2, # lambda 1, 2, 3 .03, .045, .06, # mu 1, 2, 3 .05, 0, # q12, q13 .05, .05, # q21, q23 0, .05) # q31, q32 set.seed(2) phy <- tree.musse(pars, 30, x0=1) @ The tree and its real character history are shown in figure \begin{figure} \centering <>= col <- c("#eaab00", "#004165", "#618e02") h <- history.from.sim.discrete(phy, 1:3) plot(h, phy, cex=1, col=col, no.margin=TRUE, font=1) @ \caption{Simulated MuSSE tree. Yellow is state 1, blue is state 2, and green is state 3} \label{fig:musse-tree} \end{figure} The states are numbered $1,2,3$, rather than $0,1$ in BiSSE. <<>>= states <- phy$tip.state table(states) @ Making a likelihood function is basically identical to BiSSE. The third argument needs to be the number of states. In a future version this will probably be \code{max(states)}, but there are some pitfalls about this that I am still worried about. <<>>= lik <- make.musse(phy, states, 3) @ The argument names here are in the same order as for the simulation. Just adding one more state (compared with BiSSE) has moved us up to 10 parameters. <<>>= argnames(lik) @ Rather than start with the full model, and constrain things, here I will start with a very simple model and expand. This model has all$\lambda_i$,$\mu_i$, and$q_i$the same (except for$q_{13}$and$q_{31}$, which are still zero). <<>>= lik.base <- constrain(lik, lambda2 ~ lambda1, lambda3 ~ lambda1, mu2 ~ mu1, mu3 ~ mu1, q13 ~ 0, q21 ~ q12, q23 ~ q12, q31 ~ 0, q32 ~ q12) argnames(lik.base) @ Find the ML point for this model <>= p <- starting.point.musse(phy, 3) fit.base <- find.mle(lik.base, p[argnames(lik.base)]) @ Now, allow the speciation rates to vary <>= lik.lambda <- constrain(lik, mu2 ~ mu1, mu3 ~ mu1, q13 ~ 0, q21 ~ q12, q23 ~ q12, q31 ~ 0, q32 ~ q12) fit.lambda <- find.mle(lik.lambda, p[argnames(lik.lambda)]) @ There is very little improvement here (this is a small tree) <<>>= anova(fit.base, free.lambda=fit.lambda) @ % TODO: Include a multitrait case here -- possibly just the case from % the diversitree paper. \section{Quantitative traits and diversification: QuaSSE} The QuaSSE method is by far the slowest method in diversitree. This is because to compute the likelihood, we have to solve a system of partial differential equations, which is substantially more complicated than the system of ordinary differential equations or simple algebraic expressions in the other methods. The basic interface is very similar to the other methods, but there are more options to control the behaviour of the integrator. We'll start with a simulated tree. The tree simulation differs slightly from the other methods, because there is no longer a canonical argument list (speciation and extinction rates are arbitrary functions of the character state). % Here is a set of functions; speciation rate is a sigmoidal function which ranges from$0.1$to$0.2$with an inflection point at$x=0$, extinction is constant at rate$0.03$, and the model of character evolution is Brownian motion with diffusion parameter$0.025$. <<>>= lambda <- function(x) sigmoid.x(x, 0.1, 0.2, 0, 2.5) mu <- function(x) constant.x(x, 0.03) char <- make.brownian.with.drift(0, 0.025) @ Simulate the tree: <>= set.seed(1) phy <- tree.quasse(c(lambda, mu, char), max.taxa=15, x0=0, single.lineage=FALSE) @ We need to specify the standard deviation for the states; here I will just assume that all taxa have a state standard deviation of$1/200$<<>>= states <- phy$tip.state states.sd <- 1/200 @ Then, build the likelihood as usual. The difference compared with other models is that we have to specify the speciation and extinction functions (here, \code{sigmoid.x} and \code{constant.x}, respectively). There are a number of other provided functions (see \code{?constant.x} for a list), but any function that takes \code{x} as the first argument may be used. <<>>= lik <- make.quasse(phy, states, states.sd, sigmoid.x, constant.x) @ This can be used in ML calculations as usual. There is a \code{starting.point.quasse} function that may be useful in selecting sensible starting points, but some effort is still required to convert this into a full vector as it just returns constant rate speciation, extinction, and diffusion rates. <<>>= p <- starting.point.quasse(phy, states) p @ Let's ignore drift: the argument list we need is: <<>>= lik.nodrift <- constrain(lik, drift ~ 0) argnames(lik.nodrift) @ A sensible starting point here might be <<>>= p.start <- c(p[1], p[1], mean(states), 1, p[2:3]) names(p.start) <- argnames(lik.nodrift) p.start @ Lower bounds: <<>>= lower <- c(0, 0, min(states), -Inf, 0, 0) @ Then run \code{find.mle}, as usual. The \code{control} argument here just tells the subplex algorithm to use an initial step size of 0.1 (rather than 1), which reduces the number of function evaluations somewhat. <>= fit <- find.mle(lik.nodrift, p.start, control=list(parscale=.1), lower=lower, verbose=0) @ Compare this against the constant rate speciation fit: <>= lik.constant <- constrain(lik.nodrift, l.y1 ~ l.y0, l.xmid ~ 0, l.r ~ 1) fit.constant <- find.mle(lik.constant, p.start[argnames(lik.constant)], control=list(parscale=.1), lower=0, verbose=0) @ and compare the models -- no significant difference, which is not surprising with only a 15 species tree. <<>>= anova(fit, constant=fit.constant) @ \subsection{Primate analysis} Here, I will recreate the analysis of primate diversification from \citep{FitzJohn-2010-619}, but fitting only speciation functions to keep things relatively simple. <<>>= phy <- read.nexus("data/Vos-2006.nex") d <- read.csv("data/Redding-2010.csv") mass <- log(d$mass) names(mass) <- d$tip.label @ Assume standard deviation of 1/50 for all species -- this comes from nowhere in particular, and is probably over-confident in the mass estimates for most species. <<>>= mass.sd <- 1/50 @ Starting point parameter estimates, as above: <<>>= p <- starting.point.quasse(phy, mass) p @ Create a piecewise linear'' function. This is linear in the range [\code(xr[1]), \code{xr[2]}], and flat outside this range; this satisfies the condition that the derivatives of the speciation and extinction function with respect to the character state approach zero at the edges of the modelled parameter space. <<>>= xr <- range(mass) + c(-1,1) * 20 * p["diffusion"] linear.x <- make.linear.x(xr[1], xr[2]) @ Because we are going to create a reasonable number of models, here is a function that simplifies this, requiring just speciation and extinction functions: <<>>= make.primates <- function(lambda, mu) make.quasse(phy, mass, mass.sd, lambda, mu) @ and a function that constrains drift to zero <<>>= nodrift <- function(f) constrain(f, drift ~ 0) @ Create the likelihood functions where speciation is a constant, linear, sigmoidal, or hump-shaped function of log body mass. <<>>= f.c <- make.primates(constant.x, constant.x) f.l <- make.primates(linear.x, constant.x) f.s <- make.primates(sigmoid.x, constant.x) f.h <- make.primates(noroptimal.x, constant.x) @ Start by fitting the constant model first (this will take quite a while; each function evaluation takes about $5$~s). <>= control <- list(parscale=.1, reltol=0.001) mle.c <- find.mle(nodrift(f.c), p, lower=0, control=control, verbose=0) @ Starting points for the constrained analyses based on this constrained fit. <<>>= p.c <- mle.c$par p.l <- c(p.c[1], l.m=0, p.c[2:3]) p.s <- p.h <- c(p.c[1], p.c[1], mean(xr), 1, p.c[2:3]) names(p.s) <- argnames(nodrift(f.s)) names(p.h) <- argnames(nodrift(f.h)) @ <>= mle.l <- find.mle(nodrift(f.l), p.l, control=control, verbose=0) mle.s <- find.mle(nodrift(f.s), p.s, control=control, verbose=0) mle.h <- find.mle(nodrift(f.h), p.h, control=control, verbose=0) @ \begin{figure} \centering <>= ## Work out where the "peak region" is: p.h <- coef(mle.h) y <- diff(p.h[1:2]) * .1 + p.h[1] f <- function(x) noroptimal.x(x, p.h[1], p.h[2], p.h[3], p.h[4]) - y peak <- c(uniroot(f, c(min(mass), p.h[3]))$root, uniroot(f, c(max(mass), p.h[3]))$root) ## Ranges of body sizes to plot (log scale) xr <- c(2, max(mass)) xd <- diff(xr) par(mfrow=c(1, 1), mar=c(2.5, .5, .5, .5), mgp=c(3, .5, 0)) plot(phy, cex=.9, label.offset=1.5, tip.color="white", x.lim=c(0, 107)) lastPP <- get("last_plot.phylo", envir=.PlotPhyloEnv) xx <- lastPP$xx[1:lastPP$Ntip] yy <- lastPP$yy[1:lastPP$Ntip] x0 <- max(xx) + 0.02 * diff(par("usr")[1:2]) x1 <- max(xx) + 0.15 * diff(par("usr")[1:2]) mass2 <- (mass - xr[1]) / xd mass3 <- mass2 * (x1 - x0) + x0 segments(x0, yy, mass3, yy, col=ifelse(mass < peak[1] | mass > peak[2], "black", "red")) abline(v=(peak - xr[1])/xd * (x1 - x0) + x0, lty=2) axis.at <- pretty(c(mass, xr[1], xr[2])) axis.at <- axis.at[axis.at >= xr[1] & axis.at <= xr[2]] axis(1, (axis.at - xr[1])/xd * (x1 - x0) + x0, axis.at, cex.axis=.7) axisPhylo(cex.axis=.7) dx <- diff(par("usr")[1:2]) dy <- diff(par("usr")[3:4]) text(mean(range(lastPP$xx)), min(yy) - 0.11 * dy, "Time (Ma)", xpd=NA, cex=.7) text(mean(c(x1, x0)), min(yy) - 0.11 * dy, "ln(mass)", xpd=NA, cex=.7) @ \caption{The primate tree from \citet{Vos-2006}. Log body size (in grams) is shown by the horizontal bar for each species \citep[data from][]{Redding-2010-1052}. The vertical dashed lines indicate the approximate ranges of body masses in which elevated speciation rates were inferred, and extant species whose mass falls in this range have their mass coloured red.} \label{fig:primate-tree} \end{figure} The fits can then be compared. These are all against the constant speciation rate fit (listed as full'' in the table). The support is strongest for the hump'' shaped fit. <<>>= anova(mle.c, linear=mle.l, sigmoidal=mle.s, hump=mle.h) @ Run the fits with the drift parameter added, starting from the constrained model's ML parameters: <>= mle.d.l <- find.mle(f.l, coef(mle.l, TRUE), control=control, verbose=0) mle.d.s <- find.mle(f.s, coef(mle.s, TRUE), control=control, verbose=0) mle.d.h <- find.mle(f.h, coef(mle.h, TRUE), control=control, verbose=0) @ and add these to the comparison: <<>>= anova(mle.c, linear=mle.l, sigmoidal=mle.s, hump=mle.h, drift.linear=mle.d.l, drift.sigmoidal=mle.d.s, drift.hump=mle.d.h) @ There is an improvement in the model fit when drift is added. In all cases, the drift parameter is positive, indicating an increase in mass. <<>>= c(linear=coef(mle.d.l)[["drift"]], sigmoidal=coef(mle.d.s)[["drift"]], hump=coef(mle.d.h)[["drift"]]) @ However, the precise estimate depends strongly on the speciation model assumed. Next, consider splitting the tree. MEDUSA \citep{Alfaro-2009-13410} identified a shift in diversification rates at the base of the superfamily Cercopithecoidea (old world monkeys). In this tree, this corresponds to node $153$. It's easiest to work with node names, so I'm going to add some here <<>>= phy$node.label <- paste("nd", 1:phy$Nnode, sep="") @ Then, we can construct split'' QuaSSE objects. For the constant speciation rate case: <<>>= f.cc <- make.quasse.split(phy, mass, mass.sd, constant.x, constant.x, "nd153", Inf) argnames(f.cc) @ Here, the speciation and extinction functions may either be a single function (as above), in which case all partitions get the same function for speciation and extinction. Alternatively, lists of functions can be added, in which case different partitions may have different functions. The \code{Inf} indicates that the split should go at the base of the internal edge subtending node \code{nd153}, which is consistent with MEDUSA. Passing in 0 would put it immediately prior to the node, and any other numeric value would place it at that point in time from the present. The first set of parameters refers to the background'' group, the second refers refers to the foreground'' clade rooted at \code{nd153}. To simplify things, we will constrain drift to be zero and assume that both partitions have the same diffusion coefficient. <<>>= g.cc <- constrain(f.cc, drift.1 ~ 0, drift.2 ~ 0, diffusion.2 ~ diffusion.1) argnames(g.cc) @ Generate a starting point from the single partition ML point: <<>>= p.cc <- c(p.c, p.c[1:2]) names(p.cc) <- argnames(g.cc) @ At this point, the split function should have basically the same likelihood as the single partition function: <>= mle.c$lnLik - g.cc(p.cc) @ And run the ML search <>= mle.cc <- find.mle(g.cc, p.cc, control=control, lower=0, verbose=0) @ Repeat this for linear speciation functions: <<>>= f.ll <- make.quasse.split(phy, mass, mass.sd, linear.x, constant.x, "nd153", Inf) g.ll <- constrain(f.ll, drift.1 ~ 0, drift.2 ~ 0, diffusion.2 ~ diffusion.1) g.lc <- constrain(g.ll, l.m.2 ~ 0) g.cl <- constrain(g.ll, l.m.1 ~ 0) @ Generate a starting points: start with the function where both speciation rates are linear functions. <<>>= p.cc <- coef(mle.cc) p.ll <- c(p.cc[1], 0, p.cc[2:4], 0, p.cc[5]) names(p.ll) <- argnames(g.ll) @ Run the ML searches for this model: <>= mle.ll <- find.mle(g.ll, p.ll, control=control, verbose=0) @ Then generate starting points for models with just one of the sections of the tree having a linear speciation function: <<>>= p.lc <- c(coef(mle.ll)[1:3], p.ll[c(4, 5, 7)]) p.cl <- c(p.ll[c(1, 3, 4)], coef(mle.ll)[5:7]) @ and run the ML search: <>= mle.lc <- find.mle(g.lc, p.lc, control=control, verbose=0) mle.cl <- find.mle(g.cl, p.cl, control=control, verbose=0) @ We can then compare the models again: <<>>= anova(mle.c, linear=mle.l, sigmoidal=mle.s, hump=mle.h, part.constant=mle.cc, part.linear.bg=mle.lc, part.linear.fg=mle.cl, part.linear=mle.ll) @ This supports the model with a linear foreground'' rate of speciation (lowest AIC value). Looking at the coefficients for this model: <<>>= coef(mle.cl) @ The speciation rate in the foreground clade is a negative function of body size (\code{l.m.2} is negative) -- increasing body size decreases the speciation rate. % TODO: Demonstrate how to plot the fits, as in the paper, as there is some confusion here. \subsection{Controlling the calculations} There are a number of ways to control the behaviour of the integrator: most of these are a speed/accuracy trade-off. These should be supplied in the list passed in as the argument \code{control}. The \code{make.quasse} function attempts to select defaults that result in acceptable performance and accuracy, but this is not always possible. \begin{itemize} \item \code{method}: one of \code{fftC} or \code{fftR} to switch between C (fast) and R (slow) back-ends for the integration. Both use non-adaptive fft-based convolutions. Eventually, an adaptive methods-of-lines approach will be available. \item \code{dt.max}: Maximum time step to use for the integration. By default, this will be set to 1/1000 of the tree depth. Smaller values will slow down calculations, but improve accuracy. \item \code{nx}: The number of bins into which the character space is divided (default=1024). Larger values will be slower and more accurate. For the \code{fftC} integration method, this should be an integer power of 2 (512, 2048, etc). \item \code{r}: Scaling factor that multiplies \code{nx} for a "high resolution" section at the tips of the tree (default=4, giving a high resolution character space divided into 4096 bins). This helps improve accuracy while possibly tight initial probability distributions flatten out as time progresses towards the root. Larger values will be slower and more accurate. For the \code{fftC} integration method, this should be a power of 2 (2, 4, 8, so that \code{nx*r} is a power of 2). \item \code{tc}: where in the tree to switch to the low-resolution integration (zero corresponds to the present, larger numbers moving towards the root). By default, this happens at 10\% of the tree depth. Smaller values will be faster, but less accurate. \item \code{xmid}: Mid point to centre the character space. By default this is at the mid point of the extremes of the character states. \item \code{tips.combined}: Get a modest speed-up by simultaneously integrating all tips? By default, this is \code{FALSE}, but speedups of up to 25\% are possible with this set to \code{TRUE}. \item \code{w}: Number of standard deviations of the normal distribution induced by Brownian motion to use when doing the convolutions (default=5). Probably best to leave this one alone. \end{itemize} \section{Geographic distributions and diversification: GeoSSE (by Emma Goldberg)} \label{sec:geosse} The GeoSSE model (Geographic State Speciation and Extinction) combines features of the constant-rates birth-death model with a three-state Markov model. It differs from BiSSE in parameterising the model to represent diversification and range shifts among two regions, which includes allowing widely-distributed species whose ranges may change in conjunction with a speciation event. See \citet{Goldberg-2011-451} for a full explanation of the model assumptions. \subsection{Parameters and a tree} Simulating trees under GeoSSE can be done with the function \code{tree.geosse}; a simulated tree shown in Figure \ref{fig:geosse-phy}. <>= pars <- c(1.5, 0.5, 1.0, 0.7, 0.7, 2.5, 0.5) names(pars) <- diversitree:::default.argnames.geosse() set.seed(5) phy <- tree.geosse(pars, max.t=4, x0=0) @ <>= statecols <- c("AB"="purple", "A"="blue", "B"="red") plot(phy, tip.color=statecols[phy$tip.state+1], cex=0.5) @ \begin{figure} \centering <>= statecols <- c("AB"="purple", "A"="blue", "B"="red") plot(phy, tip.color=statecols[phy$tip.state+1], cex=0.5, no.margin=TRUE) @ \caption{A GeoSSE tree simulated with {\tt params = c(1.5, 0.5, 1.0, 0.7, 0.7, 2.5, 0.5)}. The tip state colors are purple for species present in both regions (AB), blue for species only in region A, and red for species only in region B.} \label{fig:geosse-phy} \end{figure} An extremely crude starting point for parameter estimation can be obtained with: <>= p <- starting.point.geosse(phy) p @ The parameters are: speciation within region A ({\tt sA}), speciation within region B ({\tt sB}), between-region speciation ({\tt sAB}), extinction from region A ({\tt xA}), extinction from region B ({\tt xB}), dispersal from A to B (range expansion, {\tt dA}), and dispersal from B to A ({\tt dB}). \subsection{Model construction and constraining} Constructing and constraining likelihood functions works as for the other models. Here we will consider the full model, a model without between-region speciation, and a model without regional dependence of speciation or extinction rates. <>= lik1 <- make.geosse(phy, phy$tip.state) lik2 <- constrain(lik1, sAB ~ 0) lik3 <- constrain(lik1, sA ~ sB, xA ~ xB) @ \subsection{Maximum likelihood} ML parameter estimation and model comparisons: <>= ml1 <- find.mle(lik1, p) p <- coef(ml1) ml2 <- find.mle(lik2, p[argnames(lik2)]) ml3 <- find.mle(lik3, p[argnames(lik3)]) <>= round(rbind(full = coef(ml1), no.sAB = coef(ml2, TRUE), eq.div = coef(ml3, TRUE)), 3) anova(ml1, no.sAB = ml2, eq.div = ml3) @ On this tree, we reject the model of equal speciation and extinction in the two regions, concluding that there are regional differences in diversification. Including the between-region mode of speciation does not, however, significantly improve the fit. \subsection{Markov chain Monte Carlo} We will only consider the 6-parameter model here. Use the ML rate estimates as a starting point. Place a broad exponential prior on each parameter. <>= p <- coef(ml2) prior <- make.prior.exponential(1/2) @ Use a pilot run to obtain reasonable step sizes: <>= set.seed(1) tmp <- mcmc(lik2, p, nsteps=100, prior=prior, w=1, print.every=0) w <- diff(sapply(tmp[2:7], quantile, c(0.025, 0.975))) @ Now the real analysis, which will take awhile to run: <>= mcmc2 <- mcmc(lik2, p, nsteps=10000, prior=prior, w=w) @ <>= .save <- "cache/geosse-mcmc-samples.Rdata" if ( file.exists(.save) ) { load(.save) } else { <> save(mcmc2, file=.save) } @ Marginal posterior distributions are shown in Figure \ref{fig:geosse-mcmc}. We can also compare rate estimates by looking at posterior probabilities of their differences. <>= mcmc2diff <- with(mcmc2, data.frame(s.diff=sA-sB, x.diff=xA-xB, d.diff=dA-dB, div.A=sA-xA, div.B=sB-xB)) colMeans(mcmc2diff > 0) @ We correctly and confidently recover regional biases in speciation and dispersal, and positive net diversification in region A. With less confidence, we recover regional bias in extinction and negative net diversification for region B (i.e., extinction exceeds speciation in region B). \begin{figure} \centering <>= col1 <- c("red", "orange", "blue", "purple", "black", "gray") col2 <- col1[c(1,3,5)] mcmc2diff <- with(mcmc2, data.frame(s.diff=sA-sB, x.diff=xA-xB, d.diff=dA-dB)) par(mfrow=c(2,1), mar=c(3, 4, 0, 1)) profiles.plot(mcmc2[2:7], col.line=col1, xlab="", ylab="") legend("topright", argnames(lik2), col=col1, lty=1) profiles.plot(mcmc2diff, col.line=col2, xlab="", ylab="") legend("topright", colnames(mcmc2diff), col=col2, lty=1) title(xlab="rate", ylab="posterior probability density", outer=T, line=-1) @ \caption{Posterior probability distributions for the six-rate GeoSSE model, for the tree shown in Figure \ref{fig:geosse-phy}. Uncertainty is largest for dispersal and smallest for speciation. Regional differences in speciation and, to some extent, in dispersal are recovered.} \label{fig:geosse-mcmc} \end{figure} \subsection{Additional options} GeoSSE likelihood functions can be built with randomly-incomplete sampling \citep{FitzJohn-2009-595}. There is not currently support for unresolved clades. <>= p <- coef(ml1) lik1(p) lik4 <- make.geosse(phy, phy$tip.state, sampling.f=c(0.9, 0.6, 0.4)) lik4(p) @ When using the likelihood function, one can condition on survival of the clade (not done by default): <>= lik4(p, condition.surv=TRUE) @ External information about the geographic distribution of the common ancestor of the clade can be enforced by fixing the root state. For example, if you are absolutely positive that the MRCA was found only in region B: <>= lik4(p, root.p=c(0,0,1), root=ROOT.GIVEN) @ Use this procedure with caution, and only in the face of truly external data, e.g., fossil or geologic information. \subsection{Multi-clade analysis} Some applications of GeoSSE have combined multiple clades into a single analysis \citep{Anacker-2010-365,Goldberg-2011-451}. This has the advantage of providing a larger dataset and hence presumably more power, but it is important to keep in mind the assumptions that go into such an analysis. First, it treats all clades as evolving according to the same model, with the same values for the rate parameters. This can be tested by first fitting the clades individually, or it may be inherent to the hypothesis at hand. Second, it assumes that the clades are independent of each other, since their likelihoods are simply multiplied together to form a joint likelihood function. In some cases, this may be a defensible approximation, for example, when the clades being considered are only very distantly related. If you have decided that you can convince yourself and your reviewers that a multi-clade analysis is appropriate, here is one way to do it. % This could all be put into a multi.clade() function, but I am wary % of making it too easy to use. Assemble the data as lists of trees and character state vectors. Each list element is one clade. Here we will use two trees from the chaparral study, one each from the posterior sets for {\em Ceanothus} and {\em Arctostaphylos}. The trees {\tt phy.cea} and {\tt phy.arc} and the tip states {\tt chars.cea} and {\tt chars.arc} (A = chaparral, B = forest) were read in with the call to {\tt data("geosse")} above. <<>>= phy.cea <- read.tree(file="data/cea.tre") phy.arc <- read.tree(file="data/arc.tre") temp <- read.csv("data/cea.dat", header=FALSE, as.is=TRUE) chars.cea <- structure(temp[[2]], names=temp[[1]]) temp <- read.csv("data/arc.dat", header=FALSE, as.is=TRUE) chars.arc <- structure(temp[[2]], names=temp[[1]]) rm(temp) @ Create an individual likelihood function for each clade: <>= trees <- list(phy.cea, phy.arc) states <- list(chars.cea, chars.arc) sampl <- list(c(0.913, 0.941, 0.875), c(0.674, 0.533, 0.750)) liks <- mapply(make.geosse, trees, states, sampl) @ Now define the joint, multi-clade likelihood function: <<>>= lnL.multi <- combine(liks) @ Now {\tt lnL.multi()} can be used like any other likelihood function, for example, it can be constrained: <<>>= lnL.multi.constrained <- constrain(lnL.multi, sAB ~ 0) @ evaluated at arbitrary points: <>= p <- c(0.19, 0.08, 0.00, 0.29, 0.48, 1.29, 0.87) lnL.multi(p) @ or used in {\tt find.mle()} or {\tt mcmc()}. \section{Priors for MCMC} \label{sec:priors-mcmc} Above, the prior used for MCMC analysis was exponential with a mean of twice the net diversification rate for all parameters. However, any prior can be used; the only requirement is that a prior function takes as its only argument a vector of parameters, and returns as output the log prior probability. For demonstration, here is the tree from section \ref{sec:bisse}: <<>>= pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01) set.seed(4) phy <- tree.bisse(pars, max.t=30, x0=0) states <- phy$tip.state @ % $% for emacs The state-independent diversification rate can be computed as <<>>= p <- starting.point.bd(phy) r <- p[[1]] - p[[2]] r @ and we can use \code{make.prior.exponential} to make a prior with a mean of twice this rate as <<>>= prior <- make.prior.exponential(1/(2*r)) @ % Note that R's parametrisation of the exponential distribution (which diversitree uses) uses the rate of an exponential distribution, not the mean, which is why we take$1/\mathrm{mean}$. % The definition of \code{make.prior.exponential} is <<>>= make.prior.exponential @ % This is functional programming; the function \code{make.prior.exponential} takes an argument \code{r}, and \emph{returns a function} which takes an argument \code{pars}. This returned function (which we saved as \code{prior} above) uses its argument \code{pars} along with the previously specified rate \code{r} to compute the likelihood. The rate passed to \code{make.prior.exponential} need not be a scalar; it can also be the same length as \code{prior} (it could be other lengths, but I would not do this as results could be surprising). For example, suppose that we feel that the character transition rate parameters$q_{01}$and$q_{10}$(elements 5 and 6 in \code{pars}) should be lower than the speciation and extinction rates, with a mean of$r/4$, rather than$2r$. We can make a prior function that reflects this. First, make a vector of means where the first four elements ($\lambda_0$,$\lambda_1$,$\mu_0$, and$\mu_1$) have mean$2r$and the second two elements have mean$r/4$; <<>>= prior.means <- rep(c(2*r, r/4), c(4, 2)) prior.means @ then pass \code{1/prior.means} to \code{make.prior.exponential}: <<>>= prior.lowq <- make.prior.exponential(1/prior.means) @ Other distributions can be used, though you will have to write your own prior functions. For example, to use a uniform distribution on [min,max], we can use R's \code{dunif} to compute the log probability density of some parameters. The format is: <>= dunif(parameters, min, max, log=TRUE) @ % (demonstration only --- this will not run). We can make a generator function for a prior distribution where the first four parameters are exponentially distributed and the last two are uniformly distributed between \code{min} and \code{max}: <<>>= make.prior.expuniform <- function(r, min, max) { function(pars) sum(dexp(pars[1:4], r, log=TRUE)) + sum(dunif(pars[5:6], min, max, log=TRUE)) } @ Using this with the usual prior on$\lambda$and$\mu$, and letting the$q$parameters be uniformly distributed on$[0,0.01]\$ we would write: <<>>= prior.expu <- make.prior.expuniform(1/(2*r), 0, 0.01) @ \section{Topics not covered} I have not covered several included models here, but information can be found in the online documentation. \begin{itemize} % \item MEDUSA-style split'' models. Following % \citep{Alfaro-2009-13410}, these models allow different regions of % the tree to have different parameters. This is implemented for % birth-death models (\code{make.bd.split}, which is identical to % MEDUSA), BiSSE (\code{?make.bisse.split}), and MuSSE % (\code{?make.musse.split}). \item Time-dependent models. In these, time is divided into epochs'', each of which may have different parameters. This is implemented for BiSSE (\code{?make.bisse.td}) and MuSSE (\code{?make.musse.td}). Another model of time-dependence is also possible, where rates are arbitrary functions of time: implemented for plain birth-death models (\code{?make.bd.t}), BiSSE (\code{?make.bisse.td}) and MuSSE (\code{?make.musse.td}). \item Brownian motion and Ornstein-Uhlenbeck. A very simple-minded Brownian motion likelihood calculation is included. This allows estimation of the diffusion parameter of a Brownian motion process for the evolution of a single continuous trait (\code{?make.bm}). \item Stochastic character mapping. Currently implemented only for Markov models of discrete character evolution (?asr.stoch.mkn). \end{itemize} \clearpage \bibliographystyle{refstyle} \bibliography{refs} \end{document} % LocalWords: Bollback ip Gorjanc's Sweave ld Rnw bd phy lik argnames surv mle % LocalWords: nlm coef logLik lnLik anova lrt mcmc nsteps las ylab xlab abline % LocalWords: topright lty mk eb rbind tmp sapply setdiff len csv dimorph plt % LocalWords: tipmap yy loc cbind rowSums cumsum rect mfrow cex sqrt rle ym Nc % LocalWords: Chionidae Thinocoridae ifelse Dromas nex musse quasse sd nodrift % LocalWords: parscale xmid xr lc sc noroptimal reltol xmin xmax mgp lim dx sA % LocalWords: lastPP phylo envir PlotPhyloEnv Ntip uniroot usr dy xpd ln nd dt % LocalWords: Tarsiidae Hominoidea Cercopithecoidea superfamily Nnode sep GHz % LocalWords: fftC fftR fft nx tc bm diversitree Pagel Univariate Ornstein pos % LocalWords: Uhlenbeck speciation Maddison phylogenies FitzJohn MuSSE QuaSSE % LocalWords: GeoSSE BiSSEness Magnuson ClaSSE Cladogenetic classe practice sB % LocalWords: online MacPro cacheSweave taxon clades taxa subplex distr eaab % LocalWords: sim autocorrelation speciates bisse clade supertree % LocalWords: Lislevand polytomies gray polytomy monomorphic eval sigmoidal xA % LocalWords: color Vos Redding Alfaro cmp lccl AIC parameterising simtreesdd % LocalWords: loadlib geosse plottree statecols params colors parstart sAB xB % LocalWords: dA mlans eq