Analiza przeżywalności chorych z udarem mózgu
mgr Andrzej Stanisz
Collegium Medicum UJ

Udar mózgu jest trzecią co do częstości przyczyną śmierci i główną przyczyną trwałego kalectwa i braku samodzielności u osób dorosłych. Udar mózgu niedokrwienny powstaje w wyniku zwężenia lub zamknięcia przez zakrzep tętnicy mózgowej. Do tworzenia zakrzepu najczęściej dochodzi w miejscu uszkodzenia przez zmiany miażdżycowe wewnętrznej ściany naczynia. W następstwie tego procesu, do pewnego obszaru mózgu nie dochodzi krew zawierająca tlen, i komórki nerwowe w tej okolicy obumierają. W Polsce rejestruje się około 60 000 nowych zachorowań na udar rocznie. Zapadalność na tę chorobę jest podobna jak w innych krajach europejskich (około 175/100 000 mężczyzn i 125/100 000 kobiet), ale umieralność okołoudarowa i niesprawność po udarze są znacznie większe. Począwszy od lat 70., dzięki powszechnemu wprowadzeniu profilaktyki i leczenia chorób układu krążenia oraz poprawie opieki nad pacjentami z udarem, w wielu krajach świata notuje się systematyczny spadek umieralności i niesprawności związanej z udarem mózgu. W Polsce tendencji tej, niestety, się nie obserwuje. Z tego powodu problem zapobiegania i leczenia udaru w naszym kraju jest szczególnie ważny. Jedną z pierwszych inicjatyw Narodowego Programu Profilaktyki i Leczenia Udaru mózgu było opracowanie w Krakowskiej Klinice Neurologii bazy danych chorych na niedokrwienny udar mózgu. Jest to jedna z niewielu tego typu baza w Polsce. Zawiera ona informację o 1463 pacjentach (w tym 755 kobiet). Jak pokazuje to poniższy rysunek w bazie zawarte mamy wyniki wszystkich najważniejszych badań

W myśl opracowanych (w 1997 r.) ogólnych zasad postępowania lekarz w pierwszych godzinach hospitalizacji powinien przede wszystkim:

Do najważniejszych elementów tego postępowania należy ocena neurologiczna oraz wykonanie badań diagnostycznych (EKG, badanie biochemiczne i tomografia komputerowa głowy). Celem tych badań jest wykrycie zarówno czynników etiologicznych jak i wczesnych powikłań. To postępowanie terapeutyczne ma swoje odbicie w bazie. Oprócz bowiem podstawowych informacji o pacjencie zawiera ona bowiem aż pięć tabel (Biochemia_1, ...., Biochemia_5) z wynikami badań hematologicznych i biochemicznych. Wiele z nich jak widzimy to na poniższym rysunku było wykonywane wielokrotnie w trakcie hospitalizacji. W bazie spotykamy ponadto tabelę poświęconą wynikom EKG, tomografii komputerowej głowy oraz infekcjom towarzyszącym udarowi mózgu. W celach obiektywizacji oceny chorego oraz monitorowania przebiegu leczenia zalecane jest posługiwanie się odpowiednimi skalami. Tak więc wyniki badań neurologicznych i opis nasilenia poszczególnych objawów w prostych i powszechnie stosowanych skalach (Skala Skandynawska, Rankin Skale oraz Indeks Barthela) zawiera osobna tabela Stan Kliniczny.
Ze względu na zagrożenie powtórnym udarem w pierwszych miesiącach po zachorowaniu ważnym zadaniem jest zebranie informacji o czynnikach ryzyka, powikłaniach oraz wynikach leczenia. Informacje o tym zawarte są w tabelach Czynniki Ryzyka, Powikłania oraz Infekcje. Ostatnia tabela zawiera informacje o wynikach leczenia i czasie przeżycia badanych pacjentów.


Rys. 1 Fragment arkusza opisu tabeli bazy udarowej

Nie sposób opisać wszystkich analiz statystycznych przeprowadzanych na danych z bazy. Przykładowo pokażemy kilka najciekawszych związanych z różnymi analizami statystycznymi.

Przykład I - zastosowanie regresji logistycznej.

Szukamy czynników wpływających na przeżycie 90 dni po udarze mózgu. Zgodnie z wcześniejszymi badaniami i obserwacjami wzięto pod uwagę następujące zmienne:


Do analizy statystycznej użyto modelu regresji logistycznej. Przebadano 1150 pacjentów. Fragment arkusza danych widoczny jest na poniższym rysunku.

Rys. 2 Fragment arkusza danych z przykładu

Po wprowadzeniu danych i wykonaniu estymacji metodą Rosenbrocka i quasi-Newtona otrzymujemy następujący arkusz wyników:


Rys. 3 Arkusz wyników dla danych z przykładu

Wszystkie analizowane zmienne z wyjątkiem zmiennej CPK (p = 0,093) okazały się istotne statystycznie na poziomie p<0,05. Z tego wynika, że powinniśmy rozważyć prostszy model bez zmiennej NAP. Otrzymujemy wówczas arkusz wyników widoczny na poniższym rysunku.


Rys. 4 Arkusz wyników estymacji parametrów modelu

Wliczone wartości pozwalają zapisać równanie regresji logistycznej w postaci:

Ujemne szacunki parametrów odpowiadające zmiennym CUKIER_0, WIELCT, WIEK wskazują, że wzrost tych wielkości powoduje zmniejszenie się prawdopodobieństwa przeżycia 90 dni. Największy wzrost jest widoczny dla zmiennej WIELCT. Korzystając z jednostkowego ilorazu szans dla zmiennej WIELCT (0,5866) możemy zbudować tabelkę obrazującą zmniejszanie się szans na przeżycie 90 dni wraz ze wzrostem wielkości ogniska udaru w skali od 1 do 5 (przy niezmieniających się wartościach pozostałych zmiennych). Ma ona postać następującą:

Wielkość ogniska 1 2 3 4 5
Wielkość zmniejszenia           0,586 0,344 0,202 0,118

Widzimy, że osoba o wielkości ogniska "3" ma w stosunku do osoby o wielkości ogniska "1" o trzy razy mniejszą szansę przeżycia 90 pierwszych dni od momentu zachorowania. Jeszcze mniejsze szanse (0,12) ma osoba o największym ognisku udaru. Podobnie w oparciu o iloraz szans dla zmiennej WIEK możemy obliczyć, że osoby w wieku 50 lat mają w stosunku do osób 40 letnich zmniejszoną szansę przeżycia 90 dni prawie o połowę (0,567). Mamy mianowicie -

Z kolei współczynnik dodatni przy zmiennej SLALA_SK wskazuje, że wzrost tej wielkości (stanu klinicznego pacjenta) powoduje wzrost prawdopodobieństwa przeżycia 90 dni. Jest to zgodne z interpretacją medyczną. Skala stanu klinicznego jest bowiem określona od 0 do 58, przy czym 0 punktów określa sytuację najgorszą a 58 sytuację. najlepszą. Również wartość chi-kwadrat (χ²) porównująca aktualny model z modelem tylko z wyrazem wolnym jest statystycznie wysoce istotna (p<0,00000001). Model jest więc dobrze dopasowany do naszych danych. Na zakończenie klikamy na przycisku Rozkład reszt aby sprawdzić kształt rozkładu. Wykres, który uzyskaliśmy przedstawiony jest na poniższym rysunku. Widać wyraźnie, że reszty mają zasadniczo rozkład zbliżony do rozkładu normalnego.


Rys. 5. Histogram reszt

Budując model logistyczny często dołączamy lub usuwamy pewne zmienne. Tworzymy więc pewne hierarchiczne modele. Mając bowiem już pewien model, zastanawiamy się, czy dołączenie pewnej zmiennej niezależnej "poprawi" nasz model. Jest to tak zwane podejście "od szczegółu do ogółu". Podobnie czasami możemy się zastanawiać nad problemem usunięcia pewnych zmiennych z analizowanego przez nas modelu. Tego typu postępowanie jest charakterystyczne dla drugiego podejścia zwanego "od ogółu do szczegółu". W obu przypadkach, gdy usuwamy lub dodajemy jedną lub kilka zmiennych STATISTICA automatycznie wyliczy przyrostową dobroć dopasowania. Wykorzystywany do tego celu jest przyrostowy test chi-kwadrat.

Przykład II zastosowanie analizy log-liniowej

Szukamy również czynników wpływających na przeżycie 90 dni po udarze mózgu. W badaniu uwzględniono następujące zmienne:

Wyniki badań dla 50 pacjentów przedstawia poniższa tabela.

WIEK MIKKROAL HIPER P90
2 1 2 tak
2 2 1 tak
2 1 1 tak
1 1 2 tak
1 3 1 nie
1 3 1 nie
2 1 1 tak
1 1 2 tak
1 1 2 tak
2 1 2 nie
2 1 2 nie
2 2 1 nie
2 1 1 tak
2 1 2 tak
2 1 2 tak
1 1 2 tak
2 1 2 tak
1 3 1 tak
1 2 1 tak
1 1 1 tak
1 1 2 tak
1 1 2 tak
2 1 2 tak
2 2 2 nie
2 1 2 nie
WIEK MIKKROAL HIPER P90
1 1 1 tak
1 1 2 tak
2 1 2 tak
2 2 2 nie
1 1 2 tak
1 3 2 tak
1 2 2 tak
1 1 2 tak
1 1 1 tak
1 1 1 tak
1 1 2 tak
1 1 2 tak
1 1 2 tak
1 1 1 tak
1 1 2 tak
2 1 1 tak
2 1 2 nie
1 1 1 tak
1 1 2 tak
1 1 2 tak
2 2 2 nie
2 1 1 tak
1 1 2 tak
1 1 2 tak
2 1 2 nie

Do analizy statystycznej użyjemy modelu log liniowego. Po wprowadzeniu danych w panelu początkowym na liście Plik wejściowy wybieramy opcję Dane surowe. Po określeniu zmiennych klikamy OK, aby otworzyć okno Specyfikacja modelu log-liniowego. Musimy teraz określić model do testowania. W tym celu przejrzymy tabelę jednoczesnych testów na wszystkie interakcje k - czynników. Tabela ta pokazana jest na poniższym rysunku.


Rys. 6. Tabela testów na wszystkie interakcje k-czynnikowe

Jak pokazuje drugi wiersz model bez interakcji dwuczynnikowych odrzucamy (p = 0,000139). Włączenie więc wszystkich interakcji rzędu drugiego do modelu poprawia jego dopasowanie. Jednak rozszerzenie modelu o wszystkie interakcje rzędu trzeciego (k = 3 - wiersz trzeci) nie daje istotnej poprawy (p = 0,8177). Oznacza to, że najmniej złożony model, który pasuje do tabeli liczebności to model bez żadnych powiązań trójwymiarowych. Wymagane są natomiast niektóre zależności dwuwymiarowe. Jakie zależności dwuwymiarowe dołączyć do modelu, a jakie odrzucić zadecydujemy po przeglądnięciu tabeli wszystkich zależności cząstkowych i brzegowych Tabela ta widoczna jest poniżej.


Rys. 7. Testy związku brzegowego i cząstkowego

Czarna kropka zaznacza efekty dla których zależności cząstkowe i brzegowe są istotne. Przykładowo spójrzmy na efekt 24. Efekt ten prezentuje zależność lub interakcję między czynnikami 2 - MIKROAL i 2 - P90 (przeżycie 90 dni). Jeśli usuniemy ten efekt z modelu różnica w wartości chi-kwadrat wynosi 9,076 z 2 stopniami swobody. Jak widzimy wartość ta jest istotna na poziomie p = 0,011. Eliminacja tego efektu powoduje istotne pogorszenie się dopasowania modelu. Efekt 24 musimy więc pozostawić w modelu. Podobnie jest z efektem 14. Reasumując, do modelu powinniśmy włączyć następujące zależności:


Oczywiście interesują nas przede wszystkim czynniki, które są związane z przeżyciem (zmienna 4). Zakładamy więc, że nie interesują nas interakcje pomiędzy zmiennymi zależnymi np. między wielkością udaru a stężeniem glukozy. Dlatego też nie interesują nas powiązania pomiędzy zmiennymi niezależnymi, które nie wpływają na zmienną zależną. W praktyce zazwyczaj wyłącza się takie interakcje. Jednakże aby uniknąć braku dopasowania związanego z usunięciem (nie interesujących nas) interakcji pomiędzy zmiennymi niezależnymi do modelu włączamy efekt reprezentujący wszystkie interakcje pomiędzy zmiennymi niezależnymi. W naszym przykładzie jest to efekt 123.
Ostatecznie do modelu włączamy następujące zależności 14 24 123. Wartości te podajemy w oknie otwierającym się po kliknięciu przycisku Określenie modelu do testowania. Sytuacja taka przedstawiona jest na poniższym rysunku.


Rys. 8. Okno definicji modelu


Rys. 9. Okno wyników dla przykładu 2

Jak widzimy testy chi-kwadrat nie są istotne statystycznie (p = 0,924). Zatem możemy wyciągnąć wniosek, że określony przez nas model w sposób zadawalający wyjaśnia liczebności w tabeli. Reasumując, dotychczasowa analiza wykazała dwa istotne efekty, to jest zależności:


Zbadamy teraz naturę tych efektów. Wiemy, że dopasowanie modelu wiąże się z obliczaniem wartości oczekiwanych, tak aby odzwierciedlały one liczebności tabel brzegowych. Dlatego, dla głębszej interpretacji otrzymanych efektów przeglądniemy tabele brzegowe. Najpierw spójrzmy na zależność WIELK (wielkość ogniska) względem P90 (przeżycie 90 dni) przedstawioną w poniższej tabeli.


Rys. 10. Tabela brzegowa dla zależności WIELK względem P90

Analiza tej tabeli wykazuje, że stopa przeżycia (stosunek TAK do NIE) pacjentów, których ognisko jest mniejsze (WIELK = 1) jest prawie 7 razy większa niż dla pacjentów z dużym ogniskiem. Z kolei tabela MIKROAL (mikroalbuminuria) względem P90 (przeżycie 90 dni) wygląda następująco:


Rys. 11. Tabela brzegowa dla zależności MIKROAL względem P90

Widoczne jest, że stopa przeżycia najwyższa dla pacjentów z normą (4 do 1) zmniejsza się gwałtownie wraz ze wzrostem stężenia albuminów. Dla pacjentów u których zdiagnozowano białkomocz wartość ta jest najmniejsza i wynosi 1 do 1. Na podstawie analizy tej tabeli możemy stwierdzić, że głównymi czynnikami związanymi z przeżyciem przez pacjentów 90 dni jest wielkość ogniska udaru i stężenie albumin. Ten ostatni czynnik wydaje się być decydującym. Przykład zakończymy przeglądnięciem wykresu wartości obserwowanych względem dopasowanych pokazanym na poniższym rysunku.


Rys. 12. Wykres rozrzutu wartości obserwowanych względem dopasowanych

Jak widzimy, większość punktów na powyższym wykresie układa się na linii prostej. Świadczy to, że w tabeli nie ma "niedopasowanych" komórek zwanych wartościami odstającymi.

Przykład III zastosowanie analizy przeżycia (model proporcjonalnego hazardu Coxa)

Badano wpływ różnych czynników na czasy przeżycia w dniach (zmienna CZAS) pacjentów po niedokrwiennym udarze mózgu. Rozważano następujące czynniki:


Zebrane dane dla 56 pacjentów przedstawia poniższa tabela. Zmienna UCIĘTE to wskaźnik ucinania (0 - obserwacja kompletna, 1 - obserwacja ucięta).

Do analizy naszego zagadnienia użyjemy model proporcjonalnego hazardu Coxa. Oszacujemy współczynniki regresji dla naszych pięciu zmiennych niezależnych ze względu na przewidywane czasy przeżycia. Zmienne do analizy wybieramy jak na poniższym rysunku.


Rys. 13. Okno wyboru zmiennych dla przykładu

Następnie akceptując wszystkie ustawienia domyślne przechodzimy jak w poprzednim przykładzie do estymacji parametrów. Klikając OK uruchamiamy procedurę estymacji. Po wykonaniu 16 kroków, gdy program znalazł najlepsze parametry, procedura iteracyjna zatrzymała się. Następnie znów klikamy OK, aby przejść do pokazanego poniżej okna wyniki regresji.

Rys. 14. Okno wstępnych wyników modelu Coxa dla danych z przykładu

Jak widzimy ogólna wartość chi-kwadrat dla tego modelu jest istotna (p = 0,01051), możemy więc wnioskować, że co najmniej jedna zmienna niezależna jest istotnie związana z przeżyciem. Aby przekonać się która, klikamy przycisk Oceny parametrów otwierający arkusz wyników z wartościami estymatorów. Arkusz ten widoczny jest na poniższym rysunku.

Rys. 15. Arkusz wyników z wartościami estymatorów parametrów modelu

Na podstawie powyższego arkusza wyników wnioskujemy, że zmienne MIKROAL, HIPER, SK1 oraz FIBR są istotnymi predyktorami hazardu. Zmienna opisująca stężenie fibrynogenu w osoczu krwi jest z nich wszystkich najbardziej istotna. Z kolei zmienna opisująca stężenie cholesterolu w surowicy krwi okazała się nieistotna. Analizę naszą kończymy porównaniem wykresów funkcji przeżycia dla różnych poziomów zmiennej jakościowej HIPER. Krzywe te pokazane są na poniższym rysunku. Wykres po lewej stronie przedstawia przebieg funkcji przeżycia gdy HIPER = 1 (hiperglikemii nie ma). Wykres zaś po prawej stronie to krzywa przeżycia dla pacjentów z hiperglikemią (HIPER = 2). Pozostałe zmienne przyjmują wartości średnie. Jak widzimy krzywa przeżycia w przypadku wystąpienia hiperglikemii gwałtownie spada prawie do zera w ciągu pierwszych 15 dni.
Jest to w pełni zgodne z otrzymanymi wcześniej rezultatami. Na podstawie arkusza wyników (rys. 15) otrzymujemy bowiem, że współczynnik ryzyka dla zmiennej HIPER wynosi 3,33 Oznacza to, że w przypadku wystąpienia hiperglikemii powoduje ponad trzykrotne (3,33) zwiększenie ryzyka zgonu.


Rys. 16. Krzywe przeżycia bez hiperglikemii (u góry) i z hiperglikemią (na dole)

Przykład IV zastosowanie analizy wariancji z powtarzanymi pomiarami

Rehabilitacja jest jednym z ważniejszych - poza farmakologicznym - elementem. Ma ona na celu uzyskanie maksymalnej sprawności ruchowej chorego i samodzielności w czynnościach życia codziennego po udarze mózgu. Poniższy przykład opisuje porównanie dwóch metod rehabilitacji w zależności od wielkości ogniska udaru. Grupa 6 pacjentów którym zastosowano różne metody rehabilitacji była oceniana (w pewnej skali neurologicznej) pod względem braku usprawnienia. Pacjentów rehabilitowano dwiema metodami oznaczonymi w naszym arkuszu danych jako "Metoda1"i "Metoda2". W eksperymencie występują dwa czynniki powtarzanych pomiarów. jeden związany z kolejnymi pomiarami sprawności ruchowej i samodzielności co 3 dni (Okres1, Okres2, Okres3), drugi z różną wielkością ogniska (W1, W2, W3). Czynnik "METODA" jest czynnikiem międzygrupowym, ponieważ pacjenci zostali losowo przypisani do jednej z dwóch metod rehabilitacji.
Analizę statystyczną naszego zagadnienia przeprowadzimy w module ANOVA/MANOVA pakietu STATISTICA.
Po uruchomieniu modułu i określeniu zmiennych musimy określić czynniki powtarzanych pomiarów. Dokonujemy tego po kliknięciu na przycisku . Otworzy się wówczas okno dialogowe Określ czynniki powtarzanych pomiarów widoczne no poniższym rysunku.


Rys. 17. Okno określania czynników powtarzanych pomiarów

Jako pierwszy poziom wpisujemy najwolniej zmieniający się indeks. W naszym przypadku jest to zmienna czas. I tak dalej do najszybciej zmieniającego się indeksu, który powinien być wpisany jako ostatni. Aby poznać jak jest to określone w arkuszu danych wystarczy przesunąć palcem wzdłuż listy od początku i obserwować zmianę nazw zmiennych. W naszym przypadku wielkość ogniska jest najszybciej zmieniającym się indeksem. I dlatego określiliśmy go na samym końcu.
Pamiętajmy! W przypadku występowania więcej niż jednego czynnika powtarzanych pomiarów musimy poprawnie określić kolejność występowania czynników. Jeżeli źle to zrobimy określimy błędny układ danych i otrzymamy nieprawdziwe wyniki.
Następnie klikamy OK, aby rozpocząć analizę. Jej wyniki otrzymamy po kliknięciu przycisku Wszystkie efekty. Otworzy się wówczas arkusz wyników widoczny na poniższym rysunku.


Rys. 17. Okno określania czynników powtarzanych pomiarów

Jako pierwszy poziom wpisujemy najwolniej zmieniający się indeks. W naszym przypadku jest to zmienna czas. I tak dalej do najszybciej zmieniającego się indeksu, który powinien być wpisany jako ostatni. Aby poznać jak jest to określone w arkuszu danych wystarczy przesunąć palcem wzdłuż listy od początku i obserwować zmianę nazw zmiennych. W naszym przypadku wielkość ogniska jest najszybciej zmieniającym się indeksem. I dlatego określiliśmy go na samym końcu.
Pamiętajmy! W przypadku występowania więcej niż jednego czynnika powtarzanych pomiarów musimy poprawnie określić kolejność występowania czynników. Jeżeli źle to zrobimy określimy błędny układ danych i otrzymamy nieprawdziwe wyniki.
Następnie klikamy OK, aby rozpocząć analizę. Jej wyniki otrzymamy po kliknięciu przycisku Wszystkie efekty. Otworzy się wówczas arkusz wyników widoczny na poniższym rysunku.

Rys. 18. Arkusz wyników analizy wariancji dla powtarzanych pomiarów

Otrzymaliśmy trzy istotne efekty: efekt główny czynnika Czas, efekt główny czynnika Wielkość oraz efekt interakcji pomiędzy czynnikami Metoda i Czas. Aby przybliżyć efekt czynnika Wielkość, przywołujemy przyciskiem Średnie/wykresy wykres interakcji średnich, pokazany poniżej.


Rys. 19. Wykres średnich dla czynnika Wielkość

Wykres ten pokazuje, że badani z maksymalną wielkością ogniska wykazują mniejszą sprawność w stosunku do pacjentów o mniejszej wielkości ogniska. Podobnie badani o średniej wielkości ogniska wykazywali mniejsze usprawnienie niż pacjenci z minimalnym ogniskiem. Efekt ten odzwierciedla różnice pomiędzy wielkością ogniska uwzględnioną w eksperymencie. Podobnie możemy przeanalizować wykres dla efektu Czas i interakcji Metoda i Czas. Wykres dla tego ostatniego tworzymy, przyjmując domyślne ustawienie rozmieszczenia czynników. Otrzymamy wówczas poniższy rysunek.


Rys. 20. Wykres interakcji Metoda i Czas
Zauważamy, że w obydwu przypadkach sprawność pacjentów wzrasta, tzn. popełniają oni coraz mniej błędów. Jednakże dwie linie oznaczające różne metody rehabilitacji zaczynają się różnić poczynając od 2 okresu czasu. Wydaje się, że pacjenci w okresach 2 i 3 popełniają mniej błędów po rehabilitacji metodą drugą niż metodą pierwszą. Wyjątek stanowi początkowy okres1. Wyjaśnijmy dodatkowo interakcję poprzez przeprowadzenie analizy kontrastów (patrz poprzedni odcinek). W tym celu powracamy do okna Wyniki ANOVA i przyciskiem Porównania zaplanowane otwieramy okno do określania kontrastów. Określimy kontrasty oceniające różnicę pomiędzy metodami oddzielnie w momencie 2 i 3. W tym celu określamy najpierw kontrast dla czynnika Metoda, wprowadzając odpowiednio 1 i -1 (aby przeciwstawić pierwszą metodą drugiej). Następnie w kolejnym oknie dla czynnika czas wprowadzamy odpowiednio współczynniki 0, 1 i 0. Dla czynnika Wielkość pozostawiamy domyślne współczynniki 1, 1 i 1(czynnik nie stanowi na razie przedmiotu szczegółowego badania). Opisana sytuacja pokazana jest na poniższym rysunku.


Rys. 21. Okno określające kontrast

Klikamy teraz OK, aby otrzymać arkusz wyników pokazany niżej.

Rys. 22. Arkusz wyników istotności kontrastu

Widzimy, że różnica pomiędzy dwiema średnimi (metoda1 względem metod2) w momencie 2 nie jest statystycznie istotna. Analogicznie przy pomocy kontrastów możemy dokonać porównanie pomiędzy tymi dawkami w momencie 3 (dla czynnika Czas określamy współczynniki 0, 0, i 1). Jak się okaże, nie możemy mówić, że różne metody prowadzą do istotnej różnicy pomiędzy oceną braku usprawnienia w drugim i trzecim momencie. Jednakże, jak pokazaliśmy to wcześniej całkowita interakcja jest istotna. Wynika więc stąd, że za interakcję czynnika Metoda i Czas odpowiada głównie zróżnicowanie pomiędzy momentem 1 i 2 (krzyżujące się linie na rysunku 4). Również to stwierdzenie przetestujemy w oparciu o analizę kontrastów. Aby to wykonać budujemy następujące współczynniki kontrastu:

Kontrastujemy w ten sposób dwie metody leku oraz dwa pierwsze poziomy czynnika Czas (okres 1 i okres 2). Istotność tego kontrastu odczytujemy z poniższego arkusza wyników.

Rys. 23. Arkusz wyników istotności kontrastu oceny interakcji

Kontrast okazuje się istotny na poziomie p = 0,013546. Wynika stąd, że interakcję drugiego rzędu między czynnikami Metoda i Czas możemy przypisać różnym zmianom oceny niesprawności od momentu 1 do momentu 2 w zależności od stosowanej metody rehabilitacji. Pacjenci po zastosowaniu drugiej metody bardziej poprawili swoją sprawność niż pacjenci, których rehabilitowano metodą pierwszą.

Musimy wspomnieć o dodatkowym założeniu wynikających z powtarzania pomiarów. Mówiąc w skrócie, jednowymiarowa ANOVA z powtarzanymi pomiarami zakłada, że zmiany wzdłuż poziomów nie są skorelowane pomiędzy osobnikami. Założenie to jest mocno wątpliwe w większości przypadków. Jak to wygląda w przypadku naszych danych? Jest całkiem możliwe, że pacjenci, którzy poprawili się od okresu 1 do okresu 2, osiągnęli pułap swoich możliwości i w związku z tym poprawili się mniej od okresu 2 do okresu 3. Jak sprawdzić, czy tzw. założenie o sferyczności zostało naruszone?
Możemy wyliczyć współczynnik korelacji w miejscu w którym przypuszczamy, że założenie zostało naruszone. Na początku przy pomocy przycisku Porównania zaplanowane określamy kontrast 1 i -1 dla czynnika Metoda oraz równocześnie dwa niezależne kontrasty dla czynnika Czas. Wprowadzamy następujące współczynniki (*) -2, 1, 1 oraz (**) 0, 1, -1. Określiliśmy dwie niezależne hipotezy obejmujące efekt główny dla czynnika Czas. Pierwszy z kontrastów porównuje okres1 z połączonymi okresami 2 i 3. Drugi kontrast natomiast służy do sprawdzania okresu 2 względem okresu 3. Następnie przy pomocy przycisku Opcje wyników modyfikujemy sposób obliczania i wyświetlania wyników. Dla naszych potrzeb wybieramy opcje w sposób pokazany na poniższym rysunku. Zaznaczamy opcję Wyświetl macierze sum kwadratów (dla wyliczenia współczynnika korelacji) oraz opcje Test Greenhouse i Test jedno- i wielowymiarowe (dla skorygowania wyników w przypadku naruszenia założeń).


Rys. 24. Okno definiowania opcji wyników
Przy tak przyjętych opcjach podczas sprawdzania istotności kontrastu zostaną wyświetlone macierze sum kwadratów i iloczynów mieszanych. Macierz dla błędu będzie wyglądać jak poniżej. Umożliwia ona wyliczenie rzeczywistej korelacji pomiędzy analizowanymi kontrastami następującym równaniem:

Rys. 25. Macierz sumy kwadratów dla błędu

Jak widzimy, dwa zbiory kontrastów są skorelowane na poziomie -0,72. Ujemny znak współczynnika korelacji pokazuje, że większym zmianom od momentu 1 do połączonych momentów 2 i 3 towarzyszą mniejsze zmiany od momentu 2 do momentu 3. Potwierdza to nasze podejrzenia, że ci osobnicy, którzy wcześniej poprawili swoją sprawność, mają później mniejsze możliwości dalszej poprawy. Jak postąpić dalej w przypadku naruszenie tego podstawowego założenia? Mamy dwa wyjścia:


Rys. 26. Arkusz wyników skorygowanych testów jednowymiarowych

Jak widzimy, poprawka Greenhouse'a-Geissera dobrze zabezpiecza przed błędem przyjęcia interakcji jako istotnej. Rzeczywisty poziom istotności dla tego testu wynosi bowiem p = 0,0569 i jest bardzo zbliżony (jak przekonamy się o tym poniżej) do poziomu istotności dla testów wielowymiarowych. Poprawka Huynha-Feldta w naszym przypadku nie pomogła.

Rys. 27. Arkusz wyników testów wielowymiarowych

Występują różne rodzaje testu wielowymiarowego: lambda Wilksa, współczynnik R Rao, Ślad Pillai-Bartletta. Wszystkie bowiem stwierdzają, że interakcja nie jest istotna.

Reasumując, mamy do czynienia z przykładem, gdy naruszenie założenia sferyczności prowadzi do błędnej akceptacji statystycznej istotności interakcji. Stwierdzona początkowo interakcja czynników Metoda i Czas pod wpływem dalszej analizy (poważne naruszenie założeń) okazała się nieistotna. Przykład ten pokazał nam problemy z jakimi możemy się spotkać podczas analizy układów z powtarzanymi pomiarami i sposoby ich rozwiązania.

Powrót na początek
Poprzedni artykuł     Następny artykuł