# Three Rate Calibrated 'a posteriori' Time-Scaling of Paleo-Phylogenies

### Description

Time-scales an unscaled 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 time-scaling algorithm to stochastically resolve polytomies and infer potential ancestor-descendant relationships, simultaneous with the time-scaling treatment.

### Usage

1 2 3 4 5 6 7 8 9 10 11 12 | ```
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, tolerance = 1e-04,
diagnosticMode = FALSE, 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, 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 'bin'
time-scaling functions; however, see the argument |

`brRate` |
Either a single estimate of the instanteous 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 instanteous extinction rate (also known as the 'per-capita' extinction rate) or a vector of per-taxon estimates |

`sampRate` |
Either a single estimate of the instanteous sampling rate or a vector of per-taxon estimates |

`ntrees` |
Number of time-scaled trees to output |

`anc.wt` |
Weighting against inferring ancestor-descendant relationships. The argument anc.wt allows users to change the default consideration of anc-desc relationships. This value is used as a multiplier applied to the probability of choosing any node position which would infer an ancestor-descendant relationship. By default, anc.wt=1, and thus these probabilities are unaltered. if anc.wt is less than 1, the probabilities decrease and at anc.wt=0, no ancestor-descendant relationships are inferred at all. Can be a single value or a vector of per-taxon values, such as if a user wants to only allow plesiomorphic taxa to be ancestors. |

`node.mins` |
The minimum dates of internal nodes (clades) on a phylogeny can be set
using node.mins. This argument takes a vector of the same length as the number of nodes,
with dates given in the same order as nodes are ordered in the |

`dateTreatment` |
This argument controls the interpretation of timeData. The default setting
'firstLast' treats the dates in timeData as a column of precise first and last appearances.
A second option, added by great demand, is 'minMax' which
treats these dates as minimum and maximum bounds on single point dates. Under this option,
all taxa in the analysis will be treated as being point dates, such that the first appearance
is also the last. These dates will be pulled under a uniform distribution. If 'minMax' is used,
add.term becomes meaningless, and the use of it will return an error message. A third option
is 'randObs'. This assumes that the dates in the matrix are first and last appearance times,
but that the desired time of observation is unknown. Thus, this is much like 'firstLast' except
the effective time of observation (the taxon's LAD under 'firstLast') is treated an uncertain date, and is randomly
sampled between the first and last appearance times. The FAD still is treated as a fixed number, used
for dating the nodes. In previous versions of paleotree, this
was called in |

`FAD.only` |
Should the tips represent observation times at the start of the taxon ranges? If TRUE, result is similar to when terminal ranges are no added on with timePaleoPhy. If FAD.only is TRUE and dateTreatment is 'minMax' or 'randObs', the function will stop and a warning will be produced, as these combinations imply contradictory sets of times of observation. |

`adj.obs.wt` |
If the time of observation are before the LAD of a 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? Should only have an effect if time of observation ISN'T the LAD, if times of observation for a potential ancestor are earlier than the first appearance of their potential descendants and if the ancestral weights for taxa aren't set to zero (so there can be potential ancestors). |

`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 multi2di in ape rather than using the cal3 algorithm to weight the resolution of polytomies relative to sampling in the fossil record? |

`noisyDrop` |
If |

`tolerance` |
Acceptable amount of shift in tip dates from dates listed
in |

`diagnosticMode` |
If |

`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, simply give a matrix where the "site" of each FAD and LAD for every taxon is listed, as corresponding to the second matrix in timeList. If no sites matrix is given (the default), then it is assumed all fossil come from different "sites" and there is no shared temporal structure among the events. |

`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 true, dates are not stochastically pulled from uniform distributions. See below for more details. |

### 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, its sister method the sampling-rate-calibrated time-scaling method and the details of the algorithm will be given in Bapst (in prep for Paleobiology).

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 (ala 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, ala 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.

Alternatively to using `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.

By setting the argument `nonstoch.bin`

to `TRUE`

in `bin_cal3TimePaleoPhy`

, the
dates are NOT stochastically pulled from uniform bins but instead FADs are
assigned to the earliest time of whichever interval they were placed in and
LADs are placed at the most recent time in their placed interval. This
option may be useful for plotting. The sites argument becomes arbitrary if
nonstoch.bin is `TRUE`

.

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:

http://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 may not be 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`

, `binTimeData`

,
`make_durationFreqCont`

, `sProb2sRate`

, `multi2di`

### Examples

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 | ```
#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 library allows one to use rate calibrated type time-scaling methods (Bapst, in prep.)
#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); 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)
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])))))]
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... kind of low (simulated sampling rate is 0.5)
#Note: for real data, you may need to use an average int.length (no constant length)
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)
``` |