Berechnen Sie zufällige Effektvorhersagen manuell für ein lineares gemischtes Modell

Ich versuche, die Zufallseffektvorhersagen aus einem linearen gemischten Modell von Hand zu berechnen, und unter Verwendung der von Wood in Generalized Additive Models bereitgestellten Notation : eine Einführung mit R (S. 294 / S. 307 von PDF), bin ich verwirrt über die einzelnen Parameter repräsentiert.

Unten finden Sie eine Zusammenfassung von Wood.

Definieren Sie ein lineares gemischtes Modell

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

wobei b N (0, ) und N (0, )$\sim$$\psi$$\epsilon \sim$$\sigma^{2}$

Wenn b und y Zufallsvariablen mit gemeinsamer Normalverteilung sind

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

Die RE-Vorhersagen werden von berechnet

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

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

Anhand eines zufälligen Intercept-Modellbeispiels aus dem lme4R-Paket erhalte ich eine Ausgabe

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

Daher denke ich, dass = 23,51 aus und aus dem Quadrat der Residuen auf Bevölkerungsebene geschätzt werden kann .$\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)

Multiplizieren Sie diese zusammen

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

was im Vergleich zu nicht korrekt ist

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

Warum?

Zwei Probleme (ich gestehe, ich habe ungefähr 40 Minuten gebraucht, um das zweite zu erkennen):

1. Sie dürfen mit dem Quadrat der Residuen berechnen , es wird von REML als geschätzt , und es gibt keine Garantie dafür, dass die BLUPs dieselbe Varianz haben.$\sigma^2$23.51

   sig <- 23.51

   Und das ist nicht ! Welches wird geschätzt als$\psi$39.19

   psi <- 39.19

2. Die Residuen werden nicht mit, cake$angle - predict(m, re.form=NA)sondern mit erhalten residuals(m).

Etwas zusammensetzen:

> 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

das ist ähnlich wie ranef(m).

Ich verstehe wirklich nicht, was predictberechnet wird.

---

PS. Um Ihre letzte Bemerkung zu beantworten, ist der Punkt, dass wir die "Residuen" , um den Vektor wobei . Diese Matrix wird während des REML-Algorithmus berechnet. Es ist verwandt mit den BLUPs zufälliger Terme durch und PYP=V-1-V-1X(X'V - 1 X)-1X'V-1 ε =σ2PY b =ψZtPY$\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.$$

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