Oblicz ręcznie przewidywania efektów losowych dla liniowego modelu mieszanego

10

Staram się ręcznie wyliczyć przewidywania efektu losowego z liniowego modelu mieszanego i używając notacji dostarczonej przez Wooda w uogólnionych modelach addytywnych: wprowadzenie do R (str. 294 / str. 307 z pdf), mylę się co do każdego parametru reprezentuje.

Poniżej znajduje się podsumowanie z Wood.

Zdefiniuj liniowy model mieszany

Y=Xβ+Zb+ϵ

gdzie b N (0, ) i N (0, )ψ ϵ σ 2ψϵσ2)

Jeśli b i y są zmiennymi losowymi ze wspólnym rozkładem normalnym

[by]N.[[0Xβ],[ψΣbyΣybΣθσ2)]]

Prognozy RE są obliczane przez

mi[by]=ΣbyΣyy-1(y-xβ)=ΣbyΣθ-1(y-xβ)/σ2)=ψzT.Σθ-1(y-xβ)/σ2)

gdzieΣθ=ZψZT./σ2)+jan

Używając przykładowego modelu przechwytywania losowego z lme4pakietu R, otrzymuję dane wyjściowe

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

Na tej podstawie sądzę, że = 23,51, można oszacować na podstawie i z kwadratu reszt na poziomie populacji.ψ(y-Xβ)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)

Pomnożenie ich razem daje

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

co nie jest poprawne w porównaniu do

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

Dlaczego?

użytkownik2957945
źródło

Odpowiedzi:

9

Dwa problemy (przyznaję, że dostrzeżenie drugiego zajęło mi około 40 minut):

  1. Nie wolno obliczać z kwadratem reszt, jest on szacowany przez REML as i nie ma gwarancji, że BLUP będą miały tę samą wariancję.σ2)23.51

    sig <- 23.51

    A to nie jest ! Który ocenia się jakoψ39.19

    psi <- 39.19
  2. Pozostałości nie są uzyskiwane za pomocą, cake$angle - predict(m, re.form=NA)ale za pomocą residuals(m).

Składając to razem:

> 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

który jest podobny do ranef(m).

Naprawdę nie rozumiem, co predictoblicza.


ϵ^P.YP.=V.-1-V.-1X(XV.-1X)-1XV.-1

ϵ^=σ2)P.Y
b^=ψZtP.Y.

b^=ψ/σ2)Ztϵ^

Elvis
źródło
1
y-xβplot(residuals(m), cake$angle-predict(m, re.form=NULL)) ; plot(residuals(m), cake$angle-predict(m, re.form=NA))
1
Sposób, za pomocą stałego efektu, a trzecia wersja E [b | a] powyżej: z = getME(m, "Z") ; big_sig = solve(((z * psi) %*% zt ) / sig + diag(270)) ; psi/sig * zt %*% big_sig %*% (cake$angle-predict(m, re.form=NA)). Dziękujemy za wskazanie prawidłowych pozycji.
user2957945,
ΣbyΣyy
ΣybψZ
Elvis, myślałem o tym jeszcze raz (wiem, że jestem wolny). Myślę, że użycie takich reszt nie jest tak naprawdę sensowne, ponieważ wykorzystuje do obliczenia przewidywane wartości (a więc resztki) na poziomie RE, więc używamy go po obu stronach twojego równania. (więc wykorzystuje prognozy RE (E [b | y]) do prognozowania reszt, mimo że są to warunki, które staramy się przewidzieć))
2957945