equal likelihood trees

st <- sim_phyl(ct=5)
plot(st$newick,root.edge = T)
axisPhylo(backward = F,root.time = st$t[1])

First, to warm up, we plot a simple tree

obs = '((AB:1,AA:1):1,AC:2):1;'
#library(ape)
obs = read.tree(text=obs)
is.rooted(obs)

## [1] TRUE

is.ultrametric(obs)

## [1] TRUE

plot(obs,show.tip.label = T, root.edge = TRUE)
axisPhylo(backward = F,root.time = 1)

branching.times(obs)

## 4 5 
## 2 1

Then a longer one

# complete/incomplete trees
comp = '((AB:1,(AA:0.6,AD:0.3):0.4):1,((AE:0.6,AF:0.8):0.2,AC:1.5):0.5):1;'
comp= read.tree(text = comp)
is.rooted(comp)

## [1] TRUE

is.ultrametric(comp)

## [1] FALSE

plot(comp,show.tip.label = T, root.edge = TRUE)
axisPhylo(backward = F,root.time = 1)

branching.times(comp)

##   7   8   9  10  11 
## 2.0 1.0 0.6 1.5 1.3

Then we drop extinxt species

incomp = drop.fossil(comp)
write.tree(incomp)

## [1] "((AB:1,AA:1):1,AC:2):1;"

is.rooted(incomp)

## [1] TRUE

is.ultrametric(incomp)

## [1] TRUE

plot(incomp,show.tip.label = T, root.edge = TRUE)
axisPhylo(backward = F,root.time = 1)

branching.times(incomp)

## 4 5 
## 2 1

Now we would like to be able to add an extinct species

in_dmea = phylo2p(incomp)
in_dmea$wt

## [1] 1 1 1

in_dmea$E

## [1] 1 1

in_dmea$n

## [1] 1 2 3

ph2 <- update_tree(in_dmea,t_spe=0.5,t_ext = 1.5)
ph2o = p2phylo(ph2)
plot(ph2o,root.edge=T)
axisPhylo(backward = F,root.time = 0.5)

ltt.plot.coords(ph2o)

##      time N
## [1,] -3.0 1
## [2,] -2.5 1
## [3,] -2.0 2
## [4,] -1.5 3
## [5,]  0.0 4

now we would like to see equivalent trees

ph2e = p2phylo(ph2)
plot(ph2e,root.edge=T)
axisPhylo(backward = F,root.time = 0.5)

ltt.plot.coords(ph2e)

##      time N
## [1,] -3.0 1
## [2,] -2.5 1
## [3,] -2.0 2
## [4,] -1.5 3
## [5,]  0.0 4

ph2e = p2phylo(ph2)
plot(ph2e,root.edge=T)
axisPhylo(backward = F,root.time = 0.5)

ltt.plot.coords(ph2e)

##      time N
## [1,] -3.0 1
## [2,] -2.5 1
## [3,] -2.0 2
## [4,] -1.5 3
## [5,]  0.0 4

ph2e = p2phylo(ph2)
plot(ph2e,root.edge=T)
axisPhylo(backward = F,root.time = 0.5)

ltt.plot.coords(ph2e)

##      time N
## [1,] -3.0 1
## [2,] -2.5 1
## [3,] -2.0 2
## [4,] -1.5 3
## [5,]  0.0 4

plot(ph2e,root.edge=T)
axisPhylo(backward = F,root.time = 0.5)

ltt.plot.coords(ph2e)

##      time N
## [1,] -3.0 1
## [2,] -2.5 1
## [3,] -2.0 2
## [4,] -1.5 3
## [5,]  0.0 4

ph2e = p2phylo(ph2)
plot(ph2e,root.edge=T)
axisPhylo(backward = F,root.time = 0.5)

ltt.plot.coords(ph2e)

##      time N
## [1,] -3.0 1
## [2,] -2.5 1
## [3,] -2.0 2
## [4,] -1.5 3
## [5,]  0.0 4

s = sim_phyl(ct=6)
plot(s$newick)

p = phylo2p(s$newick)
p$E == s$E

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE

p$wt == s$t

## logical(0)

#p$wt
#s$t
p$n == s$n

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE

s2 = p2phylo(p)
plot(s2)

voilà!, we found two things:

This algorithm does not label species in a general way yet, but only to models where all species has same probability to speciate.
We found a nice way to ilustrate ‘equivalent’ trees.

s = sim_phyl(ct=4)
plot(s$newick)

s2=p2phylo(phylo2p(s$newick))
plot(s2)

ltt.plot.coords(s2)

##             time N
##  [1,] -4.0000000 1
##  [2,] -1.9810105 1
##  [3,] -1.4706033 2
##  [4,] -1.3859856 3
##  [5,] -1.0299193 4
##  [6,] -0.5139266 5
##  [7,] -0.2843561 6
##  [8,] -0.1934067 7
##  [9,]  0.0000000 8

nLTTstat(s2,s$newick)

## [1] 2.63678e-16

s2=p2phylo(phylo2p(s$newick))
plot(s2)

ltt.plot.coords(s2)

##             time N
##  [1,] -4.0000000 1
##  [2,] -1.9810105 1
##  [3,] -1.4706033 2
##  [4,] -1.3859856 3
##  [5,] -1.0299193 4
##  [6,] -0.5139266 5
##  [7,] -0.2843561 6
##  [8,] -0.1934067 7
##  [9,]  0.0000000 8

nLTTstat(s2,s$newick)

## [1] 1.19349e-15

s2=p2phylo(phylo2p(s$newick))
plot(s2)

ltt.plot.coords(s2)

##             time N
##  [1,] -4.0000000 1
##  [2,] -1.9810105 1
##  [3,] -1.4706033 2
##  [4,] -1.3859856 3
##  [5,] -1.0299193 4
##  [6,] -0.5139266 5
##  [7,] -0.2843561 6
##  [8,] -0.1934067 7
##  [9,]  0.0000000 8

nLTTstat(s2,s$newick)

## [1] 5.273559e-16

s2=p2phylo(phylo2p(s$newick))
plot(s2)

ltt.plot.coords(s2)

##             time N
##  [1,] -4.0000000 1
##  [2,] -1.9810105 1
##  [3,] -1.4706033 2
##  [4,] -1.3859856 3
##  [5,] -1.0299193 4
##  [6,] -0.5139266 5
##  [7,] -0.2843561 6
##  [8,] -0.1934067 7
##  [9,]  0.0000000 8

nLTTstat(s2,s$newick)

## [1] 1.151856e-15

s2=p2phylo(phylo2p(s$newick))
plot(s2)

ltt.plot.coords(s2)

##             time N
##  [1,] -4.0000000 1
##  [2,] -1.9810105 1
##  [3,] -1.4706033 2
##  [4,] -1.3859856 3
##  [5,] -1.0299193 4
##  [6,] -0.5139266 5
##  [7,] -0.2843561 6
##  [8,] -0.1934067 7
##  [9,]  0.0000000 8

nLTTstat(s2,s$newick)

## [1] 1.165734e-15

s2=p2phylo(phylo2p(s$newick))
plot(s2)

ltt.plot.coords(s2)

##             time N
##  [1,] -4.0000000 1
##  [2,] -1.9810105 1
##  [3,] -1.4706033 2
##  [4,] -1.3859856 3
##  [5,] -1.0299193 4
##  [6,] -0.5139266 5
##  [7,] -0.2843561 6
##  [8,] -0.1934067 7
##  [9,]  0.0000000 8

nLTTstat(s2,s$newick)

## [1] 1.151856e-15

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##         9        10        11        12        13        14        15 
## 1.9810105 1.4706033 1.3859856 1.0299193 0.2843561 0.5139266 0.1934067

nLTTstat(s2,s$newick)

## [1] 0

plot(s$newick)

ltt.plot.coords(s2)

##             time N
##  [1,] -4.0000000 1
##  [2,] -1.9810105 1
##  [3,] -1.4706033 2
##  [4,] -1.3859856 3
##  [5,] -1.0299193 4
##  [6,] -0.5139266 5
##  [7,] -0.2843561 6
##  [8,] -0.1934067 7
##  [9,]  0.0000000 8

Why the branching times are different ??
Why on those cases the nLTT > 0 ??