1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
# Authors: Josef Pktd and example from H Raja and rewrite from Vincent Davis
# Alexandre Gramfort <alexandre.gramfort@inria.fr>
#
# Code borrowed from statsmodels
#
# License: BSD-3-Clause
import numpy as np
def _ecdf(x):
"""No frills empirical cdf used in fdrcorrection."""
nobs = len(x)
return np.arange(1, nobs + 1) / float(nobs)
def fdr_correction(pvals, alpha=0.05, method='indep'):
"""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 : array_like
Set of p-values of the individual tests.
alpha : float
Error rate.
method : 'indep' | 'negcorr'
If 'indep' it implements Benjamini/Hochberg for independent or if
'negcorr' it corresponds to Benjamini/Yekutieli.
Returns
-------
reject : array, bool
True if a hypothesis is rejected, False if not.
pval_corrected : array
P-values adjusted for multiple hypothesis testing to limit FDR.
References
----------
.. footbibliography::
"""
pvals = np.asarray(pvals)
shape_init = pvals.shape
pvals = pvals.ravel()
pvals_sortind = np.argsort(pvals)
pvals_sorted = pvals[pvals_sortind]
sortrevind = pvals_sortind.argsort()
if method in ['i', 'indep', 'p', 'poscorr']:
ecdffactor = _ecdf(pvals_sorted)
elif method in ['n', 'negcorr']:
cm = np.sum(1. / np.arange(1, len(pvals_sorted) + 1))
ecdffactor = _ecdf(pvals_sorted) / cm
else:
raise ValueError("Method should be 'indep' and 'negcorr'")
reject = pvals_sorted < (ecdffactor * alpha)
if reject.any():
rejectmax = max(np.nonzero(reject)[0])
else:
rejectmax = 0
reject[:rejectmax] = True
pvals_corrected_raw = pvals_sorted / ecdffactor
pvals_corrected = np.minimum.accumulate(pvals_corrected_raw[::-1])[::-1]
pvals_corrected[pvals_corrected > 1.0] = 1.0
pvals_corrected = pvals_corrected[sortrevind].reshape(shape_init)
reject = reject[sortrevind].reshape(shape_init)
return reject, pvals_corrected
def bonferroni_correction(pval, alpha=0.05):
"""P-value correction with Bonferroni method.
Parameters
----------
pval : array_like
Set of p-values of the individual tests.
alpha : float
Error rate.
Returns
-------
reject : array, bool
True if a hypothesis is rejected, False if not.
pval_corrected : array
P-values adjusted for multiple hypothesis testing to limit FDR.
"""
pval = np.asarray(pval)
pval_corrected = pval * float(pval.size)
# p-values must not be larger than 1.
pval_corrected = pval_corrected.clip(max=1.)
reject = pval_corrected < alpha
return reject, pval_corrected
|
