make_inverseSurv(timeList, groups = NULL, p_cont = TRUE, q_cont = TRUE,
PA_n = "fixed", PB_1 = 0, Nb = 1, drop.extant = TRUE)modelMethods for manipulating and examining
such functions and constrainParPaleo for
constraining parameters.
The function output will take the largest number of
parameters possible with respect to groupings and
time-intervals, which means the number of parameters may
number in the hundreds. Constraining the function for
optimization is recommended except when datasets are very
large.
Parameters in the output functions are named 'p', 'q' and
'r', which are respectively the origination, extinction and
sampling parameters. If the respective arguments 'p_cont'
and 'q_cont' are TRUE, then 'p' and 'q' will represent the
instantaneous per-capita origination and extinction rates
(in units of per lineage time-units). When one of these
arguments is given as FALSE, the respective parameter (p or
q) will represent per-lineage-interval rates. For p, this
will be the per lineage-interval rate of a lineage
producing another lineage (which can exceed 1 because
diversity can more than double) and for q, this will be the
per lineage-interval 'rate' of a lineage going extinct,
which cannot be observed to exceed 1 (because the
proportion of diversity that goes extinct cannot exceed 1).
To obtain the per lineage-interval rates from a set of per
lineage-time unit rates, simply multiply the per
lineage-time-unit rate by the duration of an interval (or
divide, to do the reverse; see Foote, 2003 and 2005). 'r'
is always the instantaneous per-capita sampling rate, in
units per lineage-time units.
If PA_n or PB_1 were given as NULL in the arguments, two
additional parameters will be added, named respectively
'PA_n' and 'PB_1', and listed separately for every
additional grouping. These are the probability of a taxon
occurring before the first interval in the dataset (PB_1)
and the probability of a taxon occurring after the last
interval in a dataset (PA_n). Theses will be listed as
'PA_n.0' and 'PB_1.0' to indicate that they are not related
to any particular time-interval included in the analysis,
unlike the p, q, and r parameters (see below).
Groupings follow the parameter names, separated by periods;
by default, the parameters will be placed in groups
corresponding to the discrete intervals in the input
timeList, such that make_inverseSurv will create a function
with parameters 'p.1', 'q.1' and 'r.1' for interval 1;
'p.2', 'q.2' and 'r.2' for interval 2 and so on. Additional
groupings given by the user are listed after this first set
(e.g. 'p.1.2.2').
inverseSurv also contains an alternative mode which
takes input rates and returns the final values along with a
rudimentary plotting function. This allows users to output
per-interval and per-group parameter estimates. To select
these feature, the argument altMode must be TRUE.
This function will invisibly return the rate values for
each group and interval as a list of matrices, with each
matrix composed of the p, q and r rates for each interval,
for a specific grouping.
This plotting is extremely primitive, and most users will
probably find the invisibly returned rates to be of more
interest. The function layout is used to play plots
for different groupings in sequence, and this may lead to
plots which are either hard to read or even cause errors
(because of too many groupings, producing impossible
plots). To repress this, the argument plotPar can be
set to FALSE.
This capability means the function has more arguments that
just the usual par argument that accepts the vector
of parameters for running an optimization. The first of
these additional arguments, altMode enables this
alternative mode, instead of trying to estimate the
negative log-likelihood from the given parameters. The
other arguments augment the calculation and plotting of
rates.
To summarize, a function output by inverseSurv has the
following arguments:
[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object] }footeValues.
A similar format for likelihood models can be seen in
durationFreq.
Also see freqRat, sRate2sProb,
qsRate2Comp sProb2sRate and
qsProb2Comp.
For translating between sampling probabilities and sampling
rates, see SamplingConv.#let's simulate some taxon ranges from an imperfectly sampled fossil record
set.seed(444)
taxa <- simFossilTaxa(p=0.1,q=0.1,nruns=1,mintaxa=20,maxtaxa=30,maxtime=1000,maxExtant=0)
rangesCont <- sampleRanges(taxa,r=0.5)
#bin the ranges into discrete time intervals
rangesDisc <- binTimeData(rangesCont,int.length=5)
#apply make_inverseSurv
likFun<-make_inverseSurv(rangesDisc)
#use constrainParPaleo to make the model time-homogenous
#match.all~match.all will match parameters so only 2 pars: p=q and r
constrFun<-constrainParPaleo(likFun,match.all~match.all)
results <- optim(parInit(constrFun),constrFun,lower=parLower(constrFun),upper=parUpper(constrFun),
method="L-BFGS-B",control=list(maxit=1000000))
results
#plot the results
constrFun(results$par, altMode=TRUE)
#unconstrained function with ALL of 225 parameters!!!
# this will take forever to converge, so it isn't run
optim(parInit(likFun),likFun,lower=parLower(likFun),upper=parUpper(likFun),
method="L-BFGS-B",control=list(maxit=1000000))Run the code above in your browser using DataLab