diversitree (version 0.10-1)

make.geosse: Geographic State Speciation and Extinction Model

Description

Prepare to run GeoSSE (Geographic State Speciation and Extinction) on a phylogenetic tree and character distribution. This function creates a likelihood function that can be used in maximum likelihood or Bayesian inference.

Usage

make.geosse(tree, states, sampling.f=NULL, strict=TRUE,
    control=list())
  starting.point.geosse(tree, eps=0.5)

Arguments

tree

An ultrametric bifurcating phylogenetic tree, in ape “phylo” format.

states

A vector of character states, each of which must be 0 (in both regions/widespread; AB), 1 or 2 (endemic to one region; A or B), or NA if the state is unknown. This vector must have names that correspond to the tip labels in the phylogenetic tree (tree$tip.label).

sampling.f

Vector of length 3 with the estimated proportion of extant species in states 0, 1 and 2 that are included in the phylogeny. A value of c(0.5, 0.75, 1) means that half of species in state 0, three quarters of species in state 1, and all the species in state 2 are included in the phylogeny. By default all species are assumed to be known.

strict

The states vector is always checked to make sure that the values are 0, 1 and 2 only. If strict is TRUE (the default), then the additional check is made that every state is present. The likelihood models tend to be poorly behaved where states are missing.

control

List of control parameters for the ODE solver. See details in make.bisse.

eps

Ratio of extinction to speciation rates to be used when choosing a starting set of parameters. The procedure used is based on Magallon & Sanderson (2001).

Author

Emma E. Goldberg

Details

make.geosse returns a function of class geosse. The arguments and default values for this function are:

    f(pars, condition.surv=TRUE, root=ROOT.OBS, root.p=NULL,
      intermediates=FALSE)
  

The arguments of this function are explained in make.bisse. The parameter vector pars is ordered sA, sB, sAB, xA, xB, dA, dB.

Unresolved clade methods are not available for GeoSSE. With three states, it would rapidly become computationally infeasible.

References

FitzJohn R.G., Maddison W.P., and Otto S.P. 2009. Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies. Syst. Biol. 58:595-611.

Goldberg E.E., Lancaster L.T., and Ree R.H. 2011. Phylogenetic inference of reciprocal effects between geographic range evolution and diversification. Syst. Biol. 60:451-465.

Maddison W.P., Midford P.E., and Otto S.P. 2007. Estimating a binary character's effect on speciation and extinction. Syst. Biol. 56:701-710.

Magallon S. and Sanderson M.J. 2001. Absolute diversification rates in angiospem clades. Evol. 55:1762-1780.

See Also

constrain for making submodels, find.mle for ML parameter estimation, mcmc for MCMC integration, make.bisse for further relevant examples.

The help page for find.mle has further examples of ML searches on full and constrained BiSSE models. Things work similarly for GeoSSE, just with different parameters.

Examples

## Due to a change in sample() behaviour in newer R it is necessary to
## use an older algorithm to replicate the previous examples
if (getRversion() >= "3.6.0") {
  RNGkind(sample.kind = "Rounding")
}

## Parameter values
pars <- c(1.5, 0.5, 1.0, 0.7, 0.7, 2.5, 0.5)
names(pars) <- diversitree:::default.argnames.geosse()

## Simulate a tree
set.seed(5)
phy <- tree.geosse(pars, max.t=4, x0=0)

## See the data
statecols <- c("AB"="purple", "A"="blue", "B"="red")
plot(phy, tip.color=statecols[phy$tip.state+1], cex=0.5)

## The likelihood function
lik <- make.geosse(phy, phy$tip.state)

## With "true" parameter values
lik(pars) # -168.4791

## A guess at a starting point.
p <- starting.point.geosse(phy)

## Start an ML search from this point (takes a couple minutes to run).
if (FALSE) {
fit <- find.mle(lik, p, method="subplex")
logLik(fit) # -165.9965

## Compare with sim values.
rbind(real=pars, estimated=round(coef(fit), 2))

## A model with constraints on the dispersal rates.
lik.d <- constrain(lik, dA ~ dB)
fit.d <- find.mle(lik.d, p[-7])
logLik(fit.d) # -166.7076

## A model with constraints on the speciation rates.
lik.s <- constrain(lik, sA ~ sB, sAB ~ 0)
fit.s <- find.mle(lik.s, p[-c(2,3)])
logLik(fit.s) # -169.0123
}

## "Skeletal tree" sampling is supported.  For example, if your tree
## includes all AB species, half of A species, and a third of B species,
## create the likelihood function like this:
lik.f <- make.geosse(phy, phy$tip.state, sampling.f=c(1, 0.5, 1/3))

## If you have external evidence that the base of your tree must have
## been in state 1, say (endemic to region A), you can fix the root 
## when computing the likelihood, like this:
lik(pars, root=ROOT.GIVEN, root.p=c(0,1,0))