[R] extracting nested variances from lme4 model
Douglas Bates
bates at stat.wisc.edu
Mon Oct 16 00:59:26 CEST 2006
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