[R] extracting nested variances from lme4 model

[Spencer Graves]

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