Re: [R] Learning to do randomized block design analysis

From: Ben Bolker <bolker_at_ufl.edu>
Date: Tue, 4 Dec 2007 17:16:17 -0800 (PST)

Bert Gunter wrote:
>
> Let's be careful here. aov() treats block as a **random** error component
> of
> variance. lm() treats block as a **fixed effect**. That's a different
> kettle of fish. Perhaps both Kevin and the authors of his textbook need to
> read up on fixed versus random effects and what they mean -- and what
> sorts
> of tests make sense for each.
>
>

  Note that Kevin didn't send the last message, "T.K." did. My copies of ISwR and MASS are either at work or loaned out, so please forgive confusion in the following ... but ...

# Case Study 13.2.1, page 778

cd <- c(8, 11, 9, 16, 24)
dp <- c(2, 1, 12, 11, 19)
lm <- c(-2, 0, 6, 2, 11)
table <- data.frame(Block=LETTERS[1:5],
                    "Score changes"=c(cd, dp,
                      lm), Therapy=rep(c("Contact Desensitization",
                             "Demonstration Participation",
                             "Live Modeling"), each=5)) 

model.aov <- aov(Score.changes ~ Therapy + Error(Block), data=table) summary(model.aov)

m1 = summary(model.aov)
str(m1) ## looking at the guts

## not particularly friendly! but this will do it.
## the first element of the list is the block error info
## the second is the within-block info

blockmeansq = m1[["Error: Block"]][[1]][["Mean Sq"]] errmeansq = m1[["Error: Within"]][[1]][["Mean Sq"]][2]

fval = blockmeansq/errmeansq ## 109.5/8.55 blockdf = m1[["Error: Block"]][[1]][["Df"]] ## 4 errdf = m1[["Error: Within"]][[1]][["Df"]][2] ## 8 pf(fval,df1=blockdf,df2=errdf,lower.tail=FALSE)

  In this particular case, the fixed-effect model and the RCB design will give the same p-values:

bad.aov <- aov(Score.changes ~ Therapy + Block, data=table) summary(bad.aov)

BUT THIS IS NOT TO BE RELIED UPON IN GENERAL! You can also use lme [nlme package] or lmer [lme4 package], which can do much more complicated models.

library(nlme)
m2 = lme(Score.changes ~ Therapy, random = ~1|Block, data=table) anova(m2)

detach("package:nlme")
library(lme4)
m3 = lmer(Score.changes ~ Therapy+(1|Block), data=table) anova(m3)

-- 
View this message in context: http://www.nabble.com/Learning-to-do-randomized-block-design-analysis-tf4945409.html#a14163277
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
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.
Received on Wed 05 Dec 2007 - 01:20:16 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 Wed 05 Dec 2007 - 04:30:18 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.