liamrevell / phytools

GNU General Public License v3.0
198 stars 56 forks source link

ratebytree not converging - what should I do? #125

Closed soungalo closed 1 year ago

soungalo commented 1 year ago

I ran ratebytree on two trees and corresponding character data, using model BM. The LRT P is significant (actually 0), but I got the message: One or the other optimization may not have converged.. What should I do? Is the test result still valid? Interestingly, this only happened after I log-transformed the data. I did this following suggestions from past publications and since the ranges of values in the two trees are quite different. Maybe this is not a good idea? Also, not sure if it's relevant, but when I used fitContinuous on each of the trees separately, I didn't get any warning messages regarding optimization convergence.

Thanks!

liamrevell commented 1 year ago

Hi @soungalo. The test result will only be valid if the null (common-rates) model converged. (By valid, I mean that you can be sure that the true result, even if different, would also be significant. To see which model has converged & which hasn't, you can run the function str on the fitted model object & look at the values of convergence for each fitted model. If it shows ..$ convergence: int 0 that means this model (is thought by R to have) converged. The easiest thing you can do to try to get your models to converge is increase the optional argument maxit. Its default value is maxit=500, so try maxit=2000 and see what happens. If that doesn't work, I suggest emailing me.

soungalo commented 1 year ago

Thanks for the quick reply. I tried with maxit=2000 and maxit=100000, but still getting the same message. Here is the result of str(rbt):

List of 11
 $ multi.rate.model :List of 13
  ..$ sig2       : num [1:2] 0.00343 0.000968
  ..$ SE.sig2    : num [1:2] 5.11e-04 8.78e-05
  ..$ a          : num [1:2] 5.79 8.08
  ..$ SE.a       : num [1:2] 0.626 0.289
  ..$ alpha      : NULL
  ..$ SE.alpha   : NULL
  ..$ r          : NULL
  ..$ SE.r       : NULL
  ..$ k          : num 4
  ..$ logL       : num -20.4
  ..$ counts     : Named int [1:2] 21 21
  .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
  ..$ convergence: int 52
  ..$ message    : chr "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"
 $ common.rate.model:List of 13
  ..$ sig2       : num 0.00195
  ..$ SE.sig2    : num 0.000184
  ..$ a          : num [1:2] 5.79 8.08
  ..$ SE.a       : num [1:2] 0.472 0.41
  ..$ alpha      : NULL
  ..$ SE.alpha   : NULL
  ..$ r          : NULL
  ..$ SE.r       : NULL
  ..$ k          : num 3
  ..$ logL       : num -53.3
  ..$ counts     : Named int [1:2] 27 27
  .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
  ..$ convergence: int 0
  ..$ message    : chr "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
 $ model            : chr "BM"
 $ N                : int 2
 $ n                : int [1:2] 92 245
 $ likelihood.ratio : num 65.8
 $ regimes          : Factor w/ 2 levels "1","2": 1 2
 $ m                : int 2
 $ P.chisq          : num 4.95e-16
 $ P.sim            : NULL
 $ type             : chr "continuous"
 - attr(*, "class")= chr "ratebytree"
liamrevell commented 1 year ago

I'm guessing that this is probably fine because the common-rate model converged. Since the fitted multi-rate model should have the same values of sigma as a simple Brownian model fit to each of the characters separately, and the total log-likelihood just the sum of the individual log-likelihoods of each of the two fitted models, you could try just fitting each of those two models separately (using, say, geiger::fitContinuous) and then comparing the estimated parameters. If you can't figure out how to get all that to work, please email me your saved R workspace & a minimal code to reproduce this result and I'd be happy to help.

soungalo commented 1 year ago

Thanks, got it.