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

From: Doran, Harold <HDoran_at_air.org>
Date: Mon 31 Jul 2006 - 10:29:52 EST

OK, I see how the standardized random effects are calculated. Here is what I now see

library(nlme)

fm2 <- lme(distance~age, Orthodont)

# unstandardized
age_ranef <- ranef(fm2)[,2]

#standardized
age_Sranef <- ranef(fm2, standard=TRUE)[,2]

# We can use these to solve for the standard error, because the formula according to help for ranef.lme is

# Standardized_randomEffects = random_effects/standard error

age_ranef/age_Sranef

# OK, now note the values are exactly the same. Now, look at

VarCorr(fm2)

You can see the value used to standardize is the standard deviation of the random effect for age.

Now, the help function does say "divided by the corresponding standard error".

I've copied Doug Bates because the values in the stdDev column are the standard deviations of the variance components and not standard errors of those variance components. So, I'm not sure why the help says that the standardized random effects are divided by the corresponding SE. Maybe he can clarify if he has time.

I hope that helps
Harold

> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Doran, Harold
> Sent: Sunday, July 30, 2006 3:40 PM
> To: Spencer Graves; Dirk Enzmann
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] standardized random effects with ranef.lme()
>
>
> > > 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
>
> 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)[], "postVar")
>
> # divide the random effects by their standard error of "age"
> Sranef_lmer <- ranef(fm1)[][,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)
>
> cor(Sranef_lmer, Sranef_lme)
>  1
>
> Notice the perfect correlation. But, the actual values in
> Sranef_lme and Sranef_lmer are a bit different and I cannot
>
> Maybe somebody else sees something.
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help