I have changed the subject line because the subject has changed.
Mea culpa, here is the example that I intended.
> kiwishade$vine <- factor(rep(1:4,12))
> ## The following offers, at least in simpler cases, a check.
> unique(with(kiwishade, xtabs(~block+shade+vine)))
[1] 1
> library(nlme)
> VarCorr(lme(yield ~ shade, random=~1|block/shade/vine,
data=kiwishade))
Variance StdDev
block = pdLogChol(1)
(Intercept) 4.077830 2.019364
shade = pdLogChol(1)
(Intercept) 2.186341 1.478628
vine = pdLogChol(1)
(Intercept) 10.617337 3.258426
Residual 1.565418 1.251167
> VarCorr(lme(yield ~ shade, random=~1|block/shade, data=kiwishade))
Variance StdDev
block = pdLogChol(1)
(Intercept) 4.077868 2.019373
shade = pdLogChol(1)
(Intercept) 2.186325 1.478623
Residual 12.182756 3.490381
Notice that, with 'random=~1|block/shade/vine', the Residual
component of variances gets split in two parts. I have no idea why
this happens. I notice that lmer() throws an error if given the
formula ' ~ shade + (1|block/shade/vine)'
My earlier example doubled up on the random term for the block:shade
level. The lme() code, obviously wanting to please, obligingly split
the component of variance into two more or less equal parts. In this
instance, lmer() gives very nearly the same result.
lme4 Matrix
"0.999375-0" "0.999375-0"
> library(lme4)
Loading required package: Matrix
> lmer(yield ~ shade + (1|block/shade/plot), data=kiwishade)
. . . .
Random effects:
Groups Name Variance Std.Dev.
shade:block (Intercept) 1.0965 1.0471
plot:(shade:block) (Intercept) 1.0899 1.0440
block (Intercept) 4.0769 2.0191
Residual 12.1829 3.4904
Number of obs: 48, groups: shade:block, 12; plot:(shade:block), 12;
block, 3
You might like to compare the following, using the kiwishade data
frame from the DAAG package
VarCorr(lme(yield ~ shade, random=~1|block/shade/plot,
data=kiwishade)) ## The units are plots
VarCorr(lme(yield ~ shade, random=~1|block/shade, data=kiwishade))
The first one splits the block:shade component of variance between
block:shade and block:shade:plot This is just wrong. The alleged
block:shade:plot component has nothing to do with plots. It is
some deep mystery in the internal workings of lme() that it splits
the block:shade component of variance in this strange way rather
than throwing an error
The following is incorrect
or recognizing that block:shade:plot is just another name for the
residual, and that it should act accordingly.
The following is correct.
Doug may be able to throw light on this.
John Maindonald email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473 fax : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.