Skip to content

R-sig-mixed-models Digest, Vol 136, Issue 41

8 messages · Rune Haubo, Maarten Jung

Original
Maarten Jung
# 
Dear Rune,

am I right in thinking of Model3b as a zero-correlation-parameter
model (m_zcp) but with the variances of the operators-related effects
constrained to equality?
Specifically, is the difference between Model3b and m_zcp that m_zcp
estimates variance components for each level of the factor 'machines'
and Model3b assumes equal variances across the levels of machines and
estimates only one variance for all levels?

Model3b <- lmer(Y ~ machines + (1 | operators) + (1 | machines:operators))
m_zcp  <- lmer(Y ~ machines + (machines || operators))

Cheers,
Maarten
Rune Haubo
#
On 1 May 2018 at 15:45, Maarten Jung <Maarten.Jung at mailbox.tu-dresden.de> wrote:
Hmm, well, that only makes partial sense to me. You probably want to compare

Model2 <- lmer(Y ~ machines + (1 | machines:operators))

with m_zcp. These two models have the same number of random effects
the difference
being that Model2 assumes a single variance for all of them, while m_zcp
assumes that the vectors of random effects for each operator come from a
multivariate normal distribution the dimension of which match the number of
machines.

This is probably easier to think about if we look at a concrete example, say
using the cake data from lme4. Here recipe=machines and replicate=operators;
there are 3 levels for recipe and 15 replicates.

First, '(recipe || replicate)' is the same as/expands to '(1 |
replicate) + (0 + recipe | replicate)'
which is just an over-parameterized version of '(0 + recipe | replicate)', which
again is a re-parameterized version of '(recipe | replicate)'. These are all
representing the same model (all on 10 df though lmer() is mislead and thinks
that m_zcp has 11 df):
Groups      Name        Std.Dev. Corr
 replicate   (Intercept) 0.0000
 replicate.1 recipeA     5.0692
             recipeB     6.7300   0.951
             recipeC     7.2107   0.902 0.991
 Residual                5.3622
Groups    Name    Std.Dev. Corr
 replicate recipeA 5.0692
           recipeB 6.7300   0.951
           recipeC 7.2107   0.902 0.991
 Residual          5.3622
Data: cake
Models:
m_zcp2: angle ~ recipe + (0 + recipe | replicate)
m_zcp3: angle ~ recipe + (recipe | replicate)
m_zcp: angle ~ recipe + (recipe || replicate)
       Df  AIC    BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_zcp2 10 1741 1777.0 -860.5     1721
m_zcp3 10 1741 1777.0 -860.5     1721     0      0          1
m_zcp  11 1743 1782.6 -860.5     1721     0      1          1

If we want to enforce a diagonal covariance matrix for the random effects,
we can use the dummy function:
m_zcp4 <- lmer(angle ~ recipe +
                 (0 + dummy(recipe, "A") | replicate) +
                 (0 + dummy(recipe, "B") | replicate) +
                 (0 + dummy(recipe, "C") | replicate), data=cake)
VarCorr(m_zcp4)
 Groups      Name               Std.Dev.
 replicate   dummy(recipe, "A") 5.0429
 replicate.1 dummy(recipe, "B") 6.6476
 replicate.2 dummy(recipe, "C") 7.1727
 Residual                       5.4181

Now we have something closer to what I think you were thinking of. Here,

Model2 <- lmer(angle ~ recipe + (1 | recipe:replicate), data=cake)

and m_zcp estimate the same random effects, but Model2 assumes they have
the same variance while m_zcp says that the variance depends on the
recipe and none of the models include correlation parameters.

In this case the difference in variance between recipes is small and
the random-effect estimates are very similar (as is the test for the
fixed-effects of recipe):

head(matrix(unlist(ranef(Model2)[[1]]), ncol=3))
          [,1]      [,2]        [,3]
[1,] 10.444800 13.695174 15.22126404
[2,]  9.106993 13.397883 13.43752216
[3,]  5.242219  4.181884  3.47829667
[4,] -2.487329  1.357626  4.22152245
[5,]  2.269316  1.654916 -3.21073538
[6,] -6.054813 -2.953084 -0.08918709

head(ranef(m_zcp4)$replicate)
  dummy(recipe, "A") dummy(recipe, "B") dummy(recipe, "C")
1           9.821502          13.824869        15.58456445
2           8.563530          13.524764        13.75824831
3           4.929388           4.221487         3.56131649
4          -2.338897           1.370483         4.32228155
5           2.133894           1.670588        -3.28736906
6          -5.693489          -2.981050        -0.09131581

Cheers
Rune
Maarten Jung
# 
Sorry, I forgot that lmer() (unlike lmer_alt() from the afex package)
does not convert factors to numeric covariates when using the the
double-bar notation!
The model I was talking about would be:

m_zcp5 <- lmer_alt(angle ~ recipe  + (recipe || replicate), cake)
VarCorr(m_zcp5)
 Groups      Name        Std.Dev.
 replicate   (Intercept) 6.2359
 replicate.1 re1.recipe1 1.7034
 replicate.2 re1.recipe2 0.0000
 Residual                5.3775

This model seems to differ (and I don't really understand why) from
m_zcp6 which I think is equivalent to your m_zcp4:
m_zcp6 <- lmer_alt(angle ~ recipe  + (0 + recipe || replicate), cake)
VarCorr(m_zcp6)
 Groups      Name        Std.Dev.
 replicate   re1.recipeA 5.0429
 replicate.1 re1.recipeB 6.6476
 replicate.2 re1.recipeC 7.1727
 Residual                5.4181

anova(m_zcp6, m_zcp5, refit = FALSE)
Data: data
Models:
m_zcp6: angle ~ recipe + ((0 + re1.recipeA | replicate) + (0 + re1.recipeB |
m_zcp6:     replicate) + (0 + re1.recipeC | replicate))
m_zcp5: angle ~ recipe + ((1 | replicate) + (0 + re1.recipe1 | replicate) +
m_zcp5:     (0 + re1.recipe2 | replicate))
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
m_zcp6  7 1781.8 1807.0 -883.88   1767.8
m_zcp5  7 1742.0 1767.2 -863.98   1728.0 39.807      0  < 2.2e-16 ***

Do m_zcp5 and Model3b estimate the same random effects in this case?
If not, what is the difference between m_zcp5 and Model3b (except for
the fact that the variance depends on the
recipe in m_zcp5) and which one is the more complex model?
I would be glad if you could elaborate on this and help me and the
others understand these models.

Cheers,
Maarten
On Tue, May 1, 2018 at 9:53 PM, Rune Haubo <rune.haubo at gmail.com> wrote:
Rune Haubo
#
On 2 May 2018 at 00:27, Maarten Jung <Maarten.Jung at mailbox.tu-dresden.de> wrote:
Yes, m_zcp4 and m_zcp6 are identical.

For m_zcp5 I get:
m_zcp5 <- lmer_alt(angle ~ recipe  + (recipe || replicate), cake)
VarCorr(m_zcp5)
 Groups      Name        Std.Dev.
 replicate   (Intercept) 6.0528e+00
 replicate.1 re1.recipeB 5.8203e-07
 replicate.2 re1.recipeC 2.1303e+00
 Residual                5.4693e+00

and if we change the reference level for recipe we get yet another result:
cake2 <- cake
cake2$recipe <- relevel(cake2$recipe, "C")
m_zcp5b <- lmer_alt(angle ~ recipe  + (recipe || replicate), cake2)
VarCorr(m_zcp5b)
 Groups      Name        Std.Dev.
 replicate   (Intercept) 6.5495e+00
 replicate.1 re1.recipeA 2.5561e+00
 replicate.2 re1.recipeB 1.0259e-07
 Residual                5.4061e+00
This instability indicates that something fishy is going on...

The correlation parameters are needed in the "default" representation:
(recipe | replicate) and (0 + recipe | replicate) are equivalent
because the correlation parameters make the "appropriate adjustments",
but (recipe || replicate) is _not_ the same as (0 + recipe ||
replicate) with afex::lmer_alt. I might take it as far as to say that
(recipe | replicate) is meaningful because it is a re-parameterization
of (0 + recipe | replicate). On the other hand, while the diagonal
variance-covariance matrix parameterized by (0 + recipe || replicate)
is meaningful, a model with (recipe || replicate) using afex::lmer_alt
does _not_ make sense to me (and does not represent a diagonal
variance-covariance matrix).
Well, Model3b makes sense while m_zcp5 does not, but Model3b estimates
more random effects than the others:
Model3b <- lmerTest::lmer(angle ~ recipe + (1 | replicate) + (1 |
recipe:replicate),
                          data=cake)
length(unlist(ranef(Model3b))) # 60
length(unlist(ranef(m_zcp4))) # 45 - same for m_zcp, m_zcp2 and m_zcp6
and Model2
There is no unique 'complexity' ordering, for example, Model3b use 2
random-effect variance-covariance parameters to represent 60 random
effects, while m_zcp4 (m_zcp2) use 3 (6) random-effect
variance-covariance parameters to represent 45 random effects. But
usually the relevant 'complexity' scale is the number of parameters,
cf. likelihood ratio tests, AIC, BIC etc. There are corner-cases,
however; if x1 and x2 are continuous then (1 + x1 + x2 | group) and
'(1 + x1 | group) + (1 + x2 | group)' both use 6 random-effect
variance-covariance parameters, but the models represent different
structures and you can argue that the latter formulation is less
complex than the former since it avoids the correlation between x1 and
x2.

Cheers,
Rune
Maarten Jung
# 
Thank you for explaining this. This is *very* interesting.
As far as I understand, m_zcp5 is the model Reinhold Kliegl uses in
this RPub article[1] (actually m_zcp7 which should be identical). Also
Barr et al. (2013)[2], Bates et al. (2015)[3] and Matuschek et al.
(2017)[4] suggest similar models as the first step for model
reduction. However, their categorical independent variables have only
two levels and they work with crossed random effects.

cake3 <- cake
cake3 <- subset(cake3, recipe != "C")
cake3recipe<?factor(cake3recipe<?factor(cake3recipe)
contrasts(cake3recipe) <- c(0.5, -0.5)  # Barr and Matuschek use
effect coding  m_zcp5 <- lmer_alt(angle ~ recipe  + (recipe ||
replicate), cake3) VarCorr(m_zcp5)  Groups      Name        Std.Dev.
replicate   (Intercept) 5.8077    replicate.1 re1.recipe1 1.8601
Residual                5.5366   cake3recipe) <- c(0.5, -0.5)  # Barr
and Matuschek use effect coding  m_zcp5 <- lmer_alt(angle ~ recipe  +
(recipe || replicate), cake3) VarCorr(m_zcp5)  Groups      Name
Std.Dev.  replicate   (Intercept) 5.8077    replicate.1 re1.recipe1
1.8601    Residual                5.5366   cake3recipe_numeric <-
ifelse(cake3$recipe == "A", 0.5, -0.5)
m_zcp7 <- lmer(angle ~ recipe_numeric + (1|replicate) + (0 +
recipe_numeric|replicate), cake3)
VarCorr(m_zcp7)
 Groups      Name           Std.Dev.
 replicate   (Intercept)    5.8077
 replicate.1 recipe_numeric 1.8601
 Residual                   5.5366

Besides that, Reinhold Kliegl reduces m_zcp5 to Model3b - i.e. (recipe
|| replicate) to (1 | replicate) + (1 | recipe:replicate).
Whereas you, If I understand correctly, suggest reducing/comparing (0
+ recipe || replicate) to (1 | recipe:replicate).
Why is that? Am I missing something?

Cheers,
Maarten

[1] https://rpubs.com/Reinhold/22193
[2] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3881361/
[3] https://arxiv.org/abs/1506.04967 vignettes here:
https://github.com/dmbates/RePsychLing/tree/master/vignettes
[4] https://arxiv.org/abs/1511.01864
On Wed, May 2, 2018 at 11:56 AM, Rune Haubo <rune.haubo at gmail.com> wrote:
Maarten Jung
# 
Thank you for explaining this. This is *very* interesting.
As far as I understand, m_zcp5 is the model Reinhold Kliegl uses in
this RPub article[1] (actually m_zcp7 which should be identical). Also
Barr et al. (2013)[2], Bates et al. (2015)[3] and Matuschek et al.
(2017)[4] suggest similar models as the first step for model
reduction. However, their categorical independent variables have only
two levels and they work with crossed random effects.

cake3 <- cake
cake3 <- subset(cake3, recipe != "C")
cake3$recipe <- factor(cake3$recipe)
contrasts(cake3$recipe) <- c(0.5, -0.5)  # Barr and Matuschek use effect coding
m_zcp5 <- lmer_alt(angle ~ recipe  + (recipe || replicate), cake3)
VarCorr(m_zcp5)
 Groups      Name        Std.Dev.
 replicate   (Intercept) 5.8077
 replicate.1 re1.recipe1 1.8601
 Residual                5.5366

cake3$recipe_numeric <- ifelse(cake3$recipe == "A", 0.5, -0.5)
m_zcp7 <- lmer(angle ~ recipe_numeric + (1|replicate) + (0 +
recipe_numeric|replicate), cake3)
VarCorr(m_zcp7)
 Groups      Name           Std.Dev.
 replicate   (Intercept)    5.8077
 replicate.1 recipe_numeric 1.8601
 Residual                   5.5366

Besides that, Reinhold Kliegl reduces m_zcp5 to Model3b - i.e. (recipe
|| replicate) to (1 | replicate) + (1 | recipe:replicate).
Whereas you, If I understand correctly, suggest reducing/comparing (0
+ recipe || replicate) to (1 | recipe:replicate).
Why is that? Am I missing something?

Cheers,
Maarten

[1] https://rpubs.com/Reinhold/22193
[2] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3881361/
[3] https://arxiv.org/abs/1506.04967 vignettes here:
https://github.com/dmbates/RePsychLing/tree/master/vignettes
[4] https://arxiv.org/abs/1511.01864
On Wed, May 2, 2018 at 11:56 AM, Rune Haubo <rune.haubo at gmail.com> wrote:
1 day later
Rune Haubo
#
On 2 May 2018 at 17:50, Maarten Jung <Maarten.Jung at mailbox.tu-dresden.de> wrote:
I haven't read those articles recently in enough detail that I can comment.
So m_zcp5 and m_zcp7 are identical but I don't see how they are
meaningful in this context. Looking at the random-effect design matrix

image(getME(m_zcp5, "Z")) # identical to image(getME(m_zcp7, "Z"))

you can see that this model estimates a random main effect for
replicate (i.e. (1 | replicate)) and then a random _slope_ for recipe
at each replicate (i.e. recipe in '(recipe || replicate)' is treated
as numeric rather than factor). As far as I can tell this random slope
model is _unrelated_ to models where recipe is treated as a factor
that we have discussed previously: It is a completely different model
and I don't see how it is relevant for this design. (Notice that
'recipe' is equivalent to 'recipe_numeric' in the fixed-effects, but
not so in the random-effects!)
If anything I would say that you should look at all relevant models
and choose the one that represents the best compromise between fit to
data and complexity :-) Likelihood ratio tests can be a helpful guide,
but take care not to formally compare/test models that are not nested.

Here is an example of a set of models and sequences in which they can
be compared with LR tests:

# Random main effect of replicate, no interaction:
fm1 <- lmer(angle ~ recipe + (1 | replicate), data=cake)
# Random interaction recipe:replicate; same variance across recipes;
no main effect:
fm2 <- lmer(angle ~ recipe + (1 | recipe:replicate), data=cake)
# Random interaction with different variances across recipes; no main effect:
fm3 <- lmer(angle ~ recipe +
              (0 + dummy(recipe, "A") | replicate) +
              (0 + dummy(recipe, "B") | replicate) +
              (0 + dummy(recipe, "C") | replicate), data=cake)
# Random main effect and interaction with same variance across recipes:
fm4 <- lmer(angle ~ recipe + (1 | replicate) + (1 | recipe:replicate),
data=cake)
# Random main effect and interaction with different variances across recipes:
fm5 <- lmer(angle ~ recipe + (1 | replicate) +
              (0 + dummy(recipe, "A") | replicate) +
              (0 + dummy(recipe, "B") | replicate) +
               (0 + dummy(recipe, "C") | replicate), data=cake)
# Multivariate structure that contains both main and interaction effects with
# different variances and correlations:
fm6 <- lmer(angle ~ recipe + (0 + recipe | replicate), data=cake)
# Same model, just re-parameterized:
# fm6b <- lmer(angle ~ recipe + (recipe | replicate), data=cake)
# fm6c <- lmer(angle ~ recipe + (1 | replicate) + (0 + recipe |
replicate), data=cake)
# fm6d <- lmer(angle ~ recipe + (1 | replicate) + (recipe |
replicate), data=cake)
# fm6e <- lmer(angle ~ recipe + (1 | recipe:replicate) + (recipe |
replicate), data=cake)
# anova(fm6, fm6b, fm6c, fm6d, fm6e, refit=FALSE) # fm6 = fm6b = fm6c
= fm6d = fm6e

Note that in fm4 and fm5 the random main and interaction effects are
independent, but in fm6 they are not.

No. parameters and log-likelihood/deviance of these models:
as.data.frame(anova(fm1, fm2, fm3, fm4, fm5, fm6,
refit=FALSE))[paste0("fm", 1:6), 1:5]
    Df      AIC      BIC    logLik deviance
fm1  5 1741.019 1759.011 -865.5097 1731.019
fm2  5 1776.967 1794.959 -883.4835 1766.967
fm3  7 1779.571 1804.760 -882.7857 1765.571
fm4  6 1741.067 1762.658 -864.5337 1729.067
fm5  8 1743.437 1772.224 -863.7185 1727.437
fm6 10 1741.003 1776.987 -860.5016 1721.003

The following nesting structure indicate sequences in which models can be
compares with LR tests (arrows indicate model simplification):
fm6 -> fm5 -> fm4 -> fm2
fm6 -> fm5 -> fm4 -> fm1
fm6 -> fm5 -> fm3 -> fm2

Note that fm3 and fm4 are not nested and simply represent different structures
and so no formal LR test is available. The same is true for fm1 and
fm2 as well as
fm1 and fm3.

In addition to these models there are others which are just not as
easily fitted with lmer (to the best of my knowledge) for example a
version of fm5 where the interaction random effect is specified with a
common covariance parameter on top of the 3 variances. Theoretically
there are many options here but obtaining the fits is often not
straight forward and usually no single fit is uniquely better than the
rest.

Cheers
Rune
1 day later
Maarten Jung
# 
What you suggest makes perfect sense to me.
Many thanks for your efforts and for taking time to explain the issues that
come with these models!

That said, I wonder what other mixed model experts think about m_zcp5 when
there are only categorical independent variables (i.e. factors) and
therefore I cc'd some of the aforementioned  authors. Comments on this are
highly welcome.

Cheers,
Maarten
On Fri, May 4, 2018, 09:56 Rune Haubo <rune.haubo at gmail.com> wrote: