[R] standardized residuals (random effects) using nlme and ranef

[Dirk Enzmann]

As suggested I try another post.

First I give a reproducible example. The example data set has been 
provided by I. Plewis (1997), Statistics in Education. London: Arnold (I 
include the residuals obtained by MLWin):

# ---------------------------------------------

library(nlme)

# Example data from Plewis

BaseURL="http://www.ottersbek.de/"
ds32 = read.fwf(paste(BaseURL,'ds32.dat',sep=''),
                       widths=c(9,10,10,10,10,10),
                       col.names=c('class','pupil','zcurric',
                                   'curric','zmath1','math2'),
                       colClasses=c('factor','factor','numeric',
                                    'numeric','numeric','numeric'))

# Random slopes model (p. 44):
model.4b = lme(math2~zcurric,random=~zcurric|class,method='ML',data=ds32)

# Random effects (intercept and slope residuals) of level 1 (class):
RandEff = ranef(model.4b,level=1)

# Results of MLWin (c300 = intercept residuals of level "class",
#                   c301 = slope residuals of level "class",
#                   c304 = standardized intercept residuals,
#                   c305 = standardized slope residuals):
MLWinRes = read.fwf(paste(BaseURL,'resid_ex32.txt',sep=''),
                     widths=c(15,15,15,15),
                     col.names=c('c300','c301','c304','c305'),
                     colClasses=c('numeric','numeric','numeric','numeric'))

# They are identical to the results of MLWin (c300, c301):
round(cbind(RandEff,MLWinRes[,1:2]),5)

# However, using option "standard" the results differ considerably:
round(cbind(ranef(model.4b,level=1,standard=T),MLWinRes[,3:4]),5)

# ---------------------------------------------

 From the RSiteSearch("MLWin") there are two hits that seem to be 
relevant: One of Douglas Bates

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/53097.html

where he explains why he regards getting standard errors of the 
estimates of variance components as not being important. I think (but 
I'm not sure) that this implies that I should not use standardized 
residuals, as well.

Secondly a post of Roel de Jong

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/62280.html

However, I can't figure out how I could make use of his suggestion to 
obtain the standardized residuals I am looking for.

A part of the answer has been given in an answer to my previous post by 
Harold Doran. He showed that the so called "standardized residuals" 
obtained by ranef() using the option "standard=T" are the residuals 
devided by the standard deviation of the random effects, not divided by 
the "corresponding standard error" as stated in the ranef() help - if I 
understood him correctly.

In the multilevel list he also showed how to create a caterpillar plot 
using lmer() and the errbar() function of Hmisc (I modified his example 
to adapt it to my data):

# ------------------------------------

detach(package:nlme)
library(Matrix)
library(Hmisc)

# Fit a sample model:
fm1 = lmer(math2 ~ zcurric + (zcurric|class), ds32, method='ML', 
control=list(gradient=FALSE, niterEM=0))

.Call("mer_gradComp", fm1, PACKAGE = "Matrix")

# extract random effects:
fm1.blup = ranef(fm1)$class

#extract variances and multiple by scale factor:
fm1.post.var = fm1 at bVar$class * attr(VarCorr(fm1),"sc")**2

# posterior variance of slope on fm1:
fm1.post.slo = sqrt(fm1.post.var[2,2,])

# Slopes:
x = fm1.blup[,2]+outer(fm1.post.slo, c(-1.96,0,1.96))

# This code creates a catepillar plot using the errbar function:
x <- x[order(x[,2]) , ]
print(errbar(1:dim(x)[1], x[,2],x[,3],x[,1],ylab='Slopes', xlab='Classes'))
abline(h=0)

# ------------------------------------

Is it correct that I can obtain a respecive plot for the intercepts in a 
similar way?

# ------------------------------------

# posterior variance of intercept on fm1.post.int:
fm1.post.int = sqrt(fm1.post.var[1,1,])

# Intercepts:
x = fm1.blup[,1]+outer(fm1.post.int, c(-1.96,0,1.96))

# This code creates a catepillar plot using the errbar function
x <- x[order(x[,2]) , ]
print(errbar(1:dim(x)[1], x[,2],x[,3],x[,1],ylab='Intercepts', 
xlab='Classes'))
abline(h=0)

# ------------------------------------

To sum up, I can't figure out how MLWin calculates the standardized 
residuals. But I understand that this is not a question for the R list. 
Nevertheless, it would help if someone could point me to some arguments 
why not to use them and stick to the results obtainable by ranef().

Spencer Graves wrote:

>       Have you tried RSiteSearch("MLWin")?  I just got 29 hits.  I 
> wonder if any one of these might relate to your question?
> 
>       If you would like more help on this issue from this listserve, 
> please submit another post, preferably illustrating your question with 
> the simplest possible self-contained example that illustrates your 
> question, perhaps like the following:
> 
>       fm1.16 <- lme(distance~age, data=Orthodont[1:16,],
>            random=~age|Subject)
> 
>       Hope this helps.
>       Spencer Graves
> p.s.  PLEASE do read the posting guide 
> "www.R-project.org/posting-guide.html" and provide commented, minimal, 
> self-contained, reproducible code.
> 
> Dirk Enzmann wrote:
> 
>> Using ranef() (package nlme, version 3.1-75) with an 'lme' object I 
>> can obtain random effects for intercept and slope of a certain level 
>> (say: 1) - this corresponds to (say level 1) "residuals" in MLWin. 
>> Maybe I'm mistaken here, but the results are identical.
>>
>> However, if I try to get the standardized random effects adding the 
>> paramter "standard=T" to the specification of ranef(), the results 
>> differ considerably from the results of MLWin (although MLWin defines 
>> "standardized" in the same way as "divided by its estimated 
>> (diagnostic) standard error").
>>
>> 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)?
>>
>> Dirk
>>
>> -----------------------------
>> R version: 2.3.1 Patched (2006-06-21 r38367)

-- 
*************************************************
Dr. Dirk Enzmann
Institute of Criminal Sciences
Dept. of Criminology
Edmund-Siemers-Allee 1
D-20146 Hamburg
Germany

phone: +49-(0)40-42838.7498 (office)
        +49-(0)40-42838.4591 (Billon)
fax:   +49-(0)40-42838.2344
email: dirk.enzmann at uni-hamburg.de
www: 
http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Enzmann.html