tess.analysis: tess.analysis: Diversification rate estimation under an...

In hoehna/TESS: Diversification Rate Estimation and Fast Simulation of Reconstructed Phylogenetic Trees under Tree-Wide Time-Heterogeneous Birth-Death Processes Including Mass-Extinction Events

Description

tess.analysis estimates diversification rates under an episodic birth-death process including mass-extinction events. The method uses a reversible-jump MCMC algorithm to estimate the number, timing and magnitude of rate-shifts and mass-extinction events. It is possible to fix either number of events and provide specific values that will be used. We assume a Poison process for the number of events and a lognormal distribution with fixed, but specified, hyper-parameters for the speciation and extinction rate; and an independent Poison process for the number of mass-extinction events where each survival probability follows a Beta distribution with fixed hyper-parameters.

The MCMC algorithm can be run either for a specified number of iterations, until a time limit in seconds has been reached, or until the effective sample size (ESS) has reached a given threshold. Once the first of these requirements are met TESS will stop the analysis. Internally we use scaling and sliding proposals to change the parameter values during the MCMC and a birth-move and death-move to add/remove events (rate-shifts or mass-extinction events).

The results of the MCMC run are stored within a directory that is specified by the user. Several files will be generated containing the sampled parameter values. To summarize the output see tess.process.output(...) and tess.plot.output(...).

Usage

tess.analysis( tree,
               initialSpeciationRate,
               initialExtinctionRate,
               empiricalHyperPriors = TRUE,
               empiricalHyperPriorInflation = 10.0,
               empiricalHyperPriorForm = c("lognormal","normal","gamma"),
               speciationRatePriorMean = 0.0,
               speciationRatePriorStDev = 1.0,
               extinctionRatePriorMean = 0.0,
               extinctionRatePriorStDev = 1.0,
               initialSpeciationRateChangeTime = c(),
               initialExtinctionRateChangeTime = c(),
               estimateNumberRateChanges = TRUE,
               numExpectedRateChanges = 2,
               samplingProbability = 1,
       				 samplingStrategy = "uniform",
               missingSpecies = c(),
               timesMissingSpecies = c(),
               tInitialMassExtinction = c(),
               pInitialMassExtinction = c(),
               pMassExtinctionPriorShape1 = 5,
               pMassExtinctionPriorShape2 = 95,
               estimateMassExtinctionTimes = TRUE,
               numExpectedMassExtinctions = 2,
               estimateNumberMassExtinctions = TRUE,
               MRCA = TRUE,
               CONDITION = "survival",
               BURNIN = 10000,
               MAX_ITERATIONS = 200000,
               THINNING = 100,
               OPTIMIZATION_FREQUENCY = 500,
               CONVERGENCE_FREQUENCY = 1000,
               MAX_TIME = Inf, MIN_ESS = 500,
               ADAPTIVE = TRUE,
               dir = "" ,
               priorOnly = FALSE,
               verbose = TRUE)

Arguments

tree	The tree in 'phylo' format.
initialSpeciationRate	The initial value of the speciation rate when the MCMC is started. This can either be a single number of a vector of rates per interval.
initialExtinctionRate	The initial value of the extinction rate when the MCMC is started. This can either be a single number of a vector of rates per interval.
empiricalHyperPriors	Should we estimate the hyper-parameters empirically?
empiricalHyperPriorInflation	The scaling factor of the variance for the empirical hyperpriors.
empiricalHyperPriorForm	The possible empirical hyper prior distributions; either lognormal, normal or gamma
speciationRatePriorMean	The mean of the log-normal prior distribution for the speciation rate.
speciationRatePriorStDev	The standard deviation of the log-normal prior distribution for the speciation rate.
extinctionRatePriorMean	The mean of the log-normal prior distribution for the extinction rate.
extinctionRatePriorStDev	The standard deviation of the log-normal prior distribution for the extinction rate.
initialSpeciationRateChangeTime	The initial value of the time points when speciation rate-shifts occur. The number of time-shifts needs to be one smaller than the number of initial speciation rates.
initialExtinctionRateChangeTime	The initial value of the time points when extinction rate-shifts occur. The number of time-shifts needs to be one smaller than the number of initial extinction rates.
estimateNumberRateChanges	Do we estimate the number of rate shifts? Default is true.
numExpectedRateChanges	Expected number of rate changes which follow a Poisson process. The default gives 0.5 probability on 0 shifts.
samplingProbability	The extant taxa sampling probability at the present time.
samplingStrategy	The strategy how samples were obtained. Options are: uniform|diversified|age.
missingSpecies	The number of species missed which originated in a given time interval (empirical taxon sampling).
timesMissingSpecies	The times intervals of the missing species (empirical taxon sampling).
tInitialMassExtinction	The initial value of the vector of times of the mass-extinction events. This is used as initial values for the MCMC.
pInitialMassExtinction	The initial value of the vector of survival probabilities of the mass-extinction events. This is used as initial values for the MCMC.
pMassExtinctionPriorShape1	The alpha (first shape) parameter of the Beta prior distribution for the survival probability of a mass-extinction event.
pMassExtinctionPriorShape2	The beta (second shape) parameter of the Beta prior distribution for the survival probability of a mass-extinction event.
estimateMassExtinctionTimes	Do we estimate the times of mass-extinction events? Default is true.
numExpectedMassExtinctions	Expected number of mass-extinction events which follow a Poisson process. The default gives 0.5 probability on 0 events.
estimateNumberMassExtinctions	Do we estimate the number of mass-extinction events? Default is true.
MRCA	Does the process start with the most recent common ancestor? If not, the tree must have a root edge!
CONDITION	do we condition the process on time|survival|taxa?
BURNIN	The length of the burnin period.
MAX_ITERATIONS	The maximum number of iteration of the MCMC. The default is 200000.
THINNING	The frequency how often samples are recorded during the MCMC. The default is every 100 iterations.
OPTIMIZATION_FREQUENCY	The frequency how often the MCMC moves are optimized. The default is every 500 iterations.
CONVERGENCE_FREQUENCY	The frequency how often we check for convergence? The default is every 1000 iterations.
MAX_TIME	The maximum time the MCMC is allowed to run in seconds. The default is Inf
MIN_ESS	The minimum number of effective samples (ESS) to assume convergence. The default is 500
ADAPTIVE	Do we use auto-tuning of the MCMC moves? The default is TRUE (recommended).
dir	The subdirectory in which the output will be stored. The default is the present directoy ("")
priorOnly	Do we sample from the prior only? The default is FALSE
verbose	Do you want detailed output?

Value

There is no return value because all the results are stored into files.

Author(s)

Sebastian Hoehna

References

S. Hoehna: The time-dependent reconstructed evolutionary process with a key-role for mass-extinction events. 2015, Journal of Theoretical Biology, 380, 321-331.

S. Hoehna, MR May and BR Moore: TESS: Bayesian inference of lineage diversification rates from (incompletely sampled) molecular phylogenies in R. 2015, Bioinformatics.

MR May, S. Hoehna, and BR Moore: A Bayesian approach for detecting mass-extinction events when rates of lineage diversification vary. 2015, Systematic Biology

Examples

# we load the conifers as the test data set
data(conifers)

# for the conifers we know what the total number of species is
total <- 630
# thus, we can compute what the sampling fraction is
rho <- (conifers$Nnode+1)/total


# next, we specify the prior mean and standard deviation
# for the speciation and extinction rate
mu_lambda = 0.15
std_lambda = 0.02
mu_mu = 0.09
std_mu = 0.02

# now we can run the entire analysis.
# note that a full analyses should be run much longer
tess.analysis( tree=conifers,
               initialSpeciationRate=exp(mu_lambda),
               initialExtinctionRate=exp(mu_mu),
               empiricalHyperPriors = FALSE,
               speciationRatePriorMean = mu_lambda,
               speciationRatePriorStDev = std_lambda,
               extinctionRatePriorMean = mu_mu,
               extinctionRatePriorStDev = std_mu,
               numExpectedRateChanges = 2,
               samplingProbability = rho,
               numExpectedMassExtinctions = 2,
               BURNIN = 100,
               MAX_ITERATIONS = 200,
               THINNING = 10,
               dir = "analysis_conifer")

# You may want to look into the vignette for a more detailed description
# of the features for an analysis.
# also have a look at the functions tess.process.output and tess.plot.output

hoehna/TESS documentation built on Feb. 3, 2022, 5:59 a.m.