[R] extracting nested variances from lme4 model

[Douglas Bates]

On 10/14/06, Spencer Graves <spencer.graves at pdf.com> wrote:
>       You want to estimate a different 'cs' variance for each level of
> 'trth', as indicated by the following summary from your 'fake data set':
>
>  > tapply(dtf$x, dtf$trth, sd)
>    FALSE     TRUE
> 1.532251 8.378206
>
>       Since var(x[trth]) > var(x[!trth]), I  thought that the following
> should produce this:
>
>       mod1.<-lmer( x ~  (1|rtr)+ (trth|cs) , data=dtf)
>
>       Unfortunately, this generates an error for me:
>
> Error in lmer(x ~ (1 | rtr) + (trth | cs), data = dtf) :
>     the leading minor of order 1 is not positive definite
>
>       I do not understand this error, and I don't have other ideas about
> how to get around it using 'lmer'.

Hmm.  It's an interesting problem.  I can tell you where the error
message originates but I can't yet tell you why.

The error message originates when Cholesky decompositions of the
variance-covariance matrices for the random effects are generated at
the initial estimates.  It looks as if one of the model matrices is
not being formed correctly for some reason.  I will need to do more
debugging.


> It seems to me that 'lme' in
> library(nlme) should also be able to handle this problem, but 'lme' is
> designed for nested factors, and your 'rtr' and 'cs' are crossed.  If it
> were my problem, I think I'd study ch. 4 and possibly ch. 5 of Pinheiro
> and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer) on
> 'pdMat' and 'corrStruct' classes, because I believe those may support
> models with crossed random effects like in your example AND might also
> support estimating different variance components for different levels of
> a fixed effect like 'trth' in your example.
>
>       If we are lucky, someone else might educate us both.
>
>       I'm sorry I couldn't answer your question, but I hope these
> comments might help in some way.
>
>       Spencer Graves
>
> Frank Samuelson wrote:
> > I have a model:
> > mod1<-lmer( x ~  (1|rtr)+ trth/(1|cs) , data=dtf)  #
> >
> > Here, cs and rtr are crossed random effects.
> > cs 1-5 are of type TRUE, cs 6-10 are of type FALSE,
> > so cs is nested in trth, which is fixed.
> > So for cs I should get a fit for 1-5 and 6-10.
> >
> > This appears to be the case from the random effects:
> >  > mean( ranef(mod1)$cs[[1]][1:5] )
> > [1] -2.498002e-16
> >  > var( ranef(mod1)$cs[[1]][1:5] )
> > [1] 23.53083
> >  > mean( ranef(mod1)$cs[[1]][6:10] )
> > [1] 2.706169e-16
> >  > var( ranef(mod1)$cs[[1]][6:10] )
> > [1] 1.001065
> >
> > However VarCorr(mod1) gives me only
> > a single variance estimate for the cases.
> > How do I get the other variance estimate?
> >
> > Thanks for any help.
> >
> > -Frank
> >
> >
> >
> > My fake data set:
> > -------------------------------------------
> > data:
> > $frame
> >               x rtr  trth cs
> > 1   -4.4964808   1  TRUE  1
> > 2   -0.6297254   1  TRUE  2
> > 3    1.8062857   1  TRUE  3
> > 4    2.7273275   1  TRUE  4
> > 5   -0.6297254   1  TRUE  5
> > 6   -1.3331767   1 FALSE  6
> > 7   -0.3055796   1 FALSE  7
> > 8    1.3331767   1 FALSE  8
> > 9    0.1574314   1 FALSE  9
> > 10  -0.1574314   1 FALSE 10
> > 11  -3.0958803   2  TRUE  1
> > 12  12.4601615   2  TRUE  2
> > 13  -5.2363532   2  TRUE  3
> > 14   1.5419810   2  TRUE  4
> > 15   7.7071005   2  TRUE  5
> > 16  -0.3854953   2 FALSE  6
> > 17   0.7673098   2 FALSE  7
> > 18   2.9485984   2 FALSE  8
> > 19   0.3854953   2 FALSE  9
> > 20  -0.3716558   2 FALSE 10
> > 21  -6.4139940   3  TRUE  1
> > 22  -3.8162344   3  TRUE  2
> > 23 -11.0518554   3  TRUE  3
> > 24   2.7462631   3  TRUE  4
> > 25  16.2543990   3  TRUE  5
> > 26  -1.7353668   3 FALSE  6
> > 27   1.3521369   3 FALSE  7
> > 28   1.3521369   3 FALSE  8
> > 29   0.2054067   3 FALSE  9
> > 30  -1.7446691   3 FALSE 10
> > 31  -6.7468356   4  TRUE  1
> > 32 -11.3228243   4  TRUE  2
> > 33 -18.0088554   4  TRUE  3
> > 34   4.6264891   4  TRUE  4
> > 35   9.2973991   4  TRUE  5
> > 36  -4.1122576   4 FALSE  6
> > 37  -0.5091994   4 FALSE  7
> > 38   1.2787085   4 FALSE  8
> > 39  -1.6867089   4 FALSE  9
> > 40  -0.5091994   4 FALSE 10
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>