Dear Colleagues,
I'm attempting to use the function aovp in library lmPerm, and am getting
strange results. The package maintainer appears to be deceased, and I'm
wondering if anyone else has experience with this package.
I'm attempting a permutation test for a two-way analysis of variance
model. An
example in the aovp documentation for the lmPerm library uses the following
code:
library("lmPerm")
data(Hald17.4)
summary(aovp(Y~T+Error(block),Hald17.4))
to give a MC p-value for T of approximately .02. I compare this with the
normal theory results of
summary(aov(Y~T+Error(block),Hald17.4))
to give a p-value of 0.0221. So far so good. I then try a new example, using
data
from Higgins, Introduction to Nonparametric Statistics, p. 142:
hay<-as.data.frame(list(y=c(
1.5,2.1,1.9,2.8,1.4,1.8, 1.8,2.0,2.0,2.7,1.6,2.3,
1.9,2.5,2.5,2.6,2.1,2.4),
day=as.factor(c(rep("1",6),rep("15",6),rep("30",6))),
block=as.factor(rep(1:6,3))))
summary(aov(y~Error(block)+day,data=hay))
gives a MC p-value of about .2. Increasing the MC sample size as
summary(aovp(y~Error(block)+day,data=hay,Ca=.00001,maxIter=100000))
confirms that the p-value is .20, to two decimal paces. But the normal-theory
p-value, determined by
summary(aov(y~Error(block)+day,data=hay))
is 0.0129. Furthermore, I can calculate the permutation p-value exactly, by
examining all (3!)^6 permutations, and I obtain the p-value 0.0218. To
complicate
matters further, reordering the data by block, and rerunning the permutation
anova,
hay<-hay[order(hay$block),]
summary(aovp(y~Error(block)+day,data=hay,Ca=.00001,maxIter=100000))
gives a p-value .046.
Does anyone with experience with lmPerm have any suggestions for resolving what
looks like contradictory results? Thanks, John
______________________________________________
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.
