-
-
Notifications
You must be signed in to change notification settings - Fork 137
/
multicomp.py
494 lines (390 loc) · 17 KB
/
multicomp.py
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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
# Author: Raphael Vallat <raphaelvallat9@gmail.com>
# Date: April 2018
import numpy as np
__all__ = ["multicomp"]
##############################################################################
# INTERNAL FUNCTIONS
##############################################################################
def fdr(pvals, alpha=0.05, method='fdr_bh'):
"""P-values FDR correction with Benjamini/Hochberg and
Benjamini/Yekutieli procedure.
This covers Benjamini/Hochberg for independent or positively correlated and
Benjamini/Yekutieli for general or negatively correlated tests.
Parameters
----------
pvals : array_like
Array of p-values of the individual tests.
alpha : float
Error rate (= alpha level).
method : str
FDR correction methods ::
'fdr_bh' : Benjamini/Hochberg for independent / posit correlated tests
'fdr_by' : Benjamini/Yekutieli for negatively correlated tests
Returns
-------
reject : array, bool
True if a hypothesis is rejected, False if not
pval_corrected : array
P-values adjusted for multiple hypothesis testing using the BH or BY
correction.
See also
--------
bonf : Bonferroni correction
holm : Holm-Bonferroni correction
Notes
-----
From Wikipedia:
The **Benjamini–Hochberg** procedure (BH step-up procedure) controls the
false discovery rate (FDR) at level :math:`\\alpha`. It works as follows:
1. For a given :math:`\\alpha`, find the largest :math:`k` such that
:math:`P_{(k)}\\leq \\frac {k}{m}\\alpha.`
2. Reject the null hypothesis (i.e., declare discoveries) for all
:math:`H_{(i)}` for :math:`i = 1, \\ldots, k`.
The BH procedure is valid when the m tests are independent, and also in
various scenarios of dependence, but is not universally valid.
The **Benjamini–Yekutieli** procedure (BY) controls the FDR under arbitrary
dependence assumptions. This refinement modifies the threshold and finds
the largest :math:`k` such that:
.. math::
P_{(k)} \\leq \\frac{k}{m \\cdot c(m)} \\alpha
References
----------
- Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery
rate: a practical and powerful approach to multiple testing. Journal of
the Royal Statistical Society Series B, 57, 289–300.
- Benjamini, Y., and Yekutieli, D. (2001). The control of the false
discovery rate in multiple testing under dependency. Annals of
Statistics, 29, 1165–1188.
- https://en.wikipedia.org/wiki/False_discovery_rate
Examples
--------
FDR correction of an array of p-values
>>> import pingouin as pg
>>> pvals = [.50, .003, .32, .054, .0003]
>>> reject, pvals_corr = pg.multicomp(pvals, method='fdr_bh', alpha=.05)
>>> print(reject, pvals_corr)
[False True False False True] [0.5 0.0075 0.4 0.09 0.0015]
"""
assert method.lower() in ['fdr_bh', 'fdr_by']
# Convert to array and save original shape
pvals = np.asarray(pvals)
shape_init = pvals.shape
pvals = pvals.ravel()
num_nan = np.isnan(pvals).sum()
# Sort the (flattened) p-values
pvals_sortind = np.argsort(pvals)
pvals_sorted = pvals[pvals_sortind]
sortrevind = pvals_sortind.argsort()
ntests = pvals.size - num_nan
# Empirical CDF factor
ecdffactor = np.arange(1, ntests + 1) / float(ntests)
if method.lower() == 'fdr_by':
cm = np.sum(1. / np.arange(1, ntests + 1))
ecdffactor /= cm
# Now we adjust the p-values
pvals_corr = np.diag(pvals_sorted / ecdffactor[..., None])
pvals_corr = np.minimum.accumulate(pvals_corr[::-1])[::-1]
pvals_corr = np.clip(pvals_corr, None, 1)
# And revert to the original shape and order
pvals_corr = np.append(pvals_corr, np.full(num_nan, np.nan))
pvals_corrected = pvals_corr[sortrevind].reshape(shape_init)
with np.errstate(invalid='ignore'):
reject = np.less(pvals_corrected, alpha)
# reject = reject[sortrevind].reshape(shape_init)
return reject, pvals_corrected
def bonf(pvals, alpha=0.05):
"""P-values correction with Bonferroni method.
Parameters
----------
pvals : array_like
Array of p-values of the individual tests.
alpha : float
Error rate (= alpha level).
Returns
-------
reject : array, bool
True if a hypothesis is rejected, False if not
pval_corrected : array
P-values adjusted for multiple hypothesis testing using the Bonferroni
procedure (= multiplied by the number of tests).
See also
--------
holm : Holm-Bonferroni correction
fdr : Benjamini/Hochberg and Benjamini/Yekutieli FDR correction
Notes
-----
From Wikipedia:
Statistical hypothesis testing is based on rejecting the null hypothesis
if the likelihood of the observed data under the null hypotheses is low.
If multiple hypotheses are tested, the chance of a rare event increases,
and therefore, the likelihood of incorrectly rejecting a null hypothesis
(i.e., making a Type I error) increases.
The Bonferroni correction compensates for that increase by testing each
individual hypothesis :math:`p_i` at a significance level of
:math:`p_i = \\alpha / n` where :math:`\\alpha` is the desired overall
alpha level and :math:`n` is the number of hypotheses. For example, if a
trial is testing :math:`n=20` hypotheses with a desired
:math:`\\alpha=0.05`, then the Bonferroni correction would test each
individual hypothesis at :math:`\\alpha=0.05/20=0.0025``.
The Bonferroni adjusted p-values are defined as:
.. math::
\\widetilde {p}_{{(i)}}= n \\cdot p_{{(i)}}
The Bonferroni correction tends to be a bit too conservative.
Note that NaN values are not taken into account in the p-values correction.
References
----------
- Bonferroni, C. E. (1935). Il calcolo delle assicurazioni su gruppi
di teste. Studi in onore del professore salvatore ortu carboni, 13-60.
- https://en.wikipedia.org/wiki/Bonferroni_correction
Examples
--------
>>> import pingouin as pg
>>> pvals = [.50, .003, .32, .054, .0003]
>>> reject, pvals_corr = pg.multicomp(pvals, method='bonf', alpha=.05)
>>> print(reject, pvals_corr)
[False True False False True] [1. 0.015 1. 0.27 0.0015]
"""
pvals = np.asarray(pvals)
num_nan = np.isnan(pvals).sum()
pvals_corrected = pvals * (float(pvals.size) - num_nan)
pvals_corrected = np.clip(pvals_corrected, None, 1)
with np.errstate(invalid='ignore'):
reject = np.less(pvals_corrected, alpha)
return reject, pvals_corrected
def holm(pvals, alpha=.05):
"""P-values correction with Holm method.
Parameters
----------
pvals : array_like
Array of p-values of the individual tests.
alpha : float
Error rate (= alpha level).
Returns
-------
reject : array, bool
True if a hypothesis is rejected, False if not
pvals_corrected : array
P-values adjusted for multiple hypothesis testing using the Holm
procedure.
See also
--------
bonf : Bonferroni correction
fdr : Benjamini/Hochberg and Benjamini/Yekutieli FDR correction
Notes
-----
From Wikipedia:
In statistics, the Holm–Bonferroni method (also called the Holm method) is
used to counteract the problem of multiple comparisons. It is intended to
control the family-wise error rate and offers a simple test uniformly more
powerful than the Bonferroni correction.
The Holm adjusted p-values are the running maximum of the sorted p-values
divided by the corresponding increasing alpha level:
.. math::
\\frac{\\alpha}{n}, \\frac{\\alpha}{n-1}, ..., \\frac{\\alpha}{1}
where :math:`n` is the number of test.
The full mathematical formula is:
.. math::
\\widetilde {p}_{{(i)}}=\\max _{{j\\leq i}}\\left\\{(n-j+1)p_{{(j)}}
\\right\\}_{{1}}
Note that NaN values are not taken into account in the p-values correction.
References
----------
- Holm, S. (1979). A simple sequentially rejective multiple test procedure.
Scandinavian journal of statistics, 65-70.
- https://en.wikipedia.org/wiki/Holm%E2%80%93Bonferroni_method
Examples
--------
>>> import pingouin as pg
>>> pvals = [.50, .003, .32, .054, .0003]
>>> reject, pvals_corr = pg.multicomp(pvals, method='holm', alpha=.05)
>>> print(reject, pvals_corr)
[False True False False True] [0.64 0.012 0.64 0.162 0.0015]
"""
# Convert to array and save original shape
pvals = np.asarray(pvals)
shape_init = pvals.shape
pvals = pvals.ravel()
num_nan = np.isnan(pvals).sum()
# Sort the (flattened) p-values
pvals_sortind = np.argsort(pvals)
pvals_sorted = pvals[pvals_sortind]
sortrevind = pvals_sortind.argsort()
ntests = pvals.size - num_nan
# Now we adjust the p-values
pvals_corr = np.diag(pvals_sorted * np.arange(ntests, 0, -1)[..., None])
pvals_corr = np.maximum.accumulate(pvals_corr)
pvals_corr = np.clip(pvals_corr, None, 1)
# And revert to the original shape and order
pvals_corr = np.append(pvals_corr, np.full(num_nan, np.nan))
pvals_corrected = pvals_corr[sortrevind].reshape(shape_init)
with np.errstate(invalid='ignore'):
reject = np.less(pvals_corrected, alpha)
return reject, pvals_corrected
def sidak(pvals, alpha=0.05):
"""P-values correction with Sidak method.
Parameters
----------
pvals : array_like
Array of p-values of the individual tests.
alpha : float
Error rate (= alpha level).
Returns
-------
reject : array, bool
True if a hypothesis is rejected, False if not
pval_corrected : array
P-values adjusted for multiple hypothesis testing using the Sidak
procedure.
See also
--------
bonf, holm, fdr, multicomp
Notes
-----
The Sidak adjusted p-values are defined as:
.. math::
\\widetilde {p}_{{(i)}}= 1 - (1 - p_{{(i)}})^{n}
The Sidak correction is slightly more liberal than the Bonferroni
correction.
Note that NaN values are not taken into account in the p-values correction.
References
----------
- Šidák, Z. K. (1967). "Rectangular Confidence Regions for the Means of
Multivariate Normal Distributions". Journal of the American Statistical
Association. 62 (318): 626–633.
- https://en.wikipedia.org/wiki/%C5%A0id%C3%A1k_correction
Examples
--------
>>> import numpy as np
>>> import pingouin as pg
>>> pvals = [.50, .003, .32, .054, .0003]
>>> reject, pvals_corr = pg.multicomp(pvals, method='sidak', alpha=.05)
>>> print(reject, np.round(pvals_corr, 4))
[False True False False True] [0.9688 0.0149 0.8546 0.2424 0.0015]
"""
pvals = np.asarray(pvals)
num_nan = np.isnan(pvals).sum()
ntests = (float(pvals.size) - num_nan)
pvals_corrected = 1 - np.power((1. - pvals), ntests)
pvals_corrected = np.clip(pvals_corrected, None, 1)
with np.errstate(invalid='ignore'):
reject = np.less(pvals_corrected, alpha)
return reject, pvals_corrected
##############################################################################
# EXTERNAL FUNCTION
##############################################################################
def multicomp(pvals, alpha=0.05, method='holm'):
"""P-values correction for multiple comparisons.
Parameters
----------
pvals : array_like
Uncorrected p-values.
alpha : float
Significance level.
method : string
Method used for testing and adjustment of p-values. Can be either the
full name or initial letters. Available methods are ::
'bonf' : one-step Bonferroni correction
'sidak' : one-step Sidak correction
'holm' : step-down method using Bonferroni adjustments
'fdr_bh' : Benjamini/Hochberg FDR correction
'fdr_by' : Benjamini/Yekutieli FDR correction
'none' : pass-through option (no correction applied)
Returns
-------
reject : array, boolean
True for hypothesis that can be rejected for given alpha.
pvals_corrected : array
P-values corrected for multiple testing.
Notes
-----
This function is similar to the `p.adjust` R function.
The correction methods include the Bonferroni correction (``bonf``)
in which the p-values are multiplied by the number of comparisons.
Less conservative methods are also included such as Sidak (1967)
(``sidak``), Holm (1979) (``holm``), Benjamini & Hochberg (1995)
(``fdr_bh``), and Benjamini & Yekutieli (2001) (``fdr_by``), respectively.
The first three methods are designed to give strong control of the
family-wise error rate. Note that the Holm's method is usually preferred.
The ``fdr_bh`` and ``fdr_by`` methods control the false discovery rate,
i.e. the expected proportion of false discoveries amongst the rejected
hypotheses. The false discovery rate is a less stringent condition than
the family-wise error rate, so these methods are more powerful than the
others.
The **Bonferroni** adjusted p-values are defined as:
.. math::
\\widetilde {p}_{{(i)}}= n \\cdot p_{{(i)}}
where :math:`n` is the number of *finite* p-values (i.e. excluding NaN).
The **Sidak** adjusted p-values are defined as:
.. math::
\\widetilde {p}_{{(i)}}= 1 - (1 - p_{{(i)}})^{n}
The **Holm** adjusted p-values are the running maximum of the sorted
p-values divided by the corresponding increasing alpha level:
.. math::
\\widetilde {p}_{{(i)}}=\\max _{{j\\leq i}}\\left\\{(n-j+1)p_{{(j)}}
\\right\\}_{{1}}
The **Benjamini–Hochberg** procedure (BH step-up procedure) controls the
false discovery rate (FDR) at level :math:`\\alpha`. It works as follows:
1. For a given :math:`\\alpha`, find the largest :math:`k` such that
:math:`P_{(k)}\\leq \\frac {k}{n}\\alpha.`
2. Reject the null hypothesis for all
:math:`H_{(i)}` for :math:`i = 1, \\ldots, k`.
The BH procedure is valid when the :math:`n` tests are independent, and
also in various scenarios of dependence, but is not universally valid.
The **Benjamini–Yekutieli** procedure (BY) controls the FDR under arbitrary
dependence assumptions. This refinement modifies the threshold and finds
the largest :math:`k` such that:
.. math::
P_{(k)} \\leq \\frac{k}{n \\cdot c(n)} \\alpha
References
----------
- Bonferroni, C. E. (1935). Il calcolo delle assicurazioni su gruppi
di teste. Studi in onore del professore salvatore ortu carboni, 13-60.
- Šidák, Z. K. (1967). "Rectangular Confidence Regions for the Means of
Multivariate Normal Distributions". Journal of the American Statistical
Association. 62 (318): 626–633.
- Holm, S. (1979). A simple sequentially rejective multiple test procedure.
Scandinavian Journal of Statistics, 6, 65–70.
- Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery
rate: a practical and powerful approach to multiple testing. Journal of
the Royal Statistical Society Series B, 57, 289–300.
- Benjamini, Y., and Yekutieli, D. (2001). The control of the false
discovery rate in multiple testing under dependency. Annals of
Statistics, 29, 1165–1188.
Examples
--------
FDR correction of an array of p-values
>>> import pingouin as pg
>>> pvals = [.50, .003, .32, .054, .0003]
>>> reject, pvals_corr = pg.multicomp(pvals, method='fdr_bh')
>>> print(reject, pvals_corr)
[False True False False True] [0.5 0.0075 0.4 0.09 0.0015]
Holm correction with missing values
>>> import numpy as np
>>> pvals[2] = np.nan
>>> reject, pvals_corr = pg.multicomp(pvals, method='holm')
>>> print(reject, pvals_corr)
[False True False False True] [0.5 0.009 nan 0.108 0.0012]
"""
# Safety check
assert isinstance(pvals, (list, np.ndarray)), "pvals must be list or array"
pvals = np.squeeze(np.asarray(pvals))
assert isinstance(alpha, float), 'alpha must be a float.'
assert isinstance(method, str), 'method must be a string.'
assert 0 < alpha < 1, 'alpha must be between 0 and 1.'
if method.lower() in ['b', 'bonf', 'bonferroni']:
reject, pvals_corrected = bonf(pvals, alpha=alpha)
elif method.lower() in ['h', 'holm']:
reject, pvals_corrected = holm(pvals, alpha=alpha)
elif method.lower() in ['s', 'sidak']:
reject, pvals_corrected = sidak(pvals, alpha=alpha)
elif method.lower() in ['fdr', 'fdr_bh', 'bh']:
reject, pvals_corrected = fdr(pvals, alpha=alpha, method='fdr_bh')
elif method.lower() in ['fdr_by', 'by']:
reject, pvals_corrected = fdr(pvals, alpha=alpha, method='fdr_by')
elif method.lower() == 'none':
pvals_corrected = pvals
with np.errstate(invalid='ignore'):
reject = np.less(pvals_corrected, alpha)
else:
raise ValueError('Multiple comparison method not recognized')
return reject, pvals_corrected