P-value correction with False Discovery Rate (FDR). Correction for multiple comparison using FDR :footcite:`GenoveseEtAl2002`. This covers Benjamini/Hochberg for independent or positively correlated and Benjamini/Yekutieli for general or negatively correlated tests. Parameters
(pvals, alpha=0.05, method="indep")
| 12 | |
| 13 | |
| 14 | def fdr_correction(pvals, alpha=0.05, method="indep"): |
| 15 | """P-value correction with False Discovery Rate (FDR). |
| 16 | |
| 17 | Correction for multiple comparison using FDR :footcite:`GenoveseEtAl2002`. |
| 18 | |
| 19 | This covers Benjamini/Hochberg for independent or positively correlated and |
| 20 | Benjamini/Yekutieli for general or negatively correlated tests. |
| 21 | |
| 22 | Parameters |
| 23 | ---------- |
| 24 | pvals : array_like |
| 25 | Set of p-values of the individual tests. |
| 26 | alpha : float |
| 27 | Error rate. |
| 28 | method : 'indep' | 'negcorr' |
| 29 | If 'indep' it implements Benjamini/Hochberg for independent or if |
| 30 | 'negcorr' it corresponds to Benjamini/Yekutieli. |
| 31 | |
| 32 | Returns |
| 33 | ------- |
| 34 | reject : array, bool |
| 35 | True if a hypothesis is rejected, False if not. |
| 36 | pval_corrected : array |
| 37 | P-values adjusted for multiple hypothesis testing to limit FDR. |
| 38 | |
| 39 | References |
| 40 | ---------- |
| 41 | .. footbibliography:: |
| 42 | """ |
| 43 | pvals = np.asarray(pvals) |
| 44 | shape_init = pvals.shape |
| 45 | pvals = pvals.ravel() |
| 46 | |
| 47 | pvals_sortind = np.argsort(pvals) |
| 48 | pvals_sorted = pvals[pvals_sortind] |
| 49 | sortrevind = pvals_sortind.argsort() |
| 50 | |
| 51 | if method in ["i", "indep", "p", "poscorr"]: |
| 52 | ecdffactor = _ecdf(pvals_sorted) |
| 53 | elif method in ["n", "negcorr"]: |
| 54 | cm = np.sum(1.0 / np.arange(1, len(pvals_sorted) + 1)) |
| 55 | ecdffactor = _ecdf(pvals_sorted) / cm |
| 56 | else: |
| 57 | raise ValueError("Method should be 'indep' and 'negcorr'") |
| 58 | |
| 59 | reject = pvals_sorted < (ecdffactor * alpha) |
| 60 | if reject.any(): |
| 61 | rejectmax = max(np.nonzero(reject)[0]) |
| 62 | else: |
| 63 | rejectmax = 0 |
| 64 | reject[:rejectmax] = True |
| 65 | |
| 66 | pvals_corrected_raw = pvals_sorted / ecdffactor |
| 67 | pvals_corrected = np.minimum.accumulate(pvals_corrected_raw[::-1])[::-1] |
| 68 | pvals_corrected[pvals_corrected > 1.0] = 1.0 |
| 69 | pvals_corrected = pvals_corrected[sortrevind].reshape(shape_init) |
| 70 | reject = reject[sortrevind].reshape(shape_init) |
| 71 | return reject, pvals_corrected |