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 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 <[email protected]> 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
>
> ______________________________________________
> [email protected] 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.
>
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
