Instantly share code, notes, and snippets.

# alimuldal/dunn.py

Last active July 19, 2021 11:55
Implementation of Dunn's multiple comparison test, following a Kruskal-Wallis 1-way ANOVA
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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_dunn(groups, to_compare=None, alpha=0.05, method='bonf'): """ Kruskal-Wallis 1-way ANOVA with Dunn'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) method: string method used to adjust p-values to account for multiple corrections (see statsmodels.stats.multitest.multipletests for options) 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: --------------- Gibbons, J. D., & Chakraborti, S. (2011). Nonparametric Statistical Inference (5th ed., pp. 353-357). Boca Raton, FL: Chapman & Hall. """ # 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 upper quantiles of the standard # normal distribution p_uncorrected = stats.norm.sf(Z_pairs) * 2. # correction for multiple comparisons reject, p_corrected, alphac_sidak, alphac_bonf = multipletests( p_uncorrected, method=method ) return H, p_omnibus, Z_pairs, p_corrected, reject

### oliviac12 commented Jun 29, 2016

Hi, can you please give an example for groups?

### davidadamojr commented Aug 20, 2016

It seems as though groups is a two-dimensional list of samples e.g. [[1,2,3,4,5], [4,5,6,7,8,9]].

### jazon33y commented Sep 27, 2016

I think there might be an addition needed, although I could be completely wrong:

``````        # standardized score
ts3_ts = list(np.unique(allgroups, return_counts=True))
E_ts3_ts = sum([x**3 - x for x in ts3_ts if x>1])

if sum([x>1 for x in ts3_ts]) > 0:
warnings.warn("We see ties.")

yi = np.abs(Rmean[ii] - Rmean[jj])
theta10 = (N * (N + 1)) / 12
theta11 =  E_ts3_ts / ( 12* (N - 1) )
theta2 = (1 / n[ii] + 1 / n[jj])
theta = np.sqrt( (theta10 - theta11) * theta2 )
Zij = yi / theta
else:
Zij = (np.abs(Rmean[ii] - Rmean[jj]) /
np.sqrt((1. / 12.) * N * (N + 1) * (1. / n[ii] + 1. / n[jj])))
``````

### ricardoV94 commented May 26, 2017

@jazon33y You are correct. After adding your snippet, the output is now equal to that of R when there are ties. Thanks!

### whokan commented Sep 6, 2017

What order do the post host arrays output in? From my own testing I think for three groups it's (1 vs 2), (1 vs 3), (2 vs 3)?

### ayesop commented Nov 20, 2017 • edited

stats.chisqprob is deprecated in scipy 0.17.0; use stats.distributions.chi2.sf
`p_omnibus = stats.distributions.chi2.sf(H, df_omnibus)`

The format of 'groups' is still unclear; I get a "too many values to unpack" error when I use an array of samples as indicated.

### stuboo commented Jan 10, 2018

@whokan I think your assumption is correct. I just tested 7 groups and this was the return order:

`[(0,1), (0,2), (0,3), (0,4), (0,5), (0,6), (1,2), (1,3), (1,4), (1,5), (1,6), (2,3), (2,4), (2,5), (2,6), (3,4), (3,5), (3,6), (4,5), (4,6), (5,6)]`

### eturkes commented Jan 16, 2018 • edited

Hello there @alimuldal I would like to include your code in a GPLv3 project. Your function does not seem to be licensed, so I'd like to know if I have permission to use and modify it. And if so, how would you like me to credit you? Would this be OK?

# The following section is adapted from code written by Alistair Muldal.
# Permission for use and modification has been granted by the author.
# Source: https://gist.github.com/alimuldal/fbb19b73fa25423f02e8
# Retrieved: 2018.10.12 UTC

### jmassall commented Jan 17, 2018

@jazon33y thanks for the correction (now if you have k=2 samples the p-value is consistent with K-W test as it should). FYI for Python 2.x users that don't future import division, you will get Inf for Z-scores (and corresponding p-values of 0) because of the integer division for theta2.

### ghost commented Mar 19, 2018

Hi all,
did somebody tried the script with 3 groups or more? My p-values are wrong (compared with R, and I checked the same adjustment method was applied).
If it does work with you can you please share the script?
Thanks

### glickmac commented Nov 30, 2018 • edited

If you are getting an error with 'chisqprob' (see so question here) the Bug Fix is Here via @VincentLa14

`from scipy import stats`
`stats.chisqprob = lambda chisq, df: stats.chi2.sf(chisq, df)`