Module stikpetP.other.poho_friedman

Functions

def ph_friedman(data, levels=None, method='dunn', ties=True)

Post-Hoc Tests for a Friedman Test

A post-hoc test after a Friedman test can be used to determine which variables differ significantly.

This function provides three options: Dunn, Conover, and Nemenyi.

Parameters

data : dataframe

dataframe with a column for each variable

levels : dataframe or dictionary, optional

indication of what the levels are in order

method : {"dunn", "conover", "nemenyi"}, optional

Post-Hoc method. Default is "dunn"

ties : boolean, optional

apply a ties correction. Default is True

Returns

res : dataframe

with the following columns

• field 1,
• field 2,
• n, sample size
• *statistic", test statistic used
• df, degrees of freedom (if applicable)
• p-value, the p-value (significance)
• adj. p-value, Bonferroni adjusted p-value

Notes

Conover

Bartz-Beielstein et al. (2010, p. 319) attributes this to Conover (1999) (but also seen sites refering to Conover (1980), just different editions of the book) and uses as a formula: $$t_{1,2} = \frac{\left|R_1 - R_2\right|}{SE}$$ $$SE = \sqrt{\frac{2\times n\times\left(1-\frac{\chi_F^2}{n\times\left(k-1\right)}\right)\times\left(\sum_{i=1}^n\sum_{j=1}^k r_{i,j}^2 - \frac{n\times k\times\left(k+1\right)^2}{4}\right)}{\left(n-1\right)\times\left(k-1\right)}}$$

With: $$R_j = \sum_{i=1}^n r_{i,j}$$

In the original source it mentions $2\times k$ in $SE$ instead of $2\times n$, this was indeed an error (pers. comm. with Conover).

Gnambs (n.d.) and BrightStat (n.d.) show a different formula, that gives the same result, and is the one the function uses: $$SE = \sqrt{\frac{2\times\left(\sum_{i=1}^n\sum_{j=1}^k r_{i,j}^2 - \sum_{j=1}^k R_j^2\right)}{\left(n-1\right)\times\left(k-1\right)}}$$

The significance is then determined using: $$sig. = 2\times\left(1 - T\left(\left|t_{i,2}\right|, df\right)\right)$$

Note that in the calculation $SE$ is determined using all ranks, including those not in the comparison.

Nemenyi

Pohlert (2016, p. 15) shows the formula from Nemenyi (1963) as well as in Demšar (2006, pp. 11-12): $$q_{1,2} = \frac{\left|R_1 - R_2\right|}{\sqrt{\frac{k\times\left(k+1\right)}{6\times n}}}\times\sqrt{2}$$ $$df = n - k$$

This follows then a studentized range distribution with: $$sig. = 1 - Q\left(q_{1,2}, k, df\right)$$

Dunn

Benavoli et. al (2016, pp. 2-3) and IBM SPSS (2021, p. 814): $$z_{1,2} = \frac{\left|R_1 - R_2\right|}{SE}$$ $$SE = \sqrt{\frac{k\times\left(k+1\right)}{6\times n}}$$

This follows a standard normal distribution: $$sig. = 2\times\left(1 - \Phi\left(\left|z_{i,2}\right|\right)\right)$$

Bonferroni adjustment

The Bonferroni adjustment is done using: $$sig._{adj} = \min \left(sig. \times n_c, 1\right)$$ $$n_c = \frac{k\times\left(k-1\right)}{2}$$

Symbols Used

• $n$, the number of cases
• $k$, the number of variables
• $r_{i,j}$, the rank of case i, in variable j. The ranks are determined for each case.
• $\Phi\left(\dots\right)$, the standard normal cumulative distribution function.
• $Q\left(\dots\right)$, the studentized range distribution cumulative distribution function.
• $T\left(\dots\right)$, the Student t cumulative distribution function.
• $n_c$, the number of comparisons (pairs)

References

Benavoli, A., Corani, G., & Mangili, F. (2016). Should we really use post-hoc tests based on mean-ranks? Journal of Machine Learning Research, 17, 1–10. doi:10.48550/ARXIV.1505.02288

BrightStat. (n.d.). Friedman test. BrightStat. Retrieved November 5, 2023, from https://secure.brightstat.com/index.php?p=c&d=1&c=2&i=9

Conover, W. J. (1980). Practical nonparametric statistics (2nd ed.). Wiley.

Demšar, J. (2006). Statistical comparisons of classifiers over multiple data sets. The Journal of Machine Learning Research, 7, 1–30. doi:10.5555/1248547.1248548

Gnambs, T. (n.d.). SPSS Friedman. http://timo.gnambs.at/sites/default/files/spss_friedmanph.sps

IBM. (2021). IBM SPSS Statistics Algorithms. IBM.

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

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_friedman(data, levels=None, method = "dunn", ties=True):
    '''
    Post-Hoc Tests for a Friedman Test
    ----------------------------------
    A post-hoc test after a Friedman test can be used to determine which variables differ significantly.
    
    This function provides three options: Dunn, Conover, and Nemenyi.
    
    Parameters
    ----------
    data : dataframe
        dataframe with a column for each variable
    levels : dataframe or dictionary, optional
        indication of what the levels are in order
    method : {"dunn", "conover", "nemenyi"}, optional
        Post-Hoc method. Default is "dunn"
    ties : boolean, optional
        apply a ties correction. Default is True    
    
        
    Returns
    -------
    res : dataframe
        with the following columns
        
    * *field 1*, 
    * *field 2*, 
    * *n*, sample size
    * *statistic", test statistic used
    * *df*, degrees of freedom (if applicable)
    * *p-value*, the p-value (significance)
    * *adj. p-value*, Bonferroni adjusted p-value
    
    
    Notes
    -----
    
    **Conover**
    
    Bartz-Beielstein et al. (2010, p. 319) attributes this to Conover (1999) (but also seen sites refering to Conover (1980), just different editions of the book) and uses as a formula:
    $$t_{1,2} = \\frac{\\left|R_1 - R_2\\right|}{SE}$$
    $$SE = \\sqrt{\\frac{2\\times n\\times\\left(1-\\frac{\\chi_F^2}{n\\times\\left(k-1\\right)}\\right)\\times\\left(\\sum_{i=1}^n\\sum_{j=1}^k r_{i,j}^2 - \\frac{n\\times k\\times\\left(k+1\\right)^2}{4}\\right)}{\\left(n-1\\right)\\times\\left(k-1\\right)}}$$
    
    With:
    $$R_j = \\sum_{i=1}^n r_{i,j}$$
    
    In the original source it mentions \\(2\\times k\\) in \\(SE\\) instead of \\(2\\times n\\), this was indeed an error (pers. comm. with Conover).
    
    Gnambs (n.d.) and BrightStat (n.d.) show a different formula, that gives the same result, and is the one the function uses:
    $$SE = \\sqrt{\\frac{2\\times\\left(\\sum_{i=1}^n\\sum_{j=1}^k r_{i,j}^2 - \\sum_{j=1}^k R_j^2\\right)}{\\left(n-1\\right)\\times\\left(k-1\\right)}}$$
    
    The significance is then determined using:
    $$sig. = 2\\times\\left(1 - T\\left(\\left|t_{i,2}\\right|, df\\right)\\right)$$
    
    Note that in the calculation \\(SE\\) is determined using all ranks, including those not in the comparison.
    
    **Nemenyi**
    
    Pohlert (2016, p. 15) shows the formula from Nemenyi (1963) as well as in Demšar (2006, pp. 11-12):
    $$q_{1,2} = \\frac{\\left|R_1 - R_2\\right|}{\\sqrt{\\frac{k\\times\\left(k+1\\right)}{6\\times n}}}\\times\\sqrt{2}$$
    $$df = n - k$$
    
    This follows then a studentized range distribution with:
    $$sig. = 1 - Q\\left(q_{1,2}, k, df\\right)$$
    
    **Dunn**
    
    Benavoli et. al (2016, pp. 2-3) and IBM SPSS (2021, p. 814):
    $$z_{1,2} = \\frac{\\left|R_1 - R_2\\right|}{SE}$$
    $$SE = \\sqrt{\\frac{k\\times\\left(k+1\\right)}{6\\times n}}$$
    
    This follows a standard normal distribution:
    $$sig. = 2\\times\\left(1 - \\Phi\\left(\\left|z_{i,2}\\right|\\right)\\right)$$
    
    **Bonferroni adjustment**
    
    The Bonferroni adjustment is done using:
    $$sig._{adj} = \\min \\left(sig. \\times n_c, 1\\right)$$
    $$n_c = \\frac{k\\times\\left(k-1\\right)}{2}$$
    
    *Symbols Used*
    
    * \\(n\\), the number of cases
    * \\(k\\), the number of variables
    * \\(r_{i,j}\\), the rank of case i, in variable j. The ranks are determined for each case.
    * \\(\\Phi\\left(\\dots\\right)\\), the standard normal cumulative distribution function.
    * \\(Q\\left(\\dots\\right)\\), the studentized range distribution cumulative distribution function.
    * \\(T\\left(\\dots\\right)\\), the Student t cumulative distribution function.
    * \\(n_c\\), the number of comparisons (pairs)
    
    References
    ----------
    Benavoli, A., Corani, G., & Mangili, F. (2016). Should we really use post-hoc tests based on mean-ranks? *Journal of Machine Learning Research, 17*, 1–10. doi:10.48550/ARXIV.1505.02288
    
    BrightStat. (n.d.). Friedman test. BrightStat. Retrieved November 5, 2023, from https://secure.brightstat.com/index.php?p=c&d=1&c=2&i=9
    
    Conover, W. J. (1980). *Practical nonparametric statistics* (2nd ed.). Wiley.
    
    Demšar, J. (2006). Statistical comparisons of classifiers over multiple data sets. *The Journal of Machine Learning Research, 7*, 1–30. doi:10.5555/1248547.1248548
    
    Gnambs, T. (n.d.). SPSS Friedman. http://timo.gnambs.at/sites/default/files/spss_friedmanph.sps
    
    IBM. (2021). IBM SPSS Statistics Algorithms. IBM.
    
    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
    
    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
    
    '''
    
    
    #Remove rows with missing values and reset index
    data = data.dropna()    
    data = data.reset_index(drop=True)        
        
    nr = len(data)
    k = len(data.columns)
    
    if levels is not None:
        data = data.replace(levels)
    
    ranks = pd.DataFrame()
    for i in range(0, nr):
        rankRow = pd.Series(rankdata(data.iloc[i, :]))
        ranks[i] = rankRow
    
    ranks = ranks.transpose()
    
    rs = ranks.sum().sum()
    rm = rs/(nr*k)
    
    #Determine for each variable the average rank, and
    #the squared difference of this average with the overall average rank.
    rmj = ranks.sum()/nr
    rs2 = (ranks**2).sum().sum()
    sst = nr*((rmj - rm)**2).sum()
    sse = ranks.stack().var(ddof=0)*k/(k-1)
    
    if ties:
        qadj = sst / sse
    else:
        qadj = 12 / (nr * k * (k + 1)) * rs2 - 3 * nc * (k + 1)
    
    ncomp = k*(k - 1)/2
    
    res = pd.DataFrame([[0]*7]*int(ncomp))
    useRow=0
    for i in range(0, k-1):
        for j in range(i+1, k):
            res.iloc[useRow,0] = data.columns[i]
            res.iloc[useRow,1] = data.columns[j]
            res.iloc[useRow,2] = nr
            useRow = useRow + 1
    
    useRow=0
    if method=="dunn":                
        se = (k * (k + 1) / (6 * nr))**0.5
        df = nr - k
        for i in range(0, k-1):
            for j in range(i+1, k):
                z = (rmj[i] - rmj[j])/se
                res.iloc[useRow,3] = z
                res.iloc[useRow,4] = None
                res.iloc[useRow,5] = 2 * (1 - NormalDist().cdf(abs(z))) 
                
                useRow = useRow + 1
                
    elif method == "conover":
        r2 = ((rmj * nr)**2).sum()        
        se = (2 * (nr * rs2 - r2) / ((nr - 1) * (k - 1)))**0.5
        df = (nr - 1) * (k - 1)
        for i in range(0, k-1):
            for j in range(i+1, k):
                tVal = (rmj[i] - rmj[j]) * nr / se
                res.iloc[useRow,3] = tVal
                res.iloc[useRow,4] = df
                res.iloc[useRow,5] = 2 * (1 - t.cdf(abs(tVal), df)) 
                
                useRow = useRow + 1
                
    elif method == "nemenyi":
        se = (k * (k + 1) / (6 * nr))**0.5
        df = nr - k
        for i in range(0, k-1):
            for j in range(i+1, k):
                q = (rmj[i] - rmj[j]) / se * (2**0.5)
                res.iloc[useRow,3] = q
                res.iloc[useRow,4] = df
                res.iloc[useRow,5] = 1 - studentized_range.cdf(abs(q), k=k, df=df, scale=1)
                
                useRow = useRow + 1
                
    useRow=0
    for i in range(0, k-1):
            for j in range(i+1, k):
                res.iloc[useRow,6] = res.iloc[useRow, 5]*ncomp
                if res.iloc[useRow,6] > 1:
                    res.iloc[useRow,6] = 1
                useRow = useRow + 1
                
    res.columns = ["field 1", "field 2", "n", "statistic", "df", "p-value", "adj. p-value"]
    
    return res