Trees comparison

s <- sim_phyl()
plot(s$newick)

p <- subplex(par = c(8,0.175,0.9),fn = llik,n = s$n, E = s$E, t = s$t)
pa=c(p$par[1],p$par[2],p$par[3])
dropex <- drop.fossil(s$newick)
plot(dropex)

s2 <- phylo2p(tree=dropex)

a = data.frame(t=seq(0,15,by=0.1))
a[[2]] = approx(cumsum(s$t),s$n,xou=seq(0,15,by=0.1))$y
n_it = 100
estimations=F
if(estimations) P = matrix(nrow=n_it,ncol=4)
S = vector('list',length=n_it)
for (i in 1:n_it){
  rt = rec_tree(wt=s2$t,pars=pa)
  a[[i+2]] = approx(cumsum(rt$wt),rt$n,xou=seq(0,15,by=0.1))$y
  if(estimations){
    pars = subplex(par = c(8,0.175,0.9),fn = llik,n = rt$n, E = rt$E, t = rt$t)$par
    pars2 = c(pars[1],pars[2],pars[3],(pars[1]-pars[3])/pars[2])
    P[i,] = pars2}
  S[[i]] = rt
}

melted = melt(a, id.vars="t")
p <- ggplot(data=melted, aes(x=t, y=value, group=variable)) + geom_line() + ylab('# Lineages') + xlab('Time') + geom_line(data=melted[melted$variable=='V2',],aes(x=t, y=value),color='blue') +ggtitle('(mle)')#  + geom_line(data=data.frame(t=seq(0,15,by=0.1),variable="mean",value=rowMeans(AA, na.rm = TRUE, dims = 1)),aes(x=t, y=value),color='green')
print(p)

## Warning: Removed 828 rows containing missing values (geom_path).

## Warning: Removed 13 rows containing missing values (geom_path).

#make it a function
a = data.frame(t=seq(0,15,by=0.1))
a[[2]] = approx(cumsum(s$t),s$n,xou=seq(0,15,by=0.1))$y
n_it = 100
for (i in 1:n_it){
  s = sim_phyl()
  a[[i+2]] = approx(cumsum(s$t),s$n,xou=seq(0,15,by=0.1))$y
}
a[is.na(a)]=0
b = rowMeans(a[,-1])
c = data.frame(t=a[,1],species=b)
qplot(data=c,x=t,y=species)

another interestin way of comparison is using the nLTT package:

library(nLTT)
data(exampleTrees);
obs <- exampleTrees[[1]];
nltt_plot(obs);

# Obtain the nLTT values
dt <- 0.01
phylogenies = exampleTrees
nltt_values <- get_nltt_values(phylogenies, dt = dt)
# Plot the phylognies, where the individual nLTT values are visible
qplot(t, nltt, data = nltt_values, geom = "point",
      ylim = c(0,1),
      main = "Average nLTT plot of phylogenies", color = id, size = I(0.1)
) + stat_summary(
  fun.data = "mean_cl_boot", color = "red", geom = "smooth"
)

nLTTstat(exampleTrees[[1]], exampleTrees[[2]], distance_method = "abs")

## [1] 0.1070197