[R-sig-ME] lmer vs. nlme
Ashley Cohen
ashleyhcohen at gmail.com
Tue Jan 19 21:37:48 CET 2016
Thank you very much for your response, Ben.
I had more luck using lme to play around with the correlation structure, as
you mentioned, and was able to compare compound symmetry, AR(1) and
unstructured.
I think because I am trying to specify it (I don't know that I want to use
unstructured as we are lacking data) it is best not to 'muck around' with
the correlation matrix in lme4, as you mentioned. However, there doesn't
even seem to be much within-subject correlation anyways.
Ideally, I think an AR(1) correlation matrix works best for this data.
On Mon, Jan 18, 2016 at 5:02 PM, Ben Bolker <bbolker at gmail.com> wrote:
> I think what's going on here is that with this particular
> specification *both* models are effectively fitting _positive_
> compound symmetry models. (This is about to get a little bit
> technical, and I might not be using exactly the right vocabulary ...)
>
> if I use a term like (1|id) in either an lmer or an lme model,
> it's equivalent to adding a term b_i to the expected value for an
> observation, where i indexes the individual (id) that observation came
> from. The b_i are assumed to be exchangeable, more specifically iid
> (independent, identically distributed), and more specifically Normally
> distributed. The variance must be non-negative, so the correlation
> among values *within* a group must be non-negative too: the
> corresponding correlation matrix is block diagonal, with one block for
> the observations in each group (observations drawn from different
> groups are independent of each other), where the off-diagonal elements
> in each block are all equal to 1>= rho >= 0.
>
> When you ask for a compound symmetric model, you now have the same
> block structure, but the value of rho may be anywhere between -1 and
> 1. I'm not actually sure how to set that up for a simple
> (intercept-only) model: the only example given in Pinheiro and Bates
> 2000 is one where we have treatments (Variety) repeated across
> different blocks. From the output below it's clear that pdCompSymm()
> forces the correlation *among treatments* to be identical for all
> treatments (it also forces the variances to be identical: in principle
> one could allow heterogeneous variances with a single fixed
> correlation).
>
> It would be possible if one wanted to badly enough to muck around
> with the correlation structures in lme4, but it's not for the faint of
> heart.
>
> The fact that you got rho=0 in your output (not shown) for lme also
> indicates that you're not setting up the model you think you are.
>
> Can you say more about the statistical/subject-area model you're
> trying to fit?
>
> Ben Bolker
>
>
>
>
> set.seed(101)
> dd <- data.frame(y=rnorm(100),f=factor(rep(1:10,10)),all=1)
> library(nlme)
> data(Oats)
> m1 <- lme(yield ~ nitro, data = Oats,
> random=list(Block=pdCompSymm(~Variety-1)))
> m2 <- lme(yield ~ nitro, data = Oats,
> random=list(Block=pdLogChol(~Variety-1)))
> library(lme4)
> m3 <- lmer(yield ~ nitro + (Variety-1|Block), data = Oats)
>
> ## lme, compound symmetry
>
> VarCorr(m1)
> Block = pdCompSymm(Variety - 1)
> Variance StdDev Corr
> VarietyGolden Rain 331.5271 18.20788
> VarietyMarvellous 331.5271 18.20788 0.635
> VarietyVictory 331.5271 18.20788 0.635 0.635
> Residual 165.5585 12.86695
>
> ## lme, unstructured
>
> > VarCorr(m2)
> Block = pdLogChol(Variety - 1)
> Variance StdDev Corr
> VarietyGolden Rain 263.7242 16.23959 VrtyGR VrtyMr
> VarietyMarvellous 232.8230 15.25854 0.594
> VarietyVictory 537.9363 23.19345 0.936 0.484
> Residual 165.5585 12.86695
>
> ## lmer, unstructured
> VarCorr(m3)
> Groups Name Std.Dev. Corr
> Block VarietyGolden Rain 16.240
> VarietyMarvellous 15.259 0.594
> VarietyVictory 23.193 0.936 0.484
> Residual 12.867
>
>
>
>
> On Mon, Jan 18, 2016 at 3:48 PM, Ashley Cohen <ashleyhcohen at gmail.com>
> wrote:
> > Hi All,
> >
> > Unfortunately, as this data isn't published yet I can't give you the
> output
> > results.
> >
> > However, the models are structured as follows:
> >
> > model1 <- lmer(outcome ~ group + visit + group *visit + (1 | id), data =
> > data)
> > model2<- lme(outcome ~group + visit + group *visit ,
> > data=data,random=~1|id, na.action = na.exclude, correlation =
> corCompSymm())
> >
> > The random effects components are identical in each:
> >
> > Random effects:
> >
> > Groups Name Variance Std.Dev.
> >
> > id (Intercept) 1.280 1.131
> >
> > Residual 1.271 1.127
> >
> > Number of obs: 195, groups: id, 59
> >
> >
> > Ashley
> >
> >
> > On Mon, Jan 18, 2016 at 3:25 PM, Thierry Onkelinx <
> thierry.onkelinx at inbo.be>
> > wrote:
> >
> >> Dear Ashley,
> >>
> >> Please add some (reproducible) examples. Without that we can only guess
> >> which models you fitted and why the results are identical. Please do add
> >> the output of sessionInfo() as well.
> >>
> >> Best regards,
> >>
> >>
> >> ir. Thierry Onkelinx
> >> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and
> >> Forest
> >> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> >> Kliniekstraat 25
> >> 1070 Anderlecht
> >> Belgium
> >>
> >> To call in the statistician after the experiment is done may be no more
> >> than asking him to perform a post-mortem examination: he may be able to
> say
> >> what the experiment died of. ~ Sir Ronald Aylmer Fisher
> >> The plural of anecdote is not data. ~ Roger Brinner
> >> The combination of some data and an aching desire for an answer does not
> >> ensure that a reasonable answer can be extracted from a given body of
> data.
> >> ~ John Tukey
> >>
> >> 2016-01-18 21:20 GMT+01:00 Ashley Cohen <ashleyhcohen at gmail.com>:
> >>
> >>> Hi Mixed Models group;
> >>>
> >>>
> >>>
> >>> I was wondering if you might be able to clarify some confusion I am
> having
> >>> with results from lmer.
> >>>
> >>>
> >>>
> >>> I am running a mixed model with fixed effects for treatment, time and
> an
> >>> interaction between treatment and time. A random effect has been added
> in
> >>> for subject (for repeated measures).
> >>>
> >>>
> >>>
> >>> I was under the impression, based on the response from Ben Bolker here,
> >>> that lmer would be estimating an unstructured correlation matrix:
> >>>
> >>>
> >>>
> http://stats.stackexchange.com/questions/86958/variance-covariance-structure-for-random-effects-in-glmer
> >>>
> >>>
> >>>
> >>> However, when I ran the model with lme, specifying a compound symmetric
> >>> correlation structure and random effect, I got the exact same results
> as
> >>> from lmer.
> >>>
> >>>
> >>>
> >>> I am not sure about what is going on and was wondering if you had any
> >>> thoughts on this. I also couldn’t find a way to pull out the
> correlation
> >>> matrix from lmer.
> >>>
> >>>
> >>>
> >>> Thank you,
> >>>
> >>>
> >>> Ashley
> >>>
> >>> [[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]]
> >
> > _______________________________________________
> > 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