In franciscorichter/dmea02: Diversification and Macro-Evolution Analysis

knitr::opts_chunk$set(echo = TRUE)
library(dmea02)

tree simulation and consistency

functions

1. sim.tree to simulate trees
2. phylo2vectors to convert a phylo into the set of vectors characterizing the tree

s = sim.tree(ct=6)
plot(s$phylo)
s$tree
phylo2vectors(s$phylo)
all.equal(s$tree$wt,phylo2vectors(s$phylo)$wt)
all.equal(s$tree$E,phylo2vectors(s$phylo)$E)
all.equal(s$tree$n,phylo2vectors(s$phylo)$n)

to do:

1. phylo2vectorsdoes not give the vector S corresponding to topology. add it
2. remove the root on vector outputs, keep consistency
3. arenge the vectors2phylo function and make a lot simpler the sim.tree using it
4. add the log in an log output file

Tree manipulation

s = sim.tree(ct=2,seed=5)
plot(s$phylo)
s2=phylo2vectors(s$phylo)
s2
s$tree
up = update.tree(s2,0.5,1.5)
up

1. Doing update tree I realized that I should be more carefully in analysis when we lost crown time.
2. phylo2vectors really needs to take into account topology

more about consistency

s <- sim.tree()
plot(s$phylo)
plot(vectors2phylo(s$tree))
all.equal(s$phylo,vectors2phylo(s$tree))

The reconstruction algorithm

s = sim.tree(ct=6)
plot(s$phylo.extant)
rec = rec.tree(tree=s$tree.extant,pars=c(0.8,0.1,40))
plot(vectors2phylo(rec))
s$tree.extant
rec

let´s observe one estimation

s = sim.tree(seed = 1)
mle.tree(s$tree)
tree = s$tree.extant
st = sim.srt(tree,pars=c(0.8,0.1,40),n_trees=100)
mle.st(st)

How good is the last iteration of the MCEM algorithm?

library(dmea)
n_sim = 537
n_trees = 10
MP = matrix(nrow=n_sim,ncol=3)
RP = matrix(nrow=n_sim,ncol=3)
p = proc.time()
for(i in 1:n_sim){
  est = sim.est(n_trees=n_trees,pars=c(0.8,0.1,40),seed=i)
  RP[i,] = est$real
  MP[i,] = est$est
}
print(proc.time()-p)
par_est_vis(P=MP,par=1,PR=RP)
par_est_vis(P=MP,par=2,PR=RP)
par_est_vis(P=MP,par=3,PR=RP)

franciscorichter/dmea02 documentation built on March 21, 2020, 3:46 a.m.