shape_purebirth.md

Install libraries if you don't have them

install.packages(c("phytools","apTreeshape","ggplot2","plyr"))

Load libraries

library(phytools)
library(apTreeshape)
library(ggplot2)
library(plyr)

Define functions to generate trees and calculate metrics

simtrees_colless <- function(ntips=10, ntrees=100){
  tmp <- pbtree(n=ntips, b = 1, d=0, scale=1, nsim=ntrees)
  res <- vapply(tmp, function(z){
    colless(as.treeshape(z), norm="yule")
  }, 1)
  names(res) <- rep(ntips, length(res))
  res
}
simtrees_gamma <- function(ntips=10, ntrees=100){
  tmp <- pbtree(n=ntips, b = 1, d=0, scale=1, nsim=ntrees)
  res <- vapply(tmp, function(z) gammaStat(z), 1)
  names(res) <- rep(ntips, length(res))
  res
}

Define vector of sizes for trees

treesizes <- seq(from = 10, to = 750, by = 50)

Simulate trees and calculate colless values

out <- lapply(treesizes, simtrees_colless, ntrees=100)
names(out) <- treesizes
df <- do.call(rbind, lapply(out, function(colless) {
  tt <- data.frame(colless)
  tt$numtips <- as.numeric(names(colless)[[1]])
  tt
}))

Simulate trees and calculate gamma values

out_gamma <- lapply(treesizes, simtrees_gamma, ntrees=100)
names(out_gamma) <- treesizes
df_gamma <- do.call(rbind, lapply(out_gamma, function(gamma) {
  tt <- data.frame(gamma)
  tt$numtips <- as.numeric(names(gamma)[[1]])
  tt
}))

Plot of colless values

note: Zero is the expectation line

ggplot(df, aes(x = factor(numtips), y = colless)) +
  geom_boxplot() +
  geom_jitter(position = position_jitter(width = 0.1)) + 
  labs(x="Number of tips in phylogeny", y="Colless metric (Ic)") +
  theme_grey(base_size = 18)

Plot of gamma values

note: Zero is the expectation line

ggplot(df_gamma, aes(x = factor(numtips), y = gamma)) +
  geom_boxplot() +
  geom_jitter(position = position_jitter(width = 0.1)) + 
  labs(x="Number of tips in phylogeny", y="Gamma metric") +
  theme_grey(base_size = 18)