Cvičení 11: FMRI - t-testy

Úkolem cvičení je nalézt oblasti mozku aktivované při poslechu dvouslabičných slov pomocí t-testu z dat z funkční magnetické rezonance (fMRI).

Pro experimenty budete používat data distribuovaná k balíku SPM. (Data jsou z SPM2, dnes je k dispozici SPM8.)

Subjekt experimentu byl vystaven sluchové stimulaci pomocí dvouslabičných slov s frekvencí $60 \text{min}^{-1}$ střídané s klidem. Celkem bylo nahráno 84~objemových obrázků ($79\times 95\times 68$ s velikostí voxelu $3\times 3\times 3 \text{mm}$) v blocích po šesti (6~scanů klid, 6 scanů stimulace…). Pro měření byl použit modifikovaný systém 2T Siemens MAGNETOM Vision.

Pro účely Matlabu jsou data uložena ve 4D matici, kde první tři dimenze odpovídají prostorovým souřadnicím a čtvrtá udává pořadí v čase. (Y(:,:,:,1) je první 3D obrázek, Y(:,:,:,1:6) je první klidová fáze experimentu a Y(:,:,:,7:12) je první fáze experimentu se stimulací.)

Data jsou předregistrována aby si jednotlivé voxely časové posloupnosti odpovídaly a rozmazány 6mm Gausovským jádrem, první cyklus sekvence (klid + stimulace = 12 scanů) je odstraněn kvůli artefaktům.

Vaším úkolem bude zjistit, zda se střední hodnoty signálu změřeného při stimulaci a bez stimulace liší na dané hladině významnosti. Budeme předpokládat, že jednotlivé voxely obrázku jsou nezávislé a statistickou analýzu provedeme na každém z nich. Tam kde se budou střední hodnoty měřených signálů lišit, prohlásíme, že mozek byl aktivovaný slovní stimulací.

Provedeme Welchův test hypotézy $\mathcal{H}_{0\mu}\colon \mu_s=\mu_k$ o rovnosti středních hodnot sekvencí při a bez stimulace. Budeme testovat statistiku:

\begin{equation} T_{\mu}=\frac{\overline{x}_s-\overline{x}_k}{\sqrt{\frac{S^2_s}{m}+\frac{S^2_k}{n}}}\, , \end{equation} kde $S^2_s$ a $S^2_k$ jsou nevychýlené výběrové rozptyly měření při stimulaci, respektive v klidu definované \begin{equation} S^2=\frac{1}{n-1} ~ \sum_n \left( x_i - \overline{x} \right)^2 \, , \end{equation} $\overline{x}_s$, $\overline{x}_k$ jsou výběrové průměry a m a n je počet vzorků ve výběrech.

Tato statistika má Studentovo rozdělení s \begin{equation} \psi=\frac{\Big(\frac{S_s^2}{m}+\frac{S^2_k}{n}\Big)^2}{\frac{S^4_s}{m^3-m^2}+\frac{S^4_k}{n^3-n^2}} \end{equation} stupni volnosti, proto testujeme hodnotu statistiky $T_\mu$ s kvantily Studentova rozdělení. Nulovou hypotézu $\mathcal{H}_{0\mu}$ zamítneme, když $|T_\mu|>f_{\psi}(\alpha/2)$, kde $f_{\psi}()$ je inverzní ditribuční funkce Studentova rozdělení, $\alpha$ je hladina významnosti a $\psi$ počet stupňů volnosti.

1. Rozdělte vstupní data na dvě části, jednu pro stimulační (m vzorků na voxel) druhou pro klidovou (n vzorků na voxel) fázi experimentu (viz popis dat v Data).
2. Pro každý voxel spočtěte průměrnou hodnotu v čase pro stimulaci $\overline{x}_s(x,y,z)$ a pro klid $\overline{x}_k(x,y,z)$ (funkce mean).
3. Spočtěte nevychýlený vyběrový rozptyl dat v čase pro obě fáze měření $S^2_s(x,y,z)$ a $S^2_k(x,y,z)$ (funkce var).
4. Zjistěte, zda mají data pro klidovou a stimulační fázi experimentu stejnou střední hodnotu (testujte dle sekce Statistická analýza na hladině významnosti $\alpha=1\%$). Vytvořte matici $H_\mu$ obsahující samé logické 0 a pro voxely, kde zavrhujeme $\mathcal{H}_{0\mu}$ o rovnosti středních hodnot (použijte funkci $tinv(1-alpha/2,psi)$), dejte logickou 1.

Pro kontrolu výsledků můžete použít jednoduchý prohlížeč 3D dat projector, který napsal Jakub Krátký (syntaxe: $projector(objem,H_mu,.5);$). Své výsledky můžete porovnat s Projektorem.

Do reportu vložte časové průběhy fMRI dat pro voxel [6,48,29] (zelenou lomenou čáru pro klidovou, červenou pro stimulační část experimentu vykreslete do jednoho grafu) a tři řezy (frontální, transverzální a sagitální) bodem [41,45,29]. Do každého obrázku vykreslete řezy z prvního 3D obrázku (Y(:,:,:,1), anatomický obraz) a současně hodnoty $T_\mu$ pro oblasti kde zavrhneme $\mathcal{H}_{0\mu}$ o rovnosti středních hodnot. Návod jak takové obrázky vytvořit je v sekcích Struktura obrázků v Matlabu a Postup. Příklad, jak by takové obrázky mohly vypadat je na následujícím obrázku (vaše výsledky mohou být jiné).

Obrázky a grafy jsou v Matlabu zapouzdřeny do objektu třídy figure. Potomkem objektu figure je objekt axis, který může mít další potomky např. objekty třídy image nebo line. Přistupovat k instancím jednotlivých tříd se dá pomocí handles. Handles lze získat buď při vytváření objektů, nebo funkcí get a její speciální formy gcf a gca, které vrací handle na aktivní figure respektive axis. Pomocí funkce get můžeme získat i libovolné parametry objektu, například pole Children obsahující handles na potomky.

Funkce set pak umožňuje nastavit parametry objektů. Strukturu všech parametrů můžete zjistit pomocí helpu, nebo property inspectoru.

Popis funkcí které budete potřebovat pro vizualizaci vašich výsledků spolu s příklady jejich použití najdete zde.

1. Pomocí funkcí imagesc, colormap a squeeze vykreslete daný řez daty. Dbejte toho, aby všechny tři řezy měly stejné mapování hodnot na jasy (parametr imagesc).
2. Pomocí funkce jet vytvořte LUT (lookup table) pro mapování $T_\mu$ na pseudo-barvy a vytvořte 4-rozměrnou matici $T_{\mu 4D}$ ve které první tři souřadnice reprezentují polohu a čtvrtý rozměr bude reprezentovat RGB kanály pro visulaizaci dat.
3. Pomocí funkce h=imshow(…) vykreslete řezy maticí $T_{\mu4D}$ do odpovídajících obrázků dat. V proměnné h bude handle na obrázek (nezapomeňte na hold on a squeeze).
4. Nastavte set(h,…) parametr AlphaData objektu třídy image matice $H_{\mu}$.

Toto není jediná možnost jak vrstvit a maskovat obrázky v Matlabu, ale je to podle našich zkušeností možnost nejméně bolestivá. %Pokud byste chtěli proniknout hlouběji, můžete si pomocí funkcí gca a get nalézt handles na instance libovolného objektu.

Vyhodnocování aktivity pomocí t-testu je velmi jednoduchá, ale ne příliš robustní metoda. Jak jste si však mohli ověřit, dává na jednoduchých datech obstojné výsledky. Na druhou stranu jste si také nemohli nevšimnout velkého množství malých ostrůvků aktivity způsobené artefakty a šumem.

Pro vylepšení výsledků lze získaný binární obraz upravovat, například tak že odstraníme “aktivované” oblasti menší než daný práh. Pokud si chcete vyzkoušet efekt “ladění” výsledků, zkuste to pomocí funkce bwareaopen. Další možností jak se zbavit malých aktivovaných oblastí je filtrovat vstupní data trojrozměrnou (prostorovou) dolní propustí.

1. Pro odstranění falešných aktivací v oblastech mimo mozek si vytvořte jednoduchou 3D segmentaci mozku (např. jednoduchým pruhováním intenzit) a tu následně použijeme k maskování detekovaných aktivací.
2. Použijte Bonferroniho korekci, t.j. vyžadujte, aby pravděpodobnost falešně pozitivního výsledku byla $\alpha$ na celý obrázek. To provedete nahrazením $\alpha$ hodnotou $\frac{\alpha}{N}$, kde N je počet všech voxelů. Co pozorujete?