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)
  • $\begingroup$ One thing I noticed is also that the version of lme4 I was using was not the most recent one and therefore the simple plot(myModel.lm) did not work, perhaps this is helpful for other readers to know.. $\endgroup$ – crazjo Nov 28 '13 at 19:31

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.

  • $\begingroup$ Thanks for your answer. I would like to share my dataset and script for analysis including output to see if what I did is indeed correct. Is it possible in stack exchange? Furthermore I think I ran the wrong random factor (1|rat/task), shouldn't it just be (1|rat)? I tested 60rats (30 of each sex) on three tasks. $\endgroup$ – crazjo Nov 28 '13 at 19:30
  • 9
    $\begingroup$ I tried the code for Q4 recently and I got an error about object of type 'S4' not subsettable. Was that code intended for models fit with the lme package? What about with lme4? $\endgroup$ – emudrak Sep 22 '14 at 13:40
  • $\begingroup$ Regarding Q4, people making those plots need to keep in mind that the N for each of the plots produced will be substantially smaller than the total and therefore the plots will be much more variable. Don't expect them to look as consistently normally distributed as the overall. $\endgroup$ – John Dec 20 '15 at 23:05

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

  • $\begingroup$ Thanks for your comment. Do you suggest using the lmer over the lme4 method? And am I right in understanding the response variable does not need to be normally distributed? I will have a proper read through the Pinheiro and Bates book. $\endgroup$ – crazjo Nov 28 '13 at 11:00
  • $\begingroup$ Also, are you sure running qqnorm(resid(myModel)) on a mixed model with multiple factors works? $\endgroup$ – crazjo Nov 28 '13 at 12:07
  • $\begingroup$ The newer lmer function has more capabilities and higher performance. Have you tried qqnorm? Follow the advice at the beginning of the book on how to read it. $\endgroup$ – John Nov 28 '13 at 15:04
  • $\begingroup$ The plot I had initially looked weird, possibly because I indeed did not have the newest version of lmer. Thanks for noting this, it now works as needed. $\endgroup$ – crazjo Nov 28 '13 at 19:33

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

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.


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.


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.