1 // This file is part of PLINK 1.90, copyright (C) 2005-2020 Shaun Purcell,
2 // Christopher Chang.
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program. If not, see <http://www.gnu.org/licenses/>.
16
17
18 #include "plink_common.h"
19
20 #include "plink_stats.h"
21 #include "ipmpar.h"
22 #include "dcdflib.h"
23
24 // 2^{-40} for now, since 2^{-44} was too small on real data
25 #define FISHER_EPSILON 0.0000000000009094947017729282379150390625
26
chiprob_p(double xx,double df)27 double chiprob_p(double xx, double df) {
28 int st = 0;
29 int ww = 1;
30 double bnd = 1;
31 double pp;
32 double qq;
33 cdfchi(&ww, &pp, &qq, &xx, &df, &st, &bnd);
34 if (st) {
35 return -9;
36 }
37 return qq;
38 }
39
inverse_chiprob(double qq,double df)40 double inverse_chiprob(double qq, double df) {
41 double pp = 1 - qq;
42 int32_t st = 0;
43 int32_t ww = 2;
44 double bnd = 1;
45 double xx;
46
47 if (qq >= 1.0) {
48 return 0;
49 }
50 cdfchi(&ww, &pp, &qq, &xx, &df, &st, &bnd);
51 if (st != 0) {
52 return -9;
53 }
54 return xx;
55 }
56
calc_tprob(double tt,double df)57 double calc_tprob(double tt, double df) {
58 int32_t st = 0;
59 int32_t ww = 1;
60 double bnd = 1;
61 double pp;
62 double qq;
63 if (!realnum(tt)) {
64 return -9;
65 }
66 tt = fabs(tt);
67 cdft(&ww, &pp, &qq, &tt, &df, &st, &bnd);
68 if (st != 0) {
69 return -9;
70 }
71 return 2 * qq;
72 }
73
inverse_tprob(double dbl_qq,double df)74 double inverse_tprob(double dbl_qq, double df) {
75 double qq = dbl_qq * 0.5;
76 double pp = 1 - qq;
77 int32_t st = 0;
78 int32_t ww = 2;
79 double bnd = 1;
80 double tt;
81 cdft(&ww, &pp, &qq, &tt, &df, &st, &bnd);
82 if (st != 0) {
83 return -9;
84 }
85 return tt;
86 }
87
88 // Inverse normal distribution
89
90 //
91 // Lower tail quantile for standard normal distribution function.
92 //
93 // This function returns an approximation of the inverse cumulative
94 // standard normal distribution function. I.e., given P, it returns
95 // an approximation to the X satisfying P = Pr{Z <= X} where Z is a
96 // random variable from the standard normal distribution.
97 //
98 // The algorithm uses a minimax approximation by rational functions
99 // and the result has a relative error whose absolute value is less
100 // than 1.15e-9.
101 //
102 // Author: Peter J. Acklam
103 // Time-stamp: 2002-06-09 18:45:44 +0200
104 // E-mail: jacklam@math.uio.no
105 // WWW URL: http://www.math.uio.no/~jacklam
106 //
107 // C implementation adapted from Peter's Perl version
108
109 // Coefficients in rational approximations.
110
111 static const double ivn_a[] =
112 {
113 -3.969683028665376e+01,
114 2.209460984245205e+02,
115 -2.759285104469687e+02,
116 1.383577518672690e+02,
117 -3.066479806614716e+01,
118 2.506628277459239e+00
119 };
120
121 static const double ivn_b[] =
122 {
123 -5.447609879822406e+01,
124 1.615858368580409e+02,
125 -1.556989798598866e+02,
126 6.680131188771972e+01,
127 -1.328068155288572e+01
128 };
129
130 static const double ivn_c[] =
131 {
132 -7.784894002430293e-03,
133 -3.223964580411365e-01,
134 -2.400758277161838e+00,
135 -2.549732539343734e+00,
136 4.374664141464968e+00,
137 2.938163982698783e+00
138 };
139
140 static const double ivn_d[] =
141 {
142 7.784695709041462e-03,
143 3.224671290700398e-01,
144 2.445134137142996e+00,
145 3.754408661907416e+00
146 };
147
148 #define IVN_LOW 0.02425
149 #define IVN_HIGH 0.97575
150
ltqnorm(double p)151 double ltqnorm(double p) {
152 // assumes 0 < p < 1
153 double q, r;
154
155 if (p < IVN_LOW) {
156 // Rational approximation for lower region
157 q = sqrt(-2*log(p));
158 return (((((ivn_c[0]*q+ivn_c[1])*q+ivn_c[2])*q+ivn_c[3])*q+ivn_c[4])*q+ivn_c[5]) /
159 ((((ivn_d[0]*q+ivn_d[1])*q+ivn_d[2])*q+ivn_d[3])*q+1);
160 } else if (p > IVN_HIGH) {
161 // Rational approximation for upper region
162 q = sqrt(-2*log(1-p));
163 return -(((((ivn_c[0]*q+ivn_c[1])*q+ivn_c[2])*q+ivn_c[3])*q+ivn_c[4])*q+ivn_c[5]) /
164 ((((ivn_d[0]*q+ivn_d[1])*q+ivn_d[2])*q+ivn_d[3])*q+1);
165 } else {
166 // Rational approximation for central region
167 q = p - 0.5;
168 r = q*q;
169 return (((((ivn_a[0]*r+ivn_a[1])*r+ivn_a[2])*r+ivn_a[3])*r+ivn_a[4])*r+ivn_a[5])*q /
170 (((((ivn_b[0]*r+ivn_b[1])*r+ivn_b[2])*r+ivn_b[3])*r+ivn_b[4])*r+1);
171 }
172 }
173
SNPHWE2(int32_t obs_hets,int32_t obs_hom1,int32_t obs_hom2,uint32_t midp)174 double SNPHWE2(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, uint32_t midp) {
175 // This function implements an exact SNP test of Hardy-Weinberg
176 // Equilibrium as described in Wigginton, JE, Cutler, DJ, and
177 // Abecasis, GR (2005) A Note on Exact Tests of Hardy-Weinberg
178 // Equilibrium. American Journal of Human Genetics. 76: 000 - 000.
179 //
180 // The original version was written by Jan Wigginton.
181 //
182 // This version was written by Christopher Chang. It contains the following
183 // improvements over the original SNPHWE():
184 // - Proper handling of >64k genotypes. Previously, there was a potential
185 // integer overflow.
186 // - Detection and efficient handling of floating point overflow and
187 // underflow. E.g. instead of summing a tail all the way down, the loop
188 // stops once the latest increment underflows the partial sum's 53-bit
189 // precision; this results in a large speedup when max heterozygote count
190 // >1k.
191 // - No malloc() call. It's only necessary to keep track of a few partial
192 // sums.
193 // - Support for the mid-p variant of this test. See Graffelman J, Moreno V
194 // (2013) The mid p-value in exact tests for Hardy-Weinberg equilibrium.
195 //
196 // Note that the SNPHWE_t() function below is a lot more efficient for
197 // testing against a p-value inclusion threshold. SNPHWE2() should only be
198 // used if you need the actual p-value.
199 intptr_t obs_homc;
200 intptr_t obs_homr;
201 if (obs_hom1 < obs_hom2) {
202 obs_homc = obs_hom2;
203 obs_homr = obs_hom1;
204 } else {
205 obs_homc = obs_hom1;
206 obs_homr = obs_hom2;
207 }
208 int64_t rare_copies = 2LL * obs_homr + obs_hets;
209 int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
210 int32_t tie_ct = 1;
211 double curr_hets_t2 = obs_hets;
212 double curr_homr_t2 = obs_homr;
213 double curr_homc_t2 = obs_homc;
214 double tailp = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
215 double centerp = 0;
216 double lastp2 = tailp;
217 double lastp1 = tailp;
218 double curr_hets_t1;
219 double curr_homr_t1;
220 double curr_homc_t1;
221 double preaddp;
222 if (!genotypes2) {
223 if (midp) {
224 return 0.5;
225 } else {
226 return 1;
227 }
228 }
229
230 if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
231 // tail 1 = upper
232 while (curr_hets_t2 > 1.5) {
233 // het_probs[curr_hets] = 1
234 // het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0)
235 curr_homr_t2 += 1;
236 curr_homc_t2 += 1;
237 lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
238 curr_hets_t2 -= 2;
239 if (lastp2 < EXACT_TEST_BIAS) {
240 if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
241 tie_ct++;
242 }
243 tailp += lastp2;
244 break;
245 }
246 centerp += lastp2;
247 if (centerp == INFINITY) {
248 return 0;
249 }
250 }
251 if ((centerp == 0) && (!midp)) {
252 return 1;
253 }
254 while (curr_hets_t2 > 1.5) {
255 curr_homr_t2 += 1;
256 curr_homc_t2 += 1;
257 lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
258 curr_hets_t2 -= 2;
259 preaddp = tailp;
260 tailp += lastp2;
261 if (tailp <= preaddp) {
262 break;
263 }
264 }
265 curr_hets_t1 = obs_hets + 2;
266 curr_homr_t1 = obs_homr;
267 curr_homc_t1 = obs_homc;
268 while (curr_homr_t1 > 0.5) {
269 // het_probs[curr_hets + 2] = het_probs[curr_hets] * 4 * curr_homr * curr_homc / ((curr_hets + 2) * (curr_hets + 1))
270 lastp1 *= (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
271 preaddp = tailp;
272 tailp += lastp1;
273 if (tailp <= preaddp) {
274 break;
275 }
276 curr_hets_t1 += 2;
277 curr_homr_t1 -= 1;
278 curr_homc_t1 -= 1;
279 }
280 } else {
281 // tail 1 = lower
282 while (curr_homr_t2 > 0.5) {
283 curr_hets_t2 += 2;
284 lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
285 curr_homr_t2 -= 1;
286 curr_homc_t2 -= 1;
287 if (lastp2 < EXACT_TEST_BIAS) {
288 if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
289 tie_ct++;
290 }
291 tailp += lastp2;
292 break;
293 }
294 centerp += lastp2;
295 if (centerp == INFINITY) {
296 return 0;
297 }
298 }
299 if ((centerp == 0) && (!midp)) {
300 return 1;
301 }
302 while (curr_homr_t2 > 0.5) {
303 curr_hets_t2 += 2;
304 lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
305 curr_homr_t2 -= 1;
306 curr_homc_t2 -= 1;
307 preaddp = tailp;
308 tailp += lastp2;
309 if (tailp <= preaddp) {
310 break;
311 }
312 }
313 curr_hets_t1 = obs_hets;
314 curr_homr_t1 = obs_homr;
315 curr_homc_t1 = obs_homc;
316 while (curr_hets_t1 > 1.5) {
317 curr_homr_t1 += 1;
318 curr_homc_t1 += 1;
319 lastp1 *= (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
320 preaddp = tailp;
321 tailp += lastp1;
322 if (tailp <= preaddp) {
323 break;
324 }
325 curr_hets_t1 -= 2;
326 }
327 }
328 if (!midp) {
329 return tailp / (tailp + centerp);
330 } else {
331 return (tailp - ((1 - SMALL_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (tailp + centerp);
332 }
333 }
334
SNPHWE_t(int32_t obs_hets,int32_t obs_hom1,int32_t obs_hom2,double thresh)335 int32_t SNPHWE_t(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, double thresh) {
336 // Threshold-test-only version of SNPHWE2() which is usually able to exit
337 // from the calculation earlier. Returns 0 if these counts are close enough
338 // to Hardy-Weinberg equilibrium, 1 otherwise.
339 //
340 // Suppose, for definiteness, that the number of observed hets is no less
341 // than expectation. (Same ideas apply for the other case.) We proceed as
342 // follows:
343 // - Sum the *relative* likelihoods of more likely smaller het counts.
344 // - Determine the minimum tail mass to pass the threshold.
345 // - The majority of the time, the tail boundary elements are enough to pass
346 // the threshold; we never need to sum the remainder of the tails.
347 // - And in the case of disequilibrium, we will often be able to immediately
348 // determine that the tail sum cannot possibly pass the threshold, just by
349 // looking at the tail boundary elements and using a geometric series to
350 // upper-bound the tail sums.
351 // - Only when neither of these conditions hold do we start traveling down
352 // the tails.
353 intptr_t obs_homc;
354 intptr_t obs_homr;
355 if (obs_hom1 < obs_hom2) {
356 obs_homc = obs_hom2;
357 obs_homr = obs_hom1;
358 } else {
359 obs_homc = obs_hom1;
360 obs_homr = obs_hom2;
361 }
362 int64_t rare_copies = 2LL * obs_homr + obs_hets;
363 int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
364 double curr_hets_t2 = obs_hets; // tail 2
365 double curr_homr_t2 = obs_homr;
366 double curr_homc_t2 = obs_homc;
367
368 // Subtract epsilon from initial probability mass, so that we can compare to
369 // 1 when determining tail vs. center membership without floating point error
370 // biting us in the ass
371 double tailp1 = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
372 double centerp = 0;
373 double lastp2 = tailp1;
374 double tailp2 = 0;
375 double tail1_ceil;
376 double tail2_ceil;
377 double lastp1;
378 double curr_hets_t1;
379 double curr_homr_t1;
380 double curr_homc_t1;
381
382 // Initially, if center sum reaches this, the test can immediately fail.
383 // Once center is summed, this is recalculated, and when tail sum has reached
384 // this, we've passed.
385 double exit_thresh;
386 double exit_threshx;
387 double ratio;
388 double preaddp;
389 if (!genotypes2) {
390 return 0;
391 }
392
393 // Convert thresh into reverse odds ratio.
394 thresh = (1 - thresh) / thresh;
395
396 // Expected het count:
397 // 2 * rarefreq * (1 - rarefreq) * genotypes
398 // = 2 * (rare_copies / (2 * genotypes)) * (1 - rarefreq) * genotypes
399 // = rare_copies * (1 - (rare_copies / (2 * genotypes)))
400 // = (rare_copies * (2 * genotypes - rare_copies)) / (2 * genotypes)
401 //
402 // The computational identity is
403 // P(nhets == n) := P(nhets == n+2) * (n+2) * (n+1) /
404 // (4 * homr(n) * homc(n))
405 // where homr() and homc() are the number of homozygous rares/commons needed
406 // to maintain the same allele frequencies.
407 // This probability is always decreasing when proceeding away from the
408 // expected het count.
409
410 if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
411 // tail 1 = upper
412 if (obs_hets < 2) {
413 return 0;
414 }
415
416 // An initial upper bound on the tail sum is useful, since it lets us
417 // report test failure before summing the entire center. We use the
418 // trivial bound of 1 + floor(rare_copies / 2): that's the total number
419 // of possible het counts, and the relative probability for each count must
420 // be <= 1 if it's in the tail.
421 exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
422
423 // het_probs[curr_hets] = 1
424 // het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1) / (4 * (curr_homr + 1) * (curr_homc + 1))
425 do {
426 curr_homr_t2 += 1;
427 curr_homc_t2 += 1;
428 lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
429 curr_hets_t2 -= 2;
430 if (lastp2 < EXACT_TEST_BIAS) {
431 tailp2 = lastp2;
432 break;
433 }
434 centerp += lastp2;
435 if (centerp > exit_thresh) {
436 return 1;
437 }
438 } while (curr_hets_t2 > 1.5);
439 exit_thresh = centerp / thresh;
440 if (tailp1 + tailp2 >= exit_thresh) {
441 return 0;
442 }
443 // c + cr + cr^2 + ... = c/(1-r), which is an upper bound for the tail sum
444 ratio = (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * (curr_homr_t2 + 1) * (curr_homc_t2 + 1));
445 tail2_ceil = tailp2 / (1 - ratio);
446 curr_hets_t1 = obs_hets + 2;
447 curr_homr_t1 = obs_homr;
448 curr_homc_t1 = obs_homc;
449 // ratio for the other tail
450 lastp1 = (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
451 tail1_ceil = tailp1 / (1 - lastp1);
452 if (tail1_ceil + tail2_ceil < exit_thresh) {
453 return 1;
454 }
455 lastp1 *= tailp1;
456 tailp1 += lastp1;
457
458 if (obs_homr > 1) {
459 // het_probs[curr_hets + 2] = het_probs[curr_hets] * 4 * curr_homr * curr_homc / ((curr_hets + 2) * (curr_hets + 1))
460 exit_threshx = exit_thresh - tailp2;
461 do {
462 curr_hets_t1 += 2;
463 curr_homr_t1 -= 1;
464 curr_homc_t1 -= 1;
465 lastp1 *= (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
466 preaddp = tailp1;
467 tailp1 += lastp1;
468 if (tailp1 > exit_threshx) {
469 return 0;
470 }
471 if (tailp1 <= preaddp) {
472 break;
473 }
474 } while (curr_homr_t1 > 1.5);
475 }
476 if (tailp1 + tail2_ceil < exit_thresh) {
477 return 1;
478 }
479 exit_threshx = exit_thresh - tailp1;
480 while (curr_hets_t2 > 1) {
481 curr_homr_t2 += 1;
482 curr_homc_t2 += 1;
483 lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
484 preaddp = tailp2;
485 tailp2 += lastp2;
486 if (tailp2 >= exit_threshx) {
487 return 0;
488 }
489 if (tailp2 <= preaddp) {
490 return 1;
491 }
492 curr_hets_t2 -= 2;
493 }
494 return 1;
495 } else {
496 // tail 1 = lower
497 if (!obs_homr) {
498 return 0;
499 }
500 exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
501 do {
502 curr_hets_t2 += 2;
503 lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
504 curr_homr_t2 -= 1;
505 curr_homc_t2 -= 1;
506 if (lastp2 < EXACT_TEST_BIAS) {
507 tailp2 = lastp2;
508 break;
509 }
510 centerp += lastp2;
511 if (centerp > exit_thresh) {
512 return 1;
513 }
514 } while (curr_homr_t2 > 0.5);
515 exit_thresh = centerp / thresh;
516 if (tailp1 + tailp2 >= exit_thresh) {
517 return 0;
518 }
519 ratio = (4 * curr_homr_t2 * curr_homc_t2) / ((curr_hets_t2 + 2) * (curr_hets_t2 + 1));
520 tail2_ceil = tailp2 / (1 - ratio);
521 curr_hets_t1 = obs_hets;
522 curr_homr_t1 = obs_homr + 1;
523 curr_homc_t1 = obs_homc + 1;
524 lastp1 = (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
525 tail1_ceil = tailp1 / (1 - lastp1);
526 lastp1 *= tailp1;
527 tailp1 += lastp1;
528
529 if (tail1_ceil + tail2_ceil < exit_thresh) {
530 return 1;
531 }
532 if (obs_hets >= 4) {
533 exit_threshx = exit_thresh - tailp2;
534 do {
535 curr_hets_t1 -= 2;
536 curr_homr_t1 += 1;
537 curr_homc_t1 += 1;
538 lastp1 *= (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
539 preaddp = tailp1;
540 tailp1 += lastp1;
541 if (tailp1 > exit_threshx) {
542 return 0;
543 }
544 if (tailp1 <= preaddp) {
545 break;
546 }
547 } while (curr_hets_t1 > 3.5);
548 }
549 if (tailp1 + tail2_ceil < exit_thresh) {
550 return 1;
551 }
552 exit_threshx = exit_thresh - tailp1;
553 while (curr_homr_t2 > 0.5) {
554 curr_hets_t2 += 2;
555 lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
556 curr_homr_t2 -= 1;
557 curr_homc_t2 -= 1;
558 preaddp = tailp2;
559 tailp2 += lastp2;
560 if (tailp2 >= exit_threshx) {
561 return 0;
562 }
563 if (tailp2 <= preaddp) {
564 return 1;
565 }
566 }
567 return 1;
568 }
569 }
570
SNPHWE_midp_t(int32_t obs_hets,int32_t obs_hom1,int32_t obs_hom2,double thresh)571 int32_t SNPHWE_midp_t(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, double thresh) {
572 // Mid-p version of SNPHWE_t(). (There are enough fiddly differences that I
573 // think it's better for this to be a separate function.) Assumes threshold
574 // is smaller than 0.5.
575 intptr_t obs_homc;
576 intptr_t obs_homr;
577 if (obs_hom1 < obs_hom2) {
578 obs_homc = obs_hom2;
579 obs_homr = obs_hom1;
580 } else {
581 obs_homc = obs_hom1;
582 obs_homr = obs_hom2;
583 }
584 int64_t rare_copies = 2LL * obs_homr + obs_hets;
585 int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
586 double curr_hets_t2 = obs_hets; // tail 2
587 double curr_homr_t2 = obs_homr;
588 double curr_homc_t2 = obs_homc;
589 double tailp1 = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS * 0.5;
590 double centerp = tailp1;
591 double lastp2 = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
592 double tailp2 = 0;
593 double tail1_ceil;
594 double tail2_ceil;
595 double lastp1;
596 double curr_hets_t1;
597 double curr_homr_t1;
598 double curr_homc_t1;
599 double exit_thresh;
600 double exit_threshx;
601 double ratio;
602 double preaddp;
603 if (!genotypes2) {
604 return 0;
605 }
606 thresh = (1 - thresh) / thresh;
607 if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
608 if (obs_hets < 2) {
609 return 0;
610 }
611 exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
612 do {
613 curr_homr_t2 += 1;
614 curr_homc_t2 += 1;
615 lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
616 curr_hets_t2 -= 2;
617 if (lastp2 < EXACT_TEST_BIAS) {
618 if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
619 // tie with original contingency table, apply mid-p correction here
620 // too
621 tailp2 = tailp1;
622 centerp += tailp1;
623 } else {
624 tailp2 = lastp2;
625 }
626 break;
627 }
628 centerp += lastp2;
629 if (centerp > exit_thresh) {
630 return 1;
631 }
632 } while (curr_hets_t2 > 1.5);
633 exit_thresh = centerp / thresh;
634 if (tailp1 + tailp2 >= exit_thresh) {
635 return 0;
636 }
637 ratio = (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * (curr_homr_t2 + 1) * (curr_homc_t2 + 1));
638 // this needs to work in both the tie and no-tie cases
639 tail2_ceil = tailp2 + lastp2 * ratio / (1 - ratio);
640 curr_hets_t1 = obs_hets + 2;
641 curr_homr_t1 = obs_homr;
642 curr_homc_t1 = obs_homc;
643 lastp1 = (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
644 // always a tie here
645 tail1_ceil = tailp1 * 2 / (1 - lastp1) - tailp1;
646 if (tail1_ceil + tail2_ceil < exit_thresh) {
647 return 1;
648 }
649 lastp1 *= tailp1 * 2;
650 tailp1 += lastp1;
651
652 if (obs_homr > 1) {
653 exit_threshx = exit_thresh - tailp2;
654 do {
655 curr_hets_t1 += 2;
656 curr_homr_t1 -= 1;
657 curr_homc_t1 -= 1;
658 lastp1 *= (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
659 preaddp = tailp1;
660 tailp1 += lastp1;
661 if (tailp1 > exit_threshx) {
662 return 0;
663 }
664 if (tailp1 <= preaddp) {
665 break;
666 }
667 } while (curr_homr_t1 > 1.5);
668 }
669 if (tailp1 + tail2_ceil < exit_thresh) {
670 return 1;
671 }
672 exit_threshx = exit_thresh - tailp1;
673 while (curr_hets_t2 > 1) {
674 curr_homr_t2 += 1;
675 curr_homc_t2 += 1;
676 lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
677 preaddp = tailp2;
678 tailp2 += lastp2;
679 if (tailp2 >= exit_threshx) {
680 return 0;
681 }
682 if (tailp2 <= preaddp) {
683 return 1;
684 }
685 curr_hets_t2 -= 2;
686 }
687 return 1;
688 } else {
689 if (!obs_homr) {
690 return 0;
691 }
692 exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
693 do {
694 curr_hets_t2 += 2;
695 lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
696 curr_homr_t2 -= 1;
697 curr_homc_t2 -= 1;
698 if (lastp2 < EXACT_TEST_BIAS) {
699 if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
700 tailp2 = tailp1;
701 centerp += tailp1;
702 } else {
703 tailp2 = lastp2;
704 }
705 break;
706 }
707 centerp += lastp2;
708 if (centerp > exit_thresh) {
709 return 1;
710 }
711 } while (curr_homr_t2 > 0.5);
712 exit_thresh = centerp / thresh;
713 if (tailp1 + tailp2 >= exit_thresh) {
714 return 0;
715 }
716 ratio = (4 * curr_homr_t2 * curr_homc_t2) / ((curr_hets_t2 + 2) * (curr_hets_t2 + 1));
717 tail2_ceil = tailp2 + lastp2 * ratio / (1 - ratio);
718 curr_hets_t1 = obs_hets;
719 curr_homr_t1 = obs_homr + 1;
720 curr_homc_t1 = obs_homc + 1;
721 lastp1 = (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
722 tail1_ceil = 2 * tailp1 / (1 - lastp1) - tailp1;
723 lastp1 *= 2 * tailp1;
724 tailp1 += lastp1;
725
726 if (tail1_ceil + tail2_ceil < exit_thresh) {
727 return 1;
728 }
729 if (obs_hets >= 4) {
730 exit_threshx = exit_thresh - tailp2;
731 do {
732 curr_hets_t1 -= 2;
733 curr_homr_t1 += 1;
734 curr_homc_t1 += 1;
735 lastp1 *= (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
736 preaddp = tailp1;
737 tailp1 += lastp1;
738 if (tailp1 > exit_threshx) {
739 return 0;
740 }
741 if (tailp1 <= preaddp) {
742 break;
743 }
744 } while (curr_hets_t1 > 3.5);
745 }
746 if (tailp1 + tail2_ceil < exit_thresh) {
747 return 1;
748 }
749 exit_threshx = exit_thresh - tailp1;
750 while (curr_homr_t2 > 0.5) {
751 curr_hets_t2 += 2;
752 lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
753 curr_homr_t2 -= 1;
754 curr_homc_t2 -= 1;
755 preaddp = tailp2;
756 tailp2 += lastp2;
757 if (tailp2 >= exit_threshx) {
758 return 0;
759 }
760 if (tailp2 <= preaddp) {
761 return 1;
762 }
763 }
764 return 1;
765 }
766 }
767
fisher22(uint32_t m11,uint32_t m12,uint32_t m21,uint32_t m22,uint32_t midp)768 double fisher22(uint32_t m11, uint32_t m12, uint32_t m21, uint32_t m22, uint32_t midp) {
769 // Basic 2x2 Fisher exact test p-value calculation.
770 double tprob = (1 - FISHER_EPSILON) * EXACT_TEST_BIAS;
771 double cur_prob = tprob;
772 double cprob = 0;
773 int32_t tie_ct = 1;
774 uint32_t uii;
775 double cur11;
776 double cur12;
777 double cur21;
778 double cur22;
779 double preaddp;
780 // Ensure we are left of the distribution center, m11 <= m22, and m12 <= m21.
781 if (m12 > m21) {
782 uii = m12;
783 m12 = m21;
784 m21 = uii;
785 }
786 if (m11 > m22) {
787 uii = m11;
788 m11 = m22;
789 m22 = uii;
790 }
791 if ((((uint64_t)m11) * m22) > (((uint64_t)m12) * m21)) {
792 uii = m11;
793 m11 = m12;
794 m12 = uii;
795 uii = m21;
796 m21 = m22;
797 m22 = uii;
798 }
799 cur11 = m11;
800 cur12 = m12;
801 cur21 = m21;
802 cur22 = m22;
803 while (cur12 > 0.5) {
804 cur11 += 1;
805 cur22 += 1;
806 cur_prob *= (cur12 * cur21) / (cur11 * cur22);
807 cur12 -= 1;
808 cur21 -= 1;
809 if (cur_prob == INFINITY) {
810 return 0;
811 }
812 if (cur_prob < EXACT_TEST_BIAS) {
813 if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
814 tie_ct++;
815 }
816 tprob += cur_prob;
817 break;
818 }
819 cprob += cur_prob;
820 }
821 if ((cprob == 0) && (!midp)) {
822 return 1;
823 }
824 while (cur12 > 0.5) {
825 cur11 += 1;
826 cur22 += 1;
827 cur_prob *= (cur12 * cur21) / (cur11 * cur22);
828 cur12 -= 1;
829 cur21 -= 1;
830 preaddp = tprob;
831 tprob += cur_prob;
832 if (tprob <= preaddp) {
833 break;
834 }
835 }
836 if (m11) {
837 cur11 = m11;
838 cur12 = m12;
839 cur21 = m21;
840 cur22 = m22;
841 cur_prob = (1 - FISHER_EPSILON) * EXACT_TEST_BIAS;
842 do {
843 cur12 += 1;
844 cur21 += 1;
845 cur_prob *= (cur11 * cur22) / (cur12 * cur21);
846 cur11 -= 1;
847 cur22 -= 1;
848 preaddp = tprob;
849 tprob += cur_prob;
850 if (tprob <= preaddp) {
851 if (!midp) {
852 return preaddp / (cprob + preaddp);
853 } else {
854 return (preaddp - ((1 - FISHER_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (cprob + preaddp);
855 }
856 }
857 } while (cur11 > 0.5);
858 }
859 if (!midp) {
860 return tprob / (cprob + tprob);
861 } else {
862 return (tprob - ((1 - FISHER_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (cprob + tprob);
863 }
864 }
865
fisher22_tail_pval(uint32_t m11,uint32_t m12,uint32_t m21,uint32_t m22,int32_t right_offset,double tot_prob_recip,double right_prob,uint32_t midp,uint32_t new_m11)866 double fisher22_tail_pval(uint32_t m11, uint32_t m12, uint32_t m21, uint32_t m22, int32_t right_offset, double tot_prob_recip, double right_prob, uint32_t midp, uint32_t new_m11) {
867 // Given that the left (w.r.t. m11) reference contingency table has
868 // likelihood 1/tot_prob, the contingency table with m11 increased by
869 // right_offset has likelihood right_prob/tot_prob, and the tails (up to but
870 // not including the two references) sum to tail_sum/tot_prob, this
871 // calculates the p-value of the given m11 (which must be on one tail).
872 double left_prob = 1.0;
873 double dxx = ((intptr_t)new_m11);
874 double cur11;
875 double cur12;
876 double cur21;
877 double cur22;
878 double psum;
879 double thresh;
880 if (new_m11 < m11) {
881 cur11 = ((intptr_t)m11);
882 cur12 = ((intptr_t)m12);
883 cur21 = ((intptr_t)m21);
884 cur22 = ((intptr_t)m22);
885 dxx += 0.5; // unnecessary (53 vs. 32 bits precision), but whatever
886 do {
887 cur12 += 1;
888 cur21 += 1;
889 left_prob *= cur11 * cur22 / (cur12 * cur21);
890 cur11 -= 1;
891 cur22 -= 1;
892 } while (cur11 > dxx);
893 if (left_prob == 0) {
894 return 0;
895 }
896 if (!midp) {
897 psum = left_prob;
898 } else {
899 psum = left_prob * 0.5;
900 }
901 thresh = left_prob * (1 + FISHER_EPSILON);
902 do {
903 if (cur11 < 0.5) {
904 break;
905 }
906 cur12 += 1;
907 cur21 += 1;
908 left_prob *= cur11 * cur22 / (cur12 * cur21);
909 cur11 -= 1;
910 cur22 -= 1;
911 dxx = psum;
912 psum += left_prob;
913 } while (psum > dxx);
914 cur11 = ((intptr_t)(m11 + right_offset));
915 cur12 = ((intptr_t)(m12 - right_offset));
916 cur21 = ((intptr_t)(m21 - right_offset));
917 cur22 = ((intptr_t)(m22 + right_offset));
918 while (right_prob > thresh) {
919 cur11 += 1;
920 cur22 += 1;
921 right_prob *= cur12 * cur21 / (cur11 * cur22);
922 cur12 -= 1;
923 cur21 -= 1;
924 }
925 if (right_prob > 0) {
926 if (midp && (right_prob < thresh * (1 - 2 * FISHER_EPSILON))) {
927 psum += right_prob * 0.5;
928 } else {
929 psum += right_prob;
930 }
931 do {
932 cur11 += 1;
933 cur22 += 1;
934 right_prob *= cur12 * cur21 / (cur11 * cur22);
935 cur12 -= 1;
936 cur21 -= 1;
937 dxx = psum;
938 psum += right_prob;
939 } while (psum > dxx);
940 }
941 } else {
942 dxx -= 0.5;
943 cur11 = ((intptr_t)(m11 + right_offset));
944 cur12 = ((intptr_t)(m12 - right_offset));
945 cur21 = ((intptr_t)(m21 - right_offset));
946 cur22 = ((intptr_t)(m22 + right_offset));
947 do {
948 cur11 += 1;
949 cur22 += 1;
950 right_prob *= cur12 * cur21 / (cur11 * cur22);
951 cur12 -= 1;
952 cur21 -= 1;
953 } while (cur11 < dxx);
954 if (right_prob == 0) {
955 return 0;
956 }
957 if (!midp) {
958 psum = right_prob;
959 } else {
960 psum = right_prob * 0.5;
961 }
962 thresh = right_prob * (1 + FISHER_EPSILON);
963 do {
964 if (cur12 < 0.5) {
965 break;
966 }
967 cur11 += 1;
968 cur22 += 1;
969 right_prob *= cur12 * cur21 / (cur11 * cur22);
970 cur12 -= 1;
971 cur21 -= 1;
972 dxx = psum;
973 psum += right_prob;
974 } while (psum > dxx);
975 cur11 = ((intptr_t)m11);
976 cur12 = ((intptr_t)m12);
977 cur21 = ((intptr_t)m21);
978 cur22 = ((intptr_t)m22);
979 while (left_prob > thresh) {
980 cur12 += 1;
981 cur21 += 1;
982 left_prob *= cur11 * cur22 / (cur12 * cur21);
983 cur11 -= 1;
984 cur22 -= 1;
985 }
986 if (left_prob > 0) {
987 if (midp && (left_prob < thresh * (1 - 2 * FISHER_EPSILON))) {
988 psum += left_prob * 0.5;
989 } else {
990 psum += left_prob;
991 }
992 do {
993 cur12 += 1;
994 cur21 += 1;
995 left_prob *= cur11 * cur22 / (cur12 * cur21);
996 cur11 -= 1;
997 cur22 -= 1;
998 dxx = psum;
999 psum += left_prob;
1000 } while (psum > dxx);
1001 }
1002 }
1003 return psum * tot_prob_recip;
1004 }
1005
fisher22_precomp_pval_bounds(double pval,uint32_t midp,uint32_t row1_sum,uint32_t col1_sum,uint32_t total,uint32_t * bounds,double * tprobs)1006 void fisher22_precomp_pval_bounds(double pval, uint32_t midp, uint32_t row1_sum, uint32_t col1_sum, uint32_t total, uint32_t* bounds, double* tprobs) {
1007 // bounds[0] = m11 min
1008 // bounds[1] = m11 (max + 1)
1009 // bounds[2] = m11 min after including ties
1010 // bounds[3] = m11 (max + 1) after including ties
1011 // Treating m11 as the only variable, this returns the minimum and (maximum +
1012 // 1) values of m11 which are less extreme than the observed result, and
1013 // notes ties (2^{-40} tolerance). Also, returns the values necessary for
1014 // invoking fisher22_tail_pval() with
1015 // m11 := bounds[2] and
1016 // right_offset := bounds[3] - bounds[2] - 1
1017 // in tprobs[0], [1], and [2] (if tprobs is not NULL).
1018 //
1019 // Algorithm:
1020 // 1. Determine center.
1021 // 2. Sum unscaled probabilities in both directions to precision limit.
1022 // 3. Proceed further outwards to (pval * unscaled_psum) precision limit,
1023 // fill in the remaining return values.
1024 double tot_prob = 1.0 / EXACT_TEST_BIAS;
1025 double left_prob = tot_prob;
1026 double right_prob = tot_prob;
1027 intptr_t m11_offset = 0;
1028 double tail_prob = 0;
1029 double cmult = midp? 0.5 : 1.0;
1030 double dxx;
1031 double left11;
1032 double left12;
1033 double left21;
1034 double left22;
1035 double right11;
1036 double right12;
1037 double right21;
1038 double right22;
1039 double cur_prob;
1040 double cur11;
1041 double cur12;
1042 double cur21;
1043 double cur22;
1044 double threshold;
1045 intptr_t lii;
1046 uint32_t uii;
1047 if (!total) {
1048 // hardcode this to avoid divide-by-zero
1049 bounds[0] = 0;
1050 bounds[1] = 0;
1051 bounds[2] = 0;
1052 bounds[3] = 1;
1053 // no need to initialize the other return values, they should never be used
1054 return;
1055 } else {
1056 if (pval == 0) {
1057 if (total >= row1_sum + col1_sum) {
1058 bounds[0] = 0;
1059 bounds[1] = MINV(row1_sum, col1_sum) + 1;
1060 } else {
1061 bounds[0] = row1_sum + col1_sum - total;
1062 bounds[1] = total - MAXV(row1_sum, col1_sum) + 1;
1063 }
1064 bounds[2] = bounds[0];
1065 bounds[3] = bounds[1];
1066 return;
1067 }
1068 }
1069 // Center must be adjacent to the x which satisfies
1070 // (m11 + x)(m22 + x) = (m12 - x)(m21 - x), so
1071 // x = (m12 * m21 - m11 * m22) / (m11 + m12 + m21 + m22)
1072 if (total >= row1_sum + col1_sum) {
1073 // m11 = 0;
1074 // m12 = row1_sum;
1075 // m21 = col1_sum;
1076 // m22 = total - row1_sum - col1_sum;
1077 lii = (((uint64_t)row1_sum) * col1_sum) / total;
1078 left11 = lii;
1079 left12 = row1_sum - lii;
1080 left21 = col1_sum - lii;
1081 left22 = (total - row1_sum - col1_sum) + lii;
1082 } else {
1083 // m11 = row1_sum + col1_sum - total;
1084 // m12 = row1_sum - m11;
1085 // m21 = col1_sum - m11;
1086 // m22 = 0;
1087 lii = (((uint64_t)(total - row1_sum)) * (total - col1_sum)) / total;
1088 // Force m11 <= m22 for internal calculation, then adjust at end.
1089 m11_offset = row1_sum + col1_sum - total;
1090 left11 = lii;
1091 left12 = total - col1_sum - lii;
1092 left21 = total - row1_sum - lii;
1093 left22 = m11_offset + lii;
1094 }
1095 // We rounded x down. Should we have rounded up instead?
1096 if ((left11 + 1) * (left22 + 1) < left12 * left21) {
1097 left11 += 1;
1098 left12 -= 1;
1099 left21 -= 1;
1100 left22 += 1;
1101 }
1102 // Safe to force m12 <= m21.
1103 if (left12 > left21) {
1104 dxx = left12;
1105 left12 = left21;
1106 left21 = dxx;
1107 }
1108 // Sum right side to limit, then left.
1109 right11 = left11;
1110 right12 = left12;
1111 right21 = left21;
1112 right22 = left22;
1113 do {
1114 if (right12 < 0.5) {
1115 break;
1116 }
1117 right11 += 1;
1118 right22 += 1;
1119 right_prob *= (right12 * right21) / (right11 * right22);
1120 right12 -= 1;
1121 right21 -= 1;
1122 dxx = tot_prob;
1123 tot_prob += right_prob;
1124 } while (tot_prob > dxx);
1125 do {
1126 if (left11 < 0.5) {
1127 break;
1128 }
1129 left12 += 1;
1130 left21 += 1;
1131 left_prob *= (left11 * left22) / (left12 * left21);
1132 left11 -= 1;
1133 left22 -= 1;
1134 dxx = tot_prob;
1135 tot_prob += left_prob;
1136 } while (tot_prob > dxx);
1137 // Now traverse the tails to p-value precision limit.
1138 // Upper bound for tail sum, if current element c is included:
1139 // (c + cr + cr^2 + ...) + (c + cs + cs^2 + ...)
1140 // = c(1/(1 - r) + 1/(1 - s))
1141 // Compare this to pval * tot_prob.
1142 // I.e. compare c to pval * tot_prob * (1-r)(1-s) / (2-r-s)
1143 dxx = 1 - (left11 * left22) / ((left12 + 1) * (left21 + 1));
1144 threshold = 1 - (right12 * right21) / ((right11 + 1) * (right22 + 1));
1145 threshold = pval * tot_prob * dxx * threshold / (dxx + threshold);
1146 while (left11 > 0.5) {
1147 if (left_prob < threshold) {
1148 tail_prob = left_prob;
1149 cur11 = left11;
1150 cur12 = left12;
1151 cur21 = left21;
1152 cur22 = left22;
1153 cur_prob = left_prob;
1154 do {
1155 cur12 += 1;
1156 cur21 += 1;
1157 cur_prob *= (cur11 * cur22) / (cur12 * cur21);
1158 cur11 -= 1;
1159 cur22 -= 1;
1160 dxx = tail_prob;
1161 tail_prob += cur_prob;
1162 } while (dxx < tail_prob);
1163 left11 += 1;
1164 left22 += 1;
1165 left_prob *= (left12 * left21) / (left11 * left22);
1166 left12 -= 1;
1167 left21 -= 1;
1168 break;
1169 }
1170 left12 += 1;
1171 left21 += 1;
1172 left_prob *= (left11 * left22) / (left12 * left21);
1173 left11 -= 1;
1174 left22 -= 1;
1175 }
1176 while (right12 > 0.5) {
1177 if (right_prob < threshold) {
1178 tail_prob += right_prob;
1179 cur11 = right11;
1180 cur12 = right12;
1181 cur21 = right21;
1182 cur22 = right22;
1183 cur_prob = right_prob;
1184 do {
1185 cur11 += 1;
1186 cur22 += 1;
1187 cur_prob *= (cur12 * cur21) / (cur11 * cur22);
1188 cur12 -= 1;
1189 cur21 -= 1;
1190 dxx = tail_prob;
1191 tail_prob += cur_prob;
1192 } while (dxx < tail_prob);
1193 right12 += 1;
1194 right21 += 1;
1195 right_prob *= (right11 * right22) / (right12 * right21);
1196 right11 -= 1;
1197 right22 -= 1;
1198 break;
1199 }
1200 right11 += 1;
1201 right22 += 1;
1202 right_prob *= (right12 * right21) / (right11 * right22);
1203 right12 -= 1;
1204 right21 -= 1;
1205 }
1206 dxx = pval * tot_prob * (1 - FISHER_EPSILON / 2);
1207 threshold = pval * tot_prob * (1 + FISHER_EPSILON / 2);
1208 lii = 0;
1209 while (1) {
1210 if (left_prob < right_prob * (1 - FISHER_EPSILON / 2)) {
1211 cur_prob = tail_prob + left_prob * cmult;
1212 if (cur_prob > threshold) {
1213 break;
1214 }
1215 tail_prob += left_prob;
1216 uii = 1;
1217 } else if (right_prob < left_prob * (1 - FISHER_EPSILON / 2)) {
1218 cur_prob = tail_prob + right_prob * cmult;
1219 if (cur_prob > threshold) {
1220 break;
1221 }
1222 tail_prob += right_prob;
1223 uii = 2;
1224 } else {
1225 cur_prob = tail_prob + (left_prob + right_prob) * cmult;
1226 if (cur_prob > threshold) {
1227 if (left11 == right11) {
1228 cur_prob = tail_prob + left_prob * cmult;
1229 // center: left and right refer to same table. subcases:
1230 // 1. cur_prob > threshold: center table has less extreme pval.
1231 // lii = 0, both intervals size 1
1232 // 2. dxx < cur_prob < threshold: center table has equal pval.
1233 // lii = 1, less-than interval size 0 but leq interval size 1
1234 // 3. cur_prob < dxx: even centermost table has more extreme pval
1235 // (only possible with mid-p adj).
1236 // lii = 0, we increment left11 so both intervals size 0
1237 if (cur_prob < threshold) {
1238 if (cur_prob > dxx) {
1239 lii = 1;
1240 } else {
1241 left11++;
1242 left22++;
1243 left_prob *= (left12 * left21) / (left11 * left22);
1244 }
1245 }
1246 }
1247 break;
1248 }
1249 tail_prob += left_prob + right_prob;
1250 uii = 3;
1251 }
1252 if (cur_prob > dxx) {
1253 lii = uii;
1254 break;
1255 }
1256 // if more speed is necessary, we could use a buffer to save all unscaled
1257 // probabilities during the initial outward traversal.
1258 if (uii & 1) {
1259 left11 += 1;
1260 left22 += 1;
1261 left_prob *= (left12 * left21) / (left11 * left22);
1262 left12 -= 1;
1263 left21 -= 1;
1264 }
1265 if (uii & 2) {
1266 right12 += 1;
1267 right21 += 1;
1268 right_prob *= (right11 * right22) / (right12 * right21);
1269 right11 -= 1;
1270 right22 -= 1;
1271 }
1272 }
1273 bounds[2] = m11_offset + ((intptr_t)left11);
1274 bounds[3] = m11_offset + ((intptr_t)right11) + 1;
1275 bounds[0] = bounds[2] + (lii & 1);
1276 bounds[1] = bounds[3] - (lii >> 1);
1277 if (!tprobs) {
1278 return;
1279 }
1280 dxx = 1.0 / left_prob;
1281 tprobs[0] = left_prob / tot_prob;
1282 tprobs[1] = right_prob * dxx;
1283 /*
1284 if (lii & 1) {
1285 tail_prob -= left_prob;
1286 }
1287 if (lii >> 1) {
1288 tail_prob -= right_prob;
1289 }
1290 tprobs[2] = tail_prob * dxx;
1291 */
1292 }
1293
fisher23_tailsum(double * base_probp,double * saved12p,double * saved13p,double * saved22p,double * saved23p,double * totalp,uint32_t * tie_ctp,uint32_t right_side)1294 int32_t fisher23_tailsum(double* base_probp, double* saved12p, double* saved13p, double* saved22p, double* saved23p, double *totalp, uint32_t* tie_ctp, uint32_t right_side) {
1295 double total = 0;
1296 double cur_prob = *base_probp;
1297 double tmp12 = *saved12p;
1298 double tmp13 = *saved13p;
1299 double tmp22 = *saved22p;
1300 double tmp23 = *saved23p;
1301 double tmps12;
1302 double tmps13;
1303 double tmps22;
1304 double tmps23;
1305 double prev_prob;
1306 // identify beginning of tail
1307 if (right_side) {
1308 if (cur_prob > EXACT_TEST_BIAS) {
1309 prev_prob = tmp13 * tmp22;
1310 while (prev_prob > 0.5) {
1311 tmp12 += 1;
1312 tmp23 += 1;
1313 cur_prob *= prev_prob / (tmp12 * tmp23);
1314 tmp13 -= 1;
1315 tmp22 -= 1;
1316 if (cur_prob <= EXACT_TEST_BIAS) {
1317 break;
1318 }
1319 prev_prob = tmp13 * tmp22;
1320 }
1321 *base_probp = cur_prob;
1322 tmps12 = tmp12;
1323 tmps13 = tmp13;
1324 tmps22 = tmp22;
1325 tmps23 = tmp23;
1326 } else {
1327 tmps12 = tmp12;
1328 tmps13 = tmp13;
1329 tmps22 = tmp22;
1330 tmps23 = tmp23;
1331 while (1) {
1332 prev_prob = cur_prob;
1333 tmp13 += 1;
1334 tmp22 += 1;
1335 cur_prob *= (tmp12 * tmp23) / (tmp13 * tmp22);
1336 if (cur_prob < prev_prob) {
1337 return 1;
1338 }
1339 tmp12 -= 1;
1340 tmp23 -= 1;
1341 if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1342 // throw in extra (1 - SMALL_EPSILON) multiplier to prevent rounding
1343 // errors from causing this to keep going when the left-side test
1344 // stopped
1345 if (cur_prob > (1 - SMALL_EPSILON) * EXACT_TEST_BIAS) {
1346 break;
1347 }
1348 *tie_ctp += 1;
1349 }
1350 total += cur_prob;
1351 }
1352 prev_prob = cur_prob;
1353 cur_prob = *base_probp;
1354 *base_probp = prev_prob;
1355 }
1356 } else {
1357 if (cur_prob > EXACT_TEST_BIAS) {
1358 prev_prob = tmp12 * tmp23;
1359 while (prev_prob > 0.5) {
1360 tmp13 += 1;
1361 tmp22 += 1;
1362 cur_prob *= prev_prob / (tmp13 * tmp22);
1363 tmp12 -= 1;
1364 tmp23 -= 1;
1365 if (cur_prob <= EXACT_TEST_BIAS) {
1366 break;
1367 }
1368 prev_prob = tmp12 * tmp23;
1369 }
1370 *base_probp = cur_prob;
1371 tmps12 = tmp12;
1372 tmps13 = tmp13;
1373 tmps22 = tmp22;
1374 tmps23 = tmp23;
1375 } else {
1376 tmps12 = tmp12;
1377 tmps13 = tmp13;
1378 tmps22 = tmp22;
1379 tmps23 = tmp23;
1380 while (1) {
1381 prev_prob = cur_prob;
1382 tmp12 += 1;
1383 tmp23 += 1;
1384 cur_prob *= (tmp13 * tmp22) / (tmp12 * tmp23);
1385 if (cur_prob < prev_prob) {
1386 return 1;
1387 }
1388 tmp13 -= 1;
1389 tmp22 -= 1;
1390 if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1391 if (cur_prob > EXACT_TEST_BIAS) {
1392 break;
1393 }
1394 *tie_ctp += 1;
1395 }
1396 total += cur_prob;
1397 }
1398 prev_prob = cur_prob;
1399 cur_prob = *base_probp;
1400 *base_probp = prev_prob;
1401 }
1402 }
1403 *saved12p = tmp12;
1404 *saved13p = tmp13;
1405 *saved22p = tmp22;
1406 *saved23p = tmp23;
1407 if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1408 if (cur_prob > EXACT_TEST_BIAS) {
1409 // even most extreme table on this side is too probable
1410 *totalp = 0;
1411 return 0;
1412 }
1413 *tie_ctp += 1;
1414 }
1415 // sum tail to floating point precision limit
1416 if (right_side) {
1417 prev_prob = total;
1418 total += cur_prob;
1419 while (total > prev_prob) {
1420 tmps12 += 1;
1421 tmps23 += 1;
1422 cur_prob *= (tmps13 * tmps22) / (tmps12 * tmps23);
1423 tmps13 -= 1;
1424 tmps22 -= 1;
1425 prev_prob = total;
1426 total += cur_prob;
1427 }
1428 } else {
1429 prev_prob = total;
1430 total += cur_prob;
1431 while (total > prev_prob) {
1432 tmps13 += 1;
1433 tmps22 += 1;
1434 cur_prob *= (tmps12 * tmps23) / (tmps13 * tmps22);
1435 tmps12 -= 1;
1436 tmps23 -= 1;
1437 prev_prob = total;
1438 total += cur_prob;
1439 }
1440 }
1441 *totalp = total;
1442 return 0;
1443 }
1444
fisher23(uint32_t m11,uint32_t m12,uint32_t m13,uint32_t m21,uint32_t m22,uint32_t m23,uint32_t midp)1445 double fisher23(uint32_t m11, uint32_t m12, uint32_t m13, uint32_t m21, uint32_t m22, uint32_t m23, uint32_t midp) {
1446 // 2x3 Fisher-Freeman-Halton exact test p-value calculation.
1447 // The number of tables involved here is still small enough that the network
1448 // algorithm (and the improved variants thereof that I've seen) are
1449 // suboptimal; a 2-dimensional version of the SNPHWE2 strategy has higher
1450 // performance.
1451 // 2x4, 2x5, and 3x3 should also be practical with this method, but beyond
1452 // that I doubt it's worth the trouble.
1453 // Complexity of approach is O(n^{df/2}), where n is number of observations.
1454 double cur_prob = (1 - FISHER_EPSILON) * EXACT_TEST_BIAS;
1455 double tprob = cur_prob;
1456 double cprob = 0;
1457 double dyy = 0;
1458 uint32_t tie_ct = 1;
1459 uint32_t dir = 0; // 0 = forwards, 1 = backwards
1460 double base_probl;
1461 double base_probr;
1462 double orig_base_probl;
1463 double orig_base_probr;
1464 double orig_row_prob;
1465 double row_prob;
1466 uint32_t uii;
1467 uint32_t ujj;
1468 uint32_t ukk;
1469 double cur11;
1470 double cur21;
1471 double savedl12;
1472 double savedl13;
1473 double savedl22;
1474 double savedl23;
1475 double savedr12;
1476 double savedr13;
1477 double savedr22;
1478 double savedr23;
1479 double orig_savedl12;
1480 double orig_savedl13;
1481 double orig_savedl22;
1482 double orig_savedl23;
1483 double orig_savedr12;
1484 double orig_savedr13;
1485 double orig_savedr22;
1486 double orig_savedr23;
1487 double tmp12;
1488 double tmp13;
1489 double tmp22;
1490 double tmp23;
1491 double dxx;
1492 double preaddp;
1493 // Ensure m11 + m21 <= m12 + m22 <= m13 + m23.
1494 uii = m11 + m21;
1495 ujj = m12 + m22;
1496 if (uii > ujj) {
1497 ukk = m11;
1498 m11 = m12;
1499 m12 = ukk;
1500 ukk = m21;
1501 m21 = m22;
1502 m22 = ukk;
1503 ukk = uii;
1504 uii = ujj;
1505 ujj = ukk;
1506 }
1507 ukk = m13 + m23;
1508 if (ujj > ukk) {
1509 ujj = ukk;
1510 ukk = m12;
1511 m12 = m13;
1512 m13 = ukk;
1513 ukk = m22;
1514 m22 = m23;
1515 m23 = ukk;
1516 }
1517 if (uii > ujj) {
1518 ukk = m11;
1519 m11 = m12;
1520 m12 = ukk;
1521 ukk = m21;
1522 m21 = m22;
1523 m22 = ukk;
1524 }
1525 // Ensure majority of probability mass is in front of m11.
1526 if ((((uint64_t)m11) * (m22 + m23)) > (((uint64_t)m21) * (m12 + m13))) {
1527 ukk = m11;
1528 m11 = m21;
1529 m21 = ukk;
1530 ukk = m12;
1531 m12 = m22;
1532 m22 = ukk;
1533 ukk = m13;
1534 m13 = m23;
1535 m23 = ukk;
1536 }
1537 if ((((uint64_t)m12) * m23) > (((uint64_t)m13) * m22)) {
1538 base_probr = cur_prob;
1539 savedr12 = m12;
1540 savedr13 = m13;
1541 savedr22 = m22;
1542 savedr23 = m23;
1543 tmp12 = savedr12;
1544 tmp13 = savedr13;
1545 tmp22 = savedr22;
1546 tmp23 = savedr23;
1547 // m12 and m23 must be nonzero
1548 dxx = tmp12 * tmp23;
1549 do {
1550 tmp13 += 1;
1551 tmp22 += 1;
1552 cur_prob *= dxx / (tmp13 * tmp22);
1553 tmp12 -= 1;
1554 tmp23 -= 1;
1555 if (cur_prob <= EXACT_TEST_BIAS) {
1556 if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1557 tie_ct++;
1558 }
1559 tprob += cur_prob;
1560 break;
1561 }
1562 cprob += cur_prob;
1563 if (cprob == INFINITY) {
1564 return 0;
1565 }
1566 dxx = tmp12 * tmp23;
1567 // must enforce tmp12 >= 0 and tmp23 >= 0 since we're saving these
1568 } while (dxx > 0.5);
1569 savedl12 = tmp12;
1570 savedl13 = tmp13;
1571 savedl22 = tmp22;
1572 savedl23 = tmp23;
1573 base_probl = cur_prob;
1574 do {
1575 tmp13 += 1;
1576 tmp22 += 1;
1577 cur_prob *= (tmp12 * tmp23) / (tmp13 * tmp22);
1578 tmp12 -= 1;
1579 tmp23 -= 1;
1580 preaddp = tprob;
1581 tprob += cur_prob;
1582 } while (tprob > preaddp);
1583 tmp12 = savedr12;
1584 tmp13 = savedr13;
1585 tmp22 = savedr22;
1586 tmp23 = savedr23;
1587 cur_prob = base_probr;
1588 do {
1589 tmp12 += 1;
1590 tmp23 += 1;
1591 cur_prob *= (tmp13 * tmp22) / (tmp12 * tmp23);
1592 tmp13 -= 1;
1593 tmp22 -= 1;
1594 preaddp = tprob;
1595 tprob += cur_prob;
1596 } while (tprob > preaddp);
1597 } else {
1598 base_probl = cur_prob;
1599 savedl12 = m12;
1600 savedl13 = m13;
1601 savedl22 = m22;
1602 savedl23 = m23;
1603 if (!((((uint64_t)m12) * m23) + (((uint64_t)m13) * m22))) {
1604 base_probr = cur_prob;
1605 savedr12 = savedl12;
1606 savedr13 = savedl13;
1607 savedr22 = savedl22;
1608 savedr23 = savedl23;
1609 } else {
1610 tmp12 = savedl12;
1611 tmp13 = savedl13;
1612 tmp22 = savedl22;
1613 tmp23 = savedl23;
1614 dxx = tmp13 * tmp22;
1615 do {
1616 tmp12 += 1;
1617 tmp23 += 1;
1618 cur_prob *= dxx / (tmp12 * tmp23);
1619 tmp13 -= 1;
1620 tmp22 -= 1;
1621 if (cur_prob <= EXACT_TEST_BIAS) {
1622 if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1623 tie_ct++;
1624 }
1625 tprob += cur_prob;
1626 break;
1627 }
1628 cprob += cur_prob;
1629 if (cprob == INFINITY) {
1630 return 0;
1631 }
1632 dxx = tmp13 * tmp22;
1633 } while (dxx > 0.5);
1634 savedr12 = tmp12;
1635 savedr13 = tmp13;
1636 savedr22 = tmp22;
1637 savedr23 = tmp23;
1638 base_probr = cur_prob;
1639 do {
1640 tmp12 += 1;
1641 tmp23 += 1;
1642 cur_prob *= (tmp13 * tmp22) / (tmp12 * tmp23);
1643 tmp13 -= 1;
1644 tmp22 -= 1;
1645 preaddp = tprob;
1646 tprob += cur_prob;
1647 } while (tprob > preaddp);
1648 tmp12 = savedl12;
1649 tmp13 = savedl13;
1650 tmp22 = savedl22;
1651 tmp23 = savedl23;
1652 cur_prob = base_probl;
1653 do {
1654 tmp13 += 1;
1655 tmp22 += 1;
1656 cur_prob *= (tmp12 * tmp23) / (tmp13 * tmp22);
1657 tmp12 -= 1;
1658 tmp23 -= 1;
1659 preaddp = tprob;
1660 tprob += cur_prob;
1661 } while (tprob > preaddp);
1662 }
1663 }
1664 row_prob = tprob + cprob;
1665 orig_base_probl = base_probl;
1666 orig_base_probr = base_probr;
1667 orig_row_prob = row_prob;
1668 orig_savedl12 = savedl12;
1669 orig_savedl13 = savedl13;
1670 orig_savedl22 = savedl22;
1671 orig_savedl23 = savedl23;
1672 orig_savedr12 = savedr12;
1673 orig_savedr13 = savedr13;
1674 orig_savedr22 = savedr22;
1675 orig_savedr23 = savedr23;
1676 for (; dir < 2; dir++) {
1677 cur11 = m11;
1678 cur21 = m21;
1679 if (dir) {
1680 base_probl = orig_base_probl;
1681 base_probr = orig_base_probr;
1682 row_prob = orig_row_prob;
1683 savedl12 = orig_savedl12;
1684 savedl13 = orig_savedl13;
1685 savedl22 = orig_savedl22;
1686 savedl23 = orig_savedl23;
1687 savedr12 = orig_savedr12;
1688 savedr13 = orig_savedr13;
1689 savedr22 = orig_savedr22;
1690 savedr23 = orig_savedr23;
1691 ukk = m11;
1692 if (ukk > m22 + m23) {
1693 ukk = m22 + m23;
1694 }
1695 } else {
1696 ukk = m21;
1697 if (ukk > m12 + m13) {
1698 ukk = m12 + m13;
1699 }
1700 }
1701 ukk++;
1702 while (--ukk) {
1703 if (dir) {
1704 cur21 += 1;
1705 if (savedl23) {
1706 savedl13 += 1;
1707 row_prob *= (cur11 * (savedl22 + savedl23)) / (cur21 * (savedl12 + savedl13));
1708 base_probl *= (cur11 * savedl23) / (cur21 * savedl13);
1709 savedl23 -= 1;
1710 } else {
1711 savedl12 += 1;
1712 row_prob *= (cur11 * (savedl22 + savedl23)) / (cur21 * (savedl12 + savedl13));
1713 base_probl *= (cur11 * savedl22) / (cur21 * savedl12);
1714 savedl22 -= 1;
1715 }
1716 cur11 -= 1;
1717 } else {
1718 cur11 += 1;
1719 if (savedl12) {
1720 savedl22 += 1;
1721 row_prob *= (cur21 * (savedl12 + savedl13)) / (cur11 * (savedl22 + savedl23));
1722 base_probl *= (cur21 * savedl12) / (cur11 * savedl22);
1723 savedl12 -= 1;
1724 } else {
1725 savedl23 += 1;
1726 row_prob *= (cur21 * (savedl12 + savedl13)) / (cur11 * (savedl22 + savedl23));
1727 base_probl *= (cur21 * savedl13) / (cur11 * savedl23);
1728 savedl13 -= 1;
1729 }
1730 cur21 -= 1;
1731 }
1732 if (fisher23_tailsum(&base_probl, &savedl12, &savedl13, &savedl22, &savedl23, &dxx, &tie_ct, 0)) {
1733 break;
1734 }
1735 tprob += dxx;
1736 if (dir) {
1737 if (savedr22) {
1738 savedr12 += 1;
1739 base_probr *= ((cur11 + 1) * savedr22) / (cur21 * savedr12);
1740 savedr22 -= 1;
1741 } else {
1742 savedr13 += 1;
1743 base_probr *= ((cur11 + 1) * savedr23) / (cur21 * savedr13);
1744 savedr23 -= 1;
1745 }
1746 } else {
1747 if (savedr13) {
1748 savedr23 += 1;
1749 base_probr *= ((cur21 + 1) * savedr13) / (cur11 * savedr23);
1750 savedr13 -= 1;
1751 } else {
1752 savedr22 += 1;
1753 base_probr *= ((cur21 + 1) * savedr12) / (cur11 * savedr22);
1754 savedr12 -= 1;
1755 }
1756 }
1757 fisher23_tailsum(&base_probr, &savedr12, &savedr13, &savedr22, &savedr23, &dyy, &tie_ct, 1);
1758 tprob += dyy;
1759 cprob += row_prob - dxx - dyy;
1760 if (cprob == INFINITY) {
1761 return 0;
1762 }
1763 }
1764 if (!ukk) {
1765 continue;
1766 }
1767 savedl12 += savedl13;
1768 savedl22 += savedl23;
1769 if (dir) {
1770 while (1) {
1771 preaddp = tprob;
1772 tprob += row_prob;
1773 if (tprob <= preaddp) {
1774 break;
1775 }
1776 cur21 += 1;
1777 savedl12 += 1;
1778 row_prob *= (cur11 * savedl22) / (cur21 * savedl12);
1779 cur11 -= 1;
1780 savedl22 -= 1;
1781 }
1782 } else {
1783 while (1) {
1784 preaddp = tprob;
1785 tprob += row_prob;
1786 if (tprob <= preaddp) {
1787 break;
1788 }
1789 cur11 += 1;
1790 savedl22 += 1;
1791 row_prob *= (cur21 * savedl12) / (cur11 * savedl22);
1792 cur21 -= 1;
1793 savedl12 -= 1;
1794 }
1795 }
1796 }
1797 if (!midp) {
1798 return tprob / (tprob + cprob);
1799 } else {
1800 return (tprob - ((1 - FISHER_EPSILON) * EXACT_TEST_BIAS * 0.5) * ((int32_t)tie_ct)) / (tprob + cprob);
1801 }
1802 }
1803
chi22_get_coeffs(intptr_t row1_sum,intptr_t col1_sum,intptr_t total,double * expm11p,double * recip_sump)1804 void chi22_get_coeffs(intptr_t row1_sum, intptr_t col1_sum, intptr_t total, double* expm11p, double* recip_sump) {
1805 // chisq = (m11 - expm11)^2 * recip_sum
1806 // (see discussion for chi22_precomp_val_bounds() below.)
1807 //
1808 // expm11 = row1_sum * col1_sum / total
1809 // expm12 = row1_sum * col2_sum / total, etc.
1810 // recip_sum = 1 / expm11 + 1 / expm12 + 1 / expm21 + 1 / expm22
1811 // = total * (1 / (row1_sum * col1_sum) + 1 / (row1_sum * col2_sum) +
1812 // 1 / (row2_sum * col1_sum) + 1 / (row2_sum * col2_sum))
1813 // = total^3 / (row1_sum * col1_sum * row2_sum * col2_sum)
1814 double m11_numer = ((uint64_t)row1_sum) * ((uint64_t)col1_sum);
1815 double denom = m11_numer * (((uint64_t)(total - row1_sum)) * ((uint64_t)(total - col1_sum)));
1816 double dxx;
1817 if (denom != 0) {
1818 dxx = total;
1819 *expm11p = m11_numer / dxx;
1820 *recip_sump = dxx * dxx * dxx / denom;
1821 } else {
1822 // since an entire row or column is zero, either m11 or m22 is zero
1823 // row1_sum + col1_sum - total = m11 - m22
1824 if (row1_sum + col1_sum < total) {
1825 *expm11p = 0;
1826 } else {
1827 *expm11p = row1_sum + col1_sum - total;
1828 }
1829 *recip_sump = 0;
1830 }
1831 }
1832
chi22_eval(intptr_t m11,intptr_t row1_sum,intptr_t col1_sum,intptr_t total)1833 double chi22_eval(intptr_t m11, intptr_t row1_sum, intptr_t col1_sum, intptr_t total) {
1834 double expm11_numer = ((uint64_t)row1_sum) * ((uint64_t)col1_sum);
1835 double denom = expm11_numer * (((uint64_t)(total - row1_sum)) * ((uint64_t)(total - col1_sum)));
1836 double dxx;
1837 double dyy;
1838 if (denom != 0) {
1839 dxx = total;
1840 dyy = m11 * dxx - expm11_numer; // total * (m11 - expm11)
1841 return (dyy * dyy * dxx) / denom;
1842 } else {
1843 return 0;
1844 }
1845 }
1846
chi22_evalx(intptr_t m11,intptr_t row1_sum,intptr_t col1_sum,intptr_t total)1847 double chi22_evalx(intptr_t m11, intptr_t row1_sum, intptr_t col1_sum, intptr_t total) {
1848 // PLINK emulation. returns -9 instead of 0 if row1_sum, row2_sum, col1_sum,
1849 // or col2_sum is zero, for identical "NA" reporting.
1850 double expm11_numer = ((uint64_t)row1_sum) * ((uint64_t)col1_sum);
1851 double denom = expm11_numer * (((uint64_t)(total - row1_sum)) * ((uint64_t)(total - col1_sum)));
1852 double dxx;
1853 double dyy;
1854 if (denom != 0) {
1855 dxx = total;
1856 dyy = m11 * dxx - expm11_numer; // total * (m11 - expm11)
1857 return (dyy * dyy * dxx) / denom;
1858 } else {
1859 return -9;
1860 }
1861 }
1862
chi22_precomp_val_bounds(double chisq,intptr_t row1_sum,intptr_t col1_sum,intptr_t total,uint32_t * bounds,double * coeffs)1863 void chi22_precomp_val_bounds(double chisq, intptr_t row1_sum, intptr_t col1_sum, intptr_t total, uint32_t* bounds, double* coeffs) {
1864 // Treating m11 as the only variable, this returns the minimum and (maximum +
1865 // 1) values of m11 which produce smaller chisq statistics than given in
1866 // bounds[0] and bounds[1] respectively, and smaller-or-equal interval
1867 // bounds in bounds[2] and bounds[3].
1868 double expm11;
1869 double recip_sum;
1870 double cur11;
1871 double dxx;
1872 intptr_t ceil11;
1873 intptr_t lii;
1874 chi22_get_coeffs(row1_sum, col1_sum, total, &expm11, &recip_sum);
1875 if (coeffs) {
1876 coeffs[0] = expm11;
1877 coeffs[1] = recip_sum;
1878 }
1879 if (recip_sum == 0) {
1880 // sum-0 row or column, no freedom at all
1881 bounds[0] = (intptr_t)expm11;
1882 bounds[1] = bounds[0];
1883 bounds[2] = bounds[0];
1884 if (chisq == 0) {
1885 bounds[3] = bounds[0] + 1;
1886 } else {
1887 bounds[3] = bounds[0];
1888 }
1889 return;
1890 }
1891
1892 // double cur_stat = (cur11 - exp11) * (cur11 - exp11) * recipx11 + (cur12 - exp12) * (cur12 - exp12) * recipx12 + (cur21 - exp21) * (cur21 - exp21) * recipx21 + (cur22 - exp22) * (cur22 - exp22) * recipx22;
1893 // However, we have
1894 // cur11 - exp11 = -(cur12 - exp12) = -(cur21 - exp21) = cur22 - exp22
1895 // So the chisq statistic reduces to
1896 // (cur11 - exp11)^2 * (recipx11 + recipx12 + recipx21 + recipx22).
1897
1898 ceil11 = MINV(row1_sum, col1_sum);
1899 // chisq = (cur11 - expm11)^2 * recip_sum
1900 // -> expm11 +/- sqrt(chisq / recip_sum) = cur11
1901 recip_sum = sqrt(chisq / recip_sum);
1902 cur11 = expm11 - recip_sum;
1903 dxx = cur11 + 1 - BIG_EPSILON;
1904 if (dxx < 0) {
1905 bounds[0] = 0;
1906 bounds[2] = 0;
1907 } else {
1908 lii = (intptr_t)dxx;
1909 bounds[2] = lii;
1910 if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
1911 bounds[0] = lii + 1;
1912 } else {
1913 bounds[0] = lii;
1914 }
1915 }
1916 cur11 = expm11 + recip_sum;
1917 if (cur11 > ceil11 + BIG_EPSILON) {
1918 bounds[1] = ceil11 + 1;
1919 bounds[3] = bounds[1];
1920 } else {
1921 dxx = cur11 + 1 - BIG_EPSILON;
1922 lii = (intptr_t)dxx;
1923 bounds[1] = lii;
1924 if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
1925 bounds[3] = lii + 1;
1926 } else {
1927 bounds[3] = lii;
1928 }
1929 }
1930 }
1931
chi23_eval(intptr_t m11,intptr_t m12,intptr_t row1_sum,intptr_t col1_sum,intptr_t col2_sum,intptr_t total)1932 double chi23_eval(intptr_t m11, intptr_t m12, intptr_t row1_sum, intptr_t col1_sum, intptr_t col2_sum, intptr_t total) {
1933 // assumes no sum-zero row
1934 intptr_t m13 = row1_sum - m11 - m12;
1935 intptr_t col3_sum = total - col1_sum - col2_sum;
1936 double col1_sumd;
1937 double col2_sumd;
1938 double col3_sumd;
1939 double tot_recip;
1940 double dxx;
1941 double expect;
1942 double delta;
1943 double chisq;
1944 col1_sumd = col1_sum;
1945 col2_sumd = col2_sum;
1946 col3_sumd = col3_sum;
1947 tot_recip = 1.0 / ((double)total);
1948 dxx = row1_sum * tot_recip;
1949 expect = dxx * col1_sumd;
1950 delta = m11 - expect;
1951 chisq = delta * delta / expect;
1952 expect = dxx * col2_sumd;
1953 delta = m12 - expect;
1954 chisq += delta * delta / expect;
1955 expect = dxx * col3_sumd;
1956 delta = m13 - expect;
1957 chisq += delta * delta / expect;
1958 dxx = (total - row1_sum) * tot_recip;
1959 expect = dxx * col1_sumd;
1960 delta = (col1_sum - m11) - expect;
1961 chisq += delta * delta / expect;
1962 expect = dxx * col2_sumd;
1963 delta = (col2_sum - m12) - expect;
1964 chisq += delta * delta / expect;
1965 expect = dxx * col3_sumd;
1966 delta = (col3_sum - m13) - expect;
1967 chisq += delta * delta / expect;
1968 if (chisq < (SMALL_EPSILON * SMALL_EPSILON)) {
1969 return 0;
1970 }
1971 return chisq;
1972 }
1973
chi23_evalx(intptr_t m11,intptr_t m12,intptr_t m13,intptr_t m21,intptr_t m22,intptr_t m23,double * chip,uint32_t * dfp)1974 void chi23_evalx(intptr_t m11, intptr_t m12, intptr_t m13, intptr_t m21, intptr_t m22, intptr_t m23, double* chip, uint32_t* dfp) {
1975 // Slightly different from PLINK calculation, since it detects lone nonzero
1976 // columns.
1977 intptr_t row1_sum = m11 + m12 + m13;
1978 intptr_t row2_sum = m21 + m22 + m23;
1979 intptr_t col1_sum = m11 + m21;
1980 intptr_t col2_sum = m12 + m22;
1981 intptr_t col3_sum = m13 + m23;
1982 intptr_t total;
1983 double col1_sumd;
1984 double col2_sumd;
1985 double col3_sumd;
1986 double tot_recip;
1987 double dxx;
1988 double expect;
1989 double delta;
1990 double chisq;
1991 if ((!row1_sum) || (!row2_sum)) {
1992 *chip = -9;
1993 *dfp = 0;
1994 return;
1995 }
1996 total = row1_sum + row2_sum;
1997 if (!col1_sum) {
1998 *chip = chi22_evalx(m12, row1_sum, col2_sum, total);
1999 if (*chip != -9) {
2000 *dfp = 1;
2001 } else {
2002 *dfp = 0;
2003 }
2004 return;
2005 } else if ((!col2_sum) || (!col3_sum)) {
2006 *chip = chi22_evalx(m11, row1_sum, col1_sum, total);
2007 if (*chip != -9) {
2008 *dfp = 1;
2009 } else {
2010 *dfp = 0;
2011 }
2012 return;
2013 }
2014 col1_sumd = col1_sum;
2015 col2_sumd = col2_sum;
2016 col3_sumd = col3_sum;
2017 tot_recip = 1.0 / ((double)total);
2018 dxx = row1_sum * tot_recip;
2019 expect = dxx * col1_sumd;
2020 delta = m11 - expect;
2021 chisq = delta * delta / expect;
2022 expect = dxx * col2_sumd;
2023 delta = m12 - expect;
2024 chisq += delta * delta / expect;
2025 expect = dxx * col3_sumd;
2026 delta = m13 - expect;
2027 chisq += delta * delta / expect;
2028 dxx = row2_sum * tot_recip;
2029 expect = dxx * col1_sumd;
2030 delta = m21 - expect;
2031 chisq += delta * delta / expect;
2032 expect = dxx * col2_sumd;
2033 delta = m22 - expect;
2034 chisq += delta * delta / expect;
2035 expect = dxx * col3_sumd;
2036 delta = m23 - expect;
2037 chisq += delta * delta / expect;
2038 if (chisq < (SMALL_EPSILON * SMALL_EPSILON)) {
2039 chisq = 0;
2040 }
2041 *chip = chisq;
2042 *dfp = 2;
2043 }
2044
ca_trend_eval(intptr_t case_dom_ct,intptr_t case_ct,intptr_t het_ct,intptr_t homdom_ct,intptr_t total)2045 double ca_trend_eval(intptr_t case_dom_ct, intptr_t case_ct, intptr_t het_ct, intptr_t homdom_ct, intptr_t total) {
2046 // case_dom_ct is an allele count (2 * homa2 + het), while other inputs are
2047 // observation counts.
2048 //
2049 // If case_missing_ct is fixed,
2050 // row1_sum = case ct
2051 // col1_sum = A2 ct
2052 // case_ct * ctrl_ct * REC_ct * DOM_ct
2053 // REC_ct = 2 * obs_11 + obs_12
2054 // DOM_ct = 2 * obs_22 + obs_12
2055 // CA = (obs_U / obs_T) * (case REC ct) - (obs_A / obs_T) * (ctrl DOM ct)
2056 // = (case A2) * (obs_U / obs_T) - (obs_A / obs_T) * (DOM ct - case DOM)
2057 // = (case A2) * (obs_U / obs_T) + (case DOM) * (obs_A / obs_T) - DOM*(A/T)
2058 // = (case A2 ct) - total A2 ct * (A/T)
2059 // CAT = CA * obs_T
2060 // varCA_recip = obs_T * obs_T * obs_T /
2061 // (obs_A * obs_U * (obs_T * (obs_12 + 4 * obs_22) - DOMct * DOMct))
2062 // trend statistic = CAT * CAT * [varCA_recip / obs_T^2]
2063 double dom_ct = het_ct + 2 * homdom_ct;
2064 double totald = total;
2065 double case_ctd = case_ct;
2066 double cat = case_dom_ct * totald - dom_ct * case_ctd;
2067 double dxx = totald * (het_ct + 4 * ((int64_t)homdom_ct)) - dom_ct * dom_ct;
2068
2069 // This should never be called with dxx == 0 (which happens when two columns
2070 // are all-zero). Use ca_trend_evalx() to check for that.
2071 dxx *= case_ctd * (totald - case_ctd);
2072 return cat * cat * totald / dxx;
2073 }
2074
ca_trend_evalx(intptr_t case_dom_ct,intptr_t case_ct,intptr_t het_ct,intptr_t homdom_ct,intptr_t total)2075 double ca_trend_evalx(intptr_t case_dom_ct, intptr_t case_ct, intptr_t het_ct, intptr_t homdom_ct, intptr_t total) {
2076 double dom_ct = het_ct + 2 * homdom_ct;
2077 double totald = total;
2078 double case_ctd = case_ct;
2079 double cat = case_dom_ct * totald - dom_ct * case_ctd;
2080 double dxx = totald * (het_ct + 4 * ((int64_t)homdom_ct)) - dom_ct * dom_ct;
2081 if (dxx != 0) {
2082 dxx *= case_ctd * (totald - case_ctd);
2083 return cat * cat * totald / dxx;
2084 } else {
2085 return -9;
2086 }
2087 }
2088
ca_trend_precomp_val_bounds(double chisq,intptr_t case_ct,intptr_t het_ct,intptr_t homdom_ct,intptr_t total,uint32_t * bounds,double * coeffs)2089 void ca_trend_precomp_val_bounds(double chisq, intptr_t case_ct, intptr_t het_ct, intptr_t homdom_ct, intptr_t total, uint32_t* bounds, double* coeffs) {
2090 // If case_missing_ct is fixed,
2091 // row1_sum = case ct
2092 // col1_sum = DOM ct
2093 // case_ct * ctrl_ct * REC_ct * DOM_ct
2094 // REC_ct = 2 * obs_11 + obs_12
2095 // DOM_ct = 2 * obs_22 + obs_12
2096 // CA = (obs_U / obs_T) * (case DOM ct) - (obs_A / obs_T) * (ctrl DOM ct)
2097 // = (case DOM) * (obs_U / obs_T) - (obs_A / obs_T) * (DOM ct - case DOM)
2098 // = (case DOM) * (obs_U / obs_T) + (case DOM) * (obs_A / obs_T) - DOM*(A/T)
2099 // = (case DOM ct) - total DOM ct * (A/T)
2100 // varCA_recip = obs_T * obs_T * obs_T /
2101 // (obs_A * obs_U * (obs_T * (obs_12 + 4 * obs_22) - DOMct * DOMct))
2102 // trend statistic = CA * CA * varCA_recip
2103 intptr_t dom_ct = het_ct + 2 * homdom_ct;
2104 double dom_ctd = dom_ct;
2105 double totald = total;
2106 double case_ctd = case_ct;
2107 double tot_recip = 1.0 / totald;
2108 double expm11 = dom_ctd * case_ctd * tot_recip;
2109 double dxx = case_ctd * (totald - case_ctd) * (totald * (het_ct + 4 * ((int64_t)homdom_ct)) - dom_ctd * dom_ctd);
2110 double varca_recip;
2111 double cur11;
2112 intptr_t ceil11;
2113 intptr_t lii;
2114 if (dxx == 0) {
2115 // bounds/coefficients should never be referenced in this case
2116 return;
2117 }
2118 varca_recip = totald * totald * totald / dxx;
2119 if (coeffs) {
2120 coeffs[0] = expm11;
2121 coeffs[1] = varca_recip;
2122 }
2123
2124 // statistic: (cur11 - expm11)^2 * varca_recip
2125 ceil11 = case_ct * 2;
2126 if (dom_ct < ceil11) {
2127 ceil11 = dom_ct;
2128 }
2129 // chisq = (cur11 - expm11)^2 * varca_recip
2130 // -> expm11 +/- sqrt(chisq / varca_recip) = cur11
2131 varca_recip = sqrt(chisq / varca_recip);
2132 cur11 = expm11 - varca_recip;
2133 dxx = cur11 + 1 - BIG_EPSILON;
2134 if (dxx < 0) {
2135 bounds[0] = 0;
2136 bounds[2] = 0;
2137 } else {
2138 lii = (intptr_t)dxx;
2139 bounds[2] = lii;
2140 if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
2141 bounds[0] = lii + 1;
2142 } else {
2143 bounds[0] = lii;
2144 }
2145 }
2146 cur11 = expm11 + varca_recip;
2147 if (cur11 > ceil11 + BIG_EPSILON) {
2148 bounds[1] = ceil11 + 1;
2149 bounds[3] = bounds[1];
2150 } else {
2151 dxx = cur11 + 1 - BIG_EPSILON;
2152 lii = (intptr_t)dxx;
2153 bounds[1] = lii;
2154 if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
2155 bounds[3] = lii + 1;
2156 } else {
2157 bounds[3] = lii;
2158 }
2159 }
2160 }
2161
linear_hypothesis_chisq(uintptr_t constraint_ct,uintptr_t param_ct,double * constraints_con_major,double * coef,double * cov_matrix,double * param_df_buf,double * param_df_buf2,double * df_df_buf,MATRIX_INVERT_BUF1_TYPE * mi_buf,double * df_buf,double * chisq_ptr)2162 uint32_t linear_hypothesis_chisq(uintptr_t constraint_ct, uintptr_t param_ct, double* constraints_con_major, double* coef, double* cov_matrix, double* param_df_buf, double* param_df_buf2, double* df_df_buf, MATRIX_INVERT_BUF1_TYPE* mi_buf, double* df_buf, double* chisq_ptr) {
2163 // See PLINK model.cpp Model::linearHypothesis().
2164 //
2165 // outer = df_buf
2166 // inner = df_df_buf
2167 // tmp = param_df_buf
2168 // mi_buf only needs to be of length constraint_ct
2169 //
2170 // Since no PLINK function ever calls this with nonzero h[] values, this just
2171 // takes a df parameter for now; it's trivial to switch to the more general
2172 // interface later.
2173 double* dptr = constraints_con_major;
2174 double* dptr2;
2175 uintptr_t constraint_idx;
2176 uintptr_t constraint_idx2;
2177 uintptr_t param_idx;
2178 double dxx;
2179 double dyy;
2180 for (constraint_idx = 0; constraint_idx < constraint_ct; constraint_idx++) {
2181 dxx = 0;
2182 dptr2 = coef;
2183 for (param_idx = 0; param_idx < param_ct; param_idx++) {
2184 dxx += (*dptr++) * (*dptr2++);
2185 }
2186 df_buf[constraint_idx] = dxx;
2187 }
2188 // temporarily set param_df_buf2[][] to H-transpose
2189 transpose_copy(constraint_ct, param_ct, constraints_con_major, param_df_buf2);
2190 col_major_matrix_multiply(constraint_ct, param_ct, param_ct, param_df_buf2, cov_matrix, param_df_buf);
2191 // tmp[][] is now param-major
2192 col_major_matrix_multiply(constraint_ct, constraint_ct, param_ct, param_df_buf, constraints_con_major, df_df_buf);
2193
2194 if (invert_matrix((uint32_t)constraint_ct, df_df_buf, mi_buf, param_df_buf2)) {
2195 return 1;
2196 }
2197 dxx = 0; // result
2198 dptr = df_df_buf;
2199 for (constraint_idx = 0; constraint_idx < constraint_ct; constraint_idx++) {
2200 dyy = 0; // tmp2[c]
2201 dptr2 = df_buf;
2202 for (constraint_idx2 = 0; constraint_idx2 < constraint_ct; constraint_idx2++) {
2203 dyy += (*dptr++) * (*dptr2++);
2204 }
2205 dxx += dyy * df_buf[constraint_idx];
2206 }
2207 *chisq_ptr = dxx;
2208 return 0;
2209 }
2210
binom_2sided(uint32_t succ,uint32_t obs,uint32_t midp)2211 double binom_2sided(uint32_t succ, uint32_t obs, uint32_t midp) {
2212 // straightforward to generalize this to any success probability
2213 double cur_succ_t2 = (int32_t)succ;
2214 double cur_fail_t2 = (int32_t)(obs - succ);
2215 double tailp = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
2216 double centerp = 0;
2217 double lastp2 = tailp;
2218 double lastp1 = tailp;
2219 int32_t tie_ct = 1;
2220 // double rate_mult_incr = rate / (1 - rate);
2221 // double rate_mult_decr = (1 - rate) / rate;
2222 double cur_succ_t1;
2223 double cur_fail_t1;
2224 double preaddp;
2225 if (!obs) {
2226 return midp? 0.5 : 1;
2227 }
2228 if (obs < succ * 2) {
2229 // if (obs * rate < succ) {
2230 while (cur_succ_t2 > 0.5) {
2231 cur_fail_t2 += 1;
2232 lastp2 *= cur_succ_t2 / cur_fail_t2;
2233 // lastp2 *= rate_mult_decr * cur_succ_t2 / cur_fail_t2;
2234 cur_succ_t2 -= 1;
2235 if (lastp2 < EXACT_TEST_BIAS) {
2236 if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
2237 tie_ct++;
2238 }
2239 tailp += lastp2;
2240 break;
2241 }
2242 centerp += lastp2;
2243 if (centerp == INFINITY) {
2244 return 0;
2245 }
2246 }
2247 if ((centerp == 0) && (!midp)) {
2248 return 1;
2249 }
2250 while (cur_succ_t2 > 0.5) {
2251 cur_fail_t2 += 1;
2252 lastp2 *= cur_succ_t2 / cur_fail_t2;
2253 // lastp2 *= rate_mult_decr * cur_succ_t2 / cur_fail_t2;
2254 cur_succ_t2 -= 1;
2255 preaddp = tailp;
2256 tailp += lastp2;
2257 if (tailp <= preaddp) {
2258 break;
2259 }
2260 }
2261 cur_succ_t1 = (int32_t)(succ + 1);
2262 cur_fail_t1 = (int32_t)(obs - succ);
2263 while (cur_fail_t1 > 0.5) {
2264 lastp1 *= cur_fail_t1 / cur_succ_t1;
2265 // lastp1 *= rate_mult_incr * cur_fail_t1 / cur_succ_t1;
2266 preaddp = tailp;
2267 tailp += lastp1;
2268 if (tailp <= preaddp) {
2269 break;
2270 }
2271 cur_succ_t1 += 1;
2272 cur_fail_t1 -= 1;
2273 }
2274 } else {
2275 while (cur_fail_t2 > 0.5) {
2276 cur_succ_t2++;
2277 lastp2 *= cur_fail_t2 / cur_succ_t2;
2278 // lastp2 *= rate_mult_incr * cur_fail_t2 / cur_succ_t2;
2279 cur_fail_t2--;
2280 if (lastp2 < EXACT_TEST_BIAS) {
2281 if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
2282 tie_ct++;
2283 }
2284 tailp += lastp2;
2285 break;
2286 }
2287 centerp += lastp2;
2288 if (centerp == INFINITY) {
2289 return 0;
2290 }
2291 }
2292 if ((centerp == 0) && (!midp)) {
2293 return 1;
2294 }
2295 while (cur_fail_t2 > 0.5) {
2296 cur_succ_t2 += 1;
2297 lastp2 *= cur_fail_t2 / cur_succ_t2;
2298 // lastp2 *= rate_mult_incr * cur_fail_t2 / cur_succ_t2;
2299 cur_fail_t2 -= 1;
2300 preaddp = tailp;
2301 tailp += lastp2;
2302 if (tailp <= preaddp) {
2303 break;
2304 }
2305 }
2306 cur_succ_t1 = (int32_t)succ;
2307 cur_fail_t1 = (int32_t)(obs - succ);
2308 while (cur_succ_t1 > 0.5) {
2309 cur_fail_t1 += 1;
2310 lastp1 *= cur_succ_t1 / cur_fail_t1;
2311 // lastp1 *= rate_mult_decr * cur_succ_t1 / cur_fail_t1;
2312 preaddp = tailp;
2313 tailp += lastp1;
2314 if (tailp <= preaddp) {
2315 break;
2316 }
2317 cur_succ_t1 -= 1;
2318 }
2319 }
2320 if (!midp) {
2321 return tailp / (tailp + centerp);
2322 } else {
2323 return (tailp - ((1 - SMALL_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (tailp + centerp);
2324 }
2325 }
2326