[R-sig-ME] Equivalence of (0 + factor|group) and (1|group) + (1|group:factor) random effect specifications in case of compound symmetry
Emmanuel Curis
emmanuel.curis at parisdescartes.fr
Sat Sep 30 19:32:14 CEST 2017
Let's continue with the example given in the previous computation,
that is a 2-levels factor and the very simple case of one observation
for each of these levels for each patient (groupe). And let consider a
single group, since anyways each groups are assumed independents and
having the same random effects distribution, so the result of having
more patients will be a block-diagonal matrix, each block being the
one computed here.
I'm not completely sure of the design matrix constructed by lme4, but
let's try...
Let be Y = X Z + epsilon, with X the design matrix for the random
effects and Z the random effects themselves.
Sigma(Y) = X Sigma( Z ) tX + Sigma( epsilon )
= SigmaRE + Sigma( epsilon )
and I assume SigmaRE is what is investigated here as the "random
effects covariance matrix".
Model m1: Y ~ fixed part + (0 + factor | group)
=> there is two random effects for each group, one for each level of
factor
=> Sigma( Z ) = ( sa² r sa sb )
( r sa sb sb² )
X = ( 1 0 ) tX = X = I hence SigmaRE = Sigma( Z )
( 0 1 )
Model m2: Y ~ fixed part + (1|group) + (1|group:factor)
the (1|group) term gives a single random effect of variance sg²,
and the design matrix is X = t(1 1) because it takes the same
value for all observations in the group
=> it contributes to SigmaRE by a matrix ( sg² sg² )
( sg² sg² )
the (1|group:factor) gives a single random effect of variance sf²,
but that this time takes a different value according to the value
taken by factor in the observations of the group, and these values
are independent, so in our special example, it is equivalent to a
two-dimensional random effect of covariance matrix
SigmaGF = ( sf² 0 ) and design matrix X' = ( 1 0 )
( 0 sf² ) ( 0 1 )
so we have X SigmaGF tX = SigmaGF and finally
SigmaRE = ( sg² + sf² sg² )
( sg² sg² + sf² )
which is of compound symmetry — whereas in the first model, m1, we
had
SigmaRE = ( sa² r sa sb )
( r sa sb sb² )
which is the more general form you can imagine.
To have the two matrices equal, you must assume that 1) sb² = sa²
and 2) 0 <= r < 1 [and then sg² = r sa², sf² = sa² (1 - r)].
Since this is a special case of the model 1, you can suppose these
are nested models. But since for r it does not impose a specific
value, only an interval which is a subset of all possible values,
I'm not sure the definition of "nested models" rigorously holds,
neither that the asymptotic law is still a khi-square.
Hope this helps, and that people more familiar with nested models
theory will comment on the last point.
On Thu, Sep 28, 2017 at 04:33:41PM +0200, Maarten Jung wrote:
« @ 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
« >>
« >
« >
--
Emmanuel CURIS
emmanuel.curis at parisdescartes.fr
Page WWW: http://emmanuel.curis.online.fr/index.html
More information about the R-sig-mixed-models
mailing list