[R-sig-ME] A spurious splitting of a variance component
John Maindonald
john.maindonald at anu.edu.au
Wed Sep 12 04:02:47 CEST 2007
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.
More information about the R-sig-mixed-models
mailing list