######################################################################################################################################
######################################################################################################################################
### MiSSE -- Examines shifts in diversification in relation to hidden states ONLY
######################################################################################################################################
######################################################################################################################################
MiSSE <- function(phy, f=1, turnover=c(1,2), eps=c(1,2), fixed.eps=NULL, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL, includes.fossils=FALSE, k.samples=NULL, strat.intervals=NULL, sann=TRUE, sann.its=5000, sann.temp=5230, sann.seed=-100377, bounded.search=TRUE, max.tol=.Machine$double.eps^.50, starting.vals=NULL, turnover.upper=10000, eps.upper=3, trans.upper=100, restart.obj=NULL, ode.eps=0, dt.threads=1, expand.mode=FALSE){
misse_start_time <- Sys.time()
#This makes it easier to handle missegreedy with fixed values
if(length(fixed.eps)>0) {
if(is.na(fixed.eps)) {
fixed.eps <- NULL
}else{
eps <- numeric(length(turnover))
}
}
if(expand.mode) {
if(length(turnover)==1) {
turnover <- sequence(turnover)
}
if(length(eps)==1 & eps[1]!=0) {
eps <- sequence(eps)
}
if(length(eps)<length(turnover)) {
eps <- rep(eps, length(turnover))[sequence(length(turnover))] # i.e., c(1,2,1,2,1)
}
if(length(turnover)<length(eps)) {
turnover <- rep(turnover, length(eps))[sequence(length(eps))]
}
}
## Temporary fix for the current BUG:
if( !is.null(phy$node.label) ) phy$node.label <- NULL
if(!is.null(root.p)) {
root.p <- root.p / sum(root.p)
if(hidden.states ==TRUE){
if(length(root.p) < 26){
root.p.new <- numeric(26)
root.p.new[1:length(root.p)] <- root.p
root.p <- root.p.new
root.p <- root.p / sum(root.p)
}
}
}
setDTthreads(threads=dt.threads)
#if(!is.ultrametric(phy) & includes.fossils == FALSE){
# warning("Tree is not ultrametric. Used force.ultrametric() function to coerce the tree to be ultrametric - see note above.")
# edge_details <- GetEdgeDetails(phy, includes.intervals=FALSE, intervening.intervals=NULL)
# if(any(edge_details$type == "extinct_tip")){
# phy <- force.ultrametric(phy)
# }
#}
if(sann == FALSE & is.null(starting.vals)){
warning("You have chosen to rely on the internal starting points that generally work but does not guarantee finding the MLE.")
}
if(!root.type == "madfitz" & !root.type == "herr_als"){
stop("Check that you specified a proper root.type option. Options are 'madfitz' or 'herr_als'. See help for more details.", call.=FALSE)
}
if(length(turnover) != length(eps)){
stop("The number of turnover parameters need to match the number of extinction fraction parameters.", call.=FALSE)
}
ntips <- Ntip(phy)
param.count <- sum(c(length(unique(turnover)), length(unique(eps)), 1))
if(param.count > (ntips/10)){
warning("You might not have enough data to fit this model well", call.=FALSE, immediate.=TRUE)
}
pars <- numeric(53)
rate.cats <- hidden.states <- length(turnover)
turnover.tmp <- numeric(26)
turnover.tmp[1:length(turnover)] <- turnover
pars.tmp <- turnover
eps.tmp <- numeric(26)
eps.tmp[1:length(eps)] <- eps
eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, eps.tmp)
if(rate.cats > 1){
trans.tmp <- 1
trans.tmp <- trans.tmp + max(pars.tmp)
pars.tmp <- c(pars.tmp, trans.tmp)
}else{
trans.tmp <- 0
}
pars.tmp <- c(turnover.tmp[1], eps.tmp[1], turnover.tmp[2], eps.tmp[2], turnover.tmp[3], eps.tmp[3], turnover.tmp[4], eps.tmp[4], turnover.tmp[5], eps.tmp[5], turnover.tmp[6], eps.tmp[6], turnover.tmp[7], eps.tmp[7], turnover.tmp[8], eps.tmp[8], turnover.tmp[9], eps.tmp[9], turnover.tmp[10], eps.tmp[10], turnover.tmp[11], eps.tmp[11], turnover.tmp[12], eps.tmp[12], turnover.tmp[13], eps.tmp[13], turnover.tmp[14], eps.tmp[14], turnover.tmp[15], eps.tmp[15], turnover.tmp[16], eps.tmp[16], turnover.tmp[17], eps.tmp[17], turnover.tmp[18], eps.tmp[18], turnover.tmp[19], eps.tmp[19], turnover.tmp[20], eps.tmp[20], turnover.tmp[21], eps.tmp[21], turnover.tmp[22], eps.tmp[22], turnover.tmp[23], eps.tmp[23], turnover.tmp[24], eps.tmp[24], turnover.tmp[25], eps.tmp[25], turnover.tmp[26], eps.tmp[26], trans.tmp[1])
pars[1:length(pars.tmp)] <- pars.tmp
if(includes.fossils == TRUE){
pars <- c(pars, max(pars.tmp)+1)
}else{
pars <- c(pars, 0)
}
np <- max(pars)
pars[pars==0] <- np + 1
cat("Initializing...", "\n")
# Some new prerequisites #
if(includes.fossils == TRUE){
if(!is.null(k.samples)){
phy.og <- phy
psi.type <- "m+k"
split.times <- dateNodes(phy, rootAge=max(node.depth.edgelength(phy)))[-c(1:Ntip(phy))]
strat.cache <- NULL
k.samples <- k.samples[order(as.numeric(k.samples[,3]), decreasing=FALSE),]
phy <- AddKNodes(phy, k.samples)
fix.type <- GetKSampleMRCA(phy, k.samples)
no.k.samples <- length(k.samples[,1])
gen <- FindGenerations(phy)
dat.tab <- OrganizeDataMiSSE(phy=phy, f=f, hidden.states=hidden.states, includes.intervals=FALSE, intervening.intervals=NULL, includes.fossils=includes.fossils)
#These are all inputs for generating starting values:
edge_details <- GetEdgeDetails(phy, includes.intervals=FALSE, intervening.intervals=NULL)
fossil.taxa <- edge_details$tipward_node[which(edge_details$type == "extinct_tip")]
fossil.ages <- dat.tab$TipwardAge[which(dat.tab$DesNode %in% fossil.taxa)]
n.tax.starting <- Ntip(phy)-length(fossil.taxa)-no.k.samples
}else{
if(!is.null(strat.intervals)){
phy.og <- phy
psi.type <- "m+int"
split.times.plus.tips <- dateNodes(phy, rootAge=max(node.depth.edgelength(phy)))
split.times <- split.times.plus.tips[-c(1:Ntip(phy))]
strat.cache <- GetStratInfo(strat.intervals=strat.intervals)
k.samples <- GetIntervalToK(strat.intervals, intervening.intervals=strat.cache$intervening.intervals)
extinct.tips <- which(round(k.samples$timefrompresent,8) %in% round(split.times.plus.tips[c(1:Ntip(phy))],8))
if(length(extinct.tips > 0)){
k.samples <- k.samples[-extinct.tips,]
}
phy <- AddKNodes(phy, k.samples)
fix.type <- GetKSampleMRCA(phy, k.samples, strat.intervals=TRUE)
gen <- FindGenerations(phy)
dat.tab <- OrganizeDataMiSSE(phy=phy, f=f, hidden.states=hidden.states, includes.intervals=TRUE, intervening.intervals=strat.cache$intervening.intervals, includes.fossils=includes.fossils)
#These are all inputs for generating starting values:
edge_details <- GetEdgeDetails(phy, includes.intervals=TRUE, intervening.intervals=strat.cache$intervening.intervals)
fossil.taxa <- edge_details$tipward_node[which(edge_details$type == "extinct_tip" | edge_details$type == "k_extinct_interval")]
fossil.ages <- dat.tab$TipwardAge[which(dat.tab$DesNode %in% fossil.taxa)]
n.tax.starting <- Ntip(phy)-length(fossil.taxa)-dim(fix.type)[1]
}else{
phy.og <- phy
psi.type <- "m_only"
fix.type <- NULL
split.times <- dateNodes(phy, rootAge=max(node.depth.edgelength(phy)))[-c(1:Ntip(phy))]
strat.cache <- NULL
no.k.samples <- 0
gen <- FindGenerations(phy)
dat.tab <- OrganizeDataMiSSE(phy=phy, f=f, hidden.states=hidden.states, includes.intervals=FALSE, intervening.intervals=NULL, includes.fossils=includes.fossils)
#These are all inputs for generating starting values:
edge_details <- GetEdgeDetails(phy, includes.intervals=FALSE, intervening.intervals=NULL)
fossil.taxa <- edge_details$tipward_node[which(edge_details$type == "extinct_tip")]
fossil.ages <- dat.tab$TipwardAge[which(dat.tab$DesNode %in% fossil.taxa)]
n.tax.starting <- Ntip(phy)-length(fossil.taxa)-no.k.samples
}
}
}else{
phy.og <- phy
gen <- FindGenerations(phy)
dat.tab <- OrganizeDataMiSSE(phy=phy, f=f, hidden.states=hidden.states, includes.intervals=FALSE, intervening.intervals=NULL, includes.fossils=includes.fossils)
fossil.taxa <- NULL
fix.type <- NULL
psi.type <- NULL
strat.cache <- NULL
}
nb.tip <- Ntip(phy)
nb.node <- phy$Nnode
##########################
#This is used to scale starting values to account for sampling:
if(length(f) == 1){
samp.freq.tree <- f
}else{
if(length(f) == Ntip(phy)){
stop("This functionality has been temporarily removed.")
}else{
stop("The vector of sampling frequencies does not match the number of tips in the tree.")
}
}
if(is.null(restart.obj)){
if(sum(eps)==0){
if(includes.fossils == TRUE){
stop("You input a tree that contains fossils but you are assuming no extinction. Check your function call.")
}else{
init.pars <- starting.point.generator(phy, k=2, samp.freq.tree, yule=TRUE)
}
}else{
if(includes.fossils == TRUE){
if(!is.null(strat.intervals)){
branch.type = NULL
cols <- c("FocalNode","DesNode", "RootwardAge", "TipwardAge", "branch.type")
seg.map <- dat.tab[, cols, with=FALSE]
#remove k tips -- we do not do anything with them.
setkey(seg.map, branch.type)
#drop the k.tips because we do not do calculation on these zero length edges:
seg.map <- seg.map[branch.type != 2]
init.pars <- starting.point.generator.intervals(k=2, samp.freq.tree, n.tax=n.tax.starting, seg_map=seg.map, split.times=split.times, fossil.ages=fossil.ages, strat.cache=strat.cache)
}else{
init.pars <- starting.point.generator.fossils(n.tax=n.tax.starting, k=2, samp.freq.tree, fossil.taxa=fossil.taxa, fossil.ages=fossil.ages, no.k.samples=no.k.samples, split.times=split.times)
}
psi.start <- init.pars[length(init.pars)]
}else{
init.pars <- starting.point.generator(phy, k=2, samp.freq.tree, yule=FALSE)
}
if(init.pars[3] == 0){
init.pars[3] = 1e-6
}
}
names(init.pars) <- NULL
if(is.null(starting.vals)){
def.set.pars <- rep(NA, 2*rate.cats)
scaling=seq(from=1.2, to=0.8, length.out=rate.cats)
for(cat.index in sequence(rate.cats)) {
def.set.pars[1+(cat.index-1)*2] <- log((init.pars[1]+init.pars[3]) * ifelse(length(unique(turnover))==1, 1, scaling[cat.index]))
def.set.pars[2+(cat.index-1)*2] <- log(ifelse(length(unique(eps))==1, 1, scaling[cat.index]) * init.pars[3]/init.pars[1])
}
#def.set.pars <- rep(c(log(init.pars[1]+init.pars[2]), log(init.pars[2]/init.pars[1])), rate.cats)
#trans.start <- log(rate.cats/sum(phy$edge.length))
trans.start <- log(init.pars[5])
}else{
if(includes.fossils == TRUE){
## Check the length of the stating vector:
if( length( starting.vals ) != 4 ){
stop("Incorrect length for starting.vals vector.")
}
def.set.pars <- rep(c(log(starting.vals[1]), log(starting.vals[2])), rate.cats)
trans.start <- log(starting.vals[3])
psi.start <- log(starting.vals[4])
}else{
## Check the length of the stating vector:
if( length( starting.vals ) != 3 ){
stop("Incorrect length for starting.vals vector.")
}
def.set.pars <- rep(c(log(starting.vals[1]), log(starting.vals[2])), rate.cats)
trans.start <- log(starting.vals[3])
}
}
if(bounded.search == TRUE){
upper.full <- rep(c(log(turnover.upper), log(eps.upper)), rate.cats)
}else{
upper.full <- rep(21, length(def.set.pars))
}
if(rate.cats > 1){
if(includes.fossils == TRUE){
np.sequence <- 1:(np-2)
ip <- numeric(np-2)
upper <- numeric(np-2)
}else{
np.sequence <- 1:(np-1)
ip <- numeric(np-1)
upper <- numeric(np-1)
}
}else{
if(includes.fossils == TRUE){
np.sequence <- 1:(np-1)
ip <- numeric(np-1)
upper <- numeric(np-1)
}else{
np.sequence <- 1:np
ip <- numeric(np)
upper <- numeric(np)
}
}
for(i in np.sequence){
ip[i] <- def.set.pars[which(pars == np.sequence[i])[1]]
upper[i] <- upper.full[which(pars == np.sequence[i])[1]]
}
if(rate.cats > 1){
if(includes.fossils == TRUE){
ip <- c(ip, trans.start, log(psi.start))
upper <- c(upper, log(trans.upper), log(trans.upper))
}else{
ip <- c(ip, trans.start)
upper <- c(upper, log(trans.upper))
}
}else{
if(includes.fossils == TRUE){
ip <- c(ip, log(psi.start))
upper <- c(upper, log(trans.upper))
}
}
lower <- rep(-20, length(ip))
}else{
upper <- restart.obj$upper.bounds
lower <- restart.obj$lower.bounds
pars <- restart.obj$index.par
ip <- numeric(length(unique(restart.obj$index.par))-1)
for(k in 1:length(ip)){
ip[k] <- log(restart.obj$solution[which(restart.obj$index.par==k)][1])
}
}
if(sann == FALSE){
if(bounded.search == TRUE){
cat("Finished. Beginning bounded subplex routine...", "\n")
opts <- list("algorithm" = "NLOPT_LN_SBPLX", "maxeval" = 100000, "ftol_rel" = max.tol)
out = nloptr(x0=ip, eval_f=DevOptimizeMiSSE, ub=upper, lb=lower, opts=opts, pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, fixed.eps=fixed.eps, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type, strat.cache=strat.cache)
sann.counts <- NULL
solution <- numeric(length(pars))
solution[] <- c(exp(out$solution), 0)[pars]
loglik = -out$objective
}else{
cat("Finished. Beginning subplex routine...", "\n")
out = subplex(ip, fn=DevOptimizeMiSSE, control=list(reltol=max.tol, parscale=rep(0.1, length(ip))), pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, fixed.eps=fixed.eps, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type, strat.cache=strat.cache)
sann.counts <- NULL
solution <- numeric(length(pars))
solution[] <- c(exp(out$par), 0)[pars]
loglik = -out$value
}
}else{
cat("Finished. Beginning simulated annealing...", "\n")
out.sann = GenSA(ip, fn=DevOptimizeMiSSE, lower=log(exp(lower)+0.0000000001), upper=log(exp(upper)-0.0000000001), control=list(max.call=sann.its, temperature=sann.temp, seed=sann.seed), pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, fixed.eps=fixed.eps, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type, strat.cache=strat.cache)
sann.counts <- out.sann$counts
cat("Finished. Refining using subplex routine...", "\n")
opts <- list("algorithm" = "NLOPT_LN_SBPLX", "maxeval" = 100000, "ftol_rel" = max.tol)
out <- nloptr(x0=out.sann$par, eval_f=DevOptimizeMiSSE, ub=upper, lb=lower, opts=opts, pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, fixed.eps=fixed.eps, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type, strat.cache=strat.cache)
solution <- numeric(length(pars))
solution[] <- c(exp(out$solution), 0)[pars]
loglik = -out$objective
}
names(solution) <- c("turnover0A","eps0A", "turnover0B","eps0B", "turnover0C","eps0C", "turnover0D","eps0D", "turnover0E","eps0E", "turnover0F","eps0F", "turnover0G","eps0G", "turnover0H","eps0H", "turnover0I","eps0I", "turnover0J","eps0J", "turnover0K","eps0K", "turnover0L","eps0L", "turnover0M","eps0M", "turnover0N","eps0N", "turnover0O","eps0O", "turnover0P","eps0P", "turnover0Q","eps0Q", "turnover0R","eps0R", "turnover0S","eps0S", "turnover0T","eps0T", "turnover0U","eps0U", "turnover0V","eps0V","turnover0W","eps0W","turnover0X","eps0X", "turnover0Y","eps0Y", "turnover0Z","eps0Z", "q0", "psi")
cat("Finished. Summarizing results...", "\n")
misse_end_time <- Sys.time()
obj = list(loglik = loglik, AIC = -2*loglik+2*np, AICc = -2*loglik+(2*np*(Ntip(phy)/(Ntip(phy)-np-1))), solution=solution, index.par=pars, f=f, hidden.states=hidden.states, fixed.eps=fixed.eps, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, phy=phy.og, phy.w.k=phy, max.tol=max.tol, starting.vals=ip, upper.bounds=upper, lower.bounds=lower, ode.eps=ode.eps, turnover=turnover, eps=eps, elapsed.minutes=difftime(misse_end_time, misse_start_time, units="min"), includes.fossils=includes.fossils, k.samples=k.samples, strat.intervals=strat.intervals, fix.type=fix.type, psi.type=psi.type, sann.counts=sann.counts)
class(obj) <- append(class(obj), "misse.fit")
return(obj)
}
######################################################################################################################################
######################################################################################################################################
### MiSSEGreedy -- Automated algorithm for running MiSSE of varying complexity
######################################################################################################################################
######################################################################################################################################
# options(error = utils::recover)
# a <- hisse:::MiSSEGreedyNew(ape::rcoal(50), possible.combos=hisse:::generateMiSSEGreedyCombinations(4), n.cores=4, save.file='~/Downloads/greedy.rda')
MiSSEGreedy <- function(phy, f=1, possible.combos = generateMiSSEGreedyCombinations(shuffle.start=TRUE), stop.deltaAICc=10, save.file=NULL, n.cores=NULL, chunk.size=10, check.fits=FALSE, remove.bad=FALSE, n.tries=2, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL, includes.fossils=FALSE, k.samples=NULL, strat.intervals=NULL, sann=TRUE, sann.its=5000, sann.temp=5230, sann.seed=-100377, bounded.search=TRUE, max.tol=.Machine$double.eps^.50, starting.vals=NULL, turnover.upper=10000, eps.upper=3, trans.upper=100, restart.obj=NULL, ode.eps=0) {
misse.list <- list()
chunk.size <- ifelse(is.null(chunk.size),ifelse(is.null(n.cores),1,n.cores), chunk.size)
total.chunks <- ceiling(nrow(possible.combos)/chunk.size)
possible.combos$runorder <- NA
possible.combos$deltaAICc <- NA
possible.combos$AICc <- NA
possible.combos$predictedAICc <- NA
possible.combos$lnL <- NA
possible.combos$AIC <- NA
possible.combos$elapsedMinutes <- NA
possible.combos$predictedMinutes <- NA
final.combos <- possible.combos
all.examined <- c()
for (batch_index in sequence(total.chunks)) { # So, if we can do parallel, we do it in chunks so all cores are busy
starting.time <- Sys.time()
focal.models <- sequence(min(chunk.size, nrow(possible.combos)))
if(batch_index>1) {
#save(final.combos, possible.combos, file="stufffordebugging.Rsave")
#final.combos$predictedAICc <- stats::predict(stats::glm(AICc ~ turnover + eps + turnover*eps, data=subset(final.combos, !is.na(final.combos$AICc))), newdata=final.combos) #Idea here is to focus on the best candidate models
model.distances <- as.matrix(dist(final.combos[,c("turnover", "eps")], diag=TRUE, upper=TRUE))
unknown.indices <- which(is.na(final.combos$AICc))
for(unknown.index in unknown.indices) {
model.distances[,unknown.index] <- 1e10
}
for (model_index in sequence(nrow(final.combos))) {
best <- which(model.distances[model_index,] == min(model.distances[model_index,]))
final.combos$predictedAICc[model_index] <- min(final.combos[best, "AICc"]+model.distances[model_index, best], na.rm=TRUE)
}
best.ones <- base::order(final.combos$predictedAICc)
best.ones <- best.ones[!(best.ones %in% which(!is.na(final.combos$AICc)))]
focal.models <- best.ones[1:min(chunk.size, length(best.ones))]
}
all.examined <- append(all.examined, focal.models)
#local.combos <- possible.combos[(1+chunk.size*(batch_index-1)):min(nrow(possible.combos), chunk.size*batch_index) ,]
local.combos <- final.combos[focal.models,]
#cat("\nNow starting run with", paste(range(local.combos$turnover), collapse="-"), "turnover categories and", paste(range(local.combos$eps), collapse="-"), "extinction fraction categories", "\n")
cat("Starting at ", as.character(starting.time), "\n running on ", n.cores, " cores.", sep="")
cat("\n")
print(local.combos[,1:3])
cat("\n")
misse.list <- append(misse.list, parallel::mcmapply(
MiSSE,
eps=local.combos$eps,
turnover=local.combos$turnover,
fixed.eps=local.combos$fixed.eps,
MoreArgs=list(
phy=phy,
f=f,
condition.on.survival=condition.on.survival,
root.type=root.type,
root.p=root.p,
includes.fossils=includes.fossils,
k.samples=k.samples,
strat.intervals=strat.intervals,
sann=sann,
sann.its=sann.its,
sann.temp=sann.temp,
sann.seed=sann.seed,
bounded.search=bounded.search,
max.tol=max.tol,
starting.vals=starting.vals,
turnover.upper=turnover.upper,
eps.upper=eps.upper,
trans.upper=trans.upper,
restart.obj=restart.obj,
ode.eps=ode.eps,
expand.mode=TRUE
),
mc.cores=ifelse(is.null(n.cores),1,n.cores),
SIMPLIFY=FALSE
))
AICc <- unlist(lapply(misse.list, "[[", "AICc"))
deltaAICc <- AICc-min(AICc)
min.deltaAICc.this.chunk <- min(deltaAICc[(1+chunk.size*(batch_index-1)):min(nrow(final.combos), chunk.size*batch_index)])
# final.combos$lnL[1:min(nrow(final.combos), chunk.size*batch_index)] <- unlist(lapply(misse.list, "[[", "loglik"))
# final.combos$AIC[1:min(nrow(final.combos), chunk.size*batch_index)] <- unlist(lapply(misse.list, "[[", "AIC"))
# final.combos$AICc[1:min(nrow(final.combos), chunk.size*batch_index)] <- AICc
# final.combos$deltaAICc[1:min(nrow(final.combos), chunk.size*batch_index)] <- deltaAICc
# final.combos$elapsedMinutes[1:min(nrow(final.combos), chunk.size*batch_index)] <- unlist(lapply(misse.list, "[[", "elapsed.minutes"))
# data.for.fit <- data.frame(nparam=(final.combos$eps+possible.combos$turnover)[1:min(nrow(final.combos), chunk.size*batch_index)], logmin=log(final.combos$elapsedMinutes[1:min(nrow(final.combos), chunk.size*batch_index)]))
# data.for.prediction <- data.frame(nparam=(final.combos$eps+final.combos$turnover))
final.combos$lnL[all.examined] <- unlist(lapply(misse.list, "[[", "loglik"))
final.combos$AIC[all.examined] <- unlist(lapply(misse.list, "[[", "AIC"))
final.combos$AICc[all.examined] <- AICc
final.combos$deltaAICc[all.examined] <- deltaAICc
final.combos$elapsedMinutes[all.examined] <- unlist(lapply(misse.list, "[[", "elapsed.minutes"))
final.combos$runorder[focal.models] <- batch_index
data.for.fit <- data.frame(nparam=(final.combos$eps+final.combos$turnover)[all.examined], logmin=log(final.combos$elapsedMinutes[all.examined]))
data.for.prediction <- data.frame(nparam=(final.combos$eps+final.combos$turnover))
suppressWarnings(final.combos$predictedMinutes <- exp(predict(lm(logmin ~ nparam, data=data.for.fit), newdata=data.for.prediction)))
cat("\nResults so far\n")
final.combos[] <- lapply(final.combos, function(x) { if(!is.numeric(x)) as.numeric(x) else x })
print(round(final.combos,2))
if(!is.null(save.file)) {
save(misse.list, final.combos, file=save.file)
}
if(batch_index<total.chunks) {
if(stop.deltaAICc>min.deltaAICc.this.chunk) {
print(paste0("Best AICc in this set of parallel runs was ", round(min.deltaAICc.this.chunk,2), " which is less than the cutoff to stop running (",stop.deltaAICc,"), so starting another set of parallel runs"))
} else {
print(paste0("Best AICc in this set of parallel runs was ", round(min.deltaAICc.this.chunk,2), " which is greater than the cutoff to stop running (",stop.deltaAICc,"), so stopping here"))
break()
}
}
# print("\n")
# print(local.combos)
# misse.list <- append(misse.list, MiSSE(
# eps=local.combos$eps[1],
# turnover=local.combos$turnover[1],
# fixed.eps=local.combos$fixed.eps[1],
#
# phy=phy,
# f=f,
# condition.on.survival=condition.on.survival,
# root.type=root.type,
# root.p=root.p,
# sann=sann,
# sann.its=sann.its,
# bounded.search=bounded.search,
# max.tol=max.tol,
# starting.vals=starting.vals,
# turnover.upper=turnover.upper,
# eps.upper=eps.upper,
# trans.upper=trans.upper,
# restart.obj=restart.obj,
# ode.eps=ode.eps,
# expand.mode=TRUE
# ))
}
if(check.fits == TRUE){
cat("Checking model fits...", "\n")
misse.list.updated <- MiSSENet(misse.list=misse.list, n.tries=n.tries, remove.bad=remove.bad, dont.rerun=FALSE, save.file=save.file, n.cores=ifelse(is.null(n.cores),1,n.cores), sann=sann, sann.its=sann.its, sann.temp=sann.temp, bounded.search=bounded.search, starting.vals=starting.vals, turnover.upper=turnover.upper, eps.upper=eps.upper, trans.upper=trans.upper, restart.obj=restart.obj)
return(misse.list.updated)
}else{
if(remove.bad == TRUE){
misse.list.updated <- MiSSENet(misse.list=misse.list, n.tries=n.tries, remove.bad=remove.bad, dont.rerun=TRUE, save.file=save.file, n.cores=ifelse(is.null(n.cores),1,n.cores), sann=sann, sann.its=sann.its, sann.temp=sann.temp, bounded.search=bounded.search, starting.vals=starting.vals, turnover.upper=turnover.upper, eps.upper=eps.upper, trans.upper=trans.upper, restart.obj=restart.obj)
return(misse.list.updated)
}else{
return(misse.list)
}
}
}
SummarizeMiSSEGreedy <- function(greedy.result, min.weight=0.01, n.cores=1, recon=TRUE) {
parameters <- c("turnover", "net.div", "speciation", "extinction", "extinction.fraction")
AICc_weights <- GetAICWeights(greedy.result, criterion="AICc")
good_enough <- which(AICc_weights>min.weight)
best_model <- which.max(AICc_weights)
rates <- NA
if(recon) {
recons <- list()
rates <- array(dim=c(ape::Ntip(greedy.result[[1]]$phy) + ape::Nnode(greedy.result[[1]]$phy), length(parameters), length(good_enough)+2), dimnames=list(c(greedy.result[[1]]$phy$tip.label, (1+ape::Ntip(greedy.result[[1]]$phy)):(ape::Ntip(greedy.result[[1]]$phy) + ape::Nnode(greedy.result[[1]]$phy))), parameters, c(good_enough, "model_average", "best")))
for(i in sequence(length(good_enough))){
local_model <- greedy.result[[good_enough[i]]]
recons[[i]] <- MarginReconMiSSE(phy=local_model$phy, f=local_model$f, pars=local_model$solution[local_model$index.par], hidden.states=local_model$hidden.states, fixed.eps=local_model$fixed.eps, condition.on.survival=local_model$condition.on.survival, root.type=local_model$root.type, root.p=local_model$root.p, includes.fossils=local_model$includes.fossils, k.samples=local_model$k.samples, strat.intervals=local_model$strat.intervals, AIC=local_model$AICc, get.tips.only=FALSE, verbose=FALSE, n.cores=n.cores)
}
for(param_index in sequence(length(parameters))) {
rate.param <- parameters[param_index]
for(i in sequence(length(good_enough))){
rates.tips <- ConvertManyToRate(recons[i], rate.param, "tip.mat")
rates.internal <- ConvertManyToRate(recons[i], rate.param, "node.mat")
rates[,param_index,i] <- c(rates.tips, rates.internal)
if(good_enough[i]==best_model) {
rates[,param_index,2+length(good_enough)] <- c(rates.tips, rates.internal)
}
}
rates.tips.avg <- ConvertManyToRate(recons, rate.param, "tip.mat")
rates.internal.avg <- ConvertManyToRate(recons, rate.param, "node.mat")
rates[,param_index,1+length(good_enough)] <- c(rates.tips.avg, rates.internal.avg)
}
}
overview <- data.frame(model_number=sequence(length(greedy.result)))
overview$lnL <- unlist(lapply(greedy.result, "[[", "loglik"))
overview$AICc <- unlist(lapply(greedy.result, "[[", "AICc"))
overview$deltaAICc <- overview$AICc-min(overview$AICc)
overview$AIC <- unlist(lapply(greedy.result, "[[", "AIC"))
overview$hidden.states <- unlist(lapply(greedy.result, "[[", "hidden.states"))
overview$n.turnover <- NA
overview$n.eps <- NA
overview$n.transition_rate <- 1
overview$n.freeparam <- NA
for(i in sequence(length(greedy.result))) {
overview$n.turnover[i] <- max(greedy.result[[i]]$turnover)
overview$n.eps[i] <- max(greedy.result[[i]]$eps)
overview$n.freeparam[i] <- overview$n.turnover[i]+overview$n.eps[i]+1
}
overview$AICc_weights <- AICc_weights
overview$Included_In_Model_Average <- FALSE
overview$Included_In_Model_Average[good_enough] <- TRUE
overview$elapsed.minutes <- unlist(lapply(greedy.result, "[[", "elapsed.minutes"))
return(list(overview=overview, rates=rates))
}
generateMiSSEGreedyCombinations <- function(max.param=52, turnover.tries=sequence(26), eps.tries=sequence(26), fixed.eps.tries=NA, vary.both=TRUE, shuffle.start=TRUE) {
fixed.eps <- eps <- "CRAN wants this declared somehow"
if(vary.both) {
combos <- expand.grid(turnover=turnover.tries, eps=eps.tries, fixed.eps=fixed.eps.tries)
} else {
combos <- rbind(
expand.grid(turnover=turnover.tries, eps=1, fixed.eps=fixed.eps.tries),
expand.grid(turnover=1, eps=eps.tries, fixed.eps=fixed.eps.tries)
)
}
combos$eps[which(!is.na(combos$fixed.eps))] <- 0
rownames(combos) <- NULL
combos <- subset(combos, eps==0 | is.na(fixed.eps)) # Don't estimate multiple eps while also fixing eps
combos <- combos[!duplicated(combos),]
combos <- combos[which(combos$turnover + combos$eps <= max.param),]
combos <- combos[base::order(combos$eps, combos$turnover, decreasing=FALSE),]
if(shuffle.start) {
model_pref <- 1/(combos$turnover + 2*combos$eps)^4 #reorder with weight on simplest ones first
combos <- combos[sample.int(n=nrow(combos), prob=model_pref, replace=FALSE), ]
}
rownames(combos) <- NULL
return(combos)
}
#Original version -- now defunct but kept for posterity.
MiSSEGreedyOLD <- function(phy, f=1, turnover.tries=sequence(26), eps.constant=c(TRUE,FALSE), stop.count=2, stop.deltaAICc=10, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL, sann=FALSE, sann.its=10000, bounded.search=TRUE, max.tol=.Machine$double.eps^.50, starting.vals=NULL, turnover.upper=10000, eps.upper=3, trans.upper=100, restart.obj=NULL, ode.eps=0, n.cores=NULL) {
misse.list <- list()
first.AICc <- Inf
for (eps.index in seq_along(eps.constant)) {
best.AICc <- first.AICc #reset so we start over at one turnover parameter and work our way back up
times.since.close.enough <- 0
for (turnover.index in seq_along(turnover.tries)) {
if(turnover.index>1 | eps.index==1) { # don't do the 1 turnover, 1 eps model twice -- overcounts it
starting.time <- Sys.time()
turnover <- sequence(turnover.tries[turnover.index])
eps <- turnover
if(eps.constant[eps.index]) {
eps <- rep(1, length(turnover))
}
cat("\nNow starting run with", turnover.tries[turnover.index], "turnover categories and", length(unique(eps)), "extinction fraction categories", "\n")
cat("Starting at ", as.character(starting.time), "\n", sep="")
current.run <- MiSSE(phy, f=f, turnover=turnover, eps=eps, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, sann=sann, sann.its=sann.its, bounded.search=bounded.search, max.tol=max.tol, starting.vals=starting.vals, turnover.upper=turnover.upper, eps.upper=eps.upper, trans.upper=trans.upper, restart.obj=restart.obj, ode.eps=ode.eps)
misse.list[[length(misse.list)+1]] <- current.run
ending.time <- Sys.time()
cat("Finished at ", as.character(ending.time),"; this model took ", round(difftime(ending.time,starting.time, units="mins"),2), " minutes\n",sep="")
cat("Current AICc is ", current.run$AICc, "\n", sep="")
if(is.infinite(first.AICc)) {
first.AICc <- current.run$AICc # so when we reset, use the 1,1 AICc, not Inf
}
if(current.run$AICc < best.AICc) {
cat("Found better AICc by ", best.AICc - current.run$AICc, "\n", sep="")
best.AICc <- current.run$AICc
times.since.close.enough <- 0
} else if ((current.run$AICc - best.AICc ) < stop.deltaAICc) {
cat("Found worse AICc by ", current.run$AICc - best.AICc , ", but this is still within ", stop.deltaAICc, " of the best", "\n", sep="")
times.since.close.enough <- 0
} else {
times.since.close.enough <- times.since.close.enough + 1
cat("Found worse AICc by ", current.run$AICc - best.AICc , ", it has been ", times.since.close.enough, " models since finding one within ", stop.deltaAICc, " of the best", "\n", sep="")
if(times.since.close.enough > stop.count) {
break()
}
}
}
}
}
return(misse.list)
}
######################################################################################################################################
######################################################################################################################################
### The function used to optimize parameters:
######################################################################################################################################
######################################################################################################################################
DevOptimizeMiSSE <- function(p, pars, dat.tab, gen, hidden.states, fixed.eps, nb.tip, nb.node, condition.on.survival, root.type, root.p, np, ode.eps, fossil.taxa, fix.type, strat.cache) {
#Generates the final vector with the appropriate parameter estimates in the right place:
p.new <- exp(p)
model.vec <- numeric(length(pars))
model.vec[] <- c(p.new, 0)[pars]
cache <- ParametersToPassMiSSE(model.vec=model.vec, hidden.states=hidden.states, fixed.eps=fixed.eps, nb.tip=nb.tip, nb.node=nb.node, bad.likelihood=exp(-300), ode.eps=ode.eps)
if(!is.null(fix.type)){
if(!is.null(strat.cache)){
logl <- DownPassMisse(dat.tab=dat.tab, gen=gen, cache=cache, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, node=fix.type[,1], state=NULL, fossil.taxa=fossil.taxa, fix.type=fix.type[,2]) + (strat.cache$k*log(cache$psi)) + (cache$psi*strat.cache$l_s)
}else{
logl <- DownPassMisse(dat.tab=dat.tab, gen=gen, cache=cache, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, node=fix.type[,1], state=NULL, fossil.taxa=fossil.taxa, fix.type=fix.type[,2])
}
}else{
logl <- DownPassMisse(dat.tab=dat.tab, gen=gen, cache=cache, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, node=NULL, fossil.taxa=fossil.taxa, fix.type=NULL)
}
return(-logl)
}
######################################################################################################################################
######################################################################################################################################
### The various utility functions used
######################################################################################################################################
######################################################################################################################################
GetTreeTable <- function(phy, root.age=NULL){
if(is.null(root.age)){
node.ages <- dateNodes(phy, rootAge=max(node.depth.edgelength(phy)))
}else{
node.ages <- dateNodes(phy, rootAge=root.age)
}
max.age <- max(node.ages)
table.ages <- matrix(0, dim(phy$edge)[1], 2)
for(row.index in 1:dim(phy$edge)[1]){
table.ages[row.index,1] <- node.ages[phy$edge[row.index,1]]
table.ages[row.index,2] <- node.ages[phy$edge[row.index,2]]
}
full.table <- cbind(table.ages, phy$edge.length, phy$edge)
return(full.table)
}
GetEdgeDetails <- function(phy, includes.intervals=FALSE, intervening.intervals=NULL){
table.info <- GetTreeTable(phy, root.age=NULL)
edges_detail <- as.data.frame(table.info)
colnames(edges_detail) <- c("rootward_age", "tipward_age", "edge.length", "rootward_node", "tipward_node")
edges_detail$type <- "internal"
edges_detail$type[which(edges_detail$tipward_node<=ape::Ntip(phy))] <- "extant_tip"
for(row.index in 1:dim(table.info)[1]){
if(table.info[row.index,5]<=ape::Ntip(phy)){
if(table.info[row.index,2] > .Machine$double.eps^.50){
edges_detail$type[row.index] <- "extinct_tip"
}
}
}
edges_detail$type[which(edges_detail$tipward_node %in% which(grepl("Ksamp_",phy$tip.label)))] <- "k_tip"
if(includes.intervals == TRUE){
Ksamples_rootward_nodes <- edges_detail$rootward_node[which(edges_detail$type=="k_tip")]
Extinct_tipward_nodes <- edges_detail$tipward_node[which(edges_detail$type=="extinct_tip")]
Extant_tipward_nodes <- edges_detail$tipward_node[which(edges_detail$type=="extant_tip")]
tipwards_nodes_of_rootwards_nodes <- edges_detail$tipward_node[edges_detail$rootward_node %in% Ksamples_rootward_nodes]
for(rootward_focal in Ksamples_rootward_nodes) {
focal_rows <- which(edges_detail$rootward_node == rootward_focal)
for(row_index in focal_rows) {
###For k --> k
if(edges_detail$tipward_node[row_index] %in% Ksamples_rootward_nodes) {
edges_detail$type[row_index] <- "k_k_interval"
}
###For k --> extinct
if(edges_detail$tipward_node[row_index] %in% Extinct_tipward_nodes) {
edges_detail$type[row_index] <- "k_extinct_interval"
}
###For k --> extant
if(edges_detail$tipward_node[row_index] %in% Extant_tipward_nodes) {
edges_detail$type[row_index] <- "k_extant_interval"
}
}
}
#First go and find the intervening intervals. Should be k_k_intervals, so use intervening interval table to find and replace:
if(!is.null(intervening.intervals)){
tmp <- which(round(edges_detail$tipward_age,10) %in% round(intervening.intervals$o_i,10))
for(match.index in 1:length(tmp)){
if(edges_detail[tmp[match.index],]$type == "k_k_interval"){
edges_detail[tmp[match.index],]$type <- "intervening_interval"
}
}
}
#Next, check any tip intervals that are actually subtended by a k_k interval:
for(check.index in 1:dim(edges_detail)[1]){
if(edges_detail$type[check.index] == "k_extant_interval"){
sister.taxa <- GetSister(phy, edges_detail$rootward_node[check.index])
if(edges_detail$type[which(edges_detail$tipward_node==sister.taxa)] == "k_tip"){
rootward.of.sister.taxa <- edges_detail$rootward_node[which(edges_detail$tipward_node==sister.taxa)]
tmp <- edges_detail[which(edges_detail$rootward_node==rootward.of.sister.taxa),]
if(any(tmp$type == "k_k_interval")){
edges_detail$type[check.index] <- "extant_tip"
}
}
}
}
#Finally, check any extinct tip intervals that are actually subtended by a k_k interval:
for(check.index in 1:dim(edges_detail)[1]){
if(edges_detail$type[check.index] == "k_extinct_interval"){
sister.taxa <- GetSister(phy, edges_detail$rootward_node[check.index])
if(edges_detail$type[which(edges_detail$tipward_node==sister.taxa)] == "k_tip"){
rootward.of.sister.taxa <- edges_detail$rootward_node[which(edges_detail$tipward_node==sister.taxa)]
tmp <- edges_detail[which(edges_detail$rootward_node==rootward.of.sister.taxa),]
if(any(tmp$type == "k_k_interval")){
edges_detail$type[check.index] <- "extinct_tip"
}
}
}
}
}
return(edges_detail)
}
OrganizeDataMiSSE <- function(phy, f, hidden.states, includes.intervals=FALSE, intervening.intervals=NULL, includes.fossils=FALSE){
### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
DesNode = NULL
nb.tip <- length(phy$tip.label)
nb.node <- phy$Nnode
compD <- matrix(0, nrow=nb.tip, ncol=26)
compE <- matrix(0, nrow=nb.tip, ncol=26)
#Initializes the tip sampling and sets internal nodes to be zero:
ncols = hidden.states
if(length(f) == 1){
for(i in 1:(nb.tip)){
compD[i,1:hidden.states] <- f
compE[i,1:hidden.states] <- rep((1-f), ncols)
}
}else{
for(i in 1:(nb.tip)){
compD[i,] <- f[i]
compE[i,] <- rep((1-f[i]), ncols)
}
}
#This seems stupid but I cannot figure out how to get data.table to not make this column a factor. When a factor this is not right. For posterity, let it be known Jeremy would rather just retain the character instead of this mess:
edge_details <- GetEdgeDetails(phy, includes.intervals=includes.intervals, intervening.intervals=intervening.intervals)
if(includes.fossils == TRUE){
branch.type <- edge_details$type
branch.type[which(branch.type == "extant_tip")] <- 0
branch.type[which(branch.type == "internal")] <- 0
branch.type[which(branch.type == "extinct_tip")] <- 1
branch.type[which(branch.type == "k_tip")] <- 2
branch.type[which(branch.type == "k_k_interval")] <- 3
branch.type[which(branch.type == "k_extinct_interval")] <- 3
branch.type[which(branch.type == "k_extant_interval")] <- 3
branch.type[which(branch.type == "intervening_interval")] <- 4
}else{
branch.type <- rep(0, length(edge_details$type))
}
tmp.df <- cbind(edge_details[,1:5], 0, matrix(0, nrow(edge_details), ncol(compD)), matrix(0, nrow(edge_details), ncol(compE)), as.numeric(branch.type))
colnames(tmp.df) <- c("RootwardAge", "TipwardAge", "BranchLength", "FocalNode", "DesNode", "comp", paste("compD", 1:ncol(compD), sep="_"), paste("compE", 1:ncol(compE), sep="_"), "branch.type")
dat.tab <- as.data.table(tmp.df)
setkey(dat.tab, DesNode)
cols <- names(dat.tab)
for (j in 1:(dim(compD)[2])){
#dat.tab[data.table(c(1:nb.tip)), paste("compD", j, sep="_") := compD[,j]]
set(dat.tab, 1:nb.tip, cols[6+j], compD[,j])
#dat.tab[data.table(c(1:nb.tip)), paste("compE", j, sep="_") := compE[,j]]
set(dat.tab, 1:nb.tip, cols[32+j], compE[,j])
}
return(dat.tab)
}
SingleChildProbMiSSE <- function(cache, pars, compD, compE, start.time, end.time, branch.type){
if(any(!is.finite(c(compD, compE)))) { # something went awry at a previous step. Bail!
prob.subtree.cal <- rep(0, 26*2)
prob.subtree.cal[(1+length(prob.subtree.cal)/2):length(prob.subtree.cal)] <- cache$bad.likelihood
return(prob.subtree.cal)
}
yini <- c(E0A = compE[1], E0B = compE[2], E0C = compE[3], E0D = compE[4], E0E = compE[5], E0F = compE[6], E0G = compE[7], E0H = compE[8], E0I = compE[9], E0J = compE[10], E0K = compE[11], E0L = compE[12], E0M = compE[13], E0N = compE[14], E0O = compE[15], E0P = compE[16], E0Q = compE[17], E0R = compE[18], E0S = compE[19], E0T = compE[20], E0U = compE[21], E0V = compE[22], E0W = compE[23], E0X = compE[24], E0Y = compE[25], E0Z = compE[26], D0A = compD[1], D0B = compD[2], D0C = compD[3], D0D = compD[4], D0E = compD[5], D0F = compD[6], D0G = compD[7], D0H = compD[8], D0I = compD[9], D0J = compD[10], D0K = compD[11], D0L = compD[12], D0M = compD[13], D0N = compD[14], D0O = compD[15], D0P = compD[16], D0Q = compD[17], D0R = compD[18], D0S = compD[19], D0T = compD[20], D0U = compD[21], D0V = compD[22], D0W = compD[23], D0X = compD[24], D0Y = compD[25], D0Z = compD[26])
if(branch.type == 2){
times=c(0, end.time)
}else{
times=c(start.time, end.time)
}
runSilent <- function(branch.type) {
options(warn = -1)
on.exit(options(warn = 0))
if(branch.type == 3){
capture.output(res <- lsoda(yini, times, func = "misse_strat_derivs", pars, initfunc="initmod_misse", dllname = "hisse", rtol=1e-8, atol=1e-8))
}else{
capture.output(res <- lsoda(yini, times, func = "misse_derivs", pars, initfunc="initmod_misse", dllname = "hisse", rtol=1e-8, atol=1e-8))
}
res
}
#prob.subtree.cal.full <- lsoda(yini, times, func = "misse_derivs", pars, initfunc="initmod_misse", dll = "misse-ext-derivs", rtol=1e-8, atol=1e-8)
#prob.subtree.cal.full <- lsoda(yini, times, func = "misse_derivs", pars, initfunc="initmod_misse", dllname = "hisse", rtol=1e-8, atol=1e-8)
prob.subtree.cal.full <- runSilent(branch.type=branch.type)
######## THIS CHECKS TO ENSURE THAT THE INTEGRATION WAS SUCCESSFUL ###########
if(attributes(prob.subtree.cal.full)$istate[1] < 0){
prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
prob.subtree.cal[27:52] <- cache$bad.likelihood
return(prob.subtree.cal)
}else{
if(branch.type == 2){
prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
#So, what this does is simply maintains the original input which here should just be a bunch of ones.
prob.subtree.cal[27:52] <- compD
}else{
if(branch.type == 4){
prob.subtree.cal.wevents <- prob.subtree.cal.full[-1,-1]
prob.subtree.cal.full <- runSilent(branch.type=3)
prob.subtree.cal.noevents <- prob.subtree.cal.full[-1,-1]
unobs.spec.probs <- 1 - (prob.subtree.cal.noevents[27:52]/prob.subtree.cal.wevents[27:52])
prob.subtree.cal.wevents[27:52] <- prob.subtree.cal.wevents[27:52] * unobs.spec.probs
prob.subtree.cal.wevents[which(is.na(prob.subtree.cal.wevents))] <- 0
prob.subtree.cal <- prob.subtree.cal.wevents
}else{
prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
}
}
}
##############################################################################
if(any(is.nan(prob.subtree.cal[27:52]))){
prob.subtree.cal[27:52] <- cache$bad.likelihood
return(prob.subtree.cal)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(any(prob.subtree.cal[27:52] < 0)){
prob.subtree.cal[27:52] <- cache$bad.likelihood
return(prob.subtree.cal)
}
if(sum(prob.subtree.cal[27:52]) < cache$ode.eps){
prob.subtree.cal[27:52] <- cache$bad.likelihood
return(prob.subtree.cal)
}
return(prob.subtree.cal)
}
FocalNodeProbMiSSE <- function(cache, pars, lambdas, dat.tab, generations){
### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
DesNode = NULL
FocalNode = NULL
. = NULL
gens <- data.table(c(generations))
setkey(dat.tab, FocalNode)
CurrentGenData <- dat.tab[gens]
tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProbMiSSE(cache, pars, z[7:32], z[33:58], z[2], z[1], z[59])))
v.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),27:52] * tmp[seq(2,nrow(tmp),2),27:52], length(unique(CurrentGenData$FocalNode)), 26)
v.mat <- v.mat * matrix(lambdas, length(unique(CurrentGenData$FocalNode)), 26, byrow=TRUE)
phi.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),1:26], length(unique(CurrentGenData$FocalNode)), 26)
if(!is.null(cache$node)){
if(any(cache$node %in% generations)){
for(fix.index in 1:length(cache$node)){
if(cache$fix.type[fix.index] == "event"){
#basically we are using the node to fix the state along a branch, but we do not want to assume a true speciation event occurred here.
lambdas.check <- lambdas
lambdas.check[which(lambdas==0)] <- 1
#The initial condition for a k.sample is D(t)*psi
v.mat[which(generations == cache$node[fix.index]),] <- (v.mat[which(generations == cache$node[fix.index]),] / lambdas.check) * cache$psi
}else{
if(cache$fix.type[fix.index] == "interval"){
lambdas.check <- lambdas
lambdas.check[which(lambdas==0)] <- 1
v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] / lambdas.check
}else{
fixer = numeric(26)
fixer[cache$state] = 1
v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] * fixer
}
}
}
}
}
tmp.comp <- rowSums(v.mat)
tmp.probs <- v.mat / tmp.comp
#tmp.probs <- v.mat
setkey(dat.tab, DesNode)
rows <- dat.tab[.(generations), which=TRUE]
cols <- names(dat.tab)
for (j in 1:(dim(tmp.probs)[2])){
#dat.tab[data.table(c(generations)), paste("compD", j, sep="_") := tmp.probs[,j]]
set(dat.tab, rows, cols[6+j], tmp.probs[,j])
#dat.tab[data.table(c(generations)), paste("compE", j, sep="_") := phi.mat[,j]]
set(dat.tab, rows, cols[32+j], phi.mat[,j])
}
dat.tab[data.table(c(generations)), "comp" := tmp.comp]
return(dat.tab)
}
#Have to calculate root prob separately because it is not a descendant in our table. Could add it, but I worry about the NA that is required.
GetRootProbMiSSE <- function(cache, pars, lambdas, dat.tab, generations){
### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
FocalNode = NULL
gens <- data.table(c(generations))
setkey(dat.tab, FocalNode)
CurrentGenData <- dat.tab[gens]
tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProbMiSSE(cache, pars, z[7:32], z[33:58], z[2], z[1], z[59])))
v.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),27:52] * tmp[seq(2,nrow(tmp),2),27:52], length(unique(CurrentGenData$FocalNode)), 26)
v.mat <- v.mat * matrix(lambdas, length(unique(CurrentGenData$FocalNode)), 26, byrow=TRUE)
phi.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),1:26], length(unique(CurrentGenData$FocalNode)), 26)
if(!is.null(cache$node)){
if(any(cache$node %in% generations)){
for(fix.index in 1:length(cache$node)){
if(cache$fix.type[fix.index] == "event"){
#basically we are using the node to fix the state along a branch, but we do not want to assume a true speciation event occurred here.
lambdas.check <- lambdas
lambdas.check[which(lambdas==0)] <- 1
#The initial condition for a k.sample is D(t)*psi
v.mat[which(generations == cache$node[fix.index]),] <- (v.mat[which(generations == cache$node[fix.index]),] / lambdas.check) * cache$psi
}else{
if(cache$fix.type[fix.index] == "interval"){
#basically we are using the node to fix the state along a branch, but we do not want to assume a true speciation event occurred here.
lambdas.check <- lambdas
lambdas.check[which(lambdas==0)] <- 1
v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] / lambdas.check
}else{
fixer = numeric(26)
fixer[cache$state] = 1
v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] * fixer
}
}
}
}
}
tmp.comp <- rowSums(v.mat)
tmp.probs <- v.mat / tmp.comp
#tmp.probs <- v.mat
return(cbind(tmp.comp, phi.mat, tmp.probs))
}
#Calculates the initial conditions for fossil taxa in the tree.
GetFossilInitialsMiSSE <- function(cache, pars, lambdas, dat.tab, fossil.taxa){
### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
DesNode = NULL
. = NULL
fossils <- data.table(c(fossil.taxa))
setkey(dat.tab, DesNode)
CurrentGenData <- dat.tab[fossils]
tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProbMiSSE(cache, pars, z[7:32], z[33:58], 0, z[2], 1)))
tmp.probs <- matrix(tmp[,1:26], length(fossil.taxa), 26) * as.matrix(CurrentGenData[,7:32]) * cache$psi
phi.mat <- matrix(tmp[,1:26], length(fossil.taxa), 26)
setkey(dat.tab, DesNode)
rows <- dat.tab[.(fossils), which=TRUE]
cols <- names(dat.tab)
for (j in 1:(dim(tmp.probs)[2])){
#dat.tab[data.table(c(generations)), paste("compD", j, sep="_") := tmp.probs[,j]]
set(dat.tab, rows, cols[6+j], tmp.probs[,j])
#dat.tab[data.table(c(generations)), paste("compE", j, sep="_") := phi.mat[,j]]
set(dat.tab, rows, cols[32+j], phi.mat[,j])
}
return(dat.tab)
}
######################################################################################################################################
######################################################################################################################################
### The MiSSE type down pass that carries out the integration and returns the likelihood:
######################################################################################################################################
######################################################################################################################################
DownPassMisse <- function(dat.tab, gen, cache, condition.on.survival, root.type, root.p, get.phi=FALSE, node=NULL, state=NULL, fossil.taxa=NULL, fix.type=NULL) {
### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
DesNode = NULL
compE = NULL
. = NULL
pars <- c(cache$lambda0A, cache$mu0A, cache$lambda0B, cache$mu0B, cache$lambda0C, cache$mu0C, cache$lambda0D, cache$mu0D, cache$lambda0E, cache$mu0E, cache$lambda0F, cache$mu0F, cache$lambda0G, cache$mu0G, cache$lambda0H, cache$mu0H, cache$lambda0I, cache$mu0I, cache$lambda0J, cache$mu0J, cache$lambda0K, cache$mu0K, cache$lambda0L, cache$mu0L, cache$lambda0M, cache$mu0M, cache$lambda0N, cache$mu0N, cache$lambda0O, cache$mu0O, cache$lambda0P, cache$mu0P, cache$lambda0Q, cache$mu0Q, cache$lambda0R, cache$mu0R, cache$lambda0S, cache$mu0S, cache$lambda0T, cache$mu0T, cache$lambda0U, cache$mu0U, cache$lambda0V, cache$mu0V, cache$lambda0W, cache$mu0W, cache$lambda0X, cache$mu0X, cache$lambda0Y, cache$mu0Y, cache$lambda0Z, cache$mu0Z, cache$q0, cache$hidden.states, cache$psi)
lambda <- c(cache$lambda0A, cache$lambda0B, cache$lambda0C, cache$lambda0D, cache$lambda0E, cache$lambda0F, cache$lambda0G, cache$lambda0H, cache$lambda0I, cache$lambda0J, cache$lambda0K, cache$lambda0L, cache$lambda0M, cache$lambda0N, cache$lambda0O, cache$lambda0P, cache$lambda0Q, cache$lambda0R, cache$lambda0S, cache$lambda0T, cache$lambda0U, cache$lambda0V, cache$lambda0W, cache$lambda0X, cache$lambda0Y, cache$lambda0Z)
if(!is.null(fossil.taxa)){
#Gets initial conditions for fossil taxa:
dat.tab.copy <- copy(dat.tab)
dat.tab.copy <- GetFossilInitialsMiSSE(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, fossil.taxa=fossil.taxa)
}else{
dat.tab.copy <- copy(dat.tab)
}
TIPS <- 1:cache$nb.tip
for(i in 1:length(gen)){
if(i == length(gen)){
if(!is.null(node)){
if(any(node %in% gen[[i]])){
cache$node <- node
cache$state <- state
cache$fix.type <- fix.type
res.tmp <- GetRootProbMiSSE(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, generations=gen[[i]])
cache$node <- NULL
cache$state <- NULL
cache$fix.type <- NULL
}else{
res.tmp <- GetRootProbMiSSE(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, generations=gen[[i]])
}
}else{
res.tmp <- GetRootProbMiSSE(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, generations=gen[[i]])
}
compD.root <- res.tmp[c(28:53)]
compE.root <- res.tmp[c(2:27)]
setkey(dat.tab.copy, DesNode)
comp <- dat.tab.copy[["comp"]]
comp <- c(comp[-TIPS], res.tmp[1])
}else{
if(!is.null(node)){
if(any(node %in% gen[[i]])){
cache$node <- node
cache$state <- state
cache$fix.type <- fix.type
dat.tab.copy <- FocalNodeProbMiSSE(cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, gen[[i]])
cache$node <- NULL
cache$state <- NULL
cache$fix.type <- NULL
}else{
dat.tab.copy <- FocalNodeProbMiSSE(cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, gen[[i]])
}
}else{
dat.tab.copy <- FocalNodeProbMiSSE(cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, gen[[i]])
}
}
}
if(!is.null(fossil.taxa)){
#this overrides the post order traversal E and recalculates assuming psi=0. See pg. 400 Stadler 2010.
#if(all(fix.type == "event")){ ## not sure yet for intervals. This is unclear in Stadler et al 2018.
pars[length(pars)] <- 0
cache$psi <- 0
phi.mat <- SingleChildProbMiSSE(cache, pars, c(as.matrix(dat.tab[1,7:32])), c(as.matrix(dat.tab[1, 33:58])), 0, max(dat.tab$RootwardAge), 0)
compE.root <- matrix(phi.mat[1:26], 1, 26)
#}
}
if (is.na(sum(log(compD.root))) || is.na(log(sum(1-compE.root)))){
return(log(cache$bad.likelihood)^13)
}else{
if(root.type == "madfitz" | root.type == "herr_als"){
if(is.null(root.p)){
root.p = compD.root/sum(compD.root)
root.p[which(is.na(root.p))] = 0
}
}
if(condition.on.survival == TRUE){
if(root.type == "madfitz"){
#lambda <- c(cache$lambda0A, cache$lambda0B, cache$lambda0C, cache$lambda0D, cache$lambda0E, cache$lambda0F, cache$lambda0G, cache$lambda0H, cache$lambda0I, cache$lambda0J, cache$lambda0K, cache$lambda0L, cache$lambda0M, cache$lambda0N, cache$lambda0O, cache$lambda0P, cache$lambda0Q, cache$lambda0R, cache$lambda0S, cache$lambda0T, cache$lambda0U, cache$lambda0V, cache$lambda0W, cache$lambda0X, cache$lambda0Y, cache$lambda0Z)
compD.root <- compD.root / sum(root.p * lambda * (1 - compE.root)^2)
#Corrects for possibility that you have 0/0:
compD.root[which(is.na(compD.root))] = 0
loglik <- log(sum(compD.root * root.p)) + sum(log(comp))
#loglik <- log(sum(compD.root * root.p))
}else{
#lambda <- c(cache$lambda0A, cache$lambda0B, cache$lambda0C, cache$lambda0D, cache$lambda0E, cache$lambda0F, cache$lambda0G, cache$lambda0H, cache$lambda0I, cache$lambda0J, cache$lambda0K, cache$lambda0L, cache$lambda0M, cache$lambda0N, cache$lambda0O, cache$lambda0P, cache$lambda0Q, cache$lambda0R, cache$lambda0S, cache$lambda0T, cache$lambda0U, cache$lambda0V, cache$lambda0W, cache$lambda0X, cache$lambda0Y, cache$lambda0Z)
compD.root <- (compD.root * root.p) / (lambda * (1 - compE.root)^2)
#Corrects for possibility that you have 0/0:
compD.root[which(is.na(compD.root))] = 0
loglik <- log(sum(compD.root)) + sum(log(comp))
}
}
if(!is.finite(loglik)){
return(log(cache$bad.likelihood)^7)
}
}
if(get.phi==TRUE){
obj = NULL
obj$compD.root = compD.root/sum(compD.root)
obj$compE = compE.root
obj$root.p = root.p
return(obj)
}else{
return(loglik)
}
}
######################################################################################################################################
######################################################################################################################################
### Cache object for storing parameters that are used throughout MiSSE:
######################################################################################################################################
######################################################################################################################################
ParametersToPassMiSSE <- function(model.vec, hidden.states, fixed.eps, nb.tip, nb.node, bad.likelihood, ode.eps){
#Provides an initial object that contains all the parameters to be passed among functions. This will also be used to pass other things are we move down the tree (see DownPassGeoSSE):
obj <- NULL
obj$hidden.states <- hidden.states
obj$nb.tip <- nb.tip
obj$nb.node <- nb.node
obj$bad.likelihood <- bad.likelihood
obj$ode.eps <- ode.eps
if(is.null(fixed.eps)){
##Hidden State A
obj$lambda0A = model.vec[1] / (1 + model.vec[2])
obj$mu0A = (model.vec[2] * model.vec[1]) / (1 + model.vec[2])
##Hidden State B
obj$lambda0B = model.vec[3] / (1 + model.vec[4])
obj$mu0B = (model.vec[4] * model.vec[3]) / (1 + model.vec[4])
##Hidden State C
obj$lambda0C = model.vec[5] / (1 + model.vec[6])
obj$mu0C = (model.vec[6] * model.vec[5]) / (1 + model.vec[6])
##Hidden State D
obj$lambda0D = model.vec[7] / (1 + model.vec[8])
obj$mu0D = (model.vec[8] * model.vec[7]) / (1 + model.vec[8])
##Hidden State E
obj$lambda0E = model.vec[9] / (1 + model.vec[10])
obj$mu0E = (model.vec[10] * model.vec[9]) / (1 + model.vec[10])
##Hidden State F
obj$lambda0F = model.vec[11] / (1 + model.vec[12])
obj$mu0F = (model.vec[12] * model.vec[11]) / (1 + model.vec[12])
##Hidden State G
obj$lambda0G = model.vec[13] / (1 + model.vec[14])
obj$mu0G = (model.vec[14] * model.vec[13]) / (1 + model.vec[14])
##Hidden State H
obj$lambda0H = model.vec[15] / (1 + model.vec[16])
obj$mu0H = (model.vec[16] * model.vec[15]) / (1 + model.vec[16])
##Hidden State I
obj$lambda0I = model.vec[17] / (1 + model.vec[18])
obj$mu0I = (model.vec[18] * model.vec[17]) / (1 + model.vec[18])
##Hidden State J
obj$lambda0J = model.vec[19] / (1 + model.vec[20])
obj$mu0J = (model.vec[20] * model.vec[19]) / (1 + model.vec[20])
##Hidden State K
obj$lambda0K = model.vec[21] / (1 + model.vec[22])
obj$mu0K = (model.vec[22] * model.vec[21]) / (1 + model.vec[22])
##Hidden State L
obj$lambda0L = model.vec[23] / (1 + model.vec[24])
obj$mu0L = (model.vec[24] * model.vec[23]) / (1 + model.vec[24])
##Hidden State M
obj$lambda0M = model.vec[25] / (1 + model.vec[26])
obj$mu0M = (model.vec[26] * model.vec[25]) / (1 + model.vec[26])
##Hidden State N
obj$lambda0N = model.vec[27] / (1 + model.vec[28])
obj$mu0N = (model.vec[28] * model.vec[27]) / (1 + model.vec[28])
##Hidden State O
obj$lambda0O = model.vec[29] / (1 + model.vec[30])
obj$mu0O = (model.vec[30] * model.vec[29]) / (1 + model.vec[30])
##Hidden State P
obj$lambda0P = model.vec[31] / (1 + model.vec[32])
obj$mu0P = (model.vec[32] * model.vec[31]) / (1 + model.vec[32])
##Hidden State Q
obj$lambda0Q = model.vec[33] / (1 + model.vec[34])
obj$mu0Q = (model.vec[34] * model.vec[33]) / (1 + model.vec[34])
##Hidden State R
obj$lambda0R = model.vec[35] / (1 + model.vec[36])
obj$mu0R = (model.vec[36] * model.vec[35]) / (1 + model.vec[36])
##Hidden State S
obj$lambda0S = model.vec[37] / (1 + model.vec[38])
obj$mu0S = (model.vec[38] * model.vec[37]) / (1 + model.vec[38])
##Hidden State T
obj$lambda0T = model.vec[39] / (1 + model.vec[40])
obj$mu0T = (model.vec[40] * model.vec[39]) / (1 + model.vec[40])
##Hidden State U
obj$lambda0U = model.vec[41] / (1 + model.vec[42])
obj$mu0U = (model.vec[42] * model.vec[41]) / (1 + model.vec[42])
##Hidden State V
obj$lambda0V = model.vec[43] / (1 + model.vec[44])
obj$mu0V = (model.vec[44] * model.vec[43]) / (1 + model.vec[44])
##Hidden State W
obj$lambda0W = model.vec[45] / (1 + model.vec[46])
obj$mu0W = (model.vec[46] * model.vec[45]) / (1 + model.vec[46])
##Hidden State X
obj$lambda0X = model.vec[47] / (1 + model.vec[48])
obj$mu0X = (model.vec[48] * model.vec[47]) / (1 + model.vec[48])
##Hidden State Y
obj$lambda0Y = model.vec[49] / (1 + model.vec[50])
obj$mu0Y = (model.vec[50] * model.vec[49]) / (1 + model.vec[50])
##Hidden State Z
obj$lambda0Z = model.vec[51] / (1 + model.vec[52])
obj$mu0Z = (model.vec[52] * model.vec[51]) / (1 + model.vec[52])
}else{
##Hidden State A
obj$lambda0A = model.vec[1] / (1 + fixed.eps)
obj$mu0A = (fixed.eps * model.vec[1]) / (1 + fixed.eps)
##Hidden State B
obj$lambda0B = model.vec[3] / (1 + fixed.eps)
obj$mu0B = (fixed.eps * model.vec[3]) / (1 + fixed.eps)
##Hidden State C
obj$lambda0C = model.vec[5] / (1 + fixed.eps)
obj$mu0C = (fixed.eps * model.vec[5]) / (1 + fixed.eps)
##Hidden State D
obj$lambda0D = model.vec[7] / (1 + fixed.eps)
obj$mu0D = (fixed.eps * model.vec[7]) / (1 + fixed.eps)
##Hidden State E
obj$lambda0E = model.vec[9] / (1 + fixed.eps)
obj$mu0E = (fixed.eps * model.vec[9]) / (1 + fixed.eps)
##Hidden State F
obj$lambda0F = model.vec[11] / (1 + fixed.eps)
obj$mu0F = (fixed.eps * model.vec[11]) / (1 + fixed.eps)
##Hidden State G
obj$lambda0G = model.vec[13] / (1 + fixed.eps)
obj$mu0G = (fixed.eps * model.vec[13]) / (1 + fixed.eps)
##Hidden State H
obj$lambda0H = model.vec[15] / (1 + fixed.eps)
obj$mu0H = (fixed.eps * model.vec[15]) / (1 + fixed.eps)
##Hidden State I
obj$lambda0I = model.vec[17] / (1 + fixed.eps)
obj$mu0I = (fixed.eps * model.vec[17]) / (1 + fixed.eps)
##Hidden State J
obj$lambda0J = model.vec[19] / (1 + fixed.eps)
obj$mu0J = (fixed.eps * model.vec[19]) / (1 + fixed.eps)
##Hidden State K
obj$lambda0K = model.vec[21] / (1 + fixed.eps)
obj$mu0K = (fixed.eps * model.vec[21]) / (1 + fixed.eps)
##Hidden State L
obj$lambda0L = model.vec[23] / (1 + fixed.eps)
obj$mu0L = (fixed.eps * model.vec[23]) / (1 + fixed.eps)
##Hidden State M
obj$lambda0M = model.vec[25] / (1 + fixed.eps)
obj$mu0M = (fixed.eps * model.vec[25]) / (1 + fixed.eps)
##Hidden State N
obj$lambda0N = model.vec[27] / (1 + fixed.eps)
obj$mu0N = (fixed.eps * model.vec[27]) / (1 + fixed.eps)
##Hidden State O
obj$lambda0O = model.vec[29] / (1 + fixed.eps)
obj$mu0O = (fixed.eps * model.vec[29]) / (1 + fixed.eps)
##Hidden State P
obj$lambda0P = model.vec[31] / (1 + fixed.eps)
obj$mu0P = (fixed.eps * model.vec[31]) / (1 + fixed.eps)
##Hidden State Q
obj$lambda0Q = model.vec[33] / (1 + fixed.eps)
obj$mu0Q = (fixed.eps * model.vec[33]) / (1 + fixed.eps)
##Hidden State R
obj$lambda0R = model.vec[35] / (1 + fixed.eps)
obj$mu0R = (fixed.eps * model.vec[35]) / (1 + fixed.eps)
##Hidden State S
obj$lambda0S = model.vec[37] / (1 + fixed.eps)
obj$mu0S = (fixed.eps * model.vec[37]) / (1 + fixed.eps)
##Hidden State T
obj$lambda0T = model.vec[39] / (1 + fixed.eps)
obj$mu0T = (fixed.eps * model.vec[39]) / (1 + fixed.eps)
##Hidden State U
obj$lambda0U = model.vec[41] / (1 + fixed.eps)
obj$mu0U = (fixed.eps * model.vec[41]) / (1 + fixed.eps)
##Hidden State V
obj$lambda0V = model.vec[43] / (1 + fixed.eps)
obj$mu0V = (fixed.eps * model.vec[43]) / (1 + fixed.eps)
##Hidden State W
obj$lambda0W = model.vec[45] / (1 + fixed.eps)
obj$mu0W = (fixed.eps * model.vec[45]) / (1 + fixed.eps)
##Hidden State X
obj$lambda0X = model.vec[47] / (1 + fixed.eps)
obj$mu0X = (fixed.eps * model.vec[47]) / (1 + fixed.eps)
##Hidden State Y
obj$lambda0Y = model.vec[49] / (1 + fixed.eps)
obj$mu0Y = (fixed.eps * model.vec[49]) / (1 + fixed.eps)
##Hidden State Z
obj$lambda0Z = model.vec[51] / (1 + fixed.eps)
obj$mu0Z = (fixed.eps * model.vec[51]) / (1 + fixed.eps)
}
obj$q0 = model.vec[53]
obj$psi = model.vec[54]
return(obj)
}
print.misse.fit <- function(x,...){
## Function to print a "muhisse.fit" object.
set.zero <- max( x$index.par )
## Keep only the parameters estimated:
par.list <- x$solution[ !x$index.par == set.zero ]
ntips <- Ntip( x$phy )
nstates <- x$hidden.states
output <- c(x$loglik, x$AIC, x$AICc, ntips, nstates)
names(output) <- c("lnL", "AIC", "AICc", "n.taxa", "n.hidden.states")
cat("\n")
cat("Fit \n")
print(output)
cat("\n")
cat("Model parameters: \n")
cat("\n")
print(par.list)
cat("\n")
if(!is.null(x$fixed.eps)){
cat("Fixed eps used: \n")
cat("\n")
print(x$fixed.eps)
cat("\n")
}
if(!is.null(x$psi.type)){
if(x$psi.type == "m_only"){
cat("psi estimate reflects fossil tip sampling only. \n")
}
if(x$psi.type == "m+k"){
cat("psi estimate reflects both fossil edge and tip sampling. \n")
}
if(x$psi.type == "m+int"){
cat("psi estimate reflects stratigraphic interval sampling. \n")
}
}
}
#phy <- read.tree("../vignettes/whales_Steemanetal2009.tre")
#phy <- read.tree("whales_Slateretal2010.tre")
## print(p.new)
#gen <- hisse:::FindGenerations(phy)
#dat.tab <- hisse:::OrganizeDataMiSSE(phy=phy, f=1, hidden.states=1)
#nb.tip <- Ntip(phy)
#nb.node <- phy$Nnode
#model.vec <- c(0.103624, 5.207178e-09, rep(0,52), 1)
#cache = hisse:::ParametersToPassMiSSE(model.vec=model.vec, hidden.states=1, nb.tip=nb.tip, nb.node=nb.node, bad.likelihood=exp(-250), ode.eps=0)#
#logl <- hisse:::DownPassMisse(dat.tab=dat.tab, cache=cache, gen=gen, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL)
#right.logl <- -277.6942
#round(logl,4) == round(right.logl,4)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.