[R-sig-ME] Equivalence of (0 + factor|group) and (1|group) + (1|group:factor) random effect specifications in case of compound symmetry
Maarten Jung
Maarten.Jung at mailbox.tu-dresden.de
Thu Sep 28 16:33:41 CEST 2017
@ Emmanuel Curis Douglas Bates' statement is about the "variance-covariance
matrix for the vector-valued random effects" (slide 91) not about the
covariance matrix of the observations.
2017-09-25 21:38 GMT+02:00 Maarten Jung <Maarten.Jung at mailbox.tu-dresden.de>
:
> Thanks for the very nice explanation of compound symmetry in your
> examples.
> I understand the continuum of complexity wrt to the random-effects
> structure explained by Reinhold Kliegl, also.
> But I still don't see the link between m1 and m2 in case of compound
> symmetry - sorry for that!
> Analogous to this post http://www.rpubs.com/Reinhold/22193 mod2 should
> constrain the group-related effects to equality. What exactly does this
> have to do with compound symmetry?
>
> 2017-09-25 10:15 GMT+02:00 Emmanuel Curis <emmanuel.curis at parisdescartes
> .fr>:
>
>> Hello,
>>
>> I'm not sure if I understand correctly, but I think the "compound
>> symmetry" mention usually applies to the covariance matrix of the
>> observations, Y, in a given cluster.
>>
>> In that case, and omiting the fixed part of the model that does not
>> change the Y covariance matrix, we have
>>
>> Y = fixed effects + X U + epsilon
>>
>> with U the random effect(s), X the design matrix for it, and epsilon
>> the residual error.
>>
>> Hence, assuming random effect and epsilon independant
>>
>> Sigma(Y) = tX Sigma( U ) X + Sigma( epsilon)
>>
>> so Sigma(Y) will have compound symmetry as soon as both X Sigma(U) tX
>> and Sigma( epsilon ) have themselves compound symmetry — or if they
>> have structure that kind of "compensate" themselve. So basically, that
>> does not preclude from having correlations between levels of a fixed
>> effect factor.
>>
>> However, let consider a scholar case with two observations on the same
>> patient [cluster] that takes or not its treatment [fixed effect], with
>> a different random effect for both cases, and let consider two ways of
>> writing it.
>>
>> Let Sigma(epsilon) = ( se² rse ), Sigma(U) = ( sa² r )
>> ( rse se²' ) ( r sb² )
>>
>> (typically, rse = 0, se² = se²')
>>
>> 1) Let be U1 the random effect for patient without treatment and U2
>> therandom effect for the same patient with treatment, so
>>
>> X = ( 1 0 ) = identity.
>> ( 0 1 )
>>
>> In that case, X Sigma(U) tX = Sigma(U) so
>>
>> Sigma(Y) = ( sa² + se² r + rse )
>> ( r + rse sb² + se²' )
>>
>> and you have compound symmetry if sa² + se² = sb² + se²', that is in
>> the classical case of se² = se²', assuming equal variances for both
>> the two random effects. But there is no condition on r: you can have
>> any correlation between the levels and still have compound symmetry.
>>
>> 2) U1 for patient i and U2 for effect of treatment in patient i.
>>
>> So X = ( 1 0 ), tX = ( 1 1 )
>> ( 1 1 ) ( 0 1 )
>>
>> Then Sigma(Y) = ( sa² + se² sa² + r + rse )
>> ( sa² + r + rse sa² + sb² + 2r + se²' )
>>
>> and compound symmetry is much harder to obtain for Y, you should have
>> se² = sb² + se²' + 2r — only possible for the usual model that assumes
>> se²' = se² if r = -sb²/2, which is a special case. Other correlations
>> between levels (r) will not lead to compound symmetry.
>>
>> (Hope my computations are correct).
>>
>> This suggests that the answer also depends on the way the model is
>> coded...
>>
>> Hope this helps and really corresponds to the question asked,
>>
>> Best regards,
>>
>> On Sat, Sep 23, 2017 at 07:26:53PM +0200, Maarten Jung wrote:
>> « 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]]
>> «
>> « _______________________________________________
>> « R-sig-mixed-models at r-project.org mailing list
>> « https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>> --
>> Emmanuel CURIS
>> emmanuel.curis at parisdescartes.fr
>>
>> Page WWW: http://emmanuel.curis.online.fr/index.html
>>
>
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list