Hello, Is this a bug in the lmer routine? > library(lme4) > ### test case based on rats data from Crawley > a<-rnorm(36);b<-rep(1:3,each=12);c<-rep(1:2,each=6,3);d<-rep (1:3,each=2,6) > > ### mixed model works when c & d are numeric, lmer assumes they are factors > m <- lmer(a ~ b + (1|c/d)) > > ### but bails out when they are actually specified as factors > c<-factor(c); d<-factor(d) > m <- lmer(a ~ b + (1| c / d)) Error in lmer(a ~ b + (1 | c/d)) : entry 0 in matrix[0,0] has row 2147483647 and column 2147483647 In addition: Warning message: / not meaningful for factors in: Ops.factor(c, d) Cheers Yan
On 16 Sep 2005, at 16:57, Yan Wong wrote:> Hello, > > Is this a bug in the lmer routine? > > ... > Error in lmer(a ~ b + (1 | c/d)) : entry 0 in matrix[0,0] has row > 2147483647 and column 2147483647 > In addition: Warning message: > / not meaningful for factors in: Ops.factor(c, d)Sorry, I forgot to specify: R 2.1.1, Mac OS X 10.4.2, lme4 version 0.98-1, Matrix version 0.98-7. Yan
I think you might have confused lme code with lmer code. Why do you have c/d in the random portion? I think what you want is> lmer(a ~ b + (1 | c)+(1|d))Which gives the following using your data Linear mixed-effects model fit by REML Formula: a ~ b + (1 | c) + (1 | d) AIC BIC logLik MLdeviance REMLdeviance 108.0239 115.9415 -49.01193 94.57296 98.02386 Random effects: Groups Name Variance Std.Dev. d (Intercept) 4.2877e-10 2.0707e-05 c (Intercept) 4.2877e-10 2.0707e-05 Residual 8.5754e-01 9.2603e-01 # of obs: 36, groups: d, 3; c, 2 Fixed effects: Estimate Std. Error DF t value Pr(>|t|) (Intercept) 0.82928 0.40834 34 2.0309 0.05015 . b -0.37563 0.18903 34 -1.9872 0.05500 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Yan Wong Sent: Friday, September 16, 2005 11:58 AM To: R-help Subject: [R] Possible bug in lmer nested analysis with factors Hello, Is this a bug in the lmer routine? > library(lme4) > ### test case based on rats data from Crawley > a<-rnorm(36);b<-rep(1:3,each=12);c<-rep(1:2,each=6,3);d<-rep (1:3,each=2,6) > > ### mixed model works when c & d are numeric, lmer assumes they are factors > m <- lmer(a ~ b + (1|c/d)) > > ### but bails out when they are actually specified as factors > c<-factor(c); d<-factor(d) > m <- lmer(a ~ b + (1| c / d)) Error in lmer(a ~ b + (1 | c/d)) : entry 0 in matrix[0,0] has row 2147483647 and column 2147483647 In addition: Warning message: / not meaningful for factors in: Ops.factor(c, d) Cheers Yan ______________________________________________ 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
Doug Bates has the following article in R News. To date, it is the only
source I know of documenting the lmer function.
@Article{Rnews:Bates:2005,
author = {Douglas Bates},
title = {Fitting Linear Mixed Models in {R}},
journal = {R News},
year = 2005,
volume = 5,
number = 1,
pages = {27--30},
month = {May},
url = {http://CRAN.R-project.org/doc/Rnews/},
}
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Yan Wong
Sent: Friday, September 16, 2005 1:58 PM
To: R-help
Subject: Re: [R] Possible bug in lmer nested analysis with factors
On 16 Sep 2005, at 17:12, Doran, Harold wrote:
> I think you might have confused lme code with lmer code. Why do you
> have c/d in the random portion?
Apologies. I obviously have done something of the sort. I assumed that
the 'random' assignment in lme could just be incorporated into an lmer
call by placing it in brackets and removing the ~, so that
lme(a ~ b, random= ~ 1|c/d)
would be equivalent to
lmer(a ~ b + (1|c/d))
Is there a good guide somewhere to lmer calling conventions? I obviously
don't understand them. As you can see, I would like to nest d within c,
(and actually, c is nested in b too).
Perhaps it would be better with some explanation of the Crawley data.
There are 3 fixed drug treatments ('b') given to 2 rats (6 rats in
all: 'c'), and 3 samples ('d') are taken from each of the
rat's livers,
with some response variable recorded for each sample ('a':
here just allocated a Normal distribution for testing purposes). I.e.
c and d are random effects, with d %in% c and c %in% b.
He analyses it via
aov(a ~ b+c+d+Error(a/b/c))
I'm wondering what the lme and lmer equivalents are. I've been told that
a reasonable form of analysis using lme is
a<-rnorm(36);b<-rep(1:3,each=12);d<-rep(1:3,each=2,6)
c <- rep(1:6, each=6) #use unique labels for each rat ## I got this
wrong in my previous example
model1 <- lme(a ~ b, random= ~ 1|c/d)
Which gives what looks to be a reasonable output (but I'm new to all
this mixed modelling stuff). How would I code the same thing using lmer?
> I think what you want is
>
>> lmer(a ~ b + (1 | c)+(1|d))
>>
>
> Which gives the following using your data
I'm not sure this is what I wanted to do. But thanks for the all the
help.
Yan
______________________________________________
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
I have a couple of other comments. You can write the nested grouping factors as Sundar did without having to explicitly evaluate the interaction term and drop unused levels. The lmer function, like most modeling functions in R, uses the optional argument drop.unused.levels = TRUE when creating the model frame. John Maindonald has already suggested the use of (1|b/c) => (1|b:c) + (1|b) as "syntactic sugar" for the lmer formula and it is a reasonable request. Unfortunately, implementing this is not high on my priority list at present. (We just made a massive change in the sparse matrix implementation in the Matrix package and shaking the bugs out of that will take a while.) In any case the general recommendation for nested grouping factors is first to ensure that they are stored as factors and then to create the sequence of interaction terms. On 9/18/05, Yan Wong <h.y.wong at leeds.ac.uk> wrote:> > On 16 Sep 2005, at 17:21, Sundar Dorai-Raj wrote: > > > My guess is he wants this: > > > > c1 <- factor(c) > > d1 <- factor(d) > > m <- lmer(a ~ b + (1|c1:d1)+(1|c1)) > > > > which assumes d1 is nested within c1. > > > > Take a look at Section 3 in the "MlmSoftRev" vignette: > > > > library(mlmRev) > > vignette("MlmSoftRev") > > Ah, that vignette is extremely useful: it deserves to be more widely > known. > I mainly intended this reply to be a thank you to yourself and Harold. > > Unfortunately, there is (as always), one last thing that is still > puzzling me: > the df for fixed factors seems to vary between what I currently > understand to > be equivalent calls to lme and lmer, e.g: > > ------- > > a<-rnorm(36); > b<-factor(rep(1:3,each=12)) > c<-factor(rep(1:2,each=6,3)) > d<-factor(rep(1:3,each=2,6)) > c <- evalq(b:c)[,drop=T] #make unique factor levels > d <- evalq(c:d)[,drop=T] #make unique factor levels > > summary(lme(a ~ b, random=~1|c/d)) > # output includes estimates for fixed effects such as > # Value Std.Error DF t-value p-value > # (Intercept) 0.06908901 0.3318330 18 0.2082041 0.8374 > # b2 0.13279084 0.4692828 3 0.2829655 0.7956 > # b3 -0.26146698 0.4692828 3 -0.5571630 0.6163 > > # I understand the above lme model to be equivalent to > summary(lmer(a ~ b + (1|c) +(1|c:d)) > #but this gives fixed effects estimates with differing DF: > # Estimate Std. Error DF t value Pr(>|t|) > # (Intercept) 0.069089 0.331724 33 0.2083 0.8363 > # b2 0.132791 0.469128 33 0.2831 0.7789 > # b3 -0.261467 0.469128 33 -0.5573 0.5811 > > Again, many thanks for your help: even more so if you or anyone > else can answer this last little niggle of mine.