Quadform p-value=2

[yyang-peko]

Sometimes when running quadform.qf(), it outputs 2.0, 4, array([ 0., 0., 0., 0., 0., 0., 1000001.]) when I input, for example, chi2val=359.4848154869432 [36.94163236 24.76110077 24.74140768 18.61546473 18.61546473 18.61546473 18.61546473 18.61546473 18.58195836 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4125513 ]

Should it just be interpreted as p-value essentially = 1?

[CarlKCarlK]

YYang,

Thanks for your question and thanks for using FastLmmCLib! I should be able to look into this tomorrow.

• Carl

Carl Kadie, Ph.D. FaST-LMM & PySnpTools Teamhttps://fastlmm.github.io/ (Microsoft Research, retired) https://www.linkedin.com/in/carlk/

Join the FaST-LMM user discussion and announcement list via @.***?subject=Subscribe> (or use web sign uphttps://mail.python.org/mailman3/lists/fastlmm-user.python.org)

From: YYang @.> Sent: Thursday, September 14, 2023 1:41 PM To: fastlmm/fastlmmclib @.> Cc: Subscribed @.***> Subject: [fastlmm/fastlmmclib] Quadform p-value=2 (Issue #1) Importance: High

Sometimes when running quadform.qf(), it outputs 2.0, 4, array([ 0., 0., 0., 0., 0., 0., 1000001.]) when I input, for example, chi2val=359.4848154869432 [36.94163236 24.76110077 24.74140768 18.61546473 18.61546473 18.61546473 18.61546473 18.61546473 18.58195836 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4401378 12.4125513 ]

Should it just be interpreted as p-value essentially = 1?

- Reply to this email directly, view it on GitHubhttps://github.com/fastlmm/fastlmmclib/issues/1, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABR65P47KPTEKIUHGQVYWJLX2NTWNANCNFSM6AAAAAA4YZXIO4. You are receiving this because you are subscribed to this thread.Message ID: @.**@.>>

[CarlKCarlK]

I'm trying to reproduce the issue, but haven't been able to. Here is what I tried:

def test_sep23_qf():
    eigvals = np.array(
        [
            36.94163236,
            24.76110077,
            24.74140768,
            18.61546473,
            18.61546473,
            18.61546473,
            18.61546473,
            18.61546473,
            18.58195836,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4401378,
            12.4125513,
        ]
    )
    result = qf(359.4848154869432, eigvals)
    print(result)

And I saw this:

(0.36692794169308307, 0, array([2.70007200e+00, 4.10000000e+01, 1.00000000e+00, 3.79025439e-03,
       1.53418915e-01, 0.00000000e+00, 2.20000000e+01]))

Can you help me reproduce the issue? Thanks, Carl

[yyang-peko]

Manually running qf() on chi2vals and eigvals seem to have fixed this. I'm not sure what the cause of it is. I was calling qf() from python script running in a slurm array on an HPC, and around 6 out of 10 times my script would stop at the qf() step and cannot proceed further. Within the 40% of the times it succeeded in running it would produce the 2.0, 4, array([ 0., 0., 0., 0., 0., 0., 1000001.]) result.

Just collecting chi2vals and eigvals, then running the qf separately fixed the issue.

[CarlKCarlK]

Excellent!

From: YYang @.> Sent: Saturday, September 16, 2023 11:44 AM To: fastlmm/fastlmmclib @.> Cc: Carl Kadie @.>; Comment @.> Subject: Re: [fastlmm/fastlmmclib] Quadform p-value=2 (Issue #1) Importance: High

Manually running qf() on chi2vals and eigvals seem to have fixed this. I'm not sure what the cause of it is. I was calling qf() from python script running in a slurm array on an HPC, and around 6 out of 10 times my script would stop at the qf() step and cannot proceed further. Within the 40% of the times it succeeded in running it would produce the 2.0, 4, array([ 0., 0., 0., 0., 0., 0., 1000001.]) result.

Just collecting chi2vals and eigvals, then running the qf separately fixed the issue.

- Reply to this email directly, view it on GitHubhttps://github.com/fastlmm/fastlmmclib/issues/1#issuecomment-1722291302, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABR65PY5JE5PQDG2IGYPGB3X2XXOJANCNFSM6AAAAAA4YZXIO4. You are receiving this because you commented.Message ID: @.**@.>>

[yyang-peko]

Hi Carl,

I'm encountering a very strange issue with quadform.

I'm doing several matrix operations etc to end up with a chi2val and a numpy array of eigvals. If I save the eigvals (np.save) and load it, checking that all elements in the array is equal to the generated array, the resulting p-values and trace from qf.qf() is different when using the generated eigvals vs the save-and-loaded eigvals.

Specifically in a test case:

chi2val = 2998.1754859648086 eigvals_ = [601.283561554898 , 284.7519233469181 , 132.09058972616828 , 91.17672134789915 , 49.23921245678543 , 37.04139197122487 , 30.87515072058668 , 30.83070804169951 , 27.387885652473727 , 24.78420217167126 , 24.735588068850603 , 18.639481731694534 , 18.606664684701848 , 18.598222876777783 , 18.5877915113293 , 18.532732895132938 , 18.434416620497895 , 12.45605352114341 , 12.454616462471623 , 12.453538158305353 , 12.452084879831665 , 12.4512097640746 , 12.449710630488513 , 12.443930072105454 , 12.440381732335313 , 12.432859778920395 , 12.419744093636476 , 12.40376996354169 , 12.389340897058112 , 12.366556596613494 , 9.867422198584155]

qf.qf(chi2val, eigvals) = (0.025548898985580548, 0, array([1.19763908e+02, 7.61530000e+04, 2.00000000e+00, 2.57624534e-04, 2.08522844e+01, 2.43653829e-01, 9.50000000e+01])) the original eigvals generated via sp.sparse.linalg.eigsh() qf.qf(chi2val, eigvals_) = (0.08532721621254413, 0, array([3.36862960e+00, 2.42000000e+02, 1.00000000e+00, 2.46419103e-04, 5.93668452e-02, 0.00000000e+00, 2.20000000e+01])) saved and loaded

eigvals == eigvals_ gives 'True' for all elements.

Do you have any ideas for what could be the issue?

Thanks, Yang

On Sat, Sep 16, 2023 at 11:56 AM Carl Kadie @.***> wrote:

Excellent!

From: YYang @.> Sent: Saturday, September 16, 2023 11:44 AM To: fastlmm/fastlmmclib @.> Cc: Carl Kadie @.>; Comment @.> Subject: Re: [fastlmm/fastlmmclib] Quadform p-value=2 (Issue #1) Importance: High

Manually running qf() on chi2vals and eigvals seem to have fixed this. I'm not sure what the cause of it is. I was calling qf() from python script running in a slurm array on an HPC, and around 6 out of 10 times my script would stop at the qf() step and cannot proceed further. Within the 40% of the times it succeeded in running it would produce the 2.0, 4, array([ 0., 0., 0., 0., 0., 0., 1000001.]) result.

Just collecting chi2vals and eigvals, then running the qf separately fixed the issue.

- Reply to this email directly, view it on GitHub< https://github.com/fastlmm/fastlmmclib/issues/1#issuecomment-1722291302>, or unsubscribe< https://github.com/notifications/unsubscribe-auth/ABR65PY5JE5PQDG2IGYPGB3X2XXOJANCNFSM6AAAAAA4YZXIO4

. You are receiving this because you commented.Message ID: @.**@.>>

— Reply to this email directly, view it on GitHub https://github.com/fastlmm/fastlmmclib/issues/1#issuecomment-1722293624, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABF7K4TEHTVK6PAV7ULAD3TX2XY4FANCNFSM6AAAAAA4YZXIO4 . You are receiving this because you modified the open/close state.Message ID: @.***>

[CarlKCarlK]

Yang,

I’m happy to try to help. I’m trying to create an IPython notebook to show the problem. However, I’m having trouble creating the test case.

In your example, I see eigvals, but not eigvals. Also eigvals is a 1-D list, so it can’t be given give it to qf.qf.

Yours,

• Carl

From: YYang @.> Sent: Thursday, December 28, 2023 1:08 PM To: fastlmm/fastlmmclib @.> Cc: Carl Kadie @.>; Comment @.> Subject: Re: [fastlmm/fastlmmclib] Quadform p-value=2 (Issue #1) Importance: High

Hi Carl,

I'm encountering a very strange issue with quadform.

I'm doing several matrix operations etc to end up with a chi2val and a numpy array of eigvals. If I save the eigvals (np.save) and load it, checking that all elements in the array is equal to the generated array, the resulting p-values and trace from qf.qf() is different when using the generated eigvals vs the save-and-loaded eigvals.

Specifically in a test case:

chi2val = 2998.1754859648086 eigvals_ = [601.283561554898 , 284.7519233469181 , 132.09058972616828 , 91.17672134789915 , 49.23921245678543 , 37.04139197122487 , 30.87515072058668 , 30.83070804169951 , 27.387885652473727 , 24.78420217167126 , 24.735588068850603 , 18.639481731694534 , 18.606664684701848 , 18.598222876777783 , 18.5877915113293 , 18.532732895132938 , 18.434416620497895 , 12.45605352114341 , 12.454616462471623 , 12.453538158305353 , 12.452084879831665 , 12.4512097640746 , 12.449710630488513 , 12.443930072105454 , 12.440381732335313 , 12.432859778920395 , 12.419744093636476 , 12.40376996354169 , 12.389340897058112 , 12.366556596613494 , 9.867422198584155]

qf.qf(chi2val, eigvals) = (0.025548898985580548, 0, array([1.19763908e+02, 7.61530000e+04, 2.00000000e+00, 2.57624534e-04, 2.08522844e+01, 2.43653829e-01, 9.50000000e+01])) the original eigvals generated via sp.sparse.linalg.eigsh() qf.qf(chi2val, eigvals_) = (0.08532721621254413, 0, array([3.36862960e+00, 2.42000000e+02, 1.00000000e+00, 2.46419103e-04, 5.93668452e-02, 0.00000000e+00, 2.20000000e+01])) saved and loaded

eigvals == eigvals_ gives 'True' for all elements.

Do you have any ideas for what could be the issue?

Thanks, Yang

On Sat, Sep 16, 2023 at 11:56 AM Carl Kadie @.<mailto:@.>> wrote:

Excellent!

From: YYang @.<mailto:@.>> Sent: Saturday, September 16, 2023 11:44 AM To: fastlmm/fastlmmclib @.<mailto:@.>> Cc: Carl Kadie @.<mailto:@.>>; Comment @.<mailto:@.>> Subject: Re: [fastlmm/fastlmmclib] Quadform p-value=2 (Issue #1) Importance: High

Manually running qf() on chi2vals and eigvals seem to have fixed this. I'm not sure what the cause of it is. I was calling qf() from python script running in a slurm array on an HPC, and around 6 out of 10 times my script would stop at the qf() step and cannot proceed further. Within the 40% of the times it succeeded in running it would produce the 2.0, 4, array([ 0., 0., 0., 0., 0., 0., 1000001.]) result.

Just collecting chi2vals and eigvals, then running the qf separately fixed the issue.

- Reply to this email directly, view it on GitHub< https://github.com/fastlmm/fastlmmclib/issues/1#issuecomment-1722291302>, or unsubscribe< https://github.com/notifications/unsubscribe-auth/ABR65PY5JE5PQDG2IGYPGB3X2XXOJANCNFSM6AAAAAA4YZXIO4

. You are receiving this because you commented.Message ID: @.**@.mailto:***@***.******@***.***>>

— Reply to this email directly, view it on GitHub https://github.com/fastlmm/fastlmmclib/issues/1#issuecomment-1722293624, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABF7K4TEHTVK6PAV7ULAD3TX2XY4FANCNFSM6AAAAAA4YZXIO4 . You are receiving this because you modified the open/close state.Message ID: @.<mailto:@.>>

— Reply to this email directly, view it on GitHubhttps://github.com/fastlmm/fastlmmclib/issues/1#issuecomment-1871490179, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABR65P3PSERPKLFXCZXW77LYLXNTTAVCNFSM6AAAAAA4YZXIO6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNZRGQ4TAMJXHE. You are receiving this because you commented.Message ID: @.**@.>>

[yyang-peko]

SKS_matrix.zip

Hi Carl, I've uploaded a test data file. It's a 286 x 286 hermitian matrix stored as a numpy pickle.

SKS = np.load('SKS_matrix.npy', allow_pickle=True) eigvals = sp.sparse.linalg.eigsh(SKS, k=40, return_eigenvectors=False) eigvals = np.sort(eigvals[np.where(eigvals > 1e-8)])[::-1] Q = 585.0978079020846 qf.qf(Q, eigvals, acc=1e-7)

(1.5605687704312032e-06, 0, array([1.38445730e+02, 2.88600000e+03, 1.00000000e+00, 1.07098983e-02, 3.09010677e+01, 1.56866434e-01, 7.00000000e+01]))

However, when I do qf.qf(Q, eigvals.astype(np.float64), acc=1e-7) qf.qf(Q, eigvals.astype(float), acc=1e-7)

(0.07390013408306717, 0, array([3.73084527e+00, 1.60000000e+01, 1.00000000e+00, 6.28393232e-03, 9.16086013e-02, 0.00000000e+00, 1.80000000e+01]))

I suspect this has something to do with how sp.sparse.linalg.eigsh returns eigenvalues at machine precision, and some interaction between numpy's float64, python's float64, and C's float64_t.

for i in eigvals: print(i)

.359430815339177 24.26909591539037 22.16434645647498 20.224204984397332 20.172814719026565 19.399440077993294 17.795433563879307 17.508064982626617 16.361888020192687 14.465250735324332 14.41711515040611 14.386829364762525 14.215687314249916 13.495097410402899 13.08161070596028 10.55561990110087 10.552489150903842 10.548784209417585 10.5454783523345 10.542843626630576 10.536236974566016 10.515292888092002 10.504225556873383 10.450495164026837 10.402338819902715 10.35366667058023 9.781555167439034 8.238279080876648 8.043976264341268 2.0183786464449573 2.142120841448611e-05

for i in eigvals: print("{:.60f}".format(i))

25.359430815339177200939957401715219020843505859375000000000000 24.269095915390369810893389512784779071807861328125000000000000 22.164346456474980584516742965206503868103027343750000000000000 20.224204984397331941181619185954332351684570312500000000000000 20.172814719026565200010736589320003986358642578125000000000000 19.399440077993293840563637786544859409332275390625000000000000 17.795433563879306859689677366986870765686035156250000000000000 17.508064982626617478445041342638432979583740234375000000000000 16.361888020192687065446079941466450691223144531250000000000000 14.465250735324332254094770178198814392089843750000000000000000 14.417115150406109691516576276626437902450561523437500000000000 14.386829364762524718912573007401078939437866210937500000000000 14.215687314249915829122983268462121486663818359375000000000000 13.495097410402898674419702729210257530212402343750000000000000 13.081610705960279972259741043671965599060058593750000000000000 10.555619901100870450250113208312541246414184570312500000000000 10.552489150903841874651334364898502826690673828125000000000000 10.548784209417584989409988338593393564224243164062500000000000 10.545478352334500371512149285990744829177856445312500000000000 10.542843626630576281399953586515039205551147460937500000000000 10.536236974566016044718708144500851631164550781250000000000000 10.515292888092002243638489744625985622406005859375000000000000 10.504225556873382885214596171863377094268798828125000000000000 10.450495164026836647508389432914555072784423828125000000000000 10.402338819902714917020603024866431951522827148437500000000000 10.353666670580230046994074655231088399887084960937500000000000 9.781555167439034192966573755256831645965576171875000000000000 8.238279080876647952891289605759084224700927734375000000000000 8.043976264341267778945621103048324584960937500000000000000000 2.018378646444957347227955324342474341392517089843750000000000 0.000021421208414486109049669829151874012040934758260846138000

But, even if I create a new np.array with the above values, or save-and-load eigvals to disk, I still get the same 2nd qf() results.

[CarlKCarlK]

Do we expect SKS to contain imginary numbers?

import numpy as np
SKS = np.load('SKS_matrix.npy', allow_pickle=True)
SKS

array([[ 0.04429796-1.27082473e-09j, -0.04714658-1.32167941e-10j,
         0.01433359-3.37525097e-10j, ...,  0.02179421-4.00270757e-09j,
        -0.02024132-1.62844172e-09j, -0.03052728-7.82222069e-10j],
       [-0.04714658-3.93228427e-11j,  5.00593631+7.15638676e-10j,
      #...

[CarlKCarlK]

Also, what is phis1?

[yyang-peko]

Also, what is phis1?

Edited. Sorry.

SKS should be complex since S is the sq rt of another matrix

[CarlKCarlK]

OK. I can reproduce the problem and think I have a workaround.

I think the issue is that qf.qf expects a contiguous array, but

eigvals = np.sort(eigvals[np.where(eigvals > 1e-8)])[::-1] is noncontiguous.

The fix is to call np.ascontiguousarray (or .copy(), or as_type, etc)

qf.qf(Q, np.ascontiguousarray(eigvals), acc=1e-7), eigvals, eigvals.dtype

It then always produces the 0.07390013408306784, ... answer.

If this seems right, let me know and I'll eventually update the qf.qf code to always call np.ascontiguousarray.

• Carl p.s. The way I could see that the original eigvals was not contiguous was to print eigvals.flags. It said:

C_CONTIGUOUS : False
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False

[yyang-peko]

I see. So any modification to the generated eigvals turns it back to a contiguous array. Ah the sorting turns it non-contiguous! I thought it was the eigenvalue computation itself.

Thank you so much, Carl. Have a great New Year's.

[yyang-peko]

Hi Carl,

I discovered what was causing the original issue of this thread. When chi2val is too small or too large, the p-value can become 2.0 or a negative tiny float before rounding to 1.0 or 0.0, respectively. And where these happens depend on the acc option. Probably float under/overflows?

qf.qf(0.068, np.array([12.4]))
>>>(0.9409680834684881, 0, array([4.51497377e+00, 9.88515000e+05, 7.00000000e+00, 9.30972623e-03, 7.39576675e+05, 5.52618098e-06, 1.83000000e+02]))

qf.qf(0.067, np.array([12.4]))
>>>(2.0, 1, array([0.00000000e+00, 7.10867000e+05, 6.00000000e+00, 9.30971243e-03, 8.06810919e+05, 5.44491361e-06, 1.83000000e+02]))

qf.qf(25428, np.array([620.5, 312.0, 204.8, 137.4]))
>>>(2.4491519923230953e-13, 0, array([2.23133202e+01, 5.68300000e+03, 1.00000000e+00, 2.47007504e-04, 1.40350064e+00, 2.14053048e+00, 3.60000000e+01]))

qf.qf(25429, np.array([620.5, 312.0, 204.8, 137.4]))
>>>(-1.1501910535116622e-13, 0, array([2.23141681e+01, 5.68300000e+03, 1.00000000e+00, 2.46997786e-04, 1.40350064e+00, 2.14061910e+00, 3.60000000e+01]))

qf.qf(22961, np.array([620.5, 312.0, 204.8, 137.4]), acc=1e-7)
>>>(3.0772051573535464e-12, 0, array([2.00929994e+01, 1.77000000e+03, 1.00000000e+00, 2.73181747e-04, 4.83299118e-01, 6.07757605e+00, 3.40000000e+01]))

qf.qf(22962, np.array([620.5, 312.0, 204.8, 137.4]), acc=1e-7)
>>>(-4.4364512064021255e-13, 0, array([2.00938365e+01, 1.77000000e+03, 1.00000000e+00, 2.73169826e-04, 4.83299118e-01, 6.07785630e+00, 3.40000000e+01]))

Thanks.