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
gdzie b N (0, ) i N (0, )ψ ϵ ∼ σ 2
Jeśli b i y są zmiennymi losowymi ze wspólnym rozkładem normalnym
Prognozy RE są obliczane przez
gdzie
Używając przykładowego modelu przechwytywania losowego z lme4
pakietu 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.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?
źródło
plot(residuals(m), cake$angle-predict(m, re.form=NULL)) ; plot(residuals(m), cake$angle-predict(m, re.form=NA))
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.