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