matthewwolak / nadiv

R package that constructs (non)additive genetic relationship matrices, and their inverses, from a pedigree to be used in linear mixed effect models (A.K.A. the 'animal model').
17 stars 7 forks source link

Wrong makeA results with Matrix v1.6-1 #36

Closed Ax3man closed 1 year ago

Ax3man commented 1 year ago

Hi Matthew,

I am a very happy user of nadiv, but running into some mysterious issues. After updating some packages, my workstation install of nadiv creates obviously incorrect A matrices. E.g.:

> makeA(Mrode2)
6 x 6 sparse Matrix of class "dsCMatrix"
          1         2         3         4         5         6
1 1.8125000 0.5312500 0.5745243 0.7036456 0.4419417 0.1711633
2 0.5312500 1.7031250 0.5524272 0.2435696 0.3977476 0.4279082
3 0.5745243 0.5524272 0.6562500 0.1913664 0.3125000 0.1210307
4 0.7036456 0.2435696 0.1913664 0.9843750 0.3827328 0.1482318
5 0.4419417 0.3977476 0.3125000 0.3827328 0.6250000 0.2420615
6 0.1711633 0.4279082 0.1210307 0.1482318 0.2420615 0.4687500

(In my own data, full siblings have relatedness coefficients of 0.)

Compare with a version running on our cluster:

>  makeA(Mrode2)
6 x 6 sparse Matrix of class "dsCMatrix"
     1     2      3      4      5      6
1 1.00 .     0.5000 0.5000 0.5000 0.2500
2 .    1.000 0.5000 .      0.2500 0.6250
3 0.50 0.500 1.0000 0.2500 0.6250 0.5625
4 0.50 .     0.2500 1.0000 0.6250 0.3125
5 0.50 0.250 0.6250 0.6250 1.1250 0.6875
6 0.25 0.625 0.5625 0.3125 0.6875 1.1250

Both machines are running nadiv_2.17.2. However my workstation is using Matrix_1.6-1, while the cluster uses Matrix_1.5-3. Could that be the issue?

Any help much appreciated!

Ax3man commented 1 year ago

I rolled back Matrix using remotes::install_version('Matrix', '1.5-3'), and this indeed fixes the issue.

matthewwolak commented 1 year ago

Hi @Ax3man,

Thanks for the comment and sorry about the errors in nadiv calculations. The Matrix package maintainers found an error in their chol2inv() function that affects nadiv::makeA(). We've sorted out how the Matrix fix will affect makeA() but I just now need to update a new version of nadiv that works with the most recent (bug fixed) version of Matrix.

I will try and get a new version of nadiv to CRAN in the next 1-3 weeks, but am traveling and at a conference all of next week so that might slow things down.

In the meantime, this has been fixed in the devel branch of nadiv (see commit b822dfe). You can get this if you have the remotes package using the line to install from the devel branch:

remotes::install_github("matthewwolak/nadiv", ref = "devel")

Matthew

Ax3man commented 1 year ago

Thanks!

To clarify, does the bug in Matrix <1.6 mean that previous results of makeA() were also incorrect?

matthewwolak commented 1 year ago

Hi,

Short answer: Previous results of makeA() are correct even with the bug in Matrix<1.6. These combinations should give correct A matrices:

Long answer: The bug in Matrix<1.6 is due to Matrix::chol2inv() mixing up lower- and upper-triangle matrices along with their transposes (understandably, I always do that too). When using Matrix<1.6 to develop the last version of makeA(), I was confused by the incorrect A matrices produced by Matrix::chol2inv() and so played around (guess-and-check) with different combinations of upper- versus lower-triangle matrices and transposes to give Matrix::chol2inv() until it gave correct A matrices. What that means is now that Matrix >= 1.6 corrected these upper/lower and transpose issues the nadiv=2.17.2 makeA() currently on CRAN gives the wrong matrices to chol2inv(), but that is what was corrected in the devel branch of nadiv (i.e., 2.17.3 ).

Ax3man commented 1 year ago

Perfect, that is very clear. Thanks so much!