Test to obtain p-value of the random effect in mixed models I would like to obtain the p-value of the random effect for a simple linear mixed effect model.
To be more precise, my model is:
$$\mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\mathbf{Z}\mathbf{u}+\boldsymbol{\epsilon} $$,
where $\boldsymbol{\epsilon}\sim \mathcal{N}(0,\sigma \mathbf{I}_n)$ and $\mathbf{u}\sim\mathcal{N}(0,\tau \mathbf{I}_q)$. And I want to test whether $H_0: \tau=0$ or $H_1: \tau\neq0$. In case this is relevant, my matrix $\mathbf{Z}$ is matrix of 0 and 1 to represent the "random effect of the $q$ individuals" in my dataset of $n$ datapoints (not sure if I am clear enough).
I know this question has been directly or indirectly asked a few times in this forum however the ones I saw where pretty old so maybe there is some improvment in the field that I did not found. I also know that there is some debate whether we should even test the p-value of the random effects to maybe exclude them. And from what I understand, obtaining correct p-values for the fixed effects in mixed models is not easy/possible so trying it for the random effect might not be smart.
However I would still like to know if there is a good procedure to obtain those p-values. If it is not possible to obtain the "correct" p-values, I would like to know if there is some way to obtain "decent" p-values.
After some reasearch I have found some people that suggest using a LR test by fitting two models with lm in R and lmer function from the lme4 package:
m1: $ \mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}$ and
m2: $ \mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\mathbf{Z}\mathbf{u}+\boldsymbol{\epsilon}$,
then we compare them using anova(m1,m2). I did simulations in a simple case to check under the null hypothesis ($\tau=0$) to see whether the p-values we obtain are acceptable. And it seems those p-values are not acceptable, the p-values are not uniformly distributed at all. The p-values seem to be skewed toward 1. In this the answer of this question, someone did similar simulations as me and found similar results (although I have to admit they don't test exactly the same thing as me).
I also saw that I could use the RLRsim package, but it seems that the p-value are also based on the LR test so I haven't even bothered to try to be honest.
Do you know if there is a good test to get a "correct" p-value for $\tau=0$? Why does the LR test method seem to perform poorly?

The random intercept may be referred to as a variance component, and there is a broad literature on testing variance components. A single random intercept can be well approximated by a compound symmetry covariance matrix, like is seen in repeated measure ANOVA. To test whether a variance component is equal to 0 introduces some theoretical challenges, as the actual value of the variance component lies on the boundary of the parameter space. According to Stram and Lee, the asymptotic distribution of the likelihood ratio statistic for the nested models when the null hypothesis is true, i.e. the random-intercept variance when the data are distributed IID is given by:
$$ -2 \log LR \rightarrow_d 0.5 \chi^2_1 + 0.5 \chi^2_0$$
That is it is a mixture distribution with distribution function:

To verify this, a quick simulation may convince:
library(nlme)
n <- 100
i <- rep(1:10, 10)

lrs <- rep(0, 100)
for(j in 1:100) {
  y <- rnorm(n)
  dat <- groupedData(y ~ 1 | i)
  fm1 <- lme(y ~ 1, data=dat, correlation = corCompSymm( form = ~1|i))
  fm0 <- gls(y ~ 1 ,data=dat)
  lrs[j] <- 2*(logLik(fm1) - logLik(fm0))
}
plot(ecdf(lrs))


The p-value can be found by calculating 1 minus the DF for the sample LRT statistic. Using the above example, and a simulation with some intraclass correlation,
y <- rnorm(n) + 0.25*rnorm(n)[i]
dat <- groupedData(y ~ 1 | i)
fm1 <- lme(y ~ 1, data=dat, correlation = corCompSymm( form = ~1|i))
fm0 <- gls(y ~ 1 ,data=dat)
f <- function(q) 0.5*pchisq(q, 0 ) + 0.5*pchisq(q, 1)
lrs <- 2*(logLik(fm1) - logLik(fm0))
(1-f(lrs))*2

Gives a 2 sided p-value of:
> (1-f(lrs))*2
'log Lik.' 0.0409434 (df=4)

