====== Cvičení 11: FMRI - t-testy ====== ===== Úvod ===== Úkolem cvičení je nalézt oblasti mozku aktivované při poslechu dvouslabičných slov pomocí [[http://fmri.mchmi.com/main_index.php?strana=17#ttesty|t-testu]] z dat z funkční magnetické rezonance (fMRI). ==== Data ==== Pro experimenty budete používat [[http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/fmri-t-test/spm_data.mat|data]] distribuovaná k balíku [[http://www.fil.ion.ucl.ac.uk/spm/software/spm8/|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. {{ :courses:a6m33zsl:stimulace.jpg?600 |Schéma stimulace.}} ===== Statistická analýza ===== 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 [[http://www.wikiskripta.eu/index.php/Student%C5%AFv_t-test|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 [[https://en.wikipedia.org/wiki/Student's_t-distribution|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. {{ :courses:a6m33zsl:fmri.jpg?600 |Schéma vyhodnocování hypotézy.}} ==== Postup ==== - 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]]). - 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//). - 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//). - 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. ===== Vizualizace ===== Pro kontrolu výsledků můžete použít jednoduchý prohlížeč 3D dat //[[http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/fmri-t-test/projector.zip|projector]]//, který napsal Jakub Krátký (syntaxe: $projector(objem,H_mu,.5);$). Své výsledky můžete porovnat s Projektorem. {{ :courses:a6m33zsl:projector.jpg?600 |Vizualizace tří navzájem kolmých řezů mozku se zvýrazněnými aktivovanými částmi pomocí nástroje Projector.}} 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é). {{ :courses:a6m33zsl:fmri_t_tests_complex_rezy.png?600 |Všechny tři navzájem kolmé řezy mozkem v bodě [10,48,29] s vyznačenými aktivními oblastmi. Hodnoty T_\mu jsou kódovány barevně.}} ==== Struktura obrázků v Matlabu ==== 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. ==== Postup ==== Popis funkcí které budete potřebovat pro vizualizaci vašich výsledků spolu s příklady jejich použití najdete [[http://www.mathworks.com/company/newsletters/news_notes/2009/tips.html|zde]]. - 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//). - 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. - 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//). - 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. ==== Poznámky na závěr ==== 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í. ==== Bonus ==== - 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í. - Použijte Bonferroniho [[http://en.wikipedia.org/wiki/Bonferroni\_correction|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? {{tag>fMRI, t-tests}}