Matryce modelowe dla modeli efektów mieszanych

W lmerfunkcji w lme4w Ristnieje wezwanie do skonstruowania matrycy modelu efektów przypadkowych, $Z$ , jak opisano tutaj , na stronach 7 - 9.

Obliczanie obejmuje produkty KhatriRao i / lub Kronecker dwóch matryc, i . $Z$ $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: $J_i$ $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:

$\left[\begin{smallmatrix}1 & 1 &. &. &. &.\\.&.&1&1&.&.\\.&.&.&.&1&1 \end{smallmatrix}\right]\ast \left[\begin{smallmatrix}\,\,\,\,1 & 1 &\,\,\,\,1 &1 &\,\,\,\,1 &1\\-1&1&-1&1&-1&1 \end{smallmatrix}\right]= \left[\begin{smallmatrix}\,\,1 & 1 &.&.&.&.\\\,\,\,\,-1 &1&.&.&.&.\\ .&.&\,\,\,\,\,1 &1&.&.\\.&.&\,\,-1&1&.&.\\.&.&.&.&\,\,\,1&1\\.&.&.&.&-1&1 \end{smallmatrix}\right]$

Ale to transpozycja: $Z_i$

(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 sleepstudyw 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 sleepstudyprzykł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)
        }

wprowadź opis zdjęcia tutaj

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 1to wybór osób takich jak. Na przykład, temat 309jest włączony dla linii bazowej + dziewięć obserwacji, więc dostaje cztery 1i tak dalej. Druga część jest wyraźnie faktycznym pomiarem: 0dla wartości początkowej, 1pierwszego dnia pozbawienia snu itp.

Ale jakie są rzeczywiste i oraz odpowiadające im lub , cokolwiek jest istotny? $J_i$ $X_i$ $Z_i= (J_i^{T}∗X_i^{T})^⊤$ $Z_i= (J_i^{T}\otimes X_i^{T})^⊤$

Oto możliwość

$\left[\begin{smallmatrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 &1&1&. &. &. &. &.& . &. &. &. &. &. &.&. & . &. &. &. &. &. &.\\ .&.&.&.&. & . &. &. &. &. &1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 &1&1&.& . &. &. &. &. &. &.&.&.\\&. & . &. &. &. &. &. &.&. & . &. &. &. &. &.&.&.&.&.&1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 &1&1\end{smallmatrix}\right]\ast \left[\begin{smallmatrix}1 & 1 & 1 & 1&1&1&1 & 1 & 1 & 1\\0&1&2&3&4&5&6&7&8&9 \end{smallmatrix}\right]=$

$\left[\begin{smallmatrix}1 & 1 & 1 &1&1&1&1&1&1&1&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.\\0 & 1 & 2 &3&4&5&6&7&8&9&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.\\.&.&.&.&&.&.&.&.&.&1 &1&1&1&1&1&1&1&1&1&.&.&.&.&.&.&.&.&.&.\\ &.&.&.&.&.&.&.&.&.&0 &1&2&3&4&5&6&7&8&9&.&.&.&.&.&.&.&.&.&.\\.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&1&1&1&1&1&1&1&1&1&1\\&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&.&0&1&2&3&4&5&6&7&8&9 \end{smallmatrix}\right]$

Problem polega na tym, że nie jest to transponowana, jak się lmerwydaje, funkcja wymaga i nadal nie jest jasne, jakie reguły mają tworzyć . $X_i$

r mixed-model lme4-nlme matrix

— Antoni Parellada
źródło

Jest to o wiele łatwiejsze niż się wydaje. matrycy tutaj jest po prostu (transpozycją) produkt Kronecker z matrycy identyczny z matrycy projektu.

Z

$Z$

— Donnie

Dziękuję za podpowiedź. Nadal będę pracować nad zrozumieniem wszystkich sub-indków w liniowym szkielecie algebry tej funkcji. Jeśli zatrzaśnie się na miejscu, odpowiem na własne pytanie, ale choć wiem, że jest to bardzo proste, zgodność między nomenklaturą rusztowań matematycznych a stosowaniem do jakiegokolwiek przykładu jest myląca.

— Antoni Parellada,

Innym dobrym zasobem dla Ciebie może być ich implementacja lme4pureR , która następuje wraz z powyższą winietą i jest napisana w całości w R. Może mkZt()(poszukaj tutaj ) to dobry początek?

— alexforrence

Tworzenie matrycy wymaga wytwarzania 3 (poziom , i ) każdy z 10 pomiarów i obserwacji ( ). Zgodnie z kodem w oryginalnym linku w PO: $J_i$ 309330371nrow(sleepstudy[sleepstudy$Subject==309,]) [1] 10

f <- gl(3,10) Ji<-t(as(f,Class="sparseMatrix"))

wprowadź opis zdjęcia tutaj

Budowanie macierzy można wspomóc za pomocą funkcji jako odniesienia: $X_i$ 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")

wprowadź opis zdjęcia tutaj

Ponieważ potrzebujemy transpozycji, a obiekt Xinie jest macierzą, t(Xi)można go zbudować jako:

t_Xi <- rbind(c(rep(1,30)),c(rep(0:9,3)))

$Z_i$ oblicza się jako : $Z_i= (J_i^{T}∗X_i^{T})^⊤$

Zi<-t(KhatriRao(t_Ji,t_Xi)):

wprowadź opis zdjęcia tutaj

Odpowiada to równaniu (6) w oryginalnej pracy :

Z_{ja} = ({jot}_{ja}^{T.} * X_{ja}^{T.})^{T.} = [\begin{matrix} {jot}_{ja 1}^{T.} \otimes X_{ja 1}^{T.} \\ {jot}_{ja 2)}^{T.} \otimes X_{ja 2)}^{T.} \\ ⋮ \\ {jot}_{ja n}^{T.} \otimes X_{ja n}^{T.} \end{matrix}]

$Z_i = (J_i^T*X_i^T)^T=\begin{bmatrix}J_{i1}^T\otimes X_{i1}^T\\J_{i2}^T\otimes X_{i2}^T\\\vdots\\J_{in}^T\otimes X_{in}^T\end{bmatrix}$

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: $J_i^T$ $X_i^T$

$J_i^T=\left[\begin{smallmatrix}1&1&0&0&0&0\\0&0&1&1&0&0\\0&0&0&0&1&1\end{smallmatrix}\right]$ i . $X_i^T=\left[\begin{smallmatrix}1&1&1&1&1&1\\0&1&0&1&0&1\end{smallmatrix}\right]$

$J_i^T*X_i^T=\left[ \begin{smallmatrix} \left(\begin{smallmatrix}1\\0\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\0\end{smallmatrix}\right)&& \left(\begin{smallmatrix}1\\0\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\1\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\1\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\0\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\1\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\1\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\0\\1\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\0\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\0\\1\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\1\end{smallmatrix}\right) \end{smallmatrix} \right]$

$\small= \begin{bmatrix}J_{i1}^T\otimes X_{i1}^T && J_{i2}^T \otimes X_{i2}^T && J_{i3}^T\otimes X_{i3}^T && J_{i4}^T\otimes X_{i4}^T && J_{i5}^T\otimes X_{i5}^T && J_{i6}^T\otimes X_{i6}^T\end{bmatrix}$

$=\left[\begin{smallmatrix}1&1&0&0&0&0\\0&1&0&0&0&0\\0&0&1&1&0&0\\0&0&0&1&0&0\\0&0&0&0&1&1\\0&0&0&0&0&1\end{smallmatrix}\right]$ . Które transponowane i rozszerzone spowodują, że . $Z_i=\left[\begin{smallmatrix}1&0&0&0&0&0\\1&1&0&0&0&0\\1&2&0&0&0&0\\\vdots\\\\0&0&1&0&0&0\\0&0&1&1&0&0\\0&0&1&2&0&0\\\vdots\\\\0&0&0&0&1&0\\0&0&0&0&1&1\\0&0&0&0&1&2\\\\\vdots\end{smallmatrix}\right]$

Wyodrębnienie wektora współczynników efektów losowych można wykonać za pomocą funkcji: $b$

b <- getME(fm1,"b")

[1,] -44.1573839
[2,]  -2.4118590
[3,]  32.8633489
[4,]  -0.3998801
[5,]  11.2940350
[6,]   2.8117392

Jeśli dodamy te wartości do stałych efektów wywołania fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy), otrzymamy przechwyty:

205.3016 for 309; 282.3223 for 330; and 260.7530 for 371

i stoki:

2.407141 for 309; 4.419120 for 330; and 7.630739 for 371

wartości zgodne z:

library(lattice)
xyplot(Reaction ~ Days | Subject, groups = Subject, data = sleepstudy, 
       pch=19, lwd=2, type=c('p','r'))

wprowadź opis zdjęcia tutaj

$Zb$ można obliczyć jako as.matrix(Zi)%*%b.

— Antoni Parellada
źródło