[R-sig-ME] lmer vs. nlme
Ben Bolker
bbolker at gmail.com
Mon Jan 18 23:02:30 CET 2016
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
More information about the R-sig-mixed-models
mailing list