paleotree (version 1.8.2)

getSampProbDisc: Fit Models of Sampling Probability to Discrete-Interval Taxon Ranges

Description

Uses maximum likelihood to estimate per-interval sampling probability and extinction rate, given a set of discrete-interval taxon ranges from the fossil record. This function can fit models where there are different groupings of taxa with different parameters and different free-floating time intervals with different parameters.

Usage

getSampProbDisc(timeData, n_tbins = 1, grp1 = NA, grp2 = NA, est_only = FALSE, 
    iterations = 10000, initial = 0.5)

Arguments

timeData

A 2 column matrix with the first and last occurances of taxa given in relative time intervals. If a list of length two is given for timeData, such as would be expected if the output of binTimeData was directly input, the second element is used.

n_tbins

Number of time bins with different sampling/extinction parameters

grp1

A vector of integers or characters, the same length as the number of taxa in timeData, where each taxon-wise element gives the group ID of the taxon for the respective row of timeData

grp2

A vector of integers or characters, the same length as the number of taxa in timeData, where each taxon-wise element gives the group ID of the taxon for the respective row of timeData

est_only

If true, function will give back a matrix of ML extinction rates and sampling probabilities per species rather than usual output (see below)

iterations

Maximum number of iterations the optimizer is run for

initial

Values used as initial parameter value to initiate the optimizing function. The same starting value is used for all parameters

Value

• If est_only = TRUE, a matrix of per-taxon sampling and extinction parameters is output. If est_only = FALSE (default), then the output is a list:
• Title
  Gives details of the analysis, such as the number and type of parameters included and the number of taxa analyzed
• parameters
  Maximum-likelihood parameters of the sampling model
• log.likelihood
  The maximum support (log-likelihood) value
• AICc
  The second-order Akaike's Information Criterion, corrected for small sample sizes
• convergence
  A number indicating status of convergence; if 0 then convergence was reached; see help file for optim for the respective error list
• message
  Messages output by optim; check to make sure that model convergence occurred
• If multi-class models are fit, the element $parameters will not be present, instead there will be several different elements that describe the characteristic parameter 'components' for each class, rather than representing the parameters of actual taxa in that class. As noted in the $Title, these should not be interpreted as the actual rates/probabilities associated with any real taxa but rather as factors which must be multiplied in combination with the estimates for other classes to be meaningful. For example, for taxa of a given group in a given time bin, their extinction rate is the extinction rate component of that time bin times the extinction rate component of their group. Completeness estimates will be output with these parameters as long as classes are not overlapping, as those estimates would not otherwise refer to meaningful groups of taxa. As the limited time-window option of Foote (1997) is not implemented, any taxa listed as being in a bin with start time 0 and end time 0 (and thus being extant without question) are dropped before the model fitting it performed.

Details

This function uses maximum-likelihood solutions obtained by Foote (1997). These analyses are ideally applied to data from single stratigraphic section but potentially are applicable to regional or global datasets, although the behavior of those datasets is less well understood. getSampProbDisc allows for a considerable level of versatility in terms of the variation allowed among taxa in model parameters (extinction rate and sampling probability/rate). Essentially, taxa are divided into different (possibly overlapping) classes which have 'average' parameter values. These average parameters are multiplicatively-combined to calculate per-taxon parameters. For example, perhaps a user hypotheses that taxa which live in a particular environment have different characteristic sampling probabilities, that taxa of several different major clades have different characteristic sampling probabilities and that there may be several temporal shifts in the characteristic extinction rate. The classification IDs for the first two can be included as per-taxon vectors (of either character or integers) as the arguments 'grp1' and 'grp2' and the hypothesized number of temporal breaks included as the n_tbins argument. A model where taxa differ in parameters across time, clades and environments will thus be fit and the AIC calculated, such that this model's fit can be compared to other (probably simpler) models. By default, a simple model is fit to the range data where all taxa belong to a single class, with a single characteristic extinction rate and a single characteristic sampling probability. The timebins option allows for time intervals with different characteristic model parameters which are not defined a priori. The boundaries between time intervals with different characteristic parameters will thus be additional free-floating parameters included in the AIC calculation. If you have the prior inclination that sampling/extinction changed at a specific point in time, then seperate the taxa that originated before and after that point as two different groups and include those classifications as a group in the arguments. This function does not implement the finite window of observation modification for dealing with range data containing extant taxa and thus that continues up to the recent (Foote, 1997). This is planned for a future update, however. For now, data input into this function should be for taxa that have already gone extinct by the modern and are not presently extant. The intervals in timeData should be non-overlapping sequential intervals of roughly equal length. These should be in relative time, so the earliest interval should be 1 and it should increase as the intervals go up with age. This is so differences in interval numbers represents the same rough difference in interval timing. For example, a dataset where all taxa are listed from a set of sequential intervals of similar length, such as North American Mammal assemblage zones, microfossil faunal zones or graptolite biozones can be given as long as they are correctly numbered in sequential order in the input. As a counter example, a dataset which includes taxa resolved only to intervals as wide as the whole Jurassic and taxa resolved to biozones within the Jurassic should not be included in the same input. Drop the taxa from less poorly resolved intervals from such datasets if you want to apply this function, as long as this retains a large enough sample of taxa from the sequential intervals. Note that taxicDivDisc and the "bin_" timescaling methods do not require that intervals be truly sequential (they can be overlapping; see their helpfiles). The output from binTimeData is always sequential, at least by default. Please check the $message and $convergence elements of the output to make sure that convergence occurred. The likelihood surface can be very flat in some cases, particularly for small datasets (

References

Foote, M. 1997 Estimating Taxonomic Durations and Preservation Probability. Paleobiology 23(3):278--300. Foote, M., and D. M. Raup. 1996 Fossil preservation and the stratigraphic ranges of taxa. Paleobiology 22(2):121--140.

See Also

freqRat, getSampRateCont, sProb2sRate, qsProb2Comp

Examples

#Simulate some fossil ranges with simFossilTaxa
set.seed(444)
taxa <- simFossilTaxa(p=0.1,q=0.1,nruns=1,mintaxa=20,maxtaxa=30,maxtime=1000,maxExtant=0)
#simulate a fossil record with imperfect sampling with sampleRanges
rangesCont <- sampleRanges(taxa,r=0.5)
#Now let's use binTimeData to bin in intervals of 1 time unit
rangesDisc <- binTimeData(rangesCont,int.length=1)

#now, get an estimate of the sampling rate (we set it to 0.5 above)
#for discrete data we can estimate the sampling probability per interval (R)
    #i.e. this is not the same thing as the instantaneous sampling rate (r)
#can use sRate2sProb to see what we would expect
sRate2sProb(r=0.5)
#expect R = ~0.39

#now we can use maximum likelihood to taxon ranges to get sampling probability
SPres1 <- getSampProbDisc(rangesDisc)
print(SPres1)      #let's look at the results
sProb<-SPres1[[2]][2]
print(sProb)
#est. R = ~0.42; not too off what we would expect!
#for the rate calibrated timescaling methods, we want an estimate of the instantaneous samp rate
#we can use sProb2sRate to get the rate. We will also need to also tell it the int.length
sRate <- sProb2sRate(sProb,int.length=1)
print(sRate)
#estimates that r=0.54... Not bad!
#Note: for real data, you may need to use an average int.length (no constant length)

#this data was simulated under homogenous sampling probabilities, extinction rates
#if we fit a model with random groups and allow for multiple timebins
	#AIC should be higher (less informative models)
randomgroup <- sample(1:2,nrow(rangesDisc[[2]]),replace=TRUE)
SPres2 <- getSampProbDisc(rangesDisc,grp1=randomgroup)
SPres3 <- getSampProbDisc(rangesDisc,n_tbins=2)
print(c(SPres1$AICc,SPres2$AICc,SPres3$AICc))
#and we can see the most simple model has the lowest AICc (most informative model)

#testing temporal change in sampling probabiluty
set.seed(444)
taxa <- simFossilTaxa(p=0.1,q=0.1,nruns=1,mintaxa=100,maxtaxa=125,maxtime=1000,maxExtant=0,
  plot=T)
#let's see what the 'true' diversity curve looks like in this case
#simulate two sets of ranges at r=0.7 and r=0.1
rangesCont <- sampleRanges(taxa,r=1.1)
rangesCont2 <- sampleRanges(taxa,r=0.2)
#now make it so that taxa which originated after 850 have r=0.1
rangesCont[taxa[,3]<850,] <- rangesCont2[taxa[,3]<850,]
rangesDisc <- binTimeData(rangesCont)
#lets plot the diversity curve
taxicDivDisc(rangesDisc)
SPres1 <- getSampProbDisc(rangesDisc)
SPres2 <- getSampProbDisc(rangesDisc,n_tbins=2)
#compare the AICc of the models
print(c(SPres1$AICc,SPres2$AICc)) #model 2 looks pretty good
#when does it find the break in time intervals?
print(rangesDisc[[1]][SPres2$interval.boundaries[2],1])
#not so great: estimates 940, not 850 
    #but look at the diversity curve: most richness in bin 1 is before 940
    #might have found the right break time otherwise...
    #the parameter values it found are less great. Finds variation in q