Error message from brownie.lite and brownieREML

A user recently reported receiving the following error messages from brownieREML:

> res<-brownieREML(mtree,x)
Error in optim(rep(1, p) * runif(n = p), likelihood2, tree = tree, x = x, :
L-BFGS-B needs finite values of 'fn'

and from brownie.lite:

> res<-brownie.lite(mtree,x)
Error in solve.default(sig * C) :
Lapack routine dgesv: system is exactly singular

respectively.

The above error message is typical of trees with terminal branch lengths of zero length. This is because if two terminal taxa from an internal node have zero length branches, then the matrix vcv.phylo(tree) is singular (i.e., it has a determinant of zero and cannot be inverted). However, in the case of this user's dataset and tree, this was not the problem.

What I discovered was that the problem disappeared if I increased the scale of the data, say, by a factor of 100 or 1000. So, for instance:

> res<-brownie.lite(mtree,100*x)

worked fine.

This led me to suspect that the problem resulted from very low rates along some branches of the tree which resulted in a among species covariance matrix that was singular to numerical precision. In fact, this seems to have been the case as the fitted evolutionary rates are indeed extremely small along some of the branches of the phylogeny.

Multiplying the data by a scalar preserves the relative rates for different regimes as well as the difference in log-likelihood between different models for the scaled data; and the rates themselves are scaled by a factor of the constant squared. So, for instance, consider the following illustrative example:

> # load packages
> require(phytools); require(geiger)
> # simulate tree & data
> tree<-rescaleTree(drop.tip(birthdeath.tree(b=1,d=0, taxa.stop=101),"101"),1)
> Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
> mtree<-sim.history(tree,Q)
> sig2<-c(1,10,100); names(sig2)<-c(1,2,3)
> x<-sim.rates(mtree,sig2)
> # fit model to data on the original scale
> res1<-brownie.lite(mtree,x)
> # fit to data scaled by a factor of 10
> res10<-brownie.lite(mtree,x*10)
> # compare relative rates
> res1$sig2.multiple/res1$sig2.multiple["1"]
        2         3         1 
 9.575183 92.854992  1.000000 
> res10$sig2.multiple/res10$sig2.multiple["1"]
        2         3         1 
 9.574489 92.844736  1.000000 
> # compare absolute rates
> res10$sig2.multiple/res1$sig2.multiple
        2         3         1 
 99.99858  99.99477 100.00582
> # compare log-likelihood differences
> res1$logL.multiple-res1$logL1
[1] 24.51724
> res10$logL.multiple-res10$logL1
[1] 24.51724