# Testing the variance component in a mixed effects model

Say $y=X\beta+ Zu +\epsilon$ is our mixed effects model where $u=(u_1,..,u_r)$ and $u_{j} \stackrel{i.i.d.}{\sim} N(0, \sigma^2_{a})$ for $j=1,...,r$ and $\epsilon=(\epsilon_1,...,\epsilon_n)$ are i.i.d. $N(0, \sigma^2_{b})$, furthermore $\epsilon_j$ and $u_i$ are also assumed to be independent for all $j$'s and all $i$'s.

I am interested in testing the hypothesis $H_{0}:\sigma_{a}^2=0$ vs $H_{1}: H_{0}$ is not true. The ${\bf lmer}$ package in R does provide a bootstrap p-value for this problem, but my question is: Is there is any non-bootstrap way of obtaining this p-value (asymptotic or otherwise)?

This is usually done with maximum likelihood ratio between original model and a model omitting the variance coefficient to be estimate (random intercept/random slope/random co-variance between slope and intercept).

A good example is in these tutorials:

Sample R code:

> model1 = lmer(resp ˜ fixed1 + (1 | random1))
> model2 = lm(resp ˜ fixed1)
> chi2 = -2*logLik(model2, REML=T) +2*logLik(model1, REML=T)
> chi2
[1] 5.011
> pchisq(chi2, df=1, lower.tail=F)
[1] 0.02518675


Asymptotic test are problematic for variance parameters, because parameter space is bounded by zero. Moreover, the hypothesis you are trying to test, can't be true, as the parameter is continuous. Probability of $\sigma^2 = 0$ is exactly 0.

What you can do to make inference on the variance parameters is to switch to a Bayesian implementation, where you would get the full posterior distribution for the variance parameters. For lme4 users, the MCMCglmm package is easy to learn. You could also use JAGS or Stan. For an example, where Stan was used to compare several random effects, see [1].

[1] Schmettow, M., & Havinga, J. (2013). Are users more diverse than designs? Testing and extending a 25 years old claim . In S. Love, K. Hone, & Tom McEwan (Eds.), Proceedings of BCS HCI 2013- The Internet of Things XXVII. Uxbridge, UK: BCS Learning and Development Ltd.

• Thanks. Why cant the hypotheis be true? We test $H_{0}:\theta=$ constant vs its non constant all the time in statistics (where this $\theta$ could be the mean or any other parameter. Thanks for the Bayesioan insight But what about in the frequentist setup? I ask this because I was told that SAS provides the (asymptotic) p-values so was wondering if there was any R equivalent for that. Jul 14 '15 at 15:32