emmanuelparadis / ape

analysis of phylogenetics and evolution
http://ape-package.ird.fr/
GNU General Public License v2.0
52 stars 11 forks source link

corPhylo error - Error in Rd %*% R : non-conformable arguments #82

Closed MalharAtre closed 1 year ago

MalharAtre commented 1 year ago

Hi,

I am trying to get a phylogenetic correlation for two traits using corPhylo. The traits are being provided correctly as a single-column matrix with correct rownames for X and a list with one single column matrix with correct row names for U. Phylogeny provided is rooted, binary with branch lengths. There are no NaN values in either of trait values. The number of taxa in tree is ~2500. After running for some time, I end up with following error -

Error in Rd %*% R : non-conformable arguments

Google search suggests that this may be a matrix multiplication error associated with dimensions of matrix. However, I'm nor sure what may be triggering it.

Here is the screenshot of my command and error -

Screenshot 2023-05-03 at 12 06 23 PM

Looking forward to your reply! Thanks, Malhar

emmanuelparadis commented 1 year ago

Hi, By chance, is there an object named T in your workspace when you tried this command? You can try:

ls(pattern = "^T$")

If something is printed else than character(0), then:

get("T")
MalharAtre commented 1 year ago

Hi!

Thanks for the quick response!

There is no object named T in the workspace.

ls(pattern = "^T$")

gives-

character(0) .

There are objects named a and p, but even after removing them, the issue persists.

emmanuelparadis commented 1 year ago

Good! I noticed corphylo's code uses T in several places in place of TRUE (I'll fix that later).

It seems the error happens when executing these lines:

    ...
    L <- matrix(0, nrow = p, ncol = p)
    L[lower.tri(L, diag = T)] <- L.elements
    R <- t(L) %*% L
    Rd <- diag(diag(R)^-0.5)
    cor.matrix <- Rd %*% R %*% Rd
    ...

so all three matrices, L, R, Rd, are expected to be square. It's possible there is something in your data. You may explore this yourself after setting:

options(error = recover)

Or you can post some data here to reproduce the error so we can help.

MalharAtre commented 1 year ago

Hi! Thanks for the input.

I have tried exploring the issue and it appears that the L and R are 1x1 matrices, whereas the Rd produces a 3x3 matrix, owing to the small value of R.

Browse[3]> opt
$par
[1] 0.2602809 0.9503254

$value
[1] -2021.429

$counts
function gradient 
      51       NA 

$convergence
[1] 0

$message
NULL

Browse[3]> par
[1] 0.2602809 0.9503254
Browse[3]> LL
[1] -2021.429
Browse[3]> L.elements
[1] 0.2602809
Browse[3]> L
          [,1]
[1,] 0.2602809
Browse[3]> R
           [,1]
[1,] 0.06774617
Browse[3]> Rd
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
Browse[3]> diag(R)^-0.5
[1] 3.842002

Essentially, it appears that if the value of R drops below 0.25, ie. the optimized parameter stored in opt drops below 0.5, this error would occur at least in the case when only two traits are being compared.

Is the optimized value to not expected drop below 0.5? If so, I will go through my data accordingly to find explanations

Thank you for your help again! Looking forward to your reply,

Malhar

emmanuelparadis commented 1 year ago

It seems you didn't input the trait data correctly. From your first message above:

The traits are being provided correctly as a single-column matrix with correct rownames for X and a list with one single column matrix with correct row names for U.

So this results in p = 1 (single-column). See the help page for a description of what is the argument U. Did you input the second trait into U?