The code below may not be the most optimized, but it should do the trick.
I've implemented the permutation test, which by some nomenclatures is
distinguished from the randomization test in that the latter is exhaustive
while the former is not.
I'm eager to hear the views of the list on whether the permutation test
fares when you use a parametric stat like the F as the focal statistic. I
know that David Howell introduction to the test employs the F (
http://www.uvm.edu/~dhowell/StatPages/Resampling/Resampling.html). I
personally think that in cases of known violation of normality and/or
homogeneity of variance, it seems odd to incorporate MSerror into the focal
statistic. I look forward to the list educating me if this view turns out to
be uninformed.
#make sure the factors are treated as such
FactorA = factor(FactorA)
FactorB = factor(FactorB)
#get the observed Fs
obs.aov = aov(Resp~FactorA*FactorB)
obs.mainA.F = summary(obs.aov)[[1]]$F[1]
obs.mainB.F = summary(obs.aov)[[1]]$F[2]
obs.AbyB.F = summary(obs.aov)[[1]]$F[3]
#prepare for the permutation loop
perms = 1e4
perm.mainA.F = rep(NA,perms)
perm.mainB.F = rep(NA,perms)
perm.AbyB.F = rep(NA,perms)
#run the permutation loop
for(i in 1:perms){
Randomized.Resp = Resp[order(runif(length(Resp)))]
perm.aov = aov(Randomized.Resp~FactorA*FactorB)
perm.mainA.F[i] = summary(perm.aov)[[1]]$F[1]
perm.mainB.F[i] = summary(perm.aov)[[1]]$F[2]
perm.AbyB.F[i] = summary(perm.aov)[[1]]$F[3]
}
#compute the p-values
p.mainA = mean(perm.mainA.F>obs.mainA.F)
p.mainB = mean(perm.mainB.F>obs.mainB.F)
p.AbyB = mean(perm.AbyB.F>obs.AbyB.F)
#print the p-values
print(p.mainA)
print(p.mainB)
print(p.AbyB)
--
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University
www.thatmike.com
Looking to arrange a meeting? Do so at:
http://www.timetomeet.info/with/mike/
~ Certainty is folly... I think. ~
[[alternative HTML version deleted]]