Benedetta Cesqui
2013-Nov-25 10:13 UTC
[R] lmer specification for random effects: contradictory reults
Hi All,
I was wondering if someone could help me to solve this issue with lmer.
In order to understand the best mixed effects model to fit my data, I
compared the following options according to the procedures specified in many
papers (i.e. Baayen
<http://www.google.it/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CDsQFjAA
&url=http%3A%2F%2Fwww.ualberta.ca%2F~baayen%2Fpublications%2FbaayenDavidsonB
ates.pdf&ei=FhqTUoXuJKKV7Abds4GYBA&usg=AFQjCNFst7GT7mBX7w9lXItJTtELJSKWJg&si
g2=KGA5MHxOvEGwDxf-Gcqi6g&bvm> R.H. et al 2008)
Here, dT_purs is the response variable, T and Z are the fixed effects, and
subject is the random effect. Random and fixed effects are crossed.:
mod0 <- lmer(dT_purs ~ T + Z + (1|subject), data = x)
mod1 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject), data = x)
mod2 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod3 <- lmer(dT_purs ~ T * Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod4 <- lmer(dT_purs ~ T * Z + (1| subject), data = x)
anova(mod0, mod1,mod2, mod3, mod4)
Data: x
Models:
mod0: dT_purs ~ T + Z + (1 | subject)
mod4: dT_purs ~ T * Z + (1 | subject )
mod1: dT_purs ~ T + Z + (1 + T| subject)
mod2: dT_purs ~ T + Z + (1 + T| subject ) + (1 + Z | subject)
mod3: dT_purs ~ T * Z + (1 + T| subject) + (1 + Z | subject)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
mod0 5 -689.81 -669.46 349.91 -699.81
mod4 6 -689.57 -665.14 350.78 -701.57 1.7532 1 0.185473
mod1 7 -689.12 -660.62 351.56 -703.12 1.5504 1 0.213070
mod2 10 -695.67 -654.97 357.84 -715.67 12.5563 3 0.005701 **
mod3 11 -695.83 -651.05 358.92 -717.83 2.1580 1 0.141825
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
It turns out that mod2 has the right level of complexity for this dataset.
However when I looked at its summary, I got a correlation of -0.87 for the
random effects relative to the T effect and -1 for the random effects
relatively to the Z.
summary(mod2)
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: dT_purs ~T + Z + (1 + T | subject) + (1 + Z | subject)
Data: x
AIC BIC logLik deviance
-695.6729 -654.9655 357.8364 -715.6729
Random effects:
Groups Name Variance Std.Dev. Corr
subject (Intercept) 0.0032063 0.05662
T 0.0117204 0.10826 -0.87
subject.1 (Intercept) 0.0005673 0.02382
Z 0.0025859 0.05085 1.00
Residual 0.0104551 0.10225
Number of obs: 433, groups: soggetto, 7
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.02489 0.03833 0.650
T 0.52010 0.05905 8.808
Z -0.09019 0.02199 -4.101
Correlation of Fixed Effects:
(Intr) tempo
T -0.901
Z 0.218 -0.026
If I understand correctly what the correlation parameters reported in the
table are, the correlation of 1 means that, for the Z effects the random
intercept is perfectly collinear with the random slope. Thus, we fit the
wrong model. A random intercept only model would have sufficed.
Am I correct?
If so, should I take mod1 (mod1 <- dT_purs ~ T + Z + (1 + T | subject )
instead of mod2 to fit my data?
Why are these results contradictory?
Finally is a correlation value of -0.87 a too high or an acceptable value ?
Thanks for help me in advance!
Best
Benedetta
---
Benedetta Cesqui, Ph.D.
Laboratory of Neuromotor Physiology
IRCCS Fondazione Santa Lucia
Via Ardeatina 306
00179 Rome, Italy
tel: (+39) 06-51501485
fax:(+39) 06-51501482
E_mail: b.cesqui@hsantalucia.it
[[alternative HTML version deleted]]
ONKELINX, Thierry
2013-Nov-25 13:48 UTC
[R] lmer specification for random effects: contradictory reults
Dear Benedetta,
I think you might want (1+T+Z|subject) as random effects rather than
(1+T|subject) + (1 + Z|subject). The latter has two random intercepts per
subject: a recipe for disaster.
Follow-up posts should only go to the mixed models mailing list which I'm
cc'ing.
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx at inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more than
asking him to perform a post-mortem examination: he may be able to say what the
experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure
that a reasonable answer can be extracted from a given body of data.
~ John Tukey
-----Oorspronkelijk bericht-----
Van: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
Namens Benedetta Cesqui
Verzonden: maandag 25 november 2013 11:13
Aan: r-help at r-project.org
Onderwerp: [R] lmer specification for random effects: contradictory reults
Hi All,
I was wondering if someone could help me to solve this issue with lmer.
In order to understand the best mixed effects model to fit my data, I compared
the following options according to the procedures specified in many papers (i.e.
Baayen
<http://www.google.it/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CDsQFjAA
&url=http%3A%2F%2Fwww.ualberta.ca%2F~baayen%2Fpublications%2FbaayenDavidsonB
ates.pdf&ei=FhqTUoXuJKKV7Abds4GYBA&usg=AFQjCNFst7GT7mBX7w9lXItJTtELJSKWJg&si
g2=KGA5MHxOvEGwDxf-Gcqi6g&bvm> R.H. et al 2008) Here, dT_purs is the
response variable, T and Z are the fixed effects, and subject is the random
effect. Random and fixed effects are crossed.:
mod0 <- lmer(dT_purs ~ T + Z + (1|subject), data = x)
mod1 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject), data = x)
mod2 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod3 <- lmer(dT_purs ~ T * Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod4 <- lmer(dT_purs ~ T * Z + (1| subject), data = x)
anova(mod0, mod1,mod2, mod3, mod4)
Data: x
Models:
mod0: dT_purs ~ T + Z + (1 | subject)
mod4: dT_purs ~ T * Z + (1 | subject )
mod1: dT_purs ~ T + Z + (1 + T| subject)
mod2: dT_purs ~ T + Z + (1 + T| subject ) + (1 + Z | subject)
mod3: dT_purs ~ T * Z + (1 + T| subject) + (1 + Z | subject)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
mod0 5 -689.81 -669.46 349.91 -699.81
mod4 6 -689.57 -665.14 350.78 -701.57 1.7532 1 0.185473
mod1 7 -689.12 -660.62 351.56 -703.12 1.5504 1 0.213070
mod2 10 -695.67 -654.97 357.84 -715.67 12.5563 3 0.005701 **
mod3 11 -695.83 -651.05 358.92 -717.83 2.1580 1 0.141825
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
It turns out that mod2 has the right level of complexity for this dataset.
However when I looked at its summary, I got a correlation of -0.87 for the
random effects relative to the T effect and -1 for the random effects
relatively to the Z.
summary(mod2)
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: dT_purs ~T + Z + (1 + T | subject) + (1 + Z | subject)
Data: x
AIC BIC logLik deviance
-695.6729 -654.9655 357.8364 -715.6729
Random effects:
Groups Name Variance Std.Dev. Corr
subject (Intercept) 0.0032063 0.05662
T 0.0117204 0.10826 -0.87
subject.1 (Intercept) 0.0005673 0.02382
Z 0.0025859 0.05085 1.00
Residual 0.0104551 0.10225
Number of obs: 433, groups: soggetto, 7
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.02489 0.03833 0.650
T 0.52010 0.05905 8.808
Z -0.09019 0.02199 -4.101
Correlation of Fixed Effects:
(Intr) tempo
T -0.901
Z 0.218 -0.026
If I understand correctly what the correlation parameters reported in the
table are, the correlation of 1 means that, for the Z effects the random
intercept is perfectly collinear with the random slope. Thus, we fit the
wrong model. A random intercept only model would have sufficed.
Am I correct?
If so, should I take mod1 (mod1 <- dT_purs ~ T + Z + (1 + T | subject )
instead of mod2 to fit my data?
Why are these results contradictory?
Finally is a correlation value of -0.87 a too high or an acceptable value ?
Thanks for help me in advance!
Best
Benedetta
---
Benedetta Cesqui, Ph.D.
Laboratory of Neuromotor Physiology
IRCCS Fondazione Santa Lucia
Via Ardeatina 306
00179 Rome, Italy
tel: (+39) 06-51501485
fax:(+39) 06-51501482
E_mail: b.cesqui at hsantalucia.it
[[alternative HTML version deleted]]
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en
binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer
and may not be regarded as stating an official position of INBO, as long as the
message is not confirmed by a duly signed document.
ONKELINX, Thierry
2013-Nov-25 13:48 UTC
Re: lmer specification for random effects: contradictory reults
Dear Benedetta,
I think you might want (1+T+Z|subject) as random effects rather than
(1+T|subject) + (1 + Z|subject). The latter has two random intercepts per
subject: a recipe for disaster.
Follow-up posts should only go to the mixed models mailing list which
I''m cc''ing.
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx@inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more than
asking him to perform a post-mortem examination: he may be able to say what the
experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure
that a reasonable answer can be extracted from a given body of data.
~ John Tukey
-----Oorspronkelijk bericht-----
Van: r-help-bounces@r-project.org [mailto:r-help-bounces@r-project.org] Namens
Benedetta Cesqui
Verzonden: maandag 25 november 2013 11:13
Aan: r-help@r-project.org
Onderwerp: [R] lmer specification for random effects: contradictory reults
Hi All,
I was wondering if someone could help me to solve this issue with lmer.
In order to understand the best mixed effects model to fit my data, I compared
the following options according to the procedures specified in many papers (i.e.
Baayen
<http://www.google.it/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CDsQFjAA
&url=http%3A%2F%2Fwww.ualberta.ca%2F~baayen%2Fpublications%2FbaayenDavidsonB
ates.pdf&ei=FhqTUoXuJKKV7Abds4GYBA&usg=AFQjCNFst7GT7mBX7w9lXItJTtELJSKWJg&si
g2=KGA5MHxOvEGwDxf-Gcqi6g&bvm> R.H. et al 2008) Here, dT_purs is the
response variable, T and Z are the fixed effects, and subject is the random
effect. Random and fixed effects are crossed.:
mod0 <- lmer(dT_purs ~ T + Z + (1|subject), data = x)
mod1 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject), data = x)
mod2 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod3 <- lmer(dT_purs ~ T * Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod4 <- lmer(dT_purs ~ T * Z + (1| subject), data = x)
anova(mod0, mod1,mod2, mod3, mod4)
Data: x
Models:
mod0: dT_purs ~ T + Z + (1 | subject)
mod4: dT_purs ~ T * Z + (1 | subject )
mod1: dT_purs ~ T + Z + (1 + T| subject)
mod2: dT_purs ~ T + Z + (1 + T| subject ) + (1 + Z | subject)
mod3: dT_purs ~ T * Z + (1 + T| subject) + (1 + Z | subject)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
mod0 5 -689.81 -669.46 349.91 -699.81
mod4 6 -689.57 -665.14 350.78 -701.57 1.7532 1 0.185473
mod1 7 -689.12 -660.62 351.56 -703.12 1.5504 1 0.213070
mod2 10 -695.67 -654.97 357.84 -715.67 12.5563 3 0.005701 **
mod3 11 -695.83 -651.05 358.92 -717.83 2.1580 1 0.141825
---
Signif. codes: 0 ''***'' 0.001 ''**'' 0.01
''*'' 0.05 ''.'' 0.1 '' '' 1
It turns out that mod2 has the right level of complexity for this dataset.
However when I looked at its summary, I got a correlation of -0.87 for the
random effects relative to the T effect and -1 for the random effects
relatively to the Z.
summary(mod2)
Linear mixed model fit by maximum likelihood [''lmerMod'']
Formula: dT_purs ~T + Z + (1 + T | subject) + (1 + Z | subject)
Data: x
AIC BIC logLik deviance
-695.6729 -654.9655 357.8364 -715.6729
Random effects:
Groups Name Variance Std.Dev. Corr
subject (Intercept) 0.0032063 0.05662
T 0.0117204 0.10826 -0.87
subject.1 (Intercept) 0.0005673 0.02382
Z 0.0025859 0.05085 1.00
Residual 0.0104551 0.10225
Number of obs: 433, groups: soggetto, 7
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.02489 0.03833 0.650
T 0.52010 0.05905 8.808
Z -0.09019 0.02199 -4.101
Correlation of Fixed Effects:
(Intr) tempo
T -0.901
Z 0.218 -0.026
If I understand correctly what the correlation parameters reported in the
table are, the correlation of 1 means that, for the Z effects the random
intercept is perfectly collinear with the random slope. Thus, we fit the
wrong model. A random intercept only model would have sufficed.
Am I correct?
If so, should I take mod1 (mod1 <- dT_purs ~ T + Z + (1 + T | subject )
instead of mod2 to fit my data?
Why are these results contradictory?
Finally is a correlation value of -0.87 a too high or an acceptable value ?
Thanks for help me in advance!
Best
Benedetta
---
Benedetta Cesqui, Ph.D.
Laboratory of Neuromotor Physiology
IRCCS Fondazione Santa Lucia
Via Ardeatina 306
00179 Rome, Italy
tel: (+39) 06-51501485
fax:(+39) 06-51501482
E_mail: b.cesqui@hsantalucia.it
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en
binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer
and may not be regarded as stating an official position of INBO, as long as the
message is not confirmed by a duly signed document.
Benedetta Cesqui
2013-Nov-25 14:45 UTC
[R] R: lmer specification for random effects: contradictory reults
Dear Thierry,
thank you for the quick reply.
I have only one question about the approach you proposed.
As you suggested, imagine that the model we end up after the model selection
procedure is:
mod2.1 <- lmer(dT_purs ~ T + Z + (1 +T+Z| subject), data =x, REML= FALSE)
According to the common procedures specified in many manuals and recent
papers, if I want to compute the p_values relative to each term, I will
perform a likelihood test, in which the deviance of the (-2LL) of a model
containing the specific term is compared to another model without it.
In the case of the fixed effect terms I have no problem in the
interpretation of the results. Each comparison returns a significance
associated with the estimated coefficient of the term.
Thus in this case:
mod2.2 <- lmer(dT_purs ~ Z + (1 +T+Z|soggetto) , data = x, REML = FALSE)
mod2.3 <- lmer(dT_purs ~ T + (1 +T+Z|soggetto) , data = x, REML = FALSE)
anova(mod2.1, mod2.2)
p_valueT = 3.203e-05
anova(mod2.1, mod2.3)
p_valueZ = 0.001793
What about the p_value relative to the (1+T+Z|subject)?
One option is to compute:
mod2.4 <- lm(dT_purs ~ T + (1 +T+Z|soggetto) , data = x)
and then execute the loklikelihood test as follows:
L0 <-logLik(mod2.4)
L1 <-logLik(mod2.1)
LR <--2*(L1-L0)
pv <- pchisq(LR,2,ncp = 0, lower.tai=FALSE,log.p = FALSE)
However, what can I conclude on the random slopeif it is significant?
With the previouse approach using the model:
mod2 <- lmer(dT_purs ~ T + Z + (1 +T| subject) + (1+ Z| subject), data =x)
The comparison among the models in which the different termd were
included/excluded provided me the following results:
p_valueT = 1.269e-07;
p_valueZ =0.00322
p_valueTS = 0.4277
p_valueZS = 0.005701
I interpreted the ones relative to the random effects as if the subjects
differed not only in their overall responses, but also in the nature of
their response dT_purse values in the different T conditions, but not in the
different Z conditions.
Benedetta
-----Messaggio originale-----
Da: ONKELINX, Thierry [mailto:Thierry.ONKELINX at inbo.be]
Inviato: luned? 25 novembre 2013 14:48
A: Benedetta Cesqui; r-help at r-project.org
Cc: r-sig-mixed-models at r-project.org
Oggetto: RE: [R] lmer specification for random effects: contradictory reults
Dear Benedetta,
I think you might want (1+T+Z|subject) as random effects rather than
(1+T|subject) + (1 + Z|subject). The latter has two random intercepts per
subject: a recipe for disaster.
Follow-up posts should only go to the mixed models mailing list which I'm
cc'ing.
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest team Biometrie & Kwaliteitszorg / team Biometrics & Quality
Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx at inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more than
asking him to perform a post-mortem examination: he may be able to say what
the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
-----Oorspronkelijk bericht-----
Van: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
Namens Benedetta Cesqui
Verzonden: maandag 25 november 2013 11:13
Aan: r-help at r-project.org
Onderwerp: [R] lmer specification for random effects: contradictory reults
Hi All,
I was wondering if someone could help me to solve this issue with lmer.
In order to understand the best mixed effects model to fit my data, I
compared the following options according to the procedures specified in many
papers (i.e. Baayen
<http://www.google.it/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CDsQFjAA
&url=http%3A%2F%2Fwww.ualberta.ca%2F~baayen%2Fpublications%2FbaayenDavidsonB
ates.pdf&ei=FhqTUoXuJKKV7Abds4GYBA&usg=AFQjCNFst7GT7mBX7w9lXItJTtELJSKWJg&si
g2=KGA5MHxOvEGwDxf-Gcqi6g&bvm> R.H. et al 2008) Here, dT_purs is the
response variable, T and Z are the fixed effects, and subject is the random
effect. Random and fixed effects are crossed.:
mod0 <- lmer(dT_purs ~ T + Z + (1|subject), data = x)
mod1 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject), data = x)
mod2 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod3 <- lmer(dT_purs ~ T * Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod4 <- lmer(dT_purs ~ T * Z + (1| subject), data = x)
anova(mod0, mod1,mod2, mod3, mod4)
Data: x
Models:
mod0: dT_purs ~ T + Z + (1 | subject)
mod4: dT_purs ~ T * Z + (1 | subject )
mod1: dT_purs ~ T + Z + (1 + T| subject)
mod2: dT_purs ~ T + Z + (1 + T| subject ) + (1 + Z | subject)
mod3: dT_purs ~ T * Z + (1 + T| subject) + (1 + Z | subject)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
mod0 5 -689.81 -669.46 349.91 -699.81
mod4 6 -689.57 -665.14 350.78 -701.57 1.7532 1 0.185473
mod1 7 -689.12 -660.62 351.56 -703.12 1.5504 1 0.213070
mod2 10 -695.67 -654.97 357.84 -715.67 12.5563 3 0.005701 **
mod3 11 -695.83 -651.05 358.92 -717.83 2.1580 1 0.141825
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
It turns out that mod2 has the right level of complexity for this dataset.
However when I looked at its summary, I got a correlation of -0.87 for the
random effects relative to the T effect and -1 for the random effects
relatively to the Z.
summary(mod2)
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: dT_purs ~T + Z + (1 + T | subject) + (1 + Z | subject)
Data: x
AIC BIC logLik deviance
-695.6729 -654.9655 357.8364 -715.6729
Random effects:
Groups Name Variance Std.Dev. Corr
subject (Intercept) 0.0032063 0.05662
T 0.0117204 0.10826 -0.87
subject.1 (Intercept) 0.0005673 0.02382
Z 0.0025859 0.05085 1.00
Residual 0.0104551 0.10225
Number of obs: 433, groups: soggetto, 7
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.02489 0.03833 0.650
T 0.52010 0.05905 8.808
Z -0.09019 0.02199 -4.101
Correlation of Fixed Effects:
(Intr) tempo
T -0.901
Z 0.218 -0.026
If I understand correctly what the correlation parameters reported in the
table are, the correlation of 1 means that, for the Z effects the random
intercept is perfectly collinear with the random slope. Thus, we fit the
wrong model. A random intercept only model would have sufficed.
Am I correct?
If so, should I take mod1 (mod1 <- dT_purs ~ T + Z + (1 + T | subject )
instead of mod2 to fit my data?
Why are these results contradictory?
Finally is a correlation value of -0.87 a too high or an acceptable value ?
Thanks for help me in advance!
Best
Benedetta
---
Benedetta Cesqui, Ph.D.
Laboratory of Neuromotor Physiology
IRCCS Fondazione Santa Lucia
Via Ardeatina 306
00179 Rome, Italy
tel: (+39) 06-51501485
fax:(+39) 06-51501482
E_mail: b.cesqui at hsantalucia.it
[[alternative HTML version deleted]]
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * * Dit
bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en
binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd
is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the
writer and may not be regarded as stating an official position of INBO, as
long as the message is not confirmed by a duly signed document.
Benedetta Cesqui
2013-Nov-25 14:45 UTC
R: lmer specification for random effects: contradictory reults
Dear Thierry,
thank you for the quick reply.
I have only one question about the approach you proposed.
As you suggested, imagine that the model we end up after the model selection
procedure is:
mod2.1 <- lmer(dT_purs ~ T + Z + (1 +T+Z| subject), data =x, REML= FALSE)
According to the common procedures specified in many manuals and recent
papers, if I want to compute the p_values relative to each term, I will
perform a likelihood test, in which the deviance of the (-2LL) of a model
containing the specific term is compared to another model without it.
In the case of the fixed effect terms I have no problem in the
interpretation of the results. Each comparison returns a significance
associated with the estimated coefficient of the term.
Thus in this case:
mod2.2 <- lmer(dT_purs ~ Z + (1 +T+Z|soggetto) , data = x, REML = FALSE)
mod2.3 <- lmer(dT_purs ~ T + (1 +T+Z|soggetto) , data = x, REML = FALSE)
anova(mod2.1, mod2.2)
p_valueT = 3.203e-05
anova(mod2.1, mod2.3)
p_valueZ = 0.001793
What about the p_value relative to the (1+T+Z|subject)?
One option is to compute:
mod2.4 <- lm(dT_purs ~ T + (1 +T+Z|soggetto) , data = x)
and then execute the loklikelihood test as follows:
L0 <-logLik(mod2.4)
L1 <-logLik(mod2.1)
LR <--2*(L1-L0)
pv <- pchisq(LR,2,ncp = 0, lower.tai=FALSE,log.p = FALSE)
However, what can I conclude on the random slopeif it is significant?
With the previouse approach using the model:
mod2 <- lmer(dT_purs ~ T + Z + (1 +T| subject) + (1+ Z| subject), data =x)
The comparison among the models in which the different termd were
included/excluded provided me the following results:
p_valueT = 1.269e-07;
p_valueZ =0.00322
p_valueTS = 0.4277
p_valueZS = 0.005701
I interpreted the ones relative to the random effects as if the subjects
differed not only in their overall responses, but also in the nature of
their response dT_purse values in the different T conditions, but not in the
different Z conditions.
Benedetta
-----Messaggio originale-----
Da: ONKELINX, Thierry [mailto:Thierry.ONKELINX@inbo.be]
Inviato: lunedì 25 novembre 2013 14:48
A: Benedetta Cesqui; r-help@r-project.org
Cc: r-sig-mixed-models@r-project.org
Oggetto: RE: [R] lmer specification for random effects: contradictory reults
Dear Benedetta,
I think you might want (1+T+Z|subject) as random effects rather than
(1+T|subject) + (1 + Z|subject). The latter has two random intercepts per
subject: a recipe for disaster.
Follow-up posts should only go to the mixed models mailing list which
I''m
cc''ing.
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest team Biometrie & Kwaliteitszorg / team Biometrics & Quality
Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx@inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more than
asking him to perform a post-mortem examination: he may be able to say what
the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
-----Oorspronkelijk bericht-----
Van: r-help-bounces@r-project.org [mailto:r-help-bounces@r-project.org]
Namens Benedetta Cesqui
Verzonden: maandag 25 november 2013 11:13
Aan: r-help@r-project.org
Onderwerp: [R] lmer specification for random effects: contradictory reults
Hi All,
I was wondering if someone could help me to solve this issue with lmer.
In order to understand the best mixed effects model to fit my data, I
compared the following options according to the procedures specified in many
papers (i.e. Baayen
<http://www.google.it/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CDsQFjAA
&url=http%3A%2F%2Fwww.ualberta.ca%2F~baayen%2Fpublications%2FbaayenDavidsonB
ates.pdf&ei=FhqTUoXuJKKV7Abds4GYBA&usg=AFQjCNFst7GT7mBX7w9lXItJTtELJSKWJg&si
g2=KGA5MHxOvEGwDxf-Gcqi6g&bvm> R.H. et al 2008) Here, dT_purs is the
response variable, T and Z are the fixed effects, and subject is the random
effect. Random and fixed effects are crossed.:
mod0 <- lmer(dT_purs ~ T + Z + (1|subject), data = x)
mod1 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject), data = x)
mod2 <- lmer(dT_purs ~ T + Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod3 <- lmer(dT_purs ~ T * Z + (1 +tempo| subject) + (1+ Z| subject), data x)
mod4 <- lmer(dT_purs ~ T * Z + (1| subject), data = x)
anova(mod0, mod1,mod2, mod3, mod4)
Data: x
Models:
mod0: dT_purs ~ T + Z + (1 | subject)
mod4: dT_purs ~ T * Z + (1 | subject )
mod1: dT_purs ~ T + Z + (1 + T| subject)
mod2: dT_purs ~ T + Z + (1 + T| subject ) + (1 + Z | subject)
mod3: dT_purs ~ T * Z + (1 + T| subject) + (1 + Z | subject)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
mod0 5 -689.81 -669.46 349.91 -699.81
mod4 6 -689.57 -665.14 350.78 -701.57 1.7532 1 0.185473
mod1 7 -689.12 -660.62 351.56 -703.12 1.5504 1 0.213070
mod2 10 -695.67 -654.97 357.84 -715.67 12.5563 3 0.005701 **
mod3 11 -695.83 -651.05 358.92 -717.83 2.1580 1 0.141825
---
Signif. codes: 0 ''***'' 0.001 ''**'' 0.01
''*'' 0.05 ''.'' 0.1 '' '' 1
It turns out that mod2 has the right level of complexity for this dataset.
However when I looked at its summary, I got a correlation of -0.87 for the
random effects relative to the T effect and -1 for the random effects
relatively to the Z.
summary(mod2)
Linear mixed model fit by maximum likelihood [''lmerMod'']
Formula: dT_purs ~T + Z + (1 + T | subject) + (1 + Z | subject)
Data: x
AIC BIC logLik deviance
-695.6729 -654.9655 357.8364 -715.6729
Random effects:
Groups Name Variance Std.Dev. Corr
subject (Intercept) 0.0032063 0.05662
T 0.0117204 0.10826 -0.87
subject.1 (Intercept) 0.0005673 0.02382
Z 0.0025859 0.05085 1.00
Residual 0.0104551 0.10225
Number of obs: 433, groups: soggetto, 7
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.02489 0.03833 0.650
T 0.52010 0.05905 8.808
Z -0.09019 0.02199 -4.101
Correlation of Fixed Effects:
(Intr) tempo
T -0.901
Z 0.218 -0.026
If I understand correctly what the correlation parameters reported in the
table are, the correlation of 1 means that, for the Z effects the random
intercept is perfectly collinear with the random slope. Thus, we fit the
wrong model. A random intercept only model would have sufficed.
Am I correct?
If so, should I take mod1 (mod1 <- dT_purs ~ T + Z + (1 + T | subject )
instead of mod2 to fit my data?
Why are these results contradictory?
Finally is a correlation value of -0.87 a too high or an acceptable value ?
Thanks for help me in advance!
Best
Benedetta
---
Benedetta Cesqui, Ph.D.
Laboratory of Neuromotor Physiology
IRCCS Fondazione Santa Lucia
Via Ardeatina 306
00179 Rome, Italy
tel: (+39) 06-51501485
fax:(+39) 06-51501482
E_mail: b.cesqui@hsantalucia.it
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * * Dit
bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en
binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd
is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the
writer and may not be regarded as stating an official position of INBO, as
long as the message is not confirmed by a duly signed document.