# How to Properly Interpret the Hypotheses/Results of a Permutation Test Based on the Test Statistic?

I asked a question recently: Is This an Appropriate Application of a Permutation Test?. All information about my data, the permutation test and the results can be found there, but here is a quick summary:

• I had independent count observations, both within and between groups
• Each observation could take on a value between $$0-36$$
• Each of my 3 groups had a different number of observations $$(n_1=11, n_2=5,n_3=6)$$
• I wanted to test whether or not there was a significant difference in the means of these 3 groups
• I ran 5000 iterations of a permutation test using the observed F-statistic as my test statistic
• I obtained a p-value of $$0.131$$ (ie. the probability of obtaining an F-value from permuted data larger than the F-value I observed in my original data assuming $$H_0$$ to be true) and did not reject the null hypothesis that all group means were the same.

@Thomas Lumley Gave an excellent answer to my question and corrected my misinterpretation of the null hypothesis.

As was noted in his answer:

"Strictly speaking, the null hypothesis is that the distributions are the same, not just that they have the same means. (If they had the same means but difference variances, the test would have the wrong Type I error rate.)"

To my understanding, the F-statistic describes the variation among group means. A large F-statistic would imply that the variation among group means is larger than one might expect to see by chance. Could someone please help me understand what it is I'm not considering? How is it that this permutation test provides an overall analysis on the differences of the distributions of the populations from which we observe these 3 groups and not just their means?

Data. You have data as follows:

x1 = c(9,9,8, 11,18,4, 0,0,8, 12,12)
x2 = c(8,8,7, 0,0,4)
x3 = c(5,36,10, 10,6)


Group sample means and variances are given below (from R):

mean(x1);  sd(x1)
[1] 8.272727
[1] 5.349596
mean(x2);  sd(x2)
[1] 4.5
[1] 3.781534
mean(x3);  sd(x3)
[1] 13.4
[1] 12.83745


Possible heteroscedasticity. The last sample, which has the smallest size, is markedly more variable than the others. In these circumstances it is difficult to know, as you have mentioned, whether the F-statistic in a standard ANOVA is measuring differences among means or differences among variances. Moreover, you mentioned in discussion at your previous (linked) question that there may be extraordinary numbers of $$0$$'s in your data, which could also contribute to unequal group variances.

ANOVA without equal variance assumption. However, R implements a version oneway.test of the one-way ANOVA that does not assume equal population variances among groups. [This is accomplished by using Satterthwaite's corrected degrees of freedom for the F-statistic, somewhat as in the Welch two-sample t test, which also does not assume equal variances.]

Because sample variances seem to vary among groups, I would use the Satterthwaite F-statistic of oneway test as the metric for a permutation test. Because the degrees of freedom among F-statistics changes from one re-sample to another, it is best to use the P-values of the F-statistics for the permutation test.

The observed P-value from oneway.test is obtained as follows:

x = c(x1,x2,x3);  g = rep(1:3, c(11,6,5)
oneway.test(x ~ g)

One-way analysis of means (not assuming equal variances)

data:  x and g
F = 2.0078, num df = 2.0000, denom df = 8.5842, p-value = 0.1926

pv.obs = oneway.test(x ~ g)$p.val; pv.obs [1] 0.1925989  Permutation test using this metric. Below is a permutation test based on oneway.test. The permutation takes place by using sample(g) to put 23 observations at random into the three groups. [Without extra parameters, sample(g) randomly permutes the elements of g.] We obtain P-values for $$10\,000$$ permutations of your data. set.seed(811) pv.re = replicate(10^4, oneway.test(x~sample(g))$p.val)
mean(pv.re <= pv.obs)
[1] 0.6631


Taken at face value, this seems a very large P-value. Often a permutation test for a particular metric closely imitates the the observed P-value of a single test based on that metric, confirming that assumptions for the single test are not badly awry. The oneway.test assumes normal data.

Residuals not normal. So it seems worthwhile looking at the normality of the 21 residuals for your data. They are found by subtracting group means from each observation in the group. We subject them to a Shapiro-Wilk test, finding strong rejection of the null hypothesis of normality (P-value 0.005)---especially strong for such a small amount of data. An important difficulty is the outlying value $$36$$ in the third group.

r1 = x1 - mean(x1);  r2 = x2 - mean(x2);  r3 = x3 - mean(x3)
r = c(r1, r2, r3)
shapiro.test(r)

Shapiro-Wilk normality test

data:  r
W = 0.85958, p-value = 0.005038


An important departure from normality of residuals arises from the outlying value $$36$$ in the third group. It accounts for the point at upper-right in the normal probability plot below.

qqnorm(r);  qqline(r)


Along with a graph at the end of his Answer on the original (linked) page, @Thomas Lumley alluded in passing to the nonnormality of your data.

Stochastic domination. As a final look at your data, it may be worth noting is that Group 2 is stochastically dominated by Groups 1 and 3, as seen in the plot of their three ECDFs below. [The red ECDF lies to the left of the other two. Because of the many ties, it is not feasible to confirm what we see with formal Kolmogorov-Smirnov tests.]

plot(ecdf(x1), col="blue", main=hdr)
lines(ecdf(x2), col="red")
lines(ecdf(x3), col="darkgreen")


Satterthwaite one-way ANOVA with normal data. If data of the same sizes as yours arise from normal distributions with noticeably different means and variances, oneway.test behaves as intended. It has no problem detecting differences among the means in spite of differences among the variances. In the particular case shown below, the P-value is 0.0014.

set.seed(1234)
y1 = rnorm(11, 18, 1)
y2 = rnorm( 6, 14, 2)
y3 = rnorm( 5, 13, 4)
y = c(y1,y2,y3)
pv.oby = oneway.test(y ~ g)$p.val ; pv.oby [1] 0.001396934  Moreover, a permutation test based on oneway.test gets about the same result, indicating that the data and test are well-matched. pv.ry = replicate(10^4, oneway.test(y~sample(g))$p.val)
mean(pv.ry <= pv.oby)
[1] 0.0015


Things to try. It is difficult to give definitive reasons for the behavior of various tests with your data because there are no rejections. So we can't ponder why one test rejects and another doesn't.

(1) Perhaps you should try checking for heteroscedasticity using Levine's test which does not assume normal data.

(2) ECDF plots seem to show that Gp 2 is stochastically dominated by both of the other two groups. However, a Kruskal-Wallis test does not reject, and a permutation test using the K-W test gave about the same (non-significant) P-value as the test itself. My guess is that sample sizes are too small for these tests to have reasonable power. If you have similar datasets of larger size, it might be worthwhile investigating differences between tests that assume normality and those that don't.

kruskal.test(x ~ g)$$p.val [1] 0.135777 wilcox.test(x2,x3)$$p.val
[1] 0.1181535  # With warning message about ties

• This was an amazing answer! Thank you very much, you've addressed many of the concerns I had both on this post and my other linked question. Aug 12, 2020 at 23:31