I guess I have to retract my concern regarding the use of the F-statistic in
the randomization/permutation test. The code below demonstrates concurrently
testing three statistics when normality & heterogeneity is violated. The
statistics are: F, sum-squared deviation of each group mean from the mean
group mean ("Sm") , sum squared deviation of the group variance from the
mean group variance ("Sv"). Contrary to my apparently naive expectation, the
F and Sm statistics yield equivalent results any time you run the code. I
included Sv as a demonstration of a statistic that is not equivalent to
F/Sm, yet one that may be of interest as it directly tests the homogeneity
of variance assumption.

This is likely a very rudimentary demonstration to randomization pros, but I
include it in this thread for posterity.

get_Sm = function(dv,iv){
    means = aggregate(dv,list(iv),mean)
    Sm = sum( (means$x - mean(means$x))^2 )
    return(Sm)
}
get_Sv = function(dv,iv){
    vars = aggregate(dv,list(iv),var)
    Sv = sum( (vars$x - mean(vars$x))^2 )
    return(Sv)
}

dv = c( rnorm(10) , rnorm(10,-.25,.5) , rnorm(10)+rexp(10,1/.25) )
iv = factor(rep(c('a','b','c'),each=10))

obs.aov = aov(dv~iv)
obs.F = summary(obs.aov)[[1]]$F[1]
obs.Sm = get_Sm(dv,iv)
obs.Sv = get_Sv(dv,iv)

perms = 1e3
perm.F = rep(NA,perms)
perm.Sm = rep(NA,perms)
perm.Sv = rep(NA,perms)
for(i in 1:perms){
    reshuffled.dv = dv[order(runif(length(dv)))]
    perm.aov = aov(reshuffled.dv~iv)
    perm.F[i] = summary(perm.aov)[[1]]$F[1]
    perm.Sm[i] = get_Sm(reshuffled.dv,iv)
    perm.Sv[i] = get_Sv(reshuffled.dv,iv)
}

p.F = mean(perm.F>obs.F)
p.Sm = mean(perm.Sm>obs.Sm)
p.Sv = mean(perm.Sv>obs.Sv)

print(p.F)
print(p.Sm)
print(p.Sv)



-- 

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]]

______________________________________________
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.

Reply via email to