search

PoonLab / Kaphi

Kernel-embedded ABC-SMC for phylodynamic inference
GNU Affero General Public License v3.0
4 stars 2 forks source link

Trip is not working with the BiSSE model #137

Closed helenhe96 closed 6 years ago

helenhe96 commented 6 years ago

When I run this:

config <- load.config('tests/fixtures/test-bisse.yaml')
config <- set.model(config, 'bisse')
theta <- c(lambda0=0.2, lambda1=0.1, mu0=0.003, mu1=0.003, q01=0.03, q10=0.01)  # this is the true value
set.seed(50)
obs.tree <- speciation.model(theta, nsim=1, tips=50, model='bisse')[[1]]
obs.tree <- parse.input.tree(obs.tree, config)
theta <- c(lambda0=0.1, lambda1=0.1, mu0=0.003, mu1=0.003, q01=0.03, q10=0.01)
sim.tree <- speciation.model(theta, nsim=1, tips=50, model='bisse')[[1]]
Trip(obs.tree,sim.tree)

It produces this error:

 Error in if (nd == rootnd) break : missing value where TRUE/FALSE needed 
5.
getMRCA(y, pair) 
4.
Kaphi::Trip(x, y) at <text>#1
3.
eval(parse(text = config$dist)) 
2.
eval(parse(text = config$dist)) 
1.
distance(obs.tree, sim.tree, config)
gtng92 commented 6 years ago

Previously noted under Issue #100 that Trip/TripL did not work under BD model. Trees generated by package diversitree do not necessarily generate the same tip labels. I don't know if this was intentional or not. Trip and TripL metrics require the two given trees to have the same number and same set of tip labels.

Trees generated w/ 10 tips:

> t1 <- speciation.model(theta, nsim=2, tips=10, model='bisse')[[1]]
> t2 <- speciation.model(theta, nsim=2, tips=10, model='bisse')[[2]]
> t1$tip.label
 [1] "sp1_0"  "sp2_0"  "sp3_0"  "sp4_0"  "sp5_0"  "sp6_0"  "sp7_0"  "sp8_0"  "sp9_1"  "sp10_1"
> t2$tip.label
 [1] "sp1_0"  "sp2_0"  "sp3_0"  "sp4_0"  "sp5_0"  "sp6_0"  "sp7_0"  "sp8_0"  "sp9_0"  "sp10_0"

@ArtPoon added the '_0' and '_1' tip states to the simulated trees within diversitree-wrappers.R, but even when ignoring those appendages, trees with a higher number of tips do not necessarily have the same set of tip labels.

Trees generated w/ 50 tips:

> t1 <- speciation.model(theta, nsim=1, tips=50, model='bisse')[[1]]
> t2 <- speciation.model(theta, nsim=2, tips=50, model='bisse')[[2]]
> t1$tip.label
 [1] "sp1_1"  "sp2_1"  "sp3_1"  "sp4_0"  "sp5_1"  "sp6_0"  "sp7_1"  "sp8_1"  "sp9_0"  "sp10_0" "sp11_1" "sp12_1"
[13] "sp13_1" "sp14_1" "sp15_1" "sp16_0" "sp17_0" "sp18_0" "sp19_0" "sp20_0" "sp21_1" "sp22_0" "sp23_0" "sp24_0"
[25] "sp25_0" "sp26_0" "sp27_0" "sp28_0" "sp29_1" "sp30_0" "sp31_0" "sp32_0" "sp33_0" "sp34_0" "sp35_0" "sp36_0"
[37] "sp37_0" "sp38_0" "sp39_0" "sp40_0" "sp41_0" "sp42_0" "sp43_0" "sp44_1" "sp45_0" "sp46_0" "sp47_1" "sp48_1"
[49] "sp49_1" "sp50_1"
> t2$tip.label
 [1] "sp1_0"  "sp2_1"  "sp3_1"  "sp4_0"  "sp5_0"  "sp6_0"  "sp8_1"  "sp10_0" "sp11_0" "sp12_1" "sp13_0" "sp14_1"
[13] "sp15_0" "sp16_0" "sp17_0" "sp18_0" "sp19_0" "sp20_0" "sp21_0" "sp22_0" "sp23_0" "sp24_0" "sp25_0" "sp27_1"
[25] "sp28_1" "sp29_1" "sp30_0" "sp31_0" "sp32_0" "sp33_0" "sp34_0" "sp35_0" "sp36_0" "sp37_0" "sp38_0" "sp39_0"
[37] "sp40_0" "sp41_0" "sp42_1" "sp43_0" "sp44_0" "sp45_0" "sp46_0" "sp47_0" "sp48_0" "sp49_0" "sp50_0" "sp51_0"
[49] "sp52_0" "sp53_0"

I will be inserting a check that the tip labels for tree1 and tree2 in the Trip metric are the same so the error message is more explicit.