What is the right way to identify unique species in Paleobiology Database downloads?

(Continuing with the RMarkdown-to-Blogger experiment! I apologize for any mysterious indentations that I can’t seem to make go away! Sorry!)

So, let’s recap:

Last post, we had some fun pulling data out of the Paleobiology Database using R package paleobioDB (Varela et al., in press). We even made a nice plot of age uncertainty in occurrences sorted by species in the delightful graptolite genus Dicellograptus. Let’s see that pretty plot again, although we’ll hide some of the R code this time!

library(paleobioDB)

dicelloData<-pbdb_occurrences(limit="all", base_name="Dicellograptus", show=c("phylo","ident"))
dicelloData<-dicelloData[dicelloData$rnk==3,]   #keep only occurrences of taxa resolved to species level

plotOccPBDB(dicelloData,"Dicellograptus Species")

Now, we had a little mystery when we decided to look at species ranges from this same data using function pbdb_temp_range in package paleobioDB.

pbdb_temp_range(dicelloData, rank="species", do.plot=FALSE)

##                              max   min
## Dicellograptus alector     471.8 443.4
## Dicellograptus anceps      455.8 443.7
## Dicellograptus ornatus     453.0 443.4
## Dicellograptus complanatus 453.0 443.7

The occurrences I plotted above show something like more than twenty species. pbdb_temp_range shows only four species. What is going on?

Well, for the occurrence plot above, I used the species_name column (the column $ids, see here for its description in the paleobioDB API documentation) to distinguish different species, which we can see contains more than 20 unique identifiers:

as.character(unique(dicelloData$ids))

##  [1] "anceps"       "alector"      "gurleyi"      "mensurans"   
##  [5] "sextans"      "complanatus"  "vagus"        "intortus"    
##  [9] "divaricatus"  "caduceus"     "forchammeri"  "moffatensis" 
## [13] "intermedius"  "russonioides" "johnstrupi"   "morrisi"     
## [17] "elegans"      "smithi"       "angulatus"    "ornatus"     
## [21] "flexuosus"    "minor"        "tumidus"      "turgidus"    
## [25] "mirabilis"

However, while pbdb_temp_range uses the various ‘name’ entries for identifying independent taxa at other taxonomic levels, pbdb_temp_range uses taxon_no (column $tid) to identify unique species-level taxa. We can see this in the code in paleobioDB’s GitHub repo.

This is suggestive that taxon_no is the Paleobiology Database’s way of handling taxonomic synonymy among species. A quick survey of a few friends who actually use PBDB data on a regular basis (unlike me, I’m a PBDB-newb) reveals some… uh, lack of clarity on this point. Well, a lack of clarity often suggests lack of clear documentation. Indeed, the paleobioDB API documentation describes taxon_no simply as the “unique identifier of the identified taxonomic name” and is remarkably silent on suggesting whether it might indicate unique species-level identities. Looking through the rest of the API documentation, its never said explicitly, but one might infer that this is true, as long as taxonomic rank (column $rnk) has been been limited to species-level taxa only. I.e., all species-level taxa with unique taxon-IDs are valid species-level taxa.

I wonder though: aren’t there cases of synonymized supraspecific taxa? These must be missed if uniqueness at those levels in paleobioDB package is determined based on the various ‘names’ qualifiers. This issue seems like quite the tough nut. It seems properly sorting data out of an occurrances download from the PBDB may require directly referencing the taxonomic database of the PBDB as well.

Anyway, we can see that $tid matches the species returned by pbdb_temp_range, which finds its names for the unique tid values using match(), which just means the species name reported by pbdb_temp_range is whatever species name is listed first with that tid value in a given data table.

unique(dicelloData$tid)

## [1] 306364  33650 306226 306367

as.character(sapply(unique(dicelloData$tid),function(x) dicelloData$ids[match(x,dicelloData$tid)]))

## [1] "anceps"      "alector"     "complanatus" "ornatus"

Well, that solves that mystery, but is tid actually the right way to identify unique valid species in a dataset? I wonder, as this isn’t written down anywhere…

I decided to pursue this further and went to Fossilworks, a mirror database of the Paleobiology Database which I mentioned in the previous post. I used the download form at Fossilworks to pull the species ranges for Dicellograptus. It returns ranges for 27 separate species; here’s a summary showing just the first four columns and small selection of ranges:

genus	species	base of range (Ma)	top of range (Ma)
Dicellograptus	alector	456.1	443.7
Dicellograptus	anceps	449.5	445.6
Dicellograptus	angulatus	456.1	449.5
Dicellograptus	caduceus	460.9	456.1
Dicellograptus	complanatus	449.5	445.6
Dicellograptus	divaricatus	468.1	456.1
…	…	…	…
Dicellograptus	vagus	460.9	456.1

I also get 27 species with the download form using PBDB classic at paleobioDB (also discussed in the last post). So, these applications are clearly using species_name and not taxon_no to distinguish what species are in the occurrences data.

Okay, I thought, this is, uh… inconsistent.

Recently, the Paleobiology Database revealed a new interface, called PBDB Navigator, which uses a map and various chronological and taxonomic filters to sort through the PBDB collection data. How many species does Navigator report for Dicellograptus?

Taxon name (number of occurrences in the database)
    Dicellograptus complanatus  (1) 
    Dicellograptus ornatus  (1) 
    Dicellograptus anceps  (1) 
    Dicellograptus gravis  (1) 

Huh. Okay, so Navigator uses taxon_no to identify individual unique species, unlike the Classic PBDB which is on the same website but uses species_name.

This is making my head hurt. What’s the answer? Maybe someone could pipe up in the comments.

Anyway, until next time!

How Do We Treat Fossil Age Data? Dates, Ranges and Occurrences

(Hello all! I’m trying Rmarkdown for writing blog posts! Let’s see how this goes.)

I spend a lot of time thinking about how the datasets we (paleontologists and biologists) work with differ in their structure and information content. The obvious difference is how much paleontological data differs from the sort of data we have for living organisms in evolutionary biology (neontology… whatever…) but the one we so often overlook is how datasets from the fossil record can differ greatly from each other. I talk a lot in my book chapter last year on how much paleontological datasets for different groups of organisms or time-intervals can differ and how this impacts our intention to try phylogeny-based analyses of macroevolution (you can check out my chapter here, thanks to the publisher, including its horrid easter eggs referring to a certain popular series of Game Boy games).

That chapter barely scratched the surface in terms of ‘things to know about paleo data’. Recently, while playing with Paleobiology Database data using the recently release R package paleobioDB, I had some additional thoughts on how we think of and use chronostratigraphic data for fossil taxa.

Warning: Many of you who are paleontology are probably already familiar with a lot of what I’m about to say. In some regards, I’m very mindful that some of my biologist friends read this, so I’m trying to cover what they might not be so familiar with. Also, much of it is a retread of previous blog posts. A major goal of this blog post is an attempt to describe a set of terms for relating how different types of fossil age datasets relate to one another in terms of their information content (in a sense, their ontology), and that requires discussing some of the tiniest of minutia. I also don’t think I’m right about everything I say below–I think there is a lot more thinking we need to have about the ‘ontology’ of paleontological age data.

Four-Date Age Data (“timeList”)

I’ve thought a few ways about to talk about these issues, and the best I can do is to start with how I encountered taxonomic age data initially. So, graptoloids (the planktonic graptolites, my favorite group… why aren’t they yours too?) have a very well known fossil record. Skipping through some of the vagueries of Paleozoic time-scale making, you can generally tell based upon the species composition of a graptolite fossil assemblage which graptolite biostratigraphic zone you are in, meaning that the rocks you’re looking at come before ‘this graptolite’ showing up for the first time but after some other graptolite species shows up. The same taxa are often found across multiple continents, so there is a global correlation of such zones. The use of some complicated annealing software allows us to put absolute dates on the first appearance of graptolite taxa (Sadler et al., 2009), so if you know the graptolite zone of a rock, you know what little segment of (about) a few million years that rock is from. What nice little time-keepers… Its really too bad graptoloids are extinct, isn’t it?

Anyway, most of the data recorded in papers referring to when particular graptoloid taxa appear in the geologic record is often just which graptolite zones these species or genera first and last appear within. This is understandable, as some taxa can be found pretty continuously throughout their stratigraphic range, and all the information important for stratigraphic correlation are in these times of first and last appearance. In many cases, the data is also very well-behaved in that these graptolite zones don’t overlap and are of a rougly equal length (although maybe with more variation than we’d like; Sadler et al., 2009). This is exactly the sort of data I was first exposed to and dealt with in some undergrad research (published in Bapst et al., 2012). You can find the data here), or you can play with it in paleotree because I recently added it as example data. Let’s take a look at that sort of data.

library(paleotree)

## Loading required package: ape

data(graptDisparity)

Our item of interest is graptRanges, which is a timeList object I’ve talked about on this blog before and composed of two matrices: a matrix of the earliest and latest ages for some intervals, in this case graptolite biostratigraphic zones, which are intervals denoted by the first appearance of the taxa the zones are named after.

head(graptRanges[[1]])

##                                     start_time end_time
## Nemagraptus gracilis (Gi1)              460.86   456.35
## Orthograptus calcaratus (Gi2)           456.35   455.29
## Diplograptus lanceolatus (Ea1)          455.29   452.21
## Diplacanthograptus spiniferus (Ea2)     452.21   449.73
## Dicellograptus kirki (Ea3)              449.73   448.96
## Dicranograptus gravis (Ea4)             448.96   448.57

…and first and last intervals of appearance for some taxa.

head(graptRanges[[2]])

##                                 first_int last_int
## 'Bulmanograptus' macilentus            16       16
## 'Monograptus' arciformis               16       17
## 'Monograptus' austerus                 15       17
## 'Paramplexograptus' kiliani            12       13
## 'Paramplexograptus' paucispinus        12       13
## 'Prisitiograptus' fragilis             16       19

Information-wise, this means every taxon has four dated values associated with it: the earliest and latest dates for the intervals in which a taxon first and last appears. By separating the two into separate matrices (effectively making the data structure no longer ‘flat’) we remove some redundancy of having to list the same dates for different taxa which first or last appear in the same interval. Essentially, we’re compressing the information content of our taxon range data. Losing this redundancy is good because it minimizes needing to update the dates more than once.

Note that when we have consecutive, non-overlapping intervals (like with the graptoloid data above), its possible we could compress the information content even further; for example, we could just have the start dates of intervals listed in the interval times matrix, and their end-date is, by implication, the start date of the next interval listed. paleotree retains the second column however so that it can handle those cases where intervals do overlap or are non-consecutive. Going in the other direction, we might want more than two numbers to define dates for a given collection, if we had more information about dates than simply a minimum and maximum bound with a flat probability density inbetween (for example, if the fossils themselves can be dated using geochemistry, such as Strontium isotopes), but that sort of age information is not common.

Two-Date Age Data (“timeData”)

Particularly in discussions I had on this blog last year, I stressed that there was two other sorts of datasets that might be confused with this. If each taxon was known from only a single collection in the fossil record, such that their first appearance was their last appearance, than we would only have two dates known: the earliest and latest dates for that particular rock it was found in. Alternatively, in an extremely well resolved fossil record, we might have exceptional certainty about the dating of fossil remains, such the first and last appearance dates (FADs and LADs) of taxa were known from strata of extremely precise dating with negligble uncertainty.

Both situations are simply special cases of the four-date datasets; you could still describe them with four dates, but the values would largely be redundant. For comparability across datasets, though, it may be preferable to define taxa as

In both of these cases, each listed taxon has two dated values associated with it, but they imply very different things. I’ve encountered issues with people mistaking one for the other when it comes to providing input to paleotree functions, like timePaleoPhy, hence some recent (in the past year) changes to how timePaleoPhy handles age data (which the help file refers to as “timeData”) to make it more explicit. It doesn’t help that data reflecting age uncertainty of point observations is probably the most common type of data that people seem to have in-hand when attempting to time-scale phylogenies of fossil taxa, while I myself had in hand data of precise first and last appearances when I initially developed paleotree.

One could imagine cases where both of the above special cases are true: i.e., (1) taxa are known from singular remains and appear to not have any persistent morphotaxon duration and (2) the stratigraphic time-scale (this is sometimes called an ‘age model’) is so well-resolved that dates are known with exceptional precision; this means every fossil taxon’s age could be defined by a single dated value. Such datasets are very rare, and you won’t find anything that assumes such data as input in paleotree.

Ultimately, the commonly used functions for time-scaling phylogenies in the fossil record, my timePaleoPhy and strap’s DatePhylo use two-date age data like this. paleotree’s bin_timePaleoPhy and bin_cal3TimePaleoPhy are just wrappers which take the four-date timeList datasets, randomly sample ages within the intervals and hand their respective time-scaling functions two-date age data.

Occurrence Data

Although its possible to describe the age ranges of fossil taxa from many groups as the four-date-value type above, this format is a reduction in information content; a summary of only what is needed to know about the timing of their first and last appearances. Other than those cases where some taxon is found continuously from its first to last appearance, a given taxon may be observed to occur in a larger number of fossil collections. Thus for each of a given taxon, we might end up with two dates, reflecting the earliest and latest dates bounding the age of each collection. I’ll call this occurrence data (which I think is what everyone else calls it, but I’m not certain 100% because, well, paleontologists can be inconsistent with their terminology).

The Paleobiology Database is probably the data source most accessible to anyone who wants age data from the fossil record, and its raw data is in the form of such occurrences within collections. With the recently introduction of the new paleobioDB package (Varela et al., in press), the use of paleobioDB (aka PBDB) occurrence data is probably going to increase.

Let’s look at an example of occurrence data, using the paleobioDB package to download some PBDB records. To make this fun, we’ll look up this blog’s mascot…

library(paleobioDB)

## Loading required package: raster
## Loading required package: sp
## 
## Attaching package: 'raster'
## 
## The following objects are masked from 'package:ape':
## 
##     rotate, zoom
## 
## Loading required package: maps

nemaData<-pbdb_occurrences(limit="all", base_name="Nemagraptus", show=c("phylo","ident"))
nemaData<-nemaData[nemaData$rnk==3,]   #keep only occurrences of taxa resolved to species level
                                    
head(nemaData)

##         oid typ   cid                    tna rnk   tid               oei
## 1:4   93753 occ  6986   Nemagraptus gracilis   3 33700 Middle Ordovician
## 1:16 118184 occ  8776   Nemagraptus gracilis   3 33700           Caradoc
## 1:20 118736 occ  8887   Nemagraptus gracilis   3 33700          Actonian
## 1:22 118977 occ  9035   Nemagraptus gracilis   3 33700        Gisbornian
## 1:36 598661 occ 63445   Nemagraptus gracilis   3 33700        Gisbornian
## 1:38 601803 occ 63967 Nemagraptus ? gracilis   3 33700        Gisbornian
##        eag   lag   rid          odl   odn           cll   cln          phl
## 1:4  470.0 458.4   437 Graptoloidea 33606 Graptolithina 33534 Hemichordata
## 1:16 460.9 449.5   613 Graptoloidea 33606 Graptolithina 33534 Hemichordata
## 1:20 455.8 449.5   607 Graptoloidea 33606 Graptolithina 33534 Hemichordata
## 1:22 460.9 456.1   623 Graptoloidea 33606 Graptolithina 33534 Hemichordata
## 1:36 460.9 456.1 18179 Graptoloidea 33606 Graptolithina 33534 Hemichordata
## 1:38 460.9 456.1 18428 Graptoloidea 33606 Graptolithina 33534 Hemichordata
##        phn         idt      ids         mna mra    oli  rst  rss
## 1:4  33518 Nemagraptus gracilis Nemagraptus   5   <NA> <NA> <NA>
## 1:16 33518 Nemagraptus gracilis Nemagraptus   5   <NA> <NA> <NA>
## 1:20 33518 Nemagraptus gracilis Nemagraptus   5 Onnian <NA> <NA>
## 1:22 33518 Nemagraptus gracilis Nemagraptus   5   <NA> <NA> <NA>
## 1:36 33518 Nemagraptus gracilis Nemagraptus   5   <NA> <NA> <NA>
## 1:38 33518 Nemagraptus gracilis Nemagraptus   5   <NA>    ? <NA>

We can see each row of the obtained data-frame corresponds to a different individual occurrence of some named taxon. Myltiple occurrences are listed for the same taxon. We can simplify this by looking at just the taxon name listed for the occurrence and the earliest and latest age bounds for the associated collections.

head(nemaData)[,c("tna","eag","lag")]

##                         tna   eag   lag
## 1:4    Nemagraptus gracilis 470.0 458.4
## 1:16   Nemagraptus gracilis 460.9 449.5
## 1:20   Nemagraptus gracilis 455.8 449.5
## 1:22   Nemagraptus gracilis 460.9 456.1
## 1:36   Nemagraptus gracilis 460.9 456.1
## 1:38 Nemagraptus ? gracilis 460.9 456.1

We can already see that many of these occurrences overlap with each other and are not from intervals of the same size. Now, PBDB data is noisy and some of that is due to noise, but some of it is just due to differences in how well collections can be chronostratigraphically constrained. So, here we have a data set where taxa may be represented by hundreds of date-values, effectively two date bounds for every collection a taxon occurs in.

We could try visualizing this; first let’s group the occurrences by species…

occList<-lapply(unique(nemaData$ids),function(x) nemaData[x==nemaData$ids,c("eag","lag"), drop=FALSE])
occList<-lapply(occList,function(x) x[order(x[,1]),, drop=FALSE])
names(occList)<-unique(nemaData$ids)

print(occList)

## $gracilis
##        eag   lag
## 1:20 455.8 449.5
## 1:48 458.4 453.0
## 1:16 460.9 449.5
## 1:22 460.9 456.1
## 1:36 460.9 456.1
## 1:38 460.9 456.1
## 1:41 460.9 456.1
## 1:43 460.9 456.1
## 1:45 460.9 449.5
## 1:46 460.9 449.5
## 1:49 460.9 456.1
## 1:4  470.0 458.4
## 
## $exilis
##        eag   lag
## 1:39 460.9 456.1
## 1:40 460.9 456.1
## 1:42 460.9 456.1
## 2    460.9 449.5

And now let’s plot the bounds of the age uncertainty for each of the occurrences, color-coded for their species.

par(yaxt="n")
plot(0, 0, type="n", xlim=c(max(sapply(occList,max))*1.02, min(sapply(occList,min)))*.98,
     ylim=c(0,sum(sapply(occList,nrow))+(2*length(occList))),
     ylab = "", xlab = "Time (Ma)", main="Age Uncertainty of PBDB Occurrences for Nemagraptus Species")

occColors<-rainbow(length(occList))

count<-1
for(i in 1:length(occList)){
  for(j in 1:nrow(occList[[i]])){
    lines(occList[[i]][j,],c(count,count),lwd=3,col=occColors[i])
    count<-count+1
  }
  count<-count+2
  }

That wasn’t so impressive, as Nemagraptus is a rather species-poor graptolite genus. Let’s try its close dicranograptid cousin genus, Dicellograptus. This time, we’ll make the plotting code a function.

plotOccPBDB<-function(pbdbOccData,groupLabel=NULL,occColors=NULL,lineWidth=NULL,xlims=NULL){
  #where pbdbOccData is a matrix returned by function pbdb_occurrences in package paleobioDB
  #
  #order taxa by species (add more taxonomic levels later?)
  occList<-lapply(unique(pbdbOccData$ids),function(x) pbdbOccData[x==pbdbOccData$ids,c("eag","lag"), drop=FALSE])
  #
  #
  #order taxa by earliest occurrence
  occList<-occList[order(sapply(occList,max))]
  #order the occurrences within taxa
  occList<-lapply(occList,function(x) x[order(x[,1]),, drop=FALSE])
  names(occList)<-unique(pbdbOccData$ids)
  #
  #set xlims
  if(is.null(xlims)){
   xlims<-c(max(sapply(occList,max)), min(sapply(occList,min))) 
   xlimMod<-(xlims[1]-xlims[2])*0.01
   xlims<-c(xlims[1]+xlimMod,xlims[2]-xlimMod)
  }
  #initiate the plot with a modifiable main title
  par(yaxt="n")
  plot(0, 0, type="n", xlim=xlims,
       ylim=c(0,sum(sapply(occList,nrow))+(2*length(occList))),
       ylab = "", xlab = "Time (Ma)",
       main=ifelse(is.null(groupLabel),"Age Uncertainty of PBDB Occurrences",
             paste("Age Uncertainty of PBDB Occurrences for",groupLabel)))
  #set colors
  if(is.null(occColors)){
    occColors<-sample(rainbow(length(occList)))  #scramble colors
    }
  #set line width
  if(is.null(lineWidth)){
    lineWidth<-(-0.02*nrow(pbdbOccData))+3.2
    lineWidth<-ifelse(lineWidth>0,lineWidth,0.01)
    }
  #now plot the occurrences as lines
  count<-1
  for(i in 1:length(occList)){
    for(j in 1:nrow(occList[[i]])){
      lines(occList[[i]][j,],c(count,count),lwd=lineWidth,col=occColors[i])
      count<-count+1
    }
    count<-count+2
    }
  #return the occList as an invisible data object
  return(invisible(occList))
  }


dicelloData<-pbdb_occurrences(limit="all", base_name="Dicellograptus", show=c("phylo","ident"))
dicelloData<-dicelloData[dicelloData$rnk==3,]   #keep only occurrences of taxa resolved to species level

plotOccPBDB(dicelloData,"Dicellograptus Species")

Each color shown here is a different taxon, each line is the age uncertainty of a different occurrence. We can see for dicellograptids, the age uncertainty may potentially be a lot greater than the true sampled age ranges themselves. Many of these collections have age uncertainties that overlap and are not contiguous.

Now, the issue here is that for many statistical analyses that we apply in the fossil record, they use only the first and last appearance times, and they also often presume that these appearances are reported in contiguous, non-overlapping intervals of roughly equal length (let’s call that ‘well-behaved’ four-date data). Why? Well, most of these methods depend on the simple modeling trick that if we have contiguous intervals of roughly equal size, it becomes very easy to calculate the probability of observing a taxon X intervals into the future (or the past), given a sampling parameter. Obviously, there are considerable challenges in translating a dataset of loose occurrences to such well-behaved interval data, and that’s where my mind is right now.

Now, all that said, there are some analytical programs that take occurrence data natively. In particular, PyRate is a Python package that applies a Bayesian MCMC algorithm to estimate taxonomic origination, extinction and sampling rates from occurrence data. The method is described here and the program itself is described in MEE here. One could easily set up a functions to use [paleobioDB] package to download PBDB and use PyRate’s data preparations R functions to produce datasets for analysis in Python. Ultimately, PyRate solves the issue using a similar solution to the one I advocate via bin_timePaleoPhy (I discuss this more below): treating the age uncertainty of occurrences as uniform uncertainty bounds on the actual age of taxa and repeatedly randomly sampling ages from within those bounds and running the PyRate analysis. The sampling of ages is technically done in R; the Python component of PyRate only accepts single-date ‘point’ age data for occurrences as input. The results of multiple analyses can then be considered to see how much influence the uncertainty in taxon ages has on the resulting.

But PyRate doesn’t replace every pre-existing analysis of paleobiological data, and the ‘repeatedly sample and consider the resulting distributions’ trick isn’t an applicable solution to those analyses. So, how do we translate occurrence data into ranges within contiguous intervals? How does the manner in which we do this matter (or not) for the analyses we do? I’ve given this question a lot of thought previously with respect to time-scaling phylogenies of fossil taxa, but what about analyses where phylogenies aren’t involved? In particular, there are a number of functions in paleotree for estimating sampling rates and sampling probabilities that depend on having data from contiguous intervals (e.g., getSampProbDisc or make_durationFreqDisc or make_inverseSurv). These are questions I’ll try to address in future posts, along with a long-promised post about porting data from paleobioDB into paleotree format.

However, on a related node…

What does the Paleobiology Database report as two-date age ranges?

Some analyses only accept the two-date range data, and so people have invented algorithms for obtaining such datasets from occurrence data or four-date range data. However, such datasets must represent either precise age ranges or age uncertainty on a single occurrence, and thus information is being lost when we translate to two-dates (except in the rare cases where each taxon is known from only a single collection, and thus the two-dates are simply age uncertainty on singular finds). One can imagine some simple algorithms for obtaining two-date ranges, such as ‘maximum range’, i.e. take the earliest possible age of the oldest occurrence for a taxon and the latest possible age of the youngest occurrence. Although algorithms such as this retain some information about the geologic age of taxa, my opinion is that the information loss about age uncertainty may impact many analyses, perhaps adding an analytical bias. An alternative philosophy, which I espouse in the bin_timePaleoPhy function in paleotree, is to random draw precise FADs and LADs from four-date ranges (or occurrence data), run an analysis (such as time-scaling a phylogeny of fossil taxa), save the resulting test statistic or parameters of interest, and repeat this procedure many times to generate sampled distributions of the test statistics / parameters. Such a procedure retains (rather than loses) the uncertainty component of the ages into our results, but as I said above when discussing PyRate, this procedure isn’t applicable to every type of analysis.

Two of the various utilities or “apps” (…ugh) associated with the Paleobiology Database will return two-date age ranges for a search query. I became very curious how these ranges were calculated and so I wanted to report on my investigations.

FirstApp

Before I get into these various ways of returning ranges from occurrence data, I should mention a relavant not-range-returning application that deals with fossil taxon ages, FirstApp, a module by Matt Clapham. FirstApp takes paleobioDB data for a given taxon and displays occurrences through time as a density histogram, with occurrences placed as points at the earliest age boundary of the collections they occur within.

Ranges in Paleobiology Database Classic / Fossilworks

First, if you go to the present Paleobiology Database website, you have the option of using PBDB Classic which is a search query engine, with a download option (located here). You can access an almost identical interface at FossilWorks, which mirrors the paleobioDB data along with a number of analytical tools previously housed at the PBDB. FossilWork’s basic download option is located here. The algorithm for calculating two-date ranges used is somewhat complex but doesn’t appear to be described in any documentation, so I went straight to the “horse’s mouth” (his words), and asked John Alroy himself. I’ll paraphrase our conversation.

John has developed this “zone-of-overlap” algorithm to deal with obtaining ranges from the realities of messy occurrence data described above: age uncertainties (i.e. stratigraphic intervals) that overlap, of different lengths and don’t sort into an orderly hierarchy. Rather than returning the earliest possible date and the latest possible date like with the maximum-range algorithm described above, the zone-of-overlap algorithm “finds the oldest base that is older than at least part of all the intervals and the youngest that is younger than at least part of all the intervals” (direct quote, John Alroy).

Here’s three specific examples of how these ranges are calculated, taken from direct quotes from our email conversation. Two are from John and the third is from myself; as such, they are almost comedic examples of our taxonomic preferences:

• “…for example, if the two oldest occurrences are respectively Paleocene (not better resolved than that) and Thanetian then the base has to be the base of the Thanetian.”

• “Another example would be Hemphillian (Mio-Pliocene) and Pliocene. These work out to a range going from base of the Pliocene to the top of the Hemphillian even though the Hemphillian is ‘older’ (based on its base) and the Pliocene ‘younger.’”

• “If I follow correctly, some taxon listed in the PBDB with two occurrences, in a Katian collection (base ~456 Ma) and a Dicellograptus anceps zone (base is ~449 Ma) collection respectively, the lower age bound returned returned will be ~449 Ma.”

In other words, for calculating the first appearance datum (FAD), the algorithm looks for all occurrences that overlap with the age range of the earliest-most occurrence, obtains their earliest boundary ages and returns the latest-most earliest age boundary among these overlapping occurrences. Similarly, for calculating the latest appearance datum (LAD), the algorithm looks for all occurrences that overlap with the age range of the latest-most occurrence, obtains their latest boundary ages and returns the earliest-most latest age boundary among these overlapping occurrences. On theoretical grounds, one could probably describe the zone-of-overlap algorithm as minimizing taxonomic age ranges by assuming that overlapping occurrences probably describe a very similar FAD or LAD, and thus picks the one that extends the taxonomic range the least.

However, this does come with a downside that if these occurrences are not repeated attempts to capture the same FAD or LAD, then the zone-of-overlap algorithm isn’t an accurate depiction of the uncertainty in the ages; the age of the rocks containing the observed FADs and LADs may actually be outside of the reported range, unlike the maximum-range algorithm, which must contain the entire observed age range for a taxon as long as ourour interval dates for collections is accurate. Consider my graptolite zone example above; there is a possibility that the collection from the Katian is actually from an earlier graptolite zone than the D. anceps zone, maybe the Dicranograptus kirki zone (~453 Ma). Better stratigraphic resolution on a collection might very well provide such information in the future, but zone-of-overlap algorithm doesn’t give any weight to that possibility. On the other hand, John argues that a large number of taxa in the PBDB are listed from collections with very broadly defined intervals of time, and so some way is needed to downweight the effect of that lack of stratigraphic interval resolution. Given the nature of the data, the loss of uncertainty may be a necessary evil under the zone-of-overlap algorithm if you want to be as conservative as possible (say with providing molecular clock calibration dates for phylogenetic divergence dates for crown clades). I think the likely answer here is dependent strongly on the group we’re interested in and the reason we want range dates.

Ranges in R package paleobioDB

The paleobioDB package also has a function for obtaining and plotting two-date age ranges, the function pbdb_temp_range. Here’s a species-level example using the Dicellograptus data above, with the plotting functionality turned off.

pbdb_temp_range(dicelloData, rank="species", do.plot=FALSE)

##                              max   min
## Dicellograptus alector     471.8 443.4
## Dicellograptus anceps      455.8 443.7
## Dicellograptus ornatus     453.0 443.4
## Dicellograptus complanatus 453.0 443.7

So what are these ranges? The manual file for pbdb_temp_range just refers to them as temporal ranges or time spans of taxa.

Well, we can go to the Github repo for paleobioDB and read the code for pbdb_temp_range directly, here’s the relevant bit for species. It appears that the ranges reported by pbdb_temp_range are the maximum-ranges, the earliest and latest possible bounds on a taxon’s age range.

Wait, wait, hold on, why are there only four species listed…

So, tangentially, you’re probably wondering why so few species are reported, as we can see from the occurrences figure above that a much larger number appears to be contained within the downloaded occurrences table. Uh… I’ll address that in my next blog post.

---

PS: I just discovered how paleobioDB’s API allows you to call PhyloPic images for taxa, woah! Here, check out Nemagraptus:

http://paleobiodb.org/data1.1/taxa/thumb.png?id=281

Sweet!

Creating pretty plots of trees timescaled with paleotree via strap's geoscalePhylo

For anyone who cares, paleotree was not written for pretty plotting; other things were just more important to me. I don't have a particularly well-honed sense of aesthetics anyway; have you seen my recent Paleobiology paper, A.K.A. the grande parade of boxplots?

Over the years, I've gotten people asking me how they plot their fossil-taxon trees against, for example, a geological time-scale and I haven't been really able to tell these people what to do. But that changed recently, thanks to Mark Bell's and Graeme Lloyd's new package strap.

http://cran.r-project.org/web/packages/strap/

In addition to its own timescaling function and functions for assessing stratigraphic consistency, their package includes the extremely useful function geoscalePhylo which makes use of Mark Bell's geoscale package. This function takes a time-scaled tree and plots it against a geologic scale and, optionally, will display taxon ranges as thickened bars on the terminal branches.

Recently, I got one or two questions about using this function with trees time-scaled using paleotree, and so I decided to post an example using a dataset with a discrete interval range data; i.e. one of those datasets written for bin_timePaleoPhy and associated functions. Its mostly straightforward to do this, particularly as the 'bin_' time-scaling functions additionally outputs the taxon ranges used for a particular tree. we do need to do a minor step where we rename some columns in the range matrix for strap to accept these ranges.

So, here's a short R script that uses the example graptolite dataset from paleotree, as used in the cal3 tutorial from last year (i.e. Bates et al., 2005). See the old blog post for details about constructing the data files.

Now, you can just install strap and paleotree and give this script a whirl...
___________________________________________________________________

library(paleotree)
library(strap)

data(retiolitinae)

#need time-scaled tree
    #randomly draw FADs/LADs from within grapt biozones
timetree<-bin_timePaleoPhy(tree=retioTree,timeList=retioRanges,
    type="basic",add.term=FALSE)

#get the ranges used by the bin_ function
rangesUsed<-timetree$ranges.used
colnames(rangesUsed)<-c("FAD","LAD")

#now plot it
geoscalePhylo(tree=timetree,ages=rangesUsed,ranges=TRUE)

___________________________________________________________________________

Now, this is *much* prettier than anything paleotree makes, but we can do *much* better still.

First, this tree is a single random draw, an insignificant what-if produced by drawing the first and last appearance dates of taxa from the discrete intervals they are known to first and last occur in (here, Silurian graptolite biozones). For plotting purposes this is undesirable because you could remake this plot ten times and get ten slightly different looking plots. A better alternative, one matching plots in the literature, is to display taxa as ranging from the very start of the first interval they are found in to the very end of the last interval they occur in. This can be done easily with bin_timePaleoPhy's argument nonstoch.bin.

Secondly, the simple minimum-node-ages "basic" timescaling method almost always introduces artificial zero-length branches that, when plotted, make dichotomous nodes look confusingly like polytomies. Furthermore, this can obscure any polytomies your phylogeny contained to begin with.  So, for purposes of plotting, we can use one of the simple, arbitrary, post-hoc timescaling methods available in paleotree or strap. These methods rescale branches so there aren't any zero-length branches, but how realistic these various alternatives to "basic" are is another matter (you might want to read my most recent Paleobiology paper). Nevertheless, they are useful for showing relationships, as long as any one looking at it knows the apparent timing of divergences on the figure isn't even a best guess, just a random stab in the dark that's most probably dead wrong.

(Anyone who regularly reads paleontological papers will probably have that understanding.)

Finally, we can make the plot prettier just by playing with the graphical parameters of geoscalePhylo: make the text bigger, widen the branches, adjust the time-scale so we can see the Silurian diversification of the Retiolitinae with respect to the Ordovician and Devonian, etc.

And so...
____________________________________________________________________________

#now let's make it prettier 
    #non-random FADs and LADs
    #also use a 'prettier' time-scaling method

timetree<-bin_timePaleoPhy(tree=retioTree,timeList=retioRanges,
    type="mbl",vartime=1,nonstoch.bin=TRUE,add.term=FALSE)

rangesUsed<-timetree$ranges.used
colnames(rangesUsed)<-c("FAD","LAD")

geoscalePhylo(tree=timetree,ages=rangesUsed,ranges=TRUE,
    cex.tip=0.8,cex.ts=0.55,cex.age=0.5,x.lim=c(413,445),width=2)

_________________________________________________________________

Ta-da! Very pretty!

Enjoy!
-Dave

Raising A White Flag, paleotree v2.0 and the Reciprocal Monophyly of Various Hemichordate Groups

Hello all!

We have a lot to talk about today.

Last time, I bemoaned users who were misunderstanding my functions. Well, I have surrendered. timePaleoPhy is now happy to take taxon dates as min-max bounds. In the newest version of paleotree, version 2.0, you can find this functionality in the new argument dateTreatment for timePaleoPhy and other time-scaling functions.

In fact, a lot has changed in paleotree lately. Last night, I uploaded paleotree 2.0 to CRAN, and it should be making its way to your favorite repository as a binary of your choice soon. This, and the previous update, paleotree 1.9, have changed a lot about paleotree!

All these changes are kinda fitting, given that paleotree is now version 2.0, and the turning over of a new version number usually a rapid change in the structure of a package.  However don't really try to read too much into that, its really just a happy accident. There's a lot of very extreme opinions some people have about version numbering and I was a little naive when I started programming. So, way back in January of 2012, I just decided to go with a simple system where every public release except the most minor would be given a new increment of 0.1, and the first public release would be 1.0. I'm sure some people think that's a hideous way to do version numbering, but whatever. It ain't your package, dudes.

As always, you can look at paleotree's CHANGELOG to see what's new, but here's the last two entries for your reading pleasure.

__________________________________________________

Version - 1.9 - 01-02-14
-Changed all functions to individual R script files rather than single master script. Why didn't I do this years ago??!

-All help files converted to an roxygen2 format within function scripts

-At suggestion of Fabricio Villalobos, added option to expandTaxonTree that allows branch lengths to be retained and the added lower-taxa are connected with zero-length branches. Only useful for very specific purposes, please use with caution.

-Changed all lines which checked for class "phylo" of input objects to use is(obj,"phylo") instead of class(obj)=="phylo" per recommendation of Carl Boettiger
-Added new warning line to taxicDivDisc so people who try to pass it matrices with character strings in the matrix will get more helpful warning messages

-Upgraded probAnc and qSProb2comp based on equations in Foote (1996) for all three modes of origination, also fixed some errors in previous versions

-Added new function pqr2Ps which uses Emily King's exact derivation for the joint probability of a clade being (a) going extinct but sampled on an infinite time scale and (b) never going extinct on an infinite time scale

-Removed internal Ps function from cal3, now pqr2Ps which is exported directly to namespace

-Converted to new way to handle likelihood functions in paleotree, moving to a function-as-an-object system like diversitree.

-Following on last point, added new functions for fitting models of duration frequency data, replacing getSampProbDisc and getSampRateCont

-More models to follow in future versions of paleotree: added new function footeValues which will (eventually) support a release of a function that implements Foote's (2001, 2003, 2005) inverse survivorship models

-Modified diversitree's 'constrain' function to make a paleotree version named constrainParPaleo which is both entirely separate and fulfills needs such as constraining many similarly named parameters to a single value

-See new functions listed in ?modelMethods for manipulating functions in the new model format

-Added some not-exported hidden functions for use by the various model-handling functions

-Added a 'terrible idea' function optimPaleo which simplifies using optim with new parameter bounds functions. This function is entirely for pedagogical reasons and may be removed later.

-Added new function horizonSampRate which uses ML estimator from Solow and Smith (1997) for estimating sampling rate from precise durations in continuous time and number of sampled horizons

-Added new function perCapitaRates which estimates the per-capita origination and extinction rates from discrete interval data, following the methods from Foote (2000)


Version - 2.0 - 02-03-14
-Changed parInit to use uniform distribution to randomly draw initial parameter values between bounds, rather than take mid value between bounds

-Changed how time-scaling functions dealt with node.mins argument; can no longer use node.mins with a dataset that has unshared taxa that are to be dropped

-Altered example for use of node.mins in help files for time-scaling functions; thanks John Clarke of Oxford for the heads up! Also other modifications were made to the help descriptions, clarifying that node.mins can be used to constrain the minimum age for the root node.

-Also on a different issue brought to my attention by John Clarke, added an error message to paleotree when 'equal' is attempted by the edges leading to the root are zero-length (because 'equal' cannot run under this situation!)

-Added new function collapseNodes that collapses specific user-defined nodes, either forward or backward

-Made all lines checking for dichotomous trees check both with is.binary.tree() *and* is.rooted()

-Added new function dateNodes which returns the dates of the internal and tip nodes of a phylogeny on an absolute time-scale, with respect to the $root.time element if one exists

-On a trial basis, I have added new function inverseSurv which attempts to replicate the forward and inverse survivorship modeling applied by Foote (2001, 2003, 2005) and is useable with the newly implemented constrainParPaleo framework implemented in the previous version of paleotree. I am not yet convinced this function is a 100% faithful replicate of the original method.
-fixed error in plotTraitgram where if trait data was entered in same order as tree$tip.label, trait data was not resorted prior to running ace

-'equal' method in timePaleoPhy wasn't returning same result as 'equal' method in DatePhylo from Graeme Lloyd's original code. This turned out to be a result of differences in how we ordered nodes: Graeme ordered them by time or distance-from-the-root (using dist.nodes) and I was using node.depth, which counts number of branching events. This choice shouldn't make a considerable difference on the performance of the algorithm, but does produce some differences in the resulting time-scaled trees. For consistency, I have change timePaleoPhy to match Graeme's algorithm.

-altered timePaleoPhy and cal3timePaleoPhy to allow the point date occurrences with the first and second column of timeData interpreted as bounds instead, using the argument dateTreatment="minMax"

-related to above change, the argument rand.obs was removed from timePaleoPhy and cal3TimePaleoPhy as no longer necessary, this functionality is now available via dateTreatment="randObs"

-Although it may be strange to not report a lack of a change, but still have not added the finite time window approach for durationFreq

________________________________________________________

Hahah, and as you might guess from that last one, I've still got some more new things and changes for paleotree in store in the next few months. I included it here because I was actually partway through adding this option and had to undo those additions through commenting, since I wanted to push this version to CRAM (the undoing worked, I think, but I can't be certain, so best to add it to the CHANGELOG). 

As always, let me know what you think of the new paleotree functionality!

So what else?

Well, hemichordate phylogeny has been shaken back and forth a little lately. You might have missed them, so here's a short list:

http://onlinelibrary.wiley.com/doi/10.1002/jez.b.22510/full
http://link.springer.com/article/10.1007/s00114-013-1117-3
http://www.biolbull.org/content/225/3/194.abstract

The Stach article in particular is covered by Cambrian Mammal's blog:
http://cambrianmammal.wordpress.com/2014/01/27/the-use-of-a-larva

The Cannon et al. is really neat: with greater gene and taxon coverage than a few years ago, it looks like the different hemichordate groups really are reciprocally monophyletic and not nested, and maybe even the Rhabdopleura and Cepholdiscus groups are even reciprocally monophyletic. Big implications for what the stem deuterostome looked like... and even bigger implications for the stem graptolite.... if you're into that sort of thing. ;)

Oh, and finally, I'm also a brachiopod worker now! I've begun a remote post-doc with Sandy Carlson at UC Davis. I'm looking forward to doing some neat things with her and the rest of her lab!

-Dave

Stratigraphic Uncertainty, timePaleoPhy versus bin_timePaleoPhy, and My Own Paleontological Myopia

Hi all!

So, its been a rough week! I've gotten more emails about issues or questions about paleotree this week than I have in the past four months! Which is awesome, cause I love talking to people about using paleotree, but there's something that's been worrying me in some of these emails. These emails have given me a grave concern that a number of users don't understand the distinction between two key time-scaling functions in paleotree: timePaleoPhy and bin_timePaleoPhy. The problem is, not understanding this difference is extremely dangerous, as they could generate time-scaled trees that don't make sense without realizing it.

But I can't really blame them, because I think its my fault as the package author to explain this distinction clearly enough... the problem is I didn't understand something critical and fundamental three years ago: Our Fossil Records are Different.

See, the fossil record for your group is probably really different from the fossil record I'm working on, and that means not only how we do things is different, but the sort of assumptions we come to the table with and what we think X means is different. And that sounds kind of obvious, until you meet someone whose fossil record is radically different from yours and they don't even have a concept for this one thing you thought everybody knew about.

The problem is that people are using timePaleoPhy who probably shouldn't be using timePaleoPhy, when they should be using bin_timePaleoPhy instead. It isn't helped by the argument 'rand.obs' in timePaleoPhy, which some interpreted as representing how to bring stratigraphic uncertainty into an analysis. After explaining to several users over email why I was concerned about their use of paleotree, I felt like a specific blog post on the topic was needed.

But to understand where this misconception of some users comes from, and why I think it stems from my own misconceptions in the past that maybe fossil records were more similar than I know now, I've also got to explain my reasons for writing both functions the way I did.

So, I wrote timePaleoPhy first, when I had some stratigraphic data for graptolite ranges, with stratigraphic meters scaled to an absolute time-scale. Since graptolite species often persist over time, I had a matrix of first and last appearance dates (FADs and LADs, respectively), like so...

Normalograptus_normalis 457.9 438.3
Climacograptus_typicalis  463.2 442.8

(This is totally fake data; I don't even know if typicalis is still a Climacograptus, but maybe its the type?...I can't remember at this very moment.)

The dates here represent actual ranges; i.e. N. normalis is first found in rocks dated 457.9 Ma and is found fairly continuously until you get to rocks dated 438.3 Ma. Of course, the rocks aren't actually 'dated' so precisely, more like we're assuming continuous sedimentation rates and so we're inferring those dates using the absolute dates we do have.

This means when we time-scale the tree, unless we're allowing for ancestors (which we don't normally do unless we know something about apoomorphic and plesiomorphic taxa) the clade containing N. normalis has to be at least as old as 457.9 Ma.

And that's what timePaleoPhy does, if you hand it data like above, you get back a tree where the nodes are as old as the dates in the first column (the taxon names should be row labels), because it thinks those are first appearance dates, known very precisely.

A different question is where the tips go (or 'when', really). As I've talked about before here on this blog and in published papers, this is a difficult question. For most things, such as looking at a tree or estimating lineage diversity through time, you'd want the whole taxon range 'added' to the tree, so you should use the LAD as the tip-date. However, comparative methods for studying trait evolution assume that the date of tips represents the date of the population from which trait values were taken. Of course, maybe your trait doesn't vary much (or at all) over the duration of the taxon, but that's a story for another time. In the papers I've written, I often call these tip-dates 'times of observation'... cause I think some reviewer didn't like 'tip-date', or maybe it was a committee member. Well, reviewers seem to dislike 'times of observation' too.

(I can't freakin' win when it comes to vocabulary.)

At the time, I was most interested in comparative methods, so I assumed this would be true of everyone else. In general, it seemed like people were using the FAD primarily in the literature (or at least, I thought so at the time) as the tip-date for comparative methods, so I made the FAD the default tip-dates that timePaleoPhy spits out. To get timePaleoPhy to use the LAD as the tip-date, you have to set add.term=TRUE, which adds the range to the terminal branch the taxon sits on.

But using the FAD or the LAD wasn't really a satisfying solution: the tip-date in a comparative analysis supposed to be the time of the sampled population (hence why I call them 'times of observation'). For my long-ago purpose, I had I had size measurements from a bunch of figured specimens for graptolites but with very little data on 'when' those specimens came from. They were from 'the Late Ordovician' or 'Early Silurian' and I didn't have enough information to figure out a more specific tip-date. So that's pretty much useless information for figuring out the 'when' for these specimens, and if I want to do some PCM, well the tips are supposed to represent the date of the population where I measured the trait. I just know the Normalograptus normalis specimen I measured has a 4.2 mm long doodad. Using the FAD or the LAD wasn't really an honest approach given the uncertainty involved.

So... rand.obs is how I dealt with this uncertainty, still using the FAD for N. normalis to adjust the clade age, but allowing the tips to fall at some date randomly selected between the FAD and the LAD. Maybe I had measured doodad length from a normalis specimen from 441.3 Ma,  so, the tip would be at that date, but the clade constraining that taxon would still be constrained by the 457.9 Ma first appearance date. Note that rand.obs only works if we also set add.term=TRUE, cause the way I decided to do it (and no, I'm not changing it) was that these random observation dates were a temporary replacement of the LADs. Thankfully, timePaleoPhy returns an error if you ever try using rand.obs=TRUE while add.term=FALSE, or else people might get very confused.

Later, I needed to make use of a much less precise dataset, where taxon first and last occurrences were graptolite zones. So like, the different zones had minimum and maximum dates associated with their starting and ending, like this...

zone_A 466.2 452.4
zone_B 452.4 448.5
zone_C 448.5 441.1
zone_D 441.1 436.2
zone_E 436.2 431.0

(This is totally fake cause graptolite zones are usually shorter than that.)

To keep track of everything, I needed some way of recording the dates for the intervals and which intervals taxa first and last showed up in. Normally you see this thing displaying in papers as a range chart, where species A will have Xs in intervals 1, 2, 3, and then not 4, meaning it last appeared in interval 4. This seemed like a fairly non-straightforward way of handling the information, particularly as I was only really interested in the first and last interval each taxon appeared in.

So, what I did was I listed the first and last interval that each taxon was known from as the row of the interval matrix. Like so:

N._normalis 1 4
C._typicalis 1 3

So, if we look at this in reference to the prior matrix, this tells us that N. normalis first appears in zone A and last appears in zone D, which means it might have first shown up anywhere between 466.2 and 452.4 Ma and last appeared somewhere between 441.1 and 436.2 Ma. And so now, if we have both this 'taxon times' matrix and the 'interval times' matrix above, we have all the information on the dates associated with when these taxa were first and last occurring.

In paleotree, I decided to call a list that contained both matrices a 'timeList' object because as you have no doubt deduced, I have terrible taste when it comes to coming up with new vocabulary. A lot of functions in paleotree use these timeList objects. (It's not a real object class, cause I didn't understand S3 classes yet when I wrote paleotree.) Anyway, I went with a list of two matrices rather than two seperate matrices because the taxonTimes matrix doesn't make sense without reference to the intervalTimes matrix and so I felt I needed to keep them sorted together. To get a timeList, you can just make the two matrices by hand in your favorite spreadsheet program, read them into R and then combine as a list. The last step is as simple as typing time<-list(intervalTimes,taxonTimes) into the terminal. Note that the interval names and taxon names are rownames, so the matrices should both have only two columns each.

Now, to time-scale a tree with this sort of data, we need to somehow get precise first and last appearance dates and then make a tree from there. So I made bin_timePaleoPhy to draw dates from these intervals, under a uniform probability distribution, and then use the randomly drawn dates to make a time-scaled tree. So, maybe normalis was first seen  at 454.2 or 461.8 Ma, and last seen at 440.7 Ma, or 438.8 Ma. Or lots of other combinations! We can resample those dates and get all sorts of new FADs. and LADs and make lots of time-scaled trees. Unlike timePaleoPhy, bin_timePaleoPhy actually does something about stratigraphic uncertainty, at least as far as having data in discrete intervals goes.

But this doesn't fix the unknown times of observation issue for trait data... so bin_timePaleoPhy also has a rand.obs argument just like , so that we can let our observed doodad length for N. normalis come from 452.6 Ma or 446.3 Ma or whatever, depending on the set of FADs and LADs we plucked from the intervals in a given run. Again, just like with timePaleoPhy, the date of the nodes are still going to depend on the FADs, rand.obs is only going to impact the tip-dates, those enigmatic 'times of observation'.

So, let's be clear:

1) timePaleoPhy does *not* deal in stratigraphic uncertainty of appearance dates: it believes the numbers you are handing it are very precise first and last appearance dates, which exist for some fossil records but probably not most. It does deal with uncertainty in the times of observation of specimen measurement, though, which is what rand.obs is for.

2) bin_timePaleoPhy *does* deal with stratigraphic uncertainty, but it requires slightly more complicated input: two matrices, as a list.

And I think, somehow, I said something that is confusing people about these two facts. I'm repeatedly getting emails from people that think the matrix for timePaleoPhy is a matrix of the minimum and maximum ages for a point occurrence of a fossil. In other words, their fossil taxa only appear at a single point in time which is imprecisely known, and so their dates are the bounds on that single date, rather than two very precisely known dates that represent the first and last appearances.

This practice is very worrying, because if they put those dates in and then set rand.obs=TRUE and add.term=TRUE, they'll think their pulling their taxon occurrence dates randomly from their min and max dates... which is true for the tip-dates. But the nodes are being pushed back in time to those minimum dates, which really means that clades are 'as old as they could possibly be' or something wacky like that. And I bet people have done this without realizing that's what is happening! Its absolute nonsense and this nonsense has to stop.

I think the real issues is that when I wrote these functions, I was excessively myopic. I looked at things like a graptolite worker and I thought everyone had data like graptolite data, with taxa that typically first and last appeared at different times. I actually didn't realize that other people had data where it was more common to see things like 'species A is known from one skeleton, found somewhere in the early Cretaceous'. Some vert paleo friends of mine have kindly helped me understand that this isn't true, at least, for lots of vert datasets: they actually have taxa known from single collections that may be as poorly unconstrained as 'anywhere within a 40 million year interval', which, well, that blows my mind every single time I read this sentence.

So, given the apparent and continual misuse of timePaleoPhy, I have actually considered pulling the function entirely from the CRAN version of paleotree, or hiding the function, making it unavailable to regular users. This isn't ideal: some people, such as those who work with detailed biostratigraphic records such as for index fossils and microfossils, actually have precisely dated first and last appearances. I don't know if anyone like that uses paleotree yet, but I want to keep the door open for them! Plus, bin_timePaleoPhy is really just a dumb wrapper for timePaleoPhy. If someone really wants to write their own wrapper for timePaleoPhy, that's fine! I'm happy for them.

But honestly, the majority of people won't have data that's infinitely precise. They should be using bin_timePaleoPhy. (Or you know, bin_cal3TimePaleoPhy, but maybe I shouldn't get ahead of myself...)

So, the real question is how to convince people to want to use bin_timePaleoPhy. Well, although I just explained the function from a grapt-point-of-view, I didn't code it for that only that type of data. The timeList objects takes by bin_timePaleoPhy and other paleotree functions can actually be much less formal than my example (well, a few specific functions have very specific needs, but that isn't true of bin_timePaleoPhy).

For one thing, the intervals in the interval matrix can be completely messy. For example, you could have an interval matrix that looks like:
Lower Cumberbatchian 134.5 129.4
Cumberbatchian-Martinian 134.5 109.6
Cumberbatchian-Bloomean 134.5 99.8
Middle Upper Martinian 107.5 105.8
Martinian-Bloomean  115.4 99.8

What I've discovered in a recent collaborative project is if you get data from the PBDB, you'll probably have to have an interval matrix like this. It's okay. But if you have multiple collections for a taxa, try to get the most precise first and last intervals.

For example, say your taxa only ever appear in single interval, well than the first and last intervals are identical, so you'd just structure the second part of the input list so the columns matched, like so:

Rareosaurus 1 1
Uniqueodon 2 2

This is a little more tricky if your taxa also only ever appear in a single collection. By default, bin_timePaleoPhy assumes taxa first and last appear at different dates within an interval. For example, Rareosaurus will not be assigned matching FADs and LADs. However, if you set the argument point.occur=TRUE, then bin_timePaleoPhy will treat all your taxa as if their first and last occurrences are identical: i.e. all your taxa will be treated as point occurrences.

Of course, I have made it so that bin_timePaleoPhy has a little more functionality, to deal with those special cases that often crop up with uncertainty in appearance times; these can be usually handled with what I called the 'sites' matrix. For example, perhaps the first specimens of two taxa appear in the same fossil assemblage. That single fossil assemblage may be very poor resolved to a wide interval, but we know the first appearance date for both taxa should be the same date. If you make a site matrix as I describe in the bin_timePaleoPhy help file, and set the site number the same for the first appearance of both taxa, then bin_timePaleoPhy will always use the same randomly-drawn date for those taxa. You could do this for lots of taxa, or do this to account for a site where some taxa have their last appearance at and other taxa first appear at, or all sorts of neat modifications. I have to be honest and point out that this sites matrix idea was first suggested by Jon Mitchell over lunch about three years ago.

In reality, the point.occur argument I described is just a simplified way of modifying the site matrix so all first and last intervals are coded as having the same 'site' for each taxon. If you want to have only some taxa in a dataset constrained to be point occurrances, you can just specify a custom-made site matrix where only those taxa have their first and last interval constrained to be the same.

To be truly honest though, I don't think its a lack of advertising about the options offered by bin_timePaleoPhy that stops people from using it... really, its the timeList input format. A large number of people have suggested that this is a very obscure way of recording this data and more than a few have complained that its just not amendable to their use. Peter Smits told me that the timeList was specifically very 'ugly' to have to use two matrices. Overall, I have a feeling that I just did a shoddy job of describing timeList objects in the help files to begin with.

Well, I don't know what to do about it guys! Coming up with the two-matrix fix was the best idea I had after thinking about it for a while, and while I admit its inelegant and while I admit that the timeList structure is pretty dissimilar from the typical data structure you might see in a paleontological dataset (as opposed to a taxon by interval range chart), I just don't see a better solution, even though I can that its actively driving people to misuse the functions in paleotree in ways that I think are unpredictable and a little dangerous.

But I know that I don't know everything... so if anyone has a better idea of how to make the bin_timePaleoPhy input work better, I'd love to hear it and implement your idea in paleotree! Do you see a better, simpler way of conveying the same information? Or, is there a data format that you use, that you'd like to have a function for converting into a timeList? Let me know in the comments below.

For those of you interested, I give a real example of a very nicely behaved timeList object in my time-scaling method tutorial from this summer, which you can find here:
http://nemagraptus.blogspot.com/2013/06/a-tutorial-to-cal3-time-scaling-using.html

On a more philosophical note, I think those of us who write methods for paleontological data have to keep in mind that Fossil Records are Different. Had I been less nearsighted and graptolite-oriented in my design and help file writing in the beginning, maybe this user preference for the wrong function could have been avoided, and so for that, I really only can blame myself.

Finally, if you're reading this, and you think you might be one of those people who sent me an email this week and contributed to my concerns, I just want to thank you for emailing me and forcing me to confront the fact that there was something about my software that users were finding consistently confusing. I still love reading and responding to every paleotree email I get.

-Dave