Hitung prediksi efek acak secara manual untuk model campuran linier

Saya mencoba menghitung prediksi efek acak dari model campuran linier dengan tangan, dan menggunakan notasi yang diberikan oleh Wood dalam Generalized Additive Models: sebuah pengantar dengan R (pg 294 / pg 307 pdf), saya semakin bingung mengenai apa yang masing-masing parameter mewakili.

Di bawah ini adalah ringkasan dari Wood.

Tentukan model campuran linier

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

di mana b N (0, \ psi ), dan \ epsilon \ sim N (0, \ sigma ^ {2} )ψ ϵ ∼ σ $\sim$$\psi$$\epsilon \sim$$\sigma^{2}$

Jika b dan y adalah variabel acak dengan distribusi normal bersama

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

Prediksi RE dihitung oleh

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

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

Menggunakan contoh model intersep acak dari lme4paket R saya mendapatkan 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

Jadi dari sini, saya pikir = 23,51, dapat diperkirakan dari , dan dari kuadrat residual tingkat populasi.( y - X β $\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)

Mengalikan ini bersama memberi

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

yang tidak benar bila dibandingkan dengan

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

Mengapa?

Dua masalah (saya akui butuh waktu 40 menit untuk menemukan yang kedua):

1. Anda tidak boleh menghitung dengan kuadrat residual, diperkirakan oleh REML sebagai , dan tidak ada jaminan bahwa BLUP akan memiliki varian yang sama.$\sigma^2$23.51

   sig <- 23.51

   Dan ini bukan ! Yang diperkirakan sebagai$\psi$39.19

   psi <- 39.19

2. Sisa tidak diperoleh dengan cake$angle - predict(m, re.form=NA)tetapi dengan residuals(m).

Menyatukannya:

> 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

yang mirip dengan ranef(m).

Saya benar-benar tidak mendapatkan apa yang predictdihitung.

---

PS. Untuk menjawab komentar terakhir Anda, intinya adalah bahwa kami menggunakan "residual" sebagai cara untuk mendapatkan vektor mana . Matriks ini dihitung selama algoritma REML. Ini terkait dengan BLUP istilah acak oleh dan 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.$$

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