| cal3TimePaleoPhy {paleotree} | R Documentation |
Three Rate Calibrated a posteriori Dating of Paleontological Phylogenies
Description
Time-scales an undated cladogram of fossil taxa, using information on their ranges and estimates of the instantaneous rates of branching, extinction and sampling. The output is a sample of a posteriori time-scaled trees, as resulting from a stochastic algorithm which samples observed gaps in the fossil record with weights calculated based on the input rate estimates. This function also uses the three-rate calibrated dating algorithm to stochastically resolve polytomies and infer potential ancestor-descendant relationships, simultaneous with the time-scaling treatment.
Usage
cal3TimePaleoPhy(
tree,
timeData,
brRate,
extRate,
sampRate,
ntrees = 1,
anc.wt = 1,
node.mins = NULL,
dateTreatment = "firstLast",
FAD.only = FALSE,
adj.obs.wt = TRUE,
root.max = 200,
step.size = 0.1,
randres = FALSE,
noisyDrop = TRUE,
verboseWarnings = TRUE,
diagnosticMode = FALSE,
tolerance = 1e-04,
plot = FALSE
)
bin_cal3TimePaleoPhy(
tree,
timeList,
brRate,
extRate,
sampRate,
ntrees = 1,
anc.wt = 1,
node.mins = NULL,
dateTreatment = "firstLast",
FAD.only = FALSE,
sites = NULL,
point.occur = FALSE,
nonstoch.bin = FALSE,
adj.obs.wt = TRUE,
root.max = 200,
step.size = 0.1,
randres = FALSE,
noisyDrop = TRUE,
verboseWarnings = TRUE,
tolerance = 1e-04,
diagnosticMode = FALSE,
plot = FALSE
)
Arguments
tree |
An unscaled cladogram of fossil taxa, of class |
timeData |
Two-column matrix of first and last occurrences in absolute
continuous time, with row names as the taxon IDs used on the tree. This means the
first column is very precise FADs (first appearance dates) and the second
column is very precise LADs (last appearance dates), reflect the precise points
in time when taxa first and last appear. If there is stratigraphic uncertainty in
when taxa appear in the fossil record, it is preferable to use the |
brRate |
Either a single estimate of the instantaneous rate of branching (also known as the 'per-capita' origination rate, or speciation rate if taxonomic level of interest is species) or a vector of per-taxon estimates |
extRate |
Either a single estimate of the instantaneous extinction rate (also known as the 'per-capita' extinction rate) or a vector of per-taxon estimates |
sampRate |
Either a single estimate of the instantaneous sampling rate or a vector of per-taxon estimates |
ntrees |
Number of dated trees to output. |
anc.wt |
Weighting against inferring ancestor-descendant relationships.
The argument |
node.mins |
The minimum dates of internal nodes (clades) on a phylogeny can be set
using |
dateTreatment |
This argument controls the interpretation of A second option is A third option is With both arguments Note that |
FAD.only |
Should the tips represent observation times at the start of
the taxon ranges? |
adj.obs.wt |
If the time of observation of a taxon is before the last appearance of that taxon,
should the weight of the time of observation be adjusted to account for the
known observed history of the taxon which occurs after the time of observation?
If so, then set |
root.max |
Maximum time before the first FAD that the root can be pushed back to. |
step.size |
Step size of increments used in zipper algorithm to assign node ages. |
randres |
Should polytomies be randomly resolved using |
noisyDrop |
If |
verboseWarnings |
if |
diagnosticMode |
If |
tolerance |
Acceptable amount of shift in tip dates from dates listed
in |
plot |
If true, plots the input, "basic" time-scaled phylogeny (an intermediary step in the algorithm) and the output cal3 time-scaled phylogeny. |
timeList |
A list composed of two matrices giving interval times and
taxon appearance dates. The rownames of the second matrix should be the taxon IDs,
identical to the |
sites |
Optional two column matrix, composed of site IDs for taxon FADs
and LADs. The sites argument allows users to constrain the placement of
dates by restricting multiple fossil taxa whose FADs or LADs are from the
same very temporally restricted sites (such as fossil-rich Lagerstatten) to
always have the same date, across many iterations of time-scaled trees. To
do this, provide a |
point.occur |
If true, will automatically produce a 'sites' matrix which forces all FADs and LADs to equal each other. This should be used when all taxa are only known from single 'point occurrences', i.e. each is only recovered from a single bed/horizon, such as a Lagerstatten. |
nonstoch.bin |
If |
Details
The three-rate calibrated ("cal3") algorithm time-scales trees a posteriori by stochastically picking node divergence times relative to a probability distribution of expected waiting times between speciation and first appearance in the fossil record. This algorithm is extended to apply to resolving polytomies and designating possible ancestor-descendant relationships. The full details of this method are provided in Bapst (2013, MEE).
Briefly, cal3 time-scaling is done by examining each node separately, moving from the root upwards. Ages of branching nodes are constrained below by the ages of the nodes below them (except the root; hence the need for the root.max argument) and constrained above by the first appearance dates (FADs) of the daughter lineages. The position of the branching event within this constrained range implies different amounts of unobserved evolutionary history. cal3 considers a large number of potential positions for the branching node (the notation in the code uses the analogy of viewing the branching event as a 'zipper') and calculates the summed unobserved evolutionary history implied by each branching time. The probability density of each position is then calculated under a gamma distribution with a shape parameter of 2 (implying that it is roughly the sum of two normal waiting times under an exponential) and a rate parameter which takes into account both the probability of not observing a lineage of a certain duration and the 'twiginess' of the branch, i.e. the probability of having short-lived descendants which went extinct and never were sampled (similar to Friedman and Brazeau, 2011). These densities calculated under the gamma distribution are then used as weights to stochastically sample the possible positions for the branching node. This basic framework is extended to polytomies by allowing a branching event to fall across multiple potential lineages, adding each lineage one by one, from earliest appearing to latest appearing (the code notation refers to this as a 'parallel zipper').
As with many functions in the paleotree library, absolute time is always decreasing, i.e. the present day is zero.
These functions will intuitively drop taxa from the tree with NA for range
or that are missing from timeData.
The sampling rate used by cal3 methods is the instantaneous sampling rate,
as estimated by various other function in the paleotree package. See
make_durationFreqCont for more details.
If you have the per-time unit sampling
probability ('R' as opposed to 'r') look at the sampling parameter
conversion functions also included in this package
(e.g. sProb2sRate). Most datasets will probably use
make_durationFreqDisc and sProb2sRate
prior to using this function, as shown in an example below.
The branching and extinction rate are the 'per-capita' instantaneous
origination/extinction rates for the taxic level of the tips of the tree
being time-scaled. Any user of the cal3 time-scaling method has multiple
options for estimating these rates. One is to separately calculate the
per-capita rates (following the equations in Foote, 2001) across multiple
intervals and take the mean for each rate. A second, less preferred option,
would be to use the extinction rate calculated from the sampling rate above
(under ideal conditions, this should be very close to the mean 'per-capita'
rate calculated from by-interval FADs and LADs). The branching rate in this
case could be assumed to be very close to the extinction rate, given the
tight relationship observed in general between these two (Stanley, 1976; see
Foote et al., 1999, for a defense of this approach), and thus the extinction
rate estimate could be used also for the branching rate estimate. (This is
what is done for the examples below.) A third option for calculating all
three rates simultaneously would be to apply likelihood methods developed by
Foote (2002) to forward and reverse survivorship curves. Note that only one
of these three suggested methods is implemented in paleotree: estimating the
sampling and extinction rates from the distribution of taxon durations via
make_durationFreqCont and make_durationFreqDisc.
By default, the cal3 functions will consider that ancestor-descendant
relationships may exist among the given taxa, under a budding cladogenetic
or anagenetic modes. Which tips are designated as which is given by two
additional elements added to the output tree,
$budd.tips (taxa designated as ancestors via budding cladogenesis) and
$anag.tips (taxa designated as ancestors via anagenesis).
This can be turned off by setting anc.wt = 0. As
this function may infer anagenetic relationships during time-scaling, this
can create zero-length terminal branches in the output. Use
dropZLB to get rid of these before doing analyses of lineage
diversification.
Unlike timePaleoPhy, cal3 methods will always resolve polytomies. In
general, this is done using the rate calibrated algorithm, although if
argument randres = TRUE, polytomies will be randomly resolved with uniform
probability, similar to multi2di from ape. Also, cal3 will always add the terminal
ranges of taxa. However, because of the ability to infer potential
ancestor-descendant relationships, the length of terminal branches may be
shorter than taxon ranges themselves, as budding may have occurred during
the range of a morphologically static taxon. By resolving polytomies with
the cal3 method, this function allows for taxa to be ancestral to more than
one descendant taxon. Thus, users who believe their dataset may contain
indirect ancestors are encouraged by the package author to try cal3 methods
with their consensus trees, as opposed to using the set of most parsimonious
trees. Comparing the results of these two approaches may be very revealing.
Like timePaleoPhy, cal3TimePaleoPhy is designed for direct application to datasets
where taxon first and last appearances are precisely known in continuous time, with
no stratigraphic uncertainty. This is an uncommon form of data to have from the fossil record,
although not an impossible form (micropaleontologists often have very precise
range charts, for example). This means that most users should not use cal3TimePaleoPhy directly,
unless they have written their own code to deal with stratigraphic uncertainty. For
some groups, the more typical 'first' and 'last' dates represent the minimum
and maximum absolute ages for the fossil collections that a taxon is known
is known from. Presumably, the first and last appearances of that taxon in
the fossil record is at unknown dates within these bounds. These should not
be mistaken as the FADs and LADs desired by cal3TimePaleoPhy, as cal3TimePaleoPhy
will use the earliest dates provided to calibrate node ages, which is either
an overly conservative approach to time-scaling or fairly nonsensical.
If you have time-data in discrete intervals, consider using
bin_cal3TimePaleoPhy as an alternative to cal3TimePaleoPhy.
bin_cal3TimePaleoPhy is a wrapper of
cal3TimePaleoPhy which produces time-scaled trees for datasets which only have
interval data available. For each output tree, taxon first and last appearance
dates are placed within their listed intervals under a uniform distribution.
Thus, a large sample of time-scaled trees will approximate the uncertainty in
the actual timing of the FADs and LADs.
The input timeList object can have overlapping (i.e. non-sequential) intervals,
and intervals of uneven size. Taxa alive in the modern should be listed as last
occurring in a time interval that begins at time 0 and ends at time 0. If taxa
occur only in single collections (i.e. their first and last appearance in the
fossil record is synchronous, the argument point.occur will force all taxa
to have instantaneous durations in the fossil record. Otherwise, by default,
taxa are assumed to first and last appear in the fossil record at different points
in time, with some positive duration. The sites matrix can be used to force
only a portion of taxa to have simultaneous first and last appearances.
If timeData or the elements of timeList are actually data frames (as output
by read.csv or read.table), these will be coerced to a matrix.
A tutorial for applying the time-scaling functions in paleotree, particularly the cal3 method, along with an example using real (graptolite) data, can be found at the following link:
https://nemagraptus.blogspot.com/2013/06/a-tutorial-to-cal3-time-scaling-using.html
Value
The output of these functions is a time-scaled tree or set of
time-scaled trees, of either class phylo or multiPhylo, depending on the
argument ntrees. All trees are output with an element $root.time. This is
the time of the root on the tree and is important for comparing patterns
across trees.
Additional elements are sampledLogLike and $sumLogLike which respectively
record a vector containing
the 'log-densities' of the various node-ages selected for each tree by the 'zipper'
algorithm, and the sum of those log-densities. Although they are very similar to
log-likelihood values, they are not true likelihoods, as node ages are conditional on the other
ages selected by other nodes. However, these values may give an indication about the relative
optimality of a set of trees output by the cal3 functions.
Trees created with bin_cal3TimePaleoPhy will output with some additional
elements, in particular $ranges.used, a matrix which records the
continuous-time ranges generated for time-scaling each tree (essentially a
pseudo-timeData matrix.)
Note
Most importantly, please note the stochastic element of the three rate-calibrated time-scaling methods. These do not use traditional optimization methods, but instead draw divergence times from a distribution defined by the probability of intervals of unobserved evolutionary history. This means analyses MUST be done over many cal3 time-scaled trees for analytical rigor! No one tree is correct.
Similarly, please account for stratigraphic uncertainty in your analysis.
Unless you have exceptionally resolved data, use a wrapper with the cal3
function, either the provided bin_cal3TimePaleoPhy or code a wrapper
function of your own that accounts for stratigraphic uncertainty in
your dataset. Remember that the FADs (earliest dates) given to timePaleoPhy
will *always* be used to calibrate node ages!
Author(s)
David W. Bapst
References
Bapst, D. W. 2013. A stochastic rate-calibrated method for time-scaling phylogenies of fossil taxa. Methods in Ecology and Evolution. 4(8):724-733.
Foote, M. 2000. Origination and extinction components of taxonomic diversity: general problems. Pp. 74-102. In D. H. Erwin, and S. L. Wing, eds. Deep Time: Paleobiology's Perspective. The Paleontological Society, Lawrence, Kansas.
Foote, M. 2001. Inferring temporal patterns of preservation, origination, and extinction from taxonomic survivorship analysis. Paleobiology 27(4):602-630.
Friedman, M., and M. D. Brazeau. 2011. Sequences, stratigraphy and scenarios: what can we say about the fossil record of the earliest tetrapods? Proceedings of the Royal Society B: Biological Sciences 278(1704):432-439.
Stanley, S. M. 1979. Macroevolution: Patterns and Process. W. H. Freeman, Co., San Francisco.
See Also
timePaleoPhy,
make_durationFreqCont,
pqr2Ps,
sProb2sRate,
multi2di
Examples
# Simulate some fossil ranges with simFossilRecord
set.seed(444)
record <- simFossilRecord(p = 0.1,
q = 0.1,
nruns = 1,
nTotalTaxa = c(30,40),
nExtant = 0)
taxa <- fossilRecord2fossilTaxa(record)
# simulate a fossil record with imperfect sampling with sampleRanges
rangesCont <- sampleRanges(taxa,
r = 0.5)
# let's use taxa2cladogram to get the 'ideal' cladogram of the taxa
cladogram <- taxa2cladogram(taxa,
plot = TRUE)
# this package allows one to use
# rate calibrated type time-scaling methods (Bapst, 2014)
# to use these, we need an estimate of the sampling rate
# (we set it to 0.5 above)
likFun <- make_durationFreqCont(rangesCont)
srRes <- optim(
parInit(likFun),
likFun,
lower = parLower(likFun),
upper = parUpper(likFun),
method = "L-BFGS-B",
control = list(maxit = 1000000))
sRate <- srRes[[1]][2]
# we also need extinction rate and branching rate
# we can get extRate from getSampRateCont too
# we'll assume extRate = brRate (ala Foote et al., 1999)
# this may not always be a good assumption!
divRate <- srRes[[1]][1]
# now let's try cal3TimePaleoPhy
# which time-scales using a sampling rate to calibrate
# This can also resolve polytomies based on
# sampling rates, with some stochastic decisions
ttree <- cal3TimePaleoPhy(
cladogram,
rangesCont,
brRate = divRate,
extRate = divRate,
sampRate = sRate,
ntrees = 1,
plot = TRUE)
# notice the warning it gives!
phyloDiv(ttree)
# by default, cal3TimePaleoPhy may predict indirect ancestor-descendant relationships
# can turn this off by setting anc.wt = 0
ttree <- cal3TimePaleoPhy(
cladogram,
rangesCont,
brRate = divRate,
extRate = divRate,
sampRate = sRate,
ntrees = 1,
anc.wt = 0,
plot = TRUE)
# let's look at how three trees generated
# with very different time of obs. look
ttreeFAD <- cal3TimePaleoPhy(
cladogram,
rangesCont,
brRate = divRate,
extRate = divRate,
FAD.only = TRUE,
dateTreatment = "firstLast",
sampRate = sRate,
ntrees = 1,
plot = TRUE)
ttreeRand <- cal3TimePaleoPhy(
cladogram,
rangesCont,
brRate = divRate,
extRate = divRate,
FAD.only = FALSE,
dateTreatment = "randObs",
sampRate = sRate,
ntrees = 1,plot = TRUE)
# by default the time of observations are the LADs
ttreeLAD <- cal3TimePaleoPhy(
cladogram,
rangesCont,
brRate = divRate,
extRate = divRate,
FAD.only = FALSE,
dateTreatment = "randObs",
sampRate = sRate,
ntrees = 1,
plot = TRUE)
# and let's plot
layout(1:3)
parOrig <- par(no.readonly = TRUE)
par(mar = c(0,0,0,0))
plot(ladderize(ttreeFAD));text(5,5,
"time.obs = FAD",
cex = 1.5, pos = 4)
plot(ladderize(ttreeRand));text(5,5,
"time.obs = Random",
cex = 1.5, pos = 4)
plot(ladderize(ttreeLAD));text(5,5,
"time.obs = LAD",
cex = 1.5, pos = 4)
layout(1)
par(parOrig)
# to get a fair sample of trees
# let's increase ntrees
ttrees <- cal3TimePaleoPhy(
cladogram,
rangesCont,
brRate = divRate,
extRate = divRate,
sampRate = sRate,
ntrees = 9,
plot = FALSE)
# let's compare nine of them at once in a plot
layout(matrix(1:9,3,3))
parOrig <- par(no.readonly = TRUE)
par(mar = c(0,0,0,0))
for(i in 1:9){
plot(ladderize(ttrees[[i]]),
show.tip.label = FALSE)
}
layout(1)
par(parOrig)
# they are all a bit different!
# can plot the median diversity curve with multiDiv
multiDiv(ttrees)
# using node.mins
# let's say we have (molecular??) evidence that
# node (5) is at least 1200 time-units ago
# to use node.mins, first need to drop any unshared taxa
droppers <- cladogram$tip.label[is.na(
match(cladogram$tip.label,
names(which(!is.na(rangesCont[,1])))
)
)
]
# and then drop those taxa
cladoDrop <- drop.tip(cladogram, droppers)
# now make vector same length as number of nodes
nodeDates <- rep(NA, Nnode(cladoDrop))
nodeDates[5] <- 1200
ttree <- cal3TimePaleoPhy(cladoDrop,
rangesCont,
brRate = divRate,
extRate = divRate,
sampRate = sRate,
ntrees = 1,
node.mins = nodeDates,
plot = TRUE)
# example with time in discrete intervals
set.seed(444)
record <- simFossilRecord(p = 0.1,
q = 0.1,
nruns = 1,
nTotalTaxa = c(30,40),
nExtant = 0)
taxa <- fossilRecord2fossilTaxa(record)
# simulate a fossil record
# with imperfect sampling with sampleRanges
rangesCont <- sampleRanges(taxa,r = 0.5)
# let's use taxa2cladogram to get the 'ideal' cladogram of the taxa
cladogram <- taxa2cladogram(taxa,plot = TRUE)
# Now let's use binTimeData to bin in intervals of 1 time unit
rangesDisc <- binTimeData(rangesCont,int.length = 1)
# we can do something very similar for
# the discrete time data (can be a bit slow)
likFun <- make_durationFreqDisc(rangesDisc)
spRes <- optim(
parInit(likFun),
likFun,
lower = parLower(likFun),
upper = parUpper(likFun),
method = "L-BFGS-B",
control = list(maxit = 1000000))
sProb <- spRes[[1]][2]
# but that's the sampling PROBABILITY per bin
# NOT the instantaneous rate of change
# we can use sProb2sRate() to get the rate
# We'll need to also tell it the int.length
sRate1 <- sProb2sRate(sProb,int.length = 1)
# we also need extinction rate and branching rate (see above)
# need to divide by int.length...
divRate <- spRes[[1]][1]/1
# estimates that r = 0.3...
# that's kind of low (simulated sampling rate is 0.5)
# Note: for real data, you may need to use an average int.length
# (i.e. if intervals aren't all the same duration)
ttree <- bin_cal3TimePaleoPhy(cladogram,
rangesDisc,
brRate = divRate,
extRate = divRate,
sampRate = sRate1,
ntrees = 1,
plot = TRUE)
phyloDiv(ttree)
# can also force the appearance timings
# not to be chosen stochastically
ttree1 <- bin_cal3TimePaleoPhy(cladogram,
rangesDisc,
brRate = divRate,
extRate = divRate,
sampRate = sRate1,
ntrees = 1,
nonstoch.bin = TRUE,
plot = TRUE)
phyloDiv(ttree1)
# testing node.mins in bin_cal3TimePaleoPhy
ttree <- bin_cal3TimePaleoPhy(cladoDrop,
rangesDisc,
brRate = divRate,
extRate = divRate,
sampRate = sRate1,
ntrees = 1,
node.mins = nodeDates,
plot = TRUE)
# with randres = TRUE
ttree <- bin_cal3TimePaleoPhy(cladoDrop,
rangesDisc,
brRate = divRate,
extRate = divRate,
sampRate = sRate1,
ntrees = 1,
randres = TRUE,
node.mins = nodeDates,
plot = TRUE)
# example with multiple values of anc.wt
ancWt <- sample(0:1,
nrow(rangesDisc[[2]]),
replace = TRUE)
names(ancWt) <- rownames(rangesDisc[[2]])
ttree1 <- bin_cal3TimePaleoPhy(cladogram,
rangesDisc,
brRate = divRate,
extRate = divRate,
sampRate = sRate1,
ntrees = 1,
anc.wt = ancWt,
plot = TRUE)