Module `stikpetP.other.poho_sdcf`

Expand source code

import pandas as pd
from scipy.stats import studentized_range
from ..other.table_cross import tab_cross

def ph_sdcf(catField, ordField, categories=None, levels=None):
    '''
    Post-Hoc Steel-Dwass-Critchlow-Fligner 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. 
    
    Other post-hoc tests that could be considered are Dunn, Nemenyi, Conover, a pairwise Mann-Whitney U, or pairwise Mood-Median.
    
    Unlike the Dunn, Nemenyi and Conover-Iman test, this test re-calculates the mean ranks for each pair, using only the scores from the two categories. 
    
    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.
        
    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)
    * *statistic*, the q-value of the test
    * *std. statistic*, the standardized q value
    * *p-value*, the p-value (significance)
    
    Notes
    -----
    The formula used (Hollander & Wolfe, 1999, p. 241):
    $$q_{1,2} = \\frac{\\left| R_1 - E_1\\right|}{\\sqrt{\\sigma^2}}$$
    
    With:
    $$R_1 = \\sum_{i=1}^{n_1} r_{i,1}$$
    $$n_{1,2} = n_1 + n_2$$
    $$E_1 = \\frac{n_1\\times\\left(n_{1,2} + 1\\right)}{2}$$
    $$\\sigma^2 = \\frac{n_1\\times n_2}{12}\\times\\left(n_{1,2} + 1 - \\frac{T}{n_{1,2}-1}\\right)$$
    $$T = \\sum t_j^3 - t_j$$
    
    The p-value is then determined using (Critchlow & Fligner, 1991, p. 131):
    $$sig. = 1 - Q\\left(q_{1,2}, k, df=\\infty\\right)$$
    
    Note that while looking at the R-code for this, posted by Shigenobu (n.d.), who references Nagata and Yoshida (1997), an alternative but same result equation for the variance can be used:
    
    $$\\sigma^2 = \\frac{n_1\\times n_2}{n_{1,2}\\times\\left(n_{1,2} - 1\\right)}\\times\\left(\\sum_{i=1}^{n_1} r_{i,1}^2 + \\sum_{i=1}^{n_2} r_{i,2}^2 - \\frac{n_{1,2}\\times\\left(n_{1,2}+1\\right)^2}{4}\\right)$$.
    
    Steel (1960) and Dwass (1960) independently derived the basics for this test. Critchlow and Fligner (1991) added the case for larger samples using the Tukey Range Distribution, and in Hollander and Wolfe (1999) the version used here can be found, which includes a ties correction.
    
    *Symbols used*
    
    * \\(k\\), the number of categories
    * \\(t_j\\), the frequency of the j-th unique rank.
    * \\(n_i\\), the number of scores in category i
    * \\(r_{i,j}\\), the rank of the j-th score in category i using only the scores from the two categories in the comparison.
    * \\(Q\\left(\\dots\\right)\\), the cumulative distribution function of the standardized range distribution.
    
    References
    ----------
    Critchlow, D. E., & Fligner, M. A. (1991). On distribution-free multiple comparisons in the one-way analysis of variance. *Communications in Statistics - Theory and Methods, 20*(1), 127–139. doi:10.1080/03610929108830487
    
    Dwass, M. (1960). *Some k-sample rank-order tests*. In I. Olkin, S. G. Ghurye, W. Hoeffding, W. G. Madow, & H. B. Mann (Eds.), Contributions to probability and statistics; Essays in honor of Harold Hotelling. Stanford University Press.
    
    Hollander, M., & Wolfe, D. A. (1999). *Nonparametric statistical methods* (2nd ed.). Wiley.
    
    Nagata, Y., & Yoshida, M. (1997). *The Basics of Multiple Comparisons in Statistics*. Scientist Co.
    
    Shigenobu. (2004, July 28). Multiple comparisons using the Steel-Dwass method. http://aoki2.si.gunma-u.ac.jp/R/Steel-Dwass.html
    
    Steel, R. G. D. (1960). A rank sum test for comparing all pairs of treatments. Technometrics, 2(2), 197–207. doi:10.1080/00401706.1960.10489894
    
    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
    
    '''
    #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]
    ncomp = k * (k - 1) / 2
    res = pd.DataFrame()
    resRow = 0
    for i in range(0, k-1):
        for j in range(i+1, k):
            res.at[resRow, 0] = ct.columns[i]
            res.at[resRow, 1] = ct.columns[j]
            
            selCats = [ct.columns[i], ct.columns[j]]
            
            #create the cross table    
            ct2 = tab_cross(ordField, catField, order1=levels, order2=selCats, totals="include")
            nCat = 2
            nlvl = ct2.shape[0]-1
            
            #ties
            t = 0
            n = 0
            for i2 in range(0,nlvl):
                n = n + ct2.iloc[i2, 2]
                tf = ct2.iloc[i2, 2]
                t = t + tf**3 - tf
                
            #the ranks of the levels
            lvlRank = pd.Series(dtype="object")
            cf = 0
            for i2 in range(0,nlvl):
                lvlRank.at[i2] = (2 * cf + ct2.iloc[i2, 2] + 1) / 2
                cf = cf + ct2.iloc[i2, 2]

            #sum of ranks per category
            nEach = ct2.iloc[nlvl, 0]
            srj = pd.Series(dtype="object")
            for j2 in range(0,2):
                sr = 0
                for i2 in range(0,nlvl):
                    sr = sr + ct2.iloc[i2, j2] * lvlRank.iloc[i2]

                srj.at[j2] = sr
                
            n1 = ct2.iloc[nlvl, 0]
            n2 = ct2.iloc[nlvl, 1]
            
            e = n1 * (n + 1) / 2
            s2 = n1 * n2 / 12 * (n + 1 - t / (n * (n - 1)))
            q = (srj[0] - e) / (s2**0.5)
            
            res.at[resRow, 2] = n1
            res.at[resRow, 3] = n2
            
            res.at[resRow, 4] = srj[0]/n1
            res.at[resRow, 5] = srj[1]/n2
            
            res.at[resRow, 6] = q
            res.at[resRow, 7] = q*2**0.5
            res.at[resRow, 8] = 1 - studentized_range.cdf(abs(q*2**0.5), k, 100001)
            
            resRow = resRow + 1
            
    colNames = ["cat. 1", "cat. 2", "n1", "n2", "mean rank 1", "mean rank 2", "statistic", "std. statistic", "p-value"]
    res.columns = colNames
    
    return res

Functions

def ph_sdcf(catField, ordField, categories=None, levels=None)

Post-Hoc Steel-Dwass-Critchlow-Fligner 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.

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

Unlike the Dunn, Nemenyi and Conover-Iman test, this test re-calculates the mean ranks for each pair, using only the scores from the two categories.

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.

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)
statistic, the q-value of the test
std. statistic, the standardized q value
p-value, the p-value (significance)

Notes

The formula used (Hollander & Wolfe, 1999, p. 241): $q_{1,2} = \frac{\left| R_1 - E_1\right|}{\sqrt{\sigma^2}}$

With: $R_1 = \sum_{i=1}^{n_1} r_{i,1}$ $n_{1,2} = n_1 + n_2$ $E_1 = \frac{n_1\times\left(n_{1,2} + 1\right)}{2}$ $\sigma^2 = \frac{n_1\times n_2}{12}\times\left(n_{1,2} + 1 - \frac{T}{n_{1,2}-1}\right)$ $T = \sum t_j^3 - t_j$

The p-value is then determined using (Critchlow & Fligner, 1991, p. 131): $sig. = 1 - Q\left(q_{1,2}, k, df=\infty\right)$

Note that while looking at the R-code for this, posted by Shigenobu (n.d.), who references Nagata and Yoshida (1997), an alternative but same result equation for the variance can be used:

$\sigma^2 = \frac{n_1\times n_2}{n_{1,2}\times\left(n_{1,2} - 1\right)}\times\left(\sum_{i=1}^{n_1} r_{i,1}^2 + \sum_{i=1}^{n_2} r_{i,2}^2 - \frac{n_{1,2}\times\left(n_{1,2}+1\right)^2}{4}\right)$ .

Steel (1960) and Dwass (1960) independently derived the basics for this test. Critchlow and Fligner (1991) added the case for larger samples using the Tukey Range Distribution, and in Hollander and Wolfe (1999) the version used here can be found, which includes a ties correction.

Symbols used

$k$ , the number of categories
$t_j$ , the frequency of the j-th unique rank.
$n_i$ , the number of scores in category i
$r_{i,j}$ , the rank of the j-th score in category i using only the scores from the two categories in the comparison.
$Q\left(\dots\right)$ , the cumulative distribution function of the standardized range distribution.

References

Critchlow, D. E., & Fligner, M. A. (1991). On distribution-free multiple comparisons in the one-way analysis of variance. Communications in Statistics - Theory and Methods, 20(1), 127–139. doi:10.1080/03610929108830487

Dwass, M. (1960). Some k-sample rank-order tests. In I. Olkin, S. G. Ghurye, W. Hoeffding, W. G. Madow, & H. B. Mann (Eds.), Contributions to probability and statistics; Essays in honor of Harold Hotelling. Stanford University Press.

Hollander, M., & Wolfe, D. A. (1999). Nonparametric statistical methods (2nd ed.). Wiley.

Nagata, Y., & Yoshida, M. (1997). The Basics of Multiple Comparisons in Statistics. Scientist Co.

Shigenobu. (2004, July 28). Multiple comparisons using the Steel-Dwass method. http://aoki2.si.gunma-u.ac.jp/R/Steel-Dwass.html

Steel, R. G. D. (1960). A rank sum test for comparing all pairs of treatments. Technometrics, 2(2), 197–207. doi:10.1080/00401706.1960.10489894

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_sdcf(catField, ordField, categories=None, levels=None):
    '''
    Post-Hoc Steel-Dwass-Critchlow-Fligner 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. 
    
    Other post-hoc tests that could be considered are Dunn, Nemenyi, Conover, a pairwise Mann-Whitney U, or pairwise Mood-Median.
    
    Unlike the Dunn, Nemenyi and Conover-Iman test, this test re-calculates the mean ranks for each pair, using only the scores from the two categories. 
    
    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.
        
    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)
    * *statistic*, the q-value of the test
    * *std. statistic*, the standardized q value
    * *p-value*, the p-value (significance)
    
    Notes
    -----
    The formula used (Hollander & Wolfe, 1999, p. 241):
    $$q_{1,2} = \\frac{\\left| R_1 - E_1\\right|}{\\sqrt{\\sigma^2}}$$
    
    With:
    $$R_1 = \\sum_{i=1}^{n_1} r_{i,1}$$
    $$n_{1,2} = n_1 + n_2$$
    $$E_1 = \\frac{n_1\\times\\left(n_{1,2} + 1\\right)}{2}$$
    $$\\sigma^2 = \\frac{n_1\\times n_2}{12}\\times\\left(n_{1,2} + 1 - \\frac{T}{n_{1,2}-1}\\right)$$
    $$T = \\sum t_j^3 - t_j$$
    
    The p-value is then determined using (Critchlow & Fligner, 1991, p. 131):
    $$sig. = 1 - Q\\left(q_{1,2}, k, df=\\infty\\right)$$
    
    Note that while looking at the R-code for this, posted by Shigenobu (n.d.), who references Nagata and Yoshida (1997), an alternative but same result equation for the variance can be used:
    
    $$\\sigma^2 = \\frac{n_1\\times n_2}{n_{1,2}\\times\\left(n_{1,2} - 1\\right)}\\times\\left(\\sum_{i=1}^{n_1} r_{i,1}^2 + \\sum_{i=1}^{n_2} r_{i,2}^2 - \\frac{n_{1,2}\\times\\left(n_{1,2}+1\\right)^2}{4}\\right)$$.
    
    Steel (1960) and Dwass (1960) independently derived the basics for this test. Critchlow and Fligner (1991) added the case for larger samples using the Tukey Range Distribution, and in Hollander and Wolfe (1999) the version used here can be found, which includes a ties correction.
    
    *Symbols used*
    
    * \\(k\\), the number of categories
    * \\(t_j\\), the frequency of the j-th unique rank.
    * \\(n_i\\), the number of scores in category i
    * \\(r_{i,j}\\), the rank of the j-th score in category i using only the scores from the two categories in the comparison.
    * \\(Q\\left(\\dots\\right)\\), the cumulative distribution function of the standardized range distribution.
    
    References
    ----------
    Critchlow, D. E., & Fligner, M. A. (1991). On distribution-free multiple comparisons in the one-way analysis of variance. *Communications in Statistics - Theory and Methods, 20*(1), 127–139. doi:10.1080/03610929108830487
    
    Dwass, M. (1960). *Some k-sample rank-order tests*. In I. Olkin, S. G. Ghurye, W. Hoeffding, W. G. Madow, & H. B. Mann (Eds.), Contributions to probability and statistics; Essays in honor of Harold Hotelling. Stanford University Press.
    
    Hollander, M., & Wolfe, D. A. (1999). *Nonparametric statistical methods* (2nd ed.). Wiley.
    
    Nagata, Y., & Yoshida, M. (1997). *The Basics of Multiple Comparisons in Statistics*. Scientist Co.
    
    Shigenobu. (2004, July 28). Multiple comparisons using the Steel-Dwass method. http://aoki2.si.gunma-u.ac.jp/R/Steel-Dwass.html
    
    Steel, R. G. D. (1960). A rank sum test for comparing all pairs of treatments. Technometrics, 2(2), 197–207. doi:10.1080/00401706.1960.10489894
    
    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
    
    '''
    #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]
    ncomp = k * (k - 1) / 2
    res = pd.DataFrame()
    resRow = 0
    for i in range(0, k-1):
        for j in range(i+1, k):
            res.at[resRow, 0] = ct.columns[i]
            res.at[resRow, 1] = ct.columns[j]
            
            selCats = [ct.columns[i], ct.columns[j]]
            
            #create the cross table    
            ct2 = tab_cross(ordField, catField, order1=levels, order2=selCats, totals="include")
            nCat = 2
            nlvl = ct2.shape[0]-1
            
            #ties
            t = 0
            n = 0
            for i2 in range(0,nlvl):
                n = n + ct2.iloc[i2, 2]
                tf = ct2.iloc[i2, 2]
                t = t + tf**3 - tf
                
            #the ranks of the levels
            lvlRank = pd.Series(dtype="object")
            cf = 0
            for i2 in range(0,nlvl):
                lvlRank.at[i2] = (2 * cf + ct2.iloc[i2, 2] + 1) / 2
                cf = cf + ct2.iloc[i2, 2]

            #sum of ranks per category
            nEach = ct2.iloc[nlvl, 0]
            srj = pd.Series(dtype="object")
            for j2 in range(0,2):
                sr = 0
                for i2 in range(0,nlvl):
                    sr = sr + ct2.iloc[i2, j2] * lvlRank.iloc[i2]

                srj.at[j2] = sr
                
            n1 = ct2.iloc[nlvl, 0]
            n2 = ct2.iloc[nlvl, 1]
            
            e = n1 * (n + 1) / 2
            s2 = n1 * n2 / 12 * (n + 1 - t / (n * (n - 1)))
            q = (srj[0] - e) / (s2**0.5)
            
            res.at[resRow, 2] = n1
            res.at[resRow, 3] = n2
            
            res.at[resRow, 4] = srj[0]/n1
            res.at[resRow, 5] = srj[1]/n2
            
            res.at[resRow, 6] = q
            res.at[resRow, 7] = q*2**0.5
            res.at[resRow, 8] = 1 - studentized_range.cdf(abs(q*2**0.5), k, 100001)
            
            resRow = resRow + 1
            
    colNames = ["cat. 1", "cat. 2", "n1", "n2", "mean rank 1", "mean rank 2", "statistic", "std. statistic", "p-value"]
    res.columns = colNames
    
    return res