[R] Learning to do randomized block design analysis

[Ben Bolker]

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.