make.geosse {diversitree}  R Documentation 
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.
make.geosse(tree, states, sampling.f=NULL, strict=TRUE, control=list()) starting.point.geosse(tree, eps=0.5)
tree 
An ultrametric bifurcating phylogenetic tree, in

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 
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 
strict 
The 
control 
List of control parameters for the ODE solver. See
details in 
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). 
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.
Emma E. Goldberg
FitzJohn R.G., Maddison W.P., and Otto S.P. 2009. Estimating traitdependent speciation and extinction rates from incompletely resolved phylogenies. Syst. Biol. 58:595611.
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:451465.
Maddison W.P., Midford P.E., and Otto S.P. 2007. Estimating a binary character's effect on speciation and extinction. Syst. Biol. 56:701710.
Magallon S. and Sanderson M.J. 2001. Absolute diversification rates in angiospem clades. Evol. 55:17621780.
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.
Parameter values
pars < c(1.5, 0.5, 1, 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
## [1] 168.5
A guess at a starting point.
p < starting.point.geosse(phy)
Start an ML search from this point (takes a couple minutes to run).
[This section is not run by default]
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
[Ends not run section]
"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))
## [1] 168.3