Commit ba2039cc authored by Bedartha Goswami's avatar Bedartha Goswami
Browse files

add 'holm' function; remove save-to-disk lines

parent 59d1cd45
......@@ -6,7 +6,7 @@ A suite of functions that assist in carrying out the event-detection analysis
"""
# Created: Thu May 26, 2016 03:50PM
# Last modified: Tue Aug 15, 2017 01:25pm
# Last modified: Tue Jan 09, 2018 01:53PM
# Copyright: Bedartha Goswami <goswami@pik-potsdam.de>
import numpy as np
......@@ -14,7 +14,7 @@ from scipy.stats import percentileofscore as percentile
from lib import DATPATH
from networks import link_density
from utils import _printmsg, _savedat
from utils import _printmsg
from utils import _progressbar_start
from utils import _progressbar_update
from utils import _progressbar_finish
......@@ -142,30 +142,10 @@ def community_strength_random_model(G, time, ld, wsize, wstep, nsurr,
return Qsurr, LDsurr, tmid
def pvalue(name, dtype, ld, wsize, wstep, nsurr,
verbose=False, save=True, pbar=False):
def pvalues(Q, Qsurr, nwind, verbose=False, save=True, pbar=False):
"""
Returns the p-value of obtaining intra-community link density.
"""
_printmsg("p-VALUE OF AN EVENT FOR %s" % name.upper(), verbose)
_printmsg("=========================================", verbose)
_printmsg("Loading data...", verbose)
prefix = "%s/%s/community_strength" % (DATPATH, dtype)
suffix = "LD%.2f/WSIZE%d/%s_WSTEP%d.npz" % (ld, wsize, name, wstep)
f_data = "%s_data/%s" % (prefix, suffix)
f_cfmd = "%s_random_model/%s" % (prefix, suffix)
d_data = np.load(f_data)
d_cfmd = np.load(f_cfmd)
# check for consistency
tmid1 = d_data["tmid"]
tmid2 = d_cfmd["tmid"]
assert all(tmid1 == tmid2), "Data & random model timescales don't match!"
assert nsurr == d_cfmd["nsurr"], "Incompatible no. of surrogates!"
# initiate relevant variables
time = tmid1
Q = d_data["Q"]
Qsurr = d_cfmd["Qsurr"]
nwind = d_data["nwind"]
_printmsg("Estimating p-value...", verbose)
pvalue = np.zeros(nwind)
prog_bar = _progressbar_start(nwind, pbar)
......@@ -174,12 +154,25 @@ def pvalue(name, dtype, ld, wsize, wstep, nsurr,
pvalue[i] = 1. - percentile(Qnull, Q[i], kind="weak") / 100.
_progressbar_update(prog_bar, i)
_progressbar_finish(prog_bar)
p = "%s/%s/pvalues/" % (DATPATH, dtype)
o = "LD%.2f_WSIZE%d_%s_WSTEP%d" % (ld, wsize, name, wstep)
_savedat(p, o, verbose, save,
pvalue=pvalue,
Qsurr=Qsurr, Q=Q, time=time,
wsize=wsize, wstep=wstep,
nwind=nwind, nsurr=nsurr,
)
return None
return pvalues
def holm(pvalues, alpha=0.05, corr_type="dunn"):
"""
Returns indices of p-values using Holm's method for multiple testing.
"""
n = len(pvalues)
sortidx = np.argsort(pvalues)
p_ = pvalues[sortidx]
j = np.arange(1, n + 1)
if corr_type == "bonf":
corr_factor = alpha / (n - j + 1)
elif corr_type == "dunn":
corr_factor = 1. - (1. - alpha) ** (1. / (n - j + 1))
try:
idx = np.where(p_ <= corr_factor)[0][-1]
idx = sortidx[:idx]
except IndexError:
idx = []
return idx
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment