[R-sig-ME] Crossing Interaction and lmer

[Douglas Bates]

On Wed, Oct 21, 2009 at 4:50 AM, Johannes Schliesser <J.Schliesser at uu.nl> wrote:
> Dear mixed model experts
>
> I have troubles with a wired outcome with lmer().
>
> I was testing effects on reaction times for two fixed factors, Identity
> (levels identical "id" and non-identical ("ni") and Bonus Size ("h"igh and
> "l"ow). Subjects "subj" and "items" were Random Effects. In fact, Bonus Size
> is nested within Items, as half of the Items were of high, the other half of
> low Bonus size.
>
> The results indicate a strong crossing interaction:
>
>
>       |    h   |    l   |   ALL  |
> |id   |839.      |861  |  850   |
>
> |ni    |879    |829 |  854   |
>
> |ALL| 859    |845  |852    |
>
> Subjects tend to be about 40 ms faster for identical than for non-identical in
> the high bonus group, but about 30 ms slower for identical than for non-
> identical in the low bonus group. Grand means for identical and non-identical
> items differ only for 4 ms.
>
> So far so good.
> Classic ANOVA rejects all main effects but confirmes the interaction.
>
> The lmer analysis instead suggests a fixed effect for identity:
> ----------
> lmer.rtNA2 = lmer(rtNA~id*bonus+(1|subj)+(1|bonus/item), data = x4)

The model that you summarize below is not the model fit by this
expression and, in fact, this model would not make sense.

>  lmer.rtNA
> Linear mixed model fit by REML
> Formula: rtNA ~ id * bonus + (1 | subj) + (1 | item)
>   Data: x4
>   AIC   BIC logLik deviance REMLdev
>  24346 24385 -12166    24360   24332
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  item     (Intercept)  1926.9   43.896
>  subj     (Intercept)  9648.6   98.227
>  Residual             19026.7  137.937
> Number of obs: 1905, groups: item, 80; subj, 24
>
> Fixed effects:
>            Estimate Std. Error t value
> (Intercept)  839.804     22.132   37.95
> idni          40.264      8.942    4.50
> bonusl        22.841     13.278    1.72
> idni:bonusl  -73.346     12.644   -5.80
>
> Correlation of Fixed Effects:
>            (Intr) idni   bonusl
> idni        -0.202
> bonusl      -0.300  0.336
> idni:bonusl  0.143 -0.707 -0.476
> -----------------------------------
>
> The lmer model seems to test Identical vs. Non-Identical for the High-Bonus
> group only, but not overall (same as estimates for Bonus Size only in the "id"
> row). pvals.fnc() asssigns significance for Identity and the interaction
> Identity*Bonus Size.

You are paying too much attention to the test statistics.  The general
rule is that you check the interaction first and, if the interaction
is significant, then you do not proceed to check the main effects.  As
you have seen, the coefficient for id relates to the behavior at the
reference level of the bonus.  If there is no interaction then this
will also relate to the behavior at the other level of bonus.  When
there is an interaction the behavior with respect to id depends on the
level of bonus.  That's what an interaction means.

> I have the impression that this is missleading: we certainly cannot speak of
> slower reaction times for non-identical stimuli than for identical.

You are trying to make statements about the overall behavior of one
factor in the presence of a non-ignorable interaction with another.
Contrary to popular belief this is not possible.  Look at your
original table or, better yet, create an interaction plot.  These
coefficient estimates are just another way of expressing that table.
At the high level of bonus non-identical are slower than identical.
At the low level of bonus it is the other way around.

> Leaving out the Bonus Size as a fixed factor "lmer(rt~id+(1|subj)+(1|item),
> data = x4)" gives a realistic t-value below 1 for Identity.
>
> So, I ask:
> Is the missleading handling of crossing interactions a bug ?

No.  You have fallen into the trap of associating the result of a test
with a term when, in fact, the test is on the value of the coefficient
and the interpretation of the tests depends on how the coefficients
are defined.

In R there are many different ways of defining such coefficients
related to how the "contrasts" for a factor are assigned.  For a
two-level factorial model like this I often define the contrasts to be
contr.sum so that I get orthogonal columns in the model matrix
(assuming that the original factors are balanced).  For example

> contr.sum(2)
  [,1]
1    1
2   -1
> contr.poly(2)
             .L
[1,] -0.7071068
[2,]  0.7071068
> dat <- data.frame(iden = gl(2,2,labels = c("id","ni")),
+   bonus = gl(2,1,4, labels = c('h','l')))
> dat
  iden bonus
1   id     h
2   id     l
3   ni     h
4   ni     l
> model.matrix(~ iden * bonus, dat)
  (Intercept) idenni bonusl idenni:bonusl
1           1      0      0             0
2           1      0      1             0
3           1      1      0             0
4           1      1      1             1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$iden
[1] "contr.treatment"

attr(,"contrasts")$bonus
[1] "contr.treatment"

> options(contrasts = c("contr.sum", "contr.poly"))
> model.matrix(~ iden * bonus, dat)
  (Intercept) iden1 bonus1 iden1:bonus1
1           1     1      1            1
2           1     1     -1           -1
3           1    -1      1           -1
4           1    -1     -1            1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$iden
[1] "contr.sum"

attr(,"contrasts")$bonus
[1] "contr.sum"

> Or do I simply have to know that if there is a huge crossing interaction I
> should not generalize the Fixed Factor effects ?

Exactly.

> And more relevant: How do I get an analysis that is appropriate to the data ?

The first analysis is appropriate.  For me the main question to
address in interpretation would be why does the effect of
identical/non-identical at low bonus look so different from the effect
at high bonus.