# Re: [R] How to assess significance of random effect in lme4

From: Douglas Bates <dmbates_at_gmail.com>
Date: Sat 20 Aug 2005 - 02:44:22 EST

On 8/17/05, Shige Song <shigesong@gmail.com> wrote:
> Dear All,
>
> With kind help from several friends on the list, I am getting close.
> Now here are something interesting I just realized: for random
> effects, lmer reports standard deviation instead of standard error! Is
> there a hidden option that tells lmer to report standard error of
> random effects, like most other multilevel or mixed modeling software,
> so that we can say something like "randome effect for xxx is
> significant, while randome effect for xxx is not significant"? Thanks!

Sorry to be entering this discussion late. I know there have been several responses and a considerable amount of discussion following the original question but I think it would be worthwhile returning to it. Shige asks why standard errors of the estimates of the variance components are not reported. All other known software packages for fitting mixed-effects models report the estimates of variances and covariances of the random effects and the standard errors of these estimates. However, the summary for a fitted lmer model or a fitted lme model does not.

This is intentional. It is not an oversight.

The reason these are not reported is because they are not useful and, in fact, are quite misleading. It is useful to have a standard error of an estimate when the distribution of the estimator is approximately symmetric. Then the standard error can be used to create an approximate confidence interval or perform a hypothesis test. However, the distribution of the estimate of a variance generally looks like a scaled chi-squared.

> library(Matrix)
> Sm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)

> Ss1 <- mcmcsamp(Sm1, 10000, trans = FALSE)
> str(Ss1)

mcmc [1:10000, 1:6] 238 252 253 251 245 ...

• attr(*, "class")= chr "mcmc"
• attr(*, "mcpar")= int [1:3] 1 10000 1
• attr(*, "dimnames")=List of 2 ..\$ : NULL ..\$ : chr [1:6] "(Intercept)" "Days" "sigma^2" "Sbjc.(In)" ...

This fits an lmer model to a longitudinal data set then creates a Markov Chain Monte Carlo sample from the parameters of the model. The default is to create an MCMC sample of transformed parameters that correspond to the logarithm of any variances and Fisher's 'z' transformation of the correlations. The "trans = FALSE" optional argument allows for the sample to be on the scale of the variances and the covariances.

Check the density plots or the normal probability plots of the variance components and the covariance. The distribution is not at all normal or even symmetric. A symmetric confidence interval based upon a standard error or, even worse, a hypothesis test based on the estimate of a variance component and its standard error are clearly inappropriate.

BTW, that mcmcsamp function is frighteningly fast when applied to a linear mixed model. It takes less than 2 seconds to generate a sample of size 10000 on the machine that I use.

R-help@stat.math.ethz.ch mailing list