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
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?