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.

OK, I have worked out why the error occurs and will upload lme4_0.9975-4, which fixes this problem and also has some speedups when fitting generalized linear mixed models. Just for the record, I had assumed that the cholmod_aat function in the CHOLMOD sparse matrix library (the function cholmod_aat calculates AA' from a sparse matrix A) returned a matrix in which the columns are sorted in increasing order of the rows. That is not always correct but the situation is easily remedied. (The typical way to sort the columns in sparse matrices is to take the transpose of the transpose, which had me confused when I first saw it in some of Tim Davis's code. It turns out that a side effect of this operation is to sort the columns in increasing order of the rows.)

With this patched lme4 package the model Spencer proposed can be fit but the variance-covariance matrix of the two-dimensional random effects for the cs factor is singular because of the design (the level of trth is constant within each level of cs). I'm not sure how to specify the desired model in lmer.

> (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

Fixed effects:

            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
> > > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
> >
>

R-help@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. Received on Mon Oct 16 09:04:52 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 17 Oct 2006 - 00:30:10 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.