Merge pull request #143 from rsetienne/develop

v4.2.0

.Rbuildignore

@@ -15,3 +15,4 @@
 ^\.Rproj\.user$
 ^\.covrignore$
 ^LICENSE$
+^\.zenodo\.json$

.zenodo.json

@@ -0,0 +1,104 @@
+{
+   "creators":[
+      {
+         "orcid":"0000-0003-2142-7612",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Etienne, Rampal S."
+      },
+      {
+         "orcid":"0000-0003-4247-8785",
+         "affiliation":"Naturalis Biodiversity Center",
+         "name":"Valente, Luis"
+      },
+      {
+         "orcid":"0000-0002-6553-1553",
+         "affiliation":"Institute of Evolutionary Biology, University of Edinburgh",
+         "name":"Phillimore, Albert B."
+      },
+      {
+         "orcid":"0000-0003-2325-4727",
+         "affiliation":"Station d’Ecologie Théorique et Expérimentale, CNRS",
+         "name":"Haegeman, Bart"
+      },
+      {
+         "orcid":"0000-0001-5218-3046",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Lambert, Joshua W."
+      },
+      {
+         "orcid":"0000-0003-2561-4677",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Santos Neves, Pedro"
+      },
+      {
+         "orcid":"0000-0001-9594-946X",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Xie, Shu"
+      },
+      {
+         "orcid":"0000-0003-1107-7049",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Bilderbeek, Rich\u00e8l J.C."
+      },
+      {
+         "orcid":"0000-0002-6784-1037",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Hildenbrandt, Hanno"
+      }
+   ],
+   "contributors":[
+      {
+         "orcid":"0000-0001-5711-9457",
+         "affiliation":"Justus Liebig University",
+         "name":"Hauffe, Torsten",
+         "type":"Other"
+      },
+      {
+         "orcid":"0000-0002-2952-3345",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Laudanno, Giovanni",
+         "type":"Other"
+      },
+      {
+         "orcid":"0000-0002-9720-4581",
+         "affiliation":"National University of Singapore",
+         "name":"Kristensen, Nadiah",
+         "type":"Other"
+      },
+      {
+         "orcid":"0000-0002-1447-7630",
+         "affiliation":"Groningen Institute for Evolutionary Life Sciences, University of Groningen",
+         "name":"Scherrer, Raphael",
+         "type":"Other"
+      }
+   ],
+   "license":"GPL-3.0-only",
+   "title":"DAISIE: Dynamical Assembly of Islands by Speciation, Immigration and Extinction",
+   "description":"Simulates and computes the (maximum) likelihood of a dynamical model of island biota assembly through speciation, immigration and extinction.",
+   "access_right":"open",
+   "keywords":[
+      "island biogeography",
+      "birth-death process",
+      "maximum-likelihood estimation",
+      "R package"
+   ],
+   "related_identifiers":[
+      {
+         "scheme":"url",
+         "identifier":"https://github.com/rsetienne/DAISIE",
+         "relation":"isSupplementTo",
+         "resource_type":"software"
+      },
+      {
+         "scheme":"url",
+         "identifier":"https://cran.r-project.org/package=DAISIE",
+         "relation":"isSupplementTo",
+         "resource_type":"software"
+      },
+      {
+         "scheme":"doi",
+         "identifier":"10.5281/zenodo.4054058",
+         "relation":"isVersionOf"
+      }
+   ]
+}

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: DAISIE
 Type: Package
 Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction
-Version: 4.1.1
-Date: 2022-04-14
+Version: 4.2.0
+Date: 2022-05-24
 Depends: R (>= 3.5.0)
 biocViews:
 SystemRequirements: C++14
@@ -115,5 +115,4 @@ Encoding: UTF-8
 VignetteBuilder: knitr
 URL: https://github.com/rsetienne/DAISIE
 BugReports: https://github.com/rsetienne/DAISIE/issues
-RoxygenNote: 7.1.2
-LazyData: true
+RoxygenNote: 7.2.0

NAMESPACE

@@ -12,6 +12,7 @@ export(DAISIE_SR_ML)
 export(DAISIE_SR_ML_CS)
 export(DAISIE_SR_loglik_CS)
 export(DAISIE_SR_loglik_all)
+export(DAISIE_abm_factor)
 export(DAISIE_convertprobdist)
 export(DAISIE_dataprep)
 export(DAISIE_loglik_CS)

NEWS.md

@@ -1,3 +1,24 @@
+# DAISIE 4.2.0
+
+* Important bugfixes on estimation when data contains a lineage or a clade for 
+which (see `stac_key` vignette for details on each case):
+  * The divergence time is unknown but known to have occurred *after* a 
+  specific point in time in the island's existence.
+  * The divergence time is unknown but known to have occurred *before* a 
+  specific point in time in the island's existence (usually known when dated 
+  population level phylogenetic data is available).
+* Improve and finalise the continental island estimation scenario, when there 
+are initially species present on the island on inception. 
+* Improve tests related to the two points above.
+* Remove DAISIE time dependent estimation code for now to ease work on the remaining
+code. Development will proceed separately.
+* Rate computations in simulation are first calculated per lineage, to then be 
+calculated per island. This will be needed for time dependent estimation later on.
+* Improve zenodo metadata.
+* Add `methode = "odeint::adams_bashforth_moulton_X"` to the list of available
+numeric integrators for parameter estimation.
+
+
 # DAISIE 4.1.1
 
 * Correctly use `is.data.frame()` rather than `class(foo) == "data.frame"` to satisfy CRAN note.

R/DAISIE_ML1.R

@@ -11,18 +11,28 @@ DAISIE_loglik_all_choosepar <- function(trparsopt,
                                         CS_version = 1,
                                         abstolint = 1E-16,
                                         reltolint = 1E-10) {
-  if (sum(idparsnoshift %in% (6:10)) != 5) {
+  all_no_shift <- 6:10
+  if (max(idparsopt,-Inf) <= 6 &&
+      max(idparsfix,-Inf) <= 6 &&
+      (6 %in% idparsopt || 6 %in% idparsfix)) {
+    idparsnoshift <- 7:11
+    all_no_shift <- 7:11
+  }
+  if (sum(idparsnoshift %in% (all_no_shift)) != 5) {
     trpars1 <- rep(0, 11)
   } else {
-    trpars1 <- rep(0, 5)
-    trparsfix <- trparsfix[-which(idparsfix == 11)]
-    idparsfix <- idparsfix[-which(idparsfix == 11)]
+    trpars1 <- rep(0, 6)
+    prop_type2_present <- which(idparsfix == 11)
+    if(!is.null(prop_type2_present)) {
+       trparsfix <- trparsfix[-prop_type2_present]
+       idparsfix <- idparsfix[-prop_type2_present]
+    }
   }
   trpars1[idparsopt] <- trparsopt
   if (length(idparsfix) != 0) {
     trpars1[idparsfix] <- trparsfix
   }
-  if (sum(idparsnoshift %in% (6:10)) != 5) {
+  if (sum(idparsnoshift %in% all_no_shift) != 5) {
     trpars1[idparsnoshift] <- trpars1[idparsnoshift - 5]
   }
   if (max(trpars1) > 1 | min(trpars1) < 0) {
@@ -31,14 +41,14 @@ DAISIE_loglik_all_choosepar <- function(trparsopt,
     pars1 <- trpars1 / (1 - trpars1)
     if (pars2[6] > 0) {
       pars1 <- DAISIE_eq(datalist, pars1, pars2[-5])
-      if (sum(idparsnoshift %in% (6:10)) != 5) {
+      if (sum(idparsnoshift %in% all_no_shift) != 5) {
         pars1[idparsnoshift] <- pars1[idparsnoshift - 5]
       }
     }
     if (min(pars1) < 0) {
       loglik <- -Inf
     } else {
-      loglik <- DAISIE::DAISIE_loglik_all(
+      loglik <- DAISIE_loglik_all(
         pars1 = pars1,
         pars2 = pars2,
         datalist = datalist,
@@ -180,7 +190,35 @@ DAISIE_ML1 <- function(
     "lambda_a2",
     "prop_type2"
   )
+  all_no_shift <- 6:10
+  max_idpars <- 11
 
+  if (max(idparsopt, -Inf) <= 6 &&
+      max(idparsfix, -Inf) <= 6 &&
+      (6 %in% idparsopt || 6 %in% idparsfix)) {
+    max_idpars <- 12
+    idparsnoshift <- 7:11
+    all_no_shift <- 7:11
+    namepars <- c(
+      "lambda_c",
+      "mu",
+      "K",
+      "gamma",
+      "lambda_a",
+      "prob_init_pres",
+      "lambda_c2",
+      "mu2",
+      "K2",
+      "gamma2",
+      "lambda_a2",
+      "prop_type2"
+    )
+    nc <- NA
+    names(nc) <- "prob_init_pres"
+    out2err <- add_column_to_dataframe(df = out2err,
+                                       position = 'lambda_a',
+                                       column_to_insert = nc)
+  }
   if (length(namepars[idparsopt]) == 0) {
     optstr <- "nothing"
   } else {
@@ -194,28 +232,27 @@ DAISIE_ML1 <- function(
     fixstr <- namepars[idparsfix]
   }
   cat("You are fixing", fixstr, "\n")
-  if (sum(idparsnoshift %in% (6:10)) != 5) {
+
+  if (sum(idparsnoshift %in% (all_no_shift)) != 5) {
     noshiftstring <- namepars[idparsnoshift]
     cat("You are not shifting", noshiftstring, "\n")
   }
   idpars <- sort(c(idparsopt, idparsfix, idparsnoshift, idparseq))
-  if (!any(idpars == 11)) {
-    idpars <- c(idpars, 11)
-    idparsfix <- c(idparsfix, 11)
-    parsfix <- c(parsfix, 0)
-  }
+
   missnumspec <- unlist(lapply(datalist, function(list) {list$missing_species})) # nolint
   if (sum(missnumspec) > (res - 1)) {
     cat(
       "The number of missing species is too large relative to the
         resolution of the ODE.\n")
     return(out2err)
   }
-  if (length(idpars) != 11) {
-    cat("You have too many parameters to be optimized or fixed.\n")
+
+  if ((length(idpars) != max(idpars))) {
+    cat("The parameters to be optimized and/or fixed are incoherent.\n")
     return(out2err)
   }
-  if ((prod(idpars == (1:11)) != 1) || # nolint
+
+  if ((!all(idpars == 1:max(idpars))) || # nolint
       (length(initparsopt) != length(idparsopt)) ||
       (length(parsfix) != length(idparsfix))) {
     cat("The parameters to be optimized and/or fixed are incoherent.\n")
@@ -313,10 +350,10 @@ DAISIE_ML1 <- function(
   MLpars <- MLtrpars / (1 - MLtrpars)
   ML <- as.numeric(unlist(out$fvalues))
 
-  if (sum(idparsnoshift %in% (6:10)) != 5) {
-    MLpars1 <- rep(0, 10)
+  if (sum(idparsnoshift %in% (all_no_shift)) != 5) {
+    MLpars1 <- rep(0, 11)
   } else {
-    MLpars1 <- rep(0, 5)
+    MLpars1 <- rep(0, 6)
   }
   MLpars1[idparsopt] <- MLpars
   if (length(idparsfix) != 0) {
@@ -331,7 +368,7 @@ DAISIE_ML1 <- function(
     MLpars1[3] <- Inf
   }
 
-  if (sum(idparsnoshift %in% (6:10)) != 5) {
+  if (sum(idparsnoshift %in% (all_no_shift)) != 5) {
     if (length(idparsnoshift) != 0) {
       MLpars1[idparsnoshift] <- MLpars1[idparsnoshift - 5]
     }
@@ -355,9 +392,7 @@ DAISIE_ML1 <- function(
       conv = unlist(out$conv)
     )
     s1 <- sprintf(
-      "Maximum likelihood parameter estimates: lambda_c: %f, mu: %f, K: %f,
-      gamma: %f, lambda_a: %f, lambda_c2: %f, mu2: %f, K2: %f, gamma2: %f,
-      lambda_a2: %f, prop_type2: %f",
+      "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n lambda_c2: %f\n mu2: %f\n K2: %f\n gamma2: %f\n lambda_a2: %f\n prop_type2: %f",
       MLpars1[1],
       MLpars1[2],
       MLpars1[3],
@@ -370,6 +405,27 @@ DAISIE_ML1 <- function(
       MLpars1[10],
       MLpars1[11]
     )
+  } else if (all(all_no_shift == 7:11)) {
+    out2 <- data.frame(
+      lambda_c = MLpars1[1],
+      mu = MLpars1[2],
+      K = MLpars1[3],
+      gamma = MLpars1[4],
+      lambda_a = MLpars1[5],
+      prob_init_pres = MLpars1[6],
+      loglik = ML,
+      df = length(initparsopt),
+      conv = unlist(out$conv)
+    )
+    s1 <- sprintf(
+      "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n prob_init_pres: %f",
+      MLpars1[1],
+      MLpars1[2],
+      MLpars1[3],
+      MLpars1[4],
+      MLpars1[5],
+      MLpars1[6]
+    )
   } else {
     out2 <- data.frame(
       lambda_c = MLpars1[1],
@@ -382,8 +438,7 @@ DAISIE_ML1 <- function(
       conv = unlist(out$conv)
     )
     s1 <- sprintf(
-      "Maximum likelihood parameter estimates: lambda_c: %f, mu: %f, K: %f,
-      gamma: %f, lambda_a: %f",
+      "Maximum likelihood parameter estimates:\n lambda_c: %f\n mu: %f\n K: %f\n gamma: %f\n lambda_a: %f\n",
       MLpars1[1],
       MLpars1[2],
       MLpars1[3],