R/sim.tip.sampling.R

In FossilSim: Simulation of Fossil and Taxonomy Data

#' Include extant samples in the fossil object, with optional rho sampling.
#'
#' @param fossils Fossils object.
#' @param tree Phylo object.
#' @param taxonomy Taxonomy object.
#' @param rho Extant species sampling probability.
#' @param tol Rounding error tolerance for tip ages.
#'
#' @return An object of class fossils containing extant tip samples equal to the age of the tips (i.e. 0.0).
#'
#' @examples
#' # simulate tree
#' lambda = 0.1
#' mu = 0.05
#' tips = 8
#' t = TreeSim::sim.bd.taxa(tips, 1, lambda, mu)[[1]]
#'
#' # simulate fossils
#' f = sim.fossils.poisson(0.5, t)
#'
#' # simulate extant samples
#' f = sim.extant.samples(f, t, rho = 0.5)
#' plot(f, t)
#'
#' @export
sim.extant.samples = function(fossils, tree = NULL, taxonomy = NULL, rho = 1, tol = NULL){

  if(!"fossils" %in% class(fossils))
    stop("fossils must be an object of class \"fossils\"")

  if(is.null(tree) && is.null(taxonomy))
    stop("Specify phylo or taxonomy object")

  if(!is.null(tree) && !"phylo" %in% class(tree))
    stop("tree must be an object of class \"phylo\"")

  if(!is.null(taxonomy) && !"taxonomy" %in% class(taxonomy))
    stop("taxonomy must be an object of class \"taxonomy\"")

  if(!is.null(tree) && !is.null(taxonomy))
    warning("tree and taxonomy both defined, using taxonomy")

  if(!is.null(attr(fossils, "from.taxonomy"))){
    from.taxonomy = attr(fossils, "from.taxonomy")
    if(!is.null(taxonomy) & !from.taxonomy)
      stop("species taxonomy defined but fossils not based on taxonomy")
  }

  if(!(rho >= 0 && rho <= 1))
    stop("rho must be a probability between 0 and 1")

  if(is.null(taxonomy)){
    taxonomy = sim.taxonomy(tree, beta = 1)
    from.taxonomy = FALSE
  } else from.taxonomy = TRUE

  if(is.null(tol))
    tol = if(is.null(tree)) 1e-8 else min(min(tree$edge.length)/100, 1e-8)

  for (i in unique(taxonomy$sp)){

    end = min(taxonomy$end[which(taxonomy$sp == i)])

    if(!(end > (0 - tol) & end < (0 + tol))) next

    if(runif(1) < rho){
      # identify the edge ending zero
      edge = taxonomy$edge[which(abs(taxonomy$end) < tol & taxonomy$sp == i)]

      fossils<-rbind(fossils, data.frame(sp = i, edge = edge, hmin = 0, hmax = 0))
    }

  }
  if(!is.fossils(fossils))
    fossils = as.fossils(fossils, from.taxonomy = from.taxonomy)
  return(fossils)
  #eof
}

#' Include extant and extinct tip samples in the fossil object, with optional rho sampling.
#'
#' @param fossils Fossils object.
#' @param tree Phylo object.
#' @param taxonomy Taxonomy object.
#' @param rho Tip sampling probability.
#'
#' @return An object of class fossils containing extant or extinct tip samples equal to the age of the tips.
#'
#' @examples
#' # simulate tree
#' t = ape::rtree(6)
#'
#' # simulate fossils
#' f = sim.fossils.poisson(2, t)
#'
#' # simulate tip samples
#' f = sim.tip.samples(f, t, rho = 0.5)
#' plot(f, t)
#'
#' @export
# The tree is required to identify which edges are terminal.
sim.tip.samples<-function(fossils, tree, taxonomy = NULL, rho = 1){

  if(!"fossils" %in% class(fossils))
    stop("fossils must be an object of class \"fossils\"")

  if(!"phylo" %in% class(tree))
    stop("tree must be an object of class \"phylo\"")

  if(!is.null(taxonomy) && !"taxonomy" %in% class(taxonomy))
    stop("taxonomy must be an object of class \"taxonomy\"")

  if(!is.null(tree) && !is.null(taxonomy))
    warning("tree and taxonomy both defined, using taxonomy")

  if(!is.null(attr(fossils, "from.taxonomy"))){
    from.taxonomy = attr(fossils, "from.taxonomy")
    if(!is.null(taxonomy) & !from.taxonomy)
      stop("species taxonomy defined but fossils not based on taxonomy")
  }

  if(!(rho >= 0 && rho <= 1))
    stop("rho must be a probability between 0 and 1")

  if(is.null(taxonomy)){
    taxonomy = sim.taxonomy(tree, beta = 1)
    from.taxonomy = FALSE
  } else from.taxonomy = TRUE

  for (i in unique(taxonomy$sp)){

    # identify the terminal most edge
    end = min(taxonomy$end[which(taxonomy$sp == i)])
    edge = taxonomy$edge[which(taxonomy$sp == i & taxonomy$end == end)]

    if(is.tip(edge,tree)){

      if(runif(1) < rho){
        fossils<-rbind(fossils, data.frame(sp = i, edge = edge, hmin = end, hmax = end))
      }
    }

  }
  if(!is.fossils(fossils))
    fossils = as.fossils(fossils, from.taxonomy = from.taxonomy)
  return(fossils)
}