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

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)[])

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

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

• The same question is asked here, but a different answer is given. – filups21 May 21 '14 at 15:53

VarCorr() here provides the variance estimate of the random intercept by REML, but sd(ranef(fit1)[]) 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().