dec2bin <- function(y, ly) {
stopifnot(length(y) == 1, mode(y) == "numeric")
q1 <- (y / 2) %/% 1
r <- y - q1 * 2
res <- c(r)
while (q1 >= 1) {
q2 <- (q1 / 2) %/% 1
r <- q1 - q2 * 2
q1 <- q2
res <- c(r, res)
}
res <- c(rep(0, ly - length(res)), res)
return(res)
}
dec2binmat <- function(y) {
numrows <- 2 ^ y
res <- matrix(0, numrows, y)
for (i in 0:(numrows - 1)) {
res[i + 1, ] <- dec2bin(i, y)
}
return(res)
}
bin2dec <- function(y) {
res <- y %*% 2 ^ ((length(y) - 1):0)
return(as.numeric(res))
}
kimat <- function(dec2binmatk) {
ki <- matrix(0, dim(dec2binmatk)[1], dim(dec2binmatk)[1])
for (i in 2:dim(dec2binmatk)[1]) {
locationones <- which(dec2binmatk[i, ] == 1)
for (j in 1:length(locationones)) {
dec2binmatki <- dec2binmatk[i, ]
dec2binmatki[locationones[j]] <- 0
j2 <- 1 + bin2dec(dec2binmatki)
ki[i,j2] <- 1
}
}
return(ki)
}
kmini <- function(dec2binmatk, lxm1, lxm2, lxe, sysdim = dim(dec2binmatk)[1]) {
posc <- Matrix::rowSums(dec2binmatk)
negc <- log2(sysdim) - posc
kmin <- rep(negc, each = lxm1 * lxm2 * lxe)
dim(kmin) <- c(lxm1, lxm2, lxe, sysdim)
ki <- kimat(dec2binmatk)
res <- list(kmin = kmin, ki = ki)
return(res)
}
nndivdep <- function(lxm1, lxm2, lxe, sysdim, Kprime, M, k, l) {
nnl <- c(0, 0:(lxm1 + 1))
nnm <- c(0, 0:(lxm2 + 1))
nne <- c(0, 0, 0:(lxe + 1))
lnnl <- length(nnl)
lnnm <- length(nnm)
lnne <- length(nne)
nil2lxm1 <- 2:(lxm1 + 1)
nil2lxm2 <- 2:(lxm2 + 1)
nil2lxe <- 3:(lxe + 2)
nn <- rowSums(expand.grid(n1 = nnl, n2 = nnm, n3 = nne))
dim(nn) <- c(lnnl, lnnm, lnne)
nn <- replicate(sysdim, nn)
nilm1 <- rep(1, lxm1)
nilm2 <- rep(1, lxm2)
nile <- rep(1, lxe)
allc <- 1:sysdim
divdepfac <- pmax(array(0, dim = c(lxm1 + 3, lxm2 + 3, lxe + 4, sysdim)),
1 - (nn + k) / Kprime)
divdepfacmin1 <- pmax(array(0, dim = c(lxm1 + 3, lxm2 + 3, lxe + 4, sysdim)),
1 - (nn + k - 1) / Kprime)
divdepfac <- divdepfac[nil2lxm1, nil2lxm2, nil2lxe, allc]
divdepfacmin1 <- divdepfacmin1[nil2lxm1, nil2lxm2, nil2lxe, allc]
Mminm <- M - nn[nil2lxm1, nil2lxm2, nile, allc]
lminm1 <- l - nn[nil2lxm1, nilm2, nile, allc]
Mminlminm2 <- M - l - nn[nilm1, nil2lxm2, nile, allc]
res <- list(nn = nn,
divdepfac = divdepfac,
divdepfacmin1 = divdepfacmin1,
Mminm = Mminm,
lminm1 = lminm1,
Mminlminm2 = Mminlminm2)
return(res)
}
selectrows <- function(sysdim, order) {
mat <- NULL
for (i in seq(1, sysdim / order, by = 2)) {
group1 <- (i - 1) * order + (1:order)
group2 <- i * order + (1:order)
mat <- rbind(mat, cbind(group1, group2))
}
return(mat)
}
DAISIE_loglik_rhs_IW = function(t,x,pars)
{
lac = pars[[1]][1]
mu = pars[[1]][2]
Kprime = pars[[1]][3]
gam = pars[[1]][4]
laa = pars[[1]][5]
M = pars[[1]][6]
k = pars[[2]]
ddep = pars[[3]]
lxm1 = pars[[4]]$lxm1
lxm2 = pars[[4]]$lxm2
lxe = pars[[4]]$lxe
sysdim = pars[[4]]$sysdim
kmin = pars[[5]]$kmin
kplus = k - kmin
ki = pars[[5]]$ki
nn = pars[[6]]$nn
divdepfac = pars[[6]]$divdepfac
divdepfacmin1 = pars[[6]]$divdepfacmin1
Mminm = pars[[6]]$Mminm
Mminm[Mminm < 0] = 0
lminm1minkminplus1 = pars[[6]]$lminm1 - kmin + 1
lminm1minkminplus1[lminm1minkminplus1 < 0] = 0
Mminlminm2plus1 = pars[[6]]$Mminlminm2 + 1
Mminlminm2plus1[Mminlminm2plus1 < 0] = 0
dim(x) = c(lxm1,lxm2,lxe,sysdim)
xx = array(0,dim = c(lxm1+2,lxm2+2,lxe+3,sysdim))
xx[2:(lxm1+1),2:(lxm2+1),3:(lxe+2),1:sysdim] = x
nil2lxm1 = 2:(lxm1 + 1)
nil2lxm2 = 2:(lxm2 + 1)
nil2lxe = 3:(lxe + 2)
allc = 1:sysdim
nilm1 = rep(1,lxm1)
nile = rep(1,lxe)
nilm2 = rep(1,lxm2)
if(sysdim == 1 & lxm1 > 1)
{
dim(Mminm) = c(lxm1,lxm2,lxe)
} else
if(sysdim == 1 & lxm1 == 1)
{
dim(Mminm) = c(lxm2,lxe)
} else
if (sysdim == 1 & lxm1 == 1) {
dim(Mminm) <- c(lxm2, lxe)
} else
if (sysdim > 1 & lxm1 == 1) {
dim(Mminm) <- c(lxm2, lxe, allc)
}
dx <- gam * divdepfacmin1 * lminm1minkminplus1 * xx[nil2lxm1 - 1, nil2lxm2, nil2lxe, allc] + #immigration
gam * divdepfacmin1 * Mminlminm2plus1 * xx[nil2lxm1, nil2lxm2 - 1, nil2lxe, allc] + #immigration
mu * nn[nil2lxm1 + 1, nilm2, nile, allc] * xx[nil2lxm1 + 1, nil2lxm2, nil2lxe, allc] + #extinction non-endemics
mu * nn[nilm1, nil2lxm2 + 1, nile, allc] * xx[nil2lxm1, nil2lxm2 + 1, nil2lxe, allc] + #extinction non-endemics
mu * nn[nilm1, nilm2, nil2lxe + 1, allc] * xx[nil2lxm1, nil2lxm2, nil2lxe + 1, allc] + #extinction endemics
lac * divdepfacmin1 * nn[nil2lxm1 + 1, nilm2, nile, allc] * xx[nil2lxm1 + 1, nil2lxm2, nil2lxe - 2, allc] + #cladogenesis non-endemics
lac * divdepfacmin1 * nn[nilm1, nil2lxm2 + 1, nile, allc] * xx[nil2lxm1, nil2lxm2 + 1, nil2lxe - 2, allc] + #cladogenesis non-endemics
lac * divdepfacmin1 * nn[nilm1, nilm2, nil2lxe - 1, allc] * xx[nil2lxm1, nil2lxm2, nil2lxe - 1, allc] + #cladogenesis endemics
2 * kplus * lac * divdepfacmin1 * xx[nil2lxm1, nil2lxm2, nil2lxe - 1, allc] + #cladogenesis species in tree
laa * nn[nil2lxm1 + 1, nilm2, nile, allc] * xx[nil2lxm1 + 1, nil2lxm2, nil2lxe - 1, allc] + #anagenesis non-endemics
laa * nn[nilm1, nil2lxm2 + 1, nile, allc] * xx[nil2lxm1, nil2lxm2 + 1, nil2lxe - 1, allc] + #anagenesis non-endemics
-(laa * (nn[nil2lxm1, nil2lxm2, nile, allc] + kmin) + (gam * divdepfac * Mminm) +
(lac * divdepfac + mu) * (nn[nil2lxm1, nil2lxm2, nil2lxe, allc] + k)) * xx[nil2lxm1, nil2lxm2, nil2lxe, allc]
if (sysdim > 1) {
dx <- dx +
laa * tensor::tensor(xx[nil2lxm1, nil2lxm2, nil2lxe, allc], ki, 4, 2) + # anagenesis in colonizing lineage
2 * lac * divdepfacmin1 * tensor::tensor(xx[nil2lxm1, nil2lxm2, nil2lxe - 1, allc], ki, 4, 2) # cladogenesis in colonizing lineage
}
dim(dx) <- c(sysdim * lxm1 * lxm2 * lxe, 1)
return(list(dx))
}
#' Computes the loglikelihood of the DAISIE model with island-wide
#' diversity-dependence given data and a set of model parameters
#'
#' Computes the loglikelihood of the DAISIE model given colonization and
#' branching times for lineages on an island, and a set of model parameters for
#' the DAISIE model with island-wide diversity-dependence
#'
#' The output is a loglikelihood value
#'
#' @param pars1 Contains the model parameters: \cr \cr
#' \code{pars1[1]} corresponds to lambda^c (cladogenesis rate) \cr
#' \code{pars1[2]} corresponds to mu (extinction rate) \cr
#' \code{pars1[3]} corresponds to K (clade-level carrying capacity) \cr
#' \code{pars1[4]} corresponds to gamma (immigration rate) \cr
#' \code{pars1[5]} corresponds to lambda^a (anagenesis rate) \cr
#' \code{pars1[6]} is optional; it may contain M, the total number of species
#' on the mainland \cr \cr
#' @param pars2 Contains the model settings \cr \cr
#' \code{pars2[1]} corresponds to lx = length of ODE variable x \cr
#' \code{pars2[2]} corresponds to ddmodel = diversity-dependent model, model of diversity-dependence, which can be one
#' of\cr \cr
#' ddmodel = 0 : no diversity dependence \cr
#' ddmodel = 1 : linear dependence in speciation rate \cr
#' ddmodel = 11: linear dependence in speciation rate and in immigration rate \cr
#' ddmodel = 2 : exponential dependence in speciation rate\cr
#' ddmodel = 21: exponential dependence in speciation rate and in immigration rate\cr
#' Only ddmodel = 11 is currently implemented \cr \cr
#' \code{pars2[3]} corresponds to cond = setting of conditioning\cr \cr
#' cond = 0 : conditioning on island age \cr
#' cond = 1 : conditioning on island age and non-extinction of the island biota \cr \cr
#' \code{pars2[4]} Specifies whether intermediate output should be provided,
#' because computation may take long. Default is 0, no output. A value of 1
#' means the parameters and loglikelihood are printed. A value of 2 means also
#' intermediate progress during loglikelihood computation is shown.
#' @param datalist Data object containing information on colonisation and
#' branching times. This object can be generated using the DAISIE_dataprep
#' function, which converts a user-specified data table into a data object, but
#' the object can of course also be entered directly. It is an R list object
#' with the following elements.\cr
#' The first element of the list has two or
#' three components: \cr \cr
#' \code{$island_age} - the island age \cr
#' Then, depending on whether a distinction between types is made, we have:\cr
#' \code{$not_present} - the number of mainland lineages that are not present
#' on the island \cr The remaining elements of the list each contains
#' information on a single colonist lineage on the island and has 5
#' components:\cr \cr
#' \code{$colonist_name} - the name of the species or clade
#' that colonized the island \cr
#' \code{$branching_times} - island age and stem
#' age of the population/species in the case of Non-endemic, Non-endemic_MaxAge
#' and Endemic anagenetic species. For cladogenetic species these should be
#' island age and branching times of the radiation including the stem age of
#' the radiation.\cr
#' \code{$stac} - the status of the colonist \cr \cr
#' * Non_endemic_MaxAge: 1 \cr
#' * Endemic: 2 \cr
#' * Endemic&Non_Endemic: 3 \cr
#' * Non_endemic: 4 \cr
#' * Endemic_MaxAge: 5 \cr \cr
#' \code{$missing_species} -
#' number of island species that were not sampled for particular clade (only
#' applicable for endemic clades) \cr
#' @param methode Method of the ODE-solver. See package deSolve for details.
#' Default is "ode45"
#' @param abstolint Absolute tolerance of the integration
#' @param reltolint Relative tolerance of the integration
#' @param verbose Logical controling if progress is printed to console.
#' @return The loglikelihood
#' @author Rampal S. Etienne & Bart Haegeman
#' @seealso \code{\link{DAISIE_ML_IW}}, \code{\link{DAISIE_loglik_CS}},
#' \code{\link{DAISIE_sim_constant_rate}}
#' @references Valente, L.M., A.B. Phillimore and R.S. Etienne (2015).
#' Equilibrium and non-equilibrium dynamics simultaneously operate in the
#' Galapagos islands. Ecology Letters 18: 844-852.
#' @keywords models
#' @export DAISIE_loglik_IW
DAISIE_loglik_IW <- function(
pars1,
pars2,
datalist,
methode = "ode45",
abstolint = 1E-16,
reltolint = 1E-14,
verbose = FALSE
)
{
if(is.na(pars2[4]))
{
pars2[4] = 0
}
if (is.null(datalist[[1]]$brts_table)) {
datalist <- Add_brt_table(datalist)
}
brts = c(-abs(datalist[[1]]$brts_table[,1]),0)
clade = datalist[[1]]$brts_table[,2]
event = datalist[[1]]$brts_table[,3]
pars1 = as.numeric(pars1)
if(length(pars1) == 5)
{
np = datalist[[1]]$not_present
if(is.null(np))
{
np = datalist[[1]]$not_present_type1 + datalist[[1]]$not_present_type2
}
if(is.null(np))
{
cat('Number of species not present is misspecified.\n')
loglik = NA
return(loglik)
}
M = length(datalist) - 1 + np
} else
if(length(pars1) == 6)
{
M = pars1[6]
} else
if (length(pars1) == 6) {
M <- pars1[6]
} else {
cat("pars1 should contain 5 or 6 values.\n")
loglik <- NA
return(loglik)
}
ddep = pars2[2]
cond = pars2[3]
lac = pars1[1]
mu = pars1[2]
Kprime = pars1[3]
if(ddep == 0)
{
Kprime <- Inf
}
gam <- pars1[4]
laa <- pars1[5]
l <- 0
if (min(pars1) < 0)
{
cat("One or more parameters are negative.\n")
loglik <- -Inf
return(loglik)
}
endemic = 0
nonendemic1 = 0
nonendemic2 = 0
for(i in 2:length(datalist))
{
endemic = endemic + (datalist[[i]]$stac == 5)
nonendemic1 = nonendemic1 + (datalist[[i]]$stac == 3)
nonendemic2 = nonendemic2 + (datalist[[i]]$stac == 1)
}
if((ddep == 1 | ddep == 11) & (ceiling(Kprime) < nonendemic1 + nonendemic2 + endemic + length(brts) - 2))
{
if (verbose) {
cat('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)
}
if (ddep == 1 | ddep == 11)
{
lx <- min(1 + ceiling(Kprime), DDD::roundn(pars2[1]) )
} else {
lx <- DDD::roundn(pars2[1])
}
lxm1 = max(clade) + 1
lxm2 = min(lx,M + 1)
lxe = lx
sysdimchange = 1
sysdim = 1
totdim = lxm1 * lxm2 * lxe * sysdim
probs = rep(0,totdim)
probs[1] = 1
loglik = 0
expandvec = NULL
for(k in 0:(length(brts) - 2))
{
if(pars2[4] == 2)
{
cat(paste('k = ',k,', sysdim = ',sysdim,'\n',sep = ''))
utils::flush.console()
}
dime <- list(lxm1 = lxm1, lxm2 = lxm2, lxe = lxe, sysdim = sysdim)
if (sysdimchange == 1) {
if (sysdim > 1) {
dec2binmatk <- dec2binmat(log2(sysdim))
kmi <- kmini(dec2binmatk, lxm1, lxm2, lxe, sysdim)
} else if (sysdim == 1) {
kmi <- list(kmin = 0, ki = NULL)
}
sysdimchange <- 0
}
nndd = nndivdep(lxm1,lxm2,lxe,sysdim,Kprime,M,k,l)
parslist = list(pars = pars1,k = k,ddep = ddep,dime = dime,kmi = kmi,nndd = nndd)
y = deSolve::ode(y = probs,times = brts[(k + 1):(k + 2)],func = DAISIE_loglik_rhs_IW,parms = parslist,rtol = reltolint,atol = abstolint,method = methode)
probs = y[2,2:(totdim + 1)]
cp = checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]]
dim(probs) = c(lxm1,lxm2,lxe,sysdim)
if(k < (length(brts) - 2))
{
divdepfac = nndd$divdepfac
if(event[k + 2] == 1)
{
Mminlminm2 = nndd$Mminlminm2
Mminlminm2[Mminlminm2 < 0] = 0
probs = gam * divdepfac * Mminlminm2 * probs[,,,1:sysdim]
probs = c(probs,rep(0,totdim))
l = l + 1
sysdim = sysdim * 2
expandvec = c(expandvec,clade[k + 2])
sysdimchange = 1
} else
{
probs = lac * divdepfac * probs[,,,1:sysdim]
if(event[k + 2] == 2)
{
tocollapse = which(expandvec == clade[k + 2])
sr = selectrows(sysdim,2^(tocollapse - 1))
probs = probs[,,,sr[,1]] + probs[,,,sr[,2]]
sysdim = sysdim / 2
dim(probs) = c(lxm1,lxm2,lxe,sysdim)
expandvec = expandvec[-tocollapse]
sysdimchange = 1
}
}
cp = checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]]
totdim = lxm1 * lxm2 * lxe * sysdim
dim(probs) = c(totdim,1)
}
}
dim(probs) <- c(lxm1, lxm2, lxe, sysdim)
expandedclades <- which(pracma::histc(clade, 1:length(clade))$cnt == 1)
status <- rep(0, lexpandedclades <- length(expandedclades))
if (lexpandedclades > 0) {
for (i in lexpandedclades:1) {
if (datalist[[1 + expandedclades[i]]]$stac == 2) {
status[i] <- 1
}
}
}
if (length(status) > 0) {
decstatus <- bin2dec(rev(status))
} else {
decstatus = 0
}
#print(status)
numcol = length(datalist) - 1
loglik = loglik - (lgamma(M + 1) - lgamma(M - numcol + 1) - lgamma(nonendemic2 + 1) - lgamma(endemic + 1))
#print(loglik + log(probs[1,,1,]))
loglik = loglik + log(probs[1 + nonendemic1,1 + nonendemic2,1 + endemic,1 + decstatus])
if(cond > 0)
{
sysdim = 1
totdim = lxm1 * lxm2 * lxe * sysdim
dime = list(lxm1 = lxm1,lxm2 = lxm2,lxe = lxe,sysdim = sysdim)
probs = rep(0,totdim)
probs[1] = 1
kmi = list(kmin = 0,ki = NULL)
nndd = nndivdep(lxm2,lxe,sysdim,Kprime,M,k = 0)
parslist = list(pars = pars1,k = k,ddep = ddep,dime = dime,kmi = kmi,nndd = nndd)
y = deSolve::ode(y = probs,times = brts[(k + 1):(k + 2)],func = DAISIE_loglik_rhs_IW,parms = parslist,rtol = reltolint,atol = abstolint,method = methode)
probs = y[2,2:(totdim + 1)]
dim(probs) = c(lxm1,lxm2,lxe,sysdim)
logcond = log(1 - probs[1,1,1,1])
loglik = loglik - logcond
}
if (pars2[4] >= 1) {
s1 <- sprintf("Parameters: %f %f %f %f %f", pars1[1], pars1[2], pars1[3], pars1[4], pars1[5])
s2 <- sprintf(", Loglikelihood: %f", loglik)
cat(s1, s2, "\n", sep = "")
utils::flush.console()
}
return(as.numeric(loglik))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.