Skip to content

Instantly share code, notes, and snippets.

@josef-pkt
Forked from alimuldal/nemenyi.py
Last active August 29, 2015 14:10
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save josef-pkt/cbb72d8031a7853c854d to your computer and use it in GitHub Desktop.
Save josef-pkt/cbb72d8031a7853c854d to your computer and use it in GitHub Desktop.
import numpy as np
from scipy import stats
from itertools import combinations
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.libqsturng import psturng
import warnings
def kw_nemenyi(groups, to_compare=None, alpha=0.05):
"""
Kruskal-Wallis 1-way ANOVA with Nemenyi's multiple comparison test
Arguments:
---------------
groups: sequence
arrays corresponding to k mutually independent samples from
continuous populations
to_compare: sequence
tuples specifying the indices of pairs of groups to compare, e.g.
[(0, 1), (0, 2)] would compare group 0 with 1 & 2. by default, all
possible pairwise comparisons between groups are performed.
alpha: float
family-wise error rate used for correcting for multiple comparisons
(see statsmodels.stats.multitest.multipletests for details)
Returns:
---------------
H: float
Kruskal-Wallis H-statistic
p_omnibus: float
p-value corresponding to the global null hypothesis that the medians of
the groups are all equal
Z_pairs: float array
Z-scores computed for the absolute difference in mean ranks for each
pairwise comparison
p_corrected: float array
corrected p-values for each pairwise comparison, corresponding to the
null hypothesis that the pair of groups has equal medians. note that
these are only meaningful if the global null hypothesis is rejected.
reject: bool array
True for pairs where the null hypothesis can be rejected for the given
alpha
Reference:
---------------
"""
# omnibus test (K-W ANOVA)
# -------------------------------------------------------------------------
groups = [np.array(gg) for gg in groups]
k = len(groups)
n = np.array([len(gg) for gg in groups])
if np.any(n < 5):
warnings.warn("Sample sizes < 5 are not recommended (K-W test assumes "
"a chi square distribution)")
allgroups = np.concatenate(groups)
N = len(allgroups)
ranked = stats.rankdata(allgroups)
# correction factor for ties
T = stats.tiecorrect(ranked)
if T == 0:
raise ValueError('All numbers are identical in kruskal')
# sum of ranks for each group
j = np.insert(np.cumsum(n), 0, 0)
R = np.empty(k, dtype=np.float)
for ii in range(k):
R[ii] = ranked[j[ii]:j[ii + 1]].sum()
# the Kruskal-Wallis H-statistic
H = (12. / (N * (N + 1.))) * ((R ** 2.) / n).sum() - 3 * (N + 1)
# apply correction factor for ties
H /= T
df_omnibus = k - 1
p_omnibus = stats.chisqprob(H, df_omnibus)
# multiple comparisons
# -------------------------------------------------------------------------
# by default we compare every possible pair of groups
if to_compare is None:
to_compare = tuple(combinations(range(k), 2))
ncomp = len(to_compare)
Z_pairs = np.empty(ncomp, dtype=np.float)
p_uncorrected = np.empty(ncomp, dtype=np.float)
Rmean = R / n
for pp, (ii, jj) in enumerate(to_compare):
# standardized score
Zij = (np.abs(Rmean[ii] - Rmean[jj]) /
np.sqrt((1. / 12.) * N * (N + 1) * (1. / n[ii] + 1. / n[jj])))
Z_pairs[pp] = Zij
# corresponding p-values obtained from the upper quantiles of the
# studentized range distribution
p_corrected = psturng(Z_pairs * np.sqrt(2), ncomp, np.inf)
reject = p_corrected <= alpha
return H, p_omnibus, Z_pairs, p_corrected, reject
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment