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

[Douglas Bates]

On 8/17/05, Shige Song <shigesong at 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.

Try the following after installing version 0.98-4 or later of the Matrix package

> 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.