Laboratorium EEG/EEGLAB

Z Brain-wiki

Laboratorium_EEG/EEGLAB

Zajęcia bazują na podręczniku praktycznym użytkowania pakietu EEGLAB dostępnym na stronach:

Analiza potencjałów wywołanych

Wprowadzenie

W tej części zajęć naszym celem jest przećwiczenie pracy z danymi z eksperymentu, w którym badamy potencjały wywołane na grupie osób. Dane pochodzą z eksperymentu opisanego w artykule: Effects of Valence and Origin of Emotions in Word Processing Evidenced by Event Related Potential Correlates in a Lexical Decision Task

Skupimy się na badaniu różnic pomiędzy przetwarzaniem słów i pseudo-słów.

Proszę też zapoznać się z opisem stylu raportowania wyników w psychologii

Przebieg eksperymentu

Uczestnicy zostali poinformowani o celu eksperymentu i charakterze pomiaru EEG. Zachęciliśmy ich, aby utrzymywali wygodną postawę i kontrolowali mruganie. Protokół zapewniał 3-s przerwy na normalne mruganie co 10 prób, jak również dwie dłuższe przerwy, których czas trwania jest kontrolowany przez uczestnika, dla odpoczynku i poprawy postawy. Długie przerwy występowały co 270 prób.

Zadaniem było odczytanie bodźców, które pojawiły się na środku ekranu, i zaklasyfikowanie ich jako słów lub pseudo-słów przez naciśnięcie oznaczonych klawiszy na klawiaturze. Rejestrowano typ i opóźnienie odpowiedzi. Pojedynczy blok eksperymentalny zawierał 135 słów i 135 pseudo-słów; ten blok powtórzono trzy razy. Słowa i pseudo słowa były wyświetlane w losowej kolejności we wszystkich blokach. Próby przebiegały w następujący sposób: (1) punkt fiksacji wyświetlany przez 500 ms; (2) bodziec wyświetlany, dopóki uczestnik nie odpowie; (3) pusty ekran wyświetlany przez losowo zmienny odstęp czasu między 1000 a 1100 ms.

Przygotujmy katalogi do analizy. Przydadzą się następujące:

  • DANE
  • SKRYPTY
  • RYSUNKI

Dane dostępne są tutaj. Proszę je ściągnąć i wypakować do katalogu DANE. Powinno pojawić się 66 plików, po 2 na każdą osobę. Pliki *.set to pliki matlabowe zawierające struktury EEG, pliki *fdt zawierają właściwe sygnały.


Dwie funkcje pomocnicze, które będą potrzebne proszę zapisać w katalogu SKRYPTY:

SKRYPTY

Etapy analizy, które musimy wykonać to:

  • Filtrowanie i cięcie sygnałów
  • Dobór progów dla detekcji artefaktów na podstawie oglądania sygnału i prób z różnymi progami do automatu (abnormal values / abnorlmal trends)
  • Automatyczne usuwanie artefaktów
  • Przygotowanie danych do analiz wg. przynależności do czynnika grupującego
  • Wybór obszarów zainteresowania (grup elektrod do analizy)
  • Analiza wariancji
  • Zilustrowanie wyników analizy

Aby nasze analizy były łatwo powtarzalne przygotujemy główny skrypt np: run_me.m, do którego będziemy dodawać stopniowo skrypty i funkcje kolejnych etapów analizy. Dla czytelności kodu i możliwości uruchamiania osobno poszczególnych etapów posłużymy się technologią komórek w Matlabie (komórkę rozpoczynają dwa procenty '%%');

Zacznijmy zatem pisać run_me.m. Ustawiamy zmienne ze ścieżkami, wybieramy, które dane wchodzą do analizy (część została odrzucona ze względu na nadmiar artefaktów), częstości odcięcia filtrów dolno- i górno-przepustowego, czasy wokół momentu wyświetlenia bodźca, informacja, które kanały zawierają EEG, zakres czasu do korekty poziomu "0": 

%% USTAWIENIA:
directoryIn  = '../DANE/';                       % lokalizacja folderu z surowymi danymi
directoryOut = '../DANE/';                       % lokalizacja do zapisu przygotowanych danych
addName      = '_epoch';                        % dodatek do nazwy pliku przy zapisywaniu przygotowanych danych
subjects     = [1:6,8:11 ,15:24,26:37];         % lista badanych do wczytania (z poprawnymi danymi)
LowFreqCut   = 0.1;                             % granica odciecia filtru gornoprzepustowego
HighFreqCut  = 30;                              % granica odciecia filtru dolnoprzepustowego
startT       = -0.3;                            % start do wyciecia epok [s]
stopT        = 1.0;                             % stop do wyciecia epok [s]
eegChan      = 1:19;                            % numery kanalow EEG
baseline     = [-200, 0];                       % baseline [od,do] [ms] do odjecia


%% FILTROWANIE i CIECIE


eeglab
filtruj_i_tnij(EEG , ALLEEG , CURRENTSET, directoryIn, directoryOut, addName, subjects, LowFreqCut, HighFreqCut, startT, stopT, eegChan, baseline)

filtruj_i_tnij.m

Musimy teraz napisać funkcję filtruj_i_tnij. Musi ona kolejno:

  • wczytać nasze zbiory danych,
  • stworzyć montaż kanałów do średniej ze wszystkich kanałów, tzw. common average reference (CAR),
  • przefiltrować je kanał po kanale,
  • wyciąć interesujące nas epoki,
  • wykonać dla nich korektę poziomu "0",
  • zapisać tak przetworzone dane do nowych "*.set"-ów.

Szkielet tej funkcji wygląda następująco, proszę uzupełnić brakujące fragmenty: 

function filtruj_i_tnij(EEG, ALLEEG, CURRENTSET, directoryIn, directoryOut, addName, subjects, LowFreqCut, HighFreqCut, startT, stopT, eegChan, baseline)
    
  % zapamiętujemy stan wejściowy zmiennych globalnych eeglaba
    EEG0 = EEG; 
    ALLEEG0 = ALLEEG; 
    CURRENTSET0 = CURRENTSET;
    events = {'Pseudo','ANeg','ANeu','APos','ONeg','ONeu','OPos','RNeg','RNeu','RPos'};   % typ eventow do wyciecia
    
    for sub = subjects
        
        % tworzymy nazwę setu dla osoby z numerem 'sub'
        if sub < 10
            nameOfSubj = ['B0' num2str(sub)];
        else
            nameOfSubj = ['B' num2str(sub)];
        end

        % wczytujemy set z danymi
        EEG = pop_loadset('filename',[nameOfSubj '.set'],'filepath',directoryIn);

        % uaktualniamy zestaw danych globalnych dla eeglaba:
       [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 );
        EEG = eeg_checkset( EEG);
        
        % zostawiamy tylko potrzebne kanaly EEG , te które są wymienione w zmiennej eegChan:
        EEG.data = EEG.data(....); % wycianmy dane
        EEG.nbchan = .... ; % uaktualniamy liczbę kanałów
        EEG.chanlocs = EEG.chanlocs(1,...);  % obcianamy listę położeń kanałów
        
        % referencja - odjecie sredniej ze wszystkich kanalow
        EEG = pop_reref(.... );
        EEG = eeg_checkset(EEG);

        % filtracja danych - zastosujemy filtr Butterwartha  2-go rzędu
        [b1,a1] = butter(2, .... , 'high');
        [b2,a2] = butter(2, ... , 'low');
        [b3,a3] = butter(2, 2*[49.5 50.5]/EEG.srate, 'stop'); % filtr pasmowozaporowy  "notch" na sieć
        for channel = 1:size(EEG.data,1)
            EEG.data(channel,:) = ... % górnoprzepustowy
            EEG.data(channel,:) = ... % dolno
            EEG.data(channel,:) = ... % notch
        end
        EEG = eeg_checkset( EEG );
        
        % ekstrakcja epok: 
        % zauważmy, że podajemy w tej funkcji listę eventów wokół których mają być robione wycinki i okresy wycinków
        % przy okazji nadajemy setowi nazwę z dodanym stringiem addName: 
        EEG = pop_epoch( EEG, events, [startT, stopT], 'newname', [nameOfSubj addName], 'epochinfo', 'yes');

        [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'gui','off'); 
        EEG = eeg_checkset( EEG );
        
        % usuniecie baseline
        EEG = pop_rmbase(...);
        EEG = eeg_checkset( EEG );
        
        % zapisanie setu
        EEG = pop_saveset( EEG, 'filename',[nameOfSubj addName],'filepath',directoryOut);
        disp(['filtrowanie i cięcie ----------- dane ' nameOfSubj ' :   GOTOWE ----------------'])
        
    end

    % przywracamy  stan wejściowy zmiennych globalnych eeglaba
    EEG = EEG0; 
    ALLEEG = ALLEEG0; 
    CURRENTSET = CURRENTSET0;
end

DOBÓR PROGOW DLA ARTEFAKTÓW

Usuwanie artefaktów stanowi jeden z bardziej żmudnych ale niezbędnych etapów przygotowania danych do dalszej analizy. Jeśli mamy dostatecznie dużo danych to najbardziej bezpiecznym podejściem jest usunięcie tych epok (spośród wcześniej wyciętych), które mają artefakty. Do najczęstszych i najbardziej szkodliwych w analizie potencjałów wywołanych należą:

  • mrugnięcia
  • płynięcie sygnału
  • mięśnie

Zwykle usuwanie artefaktów wykonujemy dwuetapowo:

  • Dla ustalonych progów na wykrycie artefaktu oznaczamy automatycznie złe fragmenty
  • Dokonujemy wizualnego przeglądu tak oznaczonych danych i wykonujemy ewentualne korekty, w miejscach nietypowych zaburzeń.

Tu dla oszczędzenia czasu zrobimy tylko procedurę uproszczoną:

UWAGA: "The are several EEGLAB plugins and legacy EEGLAB menus to reject bad data and bad channels. Use menu item File → Preference and check the checkbox to If set, show all menu items from previous EEGLAB versions. Restart EEGLAB for this change to take effect. A collection of menu items for rejecting bad portions of data become available. These involve Tools → Automatic channel rejection (see the help message of the pop_rejchan.m function) and Tools → Automatic continuous rejection (see help of the pop_rejcont.m function), which were the default methods used in previous versions of EEGLAB. A collection of menus to reject bad data epochs is also described in this section of the tutorial."

Naszym celem jest takie dobranie progów

  • abnormalValue: wartosc [uV] powyżej/poniżej(na minusie) której sygnał uznawany jest za artefakt
  • abnormalTrend: maksymalne nachylenie trendu [uV/epokę], R squared limit (0-1)

Jak już zdecydujemy się na konkretne wartości to ...

Wracamy na chwilę do run_me.m: Automatyczne usuwanie artefaktów

Zaaplikujmy ustalone przed chwilą parametry do procedury automatycznego ich usuwania. W pliku run_me.m dopisujemy następujące linijki: 

%% Automatyczne usowanie artefaktow
EEG = pop_loadset('filename', 'B01_epoch.set', 'filepath', directoryIn); % wczytujemy przykładowy  set z danymi, aby struktury EEG miały sensowne wartości
abnormalValue = ...                            % wartosc [uV] powyzej/ponizej(na minusie) ktorej oznaczane sa artefakty
abnormalTrend = ...                     % [max slope [uV/epoch], R squared limit (0-1)]
oznaczArtefakty(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects, abnormalValue, abnormalTrend)

Funkcja: oznaczArtefakty.m

W pętli po zbiorach danych aplikujemy ustalone progi: 

function [] = oznaczArtefakty(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects, abnormalValue, abnormalTrend)

      % zapamiętujemy stan wejściowy zmiennych globalnych eeglaba
       ... 
    for sub = ...
        % składamy nazwę pliku
          ....

        % wczytujemy plik:
        EEG = ...
        [ALLEEG, EEG, CURRENTSET] = ...
        EEG = eeg_checkset( EEG);
        
        % automatyczne oznaczanie artefaktow
               
        EEG = pop_eegthresh(...); % abnormal value
        EEG = eeg_checkset(EEG);
        EEG = pop_rejtrend(... ); % abnormalTrend
        EEG = eeg_checkset(EEG);
        close(gcf)
        
        % zapisanie setu
        EEG = pop_saveset(EEG, 'filename', [nameOfSubj addName], 'filepath', directoryOut);
        disp(['oznaczArtefakt ----------- dane ' nameOfSubj ' :   GOTOWE ----------------'])
        
    end
    % przywracamy stan wejściowy zmiennych globalnych eeglaba
       ... 
end

Wracamy na chwilę do run_me.m: zebranie danych

Dalsze analizy będzie wygodniej prowadzić gdy dane od wszystkich osób zbierzemy w jednej macierzy o rozmiarze [kategoria x osoby x kanaly x probki]. Będziemy w niej przechowywać dane dla każdej z osób uśrednione po realizacjach należących do tej samej kategorii, tzn. 2 kategorie: (1 = Words, 2 = Pseudo).

W run_me.m dopisujemy: 

%% przygotowanie danych

data = preparujDane(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects);


Funkcja: preparujDane.m

Tym razem wczytując dane musimy iterować się po kolejnych zbiorach danych 'subIdx' abstrachując od numeru osoby w nazwie pliku. Chcemy aby w powstającej macierzy 'data' nie było zerowych wpisów z powodu odrzucenia (braku pliku) jakieś osoby. Proponuję po uzupełnieniu kodu zatrzymać się debugerm wewnątrz pętli sumującej triale i upewnić, że dokładnie rozumiemy jak są agregowane dane. 


function [data] = preparujDane(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects)
 % zapamiętujemy stan wejściowy zmiennych globalnych eeglaba
       ... 
    data    = zeros(...);       % macierz z danymi o rozmiarze [2(1.Words|2.Pseudo) x subjects x kanaly x probki]
    trialNb = zeros(...);                                         % macierz z informacja o liczbie triali na warunek [2(1.Words|2.Pseudo) x subjects]
                    

    for subIdx = 1:length(subjects)
        
        %    składamy nazwę pliku
        if subjects(subIdx) < 10
            nameOfSubj = ['B0' num2str(subjects(subIdx))];
        else
            nameOfSubj =...
        end

        % wczytujemy plik:
        EEG = ...
        [ALLEEG, EEG, CURRENTSET] = ...
        EEG = eeg_checkset( EEG);
        
        % dla danej osoby iterujemy się po trialach
    
        for trial = 1:EEG.trials
            if ~(EEG.reject.rejconst(1,trial) || EEG.reject.rejthresh(1,trial))   % wybieramy tylko triale bez artefaktow 
                type = strcmp(EEG.epoch(1,trial).eventtype,'Pseudo')+1; % określamy kategorię trialu: 
                                                         % dla triali  Pseudo, strcmp da wartość True, czyli 1 więc w wyniku dostaniemy 2,
                                                         % cała reszta to słowa i dla nich wynik porównania strcmp da fałsz , czyli 0 
                data(type,subIdx,:,:) = squeeze(data(type,subIdx,:,:)) + squeeze(EEG.data(:,:,trial)); % i sumujemy ich 
                                                         % przebiegi dla każdej z dwoch kategorii: Words/Pseudo, żeby wymiary
                                                         % dodawanych macierzy się zgadzały musimy wpierw usunąć wymiary singletowe
                trialNb(type,subIdx) = trialNb(type,subIdx) + 1;
            end
        end
        
        % dzielimy przebiegi przez zliczenia triali w kategorii aby uzyskac sredni przebieg
        data(1,subIdx,:,:) = data(1,subIdx,:,:)/squeeze(trialNb(1,subIdx));
        data(2,subIdx,:,:) = data(2,subIdx,:,:)/squeeze(trialNb(2,subIdx));
        disp(['preparujDane ----------- dane ' nameOfSubj ' :   GOTOWE ----------------'])
        
    end
    % przywracamy stan wejściowy zmiennych globalnych eeglaba
       ... 
    
end

Wybór okien czasowych dla komponentów ERP: metoda GFP

Wybór okienek do analizy ERP można oprzeć o krzywą globalnej mocy pola (ang. global field power GFP). GFP jest obliczany jako przestrzenne odchylenie standardowe i określa sumę aktywności elektrycznej wszystkich elektrod w danym punkcie czasowym. Opóźnienia maksymalnych wartości GFP wskazują na opóźnienia potencjalnych składników potencjału wywołanego (Lehmann i Skrandies, 1980; Skrandies, 1990). Intuicja, która za tym stoi mówi, że dla pewnego stanu aktywności mamy zróżnicowaną topografię amplitudy (rozkład dipolowy , kwadrupolowy itp. ) i przejście do innego stanu powinno wieść przez brak takiej różnorodności.

Do pliku run_me.m dopisujemy: 

%% GFP1

timeMarks = [...]; %wpiszmy tu czasy z artykułu
plotGFP(data,startT,stopT,EEG.srate,timeMarks,EEG.chanlocs) 

%% GFP2 - ewentualnie można  powyższe czasy można poprawić


Funkcja: plotGFP.m

Uzupełnijmy teraz funkcję plotGFP.m: 

function plotGFP(data,startT,stopT,Fs,timeMarks,chanlocs)
    


   % Najważniejsze obliczenie w tej funkcji: GFP
    GFP_dla_osob = ... ; % liczymy odchylenie standardowe po kanałach 
    GFP = ...;  % uśredniamy po osobach 
    
    
    % Poniżej są kody dotyczące rysowania GFP
    % oś czasu
    time = (startT:1/Fs:stopT);
      
    % ustawienia do rysowania
    fontSize  = 14;
    colors    = [[1 0 0];[0 0 1]];
    colorMark = [0.5 0.5 0.5];
    marginesX = 0.1;
    marginesY = 0.1;
    labels    = {'Words','Pseudo'};
    
    
   
    % obliczanie położeń poszczególnych części rysunku:
    positionPlot = [marginesX, 1-marginesY-(1-2*marginesY)*1/2, 1-2*marginesX, (1-2*marginesY)*1/2];
    positionTopo = zeros(length(timeMarks)-1,4); 
    positionTopo(:,2:4) = repmat([1-marginesY-(1-2*marginesY)*3/4, (1-2*marginesX)/length(timeMarks), (1-2*marginesY)*1/6],(length(timeMarks)-1),1);
    for tim = 1:(length(timeMarks)-1)
        positionTopo(tim,1) = marginesX+(tim-1)*(1-2*marginesX)/length(timeMarks)+(tim-1)*(1-2*marginesX)/length(timeMarks)/(length(timeMarks)-2);
    end
    positionTopo2 = positionTopo;
    positionTopo2(:,2) = 1-marginesY-(1-2*marginesY);
    positionCbar = [1-1.1*marginesX, 1-marginesY-(1-2*marginesY)*7/8, marginesX/4, (1-2*marginesY)*1/6];
    
    
    figure()
    % rysowanie krzywej GFP i pionowych linii przedzialow czasowych
    axes('position',positionPlot)
    for type=1:size(data,1)
        hold on
        plot(time,GFP(type,:)','Color',colors(type,:))
    end 
 
     % dodajemy linie pionowe oddzielające okienka czasowe
     % dobieramy pionowy zakres zmienności GFP aby wiedzieć jak wysokie mają być linie
       YLIM = [floor(min(min(GFP))), ceil(max(max(GFP)))];
        for l=1:length(timeMarks)
            line([timeMarks(l)/1000,timeMarks(l)/1000], [YLIM(1) YLIM(2)*0.9], 'Color', colorMark)
            legend off
            hold on
        end
    xlim([startT,stopT])
    ylim(YLIM)

    % jeszcze trochę kosmetyki:
    set(gca,'box','off');
    set(gca, 'FontSize', fontSize);

    drawaxis(gca,'x',0,'y',0); % ta funkcja przesuwa osie, tak aby przechodziły przez punkt 0,0

    axis off
    
    % opisy osi 
    text(stopT*0.9,-0.1*YLIM(1)-0.3,'[ms]','FontSize',fontSize,'Interpreter','Latex')
    text(startT*0.5,0.8*YLIM(2),'[\textit{$\mu$}V]','FontSize',fontSize,'Interpreter','Latex')
    legend(labels,'Box','off')
    

   % tu jeszcze dorysowujemy główki z rozkładem topograficznym uśrednionego potencjału w ramach danego odcinka czasu i po osobach badanych, 
   
    % obliczamy minimum i maksimum po wszystkich danych uśrednianych po osobach (czyli pozostaje nam zmienność po: warunkach, kanałach i czasie)
    % po to aby wszystkie główki miały ten sam zakres kolorów i były porównywalne:
    scaleTopo = [min(min(min(mean(data,2)))) max(max(max(mean(data,2))))];

    if ~isempty(timeMarks)
        % rysowanie glowek topo dla przedzialow czasowych
        for l=1:(length(timeMarks)-1)
            axes('position',positionTopo(l,:))
            idx = ... % wybierz indeksy próbek, które mieszczą się w czasie pomiędzy  timeMarks(l) a timeMarks(l+1) , zwróć uwagę aby jednostki w wektorze time i wektorze znaczników czasu były takie same
             topoAmp = .....;% uśrednij dane dla warunku 1 po osobach w zakresie czasu któremu odpowiadają idx
            topoplot( topoAmp, chanlocs, 'plotrad' ,0.6 ,'maplimits',scaleTopo);
            if l==1
                text(-1.5,0,labels{1},'FontSize',fontSize);
            end
        end
        
        for l=1:(length(timeMarks)-1)
            axes('position',positionTopo2(l,:))
               idx = ... % wybierz indeksy próbek, które mieszczą się w czasie pomiędzy  timeMarks(l) a timeMarks(l+1) , zwróć uwagę aby jednostki w wektorze time i wektorze znaczników czasu były takie same
          topoAmp = .....;% uśrednij dane dla warunku 1 po osobach w zakresie czasu któremu odpowiadają idx
            topoplot( topoAmp, chanlocs, 'plotrad' ,0.6 ,'maplimits',scaleTopo);
            if l==1
                text(-1.5,0,labels{2},'FontSize',fontSize);
            end
        end
        
        % rysowanie skali kolorow
        axes('position',positionCbar)
        c = colorbar;
        set(gca,'box','off');
        axis off
        caxis(scaleTopo);
        set(c,'FontSize',fontSize)
        text(9,0.5,'[\textit{$\mu$}V]','FontSize',fontSize,'Interpreter','Latex')
    end
    
    orient landscape
    set(findall(gcf,'-property','FontName'),'FontName','Liberation Serif')
    set(gcf,'DefaultFigurePaperType','A4');
end

Wracamy na chwilę do run_me.m: ANOVA

Prosty opis logiki i liczenia ANOVA z powtarzanymi pomiarami: [1]

Przychodzi moment na wykonanie analiz statystycznych. Chcemy zbadać czy w wybranych przez nas oknach czasowych, średnia wartość potencjału w warunkach czytania słowa różni się od odpowiadającej jej średniej w warunkach czytania pseudo-słowa. Przez średnią wartość potencjału rozumiemy, dla ustalonego okienka czasowego, dla danego ROI i ustalonej osoby, wartość uśrednioną po czasie trwania tego okienka. Do wykonania testu wykorzystamy dwu-czynnikową analizę wariancji z powtórzonymi miarami. Zmienną niezależną jest średni potencjał, czynnikami są warunek (dwa poziomy: słowo/pseudosłowo) i ROI (pięć poziomów: LF, CF, RF, LP, RP). Stosujemy wersję ANOVY z powtórzonymi miarami, bo pozwala ona wydzielić część zmienności w amplitudach średniego potencjału wynikającą z różnic pomiędzy osobami (niezależnie od pozostałych czynników) i ta część zmienności jest z punktu widzenia naszego zadania nieistotna.

Dane dla tej wersji ANOVy muszą mieć postać macierzy (w poniższym kodzie dANOVA) takiej, że:

  • w wierszach są wpisane pojedyncze obserwacje
  • w kolumnach są zmienne
amplituda osoba kod slowo/pseudo ROI
a1 sub1 0 LF
a2 sub1 1 LF
a3 sub1 0 CF
... ... ... ...
a_xxx sub36 1 RP

Macierze te, dla wygody poskładamy w jedną macierz 3-wymiarową dataANOVA, tak aby pierwszy indeks wskazywał okno czasowe.

W run_me.m dopisujemy: 

%% ANOVA

ROI.channels = [[18,13];[15,10];[19,17];[9,4];[11,6]];
ROI.labels   = {'LF','CF','RF','LP','RP'};
N_win = length(timeMarks)-1;

dataANOVA = preparujDaneDlaANOVA(data,timeMarks,ROI,startT,stopT,EEG.srate);
% dataANOVA - [5 okien czasowych x obserwacje x |1.srednia amplituda 2. subject 3.slowo(0)/pseudo(1) 4.ROI|] 


for window = 1:N_win
    dANOVA = squeeze(dataANOVA(window,:,:));
    stats  = rm_anova2(dANOVA(:,1),dANOVA(:,2),dANOVA(:,3),dANOVA(:,4),{'type','ROI'});
    stats2 = PostHocStats(stats, dANOVA, ROI, N_win);
    
    disp(['------------- Time window: ' num2str(timeMarks(window)) '-' num2str(timeMarks(window+1)) '-------------'])
    disp(stats)
    disp(stats2)
end

Funkcja: preparujDaneDlaANOVA.m

function dataANOVA = preparujDaneDlaANOVA(data,timeMarks,ROI,startT,stopT,Fs)
    
    time = (startT:1/Fs:stopT);
   % przygotujmy miejsce na macierz, która pozbiera dane potrzebne do analizy ANOVA
   % jest to macierz trzy-wymiarowa: (okna czasow )x( obserwacje )x( typ informacji)
   % pierwszy indeks numeruje okna czsowe : okien czasowych tyle co  timeMarks  -1
   % drugi numeruje obserwacje: jedna obserwacja dotyczy średniej amplitudy jaką odnotowano dla jednej osoby w konkretnym warunku , w konkretnym obszarze zaintersowania (ROI)
% trzeci indeks wskazuje jaką informację w danej komórce macierzy przechowujemy:  |1 -> srednia amplituda 2-> subject 3 ->slowo(0)/pseudo(1) 4->ROI| 

    dataANOVA = zeros(length(timeMarks)-1,size(data,1)*size(ROI.channels,1)*size(data,2),4);  

    for tim = ... % pętla po okienkach czasowych 
        licz = 1;
        for sub = ... % pętla po indeksach osób
            for type =  ... % pętla po warunkach
                for R = ...% pętla po ROI-ach jest ich tyle co size(ROI.channels,1)
                    idx = ... % wybieramy indeksy próbek należących do czasów pomiędzy timeMarks(tim) a timeMarks(tim+1)
                    tmpData = data(... ); % pobieramy wycinek danych konkretnego typ, od konkretnej osoby, z kanałów wskazywanych przez ROI.channels(R,:), i próbek wskazywanych przez idx
                    tmpData2 = mean(  tmpData ,3); % uśredniamy wybrane kanały
                    dataANOVA(tim,licz,1) = mean(squeeze( tmpData2)); % uśredniamy wybrany zakres czasów
                    dataANOVA(tim,licz,2) = sub; % odnotowujemy który to był badany
                    dataANOVA(tim,licz,3) = type-1; % kodujemy typ
                    dataANOVA(tim,licz,4) = R; % kodujemy region ROI
                    licz = licz + 1;
                end
            end
        end
    end

end


Przykładowe wyniki ANOVA mogą wyglądać tak: 


------------- Time window: 70-120-------------
    'Source'               'SS'          'df'     'MS'         'F'         'p'         
    'type'                 [  0.0039]    [  1]    [ 0.0039]    [0.1842]    [    0.6707]
    'ROI'                  [ 45.8683]    [  4]    [11.4671]    [9.5012]    [9.8407e-07]
    'type x ROI'           [  0.1136]    [  4]    [ 0.0284]    [0.6086]    [    0.6572]
    'type x Subj'          [  0.6489]    [ 31]    [ 0.0209]          []              []
    'ROI x Subj'           [149.6574]    [124]    [ 1.2069]          []              []
    'type x ROI x Subj'    [  5.7855]    [124]    [ 0.0467]          []              []

W tym przypadku widzimy tylko główny efekt ROI, ale to mówi tylko tyle, że średnie potencjały różnią się pomiędzy ROIami -> mamy "niepłaską" topografię.

Dla innych okienek czasowych możemy uzyskać istotne:

  • efekty główne dla typu: świadczyłyby one o tym żę w danym okienku czasowym średnie potencjały dla słów są różne niż dla pseudosłów
  • sprzężenie pomiędzy typem i ROIem: świadczyłyby one o tym , że w różnych ROI efekty typu są różne

Ponieważ wykonujemy ANOVe dla każdego z okienek czasowych niezależnie (np. czy mamy efekt główny typu w pięciu okienkach czasowych), to patrząc na efekty tego samego czynnika powinniśmy uwzględnić poprawkę na wielokrotne porównania, np. poprawkę Bonferroniego.

Uzyskanie istotnego efektu w analizie ANOVA mówi tylko o tym, że średnie uzyskane przy grupowaniu danych wg. tego czynnika nie są równe. Aby dowiedzieć się, które średnie różnią się od których wykoujemy tzw. testy post-hoc, np. serię testów t. W naszym przypadku ponieważ czynnik 'type' ma tylko dwa poziomy to nie ma potrzeby podążania z testami post-hoc. Jeśli wykryjemy istotny efekt sprzężenia to warto go dalej przeanalizować. Robi to kolejna funkcja.

Funkcja: PostHocStats.m

function [stat] = PostHocStats(stats,dataANOVA,ROI,N_win)

    TypeROIpval = stats(strcmp(stats(:,1),'type x ROI'),6);
    stat={};
    if TypeROIpval{1} < 0.05/N_win % poprawka Bonefroniego na liczbę okienczasowych
        stat(1,:) = {'ROI','df','t','p'};
        for R = 1:size(ROI.channels,1)
            idxR = (dataANOVA(:,4)==R);
            idxW = (dataANOVA(:,3)==0);
            idxP = (dataANOVA(:,3)==1);
            dataWords  = dataANOVA((idxR & idxW),[1,2]);
            dataPseudo = dataANOVA((idxR & idxP),[1,2]);
            [h,p,ci,stats_t] = ttest(dataWords,dataPseudo);
            if p(1)<=0.05/length(ROI.labels) % poprawka Bonefroniego na liczbę ROI
                stat(R+1,:) = {ROI.labels{R},stats_t.df(1),stats_t.tstat(1), p(1)};
            else
                stat(R+1,:) = {ROI.labels{R},[],[], []};
            end
        end
    end

end

Wracamy na chwilę do run_me.m: Rysowanie ERP

W run_me.m dopisujemy: 

%% rysowanie ERP w ROI

plotERP_ROI(data,startT,stopT,EEG.srate,timeMarks,ROI)

Funkcja: plotERP_ROI.m

Uzupełnij uśrednianie ERPa: 

function [] = plotERP_ROI(data,startT,stopT,Fs,timeMarks,ROI)
    
    position = [[1,2];[3,4];[5,6];[8,9];[10,11]];
    colors = [[0 0 1];[1 0 0]];
    time = (startT:1/Fs:stopT)*1000;
    fontSize = 12;
    font = 'Liberation Serif';
    XLIM = [-200 timeMarks(end)];
    
    % obliczanie skali Y
    YLIM = zeros(2,size(ROI.labels,2));
    for R = 1:size(ROI.labels,2)
        d = squeeze(mean(mean(data(:,:,ROI.channels(R,:),:),2),3));
        YLIM(:,R) = [min(min(d)) max(max(d))];
    end
    YLIM = [min(min(YLIM(1,:))) max(max(YLIM(2,:)))];
    
    figure()
    for R = 1:size(ROI.labels,2)
        subplot(2,6,position(R,:))
        handles = [];
               
        
        % rysowanie przebiegu ERP
        for type = 1:size(data,1)
            hold on
            erp = ... ; % sygnał dla danego typu, danego ROI, uśredniony po osobach i kanałach należących do danego ROI. 
            p = plot(time,erp,'Color',colors(type,:));
            handles = [handles p];
        end
        title(ROI.labels{R})
        xlim(XLIM)
        ylim(YLIM)
        xlabel('Czas [ms]')
        ylabel('Amplituda [uV]')
        legend(handles,{'Słowa','Pseudo'},'Box','off')
        set(gca,'box','off','TickDir','out');
         
        % rysowanie pionowych lini znacznikow czasowych
        for l=1:length(timeMarks)
            hold on
            line([timeMarks(l),timeMarks(l)], YLIM, 'Color', [0.5 0.5 0.5])
            legend off
        end
        hold on
        line([0 0],YLIM,'Color',[1 0 0])
        drawaxis(gca,'x',0,'y',0);
        axis off
    end

    set(findall(gcf,'-property','FontName'),'FontName',font)
    set(findall(gcf,'-property','FontSize'),'FontSize',fontSize)

end

Laboratorium_EEG/EEGLAB

Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=Laboratorium_EEG/EEGLAB&oldid=9163”