pvalue combining tests

Test combining pvalues with different fixed reference time using the fisher method.

import numpy, pylab
from scipy.stats import combine_pvalues


def ftop(f, ft):
   return 1 - numpy.exp(-ft * f)

def ptof(p, ft):
   return numpy.log(1.0-p) / -ft


ftime = 0.1
ifars = [0.01, 0.1, 1, 10, 100]
pvalues = numpy.logspace(-4, 0, num=50)

for ifar in ifars:
   p1 = ftop(1.0/ifar, .01 * ifar)

   cp = numpy.array([combine_pvalues([p1, p2])[1] for p2 in pvalues])

   nifar = 1.0 / ptof(cp, .01 * ifar)

   pylab.plot(pvalues, nifar / ifar, label='Original IFAR = %s' % ifar)


pylab.plot(pvalues, 1.0/pvalues, linestyle='--')
pylab.grid()
pylab.xscale('log')
pylab.yscale('log')
pylab.ylabel('ratio of ifars (combined / original)')
pylab.xlabel('added pvalue')
pylab.legend()
pylab.savefig('farchange.png')

import numpy, pylab
from scipy.stats import combine_pvalues
from numpy.random import uniform
from numpy.random import seed
seed(0)

Compare (1) scaling reference time to taking (2) taking the max of the combined pvalue and original pvalue i.e. to minimize significance loss of certain candidates and (3) fisher method with fixed reference time.

import numpy, pylab
from scipy.stats import combine_pvalues
from numpy.random import uniform
from numpy.random import seed
seed(0)

def ptof(p, ft):
   return numpy.log(1.0-p) / -ft

def ftop(f, ft):
   return 1 - numpy.exp(-ft * f)

size = 1000000

p1 = uniform(0, 1, size=size)
p2 = uniform(0, 1, size=size)

cp = [combine_pvalues([v1, v2], method='fisher')[1] for v1, v2 in zip(p1, p2)]
cp = numpy.array(cp)

ftime = .1 # Choosing units so this is a little over a month
btime = ftime * size

ifars = 1.0 / ptof(p1, ftime)
cifars = 1.0 / ptof(cp, ftime)


mifars = numpy.maximum(ifars, cifars) / 2.0
mifars.sort()

ifars.sort()
cifars.sort()

cnum = numpy.arange(size, 0, -1) / btime

p3 = ftop(cnum, .01* 1.0/cnum)
print p3
p3 = p3[0]
cp2 = [combine_pvalues([p3, v2], method='fisher')[1] for v2 in p2]
cp2 = numpy.array(cp2)
sifar = 1.0 / ptof(cp2, .01 * 1.0/cnum)
sifar.sort()

pylab.plot(sifar, cnum, label='scaling')
pylab.plot(mifars, cnum, label='max(combined, coinc)')
pylab.plot(cifars, cnum, label='combined (Fisher)')
pylab.plot(ifars, cnum, label='ifars uniform pvalues')
pylab.plot(1./cnum, cnum, linestyle='--')
pylab.legend()
pylab.xscale('log')
pylab.yscale('log')
pylab.ylabel('Cumulative Rate')
pylab.xlabel('IFAR (yrs)')
pylab.savefig('combining.png')