Develop

DESCRIPTION

@@ -3,7 +3,6 @@ Type: Package
 Title: Several Examined and Concealed States-Dependent Speciation and
  Extinction
 Version: 3.6.0
-Date: 2025-10-13
 License: GPL (>= 3) | file LICENSE
 Authors@R: c(
  person(given = "Leonel",
@@ -43,7 +42,7 @@ Description: Simultaneously infers state-dependent diversification across
 Depends: R (>= 4.2.0)
 Imports: 
  utils,
- DDD (>= 5.0),
+ DDD (>= 5.2.4),
  ape,
  geiger,
  Rcpp (>= 1.0.10),

NEWS.md

@@ -2,14 +2,10 @@
 - changed the default integration method to "odeint::runge_kutta_cash_karp54",
 although a bit slower, this method is more numerically stable.
 - fixed lack of sorting of numeric traits in the function `q_doubletrans`
-- fixed an error in preparing the state matrices when NAs were present in the 
-traits
-- ML and LL functions now optionally also return the root state, which in turn
-can be used to in `secsse_sim` as a starting point at the root / crown of the
-tree.
-- updated simulations to sample species using a binary search, instead of using
-stochastic acceptance
-
+- fixed an error in preparing the state matrices when NAs were present in the traits
+- ML and LL functions now optionally also return the root state, which in turn can be used to in `secsse_sim` as a starting point at the root / crown of the tree.
+- updated simulations to sample species using a binary search, instead of using stochastic acceptance
+- better documentation of output of ML and LL functions.
 # 3.5.0
 Version 3.5.0 uses a separate calculation for 1 - E, e.g. one minus the local
 extinction rates to avoid numerical aberrations. These are only used in the CLA

R/default_params_doc.R

@@ -65,13 +65,13 @@
 #' @param parameter list where first vector represents lambdas, the second 
 #'  mus and the third transition rates.
 #' @param setting_calculation argument used internally to speed up calculation.
-#'  It should be left blank (default : `setting_calculation = NULL`).
+#'  This should be left blank (default : `setting_calculation = NULL`).
 #' @param loglik_penalty the size of the penalty for all parameters; default is
 #'  0 (no penalty).
 #' @param num_steps number of substeps to show intermediate likelihoods
 #'  along a branch.
-#' @param see_ancestral_states Boolean for whether the ancestral states should 
-#'  be shown? Defaults to `FALSE`.
+#' @param see_ancestral_states Boolean for whether the ancestral states for each
+#'  of the internal nodes should be output. Defaults to `FALSE`.
 #' @param lambdas speciation rates, in the form of a list of matrices.
 #' @param mus extinction rates, in the form of a vector.
 #' @param qs The Q matrix, for example the result of function q_doubletrans, but

R/secsse_loglik.R

@@ -226,7 +226,14 @@ master_loglik <- function(parameter,
 
   # reset number of threads:
   RcppParallel::setThreadOptions(numThreads = 1)
+
+  if (!see_ancestral_states && !return_root_state) {
+    return(LL)
+  }
 
+
+  result <- list()
+  result$LL <- LL
   if (see_ancestral_states == TRUE) {
     states <- calcul$states
     num_tips <- ape::Ntip(phy)
@@ -235,26 +242,44 @@ master_loglik <- function(parameter,
       ancestral_states[, (1/3 * ncol(ancestral_states) + 1):(2/3 * ncol(ancestral_states))]
 
     rownames(ancestral_states) <- ances
-    return(list(ancestral_states = ancestral_states, LL = LL, states = states))
+    colnames(ancestral_states) <- names(mus)
+    colnames(states) <- c(paste("E",names(mus), sep = '_'),
+                          paste("D",names(mus), sep = '_'),
+                          paste("S",names(mus), sep = '_'))
+    result$ancestral_states <- ancestral_states
+    result$states = states
   } 
 
   if (return_root_state) {
-    return(list(LL = LL,
-                root_state = get_root_state(calcul$states,
-                                            phy,
-                                            mus,
-                                            d)))
+    root_state <- get_root_state(calcul$states,
+                                 phy,
+                                 mus,
+                                 d)
+    result$root_state <- root_state
   }
-
-  return(LL)
+  return(result)
 }
 
 #' @title Likelihood for SecSSE model
 #' Loglikelihood calculation for the SecSSE model given a set of parameters and
 #' data
 #' 
 #' @inheritParams default_params_doc
-#' @return The loglikelihood of the data given the parameter.
+#' @return A list with the following elements:
+#' $LL the loglikelihood of the data (phylogeny + tip states) given the
+#' parameters (speciation, extinction, transition rates).
+#' If see_ancestral_states = TRUE, then there will be two additional elements:
+#' $ancestral_states: a matrix with the probabilities of each state at the
+#' internal nodes
+#' $states: a matrix with the probabilities E, D (normalized) and S that are used
+#' in the calculations. The ancestral_states matrix is a submatrix of this matrix.
+#' This matrix is mostly used for package developers.
+#' If return_root_state = TRUE, then there will be one additional element:
+#' $root_state: vector with probabilities of each state at the root. This vector
+#' is the same as the top row of $ancestral_states
+#' We have used the shorthand description of "probabilities of each state",
+#' but technically, the probabilities are the normalized probabilities D of the
+#' data given each state at the internal nodes.
 #' @examples
 #' rm(list = ls(all = TRUE))
 #' library(secsse)
@@ -328,7 +353,8 @@ secsse_loglik <- function(parameter,
 #' 
 #' @inheritParams default_params_doc
 #' 
-#' @return The loglikelihood of the data given the parameters
+#' @return A List with property LL: The loglikelihood of the data given the 
+#' parameters, and potentially the root state.
 #' @examples
 #'rm(list=ls(all=TRUE))
 #'library(secsse)

R/secsse_ml.R

@@ -27,8 +27,8 @@ master_ml <- function(phy,
                       rtol = 1e-7,
                       method = "odeint::runge_kutta_cash_karp54",
                       use_normalization = TRUE,
-                      return_root_state = FALSE) {
-
+                      return_root_state = FALSE,
+                      see_ancestral_states = FALSE) {
   structure_func <- NULL
   if (!is.null(functions_defining_params)) {
     structure_func <- set_and_check_structure_func(idparsfuncdefpar,
@@ -60,7 +60,6 @@ master_ml <- function(phy,
                                             idparslist[[1]])
   }
 
-  see_ancestral_states <- FALSE
   if (!is.null(structure_func)) {
     initparsopt <- c(initparsopt, initfactors)
   }
@@ -226,12 +225,10 @@ master_ml <- function(phy,
                                               idparsfix,
                                               idparslist,
                                               structure_func)
-
-      if (!return_root_state) {
-        out2 <- list(MLpars = ml_pars1,
-                     ML = as.numeric(unlist(out$fvalues)),
-                     conv = out$conv)
-      } else {
+      out2 <- list(MLpars = ml_pars1,
+                   ML = as.numeric(unlist(out$fvalues)),
+                   conv = out$conv)
+      if (return_root_state || see_ancestral_states) {
         # we have to recover the root state by re-calculating the LL
         root_res <- master_loglik(parameter = ml_pars1,
                                   phy = phy,
@@ -262,11 +259,13 @@ master_ml <- function(phy,
         if (as.numeric(unlist(out$fvalues)) != root_res$LL) {
           warning("recalculating the root state resulted in a different LL")
         }
-
-        out2 <- list(MLpars = ml_pars1,
-                     ML = as.numeric(unlist(out$fvalues)),
-                     conv = out$conv,
-                     root_state = root_res$root_state)
+      }
+      if(return_root_state) {
+        out2$root_state <- root_res$root_state
+      }
+      if(see_ancestral_states) {
+        out2$ancestral_states <- root_res$ancestral_states
+        out2$states <- root_res$states
       }
     }
   }
@@ -276,10 +275,26 @@ master_ml <- function(phy,
 #' Maximum likehood estimation for (SecSSE)
 #' 
 #' Maximum likehood estimation under Several examined and concealed
-#' States-dependent Speciation and Extinction (SecSSE)
+#' trait States dependent Speciation and Extinction (SecSSE)
 #' @inheritParams default_params_doc
 #' 
-#' @return Parameter estimated and maximum likelihood
+#' @return A list with the following elements
+#' $MLpars: the maximum likelihood parameter estimates
+#' $ML: the maximum likelihood of the data (phylogeny + tip states) given the
+#' parameters (speciation, extinction, transition rates).
+#' $conv: whether the optimization converged or not
+#' If see_ancestral_states = TRUE, then there will be two additional elements:
+#' $ancestral_states: a matrix with the probabilities of each state at the
+#' internal nodes
+#' $states: a matrix with the probabilities E, D (normalized) and S that are used
+#' in the calculations. The ancestral_states matrix is a submatrix of this matrix.
+#' This matrix is mostly used for package developers.
+#' If return_root_state = TRUE, then there will be one additional element:
+#' $root_state: vector with probabilities of each state at the root. This vector
+#' is the same as the top row of $ancestral_states
+#' We have used the shorthand description of "probabilities of each state",
+#' but technically, the probabilities are the normalized probabilities D of the
+#' data given each state at the internal nodes.
 #' @examples
 #'# Example of how to set the arguments for a ML search.
 #'library(secsse)
@@ -357,7 +372,8 @@ secsse_ml <- function(phy,
                       rtol = 1e-7,
                       method = "odeint::runge_kutta_cash_karp54",
                       use_normalization = TRUE,
-                      return_root_state = FALSE) {
+                      return_root_state = FALSE,
+                      see_ancestral_states = FALSE) {
   master_ml(phy = phy,
             traits = traits,
             num_concealed_states = num_concealed_states,
@@ -384,7 +400,8 @@ secsse_ml <- function(phy,
             rtol = rtol,
             method = method,
             use_normalization = use_normalization,
-            return_root_state = return_root_state)
+            return_root_state = return_root_state,
+            see_ancestral_states = see_ancestral_states)
 }
 
 #' @keywords internal
@@ -446,7 +463,7 @@ secsse_loglik_choosepar <- function(trparsopt,
                             rtol = rtol,
                             display_warning = display_warning,
                             use_normalization = use_normalization,
-                            return_root_state = return_root_state)
+                            return_root_state = FALSE)
 
     if (is.nan(loglik) || is.na(loglik)) {
       warning("There are parameter values used which cause
@@ -547,7 +564,8 @@ cla_secsse_ml <- function(phy,
                           rtol = 1e-7,
                           method = "odeint::runge_kutta_cash_karp54",
                           use_normalization = TRUE,
-                          return_root_state = FALSE) {
+                          return_root_state = FALSE,
+                          see_ancestral_states = FALSE) {
   master_ml(phy = phy,
             traits = traits,
             num_concealed_states = num_concealed_states,
@@ -572,5 +590,6 @@ cla_secsse_ml <- function(phy,
             rtol = rtol,
             method = method,
             use_normalization = use_normalization,
-            return_root_state = return_root_state)
+            return_root_state = return_root_state,
+            see_ancestral_states = see_ancestral_states)
 }