######################################################################################################################################
######################################################################################################################################
### Model Average Rates -- calculates a model-averaged rates for each observed state under a set of parameters
######################################################################################################################################
######################################################################################################################################
GetModelAveRates <- function(x, AIC.weights=NULL, type=c("tips", "nodes", "both"), bound.par.matrix=cbind(c(0,-1000,0,0,0),c(1000,1000,1000,3,1000))){
hisse.results <- x
if (!type %in% c("tips", "nodes", "both")) stop("Argument 'type' needs to be one of 'tips', 'nodes', or 'both'.")
## Check if the format of bound par matrix is correct.
## This needs to be a numeric matrix with two columns and 5 rows.
## First column are minimums and second colum are maximums.
if( !inherits(bound.par.matrix, what = c("matrix","data.frame")) ){
stop(" 'bound.par.matrix' needs to be a matrix or data.frame ")
}
if( !ncol(bound.par.matrix) == 2 | !nrow(bound.par.matrix) == 5 ){
stop( " 'bound.par.matrix' needs to be a matrix with 5 rows and 2 columns. See help.")
}
## Beware that the test below might misbehave if the lengths are different.
## But since we tested the dimensions, then this seems fine.
if( !all( bound.par.matrix[,1] < bound.par.matrix[,2] ) ){
stop("Min bounds larger than max! First column of 'bound.par.matrix' are the minimum bounds and second column are the maximum.")
}
## Create flag to compute the model average for the both nodes and tips or just one of them.
if( type == "both" ){
to.mod.ave <- c("tips", "nodes")
} else{
to.mod.ave <- type
}
if( !inherits(hisse.results, what = c("list", "hisse.states", "hisse.geosse.states", "muhisse.states", "misse.states")) ) stop("x needs to be a list of model reconstructions or a single model reconstruction object of class 'hisse.states', 'hisse.geosse.states', 'muhisse.states', or 'misse.states'.")
## If hisse.results is a list of model reconstructions, then test if they have $AIC. Return error message otherwise.
## There is no need for the $AIC element if AIC.weigths argument is provided.
## AIC.weights is a substitute for AIC. See help page.
if(inherits(hisse.results, what="list")){
if( is.null( AIC.weights ) ){
empty.AIC <- sapply(hisse.results, function(x) !is.null(x$AIC) )
if(sum( empty.AIC ) != length(hisse.results)) {
for (hisse.results.index in seq_along(hisse.results)) {
if(is.null(hisse.results[[hisse.results.index]]$AIC) & !is.null(hisse.results[[hisse.results.index]]$aic)) {
hisse.results[[hisse.results.index]]$AIC <- hisse.results[[hisse.results.index]]$aic
}
}
empty.AIC <- sapply(hisse.results, function(x) !is.null(x$AIC) )
}
if( sum( empty.AIC ) != length(hisse.results) ) stop("All elements of the list need to have a '$AIC' element.")
model.class <- sapply(hisse.results, function(x) !inherits(x, what = c("hisse.states", "hisse.geosse.states", "muhisse.states", "misse.states")) )
if( as.logical( sum( model.class ) ) ) stop("x needs to be a list of model reconstruction with class 'hisse.states' or 'hisse.geosse.states' ")
} else{
if( !length(hisse.results) == length(AIC.weights) ){
stop( " AIC.weights needs to be NULL or a numeric vector with length equal to the number of models in 'x'. " )
}
}
}
if( inherits(hisse.results, what = c("hisse.states", "hisse.geosse.states", "muhisse.states", "misse.states")) ){ # we have to make a list so we can run this generally
if( !is.null( AIC.weights ) ){
stop( "If using a single model, then 'AIC.weights' needs to be NULL. " )
}
if(is.null(hisse.results$AIC)){
if(!is.null(hisse.results$aic)) {
hisse.results$AIC <- hisse.results$aic #fix capitalization
} else {
## If a user forgot to include the AIC, then we add a random value in for them
hisse.results$AIC = 42
}
}
tmp.list <- list()
tmp.list[[1]] <- hisse.results
hisse.results <- tmp.list
}
## This will compute the AIC weights just once now. ;)
if( is.null(AIC.weights) ){
AIC.weights <- GetAICWeights(hisse.results)
}
## ##################################################################################################################
## Compute the model average for the nodes:
if( "tips" %in% to.mod.ave ){
rates.tips.turnover <- ConvertManyToRate_ModelAve(hisse.results, rate.param="turnover", which.element="tip.mat")
rates.tips.net.div <- ConvertManyToRate_ModelAve(hisse.results, rate.param="net.div", which.element="tip.mat")
rates.tips.speciation <- ConvertManyToRate_ModelAve(hisse.results, rate.param="speciation", which.element="tip.mat")
rates.tips.extinct.fraction <- ConvertManyToRate_ModelAve(hisse.results, rate.param="extinction.fraction", which.element="tip.mat")
rates.tips.extinction <- ConvertManyToRate_ModelAve(hisse.results, rate.param="extinction", which.element="tip.mat")
## Here 'averaged.rates' is a list with the vector of results.
averaged.tip.rates <- CheckReconBounds(x=list(rates.tips.turnover, rates.tips.net.div, rates.tips.speciation, rates.tips.extinct.fraction, rates.tips.extinction), n.models = length(hisse.results), AIC.weights, bound.par.matrix)
## Objects will always be of list class here.
if(inherits(hisse.results[[1]], what="hisse.states")){
states.tips <- ConvertManyToBinaryState(hisse.results, which.element="tip.mat", AIC.weights=AIC.weights)
}
if(inherits(hisse.results[[1]], what="hisse.geosse.states")){
states.tips <- ConvertManyToMultiState(hisse.results, which.element="tip.mat", AIC.weights=AIC.weights)
}
if(inherits(hisse.results[[1]], what="muhisse.states")){
states.tips <- ConvertManyToMultiState(hisse.results, which.element="tip.mat", AIC.weights=AIC.weights)
}
if(inherits(hisse.results[[1]], what="misse.states")){
states.tips <- rep(0, dim(hisse.results[[1]][["tip.mat"]])[1])
}
## Note that the columns of 'averaged.tip.rates' need to be in the correct order here. They should be.
final.df.tips <- data.frame(taxon=hisse.results[[1]]$phy$tip.label, state=states.tips, turnover=averaged.tip.rates[[1]], net.div=averaged.tip.rates[[2]], speciation=averaged.tip.rates[[3]], extinct.frac=averaged.tip.rates[[4]], extinction=averaged.tip.rates[[5]])
## Check if the function need to return only the reconstruction for the tips.
if( type == "tips" ) return(final.df.tips)
}
## ##################################################################################################################
## Compute the model average for the nodes:
if( "nodes" %in% to.mod.ave ){
rates.nodes.turnover <- ConvertManyToRate_ModelAve(hisse.results, rate.param="turnover", which.element="node.mat")
rates.nodes.net.div <- ConvertManyToRate_ModelAve(hisse.results, rate.param="net.div", which.element="node.mat")
rates.nodes.speciation <- ConvertManyToRate_ModelAve(hisse.results, rate.param="speciation", which.element="node.mat")
rates.nodes.extinct.fraction <- ConvertManyToRate_ModelAve(hisse.results, rate.param="extinction.fraction", which.element="node.mat")
rates.nodes.extinction <- ConvertManyToRate_ModelAve(hisse.results, rate.param="extinction", which.element="node.mat")
## Here 'averaged.node.rates' is a list with the vector of results.
averaged.node.rates <- CheckReconBounds(x=list(rates.nodes.turnover, rates.nodes.net.div, rates.nodes.speciation, rates.nodes.extinct.fraction, rates.nodes.extinction), n.models = length(hisse.results), AIC.weights, bound.par.matrix)
## Objects will always be of list class here.
if(inherits(hisse.results[[1]], what="hisse.states")){
states.internal <- ConvertManyToBinaryState(hisse.results, which.element="node.mat", AIC.weights=AIC.weights)
}
if(inherits(hisse.results[[1]], what="hisse.geosse.states")){
states.internal <- ConvertManyToMultiState(hisse.results, which.element="node.mat", AIC.weights=AIC.weights)
}
if(inherits(hisse.results[[1]], what="muhisse.states")){
states.internal <- ConvertManyToMultiState(hisse.results, which.element="node.mat", AIC.weights=AIC.weights)
}
if(inherits(hisse.results[[1]], what="misse.states")){
states.internal <- rep(0, dim(hisse.results[[1]][["node.mat"]])[1])
}
final.df.nodes <- data.frame(id=hisse.results[[1]]$node.mat[,1], state=states.internal, turnover=averaged.node.rates[[1]], net.div=averaged.node.rates[[2]], speciation=averaged.node.rates[[3]], extinct.frac=averaged.node.rates[[4]], extinction=averaged.node.rates[[5]])
## Check if the function need to return only the reconstruction for the nodes:
if( type == "nodes" ) return(final.df.nodes)
}
## If both node and tips need to be returned:
if( type == "both" ){
final.both <- list(tips = final.df.tips, nodes = final.df.nodes)
return( final.both )
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.