En jämförelse mellan frekventistisk och Bayesiansk Dual Regression : för nätverkskartor i hjärnan vid resting-state fMRI

(1)

Kandidatuppsats i Statistik

En jämförelse mellan frekventistisk och

Bayesiansk Dual Regression

för nätverkskartor i hjärnan vid resting-state fMRI

Patrick Jonsson Jacob Welander

Avdelningen för Statistik och maskininlärning Institutionen för datavetenskap

Linköpings universitet

(2)

Handledare: Bertil Wegmann, universitetslektor Examinator: Linda Wänström, universitetslektor

(3)

Sammanfattning

Att undersöka områden i hjärnan som är aktiva utan att någon stimulans sker kan ge information om en in-divids standardnätverks basnivå. Denna basnivå kan användas för att identifiera avvikande spatiala mönster i hjärnan som associeras med sjukdomar och funktionsnedsättningar. Denna uppsats syftar till att undersöka hur skillnaderna ser ut för individspecifika nätverkskartor genom att jämföra tre olika anpassningar av Dual Regression, en frekventistisk och två Bayesianska modeller. Datamaterialet som analyseras i uppsatsen är från Cambridge-Buckner, en del av 1000 Functional Connectomes Project som innehåller fMRI-data. Från datamaterialet har även tillhörande förhandsskattade gruppvisa oberoende komponenter erhållits från 20 ut-valda individer vilket sedan används i uppsatsen för att skatta individspecifika nätverkskartor i hjärnan för tre individer från studien.

Det anpassas tre olika Dual Regressions-modeller: En frekventistisk modell med homoskedastisk varians, en Bayesiansk modell med heteroskedastisk varians med okorrelade feltermer samt en Bayesiansk modell med he-teroskedastisk varians och korrelerade feltermer. För de två Bayesianska modellerna används icke-informativa priorfördelningar. Dessa olika modeller skiljer sig åt då de kan ta hänsyn till olika mängder av information genom att ha olika komplexa kovariansstrukturer. Det observeras att den frekventistiska modellen och den Bayesianska modellen med heteroskedastisk varians och okorrelerade feltermer skattar nätverk som är i stor utsträckning lika varandra. Den Bayesianska modellen med heteroskedastisk varians och korrelerade felter-mer tenderar att skatta nätverk som är skild från de andra modellerna, där det ofta förekom skillnader i nätverkens former samt en del amplitudskillnader. I kovariansmatrisen för den Bayesianska modellen med heteroskedatisk varians och korrelerade feltermer observeras ett flertal höga korrelationer mellan feltermerna vilket indikerar på att det bör tas hänsyn till korrelerade feltermer.

Det diskuteras även om problem som förekommer hos respektive tillvägagångssätt för att skatta modellen, där frekventistiska tillvägagångssättet inte tar hänsyn till all information i data men är enkel att anpassa. Den Bayesianska modellen med heteroskedastisk varians och okorrelerade feltermer ger liknande resultat som det frekventistiska tillvägagångssättet. Den Bayesianska modellen med heteroskedastisk varians och korrelerade feltermer ger resultat som anpassar data bättre än de andra två modellerna men är mer komplex att beräkna.

(4)

Abstract

Examining regions in the brain that are active without any stimuli gives information about an individual’s default brain networks. These default mode networks can be analyzed to identify deviating spatial patterns in the brain that are associated with diseases and disabilities. This thesis aims to analyze the difference in how frequentist and Bayesian Dual Regression estimates subject specific spatial-maps. We received pre-estimated groupwise independent components from 20 individuals based off of fMRI-data from the Cambridge-Buckner dataset which is part of the 1000 Functional Connectomes Project. These are later used to create subject specific spatial-maps for 3 individuals in the study.

In this thesis 3 different types of Dual Regression models will be fitted: A frequentist Dual Regression, A Bayesian model with heteroscedastic variance and uncorrelated error terms and a Bayesian model with heteroscedastic variance and correlated error terms. Non-informative prior distributions are used for both Bayesian models. As these 3 models can account for varying amounts of information in the data due to varying complexity of the covariance structure some difference were observed in the subject specific maps. The frequentist Dual Regression and the Bayesian model with heteroscedastic variance and uncorrelated error terms often showed similar results, however the resulting networks from the Bayesian model with het-eroscedastic variance and correlated error terms often differed from the other two models. The difference was observed both in network shapes and in activation amplitude. The covariance matrix for the Bayesian model with heteroscedastic variance and correlated error terms contained a number of high correlations between the error terms, indicating that correlation among error terms should be taken into account.

Some arguments are made for respective way of fitting the model as each model has its unique advantage and disadvantage; where the frequentist model does not take into account all information from the data it is easy to fit. The Bayesian model with heteroscedastic variance and uncorrelated error terms is also relatively easy to fit and provides similar results to the frequentist model. The Bayesian model with heteroscedastic variance and correlated error terms however does account for more information and yields better results but is more computationally expensive.

(5)

Förord

Denna kandidatuppsats är skriven vid Linköpings Universitet på Institutionen för datavetenskap. Vi vill tacka Anders Eklund på Institutionen för medicinsk teknik och Institutionen för datavetenskap på Linköpings Universitet som tillhandahållit oss med datamaterialet för uppsatsen. Vi vill även särskilt tacka vår handledare och uppdragsgivare Bertil Wegmann för sitt stöd och inspiration under uppsatsens gång. Slutligen vill vi tacka våra opponenter Jakob Lindén och Sidney Rydström som efter granskning har bidragit med värdefulla åsikter på uppsatsen.

(6)

Innehåll

Centrala begrepp viii

1 Introduktion 1

1.1 Bakgrund . . . 1

1.2 Syfte . . . 2

1.2.1 Avgränsningar . . . 2

1.3 Etiska och samhälleliga aspekter . . . 2

2 Data 4 2.1 Centrering av data . . . 6

3 Metod 7 3.1 Matematiska operatörer . . . 7

3.2 Multivariat linjär regression . . . 8

3.3 Dual Regression . . . 9

3.4 Frekventistisk Dual Regression . . . 9

3.4.1 Regressionssteg för individspecifika tidsförlopp . . . 9

3.4.2 Regressionssteg för individspecifika nätverkskartor . . . 10

3.5 Sannolikhetsfördelningar . . . 10

3.5.1 Invers Wishartfördelning . . . 10

3.5.2 Matrisnormalfördelning . . . 11

3.6 Bayes sats . . . 11

3.7 Bayesiansk multivariat linjär regression . . . 12

3.8 Bayesiansk Dual Regression . . . 13

3.8.1 Regressionssteg för individspecifika tidsförlopp . . . 13

3.8.2 Regressionssteg för individspecifika nätverkskartor . . . 14

3.9 Bayesiansk Dual Regression utan korrelerade feltermer . . . 14

3.10 Skillnader mellan frekventistisk och Bayesiansk Dual Regression . . . 15

3.11 Implementering i R . . . 15

4 Analys 16 4.1 Val av priorfördelningarnas hyperparameterar . . . 16

4.2 Utvärdering av posteriordragningar . . . 16 4.2.1 Individ 1 . . . 17 4.2.2 Individ 2 . . . 19 4.2.3 Individ 3 . . . 20 4.3 Nätverkskartor . . . 22 4.3.1 Individ 1 . . . 22

(7)

4.3.3 Individ 3 . . . 28 4.4 Korrelerade feltermer . . . 31 4.4.1 Individ 1 . . . 31 4.4.2 Individ 2 . . . 32 4.4.3 Individ 3 . . . 32 5 Diskussion 33 5.1 Datamaterial . . . 33 5.2 Metod . . . 33 5.3 Resultat . . . 34 5.4 Framtida forskning . . . 34 6 Slutsatser 36 7 Bilaga I

7.1 Bilaga 1, Voxelvis utveckling av ackumulerat posteriormedelvärde . . . I 7.2 Bilaga 2, Nätverkskartor . . . VII

(8)

Figurer

2.1 Exempel på fMRI-data, framtaget i R . . . 4

2.2 Uppdelning av hjärnregioner . . . 5

2.3 GICA 1-20 . . . 6

4.1 Individ 1:Utveckling av ackumulerade posteriormedelvärde . . . 17

4.2 Individ 1: Utveckligen av voxelvisa ackumulerade posteriormedelvärdet . . . 17

4.3 Individ 1: Histogram för P(βββSM_{> 0)} _{. . . .} ₁₈

4.4 Individ 2: Utveckling av ackumulerade posteriormedelvärde . . . 19

4.6 Individ 2: Histogram för P(βββSM_{> 0)} _{. . . .} ₂₀

4.7 Individ 3: Utveckling av ackumulerade posteriormedelvärde . . . 20

4.9 Individ 3: Histogram för P(βββSM_{> 0)} _{. . . .} ₂₁

4.10 Individ 1: Nätverkskarta 1 . . . 22

4.13 Individ 1: Differensen i aktivering mellan Bayesianska modellernas nätverkskartor 1, 4 och 6 . 24 4.14 Individ 1: Standardavvikelse för nätverkskarta 1, 4 och 6 . . . 24

4.15 Individ 1: Nätverkskarta 6 för modell 3, tröskelnivå 50% och 67% . . . 25

4.19 Individ 2: Differensen i aktivering mellan Bayesianska modellernas nätverkskartor 1, 4 och 6 . 27 4.20 Individ 2: Standardavvikelse för nätverkskarta 1, 4 och 6 . . . 27

4.25 Individ 3: Differensen i aktivering mellan Bayesianska modellernas nätverkskartor 1, 4 och 6 . 30 4.26 Individ 3:Standardavvikelse för nätverkskarta 1, 4 och 6. . . 30

4.28 Individ 1: Histogram över korrelerade feltermer . . . 31

4.29 Individ 2: Histogram över korrelerade feltermer . . . 32

4.30 Individ 3: Histogram över korrelerade feltermer . . . 32 7.1 Individ 1: Utveckling i z-led för nätverkskarta 4 för modell 3 . . . I 7.2 Individ 1: Utveckling i z-led för nätverkskarta 6 för modell 3 . . . II

(9)

7.4 Individ 2: Utveckling i z-led för nätverkskarta 6 för modell 3 . . . IV 7.5 Individ 3: Utveckling i z-led för nätverkskarta 4 för modell 3 . . . V 7.6 Individ 3: Utveckling i z-led för nätverkskarta 6 för modell 3 . . . VI 7.7 Individ 1: Nätverkskarta 2 för de tre olika modellerna . . . VII 7.8 Individ 1: Nätverkskarta 3 för de tre olika modellerna . . . VII 7.9 Individ 1: Nätverkskarta 5 för de tre olika modellerna . . . VIII 7.10 Individ 2: Nätverkskarta 2 för de tre olika modellerna . . . VIII 7.11 Individ 2: Nätverkskarta 3 för de tre olika modellerna . . . IX 7.12 Individ 2: Nätverkskarta 5 för de tre olika modellerna . . . IX 7.13 Individ 3: Nätverkskarta 2 för de tre olika modellerna . . . X 7.14 Individ 3: Nätverkskarta 3 för de tre olika modellerna . . . X 7.15 Individ 3: Nätverkskarta 5 för de tre olika modellerna . . . XI

(10)

Centrala begrepp

fMRI: Funktionell magnetresonanstomografi

BOLD: Blood-oxygen-level-dependent imaging, används i fMRI för att undersöka aktiva organ RSN: Resting state networks, standardnätverk i hjärnan när den är i viloläge

RSFC: Resting state functional connectivity, områden där aktivering i nervceller inträffar fastän ingen sti-mulans sker

(11)

1. Introduktion

I detta kapitel introduceras bakgrund, syfte, etiska aspekter samt övergripande information om datamateri-alet.

1.1 Bakgrund

Functional magnetic resonance imaging (fMRI) är en icke-invasiv metod som används för att mäta magnetiska egenskaper i hjärnan genom att studera mängden syresatt och icke-syresatt blod. Vid aktivering av nervcel-lerna ökar syrekonsumtionen vilket leder till att syresatt blod transporteras till det aktiva området, således kan en signal mätas och kartlägga områden i hjärnan som är aktiva. För att kartlägga de aktiva områdena i hjärnan används Blood-oxygen-level dependent (BOLD) fMRI, där en signal ger magnituden för den indirekta aktiviteten i nervcellerna. Från BOLD-signalen kan resting-state network (RSN) tas fram för att undersöka nätverk i hjärnan som har en förändring i BOLD-signalen fastän ingen uppgift eller stimulans sker (Huettel, Song & McCarthy, 2004). Under fMRI-undersökningen kan individens aktivitet, exempelvis fysiska rörelser påverka resultaten av mätningen genom att det generar ett blodflöde till hjärnan som inte liknar blodflödet vid RSN och där med skapa brus. Genom att säkerställa att förutsättningarna vid fMRI-undersökningen är lika för alla individer som undersöks kan brus och avvikelser hanteras. Det kan även krävas databearbeting för att minska bruset (Nickerson, Smith, Öngür & Beckmann, 2017).

Independent Component Analysis (ICA) används för att hitta komponenter från BOLD-signalerna som är oberoende. Detta görs genom att anta att signalerna är centrerade runt 0 och oberoende. ICA kan använ-das för att separera oberoende spatiala aktivitetsmönster i hjärnan. Det kan även använanvän-das för att separera temporala aktiveringsmönster baserat på fMRI-mätningarna som gjorts över tid i respektive del i hjärnan. ICA har visat sig kunna separera flera olika källor av signaler i hjärnan som är oberoende. En Group-average Independent Component Analysis (GICA) kan utföras baserat på flera individer genom att skapa en samling av nätverk som är vanligt förekommande hos individerna. De spatio-temporala signaler som ofta är kopplade till RSN nätverken kan identifieras genom att spatialt kartlägga hjärnan med hjälp av GICA (Kiviniemi, Kantola, Jauhiainen, Hyvärinen & Tervonen, 2003).

Genom att undersöka hjärnans resting-state functional connectivity (RSFC) områden, vilket är områden där aktivering i nervceller inträffar fastän ingen stimulans sker, kan en basnivå för en grupp individer identifieras (Morcom & Fletcher, 2007). Denna basnivå kan i sin tur användas för att identifiera avvikande spatiala möns-ter i hjärnan för patienmöns-ter med till exempel autism, alzheimer’s, bipolär sjukdom och schizofreni (Woodward och Cascio, 2015; Calhoun, Maciejewski, Pearlson och Kiehl, 2008; Nickerson m. fl., 2017; Salman m. fl., 2019; Mejia, Nebel, Wang, Caffo och Guo, 2019). En vanlig förekommande metod för att skatta individspecifika nätverkskartor är Dual Regression. Detta är en spatio-temporal regressionsmetod som utförs i två steg för att undersöka RSFC-nätverkskartor i hjärnan. I regressionens första steg skattas de nätverksspecifika tids-serierna för varje del i hjärnan. Parametrarna från första steget används senare som förklarande variabler i regressionens andra steg som skattar de individspecifika nätverkskartorna. I tidigare forskning har Dual Regression tillsammans med GICA kunnat identifiera gruppskillnader i hjärnans funktionella konnektivitet (Nickerson m. fl., 2017). I studien kunde Dual Regression även identifiera de individspecifika tidsförloppen och nätverkskartorna där en korrekt inferens genomfördes när regressionen anpassades genom multivariat minsta kvadratmetoden (OLS).

Då det frekventistiska tillvägagångsättet anpassar en modell genom multivariat OLS vilket endast ger en punktestimering uppstår möjligheter att förbättra dessa resultat genom att istället använda sig av

(12)

Baye-sianska metoder. Detta eftersom de kan ta hänsyn till all osäkerhet för parametrarna. Vid användning av frekventistisk multivariat OLS antas det även att feltermer är temporalt och spatialt oberoende. Detta re-sulterar i en grov förenkling, eftersom aktiveringen i en av hjärnans områden antas inte vara beroende av aktivering i andra närliggande områden samt att mätningar över tid är oberoende. Detta är ett antagande som är känt att inte stämma, men som vanlig multivariat OLS skattning inte tar hänsyn till (Mejia m. fl., 2019).

Ett problem som kan uppstå vid användning av Bayesianska metoder är att modellerna är svåra att anpassa som följd av att all osäkerhet i parametrarna tas med. Detta problem uppstår på grund av att likelihoodfunk-tionen blir så komplex att den numeriska integralikelihoodfunk-tionen blir långsam, vilket leder till att modellanpassningen kan ta väldigt lång tid att utföra (Woolrich m. fl., 2009).

1.2 Syfte

Syftet med studien är att jämföra hur olika modeller skattar nätverkskartor samt att undersöka om mer komplexa modeller anpassar data bättre. Detta görs för att undersöka om det finns möjligheter att bättra inferens som görs på fMRI-data. Genom multivariat OLS anpassas en frekventistisk Dual Regression samt två olika Bayesianska Dual Regressioner. Resultaten av de tre tillvägagångsätten kommer då jämföras och visualiseras för att undersöka hur skillnaden i resultaten ser ut mellan den frekventistiska punktskattningen och de två Bayesianska modellerna som tar till hänsyn till all osäkerhet i parametrarna när de två stegen av regressionerna anpassas.

Rapportens huvudsakliga frågeställningar som kommer att besvaras är:

• Hur skiljer sig individspecifika spatiala nätverkskartor i hjärnan från anpassningar av frekventisk och Bayesiansk Dual Regression?

• Hur ser skillnaden ut för de individspecifika spatiala nätverkskartorna i Bayesiansk Dual Regression med heteroskedastisk varians när hänsyn är tagen till korrelerade feltermer?

1.2.1 Avgränsningar

Studien begränsas till att endast undersöka hjärnan med en upplösning på 6 mm3_{. Det sker även en}

begräns-ning att data på ett urval av endast 20 individer används för undersökbegräns-ningen. Av de 177 komponenterna kommer endast 20 av dessa används i analysen. Gruppvisa jämförelser kommer inte att genomföras utan endast individspecifika fall. Tre utav 20 individer från studien kommer att analyseras, där de tre olika me-toderna kommer jämföras för varje utvald individ. Dessa avgränsningar genomförs i uppsatsen på grund av arbetets omfattning och dess tidskrävande beräkningar.

1.3 Etiska och samhälleliga aspekter

Data som används i uppsatsen innehåller ingen information som kan vara röjande eller anses som känslig på individnivå. Det förekommer inga direkt kända fysiska hälsorisker när kroppen utsätts för magnetfält som används för kartläggningen av hjärnan vid en fMRI-undersökning (Shechner m. fl., 2013). Patienterna som har undersökts har inte heller utsatts för någon typ av stimuli som till exempel narkos under mätningen som kan potentiellt sett vara skadlig. Det anses därför inte finnas någon etisk känslig aspekt vid insamlingen av data.

(13)

Genom att förbättra inferensen av fMRI-data kan detta leda till bättre insikt i hur både friska och sjuka hjärnor fungerar. fMRI är en vanligt förekommande metod för mätning och fMRI-data används ofta som beslutsgrund för vilken typ av medicinering eller åtgärd som passar bäst för en mängd olika sjukdomar och hälsoproblem. Exempel på dessa kan vara upptäckten och effekten av tumörer, hjärnblödning, hjärnskador samt funktionsnedsättningar som Alzheimers och Autism m.fl. (Bobholz Julie A., 2007). Det är därför av intresse att utvärdera och potentiellt förbättra inferensen.

(14)

2. Data

fMRI-data skapas genom att en magnetkamera använder en serie av olika gradienter av elektromagnetiska fält där en så kallad magnetresonanstomografi skapas som avbildar hjärnan i form av bilder. Individen av intresse undersöks genom att magnetkameran tar bilder skiktvis under tiden bordet som individen ligger på rör sig fram och tillbaka för att skanna hela hjärnan i etapper (Huettel m. fl., 2004). Från en fMRI-mätning kan sedan BOLD-signalen erhållas genom att studera syresatt och icke-syresatt blod för att mäta den indirekta aktiviteten i hjärncellerna. Genom att transformera fMRI-mätningen till analyserbar fMRI-data sker en chif-frering, där de magnetiska egenskaperna från BOLD-signalen översätts till numeriska värden och amplituden i signalen kan mätas. En upplösning väljs vid chiffreringen vilket resulterar i ett bestämt antal voxlar som representerar utformningen av hjärnan och dess aktivitet. Voxeln är en tredimensionell motsvarighet till en pixel, som har dimensionerna x-, y- och z-led.

I figur 2.1 visualiseras ett exempel på en magnetresonanstomografi som gjorts om till analyserbar fMRI-data. Exemplets datamaterial är på 2mm och erhålls från MNITemplate (Fonov m. fl., 2011; Fonov, Evans, McKinstry, Almli och Collins, 2009; Collins, Zijdenbos, Baaré och Evans, 1999).

t z - Transversell x - Koronal y - Sagittal

Figur 2.1: Exempel på fMRI-data, framtaget i R

Första raden i figur 2.1 visualiserar data som kommer direkt från en fMRI-mätning, här är kraniet kvar innan en databerbetning har genomgåtts för att isolera hjärnan. I andra raden visualiseras individiens hjärna som nu består av voxlar som är 2mm3 _{stora, där varje voxel innehåller en tidpunkt t som består av}

fMRI-mätningens BOLD-signal. Till sist visualiseras ett exempel på tvärsnitt i x−, y− samt z-led som representerar det koronala, sagittala respektive det transversella planet i hjärnan.

(15)

Datamaterialet som analyseras i uppsatsen består av 20 stycken individers fMRI-mätningar samt tillhöran-de gruppvisa oberoentillhöran-de komponenter för tillhöran-de 20 indivitillhöran-derna. Detta har tillhandahållits av Antillhöran-ders Eklund, som i sin tur har hämtat datamaterialet Cambridge-Buckner vilket är en del i 1000 Connectome Project. Cambridge-Buckner datamaterialet består ursprungligen av 198 individer varav 75 är män och 123 är kvinnor i åldrar mellan 18 till 30 år (Biswal m. fl., 2010), data innehåller inga saknade värden. Datamaterialet har även bearbetats av Bertil Wegmann. Varje individs fMRI-data består av 4 dimensioner, där tre av dimensio-nerna skapar tredimensionella pixlar som kallas voxlar. Dessa dimensioner består av x, y och z samt den 4e dimension består av tillhörande tidserie t för varje voxel bestående av 119 stycken mätta tidpunkter. Varje voxel är en 6 mm3 _{kub och tidserien är skapad genom att ta bilder med ett tidsintervall på 3 sekunder där}

varje observation utger en BOLD-signal.

Mjukvaran FMRIB Software library (FSL) är ett program som innehåller verktyg för att analysera bland annat fMRI-data. FSL har använts av Anders Eklund för att ta fram Group Average Spatial Independent Component Analysis (GICA) som används i uppsatsen. Totalt togs 177 komponenter fram i en så kallad GICA. Av alla 177 komponenter kommer de 20 komponenterna som skattats att användas. De 20 komponen-terna som valts är baserat på enfunktion i FSL som ordnar komponenter efter de komponenter som förklarar mest unik varians.

För tolkning av aktiviteten i hjärnans nätverk kommer hjärnan delas i fem olika regioner: Pannloben, hjäss-loben, nackhjäss-loben, tinningloben samt lillhjärnan. Dessa regioner delas upp ungefärligt enligt figur 2.2.

Figur 2.2: Uppdelning av hjärnregioner

När individspecifika nätverkskartor tolkas i uppsatsen kommer detta göras på en övergripande nivå som matchar författarnas kunskapsnivå om hjärnans anatomi.

(16)

I figur 2.3 nedan visualiseras 20 komponenter av GICA där tvärsnitt från det koronala, sagittala samt det transversella planet. I uppsatsen kommer dessa 20 komponenter användas som förklarande variabler i respek-tive multivariat regressionsmodell för den temporala och spatiala delen.

Figur 2.3: GICA 1-20:Färgskalan är baserad på den centrerade amplituden i GICA där positiva värden visualiseras från 0-25. Amplitudskalan går från röd till blå färg beroende på aktivitet i hjärnan där röd indikerar mindre aktivitet och gul till blå högre aktivitet. Aktivitet under 0 analyseras inte visuellt i uppsatsen.

I figur 2.3 visualiseras de 20 komponenterna som beräknats för de 20 individerna radvis. De flesta komponenter tycks visa minst ett område i hjärnan där högre amplitud förekommer i gul eller blå färg som avviker från resterande nätverk. Detta indikerar högre aktivering inom dessa områden.

2.1 Centrering av data

Innan modellanpassningen för individspecifika nätverkskartor kan påbörjas kräver data en centrering för både tidsaspekten och den spatiala aspekten samt en skalning med hänsyn för varje tidpunkt innan Dual Regres-sionen kan beräknas. Detta genomförs för att en korrekt inferens ska kunna göras och för att ta hänsyn till de olika amplituder och rörelsebrus i fMRI-data. Centreringen och skalningen genomförs även på grund av att individer tenderar att ha olika amplituder för olika tidsförlopp samt för att kunna jämföra individer sinse-mellan. Genom denna centrering kräver modellanpassningen inget intercept och multivariat linjär regression kan därmed anpassas enklare och mer beräkningseffektivt (Mejia m. fl., 2019).

I ekvation 2.1 centreras data Y voxelvis för att sedan i ekvation 2.2 centreras data C1 för varje tidpunkt. N betecknar voxlar, T för tidsförlopp samt M komponenter.

C1N xT = YN xT − ¯YN (2.1)

C2N xT = C1N xT − ¯C1T (2.2)

Efter de två centreringarna är utförda skapas slutligen en skalning enligt:

C3N xT = A−1C2N xT (2.3)

Där en diagonalmatris A = diag {sd(YT)} där diagonalelemeten består av standardavvikelsen från varje

tidpunk T från Y som erhålls innan centreringen av C1 eller C2 genomförts. En liknande centrering och

(17)

3. Metod

I metodkapitlet ges en beskrivning av de statistiska metoderna samt utförande av dessa, med referenser till artiklar, som används i uppsatsen.

3.1 Matematiska operatörer

Matematiska operatörer presenteras som möjliggör beräkningar i matrisform: vektorisering, kroneckerpro-dukten samt spår (Henderson & Searle, 1981).

Vektorisering vec(·) antar en matris n × p och skapar en vektor med längden np enligt:

A =

a b c d

, som vid vektorisering skapar vec(A) =     a c b d     (3.1)

Kroneckerprodukten ⊗ erhåller varje möjliga kombination av produkten från elementen av två matriser. Vid A ⊗ Bdär A är m × n och B är p × q erhålls en matris av mp × nq enligt ekvation 3.2. Denna matematiska operation möjliggör multiplikation av matriser i olika dimensioner enligt:

A ⊗ B =    a11B · · · a1nB ... ... ... am1B · · · amnB    (3.2)

Spåret tr(·) antar en symmetrisk matris A av n × n där diagonalelementen summeras enligt:

tr(A) =

n

X

i=1

(18)

Den matematiska operatören pinv(·) som betecknar Moore–Penroses pseudoinvers möjliggör en generalisering för matrisinvers som kan appliceras på icke-kvadratiska matriser där pinv(A) = A+_{(Penrose, 1955). För att}

beräkna Moore-Penrose pseudoinvers för en matris A av storlek m × n används singulärvärdesuppdelning.

A = UΣV∗ (3.4)

A+= VΣ+U∗ (3.5)

Där Σ+_{skapas från Σ genom att ta inversen av varje singulärvärde som inte är noll, de singulärvärde är noll}

i Σ sätts även till noll i Σ+ _{som sedan transponeras. Där ∗ betecknar det Hermiteska konjugatet. U är av}

storleken m × m och V är av storleken n × n dessa matriser är unitära (Ben-Israel & Greville, 2003). Moo-re–Penrose pseudoinvers kan såldes användas för att lösa minsta kvadratmetoden enligt ˆβββ = A+_{Y = YA}+_.

3.2 Multivariat linjär regression

För att undersöka linjära samband mellan flera responsvariabler och förklarande variabler kan multivariat regression användas. I multivariat regression är responsvariabeln en matris enligt:

     y1,1 y1,2 . . . y1,p y2,1 y2,2 . . . y2,p ... ... ... ... yn,1 yn,2 . . . yn,p      =      x1,1 x1,2 . . . x1,k x2,1 x2,2 . . . x2,k ... ... ... ... xn,1 xn,2 . . . xn,k      ·      β1,1 β1,2 . . . β1,p β2,1 β2,2 . . . β2,p ... ... ... β3,p βk,1 βk,1 . . . βk,p      +      1,1 1,2 . . . 1,p 2,1 2,2 . . . 2,p ... ... ... 3,p n,1 n,2 . . . n,p      (3.6) I ekvation 3.6 ovan observeras matrisen för responsvariablerna Y vara en matris av n rader och p kolumner. Matrisen för förklaringsvariablerna X ∈ RN xK_{. Matrisen med regressionskoefficienterna βββ har dimensionerna}

RKxP_{. Modellens feltermer är en matris av storlek R}N xP_.

Ekvation 3.6 kan även uttryckas i matrisform:

Y = Xβββ + (3.7)

Likt antagandet i univariat regression där feltermerna ska vara normalfördelade med ett väntevärde som är 0, följer en naturlig påbyggnad i multivariat regression att matrisen för feltermerna ska vara multivariat normalfördelad med väntevärde 0 enligt:

iid∼ N (0, σ2_{· I)} _(3.8)

Modellen antar även att varje observation n i Y har samma kovariansmatris, samt att de N raderna i Y är okorrelerade med varandra. Det får existera ett beroende inom kolumnerna p för observationer i Y (Rencher & Schimek, 2003).

(19)

3.3 Dual Regression

Dual Regression används för att undersöka den funktionella konnektiviteten i hjärnan från fMRI-data där individspecifika nätverkskartor kan identifieras (Nickerson m. fl., 2017). Från varje individs fMRI-data har en gemensam GICA tagits fram som används för att identifiera de motsvarande nätverkskartorna för varje individ. GICA som betecknas ˆS bildar de förklarande variablerna för det första steget i Dual Regressionen och struktureras om till en matris (N voxlar × M komponenter). Från de 4 dimensionerna i fMRI-data skapas Y som har strukturerats om till en matris (N voxlar × T tidpunkter). T C betecknar tidsförlopp och SM betecknar spatiala nätverk.

Dual Regression är en metod som kan användas för att göra spatio-temporal regression. Regressionen har två steg, där det första steget är en temporal regression enligt:

Y

YY = ˆSβββT C + T C (3.9)

T C iid∼ NN T(0, ΣΣΣ) (3.10)

I regressionens första steg i ekvation 3.9 antar matriserna dimensionerna Y ∈ RN xT_{, ˆS ∈ R}N xM_,

β β

βT C_{∈ R}M xT_,T C _{∈ R}N xT _{samt i ekvation 3.10 har kovariansmatrisen dimensionerna Σ}_Σ_{Σ ∈ R}T xT_.

Det andra steget är en spatial regression:

YYY0= βββ0T CβββSM + SM (3.11)

SM iid∼ NT N(0, ΣΣΣ) (3.12)

I regressionens andra steg används matrisen βββT C _{från regressionens första steg som designmatris.}

Dimen-sionerna för matriserna i regressionens andra steg är Y0 _{∈ R}T xN_{, β}_β_β0T C _{∈ R}T xM_{, β}_β_βSM _{∈ R}M xN _samt

SM _{∈ R}T xN_{. Från ekvation 3.12 har kovariansmatrisen dimensionerna Σ}_Σ_{Σ ∈ R}N xN_{. Eftersom}

modellan-passning genom Dual Regression kan skattas både genom frekventistisk och Bayesiansk inferens kommer varians och kovariansstrukturen för variera mellan olika modeller.

3.4 Frekventistisk Dual Regression

I detta avsnitt anpassas Dual Regressionen på ett frekventistiskt tillvägagångssätt enligt Nickerson m.fl. 2017.

3.4.1 Regressionssteg för individspecifika tidsförlopp

I det första steget av regressionen för Dual Regression identifieras de individspecifika tidsförloppen för de gemensamma spatiala kartorna.

Y = ˆSβββT C+ T C (3.13) ˆ βˆ βˆβT C = pinv(ˆS)Y (3.14) ˆ βˆ βˆ

βT C erhålls genom multivariat OLS, där pinv(·) betecknar Moore–Penrose pseudoinvers. Y ∈ RN xT _är

in-dividens fMRI-data, ˆS ∈ RN xM _{är de gemensamma spatiala kartorna erhållna från GICA,}T C _{∈ R}N xT

är matrisen för feltermerna och ˆβˆβˆβT C ∈ RM xT _{blir således det genomsnittliga tidsförloppet beräknat över}

(20)

3.4.2 Regressionssteg för individspecifika nätverkskartor

I den andra regressionen för dual regression identifieras de individspecifika nätverkskartor. Detta görs med hjälp av ˆβˆβˆβT C _{från regressionens första steg enligt:}

Y0= ˆβββ0T CβββSM + SM (3.15)

ˆ β β

βSM = Y · pinvˆβββT C (3.16)

Från ˆβˆβˆβSM _{∈ R}N xM _{erhålls således de spatiala kartorna som motsvarar de oberoende komponenterna för en}

enskild individ. Feltermerna T C _ochSM _{antas vara multivarat normalfördelade med väntevärde 0 samt}

oberoende temporalt och spatialt. Detta resulterar i att kovariansmatriserna från ekvation 3.10 samt 3.12 antar strukturen ΣΣΣ = σ2_{· I}_{, vilket gör att modellen antar homoskedastisk varians och okorrelerade feltermer.}

3.5 Sannolikhetsfördelningar

I uppsatsen kommer särskilda sannolikhetsfördelningar antas i den Bayesianska inferensen för att paramet-rarnas posteriorfördelningar ska följa kända sannolikhetsfördelningar. Slumpgenererade tal från dessa kända sannolikhetsfördelningar kommer senare att dras för anpassning av den Bayesianska multivariata regressionen. Genom uppsatsens tidsbegränsning antas dessa följa direkt kända posteriorfördelningar för att kunna genere-ra snabba slumpmässiga dgenere-ragningar som är lösbagenere-ra analytiskt, istället för att använda metoder som MCMC där enskilda dragningarna tar längre tid att göra samt inte alltid resulterar i effektiva posteriordragningar (Kruschke, 2014).

3.5.1 Invers Wishartfördelning

Den inversa Wishartfördelningen är den naturligt konjugerande priorfördelningen för kovariansmatrisen bland variabler som följer en matrisnormalfördelning. Fördelningen har två parametrar: matrisen V som är en skalparameter samt ν som anger antalet frihetsgrader. Den inversa Wishartfördelningens täthetsfunktion är:

fx(X; V, ν) = |V|ν/2 2νp/2_Γ p ν₂ |X| −(ν+p+1)/2_e−1 2tr(VX −1₎ (3.17) Där V betecknar matrisen som är en skalparameter med dimensionerna RP xP_{. Matrisen ν betecknar antalet}

frihetsgrader samt Γpär en multivariat gamma funktion och tr(·) betecknar matrisens spår (Gupta & Nagar,

2018).

Sannolikhetsfördelningens väntevärde beräknas genom:

E[X] = V

ν − p − 1 (3.18)

I ekvation 3.18 observeras det att sannolikhetsfördelningen är restrikterad till att antalet frihetsgrader ν i fördelningen måste vara större än p + 1 för att ha ett definierat medelvärde, där p betecknar antalet parametrar. Det kan genom sannolikhetsfördelningens ekvationer för kovarians och varians observeras att när antalet frihetsgrader ν ökar blir priorfördelningen mer informativ:

V ar(Xn,k) =

(ν − p + 1)V2

n,k+ (ν − p − 1)Vn,nVk,k

(21)

Samt för variansen längs matrisens diagonal: V ar(Xn,n) =

2V2 n,n

(ν − p − 1)2_{(ν − p − 3)} (3.20)

Där det observeras att med ökande antal frihetsgrader ν i ekvationernas nämnare kommer variansen att minska.

3.5.2 Matrisnormalfördelning

En stokastisk variabel X ∈ RN xP _{som följer matrisnormalfördelning är det generella fallet av en multivariat}

normalfördelning som tillämpas för matriser där kovariansmatrisen är separerbar för till exempel temporal och spatial data (Rowe, 2002).

Matrisnormal är relaterad till multivariat normalfördelad enligt ekvationen nedan:

X ∼ MNn×p(βββ, ΛΛΛ, Σ)där (3.21)

vec(X) ∼ Nnp(vec(βββ), Σ⊗ ΛΛΛ)

Där βββ ∈ RN xP _{anger läget för alla parametrar, Λ}_Λ_{Λ ∈ R}N xN _{ger den första kovariansmatrisen för rader och}

Σ∈ RP xP ger den andra kovariansmatrisen för kolumner. Här används två stycken matematiska operatörer,

vektorisering och vec(·) enligt ekvation 3.1 samt kroneckerprodukten ⊗ enligt ekvation 3.2. Ett antagande för matrisnormalfördelning är om kroneckerprodukten är separerbar för kovarianserna så att X ∼ MN (βββ,ΛΛΛ, Σ)

antas vara den samma som vec(X) ∼ Nnp(vec(β), Σ⊗ Λ)och således kan matrisnormal tillämpas.

Täthetsfunktionen för matrisnormalfördelning är: p(X|βββ, ΛΛΛ, Σ) = exp(1₂tr[Σ−1 (X − βββ)0ΛΛΛ −1 (X − βββ]) (2π)np/2_|Σ |n/2|ΛΛΛ|p/2 (3.22)

Matrisnormalfördelningens väntevärde blir således:

E(X|βββ, ΛΛΛ, Σ) = βββ (3.23)

Och fördelningens varians:

V ar (vec (X0) |βββ, ΛΛΛ, Σ) = ΛΛΛ ⊗ Σ (3.24)

3.6 Bayes sats

För att uttrycka sannolikheten av ett utfall, givet förkunskaper kan Bayes sats användas. Bayes sats kan naturligt appliceras för användning av Bayesianska modeller:

p(θ | x1, ..., xn) =

p(x1, ..., xn) | θ) · p(θ)

p(x1, ..., xn) (3.25)

I 3.25 ovan betecknas p(θ) som priorfördelningen för parametern θ, p(x1, ..., xn)som sannolikheten för data

enligt modellen, denna normaliserar uttrycket i nämnaren vilket leder till att resultatet blir en sannolik-het.(Kruschke, 2014)

(22)

prior-p(θ | x1, ..., xn) ∝ p(x1, ..., xn) | θ) · p(θ) (3.26)

Termen p(x1, ..., xn | θ) är likelihoodfunktionen för data givet parametern θ.

Posteriorsannolikhetsfördel-ningen p(θ | x1, ..., xn)kan beräknas genom att multiplicera likelihoodfunktionen med priorfördelningen och

normalisera uttrycket med hjälp av sannolikheten för data enligt modellen (Kruschke, 2014).

3.7 Bayesiansk multivariat linjär regression

I den Bayesianska multivariata regressionen antas alla parametrar vara slumpmässiga enligt avsnitt 3.6. Det följer då att dessa parametrar följer sannolikhetsfördelningar, från ekvation 3.6 presenterat tidigare är βββ samt slumpvariabler där varje observation n i antas n vara multivariat normalfördelad med en positiv definit

kovariansmatris ΣΣΣ. Matrisen βββ antas ha följande naturlig konjugerande priorfördelning:

p(βββ, Σ) = p(Σ) · p(βββ|Σ)där (3.27)

p(Σ) ∼ W−1(V0, ν0)och

p(βββ|Σ) ∼ N (βββ0, Σ⊗ Λ−10 )

I ekvationen 3.27 ovan betecknar W−1_{den inversa Wishartfördelningen och kroneckerprodukten ⊗ för de två}

matriserna Σoch Λ−10.

Detta resulterar i följande posteriorfördelningar för parametrarna p(ΣΣΣ|Y, X)blir:

p(ΣΣΣ|Y, X) ∼ W−1(Vn, νn)där (3.28)

V

VVn= VVV0+ (Y − Xβββn)0(Y − Xβββn) + (βββn− βββ0)0ΛΛΛ0(βββn− βββ0)

νn= ν0+ n

Samt posteriorfördelningen för p(βββ|Y, X,ΣΣΣ)blir:

p(βββ|Y, X, ΣΣΣ) ∼ M Nk,p(βββn, ΛΛΛ−1n , ΣΣΣ)där (3.29)

βββn = (X0X + ΛΛΛ0)−1(X0Y + Λ0βββ0)

ΛΛΛn = X0X + Λ0

Där MNk,p betecknar en matrisnormalfördelning med parametrarna βββn, vilket är en matris som anger

för-delningens position, samt matriserna ΛΛΛ−1_n och ΣΣΣ som fungerar som fördelningens skalparametrar (Rossi,

(23)

3.8 Bayesiansk Dual Regression

Likt de generella ekvationer 3.9 samt 3.11 i avsnitt 3.3 skattas den Bayesianska Dual Regressionen i två steg. Det första steget är för de individspecifika tidsförloppen samt det andra för individspecifika spatiala nätverkskartor. I denna Dual Regression har kovariansmatriserna ΣΣΣför T C samt SM från ekvationerna 3.10 respektive 3.12 heteroskedastisk varians och tar hänsyn till korrelation bland feltermerna.

3.8.1 Regressionssteg för individspecifika tidsförlopp

Den första Bayesianska Dual Regressionsmodellen som kommer anpassas är en heteroskedastisk modell som tar hänsyn till korrelerade feltermer. I den Bayesianska Dual Regressionen dras slumpmässiga värden baserat på posteriorparametrarna presenterat i avsnitt 3.7.

I regressionens första steg:

Y = ˆSβββT C+

T C (3.30)

Anpassas den Bayesianska motsvarigheten av βββT C _{som presenterades i ekvation 3.9, denna beräknas nu genom}

tillvägagångssättet i ekvation 3.29:

βββT C_n = (ˆS0S + Λˆ ΛΛT C₀ )−1(ˆS0Y + ΛΛΛT C₀ βββT C₀ ) Där ΛΛΛT C

0 är en kvadratisk identitetsmatris som anger läget för kovariansen samt βββT C0 är en lägesmatris för

medelvärdet för βββ-parametrarna.

Posteriorfördelningen för hyperparametern ΛΛΛT C_n kan sedan beräknas genom: Λ

Λ

ΛT C_n = ˆS0S + Λˆ T C₀ Hyperparametern VVVT C

n beräknas enligt följande:

V

VVT C_n = VVVT C₀ + (Y − ˆSβββT C_n )0(Y − ˆSβββT C_n ) + (βββT C_n − βββT C₀ )0ΛΛΛT C₀ (βββT C_n − βββT C₀ ) Där VVVT C

0 är en kvadratisk identitetsmatris.

Hyperparametern νT C

n anger frihetsgraderna för inversa Wishartfördelningen enligt:

ν_nT C = ν₀T C+ nT C nT C _{anger antalet rader från Y, därefter sätts ν}T C

0 beroende på hur informativ fördelning som eftersträvas

givet att νn > p − 1. När hyperparametrarna är beräknade dras det slumpmässigt observationer från den

inversa Wishartfördelningen som sedan används till att generera observationer för posteriorfördelningen av β

β

βT C _{från matrisnormalfördelningen enligt ekvationerna 3.28 samt 3.29:}

p(ΣΣΣT C |Y, ˆS) ∼ W−1(V_nT C, ν_nT C) (3.31)

p(βββT C|Y, ˆS, ΣΣΣT C ) ∼ M Nk,p(βββT Cn , ΛΛΛ−1n

T C_{, Σ}_Σ_ΣT C

) (3.32)

Genom att generera slumpmässiga observationer från posteriorfördelningarna ovan kan de individspecifika tidsförloppen erhållas och användas till vidare användning för att skatta de individspecifika nätverkskartorna i regressionens andra steg.

(24)

3.8.2 Regressionssteg för individspecifika nätverkskartor

Samma process genomförs sedan för regressionens andra steg, där βββ0_{T C}

används som designmatris:

Y0= βββ0T CβββSM +

SM (3.33)

För βββSM

n används Moore–Penrose Psuedoinvers för att hantera skattningarna då antalet kolumner överstiger

antalet rader. Om icke-informativa priorfördelningar används där βββ0 = 0 samt ΛΛΛ0 är en nollmatris blir

ekvation 3.29 följande: β ββSMn = pinv(βββ 0_{T C} )Y Posteriorfördelningen för hyperparametern ΛΛΛSM

n kan sedan beräknas likt 3.29:

Λ Λ

ΛSM_n = βββT Cβββ0T C+ ΛSM₀

Där ΛΛΛSM0 är en kvadratisk identitetsmatris som anger läget för kovariansen.

Posteriorfördelningen för VVVSM

n beräknas likt ekvation 3.28:

V V VSM_n = VVVSM₀ + (Y − βββ0T CβββSM_n )0(Y − βββ0T Cβββ_nSM) + (βββSM_n − βββSM₀ )0ΛΛΛSM₀ (βββSM_n − βββSM₀ ) β β βSM

0 är en lägesmatris för medelvärdet för βββ-parametrarna och VVVSM0 är en kvadratisk identitetsmatris.

När dessa hyperparametrar har beräknats för andra steget kan posteriordragningar generera de individspeci-fika nätverkskartorna. Detta görs genom att göra slumpässiga dragningar från följande posteriorfördelnignar: p(ΣΣΣSM |Y, ˆS) ∼ W−1(VSM_n , ν_nSM) (3.34) p(βββ0SM|Y0, βββ0T C, ΣΣΣ) ∼ M Nk,p(βββSMn , ΛΛΛ −1 n SM_{, Σ}_Σ_ΣSM ) (3.35)

3.9 Bayesiansk Dual Regression utan korrelerade feltermer

Den andra Bayesianska Dual Regressionsmodellen som kommer anpassas är en heteroskedastisk modell där det antas att det inte finns någon korrelation mellan modellens feltermer. Jämfört mot den föregående poste-riorfördelningen för βββSM _{ekvationerna 3.32 samt 3.35 kommer den inte ha en Σ}_Σ_Σ

som följer en invers

Wishart-fördelning. Istället används en identitetsmatris med betaparametrarnas varians som diagonalelement, vilket resulterar i följande modell:

p(βββT C|Y, ˆS, σσσ2· I) ∼ M Nk,p(βββT Cn , ΛΛΛ −1 n T C_{, σ}_σ_σ2_{· I)} _(3.36) p(βββ0SM|Y, βββT C, σσσ2· I) ∼ M Nk,p(βββSMn , ΛΛΛ−1n SM_{, σ}_σ_σ2_{· I)} _(3.37)

Där matrisen med varianserna för ˆβˆβˆβT C _{skattas genom minsta kvadratmetoden:}

S(βββT C) = (Y − ˆSβββT C)0(Y − ˆSβββT C) (3.38) σσσ2_T =n − p n S(βββT C₎ n − p = 1 nS(βββ T C₎ _(3.39)

(25)

I ekvation 3.39 ovan betecknar n antalet rader och p antalet kolumner enligt ekvation 3.6. Samma process upprepas sedan för βββSM _{och σσσ}2

N, men istället för ˆS används βββ0T C likt den föregående Bayesianska Dual

Regressionsanpassningen. Kovariansmatrisen ΣΣΣi T C _{från ekvation 3.10 i regressionens första steg antar nu}

formen σσσ2

T· I, samt för kovariansmatrisen ΣΣΣfrån ekvation 3.12 i regressionens andra steg antas formen σσσ 2 N· I.

Detta resulterar i en modell med okorrelerade heteroskedastiska feltermer.

3.10 Skillnader mellan frekventistisk och Bayesiansk Dual

Regres-sion

Skillnaden mellan de två tillvägagångssätten är att den Bayesianska inferensen i avsnitt 3.8 inte kräver att feltermerna T C _{från ekvation 3.30 samt}SM _{från ekvation 3.33 ska vara okorrelerade, vilket är ett antagande}

i den frekventistiska modellen. Eftersom dessa kan tillåtas att vara korrelerade kan det antas mer komplexa strukturer på feltermerna som kan ta till hänsyn om hur information kan korrelera temporalt och spatialt. I denna uppsats kommer även den frekventistiska modellen ha homoskedastisk varians likt anpassningen i Nickerson m. fl., 2017. Detta innebär att båda Bayesianska modellerna tar hänsyn till om det finns skillnader i mängden variation mellan olika nätverkskartor, vilket den frekventistiska anpassningen inte gör.

3.11 Implementering i R

Vid anpassning av modellerna används funktioner från två olika paket i R: pracma (Borchers & Borchers, 2019) och LaplacesDemon (Hall, 2011). Från pracma används funktionen pinv() för att lösa psuedoinversen av matriserna för βββ-matriserna. Från paketet LaplacesDemon används funktionen as.symmetric.matrix() för att göra den kvadratiska matrisen i inversa Wishartfördelningens dragningar symmetrisk. Detta är ett problem som uppstår genom begränsningar i hur R lagrar decimaltal. Eftersom R använder IEEE 754 double-precision floating point (R Development Core Team m. fl., 2011) som lagringsformat är detta endast precist till 53 bits av information, vilket i decimalform blir 2−53 _{≈ 1.11 · 10}−16_{. Decimaler efter detta behandlas som brus och}

används inte precist vid beräkningar. Eftersom dragningar med hjälp av funktioner från LaplacesDemon krä-ver att matrisen ska vara exakt symmetrisk måste matrisen göras symmetrisk genom as.symmetric.matrix(). Dragningarna från den inversa Wishartfördelningen samt matrisnormalfördelningen görs sedan genom pake-tets funktioner rinvwishart() samt rmatrixnorm().

(26)

4. Analys

I analysen kommer valet av priorfördelningens hyperparametrar presenteras. Den frekventistiska modellen kommer refereras till modell 1, den Bayesiansk heteroskedastiska modellen utan korrelerade feltermer kommer att refereras till modell 2 samt den Bayesianska heteroskedastiska modellen med korrelerade feltermer refereras i analysen som modell 3. Tre stycken individer kommer analyseras för att säkerställa att de Bayesianska anpassningarna går att jämföra individer sinsemellan.

4.1 Val av priorfördelningarnas hyperparameterar

Vid anpassningen av de Bayesianska modellerna kommer priorfördelningarna som används vara icke-informativa. Detta uppnås genom att använda så lite priorinformation i fördelningarnas hyperparametrar som möjligt så att posteriorfördelningen påverkas huvudsakligen av data. Därför kommer hyperparametern ν0= 0i

regres-sionens första steg i ekvation 3.30 vilket innebär att posteriorfördelningen för νn enligt ekvation 3.28 endast

beror på data. I det första steget ger detta att för individ 1 är νn = 8349, för individ 2 är νn= 7997samt för

individ 3 är νn= 8075, en skillnad mellan individernas hyperparametervärden förekommer på grund av den

naturliga skillnaden som finns mellan de olika individernas hjärnvolym. För andra steget av regressionen i ekvation 3.33 sätts ν0= 2för att väntevärdet för fördelningen ska vara definierat enligt ekvation 3.18, vilket

resulterar i att för individ 1 är νn = 8351, för individ 2 är νn= 7999samt νn = 8077för den tredje individen.

Valet av dessa konstanter kan ha en påverkan på inferensen eftersom priorfördelningen kan bli mer informativ, men detta är oundvikligt då det är den lägsta möjliga konstanten för att kunna göra posteriordragningar med ett definierat medelvärde från den inversa Wishartfördelningen.

Matriserna VVV0 samt ΛΛΛ0från avsnitt 3.8 är identitetsmatriser samt βββ0 är en nollmatris.

4.2 Utvärdering av posteriordragningar

För att utvärdera resultaten av posteriordragningarna kommer det visualiseras ett linjediagram över det ackumulerade posteriormedelvärdet för respektive nätverkskarta. Detta görs för både modell 2 och 3. Det-ta görs för samtliga individer. Det kommer även att visualiseras varje individs voxelvisa utveckling till det ackumulerade posteriormedelvärdet för modell 3. Denna nätverkskarta visar utvecklingen över de 100 poste-riordragningarna för alla voxlar där det ackumulerade posteriormedelvärdet är positivt. Visualiseringen visas endast i z-ledets tvärsnitt på grund av begränsningar i programvarans funktioner.

Det kommer även visualiseras ett histogram över posteriordragningarna som kommer att användas som en mall för visualiseringen av nätverken som skattats genom de Bayesianska modellerna. Voxlarna som uppfyller ett kriterie där sannolikheten över 50% för posteriordragningar har ett posteriorvärde på P (βββSM _{> 0)}

klas-sificeras som aktiva regioner och kommer att synas i bilderna över nätverkskartor, med en röd till blå färg baserat på voxelns amplitud.

(27)

4.2.1 Individ 1

I figur 4.1 visualiseras utvecklingen för det ackumulerade posteriormedelvärdet över de 100 dragningarna för nätverkskarta 1 till 6 som innehåller alla voxlar för βββSM _{för individ 1. Utvecklingen som visualiseras är för}

båda Bayesianska modellerna.

5 6 3 4 1 2 0 25 50 75 100 0 25 50 75 100 −0.01 0.00 0.01 −0.01 0.00 0.01 −0.01 0.00 0.01 Antal posteriordragningar Ack. poster ior medelv ärde Modell Modell 2 Modell 3

Figur 4.1: Utveckling av ackumulerade posteriormedelvärde:Det ackumulerade posteriormedelvärdet över de 100 dragningarna för nätverkskarta 1 till 6 som innehåller alla voxlar för βββSM_{. Individ 1, modell 2 och 3.}

Det observeras att modell 2 verkar stabiliseras direkt till det gemensamma posteriormedelvärdet för alla voxlar, men att modell 3 som tar till hänsyn korrelation mellan feltermer tar lite längre tid att stabiliseras. Båda modellernas ackumulerade posteriormedelvärde verkar dock vara stabiliserat efter cirka 50 dragningar. Eftersom varje voxel enskilt bör nå ett stabilt ackumulerat posteriormedelvärde undersöktes även detta. I analysen visas endast resultatet av den första nätverkskartan för respektive individ. I bilagan går det även att se hur utvecklingarna ser ut för nätverkskarta 4 och 6. Visualiseringen av utvecklingen gjordes för de voxlar där posteriormedelvärdet var positivt. Det har skapats en video över den voxelvisa utvecklingen över posteriordragningarna från de tre individerna samt de tre utvalda nätverkskartorna: (Jonsson och Welander, 2020). Från videon erhålls stöd för att alla voxlar i z-ledets tvärsnitt når ett stabilt ackumulerat posterior-medelvärde. I figur 4.2 visualiseras nu en linje för varje voxel i z-ledets tvärsnitt, för att undersöka hur den voxelvisa utvecklas till det ackumulerade posteriormedelvärdet ser ut.

Figur 4.2: Utveckligen av voxelvisa ackumulerade posteriormedelvärdet:Det ackumulerade posteriormedelvärdet i tvärsnitt för z-axeln visualiseras för alla voxlar som har ett positivt posteriormedelvärde. Individ 1, modell 3.

(28)

Det observeras i figur 4.2 för nätverkskarta 1 att voxlarna längs z-ledets tvärsnitt verkar utvecklas till ett stabilt ackumulerat posteriormedelvärdet för varje voxel över de 100 posteriordragningarna. Det resulterande nätverkets aktivitet tolkas enligt regionerna i figur 2.2. Det verkar förekomma högst aktivering i nackloben, men nätverket finns även utspritt i delar av hjässloben samt i mindre utsträckning i främre delen av pannloben. I ekvation 4.1 summeras antalet posteriordragningar p för varje voxel v i βββSM _{som har ett värde över 0 för}

att sedan ta fram sannolikheten att posteriordragningarna för varje voxel är över 0. P(βββSM> 0) = ( P100 p=1βββ SM vp > 0) 100 (4.1)

I histogrammet i figur 4.3 visualiseras antalet posteriordragningar för P(βββSM _{> 0)} _{för modell 2 och 3 för}

nätverkskarta 1 till 6. En stor spridning i detta histogram innebär att sannolikheten att värdet på en voxel är positivt eller negativt är ungefär lika för respektive modell.

Nätverkskarta: 5 Nätverkskarta: 6 Nätverkskarta: 3 Nätverkskarta: 4 Nätverkskarta: 1 Nätverkskarta: 2 0 25 50 75 100 0 25 50 75 100 0.0 0.1 0.2 0.3 0.0 0.1 0.2 0.3 0.0 0.1 0.2 0.3 Antal posteriordragningar Sannolikhet Modell Modell 2 Modell 3 Sannolikheten för att Bsm > 0 över alla voxlar per nätverkskarta

Figur 4.3: Histogram för P(βββSM_{> 0)}_:_{Nätverkskarta 1 till 6 över antalet posteriordragningar för P(β}_β_βSM_{> 0). En} stor spridning i detta histogram innebär att sannolikheten att värdet på en voxel är positivt eller negativt är ungefär lika för respektive modell. Individ 1, modell 2 och 3.

I figur 4.3 observeras för modell två att det är två stora spikar på respektive ände av fördelningen. Detta innebär att denna modell skattar voxlarna i stor utsträckning som entydigt positiva eller negativa. För modell 3 observeras en större spridning, vilket indikerar att modellen inte skattar voxlarna lika positivt entydiga eller negativa. För modell 3 som tar hänsyn till korrelerade feltermer observeras variationen för βββSM _{vara större.}

(29)

4.2.2 Individ 2

I figur 4.4 visualiseras utvecklingen av det ackumulerade posteriormedelvärdet för 6 stycken nätverkskartorna för individ 2, visualiseringen görs för både modell 2 och 3.

Figur 4.4: Utveckling av ackumulerade posteriormedelvärde:Det ackumulerade posteriormedelvärdet över de 100 dragningarna för nätverkskarta 1 till 6 som innehåller alla voxlar för βββSM_{. Individ 2, modell 2 och 3.}

I figur 4.4 syns det att det ackumulerade posteriormedelvärdet över dragningarna för modell 2 verkar stabilt för samtliga nätverkskartor. Samma sak observeras för modell 3 men det kräver lite fler dragningar för vissa av nätverkskartorna.

I figur 4.5 visualiseras den voxelvisa utvecklingen i z-ledets tvärsnitt till det ackumulerade posteriormedel-värdet för nätverkskarta 1 för individ 2 likt föregående individ.

Figur 4.5: Utveckligen av voxelvisa ackumulerade posteriormedelvärdet:Det ackumulerade posteriormedelvärdet i tvärsnitt för z-axeln visualiseras för alla voxlar som har ett positivt posteriormedelvärde. Individ 2, modell 3. Det observeras från figur 4.5 att den voxelvisa utvecklingen av det ackumulerade posteriormedelvärdet verkar ge ett stabilt posteriormedelvärde. Denna utveckling resulterar i ett nätverk med högst amplitud i nackloben, men det förekommer även aktivitet i hjässloben och pannloben.

(30)

I histogrammet i figur 4.6 visualiseras antalet posteriordragningar för P(βββSM _{> 0)} _{för modell 2 och 3 för}

nätverkskarta 1 till 6. Detta utförs likt tidigare för att undersöka hur sannolikt modellen skattar om hur många voxlar är positiva eller negativa.

Nätverkskarta: 5 Nätverkskarta: 6 Nätverkskarta: 3 Nätverkskarta: 4 Nätverkskarta: 1 Nätverkskarta: 2 0 25 50 75 100 0 25 50 75 100 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 Antal posteriordragningar Sannolikhet Modell Modell 2 Modell 3 Sannolikheten för att Bsm > 0 över alla voxlar per nätverkskarta

Figur 4.6: Histogram för P(βββSM_{> 0)}_:_{Nätverkskarta 1 till 6 över antalet posteriordragningar för P(β}_β_βSM_{> 0). En} stor spridning i detta histogram innebär att sannolikheten att värdet på en voxel är positivt eller negativt är ungefär lika för respektive modell. Individ 2, modell 2 och 3.

I figur 4.6 observeras det att posteriordragningarna i modell 2 antingen har ett värde kring 0 eller 100 för hur många voxlar som uppfyller kriteriet, detta indikerar på att modellen skattar med väldigt hög sannolikhet att en voxel är negativ eller positiv. Vid jämförelse med modell 3 är det inte längre lika entydig om voxlarna är negativa eller positiva.

4.2.3 Individ 3

I figur 4.7 visualiseras utvecklingen för det ackumulerade posteriormedelvärdet för varje nätverkskarta för individ 3 över de 100 posteriordragningarna. Utvecklingen som visualiseras är för Bayesianska modell 2 och 3.

Figur 4.7: Utveckling av ackumulerade posteriormedelvärde:Det ackumulerade posteriormedelvärdet över de 100 dragningarna för nätverkskarta 1 till 6 som innehåller alla voxlar för βββSM. Individ 3, modell 2 och 3.

(31)

Den tredje individens voxelvisa utveckling för det första nätverket för z-ledets tvärsnitt visualiseras likt tidigare i figur 4.8.

Figur 4.8: Utveckligen av voxelvisa ackumulerade posteriormedelvärdet:Det ackumulerade posteriormedelvärdet i tvärsnitt för z-axeln visualiseras för alla voxlar som har ett positivt posteriormedelvärde. Individ 3, modell 3. Det observeras i figur 4.8 att voxlarna verkar stabilisera till sitt ackumulerade posteriormedelvärde. Aktive-ringen i nätverket verkar högst i nackloben, men det förekommer även aktivering i pannloben.

Histogrammen i figur 4.9 visar visualiseras antalet posteriordragningar för P(βββSM_{> 0)}_{. Detta görs för att}

visualisera hur sannolikt modellen skattar om en voxel är negativ eller positiv. Denna jämförelse genomförs för varje nätverkskarta. Nätverkskarta: 5 Nätverkskarta: 6 Nätverkskarta: 3 Nätverkskarta: 4 Nätverkskarta: 1 Nätverkskarta: 2 0 25 50 75 100 0 25 50 75 100 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 Antal posteriordragningar Sannolikhet Modell Modell 2 Modell 3 Sannolikheten för att Bsm > 0 över alla voxlar per nätverkskarta

Figur 4.9: Histogram för P(βββSM_{> 0)}_:

Nätverkskarta 1 till 6 över antalet posteriordragningar för P(βββSM> 0). En stor spridning i detta histogram innebär att sannolikheten att värdet på en voxel är positivt eller negativt är ungefär lika för respektive modell. Individ 3, modell 2 och 3.

Likt tidigare observeras det i figur 4.9 att resultaten för den Bayesianska modell 2 att varje voxel skattas antingen negativt eller positivt med en hög sannolikhet, vilket observeras för varje nätverkskarta. För modell 3 kan en motsatt observation göras, där det är väldigt osannolikt att voxlarnas posteriordragning är uteslutande positiv eller negativa, för denna modell är resultatet mer normalfördelat med ett medelvärde runt 50.

(32)

4.3 Nätverkskartor

För varje individ prenteras nätverkskartorna 1, 4 och 6. Figurer för övriga nätverkskartor går att se i bilaga. Det sker även jämförelser mellan hur de tre olika modellerna skattar nätverken samt nätverkens amplitud. I bilderna över nätverkskartorna presenteras dessa i tvärsnitt från tre olika led, där de visualiseras i ordningen x, y och z.

I ekvation 4.2 bestäms färgtröskeln för de Bayesianska modellerna i uppsatsen. Färgtröskeln baseras på sannolikheten för att posteriordragningar för varje voxel är över 0 enligt ekvation 4.1. För den frekventistiska modell 1 visualiseras endast värden på βββSM _{över 0. Detta görs eftersom denna modell är centrerad runt 0}

och visualiseringen blir därför baserat på voxlar som är mer aktiva än medelvärdet. VoxelFärg(voxel) =(Posteriormedelvärde om P(βββSM> 0) > 0.5

Transparent om P(βββSM_{> 0) < 0.5} (4.2)

4.3.1 Individ 1

För individ 1 så visualiseras de tre modellernas olika skattningar på nätverkskarta 1 i figur 4.10.

Figur 4.10: Nätverkskarta 1:Modell 1, 2 samt 3 visas radvis för individ 1. Färgtröskeln baseras på sannolikheten för att posteriordragningar är över 50% för de Bayesianska modellerna. Den frekventistiska modellen visar endast positiva amplituder. En hög amplitud i färgskalan indikerar en högre aktivitet i området.

I figur 4.10 tycks en hög aktivitet ske i nackloben för alla modeller, resterande nätverk är spridda utöver hjärnan med mindre aktivitet. För modell 1 och 2 tycks en skattning av hög aktivering ske i nackloben men dess amplitud är lägre i modell 3. Denna modell har inte heller lika utspridda nätverk i pannloben som modell 1 och 2. Annars observeras ingen större skillnad i nätverkets utformning.

(33)

I figur 4.11 ger alla modeller en hög aktivering för båda sidor av tinningloben där den vänstra sidan uppger en högre amplitud än den högra. Här tycks modell 3 reducera lite brus och förstärka de områden där den höga aktivering sker jämfört mot modell 1 och 2.

För de tre modellerna visualiseras nätverkskarta 6 i figur 4.12.

Det observeras från figur 4.12 hög amplitud i pannloben och i lillhjärnan hos samtliga modeller. En del aktivitet förekommer i individens hjässlob och pannlob, för modell 3 är aktiviteten i denna region inte före-kommande i samma grad som för modell 2 och 3 vilket kan observeras genom tvärsnittet i y-led.

(34)

I figur 4.13 visualiseras en differens mellan Modell 3 som tar hänsyn till korrelerade fel och Modell 2 som ej tar hänsyn till korrelerade feltermer. För individ 1 visas nätverkskarta 1,4 och 6. Detta genomförs för att se hur skillnaderna är mellan modellerna där det enda som skiljer modellerna åt är om korrelerade feltermer tas till hänsyn.

Figur 4.13: Individ 1: Differensen i aktivering mellan Bayesianska modellernas nätverkskartor 1, 4 och 6 Från figur 4.13 tycks differenserna mellan modellerna vara runt 0 för majoriteten av hjärnans områden. För tvärsnittet i y-led för nätverkskarta 1 kan en större negativ differens iakttas i individens hjässlob och nacklob, detta indikerar på att modell 3 har skattat ett värde i området som är lägre än modell 2. Även en större positiv differens på aktiveringen kan identifieras centralt i hjärnan för nätverkskarta 1 där modell 3 skattar området med mindre värde på aktivering jämfört med modell 2.

I figur 4.14 visualiseras standardavvikelsen för varje voxels posteriordragningar för modell 2 och modell 3. Detta genomförs för att få en övergripande bild om det finns voxlar och områden i hjärnan som har en hög standardavvikelse för att identifiera områden med stor variation i vardera modell.

Figur 4.14: Standardavvikelse för nätverkskarta 1, 4 och 6, individ 1:Standardavvikelsen för varje voxels poste-riordragningar för modell 2 och 3.

Från figur 4.14 varierar de olika modellernas skala på standardavvikelse där modell 2 går från 0-4 och modell 3 0-80, vilket innebär på att modell 3 har högre variation i sin anpassning än modell 2. Båda modeller tycks ha större variation i posteriordragningarna för voxlar som befinner sig centralt i y-led. Utöver de gemensamma områden med hög variation kan även ett flertal områden för modell 3 i hjärnan iakttas så som för centrala

(35)

pannloben vilket kan ses i y-led för samtliga nätverkskartor.

I figur 4.15 visualiseras nätverkskarta 6 där en ny tröskelnivå används. Denna tröskel är baserad på en 67% sannolikhet för att en voxel har positiv ett värde. Detta visualiseras tillsammans med föregående tröskelnivå som ligger på 50% sannolikhet. Detta genomförs för att undersöka vilka voxlar som har ännu högre sannolikhet att ha konsekvent positiva värden.

Figur 4.15: Nätverkskarta 6 för modell 3, individ 1:Två olika tröskelnivåer sätts på P(βββSM_{> 0). Första raden} visualiserar 50% sannolikhet och andra raden visualiserar 67% sannolikhet

Från figur 4.15 kan en skillnad iakttas mellan nätverket med tröskelnivå 50% och 67% sannolikhet. Mycket av bruset reduceras och det nya nätverket med högre tröskel visar mindre och ett mer kompakt nätverk kring den höga aktiveringen. Med den nya tröskeln är det 2 gånger så sannolikt att en voxel är positiv än negativ.

4.3.2 Individ 2

Vidare jämförs skillnaden i nätverkskarta 1 mellan de olika modellerna för individ 2 i figur 4.16.

Mellan modell 1 och modell 2 märks ingen större skillnad mellan de skattade nätverkskartorna enligt figur 4.16, varken i amplitud eller i nätverkets utseende. Generellt sett observeras högsta aktiviteten i nackloben men det finns även aktivitet utspritt i hjässloben och pannloben. Modell 3 som tar hänsyn till korrelationen

(36)

med undantag för nackloben där modell 3 har högre amplitud. Nätverkskarta 4 för individ 2 har följande utseende i figur 4.17.

Nätverkskartan i figur 4.17 visar högst aktivitet i nackloben och pannloben. Här skiljer sig modell 3 inte lika tydligt från modell 1 och 2, den visar lite högre aktivitet i de fåtalet aktiva voxlarna i lillhjärnan men i övrigt är kartorna lika.

Nätverkskarta 6 för den andra individen ser ut enligt följande i figur 4.18.

I figur 4.18 observeras högst aktivitet i pannloben, men även aktiva områden i delar av hjässloben för samtliga modeller. Modell 1 och 2 skattar mer aktivitet utspridd i nackloben jämfört mot modell 3, men denna modell verkar skatta lite högre amplitud i nätverken jämfört mot modell 1 och 2.

(37)

För den andra individen visualiseras differensen mellan posteriormedelvärdena för modell 3 och modell 2 i figur 4.19 nedan:

Figur 4.19: Individ 2: Differensen i aktivering mellan Bayesianska modellernas nätverkskartor 1, 4 och 6 För nätverkskarta 1 observeras några områden, framförallt i pannloben och nackloben där modell 3 har hög-re aktivitet än modell 2. I nätverkskartorna 4 och 6 observeras liknande nätverk för diffehög-renserna. Modell 2 visar lite högre aktivering i ett område av hjässloben, på två områden centralt i bakre pannloben och främre hjässloben samt i nackloben. Det observeras högre aktivering för modell 3 i delar av lillhjärnan, samt på y-ledets tvärsnitt observeras två områden centralt i området mellan hjässloben och tinningloben.

I figur 4.20 visualiseras standardavvikelsen för varje voxels posteriordragningar för modell 2 och modell 3 likt som för föregående individ.

Figur 4.20: Standardavvikelse för nätverkskarta 1, 4 och 6, individ 2:Standardavvikelsen för varje voxels poste-riordragningar för modell 2 och 3.

Från figur 4.20 varierar de olika modellernas skala på standardavvikelse där modell 2 går från 0-4 och modell 3 från 30-80, vilket innebär på att modell 3 har högre variation i sin anpassning än modell 2. Modell 2 tycks inte uppvisa någon högre standardavvikelse jämfört mot modell 3 som har större variation på flertal ställen för nätverkskarta 4 samt 6.

(38)

I figur 4.21 visualiseras nätverkskarta 1 där en ny tröskelnivå presenteras, 67% sannolikhet för att en voxel har positiv värde jämfört mot föregående tröskelnivå som ligger på 50%.

Figur 4.21: Nätverkskarta 1 för modell 3, individ 2:Två olika tröskelnivåer sätts på P(βββSM> 0). Första raden visualiserar 50% sannolikhet och andra raden visualiserar 67% sannolikhet

I figur 4.21 tycks den högre tröskeln reducera brus och koncentrera nätverkskartan där mest aktivitet före-kommer.

4.3.3 Individ 3

I figur 4.22 kan en hög aktivitet identifieras i en stor del av den centrala nackloben, samt en del spridd aktivitet i pannloben och hjässloben. Modell 3 visar högre amplitud jämfört mot de två andra modellerna i nackloben. I övrigt verkar samtliga modellers nätverk relativt lika.

En jämförelse mellan frekventistisk och Bayesiansk Dual Regression : för nätverkskartor i hjärnan vid resting-state fMRI

Kandidatuppsats i Statistik