Module stikpetP.other.poho_nemenyi

Functions

def ph_nemenyi(catField, ordField, categories=None, levels=None, version='auto')

Post-Hoc Nemenyi Test

This can be used as a post-hoc test for a Kruskal-Wallis test (see ts_kruskal_wallis()).

The test compares each possible pair of categories from the catField and their mean rank. The null hypothesis is that these are then equal.

Pohlert (2016) mentions the exact version should only be used if there are no ties, and suggest to use a chi-square alternative in case of ties. This is referred to by Zaiontz (n.d.-b) as the Schaich-Hamerle test (1984).

The ties correction is taken from Pohlert (2016).

Other post-hoc tests that could be considered are Dunn, Steel-Dwass, Conover, a pairwise Mann-Whitney U, or pairwise Mood-Median.

Parameters

catField : pandas series

data with categories

ordField : pandas series

data with the scores

categories : list or dictionary, optional

the categories to use from catField

levels : list or dictionary, optional

the levels or order used in ordField.

version : {"auto", "exact", "sh", "sh-ties"}, optional

version of the test to use

Returns

A dataframe with:

 

• cat. 1, one of the two categories being compared
• cat. 2, second of the two categories being compared
• n1, number of cat. 1. cases in comparison
• n2, number of cat. 2 cases in comparison
• mean rank 1, mean rank of cases in cat. 1, based on all cases (incl. categories not in comparison)
• mean rank 2, mean rank of cases in cat. 2, based on all cases (incl. categories not in comparison)
• se, the standard error used
• statistic, the q-value of the test
• p-value, the p-value (significance)

Notes

The formula used (Pohlert, 2016, p. 3): $$q_{1,2} = \frac{\bar{r}_1 - \bar{r}_2}{\sqrt{\frac{n\times\left(n+1\right)}{24}\times\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}}$$ $$sig. = 1 - Q\left(q_{1,2}, k, df=\infty\right)$$

A chi-square distribution can also be used. Zaiontz (n.d.) and BrightStat (n.d.) refer to this as the Schaich-Hamerle test.

The formula used then changes to (Sachs, 1982, p. 549): $$\chi_{1,2}^2 = \frac{\left(\bar{r}_1 - \bar{r}_2\right)^2}{\frac{n\times\left(n+1\right)}{12}\times\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}$$ $$df = k - 1$$ $$sig. = 1 - \chi^2\left(\chi_{1,2}^2, df\right)$$

A ties correction found in Pohlert (2016, p. 3) adjusts this to: $$\chi_{1,2}^2 = \frac{\left(\bar{r}_1 - \bar{r}_2\right)^2}{\left(1-T\right)\times\frac{n\times\left(n+1\right)}{12}\times\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}$$ $$T = \frac{\sum t_i^3 - t_i}{n^3 - n}$$

The original formula is most likely from Nemenyi (1963) and the Schaich and Hamerle (1984).

Symbols used

• $k$, the number of categories
• $t_j$, the frequency of the j-th unique rank.
• $n$, the total sample size, of all scores, incl. those not in the comparison
• $n_i$, the number of scores in category i
• $r_{i,j}$, the rank of the j-th score in category i, using all original scores (incl. those not in the comparison).
• $\bar{r}_i$, the average of the ranks in category i, using all original scores (incl. those not in the comparison).
• $Q\left(\dots\right)$, the cumulative distribution function of the standardized range distribution.
• $\chi^2\left(\dots\right)$, the cumulative distribution function of the chi-square distribution.

References

BrightStat. (n.d.). Kruskal-Wallis test. BrightStat. Retrieved October 25, 2023, from https://secure.brightstat.com/index.php?p=c&d=1&c=2&i=7

Nemenyi, P. (1963). Distribution-free Multiple Comparisons. Princeton University.

Pohlert, T. (2016). The pairwise multiple comparison of mean ranks package (PMCMR). https://cran.r-hub.io/web/packages/PMCMR/vignettes/PMCMR.pdf

Sachs, L. (1982). Applied statistics: A handbook of techniques. Springer-Verlag.

Schaich, E., & Hamerle, A. (1984). Verteilungsfreie statistische Prüfverfahren. Springer. doi:10.1007/978-3-642-70032-3

Zaiontz, C. (n.d.). Schaich-Hamerle Test after KW. Real Statistics Using Excel. Retrieved October 25, 2023, from https://real-statistics.com/one-way-analysis-of-variance-anova/kruskal-wallis-test/schaich-hamerle-test/

Author

Made by P. Stikker

Companion website: https://PeterStatistics.com
YouTube channel: https://www.youtube.com/stikpet
Donations: https://www.patreon.com/bePatron?u=19398076

Expand source code

def ph_nemenyi(catField, ordField, categories=None, levels=None, version="auto"):
    '''
    Post-Hoc Nemenyi Test
    ---------------------
    This can be used as a post-hoc test for a Kruskal-Wallis test (see ts_kruskal_wallis()).
    
    The test compares each possible pair of categories from the catField and their mean rank. The null hypothesis is that these are then equal. 
    
    Pohlert (2016) mentions the exact version should only be used if there are no ties, and suggest to use a chi-square alternative in case of ties. This is referred to by Zaiontz (n.d.-b) as the Schaich-Hamerle test (1984).
    
    The ties correction is taken from Pohlert (2016).
    
    Other post-hoc tests that could be considered are Dunn, Steel-Dwass, Conover, a pairwise Mann-Whitney U, or pairwise Mood-Median.
    
    Parameters
    ----------
    catField : pandas series
        data with categories
    ordField : pandas series
        data with the scores
    categories : list or dictionary, optional
        the categories to use from catField
    levels : list or dictionary, optional
        the levels or order used in ordField.
    version : {"auto", "exact", "sh", "sh-ties"}, optional
        version of the test to use
    
    Returns
    -------
    A dataframe with:
    
    * *cat. 1*, one of the two categories being compared
    * *cat. 2*, second of the two categories being compared
    * *n1*, number of cat. 1. cases in comparison
    * *n2*, number of cat. 2 cases in comparison
    * *mean rank 1*, mean rank of cases in cat. 1, based on all cases (incl. categories not in comparison)
    * *mean rank 2*, mean rank of cases in cat. 2, based on all cases (incl. categories not in comparison)
    * *se*, the standard error used
    * *statistic*, the q-value of the test
    * *p-value*, the p-value (significance)
    
    Notes
    -----
    The formula used (Pohlert, 2016, p. 3):
    $$q_{1,2} = \\frac{\\bar{r}_1 - \\bar{r}_2}{\\sqrt{\\frac{n\\times\\left(n+1\\right)}{24}\\times\\left(\\frac{1}{n_1}+\\frac{1}{n_2}\\right)}}$$
    $$sig. = 1 - Q\\left(q_{1,2}, k, df=\\infty\\right)$$
    
    A chi-square distribution can also be used. Zaiontz (n.d.) and BrightStat (n.d.) refer to this as the Schaich-Hamerle test. 
    
    The formula used then changes to (Sachs, 1982, p. 549):
    $$\\chi_{1,2}^2 = \\frac{\\left(\\bar{r}_1 - \\bar{r}_2\\right)^2}{\\frac{n\\times\\left(n+1\\right)}{12}\\times\\left(\\frac{1}{n_1}+\\frac{1}{n_2}\\right)}$$
    $$df = k - 1$$
    $$sig. = 1 - \\chi^2\\left(\\chi_{1,2}^2, df\\right)$$
    
    A ties correction found in Pohlert (2016, p. 3) adjusts this to:
    $$\\chi_{1,2}^2 = \\frac{\\left(\\bar{r}_1 - \\bar{r}_2\\right)^2}{\\left(1-T\\right)\\times\\frac{n\\times\\left(n+1\\right)}{12}\\times\\left(\\frac{1}{n_1}+\\frac{1}{n_2}\\right)}$$
    $$T = \\frac{\\sum t_i^3 - t_i}{n^3 - n}$$
    
    The original formula is most likely from Nemenyi (1963) and the Schaich and Hamerle (1984).
    
    *Symbols used*
    
    * \\(k\\), the number of categories
    * \\(t_j\\), the frequency of the j-th unique rank.
    * \\(n\\), the total sample size, of all scores, incl. those not in the comparison
    * \\(n_i\\), the number of scores in category i
    * \\(r_{i,j}\\), the rank of the j-th score in category i, using all original scores (incl. those not in the comparison).
    * \\(\\bar{r}_i\\), the average of the ranks in category i, using all original scores (incl. those not in the comparison).
    * \\(Q\\left(\\dots\\right)\\), the cumulative distribution function of the standardized range distribution.
    * \\(\\chi^2\\left(\\dots\\right)\\), the cumulative distribution function of the chi-square distribution.
    
    References
    ----------
    BrightStat. (n.d.). Kruskal-Wallis test. BrightStat. Retrieved October 25, 2023, from https://secure.brightstat.com/index.php?p=c&d=1&c=2&i=7
    
    Nemenyi, P. (1963). *Distribution-free Multiple Comparisons*. Princeton University.
    
    Pohlert, T. (2016). The pairwise multiple comparison of mean ranks package (PMCMR). https://cran.r-hub.io/web/packages/PMCMR/vignettes/PMCMR.pdf
    
    Sachs, L. (1982). *Applied statistics: A handbook of techniques*. Springer-Verlag.
    
    Schaich, E., & Hamerle, A. (1984). *Verteilungsfreie statistische Prüfverfahren*. Springer. doi:10.1007/978-3-642-70032-3
    
    Zaiontz, C. (n.d.). Schaich-Hamerle Test after KW. Real Statistics Using Excel. Retrieved October 25, 2023, from https://real-statistics.com/one-way-analysis-of-variance-anova/kruskal-wallis-test/schaich-hamerle-test/
    
    Author
    ------
    Made by P. Stikker
    
    Companion website: https://PeterStatistics.com  
    YouTube channel: https://www.youtube.com/stikpet  
    Donations: https://www.patreon.com/bePatron?u=19398076
    
    '''
    
    #set preset
    tiesPresent = False
    nSame = True
    
    #create the cross table    
    ct = tab_cross(ordField, catField, order1=levels, order2=categories, totals="include")
    
    #basic counts
    k = ct.shape[1]-1
    nlvl = ct.shape[0]-1
    n = ct.iloc[nlvl, k]
    
    #ties
    t = 0
    n = 0
    for i in range(0,nlvl):
        n = n + ct.iloc[i, k]
        tf = ct.iloc[i, k]
        t = t + tf**3 - tf
    t = t/(n**3 - n)

    #the ranks of the levels
    lvlRank = pd.Series(dtype="object")
    cf = 0
    for i in range(0,nlvl):
        lvlRank.at[i] = (2 * cf + ct.iloc[i, k] + 1) / 2
        cf = cf + ct.iloc[i, k]
    
    #sum of ranks per category
    nEach = ct.iloc[nlvl, 0]
    srj = pd.Series(dtype="object")
    for j in range(0,k):
        sr = 0
        for i in range(0,nlvl):
            sr = sr + ct.iloc[i, j] * lvlRank.iloc[i]
            
            #check for ties
            if ct.iloc[i,j] > 1:
                tiesPresent=True
                
        srj.at[j] = sr
        
        #check if sample sizes are the same
        if nEach != ct.iloc[nlvl, j]:
            nSame = False
    
    if version == "auto":
        if not(tiesPresent) and nSame:
            version = "exact"
        elif not(tiesPresent) and nSame:
            version = "sh"
        else:
            version = "sh-ties"
    
    if version == "exact":
        ff = n * (n + 1) / 24
        df = n - k
    else:
        ff = n * (n + 1) / 12
        df = k - 1
    
    ncomp = k * (k - 1) / 2
    res = pd.DataFrame()
    resRow = 0
    for i in range(0, k-1):
        for j in range(i+1, k):
            n1 = ct.iloc[nlvl, i]
            n2 = ct.iloc[nlvl, j]
            m1 = srj.iloc[i] / n1
            m2 = srj.iloc[j] / n2
            d = m1 - m2
            
            if version == "exact":
                Var = ff * (1 / n1 + 1 / n2)
            elif version == "sh-ties":
                Var = (1 - t) * ff * (1 / n1 + 1 / n2)
            else:
                Var = ff * (1 / n1 + 1 / n2)
            
            se = (Var)**0.5
            
            if version == "exact":
                stat = d / se
                p = 1 - studentized_range.cdf(abs(stat), k, 100001)
            else:
                stat = (d / se)**2
                p = chi2.sf(stat, df)
            
            res.at[resRow, 0] = ct.columns[i]
            res.at[resRow, 1] = ct.columns[j]
            res.at[resRow, 2] = n1
            res.at[resRow, 3] = n2
            res.at[resRow, 4] = m1
            res.at[resRow, 5] = m2
            res.at[resRow, 6] = se
            res.at[resRow, 7] = stat
            res.at[resRow, 8] = p
            if res.iloc[resRow, 8] > 1:
                res.at[resRow, 8] = 1            
            
            resRow = resRow + 1
    
    colNames = ["cat. 1", "cat. 2", "n1", "n2", "mean rank 1", "mean rank 2", "statistic", "se", "p-value"]
    res.columns = colNames
    return res