10
$\begingroup$ 
subject <- factor(rep(c(1,2,3,4,5,6,7),each=4, times=2))
dep <- c(5,4,9,3,4,4,2,1,10,7,8,7,1,2,1,1,5,10,1,7,3,2,1,4,3,8,7,3,1,1,2,1,15,10,20,11,2,2,1,3,11,12,9,7,2,3,1,2,11,9,8,9,3,4,2,1) 

f1 <- factor(rep(c(rep("Female",times=16),rep("Male",times=12)), times=2))
f2 <- factor(rep(c("day1","day2","day3","day4"),times=14))

data <- data.frame(sub=subject, dep=dep, f1=f1, f2=f2)

m <- lmer(dep ~ f1*f2 + (1|sub), data=data)

I'm trying to understand how I can test the assumptions for mixed models.

1) In the case of model m, should I look at homogeneity of variances for every combination of f1 and f2 like this plot(resid(m)~fitted(.)|f1:f2) or is it enough to simply do plot(resid(m))?

2) In either case my real model is showing a funnel shape, is this too problematic? What would you do in that case?

3) How can I check the linearity assumption for each categorical variable in R?

4) Apart from homogeneity of variances and normality of residuals, are there any other important assumptions that you think I should be aware of?

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ In your example, your dependent variable contains only strictly positive integers. Is it perhaps a count variable in your real data? Then you should use a GLMM. $\endgroup$
    – Roland
    Commented Nov 12, 2018 at 8:37

2 Answers 2

Reset to default
12
$\begingroup$

  1. Existence of variance: Do not need to check, in practice, it is always true.

  2. Linearity: Do not need to check, because your covariates are categorical.

  3. Homogeneity: Need to Check by plotting residuals vs predicted values.

  4. Normality of error term: need to check by histogram, QQplot of residuals, even Kolmogorov-Smirnov test.

  5. Normality of random effect: Get the estimate of random effect (in your case random intercepts), and check them as check the residual. But it is not efficient because you just have 7 random intercepts.

    Another assumption is the independent between subjects. No test, based on your judgement. Subject specific random intercept means the correlation between the response variable from the same subject are the same.

Improve this answer
$\endgroup$
10
  • 3
    $\begingroup$ Normality of the random effects is more difficult to be checked because of shrinkage (i.e., the empirical Bayes estimates of the random effects will look more normal than the original random effects were); though this assumption is in general less important. $\endgroup$
    – Dimitris Rizopoulos
    Commented Nov 10, 2018 at 6:00
  • $\begingroup$ @DimitrisRizopoulos isn't OP's bigger problem that n_cases in their example is 7 such that normality is not verifiable? At large n_cases or lots of measurement per case, shrinkage will not affect one's ability to verify, right? $\endgroup$
    – Heteroskedastic Jim
    Commented Nov 10, 2018 at 13:22
  • 1
    $\begingroup$ @HeteroskedasticJim I agree that n_cases being 7 is a bigger issue here. Nonetheless, with regard to the point I made, having many repeated measurements per case makes it even harder to verify the assumption for the random effects distribution, but also in this case a misspecification will mater less. $\endgroup$
    – Dimitris Rizopoulos
    Commented Nov 10, 2018 at 13:33
  • 1
    $\begingroup$ @IsabellaGhement You can say in that way, but it is easy to be confused (duo marginal and conditional). The exact correct way is the variance of error term ($\epsilon$). $Var(\epsilon) = Var(Y|X) < \infty$ $\endgroup$
    – user158565
    Commented Nov 10, 2018 at 15:48
  • 2
    $\begingroup$ @HeteroskedasticJim the random effects are typically estimated as the modes/means of the posterior distribution p(b_i | y_i), with b_i denoting the random effects, and y_i the response vector for case i. Now, this posterior is proportional to p(y_i | b_i) p(b_i), with p(b_i) denoting the assumed dist for the random effects. However, as n_i increases, this posterior is dominated by the first term p(y_i | b_i) as a function of b_i, which converges to a normal distribution (i.e., with enough observations the prior doesn’t matter). The rate of convergence depends on how “normal” p(y_i | b_i) is. $\endgroup$
    – Dimitris Rizopoulos
    Commented Nov 10, 2018 at 18:08
9
$\begingroup$

The commonly quoted assumptions (or "conditions" as I prefer to call some of them) of linear mixed effects models are:

  • Linearity of the predictors. This can be checked by plotting the residuals against the response and looking for any systematic shape, and by including non-linear terms (or splines) and comparing the model fit. Very often this will not be an issue, and if it is, then including non-linear terms (such as log, exp or polynomials) in the linear predictor may be sufficient. More importantly, substantive domain knowledge should inform whether the linearity condition is justified. For example, in some domains such as pharmacokinetics, we already know, based on rigorous theory and experimentation, that a linear model is not appropriate in some cases.

  • The residuals have constant variance. This can be checked with a plot of residuals against fitted values - there should be no pattern/trend.

  • The residuals are independent. This can be checked by plotting residuals against covariates - especially time-varying or spatial covariates. There should not be any systematic pattern

  • The residuals are normally distributed. This can be checked in many ways, such as a Q-Q plot and a simple histogram. Statistical tests, such as Anderson-Darling and Kolmogorov–Smirnov are also possible.

Note that "residuals" above refer to both the unit-level residuals (often called "errors") and the random effects. For the random effects, this can be problematic where only a small number of groups/clusters exist in the sample. In the simulated example given in the OP there are 7 clusters. There are lots of rules of thumb to inform a sufficient number of clusters, and 7 is generally thought to be difficult to draw any conclusions.

It is mentioned in the OP that their actual model exhibits a funneled shape plot of residuals vs fitted values. This indicates heteroskadasticity. One way to proceed is to consider a transformations of variables - ideally this should be informed by expert domain knowledge. With this in mind, Box-Cox transformations may be useful.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Thanks for the answer @Robert Long, but I don't understand the part "residuals above refer to both the unit-level residuals (often called "errors") and the random effects". Are you referring to the last assumption specifically? I checked the normality of residuals with qqnorm(resid(m)) but how can I test the normality of random effects in R? $\endgroup$
    – locus
    Commented Nov 13, 2018 at 22:58
  • 1
    $\begingroup$ @locus if using lme4 you can extract the random effects with ranef as shown here $\endgroup$
    – Robert Long
    Commented Nov 14, 2018 at 5:35
  • 1
    $\begingroup$ thank you very much! Unfortunately I cannot accept both answers, but I upvoted yours as well $\endgroup$
    – locus
    Commented Nov 14, 2018 at 23:14
  • $\begingroup$ BTW, which test do you prefer to test whether residuals are normally distributed? I know the ks.test() function but I'm not sure what the y should be when I fit a mixed model. $\endgroup$
    – locus
    Commented Nov 15, 2018 at 0:07
  • $\begingroup$ @locus the type of model in itself does not determine the test. The main point here is that you have rather few observations to test so K-S will not have much power. You can run it anyway, but I would definitely also plot a histogram and a Q-Q plot and use those too. To satisfy yourself a little more, you could run a lot of simulations, drawing 7 samples from a normal distribution and observe the many different shapes of the plots that are possible ! $\endgroup$
    – Robert Long
    Commented Nov 15, 2018 at 8:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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