Randomized Complete Block Design

This is a first attempt at an R version of the Randomized Complete Block example.

This example considers a situation where we want to test the effect of four different diets to be fed to just-weaned piglets. This is a study to simulate feeding of young infant humans, our Research Ethics Board would not allow us to carry out this early-stage trial using humans! We know/believe, from previous studies, that there is likely to be an effect of the sow, even post-weaning. So, we decide to use a Randomized Complete Block Design: we shall have 12 sows, and we shall take 4 male piglets from each sow, and assign one piglet to each of the four treatment diets. Thus, we shall need to consider both Sow and Diet as factors in our statistical model. We shall wish to consider Diet as a fixed effect; we are interested in these 4 specific diets. We shall wish to consider the sows not as fixed effects, but rather as a random effect; we want our conclusions not to be limited to these 12 sows, but to be generalisable across sows.

Diet Sow Gain
1 1 1 51.90
2 2 1 53.70
3 3 1 54.90
4 4 1 55.60
5 1 2 51.60
6 3 2 52.20
7 4 2 53.80
8 1 3 44.00
9 2 3 46.00
10 3 3 47.30
11 4 3 48.40
12 1 4 52.60
13 2 4 52.80
14 3 4 54.00
15 4 4 56.30
16 1 5 56.50
17 2 5 57.50
18 3 5 58.20
19 4 5 61.70
20 1 6 48.50
21 2 6 50.80
22 3 6 51.60
23 4 6 52.00
24 1 7 52.10
25 2 7 53.40
26 3 7 55.50
27 4 7 56.10
28 1 8 52.60
29 2 8 52.40
30 3 8 53.60
31 4 8 55.80
32 1 9 49.30
33 2 9 49.40
34 3 9 52.10
35 4 9 51.60
36 1 10 51.10
37 2 10 52.50
38 3 10 52.90
39 4 10 57.10
40 1 11 58.60
41 2 11 57.60
42 3 11 58.90
43 4 11 61.00
44 1 12 45.10
45 2 12 46.70
46 3 12 47.90
47 4 12 48.00

NOTE: the piglets are all kept in individual cages, so they are all independent of one another. Piglets were randomly assigned, within litters, one to each of the four diets. During the experiment one piglet on Diet 2 was lost, for reasons completely unrelated to the treatments (hence it is legitimate to consider it to be Missing At Random, MAR).

Thus, this is an RCBD, with sows as the 'blocks', and diets as our treatment.

The statistical model will be :

Yij = µ + Dieti + Sowj + eij

We shall read the data, from a plain ASCII text file, into R, we shall then use the as.factor() function to define Sow and Diet as factors, i.e. as qualitative levels, with no quantitative meaning. We shall use the lmerTest package, to fit a linear mixed model, and then esticon() and lsmeans() to estimate least squares means and differences.


#
# R statements, for a mixed effects 2-way model

library(doBy)
library(lmerTest)

# read data, fit Y = mu + diet + sow + e
 
ds<-read.table("c:/Users/Roger Cue/wwwstats/R/RCBD/blocks2.txt",header=TRUE)
summary(ds)
ds$Diet<-as.factor(ds$Diet)
ds$Sow<-as.factor(ds$Sow)
summary(ds)

# mixed model, Diet is a fixed effect, Sow is a Random effect
lm2 <- lmer(Gain~Diet+(1|Sow),data=ds)

# generate an ANOVA for the fixed effects (ONLY)
anova(lm2)

# test the significance of the random effect(s), via Log Likelihood and Chi-squared
rand(lm2)

lm2
str(lm2)

# show how we can output the random effects variance-covariance components

VarCorr(lm2)

(vc <- VarCorr(lm2))  ## default print method: standard dev and corr

## both variance and std.dev.
print(vc,comp=c("Variance","Std.Dev."),digits=4)
## variance only
print(vc,comp=c("Variance"))
as.data.frame(vc)
as.data.frame(vc,order="lower.tri")

lsmeans (lm2)

Results



> #
> # R statements, for a mixed effects 2-way model
> 
> library(doBy)
Loading required package: survival
> library(lmerTest)
Loading required package: Matrix
Loading required package: lme4

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

> # see below for use of lmerTest (probably better)
> 
> # read data, fit Y = mu + sow + trt + e
>  
> ds<-read.table("c:/Users/Roger Cue/wwwstats/R/RCBD/blocks2.txt",header=TRUE)
> summary(ds)
      Diet            Sow              Gain      
 Min.   :1.000   Min.   : 1.000   Min.   :44.00  
 1st Qu.:1.500   1st Qu.: 4.000   1st Qu.:50.95  
 Median :3.000   Median : 7.000   Median :52.60  
 Mean   :2.511   Mean   : 6.596   Mean   :52.83  
 3rd Qu.:3.500   3rd Qu.: 9.500   3rd Qu.:55.70  
 Max.   :4.000   Max.   :12.000   Max.   :61.70  
> ds$Diet<-as.factor(ds$Diet)
> ds$Sow<-as.factor(ds$Sow)
> summary(ds)
 Diet        Sow          Gain      
 1:12   1      : 4   Min.   :44.00  
 2:11   3      : 4   1st Qu.:50.95  
 3:12   4      : 4   Median :52.60  
 4:12   5      : 4   Mean   :52.83  
        6      : 4   3rd Qu.:55.70  
        7      : 4   Max.   :61.70  
        (Other):23                  
> 
> # mixed model
> lm2 <- lmer(Gain~Diet+(1|Sow),data=ds)
> anova(lm2)
Analysis of Variance Table of type III  with  Satterthwaite 
approximation for degrees of freedom
     Sum Sq Mean Sq NumDF  DenDF F.value   Pr(>F)    
Diet 88.461  29.487     3 32.013   47.91 6.11e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> rand(lm2)
Analysis of Random effects Table:
    Chi.sq Chi.DF p.value    
Sow   88.8      1  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> lm2
Linear mixed model fit by REML ['merModLmerTest']
Formula: Gain ~ Diet + (1 | Sow)
   Data: ds
REML criterion at convergence: 160.8686
Random effects:
 Groups   Name        Std.Dev.
 Sow      (Intercept) 3.8119  
 Residual             0.7845  
Number of obs: 47, groups:  Sow, 12
Fixed Effects:
(Intercept)        Diet2        Diet3        Diet4  
    51.1583       0.8666       2.1000       3.6250  
>

From the ANOVA for the fixed effects test we can see that the F-value is 47.91, with ndf=3 and ddf=32.013. The non-integer degrees of freedom for the denominator are because of the missing observation.

For the test of significance of the random effect of Sow, we see that the calculated Chi-squared value is 88.8, with 1 d.f., i.e. statistically significant, hence validating our choice of an RCBD.

Sow variance = 14.530

Residual variance = 0.6155


R.I. Cue ©
Department of Animal Science, McGill University
last updated : 2016 October 24th