[R-sig-ME] A spurious splitting of a variance component

[John Maindonald]

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.