Nothing
DAISIE_loglik_CS_shift <- function(pars1,
pars2,
brts,
stac,
missnumspec,
methode = "lsodes",
abstolint = 1E-16,
reltolint = 1E-10,
verbose,
tshift) {
# stac = status of the clade formed by the immigrant
# . stac == 1 : immigrant is present but has not formed an extant clade
# . stac == 2 : immigrant is not present but has formed an extant clade
# . stac == 3 : immigrant is present and has formed an extant clade
# . stac == 4 : immigrant is present but has not formed an extant clade,
# and it is known when it immigrated.
# . stac == 5 : immigrant is not present and has not formed an extant clade,
# but only an endemic species
# . stac == 6 : like 2, but with max colonization time
# . stac == 7 : like 3, but with max colonization time
# . stac == 8 : like 1, but with min colonization time
# . stac == 9 : like 5, but with min colonization time
# warn if laa becomes Inf
if (any(is.infinite(pars1)) ) {
if (verbose) {
message('One of the parameters is infinite.')
}
}
if(is.na(pars2[4]))
{
pars2[4] <- 0
}
ddep <- pars2[2]
K <- pars1[3]
# if (!is.na(pars2[5])) {
# K <- K * pars1[8]
# }
if(length(pars1) == 6) {
probability_of_init_presence <- pars1[6]
pars1 <- pars1[-6]
} else {
probability_of_init_presence <- 0
}
brts <- -sort(abs(as.numeric(brts)),decreasing = TRUE)
tshift <- -abs(tshift)
if(length(brts) == 1 & sum(brts == 0) == 1)
{
stop('The branching times contain only a 0. This means the island emerged at the present which is not allowed.');
loglik <- -Inf
return(loglik)
}
if (sum(brts == 0) == 0) {
brts[length(brts) + 1] <- 0
}
# for stac = 0, brts will contain origin of island and 0; length = 2;
# no. species should be 0
# for stac = 1, brts will contain origin of island, maximum colonization time
# (usually island age) and 0; length = 3; no. species should be 1
# for stac = 2, brts will contain origin of island, colonization event,
# branching times, 0; no. species should be no. branching times + 1
# for stac = 3, brts will contain origin of island, colonization event,
# branching times, 0; no. species should be no. branching times + 2
# for stac = 4, brts will contain origin of island, colonization event and 0;
# length = 3; no. species should be 1
# for stac = 5, brts will contain origin of island, maximum colonization time
# (usually island age), and 0; length = 2; number of species should be 1 (+ missing species)
# for stac = 6, brts will contain origin of island, maximum colonization time
# (usually island age), branching times and 0;
# number of species should be no. branching times + 1
# for stac = 7, brts will contain origin of island, maximum colonization time
# usually island age), branching times and 0;
# number of species should be no. branching times + 2
# for stac = 8, brts will contain origin of island, maximum colonization time
# usually island age), minimum colonization time and 0; length = 4;
# number of species should be 1
# for stac = 9, brts will contain origin of island, maximum colonization time
# usually island age), minimum colonization time and 0; length = 4;
# number of species should be 1
S <- 0 * (stac == 0) + (stac == 1 || stac == 4 || stac == 5 || stac == 8 || stac == 9) +
(length(brts) - 2) * (stac == 2) + (length(brts) - 1) * (stac == 3) +
(length(brts) - 2) * (stac == 6) + (length(brts) - 1) * (stac == 7)
S2 <- S - (stac == 1) - (stac == 3) - (stac == 4) - (stac == 7)
loglik <- -lgamma(S2 + missnumspec + 1) +
lgamma(S2 + 1) + lgamma(missnumspec + 1)
if (min(pars1) < 0) {
message("One or more parameters are negative.")
loglik <- -Inf
return(loglik)
}
if ((ddep == 1 | ddep == 11) & ceiling(K) < (S + missnumspec)) {
if (verbose) {
message('The proposed value of K is incompatible with the number of species
in the clade. Likelihood for this parameter set
will be set to -Inf. \n')
}
loglik <- -Inf
return(loglik)
}
# lac <- pars1[1]
# if(lac == Inf & missnumspec == 0 & length(pars1) == 5) {
# if(verbose) warning('Infinite lambda detected')
# loglik <- DAISIE_loglik_high_lambda(pars1, -brts, stac)
# } else {
if (ddep == 1 | ddep == 11) {
lx <- min(
1 + max(missnumspec, ceiling(K)),
DDD::roundn(pars2[1]) + missnumspec
)
} else {
lx <- DDD::roundn(pars2[1]) + missnumspec
}
if(loglik > -Inf)
{
# in all cases we integrate from the origin of the island to the colonization event
# (stac 2, 3, 4), the first branching point (stac = 6, 7), to the maximum colonization
# time (stac = 1, 5, 8, 9) or to the present (stac = 0)
probs <- rep(0,2 * lx + 1)
probs[1] <- 1 - probability_of_init_presence #Q^k_n
probs[lx + 1] <- probability_of_init_presence #Q^{M,k}_n
k1 <- 0
probs <- DAISIE_integrate_shift(probs,times = brts[1:2],tshift = tshift,loglik_rhs = DAISIE_loglik_rhs,pars1 = pars1,k1 = k1,ddep = ddep,rtol = reltolint,atol = abstolint,method = methode)
cp <- checkprobs2(lv = 2 * lx, loglik, probs, verbose); loglik <- cp[[1]]; probs = cp[[2]]
if(stac == 0)
{
# for stac = 0, the integration was from the origin of the island until
# the present so we can immediately evaluate the probability of no clade
# being present and no immigrant species.
loglik <- loglik + log(probs[1])
} else
{
if (stac %in% c(1, 5:9) )
{
# for stac = 1, we now integrate from the maximum colonization time
# (usually the island age + tiny time unit) until the present, where
# we set all probabilities where the immigrant is already present to 0
# and we evaluate the probability of the immigrant species being
# present, but there can be missing species.
# for stac = 5, we do exactly the same, but we evaluate the
# probability of an endemic species being present alone.
# for stac = 6 and 7, integration is from the maximum colonization
# time until the first branching time. This is the same as we did for
# stac = 1, 5.
# for stac = 8 and 9, integration is from the maximum colonization
# time until the minimum colonization time.
# In all cases we are dealing with a maximum colonization time which
# means that any colonization that took place before this maximum
# colonization time (including presence in the non-oceanic scenario)
# does not count and should be followed by another colonization.
# To allow this we introduce a third set of equations for the
# probability that colonization might have happened before but
# recolonization has not taken place yet (Q_M,n).
epss <- 1.01E-5 #We're taking the risk
if (abs(brts[2] - brts[1]) >= epss) {
probs[(2 * lx + 1):(4 * lx)] <- probs[1:(2 * lx)]
probs[1:(2 * lx)] <- 0
} else {
probs[(2 * lx + 1):(4 * lx)] <- 0
}
probs <- DAISIE_integrate_shift(probs,times = brts[2:3],tshift = tshift,loglik_rhs = DAISIE_loglik_rhs1,pars1 = pars1,k1 = k1,ddep = ddep,rtol = reltolint,atol = abstolint,method = methode)
cp <- checkprobs2(lx, loglik, probs, verbose); loglik <- cp[[1]]; probs <- cp[[2]]
if (stac %in% c(1, 5))
{
loglik <- loglik + log(probs[(stac == 1) * lx + (stac == 5) + 1 + missnumspec])
} else if (stac %in% c(6, 7, 8, 9))
{
probs2 <- rep(0, 3 * lx)
probs2[1:(lx - 1)] <- (1:(lx - 1)) * probs[2:lx]
probs2[(lx + 1):(2 * lx - 1)] <- (1:(lx - 1)) * probs[(lx + 2):(2 * lx)]
probs2[2 * lx + 1] <- probs[(lx + 1)]
probs2[(2 * lx + 2):(3 * lx)] <- 0
probs <- probs2
rm(probs2)
if (stac %in% c(8, 9))
{
k1 <- 1
probs <- DAISIE_integrate_shift(probs,times = brts[3:4],tshift = tshift,loglik_rhs = DAISIE_loglik_rhs2,pars1 = pars1,k1 = k1,ddep = ddep,rtol = reltolint,atol = abstolint,method = methode)
cp <- checkprobs2(lx, loglik, probs, verbose); loglik <- cp[[1]]; probs = cp[[2]]
loglik <- loglik + log(probs[(stac == 8) * (2 * lx + 1) + (stac == 9) + missnumspec])
}
}
} else if (stac %in% c(2, 3, 4) )
{
# for stac = 2, 3, 4, integration is then from the colonization
# event until the first branching time (stac = 2 and 3) or the present
# (stac = 4). We add a set of equations for Q_M,n, the probability
# that the process is compatible with the data, and speciation has not
# happened; during this time immigration is not allowed because it
# would alter the colonization time.
t <- brts[2]
gamvec <- divdepvec(
lac_or_gam = "gam",
pars1 = ifelse(rep(t,5) <= rep(tshift,5), pars1[1:5], pars1[6:10]),
t = t,
lx = lx,
k1 = k1,
ddep = ddep * (ddep == 11 | ddep == 21)
)
probs[(2 * lx + 1):(3 * lx)] = gamvec[1:lx] * probs[1:lx] +
gamvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)]
probs[1:(2 * lx)] = 0
k1 <- 1
probs <- DAISIE_integrate_shift(probs,times = brts[2:3],tshift = tshift,loglik_rhs = DAISIE_loglik_rhs2,pars1 = pars1,k1 = k1,ddep = ddep,rtol = reltolint,atol = abstolint,method = methode)
cp <- checkprobs2(lx,loglik,probs, verbose); loglik = cp[[1]]; probs = cp[[2]]
if (stac == 4)
# if stac = 4, we're done and we take an element from Q_M,n
{
loglik <- loglik + log(probs[2 * lx + 1 + missnumspec])
}
}
if (stac %in% c(2, 3, 6, 7) )
{
# at the first branching point all probabilities of states Q_M,n are
# transferred to probabilities where only endemics are present. Then
# go through the branching points.
S1 <- length(brts) - 1
startk <- 3
if(S1 >= startk)
{
t <- brts[startk]
lacvec <- divdepvec(
lac_or_gam = "lac",
pars1 = ifelse(rep(t,5) <= rep(tshift,5), pars1[1:5], pars1[6:10]),
t = t,
lx = lx + stac %in% c(6,7),
k1 = k1,
ddep = ddep
)
if(stac %in% c(2,3))
{
probs[1:lx] <- lacvec[1:lx] * (probs[1:lx] + probs[(2 * lx + 1):(3 * lx)])
probs[(lx + 1):(2 * lx)] <- lacvec[2:(lx + 1)] * probs[(lx + 1):(2 * lx)]
} else { # stac in c(6,7)
probs[1:(lx - 1)] <- lacvec[2:lx] *
((1:(lx - 1)) * probs[2:lx] + probs[(2 * lx + 1):(3 * lx - 1)])
probs[(lx + 1):(2 * lx - 1)] <- lacvec[3:(lx + 1)] * (1:(lx - 1)) *
probs[(lx + 2):(2 * lx)]
probs[lx] <- 0
probs[2 * lx] <- 0
}
probs <- probs[-c((2 * lx + 2):(3 * lx))]
probs[2 * lx + 1] <- 0
# After speciation, colonization is allowed again (re-immigration)
# all probabilities of states with the immigrant present are set to
# zero and all probabilities of states with endemics present are
# transported to the state with the colonist present waiting for
# speciation to happen. We also multiply by the (possibly diversity-
# dependent) immigration rate.
for (k in startk:S1)
{
k1 <- k - 1
probs <- DAISIE_integrate_shift(probs,times = brts[k:(k+1)],tshift = tshift,loglik_rhs = DAISIE_loglik_rhs,pars1 = pars1,k1 = k1,ddep = ddep,rtol = reltolint,atol = abstolint,method = methode)
cp <- checkprobs2(lx, loglik, probs, verbose); loglik <- cp[[1]]; probs = cp[[2]]
if(k < S1)
{
# speciation event
t <- brts[k + 1]
lacvec <- divdepvec(
lac_or_gam = "lac",
pars1 = ifelse(rep(t,5) <= rep(tshift,5), pars1[1:5], pars1[6:10]),
t = t,
lx = lx,
k1 = k1,
ddep = ddep
)
probs[1:(2 * lx)] <- c(lacvec[1:lx], lacvec[2:(lx + 1)]) *
probs[1:(2 * lx)]
}
}
}
# we evaluate the probability of the phylogeny with any missing species at the present without (stac = 2 or stac = 6) or with (stac = 3 or stac = 7) the immigrant species
loglik <- loglik + log(probs[(stac %in% c(3, 7)) * lx + 1 + missnumspec])
}
}
}
#}
if (length(pars1) == 11) { # CHANGE
print_parameters_and_loglik(pars = c(stac, pars1[5:10]), # should this be 6:10, or 6:11?
loglik = loglik,
verbose = pars2[4],
type = 'clade_loglik')
} else {
print_parameters_and_loglik(pars = c(stac, pars1[1:5]),
loglik = loglik,
verbose = pars2[4],
type = 'clade_loglik')
}
if (is.na(loglik)) {
message("NA in loglik encountered. Changing to -Inf.")
loglik <- -Inf
stop('stopping')
}
loglik <- as.numeric(loglik)
#testit::assert(is.numeric(loglik))
return(loglik)
}
DAISIE_integrate_shift <- function(probs,times,tshift,loglik_rhs,pars1,k1,ddep,rtol,atol,method) {
if(tshift > times[1] && tshift >= times[2]) {
probs <- DAISIE_integrate(probs,c(times[1:2]),loglik_rhs,c(pars1[1:5],k1,ddep),rtol = rtol,atol = atol,method = method)
} else
if(tshift > times[1] && tshift < times[2]) {
probs <- DAISIE_integrate(probs,c(times[1],tshift),loglik_rhs,c(pars1[1:5],k1,ddep),rtol = rtol,atol = atol,method = method)
probs <- DAISIE_integrate(probs,c(tshift,times[2]),loglik_rhs,c(pars1[6:10],k1,ddep),rtol = rtol,atol = atol,method = method)
} else
if(tshift <= times[1] && tshift < times[2]) {
probs <- DAISIE_integrate(probs,c(times[1:2]),loglik_rhs,c(pars1[6:10],k1,ddep),rtol = rtol,atol = atol,method = method)
}
return(probs)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.