[R] difficulties computing a simple anova

Simon Blomberg s.blomberg1 at uq.edu.au
Thu Jan 31 02:11:40 CET 2008


Your data setup is wrong. You have one factor (Drug) with 3 levels
(Zoloft, Naltrexone, Valium). So your data should be:

spiderdata <- data.frame(Drug=rep(c("Zoloft", "Naltrexone", "Valium"),
each=10), Response=c(9, 11, 5, 12, 15, 14, 13, 12, 7, 6, 15, 16, 12, 12,
18, 19, 23, 20, 13, 17, 9, 11, 12, 5, 13, 15, 11, 8, 6, 9))

You should be able to take it from there. (Since this is a homework
problem, I'm being intentionally vague.)

cheers,

Simon.

On Wed, 2008-01-30 at 18:47 -0600, Will Holcomb wrote:
> 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]]
> 
> ______________________________________________
> R-help at 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.
-- 
Simon Blomberg, BSc (Hons), PhD, MAppStat. 
Lecturer and Consultant Statistician 
Faculty of Biological and Chemical Sciences 
The University of Queensland 
St. Lucia Queensland 4072 
Australia
Room 320 Goddard Building (8)
T: +61 7 3365 2506
http://www.uq.edu.au/~uqsblomb
email: S.Blomberg1_at_uq.edu.au

Policies:
1.  I will NOT analyse your data for you.
2.  Your deadline is your problem.

The combination of some data and an aching desire for 
an answer does not ensure that a reasonable answer can 
be extracted from a given body of data. - John Tukey.



More information about the R-help mailing list