# Should I keep or remove random effects?

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",
data=Subjectsm1)


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?

-
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. – ocram Aug 22 '12 at 17:19
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. – Peter Flom Aug 22 '12 at 18:26
@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). – smillig Aug 23 '12 at 12:32
@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?? – Jonathan Bone Aug 23 '12 at 15:50
@Jonathan The model you want is Average.payoff ~ Game*Type + Game*Others.Type + Type*Others.Type which won't have any three-way interactions. – smillig Aug 23 '12 at 16:51
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.).
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$. – user11852 Aug 23 '12 at 22:07