is there interest in supporting this package? get_bgc_outliers big bug

[GabrieleNocchi]

is there any interest in supporting this package?
I have found a few bugs.

The get_bgc_outliers function has a very odd behaviour. The test implemented to get loci falling outside the interval does not work as intended, I think. I am trying to fix it but I haven't managed it yet.

G.

[GabrieleNocchi]

The problem is inside the get_ab_outliers function which is used by get_bgc_outliers.

The problem is here, I explain it for alpha but the same thing happens with beta:

  a.out$q <- qa.out$mean

  qn_lower <- 1 - qn

  a.out$qlb <- qnorm(qn_lower,0,sqrt(a.out$q))
  a.out$qub <- qnorm(qn,0,sqrt(a.out$q))

 a.out$outlier <- NA
  a.out$outlier[a.out$median > a.out$qub] <- 'pos'
  a.out$outlier[a.out$median < a.out$qlb] <- 'neg'

The problem is the following. The median alpha value a.out$median will always be a value between -1 and +1 , more likely between -0.6 and +0.6

For high a.out$median values (ie. positive outliers), let s say 0.3 and above you tend to get very high qa.out$mean, which becomes a.out$q (first line of code reported). What happens now is that when you you use a high value (close to 1) of a.out$q in qnorm, the resulting upper and lower interval are also large. As an example, a snippet:

         mean  median        lb        ub   qa_mean qa_median     qa_lb     qa_ub
70   0.3217946 0.32319 0.1074225 0.5245565 0.9569750   0.98187 0.7592230 0.9996700
293  0.3296822 0.33390 0.1091840 0.5309080 0.9605573   0.98487 0.7603545 0.9997295
302  0.3051092 0.30743 0.1233210 0.4896975 0.9532425   0.97712 0.7903960 0.9992290
425  0.3202411 0.31816 0.1001015 0.5455615 0.9533614   0.98099 0.7401475 0.9998000
426  0.3697829 0.37236 0.1534330 0.5769270 0.9750515   0.99200 0.8412615 0.9999100
491  0.3002155 0.30176 0.1165860 0.4871995 0.9501103   0.97474 0.7734910 0.9992595
492  0.3847032 0.38473 0.1794085 0.5868635 0.9805302   0.99367 0.8788665 0.9999300
493  0.3633194 0.36485 0.1603090 0.5520700 0.9748470   0.99128 0.8515130 0.9998300
494  0.4702170 0.47138 0.2597860 0.6755700 0.9943327   0.99889 0.9544030 0.9999900

Let s look for example at the first row, where qa$mean is 0.9569750

Those are the upper and lower bound computed with qnorm (rounding to 0.95):

qnorm(0.975,0,sqrt(0.95))
[1] 1.910337
qnorm(0.025,0,sqrt(0.95))
[1] -1.910337

You see, for a positive alpha value it is IMPOSSIBLE for the median to be outside this interval, so the crazy.a outliers that satisfy both tests are always ONLY negative alpha.

Whereas let s look at the negative values a.out$median:

     mean   median         lb         ub     qa_mean qa_median     qa_lb
1    -0.3046344 -0.30504 -0.5136270 -0.0711290 0.056005223   0.02387 0.0004110
2    -0.4718997 -0.47040 -0.6827980 -0.2757195 0.005134234   0.00113 0.0000000
32   -0.3171964 -0.31876 -0.5324780 -0.0975345 0.046875223   0.01927 0.0002900
76   -0.3008080 -0.30108 -0.4955820 -0.1087230 0.050142684   0.02586 0.0006115
118  -0.3132045 -0.31366 -0.5036815 -0.1223045 0.042734282   0.02074 0.0005205
119  -0.3006422 -0.30045 -0.4644670 -0.1328295 0.043925122   0.02562 0.0011900
418  -0.3284622 -0.33061 -0.5186505 -0.1318260 0.035739184   0.01610 0.0004300
419  -0.3361822 -0.33736 -0.5327635 -0.1266835 0.033791618   0.01406 0.0002405

IF you see, negative outliers (ie. a.out$median < -0.3) get a low qa.out$mean -- if we see the first row it is around 0.05. What happens now is:

qnorm(0.025,0,sqrt(0.05))
[1] -0.4382613

qnorm(0.975,0,sqrt(0.05))
[1] 0.4382613

Or second row:

qnorm(0.975,0,sqrt(0.005))
[1] 0.1385904
qnorm(0.025,0,sqrt(0.005))
[1] -0.1385904

You see that for a negative value, it is instead very possible to be outside this interval. But not for a positive value.

I read the paper, as well as the original, and I am unsure what to change. There is something weird in how these gamma quantiles (the qa ) is used, and I do not think is correct.

A quick proof of concept that something is odd is simple to do:

Run an analysis, and then run the same analysis but invert the order of the parental populations (ie. p0 becomes p1 and viceversa).
You will find that the actual alpha value are good, perfectly correlated (of course in opposite direction correlation of -0.9999). The SNPs whose alpha don t overlap 0 also match (of course with an opposite sign). But the SNPs that do not overlap the chosen qn interval will differ between the inverse direction analyses! That s because only negative alpha values can satisfy that interval criterium, and subsequently also crazy.a will differ. So both analyses will only report some negative alpha SNPs as outside the interval, and these of course will not match as the analyses are the same but with inverse direction.

I suggest not to use the interval outlier test but to just look at the alpha.excess (ie. alpha not overlapping 0). I would ignore alpha.outlier and crazy.a fields until this is fixed.

The same probably happens for beta.

Cheers,
G.