[R] extracting nested variances from lme4 model
Spencer Graves
spencer.graves at pdf.com
Tue Oct 17 01:57:42 CEST 2006
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 at stat.wisc.edu> wrote:
>> On 10/14/06, Spencer Graves <spencer.graves at 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 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.
>> > >
>> >
>>
More information about the R-help
mailing list