diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 17d90db..607d193 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -2,7 +2,20 @@ # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: + branches: + - main + - master + - develop + - richel + - new-ddmodels + pull_request: + branches: + - main + - master + - develop + - richel + - new-ddmodels name: R-CMD-check diff --git a/NAMESPACE b/NAMESPACE index 1f193ec..2512535 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(L2brts) export(L2phylo) export(bd_ML) export(bd_loglik) +export(both_rates_vary) export(brts2phylo) export(conv) export(dd_KI_ML) @@ -19,6 +20,7 @@ export(dd_SR_loglik) export(dd_SR_sim) export(dd_loglik) export(dd_sim) +export(get_Kprime) export(optimizer) export(phylo2L) export(rng_respecting_sample) @@ -28,5 +30,6 @@ export(simplex) export(td_sim) export(transform_pars) export(untransform_pars) +export(what_is_this_ddmodel) import(Rcpp) useDynLib(DDD) diff --git a/R/dd_KI_loglik.R b/R/dd_KI_loglik.R index 5efa169..9b4122e 100644 --- a/R/dd_KI_loglik.R +++ b/R/dd_KI_loglik.R @@ -85,11 +85,11 @@ #' @keywords models #' @examples #' -#' pars1 = c(0.25,0.12,25.51,1.0,0.16,8.61,9.8) -#' pars2 = c(200,1,0,18.8,1,2) -#' missnumspec = 0 -#' brtsM = c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) -#' brtsS = c(9.6,8.6,7.4,4.9,2.5) +#' pars1 <- c(0.25,0.12,25.51,1.0,0.16,8.61,9.8) +#' pars2 <- c(200,1,0,18.8,1,2) +#' missnumspec <- 0 +#' brtsM <- c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) +#' brtsS <- c(9.6,8.6,7.4,4.9,2.5) #' dd_KI_loglik(pars1,pars2,brtsM,brtsS,missnumspec) #' #' @export dd_KI_loglik @@ -102,46 +102,46 @@ dd_KI_loglik <- function(pars1, { if(length(pars2) == 4) { - pars2[5] = 0 - pars2[6] = 2 - pars2[7] = 1 + pars2[5] <- 0 + pars2[6] <- 2 + pars2[7] <- 1 } if(is.na(pars2[7])) { - pars2[7] = 0 + pars2[7] <- 0 } #pars2 <- as.data.frame(rbind(pars2)) - names(pars2) = c('lx','ddep','cond','t_split','verbose','soc','corr') - tinn = -abs(pars1[7]) - soc = pars2[6] + names(pars2) <- c('lx','ddep','cond','t_split','verbose','soc','corr') + tinn <- -abs(pars1[7]) + soc <- pars2[6] #pars2 <- as.vector(pars2) # order branching times - brts = -sort(abs(c(brtsM,brtsS)),decreasing = TRUE) + brts <- -sort(abs(c(brtsM,brtsS)),decreasing = TRUE) if(sum(brts == 0) == 0) { - brts[length(brts) + 1] = 0 + brts[length(brts) + 1] <- 0 } - S = length(brts) + (soc - 2) - brtsM = -sort(abs(brtsM),decreasing = TRUE) - brtsS = -sort(abs(brtsS),decreasing = TRUE) + S <- length(brts) + (soc - 2) + brtsM <- -sort(abs(brtsM),decreasing = TRUE) + brtsS <- -sort(abs(brtsS),decreasing = TRUE) tinn <- -abs(pars1[7]) if(!any(brtsM == 0)) { - brtsM = c(brtsM,0) + brtsM <- c(brtsM,0) } if(!any(brtsS == 0)) { - brtsS = c(brtsS,0) + brtsS <- c(brtsS,0) } if(!any(brtsS == tinn)) { - brtsS = c(tinn,brtsS) + brtsS <- c(tinn,brtsS) } # avoid coincidence of branching time and key innovation time - if(sum(abs(brtsM - tinn) < 1E-14) == 1) { tinn = tinn - 1E-8 } + if(sum(abs(brtsM - tinn) < 1E-14) == 1) { tinn <- tinn - 1E-8 } - ka = sum(brtsM < tinn) + ka <- sum(brtsM < tinn) brtsMtinn <- sort(c(tinn,brtsM),decreasing = FALSE) kvec <- c(soc:(soc + ka - 1), (soc + ka - 2):(soc + length(brtsMtinn) - 4), @@ -174,15 +174,15 @@ dd_KI_loglik <- function(pars1, if(pars2[5]) { - s1 = sprintf('Parameters: %f %f %f %f %f %f %f, ',pars1[1],pars1[2],pars1[3],pars1[4],pars1[5],pars1[6],pars1[7]) - s2 = sprintf('Loglikelihood: %f',loglik) + s1 <- sprintf('Parameters: %f %f %f %f %f %f %f, ',pars1[1],pars1[2],pars1[3],pars1[4],pars1[5],pars1[6],pars1[7]) + s2 <- sprintf('Loglikelihood: %f',loglik) cat(s1,s2,"\n",sep = "") utils::flush.console() } - loglik = as.numeric(loglik) + loglik <- as.numeric(loglik) if(is.nan(loglik) || is.na(loglik)) { - loglik = -Inf + loglik <- -Inf } return(loglik) } @@ -309,7 +309,7 @@ dd_KI_logliknorm <- function(brts_k_list, { if(cond == 0 || loglik == -Inf) { - logliknorm = 0 + logliknorm <- 0 } else { if(length(brts_k_list) > 2 && cond == 1) { @@ -326,101 +326,101 @@ dd_KI_logliknorm <- function(brts_k_list, tpres <- 0 # COMPUTE NORMALIZATION # compute survival probability of clade S - lx = lx_list[[2]] - nx = -1:lx + lx <- lx_list[[2]] + nx <- -1:lx lambdamu_nk <- lambdamu(nx,pars = c(pars1_list[[2]],0),ddep = ddep) lavec <- lambdamu_nk[[1]] muvec <- lambdamu_nk[[2]] - probs = rep(0,lx) # probs[1] = extinction probability - probs[2] = 1 # clade S starts with one species + probs <- rep(0,lx) # probs[1] = extinction probability + probs[2] <- 1 # clade S starts with one species if(methode != 'analytical') { - m1 = lavec[1:lx] * nx[1:lx] - m2 = muvec[3:(lx + 2)] * nx[3:(lx + 2)] - m3 = (lavec[2:(lx + 1)] + muvec[2:(lx + 1)]) * nx[2:(lx + 1)] + m1 <- lavec[1:lx] * nx[1:lx] + m2 <- muvec[3:(lx + 2)] * nx[3:(lx + 2)] + m3 <- (lavec[2:(lx + 1)] + muvec[2:(lx + 1)]) * nx[2:(lx + 1)] if (startsWith(methode, 'odeint::')) { - probs = dd_logliknorm1_odeint(probs, c(tinn,tpres), c(m1,m2,m3), abstol, reltol, methode) + probs <- dd_logliknorm1_odeint(probs, c(tinn,tpres), c(m1,m2,m3), abstol, reltol, methode) } else { - y = deSolve::ode(probs,c(tinn,tpres),dd_logliknorm_rhs1,c(m1,m2,m3),rtol = reltol,atol = abstol,method = methode) - probs = y[2,2:(lx+1)] + y <- deSolve::ode(probs,c(tinn,tpres),dd_logliknorm_rhs1,c(m1,m2,m3),rtol = reltol,atol = abstol,method = methode) + probs <- y[2,2:(lx+1)] } } else { - probs = dd_loglik_M(pars1_list[[2]],lx,0,ddep,tt = abs(tpres - tinn),probs) + probs <- dd_loglik_M(pars1_list[[2]],lx,0,ddep,tt = abs(tpres - tinn),probs) } - PS = 1 - probs[1] + PS <- 1 - probs[1] # compute survival probability of clade M - lx = lx_list[[1]] + lx <- lx_list[[1]] n <- -1:lx - nx1 = rep(n,lx + 2) - dim(nx1) = c(lx + 2,lx + 2) # row index = number of species in first group - nx2 = t(nx1) # column index = number of species in second group - nxt = nx1 + nx2 + nx1 <- rep(n,lx + 2) + dim(nx1) <- c(lx + 2,lx + 2) # row index = number of species in first group + nx2 <- t(nx1) # column index = number of species in second group + nxt <- nx1 + nx2 lambdamu_n <- lambdamu2(nxt,pars1_list[[1]],ddep) lavec <- lambdamu_n[[1]] muvec <- lambdamu_n[[2]] - probs = matrix(0,lx,lx) + probs <- matrix(0,lx,lx) # probs[1,1] = probability of extinction of both lineages # sum(probs[1:lx,1]) = probability of extinction of second lineage - probs[2,2] = 1 # clade M starts with two species + probs[2,2] <- 1 # clade M starts with two species # STEP 1: integrate from tcrown to tinn - dim(probs) = c(lx,lx) + dim(probs) <- c(lx,lx) if(methode != 'analytical') { - m1 = lavec[1:lx,2:(lx+1)] * nx1[1:lx,2:(lx+1)] - m2 = muvec[3:(lx+2),2:(lx+1)] * nx1[3:(lx+2),2:(lx+1)] - ma = lavec[2:(lx+1),2:(lx+1)] + muvec[2:(lx+1),2:(lx+1)] - m3 = ma * nx1[2:(lx+1),2:(lx+1)] - m4 = lavec[2:(lx+1),1:lx] * nx2[2:(lx+1),1:lx] - m5 = muvec[2:(lx+1),3:(lx+2)] * nx2[2:(lx+1),3:(lx+2)] - m6 = ma * nx2[2:(lx+1),2:(lx+1)] + m1 <- lavec[1:lx,2:(lx+1)] * nx1[1:lx,2:(lx+1)] + m2 <- muvec[3:(lx+2),2:(lx+1)] * nx1[3:(lx+2),2:(lx+1)] + ma <- lavec[2:(lx+1),2:(lx+1)] + muvec[2:(lx+1),2:(lx+1)] + m3 <- ma * nx1[2:(lx+1),2:(lx+1)] + m4 <- lavec[2:(lx+1),1:lx] * nx2[2:(lx+1),1:lx] + m5 <- muvec[2:(lx+1),3:(lx+2)] * nx2[2:(lx+1),3:(lx+2)] + m6 <- ma * nx2[2:(lx+1),2:(lx+1)] if (startsWith(methode, "odeint::")) { - probs = dd_logliknorm2_odeint(probs, c(tcrown,tinn), list(m1,m2,m3,m4,m5,m6), reltol, abstol, methode) + probs <- dd_logliknorm2_odeint(probs, c(tcrown,tinn), list(m1,m2,m3,m4,m5,m6), reltol, abstol, methode) } else { - dim(probs) = c(lx*lx,1) - y = deSolve::ode(probs,c(tcrown,tinn),dd_logliknorm_rhs2,list(m1,m2,m3,m4,m5,m6),rtol = reltol,atol = abstol, method = "ode45") - probs = y[2,2:(lx * lx + 1)] + dim(probs) <- c(lx*lx,1) + y <- deSolve::ode(probs,c(tcrown,tinn),dd_logliknorm_rhs2,list(m1,m2,m3,m4,m5,m6),rtol = reltol,atol = abstol, method = "ode45") + probs <- y[2,2:(lx * lx + 1)] } } else { - probs = dd_loglik_M2(pars = pars1_list[[1]],lx = lx,ddep = ddep,tt = abs(tinn - tcrown),p = probs) + probs <- dd_loglik_M2(pars = pars1_list[[1]],lx = lx,ddep = ddep,tt = abs(tinn - tcrown),p = probs) } - dim(probs) = c(lx,lx) - probs[1,1:lx] = 0 - probs[1:lx,1] = 0 + dim(probs) <- c(lx,lx) + probs[1,1:lx] <- 0 + probs[1:lx,1] <- 0 # STEP 2: transformation at tinn - nx1a = nx1[2:(lx+1),2:(lx+1)] - nx2a = nx2[2:(lx+1),2:(lx+1)] - probs = probs * nx1a/(nx1a+nx2a) - probs = rbind(probs[2:lx,1:lx], rep(0,lx)) + nx1a <- nx1[2:(lx+1),2:(lx+1)] + nx2a <- nx2[2:(lx+1),2:(lx+1)] + probs <- probs * nx1a/(nx1a+nx2a) + probs <- rbind(probs[2:lx,1:lx], rep(0,lx)) # STEP 3: integrate from tinn to tpres if(methode != 'analytical') { if (startsWith(methode, "odeint::")) { - probs = dd_logliknorm2_odeint(probs, c(tinn, tpres), list(m1,m2,m3,m4,m5,m6), reltol, abstol, 'odeint::runge_kutta_fehlberg78') + probs <- dd_logliknorm2_odeint(probs, c(tinn, tpres), list(m1,m2,m3,m4,m5,m6), reltol, abstol, 'odeint::runge_kutta_fehlberg78') } else { - dim(probs) = c(lx*lx,1) - y = deSolve::ode(probs,c(tinn, tpres),dd_logliknorm_rhs2,list(m1,m2,m3,m4,m5,m6),rtol = reltol,atol = abstol, method = "ode45") - probs = y[2,2:(lx * lx + 1)] - dim(probs) = c(lx,lx) + dim(probs) <- c(lx*lx,1) + y <- deSolve::ode(probs,c(tinn, tpres),dd_logliknorm_rhs2,list(m1,m2,m3,m4,m5,m6),rtol = reltol,atol = abstol, method = "ode45") + probs <- y[2,2:(lx * lx + 1)] + dim(probs) <- c(lx,lx) } } else { - probs = dd_loglik_M2(pars = pars1_list[[1]],lx = lx,ddep = ddep,tt = abs(tpres - tinn),p = probs) + probs <- dd_loglik_M2(pars = pars1_list[[1]],lx = lx,ddep = ddep,tt = abs(tpres - tinn),p = probs) } - dim(probs) = c(lx,lx) - PM12 = sum(probs[2:lx,2:lx]) - PM2 = sum(probs[1,2:lx]) + dim(probs) <- c(lx,lx) + PM12 <- sum(probs[2:lx,2:lx]) + PM2 <- sum(probs[1,2:lx]) #print(log(2) + log(PM12 + PS * PM2)) #print(log(2) + (log(PM12 + PM2) + log(PS))) #print(log(2) + (log(PM12) + log(PS))) - logliknorm = log(2) + (cond == 1) * log(PM12 + PS * PM2) + + logliknorm <- log(2) + (cond == 1) * log(PM12 + PS * PM2) + (cond == 4) * (log(PM12 + PM2) + log(PS)) + (cond == 5) * (log(PM12) + log(PS)) } @@ -538,7 +538,7 @@ check_for_impossible_pars <- function(pars1, #' the dynamics of main clades, potentially accompanied by a #' shift in parameters. #' @inheritParams dd_KI_loglik -#' @param pars1_list list of paramater sets one for each rate regime (subclade). +#' @param pars1_list list of parameter sets one for each rate regime (subclade). #' The parameters are: lambda (speciation rate), mu (extinction rate), and K #' (clade-level carrying capacity). #' @param brts_k_list list of matrices, one for each rate regime (subclade). Each diff --git a/R/dd_ML.R b/R/dd_ML.R index e6ebd05..6402bad 100644 --- a/R/dd_ML.R +++ b/R/dd_ML.R @@ -1,25 +1,20 @@ initparsoptdefault = function(ddmodel,brts,missnumspec) { - if(ddmodel < 5) - { - return(c(0.2,0.1,2 * (length(brts) + missnumspec)^(ddmodel != 2.3))) + if (both_rates_vary(ddmodel)) { + return(c(0.2, 0.1, 2 * (length(brts) + missnumspec), 0.01)) } else { - return(c(0.2,0.1,2 * (length(brts) + missnumspec),0.01)) + return(c(0.2, 0.1, 2 * (length(brts) + missnumspec) ^ (ddmodel != 2.3))) } } -parsfixdefault = function(ddmodel,brts,missnumspec,idparsopt) -{ - if(ddmodel < 5) - { - return(c(0.2,0.1,2*(length(brts) + missnumspec))[-idparsopt]) +parsfixdefault = function(ddmodel, brts, missnumspec, idparsopt) { + if (both_rates_vary(ddmodel)) { + return(c(0.2, 0.1, 2 * (length(brts) + missnumspec), 0)[-idparsopt]) } else { - return(c(0.2,0.1,2*(length(brts) + missnumspec),0)[-idparsopt]) + return(c(0.2, 0.1, 2 * (length(brts) + missnumspec))[-idparsopt]) } } - - #' Maximization of the loglikelihood under a diversity-dependent #' diversification model #' @@ -62,19 +57,41 @@ parsfixdefault = function(ddmodel,brts,missnumspec,idparsopt) #' \code{ddmodel == 1.5} : positive and negative dependence in speciation rate #' with parameter K' (= diversity where speciation = 0); lambda = lambda0 * #' S/K' * (1 - S/K') where S is species richness\cr -#' \code{ddmodel == 2} : exponential dependence in speciation rate with parameter +#' \code{ddmodel == 2} : exponential dependence (power function) in speciation rate with parameter #' K (= diversity where speciation = extinction)\cr -#' \code{ddmodel == 2.1} : variant of exponential dependence in speciation rate +#' \code{ddmodel == 2.1} : variant of exponential dependence (power function) in speciation rate #' with offset at infinity\cr #' \code{ddmodel == 2.2} : 1/n dependence in speciation rate\cr -#' \code{ddmodel == 2.3} : exponential dependence in speciation rate with parameter x (= +#' \code{ddmodel == 2.3} : exponential dependence (power function) in speciation rate with parameter x (= #' exponent)\cr #' \code{ddmodel == 3} : linear dependence in extinction rate \cr -#' \code{ddmodel == 4} : exponential dependence in extinction rate \cr +#' \code{ddmodel == 4} : exponential dependence (power function) in extinction rate \cr #' \code{ddmodel == 4.1} : variant of exponential dependence in extinction rate #' with offset at infinity \cr #' \code{ddmodel == 4.2} : 1/n dependence in extinction rate with offset at infinity \cr \code{ddmodel == 5} : linear #' dependence in speciation and extinction rate \cr +#' \code{ddmodel == 5} : linear dependence in speciation and +#' extinction rate \cr +#' \code{ddmodel == 6} : linear dependence in speciation rate, exponential +#' dependence (power function) in extinction rate \cr +#' \code{ddmodel == 7} : exponential dependence (power function) in speciation +#' and extinction rate \cr +#' \code{ddmodel == 8} : exponential dependence (power function) in speciation rate, +#' linear dependence in extinction rate \cr +#' \code{ddmodel == 9} : exponential dependence (exponential function) in speciation, +#' constant-rate extinction \cr +#' \code{ddmodel == 10} : constant-rate speciation, exponential dependence +#' (exponential function) in extinction \cr +#' \code{ddmodel == 11} : linear dependence in speciation, exponential +#' dependence (exponential function) in extinction\cr +#' \code{ddmodel == 12} : exponential dependence (exponential function) in +#' speciation and extinction\cr +#' \code{ddmodel == 13} : exponential dependence (exponential function) in +#' speciation, linear dependence in extinction \cr +#' \code{ddmodel == 14} : exponential dependence (exponential function) in +#' speciation, exponential dependence (power function) in extinction \cr +#' \code{ddmodel == 15} : exponential dependence (power function) in +#' speciation, exponential dependence (exponential function) in extinction \cr #' @param missnumspec The number of species that are in the clade but missing #' in the phylogeny #' @param cond Conditioning: \cr @@ -136,9 +153,9 @@ dd_ML = function( brts, initparsopt = initparsoptdefault(ddmodel,brts,missnumspec), idparsopt = 1:length(initparsopt), - idparsfix = (1:(3 + (ddmodel == 5)))[-idparsopt], + idparsfix = (1:(3 + both_rates_vary(ddmodel)))[-idparsopt], parsfix = parsfixdefault(ddmodel,brts,missnumspec,idparsopt), - res = 10*(1+length(brts)+missnumspec), + res = 10 * (1 + length(brts) + missnumspec), ddmodel = 1, missnumspec = 0, cond = 1, @@ -150,39 +167,50 @@ dd_ML = function( optimmethod = 'subplex', num_cycles = 1, methode = 'analytical', - verbose = FALSE) -{ + verbose = FALSE + ) { #options(warn = -1) - if(length(tol) != 3) - { + if(length(tol) != 3) { stop('Please specify a tolerance vector with three values') } + if (both_rates_vary(ddmodel)) { + output_error <- data.frame(lambda = -1,mu = -1,K = -1, r = -1, loglik = -1, df = -1, conv = -1) + } else { + output_error <- data.frame(lambda = -1,mu = -1,K = -1, loglik = -1, df = -1, conv = -1) + } + + if (ddmodel > 5) { + if (methode == "analytical" || cond == 3) { + stop("Sorry, ddmodel options > 5 have not been developed for method = \"analytical\" or cond = 3.") + } + } + brts = sort(abs(as.numeric(brts)),decreasing = TRUE) - if(is.numeric(brts) == FALSE) - { + if (is.numeric(brts) == FALSE) { cat("The branching times should be numeric.\n") - out2 = data.frame(lambda = -1,mu = -1,K = -1, loglik = -1, df = -1, conv = -1) - if(ddmodel == 5) {out2 = data.frame(lambda = -1,mu = -1,K = -1, r = -1, loglik = -1, df = -1, conv = -1)} + out2 <- output_error } else { idpars = sort(c(idparsopt,idparsfix)) - if((prod(idpars == (1:(3 + (ddmodel == 5)))) != 1) || (length(initparsopt) != length(idparsopt)) || (length(parsfix) != length(idparsfix))) + if (!all(idpars == (1:(3 + both_rates_vary(ddmodel)))) || (length(initparsopt) != length(idparsopt)) || (length(parsfix) != length(idparsfix))) { cat("The parameters to be optimized and/or fixed are incoherent.\n") - out2 = data.frame(lambda = -1,mu = -1,K = -1, loglik = -1, df = -1, conv = -1) - if(ddmodel == 5) {out2 = data.frame(lambda = -1,mu = -1,K = -1, r = -1, loglik = -1, df = -1, conv = -1)} + out2 <- output_error } else { - namepars = c("lambda","mu","K") - if(ddmodel == 5) {namepars = namepars = c("lambda","mu","K","r")} + if (both_rates_vary(ddmodel)) { + namepars = c("lambda","mu","K","r") + } else { + namepars = c("lambda","mu","K") + } if(length(namepars[idparsopt]) == 0) { optstr = "nothing" } else { optstr = namepars[idparsopt] } cat("You are optimizing",optstr,"\n") if(length(namepars[idparsfix]) == 0) { fixstr = "nothing" } else { fixstr = namepars[idparsfix] } cat("You are fixing",fixstr,"\n") cat("Optimizing the likelihood - this may take a while.","\n") utils::flush.console() - trparsopt = initparsopt/(1 + initparsopt) - trparsopt[which(initparsopt == Inf)] = 1 - trparsfix = parsfix/(1 + parsfix) - trparsfix[which(parsfix == Inf)] = 1 + trparsopt = initparsopt / (1 + initparsopt) + trparsopt[which(initparsopt == Inf)] <- 1 # fix NaN + trparsfix = parsfix / (1 + parsfix) + trparsfix[which(parsfix == Inf)] <- 1 # fix NaN pars2 = c(res,ddmodel,cond,btorph,verbose,soc,tol,maxiter) optimpars = c(tol,maxiter) initloglik = dd_loglik_choosepar(trparsopt = trparsopt,trparsfix = trparsfix,idparsopt = idparsopt,idparsfix = idparsfix,pars2 = pars2,brts = brts,missnumspec = missnumspec, methode = methode) @@ -191,39 +219,58 @@ dd_ML = function( if(initloglik == -Inf) { cat("The initial parameter values have a likelihood that is equal to 0 or below machine precision. Try again with different initial values.\n") - out2 = data.frame(lambda = -1,mu = -1,K = -1, loglik = -1, df = -1, conv = -1) - if(ddmodel == 5) {out2 = data.frame(lambda = -1,mu = -1,K = -1, r = -1, loglik = -1, df = -1, conv = -1)} + out2 <- output_error } else { #code up to DDD v1.6: out = optimx2(trparsopt,dd_loglik_choosepar,hess=NULL,method = "Nelder-Mead",hessian = FALSE,control = list(maximize = TRUE,abstol = pars2[8],reltol = pars2[7],trace = 0,starttests = FALSE,kkt = FALSE),trparsfix = trparsfix,idparsopt = idparsopt,idparsfix = idparsfix,brts = brts, pars2 = pars2,missnumspec = missnumspec) #out = dd_simplex(trparsopt,idparsopt,trparsfix,idparsfix,pars2,brts,missnumspec) - out = optimizer(optimmethod = optimmethod,optimpars = optimpars,fun = dd_loglik_choosepar,trparsopt = trparsopt,trparsfix = trparsfix,idparsopt = idparsopt,idparsfix = idparsfix,pars2 = pars2,brts = brts, missnumspec = missnumspec, methode = methode, num_cycles = num_cycles) - if(out$conv != 0) - { + out = optimizer( + optimmethod = optimmethod, + optimpars = optimpars, + fun = dd_loglik_choosepar, + trparsopt = trparsopt, + trparsfix = trparsfix, + idparsopt = idparsopt, + idparsfix = idparsfix, + pars2 = pars2, + brts = brts, + missnumspec = missnumspec, + methode = methode, + num_cycles = num_cycles + ) + if (out$conv != 0) { cat("Optimization has not converged. Try again with different initial values.\n") - out2 = data.frame(lambda = -1,mu = -1,K = -1, loglik = -1, df = -1, conv = unlist(out$conv)) - if(ddmodel == 5) {out2 = data.frame(lambda = -1,mu = -1,K = -1, r = -1, loglik = -1, df = -1, conv = unlist(out$conv))} + out2 <- output_error } else { MLtrpars = as.numeric(unlist(out$par)) - MLpars = MLtrpars/(1-MLtrpars) - MLpars1 = rep(0,3) - if(ddmodel == 5) {MLpars1 = rep(0,4)} + MLpars = MLtrpars / (1 - MLtrpars) + if (both_rates_vary(ddmodel)) { + MLpars1 <- rep(0,4) + } else { + MLpars1 <- rep(0,3) + } MLpars1[idparsopt] = MLpars - if(length(idparsfix) != 0) { MLpars1[idparsfix] = parsfix } - if(MLpars1[3] > 10^7){MLpars1[3] = Inf} + if (length(idparsfix) != 0) { + MLpars1[idparsfix] = parsfix + } + if (MLpars1[3] > 10 ^ 7) { + MLpars1[3] = Inf + } ML = as.numeric(unlist(out$fvalues)) - out2 = data.frame(lambda = MLpars1[1],mu = MLpars1[2],K = MLpars1[3], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) - s1 = sprintf('Maximum likelihood parameter estimates: lambda: %f, mu: %f, K: %f',MLpars1[1],MLpars1[2],MLpars1[3]) - if(ddmodel == 5) - { - s1 = sprintf('%s, r: %f',s1,MLpars1[4]) - out2 = data.frame(lambda = MLpars1[1],mu = MLpars1[2],K = MLpars1[3], r = MLpars1[4], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) + if (both_rates_vary(ddmodel)) { + s1 <- sprintf('Maximum likelihood parameter estimates: lambda: %f, mu: %f, K: %f, r: %f', MLpars1[1], MLpars1[2], MLpars1[3], MLpars1[4]) + out2 <- data.frame(lambda = MLpars1[1], mu = MLpars1[2], K = MLpars1[3], r = MLpars1[4], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) + } else { + s1 <- sprintf('Maximum likelihood parameter estimates: lambda: %f, mu: %f, K: %f', MLpars1[1], MLpars1[2], MLpars1[3]) + out2 <- data.frame(lambda = MLpars1[1], mu = MLpars1[2], K = MLpars1[3], loglik = ML, df = length(initparsopt), conv = unlist(out$conv)) + } + if(out2$conv != 0 & changeloglikifnoconv == T) { + out2$loglik = -Inf } - if(out2$conv != 0 & changeloglikifnoconv == T) { out2$loglik = -Inf } s2 = sprintf('Maximum loglikelihood: %f',ML) cat(paste("\n",s1,"\n",s2,"\n",sep = '')) } } } } - return(invisible(out2)) + return(invisible(out2)) } diff --git a/R/dd_MS_ML.R b/R/dd_MS_ML.R index 8be2fd8..ca8f83a 100644 --- a/R/dd_MS_ML.R +++ b/R/dd_MS_ML.R @@ -138,7 +138,7 @@ dd_MS_ML = function(brtsM, changeloglikifnoconv = FALSE, optimmethod = 'subplex', num_cycles = 1, - methode = 'ode45', + methode = 'odeint::runge_kutta_cash_karp54', correction = FALSE, verbose = FALSE) { diff --git a/R/dd_loglik.R b/R/dd_loglik.R index 335fa77..e9ac301 100644 --- a/R/dd_loglik.R +++ b/R/dd_loglik.R @@ -10,15 +10,35 @@ # . ddep == 1 : linear dependence in speciation rate with parameter K # . ddep == 1.3 : linear dependence in speciation rate with parameter K' # . ddep == 1.4 : positive and negative linear diversity-dependence in speciation rate with parameter K' -# . ddep == 2 : exponential dependence in speciation rate +# . ddep == 2 : exponential dependence (power function) in speciation rate # . ddep == 2.1: variant with offset at infinity # . ddep == 2.2: 1/n dependence in speciation rate -# . ddep == 2.3: exponential dependence in speciation rate with parameter x +# . ddep == 2.3: exponential dependence (power function) in speciation rate with parameter x # . ddep == 3 : linear dependence in extinction rate -# . ddep == 4 : exponential dependence in extinction rate +# . ddep == 4 : exponential dependence (power function) in extinction rate # . ddep == 4.1: variant with offset at infinity # . ddep == 4.2: 1/n dependence in speciation rate # . ddep == 5 : linear dependence in speciation and extinction rate +# . ddep == 6 : linear dependence in speciation, exponential dependence +# (power function) in extinction rate +# . ddep == 7 : exponential dependence (power function) in speciation +# and extinction rate +# . ddep == 8 : exponential dependence (power function) in speciation rate, +# linear dependence in extinction rate +# . ddep == 9 : exponential dependence (exponential function) in speciation, +# constant-rate extinction +# . ddep == 10 : constant-rate speciation, exponential dependence +# (exponential function) in extinction +# . ddep == 11 : linear dependence in speciation, exponential +# dependence (exponential function) in extinction +# . ddep == 12 : exponential dependence (exponential function) in +# speciation and extinction +# . ddep == 13 : exponential dependence (exponential function) in +# speciation, linear dependence in extinction +# . ddep == 14 : exponential dependence (exponential function) in +# speciation, exponential dependence (power function) in extinction +# . ddep == 15 : exponential dependence (power function) in +# speciation, exponential dependence (exponential function) in extinction # - pars2[3] = cond = conditioning # . cond == 0 : conditioning on stem or crown age # . cond == 1 : conditioning on stem or crown age and non-extinction of the phylogeny @@ -32,11 +52,9 @@ dd_loglik_test = function(pars1,pars2,brts,missnumspec,methode = 'analytical',rhs_func_name = 'dd_loglik_rhs') { - if(methode == 'analytical') - { + if (methode == 'analytical') { out = dd_loglik2(pars1,pars2,brts,missnumspec) - } else - { + } else { out = dd_loglik1(pars1,pars2,brts,missnumspec,methode = methode,rhs_func_name = rhs_func_name) } return(out) @@ -61,11 +79,9 @@ dd_loglik_test = function(pars1,pars2,brts,missnumspec,methode = 'analytical',rh #' \cr \cr \code{pars2[1]} sets the #' maximum number of species for which a probability must be computed. This #' must be larger than 1 + missnumspec + length(brts). -#' \cr \cr \code{pars2[2]} -#' sets the model of diversity-dependence: -#' \cr - \code{pars2[2] == 1} linear -#' dependence in speciation rate with parameter K (= diversity where speciation -#' = extinction) +#' \cr \cr \code{pars2[2]} sets the model of diversity-dependence: +#' \cr - \code{pars2[2] == 1} linear dependence in speciation rate with +#' parameter K (= diversity where speciation = extinction) #' \cr - \code{pars2[2] == 1.3} linear dependence in speciation #' rate with parameter K' (= diversity where speciation = 0) #' \cr - \code{pars2[2] == 1.4} : positive diversity-dependence in speciation rate @@ -74,28 +90,43 @@ dd_loglik_test = function(pars1,pars2,brts,missnumspec,methode = 'analytical',rh #' \cr - \code{pars2[2] == 1.5} : positive and negative diversity-dependence in #' speciation rate with parameter K' (= diversity where speciation = 0); lambda #' = lambda0 * S/K' * (1 - S/K') where S is species richness -#' \cr - \code{pars2[2] == 2} exponential dependence in speciation rate with +#' \cr - \code{pars2[2] == 2} exponential dependence (power function) in speciation rate with #' parameter K (= diversity where speciation = extinction) #' \cr - \code{pars2[2] -#' == 2.1} variant of exponential dependence in speciation rate with offset at +#' == 2.1} variant of exponential dependence (power function) in speciation rate with offset at #' infinity #' \cr - \code{pars2[2] == 2.2} 1/n dependence in speciation rate -#' \cr - \code{pars2[2] == 2.3} exponential dependence in speciation rate with +#' \cr - \code{pars2[2] == 2.3} exponential dependence (power function) in speciation rate with #' parameter x (= exponent) -#' \cr - \code{pars2[2] == 3} linear dependence in -#' extinction rate -#' \cr - \code{pars2[2] == 4} exponential dependence in -#' extinction rate -#' \cr - \code{pars2[2] == 4.1} variant of exponential -#' dependence in extinction rate with offset at infinity -#' \cr - \code{pars2[2] == -#' 4.2} 1/n dependence in extinction rate -#' \cr - \code{pars2[2] == 5} linear -#' dependence in speciation and extinction rate -#' \cr \cr \code{pars2[3]} sets -#' the conditioning: -#' \cr - \code{pars2[3] == 0} conditioning on stem or crown -#' age +#' \cr - \code{pars2[2] == 3} linear dependence in extinction rate +#' \cr - \code{pars2[2] == 4} exponential dependence (power function) in extinction rate +#' \cr - \code{pars2[2] == 4.1} variant of exponential dependence (power function) in extinction +#' rate with offset at infinity +#' \cr - \code{pars2[2] == 4.2} 1/n dependence in extinction rate +#' \cr - \code{pars2[2] == 5} linear dependence in speciation and extinction rate +#' \cr - \code{pars2[2] == 6} linear dependence in speciation rate, exponential +#' dependence (power function) in extinction rate \cr +#' \cr - \code{pars2[2] == 7} exponential dependence (power function) in speciation +#' and extinction rate \cr +#' \cr - \code{pars2[2] == 8} exponential dependence (power function) in speciation rate, +#' linear dependence in extinction rate \cr +#' \cr - \code{pars2[2] == 9} exponential dependence (exponential function) in speciation, +#' constant-rate extinction \cr +#' \cr - \code{pars2[2] == 10} constant-rate speciation, exponential dependence +#' (exponential function) in extinction \cr +#' \cr - \code{pars2[2] == 11} linear dependence in speciation, exponential +#' dependence (exponential function) in extinction\cr +#' \cr - \code{pars2[2] == 12} exponential dependence (exponential function) in +#' speciation and extinction\cr +#' \cr - \code{pars2[2] == 13} exponential dependence (exponential function) in +#' speciation, linear dependence in extinction \cr +#' \cr - \code{pars2[2] == 14} exponential dependence (exponential function) in +#' speciation, exponential dependence (power function) in extinction \cr +#' \cr - \code{pars2[2] == 15} exponential dependence (power function) in +#' speciation, exponential dependence (exponential function) in extinction \cr +#' +#' \cr \cr \code{pars2[3]} sets the conditioning: +#' \cr - \code{pars2[3] == 0} conditioning on stem or crown age #' \cr - \code{pars2[3] == 1} conditioning on stem or crown age and #' non-extinction of the phylogeny #' \cr - \code{pars2[3] == 2} conditioning on @@ -130,26 +161,241 @@ dd_loglik_test = function(pars1,pars2,brts,missnumspec,methode = 'analytical',rh #' @examples #' dd_loglik(pars1 = c(0.5,0.1,100), pars2 = c(100,1,1,1,0,2), brts = 1:10, missnumspec = 0) #' @export dd_loglik -dd_loglik = function(pars1,pars2,brts,missnumspec,methode = 'analytical') -{ - if(pars2[3] == 3) - { - rhs_func_name = 'dd_loglik_bw_rhs' - } else - { - rhs_func_name = 'dd_loglik_rhs' - } - if(methode == 'analytical') - { - out = dd_loglik2(pars1,pars2,brts,missnumspec) - } else - { - out = dd_loglik1(pars1,pars2,brts,missnumspec,methode = methode,rhs_func_name = rhs_func_name) - } - return(out) +dd_loglik = function(pars1, + pars2, + brts, + missnumspec = 0, + methode = "analytical" +) { + if(pars2[3] == 3) + { + rhs_func_name = 'dd_loglik_bw_rhs' + } else + { + rhs_func_name = 'dd_loglik_rhs' + } + if (methode == 'analytical') { + out = dd_loglik2(pars1,pars2,brts,missnumspec) + } else { + out = dd_loglik1(pars1,pars2,brts,missnumspec,methode = methode,rhs_func_name = rhs_func_name) + } + return(out) } -dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutta_cash_karp54',rhs_func_name = 'dd_loglik_rhs') +dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutta_cash_karp54',rhs_func_name = 'dd_loglik_rhs') { + # Unpack pars2 + if (length(pars2) == 4) { + pars2[5] = 0 + pars2[6] = 2 + } + ddep = pars2[2] + cond = pars2[3] + btorph = pars2[4] + verbose = pars2[5] + soc = pars2[6] + if(cond == 3) { + soc = 2 + } + # Unpack pars1 + la = pars1[1] + mu = pars1[2] + K = pars1[3] + r <- ifelse(both_rates_vary(ddep), pars1[4], 0) + Kprime <- get_Kprime(ddep, pars1) + is_speciation_linear <- ddep %in% c(1, 1.3, 5, 6, 11) + if (verbose) { + if (both_rates_vary(ddep)) { + cat("loglik for la0 =", la, "mu0 =", mu, "K =", K, "r =", r, "\n") + } else { + cat("loglik for la0 =", la, "mu0 =", mu, "K =", K, "\n") + } + } + + if ((ddep == 1) & ((mu == 0 & missnumspec == 0 & floor(K) != ceiling(K) & la > 0.05) | K == Inf)) { + loglik = bd_loglik(pars1[1:(2 + (K < Inf))],c(2*(mu == 0 & K < Inf),pars2[3:6]),brts,missnumspec) + } else { + abstol = 1e-10 + reltol = 1e-8 + brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) + if (sum(brts == 0) == 0) { + brts[length(brts) + 1] = 0 + } + S = length(brts) + (soc - 2) + if ((la == 0 | K == 0) & S > soc) { + loglik <- -Inf + return(loglik) + } + lx <- min(max(1 + missnumspec, 1 + Kprime), ceiling(pars2[1])) + ff <- lambdamu(1:lx, pars1, ddep) + lx <- min(which(ff[[2]]/ff[[1]] > 100), lx) + if (any(pars1 < 0)) { + if (verbose) cat('Model parameters cannot be negative.\n') + loglik = -Inf + } else if (la == 0) { + if (verbose) cat('la0 cannot be zero.\n') + loglik = -Inf + } else if (ddep %in% c(2, 2.1, 2.2, 4, 6, 7, 10:12, 14:15) && mu == 0) { + if (verbose) cat('mu0 cannot be exactly zero for this model.\n') + loglik = -Inf + } else if (ddep == 8 && mu == 0 && r == 0) { + if (verbose) cat('mu0 and r cannot both be zero for this model.\n') + loglik = -Inf + } else if (is_speciation_linear && Kprime < missnumspec + S) { + if (verbose) cat('K\' cannot be smaller than the nb of species in the tree.\n') + loglik = -Inf + } else if (la <= mu) { + if(verbose) cat("lambda0 cannot be smaller than mu0.\n") + loglik = -Inf + } else { + loglik = (btorph == 0) * lgamma(S) + if (cond != 3) { + qn_vec = rep(0,lx) + qn_vec[1] = 1 # change if other species at stem/crown age + for (k in 2:(S + 2 - soc)) { + k1 = k + (soc - 2) + y <- dd_integrate( + initprobs = qn_vec, + tvec = brts[(k-1):k], + rhs_func = rhs_func_name, + pars = c(pars1,k1,ddep), + rtol = reltol, + atol = abstol, + method = methode + ) + qn_vec = y[2, 2:(lx+1)] + + if (any(is.na(qn_vec)) && pars1[2] / pars1[1] < 1E-4 && missnumspec == 0) { + loglik = dd_loglik_high_lambda(pars1 = pars1, pars2 = pars2, brts = brts) + if (verbose) cat('High lambda approximation has been applied.\n') + return(loglik) + } + if (k < (S + 2 - soc)) { + qn_vec <- flavec(ddep, la, mu, K, r, lx, k1) * qn_vec # transition vector + } + cp <- check_probs(loglik, qn_vec, verbose) + loglik <- cp[[1]] + qn_vec <- cp[[2]] + if (loglik == -Inf | is.na(loglik) | is.nan(loglik)) break() + } + } else { # cond == 3 + qn_vec = rep(0,lx + 1) + qn_vec[1 + missnumspec] = 1 + for (k in (S + 2 - soc):2) { + k1 = k + (soc - 2) + y = dd_integrate( + initprobs = qn_vec, + tvec = -brts[k:(k-1)], + rhs_func = rhs_func_name, + pars = c(pars1,k1,ddep), + rtol = reltol, + atol = abstol, + method = methode + ) + qn_vec = y[2,2:(lx+2)] + if(k > soc) { + qn_vec = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * qn_vec # speciation event + } + cp <- check_probs(loglik,qn_vec[1:lx],verbose) + loglik <- cp[[1]] + qn_vec[1:lx] <- cp[[2]] + if (loglik == -Inf | is.na(loglik) | is.nan(loglik)) break() + } + } + if (qn_vec[1 + missnumspec] <= 0 | loglik == -Inf | is.na(loglik) | is.nan(loglik)) { + if(verbose) cat('Probabilities smaller than 0 or other numerical problems are encountered in final result.\n') + loglik = -Inf + } else { + loglik = loglik + (cond != 3 | soc == 1) * log(qn_vec[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) + + logliknorm = 0 + if(cond == 1 | cond == 2) { + probsn = rep(0, lx) + probsn[1] = 1 # change if other species at stem or crown age + k = soc + t1 = brts[1] + t2 = brts[S + 2 - soc] + y = dd_integrate( + initprobs = probsn, + tvec = c(t1,t2), + rhs_func = rhs_func_name, + pars = c(pars1,k,ddep), + rtol = reltol, + atol = abstol, + method = methode + ) + probsn = y[2, 2:(lx + 1)] + if(soc == 1) { + aux = 1:lx + } + if(soc == 2) { + aux = (2:(lx + 1)) * (3:(lx + 2)) / 6 + } + probsc = probsn / aux + cp <- check_probs(logliknorm,probsc,verbose) + logliknorm <- cp[[1]] + probsc <- cp[[2]] + if (cond == 1) { + logliknorm = logliknorm + log(sum(probsc)) + } + if (cond == 2) { + logliknorm = logliknorm + log(probsc[S + missnumspec - soc + 1]) + } + } + if (cond == 3) { + probsn = rep(0, lx + 1) + probsn[S + missnumspec + 1] = 1 + TT = max(1, 1 / abs(la - mu)) * 1E+10 * max(abs(brts)) # make this more efficient later + y = dd_integrate( + initprobs = probsn, + tvec = c(0, TT), + rhs_func = rhs_func_name, + pars = c(pars1, 0, ddep), + rtol = reltol, + atol = abstol, + method = methode + ) + logliknorm = log(y[2,lx + 2]) + if(soc == 2) { + probsn = rep(0,lx + 1) + probsn[1:lx] = qn_vec[1:lx] + probsn = c(flavec(ddep, la, mu, K, r, lx, 1), 1) * probsn # speciation event + y = dd_integrate( + initprobs = probsn, + tvec = c(max(abs(brts)), TT), + rhs_func = rhs_func_name, + pars = c(pars1,1,ddep), + rtol = reltol, + atol = abstol, + method = methode + ) + logliknorm = logliknorm - log(y[2,lx + 2]) + } + } + if(is.na(logliknorm) | is.nan(logliknorm) | logliknorm == Inf) { + if(verbose) cat('The normalization did not yield a number.\n') + loglik = -Inf + } else { + loglik = loglik - logliknorm + } + } + } + } + if(verbose) + { + s1 = sprintf('Parameters: %f %f %f',pars1[1],pars1[2],pars1[3]) + if(both_rates_vary(ddep)) {s1 = sprintf('%s %f',s1,pars1[4])} + s2 = sprintf(', Loglikelihood: %f',loglik) + cat(s1,s2,"\n",sep = "") + utils::flush.console() + } + loglik = as.numeric(loglik) + if(is.nan(loglik) | is.na(loglik)) { + loglik = -Inf + } + return(loglik) +} + +dd_loglik2 = function(pars1,pars2,brts,missnumspec) { if(length(pars2) == 4) { @@ -157,32 +403,42 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt pars2[6] = 2 } ddep = pars2[2] + #if (ddep > 5) { + # stop("This DD model is not implemented for the analytical method yet.") + #} cond = pars2[3] btorph = pars2[4] - verbose = pars2[5] + verbose <- pars2[5] soc = pars2[6] - if(cond == 3) { soc = 2 } + if(cond == 3) + { + soc = 2 + } la = pars1[1] mu = pars1[2] K = pars1[3] - if(ddep == 5) {r = pars1[4]} else {r = 0} + if(ddep == 5) + { + r = pars1[4] + } else + { + r = 0 + } if(ddep == 1 | ddep == 5) { lx = min(max(1 + missnumspec,1 + ceiling(la/(la - mu) * (r + 1) * K)),ceiling(pars2[1])) + } else if(ddep == 1.3) + { + lx = min(ceiling(K),ceiling(pars2[1])) } else { - if(ddep == 1.3) - { - lx = min(ceiling(K),ceiling(pars2[1])) - } else { - lx = round(pars2[1]) - } + lx = round(pars2[1]) } if((ddep == 1) & ((mu == 0 & missnumspec == 0 & floor(K) != ceiling(K) & la > 0.05) | K == Inf)) { loglik = bd_loglik(pars1[1:(2 + (K < Inf))],c(2*(mu == 0 & K < Inf),pars2[3:6]),brts,missnumspec) } else { - abstol = 1e-10 - reltol = 1e-8 + abstol = 1e-16 + reltol = 1e-10 brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) if(sum(brts == 0) == 0) { @@ -191,7 +447,6 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt S = length(brts) + (soc - 2) if(min(pars1) < 0) { - if(verbose) cat('The parameters are negative.\n') loglik = -Inf } else { if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | ddep == 4.1 | ddep == 4.2)) | (la <= mu)) @@ -199,9 +454,8 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt if(verbose) cat("These parameter values cannot satisfy lambda(N) = mu(N) for a positive and finite N.\n") loglik = -Inf } else { - if(((ddep == 1 | ddep == 5) & ceiling(la/(la - mu) * (r + 1) * K) < (S + missnumspec)) | ((ddep == 1.3) & (S + missnumspec > ceiling(K)))) + if(((ddep == 1 | ddep == 5) & ceiling(la/(la - mu) * (r + 1) * K) < (S + missnumspec)) | ((ddep == 1.3) & ((S + missnumspec) > ceiling(K)))) { - if(verbose) cat('The parameters are incompatible.\n') loglik = -Inf } else { loglik = (btorph == 0) * lgamma(S) @@ -212,19 +466,21 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt for(k in 2:(S + 2 - soc)) { k1 = k + (soc - 2) - y = dd_integrate(probs,brts[(k-1):k],rhs_func_name,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) - probs = y[2,2:(lx+1)] + #y = deSolve::ode(probs,brts[(k-1):k],rhs_func,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = "analytical") + #probs2 = y[2,2:(lx+1)] + probs = dd_loglik_M(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs) if(is.na(sum(probs)) && pars1[2]/pars1[1] < 1E-4 && missnumspec == 0) - { + { loglik = dd_loglik_high_lambda(pars1 = pars1,pars2 = pars2,brts = brts) if(verbose) cat('High lambda approximation has been applied.\n') return(loglik) } if(k < (S + 2 - soc)) { - probs = flavec(ddep,la,mu,K,r,lx,k1) * probs # speciation event + #probs = flavec(ddep,la,mu,K,r,lx,k1) * probs # speciation event + probs = lambdamu(0:(lx - 1) + k1,c(pars1[1:3],r),ddep)[[1]] * probs } - cp <- check_probs(loglik,probs,verbose); loglik <- cp[[1]]; probs <- cp[[2]]; + cp <- check_probs(loglik,probs,verbose); loglik <- cp[[1]]; probs<- cp[[2]]; } } else { probs = rep(0,lx + 1) @@ -232,20 +488,22 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt for(k in (S + 2 - soc):2) { k1 = k + (soc - 2) - y = dd_integrate(probs,-brts[k:(k-1)],rhs_func_name,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) - probs = y[2,2:(lx+2)] + #y = deSolve::ode(probs,-brts[k:(k-1)],dd_loglik_bw_rhs,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = "analytical") + #probs2 = y[2,2:(lx+2)] + probs = dd_loglik_M_bw(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs[1:lx]) + probs = c(probs,0) if(k > soc) { - probs = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * probs # speciation event - } + #probs = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * probs # speciation event + probs = c(lambdamu(0:(lx - 1) + k1 - 1,pars1,ddep)[[1]],1) * probs + } cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; } } - if(probs[1 + missnumspec] <= 0 | loglik == -Inf | is.na(loglik) | is.nan(loglik)) + if(probs[1 + (cond != 3) * missnumspec] <= 0 | loglik == -Inf) { - if(verbose) cat('Probabilities smaller than 0 or other numerical problems are encountered in final result.\n') loglik = -Inf - } else { + } else { loglik = loglik + (cond != 3 | soc == 1) * log(probs[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) logliknorm = 0 @@ -256,44 +514,51 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt k = soc t1 = brts[1] t2 = brts[S + 2 - soc] - y = dd_integrate(probsn,c(t1,t2),rhs_func_name,c(pars1,k,ddep),rtol = reltol,atol = abstol,method = methode); - probsn = y[2,2:(lx+1)] + #y = deSolve::ode(probsn,c(t1,t2),rhs_func,c(pars1,k,ddep),rtol = reltol,atol = abstol,method = "analytical") + #probsn = y[2,2:(lx+1)] + probsn = dd_loglik_M(pars1,lx,k,ddep,tt = abs(t2 - t1),probsn) if(soc == 1) { aux = 1:lx } if(soc == 2) { aux = (2:(lx+1)) * (3:(lx+2))/6 } probsc = probsn/aux - cp <- check_probs(logliknorm,probsc,verbose); logliknorm <- cp[[1]]; probsc <- cp[[2]]; - if(cond == 1) { logliknorm = logliknorm + log(sum(probsc)) } - if(cond == 2) { logliknorm = logliknorm + log(probsc[S + missnumspec - soc + 1])} + if(cond == 1) { logliknorm = log(sum(probsc)) } + if(cond == 2) { logliknorm = log(probsc[S + missnumspec - soc + 1])} } if(cond == 3) { + #probsn = rep(0,lx + 1) + #probsn[S + missnumspec + 1] = 1 #/ (S + missnumspec) + #TT = max(1,1/abs(la - mu)) * 100000000 * max(abs(brts)) # make this more efficient later + #y = deSolve::ode(probsn,c(0,TT),dd_loglik_bw_rhs,c(pars1,0,ddep),rtol = reltol,atol = abstol,method = "analytical") + #logliknorm = log(y[2,lx + 2]) probsn = rep(0,lx + 1) - probsn[S + missnumspec + 1] = 1 - - TT = 1e14 # max(1,1/abs(la - mu)) * 1E+10 * max(abs(brts)) # make this more efficient later - y = dd_integrate(probsn,c(0,TT),rhs_func_name,c(pars1,0,ddep),rtol = reltol,atol = abstol,method = methode) - logliknorm = log(y[2,lx + 2]) + probsn[2] = 1 + MM = dd_loglik_M_aux(pars1,lx + 1,k = 0,ddep) + MM = MM[-1,-1] + #probsn = SparseM::solve(-MM,probsn[2:(lx + 1)]) + MMinv = SparseM::solve(MM) + probsn = -MMinv %*% probsn[2:(lx + 1)] + logliknorm = log(probsn[S + missnumspec]) if(soc == 2) { - probsn = rep(0,lx + 1) - probsn[1:lx] = probs[1:lx] - probsn = c(flavec(ddep,la,mu,K,r,lx,1),1) * probsn # speciation event - y = dd_integrate(probsn,c(max(abs(brts)),TT),rhs_func_name,c(pars1,1,ddep),rtol = reltol,atol = abstol,method = methode) - logliknorm = logliknorm - log(y[2,lx + 2]) + #probsn = rep(0,lx + 1) + #probsn[1:lx] = probs[1:lx] + #probsn = c(flavec(ddep,la,mu,K,r,lx,1),1) * probsn # speciation event + #probsn = c(lambdamu(0:(lx - 1) + 1,pars1,ddep)[[1]],1) * probsn # speciation event + #y = deSolve::ode(probsn,c(max(abs(brts)),TT),dd_loglik_bw_rhs,c(pars1,1,ddep),rtol = reltol,atol = abstol,method = "analytical") + #logliknorm = logliknorm - log(y[2,lx + 2]) + probsn2 = rep(0,lx) + probsn2 = lambdamu(0:(lx - 1) + 1,pars1,ddep)[[1]] * probs[1:lx] + MM = dd_loglik_M_bw_aux(pars1,lx,k = 1,ddep) + MMinv = SparseM::solve(MM) + #probsn2 = SparseM::solve(-MM,probsn2[1:lx]) + probsn2 = -MMinv %*% probsn2[1:lx] + logliknorm = logliknorm - log(probsn2[1]) } } - if(is.na(logliknorm) | is.nan(logliknorm) | logliknorm == Inf) - { - if(verbose) cat('The normalization did not yield a number.\n') - loglik = -Inf - } else - { - loglik = loglik - logliknorm - } + loglik = loglik - logliknorm } } - } - } + }} if(verbose) { s1 = sprintf('Parameters: %f %f %f',pars1[1],pars1[2],pars1[3]) @@ -304,191 +569,13 @@ dd_loglik1 = function(pars1,pars2,brts,missnumspec,methode = 'odeint::runge_kutt } } loglik = as.numeric(loglik) - if(is.nan(loglik) | is.na(loglik)) + if(is.nan(loglik) | is.na(loglik) | loglik == Inf) { loglik = -Inf } return(loglik) } -dd_loglik2 = function(pars1,pars2,brts,missnumspec) -{ -if(length(pars2) == 4) -{ - pars2[5] = 0 - pars2[6] = 2 -} -ddep = pars2[2] -cond = pars2[3] -btorph = pars2[4] -verbose <- pars2[5] -soc = pars2[6] -if(cond == 3) -{ - soc = 2 -} -la = pars1[1] -mu = pars1[2] -K = pars1[3] -if(ddep == 5) -{ - r = pars1[4] -} else -{ - r = 0 -} -if(ddep == 1 | ddep == 5) -{ - lx = min(max(1 + missnumspec,1 + ceiling(la/(la - mu) * (r + 1) * K)),ceiling(pars2[1])) -} else if(ddep == 1.3) -{ - lx = min(ceiling(K),ceiling(pars2[1])) -} else { - lx = round(pars2[1]) -} -if((ddep == 1) & ((mu == 0 & missnumspec == 0 & floor(K) != ceiling(K) & la > 0.05) | K == Inf)) -{ - loglik = bd_loglik(pars1[1:(2 + (K < Inf))],c(2*(mu == 0 & K < Inf),pars2[3:6]),brts,missnumspec) -} else { -abstol = 1e-16 -reltol = 1e-10 -brts = -sort(abs(as.numeric(brts)),decreasing = TRUE) -if(sum(brts == 0) == 0) -{ - brts[length(brts) + 1] = 0 -} -S = length(brts) + (soc - 2) -if(min(pars1) < 0) -{ - loglik = -Inf -} else { -if((mu == 0 & (ddep == 2 | ddep == 2.1 | ddep == 2.2)) | (la == 0 & (ddep == 4 | ddep == 4.1 | ddep == 4.2)) | (la <= mu)) -{ - if(verbose) cat("These parameter values cannot satisfy lambda(N) = mu(N) for a positive and finite N.\n") - loglik = -Inf -} else { - if(((ddep == 1 | ddep == 5) & ceiling(la/(la - mu) * (r + 1) * K) < (S + missnumspec)) | ((ddep == 1.3) & ((S + missnumspec) > ceiling(K)))) - { - loglik = -Inf - } else { - loglik = (btorph == 0) * lgamma(S) - if(cond != 3) - { - probs = rep(0,lx) - probs[1] = 1 # change if other species at stem/crown age - for(k in 2:(S + 2 - soc)) - { - k1 = k + (soc - 2) - #y = deSolve::ode(probs,brts[(k-1):k],rhs_func,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) - #probs2 = y[2,2:(lx+1)] - probs = dd_loglik_M(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs) - if(is.na(sum(probs)) && pars1[2]/pars1[1] < 1E-4 && missnumspec == 0) - { - loglik = dd_loglik_high_lambda(pars1 = pars1,pars2 = pars2,brts = brts) - if(verbose) cat('High lambda approximation has been applied.\n') - return(loglik) - } - if(k < (S + 2 - soc)) - { - #probs = flavec(ddep,la,mu,K,r,lx,k1) * probs # speciation event - probs = lambdamu(0:(lx - 1) + k1,c(pars1[1:3],r),ddep)[[1]] * probs - } - cp <- check_probs(loglik,probs,verbose); loglik <- cp[[1]]; probs<- cp[[2]]; - } - } else { - probs = rep(0,lx + 1) - probs[1 + missnumspec] = 1 - for(k in (S + 2 - soc):2) - { - k1 = k + (soc - 2) - #y = deSolve::ode(probs,-brts[k:(k-1)],dd_loglik_bw_rhs,c(pars1,k1,ddep),rtol = reltol,atol = abstol,method = methode) - #probs2 = y[2,2:(lx+2)] - probs = dd_loglik_M_bw(pars1,lx,k1,ddep,tt = abs(brts[k] - brts[k-1]),probs[1:lx]) - probs = c(probs,0) - if(k > soc) - { - #probs = c(flavec(ddep,la,mu,K,r,lx,k1-1),1) * probs # speciation event - probs = c(lambdamu(0:(lx - 1) + k1 - 1,pars1,ddep)[[1]],1) * probs - } - cp <- check_probs(loglik,probs[1:lx],verbose); loglik <- cp[[1]]; probs[1:lx] <- cp[[2]]; - } - } - if(probs[1 + (cond != 3) * missnumspec] <= 0 | loglik == -Inf) - { - loglik = -Inf - } else { - loglik = loglik + (cond != 3 | soc == 1) * log(probs[1 + (cond != 3) * missnumspec]) - lgamma(S + missnumspec + 1) + lgamma(S + 1) + lgamma(missnumspec + 1) - - logliknorm = 0 - if(cond == 1 | cond == 2) - { - probsn = rep(0,lx) - probsn[1] = 1 # change if other species at stem or crown age - k = soc - t1 = brts[1] - t2 = brts[S + 2 - soc] - #y = deSolve::ode(probsn,c(t1,t2),rhs_func,c(pars1,k,ddep),rtol = reltol,atol = abstol,method = methode); - #probsn = y[2,2:(lx+1)] - probsn = dd_loglik_M(pars1,lx,k,ddep,tt = abs(t2 - t1),probsn) - if(soc == 1) { aux = 1:lx } - if(soc == 2) { aux = (2:(lx+1)) * (3:(lx+2))/6 } - probsc = probsn/aux - if(cond == 1) { logliknorm = log(sum(probsc)) } - if(cond == 2) { logliknorm = log(probsc[S + missnumspec - soc + 1])} - } - if(cond == 3) - { - #probsn = rep(0,lx + 1) - #probsn[S + missnumspec + 1] = 1 #/ (S + missnumspec) - #TT = max(1,1/abs(la - mu)) * 100000000 * max(abs(brts)) # make this more efficient later - #y = deSolve::ode(probsn,c(0,TT),dd_loglik_bw_rhs,c(pars1,0,ddep),rtol = reltol,atol = abstol,method = methode) - #logliknorm = log(y[2,lx + 2]) - probsn = rep(0,lx + 1) - probsn[2] = 1 - MM = dd_loglik_M_aux(pars1,lx + 1,k = 0,ddep) - MM = MM[-1,-1] - #probsn = SparseM::solve(-MM,probsn[2:(lx + 1)]) - MMinv = SparseM::solve(MM) - probsn = -MMinv %*% probsn[2:(lx + 1)] - logliknorm = log(probsn[S + missnumspec]) - if(soc == 2) - { - #probsn = rep(0,lx + 1) - #probsn[1:lx] = probs[1:lx] - #probsn = c(flavec(ddep,la,mu,K,r,lx,1),1) * probsn # speciation event - #probsn = c(lambdamu(0:(lx - 1) + 1,pars1,ddep)[[1]],1) * probsn # speciation event - #y = deSolve::ode(probsn,c(max(abs(brts)),TT),dd_loglik_bw_rhs,c(pars1,1,ddep),rtol = reltol,atol = abstol,method = methode) - #logliknorm = logliknorm - log(y[2,lx + 2]) - probsn2 = rep(0,lx) - probsn2 = lambdamu(0:(lx - 1) + 1,pars1,ddep)[[1]] * probs[1:lx] - MM = dd_loglik_M_bw_aux(pars1,lx,k = 1,ddep) - MMinv = SparseM::solve(MM) - #probsn2 = SparseM::solve(-MM,probsn2[1:lx]) - probsn2 = -MMinv %*% probsn2[1:lx] - logliknorm = logliknorm - log(probsn2[1]) - } - } - loglik = loglik - logliknorm - } - } -}} -if(verbose) -{ - s1 = sprintf('Parameters: %f %f %f',pars1[1],pars1[2],pars1[3]) - if(ddep == 5) {s1 = sprintf('%s %f',s1,pars1[4])} - s2 = sprintf(', Loglikelihood: %f',loglik) - cat(s1,s2,"\n",sep = "") - utils::flush.console() -} -} -loglik = as.numeric(loglik) -if(is.nan(loglik) | is.na(loglik) | loglik == Inf) -{ - loglik = -Inf -} -return(loglik) -} - dd_int <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) { if(method == 'analytical') @@ -508,8 +595,7 @@ dd_int <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) dd_integrate <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) { - if(method == 'analytical') - { + if(method == 'analytical') { probs <- dd_loglik_M(pars = pars[1:3], lx = length(initprobs), k = pars[4], @@ -517,11 +603,9 @@ dd_integrate <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) tt = abs(tvec[2] - tvec[1]), initprobs) y <- cbind(c(NA,NA),rbind(rep(NA,length(probs)),probs)) - } else - { - #rhs_func_name <- 'no_name' - #if(is.character(rhs_func)) - #{ + } else { + rhs_func_name <- 'no_name' + if (is.character(rhs_func)) { rhs_func_name <- rhs_func #if(rhs_func_name != 'dd_loglik_rhs' & rhs_func_name != 'dd_loglik_bw_rhs') #{ @@ -531,8 +615,7 @@ dd_integrate <- function(initprobs,tvec,rhs_func,pars,rtol,atol,method) if(rhs_func_name == 'dd_loglik_rhs' || rhs_func_name == 'dd_loglik_bw_rhs') { parsvec = c(dd_loglik_rhs_precomp(pars,initprobs),pars[length(pars) - 1]) - } else - { + } else { parsvec = pars } if (startsWith(method, "odeint::")) { diff --git a/R/dd_loglik_M.R b/R/dd_loglik_M.R index d295ce1..2184d92 100644 --- a/R/dd_loglik_M.R +++ b/R/dd_loglik_M.R @@ -1,41 +1,116 @@ lambdamu = function(n,pars,ddep) { lnn = length(n) - zeros = rep(0,lnn) - ones = rep(1,lnn) la = pars[1] mu = pars[2] K = pars[3] r = pars[4] - lavec = la * ones - muvec = mu * ones + phi <- ifelse(r == Inf, 1, r / (1 + r)) # else r/(1+r) is NaN + eq_rate <- phi * la + (1 - phi) * mu n0 = (ddep == 2 | ddep == 4) - if(ddep == 1) - { - lavec = pmax(zeros,la - (la - mu) * n / K) - } else if(ddep == 1.3) - { - lavec = pmax(zeros,la * (1 - n/K)) - } else if(ddep == 2 | ddep == 2.1 | ddep == 2.2) - { - y = -(log(la/mu)/log(K+n0))^(ddep != 2.2) - lavec = pmax(zeros,la * (n + n0)^y) - } else if(ddep == 2.3) - { - y = -K - lavec = pmax(zeros,la * (n + n0)^y) - } else if(ddep == 3) - { - lavec = la * ones - muvec = mu + (la - mu) * n/K - } else if(ddep == 4 | ddep == 4.1 | ddep == 4.2) - { - y = (log(la/mu)/log(K+n0))^(ddep != 4.2) - muvec = mu * (n + n0)^y - } else if(ddep == 5) - { - lavec = pmax(zeros,la - 1/(r + 1)*(la - mu)/K * n) - muvec = muvec = mu + r/(r + 1)*(la - mu)/K * n + if(ddep == 1) { + # linear DD on speciation (K = equilibrium diversity) + # constant-rate extinction + lavec = pmax(0, la - (la - mu) * n / K) + muvec = rep(mu, lnn) + } else if(ddep == 1.3) { + # linear DD on speciation (K = carrying capacity) + # constant-rate extinction + lavec = pmax(0, la * (1 - n / K)) + muvec = rep(mu, lnn) + } else if(ddep == 2 | ddep == 2.1 | ddep == 2.2 | ddep == 2.4) { + # "exponential" DD on speciation (power function) + # constant-rate extinction + frac <- ifelse(ddep == 2.3, 0.1, la / mu) + y = -(log( frac ) / log(K + n0)) ^ (ddep != 2.2) + lavec = pmax(0, la * (n + n0) ^ y) + muvec = rep(mu, lnn) + } else if(ddep == 2.3) { + # "exponential" DD on speciation (power function) + # constant-rate extinction + y = -K + lavec = pmax(0, la * (n + n0) ^ y) + muvec = rep(mu, lnn) + } else if(ddep == 3) { + # constant-rate speciation + # linear DD on extinction + lavec = rep(la, lnn) + muvec = mu + (la - mu) * n/K + } else if (ddep == 4 | ddep == 4.1 | ddep == 4.2) { + # constant-rate speciation + # "exponential" DD on extinction (power function) + y = (log(la / mu) / log(K + n0)) ^ (ddep != 4.2) + lavec = rep(la, lnn) + muvec = mu * (n + n0) ^ y + } else if (ddep == 5) { + # linear DD on speciation + # linear DD on extinction + lavec = pmax(0, la - (la - eq_rate) * n / K) + muvec = mu + (eq_rate - mu) * n / K + } else if (ddep == 6) { + # linear DD on speciation + # "exponential" DD on extinction (power function) + y = log(eq_rate / mu) / log(K) + lavec = pmax(0, la - (la - eq_rate) * n / K) + muvec = mu * n ^ y + } else if (ddep == 7) { + # "exponential" DD on speciation (power function) + # "exponential" DD on extinction (power function) + y1 = -log(la / eq_rate) / log(K) + y2 = log(eq_rate / mu) / log(K) + lavec = pmax(0, la * n ^ y1) + muvec = mu * n ^ y2 + } else if (ddep == 8) { + # "exponential" DD on speciation (power function) + # linear DD on extinction + y = -log(la / eq_rate) / log(K) + lavec = pmax(0, la * n ^ y) + muvec = mu + (eq_rate - mu) / K * n + } else if (ddep == 9) { + # exponential DD on speciation (exponential function) + # constant-rate extinction + y = log(la / mu) / K + lavec = pmax(0, la * exp(-n * y)) + muvec = rep(mu, lnn) + } else if (ddep == 10) { + # constant-rate speciation + # exponential DD on extinction (exponential function) + y = log(la / mu) / K + lavec = rep(la, lnn) + muvec = mu * exp(n * y) + } else if (ddep == 11) { + # linear DD on speciation + # exponential DD on extinction (exponential function) + y = log(eq_rate / mu) / K + lavec = pmax(0, la - (la - eq_rate) * n / K) + muvec = mu * exp(n * y) + } else if (ddep == 12) { + # exponential DD on speciation (exponential function) + # exponential DD on extinction (exponential function) + y1 = log(la / eq_rate) / K + y2 = log(eq_rate / mu) / K + lavec = pmax(0, la * exp(-n * y1)) + muvec = mu * exp(n * y2) + } else if (ddep == 13) { + # exponential DD on speciation (exponential function) + # linear DD on extinction + y = log(la / eq_rate) / K + lavec = pmax(0, la * exp(-n * y)) + muvec = mu + (eq_rate - mu) * n / K + } else if (ddep == 14) { + # exponential DD on speciation (exponential function) + # exponential DD on extinction (power function) + y1 = log(la / eq_rate) / K + y2 = log(eq_rate / mu) / log(K) + lavec = pmax(0, la * exp(-n * y1)) + muvec = mu * n ^ y2 + } else if (ddep == 15) { + # exponential DD on speciation (power function) + # exponential DD on extinction (exponential function) + y1 = -log(la / eq_rate) / log(K) + y2 = log(eq_rate / mu) / K + lavec = pmax(0, la * n ^ y1) + muvec = mu * exp(n * y2) } return(list(lavec,muvec)) } @@ -58,21 +133,21 @@ changepars = function(pars) { if(length(pars) <= 3) { - sel = 1:2 + sel = 1:2 } else { - sel = c(1:2,4:min(length(pars),5)) + sel = c(1:2,4:min(length(pars),5)) } if(sum(pars[sel] == Inf) > 0) { - pars[which(pars[sel] == Inf)] = 1E+10 + pars[which(pars[sel] == Inf)] = 1E+10 } if(sum(pars[sel] == 0) > 0) { - pars[which(pars[sel] == 0)] = 1E-14 + pars[which(pars[sel] == 0)] = 1E-14 } if(sum(pars[sel] == 1) > 0) { - pars[which(pars[sel] == 1)] = pars[which(pars[sel] == 1)] - 1E-14 + pars[which(pars[sel] == 1)] = pars[which(pars[sel] == 1)] - 1E-14 } return(pars) } @@ -124,9 +199,10 @@ lambdamu2 = function(nxt,pars,ddep) { lavec = pmax(matrix(0,lnn,lnn),laM * (1 - nxt/KM)) muvec = muM * matrix(1,lnn,lnn) - } else if(ddep == 2 | ddep == 2.1 | ddep == 2.2) + } else if(ddep == 2 | ddep == 2.1 | ddep == 2.2 | ddep == 2.4) { - x = -(log(laM/muM)/log(KM+n0))^(ddep != 2.2) + fracM <- ifelse(ddep == 2.4, 10, laM/muM) + x = -(log( fracM )/log(KM + n0))^(ddep != 2.2) lavec = pmax(matrix(0,lnn,lnn),laM * (nxt + n0)^x) muvec = muM * matrix(1,lnn,lnn) } else if(ddep == 2.3) diff --git a/R/dd_loglik_choosepar.R b/R/dd_loglik_choosepar.R index fdad519..a3d32b5 100644 --- a/R/dd_loglik_choosepar.R +++ b/R/dd_loglik_choosepar.R @@ -1,9 +1,9 @@ dd_loglik_choosepar = function(trparsopt,trparsfix,idparsopt,idparsfix,pars2,brts,missnumspec,methode) { - trpars1 = rep(0,3) - if(pars2[2] == 5) - { + if (both_rates_vary(pars2[2])) { trpars1 = rep(0,4) + } else { + trpars1 = rep(0,3) } trpars1[idparsopt] = trparsopt if(length(idparsfix) != 0) diff --git a/R/dd_loglik_rhs.R b/R/dd_loglik_rhs.R index b6b080b..fb2624f 100644 --- a/R/dd_loglik_rhs.R +++ b/R/dd_loglik_rhs.R @@ -4,88 +4,114 @@ dd_loglik_rhs_precomp = function(pars,x) la = pars[1] mu = pars[2] K = pars[3] - if(length(pars) < 6) - { + if (length(pars) < 6) { kk = pars[4] ddep = pars[5] } else { r = pars[4] kk = pars[5] ddep = pars[6] + phi = ifelse(r == Inf, 1, r / (1 + r)) # else r/(1+r) can be NaN + eq_rate = phi * la + (1 - phi) * mu } n0 = (ddep == 2 | ddep == 4) - nn = -1:(lx+2*kk) - lnn = length(nn) - nn = pmax(rep(0,lnn),nn) + nn <- c(0, 0:(lx + 2 * kk)) + lnn <- length(nn) - if(ddep == 1) - { - lavec = pmax(rep(0,lnn),la - (la - mu)/K * nn) - muvec = mu * rep(1,lnn) - } else { - if(ddep == 1.3) - { - lavec = pmax(rep(0,lnn),la * (1 - nn/K)) - muvec = mu * rep(1,lnn) - } else { - if(ddep == 1.4) - { - lavec = pmax(rep(0,lnn),la * nn/(nn + K)) - } else { - if(ddep == 1.5) - { - lavec = pmax(rep(0,lnn),la * nn/K * (1 - nn/K)) - muvec = mu * rep(1,lnn) - } else { - if(ddep == 2 | ddep == 2.1 | ddep == 2.2) - { - y = -(log(la/mu)/log(K+n0))^(ddep != 2.2) - lavec = pmax(rep(0,lnn),la * (nn + n0)^y) - muvec = mu * rep(1,lnn) - } else { - if(ddep == 2.3) - { - y = -K - lavec = pmax(rep(0,lnn),la * (nn + n0)^y) - muvec = mu * rep(1,lnn) - } else { - if(ddep == 3) - { - lavec = la * rep(1,lnn) - muvec = mu + (la - mu) * nn/K - } else { - if(ddep == 4 | ddep == 4.1 | ddep == 4.2) - { - lavec = la * rep(1,lnn) - y = (log(la/mu)/log(K+n0))^(ddep != 4.2) - muvec = mu * (nn + n0)^y - } else { - if(ddep == 5) - { - lavec = pmax(rep(0,lnn),la - 1/(r + 1)*(la - mu)/K * nn) - muvec = muvec = mu + r/(r + 1)*(la - mu)/K * nn - } - } - } - } - } - } - } - } + if (ddep == 1) { + lavec = pmax(0, la - (la - mu) / K * nn) + muvec = rep(mu, lnn) + } else if (ddep == 1.3) { + lavec = pmax(0, la * (1 - nn / K)) + muvec = rep(mu, lnn) + } else if (ddep == 1.4) { + lavec = pmax(0, la * nn / (nn + K)) + } else if(ddep == 1.5) { + lavec = pmax(0, la * nn / K * (1 - nn / K)) + muvec = rep(mu, lnn) + } else if (ddep == 2 | ddep == 2.1 | ddep == 2.2 | ddep == 2.4) { + frac <- ifelse(ddep == 2.4, 10, la / mu) + y = -(log( frac ) / log(K + n0)) ^ (ddep != 2.2) + lavec = pmax(0, la * (nn + n0) ^ y) + muvec = rep(mu, lnn) + } else if(ddep == 2.3) { + y = -K + lavec = pmax(0, la * (nn + n0) ^ y) + muvec = rep(mu, lnn) + } else if (ddep == 3) { + lavec = rep(la, lnn) + muvec = mu + (la - mu) * nn / K + } else if (ddep == 4 | ddep == 4.1 | ddep == 4.2) { + lavec = rep(la, lnn) + y = (log(la / mu) / log(K + n0)) ^ (ddep != 4.2) + muvec = mu * (nn + n0) ^ y + } else if (ddep == 5) { + lavec = pmax(0, la - (la - eq_rate) * nn / K) + muvec = mu + (eq_rate - mu) / K * nn + } else if (ddep == 6) { + y = log(eq_rate / mu) / log(K) + lavec = pmax(0, la - (la - eq_rate) * nn / K) + muvec = mu * nn ^ y + } else if (ddep == 7) { + y1 = -log(la / eq_rate) / log(K) + y2 = log(eq_rate / mu) / log(K) + lavec = pmax(0, la * nn ^ y1) + muvec = mu * nn ^ y2 + } else if (ddep == 8) { + y = -log(la / eq_rate) / log(K) + lavec = pmax(0, la * nn ^ y) + muvec = mu + (eq_rate - mu) / K * nn + } else if (ddep == 9) { + y = log(la / mu) / K + lavec = pmax(0, la * exp(-nn * y)) + muvec = rep(mu, lnn) + } else if (ddep == 10) { + y = log(la / mu) / K + lavec = rep(la, lnn) + muvec = mu * exp(nn * y) + } else if (ddep == 11) { + y = log(eq_rate / mu) / K + lavec = pmax(0, la - (la - eq_rate) * nn / K ) + muvec = mu * exp(nn * y) + } else if (ddep == 12) { + y1 = log(la / eq_rate) / K + y2 = log(eq_rate / mu) / K + lavec = pmax(0, la * exp(-nn * y1)) + muvec = mu * exp(nn * y2) + } else if (ddep == 13) { + y = log(la / eq_rate) / K + lavec = pmax(0, la * exp(-nn * y)) + muvec = mu + (eq_rate - mu) / K * nn + } else if (ddep == 14) { + y1 = log(la / eq_rate) / K + y2 = log(eq_rate / mu) / log(K) + lavec = pmax(0, la * exp(-nn * y1)) + muvec = mu * nn ^ y2 + } else if (ddep == 15) { + y1 = -log(la / eq_rate) / log(K) + y2 = log(eq_rate / mu) / K + lavec = pmax(0, la * nn ^ y1) + muvec = mu * exp(nn * y2) } - return(c(lavec,muvec,nn)) + return(c(lavec, muvec, nn)) } dd_loglik_rhs = function(t,x,parsvec) { + # Unpack parameter vector lv = (length(parsvec) - 1)/3 lavec = parsvec[1:lv] muvec = parsvec[(lv + 1):(2 * lv)] nn = parsvec[(2 * lv + 1):(3 * lv)] - kk = parsvec[length(parsvec)] + k = parsvec[length(parsvec)] # nb species in phylo at time t + lx = length(x) - xx = c(0,x,0) - dx = lavec[(2:(lx+1))+kk-1] * nn[(2:(lx+1))+2*kk-1] * xx[(2:(lx+1))-1] + muvec[(2:(lx+1))+kk+1] * nn[(2:(lx+1))+1] * xx[(2:(lx+1))+1] - (lavec[(2:(lx+1))+kk] + muvec[(2:(lx+1))+kk]) * nn[(2:(lx+1))+kk] * xx[2:(lx+1)] + qn_vec = c(0, x, 0) + nvec <- 2:(lx + 1) + # DDD master system + dx <- lavec[nvec + k - 1] * nn[nvec + 2 * k - 1] * qn_vec[nvec - 1] + + muvec[nvec + k + 1] * nn[nvec + 1] * qn_vec[nvec + 1] - + (lavec[nvec + k] + muvec[nvec + k]) * nn[nvec + k] * qn_vec[nvec] return(list(dx)) } \ No newline at end of file diff --git a/R/dd_sim.R b/R/dd_sim.R index b44e0a8..8c2099a 100644 --- a/R/dd_sim.R +++ b/R/dd_sim.R @@ -1,57 +1,89 @@ -dd_lamuN = function(ddmodel,pars,N) -{ +dd_lamuN = function(ddmodel,pars,N) { + testthat::expect_length(N, 1) la = pars[1] mu = pars[2] K = pars[3] n0 = (ddmodel == 2 | ddmodel == 4) - if(length(pars) == 4) - { + if (length(pars) == 4) { r = pars[4] + phi <- ifelse(r == Inf, 1, r / (1 + r)) + eq_rate <- phi * la + (1 - phi) * mu } - if(ddmodel == 1) - { + if (ddmodel == 1) { # linear dependence in speciation rate laN = max(0,la - (la - mu) * N/K) muN = mu - } - if(ddmodel == 1.3) - { + } else if (ddmodel == 1.3) { # linear dependence in speciation rate laN = max(0,la * (1 - N/K)) muN = mu - } - if(ddmodel == 2 | ddmodel == 2.1 | ddmodel == 2.2) - { + } else if (ddmodel == 2 | ddmodel == 2.1 | ddmodel == 2.2) { # exponential dependence in speciation rate al = (log(la/mu)/log(K+n0))^(ddmodel != 2.2) laN = la * (N + n0)^(-al) muN = mu - } - if(ddmodel == 2.3) - { + } else if(ddmodel == 2.3) { # exponential dependence in speciation rate al = K laN = la * (N + n0)^(-al) muN = mu - } - if(ddmodel == 3) - { + } else if(ddmodel == 3) { # linear dependence in extinction rate laN = la muN = mu + (la - mu) * N/K - } - if(ddmodel == 4 | ddmodel == 4.1 | ddmodel == 4.2) - { + } else if(ddmodel == 4 | ddmodel == 4.1 | ddmodel == 4.2) { # exponential dependence in extinction rate al = (log(la/mu)/log(K+n0))^(ddmodel != 4.2) laN = la muN = mu * (N + n0)^al - } - if(ddmodel == 5) - { + } else if (ddmodel == 5) { # linear dependence in speciation rate and extinction rate - laN = max(0,la - 1/(r+1)*(la-mu) * N/K) - muN = mu + r/(r+1)*(la-mu)/K * N + laN = max(0, la - (la - eq_rate) * N / K) + muN = mu + (eq_rate - mu) * N / K + } else if (ddmodel == 6) { + y = log(eq_rate / mu) / log(K) + laN = max(0, la - (1 - phi) * (la - mu) * N / K) + muN = mu * N ^ y + } else if (ddmodel == 7) { + y1 = -log(la / eq_rate) / log(K) + y2 = log(eq_rate / mu) / log(K) + laN = max(0, la * N ^ y1) + muN = mu * N ^ y2 + } else if (ddmodel == 8) { + y = -log(la / eq_rate) / log(K) + laN = max(0, la * N ^ y) + muN = mu + (eq_rate - mu) / K * N + } else if (ddmodel == 9) { + y = log(la / mu) / K + laN = max(0, la * exp(-N * y)) + muN = mu + } else if (ddmodel == 10) { + y = log(la / mu) / K + laN = la + muN = mu * exp(N * y) + } else if (ddmodel == 11) { + y = log(eq_rate / mu) / K + laN = max(0, la - (la - eq_rate) / K * N) + muN = mu * exp(N * y) + } else if (ddmodel == 12) { + y1 = log(la / eq_rate) / K + y2 = log(eq_rate / mu) / K + laN = max(0, la * exp(-N * y1)) + muN = mu * exp(N * y2) + } else if (ddmodel == 13) { + y = log(la / eq_rate) / K + laN = max(0, la * exp(-N * y)) + muN = mu + (eq_rate - mu) / K * N + } else if (ddmodel == 14) { + y1 = log(la / eq_rate) / K + y2 = log(eq_rate / mu) / log(K) + laN = max(0, la * exp(-N * y1)) + muN = mu * N ^ y2 + } else if (ddmodel == 15) { + y1 = -log(la / eq_rate) / log(K) + y2 = log(eq_rate / mu) / K + laN = max(0, la * N ^ y1) + muN = mu * exp(N * y2) } return(c(laN,muN)) } diff --git a/R/dd_utils.R b/R/dd_utils.R index 5815c58..c1a9e97 100644 --- a/R/dd_utils.R +++ b/R/dd_utils.R @@ -5,7 +5,6 @@ #' #' Converting a set of branching times to a phylogeny #' -#' #' @param times Set of branching times #' @param root When root is FALSE, the largest branching time will be assumed to be #' the crown age. When root is TRUE, it will be the stem age. @@ -25,11 +24,11 @@ brts2phylo <- function(times,root=FALSE,tip.label=NULL) n <- n-1 } nbr <- 2*n - 2 - + # create the data types for edges and edge-lengths edge <- matrix(NA, nbr, 2) edge.length <- numeric(nbr) - + h <- numeric(2*n - 1) # initialized with 0's pool <- 1:n # VERY VERY IMPORTANT: the root MUST have index n+1 !!! @@ -55,24 +54,24 @@ brts2phylo <- function(times,root=FALSE,tip.label=NULL) nextnode <- nextnode - 1L } } - + phy <- list(edge = edge, edge.length = edge.length) if (is.null(tip.label)) tip.label <- paste("t", 1:n, sep = "") phy$tip.label <- sample(tip.label) phy$Nnode <- n - 1L - + if ( root ) { phy$root.edge <- times[n] - times[n-1] phy$root <- times[n] - times[n-1] } - + class(phy) <- "phylo" - + phy <- ape::reorder.phylo(phy) ## to avoid crossings when converting with as.hclust: phy$edge[phy$edge[, 2] <= n, 2] <- 1:n - + return(phy) } @@ -95,60 +94,61 @@ brts2phylo <- function(times,root=FALSE,tip.label=NULL) #' @export conv conv = function(x,y) { - lx = length(x) - ly = length(y) - lxy = length(x) + length(y) - x = c(x,rep(0,lxy - lx)) - y = c(y,rep(0,lxy - ly)) - cvxy = rep(0,lxy) - for(i in 2:lxy) - { - cvxy[i] = crossprod(x[(i-1):1],y[1:(i-1)]) - } - return(cvxy[2:lxy]) + lx = length(x) + ly = length(y) + lxy = length(x) + length(y) + x = c(x,rep(0,lxy - lx)) + y = c(y,rep(0,lxy - ly)) + cvxy = rep(0,lxy) + for(i in 2:lxy) + { + cvxy[i] = crossprod(x[(i-1):1],y[1:(i-1)]) + } + return(cvxy[2:lxy]) } flavec <- function(ddep,la,mu,K,r,lx,kk) { nn <- (0:(lx - 1)) + kk - lambdamu_nk <- lambdamu(nn,c(la,mu,K,r),ddep) - return(lambdamu_nk[[1]]) + lambda_n_vec <- lambdamu(nn, c(la, mu, K, r), ddep)[[1]] + return(lambda_n_vec) } flavec2 = function(ddep,la,mu,K,r,lx,kk,n0) { - nn = (0:(lx - 1)) + kk - if(ddep == 1 | ddep == 5) - { - lavec = pmax(rep(0,lx),la - 1/(r + 1) * (la - mu)/K * nn) - } - if(ddep == 1.3) - { - lavec = pmax(rep(0,lx),la * (1 - nn/K)) - } - if(ddep == 1.4) - { - lavec = pmax(rep(0,lx),la * nn/(nn + K)) - } - if(ddep == 1.5) - { - lavec = pmax(rep(0,lx),la * nn/K * (1 - nn/K)) - } - if(ddep == 2 | ddep == 2.1 | ddep == 2.2) - { - x = -(log(la/mu)/log(K + n0))^(ddep != 2.2) - lavec = pmax(rep(0,lx),la * (nn + n0)^x) - } - if(ddep == 2.3) - { - x = -K - lavec = pmax(rep(0,lx),la * (nn + n0)^x) - } - if(ddep == 3 | ddep == 4 | ddep == 4.1 | ddep == 4.2) - { - lavec = la * rep(1,lx) - } - return(lavec) + nn = (0:(lx - 1)) + kk + if(ddep == 1 | ddep == 5) + { + lavec = pmax(rep(0,lx),la - 1/(r + 1) * (la - mu)/K * nn) + } + if(ddep == 1.3) + { + lavec = pmax(rep(0,lx),la * (1 - nn/K)) + } + if(ddep == 1.4) + { + lavec = pmax(rep(0,lx),la * nn/(nn + K)) + } + if(ddep == 1.5) + { + lavec = pmax(rep(0,lx),la * nn/K * (1 - nn/K)) + } + if(ddep == 2 | ddep == 2.1 | ddep == 2.2 | ddep == 2.4) + { + frac <- ifelse(ddep == 2.4, 10, la / mu) + x = -(log( frac )/log(K + n0))^(ddep != 2.2) + lavec = pmax(rep(0,lx),la * (nn + n0)^x) + } + if(ddep == 2.3) + { + x = -K + lavec = pmax(rep(0,lx),la * (nn + n0)^x) + } + if(ddep == 3 | ddep == 4 | ddep == 4.1 | ddep == 4.2) + { + lavec = la * rep(1,lx) + } + return(lavec) } @@ -183,51 +183,51 @@ flavec2 = function(ddep,la,mu,K,r,lx,kk,n0) #' #' @export L2phylo L2phylo = function(L,dropextinct = T) -# makes a phylogeny out of a matrix with branching times, parent and daughter species, and extinction times + # makes a phylogeny out of a matrix with branching times, parent and daughter species, and extinction times { - L = L[order(abs(L[,3])),1:4] - age = L[1,1] - L[,1] = age - L[,1] - L[1,1] = -1 - notmin1 = which(L[,4] != -1) - L[notmin1,4] = age - L[notmin1,4] - if(dropextinct == T) - { - sall = which(L[,4] == -1) - tend = age - } else { - sall = which(L[,4] >= -1) - tend = (L[,4] == -1) * age + (L[,4] > -1) * L[,4] - } - L = L[,-4] - linlist = cbind(data.frame(L[sall,]),paste("t",abs(L[sall,3]),sep = ""),tend) - linlist[,4] = as.character(linlist[,4]) - names(linlist) = 1:5 - done = 0 - while(done == 0) - { - j = which.max(linlist[,1]) - daughter = linlist[j,3] - parent = linlist[j,2] - parentj = which(parent == linlist[,3]) - parentinlist = length(parentj) - if(parentinlist == 1) - { - spec1 = paste(linlist[parentj,4],":",linlist[parentj,5] - linlist[j,1],sep = "") - spec2 = paste(linlist[j,4],":",linlist[j,5] - linlist[j,1],sep = "") - linlist[parentj,4] = paste("(",spec1,",",spec2,")",sep = "") - linlist[parentj,5] = linlist[j,1] - linlist = linlist[-j,] - } else { - #linlist[j,1:3] = L[abs(as.numeric(parent)),1:3] - linlist[j,1:3] = L[which(L[,3] == parent),1:3] - } - if(nrow(linlist) == 1) { done = 1 } - } - linlist[4] = paste(linlist[4],":",linlist[5],";",sep = "") - phy = ape::read.tree(text = linlist[1,4]) - tree = ape::as.phylo(phy) - return(tree) + L = L[order(abs(L[,3])),1:4] + age = L[1,1] + L[,1] = age - L[,1] + L[1,1] = -1 + notmin1 = which(L[,4] != -1) + L[notmin1,4] = age - L[notmin1,4] + if(dropextinct == T) + { + sall = which(L[,4] == -1) + tend = age + } else { + sall = which(L[,4] >= -1) + tend = (L[,4] == -1) * age + (L[,4] > -1) * L[,4] + } + L = L[,-4] + linlist = cbind(data.frame(L[sall,]),paste("t",abs(L[sall,3]),sep = ""),tend) + linlist[,4] = as.character(linlist[,4]) + names(linlist) = 1:5 + done = 0 + while(done == 0) + { + j = which.max(linlist[,1]) + daughter = linlist[j,3] + parent = linlist[j,2] + parentj = which(parent == linlist[,3]) + parentinlist = length(parentj) + if(parentinlist == 1) + { + spec1 = paste(linlist[parentj,4],":",linlist[parentj,5] - linlist[j,1],sep = "") + spec2 = paste(linlist[j,4],":",linlist[j,5] - linlist[j,1],sep = "") + linlist[parentj,4] = paste("(",spec1,",",spec2,")",sep = "") + linlist[parentj,5] = linlist[j,1] + linlist = linlist[-j,] + } else { + #linlist[j,1:3] = L[abs(as.numeric(parent)),1:3] + linlist[j,1:3] = L[which(L[,3] == parent),1:3] + } + if(nrow(linlist) == 1) { done = 1 } + } + linlist[4] = paste(linlist[4],":",linlist[5],";",sep = "") + phy = ape::read.tree(text = linlist[1,4]) + tree = ape::as.phylo(phy) + return(tree) } @@ -385,52 +385,52 @@ phylo2L = function(phy) #' #' @export L2brts L2brts = function(L,dropextinct = T) -# makes a phylogeny out of a matrix with branching times, parent and daughter species, and extinction times + # makes a phylogeny out of a matrix with branching times, parent and daughter species, and extinction times { - brts = NULL - L = L[order(abs(L[,3])),1:4] - age = L[1,1] - L[,1] = age - L[,1] - L[1,1] = -1 - notmin1 = which(L[,4] != -1) - L[notmin1,4] = age - L[notmin1,4] - if(dropextinct == T) - { - sall = which(L[,4] == -1) - tend = age - } else { - sall = which(L[,4] >= -1) - tend = (L[,4] == -1) * age + (L[,4] > -1) * L[,4] - } - L = L[,-4] - linlist = cbind(data.frame(L[sall,]),paste("t",abs(L[sall,3]),sep = ""),tend) - linlist[,4] = as.character(linlist[,4]) - names(linlist) = 1:5 - done = 0 - while(done == 0) - { - j = which.max(linlist[,1]) - daughter = linlist[j,3] - parent = linlist[j,2] - parentj = which(parent == linlist[,3]) - parentinlist = length(parentj) - if(parentinlist == 1) - { - spec1 = paste(linlist[parentj,4],":",linlist[parentj,5] - linlist[j,1],sep = "") - spec2 = paste(linlist[j,4],":",linlist[j,5] - linlist[j,1],sep = "") - linlist[parentj,4] = paste("(",spec1,",",spec2,")",sep = "") - linlist[parentj,5] = linlist[j,1] - brts = c(brts,linlist[j,1]) - linlist = linlist[-j,] - } else { - #linlist[j,1:3] = L[abs(as.numeric(parent)),1:3] - linlist[j,1:3] = L[which(L[,3] == parent),1:3] - } - if(nrow(linlist) == 1) { done = 1 } - } - #linlist[4] = paste(linlist[4],":",linlist[5],";",sep = "") - brts = rev(sort(age - brts)) - return(brts) + brts = NULL + L = L[order(abs(L[,3])),1:4] + age = L[1,1] + L[,1] = age - L[,1] + L[1,1] = -1 + notmin1 = which(L[,4] != -1) + L[notmin1,4] = age - L[notmin1,4] + if(dropextinct == T) + { + sall = which(L[,4] == -1) + tend = age + } else { + sall = which(L[,4] >= -1) + tend = (L[,4] == -1) * age + (L[,4] > -1) * L[,4] + } + L = L[,-4] + linlist = cbind(data.frame(L[sall,]),paste("t",abs(L[sall,3]),sep = ""),tend) + linlist[,4] = as.character(linlist[,4]) + names(linlist) = 1:5 + done = 0 + while(done == 0) + { + j = which.max(linlist[,1]) + daughter = linlist[j,3] + parent = linlist[j,2] + parentj = which(parent == linlist[,3]) + parentinlist = length(parentj) + if(parentinlist == 1) + { + spec1 = paste(linlist[parentj,4],":",linlist[parentj,5] - linlist[j,1],sep = "") + spec2 = paste(linlist[j,4],":",linlist[j,5] - linlist[j,1],sep = "") + linlist[parentj,4] = paste("(",spec1,",",spec2,")",sep = "") + linlist[parentj,5] = linlist[j,1] + brts = c(brts,linlist[j,1]) + linlist = linlist[-j,] + } else { + #linlist[j,1:3] = L[abs(as.numeric(parent)),1:3] + linlist[j,1:3] = L[which(L[,3] == parent),1:3] + } + if(nrow(linlist) == 1) { done = 1 } + } + #linlist[4] = paste(linlist[4],":",linlist[5],";",sep = "") + brts = rev(sort(age - brts)) + return(brts) } @@ -459,9 +459,9 @@ L2brts = function(L,dropextinct = T) #' @export roundn roundn = function(x, digits = 0) { - fac = 10^digits - n = trunc(fac * x + 0.5)/fac - return(n) + fac = 10^digits + n = trunc(fac * x + 0.5)/fac + return(n) } @@ -490,21 +490,21 @@ roundn = function(x, digits = 0) #' @export sample2 sample2 = function(x,size,replace = FALSE,prob = NULL) { - if(length(x) == 1) - { - x = c(x,x) - prob = c(prob,prob) - if(is.null(size)) - { - size = 1 - } - if(replace == FALSE & size > 1) - { - stop('It is not possible to sample without replacement multiple times from a single item.') - } + if(length(x) == 1) + { + x = c(x,x) + prob = c(prob,prob) + if(is.null(size)) + { + size = 1 + } + if(replace == FALSE & size > 1) + { + stop('It is not possible to sample without replacement multiple times from a single item.') } - sam = sample(x,size,replace,prob) - return(sam) + } + sam = sample(x,size,replace,prob) + return(sam) } #' Carries out optimization using a simplex algorithm (finding a minimum) @@ -536,26 +536,26 @@ simplex = function(fun,trparsopt,optimpars,...) reltolf = optimpars[2] abstolx = optimpars[3] maxiter = optimpars[4] - + ## Setting up initial simplex v = t(matrix(rep(trparsopt,each = numpar + 1),nrow = numpar + 1)) for(i in 1:numpar) { - parsoptff = 1.05 * untransform_pars(trparsopt[i]) - trparsoptff = transform_pars(parsoptff) - fac = trparsoptff/trparsopt[i] - if(v[i,i + 1] == 0) - { - v[i,i + 1] = 0.00025 - } else { - v[i,i + 1] = v[i,i + 1] * min(1.05,fac) - } + parsoptff = 1.05 * untransform_pars(trparsopt[i]) + trparsoptff = transform_pars(parsoptff) + fac = trparsoptff/trparsopt[i] + if(v[i,i + 1] == 0) + { + v[i,i + 1] = 0.00025 + } else { + v[i,i + 1] = v[i,i + 1] * min(1.05,fac) + } } fv = rep(0,numpar + 1) for(i in 1:(numpar + 1)) { - fv[i] = -fun(trparsopt = v[,i], ...) + fv[i] = -fun(trparsopt = v[,i], ...) } how = "initial" @@ -563,7 +563,7 @@ simplex = function(fun,trparsopt,optimpars,...) string = itercount for(i in 1:numpar) { - string = paste(string, untransform_pars(v[i,1]), sep = " ") + string = paste(string, untransform_pars(v[i,1]), sep = " ") } string = paste(string, -fv[1], how, "\n", sep = " ") cat(string) @@ -572,9 +572,9 @@ simplex = function(fun,trparsopt,optimpars,...) tmp = order(fv) if(numpar == 1) { - v = matrix(v[tmp],nrow = 1,ncol = 2) + v = matrix(v[tmp],nrow = 1,ncol = 2) } else { - v = v[,tmp] + v = v[,tmp] } fv = fv[tmp] @@ -588,100 +588,100 @@ simplex = function(fun,trparsopt,optimpars,...) while(itercount <= maxiter & ( ( is.nan(max(abs(fv - fv[1]))) | (max(abs(fv - fv[1])) - reltolf * abs(fv[1]) > 0) ) + ( (max(abs(v - v2) - reltolx * abs(v2)) > 0) | (max(abs(v - v2)) - abstolx > 0) ) ) ) { - ## Calculate reflection point - - if(numpar == 1) - { - xbar = v[1] - } else { - xbar = rowSums(v[,1:numpar])/numpar - } - xr = (1 + rh) * xbar - rh * v[,numpar + 1] - fxr = -fun(trparsopt = xr, ...) - - if(fxr < fv[1]) - { - ## Calculate expansion point - xe = (1 + rh * ch) * xbar - rh * ch * v[,numpar + 1] - fxe = -fun(trparsopt = xe, ...) - if(fxe < fxr) - { - v[,numpar + 1] = xe - fv[numpar + 1] = fxe - how = "expand" - } else { - v[,numpar + 1] = xr - fv[numpar + 1] = fxr - how = "reflect" - } - } else { - if(fxr < fv[numpar]) - { - v[,numpar + 1] = xr - fv[numpar + 1] = fxr - how = "reflect" - } else { - if(fxr < fv[numpar + 1]) - { - ## Calculate outside contraction point - xco = (1 + ps * rh) * xbar - ps * rh * v[,numpar + 1] - fxco = -fun(trparsopt = xco, ...) - if(fxco <= fxr) - { - v[,numpar + 1] = xco - fv[numpar + 1] = fxco - how = "contract outside" - } else { - how = "shrink" - } - } else { - ## Calculate inside contraction point - xci = (1 - ps) * xbar + ps * v[,numpar + 1] - fxci = -fun(trparsopt = xci, ...) - if(fxci < fv[numpar + 1]) - { - v[,numpar + 1] = xci - fv[numpar + 1] = fxci - how = "contract inside" - } else { - how = "shrink" - } - } - if(how == "shrink") - { - for(j in 2:(numpar + 1)) - { - - v[,j] = v[,1] + si * (v[,j] - v[,1]) - fv[j] = -fun(trparsopt = v[,j], ...) - } - } - } - } - tmp = order(fv) - if(numpar == 1) - { - v = matrix(v[tmp],nrow = 1,ncol = 2) - } else { - v = v[,tmp] - } - fv = fv[tmp] - itercount = itercount + 1 - string = itercount; - for(i in 1:numpar) - { - string = paste(string, untransform_pars(v[i,1]), sep = " ") - } - string = paste(string, -fv[1], how, "\n", sep = " ") - cat(string) - utils::flush.console() - v2 = t(matrix(rep(v[,1],each = numpar + 1),nrow = numpar + 1)) + ## Calculate reflection point + + if(numpar == 1) + { + xbar = v[1] + } else { + xbar = rowSums(v[,1:numpar])/numpar + } + xr = (1 + rh) * xbar - rh * v[,numpar + 1] + fxr = -fun(trparsopt = xr, ...) + + if(fxr < fv[1]) + { + ## Calculate expansion point + xe = (1 + rh * ch) * xbar - rh * ch * v[,numpar + 1] + fxe = -fun(trparsopt = xe, ...) + if(fxe < fxr) + { + v[,numpar + 1] = xe + fv[numpar + 1] = fxe + how = "expand" + } else { + v[,numpar + 1] = xr + fv[numpar + 1] = fxr + how = "reflect" + } + } else { + if(fxr < fv[numpar]) + { + v[,numpar + 1] = xr + fv[numpar + 1] = fxr + how = "reflect" + } else { + if(fxr < fv[numpar + 1]) + { + ## Calculate outside contraction point + xco = (1 + ps * rh) * xbar - ps * rh * v[,numpar + 1] + fxco = -fun(trparsopt = xco, ...) + if(fxco <= fxr) + { + v[,numpar + 1] = xco + fv[numpar + 1] = fxco + how = "contract outside" + } else { + how = "shrink" + } + } else { + ## Calculate inside contraction point + xci = (1 - ps) * xbar + ps * v[,numpar + 1] + fxci = -fun(trparsopt = xci, ...) + if(fxci < fv[numpar + 1]) + { + v[,numpar + 1] = xci + fv[numpar + 1] = fxci + how = "contract inside" + } else { + how = "shrink" + } + } + if(how == "shrink") + { + for(j in 2:(numpar + 1)) + { + + v[,j] = v[,1] + si * (v[,j] - v[,1]) + fv[j] = -fun(trparsopt = v[,j], ...) + } + } + } + } + tmp = order(fv) + if(numpar == 1) + { + v = matrix(v[tmp],nrow = 1,ncol = 2) + } else { + v = v[,tmp] + } + fv = fv[tmp] + itercount = itercount + 1 + string = itercount; + for(i in 1:numpar) + { + string = paste(string, untransform_pars(v[i,1]), sep = " ") + } + string = paste(string, -fv[1], how, "\n", sep = " ") + cat(string) + utils::flush.console() + v2 = t(matrix(rep(v[,1],each = numpar + 1),nrow = numpar + 1)) } if(itercount < maxiter) { - cat("Optimization has terminated successfully.","\n") + cat("Optimization has terminated successfully.","\n") } else { - cat("Maximum number of iterations has been exceeded.","\n") + cat("Maximum number of iterations has been exceeded.","\n") } out = list(par = v[,1], fvalues = -fv[1], conv = as.numeric(itercount > maxiter)) invisible(out) @@ -782,7 +782,7 @@ optimizer <- function( itermax = optimpars[4], packages = c('DDD')), fun = fun, - + ...))$optim outnew <- list(par = outnew$bestmem, fvalues = -outnew$bestval, conv = 0) } else if(substr(optimmethod,1,7) == 'optim::') @@ -864,17 +864,15 @@ untransform_pars <- function(trpars) check_probs <- function(loglik,probs,verbose) { - probs <- probs * (probs > 0) - if(is.na(sum(probs)) | is.nan(sum(probs))) - { + probs <- pmax(probs, 0) + if (any(is.na(probs)) | any(is.nan(probs))) { if(verbose) cat('NA or NaN issues encountered.\n') loglik <- -Inf - #probs <- rep(-Inf,length(probs)); this line may cause problems with subplex - } else if(sum(probs) <= 0) - { + probs <- rep(-Inf,length(probs)) + } else if (sum(probs) <= 0) { if(verbose) cat('Probabilities smaller than 0 encountered\n') loglik <- -Inf - #probs <- rep(-Inf,length(probs)); this line may cause problems with subplex + probs <- rep(-Inf, length(probs)) } else { loglik <- loglik + log(sum(probs)) probs <- probs/sum(probs) @@ -922,4 +920,110 @@ rng_respecting_sample <- function(x, size, replace, prob) { non_zero_prob <- prob[which_non_zero] non_zero_x <- x[which_non_zero] return(sample(x = non_zero_x, size = size, replace = replace, prob = non_zero_prob)) -} \ No newline at end of file +} + +#' Evaluate whether both speciation rates vary in specified DD model +#' +#' @param ddmodel a character or integer specifying a DD model, as described in +#' dd_loglik() and dd_ML() documentation +#' +#' @return TRUE if both speciation and extinction rates vary with N, FALSE otherwise +#' @author Theo Pannetier +#' @export +both_rates_vary <- function(ddmodel) { + return(ddmodel %in% c(5:8, 11:15)) +} + +#' Get carrying capacity from other parameters of the model +#' +#' Compute the carrying capacity K', the maximum diversity that is possible to reach +#' with the model, i.e. the value of N for which lambda(N) = 0. +#' This is different from K, the equilibrium diversity, i.e. the value of N for +#' which lambda(N) = mu(N). +#' +#' @param ddmodel a character or integer specifying a DD model, as described in +#' dd_loglik() and dd_ML() documentation +#' @param pars a numeric vector containing parameter values of the DD model. +#' \code{pars[1]} is the lambda0, \code{pars[2]} is mu0, \code{pars[3]} is K +#' and \code{pars[4]}, if relevant, is r. +#' +#' @export +#' @author Theo Pannetier +#' @return a numeric value, K' + +get_Kprime <- function(ddmodel, pars) { + la <- pars[1] + mu <- pars[2] + K <- pars[3] + if (ddmodel == 1) { + Kprime <- la / (la - mu) * K + } else if (ddmodel == 1.3) { + Kprime <- K + } else if (ddmodel %in% c(5, 6, 11)) { + r <- pars[4] + Kprime <- (1 + r) * la / (la - mu) * K + } else { + Kprime <- Inf + } + return(Kprime) +} + +#' Helper function to quickly find what a ddmodel contains +#' +#' Given a ddmodel code, returns a short description of the contents of the model +#' +#' @param ddmodel a integer code between 1 and 15, corresponding to one of the +#' diversity-dependent models taken as input e.g. in [dd_ML()] and [dd_loglik()]. +#' @param short logical. If FALSE (the default), returns a description of the +#' diversity-dependent functions used for speciation and extinction. If TRUE, +#' instead returns a 2-letter code (speciation + extinction) summarising these +#' functions: L for liner DD, P for a power DD function, X for an exponential DD +#' function, and C for a constant rate (no DD). +#' +#' @export +#' @author Theo Pannetier +what_is_this_ddmodel <- function(ddmodel, short = FALSE) { + if (!ddmodel %in% 1:15) { + stop("ddmodel not found") + } + if (!short) { + switch( + as.character(ddmodel), + "1" = "linear DD on speciation, constant-rate extinction", + "2" = "exponential (power function) DD on speciation, constant-rate extinction", + "3" = "constant-rate speciation, linear DD on extinction", + "4" = "constant-rate speciation, exponential (power function) DD on extinction", + "5" = "linear DD on speciation, linear DD on extinction", + "6" = "linear DD on speciation, exponential DD (power function) on extinction", + "7" = "exponential (power function) DD on speciation, exponential (power function) DD on extinction", + "8" = "exponential (power function) DD on speciation, linear DD on extinction", + "9" = "exponential (exponential function) DD on speciation, constant-rate extinction", + "10" = "constant-rate speciation, exponential (exponential function) DD on extinction", + "11" = "linear DD on speciation, exponential (exponential function) DD on extinction", + "12" = "exponential (exponential function) DD on speciation, exponential (exponential function) DD on extinction", + "13" = "exponential (exponential function) DD on speciation, linear DD on speciation", + "14" = "exponential (exponential function) DD on speciation, exponential (power function) DD on extinction", + "15" = "exponential (power function) DD on speciation, exponential (exponential function) DD on extinction" + ) + } else { + switch ( + as.character(ddmodel), + "1" = "LC", + "2" = "PC", + "3" = "CL", + "4" = "CP", + "5" = "LL", + "6" = "LP", + "7" = "PP", + "8" = "PL", + "9" = "XC", + "10" = "CX", + "11" = "LX", + "12" = "XX", + "13" = "XL", + "14" = "XP", + "15" = "PX" + ) + } +} + diff --git a/man/both_rates_vary.Rd b/man/both_rates_vary.Rd new file mode 100644 index 0000000..620bea2 --- /dev/null +++ b/man/both_rates_vary.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dd_utils.R +\name{both_rates_vary} +\alias{both_rates_vary} +\title{Evaluate whether both speciation rates vary in specified DD model} +\usage{ +both_rates_vary(ddmodel) +} +\arguments{ +\item{ddmodel}{a character or integer specifying a DD model, as described in +dd_loglik() and dd_ML() documentation} +} +\value{ +TRUE if both speciation and extinction rates vary with N, FALSE otherwise +} +\description{ +Evaluate whether both speciation rates vary in specified DD model +} +\author{ +Theo Pannetier +} diff --git a/man/dd_KI_loglik.Rd b/man/dd_KI_loglik.Rd index ab9f9a6..f2ca0cc 100644 --- a/man/dd_KI_loglik.Rd +++ b/man/dd_KI_loglik.Rd @@ -100,11 +100,11 @@ shift in parameters. } \examples{ -pars1 = c(0.25,0.12,25.51,1.0,0.16,8.61,9.8) -pars2 = c(200,1,0,18.8,1,2) -missnumspec = 0 -brtsM = c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) -brtsS = c(9.6,8.6,7.4,4.9,2.5) +pars1 <- c(0.25,0.12,25.51,1.0,0.16,8.61,9.8) +pars2 <- c(200,1,0,18.8,1,2) +missnumspec <- 0 +brtsM <- c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) +brtsS <- c(9.6,8.6,7.4,4.9,2.5) dd_KI_loglik(pars1,pars2,brtsM,brtsS,missnumspec) } diff --git a/man/dd_ML.Rd b/man/dd_ML.Rd index 5505077..f435fce 100644 --- a/man/dd_ML.Rd +++ b/man/dd_ML.Rd @@ -9,7 +9,7 @@ dd_ML( brts, initparsopt = initparsoptdefault(ddmodel, brts, missnumspec), idparsopt = 1:length(initparsopt), - idparsfix = (1:(3 + (ddmodel == 5)))[-idparsopt], + idparsfix = (1:(3 + both_rates_vary(ddmodel)))[-idparsopt], parsfix = parsfixdefault(ddmodel, brts, missnumspec, idparsopt), res = 10 * (1 + length(brts) + missnumspec), ddmodel = 1, @@ -61,19 +61,41 @@ maximum); lambda = lambda0 * S/(S + K') where S is species richness\cr \code{ddmodel == 1.5} : positive and negative dependence in speciation rate with parameter K' (= diversity where speciation = 0); lambda = lambda0 * S/K' * (1 - S/K') where S is species richness\cr -\code{ddmodel == 2} : exponential dependence in speciation rate with parameter +\code{ddmodel == 2} : exponential dependence (power function) in speciation rate with parameter K (= diversity where speciation = extinction)\cr -\code{ddmodel == 2.1} : variant of exponential dependence in speciation rate +\code{ddmodel == 2.1} : variant of exponential dependence (power function) in speciation rate with offset at infinity\cr \code{ddmodel == 2.2} : 1/n dependence in speciation rate\cr -\code{ddmodel == 2.3} : exponential dependence in speciation rate with parameter x (= +\code{ddmodel == 2.3} : exponential dependence (power function) in speciation rate with parameter x (= exponent)\cr \code{ddmodel == 3} : linear dependence in extinction rate \cr -\code{ddmodel == 4} : exponential dependence in extinction rate \cr +\code{ddmodel == 4} : exponential dependence (power function) in extinction rate \cr \code{ddmodel == 4.1} : variant of exponential dependence in extinction rate with offset at infinity \cr \code{ddmodel == 4.2} : 1/n dependence in extinction rate with offset at infinity \cr \code{ddmodel == 5} : linear -dependence in speciation and extinction rate \cr} +dependence in speciation and extinction rate \cr +\code{ddmodel == 5} : linear dependence in speciation and +extinction rate \cr +\code{ddmodel == 6} : linear dependence in speciation rate, exponential +dependence (power function) in extinction rate \cr +\code{ddmodel == 7} : exponential dependence (power function) in speciation +and extinction rate \cr +\code{ddmodel == 8} : exponential dependence (power function) in speciation rate, +linear dependence in extinction rate \cr +\code{ddmodel == 9} : exponential dependence (exponential function) in speciation, +constant-rate extinction \cr +\code{ddmodel == 10} : constant-rate speciation, exponential dependence +(exponential function) in extinction \cr +\code{ddmodel == 11} : linear dependence in speciation, exponential +dependence (exponential function) in extinction\cr +\code{ddmodel == 12} : exponential dependence (exponential function) in +speciation and extinction\cr +\code{ddmodel == 13} : exponential dependence (exponential function) in +speciation, linear dependence in extinction \cr +\code{ddmodel == 14} : exponential dependence (exponential function) in +speciation, exponential dependence (power function) in extinction \cr +\code{ddmodel == 15} : exponential dependence (power function) in +speciation, exponential dependence (exponential function) in extinction \cr} \item{missnumspec}{The number of species that are in the clade but missing in the phylogeny} diff --git a/man/dd_MS_ML.Rd b/man/dd_MS_ML.Rd index 97e95c0..8bbe676 100644 --- a/man/dd_MS_ML.Rd +++ b/man/dd_MS_ML.Rd @@ -26,7 +26,7 @@ dd_MS_ML( changeloglikifnoconv = FALSE, optimmethod = "subplex", num_cycles = 1, - methode = "ode45", + methode = "odeint::runge_kutta_cash_karp54", correction = FALSE, verbose = FALSE ) diff --git a/man/dd_loglik.Rd b/man/dd_loglik.Rd index ff85205..9945c17 100644 --- a/man/dd_loglik.Rd +++ b/man/dd_loglik.Rd @@ -4,7 +4,7 @@ \alias{dd_loglik} \title{Loglikelihood for diversity-dependent diversification models} \usage{ -dd_loglik(pars1, pars2, brts, missnumspec, methode = "analytical") +dd_loglik(pars1, pars2, brts, missnumspec = 0, methode = "analytical") } \arguments{ \item{pars1}{Vector of parameters: @@ -18,11 +18,9 @@ rate) \cr \cr \code{pars2[1]} sets the maximum number of species for which a probability must be computed. This must be larger than 1 + missnumspec + length(brts). -\cr \cr \code{pars2[2]} -sets the model of diversity-dependence: -\cr - \code{pars2[2] == 1} linear -dependence in speciation rate with parameter K (= diversity where speciation -= extinction) +\cr \cr \code{pars2[2]} sets the model of diversity-dependence: +\cr - \code{pars2[2] == 1} linear dependence in speciation rate with +parameter K (= diversity where speciation = extinction) \cr - \code{pars2[2] == 1.3} linear dependence in speciation rate with parameter K' (= diversity where speciation = 0) \cr - \code{pars2[2] == 1.4} : positive diversity-dependence in speciation rate @@ -31,28 +29,43 @@ maximum); lambda = lambda0 * S/(S + K') where S is species richness \cr - \code{pars2[2] == 1.5} : positive and negative diversity-dependence in speciation rate with parameter K' (= diversity where speciation = 0); lambda = lambda0 * S/K' * (1 - S/K') where S is species richness -\cr - \code{pars2[2] == 2} exponential dependence in speciation rate with +\cr - \code{pars2[2] == 2} exponential dependence (power function) in speciation rate with parameter K (= diversity where speciation = extinction) \cr - \code{pars2[2] -== 2.1} variant of exponential dependence in speciation rate with offset at +== 2.1} variant of exponential dependence (power function) in speciation rate with offset at infinity \cr - \code{pars2[2] == 2.2} 1/n dependence in speciation rate -\cr - \code{pars2[2] == 2.3} exponential dependence in speciation rate with +\cr - \code{pars2[2] == 2.3} exponential dependence (power function) in speciation rate with parameter x (= exponent) -\cr - \code{pars2[2] == 3} linear dependence in -extinction rate -\cr - \code{pars2[2] == 4} exponential dependence in -extinction rate -\cr - \code{pars2[2] == 4.1} variant of exponential -dependence in extinction rate with offset at infinity -\cr - \code{pars2[2] == -4.2} 1/n dependence in extinction rate -\cr - \code{pars2[2] == 5} linear -dependence in speciation and extinction rate -\cr \cr \code{pars2[3]} sets -the conditioning: -\cr - \code{pars2[3] == 0} conditioning on stem or crown -age +\cr - \code{pars2[2] == 3} linear dependence in extinction rate +\cr - \code{pars2[2] == 4} exponential dependence (power function) in extinction rate +\cr - \code{pars2[2] == 4.1} variant of exponential dependence (power function) in extinction +rate with offset at infinity +\cr - \code{pars2[2] == 4.2} 1/n dependence in extinction rate +\cr - \code{pars2[2] == 5} linear dependence in speciation and extinction rate +\cr - \code{pars2[2] == 6} linear dependence in speciation rate, exponential +dependence (power function) in extinction rate \cr +\cr - \code{pars2[2] == 7} exponential dependence (power function) in speciation +and extinction rate \cr +\cr - \code{pars2[2] == 8} exponential dependence (power function) in speciation rate, +linear dependence in extinction rate \cr +\cr - \code{pars2[2] == 9} exponential dependence (exponential function) in speciation, +constant-rate extinction \cr +\cr - \code{pars2[2] == 10} constant-rate speciation, exponential dependence +(exponential function) in extinction \cr +\cr - \code{pars2[2] == 11} linear dependence in speciation, exponential +dependence (exponential function) in extinction\cr +\cr - \code{pars2[2] == 12} exponential dependence (exponential function) in +speciation and extinction\cr +\cr - \code{pars2[2] == 13} exponential dependence (exponential function) in +speciation, linear dependence in extinction \cr +\cr - \code{pars2[2] == 14} exponential dependence (exponential function) in +speciation, exponential dependence (power function) in extinction \cr +\cr - \code{pars2[2] == 15} exponential dependence (power function) in +speciation, exponential dependence (exponential function) in extinction \cr + +\cr \cr \code{pars2[3]} sets the conditioning: +\cr - \code{pars2[3] == 0} conditioning on stem or crown age \cr - \code{pars2[3] == 1} conditioning on stem or crown age and non-extinction of the phylogeny \cr - \code{pars2[3] == 2} conditioning on diff --git a/man/dd_multiple_KI_loglik.Rd b/man/dd_multiple_KI_loglik.Rd index 1d0a8ea..8d8b983 100644 --- a/man/dd_multiple_KI_loglik.Rd +++ b/man/dd_multiple_KI_loglik.Rd @@ -16,7 +16,7 @@ dd_multiple_KI_loglik( ) } \arguments{ -\item{pars1_list}{list of paramater sets one for each rate regime (subclade). +\item{pars1_list}{list of parameter sets one for each rate regime (subclade). The parameters are: lambda (speciation rate), mu (extinction rate), and K (clade-level carrying capacity).} diff --git a/man/get_Kprime.Rd b/man/get_Kprime.Rd new file mode 100644 index 0000000..2f1ff3e --- /dev/null +++ b/man/get_Kprime.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dd_utils.R +\name{get_Kprime} +\alias{get_Kprime} +\title{Get carrying capacity from other parameters of the model} +\usage{ +get_Kprime(ddmodel, pars) +} +\arguments{ +\item{ddmodel}{a character or integer specifying a DD model, as described in +dd_loglik() and dd_ML() documentation} + +\item{pars}{a numeric vector containing parameter values of the DD model. +\code{pars[1]} is the lambda0, \code{pars[2]} is mu0, \code{pars[3]} is K +and \code{pars[4]}, if relevant, is r.} +} +\value{ +a numeric value, K' +} +\description{ +Compute the carrying capacity K', the maximum diversity that is possible to reach +with the model, i.e. the value of N for which lambda(N) = 0. +This is different from K, the equilibrium diversity, i.e. the value of N for +which lambda(N) = mu(N). +} +\author{ +Theo Pannetier +} diff --git a/man/what_is_this_ddmodel.Rd b/man/what_is_this_ddmodel.Rd new file mode 100644 index 0000000..55b8e1f --- /dev/null +++ b/man/what_is_this_ddmodel.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dd_utils.R +\name{what_is_this_ddmodel} +\alias{what_is_this_ddmodel} +\title{Helper function to quickly find what a ddmodel contains} +\usage{ +what_is_this_ddmodel(ddmodel, short = FALSE) +} +\arguments{ +\item{ddmodel}{a integer code between 1 and 15, corresponding to one of the +diversity-dependent models taken as input e.g. in [dd_ML()] and [dd_loglik()].} + +\item{short}{logical. If FALSE (the default), returns a description of the +diversity-dependent functions used for speciation and extinction. If TRUE, +instead returns a 2-letter code (speciation + extinction) summarising these +functions: L for liner DD, P for a power DD function, X for an exponential DD +function, and C for a constant rate (no DD).} +} +\description{ +Given a ddmodel code, returns a short description of the contents of the model +} +\author{ +Theo Pannetier +} diff --git a/tests/testthat/test-optimizer.R b/tests/testthat/test-optimizer.R index d99a2f1..19a9a44 100644 --- a/tests/testthat/test-optimizer.R +++ b/tests/testthat/test-optimizer.R @@ -2,15 +2,15 @@ test_that("optimizer works", { brts <- 1:10 initparsopt <- c(0.3,0.05,12) idparsopt <- 1:3 - res_simplex <- DDD::dd_ML(brts = brts, + res_simplex <- dd_ML(brts = brts, initparsopt = initparsopt, idparsopt = idparsopt, optimmethod = 'simplex') - res_subplex <- DDD::dd_ML(brts = brts, + res_subplex <- dd_ML(brts = brts, initparsopt = initparsopt, idparsopt = idparsopt, optimmethod = 'subplex') - #res_nloptr <- DDD::dd_ML(brts = brts, + #res_nloptr <- dd_ML(brts = brts, # initparsopt = initparsopt, # idparsopt = idparsopt, # optimmethod = 'NLOPT_LN_SBPLX') diff --git a/tests/testthat/test_DDD.R b/tests/testthat/test_DDD.R index 08dc7df..65ae205 100644 --- a/tests/testthat/test_DDD.R +++ b/tests/testthat/test_DDD.R @@ -17,10 +17,10 @@ test_that("DDD works", { missnumspec = 0 methode = 'analytical' - r0 <- DDD:::dd_loglik(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,methode = methode) - r1 <- DDD:::dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,rhs_func_name = 'dd_loglik_rhs',methode = methode) - r2 <- DDD:::dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,rhs_func_name = 'dd_loglik_rhs_FORTRAN',methode = methode) - r3 <- DDD:::dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,methode = 'analytical') + r0 <- dd_loglik(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,methode = methode) + r1 <- dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,rhs_func_name = 'dd_loglik_rhs',methode = methode) + r2 <- dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,rhs_func_name = 'dd_loglik_rhs_FORTRAN',methode = methode) + r3 <- dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,methode = 'analytical') testthat::expect_equal(r0,r2,tolerance = .00001) testthat::expect_equal(r1,r2,tolerance = .00001) @@ -30,6 +30,23 @@ test_that("DDD works", { r4 <- DDD::dd_SR_loglik(pars1 = c(0.2,0.1,50,0.2,0.1,70,5), pars2 = c(100,1,1,1,0,2), brts = 1:10, missnumspec = 0) testthat::expect_equal(-27.37304,r4,tolerance = .000001) + brts = 1:5 + pars2 = c(100,1,3,0,0,2) + r5 <- dd_loglik(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,methode = methode) + r6 <- dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,rhs_func_name = 'dd_loglik_bw_rhs',methode = methode) + r7 <- dd_loglik_test(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,rhs_func_name = 'dd_loglik_bw_rhs_FORTRAN',methode = 'analytical') + + testthat::expect_equal(r5,r6,tolerance = .00001) + testthat::expect_equal(r5,r7,tolerance = .01) + expect_equal_x64(-8.579058,r7,tolerance = .00001) #was -8.582413 before + + pars1 = c(0.2,0.05,1000000) + pars2 = c(1000,1,1,0,0,2) + brts = 1:10 + r8 <- dd_loglik(pars1 = pars1,pars2 = pars2,brts = brts,missnumspec = missnumspec,methode = methode) + r9 <- dd_loglik(pars1 = c(pars1[1:2],Inf),pars2 = pars2,brts = brts,missnumspec = missnumspec,methode = methode) + expect_equal_x64(r8,r9,tolerance = .00001) + #skip the following tests due to numerical issues with cond = 3 #brts = 1:5 #pars2 = c(100,1,3,0,0,2) @@ -67,7 +84,7 @@ test_that("DDD_KI works", # brts = brts, # cond = 0, # n_max = 1e3 - #); + #) high_k <- 1e7 pars1 <- c(pars[1], pars[2], high_k, pars[3], pars[4], high_k, brts[[2]][1]) pars2 <- c(500,1,0,brts[[2]][1],0,2,1.5) @@ -80,7 +97,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test,-24.4171970357049624,tolerance = .000001) ddd_test2 <- DDD::dd_KI_loglik( pars1 = pars1, @@ -89,7 +106,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'analytical' - ); + ) testthat::expect_equal(ddd_test,ddd_test2,tolerance = .000001) low_k <- 20 @@ -101,7 +118,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test,-21.1781625797899231,tolerance = .000001) ddd_test2 <- DDD::dd_KI_loglik( pars1 = pars1, @@ -110,7 +127,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'analytical' - ); + ) testthat::expect_equal(ddd_test,ddd_test2,tolerance = .000001) ddd_test3 <- DDD::dd_KI_loglik( @@ -120,7 +137,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 3, methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test3,-19.6273910107265408,tolerance = .000001) ddd_test03 <- DDD::dd_KI_loglik( @@ -130,7 +147,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = c(0,3), methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test03,-21.4981352311200595,tolerance = .000001) ddd_test12 <- DDD::dd_KI_loglik( @@ -140,7 +157,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = c(1,2), methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test12,-20.7167138427128776,tolerance = .000001) ddd_test21 <- DDD::dd_KI_loglik( @@ -150,7 +167,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = c(2,1), methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test21,-20.0405933298708874,tolerance = .000001) ddd_test30 <- DDD::dd_KI_loglik( @@ -160,7 +177,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = c(3,0), methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test30,-19.4834201422017124,tolerance = .000001) #testthat::expect_equal(ddd_test3,log(exp(ddd_test03) + exp(ddd_test12) + exp(ddd_test21) + exp(ddd_test30)),tolerance = .000001) @@ -174,7 +191,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test,-20.5299241171281643,tolerance = .000001) ddd_test2 <- DDD::dd_KI_loglik( pars1 = pars1, @@ -183,7 +200,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'analytical' - ); + ) testthat::expect_equal(ddd_test,ddd_test2,tolerance = .000001) pars2[3] <- 4 @@ -194,7 +211,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test,-20.2509115267895616,tolerance = .000001) pars2[3] <- 5 @@ -205,7 +222,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'odeint::runge_kutta_cash_karp54' - ); + ) testthat::expect_equal(ddd_test,-20.1686905596579997,tolerance = .000001) cond <- 1 @@ -216,7 +233,7 @@ test_that("DDD_KI works", # brts = brts, # cond = cond, # n_max = 1e3 - # ); + # ) t_d <- brts[[2]][1] tsplit <- min(abs(brts[[1]][abs(brts[[1]]) > t_d])) high_k <- 1e7 @@ -231,7 +248,7 @@ test_that("DDD_KI works", brtsS = brtsS, missnumspec = 0, methode = 'analytical' - ); + ) testthat::expect_equal(ddd_test1,-28.5415506633517460,tolerance = .000001) cond <- 0 pars2 <- c(200,1,cond,brts[[2]][1],0,2,1.5) @@ -267,9 +284,9 @@ test_that("DDD_KI works", test_that("conditioning_DDD_KI works", { skip_if(Sys.getenv("CI") == "", message = "Run only on CI") - ts <- seq(-9,-1,2); - p1 <- rep(0,5); - p2 <- rep(0,5); + ts <- seq(-9,-1,2) + p1 <- rep(0,5) + p2 <- rep(0,5) pars1_list <- list(c(0.5,0.4,Inf),c(0,0,Inf)) reltol <- 1e-8 abstol <- 1e-8 @@ -278,7 +295,7 @@ test_that("conditioning_DDD_KI works", for(i in 1:5) { brts_k_list <- list(rbind(c(-10,ts[i],0),c(2,1,1)),rbind(c(ts[i],0),c(1,1))) - p1[i] <- DDD:::dd_multiple_KI_logliknorm(brts_k_list = brts_k_list, + p1[i] <- dd_multiple_KI_logliknorm(brts_k_list = brts_k_list, pars1_list = pars1_list, pars2 = c(200,1,5,NA,1,2,3), loglik = 0, @@ -305,7 +322,7 @@ test_that("conditioning_DDD_KI works", pars2 <- c(500,1,5,NA,1,2,3) lx_list <- list(pars2[1],pars2[1]) brts_k_list <- list(rbind(sort(c(-10:-6,-3,-1,0)),c(2,3,4,5,6,5,6,6)),rbind(c(-3,-2,0),c(1,2,2))) - logliknorm1 <- DDD:::dd_KI_logliknorm(brts_k_list = brts_k_list, + logliknorm1 <- dd_KI_logliknorm(brts_k_list = brts_k_list, pars1_list = pars1_list, loglik = 0, cond = 5, @@ -328,7 +345,7 @@ test_that("conditioning_DDD_KI works", pars2 <- c(500,1,5,NA,1,2,3) lx_list <- list(pars2[1],pars2[1]) brts_k_list <- list(rbind(c(-10:-6,-3,-1,0),c(2,3,4,5,6,5,6,6)),rbind(c(-3,-2,0),c(1,2,2))) - logliknorm1 <- DDD:::dd_KI_logliknorm(brts_k_list = brts_k_list, + logliknorm1 <- dd_KI_logliknorm(brts_k_list = brts_k_list, pars1_list = pars1_list, loglik = 0, cond = 5, diff --git a/tests/testthat/test_ddmodels.R b/tests/testthat/test_ddmodels.R new file mode 100644 index 0000000..ef71d2b --- /dev/null +++ b/tests/testthat/test_ddmodels.R @@ -0,0 +1,214 @@ +context("test_ddmodels") + +# Case 1.: 0 < r < Inf (or 0 < phi < 1) +pars_set1 <- c( + "lambda_0" = 0.8, + "mu_0" = 0.2, + "K" = 20, + "r" = 1/3 # corresponds to phi = 1/4 +) +# Rates obtained on paper +exptd_rates_set1 <- list( + "lambda_cst" = function(N) 0.8, # not with phi != 1 + "mu_cst" = function(N) 0.2, # not with phi != 0 + "lambda_lin" = function(N) pmax(0.8 - 0.0225 * N, 0), + "mu_lin" = function(N) 0.2 + 0.0075 * N, + "lambda_pow" = function(N) pmax(0.8 * N ^ (-log(0.8 / 0.35) / log(20)), 0), + "mu_pow" = function(N) 0.2 * N ^ (log(1.75) / log(20)), + "lambda_exp" = function(N) pmax(0.8 * (7 / 16) ^ (N / 20), 0), + 'mu_exp' = function(N) 0.2 * (7 / 4) ^ (N / 20) +) + +# Case 2.: r = 0 (phi = 0) +pars_set2 <- c( + "lambda_0" = 0.8, + "mu_0" = 0.2, + "K" = 20, + "r" = 0 +) +exptd_rates_set2 <- list( + "lambda_cst" = function(N) rep(0.8, length(N)), # not with phi != 1 + "mu_cst" = function(N) rep(0.2, length(N)), # not with phi != 0 + "lambda_lin" = function(N) pmax(0.8 - 0.03 * N, 0), + "mu_lin" = function(N) rep(0.2, length(N)), + "lambda_pow" = function(N) pmax(0.8 * N ^ (-log(4) / log(20)), 0), + "mu_pow" = function(N) rep(0.2, length(N)), + "lambda_exp" = function(N) pmax(0.8 * (1 / 4) ^ (N / 20), 0), + 'mu_exp' = function(N) rep(0.2, length(N)) +) +# Case 3.: r = Inf (phi = 1) +pars_set3 <- c( + "lambda_0" = 0.8, + "mu_0" = 0.2, + "K" = 20, + "r" = Inf +) +exptd_rates_set3 <- list( + "lambda_cst" = function(N) rep(0.8, length(N)), # not with phi != 1 + "mu_cst" = function(N) rep(0.2, length(N)), # not with phi != 0 + "lambda_lin" = function(N) rep(0.8, length(N)), + "mu_lin" = function(N) 0.2 + 0.03 * N, + "lambda_pow" = function(N) rep(0.8, length(N)), + "mu_pow" = function(N) 0.2 * N ^ (log(4) / log(20)), + "lambda_exp" = function(N) rep(0.8, length(N)), + 'mu_exp' = function(N) 0.2 * 4 ^ (N / 20) +) + +# Declare test functions +## Match a ddmodel with the corresponding pair of speciation and extinction functions +match_exptd_rates <- function(ddmodel, exptd_rates_set, n_seq) { + rates_ls <- switch( + as.character(ddmodel), + "1" = list("la_N" = exptd_rates_set$lambda_lin(n_seq), "mu_N" = exptd_rates_set$mu_cst(n_seq)), + "2" = list("la_N" = exptd_rates_set$lambda_pow(n_seq), "mu_N" = exptd_rates_set$mu_cst(n_seq)), + "3" = list("la_N" = exptd_rates_set$lambda_cst(n_seq), "mu_N" = exptd_rates_set$mu_lin(n_seq)), + "4" = list("la_N" = exptd_rates_set$lambda_cst(n_seq), "mu_N" = exptd_rates_set$mu_pow(n_seq)), + "5" = list("la_N" = exptd_rates_set$lambda_lin(n_seq), "mu_N" = exptd_rates_set$mu_lin(n_seq)), + "6" = list("la_N" = exptd_rates_set$lambda_lin(n_seq), "mu_N" = exptd_rates_set$mu_pow(n_seq)), + "7" = list("la_N" = exptd_rates_set$lambda_pow(n_seq), "mu_N" = exptd_rates_set$mu_pow(n_seq)), + "8" = list("la_N" = exptd_rates_set$lambda_pow(n_seq), "mu_N" = exptd_rates_set$mu_lin(n_seq)), + "9" = list("la_N" = exptd_rates_set$lambda_exp(n_seq), "mu_N" = exptd_rates_set$mu_cst(n_seq)), + "10" = list("la_N" = exptd_rates_set$lambda_cst(n_seq), "mu_N" = exptd_rates_set$mu_exp(n_seq)), + "11" = list("la_N" = exptd_rates_set$lambda_lin(n_seq), "mu_N" = exptd_rates_set$mu_exp(n_seq)), + "12" = list("la_N" = exptd_rates_set$lambda_exp(n_seq), "mu_N" = exptd_rates_set$mu_exp(n_seq)), + "13" = list("la_N" = exptd_rates_set$lambda_exp(n_seq), "mu_N" = exptd_rates_set$mu_lin(n_seq)), + "14" = list("la_N" = exptd_rates_set$lambda_exp(n_seq), "mu_N" = exptd_rates_set$mu_pow(n_seq)), + "15" = list("la_N" = exptd_rates_set$lambda_pow(n_seq), "mu_N" = exptd_rates_set$mu_exp(n_seq)) + ) + return(rates_ls) +} + +## Test function; compare DDD output with rates obtained on paper +test_dd_loglik_rhs_precomp <- function(ddmodel, pars_set, exptd_rates_set) { + # global variables + N <- pars_set["K"] + x <- rep(NA, 10) # length of the Q_n vector + n_seq <- c(0, 0:(length(x) + 2 * N)) # based on internal code, not sure why + lnn <- length(n_seq) + + exptd_rates <- match_exptd_rates( + ddmodel = ddmodel, + exptd_rates_set = exptd_rates_set, + n_seq = n_seq + ) + ddd_output <- dd_loglik_rhs_precomp( + pars = c("pars" = pars_set, "k" = N, "ddmodel" = ddmodel), x = x + ) + names(ddd_output) <- NULL + ddd_rates <- list("la_N" = ddd_output[1:lnn], "mu_N" = ddd_output[(lnn + 1):(2 * lnn)]) + # Test + cat(paste("Testing ddmodel =", ddmodel, "\n")) + expect_equal(ddd_rates, exptd_rates) +} + +test_lambdamu <- function(ddmodel, pars_set, exptd_rates_set) { + # global variables + N <- pars_set["K"] + x <- rep(NA, 10) # length of the Q_n vector + n_seq <- c(0, 0:(length(x) + 2 * N)) # based on internal code, not sure why + lnn <- length(n_seq) + + exptd_rates <- match_exptd_rates( + ddmodel = ddmodel, + exptd_rates_set = exptd_rates_set, + n_seq = n_seq + ) + ddd_rates <- lambdamu(n = n_seq, pars = pars_set, ddep = ddmodel) + names(ddd_rates) <- c("la_N", "mu_N") + names(ddd_rates$la_N) <- NULL + names(ddd_rates$mu_N) <- NULL + # Test + cat(paste("Testing ddmodel =", ddmodel, "\n")) + expect_equal(ddd_rates, exptd_rates) +} + +test_dd_lamuN <- function(ddmodel, pars_set, exptd_rates_set) { + # global variables + N <- pars_set["K"] + x <- rep(NA, 10) # length of the Q_n vector + n_seq <- c(0, 0:(length(x) + 2 * N)) # based on internal code, not sure why + + exptd_rates <- match_exptd_rates( + ddmodel = ddmodel, + exptd_rates_set = exptd_rates_set, + n_seq = n_seq + ) + ddd_rates <- list( + "la_N" = unlist(lapply(n_seq, function(n) { + dd_lamuN(ddmodel = ddmodel, pars = pars_set, N = n)[1] + }), use.names = FALSE), + "mu_N" = unlist(lapply(n_seq, function(n) { + dd_lamuN(ddmodel = ddmodel, pars = pars_set, N = n)[2] + }), use.names = FALSE) + ) + # Test + cat(paste("Testing ddmodel =", ddmodel, "\n")) + expect_equal(ddd_rates, exptd_rates) +} + +test_that("set1", { + ddmodels <- c(5:8, 11:15) + invisible(lapply( + ddmodels, + test_dd_loglik_rhs_precomp, + pars_set = pars_set1, + exptd_rates_set = exptd_rates_set1 + )) + invisible(lapply( + ddmodels, + test_lambdamu, + pars_set = pars_set1, + exptd_rates_set = exptd_rates_set1 + )) + invisible(lapply( + ddmodels, + test_dd_lamuN, + pars_set = pars_set1, + exptd_rates_set = exptd_rates_set1 + )) +}) + +test_that("set2", { + ddmodels <- c(1, 5:9, 11:15) + invisible(lapply( + ddmodels, + test_dd_loglik_rhs_precomp, + pars_set = pars_set2, + exptd_rates_set = exptd_rates_set2 + )) + invisible(lapply( + ddmodels, + test_lambdamu, + pars_set = pars_set2, + exptd_rates_set = exptd_rates_set2 + )) + invisible(lapply( + ddmodels, + test_dd_lamuN, + pars_set = pars_set2, + exptd_rates_set = exptd_rates_set2 + )) +}) + +test_that("set3", { + ddmodels <- c(3, 5:8, 10:15) + invisible(lapply( + ddmodels, + test_dd_loglik_rhs_precomp, + pars_set = pars_set3, + exptd_rates_set = exptd_rates_set3 + )) + invisible(lapply( + ddmodels, + test_lambdamu, + pars_set = pars_set3, + exptd_rates_set = exptd_rates_set3 + )) + invisible(lapply( + ddmodels, + test_dd_lamuN, + pars_set = pars_set3, + exptd_rates_set = exptd_rates_set3 + )) +}) +