[R] Rép : ANOVA and permutation tests : beware of traps
Stéphane Adamowicz
stephane.adamowicz at avignon.inra.fr
Wed Sep 24 14:00:52 CEST 2014
Many thanks to J. Wiley and B. Ripley for their quick answers.
I suppose that many end users are aware of problems in calculation accuracy with computers. However, I would like to comment that it was not that obvious for me that the data order matters. First, I do not find any clear mention of this particular problem in the == help page, but perhaps I am experiencing difficulties with the English. Second, I do not encounter this problem neither with the piece of code I proposed to replace the dubious one, or in the following experiments :
> # experiment 1 : comparing total variances
> var(dat1$Y) == var(dat2$Y)
[1] TRUE
> # experiment 2 : comparing bilateral T tests
> abs(t.test(Y~F, dat1)$statistic) == abs(t.test(Y~F, dat2)$statistic)
t
TRUE
> # experiment 3 : applying permutations to T tests
> nperm <- 10000
> T <- abs(t.test(Y~F, dat1)$statistic)
> Tperm <- replicate(n=nperm, abs(t.test(sample(Y)~F, dat1)$statistic))
[1] 0.1018898 # that's nice !
Thus, why a naive end user as I am should expect such pitfalls with F values given by the lm function ? Furthermore, codes similar to the one I criticized can be found in the teaching documents of various Universities and thus are spreading out. I would not be surprised that some scientific papers already rely on it ...
In fact, even in R web pages, under Books ("This page gives a partially annotated list of books that are related to S or R and may be useful to the R user community"), I found only one book clearly devoted to randomization methods : "[32] Laura Chihara and Tim Hesterberg. Mathematical Statistics with Resampling and R. Wiley, 1st edition, 2011. ISBN 978-1-1180-2985-5". Looking at the author's profiles, I would say that "Beware of the trap of listening to people with no knowledge of basic numerical methods!" does not apply to them. Here is their recommended R code for one-way anova (chapter 12, but adapted to my example data):
> observed <- anova(lm(Y ~ F, dat1))$F[1]
> n <- length(dat1$Y)
> B <- 10^5 - 1
> results <- numeric(B)
> for (i in 1:B)
+ {
+ index <- sample(n)
+ Y.perm <- dat1$Y[index]
+ results[i] <- anova(lm(Y.perm ~ dat1$F))$F[1]
+ }
> (sum(results> observed) + 1) / (B + 1) # P value
[1] 0.03969
Well, it seems that I am not the only guy who does not find the trap obvious ...
Thus, my final question remains : How can we evaluate the reliability of CRAN packages that propose randomization (or bootstrap) methods ?
Cheers, Stéphane
_______________________________________________________
Stéphane Adamowicz
stephane.adamowicz at avignon.inra.fr
_______________________________________________________
UR 1115 Plantes et Systèmes de Culture Horticoles (PSH)
228, route de l'aérodrome
CS 40509
domaine St Paul, site agroparc
84914 Avignon, cedex 9
France
Tel +33 (0)4 32 72 24 35
Fax +33 (0)4 32 72 24 32
http://www.inra.fr/
https://www6.paca.inra.fr/psh
_______________________________________________________
More information about the R-help
mailing list