[R] Off topic: SS formulae for 3-way repeated measure anova (for when aov() fails)

[Mike Lawrence]

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