Рассчитать прогноз случайного эффекта вручную для линейной смешанной модели

Я пытаюсь вручную вычислить предсказания случайных эффектов на основе линейной смешанной модели и, используя обозначения, представленные Вудом в Обобщенных аддитивных моделях: введение в R (стр. 294 / стр. 307 в формате pdf), я запутываюсь в том, что каждый из параметров представляет.

Ниже приведено резюме из дерева.

Определить линейную смешанную модель

$$Y = X\beta + Zb + \epsilon$$

где b N (0, ) и N (0, )ψ ϵ ∼ σ $\sim$$\psi$$\epsilon \sim$$\sigma^{2}$

Если b и y - случайные величины с совместным нормальным распределением

$$\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*}$$

Прогнозы RE рассчитываются по

$$\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*}$$

где$\Sigma_{\theta} = Z\psi Z^{T} /\sigma^{2} + I_{n}$

Используя пример модели случайного перехвата из lme4пакета R, я получаю вывод

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

Так что с этим, я думаю = 23,51, можно оценить , и от площади остатков уровня населения.( у - Х β $\psi$$(y-X\beta)$cake$angle - predict(m, re.form=NA)sigma

th = 23.51
zt = getME(m, "Zt") 
res = cake$angle - predict(m, re.form=NA)
sig = sum(res^2) / (length(res)-1)

Умножение их вместе дает

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

что не правильно по сравнению с

> ranef(m)
$replicate
   (Intercept)
1   14.2365633
2   13.0000038
3    4.6666680
4    1.1182799
---

Почему?

Две проблемы (признаюсь, мне понадобилось около 40 минут, чтобы найти вторую):

1. Вы не должны вычислять с квадратом невязок, это оценивается REML как , и нет никакой гарантии, что BLUP будут иметь такую ​​же дисперсию.$\sigma^2$23.51

   sig <- 23.51

   И это не ! Который оценивается как$\psi$39.19

   psi <- 39.19

2. Остатки получены не с, cake$angle - predict(m, re.form=NA)а с residuals(m).

Собираем это вместе:

> 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

который похож на ranef(m).

Я действительно не понимаю, что predictвычисляет.

---

PS. Чтобы ответить на ваше последнее замечание, дело в том, что мы используем «остатки» как способ получения вектора где . Эта матрица вычисляется во время алгоритма REML. Он связан с BLUP случайных слагаемых: и РУР=V-1-V-1Х(Х'V - 1 Х)-1Х'V-1 ε =σ2РУ Ь =ψZтPY$\hat\epsilon$$PY$$P = V^{-1} - V^{-1} X \left(X' V^{-1} X\right)^{-1} X' V^{-1}$

$$\hat\epsilon = \sigma^2 PY$$ $$\hat b = \psi Z^t PY.$$

Таким образом, .$\hat b = \psi / \sigma^2 Z^t \hat\epsilon$