My grasp of R and statistics are both seriously lacking, so if this question
is completely naive, I apologize in advance. I've hunted for a couple hours
on the internet and none of the methods I've found have produced the result
I'm looking for.

I'm currently a student in a Statistics class and we are learning the ANOVA.
We had to do one by hand and then reproduce our work in SAS. I really like
the idea of understanding R however and would like to reproduce the solution
in R if possible.

Where I'm at now is this little program:
http://odin.himinbi.org/classes/psy304b/spider_analysis.r

The program calculates an anova manually (correctly, I'm pretty sure, it
agrees with the same numbers in excel). The answer that it comes up with
doesn't agree with any of the numbers I can get either the aov or anova
functions to produce.

Can anyone help me with simply the method to compute a one-way anova? Well,
specifically how to replicate the sort of anova people learn in an intro to
statistics class. All of the degrees of freedom are off from what I expect
them to be (they're all 1).

(The original problem, should it help in understanding my question, is at:
http://odin.himinbi.org/classes/psy304b/homework_1.xhtml#2 though it will
likely look pretty funky if your browser doesn't support mathml (firefox
does).)

Will

The program is as follows:

library(foreign)
# spiderdata <- read.csv("spider_data.csv")

spiderdata = data.frame(Zoloft = c(9, 11, 5, 12, 15, 14, 13, 12, 7, 6),
Naltrexone = c(15, 16, 12, 12, 18, 19, 23, 20, 13, 17),
Valium = c(9, 11, 12, 5, 13, 15, 11, 8, 6, 9))

summary(spiderdata)

# Compute a one-way ANOVA by hand

J = length(spiderdata)

sqdata <- data.frame((spiderdata[1] - mean(spiderdata[1])) ^ 2)
for(j in 2:J) {
sqdata <- cbind(sqdata, (spiderdata[j] - mean(spiderdata[j])) ^ 2)
}
sqdata

N = 0
for(j in 1:J) {
N = N + length(sqdata[[j]])
}

SSW = sum(sqdata)
MSW = SSW / (N - J)
SSB = 0
for(j in 1:(length(spiderdata))) {
SSB = SSB + length(spiderdata[[j]]) * ((mean(spiderdata[j])[[1]] -
(sum(spiderdata) / N)) ^ 2)
}
MSB = SSB / (J - 1)

F = MSB / MSW
f_prob = pf(F, J - 1, N - J)
reject_point = qf(.95, J - 1, N - J)

cat("SSW:", SSW, ", MSW:", MSW, ", SSB:", SSB, ", MSB:", MSB, ", F:", F, ",
P(F):", f_prob, ", P(", reject_point, ") = .95\n", sep = "")

anova(lm(Zoloft ~ Valium + Naltrexone, data = spiderdata))
aov(Zoloft ~ Valium + Naltrexone, data = spiderdata)

	[[alternative HTML version deleted]]

