R/sumtrees.R

Defines functions sumtree

Documented in sumtree

#' Compute node summaries for a set of chronos runs
#' 
#' @param trees.chronos object of class multPhylo, containing objects of class chronos
#' @param interval.type character specifying type of confidence interval (CI) calculation, one or more of c("CI95", "boot", "range", "HPD") or "all"
#' @return list containing tree and one or more CI marices
#' 
#' @details "HPD" only appropriate for posterior distributions
#' 
#' @examples
#' \donotrun{
#' chr.mult.s <- sumtree(chr.mult, "range")
#' phytools::plotTree.errorbars(chr.mult.s$phy, chr.mult.s$CI[,2:3])
#' }
#' @import Hmisc, LaplacesDemon

sumtree <- function(trees.chronos, interval.type  = "CI95"){
  
  require('ape')
  # requires Hmisc for boot
  require('Hmisc')
  # requires LaplacesDemon for HPD 
  require('LaplacesDemon')
  
  # test
  if(class(trees.chronos)!="multiPhylo")
    stop("trees.chronos must be of class multiPhylo")
  if(!inherits(trees.chronos[[1]], c("chronos", "phylo")))
    stop("trees.chrono must be of class chronos")
  
  # extra functions
  se <- function(x) sd(x)/length(x)
  ci95 <- function(x) { 
    s <- sd(x)
    error <- qnorm(0.975)*s/sqrt(length(x)) 
    return(error)
  }
  
  if(interval.type == "all"){
    interval.type <- c("CI95", "boot", "range", "HPD")
  }
  #### calc mean BLs and node heights
  bls <- do.call(rbind, lapply(trees.chronos, function(x) x$edge.length))
  nds <- do.call(rbind, lapply(trees.chronos, branching.times))
  mean_nds <- colMeans(nds)
  mean_bls <- colMeans(bls)
  
  # interval type.  options are ci95, boot, range, HPD
  if ("boot" %in% interval.type) {
    err <- apply(nds, 2, smean.cl.boot)
    boot <- cbind.data.frame(mean = mean_nds, lower = err[2,], upper = err[3,]) }
  
  if("CI95" %in% interval.type) {
    err <- apply(nds, 2, ci95)
    CI95 <- cbind.data.frame(mean = mean_nds, lower = mean_nds-err, upper = mean_nds+err) }
  
  if("range" %in% interval.type) {
    err <-  apply(nds, 2, range)
    range <- cbind.data.frame(mean = mean_nds, lower = err[1,], upper = err[2,]) }
  
  if ("HPD" %in% interval.type) {
    err <- apply(nds, 2, p.interval, HPD = TRUE, MM = FALSE, prob = 0.95)
    HPD <- cbind.data.frame(mean = mean_nds, lower = err[1,], upper = err[1,]) }
  
  #### allocate mean values to a random tree
  tr <- trees.chronos[[1]]
  # class(tr) <- "phylo"
  tr$edge.length <- mean_bls
  
  return(list(phy = tr, CI = mget(interval.type)))
}
FranzKrah/seta documentation built on Sept. 15, 2017, 4:50 p.m.