[R] Function aovp in library lmPerm
John Kolassa
kolassa at stat.rutgers.edu
Tue Sep 10 20:46:50 CEST 2013
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
More information about the R-help
mailing list