question on susie_get_cs

[harryyiheyang]

Dear Dr. Wang and others,
I noticed a potential bug that the output of susie_get_cs() and the main susie functions (e.g., susie_rss()) may result in different credible sets.
An example is shown below:

#####################################################
summary(fit)

Credible sets summary:

cs cs_log10bf cs_avg_r2 cs_min_r2
1 206.28338 1.0000000 1.0000000
2 11.04498 0.9044700 0.9044700
3 2.51137 0.8569526 0.5166298
#####################################################

As you can see, the output of the main susie function contained three credible sets.

However, when I used susie_get_cs() to get the credible sets, no matter how large the min_abs_corr parameter I set, it generated a suspicious 4th credible set that contains more than 3000 SNPs:
#####################################################
susie_get_cs(fit,min_abs_corr=0.999)
$cs
$cs$L1
[1] 3699

$cs$L2
[1] 2852 2857

$cs$L3
[1] 6489 6490 6491 6492 6493 6494 6495 6496 6498 6499 6500 6501 6504 6505 6506 6507 6508 6509
[19] 6510 6511 6512 6514 6516 6518 6519 6523 6525 6526 6528 6529 6530 6531 6532 6533 6537 6543
[37] 6546 6547 6548 6550 6552 6555 6556 6557 6562 6564 6566 6568 6570 6576 6577 6585 6589 6590
[55] 6591 6593 6595 6596 6598 6599 6600 6601 6602 6604 6610 6613 6620

$cs$L4
[1] 1 2 3 4 5 6 7 8 10 11 12 13 14 15 18 19 20
[18] 21 22 23 24 25 27 28 29 30 31 32 33 34 35 36 38 39
[35] 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
[52] 57 58 59 60 61 62 63 64 65 66 67 68 70 71 72 73 74
[69] 75 76 77 78 79 81 82 83 84 86 87 88 89 90 91 92 93
[86] 94 97 98 99 102 105 106 109 110 111 112 114 115 117 118 119 120
[103] 121 122 123 124 126 127 129 131 132 134 135 136 137 138 139 140 143
........................................
[ reached getOption("max.print") -- omitted 3811 entries ]

$coverage
[1] 1.0000000 0.9585776 0.9528785 0.9500147

$requested_coverage
[1] 0.95
#####################################################
I also checked the correlation between X1 and X2, which was far less than 0.999:
#####################################################
cor(X[,1],X[,2])
[1] -0.03066298
#####################################################
Is this a bug, or am I misunderstanding susie_get_cs()? It is not a big issue as I usually use the output of the main susie function to get the credible sets.

Best,
Yihe Yang

[pcarbo]

@harryyiheyang This is not a bug: min_abs_corr will have no effect unless you provide X or Xcorr, which are needed to compute the correlations among the SNPs. (We should have a warning about this when you do not provide X or Xcorr; this is a mistake in the interface.)

Also, in general, we do not recommend setting min_abs_corr to such a high number because it risks filtering out "useful" CSs.

[harryyiheyang]

Dear Dr. Carbonetto,
Thank you for your explanation. This clarifies why the results remained unchanged regardless of how large or small the min_abs_corr parameter was set.
However, I am still puzzled by the discrepancy between the number of credible sets identified by the susie_get_cs() function versus the main SuSiE function results. Specifically, susie_get_cs() typically identifies additional credible sets, and these extra credible sets, as demonstrated in my previous example, contain an unusually large number of SNPs.
Could you help me understand what might be causing this behavior and whether this is expected under certain conditions?
By the way, I should note again that this is not a big issue. I guess most users will only pay attention to the outputs of the main function, which already define the credible sets.
Best regards,
Yihe

[pcarbo]

@harryyiheyang Yes, so the reason for the discrepancy is that susie() by default applies the "purity" (minimum absolute correlation) filter, with min_abs_corr = 0.5. This would explain the discrepancy. (Again, this could be better documented in the package, so it isn't necessarily your fault for missing this.)

[harryyiheyang]

Thank you. Now I understand the discrepancy.

[pcarbo]

@harryyiheyang I'm going to leave this issue open as a reminder to clarify these points in the documentation and add a helpful warning message. Thanks.