I am currently trying to analyse data from an experiment of mine and I have done some searching for instructions on the usage of the lme() function for R, since I am looking to analyse my data with a linear mixed effects approach. However,

  1. In case lme() does work out, I am not sure how to fit a model to my data.
  2. In case there are better options than lme(), I would be glad about any help regarding the proper use of these functions.

Subjects in my experiment were randomly assigned to one of two contexts (between subjects factor B) and within this context all of the subjects completed a task which was constructed as a combination of two within subjects factors (W1 and W2). Additionally, subjects have filled a personality questionnaire, thus each subjects has a value for this, which I need to include in my analysis (Pers). The outcome variable is a continuous one (reaction time, RT). The data look like this:

B   Pers    subject W1  W2  RT
A   23      1       m   x   187
A   23      1       m   y   333
A   23      1       n   x   112
A   23      1       n   y   110
A   5       2       m   x   313
... ...     ...     ... ... ...
B   11      19      m   x   911
B   11      19      m   y   248
B   11      19      n   x   411
B   11      19      n   y   212
B   -9      20      m   x   666
... ...     ...     ... ... ...

I am interested in an analysis giving me information about the interactions of all four factor (B, Pers, W1 and W2), thus four-way interaction.

When including only B, it seems to me that I understand what to do when using lme(). My code looks like this (nb: I constructed contrasts between the different conditions before the actual lme() code using the contrasts() function):

> baselinemodel = lme(RT ~ 1, random = ~1|subject/W1/W2, data = df, method = "ML")
> W1model = update(baseline, .~. +W1)
> W2model = update(W1model, .~. +W2)
> Bmodel = update(W2model, .~. +B)
> W1W2mod = update(Bmodel, .~. +W1:W2)
> W1Bmod = update(W1W2mod, .~. +W1:B)
> W2Bmod = update(W1Bmod, .~. +W2:B)
> fullmodel = update(W2bmodel, .~. +W1:W2:B)

I.e. I set up a basline model including only the intercept and specifying W1 and W2 as nested within subjects (hence defining them as within subjects factors) and then in a sequential manner I include more predictors and their interactions in order to finally be able to compare these models using the anova() function:

> anova(baseline, W1model, W2model, Bmodel, W1W2mod, W1Bmod, W2Bmod, fullmodel)

Naturally I also take a look at my predefined contrasts using

> summary(fullmodel)

I find this approach very intuitive, as long as I do not try to include the continuous predictor Pers in addition to B. So this is where I am stuck. I would be deeply grateful if anyone could

  1. explain to me how to include Pers in addition to B, W1 and W2 (i.e. how to include a continuous predictor varying between subjects, and its interactions with the other factors into my analysis)

  2. and in case lme() is not the optimal function for my objectives how to do a comparable analysis with another function (thus broadening my still quite narrow scope of R)

I have, naturally done some searching before starting to post this question, but unfortunately others, who seemed to have related issues weren't too successful in getting answers (https://stats.stackexchange.com/questions/97669/analysis-of-a-mixed-design-with-categorical-and-continuous-variables)

Thank you in advance


1 Answer 1


I am assuming that W1 and W2 are both fixed factors, and that all subjects were exposed to all levels of W1 and W2 (i.e., subjects, W1, and W2 are all crossed, not nested). If this is the case, you can use either lme or lmer in the lme4 package:

fullmodel_lme  <-lme(RT ~ W1*W2*B2*Pers, random =~1|subject, data=df, method="ML")
fullmodel_lmer <- lmer(RT ~ W1*W2*B2*Pers + (1|subject), data=df, REML=FALSE)

By including W1*W2*B2*Pers you are testing the hypotheses that effect of personality on RT is moderated by the other predictors as well as the combinations thereof. You can use the sequential approach as you did. If you use lmer you can also use the lmerTest package to get estimates of degrees of freedom and p-values for the coefficients.

If W1 and W2 are random factors, then I suggest you use lmer since it is much easier to fit crossed models using lmer than lme.

  • $\begingroup$ Thank you very much for your answer. The code suggested by you does, in fact, look quite simple. I will immediately try this out. Thank you also for the hint regarding lmerTest(). I had heard before that lmer() does not automatically give out p-values and dfs. However, am I right in assuming that lmer() automatically differentiates between within and between subjects factors by indication of the subject variable (i.e. "subject" in my case) in the term (1|subject) and the fact that B and Pers do not vary within subjects, whereas W1 and W2 do? $\endgroup$
    – bunsenbaer
    Commented Sep 2, 2014 at 15:55
  • $\begingroup$ Correct, as long as random factors are specified and data are coded correctly, you do not have specify which is a between- or within-subject factor. $\endgroup$ Commented Sep 2, 2014 at 16:07
  • $\begingroup$ Good answer, but note that crossing is only hard when it is of random effects. lme(RT ~ W1*W2*B2*Pers, random =~1|subject, data=df, method="ML") should work fine (although lmer is probably faster for large data sets). $\endgroup$
    – Ben Bolker
    Commented Sep 2, 2014 at 21:58
  • $\begingroup$ Thanks for the answers so far. I have now reached a point at which I think that it might be useful to add a covariate (intelligence) to my model for the sake of increasing the power by reducing unexplained variance (as in ANCOVA). However, I do not quite understand how to do this. I think that it should be correct to use the lmer model proposed by Masato Nakazawa and adding intelligence in this manner lmer(RT~W1*W2*B*Pers+intelligence+(1|subject), data=df, REML =FALSE). Is this correct? What is the difference if instead of this version I added intelligence as (1+intelligence|subject)? $\endgroup$
    – bunsenbaer
    Commented Sep 9, 2014 at 17:55
  • 1
    $\begingroup$ Ok, so I have noticed that the option of including intelligence as (1+intelligence|subject)is not useful, since lmer would try to calculate a slope for intelligence for each random effect (subject). Since there is no variation in intelligence within each subject, random slopes would have to be 0 and covary perfectly with the random intercepts, which is evident in the random effectspart of the lmer-output. Also, I have found this extremely helpful thread and this thread for lmer basics. $\endgroup$
    – bunsenbaer
    Commented Sep 10, 2014 at 8:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.