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.adamow...@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
_______________________________________________________
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
