# Re: [R] extracting nested variances from lme4 model

From: Douglas Bates <bates_at_stat.wisc.edu>
Date: Sun 15 Oct 2006 - 22:59:26 GMT

On 10/15/06, Douglas Bates <bates@stat.wisc.edu> wrote:
> On 10/14/06, Spencer Graves <spencer.graves@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.

> (fm2 <- lmer(x ~(1|rtr) + (trth|cs), dtf))
Linear mixed-effects model fit by REML
Formula: x ~ (1 | rtr) + (trth | cs)

Data: dtf
AIC BIC logLik MLdeviance REMLdeviance  252.5 260.9 -121.2 244.8 242.5 Random effects:

``` Groups   Name        Variance Std.Dev. Corr
cs       (Intercept)  1.9127  1.3830
trthTRUE    24.1627  4.9156   1.000
rtr      (Intercept)  1.9128  1.3830
Residual             18.5278  4.3044
```

number of obs: 40, groups: cs, 10; rtr, 4
```            Estimate Std. Error t value
(Intercept)  -0.2128     1.2723 -0.1673
```

Warning message:
nlminb returned message false convergence (8)
in: `LMEoptimize<-`(`*tmp*`, value = list(maxIter = 200, tolerance = 1.49011611938477e-08,
```>
>

> > 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@stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help