Calculate random effect predictions manually for a linear mixed model I am trying to calculate the random effect predictions from a linear mixed model by hand, and
using notation provided by Wood in Generalized Additive Models: an introduction with R (pg 294 / pg 307 of pdf), 
I am getting confused over what each parameters represents.
Below is a summary from Wood.
Define a linear mixed model
$$ Y = X\beta + Zb + \epsilon $$
where b  $\sim$ N(0, $\psi$), and $\epsilon \sim$ N(0, $\sigma^{2}$)
If b and y are random variables with joint normal distribution
\begin{align*}
  \begin{bmatrix}
    b\\
    y
  \end{bmatrix} &\sim  N
\begin{bmatrix}
  \begin{bmatrix}
    0\\
    X\beta
  \end{bmatrix}\!\!,&
  \begin{bmatrix}
    \psi   & \Sigma_{by} \\
    \Sigma_{yb}& \Sigma_{\theta}\sigma^{2} 
  \end{bmatrix}
\end{bmatrix}\\
 \end{align*}
The RE predictions are calculated by
\begin{align*}
E[b \mid y] &= \Sigma_{by} \Sigma_{yy}^{-1} (y - x\beta)\\ 
            &= \Sigma_{by}\Sigma_{\theta}^{-1}(y - x \beta) / \sigma^{2}\\
            &= \psi z^{T}\Sigma_{\theta}^{-1} (y - x \beta) / \sigma^{2}
\end{align*}
where $\Sigma_{\theta} = Z\psi Z^{T} /\sigma^{2} + I_{n}$
Using an random intercept model example from lme4 R package I get output
library(lme4)
m = lmer(angle ~ temp + (1 | replicate), data=cake)
summary(m)

% Linear mixed model fit by REML ['lmerMod']
% Formula: angle ~ temp + (1 | replicate)
%    Data: cake
% 
% REML criterion at convergence: 1671.7
% 
% Scaled residuals: 
%      Min       1Q   Median       3Q      Max 
% -2.83605 -0.56741 -0.02306  0.54519  2.95841 
% 
% Random effects:
%  Groups    Name        Variance Std.Dev.
%  replicate (Intercept) 39.19    6.260   
%  Residual              23.51    4.849   
% Number of obs: 270, groups:  replicate, 15
% 
% Fixed effects:
%             Estimate Std. Error t value
% (Intercept)  0.51587    3.82650   0.135
% temp         0.15803    0.01728   9.146
% 
% Correlation of Fixed Effects:
%      (Intr)
% temp -0.903

So from this, I think $\psi$ = 23.51, $(y-X\beta)$ can be estimated from cake$angle - predict(m, re.form=NA), 
and sigma from the square of the population level residuals.
th = 23.51
zt = getME(m, "Zt") 
res = cake$angle - predict(m, re.form=NA)
sig = sum(res^2) / (length(res)-1)

Multiplying these together gives
th * zt %*% res / sig
         [,1]
1  103.524878
2   94.532914
3   33.934892
4    8.131864
---

which is not correct when compared to
> ranef(m)
$replicate
   (Intercept)
1   14.2365633
2   13.0000038
3    4.6666680
4    1.1182799
---

Why?
 A: Two problems (I confess it took me like 40 minutes to spot the second one):


*

*You must not compute $\sigma^2$ with the square of residuals, it is estimated by REML as 23.51, and there is no guarantee that the BLUPs will have the same variance.
sig <- 23.51

And this is not $\psi$ ! Which is estimated as 39.19
psi <- 39.19


*The residuals are not obtained with cake$angle - predict(m, re.form=NA) but with residuals(m).
Putting it together:
> psi/sig * zt %*% residuals(m)
15 x 1 Matrix of class "dgeMatrix"
         [,1]
1  14.2388572
2  13.0020985
3   4.6674200
4   1.1184601
5   0.2581062
6  -3.2908537
7  -4.6351567
8  -4.5813846
9  -4.6351567
10 -3.1833095
11 -2.1616392
12 -1.1399689
13 -0.2258429
14 -4.0974355
15 -5.3341942

which is similar to ranef(m).
I really don't get what predict computes.

PS. To answer your last remark, the point is that we use the "residuals" $\hat\epsilon$ as a way to obtain the vector $PY$ where $P = V^{-1} -  V^{-1}  X \left(X' V^{-1} X\right)^{-1} X' V^{-1}$. This matrix is computed during the REML algorithm. It is related to the BLUPs of random terms by
$$ \hat\epsilon = \sigma^2 PY $$
and
$$ \hat b = \psi Z^t PY. $$
Thus $ \hat b = \psi / \sigma^2 Z^t \hat\epsilon $.
