Maximization of the loglikelihood under the diversity-independent, possibly time-dependent diversification model

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

Usage

bd_ML(
  brts,
  initparsopt = c(0.1, 0.05 * (tdmodel <= 1) + 10 * (length(brts) + missnumspec) *
    (tdmodel > 1)),
  idparsopt = c(1, 2 + (tdmodel > 1)),
  idparsfix = (1:4)[-idparsopt],
  parsfix = rep(0, 4)[idparsfix],
  missnumspec = 0,
  tdmodel = 0,
  cond = 1,
  btorph = 1,
  soc = 2,
  tol = c(0.001, 1e-04, 1e-06),
  maxiter = 1000 * round((1.25)^length(idparsopt)),
  changeloglikifnoconv = FALSE,
  optimmethod = "simplex",
  num_cycles = 1,
  methode = "odeint::runge_kutta_cash_karp54",
  verbose = FALSE
)

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 lambda0 (speciation rate)
id == 2 corresponds to mu0 (extinction rate)
id == 3 corresponds to lamda1 (parameter controlling decline in speciation rate with time)
id == 4 corresponds to mu1 (parameter controlling decline in extinction rate with time)

idparsfix

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

parsfix

The values of the parameters that should not be optimized

missnumspec

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

tdmodel

Sets the model of time-dependence:
tdmodel == 0 : constant speciation and extinction rates
tdmodel == 1 : speciation and/or extinction rate is exponentially declining with time
tdmodel == 2 : stepwise decline in speciation rate as in diversity-dependence without extinction
tdmodel == 3 : decline in speciation rate following deterministic logistic equation for ddmodel = 1
tdmodel == 4 : decline in speciation rate such that the expected number of species matches with that of ddmodel = 1 with the same mu

cond

Conditioning:
cond == 0 : conditioning on stem or crown age
cond == 1 : conditioning on stem or crown age and non-extinction of the phylogeny
cond == 2 : conditioning on stem or crown age and on the total number of extant taxa (including missing species)
cond == 3 : conditioning on the total number of extant taxa (including missing species)

btorph

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

soc

Sets whether stem or crown age should be used (1 or 2)

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

changeloglikifnoconv

if TRUE the loglik will be set to -Inf if ML does not converge

optimmethod

Method used in optimization of the likelihood. Current default is 'simplex'. Alternative is 'subplex' (default of previous versions)

num_cycles

the number of cycles of opimization. If set at Inf, it will do as many cycles as needed to meet the tolerance set for the target function.

methode

The method used to solve the master equation under tdmodel = 4, default is 'odeint::runge_kutta_cash_karp54'.

verbose

Show the parameters and loglikelihood for every call to the loglik function

Value

A dataframe with the following elements:

lambda0

gives the maximum likelihood estimate of lambda0

mu0

gives the maximum likelihood estimate of mu0

lambda1

gives the maximum likelihood estimate of lambda1

mu1

gives the maximum likelihood estimate of mu1

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, R.S. et al. 2012, Proc. Roy. Soc. B 279: 1300-1309, doi: 10.1098/rspb.2011.1439
- Etienne, R.S. & B. Haegeman 2012. Am. Nat. 180: E75-E89, doi: 10.1086/667574

See also

bd_loglik

Author

Rampal S. Etienne & Bart Haegeman

Examples

cat("Estimating parameters for a set of branching times brts with the default settings:")
#> Estimating parameters for a set of branching times brts with the default settings:
brts = 1:20
bd_ML(brts = brts, cond = 1)
#> You are optimizing lambda0 mu0 
#> You are fixing lambda1 mu1 
#> Optimizing the likelihood - this may take a while. 
#> The loglikelihood for the initial parameter values is -67.8956450330451.
#> 1 0.1 0.05 -67.8956450330451 initial 
#> 2 0.10751708428246 0.0450354609929078 -67.6692013580063 expand 
#> 3 0.10751708428246 0.0450354609929078 -67.6692013580063 reflect 
#> 4 0.115137614678899 0.0401176470588235 -67.6605181277165 reflect 
#> 5 0.108779931584949 0.0376760563380281 -67.4968509240751 expand 
#> 6 0.108779931584949 0.0376760563380281 -67.4968509240751 reflect 
#> 7 0.10751708428246 0.0256380510440834 -67.2533235698667 expand 
#> 8 0.0919708029197083 0.029220023282887 -67.2328548115931 expand 
#> 9 0.081947148817803 0.00752136752136739 -66.6092346077597 expand 
#> 10 0.081947148817803 0.00752136752136739 -66.6092346077597 contract inside 
#> 11 0.0870118074477753 0.000537558353373736 -66.4167519746098 reflect 
#> 12 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 13 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 14 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 15 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 16 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 17 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 18 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 19 0.0870118074477753 0.000537558353373736 -66.4167519746098 contract inside 
#> 20 0.0869275940813172 9.90616347764255e-05 -66.4063799769886 reflect 
#> 21 0.0869275940813172 9.90616347764255e-05 -66.4063799769886 contract inside 
#> 22 0.0869275940813172 9.90616347764255e-05 -66.4063799769886 contract inside 
#> 23 0.0869275940813172 9.90616347764255e-05 -66.4063799769886 contract inside 
#> 24 0.0869275940813172 9.90616347764255e-05 -66.4063799769886 contract inside 
#> 25 0.0868194993048274 4.6451145110446e-05 -66.4040401804145 reflect 
#> 26 0.0868194993048274 4.6451145110446e-05 -66.4040401804145 contract inside 
#> 27 0.0868194993048274 4.6451145110446e-05 -66.4040401804145 contract inside 
#> 28 0.0868194993048274 4.6451145110446e-05 -66.4040401804145 contract inside 
#> 29 0.0866900624543031 5.79331059251881e-05 -66.4028740830437 expand 
#> 30 0.0866900624543031 5.79331059251881e-05 -66.4028740830437 contract inside 
#> 31 0.0866935194470033 1.61643430795484e-05 -66.4020061670072 reflect 
#> 32 0.0864364642522744 1.82431053102008e-05 -66.3993551222831 expand 
#> 33 0.0864364642522744 1.82431053102008e-05 -66.3993551222831 contract inside 
#> 34 0.0862090502482492 5.1388518576471e-05 -66.3978487512211 expand 
#> 35 0.0857137167375939 2.93106577196041e-05 -66.3929321539514 expand 
#> 36 0.0850122982006982 8.45651204878309e-05 -66.3890717993446 expand 
#> 37 0.083674530718775 6.8034525211653e-05 -66.3824739534804 expand 
#> 38 0.083674530718775 6.8034525211653e-05 -66.3824739534804 reflect 
#> 39 0.0824833540059418 0.00010121223160964 -66.3819361811278 contract outside 
#> 40 0.0832844158481555 7.29174544438885e-06 -66.3799884852552 expand 
#> 41 0.0832844158481555 7.29174544438885e-06 -66.3799884852552 reflect 
#> 42 0.0832844158481555 7.29174544438885e-06 -66.3799884852552 contract inside 
#> 43 0.0832844158481555 7.29174544438885e-06 -66.3799884852552 contract inside 
#> 44 0.0832844158481555 7.29174544438885e-06 -66.3799884852552 contract inside 
#> 45 0.0832844158481555 7.29174544438885e-06 -66.3799884852552 contract inside 
#> 46 0.0827775751121175 2.61625095469112e-05 -66.379905383758 contract inside 
#> 47 0.0832903834742893 1.81194330916335e-06 -66.3798568412851 reflect 
#> 48 0.0827835371553468 2.06825005957022e-05 -66.3797615164625 reflect 
#> 49 0.0827835371553468 2.06825005957022e-05 -66.3797615164625 contract inside 
#> 50 0.0830678341839042 1.07527077822123e-05 -66.3797419562324 contract inside 
#> 51 0.0829625578401348 9.74326091209499e-06 -66.3795981441118 expand 
#> 52 0.0829625578401348 9.74326091209499e-06 -66.3795981441118 contract inside 
#> 53 0.0829625578401348 9.74326091209499e-06 -66.3795981441118 reflect 
#> 54 0.0828362777077475 5.36815054069265e-06 -66.379381300286 expand 
#> 55 0.0828362777077475 5.36815054069265e-06 -66.379381300286 reflect 
#> 56 0.0828362777077475 5.36815054069265e-06 -66.379381300286 contract inside 
#> 57 0.0828362777077475 5.36815054069265e-06 -66.379381300286 contract inside 
#> 58 0.0829099885209629 1.32912255536243e-06 -66.3793276215614 reflect 
#> 59 0.0829099885209629 1.32912255536243e-06 -66.3793276215614 reflect 
#> 60 0.0828047428592901 3.19694710158561e-07 -66.3792282185678 expand 
#> 61 0.0828047428592901 3.19694710158561e-07 -66.3792282185678 contract inside 
#> 62 0.0828047428592901 3.19694710158561e-07 -66.3792282185678 reflect 
#> 63 0.0828047428592901 3.19694710158561e-07 -66.3792282185678 contract inside 
#> 64 0.0827245006717876 3.40940471567272e-07 -66.3791940531816 reflect 
#> 65 0.0827245006717876 3.40940471567272e-07 -66.3791940531816 contract inside 
#> 66 0.0827245006717876 3.40940471567272e-07 -66.3791940531816 reflect 
#> 67 0.082644270376238 3.62186233878748e-07 -66.3791777703448 reflect 
#> 68 0.082644270376238 3.62186233878748e-07 -66.3791777703448 contract inside 
#> 69 0.082644270376238 3.62186233878748e-07 -66.3791777703448 contract inside 
#> 70 0.082644270376238 3.62186233878748e-07 -66.3791777703448 reflect 
#> 71 0.082644270376238 3.62186233878748e-07 -66.3791777703448 reflect 
#> 72 0.082612956199224 1.11157429055683e-07 -66.3791692856385 reflect 
#> 73 0.082612956199224 1.11157429055683e-07 -66.3791692856385 contract outside 
#> 74 0.082612956199224 1.11157429055683e-07 -66.3791692856385 contract inside 
#> 75 0.082612956199224 1.11157429055683e-07 -66.3791692856385 reflect 
#> 76 0.082612956199224 1.11157429055683e-07 -66.3791692856385 reflect 
#> 77 0.082612956199224 1.11157429055683e-07 -66.3791692856385 contract inside 
#> 78 0.0825871516817185 2.88498843518927e-08 -66.3791676899316 reflect 
#> 79 0.0825871516817185 2.88498843518927e-08 -66.3791676899316 contract outside 
#> 80 0.0825871516817185 2.88498843518927e-08 -66.3791676899316 contract inside 
#> 81 0.0825970640107624 2.26791415952077e-08 -66.3791670652764 reflect 
#> 82 0.0825970640107624 2.26791415952077e-08 -66.3791670652764 contract inside 
#> 83 0.0825970640107624 2.26791415952077e-08 -66.3791670652764 reflect 
#> 84 0.0826069765213256 1.65083989146789e-08 -66.3791667144197 reflect 
#> 85 0.0826069765213256 1.65083989146789e-08 -66.3791667144197 contract inside 
#> 86 0.0825956544058777 7.16141905750347e-09 -66.3791666952829 reflect 
#> 87 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 reflect 
#> 88 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 contract inside 
#> 89 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 reflect 
#> 90 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 contract inside 
#> 91 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 contract inside 
#> 92 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 contract outside 
#> 93 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 contract inside 
#> 94 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 reflect 
#> 95 0.0826055668906274 9.90676568486349e-10 -66.3791663057083 contract inside 
#> 96 0.0826073456637157 6.99922452149296e-10 -66.3791662867791 reflect 
#> 97 0.0826073456637157 6.99922452149296e-10 -66.3791662867791 contract inside 
#> 98 0.0826073456637157 6.99922452149296e-10 -66.3791662867791 reflect 
#> 99 0.0826091244426492 4.09168335981318e-10 -66.3791662766609 reflect 
#> 100 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 reflect 
#> 101 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 102 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 103 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 104 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 105 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 106 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 107 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 108 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 109 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 110 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 111 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 112 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 113 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 114 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 115 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 116 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 117 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 118 0.0826072351600045 4.18852071671142e-12 -66.3791662684607 contract inside 
#> 119 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 reflect 
#> 120 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 121 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 122 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 123 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 124 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 125 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 126 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 127 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 128 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 129 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 130 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 131 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 132 0.0826072217231241 1.56066448197112e-13 -66.3791662684069 contract inside 
#> 133 0.0826072212612932 1.67128795706832e-14 -66.379166268405 reflect 
#> 134 0.0826072212612932 1.67128795706832e-14 -66.379166268405 contract inside 
#> 135 0.0826072212612932 1.67128795706832e-14 -66.379166268405 contract inside 
#> 136 0.0826072212612932 1.67128795706832e-14 -66.379166268405 contract inside 
#> 137 0.0826072212612932 1.67128795706832e-14 -66.379166268405 contract inside 
#> 138 0.0826072212612932 1.67128795706832e-14 -66.379166268405 contract inside 
#> 139 0.0826072212612932 1.67128795706832e-14 -66.379166268405 contract inside 
#> 140 0.0826072212129633 2.10471625023496e-15 -66.3791662684048 reflect 
#> 141 0.0826072212129633 2.10471625023496e-15 -66.3791662684048 contract inside 
#> 142 0.0826072212129633 2.10471625023496e-15 -66.3791662684048 contract inside 
#> 143 0.0826072212129633 2.10471625023496e-15 -66.3791662684048 contract inside 
#> 144 0.0826072212129633 2.10471625023496e-15 -66.3791662684048 contract inside 
#> 145 0.082607221207963 5.81399751544815e-16 -66.3791662684048 reflect 
#> 146 0.082607221207963 5.81399751544815e-16 -66.3791662684048 contract inside 
#> 147 0.082607221207963 5.81399751544815e-16 -66.3791662684048 contract inside 
#> 148 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 149 0.082607221207963 5.81399751544815e-16 -66.3791662684048 contract inside 
#> 150 0.082607221207963 5.81399751544815e-16 -66.3791662684048 contract outside 
#> 151 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 152 0.082607221207963 5.81399751544815e-16 -66.3791662684048 contract outside 
#> 153 0.082607221207963 5.81399751544815e-16 -66.3791662684048 contract outside 
#> 154 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 155 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 156 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 157 0.082607221207963 5.81399751544815e-16 -66.3791662684048 contract inside 
#> 158 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 159 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 160 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 161 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> 162 0.082607221207963 5.81399751544815e-16 -66.3791662684048 shrink 
#> Optimization has terminated successfully. 
#> 
#> Maximum likelihood parameter estimates: lambda0: 0.082607, mu0: 0.000000, lambda1: 0.000000, mu1: 0.000000: 
#> Maximum loglikelihood: -66.379166
#>      lambda0          mu0 lambda1 mu1    loglik df conv
#> 1 0.08260722 5.813998e-16       0   0 -66.37917  2    0