Purpose

To remove sampling from a diversity curve.

Data required

Time series of both diversity (taxonomic richness) and a sampling proxy (e.g. formations, map area) as well as dates for each time bin.

Description

Although reconstructing estimates of past diversity is a central goal of palaeobiology raw data can be misleading due to sampling and other biases. Raup (1972) made an early tabulation of such biases and set out two major ways in which they might be corrected. The first was subsampling, which has been explored thoroughly by John Alroy (Alroy et al. 2001, 2008; Alroy 2010a,b,c). The second was modelling, which has received comparatively little attention. A key paper, then, was that of Smith and McGowan (2007) who introduced a new modelling technique that assumes true diversity is constant such that observed diversity is driven purely by sampling. When this modelled estimate is subtracted from actual observed diversity then the residuals can be thought to represent a sampling-free diversity curve. Barrett et al. (2009) added a rudimentary significance test to the Smith and McGowan (2007) approach by treating the standard deviation of the residuals as a kind of confidence interval. However, this is flawed as it does not take into account how well or poorly the model fits. Most recently, Lloyd (2012) overhauled the Smith and McGowan method, by allowing non-linear relationships between the sampling proxy and diversity, implememting a more informative significance test and introducing a hinged-regression line to the resulting time series to elucidate medium-term diversity trends. This refined method has been applied to dinosaurs Lloyd (2012) and coccolithophores (Lloyd et al. 2011). Here a function that automates the Lloyd (2012) approach has been produced for use in the freely available statistical progamming language R.

Implementation

1. Install and run R (it's free!).

2. Install the earth, nlme, paleoTS and plotrix libraries from within R, making sure to include any dependencies.

3. Copy and paste the following line into R and hit enter:

    source("http://www.graemetlloyd.com/pubdata/functions_2.r")

4. Enter your diversity data in order and store it as a variable named "div" (here, as an example, I use the sauropodomorph data of Barrett et al. 2009):

    div<-c(0,0,0,0,0,1,1,1,2,9,
    9,9,5,5,5,11,11,11,10,10,
    10,12,12,12,8,9,8,5,5,5,
    7,7,8,12,12,15,15,15,15,9,
    9,9,37,37,38,31,30,30,3,3,
    3,5,5,5,4,4,5,7,6,6,
    7,9,10,14,12,13,8,6,7,5,
    2,2,3,3,3,4,4,4,9,11,
    10,13,13,14)

5. Enter your sampling proxy data in the same way and store it in a variable named "proxy" (here again, I use the dinosaur-bearing formation data of Barrett et al. 2009):

    proxy<-c(1,1,1,2,2,3,11,9,10,15,
    13,15,7,7,7,44,39,40,36,37,
    39,27,27,27,23,23,22,15,14,15,
    19,19,22,32,32,32,26,26,26,46,
    54,61,60,54,66,65,61,70,85,83,
    86,88,85,85,93,93,98,90,82,91,
    103,98,113,119,113,120,113,101,102,80,
    80,80,85,84,87,85,81,91,134,115,
    150,159,133,175)

6. Finally, enter a time value (in Ma) for each bin and store it in a variable named "time" (here again, I use the time data of Barrett et al. 2009):

    time<-c(243.67,241.00,238.33,235.50,232.50,229.50,226.08,222.25,218.42,214.35,
    210.05,205.75,202.93,201.60,200.27,199.08,198.05,197.02,195.35,193.05,
    190.75,188.50,186.30,184.10,181.77,179.30,176.83,174.93,173.60,172.27,
    170.95,169.65,168.35,167.20,166.20,165.20,164.12,162.95,161.78,160.28,
    158.45,156.62,154.88,153.25,151.62,149.92,148.15,146.38,144.62,142.85,
    141.08,139.57,138.30,137.03,135.33,133.20,131.07,129.17,127.50,125.83,
    122.83,118.50,114.17,109.93,105.80,101.67,98.58,96.55,94.52,92.80,
    91.40,90.00,88.72,87.55,86.38,85.42,84.65,83.88,81.35,77.05,
    72.75,69.75,68.05,66.35)

7. We can now run the function, storing the results in a new variable called "results":

    results<-rockmodel.predictCI(proxy,div)

8. The function returns five time series that can be accessed using '$' then the series name. So for example, to get the predicted value we can write:

    results$predicted

This is the model-predicted value based on the assumption true diversity is constant and observed diversity is driven purely by the sampling proxy. The remaining four series are the upper and lower 95% confidence intervals based on standard error:

    results$selowerCI
    results$seupperCI

and standard deviation (recommended):

    results$sdlowerCI
    results$sdupperCI

There is also information on the favoured model (y being diversity and x the sampling/rock proxy):

    results$model

9. However, it is best to plot the data to visualise it by copying and pasting the following:

    plot(time,results$predicted,type="l",xlim=c(max(time),min(time)),
        ylim=c(0,max(results$sdupperCI)),col="grey",
        xlab="Time (Ma)",ylab="Taxonomic richness")
    points(time,results$selowerCI,type="l",lty=2,col="grey")
    points(time,results$seupperCI,type="l",lty=2,col="grey")
    points(time,results$sdlowerCI,type="l",lty=4,col="grey")
    points(time,results$sdupperCI,type="l",lty=4,col="grey")
    points(time,div,type="l")

This creates a graph that shows the modelled estimate of diversity in grey (first line), adds a dashed grey line for the standard error 95% confidence interval (second and third lines), adds a dash-dot grey line for the standard deviation 95% confidence interval (fourth and fifth lines) and plots the actual data in black (last line).

10. A better way to visualise things, however, is to look at the model-detrended richness. This can be done by copying and pasting the following:

    plot(time,div-results$predicted,type="l",xlim=c(max(time),min(time)),
        xlab="Time (Ma)",ylab="Model Detrended Taxonomic Richness")
    polygon(x=c(time,rev(time)),y=c(div-results$predicted,
        rep(0,length(time))),col="grey",border=0)
    points(time,div-results$predicted,type="l")
    lines(time,rep(0,length(time)))
    lines(time,results$seupperCI-results$predicted,lty=2)
    lines(time,results$selowerCI-results$predicted,lty=2)
    lines(time,results$sdupperCI-results$predicted,lty=4)
    lines(time,results$sdlowerCI-results$predicted,lty=4)
    marsout<-mars(time,div-results$predicted)
    if(marsout$best.model != "Linear") {
      points(time,marsout$prediction,type="l")
      points(marsout$cuts,marsout$prediction[match(marsout$cuts,time)],pch=19,cex=0.5)
    }

This plots the model-detrended diversity as a grey polygon (first to fourth lines), adds a dashed line for the standard error 95% confidence interval (fifth and sixth lines), adds a dash-dot line for the standard deviation 95% confidence interval (seventh and eighth lines) and fits a hinged-regression (if optimal) to the model-detrended tme series (ninth to thirteenth lines).

References cited

Alroy, J. 2010a. The shifting balance of diversity among major marine animal groups. Science, 329, 1191-1194.

Alroy, J. 2010b. Fair sampling of taxonomic richness and unbiased estimation of origination and extinction rates. In: Alroy J, Hunt G (eds.) Quantitative Methods in Paleobiology. The Paleontological Society; New Haven, USA. 55-80.

Alroy, J. 2010c. Geographical, environmental and intrinsic biotic controls on Phanerozoic marine diversification. Palaeontology, 53, 1211-1235.

Alroy, J., Marshall, C. R., Bambach, R. K., Bezusko, K., Foote, M., Fursich, F. T., Hansen, T. A., Holland, S. M., Ivany, L. C., Jablonski, D., Jacobs, D. K., Jones, D. C., Kosnik, M. A., Lidgard, S., Low, S., Miller, A. I., Novack-Gottshall, P. M., Olszewski, T. D., Patzowsky, M. E., Raup, D. M., Roy, K., Sepkoski, J. J., Sommers, M. G., Wagner, P. J. and Webber, A. 2001. Effects of sampling standardization on estimates of Phanerozoic marine diversification. Proceedings of the National Academy of Sciences of the United States of America, 98, 6261-6266.

Alroy, J., Aberhan, M., Bottjer, D. J., Foote, M., Fursich, F. T., Harries, P. J., Hendy, A. J. W., Holland, S. M., Ivany, L. C., Kiessling, W., Kosnik, M. A., Marshall, C. R., McGowan, A. J., Miller, A. I., Olszewski, T. D., Patzowsky, M. E., Peters, S. E., Villier, L., Wagner, P. J., Bonuso, N., Borkow, P. S., Brenneis, B., Clapham, M. E., Fall, L. M., Ferguson, C. A., Hanson, V. L., Krug, A. Z., Layou, K. M., Leckey, E. H., Nurnberg, S., Powers, C. M., Sessa, J. A., Simpson, C., Tomasovych, A. and Visaggi, C. C. 2008. Phanerozoic trends in the global diversity of marine invertebrates. Science, 321, 97-100.

Barrett, P. M., McGowan, A. J. and Page, V. 2009. Dinosaur diversity and the rock record. Proceedings of the Royal Society, London B, 276, 2667-2674.

Lloyd, G. T. 2012. A refined modelling approach to assess the influence of sampling on palaeobiodiversity curves: new support for declining Cretaceous dinosaur richness. Biology Letters, 8, 123-126.

Lloyd, G. T., Smith, A. B. and Young, J. R. 2011. Quantifying the deep sea rock and fossil record bias using coccolithophores. In: Comparing the Rock and Fossil Records: Implications for Biodiversity, McGowan, A. J. and Smith, A. B. (eds.). Geological Society Special Publication 358, 167-177.

Raup, D. M. 1972. Taxonomic diversity during the Phanerozoic. Science, 177, 1065-1071.

Smith, A. B. and McGowan, A. J. 2007. The shape of the Phanerozoic diversity curve. How much can be predicted from the sedimentary rock record of Western Europe? Palaeontology, 50, 765-777.