[R] Possible bug in lmer nested analysis with factors

[Douglas Bates]

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.