W lmer
funkcji w lme4
w R
istnieje wezwanie do skonstruowania matrycy modelu efektów przypadkowych, , jak opisano tutaj , na stronach 7 - 9.
Obliczanie obejmuje produkty KhatriRao i / lub Kronecker dwóch matryc, i . J i X i
Macierz to kęs: „Macierz wskaźników wskaźników współczynnika grupowania”, ale wydaje się, że jest to rzadka macierz z fałszywym kodowaniem do wyboru, która jednostka (na przykład podmioty w powtarzalnych pomiarach) odpowiadająca wyższym poziomom hierarchicznym jest „włączona” dla każda obserwacja. matrycy wydaje się działać jako selektor pomiarów na niższy poziom hierarchiczny, tak że „połączenie obu przełączników” dawałyby macierzy formy przedstawionej na papier za pomocą następującego przykładu:X i Z i
(f<-gl(3,2))
[1] 1 1 2 2 3 3
Levels: 1 2 3
(Ji<-t(as(f,Class="sparseMatrix")))
6 x 3 sparse Matrix of class "dgCMatrix"
1 2 3
[1,] 1 . .
[2,] 1 . .
[3,] . 1 .
[4,] . 1 .
[5,] . . 1
[6,] . . 1
(Xi<-cbind(1,rep.int(c(-1,1),3L)))
[,1] [,2]
[1,] 1 -1
[2,] 1 1
[3,] 1 -1
[4,] 1 1
[5,] 1 -1
[6,] 1 1
Transponowanie każdej z tych macierzy i wykonanie mnożenia Khatri-Rao:
Ale to transpozycja:
(Zi<-t(KhatriRao(t(Ji),t(Xi))))
6 x 6 sparse Matrix of class "dgCMatrix"
[1,] 1 -1 . . . .
[2,] 1 1 . . . .
[3,] . . 1 -1 . .
[4,] . . 1 1 . .
[5,] . . . . 1 -1
[6,] . . . . 1 1
Okazuje się, że autorzy korzystają z bazy danych sleepstudy
w lme4
, ale tak naprawdę nie rozwijają matryc projektowych, ponieważ odnoszą się one do tego konkretnego badania. Staram się więc zrozumieć, w jaki sposób skompilowany kod w powyższym artykule przełożyłby się na bardziej znaczący sleepstudy
przykład.
Dla uproszczenia wizualnego ograniczyłem zestaw danych do zaledwie trzech tematów - „309”, „330” i „371”:
require(lme4)
sleepstudy <- sleepstudy[sleepstudy$Subject %in% c(309, 330, 371), ]
rownames(sleepstudy) <- NULL
Każda osoba wykazywałaby zupełnie inne przechwytywanie i nachylenie, gdyby prosta regresja OLS była rozpatrywana indywidualnie, co sugeruje potrzebę modelu mieszanego z wyższą hierarchią lub poziomem jednostek odpowiadającym badanym:
par(bg = 'peachpuff')
plot(1,type="n", xlim=c(0, 12), ylim=c(200, 360),
xlab='Days', ylab='Reaction')
for (i in sleepstudy$Subject){
fit<-lm(Reaction ~ Days, sleepstudy[sleepstudy$Subject==i,])
lines(predict(fit), col=i, lwd=3)
text(x=11, y=predict(fit, data.frame(Days=9)), cex=0.6,labels=i)
}
Wywołanie regresji z efektem mieszanym to:
fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy)
A macierz wyodrębniona z funkcji daje następujące:
parsedFormula<-lFormula(formula= Reaction~Days+(Days|Subject),data= sleepstudy)
parsedFormula$reTrms
$Ztlist
$Ztlist$`Days | Subject`
6 x 12 sparse Matrix of class "dgCMatrix"
309 1 1 1 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . .
309 0 1 2 3 4 5 6 7 8 9 . . . . . . . . . . . . . . . . . . . .
330 . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 . . . . . . . . . .
330 . . . . . . . . . . 0 1 2 3 4 5 6 7 8 9 . . . . . . . . . .
371 . . . . . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1
371 . . . . . . . . . . . . . . . . . . . . 0 1 2 3 4 5 6 7 8 9
Wydaje się to słuszne, ale jeśli tak, to czym jest algebra liniowa? Rozumiem, że wiersze 1
to wybór osób takich jak. Na przykład, temat 309
jest włączony dla linii bazowej + dziewięć obserwacji, więc dostaje cztery 1
i tak dalej. Druga część jest wyraźnie faktycznym pomiarem: 0
dla wartości początkowej, 1
pierwszego dnia pozbawienia snu itp.
Ale jakie są rzeczywiste i oraz odpowiadające im lub , cokolwiek jest istotny?X i Z i = ( J T i ∗ X T i ) ⊤ Z i = (
Oto możliwość
Problem polega na tym, że nie jest to transponowana, jak się lmer
wydaje, funkcja wymaga i nadal nie jest jasne, jakie reguły mają tworzyć .
źródło
mkZt()
(poszukaj tutaj ) to dobry początek?Odpowiedzi:
309
330
371
nrow(sleepstudy[sleepstudy$Subject==309,]) [1] 10
f <- gl(3,10)
Ji<-t(as(f,Class="sparseMatrix"))
Budowanie macierzy można wspomóc za pomocą funkcji jako odniesienia:Xja
getME
library(lme4) sleepstudy <- sleepstudy[sleepstudy$Subject %in% c(309, 330, 371), ] rownames(sleepstudy) <- NULL fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy)
Xi <- getME(fm1,"mmList")
Ponieważ potrzebujemy transpozycji, a obiekt
Xi
nie jest macierzą,t(Xi)
można go zbudować jako:t_Xi <- rbind(c(rep(1,30)),c(rep(0:9,3)))
Zi<-t(KhatriRao(t_Ji,t_Xi))
:Odpowiada to równaniu (6) w oryginalnej pracy :
Aby to zobaczyć, możemy zamiast tego bawić się obciętymi i , wyobrażając sobie, że zamiast 9 pomiarów i linii bazowej (0), jest tylko 1 pomiar (i linia bazowa). Otrzymane macierze to:jotT.ja XT.ja
X T i = [ 1 1 1 1 1 1 0 1 0 1 0 1 ]jotT.ja= [ 100100010010001001] i .XT.ja= [ 101110111011]
I
b <- getME(fm1,"b")
Jeśli dodamy te wartości do stałych efektów wywołania
fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy)
, otrzymamy przechwyty:i stoki:
wartości zgodne z:
as.matrix(Zi)%*%b
.źródło