Package 'paleotree'

Description

Analyzes, time-scales and simulates phylogenies of extinct/fossil lineages, along with calculation of diversity curves. Also fits likelihood models to estimate sampling rates from stratigraphic ranges.

Details

This package contains functions for analyzing sampling rates given ranges of fossil taxa, in both continuous and discrete time, functions for a posteriori time-scaling phylogenies of fossil taxa and functions for simulating the fossil record in both taxic and phylogenetic varieties.

Author(s)

David W. Bapst

Maintainer: David W. Bapst <[email protected]>

References

Bapst, D.W. 2012. paleotree: an R package for paleontological and phylogenetic analyses of evolution. Methods in Ecology and Evolution. 3: 803-807. doi: 10.1111/j.2041-210X.2012.00223.x

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.

Bapst, D. W. 2013. When Can Clades Be Potentially Resolved with Morphology? PLoS ONE. 8(4):e62312.

Bapst, D. W. 2014. Assessing the effect of time-scaling methods on phylogeny-based analyses in the fossil record. Paleobiology 40(3):331-351.

See Also

This package relies extensively on the phylogenetic toolkit and standards offered by the ape package, and hence lists this package as a depends, so it is loaded simultaneously.

Examples

# get the package version of paleotree
packageVersion("paleotree")

# get the citation for paleotree
citation("paleotree")

## 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)

# let's see what the 'true' diversity curve looks like in this case
   # plot the FADs and LADs with taxicDivCont()
taxicDivCont(taxa)

# simulate a fossil record with imperfect sampling with sampleRanges
rangesCont <- sampleRanges(taxa,r = 0.5)

# plot the diversity curve based on the sampled ranges
layout(1:2)
taxicDivCont(rangesCont)

# Now let's use binTimeData to bin in intervals of 10 time units
rangesDisc <- binTimeData(rangesCont,int.length = 10)

# plot with taxicDivDisc
taxicDivDisc(rangesDisc)

#compare to the continuous time diversity curve above!

layout(1)

# taxa2phylo assumes we know speciation events perfectly... what if we don't?

# first, let's use taxa2cladogram to get the 'ideal' cladogram of the taxa
cladogram <- taxa2cladogram(taxa,plot = TRUE)

# Now let's try timePaleoPhy using the continuous range data
ttree <- timePaleoPhy(cladogram,rangesCont,type = "basic",plot = TRUE)

# plot diversity curve
phyloDiv(ttree,drop.ZLB = TRUE)

# that tree lacked the terminal parts of ranges (tips stops at the taxon FADs)
# let's add those terminal ranges back on with add.term
ttree <- timePaleoPhy(
    cladogram,
    rangesCont,
    type = "basic",
    add.term = TRUE,
    plot = TRUE
    )

# plot diversity curve 
phyloDiv(ttree)

Description

Converts a matrix of simulated continuous-time first occurrences and last occurrences for fossil taxa into first and last occurrences given in some set of discrete-time intervals, either simulated or place a priori, which is output along with information of the dates of the given intervals.

Usage

binTimeData(timeData, int.length = 1, start = NA, int.times = NULL)

Arguments

Details

This function takes a simulated matrix of per-taxon first and last occurrences and, by dividing the time-scale into time intervals of non-zero length, lists taxon occurrences within those interval. By default, a set of sequential non-overlapping time-interval of equal non-zero length are used, with the length controlled by the argument int.length.

Alternatively, a two column matrix of interval start and end times to be used can be input via the argument int.times. None of these intervals can have a duration (temporal length) greater than zero. If a first or last appearance in the input range data could fit into multiple intervals (i.e. the input discrete time intervals are overlapping), then the appearance data is placed in the interval of the shortest duration. When output, the interval times matrix (see below) will be sorted from first to last.

As with many functions in the paleotree package, absolute time is always decreasing, i.e. the present day is zero. However, the numbering of intervals giving in the output increases with time, as these are numbered relative to each other, from first to last.

As of version 1.7 of paleotree, taxa which are extant as indicated in timeData as being in a time interval bounded (0, 0), unless time-bins are preset using argument int.times (prior to version 1.5 they were erroneously listed as NA).

Value

A list containing:

Note

This function is SPECIFICALLY for simulating the effect of having a discrete time-scale for analyses using simulations. This function should not be used for non-simulations uses, such as binning temporal occurrences for analyses of real data. In those case, the temporal ranges (which, in real data, will probably be given as discrete time intervals) should already be tabulated within discrete intervals prior to use in paleotree. The user should place the temporal information in a list object, as described for the output of binTimeData.

Author(s)

David W. Bapst

See Also

simFossilRecord, sampleRanges, taxicDivCont

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 via sampleRanges
rangesCont <- sampleRanges(taxa,r = 0.5)
# Now let's use binTimeData() to bin in intervals of 1 time unit
rangesDisc <- binTimeData(rangesCont,int.length = 1)
# plot with taxicDivDisc()
equalDiscInt <- taxicDivDisc(rangesDisc)

# example with pre-set intervals input (including overlapping)
presetIntervals <- cbind(
    c(1000, 990, 970, 940),
    c(980, 970, 950, 930)
    )
rangesDisc1 <- binTimeData(rangesCont,
    int.times = presetIntervals)

# plot the diversity curve with these uneven bins
taxicDivDisc(rangesDisc1)

# now let's plot the diversity from these unequal-length bins
   # with the original equal length intervals from above
taxicDivDisc(rangesDisc1, int.times = equalDiscInt[,1:2])


####################################
#example with extant taxa
set.seed(444)
record <- simFossilRecord(p = 0.1, 
                          q = 0.1, 
                          nruns = 1,
                          nTotalTaxa = c(30,40)
                          )
taxa <- fossilRecord2fossilTaxa(record)
# simulate a fossil record 
    # with imperfect sampling via sampleRanges
rangesCont <- sampleRanges(
    taxa, r = 0.5,
    modern.samp.prob = 1)
# Now let's use binTimeDat to bin into intervals of 1 time-unit
rangesDisc <- binTimeData(rangesCont,
    int.length = 1)
# plot with taxicDivDisc()
taxicDivDisc(rangesDisc)

# example with pre-set intervals input
    # (including overlapping)
presetIntervals <- cbind(
    c(40, 30, 20, 10),
    c(30, 20, 10, 0)
    )
rangesDisc1 <- binTimeData(rangesCont,
    int.times = presetIntervals)
    
taxicDivDisc(rangesDisc1)

Description

Partitions the branch lengths of a tree into several classes based on their placement.

Usage

branchClasses(tree, whichExtant = NULL, tol = 0.01)

Arguments

Details

This function will partition the internode (node to node, including internal node to terminal tip) branch lengths of a tree into four separate classes: all (all the internode branches of a tree), int (internal branches which run from one internode to another), live (terminal branches which run from an internal node to a terminal tip representing an extinction event before the present) and dead (terminal branches which run from an internal node to a terminal tip at the modern day, reflecting a still-living taxon).

The depths of the internal 'mother' node (i.e. time of origin, before the modern day) of each branch length are included as the labels of the branch length vectors.

This function is mainly of use for modeling internode branch lengths in a phylogeny including fossil taxa.

Value

The output is a list consisting of four vectors, with the labels of the vectors being their corresponding time of origin. See details.

Examples

#simulated example
set.seed(444)
record <- simFossilRecord(
    p = 0.1, 
    q = 0.1, 
    nruns = 1,
    nTotalTaxa = c(30,40), 
    nExtant = c(10,20)
    )
taxa <- fossilRecord2fossilTaxa(record)
tree <- taxa2phylo(taxa)
brlenRes <- branchClasses(tree)

#see frequency histograms of branch lengths
layout(1:4)
for(x in 1:length(brlenRes)){ 
	hist(
	    brlenRes[[x]],
	    main = "Branch Lengths",
	    xlab = names(brlenRes)[x])
	}

#see frequency histograms of branch depths
layout(1:4)
for(x in 1:length(brlenRes)){ 
	hist(
	    as.numeric(names(brlenRes[[x]])),
	    main = "Branch Depths",
	    xlab = names(brlenRes)[x])
	}

layout(1)

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

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)

Description

This function simulates trait evolution at each speciation/branching event in a matrix output from simFossilRecord, after transformation with fossilRecord2fossilTaxa.

Usage

cladogeneticTraitCont(taxa, rate = 1, meanChange = 0, rootTrait = 0)

Arguments

Details

This function simulates continuous trait evolution where change occurs under a Brownian model, but only at events that create new distinct morphotaxa (i.e. species as recognized in the fossil record), either branching events or anagenesis (pseudospeciation). These are the types of morphological differentiation which can be simulated in the function simFossilRecord. This is sometimes referred to as cladogenetic or speciation trait evolution and is related to Punctuated Equilibrium theory. Anagenetic shifts are not cladogenetic events per se (no branching!), so perhaps the best way to this of this function is it allows traits to change anytime simFossilRecord created a new 'morphotaxon' in a simulation.

Importantly, trait changes only occur at the base of 'new' species, thus allowing cladogenetic trait evolution to be asymmetrical at branching points: i.e. only one branch actually changes position in trait-space, as expected under a budding cladogenesis model. This distinction is important as converting the taxa matrix to a phylogeny and simulating the trait changes under a 'speciational' tree-transformation would assume that divergence occurred on both daughter lineages at each node. (This has been the standard approach for simulating cladogenetic trait change on trees).

Cryptic taxa generated with prop.cryptic in simFossilRecord will not differ at all in trait values. These species will all be identical.

See this link for additional details:

https://nemagraptus.blogspot.com/2012/03/simulating-budding-cladogenetictrait.html

Value

Returns a vector of trait values for each taxon, with value names being the taxa IDs (column 1 of the input) with a 't' pasted (as with rtree in the ape library).

Author(s)

David W. Bapst

See Also

simFossilRecord,

This function is similar to Brownian motion simulation functions such as rTraitCont in ape, sim.char in geiger and fastBM in phytools.

See also unitLengthTree in this package and speciationalTree in the package geiger. These are tree transformation functions; together with BM simulation functions, they would be expected to have a similar effect as this function (when cladogenesis is 'bifurcating' and not 'budding'; see above).

Examples

set.seed(444)
record <- simFossilRecord(
   p = 0.1, q = 0.1, 
   nruns = 1,
   nTotalTaxa = c(30, 1000), 
   plot = TRUE)
taxa <- fossilRecord2fossilTaxa(record)
trait <- cladogeneticTraitCont(taxa)
tree <- taxa2phylo(taxa)
plotTraitgram(trait, tree,
   conf.int = FALSE)

#with cryptic speciation
record <- simFossilRecord(
   p = 0.1, q = 0.1, 
   prop.cryptic = 0.5, 
   nruns = 1, 
   nTotalTaxa = c(30, 1000), 
   plot = TRUE)
taxa <- fossilRecord2fossilTaxa(record)
trait <- cladogeneticTraitCont(taxa)
tree <- taxa2phylo(taxa)
plotTraitgram(trait, tree,
   conf.int = FALSE)

Description

This is just a small collection of miscellaneous functions that may be useful, primarily for community ecology analyses, particularly for paleoecological data. They are here mainly for pedagogical reasons (i.e. for students) as they don't appear to be available in other ecology-focused packages.

Usage

pairwiseSpearmanRho(
  x,
  dropAbsent = "bothAbsent",
  asDistance = FALSE,
  diag = NULL,
  upper = NULL,
  na.rm = FALSE
)

HurlbertPIE(x, nAnalyze = Inf)

Arguments

Details

pairwiseSpearmanRho returns Spearman rho correlation coefficients based on the rank abundances of taxa (columns) within sites (rows) from the input matrix, by internally wrapping the function cor.test. It allows for various options that automatically allow for dropping taxa not shared between two sites (the default), as well as several other options. This allows the rho coefficient to behave like the Bray-Curtis distance, in that it is not affected by the number of taxa absent in both sites.

pairwiseSpearmanRho can also rescale the rho coefficients with (1-rho)/2 to provide a measure similar to a dissimilarity metric, bounded between 0 and 1. This function was written so several arguments would be in a similar format to the vegan library function vegdist. If used to obtain rho rescaled as a dissimilarity, the default output will be the lower triangle of a distance matrix object, just as is returned by default by vegdist. This behavior can be modified via the arguments for including the diagonal and upper triangle of the matrix. Otherwise, a full matrix is returned (by default) if the asDistance argument is not enabled.

HurlbertPIE provides the 'Probability of Interspecific Encounter' metric for relative community abundance data, a commonly used metric for evenness of community abundance data based on derivations in Hurlbert (1971). An optional argument allows users to apply Hurlbert's PIE to only a subselection of the most abundant taxa.

Value

pairwiseSpearmanRho will return either a full matrix (the default) or (if asDistance is true, a distance matrix, with only the lower triangle shown (by default). See details.

HurlbertPIE returns a named vector of PIE values for the input data.

Author(s)

David W. Bapst

References

Hurlbert, S. H. 1971. The nonconcept of species diversity: a critique and alternative parameters. Ecology 52(4):577-586.

See Also

twoWayEcologyCluster; example dataset: kanto

Examples

# let's load some example data:
# a classic dataset collected by Satoshi & Okido from the Kanto region

data(kanto)

rhoBothAbsent <- pairwiseSpearmanRho(kanto,dropAbsent = "bothAbsent")

#other dropping options
rhoEitherAbsent <- pairwiseSpearmanRho(kanto,dropAbsent = "eitherAbsent")
rhoNoDrop <- pairwiseSpearmanRho(kanto,dropAbsent = "noDrop")

#compare
layout(1:3)
lim <- c(-1,1)
plot(rhoBothAbsent, rhoEitherAbsent, xlim = lim, ylim = lim)
	abline(0,1)
plot(rhoBothAbsent, rhoNoDrop, xlim = lim, ylim = lim)
	abline(0,1)
plot(rhoEitherAbsent, rhoNoDrop, xlim = lim, ylim = lim)
	abline(0,1)
layout(1)

#using dropAbsent = "eitherAbsent" reduces the number of taxa so much that
	# the number of taxa present drops too low to be useful
#dropping none of the taxa restricts the rho measures to high coefficients
	# due to the many shared 0s for absent taxa

#############

# Try the rho coefficients as a rescaled dissimilarity
rhoDist <- pairwiseSpearmanRho(kanto,asDistance = TRUE,dropAbsent = "bothAbsent")

# What happens if we use these in typical distance matrix based analyses?

# Cluster analysis
clustRes <- hclust(rhoDist)
plot(clustRes)

# Principle Coordinates Analysis
pcoRes <- pcoa(rhoDist,correction = "lingoes")
scores <- pcoRes$vectors
#plot the PCO
plot(scores,type = "n")
text(labels = rownames(kanto),scores[,1],scores[,2],cex = 0.5)

##################################

# measuring evenness with Hurlbert's PIE

kantoPIE <- HurlbertPIE(kanto)

#histogram
hist(kantoPIE)
#evenness of the kanto data is fairly high

#barplot
parX <- par(mar = c(7,5,3,3))
barplot(kantoPIE,las = 3,cex.names = 0.7,
	ylab = "Hurlbert's PIE",ylim = c(0.5,1),xpd = FALSE)
par(parX)

#and we can see that the Tower has extremely low unevenness
	#...overly high abundance of ghosts?

# NOTE it doesn't matter whether we use absolute abundances
# or proportional (relative) abundances
kantoProp<-t(apply(kanto,1,function(x) x/sum(x)))
kantoPropPIE <- HurlbertPIE(kantoProp)
identical(kantoPIE,kantoPropPIE)

#let's look at evenness of 5 most abundant taxa
kantoPIE_5 <- HurlbertPIE(kanto,nAnalyze = 5)

#barplot
parX <- par(mar = c(7,5,3,3))
barplot(kantoPIE_5,las = 3,cex.names = 0.7,
	ylab = "Hurlbert's PIE for 5 most abundant taxa",ylim = c(0.5,1),xpd = FALSE)
par(parX)

Description

These functions take two trees and calculate the changes in node ages (for compareNodeAges) for shared clades or terminal branch lengths leading to shared tip taxa (for compareTermBranches).

Usage

compareNodeAges(tree1, tree2, dropUnshared = FALSE)

compareTermBranches(tree1, tree2)

Arguments

Details

For their most basic usage, these functions compare the time-scaling of two trees. Any taxa not-shared on both trees are dropped before analysis, based on tip labels.

As with many paleotree functions, calculations relating to time on trees are done with respect to any included $root.time elements. If these are not present, the latest tip is assumed to be at the present day (time = 0).

The function compareNodeAges calculates the changes in the clade ages among those clades shared by the two trees, relative to the first tree in absolute time. For example, a shift of +5 means the clade originates five time-units later in absolute time on the second tree, while a shift of -5 means the clade originated five time-units prior on the second tree.

For compareNodeAges, if tree2 is actually a multiPhylo object composed of multiple phylogenies, the output will be a matrix, with each row representing a different tree and each column a different clade shared between at least some subset of the trees in tree2 and the tree in tree1. values in the matrix are the changes in clade ages between from tree1 (as baseline) to tree2, with NA values representing a clade that is not contained in the tree represented by that row (but is contained in tree1 and at least one other tree in tree2). The matrix can be reduced to only those clades shared by all trees input via the argument dropUnshared. Note that this function distinguishes clades based on their shared taxa, and cannot so infer that two clades might be identical if it were not for single taxon within the crown of one considered clade, despite that such a difference should probably have no effect on compare a node divergence date. Users should consider their dataset for such scenarios prior to application of compareNodeAges, perhaps by dropping all taxa not included in all other trees to be considered (this is NOT done by this function).

compareTermBranches calculates the changes in the terminal branch lengths attached to tip taxa shared by the two trees, relative to the first tree. Thus, a shift of +5 means that this particular terminal taxon is connected to a terminal branch which is five time-units longer.

Value

compareTermBranches returns a vector of temporal shifts for terminal branches with the shared tip names as labels.

For the function compareNodeAges, if both tree1 and tree2 are single trees, outputs a vector of temporal shifts for nodes on tree2 with respect to tree1. If tree2 is multiple trees, then a matrix is output, with each row representing each tree in tree2 (and carrying the name of each tree, if any is given). The values are temporal shifts for each tree in tree2 with respect to tree1. For either case, the column names or element names (for a vector) are the sorted taxon names of the particular clade, the dates of which are given in that column. See above for more details. These names can be very long when large trees are considered.

See Also

dateNodes, taxa2phylo, phyloDiv

Examples

set.seed(444)
record <- simFossilRecord(p = 0.1, q = 0.1, nruns = 1,
	nTotalTaxa = c(30,40), nExtant = 0)
taxa <- fossilRecord2fossilTaxa(record)
#get the true tree
tree1 <- taxa2phylo(taxa)
#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 try timePaleoPhy using the continuous range data
tree2 <- timePaleoPhy(cladogram,rangesCont,type = "basic")
#let's look at the distribution of node shifts
hist(compareNodeAges(tree1,tree2))
#let's look at the distribution of terminal branch lengths
hist(compareTermBranches(tree1,tree2))

#testing ability to compare multiple trees with compareNodeAges
trees <- cal3TimePaleoPhy(cladogram,rangesCont,brRate = 0.1,extRate = 0.1,
    sampRate = 0.1,ntrees = 10)
nodeComparison <- compareNodeAges(tree1,trees)
#plot it as boxplots for each node
boxplot(nodeComparison,names = NULL);abline(h = 0)
#plot mean shift in node dates
abline(h = mean(apply(nodeComparison,2,mean,na.rm = TRUE)),lty = 2)

#just shifting a tree back in time
set.seed(444)
tree1 <- rtree(10)
tree2 <- tree1
tree1$root.time <- 10
compareNodeAges(tree1,tree2)
compareTermBranches(tree1,tree2)

Description

This function constrains a model to make submodels with fewer parameters, using a structure and syntax taken from the function constrain in Rich Fitzjohn's package diversitree.

Usage

constrainParPaleo(f, ..., formulae = NULL, names = parnames(f), extra = NULL)

Arguments

Details

This function is based on (but does not depend on) the function constrain from the package diversitree. Users should refer to this parent function for more detailed discussion of model constraint usage and detailed examples.

The parent function was forked to add functionality necessary for dealing with the high parameter count models typical to some paleontological analyses, particularly the inverse survivorship method. This necessitated that the new function be entirely separate from its parent. Names of functions involved (both exported and not) have been altered to avoid overlap in the package namespaces. Going forward, the paleotree package maintainer (Bapst) will try to be vigilant with respect to changes in constrain in the original package, diversitree.

Useful information from the diversitree manual (11/01/13):

"If f is a function that takes a vector x as its first argument, this function returns a new function that takes a shorter vector x with some elements constrained in some way; parameters can be fixed to particular values, constrained to be the same as other parameters, or arbitrary expressions of free parameters."

In general, formulae should be of the structure:

LHS ~ RHS

...where the LHS is the 'Parameter We Want to Constrain' and the RHS is whatever we are constraining the LHS to, usually another parameter. LHS and RHS are the 'left-hand side' and 'right-hand side' respectively (which I personally find obscure).

Like the original constrain function this function is based on, this function cannot remove constraints previously placed on a model object and there may be cases in which the constrained function may not make sense, leading to an error. The original function will sometimes issue nonsensical functions with an incorrect number/names of parameters if the parameters to be constrained are given in the wrong order in formulae.

Differences from the original constrain function from diversitree

This forked paleotree version of constrain has two additional features, both introduced to aid in constraining models with a high number of repetitive parameters. (I did not invent these models, so don't shoot the messenger.)

First, it allows nuanced control over the constraining of many parameters simultaneously, using the all and match descriptors. This system depends on parameters being named as such: name.group1.group2.group3 and so on. Each 'group' is reference to a system of groups, perhaps referring to a time interval, region, morphology, taxonomic group or some other discrete characterization among the data (almost all functions envisioned for paleotree are for per-taxon diversification models). The number of group systems is arbitrary, and may be from zero to a very large number; it depends on the 'make' function used and the arguments selected by the user. For example, the parameter x.1 would be for the parameter x in the first group of the first group system (generally a time interval for most paleotree functions). For a more complicated exampled, with the parameter x.1.3.1, the third group for the second group system (perhaps this taxonomic data point has a morphological feature not seen in some other taxa) and group 1 of the third group system (maybe biogeographic region 1? The possibilities are endless depending on user choices!).

The all option work like so: if x.all ~ x.1 is given as a formulae, then all x parameters will be constrained to equal x.1. For example, if there is x.1, x.2, x.3 and x.4 parameters for a model, x.all ~ x.1 would be equivalent to individually giving the formulae x.2~x.1, x.3~x.1 and x.4~x.1. This means that if there are many parameters of a particular type (say, 50 x parameters) it is easy to constrain all with a short expression. It is not necessary that the The all term can be used anywhere in the name of the parameter in a formulae, including to make all parameters for a given group the same. Furthermore, the LHS and RHS don't need to be same parameter group, and both can contain all statements, even multiple all statements. Consider these examples, each of which are legal, acceptable uses:

x.all ~ y.1

Constrains all values of the parameter x for every group to be equal to the single value for the parameter y for group 1 (note that there's only a single set of groups).

all.1 ~ x.1

Constrains all parameters for the first group to equal each other, here arbitrary named x.1. For example, if there is parameters named x.1, y.1 and z.1, all will be constrained to be equal to a single parameter value.

x.all.all ~ y.2.3

Constrains all values for x in any and all groups to equal the single value for y for group 2 (of system 1) and group 3 (of system 2).

x.all ~ y.all

Constrains all values of x for every group and y for every group to be equal to a single value, which by default will be reported as y.1

The match term is similar, allowing parameter values from the same group to be quickly matched and made equivalent. These match terms must have a matching (cue laughter) term both in the corresponding LHS and RHS of the formula. For example, consider x.match ~ y.match where there are six parameters: x.1, x.2, x.3, y.1, y.2 and y.3. This will effectively constrain x.1 ~ y.1, x.2 ~ y.2 and x.3 ~ y.3. This is efficient for cases where we have some parameters that we often treat as equal. For example, in paleontology, we sometimes make a simplifying assumption that birth and death rates are equal in multiple time intervals. Some additional legal examples are:

x.match.1 ~ y.match.1

This will constrain only parameters of x and y to to equal each other if they both belong to the same group for the first group system AND belong to group 1 of the first group.

all.match. ~ x.match

This will constrain all named parameters in each group to equal each other; for example, if there are parameters x.1, y.1, z.1, x.2, y.2 and z.2, this will constrain them such that y.1 ~ x.1, z.1 ~ x.1, y.2 ~ x.2 and z.2 ~ x.2, leaving x.1 and x.2 as the only parameters effectively.

There are two less fortunate qualities to the introduction of the above terminology.

Unfortunately, this package author apologizes that his programming skills are not good enough to allow more complex sets of constraints, as would be typical with the parent function, when all or match terms are included. For example, it would not be legal to attempt to constraint y.all ~ x.1 / 2, where the user presumably is attempting to constrain all y values to equal the x parameter to equal half of the x parameter for group 1. This will not be parsed as such and should return an error. However, there are workarounds, but they require using constrainParPaleo more than once. For the above example, a user could first use y.all ~ y.1 constraining all y values to be equal. Then a user could constrain with the formula y.1 ~ x.1 / 2 which would then constrain y.1 (and all the y values constrained to equal it) to be equal to the desired fraction.

Furthermore, this function expects that parameter names don't already have period-separated terms that are identical to all or match. No function in paleotree should produce such natively. If such were to occur, perhaps by specially replacing parameter names, constrainParPaleo would confuse these terms for the specialty terms described here.

Secondly, this altered version of constrain handles the parameter bounds included as attributes in functions output by the various 'make' functions. This means that if x.1 ~ y.1 is given, constrainParPaleo will test if the bounds on x.1 and y.1 are the same. If the bounds are not the same, constrainParPaleo will return an error. This is important, as some models in paleotree may make a parameter a rate (bounded zero to some value greater than one) or a probability (bounded zero to one), depending on user arguments. Users may not realize these differences and, in many cases, constraining a rate to equal a probability is nonsense (absolute poppycock!). If a user really wishes to constrain two parameters with different bounds to be equal (I have no idea why anyone would want to do this), they can use the parameter bound replacement functions described in modelMethods to set the parameter bounds as equal. Finally, once parameters with the same bounds are constrained, the output has updated bounds that reflect the new set of parameters for the new constrained function.

Value

Modified from the diversitree manual: This function returns a constrained function that can be passed through to the optimization functions of a user's choice, such as optim, find.mle in diversitree or mcmc. It will behave like any other function. However, it has a modified class attribute so that some methods will dispatch differently: parnames, for example, will return the names of the parameters of the constrained function and parInit will return the initial values for those same constrained set of parameters. All arguments in addition to x will be passed through to the original function f.

Additional useful information from the diversitree manual (11/01/13):

For help in designing constrained models, the returned function has an additional argument pars.only, when this is TRUE the function will return a named vector of arguments rather than evaluate the function (see Examples).

Author(s)

This function (and even this help file!) was originally written by Rich Fitzjohn for his library diversitree, and subsequently rewritten and modified by David Bapst.

References

FitzJohn, R. G. 2012. diversitree: comparative phylogenetic analyses of diversification in R. Methods in Ecology and Evolution 3(6):1084-1092.

See Also

As noted above, this function is strongly based on (but does not depend on) the function constrain from the library diversitree.

Examples

# simulation example with make_durationFreqCont, with three random groups
set.seed(444)
record <- simFossilRecord(
    p = 0.1, 
    q = 0.1, 
    nruns = 1,
    nTotalTaxa = c(30,40), 
    nExtant = 0
    )
taxa <- fossilRecord2fossilTaxa(record)
rangesCont <- sampleRanges(taxa,r = 0.5)
# create a groupings matrix
grp1 <- matrix(
    sample(1:3,nrow(taxa),replace = TRUE), , 1)   
likFun <- make_durationFreqCont(rangesCont, groups = grp1)

# can constrain both extinction rates to be equal
constrainFun <- constrainParPaleo(likFun, q.2 ~ q.1)

#see the change in parameter names and bounds
parnames(likFun)
parnames(constrainFun)
parbounds(likFun)
parbounds(constrainFun)

# some more ways to constrain stuff!

#constrain all extinction rates to be equal
constrainFun <- constrainParPaleo(likFun, q.all ~ q.1)
parnames(constrainFun)

#constrain all rates for everything to be a single parameter
constrainFun <- constrainParPaleo(likFun, r.all ~ q.all)
parnames(constrainFun)

#constrain all extinction rates to be equal & all sampling to be equal
constrainFun <- constrainParPaleo(likFun, q.all ~ q.1, r.all ~ r.1)
parnames(constrainFun)

#similarly, can use match.all to make all matching parameters equal each other
constrainFun <- constrainParPaleo(likFun, match.all ~ match.all)
parnames(constrainFun)

#Constrain rates in same group to be equal
constrainFun <- constrainParPaleo(likFun, r.match ~ q.match)
parnames(constrainFun)

Description

Takes a phylogeny in the form of an object of class phylo and outputs a set of topological constraints for MrBayes as a set of character strings, either printed in the R console or in a named text file, which can be used as commands in the MrBayes block of a NEXUS file for use with (you guessed it!) MrBayes.

Usage

createMrBayesConstraints(
  tree,
  partial = TRUE,
  file = NULL,
  includeIngroupConstraint = FALSE
)

Arguments

Details

partial = TRUE may be useful if the reason for using createMrBayesConstraints is to constrain a topology containing some of the taxa in an analysis, while allowing other taxa to freely vary. For example, Slater (2013) constrained an analysis so extant taxon relationships were held constant, using a molecular-based topology, while allowing fossil taxa to freely vary relative to their morphological character data.

Value

If argument file is NULL, then the constrain commands are output as a series of character strings.

Author(s)

David W. Bapst, with some inspiration from Graham Slater. This code was produced as part of a project funded by National Science Foundation grant EAR-1147537 to S. J. Carlson.

References

Slater, G. J. 2013. Phylogenetic evidence for a shift in the mode of mammalian body size evolution at the Cretaceous-Palaeogene boundary. Methods in Ecology and Evolution 4(8):734-744.

See Also

createMrBayesTipDatingNexus, createMrBayesTipCalibrations

Examples

set.seed(444)
tree <- rtree(10)
createMrBayesConstraints(tree)
createMrBayesConstraints(tree,partial = FALSE)

## Not run: 

createMrBayesConstraints(tree,file = "topoConstraints.txt")


## End(Not run)

Description

Takes a set of tip ages (in several possible forms, see below), and outputs a set of tip age calibrations for use with tip-dating analyses (sensu Zhang et al., 2016) in the popular phylogenetics program MrBayes. These calibrations are printed as a set of character strings, as well as a line placing an offset exponential prior on the tree age, either printed in the R console or in a named text file, which can be used as commands in the MrBayes block of a NEXUS file for use with (you guessed it!) MrBayes.

Usage

createMrBayesTipCalibrations(
  tipTimes,
  ageCalibrationType,
  whichAppearance = "first",
  treeAgeOffset,
  minTreeAge = NULL,
  collapseUniform = TRUE,
  anchorTaxon = TRUE,
  file = NULL
)

Arguments

Details

Beware: some combinations of arguments might not make sense for your data.

(But that's always true, is it not?)

Value

If argument file is NULL, then the tip age commands are output as a series of character strings.

All taxa with their ages set to fixed by the behavior of anchorTaxon or collapseUniform are returned as a list within a commented line of the returned MrBayes block.

Author(s)

David W. Bapst. This code was produced as part of a project funded by National Science Foundation grant EAR-1147537 to S. J. Carlson.

References

Zhang, C., T. Stadler, S. Klopfstein, T. A. Heath, and F. Ronquist. 2016. Total-Evidence Dating under the Fossilized Birth-Death Process. Systematic Biology 65(2):228-249.

See Also

createMrBayesConstraints, createMrBayesTipDatingNexus

Examples

# load retiolitid dataset
data(retiolitinae)

# uniform prior, with a 10 million year offset for
	# the expected tree age from the earliest first appearance

createMrBayesTipCalibrations(
   tipTimes = retioRanges, 
   whichAppearance = "first",
	  ageCalibrationType = "uniformRange", 
	  treeAgeOffset = 10)

# fixed prior, at the earliest bound for the first appearance

createMrBayesTipCalibrations(
   tipTimes = retioRanges, 
   whichAppearance = "first",
	  ageCalibrationType = "fixedDateEarlier", 
	  treeAgeOffset = 10
	  )

# fixed prior, sampled from between the bounds on the last appearance
	# you should probably never do this, fyi

createMrBayesTipCalibrations(
   tipTimes = retioRanges, 
   whichAppearance = "first",
	  ageCalibrationType = "fixedDateRandom", 
	  treeAgeOffset = 10
	  )


## Not run: 

createMrBayesTipCalibrations(
   tipTimes = retioRanges, 
   whichAppearance = "first",
	  ageCalibrationType = "uniformRange", 
	  treeAgeOffset = 10, 
	  file = "tipCalibrations.txt"
	  )


## End(Not run)

Description

This function is meant to expedite the creation of NEXUS files formatted for performing tip-dating analyses in the popular phylogenetics software MrBayes, particularly clock-less tip-dating analyses executed with 'empty' morphological matrices (i.e. where all taxa are coded for a single missing character), although a pre-existing morphological matrix can also be input by the user (see argument origNexusFile). Under some options, this pre-existing matrix may be edited by this function. The resulting full NEXUS script is output as a set of character strings either printed to the R console, or output to file which is then overwritten.

Usage

createMrBayesTipDatingNexus(
  tipTimes,
  outgroupTaxa = NULL,
  treeConstraints = NULL,
  ageCalibrationType,
  whichAppearance = "first",
  treeAgeOffset,
  minTreeAge = NULL,
  collapseUniform = TRUE,
  anchorTaxon = TRUE,
  newFile = NULL,
  origNexusFile = NULL,
  parseOriginalNexus = TRUE,
  createEmptyMorphMat = TRUE,
  orderedChars = NULL,
  morphModel = "strong",
  morphFiltered = "parsInf",
  runName = NULL,
  ngen = "100000000",
  doNotRun = FALSE,
  autoCloseMrB = FALSE,
  cleanNames = TRUE,
  printExecute = TRUE
)

Arguments

Details

Users must supply a data set of tip ages (in various formats), which are used to construct age calibrations commands on the tip taxa (via paleotree function createMrBayesTipCalibrations). The user must also supply some topological constraint: either a set of taxa designated as the outgroup, which is then converted into a command constraining the monophyly on the ingroup taxa, which is presumed to be all taxa not listed in the outgroup. Alternatively, a user may supply a tree which is then converted into a series of hard topological constraints (via function createMrBayesConstraints. Both types of topological constraints cannot be applied.

Many of the options available with createMrBayesTipCalibrations are available with this function, allowing users to choose between fixed calibrations or uniform priors that approximate stratigraphic uncertainty. In addition, the user may also supply a path to a text file presumed to be a NEXUS file containing character data formatted for use with MrBayes.

The taxa listed in tipTimes must match the taxa in treeConstraints, if such is supplied. If supplied, the taxa in outgroupTaxa must be contained within this same set of taxa. These all must have matches in the set of taxa in origNexusFile, if provided and if parseOriginalNexus is TRUE.

Note that because the same set of taxa must be contained in all inputs, relationships are constrained as 'hard' constraints, rather than 'partial' constraints, which allows some taxa to float across a partially fixed topology. See the documentation for createMrBayesConstraints, for more details.

Value

If argument newFile is NULL, then the text of the generated NEXUS script is output to the console as a series of character strings.

Note

This function allows a user to take an undated phylogenetic tree in R, and a set of age estimates for the taxa on that tree, and produce a posterior sample of dated trees using the MCMCMC in MrBayes, while treating an 'empty' morphological matrix as an uninformative set of missing characters. This 'clock-less tip-dating' approach is essentially an alternative to the cal3 method in paleotree, sharing the same fundamental theoretical model (a version of the fossilized birth-death model), but with a better algorithm that considers the whole tree simultaneously, rather than evaluating each node individually, from the root up to the tips (as cal3 does it, and which may cause artifacts). That said, cal3 still has a few advantages: tip-dating as of April 2017 still only treats OTUs as point observations, contained in a single time-point, while cal3 can consider taxa as having durations with first and last occurrences. This means it may be more straightforward to assess the extent of budding cladogenesis patterns of ancestor-descendant relationships in cal3, than in tip-dating.

Author(s)

David W. Bapst. This code was produced as part of a project funded by National Science Foundation grant EAR-1147537 to S. J. Carlson.

The basic MrBayes commands utilized in the output script are a collection of best practices taken from studying NEXUS files supplied by April Wright, William Gearty, Graham Slater, Davey Wright, and guided by the recommendations of Matzke and Wright, 2016 in Biology Letters.

References

The basic fundamentals of tip-dating, and tip-dating with the fossilized birth-death model are introduced in these two papers:

Ronquist, F., S. Klopfstein, L. Vilhelmsen, S. Schulmeister, D. L. Murray, and A. P. Rasnitsyn. 2012. A Total-Evidence Approach to Dating with Fossils, Applied to the Early Radiation of the Hymenoptera. Systematic Biology 61(6):973-999.

Zhang, C., T. Stadler, S. Klopfstein, T. A. Heath, and F. Ronquist. 2016. Total-Evidence Dating under the Fossilized Birth-Death Process. Systematic Biology 65(2):228-249.

For recommended best practices in tip-dating analyses, please see:

Matzke, N. J., and A. Wright. 2016. Inferring node dates from tip dates in fossil Canidae: the importance of tree priors. Biology Letters 12(8).

The rationale behind the two alternative morphological models are described in more detail here:

Bapst, D. W., H. A. Schreiber, and S. J. Carlson. 2018. Combined Analysis of Extant Rhynchonellida (Brachiopoda) using Morphological and Molecular Data. Systematic Biology 67(1):32-48.

See Also

This function wraps various aspects of the functions createMrBayesConstraints and the function createMrBayesTipCalibrations. In many ways, this functionality is a replacement for the probabilistic dating method represented by the cal3 dating functions.

For putting the posterior estimated trees on an absolute time scale, see functions obtainDatedPosteriorTreesMrB. Use the argument getFixedTimes = TRUE if you used a taxon with a fixed age, and function setRootAges to set the root age.

Examples

# load retiolitid dataset
data(retiolitinae)

# let's try making a NEXUS file!

# Use a uniform prior, with a 10 million year offset for
	 # the expected tree age from the earliest first appearance

# Also set average tree age to be 10 Ma earlier than first FAD

outgroupRetio <- "Rotaretiolites" 
# this taxon will now be sister to all other included taxa

# the following will create a NEXUS file 
  # with an 'empty' morph matrix
	 # where the only topological constraint is on ingroup monophyly
	 # Probably shouldn't do this: leaves too much to the FBD prior
 
# with doNotRun set to TRUE for troubleshooting

createMrBayesTipDatingNexus(
tipTimes = retioRanges,
		outgroupTaxa = outgroupRetio,
		treeConstraints = NULL,
		ageCalibrationType = "uniformRange",
		whichAppearance = "first",
		treeAgeOffset = 10,	
		newFile = NULL,	
		origNexusFile = NULL,
		createEmptyMorphMat = TRUE,
		runName = "retio_dating",
		doNotRun = TRUE
		)

# let's try it with a tree for topological constraints
     # this requires setting outgroupTaxa to NULL
# let's also set doNotRun to FALSE

createMrBayesTipDatingNexus(
   tipTimes = retioRanges,
		outgroupTaxa = NULL,
		treeConstraints = retioTree,
		ageCalibrationType = "uniformRange",
		whichAppearance = "first",
		treeAgeOffset = 10,	
		newFile = NULL,	
		origNexusFile = NULL,
		createEmptyMorphMat = TRUE,
		runName = "retio_dating",
		doNotRun = FALSE
		)

# the above is essentially cal3 with a better algorithm,
		# and no need for a priori rate estimates
# just need a tree and age estimates for the tips!

####################################################
# some more variations for testing purposes

# no morph matrix supplied or generated
	# you'll need to manually append to an existing NEXUS file
	
createMrBayesTipDatingNexus(
   tipTimes = retioRanges,
		outgroupTaxa = NULL,
		treeConstraints = retioTree,
		ageCalibrationType = "uniformRange",
		whichAppearance = "first",
		treeAgeOffset = 10,
		newFile = NULL,	
		origNexusFile = NULL,
		createEmptyMorphMat = FALSE,
		runName = "retio_dating",
		doNotRun = TRUE
		)

## Not run: 

# let's actually try writing an example with topological constraints
	# to file and see what happens

# here's my super secret MrBayes directory
file <- "D:\\dave\\workspace\\mrbayes\\exampleRetio.nex"

createMrBayesTipDatingNexus(
   tipTimes = retioRanges,
		outgroupTaxa = NULL,
		treeConstraints = retioTree,
		ageCalibrationType = "uniformRange",
		whichAppearance = "first",
		treeAgeOffset = 10,	
		newFile = file,	
		origNexusFile = NULL,
		createEmptyMorphMat = TRUE,
		runName = "retio_dating",
		doNotRun = FALSE
		)


## End(Not run)

Description

This function returns the ages of nodes (both internal and terminal tips) for a given phylogeny of class phylo. Its use is specialized for application to dated trees from paleotree, see Details below.

Usage

dateNodes(
  tree,
  rootAge = tree$root.time,
  labelDates = FALSE,
  tolerance = 0.001
)

Arguments

Details

This function is specialized for dated phylogenies, either estimated from empirical data or simulated with functions from paleotree, and thus have a $root.time element. This function will still work without such, but users should see the details for the rootAge argument.

Value

Returns a vector of length Ntip(tree) + Nnode(tree) which contains the dates for all terminal tip nodes and internal nodes for the tree, in that order, as numbered in the tree$edge matrix. These dates are always on a descending scale (i.e. time before present); see help for argument rootAge for how the present time is determined. If rootAge is so defined that some nodes may occur later than time = 0 units before present, this function may (confusingly) return negative dates and a warning message will be issued.

Author(s)

David W. Bapst, based on a function originally written by Graeme Lloyd.

See Also

compareTimescaling, nodeDates2branchLengths

Examples

#let's simulate some example data
set.seed(444)
record <- simFossilRecord(p = 0.1, q = 0.1, nruns = 1,
	nTotalTaxa = c(30,40), nExtant = 0)
taxa <- fossilRecord2fossilTaxa(record)
#get the true time-sclaed tree
tree1 <- taxa2phylo(taxa)

#now let's try dateNodes
dateNodes(tree1)

#let's ignore $root.time
dateNodes(tree1,rootAge = NULL)

#with the lengthy tip-label based labels
   #some of these will be hideously long
dateNodes(tree1,labelDates = TRUE)

Description

The function dateTaxonTreePBDB takes a input consisting of a topology, with tip and internal node labels corresponding to taxa in the Paleobiology Database, and a table of data (containing those same tip and node taxa) obtained from the taxa-list functionality of the Paleobiology Database's API, with appearance times. This function will then output a tree with nodes reflecting the ages of the respective higher taxa, based on their earliest times of appearance from the Paleobiology Database.

Usage

dateTaxonTreePBDB(
  taxaTree,
  taxaDataPBDB = taxaTree$taxaDataPBDB,
  minBranchLen = 0,
  tipTaxonDateUsed = "shallowestAge",
  dropZeroOccurrenceTaxa = TRUE,
  plotTree = FALSE,
  failIfNoInternet = TRUE
)

Arguments

Details

The dating by this function is very simplistic, representing a rather straight interpretation of what the PBDB reports. The dated trees produced should not be taken overly seriously.

Value

Returns a dated phylogeny of class phylo, with an additional element $taxaDataPBDB added containing the input taxaDataPBDB, as this might be called by other functions.

Author(s)

David W. Bapst

References

Peters, S. E., and M. McClennen. 2015. The Paleobiology Database application programming interface. Paleobiology 42(1):1-7.

The equuid tree used in the examples is from: MacFadden, B. J. 1992. Fossil horses: systematics, paleobiology, and evolution of the family Equidae. Cambridge University Press.

See Also

See getDataPBDB, makePBDBtaxonTree, and plotPhyloPicTree.

Examples

# Note that all examples here use argument 
    # failIfNoInternet = FALSE so that functions do
    # not error out but simply return NULL if internet
    # connection is not available, and thus
    # fail gracefully rather than error out (required by CRAN).
# Remove this argument or set to TRUE so functions fail
    # when internet resources (paleobiodb) is not available.
    


taxaAnimals <- c("Archaeopteryx", "Eldredgeops",
	"Corvus", "Acropora", "Velociraptor", "Gorilla", 
	"Olenellus", "Lingula", "Dunkleosteus",
	"Tyrannosaurus", "Triceratops", "Giraffa",
	"Megatheriidae", "Aedes", "Histiodella",
	"Rhynchotrema", "Pecten", "Homo", "Dimetrodon",
	"Nemagraptus", "Panthera", "Anomalocaris")

animalData <-getSpecificTaxaPBDB(taxaAnimals, 
    failIfNoInternet = FALSE)
    
if(!is.null(animalData)){

tree <- makePBDBtaxonTree(animalData, 
    rankTaxon = "genus")

#get the ranges 
timeTree <- dateTaxonTreePBDB(tree)
    
}



#####################################

## Not run: 
 
# plotting the tree with phyloPics

plotPhyloPicTree(tree = timeTree,
     depthAxisPhylo = TRUE, 
     failIfNoInternet = FALSE)


# can also plot dated tree with strap

library(strap)
#now plot it
strap::geoscalePhylo(
    tree = timeTree,
    direction = "upwards",
    ages = rangesMinMax,
    cex.tip = 0.7,
    cex.ts = 0.55,
    cex.age = 0.5,
    width = 3,
    tick.scale = 50,
    quat.rm = TRUE,
    boxes = "Period",
    arotate = 90,
    units = c("Eon","Period","Era"),
    x.lim = c(650,-20)
    )

## End(Not run)

##############################################################

## HORSES


#if(require(curl)){

# we can also use this for pre-existing trees
    # for example, this tree of equuids (horses)
    # borrowed from UCMP materials on horse evolution
    # https://evolution.berkeley.edu/evolibrary/images/HorseTree.pdf
    # (apparently from MacFadden, 1992? Citation above)

# read the tree in as Newick string
horseTree <- ape::read.tree(file=NULL, 
    text = paste0(
         "(Eohippus,(Xenicohippus,(Haplohippus,(Epihippus,",
         "(Miohippus,(((Hypohippus,Megahippus),(Anchitherium,",
         "Kalobatippus)),(Archaeohippus,(Desmatippus,(Parahippus,",
         "(Merychippus,(((Hipparion_Merychippus,(Nannippus,",
         " Cormohipparion)),(Pseudhipparion,(Neohipparion,",
         " Hipparion))),(Equine_Merychippus,((Protohippus,Calippus),",
         "(Pliohippus,(Astrohippus,(Dinohippus,Equus))))))))))))))));"
         )
    )

# note there is a message that the tree lacks node names
    # this is unexpected / atypical for taxon trees

plot(horseTree)

# now let's get data on the tip from the PBDB
    # using getSpecificTaxaPBDB
horseData <- getSpecificTaxaPBDB(horseTree$tip.label, 
    failIfNoInternet = FALSE)

if(!is.null(horseData)){

# now we can date the tree with dateTaxonTreePBDB

datedHorseTree <- dateTaxonTreePBDB(
    taxaTree = horseTree,
    taxaDataPBDB = horseData,
    minBranchLen = 1, 
    failIfNoInternet = FALSE)

# and let's try plotting it!	
plotPhyloPicTree(
    tree = datedHorseTree,
    depthAxisPhylo = TRUE, 
    failIfNoInternet = FALSE)		
	
# a fairly boring phylopic diagram
     # not many horse phylopics as of 07-16-19?

}

#}


## Not run: 

# Let's look at this horse tree with strap

library(strap)

geoscalePhylo(
    tree = datedHorseTree,
    ages = datedHorseTree$ranges.used,
    cex.tip = 0.7,
    cex.ts = 0.7,
    cex.age = 0.7,
    width = 4,
    tick.scale = 15,
    boxes = "Epoch",
    erotate = 90,
    quat.rm = TRUE,
    units = c("Period","Epoch"),
    x.lim = c(65,-10)
    )


## End(Not run)

Description

degradeTree removes a proportion of the total nodes in a tree, chosen randomly, collapsing the nodes to produce a less-resolved tree. The related function collapseNodes given a tree and a vector of nodes to collapse, removes those nodes from a tree, creating a polytomy.

Usage

degradeTree(
  tree,
  prop_collapse = NULL,
  nCollapse = NULL,
  node.depth = NA,
  leave.zlb = FALSE
)

collapseNodes(tree, nodeID, collapseType, leave.zlb = FALSE)

Arguments

Details

In the function degradeTree, the nodes are removed at random using the basic R function sample. degradeTree can be conditioned to remove nodes of a particular depth with greater probability/frequency by setting node.depth to a value between zero (favoring the removal of deep nodes close to the root) or one (shallow nodes far from the root). Depth is evaluated based on the number of descendant tips. If node.depth is not NA, the relative proportion of descendants from each node is calculated, summed to 1 and the node.depth value subtracted from this proportion. These values are then squared, normalized again to equal to 1 and then used as the probabilities for sampling nodes for removal.

By default, branch lengths are removed from the input tree prior to degradation and entirely absent from the output tree. This is changed if argument leave.zlb is TRUE.

Value

Returns the modified tree as an object of class phylo, with no edge lengths by default.

Author(s)

David W. Bapst

See Also

di2multi,timeLadderTree

Examples

set.seed(444)
tree <- rtree(100)
tree1 <- degradeTree(tree,prop_collapse = 0.5) 
tree3 <- degradeTree(tree,nCollapse = 50) 

#let's compare the input and output
layout(matrix(1:2,,2))
plot(tree,show.tip.label = FALSE,use.edge.length = FALSE)
plot(tree1,show.tip.label = FALSE,use.edge.length = FALSE)

#now with collapseNodes
tree <- rtree(10)
#collapse nodes backwards
   #let's collapse lucky node number 13!
tree1 <- collapseNodes(nodeID = 13,tree = tree,collapseType = "backward")  
#collapse nodes forwards 
tree2 <- collapseNodes(nodeID = 13,tree = tree,collapseType = "forward")
#collapse entire clade
tree3 <- collapseNodes(nodeID = 13,tree = tree,collapseType = "clade")

#let's compare
layout(1:4)
plot(tree,use.edge.length = FALSE,main = "original")
plot(tree1,use.edge.length = FALSE,main = "backward collapse")
plot(tree2,use.edge.length = FALSE,main = "forward collapse")
plot(tree3,use.edge.length = FALSE,main = "entire clade")

layout(1)

Description

Paints the edges of a phylogeny with colors relative to their depth.

Usage

depthRainbow(tree)

Arguments

Details

The only purpose of this function is to make an aesthetically-pleasing graphic of one's tree, where branches are color-coded with a rainbow palette, relative to their depth. Depth is defined relative to the number of branching nodes between the basal node of a branch and the root, not the absolute distance (i.e. branch length) to the root or the distance from the tips.

Value

No value returned, just plots a colorful phylogeny.

Examples

set.seed(444)
tree <- rtree(500)
depthRainbow(tree)

Description

An extremely simple plotting function, which plots the original taxonomic diversity versus the sampled taxonomic diversity, for use with output from the function simFossilRecord. If sampling processes were not included in the model, then it plots simply the single diversity curve.

Usage

divCurveFossilRecordSim(
  fossilRecord,
  merge.cryptic = TRUE,
  plotLegend = TRUE,
  legendPosition = "topleft",
  curveColors = c("black", "red"),
  curveLineTypes = c(1, 2)
)

Arguments

Details

This function is essentially a wrapper for paleotree function multiDiv.

Value

This function returns nothing: it just creates a plot.

Author(s)

David W. Bapst

See Also

simFossilRecord

Examples

set.seed(44)
record <- simFossilRecord(p = 0.1, q = 0.1, r = 0.1, nruns = 1,
	nTotalTaxa = c(20,30) ,nExtant = 0, plot = FALSE)

# now let's plot it
divCurveFossilRecordSim(record)

Description

Functions to plot diversity curves based on taxic range data, in both discrete and continuous time, and for phylogenies.

Usage

taxicDivCont(
  timeData,
  int.length = 1,
  int.times = NULL,
  plot = TRUE,
  plotLogRich = FALSE,
  timelims = NULL,
  drop.cryptic = FALSE
)

taxicDivDisc(
  timeList,
  int.times = NULL,
  drop.singletons = FALSE,
  plot = TRUE,
  plotLogRich = FALSE,
  timelims = NULL,
  extant.adjust = 0.001,
  split.int = TRUE
)

phyloDiv(
  tree,
  int.length = 0.1,
  int.times = NULL,
  plot = TRUE,
  plotLogRich = FALSE,
  drop.ZLB = TRUE,
  timelims = NULL
)

Arguments

Details

First, some background. Diversity curves are plots of species/taxon/lineage richness over time for a particular group of organisms. For paleontological studies, these are generally based on per-taxon range data while more recently in evolutionary biology, molecular phylogenies have been used to calculate lineage-through-time plots (LTTs). Neither of these approaches are without their particular weaknesses; reconstructing the true history of biodiversity is a difficult task no matter what data is available.

The diversity curves produced by these functions will always measure diversity within binned time intervals (and plot them as rectangular bins). For continuous-time data or phylogenies, one could decrease the int.length used to get what is essentially an 'instantaneous' estimate of diversity. This is warned against, however, as most historical diversity data will have some time-averaging or uncertain temporal resolution and thus is probably not finely-resolved enough to calculate instantaneous estimates of diversity.

As with many functions in the paleotree library, absolute time is always decreasing, i.e. the present day is zero.

As diversity is counted within binned intervals, the true standing diversity may be somewhat lower than the measured (observed) quantity, particularly if intervals are longer than the mean duration of taxa is used. This will be an issue with all diversity curve functions, but particularly the discrete-time variant. For diversity data in particularly large discrete time intervals, plotting this data in smaller bins which do not line up completely with the original intervals will create a 'spiky' diversity curve, as these smaller intersecting bins will have a large number of taxa which may have been present in either of the neighboring intervals. This will give these small bins an apparently high estimated standing diversity. This artifact is avoided with the default setting split.int = TRUE, which will split any input or calculated intervals so that they start and end at the boundaries of the discrete-time range bins.

The timeList object should be a list composed of two matrices, the first matrix giving by-interval start and end times (in absolute time), the second matrix giving the by-taxon first and last appearances in the intervals defined in the first matrix, numbered as the rows. Absolute time should be decreasing, while the intervals should be numbered so that the number increases with time. 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. See the documentation for the time-scaling function bin_timePaleoPhy and the simulation function binTimeData for more information on formatting.

Unlike some paleotree functions, such as perCapitaRates, the intervals can be overlapping or of unequal length. The diversity curve functions deal with such issues by assuming taxa occur from the base of the interval they are first found in until the end of the last interval they are occur in. Taxa in wide-ranging intervals that contain many others will be treated as occurring in all nested intervals.

phyloDiv will resolve polytomies to be dichotomous nodes separated by zero-length branches prior to calculating the diversity curve. There is no option to alter this behavior, but it should not affect the use of the function because the addition of the zero-length branches should produce an identical diversity history as a polytomy. phyloDiv will also drop zero-length terminal branches, as with the function dropZLB. This the default behavior for the function but can be turned off by setting the argument drop.zlb to FALSE.

Value

These functions will invisibly return a three-column matrix, where the first two columns are interval start and end times and the third column is the number of taxa (or lineages) counted in that interval.

Author(s)

David W. Bapst

See Also

multiDiv, timeSliceTree, binTimeData

There are several different functions for traditional LTT plots (phylogenetic diversity curves), such as the function ,ltt.plot in the package ape, the function ltt in the package phytools, the function plotLtt in the package laser and the function LTT.average.root in the package TreeSim.

Examples

# taxicDivDisc example with the retiolinae dataset
data(retiolitinae)
taxicDivDisc(retioRanges)

##################################################

# simulation examples

# 07-15-19
# note that the examples below are weird and rather old
    # the incomplete sampling can now be done
    # with the same function that simulates diversification

set.seed(444)

record <- simFossilRecord(
       p = 0.1, 
       q = 0.1, 
       nruns = 1,
       nTotalTaxa = c(30,40), 
       nExtant = 0)
taxa <- fossilRecord2fossilTaxa(record)

# let's see what the 'true' diversity curve looks like in this case
#plot the FADs and LADs with taxicDivCont
taxicDivCont(taxa)

# simulate a fossil record with imperfect sampling via sampleRanges
rangesCont <- sampleRanges(taxa, r = 0.5)

# plot the diversity curve based on the sampled ranges
layout(1:2)
taxicDivCont(rangesCont)
# Now let's use binTimeData to bin in intervals of 1 time unit
rangesDisc <- binTimeData(rangesCont,
       int.length = 1)
# plot with taxicDivDisc
taxicDivDisc(rangesDisc)
# compare to the continuous time diversity curve

layout(1)
# Now let's make a tree using taxa2phylo
tree <- taxa2phylo(taxa,obs_time = rangesCont[,2])
phyloDiv(tree)

# a simple example with phyloDiv
  # using a tree from rtree in ape
set.seed(444)
tree <- rtree(100)
phyloDiv(tree)

###########################################################

#a neat example of using phyDiv with timeSliceTree 
 #to simulate doing molecular-phylogeny studies 
 #of diversification...in the past

set.seed(444)
record <- simFossilRecord(
       p = 0.1, 
       q = 0.1, 
       nruns = 1,
       nTotalTaxa = c(30,40), 
       nExtant = 0)
taxa <- fossilRecord2fossilTaxa(record)
taxicDivCont(taxa)

#that's the whole diversity curve

#with timeSliceTree we could look at the lineage accumulation curve 
 #we'd get of species sampled at a point in time

tree <- taxa2phylo(taxa)
#use timeSliceTree to make tree of relationships up until time = 950 
tree950 <- timeSliceTree(tree,
       sliceTime = 950,
       plot = TRUE,
       drop.extinct = FALSE)

#use drop.extinct = TRUE to only get the tree of lineages extant at time = 950
tree950 <- timeSliceTree(tree,
       sliceTime = 950,
       plot = TRUE,
       drop.extinct = TRUE)

#now its an ultrametric tree with many fewer tips...
#lets plot the lineage accumulation plot on a log scale
phyloDiv(tree950, 
       plotLogRich = TRUE)

##################################################
#an example of a 'spiky' diversity curve 
       # and why split.int is a good thing
set.seed(444)
record <- simFossilRecord(
       p = 0.1, 
       q = 0.1, 
       nruns = 1,
       nTotalTaxa = c(30,40), 
       nExtant = 0)
taxa <- fossilRecord2fossilTaxa(record)

taxaDiv <- taxicDivCont(taxa)

#simulate a fossil record with imperfect sampling with sampleRanges
rangesCont <- sampleRanges(taxa, r = 0.5)
rangesDisc <- binTimeData(rangesCont,
       int.length = 10)

#now let's plot with taxicDivDisc
 # but with the intervals from taxaDiv
 # by default, split.int = TRUE

taxicDivDisc(rangesDisc,
       int.times = taxaDiv[,1:2],
       split.int = TRUE)

#look pretty!

#now let's turn off split.int
taxicDivDisc(rangesDisc,
       int.times = taxaDiv[,1:2],
       split.int = FALSE)
#looks 'spiky'!

Description

These functions construct likelihood models of the observed frequency of taxon durations, given either in discrete (make_durationFreqDisc) or continuous time (make_durationFreqCont). These models can then be constrained using functions available in this package and/or analyzed with commonly used optimizing functions.

Usage

make_durationFreqCont(
  timeData,
  groups = NULL,
  drop.extant = TRUE,
  threshold = 0.01,
  tol = 1e-04
)

make_durationFreqDisc(timeList, groups = NULL, drop.extant = TRUE)

Arguments

Details

These functions effectively replace two older functions in paleotree, now removed, getSampRateCont and getSampProbDisc. The functions here do not offer the floating time interval options of their older siblings, but do allow for greater flexibility in defining constrains on parameter values. Differences in time intervals, or any other conceivable discrete differences in parameters, can be modeled using the generic groups argument in these functions.

These functions use likelihood functions presented by Foote (1997). These analyses are ideally applied to data from single stratigraphic section but potentially are applicable to regional or global datasets, although the behavior of those datasets is less well understood.

As with many functions in the paleotree library, absolute time is always decreasing, i.e. the present day is zero and older dates are 'larger'. On the contrary, relative time is in intervals with non-zero integers that increase sequentially beginning with 1, from earliest to oldest.

For make_durationFreqDisc, the intervals in timeList should be non-overlapping sequential intervals of roughly equal length. These should be in relative time as described above, so the earliest interval should be listed as 1 and the numbering should increase as the intervals go up with age. If both previous statements are TRUE, then differences in interval numbers will represent the same rough difference in the absolute timing of those intervals. For example, a dataset where all taxa are listed from a set of sequential intervals of similar length, such as North American Mammal assemblage zones, microfossil faunal zones or graptolite biozones can be given as long as they are correctly numbered in sequential order in the input. As a counter example, a dataset which includes taxa resolved only to intervals as wide as the whole Jurassic and taxa resolved to biozones within the Jurassic should not be included in the same input. Drop taxa from less poorly resolved intervals from such datasets if you want to apply this function, as long as this retains a large enough sample of taxa listed from the sequential set of intervals.

Please check that the optimizer function you select actually converges. The likelihood surface can be very flat in some cases, particularly for small datasets (<100 taxa). If the optimizer does not converge, consider increasing iterations or changing the starting values.

Value

A function of class "paleotreeFunc", which takes a vector equal to the number of parameters and returns the *negative* log-likelihood (for use with optim and similar optimizing functions, which attempt to minimize support values). See the functions listed at modelMethods for manipulating and examining such functions and constrainParPaleo for constraining parameters.

Parameters in the output functions are named q for the instantaneous per-capita extinction rate, r for the instantaneous per-capita sampling rate and R for the per-interval taxonomic sampling probability. Groupings follow the parameter names, separated by periods; by default, the parameters will be placed in at least group '1' (of a grouping scheme containing only a single group), such that make_durationFreqCont by default creates a function with parameters named q.1 and r.1, while make_durationFreqDisc creates a function with parameters named q.1 and R.1.

Note that the q parameters estimated by make_durationFreqDisc is scaled to per lineage intervals and not to per lineage time-units. If intervals are the same length, this can be easily corrected by multiplying one by the interval length. It is unclear how to treat uneven intervals and I urge users to consider multiple strategies.

For translating these sampling probabilities and sampling rates, see SamplingConv.

Author(s)

David W. Bapst

References

Foote, M. 1997 Estimating Taxonomic Durations and Preservation Probability. Paleobiology 23(3):278–300.

Foote, M., and D. M. Raup. 1996 Fossil preservation and the stratigraphic ranges of taxa. Paleobiology 22(2):121–140.

See Also

See freqRat, sRate2sProb, qsRate2Comp sProb2sRate and qsProb2Comp.

Examples

# let's simulate some taxon ranges from 
   # an imperfectly sampled fossil record
set.seed(444)
record <- simFossilRecord(
    p = 0.1, 
    q = 0.1, 
    nruns = 1,
	   nTotalTaxa = c(30,40), 
	   nExtant = 0
	   )
taxa <- fossilRecord2fossilTaxa(record)
rangesCont <- sampleRanges(taxa,r = 0.5)
#bin the ranges into discrete time intervals
rangesDisc <- binTimeData(rangesCont,int.length = 1)
#note that we made interval lengths = 1: 
 	# thus q (per int) = q (per time) for make_durationFreqDisc

## Not run: 
#old ways of doing it (defunct as of paleotree version 2.6)
getSampRateCont(rangesCont)
getSampProbDisc(rangesDisc)

## End(Not run)

#new ways of doing it
    # we can constrain our functions
    # we can use parInit, parLower and parUpper
    # to control parameter bounds

#as opposed to getSampRateCont, we can do:
likFun <- make_durationFreqCont(rangesCont)
optim(parInit(likFun),
      likFun,
      lower = parLower(likFun),
      upper = parUpper(likFun),
      method = "L-BFGS-B",
      control = list(maxit = 1000000)
      )

#as opposed to getSampProbDisc, we can do:

likFun <- make_durationFreqDisc(rangesDisc)
optim(parInit(likFun),
      likFun,
      lower = parLower(likFun),
      upper = parUpper(likFun),
      method = "L-BFGS-B",
      control = list(maxit = 1000000)
      )

#these give the same answers (as we'd expect them to!)

#with newer functions we can constrain our functions easily
    # what if we knew the extinction rate = 0.1 a priori?
    
likFun <- make_durationFreqCont(rangesCont)
likFun <- constrainParPaleo(likFun,q.1~0.1)
optim(parInit(likFun),
      likFun,
      lower = parLower(likFun),
      upper = parUpper(likFun),
	     method = "L-BFGS-B",
	     control = list(maxit = 1000000)
	     )

#actually decreases our sampling rate estimate
   # gets further away from true generating value, r = 0.5 (geesh!)
   # but this *is* a small dataset...

Description

equation2function converts the right-hand side of an equation that can be written as a single line (like the right-hand side of an object of class formula) and creates an R function which calls the variables within as arguments and returns values consistent with the parameters of the input equation as written.

Usage

equation2function(equation, envir = parent.frame(), notName = "XXXXXXXXXXX")

Arguments

Details

This simple little function is rather 'hacky' but seems to get the job done, for a functionality that does not seem to be otherwise exist elsewhere in R.

Value

A function, with named blank (i.e. no default value) arguments.

Author(s)

David W. Bapst

Examples

# some simple examples
foo <- equation2function("x+y")
foo
foo(x = 4,y = 0.1)

foo <- equation2function("x+2*sqrt(2*y+3)^2")
foo
foo(x = 4,y = 0.1)

# what about weird long argument names and spaces
foo <- equation2function("stegosaur + 0.4 * P")
foo
foo(stegosaur = 5,P = 0.3)

Description

The following functions are for measuring and fitting various distributions to the gradual exhaustion of unexpressed character states, as originally described by Wagner (2000, Evolution).

Usage

accioExhaustionCurve(
  phyloTree,
  charData,
  charTypes = "unordered",
  outgroup = NULL,
  firstAppearances = NULL,
  missingCharValue = "?",
  inapplicableValue = "-"
)

accioBestAcquisitionModel(
  exhaustion_info,
  changesType,
  models = c("exponential", "gamma", "lognormal", "zipf")
)

charExhaustPlot(
  exhaustion_info,
  changesType,
  xlab = "Total Characters",
  ylab = NULL,
  main = NULL,
  xsize = 3
)

Arguments

Details

accioExhaustionCurve uses a Sankoff parsimony ancestral-reconstruction algorithm (written by P.J. Wagner, not the one from phangorn used elsewhere in paleotree) to calculate character changes across each branch (internode edge) of a tree, and then reports the counts of character state

accioBestAcquisitionModel takes output from accioExhaustionCurve, calculates one of two character change indices, and then fits a series of user-selected models to the dataset, returning details pertaining to the best-fit model.

charExhaustPlot is a wrapper for accioBestAcquisitionModel that produces a plot of the observed character change data against the expectation under the best-fit model.

The functions accioBestAcquisitionModel and charExhaustPlot offer users two different options for examining character change: totalAcc fits models to the total accumulated number of state changes over the phylogeny, thus using exhaustion to explore the size and distribution of character space. The other option charAlt fits models to the number of character that alter from primitive to derived over phylogeny, thus reflecting the size and distribution of state space.

accioExhaustionCurve can order its reconstruction of change by stratigraphic order of first appearances. It is unclear what this means.

Value

accioExhaustionCurve outputs a list containing two objects: first, a matrix named exhaustion consisting of three columns: "Steps", "Novel_States", and "Novel_Characters", respectively giving the counts of these respective values for each branch (internode edge). The second element of this list is named State_Derivations and is a count of how often each state, across all characters, was derived relative to the primitive position along each internode edge.

The output of accioBestAcquisitionModel is a list object containing information on the best-fit model, the parameters of that model, and the calculated probability distribution function for that model at the same intervals (for use in quantile plots).

charExhaustPlot produces a plot, and outputs no data.

Note

This family of functions presented here were originally written by Peter J. Wagner, and then modified and adapted by David W. Bapst for wider release in a CRAN-distributed package: paleotree. This makes the code presented here a very different beast than typical paleotree code, for example, there are fewer tests of correct input type, length, etc.

Author(s)

Initially written by Peter J. Wagner, with modification and documentation by David W. Bapst.

References

Wagner, P. J. 2000. Exhaustion of morphologic character states among fossil taxa. Evolution 54(2):365-386.

See Also

Also see paleotree functions minCharChange and ancPropStateMat, the latter of which is a wrapper for phangorn's function ancestral.pars.

Examples

# get data
data(SongZhangDicrano)

dicranoTree <- cladogramDicranoX13

# modify char data
charMat <- data.matrix(charMatDicrano)
charMat[is.na(charMatDicrano)] <- 0
charMat <- (charMat-1)
colnames(charMat) <- NULL
# replace missing values
charMat[is.na(charMatDicrano)] <- "?"

# the 'outgroup' is Exigraptus
   # also the first taxon listed in the matrix
exhaustionResults <- accioExhaustionCurve(
   phyloTree = dicranoTree,
   charData = charMat, charTypes = "unordered",
   outgroup = "Exigraptus_uniformis")

# fits models to exhaustion for total accumulation
accioBestAcquisitionModel(
   exhaustion_info = exhaustionResults,
   changesType = "totalAcc", 	
   models = c("exponential","gamma","lognormal","zipf")) 

# plot of exhaustion of total accumulation of character states
charExhaustPlot(exhaustion_info = exhaustionResults,
	   changesType = "totalAcc")

# plot of exhaustion of character alterations
charExhaustPlot(exhaustion_info = exhaustionResults,
	   changesType = "charAlt")

Description

This function takes a tree composed of higher-level taxa and a vector of lower-level taxa belonging to the set of higher-level taxa included in the input tree and produces a tree composed of the lower-level taxa, by treating the higher-level taxa as unresolved monophyletic polytomies. A user can also mark higher taxa as paraphyletic such that these are secondarily collapsed and do not form monophyletic clades in the output tree.

Usage

expandTaxonTree(
  taxonTree,
  taxaData,
  collapse = NULL,
  keepBrLen = FALSE,
  plot = FALSE
)

Arguments

Details

The output tree will probably be a rough unresolved view of the relationships among the taxa, due to the treatment of higher-level taxa as polytomies. This is similar to the methods used in Webb and Donoghue (2005) and Friedman (2009). Any analyses should be done by resolving this tree with multi2di in the ape package or via the various time-scaling functions found in this package (paleotree).

The taxaData vector should have one element per lower-level taxon that is to be added to the tree. The name of each element in the vector should be the names of the lower-level taxa, which will be used as new tip labels of the output lower-taxon tree. There should be no empty elements! Otherwise, expandTaxonTree won't know what to do with taxa that aren't being expanded.

By default, all higher-level taxa are treated as monophyletic clades if not otherwise specified. The collapse vector can (and probably should) be used if there is doubt about the monophyly of any higher-level taxa included in the input taxon-tree, so that such a group would be treated as a paraphyletic group in the output tree.

Also by default, the output tree will lack branch lengths and thus will not be dated, even if the input phylogeny is dated. If keepBrLen = TRUE, then the tree's edge lengths are kept and new taxa are added as zero length branches attaching to a node that represents the previous higher-taxon. This tree is probably not useful for most applications, and may even strongly bias some analyses. USE WITH CAUTION! The collapse argument, given as a vector, will cause such edges to be replaced by zero-length branches rather than fully collapsing them, which could have odd effects. If collapse is not NULL and keepBrLen = TRUE, a warning is issued that the output probably won't make much sense at all.

Value

Outputs the modified tree as an object of class phylo, with the higher-level taxa expanded into polytomies and the lower-level taxa as the tip labels.

Author(s)

David W. Bapst

References

Friedman, M. 2009 Ecomorphological selectivity among marine teleost fishes during the end-Cretaceous extinction. Proceedings of the National Academy of Sciences 106(13):5218–5223.

Webb, C. O., and M. J. Donoghue. 2005 Phylomatic: tree assembly for applied phylogenetics. Molecular Ecology Notes 5(1):181–183.

See Also

multi2di, bind.tree

Examples

set.seed(444)
# lets make our hypothetical simulated tree of higher taxa
taxtr <- rtree(10)
# taxa to place within higher taxa
taxd <- sample(taxtr$tip.label, 30, replace = TRUE)	
names(taxd) <- paste(taxd,"_x", 1:30, sep = "")
coll <- sample(taxtr$tip.label,3)		#what to collapse?
expandTaxonTree(taxonTree = taxtr, taxaData = taxd, 
                collapse = coll, plot = TRUE)

Description

Modifying a dated tree with $root.time elements often changes the actual timing of the root relative to the tips, such as when dropping tips, extending branches, or shift node ages backwards. When such modifications occur, the function fixRootTime can be used to find the correct root age. This function is mainly used as a utility function called by other tree-modifying functions discussed in the manual page for modifyTerminalBranches. This is typically performed via the function fixRootTime.

Usage

fixRootTime(
  treeOrig,
  treeNew,
  testConsistentDepth = TRUE,
  fixingMethod = "matchCladeTransferNodeAge"
)

Arguments

Value

Gives back a modified phylogeny as a phylo object, with a modified $root.time element.

Author(s)

David W. Bapst

See Also

modifyTerminalBranches, minBranchLength

Examples

#testing dropPaleoTip... and fixRootTime by extension

#simple example
tree <- read.tree(text = "(A:3,(B:2,(C:5,D:3):2):3);")
tree$root.time <- 10
plot(tree,no.margin = FALSE)
axisPhylo()

# now a series of tests, dropping various tips
(test <- dropPaleoTip(tree,"A")$root.time) #  = 7
(test[2] <- dropPaleoTip(tree,"B")$root.time) #  = 10
(test[3] <- dropPaleoTip(tree,"C")$root.time) #  = 10
(test[4] <- dropPaleoTip(tree,"D")$root.time) #  = 10
(test[5] <- dropPaleoTip(tree,c("A","B"))$root.time) #  = 5
(test[6] <- dropPaleoTip(tree,c("B","C"))$root.time) #  = 10
(test[7] <- dropPaleoTip(tree,c("A","C"))$root.time) #  = 7
(test[8] <- dropPaleoTip(tree,c("A","D"))$root.time) #  = 7

# is it all good? if not, fail so paleotree fails...
if(!identical(test,c(7,10,10,10,5,10,7,7))){stop("fixRootTime fails!")}

Description

This function calculates the intermediary values needed for fitting Foote's inverse survivorship analyses, as listed in the table of equations in Foote (2003), with the analyses themselves described further in Foote (2001) and Foote (2005).

Usage

footeValues(p, q, r, PA_n = 0, PB_1 = 0, p_cont = TRUE, q_cont = TRUE, Nb = 1)

Arguments

Details

Although most calculations in this function agree with the errata for Foote's 2003 table (see references), there were some additional corrections for the probability of D given FL (Prob(D|FL)) made as part of a personal communication in 2013 between the package author and Michael Foote.

Value

Returns a matrix with number of rows equal to the number of intervals (i.e. the length of p, q and r) and named columns representing the different values calculated by the function: "Nb", "Nbt", "NbL", "NFt", "NFL", "PD_bt", "PD_bL", "PD_Ft", "PD_FL", "PA", "PB", "Xbt", "XbL", "XFt" and "XFL".

Author(s)

David W. Bapst, with advice from Michael Foote.

References

Foote, M. 2001. Inferring temporal patterns of preservation, origination, and extinction from taxonomic survivorship analysis. Paleobiology 27(4):602-630.

Foote, M. 2003a. Origination and Extinction through the Phanerozoic: A New Approach. The Journal of Geology 111(2):125-148.

Foote, M. 2003b. Erratum: Origination and Extinction through the Phanerozoic: a New Approach. The Journal of Geology 111(6):752-753.

Foote, M. 2005. Pulsed origination and extinction in the marine realm. Paleobiology 31(1):6-20.

Examples

#very simple example with three intervals, same value for all parameters

# example rates (for the most part)
rate <- rep(0.1, 3)
                  
#all continuous
footeValues(rate,rate,rate)	

# origination pulsed
footeValues(rate,rate,rate,p_cont = FALSE)		 

# extinction pulsed
footeValues(rate,rate,rate,q_cont = FALSE) 	 

# all pulsed
footeValues(rate,rate,rate,p_cont = FALSE,q_cont = FALSE)

Description

Estimate per-interval sampling probability in the fossil record from a set of discrete-interval taxon ranges using the frequency-ratio method described by Foote and Raup (1996). Can also calculate extinction rate per interval from the same data distribution.

Usage

freqRat(timeData, calcExtinction = FALSE, plot = FALSE)

Arguments

Details

This function uses the frequency ratio ("freqRat") method of Foote and Raup (1996) to estimate the per-interval sampling rate for a set of taxa. This method assumes that intervals are of fairly similar length and that taxonomic extinction and sampling works similar to homogenous Poisson processes. These analyses are ideally applied to data from single stratigraphic section but potentially are applicable to regional or global datasets, although the behavior of those datasets is less well understood.

The frequency ratio is a simple relationship between the number of taxa observed only in a single time interval (also known as singletons), the number of taxa observed only in two time intervals and the number of taxa observed in three time intervals. These respective frequencies, respectively f1, f2 and f3 can then be related to the per-interval sampling probability with the following expression.

$Sampling.Probability = (f2^2)/(f1*f3)$

This frequency ratio is generally referred to as the 'freqRat' in paleobiological literature.

It is relatively easy to visually test if range data fits expectation that true taxon durations are exponentially distributed by plotting the frequencies of the observed ranges on a log scale: data beyond the 'singletons' category should have a linear slope, implying that durations were originally exponentially distributed. (Note, a linear scale is used for the plotting diagram made by this function when 'plot' is TRUE, so that plot cannot be used for this purpose.)

The accuracy of this method tends to be poor at small interval length and even relatively large sample sizes. A portion at the bottom of the examples in the help file examine this issue in greater detail with simulations. This package author recommends using the ML method developed in Foote (1997) instead, which is usable via the function make_durationFreqDisc.

As extant taxa should not be included in a freqRat calculation, any taxa listed as being in a bin with start time 0 and end time 0 (and thus being extant without question) are dropped before the model fitting it performed.

Value

This function returns the per-interval sampling probability as the "freqRat", and estimates

Author(s)

David W. Bapst

References

Foote, M. 1997 Estimating Taxonomic Durations and Preservation Probability. Paleobiology 23(3):278–300.

Foote, M., and D. M. Raup. 1996 Fossil preservation and the stratigraphic ranges of taxa. Paleobiology 22(2):121–140.

See Also

Model fitting methods in make_durationFreqDisc and make_durationFreqCont. Also see conversion methods in sProb2sRate, qsProb2Comp

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.1)
# Now let's use binTimeData to bin in intervals of 5 time units
rangesDisc <- binTimeData(rangesCont,int.length = 5)

# now, get an estimate of the sampling rate (we set it to 0.5 above)

# for discrete data we can estimate the sampling probability per interval (R)
    # i.e. this is not the same thing as the instantaneous sampling rate (r)
    
# can use sRate2sProb to see what we would expect
sRate2sProb(r = 0.1, int.length = 5)

# expect R = ~0.39

# now we can apply freqRat to get sampling probability
SampProb <- freqRat(rangesDisc, plot = TRUE)
SampProb

# I estimated R = ~0.25 
# Not wildly accurate, is it?

# can also calculate extinction rate per interval of time
freqRat(rangesDisc, calcExtinction = TRUE)

# est. ext rate = ~0.44 per interval
# 5 time-unit intervals, so ~0.44 / 5 = ~0.08 per time-unit
# That's pretty close to the generating value of 0.01, used in sampleRanges

## Not run: 
#################
# The following example code (which is not run by default) examines how 
	# the freqRat estimates vary with sample size, interval length
	# and compares it to using make_durationFreqDisc

# how good is the freqRat at 20 sampled taxa on avg?
set.seed(444)
r <- runif(100)
int.length = 1
# estimate R from r, assuming stuff like p = q
R <- sapply(r, sRate2sProb, int.length = 1)	
ntaxa <- freqRats <- numeric()
for(i in 1:length(r)){
	# assuming budding model
	record <- simFossilRecord(p = 0.1, 
	                          q = 0.1, 
	                          r = r[i], 
	                          nruns = 1,
	                          nSamp = c(15,25), 
	                          nExtant = 0, 
	                          plot = TRUE
	                          )
	ranges <- fossilRecord2fossilRanges(record)
	timeList <- binTimeData(ranges,int.length = int.length)
	ntaxa[i] <- nrow(timeList[[2]])
	freqRats[i] <- freqRat(timeList)
	}
plot(R,freqRats);abline(0,1)
# without the gigantic artifacts bigger than 1...
plot(R,freqRats,ylim = c(0,1));abline(0,1)
# very worrisome lookin'!

# how good is it at 100 sampled taxa on average?
set.seed(444)
r <- runif(100)
int.length = 1
R <- sapply(r,sRate2sProb,int.length = 1)
ntaxa <- freqRats <- numeric()
for(i in 1:length(r)){
    # assuming budding model
    record <- simFossilRecord(p = 0.1, 
                              q = 0.1, 
                              r = r[i], 
                              nruns = 1, 
                              nSamp = c(80,150), 
                              nExtant = 0, 
                              plot = TRUE)
	ranges <- fossilRecord2fossilRanges(record)
	timeList <- binTimeData(ranges, 
	     int.length = int.length)
	ntaxa[i] <- nrow(timeList[[2]])
	freqRats[i] <- freqRat(timeList)
	}
	
plot(R, freqRats,
     ylim = c(0,1)
     )
abline(0,1)

#not so hot, eh?

################
#LETS CHANGE THE TIME BIN LENGTH!

# how good is it at 100 sampled taxa on average, with longer time bins?
set.seed(444)
r <- runif(100)
int.length <- 10
R <- sapply(r, sRate2sProb, int.length = int.length)
ntaxa <- freqRats <- numeric()
for(i in 1:length(r)){
	   # assuming budding model
	   record <- simFossilRecord(p = 0.1, 
	                          q = 0.1, 
	                          r = r[i], 
	                          nruns = 1,	
	                          nSamp = c(80,150), 
	                          nExtant = 0, 
	                          plot = TRUE)
	   ranges <- fossilRecord2fossilRanges(record)
	   timeList <- binTimeData(ranges, int.length = int.length)
	   ntaxa[i] <- nrow(timeList[[2]])
	   freqRats[i] <- freqRat(timeList)
	   }
	
plot(R, freqRats, ylim = c(0,1))
abline(0,1)
# things get more accurate as interval length increases... odd, eh?

# how good is it at 20 sampled taxa on average, with longer time bins?
set.seed(444)
r <- runif(100)
int.length <- 10
R <- sapply(r, sRate2sProb, int.length = int.length)
ntaxa <- freqRats <- numeric()
for(i in 1:length(r)){
	# assuming budding model
	record <- simFossilRecord(p = 0.1, 
	                          q = 0.1, 
	                          r = r[i], 
	                          nruns = 1,	
	                          nSamp = c(15,25), 
	                          nExtant = 0, 
	                          plot = TRUE)
	ranges <- fossilRecord2fossilRanges(record)
	timeList <- binTimeData(ranges, int.length = int.length)
	ntaxa[i] <- nrow(timeList[[2]])
	freqRats[i] <- freqRat(timeList)
	}
plot(R, freqRats, ylim = c(0,1))
abline(0,1)
# still not so hot at low sample sizes, even with longer bins

########################
# ML METHOD

# how good is the ML method at 20 taxa, 1 time-unit bins?
set.seed(444)
r <- runif(100)
int.length <- 1
R <- sapply(r,sRate2sProb,int.length = int.length)
ntaxa <- ML_sampProb <- numeric()
for(i in 1:length(r)){
	   # assuming budding model
	   record <- simFossilRecord(p = 0.1, 
	                             q = 0.1, 
	                             r = r[i], 
	                             nruns = 1,
	                             nSamp = c(15,25), 
	                             nExtant = 0, 
	                             plot = TRUE
	                             )
	   ranges <- fossilRecord2fossilRanges(record)
	   timeList <- binTimeData(ranges, int.length = int.length)
	   ntaxa[i] <- nrow(timeList[[2]])
    likFun <- make_durationFreqDisc(timeList)
    ML_sampProb[i] <- optim(
         parInit(likFun), likFun,
		      lower = parLower(likFun),
		      upper = parUpper(likFun),
         method = "L-BFGS-B",
         control = list(maxit = 1000000)
         )[[1]][2]
	   }
	   
plot(R, ML_sampProb)
abline(0,1)

# Not so great due to likelihood surface ridges
 # but it returns values between 0-1

# how good is the ML method at 100 taxa, 1 time-unit bins?
set.seed(444)
r <- runif(100)
int.length <- 1
R <- sapply(r, sRate2sProb,
    int.length = int.length)
ntaxa <- ML_sampProb <- numeric()
for(i in 1:length(r)){
    # assuming budding model
	   record <- simFossilRecord(p = 0.1, 
	                             q = 0.1, 
	                             r = r[i], 
	                             nruns = 1,
	                             nSamp = c(80,150), 
	                             nExtant = 0, 
	                             plot = TRUE)
	   ranges <- fossilRecord2fossilRanges(record)
	   timeList <- binTimeData(ranges,int.length = int.length)
	   ntaxa[i] <- nrow(timeList[[2]])
    likFun <- make_durationFreqDisc(timeList)
    ML_sampProb[i] <- optim(parInit(likFun),
                         likFun,
		                      lower = parLower(likFun),
		                      upper = parUpper(likFun),
                         method = "L-BFGS-B",
                         control = list(maxit = 1000000)
                         )[[1]][2]
	}
	
plot(R,ML_sampProb)
abline(0,1)

# Oh, fairly nice, although still a biased uptick as R gets larger

## End(Not run)

Description

The Paleobiology Database API (link) is very easy to use, and generally any data one wishes to collect can be obtained in R through a variety of ways - the simplest being to wrap a data retrieval request to the API, specified for CSV output, with R function read.csv. The functions listed here, however, are some simple helper functions for doing tasks common to users of this package - downloading occurrence data, or taxonomic information, for particular clades, or for a list of specific taxa.

Usage

getCladeTaxaPBDB(
  taxon,
  showTaxa = c("class", "parent", "app", "img", "entname"),
  status = "accepted",
  urlOnly = FALSE,
  stopIfMissing = FALSE,
  failIfNoInternet = TRUE
)

getSpecificTaxaPBDB(
  taxa,
  showTaxa = c("class", "parent", "app", "img", "entname"),
  status = "accepted",
  urlOnly = FALSE,
  stopIfMissing = FALSE,
  failIfNoInternet = TRUE
)

getPBDBocc(
  taxa,
  showOccs = c("class", "classext", "subgenus", "ident", "entname"),
  failIfNoInternet = TRUE
)

Arguments

Details

In many cases, it might be easier to write your own query - these functions are only made to make getting data for some very specific applications in paleotree easier.

Value

These functions return a data.frame containing variables pulled for the requested taxon selection. This behavior can be modified by argument urlOnly.

Author(s)

David W. Bapst

References

Peters, S. E., and M. McClennen. 2015. The Paleobiology Database application programming interface. Paleobiology 42(1):1-7.

See Also

See makePBDBtaxonTree, makePBDBtaxonTree, and plotPhyloPicTree for functions that use taxonomic data. Occurrence data is sorted by taxon via taxonSortPBDBocc, and further utilized occData2timeList and plotOccData.

Examples

# Note that all examples here use argument 
    # failIfNoInternet = FALSE so that functions do
    # not error out but simply return NULL if internet
    # connection is not available, and thus
    # fail gracefully rather than error out (required by CRAN).
# Remove this argument or set to TRUE so functions fail
    # when internet resources (paleobiodb) is not available.

#graptolites
graptData <- getCladeTaxaPBDB("Graptolithina", 
    failIfNoInternet = FALSE)
dim(graptData)
sum(graptData$taxon_rank == "genus")

# so we can see that our call for graptolithina returned 
    # a large number of taxa, a large portion of which are
    # individual genera
# (554 and 318 respectively, as of 03-18-19)

tetrapodList<-c("Archaeopteryx", "Columba", "Ectopistes",
   "Corvus", "Velociraptor", "Baryonyx", "Bufo",
   "Rhamphorhynchus", "Quetzalcoatlus", "Natator",
   "Tyrannosaurus", "Triceratops", "Gavialis",
   "Brachiosaurus", "Pteranodon", "Crocodylus",
   "Alligator", "Giraffa", "Felis", "Ambystoma",
    "Homo", "Dimetrodon", "Coleonyx", "Equus",
   "Sphenodon", "Amblyrhynchus")

tetrapodData <-getSpecificTaxaPBDB(tetrapodList, 
    failIfNoInternet = FALSE)
dim(tetrapodData)
sum(tetrapodData$taxon_rank == "genus")
# should be 26, with all 26 as genera

#############################################
# Now let's try getting occurrence data

# getting occurrence data for a genus, sorting it
# Dicellograptus
dicelloData <- getPBDBocc("Dicellograptus", 
    failIfNoInternet = FALSE)

if(!is.null(dicelloData)){

dicelloOcc2 <- taxonSortPBDBocc(dicelloData, 
    rank = "species", onlyFormal = FALSE, 
    failIfNoInternet = FALSE)
names(dicelloOcc2)

}

Description

This dataset contains a morphological character matrix (containing both a set of 45 discrete characters, and 4 continuous characters coded as minimum and maximum range values), along with biostratigraphic range data for 183 graptoloid species-level taxa from Bapst et al. (2012, PNAS). Also includes a pre-calculated distance matrix based on the character matrix, using the algorithm applied by Bapst et al (2012). Interval dates for biostratigraphic zones is taken from Sadler et al. 2011.

Format

This dataset is composed of three objects:

graptCharMatrix

A matrix composed of mixed character data and a group code for 183 graptoloid taxa, with rows named with species names and columns named with character names.

graptRanges

A list containing two matrices: the first matrix describes the first and last interval times for 20 graptolite biozones and the second matrix contains the first and last appearances of 183 graptolite species in those same biozones. (In other words, graptRanges has the timeList format called by some paleotree functions).

graptDistMat

A 183x183 matrix of pair-wise distances (dissimilarities) for the 183 graptolite species, using the algorithm for discrete characters and min-max range values described in Bapst et al.

Details

The character matrix contains characters of two differing types with a (very) small but non-negligible amount of missing character data for some taxa. This required the use of an unconventional ad hoc distance metric for the published analysis, resulting in a (very slightly) non-Euclidean distance matrix. This breaks some assumptions of some statistical analyses or requires special corrections, such as with PCO.

Note that taxonomic data were collected only for species present within an interval defined by the base of the Uncinatus biozone (~448.57 Ma) to the end of the cyphus biozone (~439.37 Ma). Many taxa have first and last appearance dates listed in graptRanges which are outside of this interval (see examples).

Source

Source for stratigraphic ranges and character data:

Bapst, D. W., P. C. Bullock, M. J. Melchin, H. D. Sheets, and C. E. Mitchell. 2012. Graptoloid diversity and disparity became decoupled during the Ordovician mass extinction. Proceedings of the National Academy of Sciences 109(9):3428-3433.

Source for interval dates for graptolite zones:

Sadler, P. M., R. A. Cooper, and M. Melchin. 2009. High-resolution, early Paleozoic (Ordovician-Silurian) time scales. Geological Society of America Bulletin 121(5-6):887-906.

See Also

For more example graptolite datasets, see retiolitinae

This data was added mainly to provide an example dataset for nearestNeighborDist

Examples

#load data
data(graptDisparity)

#separate out two components of character matrix

#45 discrete characters
discChar <- graptCharMatrix[,1:45]

#min ranges for 4 continuous characters
cMinChar <- graptCharMatrix[,c(46,48,50,52)]
#max ranges for 4 continuous characters
cMaxChar <- graptCharMatrix[,c(47,49,51,53)]

#group (clade/paraclade) coding
groupID <- graptCharMatrix[,54]

#number of species
nspec <- nrow(graptCharMatrix)

#some plotting information from Bapst et al.'s plotting scripts
grpLabel <- c("Normalo.","Monogr.","Climaco.",
		"Dicrano.","Lasiogr.","Diplogr.","Retiol.")
grpColor <- c("red","purple",colors()[257],colors()[614],
		colors()[124],"blue",colors()[556])

##########

#plot diversity curve of taxa
taxicDivDisc(graptRanges)

#but the actual study interval for the data is much smaller
abline(v = 448.57,lwd = 3) #start of study interval
abline(v = 439.37,lwd = 3) #end of study interval

#plot diversity curve just for study interval
taxicDivDisc(graptRanges, timelims = c(448.57,439.37))

############

#distance matrix is given as graptDistMat
   #to calculate yourself, see code below in DoNotRun section

#now, is the diagonal zero? (it should be)
all(diag(graptDistMat) == 0)

#now, is the matrix symmetric? (it should be)
isSymmetric(graptDistMat)

#can apply cluster analysis
clustRes <- hclust(as.dist(graptDistMat))
plot(clustRes,labels = FALSE)

#use ape to plot with colors at the tips
dev.new(width = 15) 	# for a prettier plot
plot.phylo(as.phylo(clustRes),show.tip.label = FALSE,
		no.margin = TRUE,direction = "upwards")
tiplabels(pch = 16,col = grpColor[groupID+1])
legend("bottomright",legend = grpLabel,col = grpColor,pch = 16)
dev.set(2)

#can apply PCO (use lingoes correction to account for negative values
   #resulting from non-euclidean matrix
pco_res <- pcoa(graptDistMat,correction = "lingoes")

#relative corrected eigenvalues
rel_corr_eig <- pco_res$values$Rel_corr_eig
layout(1:2)
plot(rel_corr_eig)
#cumulative
plot(cumsum(rel_corr_eig))

#first few axes account for very little variance!!



#well let's look at those PCO axes anyway
layout(1)
pco_axes <- pco_res$vectors
plot(pco_axes[,1],pco_axes[,2],pch = 16,col = grpColor[groupID+1],
   xlab = paste("PCO Axis 1, Rel. Corr. Eigenvalue  = ",round(rel_corr_eig[1],3)),
   ylab = paste("PCO Axis 2, Rel. Corr. Eigenvalue  = ",round(rel_corr_eig[2],3)))
legend("bottomright",legend = grpLabel,col = grpColor,pch = 16,ncol = 2,cex = 0.8)


##########m##############

## Not run:  

#calculate a distance matrix (very slow!)
#Bapst et al. calculated as # char diffs / total # of chars
   #but both calculated for only non-missing characters for both taxa
#non-identical discrete states = difference for discrete traits
#non-overlapping ranges for continuous characters = difference for cont traits

distMat <- matrix(,nspec,nspec)
rownames(distMat) <- colnames(distMat) <- rownames(graptCharMatrix)
for(i in 1:nspec){ for(j in 1:nspec){    #calculate for each pair of species
   #discrete characters
   di <- discChar[i,]     #discrete character vector for species i
   dj <- discChar[j,]     #discrete character vector for species j
   #now calculate pair-wise differences for non-missing characters
   discDiff <- (di != dj)[!is.na(di)&!is.na(dj)] #logical vector
   #
   #continuous characters: need another for() loop
   contDiff <- numeric()
   for(ct in 1:4){
       #if they do not overlap, a min must be greater than a max value
       contDiff[ct] <- cMinChar[i,ct]>cMaxChar[j,ct] | cMinChar[j,ct]>cMaxChar[i,ct]
       }
   #remove NAs
   contDiff <- contDiff[!is.na(contDiff)]
   #combine
   totalDiff <- c(discDiff,contDiff)
   #divide total difference 
   distMat[i,j] <- sum(totalDiff)/length(totalDiff)
   }}

#but is it identical to the distance matrix already provided?
identical(distMat,graptDistMat)
#ehh, numerical rounding issues...

#A somewhat speeder alternative to calculate a distance matrix
distMat <- matrix(,nspec,nspec)
rownames(distMat) <- colnames(distMat) <- rownames(graptCharMatrix)
for(i in 1:(nspec-1)){ for(j in (i+1):nspec){    #calculate for each pair of species
   #now calculate pair-wise differences for non-missing characters
   discDiff <- (discChar[i,] != discChar[j,])[
      !is.na(discChar[i,])&!is.na(discChar[j,])] #logical vector
   #continuous characters: if they do not overlap, a min must be greater than a max value
      contDiff <- sapply(1:4,function(ct)
            cMinChar[i,ct]>cMaxChar[j,ct] | cMinChar[j,ct]>cMaxChar[i,ct])
   #remove NAs, combine, divide total difference 
   distMat[i,j] <- distMat[j,i] <- sum(c(discDiff,contDiff[!is.na(contDiff)]))/length(
        c(discDiff,contDiff[!is.na(contDiff)]))
   }}
diag(distMat) <- 0

#but is it identical to the distance matrix already provided?
identical(distMat,graptDistMat)
#ehh, MORE numerical rounding issues...


## End(Not run)