In richelbilderbeek/wiritttes: WIRITTE simulation

Goal

This vignette demonstrates the sampling used.

Introduction

options(warn = 2)

newick <- "(((A,B),A),A);"
phylogeny <- phytools::read.newick(text = newick)
ape::plot.phylo(phylogeny, edge.width = 3, cex = 2)

Demonstration

Load the packages we will need:

library(wiritttes)
library(ggplot2)
library(phytools)
library(nLTT)
library(plyr)

Here we will simulate a protracted birth-death process:

age <- 4
seed <- 320
set.seed(seed)

sirg <- 0.7 # species initiation rate good species
scr  <- 0.2 # speciation completion rate
siri <- 0.6 # species initiation rate incipient species
erg  <- 1.0 # extinction rate good species
eri  <- 0.6 # extinction rate incipient species

# Work on the pbd_sim output
pbd_sim_output <- PBD::pbd_sim(
  c(sirg, scr, siri, erg, eri),
  age, 
  plotit = TRUE
)
names(pbd_sim_output)
testit::assert(ribir::is_pbd_sim_output(pbd_sim_output))

Most of these trees are uninteresting for this demonstration. I will zoom in on the middle tree of the top row, which are the extant good and incipient species:

cols <- stats::setNames(c("gray", "black"), c("i", "g"))
phytools::plotSimmap(pbd_sim_output$igtree.extant, colors = cols)
ape::add.scale.bar()

From that tree, from each good species, we pick both the youngest and oldest:

stree_youngest <- pbd_sim_output$stree_youngest 
stree_oldest <- pbd_sim_output$stree_oldest
species_trees <- c(stree_youngest, stree_oldest)
for (species_tree in species_trees) {
  title <- ifelse(identical(species_tree, stree_youngest), "youngest", "oldest")
  ape::plot.phylo(species_tree, root.edge = TRUE, main = title)
}

The nLTT plots are expected to differ. Here we extract the nLTT plot is data:

df <- nLTT::get_nltt_values(species_trees, dt = 0.01)
df$id <- revalue(df$id, c("1" = "youngest", "2" = "oldest"))
names(df)

Here we plot both nLTTs:

ggplot2::ggplot(
  data = df,
  ggplot2::aes(x = t, y = nltt, colour = id)
) + ggplot2::geom_step(
  direction = "vh", size = 2, alpha = 0.5
) + ggplot2::scale_x_continuous(
  name = "Normalized time", limits = c(0,1)
) + ggplot2::scale_y_continuous(name = "Normalized number of lineages", limits = c(0,1)
) + ggplot2::ggtitle("nLTT plots")

Or we plot them as an average:

ggplot2::qplot(t, nltt, data = df, geom = "blank", ylim = c(0, 1), main = "Average nLTT plot of phylogenies") + 
  ggplot2::stat_summary(fun.data = "mean_cl_boot", color = "red", geom = "smooth")