[R-sig-ME] Random slope/intercept without correlation in lmer
Ben Bolker
bbolker at gmail.com
Tue Aug 12 18:32:27 CEST 2014
Worth pointing out (maybe) that flexLambda branch should make this easier
too. Don't have time to generate an example right now (and fL still needs
lots of work -- we really need to get our head above water to do more
development work! Pull requests welcome ...)
Ben B.
On Tue, Aug 12, 2014 at 11:45 AM, Jake Westfall <jake987722 at hotmail.com>
wrote:
> In my opinion it is best as a general rule of thumb to always manually
> code your factors into numeric objects before passing them to a model
> fitting function, whether it is lmer() or lm() or whatever. The R Gods in
> their wisdom gave us factor objects in an attempt to make life easier for
> us, but in my experience factors often just get in the way or lead to
> unexpected results. The present issues are just another example.
>
> Jake
>
> > From: gustaf.granath at gmail.com
> > Date: Tue, 12 Aug 2014 10:40:31 -0400
> > To: steve.walker at utoronto.ca; bbolker at gmail.com
> > CC: r-sig-mixed-models at r-project.org
> > Subject: Re: [R-sig-ME] Random slope/intercept without correlation in
> lmer
> >
> > Thanks. Works perfectly.
> >
> > For more than 2 levels I guess nlme is still the way to go if you want
> > to manipulate the covariances structure.
> >
> > Gustaf
> >
> > > The `dummy` function in lme4 might be useful. Here's the example from
> > > the ?dummy help file:
> > >
> > > data(Orthodont,package="nlme")
> > > lmer(distance ~ age + (age|Subject) +
> > > (0+dummy(Sex, "Female")|Subject), data = Orthodont)
> > >
> > > Here's another example:
> > >
> > > lmer(distance ~ dummy(Sex, "Female") +
> > > (dummy(Sex, "Female") || Subject),
> > > data = Orthodont)
> > >
> > > Cheers,
> > > Steve
> > >
> > >
> > >
> > > On 2014-08-11, 6:42 PM, Ben Bolker wrote:
> > >> All you need to do is set up your own dummy variable (e.g.
> > >> ntreat=as.numeric(treat)-1, or ntreat=as.numeric(treat=="1"), or
> > >> ntreat=(original treat variable before using factor(treat) and then
> use
> > >> (ntreat||group) or (1|group)+(0+ntreat|group)
> > >>
> > >> This is related but not identical to the last example in ?lmer ;
> it's
> > >> caused by an interaction between the way that R constructs model
> > >> matrices
> > >> from factors and the way lme4 uses that functionality.
> > >>
> > >>
> > >>
> > >> On Mon, Aug 11, 2014 at 2:38 PM, Gustaf Granath
> > >> <gustaf.granath at gmail.com>
> > >> wrote:
> > >>
> > >>> Hi
> > >>> I want to model random slope and intercept without a correlation
> > >>> between
> > >>> the two. Is it possible to do this in lmer when the predictor is a
> > >>> factor?
> > >>>
> > >>> For example, imagine that x has 2 levels (control and treatment). In
> > >>> nlme,
> > >>> I have been modeling uncorrelated intercept and slope like this:
> > >>> lme(y ~ x, random=list(x = pdDiag(~ group)) ) where group is a random
> > >>> factor.
> > >>> It gives me the random intercept and random slope (i.e. variation in
> > >>> treatment effect among groups).
> > >>>
> > >>> In lmer, I think the corresponding model is defined as:
> > >>> lmer( y ~ x + (x||rand) and I guess this gives me differences
> > >>> (variation
> > >>> in differences to the intercept), but it includes a covariance term.
> > >>>
> > >>> Is it possible to reproduce the above lme() model in lmer?
> > >>>
> > >>> I have a strong feeling that Im missing something here. Most of the
> > >>> literature on this subject (and R-list questions) deals with
> continuous
> > >>> variables so pls let me know if there is a good source on this topic.
> > >>>
> > >>> Below follows a small example.
> > >>>
> > >>> Cheers
> > >>>
> > >>> Gustaf
> > >>>
> > >>>
> > >>> set.seed(1)
> > >>> treat = rep(c(0, 1), each = 5, 10)
> > >>> group = rep(1:10, each = 10)
> > >>> rand.int = rep( rnorm( 10, 0, 1), each = 10)
> > >>> rand.slop = rep( rnorm(10, 0, 1), each = 10)
> > >>> e = rnorm(100, 0, 0.5)
> > >>> y = 10 + rand.int + treat + rand.slop*treat + e
> > >>> treat = factor(treat)
> > >>>
> > >>> #lmer
> > >>> library(lme4)
> > >>> # with correlation between intercept and slope
> > >>> mod = lmer(y ~ treat + (treat|group) )
> > >>> # without correlation between intercept and slope
> > >>> # gives lots of error msgs
> > >>> mod2 = lmer(y ~ treat + (treat||group) )
> > >>> summary(mod)
> > >>> summary(mod2)
> > >>> # var-covar matrix
> > >>> VarCorr(mod)$group
> > >>> VarCorr(mod2)$group.1 #still a covariance term
> > >>>
> > >>> #nlme
> > >>> # without correlation
> > >>> library(nlme)
> > >>> lme.mod <- lme(y ~ treat, random=list(group = pdDiag(~ treat)) )
> > >>> summary(lme.mod)
> > >>> getVarCov(lme.mod)
> > >>>
> > >>> --
> > >>> Gustaf Granath (PhD)
> > >>> Post doc
> > >>> McMaster University
> > >>> School of Geography & Earth Sciences
> > >>
> > >
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list