Skip to content

Learning to do randomized block design analysis

7 messages · Zembower, Kevin, T.K., Bert Gunter +1 more

Original
Zembower, Kevin
# 
We just studied randomized block design analysis in my statistics class,
and I'm trying to learn how to do them in R. I'm trying to duplicate a
case study example from my textbook [1]:
lm), Therapy=rep(c("Contact Desensitisztion", "Demonstration
Participation", "Live Modeling"), each=5))
Block Score.changes                     Therapy
1      A             8     Contact Desensitisztion
2      B            11     Contact Desensitisztion
3      C             9     Contact Desensitisztion
4      D            16     Contact Desensitisztion
5      E            24     Contact Desensitisztion
6      A             2 Demonstration Participation
7      B             1 Demonstration Participation
8      C            12 Demonstration Participation
9      D            11 Demonstration Participation
10     E            19 Demonstration Participation
11     A            -2               Live Modeling
12     B             0               Live Modeling
13     C             6               Live Modeling
14     D             2               Live Modeling
15     E            11               Live Modeling
Error: Block
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4  438.0   109.5               

Error: Within
          Df Sum Sq Mean Sq F value   Pr(>F)   
Therapy    2 260.93  130.47  15.259 0.001861 **
Residuals  8  68.40    8.55                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
I don't understand why R doesn't output a value for F and Pr for the
Error (Block) dimension, as my textbook shows 12.807 and 0.0015
respectively. All the other numbers match. Can these two values be
recovered? Also, my text shows a total line which R omits. Is this
because it's not particularly useful?

Thanks for your suggestions and advice. Also, if I'm executing this type
of problem in R inefficiently, I'd appreciate suggestions.

-Kevin

[1] An Introduction to Mathematical Statistics and Its Applications,
Larsen and Marx, fourth edition.

Kevin Zembower
Internet Services Group manager
Center for Communication Programs
Bloomberg School of Public Health
Johns Hopkins University
111 Market Place, Suite 310
Baltimore, Maryland  21202
410-659-6139
Bert Gunter
# 
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.


Bert Gunter
Genentech Nonclinical Statistics


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of T.K.
Sent: Tuesday, December 04, 2007 2:46 PM
To: Zembower, Kevin
Cc: r-help at r-project.org
Subject: Re: [R] Learning to do randomized block design analysis

This seems to work.
The trick is to use 'lm' instead of 'aov'.
data=table)
Analysis of Variance Table

Response: Score.changes
                Df Sum Sq Mean Sq F value   Pr(>F)
factor(Therapy)  2 260.93  130.47  15.259 0.001861 **
factor(Block)    4 438.00  109.50  12.807 0.001484 **
Residuals        8  68.40    8.55
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1
On Dec 4, 2007 12:15 PM, Zembower, Kevin <kzembowe at jhuccp.org> wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
Bert Gunter wrote:

            

      
      
      
    

        
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)