[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



this_aov = 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)

#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)

#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)

#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)

#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)

#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)

#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)

#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

R-help_at_r-project.org mailing list
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Sun 20 Jul 2008 - 19:08:34 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 20 Jul 2008 - 19:32:03 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive