[R-sig-ME] Equivalence of (0 + factor|group) and (1|group) + (1|group:factor) random effect specifications in case of compound symmetry

[Maarten Jung]

Thanks, this really helped me understand the models!

There is one thing I still don't understand: I think that compound symmetry
implies equal variances within the factor levels and equal
covariances/correlations between the levels. If there are equal variances
within the factor levels and no correlations between the levels I
understand the equivalence of m1, m.zcp and m2.

But what if there are correlations between the levels? you explained that
m2 forces the correlation parameters to zero. But Douglas Bates states that
the models are equivalent given compound symmetry (which I think doesn't
imply zero correlations).

Just to double check: This m1 -> m.zcp -> m2 -> m.vio is a reasonable model
reduction approach, right?

2017-09-23 13:51 GMT+02:00 Reinhold Kliegl <reinhold.kliegl at gmail.com>:

> This should always be "k(k-1)/2" correlation parameters, of course.
>
> Also perhaps you may want to read this post to the list by Douglas Bates:
> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q1/001736.html
>
> On Sat, Sep 23, 2017 at 1:41 PM, Reinhold Kliegl <
> reinhold.kliegl at gmail.com> wrote:
>
>> The models are on a continuum of complexity wrt to the random-effects
>> structure. Specifically:
>>
>> m1 estimates k variance components and k(k-1) correlation parameters for
>> the k levels of factor f
>>
>> m2 estimates 1 variance component for the intercept and a 1 variance
>> component for the k-1 contrasts defined for the k levels of factor f, that
>> it constrains the k-1 contrasts for factor f to the same value. The
>> correlation parameters are forced to zero.
>>
>> The continuum becomes transparent if one re-parameterizes m1 as a m1a
>> (see below) with 1 variance component for the intercept and k-1 variance
>> components for the k levels of factor f, and k(k-1) correlation parameters.
>>  m1 and m1a have the same number of parameters and identical deviance.
>>
>> m1a <- lmer(y ~ factor + (factor|group))
>>
>> There is an additional model specification on this continuum between m1
>> and m2.  If contrasts are converted to numeric covariates, one can force
>> correlation parameters to zero, but estimate different variance components
>> for the k-1 contrasts of factor f. We call this the zero-correlation
>> parameter model.
>>
>> cB_A <- model.matrix(m1)[,2]
>> cC_A <- model.matrix(m1)[,3]
>>
>> m.zcp <- lmer(score ~ 1 + Machine + (1 + cB_A + cC_A  ||  Worker),
>>  data=Machines, REML=FALSE)
>> print(summary(m.zcp), corr=FALSE)
>>
>> Note that m.zcp without the double-bar is equivalent to m1 and m1b.
>> m.1b <- lmer(score ~ 1 + Machine + (1 + cB_A + cC_A  |  Worker),
>>  data=Machines, REML=FALSE)
>> print(summary(m.zcp), corr=FALSE)
>>
>> Of course, there is also a simpler model than m2 - the
>> varying-intercept-only model:
>>
>> m.vio <- lmer(score ~ 1 + Machine + (1 | Worker),  data=Machines,
>> REML=FALSE)
>> print(summary(m.zcp), corr=FALSE)
>>
>> Here is some R code demonstrating all this for the Machines data.
>>
>> ####
>> library(lme4)
>> #library(RePsychLing)
>>
>> data("Machines", package = "MEMSS")
>>
>> # OP m1
>> m1 <- lmer(score ~ 1 + Machine + (0 + Machine | Worker), data=Machines,
>> REML=FALSE)
>> print(summary(m1), corr=FALSE)
>>
>> # re-parameterization of m1
>> m1a <- lmer(score ~ 1 + Machine + (1 + Machine | Worker), data=Machines,
>> REML=FALSE)
>> print(summary(m1a), corr=FALSE)
>>
>> # alternative specification of m1a
>> cB_A <- model.matrix(m1)[,2]
>> cC_A <- model.matrix(m1)[,3]
>>
>> m1b <- lmer(score ~ 1 + Machine + (1 + cB_A + cC_A  | Worker),
>>  data=Machines, REML=FALSE)
>> print(summary(m1b), corr=FALSE)
>>
>> anova(m1, m1a, m1b)
>>
>> # zero-correlation parameter LMM
>> m.zcp <- lmer(score ~ 1 + Machine + (1 + cB_A + cC_A  || Worker),
>>  data=Machines, REML=FALSE)
>> print(summary(m.zcp), corr=FALSE)
>>
>> # OP m2
>> m2 <- lmer(score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker),
>> data=Machines, REML=FALSE)
>> print(summary(m2), corr=FALSE)
>>
>> # varying-intercept-only LMM
>> m.vio <- lmer(score ~ 1 + Machine + (1 | Worker), data=Machines,
>> REML=FALSE)
>> print(summary(m.vio), corr=FALSE)
>>
>> anova(m1, m.zcp, m2, m.vio)
>>
>> sessionInfo()
>> ####
>>
>> You may also want to look this RPub:
>> http://www.rpubs.com/Reinhold/22193
>>
>> On Sat, Sep 23, 2017 at 11:43 AM, Maarten Jung <
>> Maarten.Jung at mailbox.tu-dresden.de> wrote:
>>
>>> Hello everyone,
>>>
>>> I have a question regarding the equivalence of the following models:
>>>
>>> m1 <- lmer(y ~ factor + (0 + factor|group))
>>> m2 <- lmer(y ~ factor + (1|group) + (1|group:factor))
>>>
>>> Douglas Bates states (slide 91 in this presentation [1])  that these
>>> models
>>> are equivalent in case of compound symmetry.
>>>
>>> 1. I realized that I don't really understand the random slope by factor
>>> model (m1) and espacially why it reduces to m2 given compound symmetry.
>>> Also, why is there no random intercept in m1?
>>> Can anyone explain the difference between the models and how m1 reduces
>>> to
>>> m2 in an intuitive way.
>>>
>>> 2. If m1 is a special case of m2 – this could be an interesting option
>>> for
>>> model reduction but I’ve never seen something like m2 in papers. Instead
>>> they suggest dropping the random slope and thus the interaction
>>> completely
>>> (e.g. Matuschek et al. 2017 [3]).
>>> What do you think about it?
>>>
>>> Please note that I asked the question on Stack Exchange [2] but some
>>> consider it as off-topic. So, I hope one of you can help me out.
>>>
>>>
>>> Best regards,
>>> Maarten
>>>
>>> [1] http://www.stat.wisc.edu/~bates/UseR2008/WorkshopD.pdf
>>> [2] https://stats.stackexchange.com/q/304374/136579
>>> [3] https://doi.org/10.1016/j.jml.2017.01.001
>>>
>>>         [[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]]