Randomized Complete Block Design

Dioxin Example

> dioxin <− data.frame(dioxin=scan(),

+           region=factor(rep(1:3,times=4)),

+           lab=factor(rep(1:4,each=3)))

1: 1.4 1.2 0.5 2.6 2.5 1.9 1.7 1.5 1.1 1.1 1.1 0.2

13:

Read 12 items

> dioxin

   dioxin region lab
1     1.4      1   1
2     1.2      2   1
3     0.5      3   1
4     2.6      1   2
5     2.5      2   2
6     1.9      3   2
7     1.7      1   3
8     1.5      2   3
9     1.1      3   3
10    1.1      1   4
11    1.1      2   4
12    0.2      3   4

> dioxin.lm <− lm(dioxin~., data=dioxin)

> anova(dioxin.lm)

Analysis of Variance Table

Response: dioxin
          Df Sum Sq Mean Sq F value    Pr(>F)    
region     2 1.3850  0.6925   55.40 0.0001356 ***
lab        3 4.1000  1.3667  109.33 1.260e-05 ***
Residuals  6 0.0750  0.0125                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> # Tukey's single degree of freedom test

> r <− resid(dioxin.lm)

> qij <− fitted(dioxin.lm)^2

> rstar <− resid(lm(qij~dioxin$r+dioxin$l))

> 2 * sum(r*rstar) / sum(rstar^2) # γ.hat

[1] -0.2072730

> SSna <− sum(r*rstar)^2 / sum(rstar^2)

> SSna # SSna

[1] 0.02033003

> SSrem <− sum(r^2) - SSna

> SSrem # SSrem

[1] 0.05466997

> Fgam <− SSna / (SSrem/5)

> Fgam # Fγ

[1] 1.859342

> Pgam <− pf(Fgam, 1, 5, lower.tail=FALSE)

> Pgam # Pγ

[1] 0.2308840

> # The additive model seems to be adequate since P value

> # 0.2308840 is greater than 0.05. The interaction plots

> # below confirm that factor region does not interact with

> # block factor lab.

> cols <− paste("dark",c("magenta","green","blue","orange4"),sep="")

> windows(width=6,height=3,pointsize=10)

> par(mfrow=c(1,2),mar=c(4,4,.5,2)+.1,oma=c(0,0,2,0)+.1)

> with(dioxin, {

+   interaction.plot(region,lab,dioxin,type="b",

+     col=cols[1:4],pch=1:4)

+   box(which="figure",col="grey80")

+   interaction.plot(lab,region,dioxin,type="b",

+     col=cols[1:3],pch=1:3)

+   box(which="figure",col="grey80")

+ })

> title(main="Interaction Plots", outer=T)

> dev.off()

null device 
          1

(dioxin01.png here)

> # The main effects for region are significant since

> # the P value (see ANOVA) is 0.0001356. The main effect

> # plot below confirms this.

> windows(width=4,height=3,pointsize=10)

> plot(1:3,with(dioxin,tapply(dioxin,region,mean)),

+   type="l", lwd=2, xlab="Region", ylab="mean dioxin level",

+   main="Main Effect Plot\ndioxin example",

+   xaxt="n") # suppress x-axis tick marks and labels

> axis(1,at=1:3,label=as.character(1:3)) # x-axis

> abline(h=mean(dioxin$d), lty="dotted", col="red")

> dev.off()

null device 
          1

(dioxin02.png here)

> dioxin.aov <− aov(dioxin.lm)

> require(multcomp)

Loading required package: multcomp
Loading required package: mvtnorm
[1] TRUE

> dioxin.tukey <− glht(dioxin.aov,

+   linfct=mcp(region="Tukey"))

> summary(dioxin.tukey)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = dioxin.lm)

Linear Hypotheses:
           Estimate Std. Error t value p value    
2 - 1 == 0 -0.12500    0.07906  -1.581   0.323    
3 - 1 == 0 -0.77500    0.07906  -9.803  <0.001 ***
3 - 2 == 0 -0.65000    0.07906  -8.222  <0.001 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported)

> confint(dioxin.tukey)

         Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = dioxin.lm)

Estimated Quantile = 3.068

Linear Hypotheses:
           Estimate lwr     upr    
2 - 1 == 0 -0.1250  -0.3675  0.1175
3 - 1 == 0 -0.7750  -1.0175 -0.5325
3 - 2 == 0 -0.6500  -0.8925 -0.4075

95% family-wise confidence level

> windows(width=5,height=3,pointsize=10)

> par(mar=c(5,5,4,1)+.1)

> plot(dioxin.tukey,sub="dioxin data",ylim=c(.5,3.5))

> dev.off()

null device 
          1

(dioxin03.png here)

> # It appears that only the comparison between regions

> # 1 and 2 produces insignificant result.

> #

> # Residual analysis should be carried out:

> windows(width=6,height=3,pointsize=10)

> par(mfrow=c(1,2),oma=c(0,0,2,0)+.1)

> plot(dioxin.lm, which=1:2, add.smooth=FALSE)

> dev.off()

null device 
          1

(dioxin04.png here)