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

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?

• Have you check the correlation between response and fitted value? – Kigo Chen Mar 17 '17 at 1:55
• Do you perhaps mean "heteroscedasticity", not heterogeneity? – Rose Hartman Mar 17 '17 at 3:10