[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