Error in `solve.default(C)` with `phytools::anc.ML`

[anoban]

Hi Liam,

I'm trying to fit an OU model to the continuous trait data (2 traits) that I have. The issue is I keep running into erros that originate from solve.default(C).

In the first attempt, I ran into the following error:

Error in solve.default(C): Lapack routine dgesv: system is exactly singular: U[1013,1013] = 0
Traceback:

1. anc.OU(tree, x, maxit, ...)
2. phyl.vcv(as.matrix(c(x, setNames(pp[2:tree$Nnode + 2], 2:tree$Nnode + 
 .     Ntip(tree)))), vcvPhylo(tree), lambda = 1)
3. solve(C)
4. solve.default(C)
5. .handleSimpleError(function (cnd) 
 . {
 .     watcher$capture_plot_and_output()
 .     cnd <- sanitize_call(cnd)
 .     watcher$push(cnd)
 .     switch(on_error, continue = invokeRestart("eval_continue"), 
 .         stop = invokeRestart("eval_stop"), error = NULL)
 . }, "Lapack routine dgesv: system is exactly singular: U[1013,1013] = 0", 
 .     base::quote(solve.default(C)))

I read that this happens when the phylogeny has 0 length branches. So I then replaced them with a small constant.

edge_length_adjusted_phylogeny <- phylogeny # create a copy
edge_length_adjusted_phylogeny$edge.length[edge_length_adjusted_phylogeny$edge.length <= 0] <- 1e-10 # replace zeroes with a single small positive value
phytools::anc.ML(tree = edge_length_adjusted_phylogeny, model = "OU", x = RD, maxit = 5000) # try again

But then ran into a new error:

Error in solve.default(C): system is computationally singular: reciprocal condition number = 3.13435e-17
Traceback:

1. anc.OU(tree, x, maxit, ...)
2. phyl.vcv(as.matrix(c(x, setNames(pp[2:tree$Nnode + 2], 2:tree$Nnode + 
 .     Ntip(tree)))), vcvPhylo(tree), lambda = 1)
3. solve(C)
4. solve.default(C)
5. .handleSimpleError(function (cnd) 
 . {
 .     watcher$capture_plot_and_output()
 .     cnd <- sanitize_call(cnd)
 .     watcher$push(cnd)
 .     switch(on_error, continue = invokeRestart("eval_continue"), 
 .         stop = invokeRestart("eval_stop"), error = NULL)
 . }, "system is computationally singular: reciprocal condition number = 3.13435e-17", 
 .     base::quote(solve.default(C)))

Then thinking that using a single value to replace all the 0 length branches (my phylogeny had 91 of them), I then replaced those branch lengths with small but unique random numbers.

replacement_lengths <- rnorm(n=sum(phylogeny$edge.length <= 0), mean = 1.0, sd = 0.05) * 1e-10
stopifnot(all(replacement_lengths > 0))

edge_length_adjusted_phylogeny <- phylogeny
edge_length_adjusted_phylogeny$edge.length[edge_length_adjusted_phylogeny$edge.length <= 0] <- replacement_lengths
phytools::anc.ML(tree = edge_length_adjusted_phylogeny, model = "OU", x = RD, maxit = 5000)

But ran into the system is computationally singular: reciprocal condition number = 3.0889e-17 error again. How could I resolve this? I know for a fact that my phylogeny is binary, rooted and ultrametric. phytools is the only option I have as other packages like geiger::fitContinuous does not provide the actual ACEs for internal nodes whereas phytools::fastAnc did and I hope phytools::anc.ML does too. I cannot use fastAnc because in model comparison using geiger::fitContinuous, the OU model turned out to be a better fit than the BM model. I'm including the links to the phylogeny and the trait data here.

PS. the RD variable is a named vector created as RD <- setNames(data$F00679, nm = data$binominal)
I'd really appreciate any help regarding this. Thank you!