Rate-Variable Birth Death {laser}R Documentation

Fit Rate-Variable Diversification Model With Extinction to Branching Times

Description

Fits a two-rate variant of the birth death model to branching times derived from phylogenetic data. The model assumes that a clade has diversified under a net diversification rate r1 until some time st, at which point the net diversification rate shifts to r2. The extinction fraction a is assumed constant and can take values from [0, 1). The shift time st giving the maximum log-likelihood is found by sequentially optimizing r1 and r2 across a large number of possible shift times.

Usage

rvbd(x, ai = c(0.1, 0.5, 0.95), ints = NULL, verbose = FALSE, file = "out_rvbd.txt")

Arguments

x a numeric vector of branching times
ai a vector of initial a values for optimization of the likelihood function
ints the number of intervals. See details
verbose if 'verbose = TRUE', writes likelihoods and parameter estimates for all shift points considered to the indicated file.
file a filename for output if 'verbose = TRUE'

Details

Non-linear optimization can be notoriously difficult, and the functions used in fitting this model can be susceptible to the presence of multiple optima in the likelihood surface. To guard against the possibility, the model is fitted using 3 initial values of the extinction fraction, a. You should check results obtained against those obtained using the yule2rate model. If the maximum log-likelihood under rvbd is less than that obtained using yule2rate, your analysis has been fooled by a local (sub)optimum. Repeat the analysis using additional starting values for a, e.g., 'ai = seq(0, .99. length.out = 15)'.

'ints' is used in determining the number of shift points to consider. If 'ints = NULL' (the default), the model will consider only observed branching times as possible shift points. See yule-n-rate for additional discussion of the 'ints' option.

'verbose = TRUE' will write maximum log-likelihoods and parameter estimates for each shift time under consideration to . The file can then be loaded to examine the likelihood of a rate shift at different points in time.

Value

a list with the following components:

LH the log-likelihood at the maximum
aic the Akaike Information Criterion
r1 the initial net diversification rate
r2 the final net diversification rate
a the extinction fraction
st the shift time giving the maximum log-likelihood

Note

Note that the net diversification rate r is equal to the speciation rate S minus the extinction rate E. The extinction fraction a = E/S. The model assumes a constant extinction fraction a, but this is not the same as a constant extinction rate E.

The model contains four free parameters: the initial and final net diversification rates, the extinction fraction, and the shift time.

Note that shift times, like branching times, are given in divergence units before present. Thus, if you have scaled a set of branching times to a basal divergence of 30 million years before present, you would interpret 'st1 = 19.5' as an inferred shift point 19.5 million years before present.

Large numbers of intervals or large sets of initial a values can result in high computational time

Author(s)

Dan Rabosky DLR32@cornell.edu

References

Barraclough, T. G., and A. P. Vogler. 2002. Recent diversification rates in North American tiger beetles estimated from a dated mtDNA phylogenetic tree. Mol. Biol. Evol. 19:1706-1716.

Nee, S., R. M. May, and P. H. Harvey. 1994b. The reconstructed evolutionary process. Philos. Trans. R. Soc. Lond. B 344:305-311.

Rabosky, D. L. 2006. Likelihood methods for inferring temporal shifts in diversification rates. Evolution 60:1152-1164.

See Also

yule-n-rate, fitdAICrc

Examples

 
  data(plethodon)
  result <- rvbd(plethodon, ints = 100, verbose = TRUE, file = 'temp.txt')
  
  ### now to generate a plot of log-likelihoods against shift times:
  temp <- read.table(file = 'temp.txt', header = TRUE)
  
  ### Rescaling shift times to reflect time since basal divergence
  stvec <- plethodon[1] - temp$st
  plot(temp$LH~stvec, xlab = "Time From Basal Divergence", 
        ylab = "Log-likelihood")
  
  unlink('temp.txt')  #clean-up temp file.      
  

[Package laser version 2.1 Index]