[R] fitting linear mixed model for incomplete block design

[Douglas Bates]

Caimiao Wei wrote:
> Dear R developers and users: 
>      
> I have the following data, x is the response variable, 
 >sample(individual) nested within trt, and subsample nested within nsample,
 > I want to fit trt as fixed effect, and block, nsample(trt) as
 >random effects using lme, is the following coding correct?
> 
> dat$vgrp <- getGroups(dat, form = ~ 1|trt/nsample, level = 2)
> 
> ge.lme1 <- lme(fixed=x~trt, data=dat, random=~ block + nsample|vgrp)
> 
> 
> 
>       x block trt nsample subsample 
>       -0.68984 1 1 1 1 
>       -0.2223 2 1 1 2 
>       -1.00144 3 1 2 1 
>       -2.59511 4 1 2 2 
>       2.51573 5 1 3 1 
>       -1.67577 6 1 3 2 
>       -0.31927 7 1 4 1 
>       -3.28983 8 1 4 2 
>       0.04243 9 1 5 1 
>       -0.00635 10 1 5 2 
>       -0.09925 11 1 6 1 
>       -0.73825 12 1 6 2 
>       -0.64498 13 1 7 1 
>       -2.35105 14 1 7 2 
>       1.13354 15 1 8 1 
>       -1.45652 16 1 8 2 
>       -0.05577 17 1 9 1 
>       -3.46945 18 1 9 2 
>       -0.16415 19 1 10 1 
>       -3.31346 20 1 10 2 
>       -1.67036 1 2 1 1 
>       -2.28084 2 2 1 2 
>       -0.4165 3 2 2 1 
>       -0.51967 4 2 2 2 
>       -1.20591 5 2 3 1 
>       -1.33746 6 2 3 2 
>       2.40848 7 2 4 1 
>       3.7884 8 2 4 2 
>       1.14388 9 2 5 1 
>       0.45225 10 2 5 2 
>       -0.5008 11 2 6 1 
>       -1.66627 12 2 6 2 
>       0.97082 13 2 7 1 
>       -0.14259 14 2 7 2 
>       -2.14964 15 2 8 1 
>       -1.77236 16 2 8 2 
>       -0.0305 17 2 9 1 
>       -1.23726 18 2 9 2 
>       0.43089 19 2 10 1 
>       1.07422 20 2 10 2 
> 
> 
> Thanks,
> 
> Caimiao

It looks as if there is some crossing between block and nsample(trt) so 
I think you should probably use the version of lme from the lme4 package.
 > dat$vgrp <- with(dat, nsample:trt)
 > xtabs(~block + vgrp, dat)
      vgrp
block 1:1 1:2 2:1 2:2 3:1 3:2 4:1 4:2 5:1 5:2 6:1 6:2 7:1 7:2 8:1 8:2 
9:1 9:2
    1    1   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0
    2    1   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0
    3    0   0   1   1   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0
    4    0   0   1   1   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0
    5    0   0   0   0   1   1   0   0   0   0   0   0   0   0   0   0 
  0   0
    6    0   0   0   0   1   1   0   0   0   0   0   0   0   0   0   0 
  0   0
    7    0   0   0   0   0   0   1   1   0   0   0   0   0   0   0   0 
  0   0
    8    0   0   0   0   0   0   1   1   0   0   0   0   0   0   0   0 
  0   0
    9    0   0   0   0   0   0   0   0   1   1   0   0   0   0   0   0 
  0   0
    10   0   0   0   0   0   0   0   0   1   1   0   0   0   0   0   0 
  0   0
    11   0   0   0   0   0   0   0   0   0   0   1   1   0   0   0   0 
  0   0
    12   0   0   0   0   0   0   0   0   0   0   1   1   0   0   0   0 
  0   0
    13   0   0   0   0   0   0   0   0   0   0   0   0   1   1   0   0 
  0   0
    14   0   0   0   0   0   0   0   0   0   0   0   0   1   1   0   0 
  0   0
    15   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   1 
  0   0
    16   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1   1 
  0   0
    17   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  1   1
    18   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  1   1
    19   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0
    20   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 
  0   0
      vgrp
block 10:1 10:2
    1     0    0
    2     0    0
    3     0    0
    4     0    0
    5     0    0
    6     0    0
    7     0    0
    8     0    0
    9     0    0
    10    0    0
    11    0    0
    12    0    0
    13    0    0
    14    0    0
    15    0    0
    16    0    0
    17    0    0
    18    0    0
    19    1    1
    20    1    1
 > fm1 <- lme(x ~ trt, dat, list(block = ~ 1, vgrp = ~ 1))
 > fm1
Linear mixed-effects model
Fixed: x ~ trt
  Data: dat
  log-restricted-likelihood:  -72.6419
Random effects:
  Groups   Name        Variance   Std.Dev.
  vgrp     (Intercept) 7.3585e-01 0.85781684
  block    (Intercept) 3.2395e-08 0.00017999
  Residual             1.7037e+00 1.30524021
# of obs: 40, groups: vgrp, 20; block, 20

Fixed effects:
             Estimate Std. Error DF t value Pr(>|t|)
(Intercept) -0.92005    0.39846 38 -2.3090  0.02647
trt2         0.68699    0.56350 38  1.2191  0.23030