Chapter 42 Incomplete Blocks
42.1 Wear data
I’m embarrassed to say, but I can’t really remember where I got these data. I think from one of my professors back in the dark ages. The issue is wear loss of rubber, and we have seven different treatments for reducing the loss. A block is a chunk of rubber, but a block is only large enough for four of the seven treatments. Small losses are good.
> weardata <- read.table("wear.dat.txt",header=TRUE)
> weardata <- within(weardata,{block <- factor(block);trt <- factor(trt)})
> weardata
block trt wear
1 2 1 344
2 4 1 337
3 6 1 369
4 7 1 396
5 1 2 627
6 4 2 537
7 5 2 520
8 7 2 602
9 2 3 233
10 3 3 251
11 5 3 278
12 7 3 240
13 1 4 248
14 3 4 211
15 6 4 196
16 7 4 273
17 3 5 160
18 4 5 195
19 5 5 199
20 6 5 185
21 1 6 563
22 2 6 442
23 5 6 595
24 6 6 606
25 1 7 252
26 2 7 226
27 3 7 297
28 4 7 300
The incidence matrix \(n_{ij}\). Each block contains four treatments, and each pair of treatments occurs together twice.
> with(weardata,table(block,trt))
trt
block 1 2 3 4 5 6 7
1 0 1 0 1 0 1 1
2 1 0 1 0 0 1 1
3 0 0 1 1 1 0 1
4 1 1 0 0 1 0 1
5 0 1 1 0 1 1 0
6 1 0 0 1 1 1 0
7 1 1 1 1 0 0 0
Relative efficiency is 7/8
> 7*3/(6*4)
[1] 0.875
Effective sample size is 3.5, not 4.
> .875*4
[1] 3.5
Basic intrablock analysis is treatments adjusted for blocks. It looks like there could be increasing variance.
> library(cfcdae)
Registered S3 method overwritten by 'DoE.base':
method from
factorize.factor conf.design
> library(lme4)
Loading required package: Matrix
> library(car)
Loading required package: carData
> wear.lm <- lm(wear~block+trt,data=weardata)
> plot(wear.lm,which=1)

Figure 42.1: Residual plot of wear data
Box-Cox suggests a log.
> boxCox(wear.lm)

Figure 42.2: Box-Cox plot of wear data
Refit with logs; residuals now look better.
> logwear.lm <- lm(log(wear)~block+trt,data=weardata)
> plot(logwear.lm,which=1)

Figure 42.3: Residual plot of log wear data
Residuals are a bit short tailed.
> plot(logwear.lm,which=2)

Figure 42.4: NPP of residuals of log wear data
All the usual things work. Treatment 5 looks best, but can’t really be distinguished from 4. Stay away from 2 and 6.
> anova(logwear.lm)
Analysis of Variance Table
Response: log(wear)
Df Sum Sq Mean Sq F value Pr(>F)
block 6 0.7784 0.12974 11.959 5.573e-05 ***
trt 6 3.9060 0.65100 60.006 1.230e-09 ***
Residuals 15 0.1627 0.01085
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(logwear.lm)
Call:
lm(formula = log(wear) ~ block + trt, data = weardata)
Residuals:
Min 1Q Median 3Q Max
-0.13172 -0.06245 0.01866 0.05593 0.11117
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.757176 0.019684 292.480 < 2e-16 ***
block1 0.040494 0.051545 0.786 0.444326
block2 -0.142221 0.051545 -2.759 0.014614 *
block3 -0.033005 0.051545 -0.640 0.531628
block4 0.013441 0.051545 0.261 0.797822
block5 0.050401 0.051545 0.978 0.343675
block6 -0.008792 0.051545 -0.171 0.866844
trt1 0.145530 0.051545 2.823 0.012839 *
trt2 0.542077 0.051545 10.517 2.57e-08 ***
trt3 -0.224702 0.051545 -4.359 0.000561 ***
trt4 -0.338553 0.051545 -6.568 8.91e-06 ***
trt5 -0.547229 0.051545 -10.617 2.26e-08 ***
trt6 0.562861 0.051545 10.920 1.55e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1042 on 15 degrees of freedom
Multiple R-squared: 0.9664, Adjusted R-squared: 0.9396
F-statistic: 35.98 on 12 and 15 DF, p-value: 7.699e-09
> sidelines(pairwise(logwear.lm,trt))
5 -0.547 |
4 -0.339 | |
3 -0.225 |
7 -0.140 |
1 0.146
2 0.542 |
6 0.563 |
Here we have three attempts to compute the standard error of a pairwise difference. The first is a naive application of the standard formula; this gives the incorrect answer. The second
version uses the effective sample size (3.5) and gets the correct answer. The last simply verifies the effective sample
size computation via linear.contrast()
.
> sqrt(.010849*(1/4 + (-1)^2/4))
[1] 0.07365121
> sqrt(.010849*(1/3.5 + (-1)^2/3.5))
[1] 0.07873645
> linear.contrast(logwear.lm,trt,c(1,-1,0,0,0,0,0))
estimates se t-value p-value lower-ci upper-ci
1 -0.3965467 0.07873613 -5.036401 0.0001475959 -0.5643688 -0.2287246
attr(,"class")
[1] "linear.contrast"
42.2 Interblock recovery
If you think blocks are random, you can do interblock recovery.
All you need to do is fit using lmer()
with blocks as random.
Block effects are not very large relative to either errors or treatment effects. The standard error of treatment estimates are just slightly smaller than what we achieved with intrablock analysis.
> logwear.lmer <- lmer(log(wear)~trt+(1|block),data=weardata)
> summary(logwear.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: log(wear) ~ trt + (1 | block)
Data: weardata
REML criterion at convergence: -18.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.5380 -0.7593 0.1167 0.5840 1.4015
Random effects:
Groups Name Variance Std.Dev.
block (Intercept) 0.002251 0.04744
Residual 0.010849 0.10416
Number of obs: 28, groups: block, 7
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.75718 0.02663 216.213
trt1 0.13715 0.04964 2.763
trt2 0.56873 0.04964 11.456
trt3 -0.23124 0.04964 -4.658
trt4 -0.32720 0.04964 -6.591
trt5 -0.54404 0.04964 -10.959
trt6 0.55415 0.04964 11.163
Correlation of Fixed Effects:
(Intr) trt1 trt2 trt3 trt4 trt5
trt1 0.000
trt2 0.000 -0.167
trt3 0.000 -0.167 -0.167
trt4 0.000 -0.167 -0.167 -0.167
trt5 0.000 -0.167 -0.167 -0.167 -0.167
trt6 0.000 -0.167 -0.167 -0.167 -0.167 -0.167
The KR anova yields a slightly higher F and slightly higher error df, the result of more information. And SE for pairwise contrast is slightly smaller, etc.
> anova(logwear.lm)
Analysis of Variance Table
Response: log(wear)
Df Sum Sq Mean Sq F value Pr(>F)
block 6 0.7784 0.12974 11.959 5.573e-05 ***
trt 6 3.9060 0.65100 60.006 1.230e-09 ***
Residuals 15 0.1627 0.01085
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> Anova(logwear.lmer,test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(wear)
F Df Df.res Pr(>F)
trt 61.09 6 17.337 1.036e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> linear.contrast(logwear.lm,trt,c(1,-1,0,0,0,0,0))
estimates se t-value p-value lower-ci upper-ci
1 -0.3965467 0.07873613 -5.036401 0.0001475959 -0.5643688 -0.2287246
attr(,"class")
[1] "linear.contrast"
> linear.contrast(logwear.lmer,trt,c(1,-1,0,0,0,0,0))
estimates se t-value p-value lower-ci upper-ci
1 -0.4315824 0.07854585 -5.494656 3.681049e-05 -0.5970545 -0.2661103
attr(,"class")
[1] "linear.contrast"