#' @keywords internal
master_ml <- function(phy,
traits,
num_concealed_states,
idparslist,
idparsopt,
initparsopt,
idparsfix,
parsfix,
idfactorsopt = NULL,
initfactors = NULL,
idparsfuncdefpar = NULL,
functions_defining_params = NULL,
cond = "proper_cond",
root_state_weight = "proper_weights",
sampling_fraction,
tol = c(1e-04, 1e-05, 1e-07),
maxiter = 1000 * round((1.25)^length(idparsopt)),
optimmethod = "subplex",
num_cycles = 1,
loglik_penalty = 0,
is_complete_tree = FALSE,
verbose = (optimmethod == "simplex"),
num_threads = 1,
atol = 1e-8,
rtol = 1e-7,
method = "odeint::bulirsch_stoer") {
structure_func <- NULL
if (!is.null(functions_defining_params)) {
structure_func <- set_and_check_structure_func(idparsfuncdefpar,
functions_defining_params,
idparslist,
idparsopt,
idfactorsopt,
idparsfix,
initfactors)
} else {
if (identical(as.numeric(sort(c(idparsopt, idparsfix))),
as.numeric(sort(unique(unlist(idparslist))))) == FALSE) {
stop("All elements in idparslist must be included in either
idparsopt or idparsfix ")
}
}
check_ml_conditions(traits,
idparslist,
initparsopt,
idparsopt,
idparsfix,
parsfix)
if (is.matrix(idparslist[[1]])) {
## it is a tailor case otherwise
idparslist[[1]] <- prepare_full_lambdas(traits,
num_concealed_states,
idparslist[[1]])
}
see_ancestral_states <- FALSE
if (!is.null(structure_func)) {
initparsopt <- c(initparsopt, initfactors)
}
trparsopt <- initparsopt / (1 + initparsopt)
trparsopt[which(initparsopt == Inf)] <- 1
trparsfix <- parsfix / (1 + parsfix)
trparsfix[which(parsfix == Inf)] <- 1
mus <- calc_mus(is_complete_tree,
idparslist,
idparsfix,
parsfix,
idparsopt,
initparsopt)
optimpars <- c(tol, maxiter)
num_modeled_traits <- length(idparslist[[1]]) / num_concealed_states
setting_calculation <- build_initStates_time(phy,
traits,
num_concealed_states,
sampling_fraction,
is_complete_tree,
mus,
num_modeled_traits,
traitStates =
get_trait_states(idparslist,
num_concealed_states))
initloglik <- secsse_loglik_choosepar(trparsopt = trparsopt,
trparsfix = trparsfix,
idparsopt = idparsopt,
idparsfix = idparsfix,
idparslist = idparslist,
structure_func = structure_func,
phy = phy,
traits = traits,
num_concealed_states =
num_concealed_states,
cond = cond,
root_state_weight = root_state_weight,
sampling_fraction = sampling_fraction,
setting_calculation =
setting_calculation,
see_ancestral_states =
see_ancestral_states,
loglik_penalty = loglik_penalty,
is_complete_tree = is_complete_tree,
verbose = verbose,
num_threads = num_threads,
atol = atol,
rtol = rtol,
method = method)
# Function here
print_init_ll(initloglik = initloglik, verbose = verbose)
if (initloglik == -Inf) {
stop("The initial parameter values have a likelihood that is
equal to 0 or below machine precision.
Try again with different initial values.")
} else {
out <- DDD::optimizer(optimmethod = optimmethod,
optimpars = optimpars,
fun = secsse_loglik_choosepar,
trparsopt = trparsopt,
idparsopt = idparsopt,
trparsfix = trparsfix,
idparsfix = idparsfix,
idparslist = idparslist,
structure_func = structure_func,
phy = phy,
traits = traits,
num_concealed_states = num_concealed_states,
cond = cond,
root_state_weight = root_state_weight,
sampling_fraction = sampling_fraction,
setting_calculation = setting_calculation,
see_ancestral_states = see_ancestral_states,
num_cycles = num_cycles,
loglik_penalty = loglik_penalty,
is_complete_tree = is_complete_tree,
verbose = verbose,
num_threads = num_threads,
atol = atol,
rtol = rtol,
method = method)
if (out$conv != 0) {
stop("Optimization has not converged.
Try again with different initial values.")
} else {
ml_pars1 <- secsse_transform_parameters(as.numeric(unlist(out$par)),
trparsfix,
idparsopt,
idparsfix,
idparslist,
structure_func)
out2 <- list(MLpars = ml_pars1,
ML = as.numeric(unlist(out$fvalues)),
conv = out$conv)
}
}
return(out2)
}
#' Maximum likehood estimation for (SecSSE)
#'
#' Maximum likehood estimation under Several examined and concealed
#' States-dependent Speciation and Extinction (SecSSE)
#' @inheritParams default_params_doc
#'
#' @return Parameter estimated and maximum likelihood
#' @examples
#'# Example of how to set the arguments for a ML search.
#'library(secsse)
#'library(DDD)
#'set.seed(13)
#'# lambdas for 0A and 1A and 2A are the same but need to be estimated
#'# mus are fixed to
#'# the transition rates are constrained to be equal and fixed 0.01
#'phylotree <- ape::rcoal(31, tip.label = 1:31)
#'traits <- sample(c(0,1,2), ape::Ntip(phylotree),replace=TRUE)#get some traits
#'num_concealed_states<-3
#'idparslist <- id_paramPos(traits, num_concealed_states)
#'idparslist[[1]][c(1,4,7)] <- 1
#'idparslist[[1]][c(2,5,8)] <- 2
#'idparslist[[1]][c(3,6,9)] <- 3
#'idparslist[[2]][]<-4
#'masterBlock <- matrix(5,ncol = 3,nrow = 3,byrow = TRUE)
#'diag(masterBlock) <- NA
#'diff.conceal <- FALSE
#'idparslist[[3]] <- q_doubletrans(traits,masterBlock,diff.conceal)
#'startingpoint <- DDD::bd_ML(brts = ape::branching.times(phylotree))
#'intGuessLamba <- startingpoint$lambda0
#'intGuessMu <- startingpoint$mu0
#'idparsopt <- c(1,2,3,5)
#'initparsopt <- c(rep(intGuessLamba,3),rep((intGuessLamba/5),1))
#'idparsfix <- c(0,4)
#'parsfix <- c(0,0)
#'tol <- c(1e-02, 1e-03, 1e-04)
#'maxiter <- 1000 * round((1.25)^length(idparsopt))
#'optimmethod <- 'subplex'
#'cond <- 'proper_cond'
#'root_state_weight <- 'proper_weights'
#'sampling_fraction <- c(1,1,1)
#'model<-secsse_ml(
#'phylotree,
#'traits,
#'num_concealed_states,
#'idparslist,
#'idparsopt,
#'initparsopt,
#'idparsfix,
#'parsfix,
#'cond,
#'root_state_weight,
#'sampling_fraction,
#'tol,
#'maxiter,
#'optimmethod,
#'num_cycles = 1,
#'verbose = FALSE)
#'# model$ML
#'# [1] -16.04127
#' @export
secsse_ml <- function(phy,
traits,
num_concealed_states,
idparslist,
idparsopt,
initparsopt,
idparsfix,
parsfix,
cond = "proper_cond",
root_state_weight = "proper_weights",
sampling_fraction,
tol = c(1e-04, 1e-05, 1e-07),
maxiter = 1000 * round((1.25)^length(idparsopt)),
optimmethod = "subplex",
num_cycles = 1,
loglik_penalty = 0,
is_complete_tree = FALSE,
verbose = (optimmethod == "simplex"),
num_threads = 1,
atol = 1e-8,
rtol = 1e-7,
method = "odeint::bulirsch_stoer") {
master_ml(phy = phy,
traits = traits,
num_concealed_states = num_concealed_states,
idparslist = idparslist,
idparsopt = idparsopt,
initparsopt = initparsopt,
idparsfix = idparsfix,
parsfix = parsfix,
initfactors = NULL,
idparsfuncdefpar = NULL,
cond = cond,
root_state_weight = root_state_weight,
sampling_fraction = sampling_fraction,
tol = tol,
maxiter = maxiter,
optimmethod = optimmethod,
num_cycles = num_cycles,
loglik_penalty = loglik_penalty,
is_complete_tree = is_complete_tree,
verbose = verbose,
num_threads = num_threads,
atol = atol,
rtol = rtol,
method = method)
}
#' @keywords internal
secsse_loglik_choosepar <- function(trparsopt,
trparsfix,
idparsopt,
idparsfix,
idparslist,
structure_func = structure_func,
phy = phy,
traits = traits,
num_concealed_states = num_concealed_states,
cond = cond,
root_state_weight = root_state_weight,
sampling_fraction = sampling_fraction,
setting_calculation = setting_calculation,
see_ancestral_states = see_ancestral_states,
loglik_penalty = loglik_penalty,
is_complete_tree = is_complete_tree,
verbose = verbose,
num_threads = num_threads,
atol = atol,
rtol = rtol,
method = method) {
alltrpars <- c(trparsopt, trparsfix)
if (max(alltrpars) > 1 || min(alltrpars) < 0) {
loglik <- -Inf
} else {
pars1 <- secsse_transform_parameters(trparsopt, trparsfix,
idparsopt, idparsfix,
idparslist, structure_func)
loglik <- master_loglik(parameter = pars1,
phy = phy,
traits = traits,
num_concealed_states =
num_concealed_states,
cond = cond,
root_state_weight =
root_state_weight,
sampling_fraction =
sampling_fraction,
setting_calculation =
setting_calculation,
see_ancestral_states =
see_ancestral_states,
loglik_penalty = loglik_penalty,
is_complete_tree =
is_complete_tree,
num_threads = num_threads,
method = method,
atol = atol,
rtol = rtol)
if (is.nan(loglik) || is.na(loglik)) {
warning("There are parameter values used which cause
numerical problems.")
loglik <- -Inf
}
}
if (verbose) {
out_print <- c(trparsopt / (1 - trparsopt), loglik)
message(paste(out_print, collapse = " "))
}
return(loglik)
}
#' Maximum likehood estimation for (SecSSE)
#'
#' Maximum likehood estimation under Several examined and concealed
#' States-dependent Speciation and Extinction (SecSSE) with cladogenetic option
#'
#' @inheritParams default_params_doc
#'
#' @return Parameter estimated and maximum likelihood
#' @examples
#'# Example of how to set the arguments for a ML search.
#'library(secsse)
#'library(DDD)
#'set.seed(13)
#'# Check the vignette for a better working exercise.
#'# lambdas for 0A and 1A and 2A are the same but need to be estimated
#'# (CTD model, see Syst Biol paper)
#'# mus are fixed to zero,
#'# the transition rates are constrained to be equal and fixed 0.01
#'phylotree <- ape::rcoal(31, tip.label = 1:31)
#'#get some traits
#'traits <- sample(c(0,1,2), ape::Ntip(phylotree), replace = TRUE)
#'num_concealed_states <- 3
#'idparslist <- cla_id_paramPos(traits,num_concealed_states)
#'idparslist$lambdas[1,] <- c(1,1,1,2,2,2,3,3,3)
#'idparslist[[2]][] <- 4
#'masterBlock <- matrix(5,ncol = 3,nrow = 3,byrow = TRUE)
#'diag(masterBlock) <- NA
#'diff.conceal <- FALSE
#'idparslist[[3]] <- q_doubletrans(traits,masterBlock,diff.conceal)
#'startingpoint <- bd_ML(brts = ape::branching.times(phylotree))
#'intGuessLamba <- startingpoint$lambda0
#'intGuessMu <- startingpoint$mu0
#'idparsopt <- c(1,2,3)
#'initparsopt <- c(rep(intGuessLamba,3))
#'idparsfix <- c(0,4,5)
#'parsfix <- c(0,0,0.01)
#'tol <- c(1e-04, 1e-05, 1e-07)
#'maxiter <- 1000 * round((1.25) ^ length(idparsopt))
#'optimmethod <- 'subplex'
#'cond <- 'proper_cond'
#'root_state_weight <- 'proper_weights'
#'sampling_fraction <- c(1,1,1)
#'model <- cla_secsse_ml(
#' phylotree,
#' traits,
#' num_concealed_states,
#' idparslist,
#' idparsopt,
#' initparsopt,
#' idparsfix,
#' parsfix,
#' cond,
#' root_state_weight,
#' sampling_fraction,
#' tol,
#' maxiter,
#' optimmethod,
#' num_cycles = 1,
#' num_threads = 1,
#' verbose = FALSE)
#' # [1] -90.97626
#' @export
cla_secsse_ml <- function(phy,
traits,
num_concealed_states,
idparslist,
idparsopt,
initparsopt,
idparsfix,
parsfix,
cond = "proper_cond",
root_state_weight = "proper_weights",
sampling_fraction,
tol = c(1e-04, 1e-05, 1e-07),
maxiter = 1000 * round((1.25)^length(idparsopt)),
optimmethod = "subplex",
num_cycles = 1,
loglik_penalty = 0,
is_complete_tree = FALSE,
verbose = (optimmethod == "simplex"),
num_threads = 1,
atol = 1e-8,
rtol = 1e-7,
method = "odeint::bulirsch_stoer") {
master_ml(phy = phy,
traits = traits,
num_concealed_states = num_concealed_states,
idparslist = idparslist,
idparsopt = idparsopt,
initparsopt = initparsopt,
idparsfix = idparsfix,
parsfix = parsfix,
cond = cond,
root_state_weight = root_state_weight,
sampling_fraction = sampling_fraction,
tol = tol,
maxiter = maxiter,
optimmethod = optimmethod,
num_cycles = num_cycles,
loglik_penalty = loglik_penalty,
is_complete_tree = is_complete_tree,
verbose = verbose,
num_threads = num_threads,
atol = atol,
rtol = rtol,
method = method)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.