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

From: Mike Lawrence <mla_at_dal.ca>
Date: Sun, 20 Jul 2008 16:03:56 -0300

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
Department of Psychology
Dalhousie University
www.thatmike.com

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help