R/get_marginal.r

In reslp/correlate: Character correlation with phylogenetic trees

#' Calculate the marginal probability of character occurrence.
#' 
#' This function calculates the marginal probability of observing a particular character on the provided phylogenetic tree(s). Will return a dataframe containing the marginal probabilities of character occurrence for each character.
#' @param smap List of multiPhylo objects containing stochastic character maps created with phytools.
#' @examples
#' get_marginal(smap)
get_marginal <- function(sm){
  #function get_marginal, will extract P(character) from the provided stochastic maps per tree
  #needed for later calculation of conditional
  if (class(sm) != "list") { #single trees, behaviour unclear for otherwise formatted phylo objects
    sim <- vector("list", 1)
    sim[[1]] <- sm
    sm <- sim
  }
  ntrees <- length(sm)
  nmaps <- length(sm[[1]])
  time_matrix <- vector()
  edge_length <- vector()
  time_summary <- vector()
  characters <- colnames(sm[[1]][[1]]$mapped.edge)
  df <- data.frame(matrix(0,nrow=ntrees,ncol=length(characters)))
  whichchar <- 0
  for (c in characters) {
    whichchar <- whichchar + 1
    colnames(df)[whichchar] <- c
    for (i in 1:ntrees) {
      edge_length <- sum(sm[[i]][[1]]$edge.length)
      time_matrix <- vector()
      for (j in 1:nmaps) {
        time_matrix  <- c(time_matrix,sum(as.numeric(sm[[i]][[j]]$mapped.edge[,whichchar]))) #1=character 0, 2=character 1,...
      }
      time_summary <- c(time_summary,(sum(time_matrix,na.rm=T)/nmaps)) ## maybe na.rm=T NA needed in sum
      df[i,whichchar] <- sum(time_summary/edge_length)
      time_matrix <- vector()
      time_summary <- vector()
    }

  }
  return(df)
}

iseven <- function(x) x %% 2 == 0