[R-sig-ME] [R] lmer specification for random effects: contradictory reults
ONKELINX, Thierry
Thierry.ONKELINX at inbo.be
Tue Nov 26 10:29:26 CET 2013
IMHO, if your design requires a random effect, then you add it to the model. Significance is irrelevant in that case.
Comparing mod2.1 and mod2.4 is not testing for a random slope but for the fixed effect of Z. Instead you should compare a model with (1 + T + Z|subject) versus a model with (1 + T | subject). Note that a) you must use REML for that comparison and b) you are testing on the boundary.
Zuur et al (2009) has a nice introduction on model building for mixed models.
@BOOK{
title = {{M}ixed {E}ffects {M}odels and {E}xtensions in {E}cology with {R}},
publisher = {Springer New York},
year = {2009},
author = {Zuur, Alain F. and Ieno, Elena N. and Walker, Neil J. and Saveliev,
Anatoly A. and Smith, Graham M.},
doi = {10.1007/978-0-387-87458-6}
}
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 op 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: Benedetta Cesqui [mailto:b.cesqui op hsantalucia.it]
Verzonden: maandag 25 november 2013 15:46
Aan: ONKELINX, Thierry; r-help op r-project.org
CC: r-sig-mixed-models op r-project.org
Onderwerp: 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 op inbo.be]
Inviato: lunedì 25 novembre 2013 14:48
A: Benedetta Cesqui; r-help op r-project.org
Cc: r-sig-mixed-models op 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 op 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 op r-project.org [mailto:r-help-bounces op r-project.org]
Namens Benedetta Cesqui
Verzonden: maandag 25 november 2013 11:13
Aan: r-help op 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 op hsantalucia.it
[[alternative HTML version deleted]]
______________________________________________
R-help op 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.
* * * * * * * * * * * * * 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.
More information about the R-sig-mixed-models
mailing list