cal3TimePaleoPhy: Three Rate Calibrated _a posteriori_ Dating of...

In paleotree: Paleontological and Phylogenetic Analyses of Evolution

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 phylo. Tip labels must match the taxon labels in the respective temporal data.
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_ dating functions; however, see the argument dateTreatment.
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 anc.wt allows users to alter the default consideration of ancestor-descendant 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 tree$edge matrix. Note that in tree$edge, terminal tips are given the first set of numbers (1:Ntip(tree)), so the first element of node.mins is the first internal node (the node numbered Ntip(tree)+1, which is generally the root for most phylo objects read by read.tree). Not all nodes need be given minimum dates. Nodes without minimum dates can be given as NA in node.mins, but the vector must be the same length as the number of internal nodes in tree. These are minimum date constraints, such that a node will be 'frozen' by the cal3 algorithm so that constrained nodes will always be at least as old as this date, but the final date may be even older depending on the taxon dates used, the parameters input for the cal3 algorithm and any other minimum node dates given (e.g. if a clade is given a very old minimum date, this will (of course) over-ride any minimum dates given for clades that that node is nested within). if the constrained nodes include a polytomy, this polytomy will still be resolved with respect to the cal3 algorithm, but the first divergence will be 'frozen' so that it is at least as old as the minimum age, while any additional divergences will be allowed to occur after this minimum age.
dateTreatment	This argument controls the interpretation of timeData. The default setting dateTreatment = "firstLast" treats the dates in timeData as a column of precise first and last appearances. A second option is dateTreatment = "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. Note that use of dateTreatment = "minMax" was bugged in versions paleotree <3.2.5, as the generating time-tree used for cal3 inference was generated using the input timeData as fixed first and last dates, such that the effect of dateTreatment = "minMax" was nearly identical to using dateTreatment = "firstLast" regardless of the arguments chosen by a user, with tips always being placed at the upper date constraint as if the time of observation was a fixed LAD. This is now fixed as v3.2.6, such that a new basic-time-scaled tree is generated from a randomly selected set of point occurrence times on each iteration, so that the resulting tip dates are always different. A third option is dateTreatment = "randObs", which 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 dateTreatment = "firstLast" except the effective time of observation (the taxon's LAD under dateTreatment = "firstLast") is treated as 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 cal3timePaleoPhy using the argument rand.obs, which has been removed for clarity. This temporal uncertainty in times of observation might be useful if a user is interested in applying phylogeny-based approaches to studying trait evolution, but have per-taxon measurements of traits that come from museum specimens with uncertain temporal placement. With both arguments dateTreatment = "minMax" and dateTreatment = "randObs", the time of observation of taxa is a point-occurrence with a free-floating random variable as the precise age. Thus, the option FAD.only = TRUE is incoherent with these other options for dateTreatment, and thus their use together will return an error message. Furthermore, the sampling of dates from random distributions in these approaches should compel users to produce many time-scaled trees for any given analytical purpose. Note that dateTreatment = "minMax" returns an error in 'bin_' time-scaling functions; please use points.occur instead.
FAD.only	Should the tips represent observation times at the start of the taxon ranges? FAD.only = TRUE, the resulting output is similar to when terminal ranges are no added on with timePaleoPhy. If FAD.only = TRUE and dateTreatment = "minMax" or dateTreatment = "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 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 adj.obs.wt = TRUE. This argument should only have an effect if time of observation IS NOT the LAD, if the times of observation for for a potential ancestor are earlier than the first appearance of their potential descendants, and if the ancestral weights for taxa are not 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 TRUE (the default), any taxa dropped from tree due to not having a matching entry in the time data will be listed in a system message.
verboseWarnings	if TRUE (the default), then various warnings and messages regarding best practices will be issued to the console about the analysis. If FALSE,the function will run as quietly as possible.
diagnosticMode	If TRUE, cal3timePaleoPhy will return to the console and to graphic devices an enormous number of messages, plots and ancillary information that may be useful or entirely useless to figuring out what is going wrong.
tolerance	Acceptable amount of shift in tip dates from dates listed in timeData. Shifts outside of the range of tolerance will cause a warning message to be issued that terminal tips appear to be improperly aligned.
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 tip.labels for tree. See details.
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 matrix to the sites argument 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 nonstoch.bin = TRUE (the default is FALSE, dates are not stochastically drawn from uniform distributions bounded by the upper and lower boundaries of the geologic intervals (the 'bins'), as typically occurs with 'bin_' time-scaling methods in paleotree but instead first-appearance dates are assigned to the earliest time of the interval a taxon first appears in, while last-appearance dates are placed at the youngest (the 'later-most') date in the interval that that taxon last appears in. This option may be useful for plotting. Note that if nonstoch.bin = TRUE, the sites argument becomes arbitrary and has no influence on the output.

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)
    

paleotree documentation built on Aug. 22, 2022, 9:09 a.m.