Re: [R] extracting nested variances from lme4 model

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon 16 Oct 2006 - 23:57:42 GMT

Dear Doug & Frank:

      Thanks for the reply, Doug. First a comment, then another question for you, Doug, if you have time:

  1. COMMENT: Here's an ugly hack that would seem to answer Frank's question using 'lmer':

(vt <- with(dtf, tapply(x, trth, sd)))
(vr <- vt[2]/vt[1])

mod1b <- lmer(x~trth+(1|rtr)+(I(vt[1]*(1+vr*trth))-1|cs),

              data=dtf)

      Now put this is a loop over a range of values for 'vr' and pick the one that maximizes 'logLik' (or minimizes AIC or BIC); it should be fairly close if not identical to this initial estimate.

      2. ANOTHER QUESTION: Can the desired model be fit using 'lme', with crossed random effects, and we want separate variance estimates for the random effect 'cs' for each level of the fixed effect 'trth'? I previously suggested chapters 4 and 5 in Pinheiro and Bates (2000). However, even with the scripts "ch04.R" and "ch05.R" in "~library\nlme\scripts", it might take me a while to figure out how to do it (and even answer the question of whether it can).

      Thanks again. 
      Spencer Graves

Douglas Bates wrote:
> 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 Tue Oct 17 10:13:51 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 - 01: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.