[R-sig-ME] [R] Multilevel model in lme4 and nlme

[Ben Bolker]

  [forwarding to r-sig-mixed-models list instead of r-help]

On 09/13/2011 04:54 AM, jonas garcia wrote:
> Hi Ben, thanks for your reply.
> 
> Your suggestion does not work indeed:
> 
> 
> 
> lme(y ~ x, random=list(~1|a:b, ~1|b:c), data=mydata)
> 
> Error in getGroups.data.frame(dataMix, groups) :
> 
>   Invalid formula for groups
> 
> 

  Here's a work-around, I think: define your interactions/nesting
explicitly.


mydata2 <- transform(mydata,
                     intab=interaction(a,b),
                     intabc=interaction(a,b,c),
                     intbc=interaction(b,c))

 lme(y~x,data=mydata2,random=list(intab=~1,intbc=~1))

  I'm still not sure that your nesting scheme makes sense, and note that
lme requires/assumes nesting (which lmer does not): look carefully
at the number of groups reported for each level in the model output and
convince yourself that they make sense.  (I don't have the time and
energy to puzzle through this at the moment.)

  good luck,

    Ben Bolker

> 
> Here is a reproducible example of my data:
> 
> 
> 
> set.seed(123)
> 
> library(lme4)
> 
> library(nlme)
> 
> 
> 
> y<- rnorm(30)
> 
> x<-rnorm(30)
> 
> a<- factor(sort(rep(c("alpha", "beta", "charlie"), 10)))
> 
> b<- factor(rep(c("rho", "epsilon", "lambda"), 10))
> 
> c<- factor(c(sort(rep(1:2, 5)), sort(rep(3:4, 5)), sort(rep(5:6, 5))))
> 
> mydata<- data.frame(y,x,a,b,c)
> 
> 
> 
> 
> 
> mod1<- lmer(y ~ x + (1|a:b) + (1|b:c), data=mydata)
> 
> 
> 
> mod2.lme<- lme(y ~ x, random=list(a=~1, b=~1, c=~1), data=mydata)
> 
> mod2.lmer<- lmer(y ~ x + (1|a) + (1|a:b) + (1|a:b:c), data=mydata)
> 
> 
> 
> My objective is to specify mod1 using function lme.
> 
> Anyone knows how to do it?
> 
> Thanks
> 
> 
> 
> J
> 
> 
> 
> 
> On Mon, Sep 12, 2011 at 9:43 PM, Ben Bolker <bbolker at gmail.com> wrote:
> 
>> jonas garcia <garcia.jonas80 <at> googlemail.com> writes:
>>
>>> I am trying to fit some mixed models using packages lme4 and nlme.
>>>
>>> I did the model selection using lmer but I suspect that I may have some
>>> autocorrelation going on in my data so I would like to have a look using
>> the
>>> handy correlation structures available in nlme.
>>>
>>> The problem is that I cannot translate my lmer model to lme:
>>>
>>> mod1<- lmer(y~x + (1|a:b) + (1|b:c), data=mydata)
>>>
>>> "a", "b" and "c" are factors with "c" nested in "b" and "b" nested in "a"
>>>
>>> The best I can do with lme is:
>>>
>>> mod2<- lme(y~x, random=list(a=~1, b=~1, c=~1), data=mydata)
>>>
>>> which is the same as:
>>>
>>> lmer(y~x + (1|a) + (1|a:b) + (1|a:b:c), data=mydata)
>>>
>>> I am not at all interested in random effects (1|a) and (1|a:b:c) as they
>> are
>>> not significant. I just need two random intercepts as specified in mod1.
>> How
>>> can I translate mod1 into lme language?
>>>
>>> Any help on this would be much appreciated.
>>
>>  This would probably be better on the r-sig-mixed-models list.
>>
>>  Does random=list(~1|a:b,~1|b:c) work?
>>
>>  I would be a little bit careful throwing out ~1|a (non-significance
>> is not necessarily sufficient reason to discard a term from the model --
>> it depends a lot on your procedure), and with the interpretation of
>> your nesting.  If b is only explicitly and not implicitly nested in a
>> (i.e. if there a levels of 'b' that occur in more than one level of 'a',
>> for example if a corresponded to families, b corresponded to individuals,
>> and you labeled individuals 1..N_b_i in each family) then I'm not
>> sure how you would actually interpret b:c, as it would be crossed
>> rather than nested.  But assuming that your model specification in
>> lmer is correct and sensible, I think my suggestion above should (?)
>> work to get the equivalent in lme.
>>>
>>> Jonas
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>