[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
Mon Sep 25 21:38:07 CEST 2017
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