A flexible Bayesian method for estimating stratigraphic intervals and their co-occurrence in time

A flexible Bayesian method for estimating stratigraphic intervals and their co-occurrence in time | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results A flexible Bayesian method for estimating stratigraphic intervals and their co-occurrence in time View ORCID Profile Gustavo A. Ballen doi: https://doi.org/10.1101/2025.02.13.638199 Gustavo A. Ballen 1 Instituto de Biociências de Botucatu, Universidade Estadual Paulista “Júlio de Mesquita Filho” , Botucatu, SP, Brazil 2 School of Biological and Behavioural Sciences, Queen Mary University of London , London, United Kingdom 3 Center for Tropical Paleoecology and Archaeology, Smithsonian Tropical Research Institute , Ancón, Balboa, Panamá 4 Museu de Zoologia da Universidade de São Paulo, Universidade de São Paulo , São Paulo, SP, Brazil Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Gustavo A. Ballen For correspondence: gustavo.a.ballen{at}gmail.com gustavo.ballen{at}unesp.br Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract The fossil record is both incomplete and full of biases, which hinders making inferences about the past of life on earth. Stratigraphic intervals help us estimate two important quantities from the fossil record of a given lineage, its origination, and extinction times. Many models are available in the literature, although specialising in implementation and therefore resulting in limitations of application. This toolset is extended herein by proposing a flexible method which allows one to estimate origination and extinction times, as well as model preservation potential in flexible ways. The present method is more general than previous alternatives, many of which constitute special cases. A method is presented to represent the stratigraphic interval as a continuous distribution by integrating over parameter uncertainty through the use of the posterior predictive distribution. It is possible then to combine these distributions using conflation, to build a time probabilistic model for the co-occurrence of these intervals. Empirical examples using the method for estimating the origination and extinction times of Palynomorphs and the co-occurrence time for an assemblage of unknown age, as well as the origination time of the marine Barracudas, are used for illustrating the potential of the method. It is implemented in the StratIntervals.jl Julia package, which allows a large set of possible prior distributions and MCMC samplers. 1 Introduction There is a powerful, yet blurry window into the past of life on earth, which is the fossil record. Many insights into how the life came to be in its present form, and caveats on the possible directions that we can go regarding climate change and the biodiversity crisis are a few of the many applications of the knowledge we can gather from studying fossils ( Smith et al., 2025 ). It has long been recognised that the fossil record is very incomplete, and multiple sources of bias hinder our ability to use it for inferences about the environments in deep time ( Signor and Lipps, 1982 ). However, statistical tools have proved to be effective in helping us deal with the imperfections of this record, and gather insights from it through quantitative palaeobiology and its integration with evolutionary biology ( Holland et al., 2024 ). One of the most difficult questions to ask about any extinct lineage are when they appeared, and then when they disappeared, that is, to estimate origination and extinction times. Stratigraphic intervals have been proposed to reach this goal using information from the set time occurrences of a given lineage preserved in the fossil record, and statistical methods which describe the relationship between the pattern of such occurrences and the endpoints (i.e., origination and extinction times) of the interval ( Strauss and Sadler, 1989 ). Although these methods were first developed with palaeobiological motivations, they have also been applied and extended in the realm of conservation biology by using collections of sightings in time as occurrences ( Rivadeneira et al., 2009 ). In general, any system where occurrences in time are recorded, and which has start and end times could be studied using stratigraphic intervals (e.g., archaeology and cultural evolution). Furthermore, they are still useful and applicable when only one of the times is of interest, for instance, the origination time but not the extinction time, or vice versa. There are several methods available for inference on stratigraphic intervals, with a particular focus on confidence intervals on endpoint parameters. These methods vary in assumptions (e.g., whether preservation potential is uniform), statistical paradigm (e.g., whether frequentist or Bayesian), and motivation for development (e.g., applied to palaeobiology or conservation biology). Comprehensive reviews on this subject are available in the literature ( Marshall, 1990 ; Solow, 2005 ; Rivadeneira et al., 2009 ; Marshall, 2010 ; Boakes et al., 2015 ; McCrea et al., 2024 ). Rivadeneira et al. (2009) arranged the methods known at that time into three classes: Class 1 methods, which assume constant preservation potential (e.g., Strauss and Sadler, 1989 ; Solow, 1993a ; Ferraz et al., 2003 ; McInerny et al., 2006 ), Class 2 methods, which relax that assumption (e.g., Solow, 1993b ; Marshall, 1997 ; McCarthy, 1998 ; Ferraz et al., 2003 ), and Class 3 methods, which are distribution-free (e.g., Marshall, 1994 ; Solow and Roberts, 2003 ). Later methods, especially that of Wang et al. (2016) and the one presented here fall in Class 2 because preservation is not assumed to be constant but still rely on distributions. Bayesian approaches have been proposed since the first attempts at estimation on stratigraphic intervals. Strauss and Sadler (1989) already recognised the advantages of Bayesian formulations of their method and provided a specific one using a uniform prior and a likelihood function constrained to the interval [0, 1]. Weiss and Marshall (1999) assume a constant preservation potential and were the first to consider discrete sampling and therefore a non-continuous likelihood function. Later, Weiss et al. (2004) extended this framework for dealing with non-constant preservation by using abundance data. Both methods are different from any other Bayesian one in that the likelihood function requires discrete sampling. Ferraz et al. (2003) provided two other Bayesian formulations, one assuming linearly decreasing preservation potential, and one uniform, in which suddenly extinction occurs. They favoured the former model, claiming that it is more realistic. Wang et al. (2016) proposed a flexible method in which different types of preservation potentials can be fitted, although they used fixed priors for the free endpoint and the preservation potential parameter. Then they used analytical formulae to calculate exactly the posterior distributions of the model parameters. However, the method conditions on one of the endpoints being zero. Alroy and Solow discussed in detail several aspects of Bayesian methods starting with a pair of papers by Alroy (2014 , 2015 ) applied to the estimation of extinction, including the specification of priors and justification for likelihood functions ( Alroy, 2016a , b ; Solow, 2016a , b ). Finally, Kodikara et al. (2020) proposed a method for estimating the year of extinction using sighting data and implemented the hierarchical Bayesian model to sample the posterior parameter distributions. Wang (2005) examined the correlation between parameters in Beta models and concluded that endpoint parameters tend to be correlated, which is an issue for inference using Metropolis-Hasting samplers. It is clear from the diversity of approaches using Bayesian techniques that both flexible models and their implementation in computational tools are necessary to explore and extend these Bayesian methods. The goal of the present paper is to provide a general and flexible implementation of stratigraphic intervals which allows the estimation of both origination and extinction, as well as their co-occurrence in time. 2 Materials and methods 2.1 Implementation The Bayesian implementation of the three-parameter Beta distribution as a model for stratigraphic intervals has been implemented in the StratIntervals.jl package, written in pure Julia ( Bezanson et al., 2017 ). The package uses several dependencies ( DataFrames , Distributions , Interpolations , KernelDensity , QuadGK , Random , SpecialFunctions , and StatsPlots ; Bouchet-Valat and Kamiński, 2023 ; Besaņcon et al., 2021 ; Johnson, 2013 ) which provide data structures, plotting engines, or implement specific mathematical and statistical procedures not available in the Base Julia library. In particular, StratIntervals uses Turing ( Ge et al., 2018 ), a package providing a probabilistic modelling language which allows for flexible implementation of MCMC algorithms. The stable version of StratIntervals can be found in the official Julia Registry and installed with: # in the Julia command line using Pkg # install StratIntervals Pkg.add("StratIntervals") # then load the package using StratIntervals The development version can be found at https://github.com/gaballench/StratIntervals.jl . A static website with documentation and vignettes can be found at https://gaballench.github.io/StratIntervals.jl . 2.2 Implementation evaluation Three sets of simulations were performed to verify the properties of the implementation under increasing sample size and different prior specifications. For all simulations, a true interval with θ 1 = 150.0 Ma and θ 2 = 100.0 Ma was chosen, and a true preservation parameter λ = 0.0. All simulations were implemented in pure Julia code, using GNU parallel ( Tange, 2011 , 2024 ) when appropriate, and run in the server gymnotus at the Instituto de Biociências, Universidade Estadual Paulista, Botucatu. All the simulation code is available at url https://github.com/gaballench/bayesian_collection_stratintervals as well as on zenodo (DOI: 10.5281/zenodo.14868381). 2.2.1 Maximum likelihood estimation The model herein implemented has three parameters (see below under Results) which are mentioned here for completeness. Endpoint parameters are θ 1 , θ 2 and the preservation parameter λ . Increasing datasets of 10, 50, 100, 1000, and 10000 occurrences were simulated to plot the joint likelihood surface for combination of ( θ 1 , λ ), ( θ 2 , λ ), and ( θ 1 , θ 2 ), while keeping the remaining parameter fixed to the true value. The true value of each parameter was then compared to the likelihood surface to see if parameter estimation is consistent. 2.2.2 Bayesian inference estimation Increasing sample size of 10, 50, 100, 200, and 500 occurrences were simulated under the true stratigraphic interval. Then, MCMC sampling was used to calculate the posterior distributions of all three parameters. The no-U-turn sampler (NUTS) and the Hamiltonian Monte Carlo sampler show advantages in situations where parameters may be correlated as appears to be the case of Beta models ( Wang, 2005 ), therefore we sampled through 10000 iterations using the NUTS ( Hoffman and Gelman, 2014 ) as it required to set less quantities than the Hamiltonian Monte Carlo. We are here mostly interested in the estimation of endpoint parameters θ 1 , θ 2 , and following the results suggesting an appropriate behaviour of λ estimation during simulations and maximum likelihood estimation (see below), it was fixed to its true value instead. Normal priors N ( µ, σ ) were set in four different combinations: Correct–incorrect, and informative– uninformative. A correct prior is defined as one centred at the true parameter value, whereas an incorrect prior is centred elsewhere and outside the stratigraphic interval (e.g. 20 units away from the true value). An informative prior (i.e., one with low variance) was set to σ = 1.0, whereas an uninformative one to σ = 20.0. Thus, for θ 1 the four combinations of priors are N (150.0, 1.0), N (150.0, 20.0), N (150.0 + 20.0, 1.0), and N (150.0 + 20.0, 20.0), whereas for θ 2 these are N (100.0, 1.0), N (100.0, 20.0), N (100.0 − 20.0, 1.0), and N (100.0 − 20.0, 20.0). The effect of co-estimating θ 1 and θ 2 was assessed by running the simulations with priors set as described above, whereas the effect of estimating just one at a time was measured by fixing the other parameter to its true value and sampling just one at a time. A total of 1000 simulations were run for each of these eight combinations: prior correctness, informativeness, and co-estimation. A separate set of simulations was carried out for assessing the impact of the value of λ on the estimation of endpoint θ parameters. Here, a dataset of 75 occurrences was simulated from the Three-Parameter Beta distribution with the true value of λ varying (-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5), whereas the true value of θ 1 and θ 2 were set to 550.0 and 500.0 Ma respectively, for a total of 7 simulated datasets. Priors on model parameters were set to θ 1 ∼ N ( µ = 550.0 , σ = 10.0), θ 2 ∼ N ( µ = 500.0 , σ = 10.0) and λ ∼ N (0.0, 2.0). The posterior distributions of model parameters were sampled using 100000 iterations and the NUTS. A generalisation of th Exponential distribution, called herein the Reflected-Offset Exponential (RefOffExp or RefOffExponential) with parameters s, ρ, o where s is the scale, ρ the reflection constant (either 1 or -1) and o is the offset, has been implemented in order to allow priors which start from a point different from zero, as well as a reflection parameter which allows it to decrease towards the left instead of the right in the standard Exponential distribution. 2.3 Empirical examples 2.3.1 Empirical example 1: The Palynomorphs of the Cerrejón formation The Cerrejón formation records the sedimentary dynamics during the Palaeocene in northern South America. The nearly 1000 m stratigraphic section is rich in plant macro and microfossils, preserved in a coastal forest environment. Also, it is very important in economic terms, because it preserves multiple coal seams of very high quality, which has been extensively exploited in the Cerrejón open-pit coal mine. Jaramillo et al. (2007) studied the Palynomorphs of a long section from the Cerrejón coal mine spanning the Manantial, Cerrejón, Tabaco, and Palmito formations from base to top. About 195 samples are available in the combined 4752 and 4774 sections, which were identified, counted, and subject to palaeoecological analyses. The data on the 216 different Palynomorphs found were herein analysed to illustrate the application of stratigraphic interval estimation, as well as to estimate a distribution for a co-occurrence of stratigraphic intervals. Note that in biostratigraphy, time in the stratigraphic sections (frequently cores) is measured in depth from either top or base. Accordingly, time in this application will represent position in the stratigraphic column instead of absolute time. Stratigraphic intervals of assorted Palynomorphs Three Palynomorphs were chosen for estimating the stratigraphic interval parameters using depth occurrences. Most of the 216 taxa were discarded because they (1) spanned the whole sequence, or (2) had very few occurrences. From the selected Palynomorphs, three were chosen which range from 26 to 46 occurrences, and varying in the presence of occurrences close to the top. In two of them, we co-estimated both θ 1 , θ 2, whereas in a third one we fixed θ 2 instead of co-estimating both endpoint parameters, because the latest occurrences were very close to the top of the stratigraphic column. The priors on θ 1 ∼ N ( µ = max(occurrences) , σ = 140.0), θ 2 ∼ Exp ( θ = min(occurrences) / 2), and λ ∼ N ( µ = 0.0 , σ = 2.0) were used for each stratigraphic interval. These priors are rather uninformative, following the results herein on the possible sensitivity to the prior in the presence of small sample sizes. The parameter θ 2 was fixed to 0.0 in Ischyosporites problematicus because the latest occurrences were very close to the top of the section, thus estimating only θ 1 and λ . The posterior predictive distribution was not calculated in this setting, only the posterior parameter distribution. MCMC sampling was carried out with the NUTS, with 10000 iterations. All parameter traces attained effective sample sizes above 600 ( Table 1 ). View this table: View inline View popup Download powerpoint Table 1: MCMC efficiency and convergence measures for the analyses of the Cerrejón formation Palynomorphs. Abbreviations are Monte-Carlo standard error (MCSE), effective sample size (ESS), effective sample size per second (ESS s-1). R̂ following Vehtari et al. (2021) . The lost sample of Coal Seam 115 Let us imagine a fictitious situation: While processing samples collected at the Cerrejón coal mine in northern Colombia, an analyst found out that one important sample had missing information for depth and stratigraphic section in the label, with “Coal Seam 1-something” as the only information available. We do not have first-hand stratigraphic information for the sample, but still want to figure out where in the generalised column the sample came from, as important fossil leaves were also found just at the top of that stratigraphic level and therefore a finer provenance information is key. What is the best estimate of stratigraphic position we can have? After processing the samples for Palynomorphs, six of these were found in Coal Seam 115 (true depth 285.99 m which we do not know in this fictitious scenario): Psilatricolporites sp. ( n = 70), Retitricolporites sp. ( n = 61), Psilamonocolpites grandis ( n = 44), Echinatisporis minutus ( n = 41), Retitricolpites sp. ( n = 23), and Matonisporites sp ( n = 23). A subset of only occurrences for these six Palynomorphs was selected, and the occurrence corresponding to the level 258.88 was removed in order to mimic the absence of information for the position of that level, which reduced in 1 the number of occurrences. The goal here is to first build an age model for each Palynomorph using the posterior predictive using the preserved occurrences and the combine these six distributions. The priors were set as in the previous section. The parameter θ 2 was again fixed to 0.0 in Psilatricolporites sp., Retitricolporites sp., and Matonisporites sp., for the reasons discussed above. The posterior predictive distribution was calculated here in order to conflate such distribution for each stratigraphic interval into a combined distribution for Coal Seam 115. MCMC sampling was carried out again with the NUTS, but this time using 100000 iterations because 10000 did show issues of convergence for three out of the six species. Afterwards, parameter traces attained effective sample sizes above 600 again. Quantiles were calculated with the area under the curve of the conflated PDF and then compared with the true position for Coal Seam 115. 2.3.2 Empirical example 2: The time of origin of Barracudas The origination time of the Barracudas (family Sphyraenidae) will be estimated using a curated dataset of fossil occurrences of the family ( Ballen, 2020 ). Because the group is still living, the extinction time θ 2 will be fixed to zero during Bayesian posterior parameter estimation of θ 1 and λ . The effect of using different priors was assessed by setting θ 1 ∼ N ( µ = 100.0 , σ = 25.0), RefOffExp( θ = 50.0 , o = 50.0 , ρ = 1), and U (50.0, 300.0). The prior λ ∼ N ( µ = 0.0 , σ = 2.0) was used in all cases. Also, comparisons between removing and leaving duplicate time occurrences were carried out. It is important to note that these duplicates are occurrences from different localities, not redundancies form the same locality, and therefore have the same age just because of uncertainty in dating, and thus are not true duplicates. The effect of this remotion is to reduce sample size. The posterior predictive distribution was not calculated in this setting, only the posterior parameter distribution. MCMC sampling was carried out with the NUTS, with 10000 iterations. All parameter traces attained effective sample sizes above 600 ( Table 2 ). Simulations were carried out in parallel for efficiency using GNU parallel ( Tange, 2011 , 2024 ). View this table: View inline View popup Download powerpoint Table 2: MCMC sampling efficiency measures for the analyses of the origination time of Barracudas. Abbreviations are Monte-Carlo standard error (MCSE), effective sample size (ESS), effective sample size per second (ESS s-1). R̂ following Vehtari et al. (2021) . Clean represents whether duplicates were removed or not. 3 Results Each stratigraphic interval represents our known occurrences of a given lineage in geologic time. These lineages have unobserved origination and extinction times θ. The occurrences of the lineage are preserved in the fossil record in a way which has been described by different models in the literature, from uniform to distribution-free ways. It is reasonable to consider that preservation is most likely non-uniform across geologic time for a given lineage. This can be accomplished by using a generalised beta distribution which has a very flexible set of shapes governed by the parameters α and β . Together these can describe a nonuniform preservation pattern through geologic time. A reparametrisation in order to describe an equally rich set of preservation possibilities is applied here in order to lower the free parameters in the model while retaining flexibility. The data are the known sets of fossil occurrence for a given fossil assemblage, that is, a fossil community composed of different species which coexist in space and time. Each of these species is expected to be also known from other points in time (and maybe space) and thus each of the species can be termed a stratigraphic series τ composed of all the points in time at which we know the lineage to have existed (e.g., the shells representing fossilised time occurrences in Figure 1 ). Download figure Open in new tab Figure 1: The model of a stratigraphic interval of time occurrences given the parameters. A) The fossil lineage (exemplified by the Ammonite) spans the time interval [ θ 1 , θ 2 ], and then fossilised occurrences are preserved at time points shown by the fossil figures. These define the stratigraphic interval with time occurrences [ τ 1 , …, τ 6 ] which are realisations of the underlying distribution f ( τ | θ 1 , θ 2 , λ ) (uniform in this example), whose probability is defined by the recovery potential governed by λ . That last parameter describe the relative frequency of time occurrences along the stratigraphic interval. B) When λ = 0.0, the distribution of occurrences is uniform through time. C) When it is negative the trend is decreasing in time, and D) the opposite when it is positive. The independent variable is time, depicted by τ in all figures, whereas the dependent variable is the probability density function f ( τ | θ 1 , θ 2 , λ ). E) Assemblage generating process. Each lineage has an origination and extinction time designed by θ. . Coloured ammonites represent occurrences in time, and we assume that all except those inside the blue dashed region have their age given, so that we are interested in estimating their age in that region, which represent the coexistence time of all five lineages. Such time is given by the conflation (blue distribution below) of all the probability density functions (black lines in top four panels) representing the stratigraphic interval of each species. Ammonite reconstructions available under Creative Commons license (in A) thanks to Kouhei Futaka and Izolende, and thanks to Jorge W. Moreno-Bernal (in E). Besides providing a flexible Bayesian method for stratigraphic intervals, the goal of the present method is to construct a distribution for the time of co-occurrence of stratigraphic intervals. We can think of these cooccurrences as assemblages of species as being sampled at some point in time co-occurring in time and space. We want to use their time model, that is, the stratigraphic interval for each component of the assemblage, to construct the co-occurrence time. 3.1 The model: Likelihood function of a stratigraphic interval Let τ be a random variable representing an occurrence event in time. A stratigraphic interval ( Strauss and Sadler 1989 ) is a segment of time bound by an origination time θ 1 and an extinction time θ 2 , where a lineage is existing on earth. Although in principle such interval in continuous, we can only access or sample it through point time occurrences inside the bound interval, whose distribution is governed by a parameter or set of parameters describing the potential of preservation. We construct a function that describes the likelihood of observing a given time occurrence with the four-parameter Beta distribution (e.g., Wang, 2005 ) defined over x ∈ [ a, c ] as where B ( α, β ) is the Beta function. We have renamed here the variables y = τ , a = θ 1 , and c = θ 2 in order to describe the process above. Now let us reparametrise the function so that it depends only on a parameter λ instead of α and β in a similar way as ( Wang et al., 2016 ). This has the advantage of requiring less data while still giving the likelihood function enough flexibility to describe different preservation scenarios. We shall call this function the three-parameter Beta distribution. Note however that here we do not fix any of the interval parameters, which gives the likelihood function: and which gives and Assuming that time occurrences are independent and identically distributed, we can also calculate the likelihood of observing a vector representing the collection of N time occurrences given a stratigraphic interval by the product of likelihoods of individual time occurrences: 3.2 The posterior distribution of parameters and posterior predictive of τ We now apply the Bayes’ theorem in order to construct the posterior densities for these parameters given the observed vector of occurrences in time: Depending on the specific choice of priors for θ 1 , θ 2 , and λ , this expression may have a closed form (e.g., using conjugate priors). However, for several prior functions this expression may not be closed and therefore we can resort to MCMC sampling in order approximate the integral in the denominator. Now that we have a posterior of parameters, we can sample from this conditional distribution in order to build a posterior predictive distribution of τ by integrating over the parameter uncertainty. This posterior predictive distribution represents the probabilistic model for the occurrences through time while accounting for uncertainty in endpoint and preservation parameter estimation. 3.3 The distribution of the co-occurrence of multiple stratigraphic intervals Combining different distributions is not straightforward ( Genest and Zidek, 1986 ). However, the conflation of probability density functions is a useful procedure which combines them provided that each of them is independent ( Hill, 2011 ; Hill and Miller, 2011 ). Assuming that each distribution describing the posterior predictive distribution of each interval is independent, and simplifying notation so that τ̃ = τ̃ | τ , we can define the composite distribution of τ̃ for M intervals as the conflation of individual posterior predictive distributions: Such distribution is useful e.g. when we want to build credible intervals for the time τ of co-occurrence of stratigraphic intervals. The conceptual interpretation is as follows: M different lineages as represented by their stratigraphic intervals should coexist at most for some time interval when they all were alive. As the conflation of densities is a density itself, it can be used for asking questions on the probability of co-existence of lineages during some arbitrary time interval given the distribution. 3.4 Implementation and parameter identifiability 3.4.1 Asymptotic properties and parameter estimation under maximum likelihood and Bayesian inference The likelihood surface was found to be flat with respect to the endpoint parameters, whereas it it well-behaved for the preservation parameter. The endpoint parameter estimators are biased, which means that the first appearance datum (FAD) and last appearance datum (LAD), instead of the true value, always maximise the likelihood surface. For small sample size, they are also unidentifiable. However, as n rises and the LAD and FAD tend to approach the true endpoint values, and therefore the MLE converges to the true value ( Figure 2 ). Download figure Open in new tab Figure 2: Log-likelihood surface of different pairs of parameters from the ThreeParBeta( τ ; θ 1 , θ 2 , λ ) distribution as a function of a fixed dataset of varying size ( n = 10 or n = 10000). Left-most columns are 2-dimensional representations of the 3-dimensional surfaces to the right. Isoclines represent log-likelihood, the brighter the higher. The x and y axes are the parameter pairs ( θ 1 , λ ), ( θ 2 , λ ), and ( θ 1 , θ 2 ) from top to bottom. In all cases, the remaining parameter was fixed to its true value. The pair of true values is plotted as a dot, which changes in colour just to aid in visualisation. The estimation of λ is unbiased, whereas the estimation of both θ 1 and θ 2 is biased and only converges to the true value asymptotically. The Bayesian version of the method, however, does not show bias in estimation and is better behaved in terms of estimation accuracy ( Figure 3 ). The cost of co-estimating the both endpoint parameters in simulations as compared to fixing one of them to the true value is negligible (Supplementary Figure S1) and therefore we will discuss here only simulations estimating one endpoint parameter while fixing the other to the true value. Download figure Open in new tab Figure 3: Posterior distribution of a single simulated dataset of varying sample size from the ThreeParBeta( τ ; θ 1 , θ 2 , λ ) distribution to illustrate the effect of sample size ( x axis) and prior variance on θ 1 and θ 2 ( y axis, in Ma). A) Posterior distribution of θ 1 with prior N ( µ = 150.0 , σ = 1.0). B) Posterior distribution of θ 2 with prior N ( µ = 100.0 , σ = 1.0). C) Posterior distribution of θ 1 with prior N ( µ = 150.0 , σ = 20.0). D) Posterior distribution of θ 2 with prior N ( µ = 100.0 , σ = 20.0). The dashed line represents the true value of the parameter in all plots. When simulating multiple datasets (1000) of varying sample size, we find that the distribution of posterior medians concentrate around the true value ( Figure 4 ), and the variance of this sample of medians decreases as we increase in sample size. This measures the precision of estimation of the model in general, over multiple simulated datasets. Download figure Open in new tab Figure 4: Distribution of posterior medians ( y axis, in Ma) across 1000 simulations for combinations of correct and incorrect prior, and informative and uninformative prior, as a function of sample size ( x axis). A) Correct and informative prior ( θ 1 ∼ N (150.0, 1.0)). B) Incorrect and informative prior ( θ 1 ∼ N (150.0 + 20.0, 1.0)). C) Correct and informative prior ( θ 2 ∼ N (100.0, 1.0)). D) Incorrect and informative prior ( θ 2 ∼ N (100.0 −20.0, 1.0)). E) Correct and uninformative prior ( θ 1 ∼ N (150.0, 20.0)). F) Incorrect and uninformative prior ( θ 1 ∼ N (150.0 + 20.0, 20.0)). G) Correct and uninformative prior ( θ 2 ∼ N (100.0, 20.0)). H) Incorrect and uninformative prior ( θ 2 ∼ N (100.0 − 20.0, 20.0)). The dashed line represents the true value in all cases. The inter-quartile range (IQR) of the posterior distribution, which measures better the spread in asymmetrical ones, also reduces as we grow in sample size ( Figure 5 ), therefore showing a gain in precision within each posterior estimation. This measures the precision of estimation conditional to a specific dataset. However, when a prior is incorrect and informative, the precision actually decreases with sample size ( Figures 5B,D ). This is not the case for incorrect–uninformative, correct–informative, and correct-uninformative priors, all of which show an increase in precision with increase in sample size. Download figure Open in new tab Figure 5: Distribution of posterior inter-quartile range (IQR, y axis) across 1000 simulations for combinations of correct and incorrect prior, and informative and uninformative prior, as a function of sample size ( x axis). A) Correct and informative prior ( θ 1 ∼ N (150.0, 1.0)). B) Incorrect and informative prior ( θ 1 ∼ N (150.0 + 20.0, 1.0)). C) Correct and informative prior ( θ 2 ∼ N (100.0, 1.0)). D) Incorrect and informative prior ( θ 2 ∼ N (100.0 − 20.0, 1.0)). E) Correct and uninformative prior ( θ 1 ∼ N (150.0, 20.0)). F) Incorrect and uninformative prior ( θ 1 ∼ N (150.0 + 20.0, 20.0)). G) Correct and uninformative prior ( θ 2 ∼ N (100.0, 20.0)). H) Incorrect and uninformative prior ( θ 2 ∼ N (100.0 − 20.0, 20.0)). The mean squared error (MSE) lowers by raising sample size ( Figure 6 ), which is consistent with the results of the IQR, and highlights both the increase in accuracy and precision during estimation across multiple simulated datasets from the same model. The decrease in MSE in consistent regardless of the presence of incorrect and/or uninformative priors. Download figure Open in new tab Figure 6: Mean squared error (MSE, y axis) across 1000 simulations for combinations of correct and incorrect prior, and informative and uninformative prior, as a function of sample size ( x axis). A) Correct and informative prior ( θ 1 ∼ N (150.0, 1.0)). B) Incorrect and informative prior ( θ 1 ∼ N (150.0 + 20.0, 1.0)). C) Correct and informative prior ( θ 2 ∼ N (100.0, 1.0)). D) Incorrect and informative prior ( θ 2 ∼ N (100.0 − 20.0, 1.0)). E) Correct and uninformative prior ( θ 1 ∼ N (150.0, 20.0)). F) Incorrect and uninformative prior ( θ 1 ∼ N (150.0 + 20.0, 20.0)). G) Correct and uninformative prior ( θ 2 ∼ N (100.0, 20.0)). H) Incorrect and uninformative prior ( θ 2 ∼ N (100.0 − 20.0, 20.0)). The value of λ and therefore the shape of the Three-Parameter Beta distribution have an influence on the posterior variance of the θ endpoint parameters (also see Supplementary Figure S2). When λ = 0.0, and therefore the model has a uniform shape, the posterior variance of both θ 1 and θ 2 is the same ( Figure 7A ). However, when λ is positive (e.g., 1.5), the tail goes to the present ( Figure 7B ), whereas if it is negative (e.g., = -1.5), there is a tail towards the past ( Figure 7C ). Download figure Open in new tab Figure 7: Effect of λ on θ 1 and θ 2 for a sample size of 75 occurrences sampled from distributions with different values of λ . A) λ = 0.0. B) λ = 1.5. C) λ = −1.5. In all plots the green line represents the ThreeParBeta( τ ; θ 1 = 550.0 , θ 2 = 500.0 , λ ) from where the 75 occurrences were sampled. Blue lines represent the rather uninformative prior centred at the true value on each ( θ 1 ∼ N ( µ = 550.0 , σ = 10.0) and θ 2 ∼ N ( µ = 500.0 , σ = 10.0)). The black lines represent the posterior distribution. In all plots the x axis is time (Ma) whereas the y axis is density. 3.4.2 The effect of priors The bias and semi-identifiability of endpoint parameters in maximum likelihood estimation is not seen when using Bayesian inference. By applying priors to the θ. parameters we end up with a posterior surface which allows proper statistical inference, that is, parameter identifiability ( Figure 3 ). Because of this important role of the priors, we explore now the behaviour of the Bayesian model under two situations: Priors for θ 1 , θ 2 with relative low and high variance, and priors which are consistent or inconsistent with the true parameter value, that is, whether the true value is near the distribution mean, or very off from it. This allows to verify how sensitive to the prior it is as a function of sample size. When applying priors that are centred on the true parameter values, both low- and high-variance priors are capable of estimating the posterior distributions that include the true value. The effect of the prior variance on the posterior variance however reduces as we grow in sample size ( Figures 3 , 5 ). If an incorrect prior is chosen, one could potentially run into issues where the posterior do not include the true value even for large sample sizes. By applying a prior which is incorrect in the sense that it is centred away from the true values, we can assess the effect of informative and uninformative priors on the posterior distriutions. We find that a low-variance prior which are centred away from the true values have difficulties in converging to the true values even for large sample sizes. This effect is however not seen in high-variance priors, which quickly are capable of converging to the region around the true value as sample size grows ( Figures 4 , 5 , 6 , see panels B, D, F, H for incorrect priors). Overall the method is well behaved as regardless of whether priors are correct or incorrect, informative or uninformative, the mean squared error decreases with an increase in sample size ( Figure 6 ). The interquartile ranges tend to increase with sample size only for priors which are incorrect and informative, whereas they always decrease for correct informative and uninformative priors, as well as for incorrect uninformative priors ( Figure 5A,C,E–H ). 3.5 Empirical examples 3.5.1 Empirical example 1: The Palynomorphs of the Cerrejón formation Stratigraphic intervals of select Palynomorphs Posterior sampling showed good efficiency and convergence ( Table 1 ), and recovered behaviours also seen in simulations. First, stratigraphic intervals where λ created a tail towards θ 1 by decreasing the number of occurrences towards the present also made the inference of this parameter to show higher variance than in cases where λ was closer to 0.0 (compare Figure 8A and 8B ). All posterior distributions of θ 1 were mostly contained into the measured sections ( Figure 8 ). Posterior distribution of θ 2 was bound between 0.0 and the latest occurrence, however, it did not tend to follow the shape of the prior distribution for Arecipites regio , thus suggesting a stronger influence of more data (46 occurrences, Table 3 , Figure 8A ). In contrast, the posterior was similar to the prior in Scabratriporites triangularis , presumably due to the smaller amount of data (26 occurrences, Table 3 , Figure 8C ). The posterior distribution of θ 2 was consistently different from the prior for all species, regardless of sample size. The posterior distribution of λ were always quite different from 0.0, which means that the stratigraphic interval model is different from uniform. Download figure Open in new tab Figure 8: Stratigraphic intervals of A) Arecipites regio ( n = 46), B) Ischyosporites problematicus ( n = 33), and C) Scabratriporites triangularis ( n = 26). The parameter θ 2 was fixed for Ischyosporites problematicus as its latest occurrences were very close to the uppermost samples in the stratigraphic column. Offset left is the generalised stratigraphic column of the Cerrejón formation starting from the top of core WVR-0475 and finishing below core WVR-4774. Subplots A–C are aligned with the generalised column on the depth (m) axis. Lowermost densities correspond to θ 1 , uppermost ones to θ 2 , and inlay plots to λ . Prior distributions are blue, whereas posterior distributions are black. Data are the grey bars in the histogram. Dashed grey lines are the posterior quantiles 0.025, 0.5, and 0.975. View this table: View inline View popup Download powerpoint Table 3: Summary statistics for the posterior sampling of three Palynomorphs from the Cerrejón formation. Percentage columns are quantiles at that given probability value. Abbreviations are standard deviation (SD). The lost sample from Coal Seam 115 The use of six Palynomorphs for building a conflated distribution for their co-occurrence restricts the possible depth of the mystery sample to the 95% interval 341.8m–54.4m, with a width of 287.4 m. The median of the conflated distribution was 193.8 m. This interval includes the true depth of the Coal seam 115 which is 285.99 m ( Figure 9 ). Four coal seams are found in this interval, 110, 115, 150, and 155. The 95% interval of the conflated distribution was in general restricted to the stratigraphic span of section WRV-4752. Download figure Open in new tab Figure 9: The most probable depth for a mysterious sample labelled Coal Seam 1–something as inferred from the conflation of the posterior predictive distributions of the Palynomorphs recovered from that sample. Offset left is the generalised stratigraphic column of the Cerrejón formation starting from the top of core WVR-0475 and finishing below core WVR-4774. The true depth (Coal Seamn 115) falls within the highest density interval of the conflated distribution, although its median is not aligned to the true value. Most of the Palynomorph age models, represented by their posterior predictive distributions, agreed more or less on the interval which determined in great part the shape of the conflation. However, Psilatricolporites sp. showed a wider distribution, not contributing much to the final shape because of its greater variance when compared to the posterior predictive distributions of the remaining Palynomorphs. 3.5.2 Empirical example 2: The time of origin of Barracudas Leaving the original dataset with time duplicates (i.e., occurrences with the same time but from different localities, therefore not complete duplicates) results in 75 occurrences. In this scenario the shape of the posterior distribution is quite insensitive to the choice of the prior. Removing duplicates leaves 33 occurrences, for which the choice of the prior impacts somewhat the shape of the posterior distribution. The most important factors influencing the sensitivity of the posterior seem to be both the prior mean and variance ( Fig. 10 ). The Normal prior tends to have higher density around mean ( σ 2 = 625.0), and therefore tends to push more the posterior towards its mean, in this case, into the past. In contrast, the priors RefOffExponential and Uniform, with more variance ( σ 2 = 2500.0 and σ 2 = 5208.33 respectively), tend to have a slighter impact on the posterior distribution. Download figure Open in new tab Figure 10: Origination time of the Barracudas. Left column are datasets without removing time duplicates. Right column are datasets removing time duplicates. Top row is a prior θ 1 ∼ N (100.0, 25.0), middle row is a prior θ 1 ∼ RefOffExp(50.0, 50.0, 1), and bottom row is a prior θ 1 ∼ U (50.0, 300.0). Axes represent time (in Ma) on the horizontal, and density on the vertical. Prior distribution in blue, posterior distribution in black. Data are the grey bars in the histogram. The posterior distribution of θ 1 , the origination time, has a highest posterior density interval between 54.3–53.1 Ma (mean age 54.2–54.4 Ma depending on the prior) when using the complete dataset, regardless of the prior used ( Table 4 ). Such interval grows to 113.8–53.8 Ma (narrowest, RefOffExponential prior) and even to 146.4–54.1 Ma (widest, Uniform prior) when using the reduced dataset, this time showing more influence of the prior. In this case, the posterior mean was 70.1 Ma (RefOffExponential prior), 79.1 Ma (Uniform prior), or 83.5 Ma (Normal prior). View this table: View inline View popup Download powerpoint Table 4: Summary statistics for the posterior sampling of the origination time of Barracudas. Percentage columns are quantiles at that given probability value. Clean represents whether duplicates were removed or not. Abbreviations are standard deviation (SD). The posterior distribution of λ included zero but showed a prevalence of values on the negative side of the axis (e.g., -0.4–0.1) for the complete dataset, and was consistently restricted to the negative side of the axis for the reduced dataset (e.g., -3.4–-1.2). This means that the distribution of occurrences tended to have a tail towards the past, which makes the posterior of θ 1 wider than if it was uniform (cf. Figure 7 ). 4 Discussion This is one of the first Bayesian methods for stratigraphic intervals to be available in code, and the first Julia package dedicated to estimation of stratigraphic intervals. It is the first to allow flexible specification of a large number of prior distributions, both via the Distributions.jl package or implementing custom ones. Also, the possibility of specifying custom MCMC samplers helps attain convergence very easily. The choice of the NUTS has proven to be quite useful in general cases, and more efficient than the Hamiltonian Monte Carlo, which requires to set more arguments and may run into issues of convergence. The Turing.jl package also provides other samplers including Random-walk Metropolis-Hastings, Gibbs, and compositional Gibbs, thus extending the possibilities of fine-grain control over MCMC sampling. Finally, allowing to fix parameters to given values allows to use stratigraphic intervals in situations where e.g. the lineages in question are extant, but not restricting the user to this as happens with the Adaptive Beta of Wang et al. (2016) . Maximum likelihood and Bayesian estimation The bias in estimation where the FAD and LAD are the maximum likelihood estimates is an issue present in other methods for stratigraphic intervals, not just the present one ( Strauss and Sadler, 1989 ; Marshall, 1990 ; Rivadeneira et al., 2009 ; Wang et al., 2016 , T. Quental, pers. comm.). This makes the use of Bayesian methods more convenient for avoiding biased estimates of origination and extinction times. Wang (2005) studied the maximum likelihood properties of the Four-Parameter Beta distribution and found similar issues with parameter estimation, also advocating for constraining the likelihood surface using uninformative priors on the endpoint parameters θ 1 , θ 2 . He discussed the conditions under which parameter estimation is possible, through the use of Bayesian methods. The approach herein is in line with his findings, but is of more general application as it does not enforce the use of a specific prior on endpoint parameters. Wang (2005) examined the correlation between parameters and concluded that endpoints tend to be correlated. This is an issue for inference using Metropolis-Hasting samplers, and it was unnecessary to examine it here because both the Hamiltonian and NUTS are very efficient traversing surfaces in these situations. Priors The behaviour of different priors was in general terms similar in that given sufficient data, the likelihood function is capable of dominating the information content, thus determining to a great extent the parameter posterior distributions. However, it is clear that even fairly uninformative priors can exert some influence when sample size is small ( Figure 10 ). It is important to note that priors which are very off form the true values (herein called incorrect priors) have an undesirable influence on the posterior distributions, which is proportional to the information content in the prior distribution, and therefore inversely proportional to the variance. Thus, even if a prior is incorrect because its mean is off from the true value, we can make it behave more properly by setting a larger variance for it, that is, making it fairly uninformative. These findings suggest that the method is robust to prior specification, and it is recommended that empirical applications with small sample size (around an order to 50 occurrences or less) specify wide priors in order to avoid any undesirable influence of the prior on the posterior. Although priors which are both incorrect and informative behave very badly in the sense that they do not recover the true value in simulations, they tend to approach the true values which may eventually converge if letting simulations to run for very long and with very large sample sizes (above 500 in this study). Comparison with other methods The Three-Parameter Beta model, in the Bayesian formulation, shares features with several methods already available. On one side, it assumes continuous sampling, as is common to available tools, but in contrast with methods for discrete sampling ( Weiss and Marshall, 1999 ; Weiss et al., 2004 ). It can accommodate both uniform ( Strauss and Sadler, 1989 ; Solow, 1993a ; Ferraz et al., 2003 ; Weiss and Marshall, 1999 ; McInerny et al., 2006 ; Wang et al., 2016 ) and non-uniform ( Solow, 1993b ; Marshall, 1997 ; McCarthy, 1998 ; Weiss et al., 2004 ; Ferraz et al., 2003 ; Wang et al., 2016 ) preservation potentials, thus being similar to many Class 1 and Class 2 methods, and in contrast with non-parametric or distribution free methods ( Marshall, 1994 ; Solow and Roberts, 2003 ). The method can condition on a fixed endpoint parameter like many do (e.g., Wang et al., 2016 ), in particular those tailored for sighting data in conservation biology, which assume that the oldest occurrence is an arbitrary zero point, and then estimate the extinction time ( Ferraz et al., 2003 ; Solow, 2005 ; Kodikara et al., 2020 ), but it is not forced to such conditioning and therefore can estimate both the origination and extinction times like some methods from palaeobiology ( Strauss and Sadler, 1989 ; Marshall, 1994 ). It is unique however in allowing the user to specify any prior they want, which is supported by the Distributions.jl package or implemented in StratIntervals.jl (e.g., the RefOffExponential distribution). Also, it is unique among many Bayesian methods in that it is fully implemented in a software package which can be easily integrated with others from the Julia language ecosystem. This helps in allowing the user to extend the models available into three directions: Choice of priors, which multiple works have found to matter ( Alroy, 2016a , b ; Solow, 2016a , b ), choice of samplers, therefore allowing to improve sampling efficiency, and the possibility to extend methods by using different likelihood functions. Empirical applications The combined effect of the tail of the distribution towards the past, as well as the lower sample size in the reduced dataset explains why the posterior distribution of θ 1 was considerably wider, and should be considered when assessing the results of the application of the present method to empirical datasets. Although the choice of the prior for θ 1 has some influence on the posterior distribution, this is greatly reduced by sample size, under the same inferential settings. If the distr ibution of occurrences is concentrate towards any of the endpoints, the posterior estimation of the opposite endpoint will suffer from wider variance. This effect is seen both in the Palynomorph and the Barracudas datasets. This is a feature of the data and the distribution rather than a pathological behaviour of the method. It is comparable to trying to fit a curve to an asymmetric distribution, variations on the tailed side can vary to a larger extent than the part attached to the other side. This makes sense as more of the density of the information in the data are concentrated to the other side of the distribution, and therefore we have less information on the tailed side ( Figures 8 and 10 ). The conflation of distributions showed desirable properties in real-world applications. The fact that we recovered a reasonable position for the Coal seam 115 just by observing the occurrences of only six Palynomorphs present in the sample is exciting, as normally much more taxa are available in real-world samples for biostratigraphic analyses ( Jaramillo et al., 2007 ). Its potential however is not restricted to micropalaeontology but instead to any application where combinations of distributions is of interest. For instance, Ballen and Reinales (2024) showed its application to the construction of secondary calibrations in divergence time estimation. Although methods for stratigraphic intervals were born from the field of quantitative palaeobiology, their application to questions in evolutionary biology and related fields can be much wider, especially in the field of divergence time estimation. Calculating the co-occurrence of multiple lineages in time has shown to be useful in estimating the time of co-occurrence of collections of objects in time. This method does not require the intervals to represent living organisms, and therefore can be easily applied to historical objects such as musical pieces which are preserved in multiples sources, e.g. Gregorian chant melodies in manuscripts ( Ballen et al., 2024 ). These scenarios are quite common in the study of cultural evolution and digital humanities, which are interesting fields where data sometimes correspond to series of occurrences in time, the methods presented herein are promising. Conflict of interest The author declares no conflict of interest. Author’s contribution The author contributed to conceptualisation, development, and writing. Acknowledgments This research was supported by FAPESP through a doctoral scholarship (2014/11558-5), an internship abroad (2016/02253-1), and a postdoctoral fellowship (2023/07838-1), as well as by the BBSRC (grant BB/T01282X/1 awarded to Mario dos Reis). I thank Mario dos Reis and Ziheng Yang for valuable feedback during early phases of this project, I could not have finished it without their guidance and criticism in methods development. I thank Phil Donoghue, Davide Pisani and the Molecular Palaeobiology Group at University of Bristol for inviting me to deliver a talk on the topic of this paper when it was still in early phases of development. I thank Tiago Quental and Gustavo Burin for discussions on methods for palaeobiology and macroevolution throughout the years, and for early discussions on methods for stratigraphic intervals since my first year of doctoral studies. I thank Carlos Jaramillo for supporting my career in palaeontology, and for encouraging me and several other Colombian students to build analytical and computational skills in palaeobiology. Discussions on palaeontology and geology with many colleagues (Javier Luque, Jorge W. Moreno, Camila Martinez, Maria Camila Vallejo, Federico Moreno, Andŕes Ćardenas, Edwin Cadena, Camilo Montes, Catalina Suarez, and Sandra Reinales) have been a valuable source of inspiration. Living ammonite reconstruction in Figure 1A thanks to Kouhei Futaka, available under Creative Commons license ( https://doi.org/10.7875/togopic.2020.08 ), and ammonite drawing in Figure 1A thanks to Izolende, also available under Creative Commons ( https://commons.wikimedia.org/wiki/User:Ilzolende ). Special thanks to Jorge W. Moreno-Bernal for the ammonite reconstructions in Figure 1E . Finally, I thank Sandra Reinales for supporting me personally and professionally in this and many other projects throughout the years. Funder Information Declared FAPESP Footnotes Numerous typos are corrected, and some reference errors among figure legends. Other than that the content remains the same. References ↵ Alroy , J . ( 2014 ). A simple Bayesian method of inferring extinction . Paleobiology , 40 ( 4 ): 584 – 607 . OpenUrl Abstract / FREE Full Text ↵ Alroy , J . ( 2015 ). Current extinction rates of reptiles and amphibians . Proceedings of the National Academy of Sciences , 112 ( 42 ): 13003 – 13008 . OpenUrl Abstract / FREE Full Text ↵ Alroy , J . ( 2016a ). A simple Bayesian method of inferring extinction: Reply . Ecology , 97 ( 3 ). ↵ Alroy , J . ( 2016b ). Reply to Solow: Sense and nonsense in the choice of extinction priors . Proceedings of the National Academy of Sciences , 113 ( 9 ): E1133 – E1133 . OpenUrl FREE Full Text ↵ Ballen , G. A . ( 2020 ). New records of the genus Sphyraena (Teleostei: Sphyraenidae) from the Caribbean with comments on dental characters in the genus . Journal of Vertebrate Paleontology , 40 ( 6 ): e1849246 . OpenUrl CrossRef ↵ Ballen , G. A. , Mühlová , K. H. , and Hajič jr. , J. ( 2024 ). What did the dove sing to Pope Gregory? Ancestral melody reconstruction in Gregorian chant using Bayesian phylogenetics . BioRxiv , pages 1 – 17 . ↵ Ballen , G. A. and Reinales , S . ( 2024 ). tbea: Tools for pre-and post-processing in Bayesian evolutionary analyses . bioRxiv, pages 2024 – 06 . ↵ Besaņcon , M. , Papamarkou , T. , Anthoff , D. , Arslan , A. , Byrne , S. , Lin , D. , and Pearson , J . ( 2021 ). Distributions.jl: Definition and modeling of probability distributions in the JuliaStats Ecosystem . Journal of Statistical Software , 98 ( 16 ): 1 – 30 . OpenUrl ↵ Bezanson , J. , Edelman , A. , Karpinski , S. , and Shah , V. B . ( 2017 ). Julia: A fresh approach to numerical computing . SIAM Review , 59 ( 1 ): 65 – 98 . OpenUrl CrossRef ↵ Boakes , E. H. , Rout , T. M. , and Collen , B . ( 2015 ). Inferring species extinction: The use of sighting records . Methods in Ecology and Evolution , 6 ( 6 ): 678 – 687 . OpenUrl ↵ Bouchet-Valat , M. and Kamiński , B. ( 2023 ). DataFrames.jl: Flexible and fast tabular data in Julia . Journal of Statistical Software , 107 ( 4 ): 1 – 32 . OpenUrl ↵ Ferraz , G. , Russell , G. J. , Stouffer , P. C. , Bierregaard Jr , R. O. , Pimm , S. L. , and Lovejoy , T. E . ( 2003 ). Rates of species loss from Amazonian forest fragments . Proceedings of the National Academy of Sciences , 100 ( 24 ): 14069 – 14073 . OpenUrl Abstract / FREE Full Text ↵ Ge , H. , Xu , K. , and Ghahramani , Z . ( 2018 ). Turing: A language for flexible probabilistic inference . In International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain , pages 1682 – 1690 . ↵ Genest , C. and Zidek , J. V . ( 1986 ). Combining probability distributions: A critique and an annotated bibliography . Statistical Science , 1 ( 1 ): 114 – 135 . OpenUrl CrossRef ↵ Hill , T. P . ( 2011 ). Conflations of probability distributions . Transactions of the American Mathematical Society , 363 ( 6 ): 3351 – 3372 . OpenUrl CrossRef ↵ Hill , T. P. and Miller , J . ( 2011 ). How to combine independent data sets for the same quantity . Chaos: An Interdisciplinary Journal of Nonlinear Science , 21 ( 3 ): 033102 . OpenUrl ↵ Hoffman , M. D. and Gelman , A . ( 2014 ). The No-U-Turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo . Journal of Machine Learning Research , 15 ( 1 ): 1593 – 1623 . OpenUrl ↵ Holland , S. M. , Patzkowsky , M. E. , and Loughney , K. M . ( 2024 ). Stratigraphic Paleobiology . Paleobiology, page 1 – 18 . ↵ Jaramillo , C. A. , Bayona , G. , Pardo-Trujillo , A. , Rueda , M. , Torres , V. , Harrington , G. J. , and Mora , G . ( 2007 ). The palynology of the Cerrejón Formation (upper Paleocene) of northern Colombia . Palynology , 31 ( 1 ): 153 – 189 . OpenUrl Abstract / FREE Full Text ↵ Johnson , S. G. ( 2013 ). QuadGK.jl: Gauss–Kronrod integration in Julia . https://github.com/JuliaMath/QuadGK.jl . ↵ Kodikara , S. , Demirhan , H. , Wang , Y. , Solow , A. , and Stone , L . ( 2020 ). Inferring extinction year using a Bayesian approach . Methods in Ecology and Evolution , 11 ( 8 ): 964 – 973 . OpenUrl ↵ Marshall , C. R . ( 1990 ). Confidence intervals on stratigraphic ranges . Paleobiology , 16 ( 1 ): 1 – 10 . OpenUrl Abstract ↵ Marshall , C. R . ( 1994 ). Confidence intervals on stratigraphic ranges: Partial relaxation of the assumption of randomly distributed fossil horizons . Paleobiology , 20 ( 4 ): 459 – 469 . OpenUrl Abstract ↵ Marshall , C. R . ( 1997 ). Confidence intervals on stratigraphic ranges with nonrandom distributions of fossil horizons . Paleobiology , 23 ( 2 ): 165 – 173 . OpenUrl Abstract ↵ Marshall , C. R . ( 2010 ). Using confidence intervals to quantify the uncertainty in the end-points of stratigraphic ranges . Quantitative Methods in Paleobiology , 16 : 291 – 316 . OpenUrl ↵ McCarthy , M. A . ( 1998 ). Identifying declining and threatened species with museum data . Biological conservation , 83 ( 1 ): 9 – 17 . OpenUrl CrossRef Web of Science ↵ McCrea , R. S. , Cheale , T. , Campillo-Funollet , E. , and Roberts , D. L . ( 2024 ). Inferring species extinction from sighting data . Cambridge Prisms: Extinction , 2 : e19 . OpenUrl ↵ McInerny , G. J. , Roberts , D. L. , Davy , A. J. , and Cribb , P. J . ( 2006 ). Significance of sighting rate in inferring extinction and threat . Conservation biology , 20 ( 2 ): 562 – 567 . OpenUrl CrossRef PubMed ↵ Rivadeneira , M. M. , Hunt , G. , and Roy , K . ( 2009 ). The use of sighting records to infer species extinctions: An evaluation of different methods . Ecology , 90 ( 5 ): 1291 – 1300 . OpenUrl CrossRef PubMed ↵ Signor , P. W. and Lipps , J. H . ( 1982 ). Sampling bias, gradual extinction patterns and catastrophes in the fossil record . Geological Society of America Special Paper , 190 : 291 – 296 . OpenUrl CrossRef ↵ Smith , J. , Dowding , E. , Abdelhady , A. , Abondio , P. , Araujo , R. , Aze , T. , Balisi , M. , Buatois , L. , Carvajal-Chitty , H. , Chattopadhyay , D. , Coiro , M. , Dietl , G. , Gonzalez Arango , C. , Kevrekidis , C. , Kimmig , J. , Mychajliw , A. , Pimiento , C. , Regalado Fernandez , O. , K.M., S., Warnock , R. , Yang , T. , Yasuhara , M. , Akita , L. , Allen , B. , Anderson , B. , Andreoletti , J. , Archuby , F. , Ballen , G. , Bari , M. , Benton , M. , Bergh , E. , Brambilla , L. , Brombacher , A. , Chan , Y. , Chiarenza , A. , Chinzorig , T. , Coates , K. , Cordie , D. , Cortes-Sanchez , M. , Cruz-Vega , E. , Cybulski , J. , De Baets , K. , De Entrambasaguas , J. , Dillon , E. , Du , A. , Dunhill , A. , Erlandson , J. , Forel , M. , Foster , W. , Gates , T. , Gavryushkina , A. , Grace , M. , Grossart , H. , Hansel , P. , Harnik , P. , Hopkins , M. , Hopkins , S. , Hu , K. , Huang , H. , Irmis , R. , Jaques , V. , Jenkins , X. , Jukar , A. , Kelley , P. , Kihn , R. , Klompmaker , A. , Kocsis , A. , Kriwet , J. , Lazarus , D. , Liao , C. , Lin , C. , Louys , J. , Lozano-Fernandez , J. , Lozano-Francisco , M. , Lueders-Dumont , J. , Malve , M. , Martindale , R. , Mazzini , I. , Modenini , G. , Mondal , S. , Mondini , M. , Monferran , M. , Mulvey , L. , Nanglu , K. , Nguyen , J. , Norris , R. , O’Dea , A. , Ollendorf , A. , Orihuela , J. , Pandolfi , J. , Pereira , T. , Piro , A. , Plotnick , R. , Plaza-Torres , S. , Porto , A. , Prieto-Marquez , A. , Punyasena , S. , Quental , T. , Raja , N. , Ranaivosoa , V. , Ribas-Deulofeu , L. , Rivals , F. , Roden , V. , Rosso , A. , Saleh , F. , Salvador , R. , Saupe , E. , Schneider , S. , Sclafani , J. , Smith , M. , Souron , A. , Steinbauer , M. , Stewart , M. , Tambussi , C. , Thomas , E. , Tschopp , E. , Tutken , T. , Varela , S. , Vezzosi , R. , Villasenor , A. , Weinkauf , M. , Zanno , L. , Zhang , C. , Zhao , Q. , and Kiessling , W. ( 2025 ). Big questions in paleontology: A community-driven project to motivate new insights about the history of life on Earth . Paleobiology (In press) , XX : XX – XX . OpenUrl ↵ Solow , A. R . ( 1993a ). Inferring extinction from sighting data . Ecology , 74 ( 3 ): 962 – 964 . OpenUrl CrossRef Web of Science ↵ Solow , A. R . ( 1993b ). Inferring extinction in a declining population . Journal of Mathematical Biology , 32 : 79 – 82 . OpenUrl CrossRef ↵ Solow , A. R . ( 2005 ). Inferring extinction from a sighting record . Mathematical biosciences , 195 ( 1 ): 47 – 55 . OpenUrl CrossRef PubMed Web of Science ↵ Solow , A. R . ( 2016a ). A simple Bayesian method of inferring extinction: Comment . Ecology , 97 ( 3 ). ↵ Solow , A. R . ( 2016b ). On Bayesian inference about extinction . Proceedings of the National Academy of Sciences , 113 ( 9 ): E1132 – E1132 . OpenUrl FREE Full Text ↵ Solow , A. R. and Roberts , D. L . ( 2003 ). A nonparametric test for extinction based on a sighting record . Ecology , 84 ( 5 ): 1329 – 1332 . OpenUrl CrossRef ↵ Strauss , D. and Sadler , P. M . ( 1989 ). Classical confidence intervals and Bayesian probability estimates for ends of local taxon ranges . Mathematical Geology , 21 ( 4 ): 411 – 427 . OpenUrl CrossRef GeoRef PubMed Web of Science ↵ Tange , O . ( 2011 ). Gnu parallel–the command-line power tool . USENIX Mag , 36 ( 1 ): 42 . OpenUrl ↵ Tange , O . ( 2024 ). GNU Parallel 20240522 (’Tbilisi’) . doi: 10.5281/zenodo.11247979 . OpenUrl CrossRef ↵ Vehtari , A. , Gelman , A. , Simpson , D. , Carpenter , B. , and Bürkner , P.-C. ( 2021 ). Rank-normalization, folding, and localization: An improved R̂ for assessing convergence of MCMC . Bayesian analysis , 16 ( 2 ): 667 – 718 . OpenUrl ↵ Wang , J. Z . ( 2005 ). A Note on Estimation in the Four-Parameter Beta Distribution . Communications in Statistics - Simulation and Computation , 34 ( 3 ): 495 – 501 . OpenUrl CrossRef ↵ Wang , S. C. , Everson , P. J. , Zhou , H. J. , Park , D. , and Chudzicki , D. J . ( 2016 ). Adaptive credible intervals on stratigraphic ranges when recovery potential is unknown . Paleobiology , 42 ( 2 ): 240 – 256 . OpenUrl Abstract / FREE Full Text ↵ Weiss , R. E. , Basu , S. , and Marshall , C. R . ( 2004 ). A framework for analysing fossil record data . In Tools for constructing chronologies: crossing disciplinary boundaries, pages 213 – 230 . Springer . ↵ Weiss , R. E. and Marshall , C. R . ( 1999 ). The uncertainty in the true end point of a fossil’s stratigraphic range when stratigraphic sections are sampled discretely . Mathematical Geology , 31 : 435 – 453 . OpenUrl CrossRef GeoRef Web of Science View the discussion thread. Back to top Previous Next Posted August 25, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following A flexible Bayesian method for estimating stratigraphic intervals and their co-occurrence in time Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share A flexible Bayesian method for estimating stratigraphic intervals and their co-occurrence in time Gustavo A. Ballen bioRxiv 2025.02.13.638199; doi: https://doi.org/10.1101/2025.02.13.638199 Share This Article: Copy Citation Tools A flexible Bayesian method for estimating stratigraphic intervals and their co-occurrence in time Gustavo A. Ballen bioRxiv 2025.02.13.638199; doi: https://doi.org/10.1101/2025.02.13.638199 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7624) Biochemistry (17651) Bioengineering (13873) Bioinformatics (41887) Biophysics (21424) Cancer Biology (18566) Cell Biology (25465) Clinical Trials (138) Developmental Biology (13365) Ecology (19871) Epidemiology (2067) Evolutionary Biology (24293) Genetics (15591) Genomics (22478) Immunology (17715) Microbiology (40331) Molecular Biology (17150) Neuroscience (88492) Paleontology (666) Pathology (2828) Pharmacology and Toxicology (4817) Physiology (7635) Plant Biology (15114) Scientific Communication and Education (2044) Synthetic Biology (4286) Systems Biology (9817) Zoology (2268)