Regresja logistyczna mieszana (uogólniony model mieszany)

Kiedy stosować:

Regresja logistyczna mieszana (uogólniony liniowy model mieszany, ang. generalized linear mixed model, GLMM) jest rozszerzeniem regresji logistycznej na dane, w których obserwacje tworzą grupy: powtarzane oceny tego samego pacjenta (np. obecność objawu na kolejnych wizytach), pacjenci z tego samego ośrodka, zwierzęta z tego samego stada, uczniowie z jednej klasy (zob. wprowadzenie do modeli mieszanych). Wynik binarny (0/1) obserwacji z jednej grupy jest skorelowany, a zwykła regresja logistyczna, która tę korelację pomija, zaniża błędy standardowe i zawyża istotność.

Model opisuje logit prawdopodobieństwa zdarzenia:

\[ \ln\frac{P(Y_{ij}=1)}{1-P(Y_{ij}=1)}=\beta_0+\beta_1X_{1ij}+\ldots+\beta_kX_{kij}+u_{0j}+u_{1j}X_{1ij}, \]

gdzie:
\(Y_{ij}\) – wynik \(i\)-tej obserwacji w grupie \(j\) (0 lub 1),
\(\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 efekty losowe mają rozkład normalny o średniej 0 i macierzy kowariancji \(\Sigma\) na skali logitowej.

Losowy wyraz wolny oznacza, że każda grupa ma własne „bazowe” szanse zdarzenia: grupa o \(u_{0j}=1\) ma szanse \(e^{1}\approx2.7\) razy większe niż grupa przeciętna, przy tych samych wartościach zmiennych. W przeciwieństwie do regresji wielorakiej mieszanej model nie ma wariancji reszt: na skali logitowej odpowiada jej stała \(\pi^2/3\approx3.29\).

Model ten zastępuje ANOVA Q-Cochrana, gdy pomiary powtarzane są niekompletne lub gdy oprócz czynnika powtarzanego w modelu mają być zmienne ciągłe i zmienne opisujące grupę.

Warunki stosowania

zmienna zależna dwuwartościowa,
efekty losowe o rozkładzie normalnym na skali logitowej; grupy niezależne od siebie,
co najmniej 3 grupy, a do wiarygodnego oszacowania wariancji zwykle kilkanaście lub więcej; przy losowych nachyleniach potrzeba też kilku obserwacji w każdej grupie,
wystarczająca liczba zdarzeń (jedynek) i braków zdarzeń (zer) – jak w regresji logistycznej,
braki danych losowe.

Uwaga! Dane muszą być w układzie długim: każdy wiersz to jedna obserwacja (jeden wynik 0/1), a osobna kolumna wskazuje grupę. Dane zagregowane (liczba zdarzeń i liczba prób w grupie) należy najpierw rozpisać na pojedyncze wiersze.

Definicje miar

Warunkowy iloraz szans \(OR=e^{b}\) – mówi, ile razy zmieniają się szanse zdarzenia, gdy \(X\) rośnie o jednostkę w obrębie tej samej grupy (u tego samego pacjenta, w tym samym ośrodku). Zwykła regresja logistyczna daje iloraz szans populacyjny, uśredniony po wszystkich grupach. Warunkowy OR jest co do wartości bezwzględnej zawsze większy (dalej od 1) niż populacyjny, a różnica rośnie wraz z wariancją efektów losowych; przy wariancji \(\tau^2\) populacyjny współczynnik wynosi w przybliżeniu \(b/\sqrt{1+0.346\,\tau^2}\). Oba są poprawne, ale odpowiadają na różne pytania i nie należy ich porównywać wprost.

Wariancja efektów losowych \(\tau^2\) na skali logitowej i jej pierwiastek – odchylenie standardowe. Odchylenie standardowe losowego wyrazu wolnego równe 1 oznacza, że grupa o poziomie wyższym o jedno odchylenie ma szanse zdarzenia \(e^{1}\approx2.7\) razy większe niż grupa przeciętna, a grupy skrajne (o \(\pm2\) odchylenia) różnią się szansami \(e^{4}\approx55\) razy.

ICC na skali logitowej – z modelu pustego:

\[ ICC=\frac{\tau^2}{\tau^2+\pi^2/3}, \]

interpretowany jako korelacja utajonych skłonności do zdarzenia w obrębie grupy (udział zmienności między grupami w zmienności utajonej zmiennej ciągłej stojącej za wynikiem 0/1).

\(R^2\) brzegowe i warunkowe (Nakagawa) – na skali logitowej, z \(\pi^2/3\) jako wariancją reszt: udział wariancji wyjaśnianej przez efekty stałe oraz przez efekty stałe i losowe razem.

Aproksymacja wiarygodności – wiarygodność modelu wymaga scałkowania po efektach losowych, czego nie da się zrobić dokładnie. Laplace przybliża całkę jedną funkcją wokół jej maksimum – szybko i dostatecznie dokładnie przy większych grupach. Kwadratura Gaussa-Hermite'a (adaptacyjna, 7, 15 lub 25 węzłów) liczy całkę numerycznie w kilku punktach i jest dokładniejsza, zwłaszcza przy małych grupach (kilka obserwacji na grupę) i dużej wariancji; kosztem jest czas obliczeń. Kwadratura jest dostępna tylko dla modelu z samym losowym wyrazem wolnym; przy losowych nachyleniach program stosuje aproksymację Laplace'a. Wyniki różnych aproksymacji różnią się zwykle na trzecim miejscu znaczącym.

Test LR efektów losowych, BLUP, dopasowanie osobliwe – jak w regresji wielorakiej mieszanej; BLUP są odchyleniami grup na skali logitowej.

Predykcja warunkowa i populacyjna – prawdopodobieństwo przewidywane z efektami losowymi grupy („dla tego pacjenta”) oraz z samych efektów stałych („dla przeciętnego pacjenta”). Klasyfikacja i krzywa ROC mogą być policzone dla obu; predykcja warunkowa jest zawsze dokładniejsza, bo korzysta z informacji o grupie, ale nie da się jej użyć dla nowych, nieznanych grup.

Jak ustawić analizę

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

W oknie wskazujemy Zmienną zależną Y i w polu Grupa (1) wartość oznaczającą zdarzenie, Zmienną grupującą (ID) oraz Zmienne niezależne X1, X2, …; zmienne fikcyjne i interakcje przygotowuje się jak w regresji logistycznej. Opcje modelu mieszanego:

Losowe nachylenia (z X) – zmienne spośród \(X\), których efekt ma się różnić między grupami; muszą zmieniać się wewnątrz grup. Losowy wyraz wolny jest zawsze w modelu. Losowe nachylenia w modelu logistycznym wymagają większej liczby obserwacji na grupę niż w modelu liniowym, bo każda obserwacja niesie tylko jeden bit informacji,
Kowariancja efektów losowych – pełna z korelacją lub niezależne; ta druga jest stabilniejsza przy małej liczbie grup,
Aproksymacja wiarygodności – Laplace (domyślnie, zawsze dostępna) albo kwadratura G-H z 7, 15 lub 25 węzłami (tylko losowy wyraz wolny). Zmiana aproksymacji zmienia nieznacznie wszystkie liczby w raporcie; do publikacji, przy małych grupach, warto użyć kwadratury z 15–25 węzłami,
Model pusty i ICC, Test efektów losowych (LR), Efekty losowe grup (BLUP), M. korelacji i kowariancji, Analiza reszt – jak w regresji wielorakiej mieszanej,
Klasyfikacja i krzywa ROC – czułość, swoistość, udział poprawnych klasyfikacji i AUC dla wskazanej linii odcięcia, osobno dla predykcji warunkowej i populacyjnej,
Dołącz wykresy – wykres caterpillar efektów losowych grup.

Wyniki i interpretacja

Liczności – liczba obserwacji, liczność \(Y=1\), liczba grup, minimum, mediana i maksimum obserwacji na grupę. Przy medianie poniżej 3 obserwacji na grupę aproksymacja Laplace'a jest najmniej dokładna – warto wtedy użyć kwadratury.

Aproksymacja i zbieżność – metoda (Laplace lub Gauss-Hermite z liczbą węzłów), liczba iteracji, zbieżność, dopasowanie osobliwe.

Dopasowanie modelu – logarytm wiarygodności, \(-2LL\), AIC, AICc, BIC (estymacja ML), \(R^2\) brzegowe i warunkowe oraz – dla porównania – \(-2LL\) zwykłej regresji logistycznej (bez efektów losowych) i \(-2LL\) modelu pustego.

Test LR efektów losowych (vs regresja logistyczna) – czy struktura grupowa jest potrzebna; wartość \(p\) skorygowana dla hipotezy na granicy.

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

Efekty stałe (warunkowe) – \(b\), błąd, przedział ufności, statystyka \(z\), wartość \(p\) oraz iloraz szans (warunkowy) z przedziałem ufności. Test jest testem \(z\) Walda (bez stopni swobody), tak jak w regresji logistycznej.

Efekty losowe – wariancja i odchylenie standardowe każdego składnika na skali logitowej, testy LR potrzeby składników; ostatni wiersz to stała \(\pi^2/3\). Przy więcej niż jednym składniku – macierz korelacji efektów losowych.

Efekty losowe grup (BLUP) – odchylenia grup od populacji na skali logitowej z błędem. Grupa o BLUP \(=0.5\) ma szanse zdarzenia \(e^{0.5}\approx1.65\) razy większe niż grupa przeciętna. Grupy z małą liczbą obserwacji lub bez zdarzeń mają BLUP ściągnięte w stronę zera.

Klasyfikacja i krzywa ROC – dla predykcji warunkowej (z BLUP) i populacyjnej: czułość, swoistość, udział poprawnych, AUC z błędem i przedziałem ufności. Różnica AUC między nimi pokazuje, ile zdolności rozróżniania wnosi znajomość grupy.

Analiza reszt – prawdopodobieństwo przewidywane warunkowe i populacyjne oraz reszty Pearsona dla każdej obserwacji.

Wykonanie i interpretacja

Przykład 1. (plik stadaCBPP.pqs)
W badaniu zakaźnego zapalenia płuc i opłucnej bydła (CBPP) w Etiopii (Lesnoff i in. 2004[189]) obserwowano 15 stad zebu przez cztery kolejne okresy. Dla każdego zwierzęcia w każdym okresie zapisano, czy pojawił się u niego nowy przypadek choroby. Dane są w układzie długim: Stado (1–15), Okres (1–4) i Zachorowanie (1 – nowy przypadek, 0 – brak); 842 obserwacje, w tym 99 zachorowań, od 26 do 96 obserwacji na stado. Pytamy, czy zapadalność zmienia się w kolejnych okresach epidemii i jak bardzo stada różnią się między sobą.

Zwierzęta z jednego stada dzielą warunki, kontakty i odporność, więc ich wyniki nie są niezależne – stado jest naturalną zmienną grupującą. Budujemy model: Zachorowanie jako \(Y\) z wartością 1 w polu Grupa (1), Stado jako zmienną grupującą, Okres jako zmienną niezależną. Okres jest zmienną jakościową o czterech kategoriach, więc przyciskiem Zm. fikcyjne zamieniamy go na trzy zmienne fikcyjne z okresem 1 jako kategorią odniesienia. Zostawiamy sam losowy wyraz wolny (okres nie ma losowego nachylenia), aproksymację Laplace, włączamy Model pusty i ICC, Test efektów losowych (LR), Efekty losowe grup (BLUP), Klasyfikacja i krzywa ROC oraz Dołącz wykresy.

Krok 1 – dopasowanie. 842 obserwacje, 99 zdarzeń, 15 grup (26/61/96 obserwacji na grupę), zbieżność osiągnięta, bez dopasowania osobliwego. \(-2LL\) modelu wynosi 555.00 (AIC \(=565.00\)), podczas gdy zwykła regresja logistyczna z tymi samymi zmiennymi ma \(-2LL=569.01\). Test LR efektów losowych: \(\chi^2=14.01\), \(df=1\), \(p=0.00009\) (z korektą granicy) – stada istotnie różnią się zapadalnością i model bez efektu losowego byłby niepoprawny.

Krok 2 – ICC. W modelu pustym wariancja między stadami wynosi 0.659, stąd ICC \(=0.659/(0.659+3.29)=0.167\) (95% CI: 0.037–0.297): około 17% utajonej zmienności skłonności do zachorowania to różnice między stadami.

Krok 3 – efekty stałe. Wyraz wolny \(b_0=-1.399\) odpowiada prawdopodobieństwu zachorowania w okresie 1 w przeciętnym stadzie równemu \(1/(1+e^{1.399})=0.198\). Warunkowe ilorazy szans dla kolejnych okresów względem okresu 1: okres 2 – \(OR=0.371\) (95% CI: 0.203–0.676, \(p=0.0012\)), okres 3 – \(OR=0.324\) (0.171–0.614, \(p=0.0005\)), okres 4 – \(OR=0.206\) (0.089–0.476, \(p=0.0002\)). W tym samym stadzie szanse nowego zachorowania w okresie 2 są o 63% mniejsze niż w okresie 1, a w okresie 4 – o 79% mniejsze: epidemia wygasa. Słowo „warunkowy” przypomina, że porównujemy okresy wewnątrz stada; iloraz szans uśredniony po wszystkich stadach (ze zwykłej regresji logistycznej) byłby nieco bliższy 1.

Krok 4 – efekty losowe. Wariancja losowego wyrazu wolnego wynosi 0.413, odchylenie standardowe 0.642 na skali logitowej. Stado o poziomie wyższym o jedno odchylenie ma szanse zachorowania \(e^{0.642}\approx1.9\) razy większe niż stado przeciętne, a różnica między stadami skrajnymi (\(\pm2\) odchylenia) to \(e^{2.57}\approx13\)-krotność szans. \(R^2\) brzegowe \(=0.088\), warunkowe \(=0.190\): okres wyjaśnia niecałe 9% zmienności na skali logitowej, a razem ze stadem – 19%.

Krok 5 – stada nietypowe. BLUP stada 1 wynosi \(+0.590\) (błąd 0.348): jego szanse zachorowania są \(e^{0.59}\approx1.8\) razy większe niż przeciętnego stada; w okresie 1 daje to prawdopodobieństwo \(1/(1+e^{1.399-0.590})=0.308\) zamiast 0.198. Stado 2 ma BLUP \(-0.299\) (szanse o 26% mniejsze). Wykres caterpillar porządkuje stada od najmniej do najbardziej zagrożonego; stada, których przedział nie obejmuje zera, wyraźnie odbiegają od przeciętnej. Warto pamiętać, że stada z małą liczbą zwierząt mają BLUP ściągnięte w stronę zera.

Krok 6 – klasyfikacja. Tabela klasyfikacji i ROC pokazuje dwie wersje: predykcję warunkową (z BLUP stada) i populacyjną (sam okres). AUC predykcji warunkowej jest wyższe – znajomość stada poprawia rozróżnianie zwierząt, które zachorują, od tych, które nie zachorują. Wersja populacyjna mówi, jak model sprawdzi się w nowym, nieznanym stadzie.

OR warunkowe: okres 2 – 0.37, okres 3 – 0.32, okres 4 – 0.21 (wszystkie \(p<0.002\)); \(SD\) stad \(=0.64\) (LR \(p=0.00009\)), ICC \(=0.17\)

Zapadalność maleje w kolejnych okresach epidemii, a stada istotnie różnią się bazowym ryzykiem.

Co zmienia aproksymacja. Ten sam model z kwadraturą G-H, 25 węzłów daje \(-2LL=554.92\), wariancję stad 0.419 zamiast 0.413 i współczynniki różniące się na trzecim miejscu po przecinku (np. okres 4: \(b=-1.579\) zamiast \(-1.580\)). Przy stadach liczących 26–96 zwierząt aproksymacja Laplace'a jest w pełni wystarczająca; gdyby grupy liczyły po 2–3 obserwacje, różnice byłyby większe i kwadratura stałaby się metodą z wyboru.

DLA ZAINTERESOWANYCH

Estymacja: Wiarygodność brzegowa wymaga scałkowania po efektach losowych, które nie ma postaci zamkniętej. Program stosuje aproksymację Laplace'a (Bates i in. 2015[184]): dla każdej grupy znajduje modę \(\hat u_j\) penalizowanej log-wiarygodności

\[ f_j(u)=-\sum_{i\in j}\left[y_{ij}\eta_{ij}-\ln(1+e^{\eta_{ij}})\right]+\frac{1}{2}u^Tu,\qquad \eta_{ij}=x_{ij}^T\beta+z_{ij}^T\Lambda u, \]

i przybliża dewiancję przez \(\sum_j\left[2f_j(\hat u_j)+\ln\det(I+\Lambda^TZ_j^TW_jZ_j\Lambda)\right]\), gdzie \(W_j\) jest macierzą wag \(p(1-p)\), a \(\Sigma=\Lambda\Lambda^T\). Dla jednego składnika losowego całka może być liczona adaptacyjną kwadraturą Gaussa-Hermite'a wokół mody, z węzłami przesuniętymi i przeskalowanymi według mody i krzywizny każdej grupy. Parametry \((\theta,\beta)\) są optymalizowane łącznie metodą Neldera-Meada ze startem z regresji logistycznej.

Błędy standardowe: z odwrotności hesjanu dewiancji po wszystkich parametrach \((\theta,\beta)\), liczonego numerycznie w optimum, jak w pakiecie lme4.

ICC: \(\tau^2/(\tau^2+\pi^2/3)\) z modelu pustego, z przedziałem ufności metodą delta.

\(R^2\) Nakagawy: jak w regresji wielorakiej mieszanej, z \(\pi^2/3\) w miejsce wariancji reszt.

Test na granicy: statystyka \(2(\ell_1-\ell_0)\) dla hipotezy o zerowej wariancji ma rozkład będący mieszanką 50:50 rozkładów \(\chi^2_{q-1}\) i \(\chi^2_q\) (Self i Liang 1987[193]); dla jednego parametru wartość \(p\) z \(\chi^2_1\) dzieli się przez dwa.

Efekt populacyjny: przy losowym wyrazie wolnym o wariancji \(\tau^2\) współczynnik populacyjny (brzegowy) wynosi w przybliżeniu \(\beta/\sqrt{1+0.346\,\tau^2}\) (Zeger i in. 1988[195]).