How to correctly apply the Nemenyi post-hoc test after the Friedman test

Question

I'm comparing the performance of multiple algorithms on multiple data sets. Since those performance measurements are not guaranteed to be normally distributed, I chose the Friedman Test with the Nemenyi post-hoc test based on Demšar (2006).

I then found another paper that, aside from suggesting other methods like the Quade test with subsequent Shaffer post-hoc test, they apply the Nemenyi test differently.

How do I apply the Nemenyi post-hoc test correctly?

1. Using the Studentized range statistic?

In Demšar's paper it says to reject the null hypothesis (no performance difference of two algorithms) if the average rank difference is greater than the critical distance CD with $$ CD = q_{\alpha}\sqrt{{k(k+1)}\over{6N}} $$

"where critical values qα are based on the Studentized range statistic divided by $\sqrt{2}.$"

After some digging I've found that you those "critical values" can be looked up for certain alphas, for example in a table for $\alpha = 0.05$, for infinite degrees of freedom (at the bottom of each table).

2. or using the normal distribution?

Just when I thought I knew what to do, I found another paper that confused me again, because they were only using the normal distribution. Demšar is stating a similar thing at page 12:

The test statistics for comparing the i-th and j-th classifier using these methods is $$ z = {{(R_i − R_j)}\over{\sqrt{{k(k +1)}\over{6N}}}} $$ The z value is used to find the corresponding probability from the table of normal distribution, which is then compared with an appropriate $\alpha$. The tests differ in the way they adjust the value of $\alpha$ to compensate for multiple comparisons.

At this paragraph he was talking about comparing all algorithms to a control algorithm, but the remark "differ in the way they adjust ... to compensate for multiple comparisons" suggests that this should also hold for the Nemenyi test.

So what seems logical to me is to calculate the p-value based on the test statistic $z$, which is normally distributed, and correct that one by dividing through $k(k-1)/2$.

However, that yields completely different rank differences at which to reject the null hypothesis. And now I'm stuck and don't know which method to apply. I'm strongly leaning towards the one using the normal distribution, because it is simpler and more logical to me. I also don't need to look up values in tables and I'm not bound to certain significance values.

Then again, I've never worked with the studentized range statistic and I don't understand it.

Answer 1

I also just started to look at this question.

As mentioned before, when we use the normal distribution to calculate p-values for each test, then these p-values do not take multiple testing into account. To correct for it and control the family-wise error rate, we need some adjustments. Bonferonni, i.e. dividing the significance level or multiplying the raw p-values by the number of tests, is only one possible correction. There are a large number of other multiple testing p-value corrections that are in many cases less conservative.

These p-value corrections do not take the specific structure of the hypothesis tests into account.

I am more familiar with the pairwise comparison of the original data instead of the rank transformed data as in Kruskal-Wallis or Friedman tests. In that case, which is the Tukey HSD test, the test statistic for the multiple comparison is distributed according to the studentized range distribution, which is the distribution for all pairwise comparisons under the assumption of independent samples. It is based on probabilities of multivariate normal distribution which could be calculated by numerical integration but are usually used from tables.

My guess, since I don't know the theory, is that the studentized range distribution can be applied to the case of rank tests in a similar way as in the Tukey HSD pairwise comparisons.

So, using (2) normal distribution plus multiple testing p-value corrections and using (1) studentized range distributions are two different ways of getting an approximate distribution of the test statistics. However, if the assumptions for the use of the studentized range distribution are satisfied, then it should provide a better approximation since it is designed for the specific problem of all pairwise comparisons.

Answer 2

As far as I know, when comparing only 2 algorithms, Demšar suggests the Wilcoxon signed rank test rather than Friedman + posthoc. I am, sadly, just as confuzed as you when it comes to decyphering what demšar's dividing by k-1 is supposed to mean.

Answer 3

I think the first approach is more correct, likewise, the Nemenyi test (also known as the Wilcoxon–Nemenyi–McDonald–Thompson test) suggested in the book "Nonparametric statistical methods" (Hollander, Wolfe, & Chicken, 2014) and implemented in the R NSM3 package (i.e., pWNMTfunction). For each $ij$ pair of ranks sum, $b$ blocks, $k$ groups/treatments:

$Q=\frac{|R_i-R_j|}{\sqrt{\frac{bk(k+1)}{12}}} \sim q_{(k, df)}$

Thus, reject $H_0$ if:

$|R_i-R_j| \geq q_{crit}\times{\sqrt{\frac{bk(k+1)}{12}}}$

The function calculates the p-value and difference between the sums of ranks. Example using the Monte Carlo (with 10000 Iterations):

ds <- data.frame(block = factor(rep(1:8, rep(3 ,8))),
                 group = factor(rep( c( "A" , "B" , "C" ), 8 )),
                 score = c(98,109,105,
                           65,69,80,
                           90,66,105,
                           90,92,97,
                           72,93,81,
                           94,126,101,
                           96,90,103,
                           100,116,105))

NSM3::pWNMT(x = ds$score,
        b = ds$block,
            trt = ds$group |> as.numeric(), n.mc=10000)
#> Number of blocks: n=8
#> Number of treatments: k=3
#> Using the Monte Carlo (with 10000 Iterations) method: 
#>  
#> For treatments 1 - 2, the Wilcoxon, Nemenyi, McDonald-Thompson R Statistic is 8. 
#> The smallest experimentwise error rate leading to rejection is 0.1476.
#>   
#> For treatments 1 - 3, the Wilcoxon, Nemenyi, McDonald-Thompson R Statistic is 10. 
#> The smallest experimentwise error rate leading to rejection is 0.0403.
#>   
#> For treatments 2 - 3, the Wilcoxon, Nemenyi, McDonald-Thompson R Statistic is 2. 
#> The smallest experimentwise error rate leading to rejection is 0.9663.
#> 

The $q_{crt}$ can be obtained via the function cRangeNor(α, k) (see note), which will provide critical values almost equal to the ones on qtukey function from the stats package. Example:

k <- 3
alpha <- .05
df <- Inf

NSM3::cRangeNor(alpha = alpha, k = k)
#> [1] 3.315
stats::qtukey(p = alpha, nmeans = k, df = df, lower.tail = F)
#> [1] 3.314493

The test statistic presented above is equivalent to the one provided by Pohlert (2014):

$Q=\frac{|\bar R_i - \bar R_j|}{\sqrt{\frac{k(k+1)}{6b}}} \times \sqrt{2}$

Thus, reject $H_0$ if:

$|\bar R_i - \bar R_j| \geq \frac{q_{crit}}{\sqrt{2}} \times {\sqrt{\frac{k(k+1)}{6b}}}$

(be aware that you did not use the mean of the ranks (i.e, $\bar R_i$) in your question, you just used ($R_i$).

The function posthoc.friedman.nemenyi.test from this PMCMR package (currently frdAllPairsNemenyiTest function in the PMCMRpluspacakge) will provide (default) Asymptotic p-values. While the function pWNMT from the NSM3 package will provide p-values using the Monte Carlo (default). In this function, you can set method = "Exact" to obtain exact p-values if permutations are $10000$ or less. The exact method can take ages to present results...

^{Note: that $df = Inf$ is assumed, some sources assume $df = b-k$.}

Answer 4

I also stumbled across the questio whether to compute the p-value from a normal or studentized t-distribution. Unfortunately, I still can't answer it, because different papers communicate different methods.

Nevertheless, for calculating adjusted p-values, you have to multiply the uncorrected p-value with the adjustment factor, e.g. p*(k-1) in case of comparisons against one control method or p*((k*(k-1))/2) for n x n comparisons.

What you should divide by the adjustment factor is the alpha value, if compared with unadjusted p's.