Commits

Davide Cittaro committed 95c2d18

added FDR

Comments (0)

Files changed (1)

+import numpy as np
+def p_adjust(p, method = 'BH', n = None):
+  if not n:
+    n = len(p)
+  p = np.array(p)
+  p0 = p
+  nna = np.isnan(p) ^ True
+  p = p[nna]
+  lp = len(p)
+  if n <= 1:
+    return p0
+  if method == 'bonferroni':
+    p0[nna] = np.minimum(1, n * p)
+  elif method == 'holm':
+    i = 1 + np.arange(lp)
+    o = np.argsort(p)
+    ro = np.argsort(o)
+    p0[nna] = np.minimum(1, np.maximum.accumulate((n - i + 1.) * p[o]))[ro]
+  elif method == 'hochberg':
+    i = np.arange(lp, 0, -1.)
+    o = np.argsort(p)
+    o = o[::-1]
+    ro = np.argsort(o)
+    p0[nna] = np.minimum(1, np.minimum.accumulate((n - i + 1.) * p[o] ))[ro]
+  elif method == 'BH':
+    i = np.arange(lp, 0, -1.)
+    o = np.argsort(p)
+    o = o[::-1]
+    ro = np.argsort(o)
+    p0[nna] = np.minimum(1, np.minimum.accumulate(n / i * p[o]))[ro]
+  elif method == 'BY':
+    i = np.arange(lp, 0, -1.)
+    o = np.argsort(p)
+    o = o[::-1]
+    ro = np.argsort(o)
+    q = np.sum(1. / np.arange(1, n + 1))
+    p0[nna] = np.minimum(1, np.minimum.accumulate(q * n / i * p[o]))[ro]
+  elif method == 'hommel':
+    if n > lp:
+      p = np.append(p, np.ones(n - lp))
+    i = 1 + np.arange(n)
+    o = np.argsort(p)
+    ro = np.argsort(o)
+    p = p[o]
+    q = pa = np.min(n * p / i) * np.ones(n)
+    for j in range(n , 2, -1):
+      ij = np.arange( n - j + 1)
+      i2 = np.arange(n - j + 1, n)    
+      q1 = np.min(j * p[i2]/np.arange(2, j + 1))
+      q[ij] = np.minimum(j * p[ij], q1)
+      q[i2] = q[ n - j ]
+      pa = np.maximum(pa, q)
+    if lp < n:
+      mask = ro[:lp]
+    else:
+      mask = ro    
+    p0[nna] = np.maximum(pa, p)[ro]
+  return p0  
+      
+