Develop

DESCRIPTION

@@ -2,7 +2,7 @@ Package: secsse
 Type: Package
 Title: Several Examined and Concealed States-Dependent Speciation and
  Extinction
-Version: 3.6.0
+Version: 3.6.1
 License: GPL (>= 3) | file LICENSE
 Authors@R: c(
  person(given = "Leonel",
@@ -50,7 +50,8 @@ Imports:
  ggplot2,
  tibble,
  rlang,
- treestats
+ treestats,
+ pracma
 Suggests: 
  diversitree,
  phytools,

NEWS.md

@@ -1,10 +1,18 @@
+# 3.6.1
+- added support for having a different number of concealed states compared to
+the number of observed states when creating a q_matrix using `q_doubletrans`.
+
 # 3.6.0
 - 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

R/default_params_doc.R

@@ -33,7 +33,8 @@
 #' @param cond condition on the existence of a node root: `"maddison_cond"`,
 #'  `"proper_cond"` (default). For details, see vignette.
 #' @param root_state_weight the method to weigh the states:
-#'  `"maddison_weights"`, `"proper_weights"` (default) or `"equal_weights"`.
+#'  `"maddison_weights"`, `"proper_weights"` (default), `"equal_weights"'.
+#'  or `"stationary_weights"`
 #'  It can also be specified for the root state: the vector `c(1, 0, 0)` 
 #'  indicates state 1 was the root state. When
 #'  using a `multiPhylo` object, root_state_weight should be list where each

R/seccse_multi_loglik.R

@@ -45,31 +45,31 @@ multi_loglik <- function(parameter,
     }
 
     if (length(phy[[i]]$tip.label) == 1) {
-      local_res <- secsse::secsse_single_branch_loglik(parameter = parameter,
-                                                      phy = phy[[i]],
-                                                      traits = traits[[i]],
-                                                      num_concealed_states =
-                                                        num_concealed_states,
-                                                      cond = cond,
-                                                      root_state_weight = 
-                                                        focal_root_state_weight,
-                                                      sampling_fraction = 
-                                                        focal_sampling_fraction,
-                                                      setting_calculation = 
-                                                        focal_setting_calculation,
-                                                      see_ancestral_states = FALSE,
-                                                      loglik_penalty = loglik_penalty,
-                                                      is_complete_tree = 
-                                                        is_complete_tree,
-                                                      take_into_account_root_edge = 
-                                                        take_into_account_root_edge,
-                                                      num_threads = num_threads,
-                                                      atol = atol,
-                                                      rtol = rtol,
-                                                      method = method,
-                                                      display_warning = display_warning,
-                                                      use_normalization = use_normalization,
-                                                      return_root_state = return_root_state)
+      local_res <- secsse_single_branch_loglik(parameter = parameter,
+                                               phy = phy[[i]],
+                                               traits = traits[[i]],
+                                               num_concealed_states =
+                                                 num_concealed_states,
+                                               cond = cond,
+                                               root_state_weight = 
+                                                 focal_root_state_weight,
+                                               sampling_fraction = 
+                                                 focal_sampling_fraction,
+                                               setting_calculation = 
+                                                 focal_setting_calculation,
+                                               see_ancestral_states = FALSE,
+                                               loglik_penalty = loglik_penalty,
+                                               is_complete_tree = 
+                                                 is_complete_tree,
+                                               take_into_account_root_edge = 
+                                                 take_into_account_root_edge,
+                                               num_threads = num_threads,
+                                               atol = atol,
+                                               rtol = rtol,
+                                               method = method,
+                                               display_warning = display_warning,
+                                               use_normalization = use_normalization,
+                                               return_root_state = return_root_state)
       if (return_root_state) {
         root_states[[i]] <- local_res$root_state
       } 
@@ -107,6 +107,8 @@ multi_loglik <- function(parameter,
 
   ll <- do.call(sum, res)
 
+
+
   if (return_root_state) {
     return(list(LL = ll,
                 root_state = root_states))

R/secsse_loglik.R

@@ -54,32 +54,32 @@ master_loglik <- function(parameter,
   }
 
   if (length(phy$tip.label) == 1) {
-    return(secsse::secsse_single_branch_loglik(
-                       parameter = parameter,
-                       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,
-                       take_into_account_root_edge = 
-                         take_into_account_root_edge,
-                       num_threads = num_threads,
-                       atol = atol,
-                       rtol = rtol,
-                       method = method,
-                       display_warning = display_warning,
-                       use_normalization = use_normalization,
-                       return_root_state = return_root_state))
+    return(secsse_single_branch_loglik(
+      parameter = parameter,
+      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,
+      take_into_account_root_edge = 
+        take_into_account_root_edge,
+      num_threads = num_threads,
+      atol = atol,
+      rtol = rtol,
+      method = method,
+      display_warning = display_warning,
+      use_normalization = use_normalization,
+      return_root_state = return_root_state))
   }
 
   lambdas <- parameter[[1]]
@@ -113,7 +113,6 @@ master_loglik <- function(parameter,
   states <- setting_calculation$states
   forTime <- setting_calculation$forTime
   ances <- setting_calculation$ances
-
   d <- ncol(states) / 3
 
   # with a complete tree, we need to re-calculate the states every time we
@@ -193,7 +192,8 @@ master_loglik <- function(parameter,
                                      lambdas,
                                      nodeM,
                                      d,
-                                     is_cla = using_cla)
+                                     is_cla = using_cla,
+                                     Q = q_matrix)
 
   if (is_complete_tree) {
     nodeM <- update_complete_tree(phy,
@@ -325,25 +325,26 @@ secsse_loglik <- function(parameter,
                           display_warning = TRUE,
                           use_normalization = TRUE,
                           return_root_state = FALSE) {
-  master_loglik(parameter = parameter,
-                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,
-                take_into_account_root_edge = take_into_account_root_edge,
-                num_threads = num_threads,
-                atol = atol,
-                rtol = rtol,
-                method = method,
-                display_warning = display_warning,
-                use_normalization = use_normalization,
-                return_root_state = return_root_state)
+  ll <- master_loglik(parameter = parameter,
+                      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,
+                      take_into_account_root_edge = take_into_account_root_edge,
+                      num_threads = num_threads,
+                      atol = atol,
+                      rtol = rtol,
+                      method = method,
+                      display_warning = display_warning,
+                      use_normalization = use_normalization,
+                      return_root_state = return_root_state)
+  return(ll)
 }
 
 

R/secsse_ml.R

@@ -29,6 +29,7 @@ master_ml <- function(phy,
                       use_normalization = TRUE,
                       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,
@@ -81,16 +82,6 @@ master_ml <- function(phy,
   num_modeled_traits <- length(idparslist[[1]]) / num_concealed_states
 
   if (!is.list(traits)) {
-
-    if (length(phy$tip.label) == 1) {
-      fake_phy <- ape::rphylo(n = 2, birth = 1, death = 0)
-      fake_phy$edge.length[1:2] <- phy$edge.length[1]
-      old_phy <- phy
-      old_traits <- traits
-      phy <- fake_phy
-      traits <- c(traits, traits)
-    }
-
     setting_calculation <- build_initStates_time(phy,
                                                  traits,
                                                  num_concealed_states,
@@ -100,17 +91,11 @@ master_ml <- function(phy,
                                                  num_modeled_traits,
                                                  traitStates = 
                                                    get_trait_states(idparslist,
-                                                                    num_concealed_states, FALSE))
-    if (length(phy$tip.label) == 1) {
-      setting_calculation$states <- setting_calculation$states[-2, ]
-      setting_calculation$forTime <- setting_calculation$forTime[-2, ]
-      phy <- old_phy
-      traits <- old_traits
-    }
+                                                                    num_concealed_states, verbose))
+
   } else {
     setting_calculation <- list()
     for (i in 1:length(phy)) {
-
       input_phy <- phy[[i]]
       input_traits <- traits[[i]]
 
@@ -119,14 +104,7 @@ master_ml <- function(phy,
       } else {
         input_sampling_fraction <- sampling_fraction
       }
-
-      if (length(input_phy$tip.label) == 1) {
-        fake_phy <- ape::rphylo(n = 2, birth = 1, death = 0)
-        fake_phy$edge.length[1:2] <- input_phy$edge.length[1]
-        input_phy <- fake_phy
-        input_traits <- c(input_traits, input_traits)
-      }
-
+
       setting_calculation[[i]] <- build_initStates_time(phy = input_phy,
                                                         traits = input_traits,
                                                         num_concealed_states =
@@ -141,7 +119,7 @@ master_ml <- function(phy,
                                                         first_time = FALSE,
                                                         traitStates = 
                                                           get_trait_states(idparslist,
-                                                                           num_concealed_states, FALSE))
+                                                                           num_concealed_states, verbose))
     }
   }