After attempting to produce a linear mixed model, I was left with a great deal of heterogeneity.

lme1 <- lme(Average.payoff ~ Game + Type + Others.Type + Game:Type + 
  Game:Others.Type + Type:Others.Type, random=~1|Subjects, method="REML", 

The response term Average.payoff is continuous, whilst all explanatory variables are binary.

When I look to validate, I can clearly see that the spread of the residuals decreases with the larger fitted values.

I decided alternatively to see what would happen if I just fit using gls().

gls1 <- gls(Average.payoff ~ Game + Type + Others.Type + Game:Type + 
  Game:Others.Type + Type:Others.Type, method="REML", data=Subjectsm1)
anova(lme1, gls1)
#         Model df      AIC      BIC    logLik   Test   L.Ratio p-value
#    lme1     1  9 67.81805 81.28662 -24.90902                         
#    gls1     2  8 66.14661 78.11867 -25.07330 1 vs 2 0.3285588  0.5665

As you can see the gls() model gives a better fit than the lme() model. Even though I know from the experimental design that I have random effects, am I justified to remove them from the model if the fit is better without them?

  • 4
    $\begingroup$ How do you conclude that the "gls model" has a better fit? If you base your decision on AIC, then be aware that this MIGHT be dangerous. Indeed, "lme" and "gls" do not necessarily use the same log-likelihood: it may be that one of them left all additive constants out of the log-likelihood while the other include them; and this would compromise the comparison of AIC's. $\endgroup$
    – ocram
    Aug 22, 2012 at 17:19
  • 1
    $\begingroup$ No, you aren't, unless you can show that the assumptions of gls are met. Usually, the reason for adapting a mixed model is that the assumptions of other models are not met. $\endgroup$
    – Peter Flom
    Aug 22, 2012 at 18:26
  • $\begingroup$ @Jonathan as an aside, you don't need all those extra terms in the model since R will automatically expand interactions specified with a * (eg. lm(y ~ x*z) is the same as lm(y ~ x + z + x:z) and helps make the code a little easier to follow). $\endgroup$
    – smillig
    Aug 23, 2012 at 12:32
  • $\begingroup$ @ocram I used both AIC and loglik. I have read elsewhere that the gls() and lme() likelihoods are the same (but they are not for lm()).Bearing in mind the experiment involved 20 subjects, for each subject I have two samples, one where Others.Type=0 and one where Others.Type=1, do you have any further thoughts?? $\endgroup$ Aug 23, 2012 at 15:50
  • 1
    $\begingroup$ @Jonathan The model you want is Average.payoff ~ Game*Type + Game*Others.Type + Type*Others.Type which won't have any three-way interactions. $\endgroup$
    – smillig
    Aug 23, 2012 at 16:51

1 Answer 1


An initial issue: when estimating the final parameters of a ME model, the restricted maximum likelihood (ReML) is used; as you have already done. ReML tries practically to find linear combinations of the responses, $k$, such that $k^TX= 0$ and thus to exclude any fixed terms parameters from the likelihood function ($X$ : $n\times p$ model matrix). However, ML should be used for the model selection procedure as the theory for model comparisons is based on ML estimation. As ReML will try to transform the fixed effect responses in the manner described above, this would lead to a series of different transformations for each model setting, making them incomparable. Therefore one should use ML estimators if likelihood ratio tests are to be employed (AIC etc.).

(Chapt.4 in Applied Longitudinal Analysis by Fitzmaurice et al. (towards the end of section 4.5 - I have the book but I just checked it and that section seems available in Google books also.))

  • $\begingroup$ thanks, strangely the book I am heavily using (Mixed Effects Models and Extensions in Ecology with R by Alain Zuur) says that I SHOULD be using REML for log likelihood ratio testing. I have though read other sources that agree with your suggestion. Is this a debated topic?? $\endgroup$ Aug 23, 2012 at 17:07
  • 1
    $\begingroup$ Em... not really debated... Putted simply (and paraphrasing Julian Faraway's "Extending the Linear Model with R", section 8.4) : If you compare two LME models with different fixed effects, the linear combinations, $k$, of them will be different and the obtained maximum likelihoods will be incomparable. I haven't seen Zuur's book but are you sure he is referring to model comparison? In your case it is a bit shady cause you have actually the same $X$ matrix in the LME and the GLS model but even then you assume uniqueness of $k$. $\endgroup$
    – usεr11852
    Aug 23, 2012 at 22:07
  • 1
    $\begingroup$ I too am using the ecology book and I think you should check as it says for testing fixed effects using use ML specifically when using the log likelihood test. As mine is an e-book I can not give you the page number. The REML should be run when you are happy with your fixed effects model. $\endgroup$
    – user20925
    Feb 17, 2013 at 2:59

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.