Nodes getting different labels in get.slide.freq.table for make plot.asdsf?

[lukejkelly]

I am getting some unexpected behaviour when plotting multiple chains which have mostly converged to the same tree, after discarding lots of burn-in the cumulative ASDSF is still around 10% so I am wondering if it is being caused by the following behaviour.

In using the debugger, it seems like ape:::prop.part uses the first tree in a window for ordering the labels whereas rwty:::get.slide.freq.table uses the first in the entire list of trees so node indices are not necessarily the same between windows.

This example produces similar behaviour. The tree t2 has leaves b and d swapped and appears once in each window in each chain.

library("ape")
library("rwty")

t1 <- read.tree(text = "(a:2,(b:1,c:1):1,(d:1,e:1):1);", keep.multi = TRUE)
t2 <- read.tree(text = "(a:2,(d:1,c:1):1,(b:1,e:1):1);", keep.multi = TRUE)

f1 <- list(trees = c(t1, t1, t2, t1, t1, t2, t1, t1), ptable =  NULL, gens.per.tree = 1)
f2 <- list(trees = c(t1, t1, t1, t2, t2, t1, t1, t1), ptable =  NULL, gens.per.tree = 1)

class(f1) <- "rwty.chain"
class(f2) <- "rwty.chain"

debug(rwty:::get.slide.freq.table)
makeplot.asdsf(list(f1, f2), 0, 4)

For processing the first chain, everything is fine. For the second chain, the clade tables inside clade.freq for the windows are as follows:

1: a
2: b
3: c
4: d
5: e

==> 4 time(s):[1] 1 2 3 4 5
==> 3 time(s):[1] 1 4 5
==> 3 time(s):[1] 1 2 3
==> 1 time(s):[1] 1 2 5
==> 1 time(s):[1] 1 3 4

and

1: a
2: d
3: c
4: b
5: e

==> 4 time(s):[1] 1 2 3 4 5
==> 1 time(s):[1] 1 4 5
==> 1 time(s):[1] 1 2 3
==> 3 time(s):[1] 1 2 5
==> 3 time(s):[1] 1 3 4

Each clade appeared the same number of times but the mapping between indices and leaf labels is different.

The output clade.freq.list only has the leaf indices

  cladenames cladefreqs
1  1 2 3 4 5       1.00
2      1 4 5       0.75
3      1 2 3       0.75
4      1 2 5       0.25
5      1 3 4       0.25

$`5e+00`
  cladenames cladefreqs
1  1 2 3 4 5       1.00
2      1 4 5       0.25
3      1 2 3       0.25
4      1 2 5       0.75
5      1 3 4       0.75

In producing the output for get.slide.freq.table, the clade.freq.list seems to get merged by the node indices even though they weren't the same. The tip names in the translation table also seem to be based on the first tree in the list rather than the first in each window.

$slide.table
  1e+00 5e+00        sd mean ess
1  0.75  0.25 0.3535534  0.5   0
3  0.25  0.75 0.3535534  0.5   0
4  0.25  0.75 0.3535534  0.5   0
5  0.75  0.25 0.3535534  0.5   0
2  1.00  1.00 0.0000000  1.0   0

$translation
     Clade number Tip numbers Tip names      
[1,] "1"          "1 2 3"     "a, b, c"      
[2,] "3"          "1 2 5"     "a, b, e"      
[3,] "4"          "1 3 4"     "a, c, d"      
[4,] "5"          "1 4 5"     "a, d, e"      
[5,] "2"          "1 2 3 4 5" "a, b, c, d, e"

If I instead have a chain

f3 <- list(trees = c(t1, t1, t1, t2, t1, t1, t2, t1), ptable =  NULL, gens.per.tree = 1)

then the outputs of get.slide.freq.table are identical for the two windows.

[danlwarren]

Hi there! I'm on holiday until the 25th but I'll have a good look at this when I get back if Rob doesn't get to it first. Thanks for bringing it to our attention, and even more for the pull request!

[lukejkelly]

Hi Dan, thanks a lot. I stepped through the code again this morning but will try some more examples when I get the chance. I have been working on diagnosing convergence in phylogenetic problems using pairs of chains run in parallel and trying to get them to meet, the preprint will hopefully appear next week.

[danlwarren]

That sounds awesome! I do have one request though: in the future could you submit pull requests to develop instead of master? I typically only update the master branch on github after a bit of testing, usually right before submitting a new version to CRAN.

Thanks again!