[R] R: lmer specification for random effects: contradictory reults

Mon Nov 25 15:45:30 CET 2013

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.