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

On 1 May 2018 at 15:45, Maarten Jung <Maarten.Jung at mailbox.tu-dresden.de> wrote:

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

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):

library(lmerTest)
m_zcp <- lmer(angle ~ recipe  + (recipe || replicate), cake)
VarCorr(m_zcp)

 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

m_zcp2 <- lmer(angle ~ recipe + (0 + recipe | replicate), cake)
VarCorr(m_zcp2)

 Groups    Name    Std.Dev. Corr
 replicate recipeA 5.0692
           recipeB 6.7300   0.951
           recipeC 7.2107   0.902 0.991
 Residual          5.3622

m_zcp3 <- lmer(angle ~ recipe + (recipe | replicate), cake)
anova(m_zcp, m_zcp2, m_zcp3, refit=FALSE)

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