==== Lab 12 - fMRI ====

V dnešním cvičení se podíváme na data z funkční magnetické rezonance (fMRI) a bude nás zajímat detekce regionů v mozku, které jsou aktivovány při poslouchání dvouslabičných slov. Aktivita mozku je měřena pomocí fMRI, konkrétněji měříme BOLD signál při resting-state fMRI sekvenci. Statisticky významné rozdíly mezi běžnou aktivitou (//resting state//) a aktivitou v průběhu stimulace (//active state//) budeme určovat na základě $t$-testu.

==== Domácí cvičení ====

Domácí cvičení má dvě části:

**[[#Statistická analýza]]**
  - Rozdělte vstupní data na části odpovídající klidovým a aktivním blokům (resting/active stav)
  - Spočítejte výběrový průměr $\bar{x}_a, \bar{x}_r$, v každém voxelu $\bar{x}$
  - Spočítejte výběrové rozptyly $S^2_a$ a $S^2_r$
  - Najdětě voxely $\bar{x}$ se signifikatními rozdíly mezi resting a active stavy (na hladině významnosti $\alpha=0.01$).
  
  
**[[#Vizualizace]]**
  - vykreslete signál aktivity v aktivním (plná čára) a v klidovém stavu (čárkované) pro pozice: [47,30,28] and [47,32,20]
  - vykreslete statistiku //T// jako semi-transparentní barevnou masku přes první snímek fMRI sekvence. Použije jet-colormap
     * vykreslete všechny tři standardní roviny procházející voxelem [47,30,28]
     * pro lepší orientaci vykreslete v každém řezu barevné čáry, které znázorňují postavení ostatních rovin (jak to znáte např. z MITK)
  

==== Data ====

Použíjeme experimentální data z balíčku SPM. Design snímaného experimentu je následující: 
  * subjekt je bez cílených vnějších vlivů (//resting state//)
  * subjekt poslouchá dvouslabičná slova (//active state//)
  
Sekvence obsahuje dohromady 84 volumetrických snímků (dimenze $53 \times 63 \times 52$ s izotropními voxely o rozměru $3~mm$). Data pochází z přístroje 2$T$ Siemens MAGENTOM Vision. Jako první je 6 skenů v //resting state//, pak následuje 6 snímků v //active state// a takto alternují až do konce datové řady (snímek 84). 

Snímky jsou k dispozici{{ :courses:zsl:fmri_auditory.mat |zde}}. Data jsou uložena (pro zpracování v Matlabi) jako 4D-matice $Y$, kde první tři komponenty odpovídají prostorovým dimenzím a časová osa je v poslední komponentě, tedy $Y(:,:,:,1:6)$ je výběr 3D snímků prvního //resting state// bloku, $Y(:,:,:,7:12)$ pak první //active state// sekvence.

==== Statistická analýza ====

Hypotáze je, že se //střední hodnoty nahrané v průběhu stimulace liší od stř. hodnot nahraných v klidové fázi//, testovat budeme na dané hladině významnosti $\alpha$. Pro jednoduchost budeme předpokládat, že jsou aktivity v jednotlivých voxelech navzájem nezávislé a budeme testovat zvlášť pro každý voxel. Místa, ve kterých detekujeme signifikatní rozdíly by měla být do určité míry provázána s daným úkolem.

=== Testování ===

Použijeme Welchúv test abychom ověřili hypotézu o rovnosti středních hodnot $$\mathcal{H}_{0,\mu}: \mu_a = \mu_r$$ v průběhu activation ($\mu_a$) a v klidovém stavu ($\mu_r$). Použijeme testovací statistiku:

\begin{equation}
T_{\mu}=\frac{\overline{x}_s-\overline{x}_k}{\sqrt{\frac{S^2_a}{m}+\frac{S^2_r}{n}}}\, ,
\end{equation}
s odhadem směrodatné odchylky $S^2_a$ a $S^2_r$, definovaným následovně:
\begin{equation}
	S^2=\frac{1}{n-1} ~ \sum_n \left( x_i - \overline{x} \right)^2
\end{equation}

$\overline{x}_a$, $\overline{x}_r$ jsou výběrové průmery, $m$ a $n$ počty vzorků ve výběru.

Statistika $T_{\mu}$ má Studentovo $t-$[[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. Použijeme tedy kvantily Studentova rozdělení abychom otestovali hodnotu $T_\mu$. Zamítneme //nulovou hypotézu// pokud $|T_\mu|>f_{\psi}(\alpha/2)$. $f_{\psi}()$ je inverzní funkce ke Studentovu rozdělení (MATLAB: //tinv//), $\alpha$ hladina významnosti a $\psi$ počet stupnů volnosti. 

Výsledkem téte analýzy je 3D matice $\mathbf{H}$ o stejné prostorové dimenzi jako vstupní data, s hodnotou 1 ve voxelech, kde pro danou hladinu alpha zamítáme nulovou hypotézu a s hodnotou 0 jinak. 

/*
** Bonus (Bonferroniho korekce)**

Cílem je dosáhnout hladiny významnosti $\alpha$ pro celý obrázek, ne pro jednotlivé voxely, jako to máme ve zjednodušeném modelu výše. Nejsnazší, ale zároveň nejvíce restriktivní transformací je //Bonfeerroniho korekce//, která namísto hladiny $\alpha$ testuje na hladíně $\frac{\alpha}{N}$. Abychom touto korekčí nepřišli o veškerou detekovanou aktivitu, je dobré dosáhnout co nejmenších hodnot $N$ -- před výpočtem korigované statistiky si spočteme binární masku pro mozkovou tkáň (např. threshold aplikovaný na první volume) a následně budeme vyhodnocovat pouze pro voxely unitř masky.
*/

==== Vizualizace ====

Abychom dostali náhled na jednotlivé roviny, použijeme běžný //array-slicing//. Axiální rovinu ($z = 20$) vstupního datasetu $Y$ pro čas $t=1$ získme přímo voláním $Y(:,:,20,1)$. Aplikací funkce //squeeze// po výběru roviny odstaníme //prázdné// dimenze, to je potřeba pro vykreslování pomocí imshow nebo imagesc. Aby souhlasila barevná škála ve všech třech řezech, použijte identické rozpětí intenzit ''[v_min, v_max]'' v každé volání funkce //imagesc//. 

Překrývající vrstvy lze vykreslovat v režimu ''hold on''. Jako první vykreslíme obrázek na pozadí, pak vykreslíme hodnotu statistiky (zde si musíme zapamatovat handle //h//) použitím barevné mapy ''jet''. Handle //h// potřebujeme, abychom pomocí masky signifikantních voxelů nastavili masku průhlednosti<code matlab>
...
hold on;
h = imshow( T_image ... );
set(h,'AlphaData', H_mask, ...);
</code>

Ukázka jedné z rovin je na dalším obrázku:

{{:courses:zsl:fmri_vis.png?600|}}