Why computing the variance of the extracted random effect (using ranef) is not the same as the output from lme?

Question

I am applying a mixed model to predict tumor progression (y) using tumor volume (x) as the fixed effect and center ($i=1,...10$) as random intercept. The model can be written as: $$y_{ij}=\alpha+\beta x_{ij}+b_{i}+\epsilon_{ij}$$

I used the lme() function in R:

fit1 <- lme(PD ~ log(Volume), random = ~1 | CenterID, data=Data)

The result shows that the standard deviation of the random intercept is 0.079. Thus $b_{i}$ follows a normal distribution $N(0,0.079)$.

In the meantime, I can extract the random intercept by applying ranef(fit1). This gives a list of $b_{i}$ corresponding to each center. Then I compute the standard of this vector.

sd(ranef(fit1)[[1]])

I would expect that it gives similar result as 0.079. However, it is far different.

Can someone tell me why sd(ranef(fit1)[[1]]) gives different result than the model output VarCorr(fit1)? What is exactly the relation between ranef(fit1) and VarCorr(fit1)?

Answer 1

VarCorr() here provides the variance estimate of the random intercept by REML, but sd(ranef(fit1)[[1]]) shows the standard deviation of the empirical Bayesian estimates of the random intercept. In other words, say we have parameter $\sigma$, VarCorr() gives us the estimate $\hat \sigma$; we can generate a random sample based on $\hat \sigma$, and then calculate the standard deviation of the random sample by sd().

The two would be different since the result by sd() is an estimate of the result by VarCorr(). The similarity depends on the data set, and the two are similar when tested with the Orthodont example in lme().

> summary(fit1 <- lme(distance ~ age, random = ~1 | Subject, data=Orthodont))
Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    2.114724 1.431592

> sd(ranef(fit1)[[1]])
[1] 2.003087