R/result_analysis.R

In Neves-P/utilSIE: Utility Functions for DAISIE

#' Calculate means and medians of parameter estimates
#'
#' @return list with median number of independent clades
#' @export
#'
get_summary_stats <- function() {
  project_name <- utilSIE::get_project_name()
  project_folder <- utilSIE::get_project_folder(project_name)
  platform <- .Platform$OS.type

  if (platform == "windows") {
    data_folder <- local_data_folder <- file.path(project_folder, "data")
    testit::assert(dir.exists(data_folder))
  } else {
    data_folder <- file.path(utilSIE::get_project_name(), "data")
  }
  out <- NULL # nolint; keep R CMD CHECK happy
  files <- list.files(data_folder)
  for (file in seq_along(files)) {
    load(
      file = file.path(local_data_folder, files[file])
    )
    assign(paste0("sim_", file), out) # nolint
    assign(paste0("args_", file), args) # nolint
  }

  # Initialize variables
  datasets <- ls(pattern = "sim")
  n_indep_clade <- c()
  medians_indep_clades <- c()
  medians_n_spec_island <- c()
  n_spec_island <- c()

  for (dataset_id in seq_along(datasets)) {
    # Calc median of indep clades
    for (repl in seq_along(get(datasets[dataset_id]))) {
      n_indep_clade[repl] <- # Count number of indep clades using last line
        get(datasets[dataset_id])[[repl]][[1]]$stt_all[
          nrow(get(datasets[dataset_id])[[repl]][[1]]$stt_all), 5]
    }

    # Calc median of number of species
    for (repl in seq_along(get(datasets[dataset_id]))) {
      n_spec_island[repl] <- sum(get(datasets[dataset_id])[[repl]][[1]]$stt_all[
        nrow(get(datasets[dataset_id])[[repl]][[1]]$stt_all), 2:4])
    }
    medians_indep_clades[dataset_id] <- stats::median(n_indep_clade)
    medians_n_spec_island[dataset_id] <- stats::median(n_spec_island)
  }
  return(list(medians_indep_clades, medians_n_spec_island))
}


#' Trims and adds ids to MLE results given a lac threshold
#'
#' @inheritParams run_simulations_from_MLE
#'
#' @return Same data frame as \code{results_ML} but with extra column with the
#' ID of each simulation ran (so that these can be matched to the generating
#' simulations, e.g. ontogeny simulations).
#' @export
get_simulations_to_run <- function(results_ML, lac_treshold = 10) {
  assertthat::is.number(lac_treshold)
  assertthat::assert_that(lac_treshold < 20 && lac_treshold >= 0)

  normal_lac_res <- results_ML[which(results_ML$lambda_c < lac_treshold), ]
  normal_lac_res$id <- which(results_ML$lambda_c < lac_treshold)
  normal_lac_res
}


#' Start standard simulations from set of MLE results
#'
#' @description Useful function for runing robustness analyses. Runs a set of simulations
#' per result of MLE obtained previously.
#'
#' @param results_ML Data frame with MLE as returned
#' by \code{\link[DAISIE]{DAISIE_ML}}.
#' @param time Numeric containing lenght of simulation in time units.
#' @param replicates Number of replicates to be run per each set of MLE
#' (i.e., per line of the \code{results_ML} object).
#' @param M Numeric specifying number of mainland species.
#' @param lac_treshold Numeric with cutoff value for cladogenesis estimates.
#' Lines of \code{results_ML} with values above this threshold are not included
#' in the simulations.
#'
#' @return Same data frame as \code{results_ML} but with extra column with the
#' ID of each simulation ran (so that these can be matched to the generating
#' simulations, e.g. ontogeny simulations).
#' @export
run_simulations_from_MLE <- function(results_ML,
                                     time,
                                     replicates,
                                     M,
                                     lac_treshold = 10) {


  normal_lac_res <- get_simulations_to_run(results_ML, lac_treshold)



  for (i in 1:nrow(normal_lac_res)) {
    execute_next_setup(
      time = time,
      M = M,
      lac = normal_lac_res$lambda_c[i],
      mu = normal_lac_res$mu[i],
      K = normal_lac_res$K[i],
      gam = normal_lac_res$gamma[i],
      laa = normal_lac_res$lambda_a[i],
      island_ontogeny = FALSE,
      branch = "@develop",
      replicates = replicates,
      complete_analysis = FALSE
    )
  }
  return(normal_lac_res)
}