Fix buggy brts_table for full_table = TRUE and change tests according…

…ly. Also add test on IW loglik.

R/DAISIE_format_IW.R

@@ -117,7 +117,7 @@ DAISIE_format_IW_trait <- function(island_replicates,
   return(several_islands)
 }
 
-add_brt_table <- function(island, full_table = FALSE) {
+add_brt_table_old <- function(island, full_table = FALSE) {
   island_age <- island[[1]]$island_age
   island_top <- island[[1]]
   if (length(island) == 1) {
@@ -204,8 +204,7 @@ add_brt_table <- function(island, full_table = FALSE) {
   return(island)
 }
 
-#The function below makes more sense, but DAISIE_loglik_IW needs to be adjusted to it
-add_brt_table2 <- function(island, full_table = FALSE) {
+add_brt_table <- function(island, full_table = TRUE) {
   island_age <- island[[1]]$island_age
   island_top <- island[[1]]
   if (length(island) == 1) {
@@ -247,7 +246,7 @@ add_brt_table2 <- function(island, full_table = FALSE) {
       for (i in 1:length(btimes)) {
         the_stac <- island_no_stac1or5[[i]]$stac
         if(!is.null(island[[i]]$all_colonisations) & full_table == TRUE) {
-          if(length(island[[i]]$all_colonisations) > 0) print(i)
+          #if(length(island[[i]]$all_colonisations) > 0) print(i)
           for (k in 1:length(island[[i]]$all_colonisations)) {
             the_brts <- island[[i]]$all_colonisations[[k]]$event_times[-1]
             pos2 <- pos1 + length(the_brts)

R/DAISIE_loglik_IW.R

@@ -427,7 +427,12 @@ DAISIE_loglik_IW <- function(
           test_for_colonization <- TRUE
         } else
         {
-          test_for_colonization <- (max(event[which(clade == col[k + 2])]) > 1)
+          #test_for_colonization <- (max(event[which(clade == col[k + 2])]) > 1)
+          # this tests whether the original clade has diversified, but this does not need to be the case
+          # it could also be a case of anagenesis followed by colonization, so this test seems inappropriate
+          # also, one should expect that colonizations in the data are all real colonizations; a colonization
+          # which is not followed by speciation but by recolonization should only occur once (the last colonization)
+          test_for_colonization <- TRUE
         }
         if(test_for_colonization) # new colonization or recolonization after speciation
         {

tests/testthat/test-DAISIE_format_IW.R

@@ -94,11 +94,11 @@ test_that("silent with non-empty island with correct output", {
   brts_table <- matrix(ncol = 5, nrow = 6)
   colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col")
   brts_table[1, ] <- c(1, 0, 0, NA, NA)
-  brts_table[2, ] <- c(0.9244818166871660, 1, 1, 1, NA)
-  brts_table[3, ] <- c(0.9105856673960619, 1, 2, 1, NA)
-  brts_table[4, ] <- c(0.5557734125062590, 2, 1, 0, NA)
-  brts_table[5, ] <- c(0.5288428248966160, 3, 1, 0, NA)
-  brts_table[6, ] <- c(0.3146835586399670, 1, 3, 1, NA)
+  brts_table[2, ] <- c(0.9244818166871660, 1, 1, 1, 1)
+  brts_table[3, ] <- c(0.9105856673960619, 1, 2, 1, 1)
+  brts_table[4, ] <- c(0.5557734125062590, 2, 1, 0, 1)
+  brts_table[5, ] <- c(0.5288428248966160, 3, 1, 0, 1)
+  brts_table[6, ] <- c(0.3146835586399670, 1, 3, 1, 1)
   expected_IW_format[[1]][[1]] <- list(island_age = 1,
                                        not_present = 7,
                                        stt_all = stt_all,
@@ -255,8 +255,8 @@ test_that("silent with non-empty nonoceanic island with
             brts_table <- matrix(ncol = 5, nrow = 3)
             colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col")
             brts_table[1, ] <- c(1, 0, 0, NA, NA)
-            brts_table[2, ] <- c(1, 2, 1, 1, NA)
-            brts_table[3, ] <- c(1, 1, 1, 1, NA)
+            brts_table[2, ] <- c(1, 2, 1, 1, 1)
+            brts_table[3, ] <- c(1, 1, 1, 1, 1)
             expected_IW_format[[1]][[1]] <- list(island_age = 1,
                                                  not_present = 8,
                                                  stt_all = stt_all,
@@ -316,8 +316,8 @@ test_that("silent with non-empty nonoceanic island with
             brts_table <- matrix(ncol = 5, nrow = 3)
             colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col")
             brts_table[1, ] <- c(1, 0, 0, NA, NA)
-            brts_table[2, ] <- c(1, 2, 1, 1, NA)
-            brts_table[3, ] <- c(1, 1, 1, 1, NA)
+            brts_table[2, ] <- c(1, 2, 1, 1, 1)
+            brts_table[3, ] <- c(1, 1, 1, 1, 1)
             expected_IW_format[[1]][[1]] <- list(island_age = 1,
                                                  not_present = 8,
                                                  stt_all = stt_all,
@@ -385,10 +385,10 @@ test_that("add_brt_table output is correct when length(island) != 1", {
   brt_table <- matrix(ncol = 5, nrow = 5)
   colnames(brt_table) <- c("brt", "clade", "event", "endemic", "col")
   brt_table[1, ] <- c(1, 0, 0, NA, NA)
-  brt_table[2, ] <- c(0.9244818, 1, 1, 1, NA)
-  brt_table[3, ] <- c(0.9105857, 1, 2, 1, NA)
-  brt_table[4, ] <- c(0.5557734, 2, 1, 0, NA)
-  brt_table[5, ] <- c(0.3146836, 1, 3, 1, NA)
+  brt_table[2, ] <- c(0.9244818, 1, 1, 1, 1)
+  brt_table[3, ] <- c(0.9105857, 1, 2, 1, 1)
+  brt_table[4, ] <- c(0.5557734, 2, 1, 0, 1)
+  brt_table[5, ] <- c(0.3146836, 1, 3, 1, 1)
   expected_brt <- list()
   expected_brt[[1]] <- list(island_age = 1,
                             not_present = 3,
@@ -540,9 +540,9 @@ test_that("silent when species with two trait states with
             brts_table <- matrix(ncol = 5, nrow = 4)
             colnames(brts_table) <- c("brt", "clade", "event", "endemic", "col")
             brts_table[1, ] <- c(5.00000000000000, 0, 0, NA, NA)
-            brts_table[2, ] <- c(3.10261367452990, 1, 1, 1, NA)
-            brts_table[3, ] <- c(1.50562999775257, 2, 1, 1, NA)
-            brts_table[4, ] <- c(1.26245655913561, 2, 2, 1, NA)
+            brts_table[2, ] <- c(3.10261367452990, 1, 1, 1, 1)
+            brts_table[3, ] <- c(1.50562999775257, 2, 1, 1, 1)
+            brts_table[4, ] <- c(1.26245655913561, 2, 2, 1, 1)
             expected_IW_format[[1]][[1]] <- list(island_age = 5,
                                                  not_present = 13,
                                                  stt_all = stt_all,

tests/testthat/test-DAISIE_format_IW_full_stt.R

@@ -67,7 +67,7 @@ test_that("complete stt, 1 type, no geodynamics, oceanic island, one trait state
     c(Time = 0.0, nI = 0.0, nA = 1.0, nC = 0.0)
   )
   expected_brts_table <- matrix(
-    c(1.0, 0.24481816687165, 0, 1, 0, 1, NA, 1, NA, NA),
+    c(1.0, 0.24481816687165, 0, 1, 0, 1, NA, 1, NA, 1),
     nrow = 2
   )
   colnames(expected_brts_table) <- c("brt", "clade", "event", "endemic", "col")
@@ -109,7 +109,7 @@ test_that("complete stt, 1 type, no geodynamics, oceanic island, one trait state
   )
 
   expected_brts_table <- matrix(
-    c(1.0, 0.741771912202239, 0, 1, 0, 1, NA, 0, NA, NA),
+    c(1.0, 0.741771912202239, 0, 1, 0, 1, NA, 0, NA, 1),
     nrow = 2
   )
   colnames(expected_brts_table) <- c("brt", "clade", "event", "endemic", "col")
@@ -245,7 +245,7 @@ test_that("complete stt, 1 type, geodynamics, oceanic island, one trait state
   )
   testthat::expect_equal(
     formatted_IW_sim[[1]][[1]]$brts_table[5, ],
-    c(brt = 0.83094531417507, clade = 4, event = 1, endemic = 1, col = NA)
+    c(brt = 0.83094531417507, clade = 4, event = 1, endemic = 1, col = 1)
   )
 
 })

tests/testthat/test-DAISIE_format_IW_sampled_stt.R

@@ -59,9 +59,9 @@ test_that("sampled stt, 1 type, no geodynamics, oceanic island (same arguments a
   expected_brts_table <- matrix(
     data = c(
       1.0, 0.0, 0.0, NA, NA,
-      0.244818166871655, 1, 1, 1, NA,
-      0.173128288990374, 1, 2, 1, NA,
-      0.029668240213840, 1, 3, 1, NA
+      0.244818166871655, 1, 1, 1, 1,
+      0.173128288990374, 1, 2, 1, 1,
+      0.029668240213840, 1, 3, 1, 1
     ),
     ncol = 5,
     byrow = TRUE

tests/testthat/test-integration_DAISIE.R

@@ -100,6 +100,44 @@ test_that("IW and CS loglik is same when K = Inf", {
   testthat::expect_equal(loglik_IW, loglik_CS, tol = 5E-6)
 })
 
+test_that("IW and CS loglik is same when K = Inf", {
+  skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"),
+          message = "Run only on CI")
+  skip_on_cran()
+  M <- 1000
+  data(frogs_datalist, package = "DAISIE")
+  frogs_datalist[[1]]$not_present <- frogs_datalist[[1]]$not_present + (M - 300)
+  ddmodel <- 11
+  initparsopt <- c(4.012298e-01,1.699521e-01,1.319595e+02,3.487955e-04)
+  idparsopt <- c(1,2,3,4)
+  parsfix <- 0
+  idparsfix <- 5
+  verbose <- 1
+  cond <- 1
+  res <- 200
+  methode <- 'odeint::runge_kutta_fehlberg78'
+  tolint <- c(1E-16, 1E-14)
+
+  trparsopt <- initparsopt / (1 + initparsopt)
+  trparsopt[which(initparsopt == Inf)] <- 1
+  trparsfix <- parsfix / (1 + parsfix)
+  trparsfix[which(parsfix == Inf)] <- 1
+  pars2 <- c(res, ddmodel, cond, verbose)
+  initloglik_IW <- DAISIE_loglik_IW_choosepar(trparsopt = trparsopt,
+                                              trparsfix = trparsfix,
+                                              idparsopt = idparsopt,
+                                              idparsfix = idparsfix,
+                                              M = M,
+                                              pars2 = pars2,
+                                              datalist = frogs_datalist,
+                                              methode = methode,
+                                              abstolint = tolint[1],
+                                              reltolint = tolint[2])
+
+  testthat::expect_equal(initloglik_IW, -215.097677998973, tol = 1E-6)
+})
+
+
 test_that("DAISIE_ML simple case works", {
   skip_if(Sys.getenv("CI") == "" && !(Sys.getenv("USERNAME") == "rampa"),
           message = "Run only on CI")