Skip to content

Commit

Permalink
Merge pull request #25 from tece-lab/develop
Browse files Browse the repository at this point in the history
Adding 2-type and relaxed-rate functionality
  • Loading branch information
Neves-P authored Mar 1, 2023
2 parents 4ced281 + bc6a8fe commit 7b822b7
Show file tree
Hide file tree
Showing 28 changed files with 1,565 additions and 32 deletions.
2 changes: 0 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -50,5 +50,3 @@ Imports:
DAISIE (>= 4.2.1),
ggplot2,
cowplot
Remotes:
rsetienne/DAISIE
2 changes: 2 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,11 @@ export(is_on_cluster)
export(plot_bootstrap_results)
export(print_metadata)
export(read_model_results)
export(run_daisie_2type_ml)
export(run_daisie_ml)
export(run_sim)
export(sensitivity)
export(setup_2type_model)
export(setup_model)
export(setup_std_pars2)
export(summarize_bootstrap_results)
9 changes: 8 additions & 1 deletion R/calc_bic.R
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,14 @@
#' @author Joshua W. Lambert, Pedro Santos Neves
calc_bic <- function(results, daisie_data) {
k <- results$df
m <- daisie_data[[1]]$not_present + length(daisie_data) - 1
if (is.null(daisie_data[[1]]$not_present)) {
m <- daisie_data[[1]]$not_present_type1 +
daisie_data[[1]]$not_present_type2 + length(daisie_data) - 1
} else {
m <- daisie_data[[1]]$not_present + length(daisie_data) - 1
}
bic <- k * (log(m) + log(2 * pi)) - 2 * results$loglik
testit::assert(is.numeric(bic))
testit::assert(identical(length(bic), 1L))
bic
}
26 changes: 25 additions & 1 deletion R/default_params_doc.R
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,10 @@
#' extinction. All other parameters free.
#' * `"rr_k"` Clade specific model - diversity dependent and relaxed carrying
#' capacity. All other parameters free.
#' * `"rr_gam_di"` Clade specific model - diversity independent (K = Inf)
#' and relaxed colonisation. All other parameters free.
#' * `"rr_gam_dd"` Clade specific model - diversity dependent and relaxed
#' colonisation. All other parameters free.
#' * `"rr_laa_di"` Clade specific model - diversity independent (K = Inf)
#' and relaxed anagenesis. All other parameters free.
#' * `"rr_laa_dd"` Clade specific model - diversity dependent and relaxed
Expand All @@ -36,6 +40,12 @@
#' * `"rr_k_0lac"` Clade specific model - diversity dependent, relaxed
#' carrying capacity, and no cladogenesis (lac fixed to zero). All other
#' parameters free.
#' * `"rr_gam_di_0lac"` Clade specific model - diversity independent,
#' relaxed colonisation, and no cladogenesis (lac fixed to zero). All other
#' parameters free.
#' * `"r_gam_dd_0lac"` Clade specific model - diversity dependent, relaxed
#' colonisation, and no cladogenesis (lac fixed to zero). All other
#' parameters free.
#' * `"rr_laa_di_0lac"` Clade specific model - diversity independent,
#' relaxed anagenesis, and no cladogenesis (lac fixed to zero). All other
#' parameters free.
Expand All @@ -48,6 +58,12 @@
#' * `"rr_mu_di_0laa"` Clade specific model - diversity independent
#' * `"rr_mu_dd_0laa"` Clade specific model - diversity dependent
#' * `"rr_k_0laa"` Clade specific model - diversity dependent
#' * `"rr_gam_di_0laa"` Clade specific model - diversity independent, relaxed
#' colonisation, and no anagenesis (laa fixed to zero). All other parameters
#' free
#' * `"rr_gam_dd_0laa"` Clade specific model - diversity dependent, relaxed
#' colonisation, and no anagenesis (laa fixed to zero). All other parameters
#' free
#' @param data_name String. Will be used for the name of the created output
#' folder.
#' @param results_root_folder Character. A path to the root folder containing
Expand Down Expand Up @@ -184,6 +200,12 @@
#' calculations, it sets the limit for the maximum number of species for
#' which a probability must be computed, which must be larger than the size
#' of the largest clade.
#' @param prop_type2_pool A numeric determining the proportion of the mainland
#' species pool that is composed on type 2 species.
#' @param par_upper_bound A numeric defining the upper limit of the integration
#' of a parameter when fitting the relaxed-rate DAISIE model. If the DAISIE
#' model being applied is not the relaxed-rate model, this parameter can be
#' ignored and left as its default as it does not influence the model.
#'
#' @return Nothing
#' @keywords internal
Expand Down Expand Up @@ -225,6 +247,8 @@ default_params_doc <- function(model,
mainland_n,
low_rates,
rep_index,
res) {
res,
prop_type2_pool,
par_upper_bound) {
# Nothing
}
14 changes: 10 additions & 4 deletions R/get_available_models.R
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,19 @@ get_available_models <- function() {
"cr_dd", "cr_di", "cr_dd_0laa", "cr_di_0laa",
"cr_di_0lac", "cr_dd_0lac", "rr_lac_di",
"rr_lac_dd", "rr_mu_di", "rr_mu_dd", "rr_k",
"rr_gam_di", "rr_gam_dd",
"rr_laa_di", "rr_laa_dd", "rr_mu_di_0lac",
"rr_mu_dd_0lac", "rr_k_0lac", "rr_laa_di_0lac",
"rr_laa_dd_0lac", "rr_lac_di_0laa",
"rr_mu_dd_0lac", "rr_k_0lac",
"rr_gam_di_0lac", "rr_gam_dd_0lac",
"rr_laa_di_0lac", "rr_laa_dd_0lac", "rr_lac_di_0laa",
"rr_lac_dd_0laa", "rr_mu_di_0laa", "rr_mu_dd_0laa",
"rr_k_0laa", "nonoceanic_cr_dd", "nonoceanic_cr_di",
"rr_k_0laa", "rr_gam_di_0laa", "rr_gam_dd_0laa",
"nonoceanic_cr_dd", "nonoceanic_cr_di",
"nonoceanic_cr_dd_0laa", "nonoceanic_cr_di_0laa",
"nonoceanic_cr_di_0lac", "nonoceanic_cr_dd_0lac"
"nonoceanic_cr_di_0lac", "nonoceanic_cr_dd_0lac",
"cr_dd_2type_lac", "cr_dd_2type_mu", "cr_dd_2type_k",
"cr_dd_2type_lac_mu", "cr_dd_2type_lac_k", "cr_dd_2type_mu_k",
"cr_dd_2type_lac_mu_k"
)
available_models
}
134 changes: 134 additions & 0 deletions R/run_daisie_2type_ml.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,134 @@
#' Run 2 type DAISIE analysis
#'
#' @inheritParams default_params_doc
#'
#' @return Nothing. Writes [`DAISIE`] analysis results to a `.rds` file. This
#' file is stored in `file_path`. The directory in of `file_path` is created
#' if it doesn't exist.
#' @export
#'
#' @examples
#' \dontrun{
#' data(Galapagos_datalist, package = "DAISIE")
#' run_daisie_2type_ml(
#' daisie_data = Galapagos_datalist,
#' data_name = "Galapagos_datalist",
#' model = "cr_dd",
#' array_index = 1,
#' cond = 1
#' )
#' }
#' @author Pedro Santos Neves, Joshua W. Lambert, Luis Valente
run_daisie_2type_ml <- function(daisie_data,
data_name,
model,
array_index,
cond,
methode = "lsodes",
optimmethod = "subplex",
results_dir = NULL,
low_rates = FALSE,
rep_index = "NULL",
res = 100,
prop_type2_pool,
test = FALSE) {
if (test) {
seed <- array_index
} else {
seed <- as.integer(Sys.time()) %% 10000L * array_index
}

set.seed(
seed,
kind = "Mersenne-Twister",
normal.kind = "Inversion",
sample.kind = "Rejection"
)

print_metadata(
data_name = data_name,
model = model,
array_index = array_index,
seed = seed,
methode = methode,
optimmethod = optimmethod
)

output_folder_path <- create_output_folder(
data_name = data_name,
results_dir = results_dir
)

testit::assert(is.numeric(array_index) && is.finite(array_index))
testit::assert(is.numeric(cond) && is.finite(cond))

model_arguments <- setup_2type_model(
model = model,
prop_type2_pool = prop_type2_pool,
low_rates = low_rates
)

initparsopt <- model_arguments$initparsopt
idparsnoshift <- model_arguments$idparsnoshift
idparsopt <- model_arguments$idparsopt
parsfix <- model_arguments$parsfix
idparsfix <- model_arguments$idparsfix
ddmodel <- model_arguments$ddmodel
cs_version <- model_arguments$cs_version

##### ML Optimization ####
lik_res <- DAISIE::DAISIE_ML(
datalist = daisie_data,
initparsopt = initparsopt,
idparsnoshift = idparsnoshift,
idparsopt = idparsopt,
parsfix = parsfix,
idparsfix = idparsfix,
res = res,
ddmodel = ddmodel,
cond = cond,
methode = methode,
optimmethod = optimmethod,
CS_version = cs_version
)

bic <- calc_bic(results = lik_res, daisie_data = daisie_data)
lik_res <- cbind(lik_res, bic)

if (identical(rep_index, "NULL")) {
output_path <- file.path(
output_folder_path,
paste0(
data_name,
"_",
model,
"_",
array_index,
".rds"
)
)
} else {
output_path <- file.path(
output_folder_path,
paste0(
data_name,
"_",
model,
"_",
array_index,
"_",
rep_index,
".rds"
)
)
}

if (is.na(output_path)) {
return(lik_res)
} else {
saveRDS(
lik_res,
file = output_path
)
}
}
4 changes: 3 additions & 1 deletion R/run_daisie_ml.R
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ run_daisie_ml <- function(daisie_data,
low_rates = FALSE,
rep_index = "NULL",
res = 100,
par_upper_bound = Inf,
test = FALSE) {
if (test) {
seed <- array_index
Expand Down Expand Up @@ -63,7 +64,8 @@ run_daisie_ml <- function(daisie_data,

model_arguments <- setup_model(
model = model,
low_rates = low_rates
low_rates = low_rates,
par_upper_bound = par_upper_bound
)

initparsopt <- model_arguments$initparsopt
Expand Down
Loading

0 comments on commit 7b822b7

Please sign in to comment.