DDD (version 1.9)

dd_ML: Maximization of the loglikelihood under a diversity-dependent diversification model

Description

This function computes the maximum likelihood estimates of the parameters of a diversity-dependent diversification model for a given set of phylogenetic branching times. It also outputs the corresponding loglikelihood that can be used in model comparisons.

Usage

dd_ML(brts, initparsopt = if (ddmodel < 5) { c(0.1 + bd(sort(as.numeric(brts),
   decreasing = TRUE))$r1/(1 - bd(sort(as.numeric(brts),decreasing = TRUE))$a),0.1,
   2 * (length(brts) + missnumspec))} else { c(0.1 + bd(sort(as.numeric(brts),
   decreasing = TRUE))$r1/(1 - bd(sort(as.numeric(brts),decreasing = TRUE))$a), 0.1,
   2 * (length(brts) + missnumspec), 0.01)}, idparsopt = 1:length(initparsopt),
   idparsfix = (1:(3 + (ddmodel == 5)))[-idparsopt], parsfix = (ddmodel < 5) * 
   c(0.2,0.1, 2 * (length(brts) + missnumspec))[-idparsopt] + (ddmodel == 5) * 
   c(0.2, 0.1, 2 * (length(brts) + missnumspec), 0)[-idparsopt], res = 10 * (1 +
   length(brts) + missnumspec), ddmodel = 1, missnumspec = 0, cond = TRUE, btorph = 1,
   tol = c(1E-3,1E-4,1E-6), maxiter = 1000 * round((1.25)^length(idparsopt)))

Arguments

brts

A set of branching times of a phylogeny, all positive

initparsopt

The initial values of the parameters that must be optimized

idparsopt

The ids of the parameters that must be optimized, e.g. 1:3 for intrinsic speciation rate, extinction rate and carrying capacity. The ids are defined as follows: id == 1 corresponds to lambda (speciation rate) id == 2 corresponds to mu (extinction ra

idparsfix

The ids of the parameters that should not be optimized, e.g. c(1,3) if lambda and K should not be optimized, but only mu. In that case idparsopt must be 2. The default is to fix all parameters not specified in idparsopt.

parsfix

The values of the parameters that should not be optimized

res

Sets the maximum number of species for which a probability must be computed, must be larger than 1 + length(brts)

ddmodel

Sets the model of diversity-dependence: ddmodel == 1 : linear dependence in speciation rate ddmodel == 2 : exponential dependence in speciation rate ddmodel == 3 : linear dependence in extinction rate ddmodel == 4 : exponential dependence in extinctio

missnumspec

The number of species that are in the clade but missing in the phylogeny

cond

Conditioning on non-extinction of the phylogeny (TRUE or FALSE)

btorph

Sets whether the likelihood is for the branching times (0) or the phylogeny (1)

tol

Sets the tolerances in the optimization. Consists of: reltolx = relative tolerance of parameter values in optimization reltolf = relative tolerance of function value in optimization abstolx = absolute tolerance of parameter values in optimization

maxiter

Sets the maximum number of iterations in the optimization

Value

• lambda
  gives the maximum likelihood estimate of lambda
• mu
  gives the maximum likelihood estimate of mu
• K
  gives the maximum likelihood estimate of K
• r
  (only if ddmodel == 5) gives the ratio of linear dependencies in speciation and extinction rates
• loglik
  gives the maximum loglikelihood
• df
  gives the number of estimated parameters, i.e. degrees of feedom
• conv
  gives a message on convergence of optimization; conv = 0 means convergence

Details

The output is a dataframe containing estimated parameters and maximum loglikelihood. The computed loglikelihood contains the factor q! m! / (q + m)! where q is the number of species in the phylogeny and m is the number of missing species, as explained in the supplementary material to Etienne et al. 2012.

References

Etienne et al. 2011, Proc. Roy. Soc. B, doi: 10.1098/rspb.2011.1439

See Also

dd_loglik, dd_SR_ML, dd_KI_ML,

Examples

cat("Estimating parameters for a set of branching times brts with the default settings:")
brts = 1:20
dd_ML(brts)

cat("Estimating the intrinsic speciation rate lambda and the carrying capacity K")
cat("for a fixed extinction rate of 0.1:")
brts = 1:20
dd_ML(brts,initparsopt = c(0.2,100), idparsopt = c(1,3), parsfix = 0.1)

cat("This will estimate all parameters under the linear diversity-dependent")
cat("extinction model:")
brts = 1:20
dd_ML(brts, ddmodel = 3)

cat("This will estimate all parameters under the default model")
cat("when there are 2 missing species:")
brts = 1:20
dd_ML(brts, missnumspec = 2)
dd_ML(1:3,c(1.14285,0.000001,4))