[R-sig-ME] test for significance of interaction in linear mixed models in nlme package in R

[Ben Bolker]

Gu Mi <neo.migu at ...> writes:

> 
> Dear All:
> 
> I use lme function in the nlme R package to test if levels of factor items
> has significant interaction with levels of factor condition. The factor
> condition has two levels: Control and Treatment, and the factor items has 3
> levels: E1,E2,E3. I use the following code:
> 
> f.lme.1 = lme(response ~ 0 + factor(condition) * factor(items), 
>   random = ~1|subject)

  Minor points:
  * is generally better practice to pull your input variables from
a data frame rather than leaving them lying around the workspace
(or, worse, using attach())
  * your results will be cleaner if you convert your factors
beforehand

  d <- data.frame(response,condition=factor(condition),
                   items=factor(items),subject)
 lme(response ~ 0 + condition*items,
      random = ~1|subject, data=d)

> where subject is the random effect. In this way, when I run:
> summary(f.lme.1)$tTable, I will get the following output:

[snip]
 
> together with Value, Std.Error, DF, t-value, p-value columns. I have two
> questions:
> 
> (1) If I want to compare Control vs. Treatment, shall I follow the steps:
> (a) run a smaller model
> 
> f.lme.0 =  lme(response ~ 0 + factor(condition) + factor(items), random =
> ~1|subject)
> 
> and (b) anova(f.lme.0, f.lme.1), and (c) if fail to reject, just use
> estimable() function in gmodels and make a contrast of (-1,1,0,0) in the
> smaller model?

  That is typical.

> (2) I am actually more interested in whether levels of items, i.e. E1, E2,
> E3 are different across condition, i.e. whether the interaction terms are
> significant (by just checking the p-value column?):
> 
> factor(condition)Treatment:factor(items)E2
> factor(condition)Treatment:factor(items)E3
> 
> However, how can I tell if the first interaction term:
> 
> factor(condition)Treatment:factor(items)E1
> 
> is significant or not? This term is not shown in the summary output and I
> think it has something to do with the contrast option used in R... How can
> I get the p-values for all interactions directly? Thanks a lot!

  You can fit

 lme(response ~ 0 + condition:items,
      random = ~1|subject, data=d)

If any of your conditions or items are repeated within subjects
you should consider testing an interaction with subjects (e.g.
see Schielzeth and Forstmeier 2009 _Behavioural Ecology_)