R/geohisse.old.R
In thej022214/hisse: Hidden State Speciation and Extinction
Defines functions
print.geohisse.old.fit
ParametersToPassMuSSE
ParametersToPassGeoHiSSE
DownPassMusse
DownPassGeoHisse
starting.point.geosse
DevOptimizeGeoHiSSE
GeoHiSSE.old
Documented in
GeoHiSSE.old
######################################################################################################################################
######################################################################################################################################
### GeoHiSSE -- Expanded set of GeoSSE models for examining diversification in relation to geographic range evolution
######################################################################################################################################
######################################################################################################################################
GeoHiSSE.old <- function(phy, data, f=c(1,1,1), speciation=c(1,2,3), extirpation=c(1,2), hidden.areas=FALSE, trans.rate=NULL, assume.cladogenetic=TRUE, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL, sann=TRUE, sann.its=1000, bounded.search=TRUE, max.tol=.Machine$double.eps^.50, mag.san.start=0.5, starting.vals=NULL, speciation.upper=1000, extirpation.upper=1000, trans.upper=100, ode.eps=0){
## Temporary fix for the current BUG:
if( !is.null(phy$node.label) ) phy$node.label <- NULL
if(!is.null(root.p)) {
if(hidden.areas ==TRUE){
if( length( root.p ) == 3 ){
root.p <- rep(root.p, 5)
root.p <- root.p / sum(root.p)
warning("For hidden areas, you need to specify the root.p for all hidden areas. We have adjusted it so that there's equal chance for among all hidden areas.")
} else{
root.p.new <- numeric(15)
root.p.new[1:length(root.p)] <- root.p
root.p <- root.p.new
root.p <- root.p / sum(root.p)
}
}else{
## All good:
root.p <- root.p / sum(root.p)
}
}
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(is.null(trans.rate)){
stop("Rate matrix needed. See TransMatMakerGeoHiSSE() to create one.")
}
if(hidden.areas == TRUE & dim(trans.rate)[1]<3){
stop("You chose a hidden state but this is not reflected in the transition matrix")
}
if(assume.cladogenetic == FALSE){
if(dim(trans.rate)[1] != length(extirpation)){
stop("You have not specified enough extinction parameters.")
}
}
## Return error message if the data is not in the correct format.
if( !inherits(data, what = c("matrix","data.frame")) ){
stop("'data' needs to be a matrix or data.frame with 2 columns. See help.")
}
if( !ncol( data ) == 2 ){
stop("'data' needs to be a matrix or data.frame with 2 columns. See help.")
}
## Check if the states are %in% 0:2:
states.check <- all( as.numeric(data[,2]) %in% 0:2 )
if( !states.check ){
stop("states need to be one of 0, 1, or 2. See help.")
}
## Check if 'hidden.areas' parameter is congruent with the speciation vector:
if( length(speciation) > 3 & !hidden.areas ){
stop("Speciation has more than 3 elements but 'hidden.areas' was set to FALSE. Please set 'hidden.areas' to TRUE if the model include more than one rate class.")
}
if( length(speciation) == 3 & hidden.areas ){
stop("Speciation has only 3 elements but 'hidden.areas' was set to TRUE. Please set 'hidden.areas' to FALSE if the model does not include hidden rate classes.")
}
if(assume.cladogenetic == TRUE){
pars <- numeric(115)
if(dim(trans.rate)[2]==3){
rate.cats <- 1
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
trans.tmp <- c(trans.rate["(0)", "(1)"], trans.rate["(0)", "(01)"], trans.rate["(1)", "(0)"], trans.rate["(1)", "(01)"], trans.rate["(01)", "(0)"], trans.rate["(01)", "(1)"])
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
category.rates.unique <- 0
pars.tmp <- c(pars.tmp, trans.tmp)
pars[1:11] <- pars.tmp
}
if(dim(trans.rate)[2]==6){
rate.cats <- 2
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 1), rep("(1A)", 1), rep("(01A)", 1), rep("(0B)", 1), rep("(1B)", 1), rep("(01B)", 1))
cols <- c("(0B)", "(1B)", "(01B)", "(0A)", "(1A)", "(01A)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- c(category.rate.shift[1], rep(0,3), category.rate.shift[2], rep(0,3), category.rate.shift[3], rep(0,3))
category.rate.shiftB <- c(category.rate.shift[4], rep(0,3), category.rate.shift[5], rep(0,3), category.rate.shift[6], rep(0,3))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==9){
rate.cats <- 3
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 2), rep("(1A)", 2), rep("(01A)", 2), rep("(0B)", 2), rep("(1B)", 2), rep("(01B)", 2), rep("(0C)", 2), rep("(1C)", 2), rep("(01C)", 2))
cols <- c("(0B)", "(0C)", "(1B)", "(1C)", "(01B)", "(01C)", "(0A)", "(0C)", "(1A)", "(1C)", "(01A)", "(01C)", "(0A)", "(0B)", "(1A)", "(1B)", "(01A)", "(01B)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- c(category.rate.shift[1:2], rep(0,2), category.rate.shift[3:4], rep(0,2), category.rate.shift[5:6], rep(0,2))
category.rate.shiftB <- c(category.rate.shift[7:8], rep(0,2), category.rate.shift[9:10], rep(0,2), category.rate.shift[11:12], rep(0,2))
category.rate.shiftC <- c(category.rate.shift[13:14], rep(0,2), category.rate.shift[15:16], rep(0,2), category.rate.shift[17:18], rep(0,2))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[5:6], trans.tmp[13:18], category.rate.shiftC)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==12){
rate.cats <- 4
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 3), rep("(1A)", 3), rep("(01A)", 3), rep("(0B)", 3), rep("(1B)", 3), rep("(01B)", 3), rep("(0C)", 3), rep("(1C)", 3), rep("(01C)", 3), rep("(0D)", 3), rep("(1D)", 3), rep("(01D)", 3))
cols <- c("(0B)", "(0C)", "(0D)", "(1B)", "(1C)", "(1D)", "(01B)", "(01C)", "(01D)", "(0A)", "(0C)", "(0D)", "(1A)", "(1C)", "(1D)", "(01A)", "(01C)", "(01D)", "(0A)", "(0B)", "(0D)", "(1A)", "(1B)", "(1D)", "(01A)", "(01B)", "(01D)", "(0A)", "(0B)", "(0C)", "(1A)", "(1B)", "(1C)", "(01A)", "(01B)", "(01C)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- c(category.rate.shift[1:3], rep(0,1), category.rate.shift[4:6], rep(0,1), category.rate.shift[7:9], rep(0,1))
category.rate.shiftB <- c(category.rate.shift[10:12], rep(0,1), category.rate.shift[13:15], rep(0,1), category.rate.shift[16:18], rep(0,1))
category.rate.shiftC <- c(category.rate.shift[19:21], rep(0,1), category.rate.shift[22:24], rep(0,1), category.rate.shift[25:27], rep(0,1))
category.rate.shiftD <- c(category.rate.shift[28:30], rep(0,1), category.rate.shift[31:33], rep(0,1), category.rate.shift[34:36], rep(0,1))
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[5:6], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[7:8], trans.tmp[19:24], category.rate.shiftD)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==15){
rate.cats <- 5
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)", "(0E)", "(0E)", "(1E)", "(1E)", "(01E)", "(01E)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)", "(1E)", "(01E)", "(0E)", "(01E)", "(0E)", "(1E)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 4), rep("(1A)", 4), rep("(01A)", 4), rep("(0B)", 4), rep("(1B)", 4), rep("(01B)", 4), rep("(0C)", 4), rep("(1C)", 4), rep("(01C)", 4), rep("(0D)", 4), rep("(1D)", 4), rep("(01D)", 4), rep("(0E)", 4), rep("(1E)", 4), rep("(01E)", 4))
cols <- c("(0B)", "(0C)", "(0D)", "(0E)", "(1B)", "(1C)", "(1D)", "(1E)", "(01B)", "(01C)", "(01D)", "(01E)", "(0A)", "(0C)", "(0D)", "(0E)", "(1A)", "(1C)", "(1D)", "(1E)", "(01A)", "(01C)", "(01D)", "(01E)", "(0A)", "(0B)", "(0D)", "(0E)", "(1A)", "(1B)", "(1D)", "(1E)", "(01A)", "(01B)", "(01D)", "(01E)", "(0A)", "(0B)", "(0C)", "(0E)", "(1A)", "(1B)", "(1C)", "(1E)", "(01A)", "(01B)", "(01C)", "(01E)", "(0A)", "(0B)", "(0C)", "(0D)", "(1A)", "(1B)", "(1C)", "(1D)", "(01A)", "(01B)", "(01C)", "(01D)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- category.rate.shift[1:12]
category.rate.shiftB <- category.rate.shift[13:24]
category.rate.shiftC <- category.rate.shift[25:36]
category.rate.shiftD <- category.rate.shift[37:48]
category.rate.shiftE <- category.rate.shift[49:60]
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD, category.rate.shiftE)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[5:6], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[7:8], trans.tmp[19:24], category.rate.shiftD, speciation[13:15], extirpation.tmp[9:10], trans.tmp[25:30], category.rate.shiftE)
pars[1:length(pars.tmp)] <- pars.tmp
}
np <- max(pars)
pars[pars==0] <- np+1
cat("Initializing...", "\n")
data.new <- data.frame(data[,2], data[,2], row.names=data[,1])
data.new <- data.new[phy$tip.label,]
#This is used to scale starting values to account for sampling:
if(length(f) == 3){
freqs <- table(data.new[,1])
if(length(freqs == 2)){
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}else{
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}
}else{
if(length(f) == Ntip(phy)){
stop("This is functionality has been temporarily removed.")
#samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / mean(f[as.numeric(names(freqs))+1]))
}else{
stop("The vector of sampling frequencies does not match the number of tips in the tree.")
}
}
if(sum(extirpation)==0){
init.pars <- starting.point.geosse(phy, eps=0, samp.freq.tree=samp.freq.tree)
}else{
init.pars <- starting.point.geosse(phy, eps=mag.san.start, samp.freq.tree=samp.freq.tree)
}
names(init.pars) <- NULL
if(is.null(starting.vals)){
def.set.pars <- rep(c(log(init.pars[1:3]), log(init.pars[4:5]), rep(log(init.pars[6:7]),3), rep(log(0.01), 12)), rate.cats)
}else{
## Check if 'starting.vals' has the correct format.
if( !length(starting.vals) %in% c(3,7) ){
stop("Wrong length of starting.vals")
}
if( length(starting.vals) == 7 ){
cat("Using developer mode for starting.vals.", "\n")
def.set.pars <- rep(c(log(starting.vals[1:3]), log(starting.vals[4:5]), rep(log(starting.vals[6:7]),3), rep(log(0.01), 12)), rate.cats)
} else{
def.set.pars <- rep(c(log( rep(starting.vals[1],3) ), log( rep(starting.vals[2],2) ), rep(log(starting.vals[3]),6), rep(log(0.01), 12)), rate.cats)
}
}
if(bounded.search == TRUE){
upper.full <- rep(c(rep(log(speciation.upper),3), rep(log(extirpation.upper),2), rep(log(trans.upper),2*3), rep(log(10), 12)), rate.cats)
}else{
upper.full <- rep(21,length(def.set.pars))
}
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]]
}
lower <- rep(-20, length(ip))
}else{
pars <- numeric(120)
if(dim(trans.rate)[2]==3){
rate.cats <- 1
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
trans.tmp <- c(trans.rate["(0)", "(1)"], trans.rate["(0)", "(01)"], trans.rate["(1)", "(0)"], trans.rate["(1)", "(01)"], trans.rate["(01)", "(0)"], trans.rate["(01)", "(1)"])
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
category.rates.unique <- 0
pars.tmp <- c(pars.tmp, trans.tmp)
pars[1:12] <- pars.tmp
}
if(dim(trans.rate)[2]==6){
rate.cats <- 2
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- c(category.rate.shift[1], rep(0,3), category.rate.shift[2], rep(0,3), category.rate.shift[3], rep(0,3))
category.rate.shiftB <- c(category.rate.shift[4], rep(0,3), category.rate.shift[5], rep(0,3), category.rate.shift[6], rep(0,3))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==9){
rate.cats <- 3
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- c(category.rate.shift[1:2], rep(0,2), category.rate.shift[3:4], rep(0,2), category.rate.shift[5:6], rep(0,2))
category.rate.shiftB <- c(category.rate.shift[7:8], rep(0,2), category.rate.shift[9:10], rep(0,2), category.rate.shift[11:12], rep(0,2))
category.rate.shiftC <- c(category.rate.shift[13:14], rep(0,2), category.rate.shift[15:16], rep(0,2), category.rate.shift[17:18], rep(0,2))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[7:9], trans.tmp[13:18], category.rate.shiftC)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==12){
rate.cats <- 4
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- c(category.rate.shift[1:3], rep(0,1), category.rate.shift[4:6], rep(0,1), category.rate.shift[7:9], rep(0,1))
category.rate.shiftB <- c(category.rate.shift[10:12], rep(0,1), category.rate.shift[13:15], rep(0,1), category.rate.shift[16:18], rep(0,1))
category.rate.shiftC <- c(category.rate.shift[19:21], rep(0,1), category.rate.shift[22:24], rep(0,1), category.rate.shift[25:27], rep(0,1))
category.rate.shiftD <- c(category.rate.shift[28:30], rep(0,1), category.rate.shift[31:33], rep(0,1), category.rate.shift[34:36], rep(0,1))
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[7:9], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[10:12], trans.tmp[19:24], category.rate.shiftD)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==15){
rate.cats <- 5
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)", "(0E)", "(0E)", "(1E)", "(1E)", "(01E)", "(01E)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)", "(1E)", "(01E)", "(0E)", "(01E)", "(0E)", "(1E)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- category.rate.shift[1:12]
category.rate.shiftB <- category.rate.shift[13:24]
category.rate.shiftC <- category.rate.shift[25:36]
category.rate.shiftD <- category.rate.shift[37:48]
category.rate.shiftE <- category.rate.shift[49:60]
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD, category.rate.shiftE)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[7:9], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[10:12], trans.tmp[19:24], category.rate.shiftD, speciation[13:15], extirpation.tmp[13:15], trans.tmp[25:30], category.rate.shiftE)
pars[1:length(pars.tmp)] <- pars.tmp
}
np <- max(pars)
pars[pars==0] <- np+1
cat("Initializing...", "\n")
data.new <- data.frame(data[,2], data[,2], row.names=data[,1])
data.new <- data.new[phy$tip.label,]
#This is used to scale starting values to account for sampling:
if(length(f) == 3){
freqs <- table(data.new[,1])
if(length(freqs == 2)){
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}else{
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}
}else{
if(length(f) == Ntip(phy)){
stop("This functionality has been temporarily removed.")
#samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / mean(f[as.numeric(names(freqs))+1]))
}else{
stop("The vector of sampling frequencies does not match the number of tips in the tree.")
}
}
if(sum(extirpation)==0){
init.pars <- starting.point.generator(phy, 3, samp.freq.tree, yule=TRUE)
}else{
init.pars <- starting.point.generator(phy, 3, samp.freq.tree, yule=FALSE)
if(init.pars[4] == 0){
init.pars[4:6] = 1e-6
}
}
names(init.pars) <- NULL
if(is.null(starting.vals)){
def.set.pars <- rep(c(log(init.pars[1:3]), log(init.pars[4:6]), log(init.pars[7:12]), rep(log(.01), 12)), rate.cats)
}else{
## Check if 'starting.vals' has the correct format.
if( !length(starting.vals) %in% c(3,12) ){
stop("Wrong length of starting.vals")
}
if( length(starting.vals) == 12 ){
cat("Using developer mode for starting.vals.", "\n")
def.set.pars <- rep(c(log(starting.vals[1:3]), log(starting.vals[4:6]), log(starting.vals[7:12]), rep(log(0.01), 12)), rate.cats)
} else{
def.set.pars <- rep(c(log( rep(starting.vals[1],3) ), log( rep(starting.vals[2],2) ), rep(log(starting.vals[3]),6), rep(log(0.01), 12)), rate.cats)
}
}
if(bounded.search == TRUE){
upper.full <- rep(c(rep(log(speciation.upper),3), rep(log(extirpation.upper),3), rep(log(trans.upper),6), rep(log(10), 12)), rate.cats)
}else{
upper.full <- rep(21,length(def.set.pars))
}
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]]
}
lower <- rep(-20, length(ip))
}
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=DevOptimizeGeoHiSSE, ub=upper, lb=lower, opts=opts, pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
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=DevOptimizeGeoHiSSE, control=list(reltol=max.tol, parscale=rep(0.1, length(ip))), pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
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=DevOptimizeGeoHiSSE, lower=lower, upper=upper, control=list(max.call=sann.its), pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
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=DevOptimizeGeoHiSSE, ub=upper, lb=lower, opts=opts, pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
solution <- numeric(length(pars))
solution[] <- c(exp(out$solution), 0)[pars]
loglik = -out$objective
}
if(assume.cladogenetic == TRUE){
names(solution) <- c("s0A", "s1A", "s01A", "x0A", "x1A", "d0A_1A ", "d0A_01A", "d1A_0A", "d1A_01A", "d01A_0A", "d01A_1A", "d0A_0B", "d0A_0C", "d0A_0D", "d0A_0E", "d1A_1B", "d1A_1C", "d1A_1D", "d1A_1E", "d01A_01B", "d01A_01C", "d01A_01D", "d01A_01E", "s0B", "s1B", "s01B", "x0B", "x1B", "d0B_1B ", "d0B_01B", "d1B_0B", "d1B_01B", "d01B_0B", "d01B_1B", "d0B_0A", "d0B_0C", "d0B_0D", "d0B_0E", "d1B_1A", "d1B_1C", "d1B_1D", "d1B_1E", "d01B_01A", "d01B_01C", "d01B_01D", "d01B_01E", "s0C", "s1C", "s01C", "x0C", "x1C", "d0C_1C ", "d0C_01C", "d1C_0C", "d1C_01C", "d01C_0C", "d01C_1C", "d0C_0A", "d0C_0B", "d0C_0D", "d0C_0E", "d1C_1A", "d1C_1B", "d1C_1D", "d1C_1E", "d01C_01A", "d01C_01B", "d01C_01D", "d01C_01E", "s0D", "s1D", "s01D", "x0D", "x1D", "d0D_1D ", "d0D_01D", "d1D_0D", "d1D_01D", "d01D_0D", "d01D_1D", "d0D_0A", "d0D_0B", "d0D_0C", "d0D_0E", "d1D_1A", "d1D_1B", "d1D_1C", "d1D_1E", "d01D_01A", "d01D_01B", "d01D_01C", "d01D_01E", "s0E", "s1E", "s01E", "x0E", "x1E", "d0E_1E ", "d0E_01E", "d1E_0E", "d1E_01E", "d01E_0E", "d01E_1E", "d0E_0A", "d0E_0B", "d0E_0C", "d0E_0D", "d1E_1A", "d1E_1B", "d1E_1C", "d1E_1D", "d01E_01A", "d01E_01B", "d01E_01C", "d01E_01D")
}else{
names(solution) <- c("s0A", "s1A", "s01A", "x0A", "x1A", "x01A", "d0A_1A ", "d0A_01A", "d1A_0A", "d1A_01A", "d01A_0A", "d01A_1A", "d0A_0B", "d0A_0C", "d0A_0D", "d0A_0E", "d1A_1B", "d1A_1C", "d1A_1D", "d1A_1E", "d01A_01B", "d01A_01C", "d01A_01D", "d01A_01E", "s0B", "s1B", "s01B", "x0B", "x1B", "x01B", "d0B_1B ", "d0B_01B", "d1B_0B", "d1B_01B", "d01B_0B", "d01B_1B", "d0B_0A", "d0B_0C", "d0B_0D", "d0B_0E", "d1B_1A", "d1B_1C", "d1B_1D", "d1B_1E", "d01B_01A", "d01B_01C", "d01B_01D", "d01B_01E", "s0C", "s1C", "s01C", "x0C", "x1C", "x01C", "d0C_1C ", "d0C_01C", "d1C_0C", "d1C_01C", "d01C_0C", "d01C_1C", "d0C_0A", "d0C_0B", "d0C_0D", "d0C_0E", "d1C_1A", "d1C_1B", "d1C_1D", "d1C_1E", "d01C_01A", "d01C_01B", "d01C_01D", "d01C_01E", "s0D", "s1D", "s01D", "x0D", "x1D", "x01D", "d0D_1D ", "d0D_01D", "d1D_0D", "d1D_01D", "d01D_0D", "d01D_1D", "d0D_0A", "d0D_0B", "d0D_0C", "d0D_0E", "d1D_1A", "d1D_1B", "d1D_1C", "d1D_1E", "d01D_01A", "d01D_01B", "d01D_01C", "d01D_01E", "s0E", "s1E", "s01E", "x0E", "x1E", "x01E", "d0E_1E ", "d0E_01E", "d1E_0E", "d1E_01E", "d01E_0E", "d01E_1E", "d0E_0A", "d0E_0B", "d0E_0C", "d0E_0D", "d1E_1A", "d1E_1B", "d1E_1C", "d1E_1D", "d01E_01A", "d01E_01B", "d01E_01C", "d01E_01D")
}
cat("Finished. Summarizing results...", "\n")
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.areas=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, phy=phy, data=data, trans.matrix=trans.rate, max.tol=max.tol, starting.vals=ip, upper.bounds=upper, lower.bounds=lower, ode.eps=ode.eps)
class(obj) <- append(class(obj), "geohisse.old.fit")
return(obj)
}
######################################################################################################################################
######################################################################################################################################
### The function used to optimize parameters:
######################################################################################################################################
######################################################################################################################################
DevOptimizeGeoHiSSE <- function(p, pars, phy, data, f, hidden.states, assume.cladogenetic, condition.on.survival, root.type, root.p, np, ode.eps) {
#Generates the final vector with the appropriate parameter estimates in the right place:
p.new <- exp(p)
## print(p.new)
model.vec <- numeric(length(pars))
model.vec[] <- c(p.new, 0)[pars]
if(assume.cladogenetic == TRUE){
cache = ParametersToPassGeoHiSSE(phy=phy, data=data, f=f, model.vec=model.vec, hidden.states=hidden.states)
logl <- DownPassGeoHisse(phy, cache, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, ode.eps=ode.eps)
}else{
cache = ParametersToPassMuSSE(phy=phy, data=data, f=f, model.vec=model.vec, hidden.states=hidden.states)
logl <- DownPassMusse(phy, cache, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, ode.eps=ode.eps)
}
return(-logl)
}
#Taken from the GeoSSE code modified to account for sampling -- credit goes to Emma Goldberg:
starting.point.geosse <- function(tree, eps=0.5, samp.freq.tree) {
if (eps == 0) {
s <- (log(Ntip(tree)) - log(2)) / max(branching.times(tree))
s <- s/samp.freq.tree
x <- 0
d <- s/10
} else {
n <- Ntip(tree)
r <- ( log( (n/2) * (1 - eps*eps) + 2*eps + (1 - eps)/2 *
sqrt( n * (n*eps*eps - 8*eps + 2*n*eps + n))) - log(2)
) / max(branching.times(tree))
s <- r / (1 - eps)
x <- s * eps
s <- s/samp.freq.tree
x <- x - ((s * samp.freq.tree) * (1 - 1/samp.freq.tree))
d <- x
}
p <- c(s, s, s, x, x, d, d)
names(p) <- c("sA", "sB", "sAB", "xA" , "xB" ,"dA" ,"dB")
p
}
######################################################################################################################################
######################################################################################################################################
### The GeoHiSSE down pass that carries out the integration and returns the likelihood:
######################################################################################################################################
######################################################################################################################################
DownPassGeoHisse <- function(phy, cache, hidden.states, bad.likelihood=-10000000, condition.on.survival, root.type, root.p, get.phi=FALSE, node=NULL, state=NULL, ode.eps=0) {
#Some preliminaries:
nb.tip <- length(phy$tip.label)
nb.node <- phy$Nnode
phy <- reorder(phy, "pruningwise")
anc <- unique(phy$edge[,1])
TIPS <- 1:nb.tip
if(hidden.states == FALSE){
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
}else{
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
if(!is.null(root.p)){
root.p.new <- numeric(15)
root.p.new[1:length(root.p)] <- root.p
root.p <- root.p.new
}
}
#Initializes the tip sampling and sets internal nodes to be zero:
ncols = dim(compD)[2]
if(length(cache$f) == 3){
for(i in 1:(nb.tip)){
compD[i,] <- cache$f * cache$states[i,]
compE[i,] <- rep((1-cache$f), ncols/3)
}
}else{
for(i in 1:(nb.tip)){
compD[i,] <- cache$f[i] * cache$states[i,]
compE[i,] <- rep((1-cache$f[i]), ncols/3)
}
}
logcomp <- c()
#Start the postorder traversal indexing lists by node number:
for (i in seq(from = 1, length.out = nb.node)) {
#A vector of all the internal nodes:
focal <- anc[i]
desRows <- which(phy$edge[,1]==focal)
desNodes <- phy$edge[desRows,2]
#Note: when the tree has been reordered branching.times are no longer valid. Fortunately, we extract this information in the initial cache setup. Also focal is the rootward node, whereas desNodes represent a vector of all descendant nodes:
cache$rootward.age <- cache$split.times[which(names(cache$split.times)==focal)]
v <- c()
phi <- c()
for (desIndex in sequence(length(desRows))){
cache$focal.edge.length <- phy$edge.length[desRows[desIndex]]
cache$tipward.age <- cache$rootward.age - cache$focal.edge.length
#Strange rounding errors. A tip age should be zero. This ensures that:
if(cache$tipward.age < .Machine$double.eps^0.5){
cache$tipward.age = 0
}
cache$node.D <- compD[desNodes[desIndex],]
cache$node.E <- compE[desNodes[desIndex],]
##Call to lsoda that utilizes C code. Requires a lot of inputs. Note that for now we hardcode the NUMELEMENTS arguments. The reason for this is because with lsoda we can only pass a vector of parameters.
if(hidden.states == FALSE){
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A)
NUMELEMENTS <- 11 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <-c(E_0=cache$node.E[1], E_1=cache$node.E[2], E_01=cache$node.E[3], D_N0=cache$node.D[1], D_N1=cache$node.D[2], D_N2=cache$node.D[3])
times=c(cache$tipward.age, cache$rootward.age)
runSilent <- function() {
options(warn = -1)
on.exit(options(warn = 0))
capture.output(res <- lsoda(yini, times, func = "classe_geosse_equivalent_derivs", padded.pars, initfunc="initmod_geosse", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "classe_geosse_equivalent_derivs", padded.pars, initfunc="initmod_geosse", dll = "canonical_geosse-ext-derivs", rtol=1e-8, atol=1e-8))
res
}
prob.subtree.cal.full <- runSilent()
}else{
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A, cache$d0A_0B, cache$d0A_0C, cache$d0A_0D, cache$d0A_0E, cache$d1A_1B, cache$d1A_1C, cache$d1A_1D, cache$d1A_1E, cache$d01A_01B, cache$d01A_01C, cache$d01A_01D, cache$d01A_01E, cache$s0B, cache$s1B, cache$s01B, cache$x0B, cache$x1B, cache$d0B_1B , cache$d0B_01B, cache$d1B_0B, cache$d1B_01B, cache$d01B_0B, cache$d01B_1B, cache$d0B_0A, cache$d0B_0C, cache$d0B_0D, cache$d0B_0E, cache$d1B_1A, cache$d1B_1C, cache$d1B_1D, cache$d1B_1E, cache$d01B_01A, cache$d01B_01C, cache$d01B_01D, cache$d01B_01E, cache$s0C, cache$s1C, cache$s01C, cache$x0C, cache$x1C, cache$d0C_1C , cache$d0C_01C, cache$d1C_0C, cache$d1C_01C, cache$d01C_0C, cache$d01C_1C, cache$d0C_0A, cache$d0C_0B, cache$d0C_0D, cache$d0C_0E, cache$d1C_1A, cache$d1C_1B, cache$d1C_1D, cache$d1C_1E, cache$d01C_01A, cache$d01C_01B, cache$d01C_01D, cache$d01C_01E, cache$s0D, cache$s1D, cache$s01D, cache$x0D, cache$x1D, cache$d0D_1D , cache$d0D_01D, cache$d1D_0D, cache$d1D_01D, cache$d01D_0D, cache$d01D_1D, cache$d0D_0A, cache$d0D_0B, cache$d0D_0C, cache$d0D_0E, cache$d1D_1A, cache$d1D_1B, cache$d1D_1C, cache$d1D_1E, cache$d01D_01A, cache$d01D_01B, cache$d01D_01C, cache$d01D_01E, cache$s0E, cache$s1E, cache$s01E, cache$x0E, cache$x1E, cache$d0E_1E , cache$d0E_01E, cache$d1E_0E, cache$d1E_01E, cache$d01E_0E, cache$d01E_1E, cache$d0E_0A, cache$d0E_0B, cache$d0E_0C, cache$d0E_0D, cache$d1E_1A, cache$d1E_1B, cache$d1E_1C, cache$d1E_1D, cache$d01E_01A, cache$d01E_01B, cache$d01E_01C, cache$d01E_01D)
NUMELEMENTS <- 115 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <- c(E_0A=cache$node.E[1], E_1A=cache$node.E[2], E_01A=cache$node.E[3], E_0B=cache$node.E[4], E_1B=cache$node.E[5], E_01B=cache$node.E[6], E_0C=cache$node.E[7], E_1C=cache$node.E[8], E_01C=cache$node.E[9], E_0D=cache$node.E[10], E_1D=cache$node.E[11], E_01D=cache$node.E[12], E_0E=cache$node.E[13], E_1E=cache$node.E[14], E_01E=cache$node.E[15], D_N0A=cache$node.D[1], D_N1A=cache$node.D[2], D_N01A=cache$node.D[3], D_N0B=cache$node.D[4], D_N1B=cache$node.D[5], D_N01B=cache$node.D[6], D_N0C=cache$node.D[7], D_N1C=cache$node.D[8], D_N01C=cache$node.D[9], D_N0D=cache$node.D[10], D_N1D=cache$node.D[11], D_N01D=cache$node.D[12], D_N0E=cache$node.D[13], D_N1E=cache$node.D[14], D_N01E=cache$node.D[15])
times=c(cache$tipward.age, cache$rootward.age)
runSilent <- function() {
options(warn = -1)
on.exit(options(warn = 0))
capture.output(res <- lsoda(yini, times, func = "geohisse_derivs", padded.pars, initfunc="initmod_geohisse", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "geohisse_derivs", padded.pars, initfunc="initmod_geohisse", dll = "geohisse-ext-derivs", rtol=1e-8, atol=1e-8))
res
}
prob.subtree.cal.full <- runSilent()
}
######## THIS CHECKS TO ENSURE THAT THE INTEGRATION WAS SUCCESSFUL ###########
if(attributes(prob.subtree.cal.full)$istate[1] < 0){
return(bad.likelihood)
}else{
prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
}
##############################################################################
if(hidden.states == FALSE){
if(is.nan(prob.subtree.cal[3]) | is.nan(prob.subtree.cal[4]) | is.nan(prob.subtree.cal[5])){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(prob.subtree.cal[4]<0 | prob.subtree.cal[5]<0 | prob.subtree.cal[6]<0){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[4:6]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:3])
v <- rbind(v, prob.subtree.cal[4:6])
}else{
if(any(is.nan(prob.subtree.cal[16:30]))){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(any(prob.subtree.cal[16:30] < 0)){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[16:30]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:15])
v <- rbind(v, prob.subtree.cal[16:30])
}
}
if(hidden.states == TRUE){
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0A + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1A + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01A
v <- v[,-c(1:3)]
compD[focal,4] <- v[1,1] * v[2,1] * cache$s0B
compD[focal,5] <- v[1,2] * v[2,2] * cache$s1B
compD[focal,6] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0B + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1B + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01B
v <- v[,-c(1:3)]
compD[focal,7] <- v[1,1] * v[2,1] * cache$s0C
compD[focal,8] <- v[1,2] * v[2,2] * cache$s1C
compD[focal,9] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0C + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1C + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01C
v <- v[,-c(1:3)]
compD[focal,10] <- v[1,1] * v[2,1] * cache$s0D
compD[focal,11] <- v[1,2] * v[2,2] * cache$s1D
compD[focal,12] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0D + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1D + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01D
v <- v[,-c(1:3)]
compD[focal,13] <- v[1,1] * v[2,1] * cache$s0E
compD[focal,14] <- v[1,2] * v[2,2] * cache$s1E
compD[focal,15] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0E + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1E + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01E
compE[focal,] <- phi[1:15]
if(!is.null(node)){
fixer = numeric(15)
fixer[state] = 1
if(node == focal){
compD[focal,] <- compD[focal,] * fixer
}
#compE[focal,] <- compE[focal,] * fixer
}
}else{
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0A + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1A + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01A
compE[focal,] <- phi[1:3]
if(!is.null(node)){
if(node == focal){
fixer = c(0,0,0)
fixer[state] = 1
compD[focal,] <- compD[focal,] * fixer
}
}
}
###########################
#Logcompensation bit for dealing with underflow issues. Need to give a necessary shoutout to Rich FitzJohn -- we follow his diversitree approach. VERIFIED that it works properly:
tmp <- sum(compD[focal,])
compD[focal,] <- compD[focal,] / tmp
logcomp <- c(logcomp, log(tmp))
}
root.node <- nb.tip + 1L
if (is.na(sum(log(compD[root.node,]))) || is.na(log(sum(1-compE[root.node,])))){
return(bad.likelihood)
}else{
if(root.type == "madfitz" | root.type == "herr_als"){
if(is.null(root.p)){
root.p = compD[root.node,]/sum(compD[root.node,])
root.p[which(is.na(root.p))] = 0
}
}
if(condition.on.survival == TRUE){
if(hidden.states == FALSE){
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)))
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)))
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}else{
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)), cache$s0B, cache$s1B, sum(c(cache$s0B, cache$s1B, cache$s01B)), cache$s0C, cache$s1C, sum(c(cache$s0C, cache$s1C, cache$s01C)), cache$s0D, cache$s1D, sum(c(cache$s0D, cache$s1D, cache$s01D)), cache$s0E, cache$s1E, sum(c(cache$s0E, cache$s1E, cache$s01E)))
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)), cache$s0B, cache$s1B, sum(c(cache$s0B, cache$s1B, cache$s01B)), cache$s0C, cache$s1C, sum(c(cache$s0C, cache$s1C, cache$s01C)), cache$s0D, cache$s1D, sum(c(cache$s0D, cache$s1D, cache$s01D)), cache$s0E, cache$s1E, sum(c(cache$s0E, cache$s1E, cache$s01E)))
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}
}
}
if(get.phi==TRUE){
obj = NULL
obj$compD.root = compD[root.node,]/sum(compD[root.node,])
obj$compE = compE
obj$root.p = root.p
return(obj)
}else{
return(loglik)
}
}
######################################################################################################################################
######################################################################################################################################
### The MuSSE type down pass that carries out the integration and returns the likelihood:
######################################################################################################################################
######################################################################################################################################
DownPassMusse <- function(phy, cache, hidden.states, bad.likelihood=-10000000, condition.on.survival, root.type, root.p, get.phi=FALSE, node=NULL, state=NULL, ode.eps=0) {
#Some preliminaries:
nb.tip <- length(phy$tip.label)
nb.node <- phy$Nnode
phy <- reorder(phy, "pruningwise")
anc <- unique(phy$edge[,1])
TIPS <- 1:nb.tip
if(hidden.states == FALSE){
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
}else{
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
if(!is.null(root.p)){
root.p.new <- numeric(15)
root.p.new[1:length(root.p)] <- root.p
root.p <- root.p.new
}
}
#Initializes the tip sampling and sets internal nodes to be zero:
ncols = dim(compD)[2]
if(length(cache$f) == 3){
for(i in 1:(nb.tip)){
compD[i,] <- cache$f * cache$states[i,]
compE[i,] <- rep((1-cache$f), ncols/3)
}
}else{
for(i in 1:(nb.tip)){
compD[i,] <- cache$f[i] * cache$states[i,]
compE[i,] <- rep((1-cache$f[i]), ncols/3)
}
}
logcomp <- c()
#Start the postorder traversal indexing lists by node number:
for (i in seq(from = 1, length.out = nb.node)) {
#A vector of all the internal nodes:
focal <- anc[i]
desRows <- which(phy$edge[,1]==focal)
desNodes <- phy$edge[desRows,2]
#Note: when the tree has been reordered branching.times are no longer valid. Fortunately, we extract this information in the initial cache setup. Also focal is the rootward node, whereas desNodes represent a vector of all descendant nodes:
cache$rootward.age <- cache$split.times[which(names(cache$split.times)==focal)]
v <- c()
phi <- c()
for (desIndex in sequence(length(desRows))){
cache$focal.edge.length <- phy$edge.length[desRows[desIndex]]
cache$tipward.age <- cache$rootward.age - cache$focal.edge.length
#Strange rounding errors. A tip age should be zero. This ensures that:
if(cache$tipward.age < .Machine$double.eps^0.5){
cache$tipward.age = 0
}
cache$node.D <- compD[desNodes[desIndex],]
cache$node.E <- compE[desNodes[desIndex],]
##Call to lsoda that utilizes C code. Requires a lot of inputs. Note that for now we hardcode the NUMELEMENTS arguments. The reason for this is because with lsoda we can only pass a vector of parameters.
if(hidden.states == FALSE){
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$x01A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A)
NUMELEMENTS <- 12 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <-c(E_0=cache$node.E[1], E_1=cache$node.E[2], E_01=cache$node.E[3], D_N0=cache$node.D[1], D_N1=cache$node.D[2], D_N2=cache$node.D[3])
times=c(cache$tipward.age, cache$rootward.age)
#runSilent <- function() {
#options(warn = -1)
#on.exit(options(warn = 0))
#capture.output(res <- lsoda(yini, times, func = "notclasse_derivs", padded.pars, initfunc="initmod_musse", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "notclasse_derivs", padded.pars, initfunc="initmod_musse", dll = "notclasse-ext-derivs", rtol=1e-8, atol=1e-8))
#res
#}
#prob.subtree.cal.full <- runSilent()
prob.subtree.cal.full <- lsoda(yini, times, func = "notclasse_derivs", padded.pars, initfunc="initmod_noclass", dllname = "hisse", rtol=1e-8, atol=1e-8)
}else{
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$x01A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A, cache$d0A_0B, cache$d0A_0C, cache$d0A_0D, cache$d0A_0E, cache$d1A_1B, cache$d1A_1C, cache$d1A_1D, cache$d1A_1E, cache$d01A_01B, cache$d01A_01C, cache$d01A_01D, cache$d01A_01E, cache$s0B, cache$s1B, cache$s01B, cache$x0B, cache$x1B, cache$x01B, cache$d0B_1B , cache$d0B_01B, cache$d1B_0B, cache$d1B_01B, cache$d01B_0B, cache$d01B_1B, cache$d0B_0A, cache$d0B_0C, cache$d0B_0D, cache$d0B_0E, cache$d1B_1A, cache$d1B_1C, cache$d1B_1D, cache$d1B_1E, cache$d01B_01A, cache$d01B_01C, cache$d01B_01D, cache$d01B_01E, cache$s0C, cache$s1C, cache$s01C, cache$x0C, cache$x1C, cache$x01C, cache$d0C_1C , cache$d0C_01C, cache$d1C_0C, cache$d1C_01C, cache$d01C_0C, cache$d01C_1C, cache$d0C_0A, cache$d0C_0B, cache$d0C_0D, cache$d0C_0E, cache$d1C_1A, cache$d1C_1B, cache$d1C_1D, cache$d1C_1E, cache$d01C_01A, cache$d01C_01B, cache$d01C_01D, cache$d01C_01E, cache$s0D, cache$s1D, cache$s01D, cache$x0D, cache$x1D, cache$x01D, cache$d0D_1D , cache$d0D_01D, cache$d1D_0D, cache$d1D_01D, cache$d01D_0D, cache$d01D_1D, cache$d0D_0A, cache$d0D_0B, cache$d0D_0C, cache$d0D_0E, cache$d1D_1A, cache$d1D_1B, cache$d1D_1C, cache$d1D_1E, cache$d01D_01A, cache$d01D_01B, cache$d01D_01C, cache$d01D_01E, cache$s0E, cache$s1E, cache$s01E, cache$x0E, cache$x1E, cache$x01E, cache$d0E_1E, cache$d0E_01E, cache$d1E_0E, cache$d1E_01E, cache$d01E_0E, cache$d01E_1E, cache$d0E_0A, cache$d0E_0B, cache$d0E_0C, cache$d0E_0D, cache$d1E_1A, cache$d1E_1B, cache$d1E_1C, cache$d1E_1D, cache$d01E_01A, cache$d01E_01B, cache$d01E_01C, cache$d01E_01D)
NUMELEMENTS <- 120 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <- c(E_0A=cache$node.E[1], E_1A=cache$node.E[2], E_01A=cache$node.E[3], E_0B=cache$node.E[4], E_1B=cache$node.E[5], E_01B=cache$node.E[6], E_0C=cache$node.E[7], E_1C=cache$node.E[8], E_01C=cache$node.E[9], E_0D=cache$node.E[10], E_1D=cache$node.E[11], E_01D=cache$node.E[12], E_0E=cache$node.E[13], E_1E=cache$node.E[14], E_01E=cache$node.E[15], D_N0A=cache$node.D[1], D_N1A=cache$node.D[2], D_N01A=cache$node.D[3], D_N0B=cache$node.D[4], D_N1B=cache$node.D[5], D_N01B=cache$node.D[6], D_N0C=cache$node.D[7], D_N1C=cache$node.D[8], D_N01C=cache$node.D[9], D_N0D=cache$node.D[10], D_N1D=cache$node.D[11], D_N01D=cache$node.D[12], D_N0E=cache$node.D[13], D_N1E=cache$node.D[14], D_N01E=cache$node.D[15])
times=c(cache$tipward.age, cache$rootward.age)
#runSilent <- function() {
# options(warn = -1)
# on.exit(options(warn = 0))
# capture.output(res <- lsoda(yini, times, func = "notclasse_more_derivs", padded.pars, initfunc="initmod_mussem", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "notclasse_more_derivs", padded.pars, initfunc="initmod_mussem", dll = "notclasse-more-ext-derivs", rtol=1e-8, atol=1e-8))
# res
#}
#prob.subtree.cal.full <- runSilent()
prob.subtree.cal.full <- lsoda(yini, times, func = "notclasse_more_derivs", padded.pars, initfunc="initmod_hinoclass", dllname = "hisse", rtol=1e-8, atol=1e-8)
}
######## THIS CHECKS TO ENSURE THAT THE INTEGRATION WAS SUCCESSFUL ###########
if(attributes(prob.subtree.cal.full)$istate[1] < 0){
return(bad.likelihood)
}else{
prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
}
##############################################################################
if(hidden.states == FALSE){
if(is.nan(prob.subtree.cal[3]) | is.nan(prob.subtree.cal[4]) | is.nan(prob.subtree.cal[5])){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(prob.subtree.cal[4]<0 | prob.subtree.cal[5]<0 | prob.subtree.cal[6]<0){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[4:6]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:3])
v <- rbind(v, prob.subtree.cal[4:6])
}else{
if(any(is.nan(prob.subtree.cal[16:30]))){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(any(prob.subtree.cal[16:30] < 0)){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[10:18]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:15])
v <- rbind(v, prob.subtree.cal[16:30])
}
}
##Allows for a 3-state MuSSE type model. Important for cases where you have a random tree and a random trait.
if(hidden.states == TRUE){
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- v[1,3] * v[2,3] * cache$s01A
v <- v[,-c(1:3)]
compD[focal,4] <- v[1,1] * v[2,1] * cache$s0B
compD[focal,5] <- v[1,2] * v[2,2] * cache$s1B
compD[focal,6] <- v[1,3] * v[2,3] * cache$s01B
v <- v[,-c(1:3)]
compD[focal,7] <- v[1,1] * v[2,1] * cache$s0C
compD[focal,8] <- v[1,2] * v[2,2] * cache$s1C
compD[focal,9] <- v[1,3] * v[2,3] * cache$s01C
v <- v[,-c(1:3)]
compD[focal,10] <- v[1,1] * v[2,1] * cache$s0D
compD[focal,11] <- v[1,2] * v[2,2] * cache$s1D
compD[focal,12] <- v[1,3] * v[2,3] * cache$s01D
v <- v[,-c(1:3)]
compD[focal,13] <- v[1,1] * v[2,1] * cache$s0E
compD[focal,14] <- v[1,2] * v[2,2] * cache$s1E
compD[focal,15] <- v[1,3] * v[2,3] * cache$s01E
compE[focal,] <- phi[1:15]
if(!is.null(node)){
fixer = numeric(15)
fixer[state] = 1
if(node == focal){
compD[focal,] <- compD[focal,] * fixer
}
#compE[focal,] <- compE[focal,] * fixer
}
}else{
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- v[1,3] * v[2,3] * cache$s01A
compE[focal,] <- phi[1:3]
if(!is.null(node)){
if(node == focal){
fixer = c(0,0,0)
fixer[state] = 1
compD[focal,] <- compD[focal,] * fixer
}
}
}
###########################
#Logcompensation bit for dealing with underflow issues. Need to give a necessary shoutout to Rich FitzJohn -- we follow his diversitree approach. VERIFIED that it works properly:
tmp <- sum(compD[focal,])
compD[focal,] <- compD[focal,] / tmp
logcomp <- c(logcomp, log(tmp))
}
root.node <- nb.tip + 1L
if (is.na(sum(log(compD[root.node,]))) || is.na(log(sum(1-compE[root.node,])))){
return(bad.likelihood)
}else{
if(root.type == "madfitz" | root.type == "herr_als"){
if(is.null(root.p)){
root.p = compD[root.node,]/sum(compD[root.node,])
root.p[which(is.na(root.p))] = 0
}
}
if(condition.on.survival == TRUE){
if(hidden.states == FALSE){
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, cache$s01A)
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, cache$s01A)
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}else{
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, cache$s01A, cache$s0B, cache$s1B, cache$s01B, cache$s0C, cache$s1C, cache$s01C, cache$s0D, cache$s1D, cache$s01D, cache$s0E, cache$s1E, cache$s01E)
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, cache$s01A, cache$s0B, cache$s1B, cache$s01B, cache$s0C, cache$s1C, cache$s01C, cache$s0D, cache$s1D, cache$s01D, cache$s0E, cache$s1E, cache$s01E)
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}
}
}
if(get.phi==TRUE){
obj = NULL
obj$compD.root = compD[root.node,]/sum(compD[root.node,])
obj$compE = compE
obj$root.p = root.p
return(obj)
}else{
return(loglik)
}
}
######################################################################################################################################
######################################################################################################################################
### Cache object for storing parameters that are used throughout GeoHiSSE:
######################################################################################################################################
######################################################################################################################################
ParametersToPassGeoHiSSE <- function(phy, data, f, model.vec, hidden.states){
#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$phy = phy
if(hidden.states == FALSE){
states = matrix(0,Ntip(phy),3)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==0){states[i,3]=1}
}
}
if(hidden.states == TRUE){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,c(1,4,7,10,13)]=1}
if(data[i]==2){states[i,c(2,5,8,11,14)]=1}
if(data[i]==0){states[i,c(3,6,9,12,15)]=1}
}
}
if(hidden.states == "TEST"){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==0){states[i,3]=1}
if(data[i]==4){states[i,4]=1}
if(data[i]==5){states[i,5]=1}
if(data[i]==3){states[i,6]=1}
}
}
obj$states = states
obj$tot_time = max(branching.times(phy))
obj$f = f
obj$s0A = model.vec[1]
obj$s1A = model.vec[2]
obj$s01A = model.vec[3]
obj$x0A = model.vec[4]
obj$x1A = model.vec[5]
obj$d0A_1A = model.vec[6]
obj$d0A_01A = model.vec[7]
obj$d1A_0A = model.vec[8]
obj$d1A_01A = model.vec[9]
if(model.vec[10]==0){
obj$d01A_0A = model.vec[5]
}else{
obj$d01A_0A = model.vec[10]
}
if(model.vec[11]==0){
obj$d01A_1A = model.vec[4]
}else{
obj$d01A_1A = model.vec[11]
}
obj$d0A_0B = model.vec[12]
obj$d0A_0C = model.vec[13]
obj$d0A_0D = model.vec[14]
obj$d0A_0E = model.vec[15]
obj$d1A_1B = model.vec[16]
obj$d1A_1C = model.vec[17]
obj$d1A_1D = model.vec[18]
obj$d1A_1E = model.vec[19]
obj$d01A_01B = model.vec[20]
obj$d01A_01C = model.vec[21]
obj$d01A_01D = model.vec[22]
obj$d01A_01E = model.vec[23]
obj$s0B = model.vec[24]
obj$s1B = model.vec[25]
obj$s01B = model.vec[26]
obj$x0B = model.vec[27]
obj$x1B = model.vec[28]
obj$d0B_1B = model.vec[29]
obj$d0B_01B = model.vec[30]
obj$d1B_0B = model.vec[31]
obj$d1B_01B = model.vec[32]
if(model.vec[33] == 0){
obj$d01B_0B = model.vec[28]
}else{
obj$d01B_0B = model.vec[33]
}
if(model.vec[34] == 0){
obj$d01B_1B = model.vec[27]
}else{
obj$d01B_1B = model.vec[34]
}
obj$d0B_0A = model.vec[35]
obj$d0B_0C = model.vec[36]
obj$d0B_0D = model.vec[37]
obj$d0B_0E = model.vec[38]
obj$d1B_1A = model.vec[39]
obj$d1B_1C = model.vec[40]
obj$d1B_1D = model.vec[41]
obj$d1B_1E = model.vec[42]
obj$d01B_01A = model.vec[43]
obj$d01B_01C = model.vec[44]
obj$d01B_01D = model.vec[45]
obj$d01B_01E = model.vec[46]
obj$s0C = model.vec[47]
obj$s1C = model.vec[48]
obj$s01C = model.vec[49]
obj$x0C = model.vec[50]
obj$x1C = model.vec[51]
obj$d0C_1C = model.vec[52]
obj$d0C_01C = model.vec[53]
obj$d1C_0C = model.vec[54]
obj$d1C_01C = model.vec[55]
if(model.vec[56] == 0){
obj$d01C_0C = model.vec[51]
}else{
obj$d01C_0C = model.vec[56]
}
if(model.vec[57] == 0){
obj$d01C_1C = model.vec[50]
}else{
obj$d01C_1C = model.vec[57]
}
obj$d0C_0A = model.vec[58]
obj$d0C_0B = model.vec[59]
obj$d0C_0D = model.vec[60]
obj$d0C_0E = model.vec[61]
obj$d1C_1A = model.vec[62]
obj$d1C_1B = model.vec[63]
obj$d1C_1D = model.vec[64]
obj$d1C_1E = model.vec[65]
obj$d01C_01A = model.vec[66]
obj$d01C_01B = model.vec[67]
obj$d01C_01D = model.vec[68]
obj$d01C_01E = model.vec[69]
obj$s0D = model.vec[70]
obj$s1D = model.vec[71]
obj$s01D = model.vec[72]
obj$x0D = model.vec[73]
obj$x1D = model.vec[74]
obj$d0D_1D = model.vec[75]
obj$d0D_01D = model.vec[76]
obj$d1D_0D = model.vec[77]
obj$d1D_01D = model.vec[78]
if(model.vec[79] == 0){
obj$d01D_0D = model.vec[74]
}else{
obj$d01D_0D = model.vec[79]
}
if(model.vec[80] == 0){
obj$d01D_1D = model.vec[73]
}else{
obj$d01D_1D = model.vec[80]
}
obj$d0D_0A = model.vec[81]
obj$d0D_0B = model.vec[82]
obj$d0D_0C = model.vec[83]
obj$d0D_0E = model.vec[84]
obj$d1D_1A = model.vec[85]
obj$d1D_1B = model.vec[86]
obj$d1D_1C = model.vec[87]
obj$d1D_1E = model.vec[88]
obj$d01D_01A = model.vec[89]
obj$d01D_01B = model.vec[90]
obj$d01D_01C = model.vec[91]
obj$d01D_01E = model.vec[92]
obj$s0E = model.vec[93]
obj$s1E = model.vec[94]
obj$s01E = model.vec[95]
obj$x0E = model.vec[96]
obj$x1E = model.vec[97]
obj$d0E_1E = model.vec[98]
obj$d0E_01E = model.vec[99]
obj$d1E_0E = model.vec[100]
obj$d1E_01E = model.vec[101]
if(model.vec[102] == 0){
obj$d01E_0E = model.vec[97]
}else{
obj$d01E_0E = model.vec[102]
}
if(model.vec[103] == 0){
obj$d01E_1E = model.vec[96]
}else{
obj$d01E_1E = model.vec[103]
}
obj$d0E_0A = model.vec[104]
obj$d0E_0B = model.vec[105]
obj$d0E_0C = model.vec[106]
obj$d0E_0D = model.vec[107]
obj$d1E_1A = model.vec[108]
obj$d1E_1B = model.vec[109]
obj$d1E_1C = model.vec[110]
obj$d1E_1D = model.vec[111]
obj$d01E_01A = model.vec[112]
obj$d01E_01B = model.vec[113]
obj$d01E_01C = model.vec[114]
obj$d01E_01D = model.vec[115]
obj$split.times = sort(branching.times(phy), decreasing=TRUE)
return(obj)
}
######################################################################################################################################
######################################################################################################################################
### Cache object for storing parameters that are used throughout MuSSE:
######################################################################################################################################
######################################################################################################################################
ParametersToPassMuSSE <- function(phy, data, f, model.vec, hidden.states){
#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$phy = phy
if(hidden.states == FALSE){
states = matrix(0,Ntip(phy),3)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==0){states[i,3]=1}
}
}
if(hidden.states == TRUE){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,c(1,4,7,10,13)]=1}
if(data[i]==2){states[i,c(2,5,8,11,14)]=1}
if(data[i]==0){states[i,c(3,6,9,12,15)]=1}
}
}
if(hidden.states == "TEST1"){
states = matrix(0,Ntip(phy),3)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==4){states[i,3]=1}
}
}
if(hidden.states == "TEST2"){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,4]=1}
if(data[i]==2){states[i,5]=1}
if(data[i]==4){states[i,6]=1}
}
}
obj$states = states
obj$tot_time = max(branching.times(phy))
obj$f = f
obj$s0A = model.vec[1]
obj$s1A = model.vec[2]
obj$s01A = model.vec[3]
obj$x0A = model.vec[4]
obj$x1A = model.vec[5]
obj$x01A = model.vec[6]
obj$d0A_1A = model.vec[7]
obj$d0A_01A = model.vec[8]
obj$d1A_0A = model.vec[9]
obj$d1A_01A = model.vec[10]
if(model.vec[11] == 0){
obj$d01A_0A = model.vec[5]
}else{
obj$d01A_0A = model.vec[11]
}
if(model.vec[12] == 0){
obj$d01A_1A = model.vec[4]
}else{
obj$d01A_1A = model.vec[12]
}
obj$d0A_0B = model.vec[13]
obj$d0A_0C = model.vec[14]
obj$d0A_0D = model.vec[15]
obj$d0A_0E = model.vec[16]
obj$d1A_1B = model.vec[17]
obj$d1A_1C = model.vec[18]
obj$d1A_1D = model.vec[19]
obj$d1A_1E = model.vec[20]
obj$d01A_01B = model.vec[21]
obj$d01A_01C = model.vec[22]
obj$d01A_01D = model.vec[23]
obj$d01A_01E = model.vec[24]
obj$s0B = model.vec[25]
obj$s1B = model.vec[26]
obj$s01B = model.vec[27]
obj$x0B = model.vec[28]
obj$x1B = model.vec[29]
obj$x01B = model.vec[30]
obj$d0B_1B = model.vec[31]
obj$d0B_01B = model.vec[32]
obj$d1B_0B = model.vec[33]
obj$d1B_01B = model.vec[34]
if(model.vec[35] == 0){
obj$d01B_0B = model.vec[29]
}else{
obj$d01B_0B = model.vec[35]
}
if(model.vec[36] == 0){
obj$d01B_1B = model.vec[28]
}else{
obj$d01B_1B = model.vec[36]
}
obj$d0B_0A = model.vec[37]
obj$d0B_0C = model.vec[38]
obj$d0B_0D = model.vec[39]
obj$d0B_0E = model.vec[40]
obj$d1B_1A = model.vec[41]
obj$d1B_1C = model.vec[42]
obj$d1B_1D = model.vec[43]
obj$d1B_1E = model.vec[44]
obj$d01B_01A = model.vec[45]
obj$d01B_01C = model.vec[46]
obj$d01B_01D = model.vec[47]
obj$d01B_01E = model.vec[48]
obj$s0C = model.vec[49]
obj$s1C = model.vec[50]
obj$s01C = model.vec[51]
obj$x0C = model.vec[52]
obj$x1C = model.vec[53]
obj$x01C = model.vec[54]
obj$d0C_1C = model.vec[55]
obj$d0C_01C = model.vec[56]
obj$d1C_0C = model.vec[57]
obj$d1C_01C = model.vec[58]
if(model.vec[59] == 0){
obj$d01C_0C = model.vec[53]
}else{
obj$d01C_0C = model.vec[59]
}
if(model.vec[60] == 0){
obj$d01C_1C = model.vec[52]
}else{
obj$d01C_1C = model.vec[60]
}
obj$d0C_0A = model.vec[61]
obj$d0C_0B = model.vec[62]
obj$d0C_0D = model.vec[63]
obj$d0C_0E = model.vec[64]
obj$d1C_1A = model.vec[65]
obj$d1C_1B = model.vec[66]
obj$d1C_1D = model.vec[67]
obj$d1C_1E = model.vec[68]
obj$d01C_01A = model.vec[69]
obj$d01C_01B = model.vec[70]
obj$d01C_01D = model.vec[71]
obj$d01C_01E = model.vec[72]
obj$s0D = model.vec[73]
obj$s1D = model.vec[74]
obj$s01D = model.vec[75]
obj$x0D = model.vec[76]
obj$x1D = model.vec[77]
obj$x01D = model.vec[78]
obj$d0D_1D = model.vec[79]
obj$d0D_01D = model.vec[80]
obj$d1D_0D = model.vec[81]
obj$d1D_01D = model.vec[82]
if(model.vec[83] == 0){
obj$d01D_0D = model.vec[77]
}else{
obj$d01D_0D = model.vec[83]
}
if(model.vec[84] == 0){
obj$d01D_1D = model.vec[76]
}else{
obj$d01D_1D = model.vec[84]
}
obj$d0D_0A = model.vec[85]
obj$d0D_0B = model.vec[86]
obj$d0D_0C = model.vec[87]
obj$d0D_0E = model.vec[88]
obj$d1D_1A = model.vec[89]
obj$d1D_1B = model.vec[90]
obj$d1D_1C = model.vec[91]
obj$d1D_1E = model.vec[92]
obj$d01D_01A = model.vec[93]
obj$d01D_01B = model.vec[94]
obj$d01D_01C = model.vec[95]
obj$d01D_01E = model.vec[96]
obj$s0E = model.vec[97]
obj$s1E = model.vec[98]
obj$s01E = model.vec[99]
obj$x0E = model.vec[100]
obj$x1E = model.vec[101]
obj$x01E = model.vec[102]
obj$d0E_1E = model.vec[103]
obj$d0E_01E = model.vec[104]
obj$d1E_0E = model.vec[105]
obj$d1E_01E = model.vec[106]
if(model.vec[107] == 0){
obj$d01E_0E = model.vec[101]
}else{
obj$d01E_0E = model.vec[107]
}
if(model.vec[108] == 0){
obj$d01E_1E = model.vec[100]
}else{
obj$d01E_1E = model.vec[108]
}
obj$d0E_0A = model.vec[109]
obj$d0E_0B = model.vec[110]
obj$d0E_0C = model.vec[111]
obj$d0E_0D = model.vec[112]
obj$d1E_1A = model.vec[113]
obj$d1E_1B = model.vec[114]
obj$d1E_1C = model.vec[115]
obj$d1E_1D = model.vec[116]
obj$d01E_01A = model.vec[117]
obj$d01E_01B = model.vec[118]
obj$d01E_01C = model.vec[119]
obj$d01E_01D = model.vec[120]
obj$split.times = sort(branching.times(phy), decreasing=TRUE)
return(obj)
}
######################################################################################################################################
######################################################################################################################################
### Print function for our diversity class:
######################################################################################################################################
######################################################################################################################################
print.geohisse.old.fit <- function(x,...){
## Function to print a "geohisse.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 )
nareas <- ncol( x$trans.matrix )/3
output <- c(x$loglik, x$AIC, x$AICc, ntips, nareas)
names(output) <- c("lnL", "AIC", "AICc", "n.taxa", "n.hidden.classes")
cat("\n")
cat("Fit \n")
print(output)
cat("\n")
cat("Model parameters: \n")
cat("\n")
print(par.list)
cat("\n")
}
######################################################################################################################################
######################################################################################################################################
### Code for testing purposes:
######################################################################################################################################
######################################################################################################################################
#pars <- c(1.5, 0.5, 1.0, 0.7, 0.7, 2.5, 0.5)
#names(pars) <- diversitree:::default.argnames.geosse()
## Simulate a tree
#set.seed(5)
#phy <- tree.geosse(pars, max.t=4, x0=0)
#data <- data.frame(g_s = names(phy$tip.state), states = phy$tip.state)
## See the data
#statecols <- c("AB"="purple", "A"="blue", "B"="red")
## The likelihood function
#lik <- make.geosse(phy, phy$tip.state)
#lik(pars)
#states <- data.frame(phy$tip.state, phy$tip.state, row.names=names(phy$tip.state))
#states <- states[phy$tip.label,]
#names(pars) <- NULL
#model.vec <- numeric(95)
#tests against GeoSSE
#model.vec <- rep(c(pars[1:3], pars[4:5], pars[6:7], rep(0,12)), 5)
#model.vec <- numeric(95)
#model.vec[1:7] <- c(pars[1:3], pars[4:5], pars[6:7])
#phy$node.label <- NULL
#cache <- ParametersToPassGeoHiSSE(phy, states[,1], f=c(1,1,1), model.vec, hidden.states=FALSE)
#ll.geosse <- DownPassGeoHisse(phy=phy, cache=cache, hidden.states=FALSE, bad.likelihood=-10000000000, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL)
#cache <- ParametersToPassGeoHiSSE(phy, states[,1], f=c(1,1,1), model.vec, hidden.states=TRUE)
#ll.geohisse <- DownPassGeoHisse(phy=phy, cache=cache, hidden.states=TRUE, bad.likelihood=-10000000000, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL)
#tests against NULL
#lambda_A <- pars[1]
#lambda_B <- pars[2]
#lambda_C <- pars[3]
#model.vec <- c(rep(lambda_A, 3), pars[4:5], pars[6:7], rep(1,6), rep(lambda_B, 3), pars[4:5], pars[6:7], rep(lambda_C, 3), pars[1:3], pars[4:5], pars[6:7], rep(1,6))
#phy$node.label <- NULL
#cache <- ParametersToPassGeoSSE(phy, states[,1], f=c(1,1,1), model.vec, hidden.states=TRUE)
#ll <- DownPassGeosse(phy=phy, cache=cache, hidden.states=TRUE, bad.likelihood=-10000000000, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL)
thej022214/hisse documentation built on Sept. 20, 2023, 12:40 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.geohisse.old.fit ParametersToPassMuSSE ParametersToPassGeoHiSSE DownPassMusse DownPassGeoHisse starting.point.geosse DevOptimizeGeoHiSSE GeoHiSSE.old
Documented in GeoHiSSE.old
######################################################################################################################################
######################################################################################################################################
### GeoHiSSE -- Expanded set of GeoSSE models for examining diversification in relation to geographic range evolution
######################################################################################################################################
######################################################################################################################################
GeoHiSSE.old <- function(phy, data, f=c(1,1,1), speciation=c(1,2,3), extirpation=c(1,2), hidden.areas=FALSE, trans.rate=NULL, assume.cladogenetic=TRUE, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL, sann=TRUE, sann.its=1000, bounded.search=TRUE, max.tol=.Machine$double.eps^.50, mag.san.start=0.5, starting.vals=NULL, speciation.upper=1000, extirpation.upper=1000, trans.upper=100, ode.eps=0){
## Temporary fix for the current BUG:
if( !is.null(phy$node.label) ) phy$node.label <- NULL
if(!is.null(root.p)) {
if(hidden.areas ==TRUE){
if( length( root.p ) == 3 ){
root.p <- rep(root.p, 5)
root.p <- root.p / sum(root.p)
warning("For hidden areas, you need to specify the root.p for all hidden areas. We have adjusted it so that there's equal chance for among all hidden areas.")
} else{
root.p.new <- numeric(15)
root.p.new[1:length(root.p)] <- root.p
root.p <- root.p.new
root.p <- root.p / sum(root.p)
}
}else{
## All good:
root.p <- root.p / sum(root.p)
}
}
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(is.null(trans.rate)){
stop("Rate matrix needed. See TransMatMakerGeoHiSSE() to create one.")
}
if(hidden.areas == TRUE & dim(trans.rate)[1]<3){
stop("You chose a hidden state but this is not reflected in the transition matrix")
}
if(assume.cladogenetic == FALSE){
if(dim(trans.rate)[1] != length(extirpation)){
stop("You have not specified enough extinction parameters.")
}
}
## Return error message if the data is not in the correct format.
if( !inherits(data, what = c("matrix","data.frame")) ){
stop("'data' needs to be a matrix or data.frame with 2 columns. See help.")
}
if( !ncol( data ) == 2 ){
stop("'data' needs to be a matrix or data.frame with 2 columns. See help.")
}
## Check if the states are %in% 0:2:
states.check <- all( as.numeric(data[,2]) %in% 0:2 )
if( !states.check ){
stop("states need to be one of 0, 1, or 2. See help.")
}
## Check if 'hidden.areas' parameter is congruent with the speciation vector:
if( length(speciation) > 3 & !hidden.areas ){
stop("Speciation has more than 3 elements but 'hidden.areas' was set to FALSE. Please set 'hidden.areas' to TRUE if the model include more than one rate class.")
}
if( length(speciation) == 3 & hidden.areas ){
stop("Speciation has only 3 elements but 'hidden.areas' was set to TRUE. Please set 'hidden.areas' to FALSE if the model does not include hidden rate classes.")
}
if(assume.cladogenetic == TRUE){
pars <- numeric(115)
if(dim(trans.rate)[2]==3){
rate.cats <- 1
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
trans.tmp <- c(trans.rate["(0)", "(1)"], trans.rate["(0)", "(01)"], trans.rate["(1)", "(0)"], trans.rate["(1)", "(01)"], trans.rate["(01)", "(0)"], trans.rate["(01)", "(1)"])
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
category.rates.unique <- 0
pars.tmp <- c(pars.tmp, trans.tmp)
pars[1:11] <- pars.tmp
}
if(dim(trans.rate)[2]==6){
rate.cats <- 2
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 1), rep("(1A)", 1), rep("(01A)", 1), rep("(0B)", 1), rep("(1B)", 1), rep("(01B)", 1))
cols <- c("(0B)", "(1B)", "(01B)", "(0A)", "(1A)", "(01A)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- c(category.rate.shift[1], rep(0,3), category.rate.shift[2], rep(0,3), category.rate.shift[3], rep(0,3))
category.rate.shiftB <- c(category.rate.shift[4], rep(0,3), category.rate.shift[5], rep(0,3), category.rate.shift[6], rep(0,3))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==9){
rate.cats <- 3
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 2), rep("(1A)", 2), rep("(01A)", 2), rep("(0B)", 2), rep("(1B)", 2), rep("(01B)", 2), rep("(0C)", 2), rep("(1C)", 2), rep("(01C)", 2))
cols <- c("(0B)", "(0C)", "(1B)", "(1C)", "(01B)", "(01C)", "(0A)", "(0C)", "(1A)", "(1C)", "(01A)", "(01C)", "(0A)", "(0B)", "(1A)", "(1B)", "(01A)", "(01B)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- c(category.rate.shift[1:2], rep(0,2), category.rate.shift[3:4], rep(0,2), category.rate.shift[5:6], rep(0,2))
category.rate.shiftB <- c(category.rate.shift[7:8], rep(0,2), category.rate.shift[9:10], rep(0,2), category.rate.shift[11:12], rep(0,2))
category.rate.shiftC <- c(category.rate.shift[13:14], rep(0,2), category.rate.shift[15:16], rep(0,2), category.rate.shift[17:18], rep(0,2))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[5:6], trans.tmp[13:18], category.rate.shiftC)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==12){
rate.cats <- 4
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 3), rep("(1A)", 3), rep("(01A)", 3), rep("(0B)", 3), rep("(1B)", 3), rep("(01B)", 3), rep("(0C)", 3), rep("(1C)", 3), rep("(01C)", 3), rep("(0D)", 3), rep("(1D)", 3), rep("(01D)", 3))
cols <- c("(0B)", "(0C)", "(0D)", "(1B)", "(1C)", "(1D)", "(01B)", "(01C)", "(01D)", "(0A)", "(0C)", "(0D)", "(1A)", "(1C)", "(1D)", "(01A)", "(01C)", "(01D)", "(0A)", "(0B)", "(0D)", "(1A)", "(1B)", "(1D)", "(01A)", "(01B)", "(01D)", "(0A)", "(0B)", "(0C)", "(1A)", "(1B)", "(1C)", "(01A)", "(01B)", "(01C)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- c(category.rate.shift[1:3], rep(0,1), category.rate.shift[4:6], rep(0,1), category.rate.shift[7:9], rep(0,1))
category.rate.shiftB <- c(category.rate.shift[10:12], rep(0,1), category.rate.shift[13:15], rep(0,1), category.rate.shift[16:18], rep(0,1))
category.rate.shiftC <- c(category.rate.shift[19:21], rep(0,1), category.rate.shift[22:24], rep(0,1), category.rate.shift[25:27], rep(0,1))
category.rate.shiftD <- c(category.rate.shift[28:30], rep(0,1), category.rate.shift[31:33], rep(0,1), category.rate.shift[34:36], rep(0,1))
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[5:6], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[7:8], trans.tmp[19:24], category.rate.shiftD)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==15){
rate.cats <- 5
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
for.late.adjust <- max(pars.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)", "(0E)", "(0E)", "(1E)", "(1E)", "(01E)", "(01E)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)", "(1E)", "(01E)", "(0E)", "(01E)", "(0E)", "(1E)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
rows <- c(rep("(0A)", 4), rep("(1A)", 4), rep("(01A)", 4), rep("(0B)", 4), rep("(1B)", 4), rep("(01B)", 4), rep("(0C)", 4), rep("(1C)", 4), rep("(01C)", 4), rep("(0D)", 4), rep("(1D)", 4), rep("(01D)", 4), rep("(0E)", 4), rep("(1E)", 4), rep("(01E)", 4))
cols <- c("(0B)", "(0C)", "(0D)", "(0E)", "(1B)", "(1C)", "(1D)", "(1E)", "(01B)", "(01C)", "(01D)", "(01E)", "(0A)", "(0C)", "(0D)", "(0E)", "(1A)", "(1C)", "(1D)", "(1E)", "(01A)", "(01C)", "(01D)", "(01E)", "(0A)", "(0B)", "(0D)", "(0E)", "(1A)", "(1B)", "(1D)", "(1E)", "(01A)", "(01B)", "(01D)", "(01E)", "(0A)", "(0B)", "(0C)", "(0E)", "(1A)", "(1B)", "(1C)", "(1E)", "(01A)", "(01B)", "(01C)", "(01E)", "(0A)", "(0B)", "(0C)", "(0D)", "(1A)", "(1B)", "(1C)", "(1D)", "(01A)", "(01B)", "(01C)", "(01D)")
category.tmp <- trans.rate[cbind(rows,cols)]
category.rate.shift <- category.tmp + for.late.adjust
category.rate.shift[is.na(category.rate.shift)] <- 0
category.rate.shiftA <- category.rate.shift[1:12]
category.rate.shiftB <- category.rate.shift[13:24]
category.rate.shiftC <- category.rate.shift[25:36]
category.rate.shiftD <- category.rate.shift[37:48]
category.rate.shiftE <- category.rate.shift[49:60]
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD, category.rate.shiftE)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:2], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[3:4], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[5:6], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[7:8], trans.tmp[19:24], category.rate.shiftD, speciation[13:15], extirpation.tmp[9:10], trans.tmp[25:30], category.rate.shiftE)
pars[1:length(pars.tmp)] <- pars.tmp
}
np <- max(pars)
pars[pars==0] <- np+1
cat("Initializing...", "\n")
data.new <- data.frame(data[,2], data[,2], row.names=data[,1])
data.new <- data.new[phy$tip.label,]
#This is used to scale starting values to account for sampling:
if(length(f) == 3){
freqs <- table(data.new[,1])
if(length(freqs == 2)){
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}else{
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}
}else{
if(length(f) == Ntip(phy)){
stop("This is functionality has been temporarily removed.")
#samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / mean(f[as.numeric(names(freqs))+1]))
}else{
stop("The vector of sampling frequencies does not match the number of tips in the tree.")
}
}
if(sum(extirpation)==0){
init.pars <- starting.point.geosse(phy, eps=0, samp.freq.tree=samp.freq.tree)
}else{
init.pars <- starting.point.geosse(phy, eps=mag.san.start, samp.freq.tree=samp.freq.tree)
}
names(init.pars) <- NULL
if(is.null(starting.vals)){
def.set.pars <- rep(c(log(init.pars[1:3]), log(init.pars[4:5]), rep(log(init.pars[6:7]),3), rep(log(0.01), 12)), rate.cats)
}else{
## Check if 'starting.vals' has the correct format.
if( !length(starting.vals) %in% c(3,7) ){
stop("Wrong length of starting.vals")
}
if( length(starting.vals) == 7 ){
cat("Using developer mode for starting.vals.", "\n")
def.set.pars <- rep(c(log(starting.vals[1:3]), log(starting.vals[4:5]), rep(log(starting.vals[6:7]),3), rep(log(0.01), 12)), rate.cats)
} else{
def.set.pars <- rep(c(log( rep(starting.vals[1],3) ), log( rep(starting.vals[2],2) ), rep(log(starting.vals[3]),6), rep(log(0.01), 12)), rate.cats)
}
}
if(bounded.search == TRUE){
upper.full <- rep(c(rep(log(speciation.upper),3), rep(log(extirpation.upper),2), rep(log(trans.upper),2*3), rep(log(10), 12)), rate.cats)
}else{
upper.full <- rep(21,length(def.set.pars))
}
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]]
}
lower <- rep(-20, length(ip))
}else{
pars <- numeric(120)
if(dim(trans.rate)[2]==3){
rate.cats <- 1
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
trans.tmp <- c(trans.rate["(0)", "(1)"], trans.rate["(0)", "(01)"], trans.rate["(1)", "(0)"], trans.rate["(1)", "(01)"], trans.rate["(01)", "(0)"], trans.rate["(01)", "(1)"])
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
category.rates.unique <- 0
pars.tmp <- c(pars.tmp, trans.tmp)
pars[1:12] <- pars.tmp
}
if(dim(trans.rate)[2]==6){
rate.cats <- 2
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- c(category.rate.shift[1], rep(0,3), category.rate.shift[2], rep(0,3), category.rate.shift[3], rep(0,3))
category.rate.shiftB <- c(category.rate.shift[4], rep(0,3), category.rate.shift[5], rep(0,3), category.rate.shift[6], rep(0,3))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==9){
rate.cats <- 3
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- c(category.rate.shift[1:2], rep(0,2), category.rate.shift[3:4], rep(0,2), category.rate.shift[5:6], rep(0,2))
category.rate.shiftB <- c(category.rate.shift[7:8], rep(0,2), category.rate.shift[9:10], rep(0,2), category.rate.shift[11:12], rep(0,2))
category.rate.shiftC <- c(category.rate.shift[13:14], rep(0,2), category.rate.shift[15:16], rep(0,2), category.rate.shift[17:18], rep(0,2))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[7:9], trans.tmp[13:18], category.rate.shiftC)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==12){
rate.cats <- 4
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- c(category.rate.shift[1:3], rep(0,1), category.rate.shift[4:6], rep(0,1), category.rate.shift[7:9], rep(0,1))
category.rate.shiftB <- c(category.rate.shift[10:12], rep(0,1), category.rate.shift[13:15], rep(0,1), category.rate.shift[16:18], rep(0,1))
category.rate.shiftC <- c(category.rate.shift[19:21], rep(0,1), category.rate.shift[22:24], rep(0,1), category.rate.shift[25:27], rep(0,1))
category.rate.shiftD <- c(category.rate.shift[28:30], rep(0,1), category.rate.shift[31:33], rep(0,1), category.rate.shift[34:36], rep(0,1))
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[7:9], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[10:12], trans.tmp[19:24], category.rate.shiftD)
pars[1:length(pars.tmp)] <- pars.tmp
}
if(dim(trans.rate)[2]==15){
rate.cats <- 5
pars.tmp <- speciation
extirpation.tmp <- extirpation
extirpation.tmp[which(extirpation.tmp > 0)] = (extirpation.tmp[which( extirpation.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, extirpation.tmp)
rows <- c("(0A)", "(0A)", "(1A)", "(1A)", "(01A)", "(01A)", "(0B)", "(0B)", "(1B)", "(1B)", "(01B)", "(01B)", "(0C)", "(0C)", "(1C)", "(1C)", "(01C)", "(01C)", "(0D)", "(0D)", "(1D)", "(1D)", "(01D)", "(01D)", "(0E)", "(0E)", "(1E)", "(1E)", "(01E)", "(01E)")
cols <- c("(1A)", "(01A)", "(0A)", "(01A)", "(0A)", "(1A)", "(1B)", "(01B)", "(0B)", "(01B)", "(0B)", "(1B)", "(1C)", "(01C)", "(0C)", "(01C)", "(0C)", "(1C)", "(1D)", "(01D)", "(0D)", "(01D)", "(0D)", "(1D)", "(1E)", "(01E)", "(0E)", "(01E)", "(0E)", "(1E)")
trans.tmp <- trans.rate[cbind(rows,cols)]
trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
pars.tmp <- c(pars.tmp, trans.tmp)
category.tmp <- trans.rate[which(trans.rate==max(trans.rate, na.rm=TRUE))]
category.rate.shift <- rep(max(pars.tmp)+1, length(category.tmp))
category.rate.shiftA <- category.rate.shift[1:12]
category.rate.shiftB <- category.rate.shift[13:24]
category.rate.shiftC <- category.rate.shift[25:36]
category.rate.shiftD <- category.rate.shift[37:48]
category.rate.shiftE <- category.rate.shift[49:60]
category.rates.all <- c(category.rate.shiftA, category.rate.shiftB, category.rate.shiftC, category.rate.shiftD, category.rate.shiftE)
category.rates.unique <- length(unique(category.rates.all[category.rates.all>0]))
pars.tmp <- c(speciation[1:3], extirpation.tmp[1:3], trans.tmp[1:6], category.rate.shiftA, speciation[4:6], extirpation.tmp[4:6], trans.tmp[7:12], category.rate.shiftB, speciation[7:9], extirpation.tmp[7:9], trans.tmp[13:18], category.rate.shiftC, speciation[10:12], extirpation.tmp[10:12], trans.tmp[19:24], category.rate.shiftD, speciation[13:15], extirpation.tmp[13:15], trans.tmp[25:30], category.rate.shiftE)
pars[1:length(pars.tmp)] <- pars.tmp
}
np <- max(pars)
pars[pars==0] <- np+1
cat("Initializing...", "\n")
data.new <- data.frame(data[,2], data[,2], row.names=data[,1])
data.new <- data.new[phy$tip.label,]
#This is used to scale starting values to account for sampling:
if(length(f) == 3){
freqs <- table(data.new[,1])
if(length(freqs == 2)){
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}else{
samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / f[as.numeric(names(freqs))+1])
}
}else{
if(length(f) == Ntip(phy)){
stop("This functionality has been temporarily removed.")
#samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / mean(f[as.numeric(names(freqs))+1]))
}else{
stop("The vector of sampling frequencies does not match the number of tips in the tree.")
}
}
if(sum(extirpation)==0){
init.pars <- starting.point.generator(phy, 3, samp.freq.tree, yule=TRUE)
}else{
init.pars <- starting.point.generator(phy, 3, samp.freq.tree, yule=FALSE)
if(init.pars[4] == 0){
init.pars[4:6] = 1e-6
}
}
names(init.pars) <- NULL
if(is.null(starting.vals)){
def.set.pars <- rep(c(log(init.pars[1:3]), log(init.pars[4:6]), log(init.pars[7:12]), rep(log(.01), 12)), rate.cats)
}else{
## Check if 'starting.vals' has the correct format.
if( !length(starting.vals) %in% c(3,12) ){
stop("Wrong length of starting.vals")
}
if( length(starting.vals) == 12 ){
cat("Using developer mode for starting.vals.", "\n")
def.set.pars <- rep(c(log(starting.vals[1:3]), log(starting.vals[4:6]), log(starting.vals[7:12]), rep(log(0.01), 12)), rate.cats)
} else{
def.set.pars <- rep(c(log( rep(starting.vals[1],3) ), log( rep(starting.vals[2],2) ), rep(log(starting.vals[3]),6), rep(log(0.01), 12)), rate.cats)
}
}
if(bounded.search == TRUE){
upper.full <- rep(c(rep(log(speciation.upper),3), rep(log(extirpation.upper),3), rep(log(trans.upper),6), rep(log(10), 12)), rate.cats)
}else{
upper.full <- rep(21,length(def.set.pars))
}
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]]
}
lower <- rep(-20, length(ip))
}
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=DevOptimizeGeoHiSSE, ub=upper, lb=lower, opts=opts, pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
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=DevOptimizeGeoHiSSE, control=list(reltol=max.tol, parscale=rep(0.1, length(ip))), pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
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=DevOptimizeGeoHiSSE, lower=lower, upper=upper, control=list(max.call=sann.its), pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
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=DevOptimizeGeoHiSSE, ub=upper, lb=lower, opts=opts, pars=pars, phy=phy, data=data.new[,1], f=f, hidden.states=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, ode.eps=ode.eps)
solution <- numeric(length(pars))
solution[] <- c(exp(out$solution), 0)[pars]
loglik = -out$objective
}
if(assume.cladogenetic == TRUE){
names(solution) <- c("s0A", "s1A", "s01A", "x0A", "x1A", "d0A_1A ", "d0A_01A", "d1A_0A", "d1A_01A", "d01A_0A", "d01A_1A", "d0A_0B", "d0A_0C", "d0A_0D", "d0A_0E", "d1A_1B", "d1A_1C", "d1A_1D", "d1A_1E", "d01A_01B", "d01A_01C", "d01A_01D", "d01A_01E", "s0B", "s1B", "s01B", "x0B", "x1B", "d0B_1B ", "d0B_01B", "d1B_0B", "d1B_01B", "d01B_0B", "d01B_1B", "d0B_0A", "d0B_0C", "d0B_0D", "d0B_0E", "d1B_1A", "d1B_1C", "d1B_1D", "d1B_1E", "d01B_01A", "d01B_01C", "d01B_01D", "d01B_01E", "s0C", "s1C", "s01C", "x0C", "x1C", "d0C_1C ", "d0C_01C", "d1C_0C", "d1C_01C", "d01C_0C", "d01C_1C", "d0C_0A", "d0C_0B", "d0C_0D", "d0C_0E", "d1C_1A", "d1C_1B", "d1C_1D", "d1C_1E", "d01C_01A", "d01C_01B", "d01C_01D", "d01C_01E", "s0D", "s1D", "s01D", "x0D", "x1D", "d0D_1D ", "d0D_01D", "d1D_0D", "d1D_01D", "d01D_0D", "d01D_1D", "d0D_0A", "d0D_0B", "d0D_0C", "d0D_0E", "d1D_1A", "d1D_1B", "d1D_1C", "d1D_1E", "d01D_01A", "d01D_01B", "d01D_01C", "d01D_01E", "s0E", "s1E", "s01E", "x0E", "x1E", "d0E_1E ", "d0E_01E", "d1E_0E", "d1E_01E", "d01E_0E", "d01E_1E", "d0E_0A", "d0E_0B", "d0E_0C", "d0E_0D", "d1E_1A", "d1E_1B", "d1E_1C", "d1E_1D", "d01E_01A", "d01E_01B", "d01E_01C", "d01E_01D")
}else{
names(solution) <- c("s0A", "s1A", "s01A", "x0A", "x1A", "x01A", "d0A_1A ", "d0A_01A", "d1A_0A", "d1A_01A", "d01A_0A", "d01A_1A", "d0A_0B", "d0A_0C", "d0A_0D", "d0A_0E", "d1A_1B", "d1A_1C", "d1A_1D", "d1A_1E", "d01A_01B", "d01A_01C", "d01A_01D", "d01A_01E", "s0B", "s1B", "s01B", "x0B", "x1B", "x01B", "d0B_1B ", "d0B_01B", "d1B_0B", "d1B_01B", "d01B_0B", "d01B_1B", "d0B_0A", "d0B_0C", "d0B_0D", "d0B_0E", "d1B_1A", "d1B_1C", "d1B_1D", "d1B_1E", "d01B_01A", "d01B_01C", "d01B_01D", "d01B_01E", "s0C", "s1C", "s01C", "x0C", "x1C", "x01C", "d0C_1C ", "d0C_01C", "d1C_0C", "d1C_01C", "d01C_0C", "d01C_1C", "d0C_0A", "d0C_0B", "d0C_0D", "d0C_0E", "d1C_1A", "d1C_1B", "d1C_1D", "d1C_1E", "d01C_01A", "d01C_01B", "d01C_01D", "d01C_01E", "s0D", "s1D", "s01D", "x0D", "x1D", "x01D", "d0D_1D ", "d0D_01D", "d1D_0D", "d1D_01D", "d01D_0D", "d01D_1D", "d0D_0A", "d0D_0B", "d0D_0C", "d0D_0E", "d1D_1A", "d1D_1B", "d1D_1C", "d1D_1E", "d01D_01A", "d01D_01B", "d01D_01C", "d01D_01E", "s0E", "s1E", "s01E", "x0E", "x1E", "x01E", "d0E_1E ", "d0E_01E", "d1E_0E", "d1E_01E", "d01E_0E", "d01E_1E", "d0E_0A", "d0E_0B", "d0E_0C", "d0E_0D", "d1E_1A", "d1E_1B", "d1E_1C", "d1E_1D", "d01E_01A", "d01E_01B", "d01E_01C", "d01E_01D")
}
cat("Finished. Summarizing results...", "\n")
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.areas=hidden.areas, assume.cladogenetic=assume.cladogenetic, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, phy=phy, data=data, trans.matrix=trans.rate, max.tol=max.tol, starting.vals=ip, upper.bounds=upper, lower.bounds=lower, ode.eps=ode.eps)
class(obj) <- append(class(obj), "geohisse.old.fit")
return(obj)
}
######################################################################################################################################
######################################################################################################################################
### The function used to optimize parameters:
######################################################################################################################################
######################################################################################################################################
DevOptimizeGeoHiSSE <- function(p, pars, phy, data, f, hidden.states, assume.cladogenetic, condition.on.survival, root.type, root.p, np, ode.eps) {
#Generates the final vector with the appropriate parameter estimates in the right place:
p.new <- exp(p)
## print(p.new)
model.vec <- numeric(length(pars))
model.vec[] <- c(p.new, 0)[pars]
if(assume.cladogenetic == TRUE){
cache = ParametersToPassGeoHiSSE(phy=phy, data=data, f=f, model.vec=model.vec, hidden.states=hidden.states)
logl <- DownPassGeoHisse(phy, cache, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, ode.eps=ode.eps)
}else{
cache = ParametersToPassMuSSE(phy=phy, data=data, f=f, model.vec=model.vec, hidden.states=hidden.states)
logl <- DownPassMusse(phy, cache, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, ode.eps=ode.eps)
}
return(-logl)
}
#Taken from the GeoSSE code modified to account for sampling -- credit goes to Emma Goldberg:
starting.point.geosse <- function(tree, eps=0.5, samp.freq.tree) {
if (eps == 0) {
s <- (log(Ntip(tree)) - log(2)) / max(branching.times(tree))
s <- s/samp.freq.tree
x <- 0
d <- s/10
} else {
n <- Ntip(tree)
r <- ( log( (n/2) * (1 - eps*eps) + 2*eps + (1 - eps)/2 *
sqrt( n * (n*eps*eps - 8*eps + 2*n*eps + n))) - log(2)
) / max(branching.times(tree))
s <- r / (1 - eps)
x <- s * eps
s <- s/samp.freq.tree
x <- x - ((s * samp.freq.tree) * (1 - 1/samp.freq.tree))
d <- x
}
p <- c(s, s, s, x, x, d, d)
names(p) <- c("sA", "sB", "sAB", "xA" , "xB" ,"dA" ,"dB")
p
}
######################################################################################################################################
######################################################################################################################################
### The GeoHiSSE down pass that carries out the integration and returns the likelihood:
######################################################################################################################################
######################################################################################################################################
DownPassGeoHisse <- function(phy, cache, hidden.states, bad.likelihood=-10000000, condition.on.survival, root.type, root.p, get.phi=FALSE, node=NULL, state=NULL, ode.eps=0) {
#Some preliminaries:
nb.tip <- length(phy$tip.label)
nb.node <- phy$Nnode
phy <- reorder(phy, "pruningwise")
anc <- unique(phy$edge[,1])
TIPS <- 1:nb.tip
if(hidden.states == FALSE){
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
}else{
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
if(!is.null(root.p)){
root.p.new <- numeric(15)
root.p.new[1:length(root.p)] <- root.p
root.p <- root.p.new
}
}
#Initializes the tip sampling and sets internal nodes to be zero:
ncols = dim(compD)[2]
if(length(cache$f) == 3){
for(i in 1:(nb.tip)){
compD[i,] <- cache$f * cache$states[i,]
compE[i,] <- rep((1-cache$f), ncols/3)
}
}else{
for(i in 1:(nb.tip)){
compD[i,] <- cache$f[i] * cache$states[i,]
compE[i,] <- rep((1-cache$f[i]), ncols/3)
}
}
logcomp <- c()
#Start the postorder traversal indexing lists by node number:
for (i in seq(from = 1, length.out = nb.node)) {
#A vector of all the internal nodes:
focal <- anc[i]
desRows <- which(phy$edge[,1]==focal)
desNodes <- phy$edge[desRows,2]
#Note: when the tree has been reordered branching.times are no longer valid. Fortunately, we extract this information in the initial cache setup. Also focal is the rootward node, whereas desNodes represent a vector of all descendant nodes:
cache$rootward.age <- cache$split.times[which(names(cache$split.times)==focal)]
v <- c()
phi <- c()
for (desIndex in sequence(length(desRows))){
cache$focal.edge.length <- phy$edge.length[desRows[desIndex]]
cache$tipward.age <- cache$rootward.age - cache$focal.edge.length
#Strange rounding errors. A tip age should be zero. This ensures that:
if(cache$tipward.age < .Machine$double.eps^0.5){
cache$tipward.age = 0
}
cache$node.D <- compD[desNodes[desIndex],]
cache$node.E <- compE[desNodes[desIndex],]
##Call to lsoda that utilizes C code. Requires a lot of inputs. Note that for now we hardcode the NUMELEMENTS arguments. The reason for this is because with lsoda we can only pass a vector of parameters.
if(hidden.states == FALSE){
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A)
NUMELEMENTS <- 11 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <-c(E_0=cache$node.E[1], E_1=cache$node.E[2], E_01=cache$node.E[3], D_N0=cache$node.D[1], D_N1=cache$node.D[2], D_N2=cache$node.D[3])
times=c(cache$tipward.age, cache$rootward.age)
runSilent <- function() {
options(warn = -1)
on.exit(options(warn = 0))
capture.output(res <- lsoda(yini, times, func = "classe_geosse_equivalent_derivs", padded.pars, initfunc="initmod_geosse", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "classe_geosse_equivalent_derivs", padded.pars, initfunc="initmod_geosse", dll = "canonical_geosse-ext-derivs", rtol=1e-8, atol=1e-8))
res
}
prob.subtree.cal.full <- runSilent()
}else{
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A, cache$d0A_0B, cache$d0A_0C, cache$d0A_0D, cache$d0A_0E, cache$d1A_1B, cache$d1A_1C, cache$d1A_1D, cache$d1A_1E, cache$d01A_01B, cache$d01A_01C, cache$d01A_01D, cache$d01A_01E, cache$s0B, cache$s1B, cache$s01B, cache$x0B, cache$x1B, cache$d0B_1B , cache$d0B_01B, cache$d1B_0B, cache$d1B_01B, cache$d01B_0B, cache$d01B_1B, cache$d0B_0A, cache$d0B_0C, cache$d0B_0D, cache$d0B_0E, cache$d1B_1A, cache$d1B_1C, cache$d1B_1D, cache$d1B_1E, cache$d01B_01A, cache$d01B_01C, cache$d01B_01D, cache$d01B_01E, cache$s0C, cache$s1C, cache$s01C, cache$x0C, cache$x1C, cache$d0C_1C , cache$d0C_01C, cache$d1C_0C, cache$d1C_01C, cache$d01C_0C, cache$d01C_1C, cache$d0C_0A, cache$d0C_0B, cache$d0C_0D, cache$d0C_0E, cache$d1C_1A, cache$d1C_1B, cache$d1C_1D, cache$d1C_1E, cache$d01C_01A, cache$d01C_01B, cache$d01C_01D, cache$d01C_01E, cache$s0D, cache$s1D, cache$s01D, cache$x0D, cache$x1D, cache$d0D_1D , cache$d0D_01D, cache$d1D_0D, cache$d1D_01D, cache$d01D_0D, cache$d01D_1D, cache$d0D_0A, cache$d0D_0B, cache$d0D_0C, cache$d0D_0E, cache$d1D_1A, cache$d1D_1B, cache$d1D_1C, cache$d1D_1E, cache$d01D_01A, cache$d01D_01B, cache$d01D_01C, cache$d01D_01E, cache$s0E, cache$s1E, cache$s01E, cache$x0E, cache$x1E, cache$d0E_1E , cache$d0E_01E, cache$d1E_0E, cache$d1E_01E, cache$d01E_0E, cache$d01E_1E, cache$d0E_0A, cache$d0E_0B, cache$d0E_0C, cache$d0E_0D, cache$d1E_1A, cache$d1E_1B, cache$d1E_1C, cache$d1E_1D, cache$d01E_01A, cache$d01E_01B, cache$d01E_01C, cache$d01E_01D)
NUMELEMENTS <- 115 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <- c(E_0A=cache$node.E[1], E_1A=cache$node.E[2], E_01A=cache$node.E[3], E_0B=cache$node.E[4], E_1B=cache$node.E[5], E_01B=cache$node.E[6], E_0C=cache$node.E[7], E_1C=cache$node.E[8], E_01C=cache$node.E[9], E_0D=cache$node.E[10], E_1D=cache$node.E[11], E_01D=cache$node.E[12], E_0E=cache$node.E[13], E_1E=cache$node.E[14], E_01E=cache$node.E[15], D_N0A=cache$node.D[1], D_N1A=cache$node.D[2], D_N01A=cache$node.D[3], D_N0B=cache$node.D[4], D_N1B=cache$node.D[5], D_N01B=cache$node.D[6], D_N0C=cache$node.D[7], D_N1C=cache$node.D[8], D_N01C=cache$node.D[9], D_N0D=cache$node.D[10], D_N1D=cache$node.D[11], D_N01D=cache$node.D[12], D_N0E=cache$node.D[13], D_N1E=cache$node.D[14], D_N01E=cache$node.D[15])
times=c(cache$tipward.age, cache$rootward.age)
runSilent <- function() {
options(warn = -1)
on.exit(options(warn = 0))
capture.output(res <- lsoda(yini, times, func = "geohisse_derivs", padded.pars, initfunc="initmod_geohisse", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "geohisse_derivs", padded.pars, initfunc="initmod_geohisse", dll = "geohisse-ext-derivs", rtol=1e-8, atol=1e-8))
res
}
prob.subtree.cal.full <- runSilent()
}
######## THIS CHECKS TO ENSURE THAT THE INTEGRATION WAS SUCCESSFUL ###########
if(attributes(prob.subtree.cal.full)$istate[1] < 0){
return(bad.likelihood)
}else{
prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
}
##############################################################################
if(hidden.states == FALSE){
if(is.nan(prob.subtree.cal[3]) | is.nan(prob.subtree.cal[4]) | is.nan(prob.subtree.cal[5])){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(prob.subtree.cal[4]<0 | prob.subtree.cal[5]<0 | prob.subtree.cal[6]<0){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[4:6]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:3])
v <- rbind(v, prob.subtree.cal[4:6])
}else{
if(any(is.nan(prob.subtree.cal[16:30]))){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(any(prob.subtree.cal[16:30] < 0)){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[16:30]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:15])
v <- rbind(v, prob.subtree.cal[16:30])
}
}
if(hidden.states == TRUE){
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0A + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1A + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01A
v <- v[,-c(1:3)]
compD[focal,4] <- v[1,1] * v[2,1] * cache$s0B
compD[focal,5] <- v[1,2] * v[2,2] * cache$s1B
compD[focal,6] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0B + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1B + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01B
v <- v[,-c(1:3)]
compD[focal,7] <- v[1,1] * v[2,1] * cache$s0C
compD[focal,8] <- v[1,2] * v[2,2] * cache$s1C
compD[focal,9] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0C + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1C + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01C
v <- v[,-c(1:3)]
compD[focal,10] <- v[1,1] * v[2,1] * cache$s0D
compD[focal,11] <- v[1,2] * v[2,2] * cache$s1D
compD[focal,12] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0D + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1D + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01D
v <- v[,-c(1:3)]
compD[focal,13] <- v[1,1] * v[2,1] * cache$s0E
compD[focal,14] <- v[1,2] * v[2,2] * cache$s1E
compD[focal,15] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0E + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1E + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01E
compE[focal,] <- phi[1:15]
if(!is.null(node)){
fixer = numeric(15)
fixer[state] = 1
if(node == focal){
compD[focal,] <- compD[focal,] * fixer
}
#compE[focal,] <- compE[focal,] * fixer
}
}else{
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- 0.5 * (v[1,3] * v[2,1] + v[1,1] * v[2,3]) * cache$s0A + 0.5 * (v[1,3] * v[2,2] + v[1,2] * v[2,3]) * cache$s1A + 0.5 * (v[1,1] * v[2,2] + v[1,2] * v[2,1]) * cache$s01A
compE[focal,] <- phi[1:3]
if(!is.null(node)){
if(node == focal){
fixer = c(0,0,0)
fixer[state] = 1
compD[focal,] <- compD[focal,] * fixer
}
}
}
###########################
#Logcompensation bit for dealing with underflow issues. Need to give a necessary shoutout to Rich FitzJohn -- we follow his diversitree approach. VERIFIED that it works properly:
tmp <- sum(compD[focal,])
compD[focal,] <- compD[focal,] / tmp
logcomp <- c(logcomp, log(tmp))
}
root.node <- nb.tip + 1L
if (is.na(sum(log(compD[root.node,]))) || is.na(log(sum(1-compE[root.node,])))){
return(bad.likelihood)
}else{
if(root.type == "madfitz" | root.type == "herr_als"){
if(is.null(root.p)){
root.p = compD[root.node,]/sum(compD[root.node,])
root.p[which(is.na(root.p))] = 0
}
}
if(condition.on.survival == TRUE){
if(hidden.states == FALSE){
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)))
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)))
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}else{
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)), cache$s0B, cache$s1B, sum(c(cache$s0B, cache$s1B, cache$s01B)), cache$s0C, cache$s1C, sum(c(cache$s0C, cache$s1C, cache$s01C)), cache$s0D, cache$s1D, sum(c(cache$s0D, cache$s1D, cache$s01D)), cache$s0E, cache$s1E, sum(c(cache$s0E, cache$s1E, cache$s01E)))
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, sum(c(cache$s0A, cache$s1A, cache$s01A)), cache$s0B, cache$s1B, sum(c(cache$s0B, cache$s1B, cache$s01B)), cache$s0C, cache$s1C, sum(c(cache$s0C, cache$s1C, cache$s01C)), cache$s0D, cache$s1D, sum(c(cache$s0D, cache$s1D, cache$s01D)), cache$s0E, cache$s1E, sum(c(cache$s0E, cache$s1E, cache$s01E)))
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}
}
}
if(get.phi==TRUE){
obj = NULL
obj$compD.root = compD[root.node,]/sum(compD[root.node,])
obj$compE = compE
obj$root.p = root.p
return(obj)
}else{
return(loglik)
}
}
######################################################################################################################################
######################################################################################################################################
### The MuSSE type down pass that carries out the integration and returns the likelihood:
######################################################################################################################################
######################################################################################################################################
DownPassMusse <- function(phy, cache, hidden.states, bad.likelihood=-10000000, condition.on.survival, root.type, root.p, get.phi=FALSE, node=NULL, state=NULL, ode.eps=0) {
#Some preliminaries:
nb.tip <- length(phy$tip.label)
nb.node <- phy$Nnode
phy <- reorder(phy, "pruningwise")
anc <- unique(phy$edge[,1])
TIPS <- 1:nb.tip
if(hidden.states == FALSE){
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=3)
}else{
compD <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
compE <- matrix(0, nrow=nb.tip + nb.node, ncol=15)
if(!is.null(root.p)){
root.p.new <- numeric(15)
root.p.new[1:length(root.p)] <- root.p
root.p <- root.p.new
}
}
#Initializes the tip sampling and sets internal nodes to be zero:
ncols = dim(compD)[2]
if(length(cache$f) == 3){
for(i in 1:(nb.tip)){
compD[i,] <- cache$f * cache$states[i,]
compE[i,] <- rep((1-cache$f), ncols/3)
}
}else{
for(i in 1:(nb.tip)){
compD[i,] <- cache$f[i] * cache$states[i,]
compE[i,] <- rep((1-cache$f[i]), ncols/3)
}
}
logcomp <- c()
#Start the postorder traversal indexing lists by node number:
for (i in seq(from = 1, length.out = nb.node)) {
#A vector of all the internal nodes:
focal <- anc[i]
desRows <- which(phy$edge[,1]==focal)
desNodes <- phy$edge[desRows,2]
#Note: when the tree has been reordered branching.times are no longer valid. Fortunately, we extract this information in the initial cache setup. Also focal is the rootward node, whereas desNodes represent a vector of all descendant nodes:
cache$rootward.age <- cache$split.times[which(names(cache$split.times)==focal)]
v <- c()
phi <- c()
for (desIndex in sequence(length(desRows))){
cache$focal.edge.length <- phy$edge.length[desRows[desIndex]]
cache$tipward.age <- cache$rootward.age - cache$focal.edge.length
#Strange rounding errors. A tip age should be zero. This ensures that:
if(cache$tipward.age < .Machine$double.eps^0.5){
cache$tipward.age = 0
}
cache$node.D <- compD[desNodes[desIndex],]
cache$node.E <- compE[desNodes[desIndex],]
##Call to lsoda that utilizes C code. Requires a lot of inputs. Note that for now we hardcode the NUMELEMENTS arguments. The reason for this is because with lsoda we can only pass a vector of parameters.
if(hidden.states == FALSE){
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$x01A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A)
NUMELEMENTS <- 12 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <-c(E_0=cache$node.E[1], E_1=cache$node.E[2], E_01=cache$node.E[3], D_N0=cache$node.D[1], D_N1=cache$node.D[2], D_N2=cache$node.D[3])
times=c(cache$tipward.age, cache$rootward.age)
#runSilent <- function() {
#options(warn = -1)
#on.exit(options(warn = 0))
#capture.output(res <- lsoda(yini, times, func = "notclasse_derivs", padded.pars, initfunc="initmod_musse", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "notclasse_derivs", padded.pars, initfunc="initmod_musse", dll = "notclasse-ext-derivs", rtol=1e-8, atol=1e-8))
#res
#}
#prob.subtree.cal.full <- runSilent()
prob.subtree.cal.full <- lsoda(yini, times, func = "notclasse_derivs", padded.pars, initfunc="initmod_noclass", dllname = "hisse", rtol=1e-8, atol=1e-8)
}else{
pars <- list(cache$s0A, cache$s1A, cache$s01A, cache$x0A, cache$x1A, cache$x01A, cache$d0A_1A, cache$d0A_01A, cache$d1A_0A, cache$d1A_01A, cache$d01A_0A, cache$d01A_1A, cache$d0A_0B, cache$d0A_0C, cache$d0A_0D, cache$d0A_0E, cache$d1A_1B, cache$d1A_1C, cache$d1A_1D, cache$d1A_1E, cache$d01A_01B, cache$d01A_01C, cache$d01A_01D, cache$d01A_01E, cache$s0B, cache$s1B, cache$s01B, cache$x0B, cache$x1B, cache$x01B, cache$d0B_1B , cache$d0B_01B, cache$d1B_0B, cache$d1B_01B, cache$d01B_0B, cache$d01B_1B, cache$d0B_0A, cache$d0B_0C, cache$d0B_0D, cache$d0B_0E, cache$d1B_1A, cache$d1B_1C, cache$d1B_1D, cache$d1B_1E, cache$d01B_01A, cache$d01B_01C, cache$d01B_01D, cache$d01B_01E, cache$s0C, cache$s1C, cache$s01C, cache$x0C, cache$x1C, cache$x01C, cache$d0C_1C , cache$d0C_01C, cache$d1C_0C, cache$d1C_01C, cache$d01C_0C, cache$d01C_1C, cache$d0C_0A, cache$d0C_0B, cache$d0C_0D, cache$d0C_0E, cache$d1C_1A, cache$d1C_1B, cache$d1C_1D, cache$d1C_1E, cache$d01C_01A, cache$d01C_01B, cache$d01C_01D, cache$d01C_01E, cache$s0D, cache$s1D, cache$s01D, cache$x0D, cache$x1D, cache$x01D, cache$d0D_1D , cache$d0D_01D, cache$d1D_0D, cache$d1D_01D, cache$d01D_0D, cache$d01D_1D, cache$d0D_0A, cache$d0D_0B, cache$d0D_0C, cache$d0D_0E, cache$d1D_1A, cache$d1D_1B, cache$d1D_1C, cache$d1D_1E, cache$d01D_01A, cache$d01D_01B, cache$d01D_01C, cache$d01D_01E, cache$s0E, cache$s1E, cache$s01E, cache$x0E, cache$x1E, cache$x01E, cache$d0E_1E, cache$d0E_01E, cache$d1E_0E, cache$d1E_01E, cache$d01E_0E, cache$d01E_1E, cache$d0E_0A, cache$d0E_0B, cache$d0E_0C, cache$d0E_0D, cache$d1E_1A, cache$d1E_1B, cache$d1E_1C, cache$d1E_1D, cache$d01E_01A, cache$d01E_01B, cache$d01E_01C, cache$d01E_01D)
NUMELEMENTS <- 120 #needed for passing in vector to C
padded.pars <- rep(0, NUMELEMENTS)
pars <- c(unlist(pars))
stopifnot(length(padded.pars)<=NUMELEMENTS)
padded.pars[sequence(length(pars))]<-pars
yini <- c(E_0A=cache$node.E[1], E_1A=cache$node.E[2], E_01A=cache$node.E[3], E_0B=cache$node.E[4], E_1B=cache$node.E[5], E_01B=cache$node.E[6], E_0C=cache$node.E[7], E_1C=cache$node.E[8], E_01C=cache$node.E[9], E_0D=cache$node.E[10], E_1D=cache$node.E[11], E_01D=cache$node.E[12], E_0E=cache$node.E[13], E_1E=cache$node.E[14], E_01E=cache$node.E[15], D_N0A=cache$node.D[1], D_N1A=cache$node.D[2], D_N01A=cache$node.D[3], D_N0B=cache$node.D[4], D_N1B=cache$node.D[5], D_N01B=cache$node.D[6], D_N0C=cache$node.D[7], D_N1C=cache$node.D[8], D_N01C=cache$node.D[9], D_N0D=cache$node.D[10], D_N1D=cache$node.D[11], D_N01D=cache$node.D[12], D_N0E=cache$node.D[13], D_N1E=cache$node.D[14], D_N01E=cache$node.D[15])
times=c(cache$tipward.age, cache$rootward.age)
#runSilent <- function() {
# options(warn = -1)
# on.exit(options(warn = 0))
# capture.output(res <- lsoda(yini, times, func = "notclasse_more_derivs", padded.pars, initfunc="initmod_mussem", dllname = "hisse", rtol=1e-8, atol=1e-8))
#capture.output(res <- lsoda(yini, times, func = "notclasse_more_derivs", padded.pars, initfunc="initmod_mussem", dll = "notclasse-more-ext-derivs", rtol=1e-8, atol=1e-8))
# res
#}
#prob.subtree.cal.full <- runSilent()
prob.subtree.cal.full <- lsoda(yini, times, func = "notclasse_more_derivs", padded.pars, initfunc="initmod_hinoclass", dllname = "hisse", rtol=1e-8, atol=1e-8)
}
######## THIS CHECKS TO ENSURE THAT THE INTEGRATION WAS SUCCESSFUL ###########
if(attributes(prob.subtree.cal.full)$istate[1] < 0){
return(bad.likelihood)
}else{
prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
}
##############################################################################
if(hidden.states == FALSE){
if(is.nan(prob.subtree.cal[3]) | is.nan(prob.subtree.cal[4]) | is.nan(prob.subtree.cal[5])){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(prob.subtree.cal[4]<0 | prob.subtree.cal[5]<0 | prob.subtree.cal[6]<0){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[4:6]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:3])
v <- rbind(v, prob.subtree.cal[4:6])
}else{
if(any(is.nan(prob.subtree.cal[16:30]))){
return(bad.likelihood)
}
#This is default and cannot change, but if we get a negative probability, discard the results:
if(any(prob.subtree.cal[16:30] < 0)){
return(bad.likelihood)
}
#This can be modified at the input, but if the sum of the D's at the end of a branch are less than some value, then discard the results. A little more stringent than diversitree, but with difficult problems, this stabilizes things immensely.
if(sum(prob.subtree.cal[10:18]) < ode.eps){
return(bad.likelihood)
}
#Designating phi here because of its relation to Morlon et al (2011) and using "e" would be confusing:
phi <- c(phi, prob.subtree.cal[1:15])
v <- rbind(v, prob.subtree.cal[16:30])
}
}
##Allows for a 3-state MuSSE type model. Important for cases where you have a random tree and a random trait.
if(hidden.states == TRUE){
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- v[1,3] * v[2,3] * cache$s01A
v <- v[,-c(1:3)]
compD[focal,4] <- v[1,1] * v[2,1] * cache$s0B
compD[focal,5] <- v[1,2] * v[2,2] * cache$s1B
compD[focal,6] <- v[1,3] * v[2,3] * cache$s01B
v <- v[,-c(1:3)]
compD[focal,7] <- v[1,1] * v[2,1] * cache$s0C
compD[focal,8] <- v[1,2] * v[2,2] * cache$s1C
compD[focal,9] <- v[1,3] * v[2,3] * cache$s01C
v <- v[,-c(1:3)]
compD[focal,10] <- v[1,1] * v[2,1] * cache$s0D
compD[focal,11] <- v[1,2] * v[2,2] * cache$s1D
compD[focal,12] <- v[1,3] * v[2,3] * cache$s01D
v <- v[,-c(1:3)]
compD[focal,13] <- v[1,1] * v[2,1] * cache$s0E
compD[focal,14] <- v[1,2] * v[2,2] * cache$s1E
compD[focal,15] <- v[1,3] * v[2,3] * cache$s01E
compE[focal,] <- phi[1:15]
if(!is.null(node)){
fixer = numeric(15)
fixer[state] = 1
if(node == focal){
compD[focal,] <- compD[focal,] * fixer
}
#compE[focal,] <- compE[focal,] * fixer
}
}else{
compD[focal,1] <- v[1,1] * v[2,1] * cache$s0A
compD[focal,2] <- v[1,2] * v[2,2] * cache$s1A
compD[focal,3] <- v[1,3] * v[2,3] * cache$s01A
compE[focal,] <- phi[1:3]
if(!is.null(node)){
if(node == focal){
fixer = c(0,0,0)
fixer[state] = 1
compD[focal,] <- compD[focal,] * fixer
}
}
}
###########################
#Logcompensation bit for dealing with underflow issues. Need to give a necessary shoutout to Rich FitzJohn -- we follow his diversitree approach. VERIFIED that it works properly:
tmp <- sum(compD[focal,])
compD[focal,] <- compD[focal,] / tmp
logcomp <- c(logcomp, log(tmp))
}
root.node <- nb.tip + 1L
if (is.na(sum(log(compD[root.node,]))) || is.na(log(sum(1-compE[root.node,])))){
return(bad.likelihood)
}else{
if(root.type == "madfitz" | root.type == "herr_als"){
if(is.null(root.p)){
root.p = compD[root.node,]/sum(compD[root.node,])
root.p[which(is.na(root.p))] = 0
}
}
if(condition.on.survival == TRUE){
if(hidden.states == FALSE){
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, cache$s01A)
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, cache$s01A)
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}else{
if(root.type == "madfitz"){
lambda <- c(cache$s0A, cache$s1A, cache$s01A, cache$s0B, cache$s1B, cache$s01B, cache$s0C, cache$s1C, cache$s01C, cache$s0D, cache$s1D, cache$s01D, cache$s0E, cache$s1E, cache$s01E)
compD[root.node,] <- compD[root.node,] / sum(root.p * lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,] * root.p)) + sum(logcomp)
}else{
lambda <- c(cache$s0A, cache$s1A, cache$s01A, cache$s0B, cache$s1B, cache$s01B, cache$s0C, cache$s1C, cache$s01C, cache$s0D, cache$s1D, cache$s01D, cache$s0E, cache$s1E, cache$s01E)
compD[root.node,] <- (compD[root.node,] * root.p) / (lambda * (1 - compE[root.node,])^2)
#Corrects for possibility that you have 0/0:
compD[root.node,which(is.na(compD[root.node,]))] = 0
loglik <- log(sum(compD[root.node,])) + sum(logcomp)
}
}
}
}
if(get.phi==TRUE){
obj = NULL
obj$compD.root = compD[root.node,]/sum(compD[root.node,])
obj$compE = compE
obj$root.p = root.p
return(obj)
}else{
return(loglik)
}
}
######################################################################################################################################
######################################################################################################################################
### Cache object for storing parameters that are used throughout GeoHiSSE:
######################################################################################################################################
######################################################################################################################################
ParametersToPassGeoHiSSE <- function(phy, data, f, model.vec, hidden.states){
#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$phy = phy
if(hidden.states == FALSE){
states = matrix(0,Ntip(phy),3)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==0){states[i,3]=1}
}
}
if(hidden.states == TRUE){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,c(1,4,7,10,13)]=1}
if(data[i]==2){states[i,c(2,5,8,11,14)]=1}
if(data[i]==0){states[i,c(3,6,9,12,15)]=1}
}
}
if(hidden.states == "TEST"){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==0){states[i,3]=1}
if(data[i]==4){states[i,4]=1}
if(data[i]==5){states[i,5]=1}
if(data[i]==3){states[i,6]=1}
}
}
obj$states = states
obj$tot_time = max(branching.times(phy))
obj$f = f
obj$s0A = model.vec[1]
obj$s1A = model.vec[2]
obj$s01A = model.vec[3]
obj$x0A = model.vec[4]
obj$x1A = model.vec[5]
obj$d0A_1A = model.vec[6]
obj$d0A_01A = model.vec[7]
obj$d1A_0A = model.vec[8]
obj$d1A_01A = model.vec[9]
if(model.vec[10]==0){
obj$d01A_0A = model.vec[5]
}else{
obj$d01A_0A = model.vec[10]
}
if(model.vec[11]==0){
obj$d01A_1A = model.vec[4]
}else{
obj$d01A_1A = model.vec[11]
}
obj$d0A_0B = model.vec[12]
obj$d0A_0C = model.vec[13]
obj$d0A_0D = model.vec[14]
obj$d0A_0E = model.vec[15]
obj$d1A_1B = model.vec[16]
obj$d1A_1C = model.vec[17]
obj$d1A_1D = model.vec[18]
obj$d1A_1E = model.vec[19]
obj$d01A_01B = model.vec[20]
obj$d01A_01C = model.vec[21]
obj$d01A_01D = model.vec[22]
obj$d01A_01E = model.vec[23]
obj$s0B = model.vec[24]
obj$s1B = model.vec[25]
obj$s01B = model.vec[26]
obj$x0B = model.vec[27]
obj$x1B = model.vec[28]
obj$d0B_1B = model.vec[29]
obj$d0B_01B = model.vec[30]
obj$d1B_0B = model.vec[31]
obj$d1B_01B = model.vec[32]
if(model.vec[33] == 0){
obj$d01B_0B = model.vec[28]
}else{
obj$d01B_0B = model.vec[33]
}
if(model.vec[34] == 0){
obj$d01B_1B = model.vec[27]
}else{
obj$d01B_1B = model.vec[34]
}
obj$d0B_0A = model.vec[35]
obj$d0B_0C = model.vec[36]
obj$d0B_0D = model.vec[37]
obj$d0B_0E = model.vec[38]
obj$d1B_1A = model.vec[39]
obj$d1B_1C = model.vec[40]
obj$d1B_1D = model.vec[41]
obj$d1B_1E = model.vec[42]
obj$d01B_01A = model.vec[43]
obj$d01B_01C = model.vec[44]
obj$d01B_01D = model.vec[45]
obj$d01B_01E = model.vec[46]
obj$s0C = model.vec[47]
obj$s1C = model.vec[48]
obj$s01C = model.vec[49]
obj$x0C = model.vec[50]
obj$x1C = model.vec[51]
obj$d0C_1C = model.vec[52]
obj$d0C_01C = model.vec[53]
obj$d1C_0C = model.vec[54]
obj$d1C_01C = model.vec[55]
if(model.vec[56] == 0){
obj$d01C_0C = model.vec[51]
}else{
obj$d01C_0C = model.vec[56]
}
if(model.vec[57] == 0){
obj$d01C_1C = model.vec[50]
}else{
obj$d01C_1C = model.vec[57]
}
obj$d0C_0A = model.vec[58]
obj$d0C_0B = model.vec[59]
obj$d0C_0D = model.vec[60]
obj$d0C_0E = model.vec[61]
obj$d1C_1A = model.vec[62]
obj$d1C_1B = model.vec[63]
obj$d1C_1D = model.vec[64]
obj$d1C_1E = model.vec[65]
obj$d01C_01A = model.vec[66]
obj$d01C_01B = model.vec[67]
obj$d01C_01D = model.vec[68]
obj$d01C_01E = model.vec[69]
obj$s0D = model.vec[70]
obj$s1D = model.vec[71]
obj$s01D = model.vec[72]
obj$x0D = model.vec[73]
obj$x1D = model.vec[74]
obj$d0D_1D = model.vec[75]
obj$d0D_01D = model.vec[76]
obj$d1D_0D = model.vec[77]
obj$d1D_01D = model.vec[78]
if(model.vec[79] == 0){
obj$d01D_0D = model.vec[74]
}else{
obj$d01D_0D = model.vec[79]
}
if(model.vec[80] == 0){
obj$d01D_1D = model.vec[73]
}else{
obj$d01D_1D = model.vec[80]
}
obj$d0D_0A = model.vec[81]
obj$d0D_0B = model.vec[82]
obj$d0D_0C = model.vec[83]
obj$d0D_0E = model.vec[84]
obj$d1D_1A = model.vec[85]
obj$d1D_1B = model.vec[86]
obj$d1D_1C = model.vec[87]
obj$d1D_1E = model.vec[88]
obj$d01D_01A = model.vec[89]
obj$d01D_01B = model.vec[90]
obj$d01D_01C = model.vec[91]
obj$d01D_01E = model.vec[92]
obj$s0E = model.vec[93]
obj$s1E = model.vec[94]
obj$s01E = model.vec[95]
obj$x0E = model.vec[96]
obj$x1E = model.vec[97]
obj$d0E_1E = model.vec[98]
obj$d0E_01E = model.vec[99]
obj$d1E_0E = model.vec[100]
obj$d1E_01E = model.vec[101]
if(model.vec[102] == 0){
obj$d01E_0E = model.vec[97]
}else{
obj$d01E_0E = model.vec[102]
}
if(model.vec[103] == 0){
obj$d01E_1E = model.vec[96]
}else{
obj$d01E_1E = model.vec[103]
}
obj$d0E_0A = model.vec[104]
obj$d0E_0B = model.vec[105]
obj$d0E_0C = model.vec[106]
obj$d0E_0D = model.vec[107]
obj$d1E_1A = model.vec[108]
obj$d1E_1B = model.vec[109]
obj$d1E_1C = model.vec[110]
obj$d1E_1D = model.vec[111]
obj$d01E_01A = model.vec[112]
obj$d01E_01B = model.vec[113]
obj$d01E_01C = model.vec[114]
obj$d01E_01D = model.vec[115]
obj$split.times = sort(branching.times(phy), decreasing=TRUE)
return(obj)
}
######################################################################################################################################
######################################################################################################################################
### Cache object for storing parameters that are used throughout MuSSE:
######################################################################################################################################
######################################################################################################################################
ParametersToPassMuSSE <- function(phy, data, f, model.vec, hidden.states){
#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$phy = phy
if(hidden.states == FALSE){
states = matrix(0,Ntip(phy),3)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==0){states[i,3]=1}
}
}
if(hidden.states == TRUE){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,c(1,4,7,10,13)]=1}
if(data[i]==2){states[i,c(2,5,8,11,14)]=1}
if(data[i]==0){states[i,c(3,6,9,12,15)]=1}
}
}
if(hidden.states == "TEST1"){
states = matrix(0,Ntip(phy),3)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,1]=1}
if(data[i]==2){states[i,2]=1}
if(data[i]==4){states[i,3]=1}
}
}
if(hidden.states == "TEST2"){
states = matrix(0,Ntip(phy),15)
for(i in 1:Ntip(phy)){
if(data[i]==1){states[i,4]=1}
if(data[i]==2){states[i,5]=1}
if(data[i]==4){states[i,6]=1}
}
}
obj$states = states
obj$tot_time = max(branching.times(phy))
obj$f = f
obj$s0A = model.vec[1]
obj$s1A = model.vec[2]
obj$s01A = model.vec[3]
obj$x0A = model.vec[4]
obj$x1A = model.vec[5]
obj$x01A = model.vec[6]
obj$d0A_1A = model.vec[7]
obj$d0A_01A = model.vec[8]
obj$d1A_0A = model.vec[9]
obj$d1A_01A = model.vec[10]
if(model.vec[11] == 0){
obj$d01A_0A = model.vec[5]
}else{
obj$d01A_0A = model.vec[11]
}
if(model.vec[12] == 0){
obj$d01A_1A = model.vec[4]
}else{
obj$d01A_1A = model.vec[12]
}
obj$d0A_0B = model.vec[13]
obj$d0A_0C = model.vec[14]
obj$d0A_0D = model.vec[15]
obj$d0A_0E = model.vec[16]
obj$d1A_1B = model.vec[17]
obj$d1A_1C = model.vec[18]
obj$d1A_1D = model.vec[19]
obj$d1A_1E = model.vec[20]
obj$d01A_01B = model.vec[21]
obj$d01A_01C = model.vec[22]
obj$d01A_01D = model.vec[23]
obj$d01A_01E = model.vec[24]
obj$s0B = model.vec[25]
obj$s1B = model.vec[26]
obj$s01B = model.vec[27]
obj$x0B = model.vec[28]
obj$x1B = model.vec[29]
obj$x01B = model.vec[30]
obj$d0B_1B = model.vec[31]
obj$d0B_01B = model.vec[32]
obj$d1B_0B = model.vec[33]
obj$d1B_01B = model.vec[34]
if(model.vec[35] == 0){
obj$d01B_0B = model.vec[29]
}else{
obj$d01B_0B = model.vec[35]
}
if(model.vec[36] == 0){
obj$d01B_1B = model.vec[28]
}else{
obj$d01B_1B = model.vec[36]
}
obj$d0B_0A = model.vec[37]
obj$d0B_0C = model.vec[38]
obj$d0B_0D = model.vec[39]
obj$d0B_0E = model.vec[40]
obj$d1B_1A = model.vec[41]
obj$d1B_1C = model.vec[42]
obj$d1B_1D = model.vec[43]
obj$d1B_1E = model.vec[44]
obj$d01B_01A = model.vec[45]
obj$d01B_01C = model.vec[46]
obj$d01B_01D = model.vec[47]
obj$d01B_01E = model.vec[48]
obj$s0C = model.vec[49]
obj$s1C = model.vec[50]
obj$s01C = model.vec[51]
obj$x0C = model.vec[52]
obj$x1C = model.vec[53]
obj$x01C = model.vec[54]
obj$d0C_1C = model.vec[55]
obj$d0C_01C = model.vec[56]
obj$d1C_0C = model.vec[57]
obj$d1C_01C = model.vec[58]
if(model.vec[59] == 0){
obj$d01C_0C = model.vec[53]
}else{
obj$d01C_0C = model.vec[59]
}
if(model.vec[60] == 0){
obj$d01C_1C = model.vec[52]
}else{
obj$d01C_1C = model.vec[60]
}
obj$d0C_0A = model.vec[61]
obj$d0C_0B = model.vec[62]
obj$d0C_0D = model.vec[63]
obj$d0C_0E = model.vec[64]
obj$d1C_1A = model.vec[65]
obj$d1C_1B = model.vec[66]
obj$d1C_1D = model.vec[67]
obj$d1C_1E = model.vec[68]
obj$d01C_01A = model.vec[69]
obj$d01C_01B = model.vec[70]
obj$d01C_01D = model.vec[71]
obj$d01C_01E = model.vec[72]
obj$s0D = model.vec[73]
obj$s1D = model.vec[74]
obj$s01D = model.vec[75]
obj$x0D = model.vec[76]
obj$x1D = model.vec[77]
obj$x01D = model.vec[78]
obj$d0D_1D = model.vec[79]
obj$d0D_01D = model.vec[80]
obj$d1D_0D = model.vec[81]
obj$d1D_01D = model.vec[82]
if(model.vec[83] == 0){
obj$d01D_0D = model.vec[77]
}else{
obj$d01D_0D = model.vec[83]
}
if(model.vec[84] == 0){
obj$d01D_1D = model.vec[76]
}else{
obj$d01D_1D = model.vec[84]
}
obj$d0D_0A = model.vec[85]
obj$d0D_0B = model.vec[86]
obj$d0D_0C = model.vec[87]
obj$d0D_0E = model.vec[88]
obj$d1D_1A = model.vec[89]
obj$d1D_1B = model.vec[90]
obj$d1D_1C = model.vec[91]
obj$d1D_1E = model.vec[92]
obj$d01D_01A = model.vec[93]
obj$d01D_01B = model.vec[94]
obj$d01D_01C = model.vec[95]
obj$d01D_01E = model.vec[96]
obj$s0E = model.vec[97]
obj$s1E = model.vec[98]
obj$s01E = model.vec[99]
obj$x0E = model.vec[100]
obj$x1E = model.vec[101]
obj$x01E = model.vec[102]
obj$d0E_1E = model.vec[103]
obj$d0E_01E = model.vec[104]
obj$d1E_0E = model.vec[105]
obj$d1E_01E = model.vec[106]
if(model.vec[107] == 0){
obj$d01E_0E = model.vec[101]
}else{
obj$d01E_0E = model.vec[107]
}
if(model.vec[108] == 0){
obj$d01E_1E = model.vec[100]
}else{
obj$d01E_1E = model.vec[108]
}
obj$d0E_0A = model.vec[109]
obj$d0E_0B = model.vec[110]
obj$d0E_0C = model.vec[111]
obj$d0E_0D = model.vec[112]
obj$d1E_1A = model.vec[113]
obj$d1E_1B = model.vec[114]
obj$d1E_1C = model.vec[115]
obj$d1E_1D = model.vec[116]
obj$d01E_01A = model.vec[117]
obj$d01E_01B = model.vec[118]
obj$d01E_01C = model.vec[119]
obj$d01E_01D = model.vec[120]
obj$split.times = sort(branching.times(phy), decreasing=TRUE)
return(obj)
}
######################################################################################################################################
######################################################################################################################################
### Print function for our diversity class:
######################################################################################################################################
######################################################################################################################################
print.geohisse.old.fit <- function(x,...){
## Function to print a "geohisse.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 )
nareas <- ncol( x$trans.matrix )/3
output <- c(x$loglik, x$AIC, x$AICc, ntips, nareas)
names(output) <- c("lnL", "AIC", "AICc", "n.taxa", "n.hidden.classes")
cat("\n")
cat("Fit \n")
print(output)
cat("\n")
cat("Model parameters: \n")
cat("\n")
print(par.list)
cat("\n")
}
######################################################################################################################################
######################################################################################################################################
### Code for testing purposes:
######################################################################################################################################
######################################################################################################################################
#pars <- c(1.5, 0.5, 1.0, 0.7, 0.7, 2.5, 0.5)
#names(pars) <- diversitree:::default.argnames.geosse()
## Simulate a tree
#set.seed(5)
#phy <- tree.geosse(pars, max.t=4, x0=0)
#data <- data.frame(g_s = names(phy$tip.state), states = phy$tip.state)
## See the data
#statecols <- c("AB"="purple", "A"="blue", "B"="red")
## The likelihood function
#lik <- make.geosse(phy, phy$tip.state)
#lik(pars)
#states <- data.frame(phy$tip.state, phy$tip.state, row.names=names(phy$tip.state))
#states <- states[phy$tip.label,]
#names(pars) <- NULL
#model.vec <- numeric(95)
#tests against GeoSSE
#model.vec <- rep(c(pars[1:3], pars[4:5], pars[6:7], rep(0,12)), 5)
#model.vec <- numeric(95)
#model.vec[1:7] <- c(pars[1:3], pars[4:5], pars[6:7])
#phy$node.label <- NULL
#cache <- ParametersToPassGeoHiSSE(phy, states[,1], f=c(1,1,1), model.vec, hidden.states=FALSE)
#ll.geosse <- DownPassGeoHisse(phy=phy, cache=cache, hidden.states=FALSE, bad.likelihood=-10000000000, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL)
#cache <- ParametersToPassGeoHiSSE(phy, states[,1], f=c(1,1,1), model.vec, hidden.states=TRUE)
#ll.geohisse <- DownPassGeoHisse(phy=phy, cache=cache, hidden.states=TRUE, bad.likelihood=-10000000000, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL)
#tests against NULL
#lambda_A <- pars[1]
#lambda_B <- pars[2]
#lambda_C <- pars[3]
#model.vec <- c(rep(lambda_A, 3), pars[4:5], pars[6:7], rep(1,6), rep(lambda_B, 3), pars[4:5], pars[6:7], rep(lambda_C, 3), pars[1:3], pars[4:5], pars[6:7], rep(1,6))
#phy$node.label <- NULL
#cache <- ParametersToPassGeoSSE(phy, states[,1], f=c(1,1,1), model.vec, hidden.states=TRUE)
#ll <- DownPassGeosse(phy=phy, cache=cache, hidden.states=TRUE, bad.likelihood=-10000000000, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.