Regresja wieloraka mieszana (liniowy model mieszany)

Kiedy stosować:

Regresja wieloraka mieszana (liniowy model mieszany, ang. linear mixed model, LMM) jest rozszerzeniem liniowej regresji wielorakiej na dane, w których obserwacje nie są niezależne, lecz tworzą grupy: wielokrotne pomiary tego samego pacjenta, pacjenci leczeni w tym samym ośrodku, uczniowie z jednej klasy, zwierzęta z jednego miotu (zob. wprowadzenie do modeli mieszanych). Obserwacje z tej samej grupy są do siebie bardziej podobne niż obserwacje z różnych grup. Zwykła regresja tę zależność ignoruje, przez co zaniża błędy standardowe współczynników i zawyża istotność wyników. Model mieszany rozbija zmienność zmiennej zależnej na część między grupami i część wewnątrz grup, a każdy współczynnik otrzymuje właściwy błąd.

Model zawiera dwa rodzaje efektów:

• efekty stałe – zwykłe współczynniki regresji \(\beta_0,\beta_1,\ldots,\beta_k\) opisujące populację,
• efekty losowe – odchylenia poszczególnych grup od populacji: losowy wyraz wolny (każda grupa ma własny poziom \(Y\)) i ewentualnie losowe nachylenia (efekt wybranej zmiennej \(X\) różni się między grupami).

Model z losowym wyrazem wolnym i losowym nachyleniem zmiennej \(X_1\) ma postać:

gdzie:
\(Y_{ij}\) – \(i\)-ta obserwacja w \(j\)-tej grupie,
\(\beta_0,\beta_1,\ldots,\beta_k\) – efekty stałe,
\(u_{0j}\) – losowy wyraz wolny grupy \(j\), \(u_{1j}\) – losowe nachylenie grupy \(j\); zakłada się, że \((u_{0j},u_{1j})\) mają rozkład normalny o średniej 0 i macierzy kowariancji \(\Sigma\) (wariancje \(\tau_0^2\), \(\tau_1^2\) i korelacja \(\rho_{01}\)),
\(\epsilon_{ij}\) – reszta, o rozkładzie normalnym ze średnią 0 i wariancją \(\sigma^2\).

Dla grupy \(j\) model opisuje więc własną prostą: jej wyraz wolny to \(\beta_0+u_{0j}\), a nachylenie \(\beta_1+u_{1j}\). Efekty stałe opisują prostą „przeciętnej” grupy, a wariancje efektów losowych – jak bardzo proste poszczególnych grup są rozrzucone wokół niej.

Model ten zastępuje ANOVA dla grup zależnych, gdy pomiary powtarzane są niekompletne lub nierówno rozmieszczone w czasie, gdy liczba pomiarów różni się między osobami, albo gdy oprócz czynnika powtarzanego w modelu mają się znaleźć zmienne ciągłe (np. wiek, dawka) lub zmienne opisujące grupę (np. płeć, ośrodek).

Warunki stosowania

Definicje miar

Wariancje efektów losowych – \(\tau_0^2\) (losowy wyraz wolny) mówi, jak bardzo grupy różnią się poziomem \(Y\); \(\tau_1^2\) (losowe nachylenie) – jak bardzo różnią się efektem zmiennej \(X_1\). Wygodniej czyta się ich pierwiastki, czyli odchylenia standardowe, bo są wyrażone w jednostkach \(Y\) (lub \(Y\) na jednostkę \(X_1\)): typowa grupa odbiega od populacji o mniej więcej jedno odchylenie standardowe, a około 95% grup mieści się w przedziale \(\pm2\) odchylenia.

Korelacja efektów losowych \(\rho_{01}\) – związek między poziomem grupy a jej nachyleniem. Korelacja ujemna oznacza, że grupy startujące wyżej zmieniają się wolniej. Uwaga: jej wartość zależy od tego, gdzie leży zero zmiennej \(X_1\) – wyraz wolny to poziom przy \(X_1=0\). Gdy \(X_1\) to np. wiek 8–14 lat, warto go wycentrować (odjąć wartość środkową), aby korelacja opisywała poziom w środku zakresu, a nie hipotetyczny poziom w wieku 0.

Współczynnik korelacji wewnątrzklasowej (ICC) – liczony z modelu pustego (bez zmiennych niezależnych, tylko losowy wyraz wolny):

gdzie \(\tau^2\) to wariancja między grupami, a \(\sigma^2\) wariancja wewnątrz grup (reszt). ICC to udział zmienności między grupami w zmienności całkowitej, a zarazem korelacja dwóch obserwacji z tej samej grupy. Przedział ufności ICC jest orientacyjny (metoda delta).

\(R^2\) brzegowe i warunkowe (Nakagawa i Schielzeth 2013[191]) – \(R^2\) brzegowe to udział wariancji wyjaśnianej przez same efekty stałe, \(R^2\) warunkowe – przez efekty stałe i losowe razem. Różnica między nimi mówi, ile zmienności „tłumaczy” przynależność do grupy.

Test LR efektów losowych – test ilorazu wiarygodności porównujący model z danym efektem losowym i bez niego. Hipoteza zerowa (wariancja równa 0) leży na granicy przestrzeni parametrów, więc wartość \(p\) z rozkładu \(\chi^2\) jest korygowana (mieszanka rozkładów; dla jednego parametru \(p\) dzieli się przez 2). Program zawsze liczy ten test metodą ML, niezależnie od metody estymacji wybranej dla modelu głównego.

Efekty losowe grup (BLUP) – przewidywane wartości \(u_{0j}\), \(u_{1j}\) dla każdej grupy (ang. best linear unbiased predictor). Są to oszacowania „ściągnięte” w stronę zera: grupa z małą liczbą obserwacji lub o dużym rozrzucie dostaje efekt bliższy 0, niż wynikałoby ze zwykłej średniej jej obserwacji. Służą do wskazania grup nietypowych i do predykcji dla konkretnej grupy, nie do testowania istotności.

Dopasowanie osobliwe (ang. singular fit) – sytuacja, gdy oszacowana wariancja któregoś efektu losowego jest praktycznie równa 0 albo korelacja efektów losowych wynosi \(\pm1\). Model jest wtedy zbyt złożony w stosunku do danych i należy go uprościć (usunąć losowe nachylenie lub wybrać kowariancję niezależne).

Stopnie swobody Satterthwaite'a – w modelu mieszanym nie ma jednej prostej liczby stopni swobody dla testu efektu stałego. Zmienna opisująca grupę (np. płeć pacjenta) ma w istocie tyle niezależnych obserwacji, ile jest grup, a zmienna zmieniająca się wewnątrz grup – więcej. Metoda Satterthwaite'a wyznacza dla każdego współczynnika osobne, zwykle niecałkowite stopnie swobody, dzięki czemu test \(t\) jest wiarygodny także przy małej liczbie grup.

Struktura reszt w grupie – domyślnie reszty \(\epsilon_{ij}\) są niezależne, a cała korelacja wewnątrz grupy wynika z efektów losowych. Dla pomiarów w czasie można dodatkowo założyć, że reszty sąsiednich pomiarów są skorelowane: AR(1) (korelacja \(\rho\) maleje wykładniczo z odległością w kolejności wierszy grupy: \(\rho, \rho^2, \rho^3, \ldots\)) lub wymienność (CS) (jednakowa korelacja wszystkich par reszt w grupie).

Poziom 2 – drugi, nadrzędny poziom grupowania, np. pacjent (poziom 1) zagnieżdżony w ośrodku (poziom 2). Model dostaje wtedy dwa losowe wyrazy wolne, a raport – wariancję i ICC osobno dla każdego poziomu.

Jak ustawić analizę

Statystyka \(\to\) Modele mieszane \(\to\) Regresja wieloraka mieszana.

Okno analizy ma układ znany z regresji wielorakiej. W jego lewej części wskazujemy zmienne:

• Zmienna zależna Y – zmienna ciągła, która ma być wyjaśniana,
• Zmienna grupująca (ID) – kolumna identyfikująca grupę (pacjenta, ośrodek, miot); może być liczbowa lub tekstowa,
• Zmienne niezależne X1, X2, … – predyktory; zmienne fikcyjne i interakcje przygotowuje się przyciskami Zm. fikcyjne i Interakcje tak jak w regresji wielorakiej (rozdział Przygotowanie zmiennych do analizy w modelach wielowymiarowych).

Prawa część okna zawiera opcje modelu mieszanego. Poniżej opisano, co zmienia każda z nich.

Efekty losowe

• Losowe nachylenia (z X) – lista wypełnia się zmiennymi wybranymi jako \(X\). Zaznaczenie zmiennej oznacza, że jej efekt ma się różnić między grupami (np. tempo zmian w czasie różne dla każdego pacjenta). Losowy wyraz wolny jest zawsze w modelu i nie trzeba go wybierać. Zmienna z losowym nachyleniem musi zmieniać się wewnątrz grup; zmienna stała w grupie (np. płeć pacjenta) nie może mieć losowego nachylenia – program zgłosi wtedy błąd. Każde losowe nachylenie dodaje do modelu wariancję i (przy kowariancji pełnej) korelacje, więc wymaga większej liczby grup; rozsądnie jest zaczynać od jednego nachylenia, najczęściej dla czasu.
• Kowariancja efektów losowych – pełna z korelacją szacuje także korelację między losowym wyrazem wolnym a nachyleniem (oraz między nachyleniami); niezależne przyjmuje te korelacje równe 0. Wariant niezależne ma mniej parametrów, jest stabilniejszy przy małej liczbie grup i zalecany, gdy korelacja w modelu pełnym wychodzi bliska 0 albo \(\pm1\) (dopasowanie osobliwe). Przy braku losowych nachyleń opcja nie ma znaczenia.

Estymacja

• Estymacja – REML (ograniczona największa wiarygodność, domyślnie) daje nieobciążone oszacowania wariancji i jest zalecana dla modelu końcowego; ML (największa wiarygodność) zaniża nieco wariancje, ale tylko wiarygodność ML pozwala porównywać kryteriami informacyjnymi modele różniące się efektami stałymi. Wybór metody nie zmienia w zauważalny sposób współczynników efektów stałych, zmienia natomiast wariancje efektów losowych, ich błędy i wartości \(-2LL\), AIC, BIC. Zasada: dobór zmiennych – ML, raportowanie modelu – REML.
• Stopnie swobody efektów stałych – Satterthwaite (domyślnie) wyznacza dla każdego współczynnika własne stopnie swobody i test \(t\); zalecane zawsze, a szczególnie przy małej liczbie grup. Test z Walda traktuje statystykę jak w dużej próbie (rozkład normalny, bez stopni swobody); daje mniejsze wartości \(p\) i jest odpowiedni tylko przy dużej liczbie grup.
• Struktura reszt w grupie – niezależne (domyślnie), AR(1) (kolejność wierszy) lub wymienność (CS). AR(1) wymaga, by wiersze każdej grupy były w arkuszu uporządkowane w czasie i mniej więcej równoodległe. W raporcie pojawia się parametr \(\rho\). Wymienność razem z losowym wyrazem wolnym jest nieidentyfikowalna (obie opisują tę samą stałą korelację w grupie), więc ma sens głównie jako porównanie z modelem AR(1).
• Poziom 2 – grupa nadrzędna (opcjonalnie) – kolumna wskazująca grupę wyższego poziomu, w której zagnieżdżona jest zmienna grupująca (np. ośrodek dla pacjentów). Dodaje losowy wyraz wolny poziomu 2, a w raporcie osobną wariancję, ICC i tabelę BLUP tego poziomu.

Elementy raportu

• Model pusty i ICC – dodatkowo liczony jest model bez zmiennych niezależnych, z którego wyznacza się ICC z przedziałem ufności; warto zostawić włączone, bo ICC to pierwsza liczba, na którą należy spojrzeć.
• Test efektów losowych (LR) – testy ilorazu wiarygodności dla każdego składnika losowego osobno oraz dla wszystkich łącznie (model mieszany wobec zwykłej regresji wielorakiej). Wydłuża obliczenia, bo wymaga dopasowania dodatkowych modeli.
• Efekty losowe grup (BLUP) – dołącza tabelę przewidywanych efektów każdej grupy z błędami; przy setkach grup tabela jest długa, a te same wartości pokazuje wykres caterpillar.
• M. korelacji i kowariancji – macierze korelacji i kowariancji ocen efektów stałych.
• Analiza reszt – wartości przewidywane (warunkowe i populacyjne), reszty i reszty standaryzowane dla każdej obserwacji.
• Dołącz wykresy – wykres caterpillar efektów losowych grup oraz (przy losowym nachyleniu) trajektorie grup na tle prostej populacyjnej.

Pozostałe elementy okna (poziom istotności, filtr danych, łączenie raportów) działają tak jak w innych analizach.

Wyniki i interpretacja

Liczności – liczba obserwacji, liczba grup oraz minimalna, medianowa i maksymalna liczba obserwacji na grupę. Przy mniej niż dziesięciu grupach raport ostrzega, że oszacowania wariancji efektów losowych są niepewne.

Estymacja – metoda (REML/ML), liczba iteracji, informacja o zbieżności oraz o dopasowaniu osobliwym. Brak zbieżności lub dopasowanie osobliwe oznacza, że model należy uprościć, zanim zaczniemy interpretować pozostałe wyniki.

Dopasowanie modelu – logarytm wiarygodności, \(-2LL\), kryteria informacyjne AIC, AICc i BIC (im mniejsze, tym lepszy model) oraz \(R^2\) brzegowe i warunkowe. Kryteria informacyjne wolno porównywać tylko między modelami liczonymi tą samą metodą; modele różniące się efektami stałymi – tylko przy ML.

Test LR efektów losowych (vs regresja wieloraka, ML) – czy uwzględnienie struktury grupowej jest w ogóle potrzebne. Istotny wynik oznacza, że zwykła regresja wieloraka nie jest dla tych danych poprawna. Wynik nieistotny przy hierarchicznym układzie badania nie zwalnia jednak z użycia modelu mieszanego.

Model pusty – wariancja między grupami \(\tau^2\), wariancja reszt \(\sigma^2\) i ICC z przedziałem ufności.

Efekty stałe (warunkowe) – współczynnik \(b\), jego błąd, przedział ufności, statystyka \(t\) ze stopniami swobody Satterthwaite'a (mogą być niecałkowite) lub statystyka \(z\) oraz wartość \(p\). Interpretacja jak w regresji wielorakiej: \(b\) to zmiana \(Y\) na jednostkę \(X\) przy pozostałych zmiennych ustalonych, w obrębie tej samej grupy. Stopnie swobody podpowiadają, na jakim poziomie „działa” zmienna: liczba bliska liczbie grup oznacza zmienną opisującą grupy, liczba bliska liczbie obserwacji – zmienną zmieniającą się wewnątrz grup.

Efekty losowe – dla każdego składnika wariancja i odchylenie standardowe. Kolumny LR \(\chi^2\), df i \(p\) zawierają test, czy dany składnik jest potrzebny (model z nim wobec modelu bez niego, ML, z korektą granicy). Ostatni wiersz to wariancja reszt. Przy strukturze AR(1) lub CS raport podaje parametr \(\rho\), a przy poziomie 2 – wariancję i ICC tego poziomu.

Korelacje efektów losowych – przy więcej niż jednym składniku macierz korelacji między losowym wyrazem wolnym a nachyleniami. Wartość bliska \(\pm1\) sygnalizuje dopasowanie osobliwe.

Efekty losowe grup (BLUP) – przewidywane odchylenie każdej grupy od populacji wraz z błędem (dla wyrazu wolnego i każdego losowego nachylenia). Suma efektu stałego i BLUP daje własny współczynnik grupy. Na wykresie caterpillar grupy, których przedział nie obejmuje zera, wyraźnie odbiegają od przeciętnej.

Analiza reszt – wartości przewidywane warunkowe (z efektami losowymi grupy, czyli „dla tego pacjenta”) i populacyjne (same efekty stałe, czyli „dla przeciętnego pacjenta”), reszty oraz reszty standaryzowane; obserwacje o \(|reszta|\ge3\) są wyróżnione jako odstające.

Wykresy – Efekty losowe grup (BLUP \(\pm\) błąd): grupy uporządkowane według wartości efektu, z przedziałem błędu i linią zera; osobny panel dla wyrazu wolnego i każdego nachylenia. Trajektorie grup: dla modelu z losowym nachyleniem cienkie proste dopasowane dla poszczególnych grup na tle grubej prostej populacyjnej – najbardziej obrazowe podsumowanie pomiarów powtarzanych.

Wykonanie i interpretacja