# Test for effect of groups in a mixed effects model

This model is a simple linear regression:

mtcars_lm <- lm(mpg ~ wt, mtcars)


And this model adds cyl as a random effect:

library(lme4)
mtcars_mixed_effects <- lmer(mpg ~ wt + (1 | cyl), mtcars)


Is there a way to test whether adding cyl as random effect is worthwhile? I tried this but it threw an error:

anova(mtcars_mixed_effects, mtcars_lm)


(please disregard the fact that cyl only has three groups, I'm just using one of R's built in datasets to make question reprodicible).

Yes there is, in the package RLRsim:

install.packages("RLRsim")
library("RLRsim")

set.seed(1)
exactRLRT(mtcars_mixed_effects)

simulated finite sample distribution of RLRT.

(p-value based on 10000 simulated values)

data:
RLRT = 6.0969, p-value = 0.0034


... indicating (p-value = 0.0034) that the variance of the random effect is unequal zero so that the inclusion is worthwile

The R package nlme will let you do a simple comparison with the anova command, but will report the wrong p-values.

The reason is that when testing for the inclusion of a random effect, you are testing if some variance is equal to zero, i.e $H_0:\sigma_A^2=0$ which is on the boundary of the parameter space (as variances must be larger than zero). Asymptotic theory for LR testing requires the hypothesized value to be on the inside of the parameter space.

The correct p-value is computed from a mixture of two $\chi^2$-distributions. When including only a random intercept the p-value is computed from a 50/50 mixture of the $\chi^2_0$ (all weight on zero) and $\chi^2_1$ distributions.

In your case, this code will do what you want:

library(nlme)
fit1 = gls(mpg ~ wt, data=mtcars,method="REML")
fit2 = lme(mpg ~ wt,random=~1|cyl,data=mtcars,method="REML" )

> anova(fit1,fit2)
Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit1     1  3 164.8199 169.0235 -79.40996
fit2     2  4 160.7230 166.3278 -76.36149 1 vs 2 6.096934  0.0135


I use the gls command to estimate the normal linear model using REML, which makes the models comparable by LR. Like I said though, the p-value is incorrect. A correct p-value is found by:

p_value = 0.5*(pchisq(6.096934,df=0,lower.tail=F)+pchisq(6.096934,df=1,lower.tail=F))
> p_value
[1] 0.006770832


The low p-value indicating that the inclusion is worthwhile, the variance in question is larger than zero.

For a general inclusion of mixed effects, the mixture of distributions will change, but I'm sure there are packages that will do this automatically for you.

The approach popularized by many authors is to judge if to include a random effect or not by using intraclass correlation. To compute it you use variances of all the random effects ($\sigma^2_0$, $\sigma^2_1$, ...) and variance of residuals ($\sigma^2_r$):

$$\rho_0 = \frac{\sigma^2_0}{\sigma^2_0 + \sigma^2_1 + ... + \sigma^2_r}$$

However, this approach has its flaws and may give you misleading results so is discouraged by others - check this question.

In most cases I know you include random effects because you are interested in estimating their influence - so it is not a question about "if" to include them, but rather "how big" they are. If you re-frame the question like this, then you can decide to drop the random effects if their values are (relatively) very small - that is, including them in the model does not change much in the estimates.