Nothing
#' @export
pbd_ML = function(
brts,
initparsopt = c(0.2,0.1,1),
idparsopt = 1:length(initparsopt),
idparsfix = NULL,
parsfix = NULL,
exteq = 1,
parsfunc = c(function(t,pars) {pars[1]},function(t,pars) {pars[2]},function(t,pars) {pars[3]},function(t,pars) {pars[4]}),
missnumspec = 0,
cond = 1,
btorph = 1,
soc = 2,
methode = "lsoda",
n_low = 0,
n_up = 0,
tol = c(1E-6, 1E-6, 1E-6),
maxiter = 1000 * round((1.25)^length(idparsopt)),
optimmethod = 'subplex',
verbose = TRUE)
{
# brts = branching times (positive, from present to past)
# - max(brts) = crown age
# - min(brts) = most recent branching time
# initparsopt contains initial parameter values
# - initparsopt[1] = b (= la_1 in ER2012) = speciation initiation rate
# - initparsopt[2] = mu_1 (= mu_g in ER2012) = extinction rate of good species
# - initparsopt[3] = la_1 (= la_2 in ER2012) = speciation completion rate
# - initparsopt[4] = mu_2 (= mu_i in ER2012) = extinction rate of incipient species
# exteq = incipient species have the same (1) or different (0) extinction rate as good species
# parsfunc = functions of parameters
# res = resolution of the method; res should be larger than the total number of species
# missnumspec = number of missing species
# cond = conditioning
# . cond = 0 conditioning on stem or clade age
# . cond = 1 conditioning on age and non-extinction of the phylogeny
# . cond = 2 conditioning on age and on number of extant taxa
# btorph = likelihood of branching times (0) or phylogeny (1), differ by a factor (S - 1)! where S is the number of extant species
# soc = stem (1) or crown (2) age
# methode = method of the numerical integration; see package deSolve for details
# n_low = lower bound on number of species (cond = 2)
# n_up = upper bound on number of species (cond = 2)
# tol = tolerance in optimization
# - 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 = the maximum number of iterations in the optimization
# optimmethod = 'subplex' (current default) or 'simplex' (default of previous versions)
# verbose = whether intermediate output should be shown
options(warn=-1)
brts = sort(abs(as.numeric(brts)),decreasing = TRUE)
if(is.numeric(brts) == FALSE)
{
cat("The branching times should be numeric.\n")
out2 = data.frame(b = -1, mu_1 = -1, lambda_1 = -1,mu_2 = -1, loglik = -1, df = -1, conv = -1)
} else {
if(exteq == 1){ idexteq = 4 } else { idexteq = NULL }
idpars = sort(c(idparsopt,idparsfix,idexteq))
if((prod(idpars == (1:4)) != 1) || (length(initparsopt) != length(idparsopt)) || (length(parsfix) != length(idparsfix)))
{
cat("The arguments should be coherent.\n")
out2 = data.frame(b = -1, mu_1 = -1, lambda_1 = -1,mu_2 = -1, loglik = -1, df = -1, conv = -1)
} else {
namepars = c("b", "mu_1", "lambda_1", "mu_2")
if(length(namepars[idparsopt]) == 0) { optstr = "nothing" } else { optstr = namepars[idparsopt] }
if (verbose) { cat("You are optimizing",optstr,"\n") }
if(length(namepars[idparsfix]) == 0) { fixstr = "nothing" } else { fixstr = namepars[idparsfix] }
if (verbose) { cat("You are fixing",fixstr,"\n") }
if(exteq == 1) { fixstr = "exactly" } else { fixstr = "not" }
if (verbose) { cat("Extinction rate of incipient species is",fixstr,"the same as for good species.\n") }
trparsopt = initparsopt/(1 + initparsopt)
trparsfix = parsfix/(1 + parsfix)
trparsfix[parsfix == Inf] = 1
pars2 = c(cond,btorph,soc,0,methode,n_low,n_up)
utils::flush.console()
initloglik = pbd_loglik_choosepar(trparsopt = trparsopt,trparsfix = trparsfix,idparsopt = idparsopt,idparsfix = idparsfix,exteq = exteq,parsfunc = parsfunc,pars2 = pars2,brts = brts,missnumspec = missnumspec)
if (verbose) { cat("The likelihood for the initial parameter values is",initloglik,"\n") }
utils::flush.console()
if(initloglik == -Inf)
{
cat("The initial parameter values have a likelihood that is equal to 0 or below machine precision. Try again with different initial values.\n")
out2 = data.frame(b = -1, mu_1 = -1, lambda_1 = -1,mu_2 = -1, loglik = -1, df = -1, conv = -1)
} else {
if (verbose) { cat("Optimizing the likelihood - this may take a while.","\n") }
utils::flush.console()
optimpars = c(tol,maxiter)
out = DDD::optimizer(optimmethod = optimmethod,optimpars = optimpars,fun = pbd_loglik_choosepar,trparsopt = trparsopt,trparsfix = trparsfix,idparsopt = idparsopt,idparsfix = idparsfix,exteq = exteq, parsfunc = parsfunc, pars2 = pars2,brts = brts, missnumspec = missnumspec)
if(out$conv > 0)
{
cat("Optimization has not converged. Try again with different initial values.\n")
out2 = data.frame(b = -1, mu_1 = -1, lambda_1 = -1,mu_2 = -1, loglik = -1, df = -1, conv = -1)
} else {
MLtrpars = as.numeric(unlist(out$par))
MLpars = MLtrpars/(1-MLtrpars)
MLpars1 = rep(0,4)
MLpars1[idparsopt] = MLpars
if(length(idparsfix) != 0) { MLpars1[idparsfix] = parsfix }
if(exteq == 1) { MLpars1[4] = MLpars[2] }
if(MLpars1[3] > 10^7){MLpars1[3] = Inf}
ML = as.numeric(unlist(out$fvalues))
out2 = data.frame(b = MLpars1[1],mu_1 = MLpars1[2],lambda_1 = MLpars1[3], mu_2 = MLpars1[4], loglik = ML, df = length(initparsopt), conv = unlist(out$conv))
if(verbose)
{
s1 = sprintf('Maximum likelihood parameter estimates: b: %f, mu_1: %f, lambda_1: %f, mu_2: %f',MLpars1[1],MLpars1[2],MLpars1[3],MLpars1[4])
s2 = sprintf('Maximum loglikelihood: %f',ML)
s3 = sprintf('The expected duration of speciation for these parameters is: %f',pbd_durspec_mean(c(MLpars1[1],MLpars1[3],MLpars1[4])))
s4 = sprintf('The median duration of speciation for these parameters is: %f',pbd_durspec_quantile(c(MLpars1[1],MLpars1[3],MLpars1[4]),0.5))
cat("\n",s1,"\n",s2,"\n",s3,"\n",s4,"\n")
utils::flush.console()
}
}
}
}
}
invisible(out2)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.