R: Maximization of the loglikelihood under a diversity-dependent...
dd_SR_ML
R Documentation
Maximization of the loglikelihood under a diversity-dependent diversification model with a shift in the parameters
Description
This function computes the maximum likelihood estimates of the parameters of a
diversity-dependent diversification model with shifting parameters at time t = tshift for a given set of phylogenetic branching times.
It also outputs the corresponding loglikelihood that can be used in model comparisons.
A set of branching times of a phylogeny, all positive
initparsopt
The initial values of the parameters that must be optimized
parsfix
The values of the parameters that should not be optimized
idparsopt
The ids of the parameters that must be optimized, e.g. 1:7 for all parameters.
The ids are defined as follows:
id == 1 corresponds to lambda (speciation rate) before the shift
id == 2 corresponds to mu (extinction rate) before the shift
id == 3 corresponds to K (clade-level carrying capacity) before the shift
id == 4 corresponds to lambda (speciation rate) after the shift
id == 5 corresponds to mu (extinction rate) after the shift
id == 6 corresponds to K (clade-level carrying capacity) after the shift
id == 7 corresponds to tshift (the time of shift)
idparsfix
The ids of the parameters that should not be optimized, e.g. c(1,3,4,6) if lambda and K
should not be optimized, but only mu. In that case idparsopt must be c(2,5,7). The default is to fix all parameters
not specified in idparsopt.
idparsnoshift
The ids of the parameters that should not shift;
This can only apply to ids 4, 5 and 6, e.g. idparsnoshift = c(4,5) means that lambda and mu have the same values before and after tshift
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 == 2.1 : variant of exponential dependence in speciation rate with offset at infinity
ddmodel == 2.2 : 1/n dependence in speciation rate
ddmodel == 3 : linear dependence in extinction rate
ddmodel == 4 : exponential dependence in extinction rate
ddmodel == 4.1 : variant of exponential dependence in extinction rate with offset at infinity
ddmodel == 4.2 : 1/n dependence in extinction rate with offset at infinity
missnumspec
The number of species that are in the clade but missing in the phylogeny
cond
Conditioning:
cond == 0 : no conditioning
cond == 1 : conditioning on non-extinction of the phylogeny
cond == 2 : conditioning on non-extinction of the phylogeny and on the total number of extant taxa (including missing species)
cond == 3 : conditioning on the total number of extant taxa (including missing species)
(including missing species)
Note: cond == 3 assumes a uniform prior on stem age, as is the standard in constant-rate birth-death models, see
e.g. D. Aldous & L. Popovic 2004. Adv. Appl. Prob. 37: 1094-1115 and T. Stadler 2009. J. Theor. Biol. 261: 58-66.
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)
allbp
Sets whether a search should be done with various initial conditions, with tshift at each of the branching points (TRUE/FALSE)
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 'subplex'. Alternative is 'simplex' (default of previous versions)
methode
The method used to solve the master equation, default is 'analytical' which uses matrix exponentiation; alternatively numerical ODE solvers can be used, such as 'lsoda' or 'ode45'. These were used in the package before version 3.1.
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.
Value
lambda_1
gives the maximum likelihood estimate of lambda before the shift
mu_1
gives the maximum likelihood estimate of mu before the shift
K_1
gives the maximum likelihood estimate of K before the shift
lambda_2
gives the maximum likelihood estimate of lambda after the shift
mu_2
gives the maximum likelihood estimate of mu after the shift
K_2
gives the maximum likelihood estimate of K after the shift
t_shift
gives the time of the shift
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
Note
The optimization may get trapped in local optima. Try different starting values to search for the global optimum.
Author(s)
Rampal S. Etienne & Bart Haegeman
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
dd_SR_loglik,
dd_ML,
dd_KI_ML,
Examples
cat("This will estimate parameters for a sets of branching times brts without conditioning.\n")
cat("The tolerance of the optimization is set ridiculously high to make runtime fast.\n")
cat("In real applications, use the default or more stringent settings for tol.\n")
brts = 1:10
dd_SR_ML(brts = brts, initparsopt = c(0.4581, 1E-6, 17.69, 11.09, 8.9999), idparsopt = c(1:3,6,7),
idparsfix = NULL, parsfix = NULL, idparsnoshift = c(4,5), cond = 0,
tol = c(1E-1,1E-1,1E-1),optimmethod = 'simplex'
)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(DDD)
Loading required package: deSolve
Attaching package: 'deSolve'
The following object is masked from 'package:graphics':
matplot
Loading required package: ape
Loading required package: ade4
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/DDD/dd_SR_ML.Rd_%03d_medium.png", width=480, height=480)
> ### Name: dd_SR_ML
> ### Title: Maximization of the loglikelihood under a diversity-dependent
> ### diversification model with a shift in the parameters
> ### Aliases: dd_SR_ML
> ### Keywords: models
>
> ### ** Examples
>
> cat("This will estimate parameters for a sets of branching times brts without conditioning.\n")
This will estimate parameters for a sets of branching times brts without conditioning.
> cat("The tolerance of the optimization is set ridiculously high to make runtime fast.\n")
The tolerance of the optimization is set ridiculously high to make runtime fast.
> cat("In real applications, use the default or more stringent settings for tol.\n")
In real applications, use the default or more stringent settings for tol.
> brts = 1:10
> dd_SR_ML(brts = brts, initparsopt = c(0.4581, 1E-6, 17.69, 11.09, 8.9999), idparsopt = c(1:3,6,7),
+ idparsfix = NULL, parsfix = NULL, idparsnoshift = c(4,5), cond = 0,
+ tol = c(1E-1,1E-1,1E-1),optimmethod = 'simplex'
+ )
You are optimizing la mu K K2 tshift
You are fixing nothing
You are not shifting la2 mu2
Optimizing the likelihood - this may take a while.
The loglikelihood for the initial parameter values is -24.52893
1 0.4581 1e-06 17.69 11.09 8.9999 -24.5289252582909 initial
Optimization has terminated successfully.
Maximum likelihood parameter estimates: 0.458100 0.000001 18.574500 0.458100 0.000001 11.090000 8.999900
Maximum loglikelihood: -24.527807
>
>
>
>
>
> dev.off()
null device
1
>