Multiple Pairwise Comparisons (Post Hoc) Tests in Python

A quick fix for Issue #51 is by adding norm=colors.NoNorm() to ColorbarBase. Probably should be double-checked to ensure the behaviour is correct for all cases.

import matplotlib.pyplot as plt
import scikit_posthocs as sp
import statsmodels.api as sa

x = sa.datasets.get_rdataset('iris').data
x.columns = x.columns.str.replace('.', '')
pc = sp.posthoc_ttest(x, val_col='SepalWidth', group_col='Species', p_adjust='holm')
heatmap_args = {'linewidths': 0.25, 'linecolor': '0.5', 'clip_on': False, 'square': True,
                'cbar_ax_bbox': [0.82, 0.35, 0.04, 0.3]}
sp.sign_plot(pc, **heatmap_args)
plt.show()

Describe the bug outliers_gesd has a bug ,when the outliers increase and abs_d numpy array size decrease;

in file '_outliers.py' line

      # Masked values
        lms = ms[-1] if len(ms) > 0 else []
        ms.append(lms + [np.argmax(abs_d)])

the abs_d's size is not data's size any more, so the np.argmax(abs_d) is not the true outlier index in data numpy array.

First off, congrats on the great idea of porting this to Python. I just ran it on R's InsectSprays data (with no p_adjust method) and I'm getting slightly different p-values than those from PMCMR for a couple of the group comparisons (B & C and F & C). Do you know why that might be? Thanks!

in definition posthoc_dscf, ni and nj have been replaced. CURRENT

def posthoc_dscf(a, val_col=None, group_col=None, sort=False):
....
　　def compare(i, j):
　　　　....
   　　　　     u = np.array([nj * ni + (nj * (nj + 1) / 2), 
　　　　　　　　　　　　nj * ni + (ni * (ni + 1) / 2)]) - r
　　　　....

TRUTH

 u = np.array([nj * ni + (ni * (ni + 1) / 2), 
　　　　　　nj * ni + (nj* (nj + 1) / 2)]) - r

Hi there!

Thanks for the nice library!

A small thing to suggest on sign_plot method, I believe it's better to put legend in the order to ['NS', 'p<0.05', 'p<0.01', 'p<0.001'] rather than the current version with 'NS' at the end. Because 'NS' is the situation where p>0.05, and it's more logical to sort the colormap on either descending or ascending order.

Best

Hi,

I cannot use post-hocs test for dataframes with melted = True and group_col != 'groups', block_col != 'blocks' and y_col != 'y'. Basically, anything which deviates from the example

sp.posthoc_nemenyi_friedman(data, y_col='y', block_col='blocks', group_col='groups', melted=True)

breaks the code. The error is likely due to __convert_to_block_df (https://github.com/maximtrp/scikit-posthocs/blob/master/scikit_posthocs/_posthocs.py) which returns the old y_col, group_col, block_col values but assigns the column names "groups" / "blocks" / "y"

def __convert_to_block_df(a, y_col=None, group_col=None, block_col=None, melted=False):
    # ...
    elif isinstance(a, DataFrame) and melted:
        x = DataFrame.from_dict({'groups': a[group_col],
                                 'blocks': a[block_col],
                                 'y': a[y_col]})**
    # ...
    return x, y_col, group_col, block_col

On a somewhat related note: I wanted to implement / use these tests to plot CD diagrams as suggested in "J. Demsar (2006), Statistical comparisons of classifiers over multiple data sets, Journal of Machine Learning Research, 7, 1-30." which you also cite in the documentation. However, I have a difficult time to understand what "block", "groups", and "y" mean in this context. More specifically, are blocks (or groups?) different classifiers or datasets and is y the ranks or the accuracies? You dont happen to have some example code and or explanation how to plot CD diagrams?

Thank

Describe the bug A clear and concise description of what the bug is.

Dataset Please provide a link to the dataset you get the bug with.

To Reproduce Steps to reproduce the behavior:

1. Go to '...'
2. Click on '....'
3. Scroll down to '....'
4. See error

Expected behavior A clear and concise description of what you expected to happen.

System and package information (please complete the following information):

• OS: (e.g. Linux 4.20.0-arch1-1-ARCH x86_64 GNU/Linux)
• Package version: (e.g. 0.4.0)

Additional context Add any other context about the problem here.

Hi

When I try to run a Conover posthoc on my pandas dataframe I get the following error: /opt/conda/lib/python3.6/site-packages/scikit_posthocs/_posthocs.py:705: RuntimeWarning: invalid value encountered in sqrt tval = dif / np.sqrt(A / B)

I also suspect the Nemenyi to be affected as I am getting values of 0.001 exactly for every conditions.

Best,

Hi there,

Thanks for the package! I am using friedman tests paired with the Nemenyi test, and this works very nicely for looking at all pairwise comparisons.

I was wondering if there was a way to compare vs a control method using say the Holm method? I can see there is a Holm option with both conover and siegel, but I do not believe these compare against a control (correct me if Im wrong).

Thank you

I'm really excited to use the new posthocs package, and have been trying to run the Dunn test (posthoc_dunn) on results from a survey over the last few days (I have about 1100 respondents). I have no problem when I run it on results that represent the difference between two feeling thermometers (a variable that ranges from -100 to 100). But every time I try to run it on a Likert-type scale item that takes values of 1 through 3 or 1 through 5, it returns a table full of null results (NaN in all cells except the diagonal). This comes along with a series of warnings as follows:

1. _posthocs.py:191: RuntimeWarning: invalid value encountered in sqrt: z_value = diff / np.sqrt((A - x_ties) * B)
2. multitest.py:176: RuntimeWarning: invalid value encountered in greater notreject = pvals > alphaf / np.arange(ntests, 0, -1)
3. multitest.py:251: RuntimeWarning: invalid value encountered in greater pvals_corrected[pvals_corrected>1] = 1

I am not a programming expert, but my impression is that what is happening here is that the compare_dunn function (lines 187-193 in posthoc.py) is not returning valid p-values, and I am guessing that this is because (A - x_ties) is negative for some reason and so the np.sqrt function isn't computing a value for the z_value.

I played around with some groups of small arrays involving combinations of values ranging from 1 to 3 and 1 to 5, on the same scale as my data. Sometimes these had no problem returning valid results and other times they yielded the same NaNs that I get with my full dataset. I'm wondering if the issue has something to do with the total number or overall proportion of ties in the data. Obviously with Likert-type scale items there are a lot of ties. I'd love your thoughts on whether it's something that can be fixed to make analysis on this type of data possible. Thanks!!

Hi,

I have been using scikit_posthocs(Version: 0.6.6) in Python 3.7.1 for a while. Today, when I tried to import scikit_posthocs, I got an error as ModuleNotFoundError: No module named 'statsmodels.stats.libqsturng'

I checked with pip show statsmodels and Version: 0.12.0 was shown. I reinstalled by 'pip install statsmodels' but got the same error message. I reinstalled by 'pip install scikit-posthocs' but problem remained. Running out of ideas. What could be the reason?

Thanks!

It seems there are many post-hoc tests implemented and will be really good if there's an easy interface to plot critical difference diagrams with any chosen pairwise tests.

In Orange framework, there's only few types of tests implemented and the interface isn't very intuitive.

Hi,

I often use your _posthoc_mannwhitneyu and I get ValueError 'All numbers are identical in mannwhitneyu' when two groups are composed from idential numbers. But I thought we should adjust p-values including the p-value(=1.0) from those comparisons, so I modified the code in _pothoc.py like this.

def _posthoc_mannwhitney(
        a: Union[list, np.ndarray, DataFrame],
        val_col: str = None,
        group_col: str = None,
        use_continuity: bool = True,
        alternative: str = 'two-sided',
        p_adjust: str = None,
        sort: bool = True) -> DataFrame:
    '''Pairwise comparisons with Mann-Whitney rank test.

    Parameters
    ----------
    a : array_like or pandas DataFrame object
        An array, any object exposing the array interface or a pandas
        DataFrame. Array must be two-dimensional.

    val_col : str, optional
        Name of a DataFrame column that contains dependent variable values (test
        or response variable). Values should have a non-nominal scale. Must be
        specified if `a` is a pandas DataFrame object.

    group_col : str, optional
        Name of a DataFrame column that contains independent variable values
        (grouping or predictor variable). Values should have a nominal scale
        (categorical). Must be specified if `a` is a pandas DataFrame object.

    use_continuity : bool, optional
        Whether a continuity correction (1/2.) should be taken into account.
        Default is True.

    alternative : ['two-sided', 'less', or 'greater'], optional
        Whether to get the p-value for the one-sided hypothesis
        ('less' or 'greater') or for the two-sided hypothesis ('two-sided').
        Defaults to 'two-sided'.

    p_adjust : str, optional
        Method for adjusting p values.
        See statsmodels.sandbox.stats.multicomp for details.
        Available methods are:
        'bonferroni' : one-step correction
        'sidak' : one-step correction
        'holm-sidak' : step-down method using Sidak adjustments
        'holm' : step-down method using Bonferroni adjustments
        'simes-hochberg' : step-up method  (independent)
        'hommel' : closed method based on Simes tests (non-negative)
        'fdr_bh' : Benjamini/Hochberg  (non-negative)
        'fdr_by' : Benjamini/Yekutieli (negative)
        'fdr_tsbh' : two stage fdr correction (non-negative)
        'fdr_tsbky' : two stage fdr correction (non-negative)

    sort : bool, optional
        Specifies whether to sort DataFrame by group_col or not. Recommended
        unless you sort your data manually.

    Returns
    -------
    result : pandas.DataFrame
        P values.

    Notes
    -----
    Refer to `scipy.stats.mannwhitneyu` reference page for further details.

    Examples
    --------
    >>> x = [[1,2,3,4,5], [35,31,75,40,21], [10,6,9,6,1]]
    >>> sp.posthoc_mannwhitney(x, p_adjust = 'holm')
    '''
    x, _val_col, _group_col = __convert_to_df(a, val_col, group_col)
    x = x.sort_values(by=[_group_col, _val_col], ascending=True) if sort else x

    groups = x[_group_col].unique()
    x_len = groups.size
    vs = np.zeros((x_len, x_len))
    xg = x.groupby(_group_col)[_val_col]
    tri_upper = np.triu_indices(vs.shape[0], 1)
    tri_lower = np.tril_indices(vs.shape[0], -1)
    vs[:, :] = 0

    combs = it.combinations(range(x_len), 2)

    for i, j in combs: ##I modified this section##
        try:
            vs[i, j] = ss.mannwhitneyu(
                xg.get_group(groups[i]),
                xg.get_group(groups[j]),
                use_continuity=use_continuity,
                alternative=alternative)[1]
        except ValueError as e:
            if str(e)=="All numbers are identical in mannwhitneyu":
                vs[i, j] =1.0
            else:
                raise e

    if p_adjust:
        vs[tri_upper] = multipletests(vs[tri_upper], method=p_adjust)[1]

    vs[tri_lower] = np.transpose(vs)[tri_lower]
    np.fill_diagonal(vs, 1)
    return DataFrame(vs, index=groups, columns=groups)

Is this a right solution?

I'm not sure but this error may not occur with other versions of scipy.stats.

Hi! Thanks a lot for creating this analysis tool.

I would like to check if it is normal that a post-hoc analysis using Dunn test, after Kruskal-Wallis, returns no significant result at all between the pairwise comparisons?

Another question, does Dunn test require multiple comparison correction? Either way (with or without correction), I don't get any significant even though Kruskal-Wallis test rejects the null hypothesis.

Hi, I was wondering if there is any chance to include a feature where post hoc results are grouped according to their relationship. I know that in R there are the packages multcompLetters and multcompview, which offer such feature. I could find some people looking for a feature like this, but no feasible was found.

Example: https://stackoverflow.com/questions/48841650/python-algorithm-on-letter-based-representation-of-all-pairwise-comparisons

There is a solution attempt at those topics, but I could not reproduce them: https://stackoverflow.com/questions/43987651/tukey-test-grouping-and-plotting-in-scipy https://stackoverflow.com/questions/49963138/label-groups-from-tuekys-test-results-according-to-significant

It looks like there is a paper describing the algorithm for implementing this: Hans-Peter Piepho (2004) An Algorithm for a Letter-Based Representation of All-Pairwise Comparisons, Journal of Computational and Graphical Statistics, 13:2, 456-466, DOI: 10.1198/1061860043515

By the way, thanks for this project, it is awesome!

Hi Maksim!

Nice package, thanks for sharing! I am trying to use it and hoping to contribute.

Before I use it, I'd like to understand how it is tested. For example, if you look at https://github.com/maximtrp/scikit-posthocs/blob/4709cf2821ce98dffef542b9b916f7c4a5f00ff4/tests/test_posthocs.py#L27

or any other tests, you seem have to pre-selected the outputs/results to compare.. Where do they come from? From the PMCMR R-package? I tried to look into the PMCMR package to look at their tests - I can't see them in the package distributed, do you know if it is tested?

What are the other principles behind the tests you wrote?