How does Stata calculate the "predict varname, u" after xtreg random-effect? What "predict varname, u" after xtreg with random-effects really do in Stata? How it works? I mean, how the ("individual") random-error component u_i is extracted from the overall e_it error component?
Is there a way doing this in R?
 A: Let me expand my comment above. Assume your mixed-effects model is $$\boldsymbol y_i = \boldsymbol X_i \boldsymbol \beta + \boldsymbol Z_i \boldsymbol b_i + \boldsymbol \epsilon_i,$$ where $i$ denotes cluster $i$, random effects $\boldsymbol b_i \sim N(0,\boldsymbol D)$. Without loss of generality, assume $\boldsymbol \epsilon_i \sim N(0,\sigma^2 \boldsymbol I)$. Then the formula to calculate random effects is $$\hat{\boldsymbol b}_i=(\boldsymbol Z_i'\boldsymbol Z_i+\sigma^2 \boldsymbol D^{-1})^{-1}\boldsymbol Z_i'(\boldsymbol y_i - \boldsymbol X_i \boldsymbol \beta),$$ which can be generalized if you have a more general structure for $\mathrm{var}(\boldsymbol \epsilon_i).$ 
As far as I know, we can think the estimate $\hat{\boldsymbol b}_i$ in at least three ways:


*

*the best linear unbiased prediction (BLUP).

*the conditional mean/mode of $\boldsymbol b_i | \boldsymbol
   y_i$. If we use EM algorithm for the estimation, we will see the
above formula in the E step.

*the empirical Bayes estimator, since we do not specify a prior for $\boldsymbol b_i$ as in full Bayesian analysis. In some areas, it is also called "expected a posteriori".

