Mixed-Model

手動計算線性混合模型的隨機效應預測

  • October 27, 2016

我正在嘗試手動計算來自線性混合模型的隨機效應預測,並使用 Wood 在Generalized Additive Models: an Introduction with R (pg 294 / pg 307 of pdf) 中提供的符號,我對每個參數的含義感到困惑代表。

以下是伍德的摘要。

定義線性混合模型

哪裡 b N(0,), 和N(0,)

如果 b 和 y 是具有聯合正態分佈的隨機變量

RE預測由下式計算

在哪裡

使用lme4R 包中的隨機截距模型示例,我得到輸出

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 分鐘才發現第二個問題):

  1. 你不能計算使用殘差的平方,REML 將其估計為23.51,並且不能保證 BLUP 具有相同的方差。
sig <- 23.51

這不是!估計為39.19

psi <- 39.19

  1. 殘差不是用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 有關

因此.

引用自:https://stats.stackexchange.com/questions/242759

comments powered by Disqus