Mixed-Model
手動計算線性混合模型的隨機效應預測
我正在嘗試手動計算來自線性混合模型的隨機效應預測,並使用 Wood 在Generalized Additive Models: an Introduction with R (pg 294 / pg 307 of pdf) 中提供的符號,我對每個參數的含義感到困惑代表。
以下是伍德的摘要。
定義線性混合模型
哪裡 b N(0,), 和N(0,)
如果 b 和 y 是具有聯合正態分佈的隨機變量
RE預測由下式計算
在哪裡
使用
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,可以從
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 分鐘才發現第二個問題):
- 你不能計算使用殘差的平方,REML 將其估計為
23.51
,並且不能保證 BLUP 具有相同的方差。sig <- 23.51
這不是!估計為
39.19
psi <- 39.19
- 殘差不是用
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 有關
和
因此.