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 :

Y_{ij} = µ + Diet_{i} + Sow_{j} +
e_{ij}

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 24