[R] standardized random effects with ranef.lme()

[Doran, Harold]

> > Why do the results differ although the estimates (random 
> effects and 
> > thus their variances) are almost identical? I noticed that 
> lme() does 
> > not compute the standard errors of the variances of the 
> random effects 
> > - for several reasons, but if this is true, how does 
> ranef() calculate 
> > the standardized random effects (the help says: 
> '"standardized" (i.e.
> > divided by the corresponding estimated standard error)').
> > 
> > Is there a way to obtain similar results as in MLWin (or: should I 
> > prefer the results of ranef() for certain reasons)?

I think there are two different issues here. The lme function does not
produce a standard error of the variance component as some other
multilevel packages do. It is often recommended in the multilevel
literature to consider the p-value of the variance components and "fix"
or retain the variance if p < .05. There are good reasons not to follow
this practice.

If you were using lmer(), you still wouldn't get this statistic, but you
could use the MCMCsamp() function to examine the distribution of the
random effects.

The second issue (I think) is that the conditional variance of the
random effect is not the same as the standard error of the variance of
the random effect. From the definition in the help of ranef.lme, I
believe it is the random effect divided by its conditional standard
error. I don't know how to get the posterior variance of the random
effects in lme, but I do in lmer, so we can experiment a bit. It has
been a while since I have really used lme and I do not think there is an
extractor function to get these and I didn't see the variances in the
model object.

Let's work through an example to see what we get. Here is what I see

library(nlme)
data(Orthodont)
detach(package:nlme)

library(Matrix)
fm1 <- lmer(distance ~ age + (age|Subject), data = Orthodont)
# equivalent to
# fm1 <- lme(distance ~ age, data=Orthodont)

# Extract the variances of the random effects
qq <- attr(ranef(fm1, postVar = TRUE)[[1]], "postVar")

# divide the random effects by their standard error of "age"
Sranef_lmer <- ranef(fm1)[[1]][,2]/ sqrt(qq[2,2,])

library(nlme)
# Now, run the lme model
fm2 <- lme(distance ~ age, Orthodont)

# get the standardized random effects from lme
Sranef_lme <- ranef(fm2, standard=T)[2]

cor(Sranef_lmer, Sranef_lme)
[1] 1

Notice the perfect correlation. But, the actual values in Sranef_lme and
Sranef_lmer are a bit different and I cannot see why just yet. I need to
go eat lunch, but I'll think about this. 

Maybe somebody else sees something.