Correlation between standardized residuals and fitted values in a linear mixed effect model: Course of action?

Question

I am fitting a linear mixed effect model in R with lme from nlmer, using the approach described in Zuur et al. "Mixed Effects Models and Extensions in Ecology with R".

As a first step, I use gls to fit a linear model and look for evidence of heterogeneity, starting with a full model.

require(nlme)
M0 <- gls(ptltcued ~ facesex*emo*agegroup, data=data, na.action = na.exclude)
plot(M0)

I find no (?) evidence of heterogeneity.

Next, I use AIC and REML fits to choose between models with no random effect, a random intercept, or random intercept and slope. I have 2 fixed between-subject factors (agegroup, emo), 1 fixed within-subject factor (facesex), and one random subject factor (numsubj).

M1 <- gls(ptltcued ~ facesex*emo*agegroup, method="REML", na.action= na.exclude, data=data)
M2 <- lme(ptltcued ~ facesex*emo*agegroup, random=~1|numsubj,method="REML", na.action= na.exclude, data=data) 
M3 <- lme(ptltcued ~ facesex*emo*agegroup, random=~1+facesex|numsubj, method="REML", na.action= na.exclude, data=data) 
AIC(M1,M2,M3)

AIC gives the following result:

• M1 df 13 AIC -215.2172
• M2 df 14 AIC -213.2172
• M3 df 16 AIC -221.1735

Based on AIC, I decide for a model with random intercept and slope. Next, I validate M3 before going on to select fixed effects.

Normality checks look good (?).

hist(resid(M3))

qqnorm(resid(M3))
qqline(resid(M3))

Independence checks look good (?), for example here for facesex:

plot(data$facesex,resid(M3))

Heterogeneity check looks... well, hum. I have never seen this kind of pattern in a plot of standardized residuals versus fitted values. There is no change in spread along the fitted values, but there is an obvious correlation between residuals and fitted values.

plot(fitted(M3),resid(M3))
abline(h=0,col="grey")
lines(lowess(fitted(M3)[is.finite(fitted(M3))],resid(M3)[is.finite(fitted(M3))]),col="red")

The pattern is absent in M1 and M2. For example, here is the plot for M2:

I don't understand the reason behind this. It seems to me that I should abandon M3 and go for M1, the next best in terms of AIC. However, I have a feeling that this kind of obvious relationship probably has an obvious cause such as a formula mistake or something else which I am unaware of.

So: What should be the course of action here? Is there an obvious reason for this pattern?

Answer 1

OK based on this question/answer which states that

A strong correlation [between residuals and fitted values] is not necessarily cause for alarm. This may simply means the underlying process is noisy. However, a low R2 (and hence high correlation between error and dependent) may be due to model misspecification.

and because I can't see any misspecification in my formula I decided in favor of M3. Fortunately M1, M2 and M3 give the same conclusions in terms of which fixed effects are statistically significant.

It seems that the high correlation is simply due to the low prediction value of my factor variables. Not sure exactly why this wasn't evident in M0, M1 and M2 but then again the range of fitted values was also much lower.