'one-way anova' randomization test in R

Question

I want to implement basic randomization tests as described in Ernst, 2004 Permutation methods for exact inference. I have observations of a single continuous variable for n individuals falling into k groups. I want to know if differences in mean score per group are greater than would be expected to occur by chance.

Consider the following data (n=106 observations, k=8 groups):

score<-c(71,57,52,44,67,53,40,59,46,64,58,57,54,74,48,54,57,59,72,59,73,66,67,70,66,44,65,57,54,41,57,50,51,51,50,58,66,69,54,63,44,73,63,56,57,58,50,64,61,56,55,65,60,58,46,64,66,62,63,71,56,53,39,66,35,46,62,39,64,65,43,61,67,54,54,56,61,67,59,56,61,45,44,69,60,53,64,51,61,68,60,59,53,63,52,65,68,59,71,63,62,64,61,47,50,69)
group<-c(5,2,5,6,4,6,5,5,5,3,4,5,7,3,6,5,5,7,2,6,4,4,8,5,4,6,4,5,2,6,6,5,6,6,5,5,5,5,5,2,6,5,5,6,3,6,5,5,4,4,6,3,5,6,6,5,5,4,5,4,5,6,5,5,6,5,5,3,8,2,2,4,6,6,5,6,5,5,4,8,5,3,6,5,5,4,6,5,5,6,8,5,6,5,6,6,5,4,5,5,5,5,5,4,6,5)

Following Ernst, 2004: 5.1 a test for differences between means is:

obs_means<-aggregate(values~group,FUN=mean)
group_size<-table(group)
T_stat<-sum((obs_means$values^2)*group_size) #the test statistic
#shuffle the data and recalculate
results<-numeric(10000)
for (i in 1:length(results)){
  draw<-sample(values,length(values),replace=FALSE) #no replacement!
  shuffled<-data.frame(group,draw)
  shuffled_means<-aggregate(draw~group,data=shuffled,FUN=mean)
  shuffled_T_stat<-sum((shuffled_means$draw^2)*group_size)
  results[i]<-shuffled_T_stat
}
crit_val<-quantile(results,0.95)
T_stat>crit_val #are there significant differences between means?

The test statistic here ("T_stat") is apparently equivalent to the F-statistic in a one-way ANOVA, so these results indicate that there is a significant difference between means at p=0.05.

Now, we want to know where that difference lies. Following Ernst, 2004: 5.2 the (approximate!) test is:

#mean scores per group
obs_means<-aggregate(values~group,FUN=mean)
#get pairwise differences
pairs<-combn(obs_means[,2],2)
pairnames<-combn(obs_means[,1],2)
diffs<-abs(pairs[1,]-pairs[2,])
names(diffs)<-paste(pairnames[1,],"vs",pairnames[2,]) # add names

#randomization
results<-numeric(10000)
for (i in 1:length(results)){
  draw<-sample(values,length(values),replace=FALSE)
  shuffled<-data.frame(group,draw)
  shuffled_means<-aggregate(draw~group,data=shuffled,FUN=mean)
  shuffled_pairs<-combn(shuffled_means[,2],2)
  shuffled_differences<-abs(shuffled_pairs[1,]-shuffled_pairs[2,])
  results[i]<-max(shuffled_differences)
}
crit_val<-quantile(results,0.95)
#print pairs with differences greater than the empirical critical value (there aren't any)
diffs[diffs>crit_val]

Assuming I've done that right, we find no differences between means when making pairwise comparisons. Since the test statistic we calculated in the first part for the overall difference between means was only slightly above the 'critical value' (the 95th percentile of the empirical distribution of test statistics that we simulated) I would chalk this up to this not being an exact test and conclude there are no differences between groups.

My questions are (i) is this computation and interpretation correct? (ii) is there an existing function in R that does this, perhaps better? For example, I would use the XNomial package if these data were categorical.

Answer 1

After reading the section on anova, my guess is that the T statistic is a quick approximation for the F statistic. It does not account for explained variance so I guess if the groups have different spreads, then it might be problematic.

In your example it works ok. You can see how it compares with a real F statistic further into my answer.

My suggestion here is that you can compute a lot of these using built-in R functions, for example we can calculate the Tstat like this:

obs_means<-aggregate(values~group,FUN=mean)
group_size<-table(group)
T_stat<-sum((obs_means$values^2)*group_size)
> T_stat
[1] 358402

obs_T = sum(tapply(values,group,function(i)length(i)*mean(i)^2))
>obs_T
[1] 358402

We can wrap this up and calculate all at one go:

set.seed(111)
Perm_Tstat = replicate(1000,
sum(tapply(values,sample(group),function(i)length(i)*mean(i)^2))
)
group = factor(group)
Perm_Fstat = replicate(1000, anova(lm(values ~ sample(group)))[1,4])
obs_F = anova(lm(values ~ group))[1,4]

par(mfrow=c(1,2))
hist(Perm_Tstat,br=50,main="Perm Tstat");abline(v=obs_T,col="blue")
hist(Perm_Fstat,br=50,main="Perm Fstat");abline(v=obs_F,col="blue")

For the second part, according to the paper:

We determine Cα by calculating the largest difference in treatment means for each randomization and finding the 1 − α quantile of this distribution.

So the statistic is the biggest absolute difference between two pairs, effectively maximum - minimum:

Perm_diff <- replicate(1000,diff(range(tapply(values,sample(group),mean))))

quantile(Perm_diff,0.95)
 95% 
16.5

any(diffs>16.5)
[1] FALSE

And I get the same result no difference between means. This method does not take into account the variance. If you look at the data:

boxplot(values ~ group)

The spread in groups 4,6,7,8 are much smaller and these are the ones giving you the significant difference. We can try a pairwise.t.test:

pairwise.t.test(values,group)

    Pairwise comparisons using t tests with pooled SD 

data:  values and group 

  2     3     4     5     6     7    
3 1.000 -     -     -     -     -    
4 1.000 1.000 -     -     -     -    
5 1.000 1.000 1.000 -     -     -    
6 1.000 1.000 0.025 0.017 -     -    
7 1.000 1.000 1.000 1.000 1.000 -    
8 1.000 1.000 1.000 1.000 0.860 1.000

P value adjustment method: holm 

You can see there is a difference and this agrees with the anova. In conclusion I think you need to think more about whether the difference method is suitable for you are doing, and whether there's a better way to test your hypothesis

Answer 2

re (ii): the coin package offers oneway_test as general test of the equality of distributions between groups. In this case running

oneway_test(score~group,distribution=approximate(nresample=10000))

gives the Fisher-Pitman permutation test based on 10k iterations. But as alluded to by StupidWolf, there's a problem here given that the goal is to test for differences between means. The Fisher-Pitman test compares both the mean and the variance simultaneously; if either is different the test rejects the null of equal distributions, and you don't know which one it is. In these data, a significant result could just be picking up on unequal variances. See doi: 10.2466/pr0.94.1.189-194 for more.

However, we can get a test of different central tendencies using 'median_test' from the same package (which uses the same underlying function). This gives the Brown-Mood Median test, which I haven't found much documentation on. But it does test differences in medians only, so this seems like the right test to use from coin. To determine post-hoc differences, we can use the pairwiseMedianTest function from the rcompanion package. See here for an example by the package creator, and don't forget the correction for multiple testing.

Pleasingly, you get the same result (overall significant differences, located by post-hoc testing between one group and two others) with both the Fisher-Pitman and the Brown-Mood. This suggests that Ernst's 'T-statistic', if not incorrect, is at least lower power.

To work out why this might be so we need to define the relationship between the F test, Ernst's 'T statistic', and the permutation tests available in coin. If anyone feels like posting that, we're listening!