Tuesday, November 13, 2012

Bug fix for phyl.RMA

I recently received the following bug report for the phytools function for reduced major axis regression (phyl.RMA):

I've been encountering an error when trying to run phyl.RMA:
Error in D %*% t(a) : non-conformable arguments
In addition: Warning message:
In if (lambda == 1) return(C) else { :
the condition has length > 1 and only the first element will be used
What's odd is that with R 2.15.1 and the previous version of phytools, it ran perfectly. But today I updated to R 2.15.2 and the latest version of phytools, and it's a no-go. Any idea what's going on here?


Indeed - I can reproduce this error easily:

> require(phytools)
Loading required package: phytools
...
> tree<-pbtree(n=100)
> x<-fastBM(tree)
> y<-fastBM(tree)
> args(phyl.RMA)
function (x, y, tree, method = "BM", lambda = NULL, fixed = FALSE,
    h0 = 1)
NULL
> phyl.RMA(x,y,tree)
Error in D %*% t(a) : non-conformable arguments
In addition: Warning message:
In if (lambda == 1) return(C) else { :
  the condition has length > 1 and only the first element will be used

Turns out, though, that the association with R 2.15.2 is completely spurious. The real problem came when I moved the code for various utility functions out of the functions that used them, such as phyl.pca, phyl.cca, and, lo & behold, phyl.RMA. As it turned out, some of the argument list order differed between different iterations of internally called functions (the culprit in this case was lambda.transform).

I have fixed the problem, and the revised function (here) seems to work:

> source("phyl.RMA.R")
> phyl.RMA(x,y,tree)
$RMA.beta
[1] 0.06725208 0.86425086
$V
           x          y
x 0.93872556 0.06935903
y 0.06935903 0.70116186
$lambda
[1] 1
$logL
[1] -264.4581
$test
        r2           T          df           P
0.00730885  1.44956623 99.64317035  0.15031966
$resid
              [,1]
t39  -0.2958474582
t40   0.5487124406
t9    ...

Good.

No comments:

Post a Comment