Checking assumptions lmer/lme mixed models in R

Question

I ran a repeated design whereby I tested 30 males and 30 females across three different tasks. I want to understand how the behaviour of males and females is different and how that depends on the task. I used both the lmer and lme4 package to investigate this, however, I am stuck with trying to check assumptions for either method. The code I run is

lm.full <- lmer(behaviour ~ task*sex + (1|ID/task), REML=FALSE, data=dat)
lm.full2 <-lme(behaviour ~ task*sex, random = ~ 1|ID/task, method="ML", data=dat)

I checked if the interaction was the best model by comparing it with the simpler model without the interaction and running an anova:

lm.base1 <- lmer(behaviour ~ task+sex+(1|ID/task), REML=FALSE, data=dat)
lm.base2 <- lme(behaviour ~ task+sex, random= ~1|ID/task), method="ML", data=dat)
anova(lm.base1, lm.full)
anova(lm.base2, lm.full2)

Q1: Is it ok to use these categorical predictors in a linear mixed model?
Q2: Do I understand correctly it is fine the outcome variable ("behaviour") does not need to be normally distributed itself (across sex/tasks)?
Q3: How can I check homogeneity of variance? For a simple linear model I use plot(LM$fitted.values,rstandard(LM)). Is using plot(reside(lm.base1)) sufficient?
Q4: To check for normality is using the following code ok?

hist((resid(lm.base1) - mean(resid(lm.base1))) / sd(resid(lm.base1)), freq = FALSE); curve(dnorm, add = TRUE)

Answer 1

Q1: Yes - just like any regression model.

Q2: Just like general linear models, your outcome variable does not need to be normally distributed as a univariate variable. However, LME models assume that the residuals of the model are normally distributed. So a transformation or adding weights to the model would be a way of taking care of this (and checking with diagnostic plots, of course).

Q3: plot(myModel.lme)

Q4: qqnorm(myModel.lme, ~ranef(., level=2)). This code will allow you to make QQ plots for each level of the random effects. LME models assume that not only the within-cluster residuals are normally distributed, but that each level of the random effects are as well. Vary the level from 0, 1, to 2 so that you can check the rat, task, and within-subject residuals.

EDIT: I should also add that while normality is assumed and that transformation likely helps reduce problems with non-normal errors/random effects, it's not clear that all problems are actually resolved or that bias isn't introduced. If your data requires a transformation, then be cautious about estimation of the random effects. Here's a paper addressing this.

Answer 2

You seem quite mislead about the assumptions surrounding multi-level models. There is not an assumption of homogeneity of variance in the data, just that the residuals should be approximately normally distributed. And categorical predictors are used in regression all of the time (the underlying function in R that runs an ANOVA is the linear regression command).

For details on examining assumptions check out the Pinheiro and Bates book (p. 174, section 4.3.1). Also, if you plan to use lme4 (which the book isn't written around) you can replicate their plots using plot with an lmer model (?plot.merMod).

To quickly check normality it would just be qqnorm(resid(myModel)).

Answer 3

Regarding Q2:

According to Pinheiro and Bates' book you may use the following approach:

"The lme function allow the modeling of heteroscesdasticity of the within-error group via a weights argument. This topic will be covered in detail in § 5.2, but, for now, it suffices to know that the varIdent variance function structure allows different variances for each level of a factor and can be used to fit the heteroscedastic model [...]"

Pinheiro and Bates, p. 177

If you would like to check for equal variances between sex you may use this approach:

plot( lm.base2, resid(., type = "p") ~ fitted(.) | sex,
  id = 0.05, adj = -0.3 )

If variances are different, you can update your model in the following manner:

lm.base2u <- update( lm.base2, weights = varIdent(form = ~ 1 | sex) )
summary(lm.base2u)

Further more, you may have a look at the robustlmm package which also uses a weighing approach. Koller's PhD thesis about this concept is available as open access ("Robust Estimation of Linear Mixed Models"). The abstract states:

"A new scale estimate, the Design Adaptive Scale estimate, is developed with the aim to provide a sound basis for subsequent robust tests. It does so by equalize the natural heteroskedasticity of the residuals and to adjust for the robust estimating equation for the scale itself. These design adaptive corrections are crucial in small sample settings, where the number of observations might be merely five times the number of parameters to be estimated or less."

---

---

I do not have enough points for comments. I see however the necessity to clarify some aspect of @John 's answer above. Pinheiro and Bates state on p. 174:

Assumption 1 - the within-group errors are independent and identically normally distributed, with mean zero and variance σ2, and they are independent of the random effects.

This statement is indeed not clear about homogeneous variances and I am not deep enough into statistics to know all the maths behind the LME concept. However, on p. 175, §4.3.1, the section dealing with Assumption 1 they write:

In this section, we concentrate on methods for assessing the assumption that the within-group errors are normally distributed, are centered at zero, and have constant variance.

Also, in the following examples "constant variances" are indeed important. Thus, one may speculate whether they imply homogeneous variances when they write "identically normally distributed" on p. 174 without addressing it more directly.

Answer 4

Q1: Yes, why not?

Q2: I think the requirement is that the errors are normally distributed.

Q3: Can be tested with Leven's test for example.