Re: [R] extracting nested variances from lme4 model

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Sun 15 Oct 2006 - 02:00:06 GMT

      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'. 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 Sun Oct 15 12:07:15 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 Sun 15 Oct 2006 - 17: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.