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)
> library(lme4)
> library(car)
> wear.lm <- lm(wear~block+trt,data=weardata)
> plot(wear.lm,which=1)
Residual plot of wear data

Figure 42.1: Residual plot of wear data

Box-Cox suggests a log.

> boxCox(wear.lm)
Box-Cox plot of wear data

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)
Residual plot of log wear data

Figure 42.3: Residual plot of log wear data

Residuals are a bit short tailed.

> plot(logwear.lm,which=2)
NPP of residuals of log wear data

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"