Mike Lawrence
2008-Jul-20 19:03 UTC
[R] Off topic: SS formulae for 3-way repeated measure anova (for when aov() fails)
Pursuant to a prior "on topic" thread (http://tolstoy.newcastle.edu.au/R/e4/help/08/07/17192.html ) where I found I could not use AOV to perform an anova on my large data set, I'm now trying to code the analysis "by hand" so to speak. However, as demonstrated below, when comparing my attempt to aov() using a smaller data set, I seem to betray some sort of misunderstanding when I try to compute SSerr for the first interaction. Obviously I have missed something and although I have looked around for the explicit SSerr formulas for this design (my work thus far was extrapolated from understanding of a 2-way design), I can't seem to find any. Any help would be much obliged. N = 20 levs.a = 2 levs.b = 2 levs.d = 10 temp.sub = factor(1:N) temp.a = factor(1:levs.a) temp.b = factor(1:levs.b) temp.d = factor(1:levs.d) temp = expand.grid(sub=temp.sub, a=temp.a, b=temp.b, d=temp.d) temp$x = rnorm(length(temp[, 1])) #generate some random DV data sub=temp$sub a=temp$a b=temp$b d=temp$d x=temp$x this_aov = aov( x~a*b*d+Error(sub/(a*b*d)) ) summary(this_aov) #now let's try by hand, checking each sum-of-squares # ss against the analogous aov() produced ss (rounding # each to avoid small computational differences) #Get ss.subs sub.means = aggregate(x,list(sub=sub), mean) grand.mean = mean(sub.means$x) ss.total = sum((x-grand.mean)^2) ss.subs = levs.a*levs.b*levs.d*sum((sub.means$x-grand.mean)^2) round(ss.subs, 10)==round(summary(this_aov)[[1]][[1]]$Sum, 10) #Get ss.a a.means = aggregate(x, list(a=a), mean) ss.a = N*levs.b*levs.d*sum((a.means$x-grand.mean)^2) round(ss.a, 10)==round(summary(this_aov)[[2]][[1]]$Sum[1], 10) #ok! #Get ss.a.error a.cells = aggregate(x, list(a=a, sub=sub), mean) ss.a.cells = levs.b*levs.d*sum((a.cells$x-grand.mean)^2) ss.a.error = ss.a.cells - ss.a - ss.subs round(ss.a.error, 10)==round(summary(this_aov)[[2]][[1]]$Sum[2], 10) #ok! #Get ss.b b.means = aggregate(x, list(b=b), mean) ss.b = N*levs.a*levs.d*sum((b.means$x-grand.mean)^2) round(ss.b, 10)==round(summary(this_aov)[[3]][[1]]$Sum[1], 10) #ok! #Get ss.b.error b.cells = aggregate(x, list(b=b, sub=sub), mean) ss.b.cells = levs.a*levs.d*sum((b.cells$x-grand.mean)^2) ss.b.error = ss.b.cells - ss.b - ss.subs round(ss.b.error, 10)==round(summary(this_aov)[[3]][[1]]$Sum[2], 10) #ok! #Get ss.d d.means = aggregate(x, list(d=d), mean) ss.d = N*levs.a*levs.b*sum((d.means$x-grand.mean)^2) round(ss.d, 10)==round(summary(this_aov)[[4]][[1]]$Sum[1], 10) #ok! #Get ss.d.error d.cells = aggregate(x, list(d=d, sub=sub), mean) ss.d.cells = levs.a*levs.b*sum((d.cells$x-grand.mean)^2) ss.d.error = ss.d.cells - ss.d - ss.subs round(ss.d.error, 10)==round(summary(this_aov)[[4]][[1]]$Sum[2], 10) #ok! #Get ss.aBYb aBYb.means = aggregate(x, list(a=a, b=b), mean) ss.aBYb = N*levs.d*sum((aBYb.means$x-grand.mean)^2) - ss.a - ss.b round(ss.aBYb, 10)==round(summary(this_aov)[[5]][[1]]$Sum[1], 10) #ok! #Get ss.aBYb.error aBYb.cells = aggregate(x, list(a=a, b=b, sub=sub), mean) ss.aBYb.cells = levs.d*sum((aBYb.cells$x-grand.mean)^2) ss.aBYb.error = ss.aBYb.cells - ss.aBYb - ss.subs round(ss.aBYb.error, 10)==round(summary(this_aov)[[5]][[1]]$Sum[2], 10) #not ok :( -- Mike Lawrence Graduate Student Department of Psychology Dalhousie University www.thatmike.com