[R-sig-ME] R-sig-mixed-models Digest, Vol 136, Issue 41

[Nicolas Flaibani]

Hi everyone!


I would like to join the discussion of the second topic (Specifying Multiple Random Effects in NLME by Dan) but in my case I’m going to talk about the lme4 package, instead of the nmle. I think that the question is still relevant.

I’ve had a doubt for a long time about how to investigate the interactions between random and fixed effects. I’ve read a lot of forums, papers and help’s packages and I’ve always concluded that the correct form of testing the interaction between a random and a fixed variable was:


Model1 <- lmer (Y ~ X1 + X2 + (X1 | Random Variable)


However, I found in some forums and personal communications from several researchers that there is another way to investigate the interaction between random and fixed variables and has the following syntax:


Model2 <- lmer (Y ~ X1 + X2 + (1 | Random Variable: X1)


I understand that this syntax


(1|Random Variable/X1) = (1|Random Variable)+(1|Random Variable:X1)


specify a nested structure between the variables and this is not the case of interest.

My particular question is whether the syntax of the Model 2 is correct to test interactions between random and fixed variables. If this model is correct, which are the differences with the syntax of Model 1, since the resulting models are clearly different? Besides, in coincidence with the question of Dan (“Are there cases where one vs. the other formulation should absolutely be used? My understanding that for continuous variables, e.g., multiple measurements across multiple days, Days|Object would be the correct syntax. But here we're talking about a factor variable”), I ask if one type of syntax should be used if the fixed variables are continuous or there are factors.

If I compare the summary from a model with the latter syntax (model 2), with the summary of the same analysis made with a statistic program (like Statistica), the results are very similar. That’s not the case with the model 1.

For example, if I analyze a morphological trait with the syntax


M2 <- lmer (Wing ~ Temperature * Sex + Temperature + Sex + (1 | Line) + (1 | Line:Sex:Temperature) + (1 | Line:Sex) + (1 | Line:Temperature))


the summary is the following:


Random effects:

Groups Name Variance Std.Dev.

Line:Sex:Temperature (Intercept) 14.6231 3.8240

Line:Temperature (Intercept) 154.7685 12.4406

Line:Sex (Intercept) 0.6947 0.8335

Line (Intercept) 72.5945 8.5202

Residual 180.0664 13.4189


Fixed effects:

Estimate Std. Error df t value Pr(>|t|)

(Intercept) 501.141 2.268 96.940 221.009 < 2e-16 ***

Temperature25 -57.960 2.699 54.800 -21.473 < 2e-16 ***

SexM -53.639 1.001 96.260 -53.584 < 2e-16 ***

Temperature25:SexM -6.488 1.391 48.300 -4.663 2.49e-05 ***


I found that the function rand() from the lmerTest package gives me the p values of the random effects if I write the model like this:

> rand(M2)

Analysis of Random effects Table:


Chi.sq Chi.DF p.value

Line 4.6152 1 0.03 *

Line:Sex:Temperature 30.8130 1 3e-08 ***

Line:Sex 0.0391 1 0.84

Line:Temperature 112.1539 1 <2e-16 ***


I don’t know if this function is reliable because it is not mentioned for testing the significance of the random effects in the page https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html

The summary of the same analysis made with the Statistica is:




Effect

SS

Degr. of

MS

Den.Syn.

Den.Syn.

F

p

Intercept

Fixed

764579151

1

764579151

48,001

12655,91

60412,83

0,000000

Line

Random

608181

48

12670

47,800

6489,64

1,95

0,011254

Sex

Fixed

3138661

1

3138661

48,038

495,62

6332,81

0,000000

Temperature

Fixed

3686660

1

3686660

48,003

6459,62

570,72

0,000000

Line*Sex

Random

23808

48

496

48,000

473,30

1,05

0,435866

Line*Temperature

Random

310413

48

6467

48,000

473,30

13,66

0,000000

Sex*Temperature

Fixed

10075

1

10075

48,040

472,94

21,30

0,000029

Line*Sex*Temperature

Random

22718

48

473

3696,000

167,33

2,83

0,000000

Error

618467

3696

167




But if I write the model with the other syntax:

M1 <- lmer(Wing ~ Temperature * Sex + (Temperature * Sex | Line))



the summary is the following:



REML criterion at convergence: 31440.9



Random effects:

Groups Name Variance Std.Dev. Corr

Line (Intercept) 266.78 16.333

Temperature25 398.27 19.957 -0.60

SexM 41.54 6.446 -0.56 0.46

Temperature25:SexM 61.34 7.832 0.56 -0.61 -0.80

Residual 167.33 12.936



Fixed effects:

Estimate Std. Error df t value Pr(>|t|)

(Intercept) 501.603 2.371 48.046 211.586 < 2e-16 ***

Temperature25 -58.423 2.911 48.027 -20.070 < 2e-16 ***

SexM -53.659 1.095 47.964 -49.023 < 2e-16 ***

Temperature 25:SexM -6.470 1.393 48.278 -4.644 2.66e-05 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



In addition, if I apply the ‘rand functionŽ for this syntax (M1), it doesn’t retourn the whole p-values of the random effects (Do not give me p value for line and line*Temperature*Sex)

Analysis of Random effects Table:

Chi.sq Chi.DF p.value

Temperatura:Línea 0.00e+00 0 1

Sexo:Línea 1.46e-10 0 <2e-16 ***


I really appreciate your time and dedication for answering this questions. Thank you for trying to help us understand a little more about the syntax of these complex models and thus better understand their correct approach.

Thank you very much for your time everyone.



Greetings,

Nicolas


----------------------------------------------------------------------------------------
Message: 2
Date: Wed, 25 Apr 2018 16:11:38 -0400
From: Dan <***@gmail.com>
To: "R-SIG-Mixed-***@R-project.org"
<R-sig-mixed-***@r-project.org>, Ben Bolker <***@gmail.com>
Subject: [R-sig-ME] Specifying Multiple Random Effects in NLME
Message-ID:
<CAET4i1f-SH5zA6xcCxNXc091HCk+snMv+***@mail.gmail.com>
Content-Type: text/plain; charset="utf-8"

Hi all:

I am curating an off-list thread about specifying multiple random effects
in NLME.

1. If it's (1|Object) + (1|Object:Coating) that you want then
you should be able to use a nested specification (which nlme *can*
handle relatively easily), i.e. something like

random=a+b+c~1|Object/Coating


Although (Coating|Object) and (1|Object:Coating) both in some sense
represent "interactions" the latter is *much* simpler/more parsimonious.

If you're OK with 1|Object:Coating rather than Coating|Object it
should be *much* faster. If you don't understand the distinction (which
would be absolutely fine and understandable) can we resume the
discussion on r-sig-mixed-models ... ?
-----------

So:
random=a+b+c~Coating|Object
does not fit.

But:
random=a+b+c~Object/Coating
fits.

Can you better explain the distinction here? I have sometimes used the
1|Object:Coating + 1|Object syntax and sometimes the Coating|Object syntax
in other models. My experience/understanding is that the former syntax with
multiple "within subject" variables produces exactly matching output to the
standard "repeated measures ANOVA" with the lmer assumption of compound
symmetry.

Are there cases where one vs. the other formulation should absolutely be
used? My understanding that for continuous variables, e.g., multiple
measurements across multiple days, Days|Object would be the correct syntax.
But here we're talking about a factor variable.


2. I'm trying to read the "random" section for nlme right now but it's
kind of making my head explode (and I think there's a typo: " the same
as the order of the order of the elements in the list"). It *sounds*
like (1) explicitly creating an interaction
ObjCoating=interaction(Object,Coating) and (2) using something like

list(ObjCoating=a~1,Object=b~1,Object=c~1)

should work (grouping factors as names, then [right-hand-side variable
name]~[random effects model], but I'm worried about the phrase "The
order of nesting will be assumed the same as the order of the elements
in the list": what nesting?
-----------

I think that formulation is explicitly in order. I replaced your first
ObjCoating with simply Object, just to test what would happen:

Random effects:
Formula: a ~ 1 | Object
a.(Intercept)
StdDev: 1.305816

Formula: b ~ 1 | Object %in% Object
b.(Intercept)
StdDev: 0.01576521

Formula: c ~ 1 | Object %in% Object %in% Object
c.(Intercept) Residual
StdDev: 2.677883 2.219676

[[alternative HTML version deleted]]




*****

[[alternative HTML version deleted]]

[Rune Haubo]

On 27 April 2018 at 22:33, Nicolas Flaibani <***@hotmail.com> wrote:
>
> Hi everyone!
>
>
> I would like to join the discussion of the second topic (Specifying Multiple Random Effects in NLME by Dan) but in my case I’m going to talk about the lme4 package, instead of the nmle. I think that the question is still relevant.
>
> I’ve had a doubt for a long time about how to investigate the interactions between random and fixed effects. I’ve read a lot of forums, papers and help’s packages and I’ve always concluded that the correct form of testing the interaction between a random and a fixed variable was:
>
>
> Model1 <- lmer (Y ~ X1 + X2 + (X1 | Random Variable)
>
>
> However, I found in some forums and personal communications from several researchers that there is another way to investigate the interaction between random and fixed variables and has the following syntax:
>
>
> Model2 <- lmer (Y ~ X1 + X2 + (1 | Random Variable: X1)
>
The 'classical' (think DOE and 'variance-components') interaction is
Model2 if X1 is categorical/factor and Model1 if X1 is continuous
(then Model1 is sometimes called the 'random coefficient model').

If X1 is continuous Model2 doesn't make sense - on the other hand
Model1 can be fitted if X1 is a factor. In that case Model1 is rather
complex and the number of parameters grows rapidly with the number of
levels of X1 - in comparison the random term in Model2 uses just one
(variance) parameter. While some people often favour 'Model1 with
factor X1'- construction, I often think it is difficult to explain
what this model actually means; it is also very often overfitting
since it requires a lot of data and special data-generating mechanism
to support models of that form.

>
> I understand that this syntax
>
>
> (1|Random Variable/X1) = (1|Random Variable)+(1|Random Variable:X1)
>
>
> specify a nested structure between the variables and this is not the case of interest.

This construction serves two purposes: one is when the random-effect
variables have a 'truly' nested structure such as pupils in classes
(in schools, in districts, etc). The other purpose often applies to
designed experiments where you might have machines (fixed) and
operators (random). The main effects model is then

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

and a model that includes the interaction reads

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

Technically the 'machines:operators' combinations are nested in
'operators' but we usually don't think of it that way.

The point here is that we need to also consider Model3b as an
alternative to Model1 and Model2. Of these models, Model3 _and_ Model2
are the simplest while Model1 is the most complex with Model3b in the
middle often serving as an appropriate compromise.

>
> My particular question is whether the syntax of the Model 2 is correct to test interactions between random and fixed variables. If this model is correct, which are the differences with the syntax of Model 1, since the resulting models are clearly different? Besides, in coincidence with the question of Dan (“Are there cases where one vs. the other formulation should absolutely be used? My understanding that for continuous variables, e.g., multiple measurements across multiple days, Days|Object would be the correct syntax. But here we're talking about a factor variable”), I ask if one type of syntax should be used if the fixed variables are continuous or there are factors.
>
> If I compare the summary from a model with the latter syntax (model 2), with the summary of the same analysis made with a statistic program (like Statistica), the results are very similar. That’s not the case with the model 1.
>
> For example, if I analyze a morphological trait with the syntax
>
>
> M2 <- lmer (Wing ~ Temperature * Sex + Temperature + Sex + (1 | Line) + (1 | Line:Sex:Temperature) + (1 | Line:Sex) + (1 | Line:Temperature))
>
>
> the summary is the following:
>
>
> Random effects:
>
> Groups Name Variance Std.Dev.
>
> Line:Sex:Temperature (Intercept) 14.6231 3.8240
>
> Line:Temperature (Intercept) 154.7685 12.4406
>
> Line:Sex (Intercept) 0.6947 0.8335
>
> Line (Intercept) 72.5945 8.5202
>
> Residual 180.0664 13.4189
>
>
> Fixed effects:
>
> Estimate Std. Error df t value Pr(>|t|)
>
> (Intercept) 501.141 2.268 96.940 221.009 < 2e-16 ***
>
> Temperature25 -57.960 2.699 54.800 -21.473 < 2e-16 ***
>
> SexM -53.639 1.001 96.260 -53.584 < 2e-16 ***
>
> Temperature25:SexM -6.488 1.391 48.300 -4.663 2.49e-05 ***
>
>
> I found that the function rand() from the lmerTest package gives me the p values of the random effects if I write the model like this:
>
>> rand(M2)
>
> Analysis of Random effects Table:
>
>
> Chi.sq Chi.DF p.value
>
> Line 4.6152 1 0.03 *
>
> Line:Sex:Temperature 30.8130 1 3e-08 ***
>
> Line:Sex 0.0391 1 0.84
>
> Line:Temperature 112.1539 1 <2e-16 ***
>
>
> I don’t know if this function is reliable because it is not mentioned for testing the significance of the random effects in the page https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html

Having written rand (=ranova) my views may be biased but I should say
that it _is_ reliable. That is, I haven't seen cases where it's not
doing what was intended ;-) rand() and ranova() simply compares two
models using anova(m1, m2, refit=FALSE) where m2 differs from m1 by
having one of its random-effect terms reduced or removed.

>
> The summary of the same analysis made with the Statistica is:
>
>
>
>
> Effect
>
> SS
>
> Degr. of
>
> MS
>
> Den.Syn.
>
> Den.Syn.
>
> F
>
> p
>
> Intercept
>
> Fixed
>
> 764579151
>
> 1
>
> 764579151
>
> 48,001
>
> 12655,91
>
> 60412,83
>
> 0,000000
>
> Line
>
> Random
>
> 608181
>
> 48
>
> 12670
>
> 47,800
>
> 6489,64
>
> 1,95
>
> 0,011254
>
> Sex
>
> Fixed
>
> 3138661
>
> 1
>
> 3138661
>
> 48,038
>
> 495,62
>
> 6332,81
>
> 0,000000
>
> Temperature
>
> Fixed
>
> 3686660
>
> 1
>
> 3686660
>
> 48,003
>
> 6459,62
>
> 570,72
>
> 0,000000
>
> Line*Sex
>
> Random
>
> 23808
>
> 48
>
> 496
>
> 48,000
>
> 473,30
>
> 1,05
>
> 0,435866
>
> Line*Temperature
>
> Random
>
> 310413
>
> 48
>
> 6467
>
> 48,000
>
> 473,30
>
> 13,66
>
> 0,000000
>
> Sex*Temperature
>
> Fixed
>
> 10075
>
> 1
>
> 10075
>
> 48,040
>
> 472,94
>
> 21,30
>
> 0,000029
>
> Line*Sex*Temperature
>
> Random
>
> 22718
>
> 48
>
> 473
>
> 3696,000
>
> 167,33
>
> 2,83
>
> 0,000000
>
> Error
>
> 618467
>
> 3696
>
> 167
>
>
>
>
> But if I write the model with the other syntax:
>
> M1 <- lmer(Wing ~ Temperature * Sex + (Temperature * Sex | Line))
>
Writing the model as
M1 <- lmer(Wing ~ Temperature * Sex + (0 + Temperature:Sex | Line))
is often preferable as it makes the random-effect variance-covariance
matrix easier to interpret.

>
>
> the summary is the following:
>
>
>
> REML criterion at convergence: 31440.9
>
>
>
> Random effects:
>
> Groups Name Variance Std.Dev. Corr
>
> Line (Intercept) 266.78 16.333
>
> Temperature25 398.27 19.957 -0.60
>
> SexM 41.54 6.446 -0.56 0.46
>
> Temperature25:SexM 61.34 7.832 0.56 -0.61 -0.80
>
> Residual 167.33 12.936
>
>
>
> Fixed effects:
>
> Estimate Std. Error df t value Pr(>|t|)
>
> (Intercept) 501.603 2.371 48.046 211.586 < 2e-16 ***
>
> Temperature25 -58.423 2.911 48.027 -20.070 < 2e-16 ***
>
> SexM -53.659 1.095 47.964 -49.023 < 2e-16 ***
>
> Temperature 25:SexM -6.470 1.393 48.278 -4.644 2.66e-05 ***
>
> ---
>
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
> In addition, if I apply the ‘rand function´ for this syntax (M1), it doesn’t retourn the whole p-values of the random effects (Do not give me p value for line and line*Temperature*Sex)

Use lmerTest::ranova(M1, reduce.terms=FALSE) to achieve a test for the
entire random-effect term (and make sure you have lmerTest version >=
3.0-0 installed).

>
> Analysis of Random effects Table:
>
> Chi.sq Chi.DF p.value
>
> Temperatura:Línea 0.00e+00 0 1
>
> Sexo:Línea 1.46e-10 0 <2e-16 ***
>
>
> I really appreciate your time and dedication for answering this questions. Thank you for trying to help us understand a little more about the syntax of these complex models and thus better understand their correct approach.

You are not alone if you think this is complicated. My students are
for the most part challenged in simply writing up the mathematical
formula that correctly represents models such as Model1 -
transitioning from scalar random-effect terms (e.g. Model3b) to
vector-valued random-effect terms (e.g. Model1) often takes time and
dedication.

Cheers
Rune

>
> Thank you very much for your time everyone.
>
>
>
> Greetings,
>
> Nicolas
>
>
> ----------------------------------------------------------------------------------------
> Message: 2
> Date: Wed, 25 Apr 2018 16:11:38 -0400
> From: Dan <***@gmail.com>
> To: "R-SIG-Mixed-***@R-project.org"
> <R-sig-mixed-***@r-project.org>, Ben Bolker <***@gmail.com>
> Subject: [R-sig-ME] Specifying Multiple Random Effects in NLME
> Message-ID:
> <CAET4i1f-SH5zA6xcCxNXc091HCk+snMv+***@mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
>
> Hi all:
>
> I am curating an off-list thread about specifying multiple random effects
> in NLME.
>
> 1. If it's (1|Object) + (1|Object:Coating) that you want then
> you should be able to use a nested specification (which nlme *can*
> handle relatively easily), i.e. something like
>
> random=a+b+c~1|Object/Coating
>
>
> Although (Coating|Object) and (1|Object:Coating) both in some sense
> represent "interactions" the latter is *much* simpler/more parsimonious.
>
> If you're OK with 1|Object:Coating rather than Coating|Object it
> should be *much* faster. If you don't understand the distinction (which
> would be absolutely fine and understandable) can we resume the
> discussion on r-sig-mixed-models ... ?
> -----------
>
> So:
> random=a+b+c~Coating|Object
> does not fit.
>
> But:
> random=a+b+c~Object/Coating
> fits.
>
> Can you better explain the distinction here? I have sometimes used the
> 1|Object:Coating + 1|Object syntax and sometimes the Coating|Object syntax
> in other models. My experience/understanding is that the former syntax with
> multiple "within subject" variables produces exactly matching output to the
> standard "repeated measures ANOVA" with the lmer assumption of compound
> symmetry.
>
> Are there cases where one vs. the other formulation should absolutely be
> used? My understanding that for continuous variables, e.g., multiple
> measurements across multiple days, Days|Object would be the correct syntax.
> But here we're talking about a factor variable.
>
>
> 2. I'm trying to read the "random" section for nlme right now but it's
> kind of making my head explode (and I think there's a typo: " the same
> as the order of the order of the elements in the list"). It *sounds*
> like (1) explicitly creating an interaction
> ObjCoating=interaction(Object,Coating) and (2) using something like
>
> list(ObjCoating=a~1,Object=b~1,Object=c~1)
>
> should work (grouping factors as names, then [right-hand-side variable
> name]~[random effects model], but I'm worried about the phrase "The
> order of nesting will be assumed the same as the order of the elements
> in the list": what nesting?
> -----------
>
> I think that formulation is explicitly in order. I replaced your first
> ObjCoating with simply Object, just to test what would happen:
>
> Random effects:
> Formula: a ~ 1 | Object
> a.(Intercept)
> StdDev: 1.305816
>
> Formula: b ~ 1 | Object %in% Object
> b.(Intercept)
> StdDev: 0.01576521
>
> Formula: c ~ 1 | Object %in% Object %in% Object
> c.(Intercept) Residual
> StdDev: 2.677883 2.219676
>
> [[alternative HTML version deleted]]
>
>
>
>
> *****
>
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-mixed-***@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

[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 <***@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

[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 <***@gmail.com> wrote:
> On 1 May 2018 at 15:45, Maarten Jung <***@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

[Rune Haubo]

On 2 May 2018 at 00:27, Maarten Jung <***@mailbox.tu-dresden.de> wrote:
> 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 ***
>

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

> Do m_zcp5 and Model3b estimate the same random effects in this case?

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

> 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?

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

> I would be glad if you could elaborate on this and help me and the
> others understand these models.
>
> Cheers,
> Maarten
>

[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 <***@gmail.com> wrote:
> On 2 May 2018 at 00:27, Maarten Jung <***@mailbox.tu-dresden.de> wrote:
>> 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 ***
>>
>
> 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).
>
>> Do m_zcp5 and Model3b estimate the same random effects in this case?
>
> 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
>
>> 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?
>
> 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
>
>> I would be glad if you could elaborate on this and help me and the
>> others understand these models.
>>
>> Cheers,
>> Maarten
>>

[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 <***@gmail.com> wrote:
> On 2 May 2018 at 00:27, Maarten Jung <***@mailbox.tu-dresden.de> wrote:
>> 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 ***
>>
>
> 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).
>
>> Do m_zcp5 and Model3b estimate the same random effects in this case?
>
> 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
>
>> 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?
>
> 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
>
>> I would be glad if you could elaborate on this and help me and the
>> others understand these models.
>>
>> Cheers,
>> Maarten
>>

[Rune Haubo]

On 2 May 2018 at 17:50, Maarten Jung <***@mailbox.tu-dresden.de> wrote:
> 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.

I haven't read those articles recently in enough detail that I can comment.

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

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

>
> 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?

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

>
> 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 <***@gmail.com> wrote:
>> On 2 May 2018 at 00:27, Maarten Jung <***@mailbox.tu-dresden.de> wrote:
>>> 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 ***
>>>
>>
>> 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).
>>
>>> Do m_zcp5 and Model3b estimate the same random effects in this case?
>>
>> 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
>>
>>> 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?
>>
>> 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
>>
>>> I would be glad if you could elaborate on this and help me and the
>>> others understand these models.
>>>
>>> Cheers,
>>> Maarten
>>>

[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 <***@gmail.com> wrote:

> On 2 May 2018 at 17:50, Maarten Jung <***@mailbox.tu-dresden.de>
> wrote:
> > 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.
>
> I haven't read those articles recently in enough detail that I can comment.
>
> >
> > 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
>
> 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!)
>
> >
> > 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?
>
> 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
>
> >
> > 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 <***@gmail.com>
> wrote:
> >> On 2 May 2018 at 00:27, Maarten Jung <
> ***@mailbox.tu-dresden.de> wrote:
> >>> 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 ***
> >>>
> >>
> >> 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).
> >>
> >>> Do m_zcp5 and Model3b estimate the same random effects in this case?
> >>
> >> 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
> >>
> >>> 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?
> >>
> >> 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
> >>
> >>> I would be glad if you could elaborate on this and help me and the
> >>> others understand these models.
> >>>
> >>> Cheers,
> >>> Maarten
> >>>
>

[[alternative HTML version deleted]]

[Nicolas Flaibani]

Hi, Rune,
Thank you, a lot, for answering so fast to my question (R-sig-mixed-models Digest, Vol 136, Issue 41). These days I’ve been reading the mails that you have been sending about how to investigate the interaction between random and fixed effects and how to estimate the variance components in these models.
To put the question and examples in context, I want to briefly comment the research subject of my laboratory. The working line that is developed here study different characters (morphological, physiological and behavioral) in several species of Drosophila. Because of its easy way of maintenance, Drosophila is perfect to generate isolines and have an accurate estimate of genetics components. For example, one of the working lines involucrates the study of the thorax or wing length (like estimators of body size) in fly raised at different temperatures for diverse genetic lines. One of the principal interests is the study of the interaction between the fixed effects variables (i. e. Temperature) and the random effect variable (Line) to study the genotype-phenotype interaction. In this sense, to further compare different populations we are interested in estimate the variance components (by percentages of variance explained by Line and the interaction Line*Temperature, qualitatively).

As you said “…Model1 can be fitted if X1 is a factor. In that case Model1 is rather complex and the number of parameters grows rapidly with the number of levels of X1 - in comparison the random term in Model2 uses just one (variance) parameter”, in the first place I am going to simplify the model to try to better understand.
So, as the first option, I decided to try a model like the Model 1:
mt <-lmer(Torax ~Temperature+Sex+ (Temperature|Line), data)
Random effects:
Groups Name Variance Std.Dev. Corr
Line (Intercept) 49.27 7.019
Temperature25 52.80 7.266 -0.55
Residual 19.21 4.383
Number of obs: 780, groups: Linea, 49

In addition, trying to resume my interest in knowing which is the variance explained by the interaction Line*Temperature and by Line, I’ve arrived to these models:
mtb <- lmer(Torax ~Temperature +Sex+ (0+Temperature|Line), data) AIC=4814,43
Random effects:
Groups Name Variance Std.Dev. Corr
Line Temperature17 49.27 7.019
Temperature25 46.38 6.811 0.45
Residual 19.21 4.383
Number of obs: 780, groups: Linea, 49

mtb2 <- lmer(Torax ~Temperature +Sex+ (Temperature||Line), data) AIC=4816,43
Random effects:
Groups Name Variance Std.Dev. Corr
Line (Intercept) 22.41 4.734
Line1 Temperature17 26.86 5.182
Temperature25 23.97 4.896 -0.04
Residual 19.21 4.383
Number of obs: 780, groups: Linea, 49
> AIC(mt,mtb,mtb2)
df AIC
mt 7 4814.428
mtb 7 4814.428
mtb2 8 4816.428

They have the same fixed effects, but the difference resides in how they expressed the variance components. So, if I’m interested in obtain the whole variance components (Line and Line*Temperature), do I need to run the mtb2 model? The other models aren’t useful? The difference between the three models is how we interpret and how they show us the variance?
On one hand, the difference between mtb and mtb2 is the existence of one parameter of variance related to Line? (If I want to estimate the variance for each component it will be Line=22.41 and Line:Temperature(26.86+23.97) or this last sum doesnŽt make sense?)
It's very difficult for me to understand how to interpret the variance components in the models mt and mtb. In addition, in your first answer to Jung, you said “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)…” I don’t understand why you said that the model (recipe || replicate) is over-parameterized, if it has the same parameters that (0 + recipe | replicate). Even if I notice that in the variance components I have different amounts of estimated variances, however I would not be seeing that this has an impact on the model worsening due to excess of estimated parameters. I do not know if I explain myself…
On the other hand, going back to your initial advice “The point here is that we need to also consider Model3b as an alternative to Model1 and Model2. Of these models, Model3 _and_ Model2 are the simplest while Model1 is the most complex with Model3b in the middle often serving as an appropriate compromise..” I run a model as
mt3 <- lmer(Torax ~ Sex +Temperature +(1|Line:Temperature),
AIC=4820.112
Random effects:
Groups Name Variance Std.Dev.
Line:Temperature (Intercept) 47.81 6.915
Residual 19.21 4.383
Number of obs: 780, groups: Line:Temperature, 98

However, this syntax doesn’t let me see all variance components (doesnŽt show Line). So, I decided to run another syntax.
mt3b <- lmer(Torax ~Sex +Temperature+ (1 | Line) +(1 | Temperature:Line),data)
Random effects:
Groups Name Variance Std.Dev.
Temperature:Line (Intercept) 26.40 5.138
Line (Intercept) 21.42 4.628
Residual 19.21 4.383
Number of obs: 780, groups: Temperature:Line, 98; Linea, 49

If I compare all models:
df AIC
mt 7 4814.428
mtb 7 4814.428
mtb2 8 4816.428
mt3 5 4820.112
mt3b 6 4812.476

and only I focus in those that show me all variance components (Line and Line*Temperature, mtb2 and mt3b), can I choose the best model in function of AIC?
Going back to the general objective of estimate the variance components and evaluate the interaction between the random variables, my question is: the difference between the last models is that mtb2 is estimating 2 parameters more which correspond, on one hand, to the discrimination of the interaction’s variance of Line*Temperature in two components (Line*Temperature17 and Line*Temperature25) and on the other hand, to the correlation between those components. If this is correct, according to my objective, is more parsimonious mt3b? In adittion, the sum of the variance explained for the interaction Line*Temperature17 and Line*Temperature25 (if this make any sense) shouldn’t be similar or identical to the variance explained for Temperature:Line? Because if I compare those values I don’t see that happen.

Finally, in your last mail you enumerate different models that can be compared. I don’t understand the difference between the model fm5 and the model fm6c, I don’t know if you can clarify me a bit.

Thank you so much again for your patience and dedication.
Cheers,
Nicolás.


________________________________
De: Rune Haubo <***@gmail.com>
Enviado: sábado, 28 de abril de 2018 06:32
Para: Nicolas Flaibani
Cc: r-sig-mixed-***@r-project.org
Asunto: Re: [R-sig-ME] R-sig-mixed-models Digest, Vol 136, Issue 41

On 27 April 2018 at 22:33, Nicolas Flaibani <***@hotmail.com> wrote:
>
> Hi everyone!
>
>
> I would like to join the discussion of the second topic (Specifying Multiple Random Effects in NLME by Dan) but in my case I’m going to talk about the lme4 package, instead of the nmle. I think that the question is still relevant.
>
> I’ve had a doubt for a long time about how to investigate the interactions between random and fixed effects. I’ve read a lot of forums, papers and help’s packages and I’ve always concluded that the correct form of testing the interaction between a random and a fixed variable was:
>
>
> Model1 <- lmer (Y ~ X1 + X2 + (X1 | Random Variable)
>
>
> However, I found in some forums and personal communications from several researchers that there is another way to investigate the interaction between random and fixed variables and has the following syntax:
>
>
> Model2 <- lmer (Y ~ X1 + X2 + (1 | Random Variable: X1)
>
The 'classical' (think DOE and 'variance-components') interaction is
Model2 if X1 is categorical/factor and Model1 if X1 is continuous
(then Model1 is sometimes called the 'random coefficient model').

If X1 is continuous Model2 doesn't make sense - on the other hand
Model1 can be fitted if X1 is a factor. In that case Model1 is rather
complex and the number of parameters grows rapidly with the number of
levels of X1 - in comparison the random term in Model2 uses just one
(variance) parameter. While some people often favour 'Model1 with
factor X1'- construction, I often think it is difficult to explain
what this model actually means; it is also very often overfitting
since it requires a lot of data and special data-generating mechanism
to support models of that form.

>
> I understand that this syntax
>
>
> (1|Random Variable/X1) = (1|Random Variable)+(1|Random Variable:X1)
>
>
> specify a nested structure between the variables and this is not the case of interest.

This construction serves two purposes: one is when the random-effect
variables have a 'truly' nested structure such as pupils in classes
(in schools, in districts, etc). The other purpose often applies to
designed experiments where you might have machines (fixed) and
operators (random). The main effects model is then

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

and a model that includes the interaction reads

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

Technically the 'machines:operators' combinations are nested in
'operators' but we usually don't think of it that way.

The point here is that we need to also consider Model3b as an
alternative to Model1 and Model2. Of these models, Model3 _and_ Model2
are the simplest while Model1 is the most complex with Model3b in the
middle often serving as an appropriate compromise.

>
> My particular question is whether the syntax of the Model 2 is correct to test interactions between random and fixed variables. If this model is correct, which are the differences with the syntax of Model 1, since the resulting models are clearly different? Besides, in coincidence with the question of Dan (“Are there cases where one vs. the other formulation should absolutely be used? My understanding that for continuous variables, e.g., multiple measurements across multiple days, Days|Object would be the correct syntax. But here we're talking about a factor variable”), I ask if one type of syntax should be used if the fixed variables are continuous or there are factors.
>
> If I compare the summary from a model with the latter syntax (model 2), with the summary of the same analysis made with a statistic program (like Statistica), the results are very similar. That’s not the case with the model 1.
>
> For example, if I analyze a morphological trait with the syntax
>
>
> M2 <- lmer (Wing ~ Temperature * Sex + Temperature + Sex + (1 | Line) + (1 | Line:Sex:Temperature) + (1 | Line:Sex) + (1 | Line:Temperature))
>
>
> the summary is the following:
>
>
> Random effects:
>
> Groups Name Variance Std.Dev.
>
> Line:Sex:Temperature (Intercept) 14.6231 3.8240
>
> Line:Temperature (Intercept) 154.7685 12.4406
>
> Line:Sex (Intercept) 0.6947 0.8335
>
> Line (Intercept) 72.5945 8.5202
>
> Residual 180.0664 13.4189
>
>
> Fixed effects:
>
> Estimate Std. Error df t value Pr(>|t|)
>
> (Intercept) 501.141 2.268 96.940 221.009 < 2e-16 ***
>
> Temperature25 -57.960 2.699 54.800 -21.473 < 2e-16 ***
>
> SexM -53.639 1.001 96.260 -53.584 < 2e-16 ***
>
> Temperature25:SexM -6.488 1.391 48.300 -4.663 2.49e-05 ***
>
>
> I found that the function rand() from the lmerTest package gives me the p values of the random effects if I write the model like this:
>
>> rand(M2)
>
> Analysis of Random effects Table:
>
>
> Chi.sq Chi.DF p.value
>
> Line 4.6152 1 0.03 *
>
> Line:Sex:Temperature 30.8130 1 3e-08 ***
>
> Line:Sex 0.0391 1 0.84
>
> Line:Temperature 112.1539 1 <2e-16 ***
>
>
> I don’t know if this function is reliable because it is not mentioned for testing the significance of the random effects in the page https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html
GLMM FAQ - GitHub Pages<https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html>
bbolker.github.io
Introduction. This is an informal FAQ list for the r-sig-mixed-models mailing list.. The most commonly used functions for mixed modeling in R are. linear mixed models: aov(), nlme::lme 1, lme4::lmer;




Having written rand (=ranova) my views may be biased but I should say
that it _is_ reliable. That is, I haven't seen cases where it's not
doing what was intended ;-) rand() and ranova() simply compares two
models using anova(m1, m2, refit=FALSE) where m2 differs from m1 by
having one of its random-effect terms reduced or removed.

>
> The summary of the same analysis made with the Statistica is:
>
>
>
>
> Effect
>
> SS
>
> Degr. of
>
> MS
>
> Den.Syn.
>
> Den.Syn.
>
> F
>
> p
>
> Intercept
>
> Fixed
>
> 764579151
>
> 1
>
> 764579151
>
> 48,001
>
> 12655,91
>
> 60412,83
>
> 0,000000
>
> Line
>
> Random
>
> 608181
>
> 48
>
> 12670
>
> 47,800
>
> 6489,64
>
> 1,95
>
> 0,011254
>
> Sex
>
> Fixed
>
> 3138661
>
> 1
>
> 3138661
>
> 48,038
>
> 495,62
>
> 6332,81
>
> 0,000000
>
> Temperature
>
> Fixed
>
> 3686660
>
> 1
>
> 3686660
>
> 48,003
>
> 6459,62
>
> 570,72
>
> 0,000000
>
> Line*Sex
>
> Random
>
> 23808
>
> 48
>
> 496
>
> 48,000
>
> 473,30
>
> 1,05
>
> 0,435866
>
> Line*Temperature
>
> Random
>
> 310413
>
> 48
>
> 6467
>
> 48,000
>
> 473,30
>
> 13,66
>
> 0,000000
>
> Sex*Temperature
>
> Fixed
>
> 10075
>
> 1
>
> 10075
>
> 48,040
>
> 472,94
>
> 21,30
>
> 0,000029
>
> Line*Sex*Temperature
>
> Random
>
> 22718
>
> 48
>
> 473
>
> 3696,000
>
> 167,33
>
> 2,83
>
> 0,000000
>
> Error
>
> 618467
>
> 3696
>
> 167
>
>
>
>
> But if I write the model with the other syntax:
>
> M1 <- lmer(Wing ~ Temperature * Sex + (Temperature * Sex | Line))
>
Writing the model as
M1 <- lmer(Wing ~ Temperature * Sex + (0 + Temperature:Sex | Line))
is often preferable as it makes the random-effect variance-covariance
matrix easier to interpret.

>
>
> the summary is the following:
>
>
>
> REML criterion at convergence: 31440.9
>
>
>
> Random effects:
>
> Groups Name Variance Std.Dev. Corr
>
> Line (Intercept) 266.78 16.333
>
> Temperature25 398.27 19.957 -0.60
>
> SexM 41.54 6.446 -0.56 0.46
>
> Temperature25:SexM 61.34 7.832 0.56 -0.61 -0.80
>
> Residual 167.33 12.936
>
>
>
> Fixed effects:
>
> Estimate Std. Error df t value Pr(>|t|)
>
> (Intercept) 501.603 2.371 48.046 211.586 < 2e-16 ***
>
> Temperature25 -58.423 2.911 48.027 -20.070 < 2e-16 ***
>
> SexM -53.659 1.095 47.964 -49.023 < 2e-16 ***
>
> Temperature 25:SexM -6.470 1.393 48.278 -4.644 2.66e-05 ***
>
> ---
>
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
> In addition, if I apply the ‘rand functionŽ for this syntax (M1), it doesn’t retourn the whole p-values of the random effects (Do not give me p value for line and line*Temperature*Sex)

Use lmerTest::ranova(M1, reduce.terms=FALSE) to achieve a test for the
entire random-effect term (and make sure you have lmerTest version >=
3.0-0 installed).

>
> Analysis of Random effects Table:
>
> Chi.sq Chi.DF p.value
>
> Temperatura:Línea 0.00e+00 0 1
>
> Sexo:Línea 1.46e-10 0 <2e-16 ***
>
>
> I really appreciate your time and dedication for answering this questions. Thank you for trying to help us understand a little more about the syntax of these complex models and thus better understand their correct approach.

You are not alone if you think this is complicated. My students are
for the most part challenged in simply writing up the mathematical
formula that correctly represents models such as Model1 -
transitioning from scalar random-effect terms (e.g. Model3b) to
vector-valued random-effect terms (e.g. Model1) often takes time and
dedication.

Cheers
Rune

>
> Thank you very much for your time everyone.
>
>
>
> Greetings,
>
> Nicolas
>
>
> ----------------------------------------------------------------------------------------
> Message: 2
> Date: Wed, 25 Apr 2018 16:11:38 -0400
> From: Dan <***@gmail.com>
> To: "R-SIG-Mixed-***@R-project.org"
> <R-sig-mixed-***@r-project.org>, Ben Bolker <***@gmail.com>
> Subject: [R-sig-ME] Specifying Multiple Random Effects in NLME
> Message-ID:
> <CAET4i1f-SH5zA6xcCxNXc091HCk+snMv+***@mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
>
> Hi all:
>
> I am curating an off-list thread about specifying multiple random effects
> in NLME.
>
> 1. If it's (1|Object) + (1|Object:Coating) that you want then
> you should be able to use a nested specification (which nlme *can*
> handle relatively easily), i.e. something like
>
> random=a+b+c~1|Object/Coating
>
>
> Although (Coating|Object) and (1|Object:Coating) both in some sense
> represent "interactions" the latter is *much* simpler/more parsimonious.
>
> If you're OK with 1|Object:Coating rather than Coating|Object it
> should be *much* faster. If you don't understand the distinction (which
> would be absolutely fine and understandable) can we resume the
> discussion on r-sig-mixed-models ... ?
> -----------
>
> So:
> random=a+b+c~Coating|Object
> does not fit.
>
> But:
> random=a+b+c~Object/Coating
> fits.
>
> Can you better explain the distinction here? I have sometimes used the
> 1|Object:Coating + 1|Object syntax and sometimes the Coating|Object syntax
> in other models. My experience/understanding is that the former syntax with
> multiple "within subject" variables produces exactly matching output to the
> standard "repeated measures ANOVA" with the lmer assumption of compound
> symmetry.
>
> Are there cases where one vs. the other formulation should absolutely be
> used? My understanding that for continuous variables, e.g., multiple
> measurements across multiple days, Days|Object would be the correct syntax.
> But here we're talking about a factor variable.
>
>
> 2. I'm trying to read the "random" section for nlme right now but it's
> kind of making my head explode (and I think there's a typo: " the same
> as the order of the order of the elements in the list"). It *sounds*
> like (1) explicitly creating an interaction
> ObjCoating=interaction(Object,Coating) and (2) using something like
>
> list(ObjCoating=a~1,Object=b~1,Object=c~1)
>
> should work (grouping factors as names, then [right-hand-side variable
> name]~[random effects model], but I'm worried about the phrase "The
> order of nesting will be assumed the same as the order of the elements
> in the list": what nesting?
> -----------
>
> I think that formulation is explicitly in order. I replaced your first
> ObjCoating with simply Object, just to test what would happen:
>
> Random effects:
> Formula: a ~ 1 | Object
> a.(Intercept)
> StdDev: 1.305816
>
> Formula: b ~ 1 | Object %in% Object
> b.(Intercept)
> StdDev: 0.01576521
>
> Formula: c ~ 1 | Object %in% Object %in% Object
> c.(Intercept) Residual
> StdDev: 2.677883 2.219676
>
> [[alternative HTML version deleted]]
>
>
>
>
> *****
>
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-mixed-***@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

[[alternative HTML version deleted]]