Laboratorium EEG/Analiza zjawiska ERD/ERS

Laboratorium_EEG/Analiza zjawiska ERD/ERS

Problem wyboru optymalnego estymatora gęstości energii w przestrzeni czas-częstość dla sygnałów EEG

Dwa aspekty, które trzeba wziąć pod uwagę w praktyce analizy sygnałów EEG to:

• rozdzielczość czas-częstość oferowaną przez daną metodę;
• łatwość interpretacji wyników dla sygnałów złożonych z wielu, a priori nieznanych struktur — takim niezwykle złożonym sygnałem jest EEG.

Jeśli chodzi o rozdzielczość czasowo-częstościową estymatorów gęstości energii sygnałów to wiadomo, że jest ona ograniczona przez zasadę nieoznaczoności Heisenberga w analizie sygnałów.

Minimum iloczynu wariancji w dziedzinie czasu i wariancji w dziedzinie częstości, osiągane jest przez dystrybucję Wigner-Villa dla sygnału będącego funkcją Gabora. Dystrybucja Wigner-Villa jest formą kwadratową. W przypadku bardziej złożonych sygnałów własność ta prowadzi do powstawania wyrazów mieszanych, które utrudniają prawidłową interpretację wyników. Dystrybucje wywodzące się z dystrybucji Wigner-Villa tworzą klasę dystrybucji Cohena. W klasie tych dystrybucji do tłumienia wyrazów mieszanych stosowane jest jądro, będące efektywnie czasowo-częstościowym filtrem dolnoprzepustowym. Często stosowanym w czasowo-częstościowej analizie sygnałów reprezentantem klasy dystrybucji Cohena jest dystrybucja Choi-Williamsa, gdzie jądrem filtru jest dwuwymiarowa funkcja Gaussa. Odpowiedni dobór własności jądra umożliwia stłumienie niektórych wyrazów mieszanych --- metoda ta może być szczególnie efektywna jeśli znane są a priori czasowo-częstościowe własności struktur obecnych w sygnale. Należy także pamiętać, że filtrowanie prowadzi do zmniejszenia rozdzielczości czas-częstość. Alternatywą dla dystrybucji czas-częstość są metody oparte o liniowe reprezentacje czas-częstość takie jak krótkoczasowa transformata Fouriera (spektrogram) czy transformata falkowa (skalogram). W metodach tych sygnał analizowany jest przy pomocy funkcji zlokalizowanych w czasie i częstości. Zastosowanie funkcji zlokalizowanych częściowo zmniejsza problem wyrazów mieszanych, jednak ich rozdzielczość czas-częstość jest związana z rozciągłością czasowo-częstościową stosowanych funkcji. Metody te nie zapewniają dopasowania do struktur obecnych w sygnale.

W pracy H1 zademonstrowano, że najbardziej optymalne dla analizy sygnałów EEG własności posiada metoda oparta o dekompozycję sygnału algorytmem Matching Pursuit (MP). Ten estymator gęstości energii konstruowany jest jako suma dystrybucji Wigner-Villa poszczególnych składowych wyodrębnionych z sygnału w trakcie dekompozycji. W sposób wynikający bezpośrednio z konstrukcji, estymator ten nie posiada wyrazów mieszanych oferując jednocześnie najwyższą możliwą rozdzielczość czasowo-częstościową. Kompromis pomiędzy rozdzielczością w czasie a rozdzielczością w częstości jest w każdej iteracji algorytmu MP optymalizowany ze względu na struktury występujące w sygnale. W tym sensie metoda ta jest lokalnie adaptywną reprezentacją czas-częstość, w odróżnieniu od krótkoczasowej transformaty Fouriera czy transformat falkowych.

Klasyczna technika pomiaru efektów ERD/ERS opisana w naszym podręczniku bazuje na uśrednianiu energii sygnałów przefiltrowanych w wybranym paśmie częstości. Wymaga ona wyboru reaktywnych (dających największe zmiany mocy) pasm częstości. Pasm takich zwykle poszukiwano metodą prób i błędów. Ze względu na dużą zmienność osobniczą pasm reaktywnych, znacznie lepszym podejściem, które umożliwia całościowe spojrzenie na efekty modyfikacji widma mocy sygnału EEG, jest analiza czasowo-częstościowa. Polega ona na estymowaniu gęstości energii w przestrzeni czas-częstość dla każdego powtórzenia osobno. Następnie dla każdej częstości, obliczana jest miara ERD/ERS, przy pomocy równania

[math]\mathrm{ERD/ERS}_f(t) = \frac{P_f(t) - R_f}{R_f} 100\% [/math]

gdzie:

• [math]P_f(t)[/math] - średnia moc chwilowa w paśmie częstości [math]f[/math]. Uśrednianie przebiega po realizacjach.
• [math]R_f[/math] - średnia moc w paśmie [math]f[/math] w okresie referencyjnym. Uśrednianie przebiega po realizacjach i po czasie trwania okresu referencyjnego.

W efekcie uzyskiwana jest mapa prezentująca zjawiska ERD/ERS w przestrzeni czas-częstość.

Najczęściej stosowanymi do szacowania ERD/ERS estymatorami gęstości energii są: spektrogram (Makeig, 1993), skalogram (Tallon, 1996), dystrybucje klasy Cohena (Lachaux, 2000) lub filtrowanie pasmowe w zachodzących na siebie pasmach częstości (Graimann, 2002). W sytuacjach gdy istotna jest wysoka rozdzielczość czasowo-częstościowa zaproponowana została metoda Matching Pursuit (Durka, 2001). Raportowane wyniki zależą w pewnym stopniu od własności sygnału, ale również w pewnym stopniu od zastosowanego estymatora i jego parametrów.

Prace innych autorów zwykle skupiały się na zastosowaniach jednego konkretnego estymatora. Nie dyskutowano w nich jaki wpływ na identyfikację efektów reakcji na bodziec ma wybór metody estymacji gęstości energii. Problem ten został zbadany w pracach H2 i H3. Wszystkie badane estymatory zostały zastosowane do tych samych zestawów danych co umożliwiło praktyczne porównanie wpływu estymatora na czytelność i jednoznaczność wyników. W pracy H2 pokazano, że algorytm MP daje estymaty o wyższej rozdzielczości niż spektrogram, niemniej w kontekście identyfikacji obszarów przestrzeni czas-częstość, w których przejawia się reakcja związana z bodźcem, wyniki otrzymane przy pomocy obu estymatorów są podobne i spójne. W pracy H3 zagadnienie badania optymalności estymatora zostało rozszerzone w dwóch aspektach. Do porównania włączono skalogram zaś porównanie metod przeprowadzono zarówno dla danych EEG jak i ECoG. Wykazane zostało, że wszystkie wymienione powyżej metody estymacji gęstości energii w dziedzinie czas-częstość dają spójne wyniki. Wykazano, że estymator oparty o dekompozycję MP umożliwia badanie mikrostruktury zmian gęstości energii dzięki najwyższej rozdzielczości czas-częstość. Jednocześnie, dzięki lokalnie adaptywnemu kompromisowi pomiędzy rozdzielczością w czasie a rozdzielczością w częstości daje on wyniki, w których eksponowana jest struktura sygnału nie obarczona w sposób systematyczny własnościami rozdzielczości czasowo-częstościowej metody. W przypadku pozostałych metod, przy interpretacji otrzymanych map czas-częstość należy brać pod uwagę specyficzną dla danej metody strukturę wyrazów mieszanych i zależność pomiędzy rozdzielczością w czasie a rozdzielczością w częstości (porównaj H3, Rys.3).

Problem istotności statystycznej efektów ERD i ERS w przestrzeni czas-częstość

W badaniach naukowych bardzo ważnym elementem jest ocena istotności statystycznej obserwowanych zjawisk. Dążeniem naszym jest aby identyfikowane zmiany mocy sygnału był efektem związanym z bodźcem, a nie jedynie artefaktem fluktuacji statystycznych. Problem prawidłowej oceny istotności statystycznej zmian w rozkładach gęstości energii w dziedzinie czas-częstość został opisany w serii dwóch artykułów H2 i H3. Zaproponowane rozwiązanie opiera się na masowym zastosowaniu testów jednowymiarowych i korygowaniu wyników uwzględniając wielokrotność testowania. Zbadano dwa zagadnienia:

• wybór statystyki
• korekty poziomu istotności testów ze względu na wielokrotność porównań.

Wcześniejsze prace innych autorów nie wskazywały poprawnego metodologicznie rozwiązania zagadnienia wielokrotnych porównań, które stanowi tu główny problem. Przykładowo Makeig (1993) jako rozwiązanie problemu wielokrotnych porównań proponuje przyjęcie poziomu istotności pojedynczego testu jako 0.001. W pracy (Lachaux, 2000) istotność wzrostu mocy była szacowana dla pojedynczego pasma częstości, problem wielokrotnych porównań był wspomniany, ale nie rozwiązany w przypadku porównywania aktywności wielu pasmach jednocześnie. W pracy (Graimann, 2002) szacowano istotność statystyczną w poszczególnych pasmach niezależnie, a następnie prezentowano mapę stworzoną z tych pasm nie korygując poziomów istotności.

Przyczyny takiego stanu rzeczy można upatrywać w fakcie, że powszechnie stosowaną poprawką na poziom błędu I rodzaju ze względu na wielokrotność porównań jest poprawka Bonferroniego. Poprawka ta daje prawidłową korektę poziomu istotności dla testów niezależnych. W przypadku testów skorelowanych poprawka ta jest zbyt konserwatywna. W analizie czasowo-częstościowej mamy do czynienia z bardzo dużą ilością testów, które są ze sobą skorelowane w a priori nieznanym stopniu. Zastosowanie tu poprawki Bonferroniego prowadzi do testów skrajnie konserwatywnych. Jako rozwiązanie tego problemu w pracy H2 została zaproponowana rezygnacja z kontroli błędu I rodzaju na rzecz kontroli frakcji fałszywych odkryć (false discovery rate, FDR) (Benjamini & Yekutieli, 2001). W kolejnych paragrafach omówione zostaną istotne elementy rozwiązań zaproponowanych w pracach H2 i H3.

Podział przestrzeni czas-częstość.

W pracy H2 stwierdzono, że aby móc wygodnie sformułować hipotezy i poddać je testom należy podzielić przestrzeń czas-częstość na elementy rozdzielczościowe (ang. resolution element---resel). Iloczyn rozciągłości czasowej i częstościowej jest od dołu ograniczony przez zasadę nieoznaczoności. W praktyce, co zostało zbadane w pracy H2, rozmiar resela musi być większy ze względu na wariancję estymatorów i wariancję danych pomiędzy realizacjami.

Hipotezy.

Stosując opisany powyżej podział przestrzeni czas-częstość na resele można sformułować następującą rodzinę hipotez: dla każdego resela hipoteza zerowa stwierdza, że średnia energia rozważanego resela równa jest średniej energii reseli o tej samej częstości należących do okresu referencyjnego, hipoteza alternatywna stwierdza, że średnia energia w okresie referencyjnym i w rozważanym reselu są różne. Uśrednianie przebiega po realizacjach.

Wybór statystyki.

W pracy H2 zademonstrowano, że rozkłady energii w reselach są zależne od zastosowanego estymatora gęstości energii, ale w większości przypadków są dalekie od rozkładów normalnych. Wykazano, że zastosowanie intensywnych obliczeniowo testów nieparametrycznych prowadzi do poprawnych i nieobciążonych wyników. Spośród dwóch zbadanych testów bardziej efektywnym okazał się test oparty na estymowaniu empirycznego rozkładu statystyki pseudo-t w okresie referencyjnym.

W pracy H3 wykazano, że w wielu praktycznych przypadkach poprawne rezultaty można uzyskać stosując odpowiednio testy parametryczne. Pierwszym koniecznym krokiem jest normalizacja danych, czyli poddanie energii reseli takiej transformacji, że rozkład otrzymanych danych jest dobrze przybliżony rozkładem normalnym. Zaproponowano w tym celu zastosowanie transformacji Boxa-Coxa (Box & Cox, 1964). Kolejną cechą rozważanych testów jest niejednorodność wariancji w okresie referencyjnym i pobodźcowym. Teoretycznym rozwiązaniem tego problemu jest test t ze zmodyfikowaną ilością stopni swobody zgodnie z poprawką Welcha (Welch, 1938). W pracy H3 pokazano, że dla danych testowych zarówno EEG (Rys. 1) jak i ECoG (Rys. 2) standardowy test t jak i test t z poprawką Welcha daje wyniki niemal identyczne z testem nieparametrycznym zaproponowanym w H2. Zastosowanie testów parametrycznych znacząco redukuje złożoność obliczeniową problemu oceny istotności statystycznej zmian mocy sygnału w dziedzinie czas-częstość.

Problem wielokrotnych porównań.

Ocena istotności statystycznej dla wszystkich analizowanych reseli w dziedzinie czas-częstość prowadzi w naturalny sposób do powstania problemu wielokrotnych porównań. Jak wspomniano we wstępie do tej sekcji, we wcześniejszych pracach dotyczących analizy odpowiedzi fazowo nie związanej z bodźcem w dziedzinie czas-częstość nie zaproponowano poprawnego rozwiązania problemu wielokrotnych porównań. W pracy H2 jako efektywne i ścisłe rozwiązanie problemu wielokrotnych porównań została zaproponowana metoda kontroli frakcji fałszywych odkryć (FDR). Rozważono założenia i porównano wyniki otrzymywane dla klasycznej metody kontrolowania błędu I rodzaju dla problemu wielokrotnych porównań Bonferroniego-Holmsa z metodą opartą o kontrolę FDR. Aby praktycznie zilustrować efektywność i skuteczność metody kontroli błędów w oparciu o podejście FDR zaprojektowano i przeprowadzono eksperyment, w którym celowo wydłużono czas rejestracji sygnału EEG, tak aby zawierał on znaczący udział zapisu spoczynkowego, w którym a priori nie powinny występować efekty związane z bodźcem. Otrzymano poprawną kontrolę ilości fałszywych detekcji zachowując jednocześnie moc testów wystarczającą do potwierdzenia istotności znanych efektów ERD/ERS. Ilustracją tego wyniku są Rys. 4 i 5 z pracy H2. Z przedstawionych w pracy H2 wyników wyciągnięto wniosek, że w przypadku wielokrotnych testów wykonywanych dla reseli w dziedzinie czas-częstość procedura Bonferroniego-Holmsa jest zbyt konserwatywna, natomiast procedura FDR zachowuje wystarczającą moc aby być w praktyce przydatną.

Zadania

• Proszę zaimplementować technikę szacowania istotności statystycznej zmian ERD/ERS opisaną w pracy H3 :
  • wersja z testem t i poprawkami Welcha na nierówność wariancji w okresie referencyjnym i badanym, oraz FDR na wielokrotność porównań
  • zastosować ją do danych z zadania zadania 3 z poprzedniego działu. Mapy czas-częstość proszę szacować techniką spektrogramu o odpowiednio dobranym okienku. Przed uśrednianiem mapy proszę zlogarytmować.

Literatura

H1

Blinowska, K. J., Durka, P. J., Żygierewicz, J., 2004. Time-frequency analysis of brain electrical activity–adaptive approximations. Methods Inf Med 43, 70–73.

H2

Durka, P. J., Żygierewicz, J., Klekowicz, H., Ginter, J., Blinowska, K. J., 2004. On the statistical significance of event-related EEG desynchroniza- tion and synchronization in the time-frequency plane. IEEE Trans Biomed Eng 51, 1167–1175.

H3

Żygierewicz, J., Durka, P. J., Klekowicz, H., Franaszczuk, P. J., Crone, N. E., 2005. Computationally efficient approaches to calculating significant ERD/ERS changes in the time-frequency plane. J Neurosci Methods 145, 267–276. [[1]]

Benjamini, Y., Yekutieli, Y., 2001. The control of the false discovery rate under dependency. Annals of Statistics 29, 1165–1188.

Box, G. E. P., Cox, D. R., 1964. An analysis of transformations. Journal of the Royal Statistical Society 2, 211–252.

Durka, P. J., Ircha, D., Neuper, C., Pfurtscheller, G., 2001. Time-frequency mi- crostructure of event-related desynchronization and synchronization. Medical and Biological Engineering and Computing 39 (3), 315–321.

Graimann, B., Huggins, J. E., Levine, S. P., Pfurtscheller, G., 2002. Visualiza- tion of significant ERD/ERS patterns in multichannel EEG and ECoG data. Clinical Neurophysiology 113, 43–47.

Lachaux, J.-P., Rodriguez, E., Martinerie, J., Adam, C., Hasboun, D., Varela, F. J., 2000. A quantitative study of gamma-band activity in human intracra- nial recordings triggered by visual stimuli. European Journal of Neuroscience 12, 2608–2622.

Makeig, S., 1993. Auditory event-related dynamics of the EEG spectrum and effects of exposure to tones. Electroencephalography and Clinical Neurophy- siology 86, 283–293.

Tallon-Baudry, C., Bertrand, O., Delpuech, C., Pernier, J., Jul 1996. Stimulus specificity of phase-locked and non-phase-locked 40 Hz visual responses in human. The Journal of Neuroscience 16, 4240–4249.

Welch, B., 1938. The significance of the difference between two means when the population variances are unequal. Biometrika 29, 350–362.

Laboratorium_EEG/Analiza zjawiska ERD/ERS