On 10/15/06, Douglas Bates <[EMAIL PROTECTED]> wrote: > On 10/14/06, Spencer Graves <[EMAIL PROTECTED]> 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.