brj1 / node.dating Goto Github PK

Hello!
I've used node.dating for my research and encountered weird behavior - some nodes received unrealistically early date. So I've taken a closer look inside the estimate.dates function and I am convinced that I've found a bug. Which I'd like to report. Here I go.
The function calc.Like takes as an input ch.node, ch.edge and x. Inside this function the vector t$edge.length[ch.edge] is calcultaed based on ch.edge value - the length of the edges leading to children nodes, based on the children nodes IDs. But inside the opt.fun function the function calc.Like is called with arguments dates[ch], t$edge.length[ch.edge], x. So, instead of passing ch.edge to calc.Like function, leaving it to calculate t$edge.length[ch.edge] itself, it passes t$edge.length[ch.edge]. There is another call of calc.Like inside opt.fun, and arguments there also need to be fixed in the same manner.
I guess with edge length being expressed as floating point numbers such operation should result in straight up error, but in my case edge lengths were expressed in integers. Because of that opt.fun calculated the lengths of edges, leading to children, used these as IDs for children nodes and inside calc.Like function the length of edges corresponding to these IDs were taken. Needless to say, it's not what I wanted)
Dmitry

@emvolz-phylodynamics is getting discrepant results from running test script tester.get.rmse.R on the simulated data in the repo. Specifically, he is getting a mean RMSE of about 420, whereas we reported a value of about 29 in the manuscript. I attempted to reproduce the issue on my Mac and got an RMSE of 29.21 as expected. However, when I ran the same script and data on an Ubuntu workstation, I got an RMSE of about 426.

Minimal case:

> require(ape)
> vect <- c(1997.499, 2007.262, 1999.498, 2010.498)
> tr <- read.tree(text="((A:0.084018109,B:0.074898812):0.1541108577,(D:0.039097599,E:0.0811621)1.000:0.1864050113)NA;")
> estimate.mu(tr, vect)
Error in pchisq(AIC(null.g) - AIC(g)) : 
  argument "df" is missing, with no default

The problem seems to be on line 8 of estimate.mu:

        warning(paste("Cannot reject null hypothesis (p=", (1 - 
            pchisq(AIC(null.g) - AIC(g))), ")"))

which should be replaced with:

        warning(paste("Cannot reject null hypothesis (p=", (1 - 
            pchisq(AIC(null.g) - AIC(g), df=1)), ")"))

We've determined that the problem arises when calling solve.cube with the following arguments:

> c(m, min(dates[ch]))
[1] 2007.183 2007.183
> dates[ch]
[1] 2007.202 2007.183
> ch.edge
[1] 5864 5865
> dates[p]
[1] 2007.183
> p.edge
[1] 5863
> solve.cube(c(m, min(dates[ch])), dates[ch], 
+                      ch.edge, dates[p], p.edge)
Warning in log(tim) : NaNs produced
Warning in log(tim) : NaNs produced
numeric(0)