Caimiao Wei
2004-Oct-29 18:14 UTC
[R] fitting linear mixed model for incomplete block design
Dear R developers and users:
I have the following data, x is the response vaiable, nsample(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
[[alternative HTML version deleted]]
Douglas Bates
2004-Oct-29 19:47 UTC
[R] fitting linear mixed model for incomplete block design
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, > > CaimiaoIt 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