2
$\begingroup$

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)

residual versus fitted for 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))

Histogram of residuals

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

QQ plot of residuals Independence checks look good (?), for example here for facesex:

plot(data$facesex,resid(M3))

Residuals versus Xi for factor facesex

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")

Resid vs fitted for M3

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

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

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

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.