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

[Doran, Harold]

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 at stat.math.ethz.ch 
> [mailto:r-help-bounces at 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 at 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 
> 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.
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>